| Geant4 Cross Reference |
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
992