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