Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CexmcEventAction.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /*
27  * ============================================================================
28  *
29  * Filename: CexmcEventAction.cc
30  *
31  * Description: event action
32  *
33  * Version: 1.0
34  * Created: 27.10.2009 22:48:08
35  * Revision: none
36  * Compiler: gcc
37  *
38  * Author: Alexey Radkov (),
39  * Company: PNPI
40  *
41  * ============================================================================
42  */
43 
44 #include <cmath>
45 #ifdef CEXMC_USE_PERSISTENCY
46 #include <boost/archive/binary_oarchive.hpp>
47 #endif
48 #include <G4DigiManager.hh>
49 #include <G4Event.hh>
50 #include <G4Circle.hh>
51 #include <G4VisAttributes.hh>
52 #include <G4VisManager.hh>
53 #include <G4VTrajectory.hh>
54 #include <G4TrajectoryContainer.hh>
55 #include <G4Colour.hh>
56 #include "CexmcEventAction.hh"
57 #include "CexmcEventActionMessenger.hh"
58 #include "CexmcEventInfo.hh"
59 #include "CexmcEventSObject.hh"
60 #include "CexmcEventFastSObject.hh"
61 #include "CexmcTrackingAction.hh"
62 #include "CexmcChargeExchangeReconstructor.hh"
63 #include "CexmcRunManager.hh"
64 #include "CexmcHistoManager.hh"
65 #include "CexmcRun.hh"
66 #include "CexmcPhysicsManager.hh"
67 #include "CexmcProductionModel.hh"
68 #include "CexmcProductionModelData.hh"
69 #include "CexmcEnergyDepositDigitizer.hh"
70 #include "CexmcEnergyDepositStore.hh"
71 #include "CexmcTrackPointsDigitizer.hh"
72 #include "CexmcTrackPointsStore.hh"
73 #include "CexmcTrackPointInfo.hh"
74 #include "CexmcException.hh"
75 #include "CexmcCommon.hh"
76 
77 
78 namespace
79 {
80  G4double CexmcSmallCircleScreenSize( 5.0 );
81  G4double CexmcBigCircleScreenSize( 10.0 );
82  G4Colour CexmcTrackPointsMarkerColour( 0.0, 1.0, 0.4 );
83  G4Colour CexmcRecTrackPointsMarkerColour( 1.0, 0.4, 0.0 );
84 
85 
86  inline G4double CexmcGetKinEnergy( G4double momentumAmp, G4double mass )
87  {
88  return std::sqrt( momentumAmp * momentumAmp + mass * mass ) - mass;
89  }
90 }
91 
92 
93 CexmcEventAction::CexmcEventAction( CexmcPhysicsManager * physicsManager,
94  G4int verbose ) :
95  physicsManager( physicsManager ), reconstructor( NULL ), opKinEnergy( 0. ),
96  verbose( verbose ), verboseDraw( 4 ), messenger( NULL )
97 {
98  G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
99  digiManager->AddNewModule( new CexmcEnergyDepositDigitizer(
100  CexmcEDDigitizerName ) );
101  digiManager->AddNewModule( new CexmcTrackPointsDigitizer(
102  CexmcTPDigitizerName ) );
103  reconstructor = new CexmcChargeExchangeReconstructor(
104  physicsManager->GetProductionModel() );
105  messenger = new CexmcEventActionMessenger( this );
106 }
107 
108 
109 CexmcEventAction::~CexmcEventAction()
110 {
111  delete reconstructor;
112  delete messenger;
113 }
114 
115 
116 void CexmcEventAction::BeamParticleChangeHook( void )
117 {
118  reconstructor->SetupBeamParticle();
119 }
120 
121 
122 void CexmcEventAction::BeginOfEventAction( const G4Event * )
123 {
124  G4RunManager * runManager( G4RunManager::GetRunManager() );
125  CexmcTrackingAction * trackingAction
126  ( static_cast< CexmcTrackingAction * >(
127  const_cast< G4UserTrackingAction * >(
128  runManager->GetUserTrackingAction() ) ) );
129  trackingAction->BeginOfEventAction();
130 
131  physicsManager->ResetNumberOfTriggeredStudiedInteractions();
132 }
133 
134 
135 CexmcEnergyDepositStore * CexmcEventAction::MakeEnergyDepositStore(
136  const CexmcEnergyDepositDigitizer * digitizer )
137 {
138  G4double monitorED( digitizer->GetMonitorED() );
139  G4double vetoCounterEDLeft( digitizer->GetVetoCounterEDLeft() );
140  G4double vetoCounterEDRight( digitizer->GetVetoCounterEDRight() );
141  G4double calorimeterEDLeft( digitizer->GetCalorimeterEDLeft() );
142  G4double calorimeterEDRight( digitizer->GetCalorimeterEDRight() );
143  G4int calorimeterEDLeftMaxX( digitizer->GetCalorimeterEDLeftMaxX() );
144  G4int calorimeterEDLeftMaxY( digitizer->GetCalorimeterEDLeftMaxY() );
145  G4int calorimeterEDRightMaxX( digitizer->GetCalorimeterEDRightMaxX() );
146  G4int calorimeterEDRightMaxY( digitizer->GetCalorimeterEDRightMaxY() );
147 
148  const CexmcEnergyDepositCalorimeterCollection &
149  calorimeterEDLeftCollection(
150  digitizer->GetCalorimeterEDLeftCollection() );
151  const CexmcEnergyDepositCalorimeterCollection &
152  calorimeterEDRightCollection(
153  digitizer->GetCalorimeterEDRightCollection() );
154 
155  /* ATTENTION: return object in heap - must be freed by caller! */
156  return new CexmcEnergyDepositStore( monitorED, vetoCounterEDLeft,
157  vetoCounterEDRight, calorimeterEDLeft, calorimeterEDRight,
158  calorimeterEDLeftMaxX, calorimeterEDLeftMaxY,
159  calorimeterEDRightMaxX, calorimeterEDRightMaxY,
160  calorimeterEDLeftCollection, calorimeterEDRightCollection );
161 }
162 
163 
164 CexmcTrackPointsStore * CexmcEventAction::MakeTrackPointsStore(
165  const CexmcTrackPointsDigitizer * digitizer )
166 {
167  const CexmcTrackPointInfo &
168  monitorTP( digitizer->GetMonitorTP() );
169  const CexmcTrackPointInfo &
170  targetTPBeamParticle(
171  digitizer->GetTargetTPBeamParticle() );
172  const CexmcTrackPointInfo &
173  targetTPOutputParticle(
174  digitizer->GetTargetTPOutputParticle() );
175  const CexmcTrackPointInfo &
176  targetTPNucleusParticle(
177  digitizer->GetTargetTPNucleusParticle() );
178  const CexmcTrackPointInfo &
179  targetTPOutputParticleDecayProductParticle1(
180  digitizer->
181  GetTargetTPOutputParticleDecayProductParticle( 0 ) );
182  const CexmcTrackPointInfo &
183  targetTPOutputParticleDecayProductParticle2(
184  digitizer->
185  GetTargetTPOutputParticleDecayProductParticle( 1 ) );
186  const CexmcTrackPointInfo &
187  vetoCounterTPLeft(
188  digitizer->GetVetoCounterTPLeft() );
189  const CexmcTrackPointInfo &
190  vetoCounterTPRight(
191  digitizer->GetVetoCounterTPRight() );
192  const CexmcTrackPointInfo &
193  calorimeterTPLeft(
194  digitizer->GetCalorimeterTPLeft() );
195  const CexmcTrackPointInfo &
196  calorimeterTPRight(
197  digitizer->GetCalorimeterTPRight() );
198 
199  /* ATTENTION: return object in heap - must be freed by caller! */
200  return new CexmcTrackPointsStore( monitorTP, targetTPBeamParticle,
201  targetTPOutputParticle, targetTPNucleusParticle,
202  targetTPOutputParticleDecayProductParticle1,
203  targetTPOutputParticleDecayProductParticle2,
204  vetoCounterTPLeft, vetoCounterTPRight,
205  calorimeterTPLeft, calorimeterTPRight );
206 }
207 
208 
209 void CexmcEventAction::PrintEnergyDeposit(
210  const CexmcEnergyDepositStore * edStore )
211 {
212  G4cout << " --- Energy Deposit" << G4endl;
213  G4cout << " monitor : " <<
214  G4BestUnit( edStore->monitorED, "Energy" ) << G4endl;
215  G4cout << " vc (l) : " <<
216  G4BestUnit( edStore->vetoCounterEDLeft, "Energy" ) << G4endl;
217  G4cout << " vc (r) : " <<
218  G4BestUnit( edStore->vetoCounterEDRight, "Energy" ) << G4endl;
219  G4cout << " cal (l) : " <<
220  G4BestUnit( edStore->calorimeterEDLeft, "Energy" );
221  G4cout << edStore->calorimeterEDLeftCollection;
222  G4cout << " cal (r) : " <<
223  G4BestUnit( edStore->calorimeterEDRight, "Energy" );
224  G4cout << edStore->calorimeterEDRightCollection;
225 }
226 
227 
228 void CexmcEventAction::PrintTrackPoints(
229  const CexmcTrackPointsStore * tpStore )
230 {
231  if ( ! tpStore )
232  return;
233 
234  G4cout << " --- Track Points" << G4endl;
235  G4cout << " monitor : " << tpStore->monitorTP << G4endl;
236  G4cout << " target : " << tpStore->targetTPBeamParticle << G4endl;
237  G4cout << " : " << tpStore->targetTPOutputParticle << G4endl;
238  G4cout << " : " << tpStore->targetTPNucleusParticle << G4endl;
239  G4cout << " : " <<
240  tpStore->targetTPOutputParticleDecayProductParticle1 << G4endl;
241  G4cout << " : " <<
242  tpStore->targetTPOutputParticleDecayProductParticle2 << G4endl;
243  G4cout << " vc (l) : " << tpStore->vetoCounterTPLeft << G4endl;
244  G4cout << " vc (r) : " << tpStore->vetoCounterTPRight << G4endl;
245  G4cout << " cal (l) : " << tpStore->calorimeterTPLeft << G4endl;
246  G4cout << " cal (r) : " << tpStore->calorimeterTPRight << G4endl;
247  G4cout << " ---" << G4endl;
248  G4cout << " angle between the " <<
249  tpStore->targetTPOutputParticle.particle->GetParticleName() <<
250  " decay products : " <<
251  tpStore->targetTPOutputParticleDecayProductParticle1.directionWorld.
252  angle(
253  tpStore->targetTPOutputParticleDecayProductParticle2.directionWorld ) /
254  deg << " deg" << G4endl;
255 }
256 
257 
258 void CexmcEventAction::PrintProductionModelData(
259  const CexmcAngularRangeList & angularRanges,
260  const CexmcProductionModelData & pmData )
261 {
262  G4cout << " --- Triggered angular ranges: " << angularRanges;
263  G4cout << " --- Production model data: " << pmData;
264 }
265 
266 
267 void CexmcEventAction::PrintReconstructedData(
268  const CexmcAngularRangeList & triggeredRecAngularRanges,
269  const CexmcAngularRange & angularGap ) const
270 {
271  G4cout << " --- Reconstructed data: " << G4endl;
272  G4cout << " -- entry points:" << G4endl;
273  G4cout << " left: " << G4BestUnit(
274  reconstructor->GetCalorimeterEPLeftPosition(), "Length" ) << G4endl;
275  G4cout << " right: " << G4BestUnit(
276  reconstructor->GetCalorimeterEPRightPosition(), "Length" ) << G4endl;
277  G4cout << " target: " << G4BestUnit(
278  reconstructor->GetTargetEPPosition(), "Length" ) << G4endl;
279  G4cout << " -- the angle: " << reconstructor->GetTheAngle() / deg <<
280  " deg" << G4endl;
281  G4cout << " -- mass of the output particle: " << G4BestUnit(
282  reconstructor->GetOutputParticleMass(), "Energy" ) << G4endl;
283  G4cout << " -- mass of the nucleus output particle: " << G4BestUnit(
284  reconstructor->GetNucleusOutputParticleMass(), "Energy" ) << G4endl;
285  if ( reconstructor->IsMassCutUsed() )
286  {
287  if ( reconstructor->HasMassCutTriggered() )
288  G4cout << " < mass cut passed >" << G4endl;
289  else
290  G4cout << " < mass cut failed >" << G4endl;
291  }
292  if ( reconstructor->IsAbsorbedEnergyCutUsed() )
293  {
294  if ( reconstructor->HasAbsorbedEnergyCutTriggered() )
295  G4cout << " < absorbed energy cut passed >" << G4endl;
296  else
297  G4cout << " < absorbed energy cut failed >" << G4endl;
298  }
299  const CexmcProductionModelData & pmData(
300  reconstructor->GetProductionModelData() );
301  G4cout << " -- production model data: " << pmData;
302  G4cout << " -- triggered angular ranges: ";
303  if ( triggeredRecAngularRanges.empty() )
304  G4cout << "< orphan detected, gap " << angularGap << " >" << G4endl;
305  else
306  G4cout << triggeredRecAngularRanges;
307 }
308 
309 
310 #ifdef CEXMC_USE_ROOT
311 
312 void CexmcEventAction::FillEDTHistos( const CexmcEnergyDepositStore * edStore,
313  const CexmcAngularRangeList & triggeredAngularRanges ) const
314 {
315  CexmcHistoManager * histoManager( CexmcHistoManager::Instance() );
316 
317  histoManager->Add( CexmcAbsorbedEnergy_EDT_Histo, 0,
318  edStore->calorimeterEDLeft,
319  edStore->calorimeterEDRight );
320 
321  for ( CexmcAngularRangeList::const_iterator
322  k( triggeredAngularRanges.begin() );
323  k != triggeredAngularRanges.end(); ++k )
324  {
325  histoManager->Add( CexmcAbsEnInLeftCalorimeter_ARReal_EDT_Histo,
326  k->index, edStore->calorimeterEDLeft );
327  histoManager->Add( CexmcAbsEnInRightCalorimeter_ARReal_EDT_Histo,
328  k->index, edStore->calorimeterEDRight );
329  }
330 }
331 
332 
333 void CexmcEventAction::FillTPTHistos( const CexmcTrackPointsStore * tpStore,
334  const CexmcProductionModelData & pmData,
335  const CexmcAngularRangeList & triggeredAngularRanges ) const
336 {
337  CexmcHistoManager * histoManager( CexmcHistoManager::Instance() );
338 
339  if ( tpStore->monitorTP.IsValid() )
340  {
341  histoManager->Add( CexmcMomentumBP_TPT_Histo, 0,
342  tpStore->monitorTP.momentumAmp );
343  histoManager->Add( CexmcTPInMonitor_TPT_Histo, 0,
344  tpStore->monitorTP.positionLocal.x(),
345  tpStore->monitorTP.positionLocal.y() );
346  }
347 
348  if ( tpStore->targetTPOutputParticle.IsValid() )
349  {
350  histoManager->Add( CexmcTPInTarget_TPT_Histo, 0,
351  tpStore->targetTPOutputParticle.positionLocal.x(),
352  tpStore->targetTPOutputParticle.positionLocal.y(),
353  tpStore->targetTPOutputParticle.positionLocal.z() );
354  if ( histoManager->GetVerboseLevel() > 0 )
355  {
356  histoManager->Add( CexmcMomentumIP_TPT_Histo, 0,
357  pmData.incidentParticleLAB.rho() );
358  }
359  }
360 
361  for ( CexmcAngularRangeList::const_iterator
362  k( triggeredAngularRanges.begin() );
363  k != triggeredAngularRanges.end(); ++k )
364  {
365  if ( tpStore->calorimeterTPLeft.IsValid() )
366  {
367  /* kinetic energy and momentum of gamma are equal */
368  G4double kinEnergy( tpStore->calorimeterTPLeft.momentumAmp );
369  histoManager->Add( CexmcKinEnAtLeftCalorimeter_ARReal_TPT_Histo,
370  k->index, kinEnergy );
371  }
372  if ( tpStore->calorimeterTPRight.IsValid() )
373  {
374  G4double kinEnergy( tpStore->calorimeterTPRight.momentumAmp );
375  histoManager->Add( CexmcKinEnAtRightCalorimeter_ARReal_TPT_Histo,
376  k->index, kinEnergy );
377  }
378  if ( tpStore->targetTPOutputParticle.IsValid() )
379  {
380  histoManager->Add( CexmcTPInTarget_ARReal_TPT_Histo, k->index,
381  tpStore->targetTPOutputParticle.positionLocal.x(),
382  tpStore->targetTPOutputParticle.positionLocal.y(),
383  tpStore->targetTPOutputParticle.positionLocal.z() );
384  histoManager->Add( CexmcKinEnOP_LAB_ARReal_TPT_Histo, k->index,
385  opKinEnergy );
386  histoManager->Add( CexmcAngleOP_SCM_ARReal_TPT_Histo, k->index,
387  pmData.outputParticleSCM.cosTheta() );
388  }
389  if ( tpStore->targetTPOutputParticleDecayProductParticle1.IsValid() &&
390  tpStore->targetTPOutputParticleDecayProductParticle2.IsValid() )
391  {
392  G4double openAngle(
393  tpStore->targetTPOutputParticleDecayProductParticle1.
394  directionWorld.angle( tpStore->
395  targetTPOutputParticleDecayProductParticle2.
396  directionWorld ) / deg );
397  histoManager->Add( CexmcOpenAngle_ARReal_TPT_Histo, k->index,
398  openAngle );
399  }
400  }
401 }
402 
403 
404 void CexmcEventAction::FillRTHistos( G4bool reconstructorHasFullTrigger,
405  const CexmcEnergyDepositStore * edStore,
406  const CexmcTrackPointsStore * tpStore,
407  const CexmcProductionModelData & pmData,
408  const CexmcAngularRangeList & triggeredAngularRanges ) const
409 {
410  CexmcHistoManager * histoManager( CexmcHistoManager::Instance() );
411 
412  G4double opMass( reconstructor->GetOutputParticleMass() );
413  G4double nopMass( reconstructor->GetNucleusOutputParticleMass() );
414 
415  histoManager->Add( CexmcRecMasses_EDT_Histo, 0, opMass, nopMass );
416 
417  for ( CexmcAngularRangeList::const_iterator
418  k( triggeredAngularRanges.begin() );
419  k != triggeredAngularRanges.end(); ++k )
420  {
421  if ( tpStore->calorimeterTPLeft.IsValid() )
422  {
423  histoManager->Add( CexmcOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
424  k->index,
425  tpStore->calorimeterTPLeft.positionLocal.x(),
426  tpStore->calorimeterTPLeft.positionLocal.y() );
427  }
428  if ( tpStore->calorimeterTPRight.IsValid() )
429  {
430  histoManager->Add( CexmcOPDPAtRightCalorimeter_ARReal_EDT_Histo,
431  k->index,
432  tpStore->calorimeterTPRight.positionLocal.x(),
433  tpStore->calorimeterTPRight.positionLocal.y() );
434  }
435  histoManager->Add( CexmcRecOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
436  k->index,
437  reconstructor->GetCalorimeterEPLeftPosition().x(),
438  reconstructor->GetCalorimeterEPLeftPosition().y() );
439  histoManager->Add( CexmcRecOPDPAtRightCalorimeter_ARReal_EDT_Histo,
440  k->index,
441  reconstructor->GetCalorimeterEPRightPosition().x(),
442  reconstructor->GetCalorimeterEPRightPosition().y() );
443  }
444 
445  if ( ! reconstructorHasFullTrigger )
446  return;
447 
448  if ( tpStore->monitorTP.IsValid() )
449  {
450  histoManager->Add( CexmcMomentumBP_RT_Histo, 0,
451  tpStore->monitorTP.momentumAmp );
452  }
453 
454  if ( tpStore->targetTPOutputParticle.IsValid() )
455  {
456  histoManager->Add( CexmcTPInTarget_RT_Histo, 0,
457  tpStore->targetTPOutputParticle.positionLocal.x(),
458  tpStore->targetTPOutputParticle.positionLocal.y(),
459  tpStore->targetTPOutputParticle.positionLocal.z() );
460  }
461 
462  histoManager->Add( CexmcRecMasses_RT_Histo, 0,
463  reconstructor->GetOutputParticleMass(),
464  reconstructor->GetNucleusOutputParticleMass() );
465 
466  histoManager->Add( CexmcAbsorbedEnergy_RT_Histo, 0,
467  edStore->calorimeterEDLeft,
468  edStore->calorimeterEDRight );
469 
470  G4double recCosTheta( reconstructor->GetProductionModelData().
471  outputParticleSCM.cosTheta());
472 
473  for ( CexmcAngularRangeList::const_iterator
474  k( triggeredAngularRanges.begin() );
475  k != triggeredAngularRanges.end(); ++k )
476  {
477  histoManager->Add( CexmcRecMassOP_ARReal_RT_Histo, k->index, opMass );
478  histoManager->Add( CexmcRecMassNOP_ARReal_RT_Histo, k->index, nopMass );
479  if ( tpStore->calorimeterTPLeft.IsValid() )
480  {
481  G4double kinEnergy( tpStore->calorimeterTPLeft.momentumAmp );
482  histoManager->Add( CexmcKinEnAtLeftCalorimeter_ARReal_RT_Histo,
483  k->index, kinEnergy );
484  histoManager->Add( CexmcMissEnFromLeftCalorimeter_ARReal_RT_Histo,
485  k->index,
486  kinEnergy - edStore->calorimeterEDLeft );
487  }
488  if ( tpStore->calorimeterTPRight.IsValid() )
489  {
490  G4double kinEnergy( tpStore->calorimeterTPRight.momentumAmp );
491  histoManager->Add( CexmcKinEnAtRightCalorimeter_ARReal_RT_Histo,
492  k->index, kinEnergy );
493  histoManager->Add( CexmcMissEnFromRightCalorimeter_ARReal_RT_Histo,
494  k->index,
495  kinEnergy - edStore->calorimeterEDRight );
496  }
497  if ( tpStore->targetTPOutputParticle.IsValid() )
498  {
499  histoManager->Add( CexmcTPInTarget_ARReal_RT_Histo, k->index,
500  tpStore->targetTPOutputParticle.positionLocal.x(),
501  tpStore->targetTPOutputParticle.positionLocal.y(),
502  tpStore->targetTPOutputParticle.positionLocal.z() );
503  histoManager->Add( CexmcKinEnOP_LAB_ARReal_RT_Histo, k->index,
504  opKinEnergy );
505  histoManager->Add( CexmcAngleOP_SCM_ARReal_RT_Histo, k->index,
506  pmData.outputParticleSCM.cosTheta() );
507  G4double diffCosTheta( pmData.outputParticleSCM.cosTheta() -
508  recCosTheta );
509  histoManager->Add( CexmcDiffAngleOP_SCM_ARReal_RT_Histo, k->index,
510  diffCosTheta );
511  }
512  if ( tpStore->targetTPOutputParticleDecayProductParticle1.IsValid() &&
513  tpStore->targetTPOutputParticleDecayProductParticle2.IsValid() )
514  {
515  G4double openAngle(
516  tpStore->targetTPOutputParticleDecayProductParticle1.
517  directionWorld.angle( tpStore->
518  targetTPOutputParticleDecayProductParticle2.
519  directionWorld ) / deg );
520  histoManager->Add( CexmcOpenAngle_ARReal_RT_Histo, k->index,
521  openAngle );
522  G4double diffOpenAngle( openAngle - reconstructor->GetTheAngle() /
523  deg );
524  histoManager->Add( CexmcDiffOpenAngle_ARReal_RT_Histo, k->index,
525  diffOpenAngle );
526  }
527  if ( tpStore->calorimeterTPLeft.IsValid() )
528  {
529  histoManager->Add( CexmcOPDPAtLeftCalorimeter_ARReal_RT_Histo,
530  k->index,
531  tpStore->calorimeterTPLeft.positionLocal.x(),
532  tpStore->calorimeterTPLeft.positionLocal.y() );
533  }
534  if ( tpStore->calorimeterTPRight.IsValid() )
535  {
536  histoManager->Add( CexmcOPDPAtRightCalorimeter_ARReal_RT_Histo,
537  k->index,
538  tpStore->calorimeterTPRight.positionLocal.x(),
539  tpStore->calorimeterTPRight.positionLocal.y() );
540  }
541  histoManager->Add( CexmcAbsEnInLeftCalorimeter_ARReal_RT_Histo,
542  k->index, edStore->calorimeterEDLeft );
543  histoManager->Add( CexmcAbsEnInRightCalorimeter_ARReal_RT_Histo,
544  k->index, edStore->calorimeterEDRight );
545  histoManager->Add( CexmcRecAngleOP_SCM_ARReal_RT_Histo,
546  k->index, recCosTheta );
547  histoManager->Add( CexmcRecOpenAngle_ARReal_RT_Histo,
548  k->index, reconstructor->GetTheAngle() / deg );
549  histoManager->Add( CexmcRecOPDPAtLeftCalorimeter_ARReal_RT_Histo,
550  k->index,
551  reconstructor->GetCalorimeterEPLeftPosition().x(),
552  reconstructor->GetCalorimeterEPLeftPosition().y() );
553  histoManager->Add( CexmcRecOPDPAtRightCalorimeter_ARReal_RT_Histo,
554  k->index,
555  reconstructor->GetCalorimeterEPRightPosition().x(),
556  reconstructor->GetCalorimeterEPRightPosition().y() );
557  }
558 }
559 
560 #endif
561 
562 
563 void CexmcEventAction::DrawTrajectories( const G4Event * event )
564 {
565  G4VisManager * visManager( static_cast< G4VisManager * >(
566  G4VVisManager::GetConcreteInstance() ) );
567  if ( ! visManager || ! visManager->GetCurrentGraphicsSystem() )
568  return;
569 
570  G4int nTraj( 0 );
571  G4TrajectoryContainer * trajContainer( event->GetTrajectoryContainer() );
572 
573  if ( ! trajContainer )
574  return;
575 
576  nTraj = trajContainer->entries();
577 
578  for ( int i( 0 ); i < nTraj; ++i )
579  {
580  G4VTrajectory * traj( ( *trajContainer )[ i ] );
581  traj->DrawTrajectory();
582  }
583 }
584 
585 
586 void CexmcEventAction::DrawTrackPoints(
587  const CexmcTrackPointsStore * tpStore ) const
588 {
589  G4VisManager * visManager( static_cast< G4VisManager * >(
590  G4VVisManager::GetConcreteInstance() ) );
591  if ( ! visManager || ! visManager->GetCurrentGraphicsSystem() )
592  return;
593 
594  G4Circle circle;
595  G4VisAttributes visAttributes( CexmcTrackPointsMarkerColour );
596  circle.SetScreenSize( CexmcSmallCircleScreenSize );
597  circle.SetFillStyle( G4Circle::filled );
598  circle.SetVisAttributes( visAttributes );
599 
600  if ( tpStore->monitorTP.IsValid() )
601  {
602  circle.SetPosition( tpStore->monitorTP.positionWorld );
603  visManager->Draw( circle );
604  }
605 
606  if ( tpStore->targetTPBeamParticle.IsValid() )
607  {
608  circle.SetPosition( tpStore->targetTPBeamParticle.positionWorld );
609  visManager->Draw( circle );
610  }
611 
612  if ( tpStore->targetTPOutputParticle.IsValid() )
613  {
614  circle.SetPosition( tpStore->targetTPOutputParticle.positionWorld );
615  visManager->Draw( circle );
616  }
617 
618  if ( tpStore->vetoCounterTPLeft.IsValid() )
619  {
620  circle.SetPosition( tpStore->vetoCounterTPLeft.positionWorld );
621  visManager->Draw( circle );
622  }
623 
624  if ( tpStore->vetoCounterTPRight.IsValid() )
625  {
626  circle.SetPosition( tpStore->vetoCounterTPRight.positionWorld );
627  visManager->Draw( circle );
628  }
629 
630  if ( tpStore->calorimeterTPLeft.IsValid() )
631  {
632  circle.SetPosition( tpStore->calorimeterTPLeft.positionWorld );
633  visManager->Draw( circle );
634  }
635 
636  if ( tpStore->calorimeterTPRight.IsValid() )
637  {
638  circle.SetPosition( tpStore->calorimeterTPRight.positionWorld );
639  visManager->Draw( circle );
640  }
641 }
642 
643 
644 void CexmcEventAction::DrawReconstructionData( void )
645 {
646  G4VisManager * visManager( static_cast< G4VisManager * >(
647  G4VVisManager::GetConcreteInstance() ) );
648  if ( ! visManager || ! visManager->GetCurrentGraphicsSystem() )
649  return;
650 
651  G4Circle circle( reconstructor->GetTargetEPWorldPosition() );
652  circle.SetScreenSize( CexmcSmallCircleScreenSize );
653  circle.SetFillStyle( G4Circle::filled );
654  G4VisAttributes visAttributes( CexmcRecTrackPointsMarkerColour );
655  circle.SetVisAttributes( visAttributes );
656  visManager->Draw( circle );
657 
658  circle.SetScreenSize( CexmcBigCircleScreenSize );
659  circle.SetPosition( reconstructor->GetCalorimeterEPLeftWorldPosition() );
660  visManager->Draw( circle );
661 
662  circle.SetPosition( reconstructor->GetCalorimeterEPRightWorldPosition() );
663  visManager->Draw( circle );
664 }
665 
666 
667 void CexmcEventAction::UpdateRunHits(
668  const CexmcAngularRangeList & aRangesReal,
669  const CexmcAngularRangeList & aRangesRec,
670  G4bool tpDigitizerHasTriggered,
671  G4bool edDigitizerHasTriggered,
672  G4bool edDigitizerMonitorHasTriggered,
673  G4bool reconstructorHasFullTrigger,
674  const CexmcAngularRange & aGap )
675 {
676  G4RunManager * runManager( G4RunManager::GetRunManager() );
677  const CexmcRun * run( static_cast< const CexmcRun * >(
678  runManager->GetCurrentRun() ) );
679  CexmcRun * theRun( const_cast< CexmcRun * >( run ) );
680 
681  if ( tpDigitizerHasTriggered )
682  {
683  for ( CexmcAngularRangeList::const_iterator k( aRangesReal.begin() );
684  k != aRangesReal.end(); ++k )
685  {
686  theRun->IncrementNmbOfHitsSampledFull( k->index );
687  if ( edDigitizerMonitorHasTriggered )
688  theRun->IncrementNmbOfHitsSampled( k->index );
689  if ( reconstructorHasFullTrigger )
690  theRun->IncrementNmbOfHitsTriggeredRealRange( k->index );
691  }
692  if ( reconstructorHasFullTrigger )
693  {
694  if ( aRangesRec.empty() )
695  {
696  theRun->IncrementNmbOfOrphanHits( aGap.index );
697  }
698  else
699  {
700  for ( CexmcAngularRangeList::const_iterator
701  k( aRangesRec.begin() ); k != aRangesRec.end(); ++k )
702  {
703  theRun->IncrementNmbOfHitsTriggeredRecRange( k->index );
704  }
705  }
706  }
707  }
708  else
709  {
710  if ( edDigitizerHasTriggered )
711  theRun->IncrementNmbOfFalseHitsTriggeredEDT();
712  if ( reconstructorHasFullTrigger )
713  theRun->IncrementNmbOfFalseHitsTriggeredRec();
714  }
715 }
716 
717 
718 #ifdef CEXMC_USE_PERSISTENCY
719 
720 void CexmcEventAction::SaveEvent( const G4Event * event,
721  G4bool edDigitizerHasTriggered,
722  const CexmcEnergyDepositStore * edStore,
723  const CexmcTrackPointsStore * tpStore,
724  const CexmcProductionModelData & pmData )
725 {
726  CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
727  G4RunManager::GetRunManager() ) );
728  if ( ! runManager->ProjectIsSaved() )
729  return;
730 
731  if ( runManager->GetEventDataVerboseLevel() == CexmcWriteNoEventData )
732  return;
733 
734  if ( ! edDigitizerHasTriggered && runManager->GetEventDataVerboseLevel() !=
735  CexmcWriteEventDataOnEveryTPT )
736  return;
737 
738  boost::archive::binary_oarchive * archive(
739  runManager->GetEventsArchive() );
740  if ( archive )
741  {
742  CexmcEventSObject sObject = { event->GetEventID(),
743  edDigitizerHasTriggered, edStore->monitorED,
744  edStore->vetoCounterEDLeft, edStore->vetoCounterEDRight,
745  edStore->calorimeterEDLeft, edStore->calorimeterEDRight,
746  edStore->calorimeterEDLeftCollection,
747  edStore->calorimeterEDRightCollection,
748  tpStore->monitorTP, tpStore->targetTPBeamParticle,
749  tpStore->targetTPOutputParticle, tpStore->targetTPNucleusParticle,
750  tpStore->targetTPOutputParticleDecayProductParticle1,
751  tpStore->targetTPOutputParticleDecayProductParticle2,
752  tpStore->vetoCounterTPLeft, tpStore->vetoCounterTPRight,
753  tpStore->calorimeterTPLeft, tpStore->calorimeterTPRight, pmData };
754  archive->operator<<( sObject );
755  const CexmcRun * run( static_cast< const CexmcRun * >(
756  runManager->GetCurrentRun() ) );
757  CexmcRun * theRun( const_cast< CexmcRun * >( run ) );
758  theRun->IncrementNmbOfSavedEvents();
759  }
760 }
761 
762 
763 void CexmcEventAction::SaveEventFast( const G4Event * event,
764  G4bool tpDigitizerHasTriggered,
765  G4bool edDigitizerHasTriggered,
766  G4bool edDigitizerMonitorHasTriggered,
767  G4double opCosThetaSCM )
768 {
769  CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
770  G4RunManager::GetRunManager() ) );
771  if ( ! runManager->ProjectIsSaved() )
772  return;
773 
774  if ( runManager->GetEventDataVerboseLevel() == CexmcWriteNoEventData )
775  return;
776 
777  boost::archive::binary_oarchive * archive(
778  runManager->GetFastEventsArchive() );
779  if ( archive )
780  {
781  if ( ! tpDigitizerHasTriggered )
782  opCosThetaSCM = CexmcInvalidCosTheta;
783 
784  CexmcEventFastSObject sObject = { event->GetEventID(), opCosThetaSCM,
785  edDigitizerHasTriggered,
786  edDigitizerMonitorHasTriggered };
787  archive->operator<<( sObject );
788  const CexmcRun * run( static_cast< const CexmcRun * >(
789  runManager->GetCurrentRun() ) );
790  CexmcRun * theRun( const_cast< CexmcRun * >( run ) );
791  theRun->IncrementNmbOfSavedFastEvents();
792  }
793 }
794 
795 #endif
796 
797 
798 void CexmcEventAction::EndOfEventAction( const G4Event * event )
799 {
800  G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
801  CexmcEnergyDepositDigitizer * energyDepositDigitizer(
802  static_cast< CexmcEnergyDepositDigitizer * >( digiManager->
803  FindDigitizerModule( CexmcEDDigitizerName ) ) );
804  CexmcTrackPointsDigitizer * trackPointsDigitizer(
805  static_cast< CexmcTrackPointsDigitizer * >( digiManager->
806  FindDigitizerModule( CexmcTPDigitizerName ) ) );
807 
808  energyDepositDigitizer->Digitize();
809  trackPointsDigitizer->Digitize();
810 
811  G4bool edDigitizerMonitorHasTriggered(
812  energyDepositDigitizer->MonitorHasTriggered() );
813  G4bool edDigitizerHasTriggered( false );
814 
815  CexmcEventInfo * eventInfo( static_cast< CexmcEventInfo * >(
816  event->GetUserInformation() ) );
817  if ( ! eventInfo || eventInfo->EdTriggerIsOk() )
818  edDigitizerHasTriggered = energyDepositDigitizer->HasTriggered();
819 
820  G4bool tpDigitizerHasTriggered( trackPointsDigitizer->HasTriggered() );
821  G4bool reconstructorHasBasicTrigger( false );
822  G4bool reconstructorHasFullTrigger( false );
823 
824  CexmcEnergyDepositStore * edStore( MakeEnergyDepositStore(
825  energyDepositDigitizer ) );
826  CexmcTrackPointsStore * tpStore( MakeTrackPointsStore(
827  trackPointsDigitizer ) );
828 
829  try
830  {
831  CexmcProductionModel * productionModel(
832  physicsManager->GetProductionModel() );
833 
834  if ( ! productionModel )
835  throw CexmcException( CexmcWeirdException );
836 
837  const CexmcAngularRangeList & angularRanges(
838  productionModel->GetAngularRanges() );
839  const CexmcAngularRangeList & triggeredAngularRanges(
840  productionModel->GetTriggeredAngularRanges() );
841  const CexmcProductionModelData & pmData(
842  productionModel->GetProductionModelData() );
843 
844  if ( edDigitizerHasTriggered )
845  {
846  reconstructor->Reconstruct( edStore );
847  reconstructorHasBasicTrigger = reconstructor->HasBasicTrigger();
848  reconstructorHasFullTrigger = reconstructor->HasFullTrigger();
849  }
850 
851  CexmcAngularRangeList triggeredRecAngularRanges;
852 
853  if ( reconstructorHasBasicTrigger )
854  {
855  for ( CexmcAngularRangeList::const_iterator
856  k( angularRanges.begin() ); k != angularRanges.end(); ++k )
857  {
858  G4double cosTheta( reconstructor->GetProductionModelData().
859  outputParticleSCM.cosTheta() );
860  if ( cosTheta <= k->top && cosTheta > k->bottom )
861  triggeredRecAngularRanges.push_back( CexmcAngularRange(
862  k->top, k->bottom, k->index ) );
863  }
864  }
865 
866  CexmcAngularRange angularGap( 0.0, 0.0, 0 );
867  if ( triggeredRecAngularRanges.empty() )
868  {
869  CexmcAngularRangeList angularGaps;
870  GetAngularGaps( angularRanges, angularGaps );
871  for ( CexmcAngularRangeList::const_iterator
872  k( angularGaps.begin() ); k != angularGaps.end(); ++k )
873  {
874  G4double cosTheta( reconstructor->GetProductionModelData().
875  outputParticleSCM.cosTheta() );
876  if ( cosTheta <= k->top && cosTheta > k->bottom )
877  {
878  angularGap = *k;
879  break;
880  }
881  }
882  }
883 
884  UpdateRunHits( triggeredAngularRanges, triggeredRecAngularRanges,
885  tpDigitizerHasTriggered, edDigitizerHasTriggered,
886  edDigitizerMonitorHasTriggered,
887  reconstructorHasFullTrigger, angularGap );
888 
889  if ( verbose > 0 )
890  {
891  G4bool printMessages( verbose > 3 ||
892  ( ( verbose == 1 ) && tpDigitizerHasTriggered ) ||
893  ( ( verbose == 2 ) && edDigitizerHasTriggered ) ||
894  ( ( verbose == 3 ) && ( tpDigitizerHasTriggered ||
895  edDigitizerHasTriggered ) ) );
896  if ( printMessages )
897  {
898  G4cout << "Event " << event->GetEventID() << G4endl;
899  if ( tpDigitizerHasTriggered )
900  {
901  PrintTrackPoints( tpStore );
902  PrintProductionModelData( triggeredAngularRanges, pmData );
903  }
904  if ( reconstructorHasBasicTrigger )
905  PrintReconstructedData( triggeredRecAngularRanges,
906  angularGap );
907  if ( edDigitizerHasTriggered )
908  PrintEnergyDeposit( edStore );
909  }
910  }
911 
912  if ( verboseDraw > 0 )
913  {
914  G4bool drawTrajectories( verboseDraw > 3 ||
915  ( ( verboseDraw == 1 ) && tpDigitizerHasTriggered ) ||
916  ( ( verboseDraw == 2 ) && edDigitizerHasTriggered ) ||
917  ( ( verboseDraw == 3 ) && ( tpDigitizerHasTriggered ||
918  edDigitizerHasTriggered ) ) );
919  if ( drawTrajectories )
920  {
921  DrawTrajectories( event );
922  if ( tpDigitizerHasTriggered )
923  DrawTrackPoints( tpStore );
924  if ( reconstructorHasBasicTrigger )
925  DrawReconstructionData();
926  }
927  }
928 
929 #ifdef CEXMC_USE_PERSISTENCY
930  if ( edDigitizerHasTriggered || tpDigitizerHasTriggered )
931  {
932  SaveEventFast( event, tpDigitizerHasTriggered,
933  edDigitizerHasTriggered,
934  edDigitizerMonitorHasTriggered,
935  pmData.outputParticleSCM.cosTheta() );
936  SaveEvent( event, edDigitizerHasTriggered, edStore, tpStore,
937  pmData );
938  }
939 #endif
940 
941 #ifdef CEXMC_USE_ROOT
942  /* opKinEnergy will be used in several histos */
943  if ( tpStore->targetTPOutputParticle.IsValid() )
944  {
945  opKinEnergy = CexmcGetKinEnergy(
946  tpStore->targetTPOutputParticle.momentumAmp,
947  tpStore->targetTPOutputParticle.particle->GetPDGMass() );
948  }
949 
950  if ( edDigitizerHasTriggered )
951  FillEDTHistos( edStore, triggeredAngularRanges );
952 
953  /* fill TPT histos only when the monitor has triggered because events
954  * when it was missed have less value for us */
955  if ( tpDigitizerHasTriggered && edDigitizerMonitorHasTriggered )
956  FillTPTHistos( tpStore, pmData, triggeredAngularRanges );
957 
958  if ( reconstructorHasBasicTrigger )
959  FillRTHistos( reconstructorHasFullTrigger, edStore, tpStore,
960  pmData, triggeredAngularRanges );
961 #endif
962 
963  G4Event * theEvent( const_cast< G4Event * >( event ) );
964  if ( eventInfo )
965  {
966  delete eventInfo;
967  theEvent->SetUserInformation( NULL );
968  }
969  theEvent->SetUserInformation( new CexmcEventInfo(
970  edDigitizerHasTriggered,
971  tpDigitizerHasTriggered,
972  reconstructorHasFullTrigger ) );
973  }
974  catch ( CexmcException & e )
975  {
976  G4cout << e.what() << G4endl;
977  }
978  catch ( ... )
979  {
980  G4cout << "Unknown exception caught" << G4endl;
981  }
982 
983  delete edStore;
984  delete tpStore;
985 }
986 
Generated on Sat May 25 2013 14:32:06 for Geant4 by   doxygen 1.8.4