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