File indexing completed on 2025-01-31 09:21:52
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044 #include <cmath>
0045 #ifdef CEXMC_USE_PERSISTENCY
0046 #include <boost/archive/binary_oarchive.hpp>
0047 #endif
0048 #include <G4DigiManager.hh>
0049 #include <G4Event.hh>
0050 #include <G4Circle.hh>
0051 #include <G4VisAttributes.hh>
0052 #include <G4VisManager.hh>
0053 #include <G4VTrajectory.hh>
0054 #include <G4TrajectoryContainer.hh>
0055 #include <G4Colour.hh>
0056 #include <G4SystemOfUnits.hh>
0057 #include "CexmcEventAction.hh"
0058 #include "CexmcEventActionMessenger.hh"
0059 #include "CexmcEventInfo.hh"
0060 #include "CexmcEventSObject.hh"
0061 #include "CexmcEventFastSObject.hh"
0062 #include "CexmcTrackingAction.hh"
0063 #include "CexmcChargeExchangeReconstructor.hh"
0064 #include "CexmcRunManager.hh"
0065 #include "CexmcHistoManager.hh"
0066 #include "CexmcRun.hh"
0067 #include "CexmcPhysicsManager.hh"
0068 #include "CexmcProductionModel.hh"
0069 #include "CexmcProductionModelData.hh"
0070 #include "CexmcEnergyDepositDigitizer.hh"
0071 #include "CexmcEnergyDepositStore.hh"
0072 #include "CexmcTrackPointsDigitizer.hh"
0073 #include "CexmcTrackPointsStore.hh"
0074 #include "CexmcTrackPointInfo.hh"
0075 #include "CexmcException.hh"
0076 #include "CexmcCommon.hh"
0077
0078
0079 namespace
0080 {
0081 G4double CexmcSmallCircleScreenSize( 5.0 );
0082 G4double CexmcBigCircleScreenSize( 10.0 );
0083 G4Colour CexmcTrackPointsMarkerColour( 0.0, 1.0, 0.4 );
0084 G4Colour CexmcRecTrackPointsMarkerColour( 1.0, 0.4, 0.0 );
0085
0086 #ifdef CEXMC_USE_ROOT
0087 inline G4double CexmcGetKinEnergy( G4double momentumAmp, G4double mass )
0088 {
0089 return std::sqrt( momentumAmp * momentumAmp + mass * mass ) - mass;
0090 }
0091 #endif
0092 }
0093
0094
0095 CexmcEventAction::CexmcEventAction( CexmcPhysicsManager * physicsManager_,
0096 G4int verbose_ ) :
0097 physicsManager( physicsManager_ ), reconstructor( NULL ),
0098 #ifdef CEXMC_USE_ROOT
0099 opKinEnergy( 0. ),
0100 #endif
0101 verbose( verbose_ ), verboseDraw( 4 ), messenger( NULL )
0102 {
0103 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
0104 digiManager->AddNewModule( new CexmcEnergyDepositDigitizer(
0105 CexmcEDDigitizerName ) );
0106 digiManager->AddNewModule( new CexmcTrackPointsDigitizer(
0107 CexmcTPDigitizerName ) );
0108 reconstructor = new CexmcChargeExchangeReconstructor(
0109 physicsManager->GetProductionModel() );
0110 messenger = new CexmcEventActionMessenger( this );
0111 }
0112
0113
0114 CexmcEventAction::~CexmcEventAction()
0115 {
0116 delete reconstructor;
0117 delete messenger;
0118 }
0119
0120
0121 void CexmcEventAction::BeamParticleChangeHook( void )
0122 {
0123 reconstructor->SetupBeamParticle();
0124 }
0125
0126
0127 void CexmcEventAction::BeginOfEventAction( const G4Event * )
0128 {
0129 G4RunManager * runManager( G4RunManager::GetRunManager() );
0130 CexmcTrackingAction * trackingAction
0131 ( static_cast< CexmcTrackingAction * >(
0132 const_cast< G4UserTrackingAction * >(
0133 runManager->GetUserTrackingAction() ) ) );
0134 trackingAction->BeginOfEventAction();
0135
0136 physicsManager->ResetNumberOfTriggeredStudiedInteractions();
0137 }
0138
0139
0140 CexmcEnergyDepositStore * CexmcEventAction::MakeEnergyDepositStore(
0141 const CexmcEnergyDepositDigitizer * digitizer )
0142 {
0143 G4double monitorED( digitizer->GetMonitorED() );
0144 G4double vetoCounterEDLeft( digitizer->GetVetoCounterEDLeft() );
0145 G4double vetoCounterEDRight( digitizer->GetVetoCounterEDRight() );
0146 G4double calorimeterEDLeft( digitizer->GetCalorimeterEDLeft() );
0147 G4double calorimeterEDRight( digitizer->GetCalorimeterEDRight() );
0148 G4int calorimeterEDLeftMaxX( digitizer->GetCalorimeterEDLeftMaxX() );
0149 G4int calorimeterEDLeftMaxY( digitizer->GetCalorimeterEDLeftMaxY() );
0150 G4int calorimeterEDRightMaxX( digitizer->GetCalorimeterEDRightMaxX() );
0151 G4int calorimeterEDRightMaxY( digitizer->GetCalorimeterEDRightMaxY() );
0152
0153 const CexmcEnergyDepositCalorimeterCollection &
0154 calorimeterEDLeftCollection(
0155 digitizer->GetCalorimeterEDLeftCollection() );
0156 const CexmcEnergyDepositCalorimeterCollection &
0157 calorimeterEDRightCollection(
0158 digitizer->GetCalorimeterEDRightCollection() );
0159
0160
0161 return new CexmcEnergyDepositStore( monitorED, vetoCounterEDLeft,
0162 vetoCounterEDRight, calorimeterEDLeft, calorimeterEDRight,
0163 calorimeterEDLeftMaxX, calorimeterEDLeftMaxY,
0164 calorimeterEDRightMaxX, calorimeterEDRightMaxY,
0165 calorimeterEDLeftCollection, calorimeterEDRightCollection );
0166 }
0167
0168
0169 CexmcTrackPointsStore * CexmcEventAction::MakeTrackPointsStore(
0170 const CexmcTrackPointsDigitizer * digitizer )
0171 {
0172 const CexmcTrackPointInfo &
0173 monitorTP( digitizer->GetMonitorTP() );
0174 const CexmcTrackPointInfo &
0175 targetTPBeamParticle(
0176 digitizer->GetTargetTPBeamParticle() );
0177 const CexmcTrackPointInfo &
0178 targetTPOutputParticle(
0179 digitizer->GetTargetTPOutputParticle() );
0180 const CexmcTrackPointInfo &
0181 targetTPNucleusParticle(
0182 digitizer->GetTargetTPNucleusParticle() );
0183 const CexmcTrackPointInfo &
0184 targetTPOutputParticleDecayProductParticle1(
0185 digitizer->
0186 GetTargetTPOutputParticleDecayProductParticle( 0 ) );
0187 const CexmcTrackPointInfo &
0188 targetTPOutputParticleDecayProductParticle2(
0189 digitizer->
0190 GetTargetTPOutputParticleDecayProductParticle( 1 ) );
0191 const CexmcTrackPointInfo &
0192 vetoCounterTPLeft(
0193 digitizer->GetVetoCounterTPLeft() );
0194 const CexmcTrackPointInfo &
0195 vetoCounterTPRight(
0196 digitizer->GetVetoCounterTPRight() );
0197 const CexmcTrackPointInfo &
0198 calorimeterTPLeft(
0199 digitizer->GetCalorimeterTPLeft() );
0200 const CexmcTrackPointInfo &
0201 calorimeterTPRight(
0202 digitizer->GetCalorimeterTPRight() );
0203
0204
0205 return new CexmcTrackPointsStore( monitorTP, targetTPBeamParticle,
0206 targetTPOutputParticle, targetTPNucleusParticle,
0207 targetTPOutputParticleDecayProductParticle1,
0208 targetTPOutputParticleDecayProductParticle2,
0209 vetoCounterTPLeft, vetoCounterTPRight,
0210 calorimeterTPLeft, calorimeterTPRight );
0211 }
0212
0213
0214 void CexmcEventAction::PrintEnergyDeposit(
0215 const CexmcEnergyDepositStore * edStore )
0216 {
0217 G4cout << " --- Energy Deposit" << G4endl;
0218 G4cout << " monitor : " <<
0219 G4BestUnit( edStore->monitorED, "Energy" ) << G4endl;
0220 G4cout << " vc (l) : " <<
0221 G4BestUnit( edStore->vetoCounterEDLeft, "Energy" ) << G4endl;
0222 G4cout << " vc (r) : " <<
0223 G4BestUnit( edStore->vetoCounterEDRight, "Energy" ) << G4endl;
0224 G4cout << " cal (l) : " <<
0225 G4BestUnit( edStore->calorimeterEDLeft, "Energy" );
0226 G4cout << edStore->calorimeterEDLeftCollection;
0227 G4cout << " cal (r) : " <<
0228 G4BestUnit( edStore->calorimeterEDRight, "Energy" );
0229 G4cout << edStore->calorimeterEDRightCollection;
0230 }
0231
0232
0233 void CexmcEventAction::PrintTrackPoints(
0234 const CexmcTrackPointsStore * tpStore )
0235 {
0236 if ( ! tpStore )
0237 return;
0238
0239 G4cout << " --- Track Points" << G4endl;
0240 G4cout << " monitor : " << tpStore->monitorTP << G4endl;
0241 G4cout << " target : " << tpStore->targetTPBeamParticle << G4endl;
0242 G4cout << " : " << tpStore->targetTPOutputParticle << G4endl;
0243 G4cout << " : " << tpStore->targetTPNucleusParticle << G4endl;
0244 G4cout << " : " <<
0245 tpStore->targetTPOutputParticleDecayProductParticle1 << G4endl;
0246 G4cout << " : " <<
0247 tpStore->targetTPOutputParticleDecayProductParticle2 << G4endl;
0248 G4cout << " vc (l) : " << tpStore->vetoCounterTPLeft << G4endl;
0249 G4cout << " vc (r) : " << tpStore->vetoCounterTPRight << G4endl;
0250 G4cout << " cal (l) : " << tpStore->calorimeterTPLeft << G4endl;
0251 G4cout << " cal (r) : " << tpStore->calorimeterTPRight << G4endl;
0252 G4cout << " ---" << G4endl;
0253 G4cout << " angle between the " <<
0254 tpStore->targetTPOutputParticle.particle->GetParticleName() <<
0255 " decay products : " <<
0256 tpStore->targetTPOutputParticleDecayProductParticle1.directionWorld.
0257 angle(
0258 tpStore->targetTPOutputParticleDecayProductParticle2.directionWorld ) /
0259 deg << " deg" << G4endl;
0260 }
0261
0262
0263 void CexmcEventAction::PrintProductionModelData(
0264 const CexmcAngularRangeList & angularRanges,
0265 const CexmcProductionModelData & pmData )
0266 {
0267 G4cout << " --- Triggered angular ranges: " << angularRanges;
0268 G4cout << " --- Production model data: " << pmData;
0269 }
0270
0271
0272 void CexmcEventAction::PrintReconstructedData(
0273 const CexmcAngularRangeList & triggeredRecAngularRanges,
0274 const CexmcAngularRange & angularGap ) const
0275 {
0276 G4cout << " --- Reconstructed data: " << G4endl;
0277 G4cout << " -- entry points:" << G4endl;
0278 G4cout << " left: " << G4BestUnit(
0279 reconstructor->GetCalorimeterEPLeftPosition(), "Length" ) << G4endl;
0280 G4cout << " right: " << G4BestUnit(
0281 reconstructor->GetCalorimeterEPRightPosition(), "Length" ) << G4endl;
0282 G4cout << " target: " << G4BestUnit(
0283 reconstructor->GetTargetEPPosition(), "Length" ) << G4endl;
0284 G4cout << " -- the angle: " << reconstructor->GetTheAngle() / deg <<
0285 " deg" << G4endl;
0286 G4cout << " -- mass of the output particle: " << G4BestUnit(
0287 reconstructor->GetOutputParticleMass(), "Energy" ) << G4endl;
0288 G4cout << " -- mass of the nucleus output particle: " << G4BestUnit(
0289 reconstructor->GetNucleusOutputParticleMass(), "Energy" ) << G4endl;
0290 if ( reconstructor->IsMassCutUsed() )
0291 {
0292 if ( reconstructor->HasMassCutTriggered() )
0293 G4cout << " < mass cut passed >" << G4endl;
0294 else
0295 G4cout << " < mass cut failed >" << G4endl;
0296 }
0297 if ( reconstructor->IsAbsorbedEnergyCutUsed() )
0298 {
0299 if ( reconstructor->HasAbsorbedEnergyCutTriggered() )
0300 G4cout << " < absorbed energy cut passed >" << G4endl;
0301 else
0302 G4cout << " < absorbed energy cut failed >" << G4endl;
0303 }
0304 const CexmcProductionModelData & pmData(
0305 reconstructor->GetProductionModelData() );
0306 G4cout << " -- production model data: " << pmData;
0307 G4cout << " -- triggered angular ranges: ";
0308 if ( triggeredRecAngularRanges.empty() )
0309 G4cout << "< orphan detected, gap " << angularGap << " >" << G4endl;
0310 else
0311 G4cout << triggeredRecAngularRanges;
0312 }
0313
0314
0315 #ifdef CEXMC_USE_ROOT
0316
0317 void CexmcEventAction::FillEDTHistos( const CexmcEnergyDepositStore * edStore,
0318 const CexmcAngularRangeList & triggeredAngularRanges ) const
0319 {
0320 CexmcHistoManager * histoManager( CexmcHistoManager::Instance() );
0321
0322 histoManager->Add( CexmcAbsorbedEnergy_EDT_Histo, 0,
0323 edStore->calorimeterEDLeft,
0324 edStore->calorimeterEDRight );
0325
0326 for ( CexmcAngularRangeList::const_iterator
0327 k( triggeredAngularRanges.begin() );
0328 k != triggeredAngularRanges.end(); ++k )
0329 {
0330 histoManager->Add( CexmcAbsEnInLeftCalorimeter_ARReal_EDT_Histo,
0331 k->index, edStore->calorimeterEDLeft );
0332 histoManager->Add( CexmcAbsEnInRightCalorimeter_ARReal_EDT_Histo,
0333 k->index, edStore->calorimeterEDRight );
0334 }
0335 }
0336
0337
0338 void CexmcEventAction::FillTPTHistos( const CexmcTrackPointsStore * tpStore,
0339 const CexmcProductionModelData & pmData,
0340 const CexmcAngularRangeList & triggeredAngularRanges ) const
0341 {
0342 CexmcHistoManager * histoManager( CexmcHistoManager::Instance() );
0343
0344 if ( tpStore->monitorTP.IsValid() )
0345 {
0346 histoManager->Add( CexmcMomentumBP_TPT_Histo, 0,
0347 tpStore->monitorTP.momentumAmp );
0348 histoManager->Add( CexmcTPInMonitor_TPT_Histo, 0,
0349 tpStore->monitorTP.positionLocal.x(),
0350 tpStore->monitorTP.positionLocal.y() );
0351 }
0352
0353 if ( tpStore->targetTPOutputParticle.IsValid() )
0354 {
0355 histoManager->Add( CexmcTPInTarget_TPT_Histo, 0,
0356 tpStore->targetTPOutputParticle.positionLocal.x(),
0357 tpStore->targetTPOutputParticle.positionLocal.y(),
0358 tpStore->targetTPOutputParticle.positionLocal.z() );
0359 if ( histoManager->GetVerboseLevel() > 0 )
0360 {
0361 histoManager->Add( CexmcMomentumIP_TPT_Histo, 0,
0362 pmData.incidentParticleLAB.rho() );
0363 }
0364 }
0365
0366 for ( CexmcAngularRangeList::const_iterator
0367 k( triggeredAngularRanges.begin() );
0368 k != triggeredAngularRanges.end(); ++k )
0369 {
0370 if ( tpStore->calorimeterTPLeft.IsValid() )
0371 {
0372
0373 G4double kinEnergy( tpStore->calorimeterTPLeft.momentumAmp );
0374 histoManager->Add( CexmcKinEnAtLeftCalorimeter_ARReal_TPT_Histo,
0375 k->index, kinEnergy );
0376 }
0377 if ( tpStore->calorimeterTPRight.IsValid() )
0378 {
0379 G4double kinEnergy( tpStore->calorimeterTPRight.momentumAmp );
0380 histoManager->Add( CexmcKinEnAtRightCalorimeter_ARReal_TPT_Histo,
0381 k->index, kinEnergy );
0382 }
0383 if ( tpStore->targetTPOutputParticle.IsValid() )
0384 {
0385 histoManager->Add( CexmcTPInTarget_ARReal_TPT_Histo, k->index,
0386 tpStore->targetTPOutputParticle.positionLocal.x(),
0387 tpStore->targetTPOutputParticle.positionLocal.y(),
0388 tpStore->targetTPOutputParticle.positionLocal.z() );
0389 histoManager->Add( CexmcKinEnOP_LAB_ARReal_TPT_Histo, k->index,
0390 opKinEnergy );
0391 histoManager->Add( CexmcAngleOP_SCM_ARReal_TPT_Histo, k->index,
0392 pmData.outputParticleSCM.cosTheta() );
0393 }
0394 if ( tpStore->targetTPOutputParticleDecayProductParticle1.IsValid() &&
0395 tpStore->targetTPOutputParticleDecayProductParticle2.IsValid() )
0396 {
0397 G4double openAngle(
0398 tpStore->targetTPOutputParticleDecayProductParticle1.
0399 directionWorld.angle( tpStore->
0400 targetTPOutputParticleDecayProductParticle2.
0401 directionWorld ) / deg );
0402 histoManager->Add( CexmcOpenAngle_ARReal_TPT_Histo, k->index,
0403 openAngle );
0404 }
0405 }
0406 }
0407
0408
0409 void CexmcEventAction::FillRTHistos( G4bool reconstructorHasFullTrigger,
0410 const CexmcEnergyDepositStore * edStore,
0411 const CexmcTrackPointsStore * tpStore,
0412 const CexmcProductionModelData & pmData,
0413 const CexmcAngularRangeList & triggeredAngularRanges ) const
0414 {
0415 CexmcHistoManager * histoManager( CexmcHistoManager::Instance() );
0416
0417 G4double opMass( reconstructor->GetOutputParticleMass() );
0418 G4double nopMass( reconstructor->GetNucleusOutputParticleMass() );
0419
0420 histoManager->Add( CexmcRecMasses_EDT_Histo, 0, opMass, nopMass );
0421
0422 for ( CexmcAngularRangeList::const_iterator
0423 k( triggeredAngularRanges.begin() );
0424 k != triggeredAngularRanges.end(); ++k )
0425 {
0426 if ( tpStore->calorimeterTPLeft.IsValid() )
0427 {
0428 histoManager->Add( CexmcOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
0429 k->index,
0430 tpStore->calorimeterTPLeft.positionLocal.x(),
0431 tpStore->calorimeterTPLeft.positionLocal.y() );
0432 }
0433 if ( tpStore->calorimeterTPRight.IsValid() )
0434 {
0435 histoManager->Add( CexmcOPDPAtRightCalorimeter_ARReal_EDT_Histo,
0436 k->index,
0437 tpStore->calorimeterTPRight.positionLocal.x(),
0438 tpStore->calorimeterTPRight.positionLocal.y() );
0439 }
0440 histoManager->Add( CexmcRecOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
0441 k->index,
0442 reconstructor->GetCalorimeterEPLeftPosition().x(),
0443 reconstructor->GetCalorimeterEPLeftPosition().y() );
0444 histoManager->Add( CexmcRecOPDPAtRightCalorimeter_ARReal_EDT_Histo,
0445 k->index,
0446 reconstructor->GetCalorimeterEPRightPosition().x(),
0447 reconstructor->GetCalorimeterEPRightPosition().y() );
0448 }
0449
0450 if ( ! reconstructorHasFullTrigger )
0451 return;
0452
0453 if ( tpStore->monitorTP.IsValid() )
0454 {
0455 histoManager->Add( CexmcMomentumBP_RT_Histo, 0,
0456 tpStore->monitorTP.momentumAmp );
0457 }
0458
0459 if ( tpStore->targetTPOutputParticle.IsValid() )
0460 {
0461 histoManager->Add( CexmcTPInTarget_RT_Histo, 0,
0462 tpStore->targetTPOutputParticle.positionLocal.x(),
0463 tpStore->targetTPOutputParticle.positionLocal.y(),
0464 tpStore->targetTPOutputParticle.positionLocal.z() );
0465 }
0466
0467 histoManager->Add( CexmcRecMasses_RT_Histo, 0,
0468 reconstructor->GetOutputParticleMass(),
0469 reconstructor->GetNucleusOutputParticleMass() );
0470
0471 histoManager->Add( CexmcAbsorbedEnergy_RT_Histo, 0,
0472 edStore->calorimeterEDLeft,
0473 edStore->calorimeterEDRight );
0474
0475 G4double recCosTheta( reconstructor->GetProductionModelData().
0476 outputParticleSCM.cosTheta());
0477
0478 for ( CexmcAngularRangeList::const_iterator
0479 k( triggeredAngularRanges.begin() );
0480 k != triggeredAngularRanges.end(); ++k )
0481 {
0482 histoManager->Add( CexmcRecMassOP_ARReal_RT_Histo, k->index, opMass );
0483 histoManager->Add( CexmcRecMassNOP_ARReal_RT_Histo, k->index, nopMass );
0484 if ( tpStore->calorimeterTPLeft.IsValid() )
0485 {
0486 G4double kinEnergy( tpStore->calorimeterTPLeft.momentumAmp );
0487 histoManager->Add( CexmcKinEnAtLeftCalorimeter_ARReal_RT_Histo,
0488 k->index, kinEnergy );
0489 histoManager->Add( CexmcMissEnFromLeftCalorimeter_ARReal_RT_Histo,
0490 k->index,
0491 kinEnergy - edStore->calorimeterEDLeft );
0492 }
0493 if ( tpStore->calorimeterTPRight.IsValid() )
0494 {
0495 G4double kinEnergy( tpStore->calorimeterTPRight.momentumAmp );
0496 histoManager->Add( CexmcKinEnAtRightCalorimeter_ARReal_RT_Histo,
0497 k->index, kinEnergy );
0498 histoManager->Add( CexmcMissEnFromRightCalorimeter_ARReal_RT_Histo,
0499 k->index,
0500 kinEnergy - edStore->calorimeterEDRight );
0501 }
0502 if ( tpStore->targetTPOutputParticle.IsValid() )
0503 {
0504 histoManager->Add( CexmcTPInTarget_ARReal_RT_Histo, k->index,
0505 tpStore->targetTPOutputParticle.positionLocal.x(),
0506 tpStore->targetTPOutputParticle.positionLocal.y(),
0507 tpStore->targetTPOutputParticle.positionLocal.z() );
0508 histoManager->Add( CexmcKinEnOP_LAB_ARReal_RT_Histo, k->index,
0509 opKinEnergy );
0510 histoManager->Add( CexmcAngleOP_SCM_ARReal_RT_Histo, k->index,
0511 pmData.outputParticleSCM.cosTheta() );
0512 G4double diffCosTheta( pmData.outputParticleSCM.cosTheta() -
0513 recCosTheta );
0514 histoManager->Add( CexmcDiffAngleOP_SCM_ARReal_RT_Histo, k->index,
0515 diffCosTheta );
0516 }
0517 if ( tpStore->targetTPOutputParticleDecayProductParticle1.IsValid() &&
0518 tpStore->targetTPOutputParticleDecayProductParticle2.IsValid() )
0519 {
0520 G4double openAngle(
0521 tpStore->targetTPOutputParticleDecayProductParticle1.
0522 directionWorld.angle( tpStore->
0523 targetTPOutputParticleDecayProductParticle2.
0524 directionWorld ) / deg );
0525 histoManager->Add( CexmcOpenAngle_ARReal_RT_Histo, k->index,
0526 openAngle );
0527 G4double diffOpenAngle( openAngle - reconstructor->GetTheAngle() /
0528 deg );
0529 histoManager->Add( CexmcDiffOpenAngle_ARReal_RT_Histo, k->index,
0530 diffOpenAngle );
0531 }
0532 if ( tpStore->calorimeterTPLeft.IsValid() )
0533 {
0534 histoManager->Add( CexmcOPDPAtLeftCalorimeter_ARReal_RT_Histo,
0535 k->index,
0536 tpStore->calorimeterTPLeft.positionLocal.x(),
0537 tpStore->calorimeterTPLeft.positionLocal.y() );
0538 }
0539 if ( tpStore->calorimeterTPRight.IsValid() )
0540 {
0541 histoManager->Add( CexmcOPDPAtRightCalorimeter_ARReal_RT_Histo,
0542 k->index,
0543 tpStore->calorimeterTPRight.positionLocal.x(),
0544 tpStore->calorimeterTPRight.positionLocal.y() );
0545 }
0546 histoManager->Add( CexmcAbsEnInLeftCalorimeter_ARReal_RT_Histo,
0547 k->index, edStore->calorimeterEDLeft );
0548 histoManager->Add( CexmcAbsEnInRightCalorimeter_ARReal_RT_Histo,
0549 k->index, edStore->calorimeterEDRight );
0550 histoManager->Add( CexmcRecAngleOP_SCM_ARReal_RT_Histo,
0551 k->index, recCosTheta );
0552 histoManager->Add( CexmcRecOpenAngle_ARReal_RT_Histo,
0553 k->index, reconstructor->GetTheAngle() / deg );
0554 histoManager->Add( CexmcRecOPDPAtLeftCalorimeter_ARReal_RT_Histo,
0555 k->index,
0556 reconstructor->GetCalorimeterEPLeftPosition().x(),
0557 reconstructor->GetCalorimeterEPLeftPosition().y() );
0558 histoManager->Add( CexmcRecOPDPAtRightCalorimeter_ARReal_RT_Histo,
0559 k->index,
0560 reconstructor->GetCalorimeterEPRightPosition().x(),
0561 reconstructor->GetCalorimeterEPRightPosition().y() );
0562 }
0563 }
0564
0565 #endif
0566
0567
0568 void CexmcEventAction::DrawTrajectories( const G4Event * event )
0569 {
0570 G4VisManager * visManager( static_cast< G4VisManager * >(
0571 G4VVisManager::GetConcreteInstance() ) );
0572 if ( ! visManager || ! visManager->GetCurrentGraphicsSystem() )
0573 return;
0574
0575 G4int nTraj( 0 );
0576 G4TrajectoryContainer * trajContainer( event->GetTrajectoryContainer() );
0577
0578 if ( ! trajContainer )
0579 return;
0580
0581 nTraj = trajContainer->entries();
0582
0583 for ( int i( 0 ); i < nTraj; ++i )
0584 {
0585 G4VTrajectory * traj( ( *trajContainer )[ i ] );
0586 traj->DrawTrajectory();
0587 }
0588 }
0589
0590
0591 void CexmcEventAction::DrawTrackPoints(
0592 const CexmcTrackPointsStore * tpStore ) const
0593 {
0594 G4VisManager * visManager( static_cast< G4VisManager * >(
0595 G4VVisManager::GetConcreteInstance() ) );
0596 if ( ! visManager || ! visManager->GetCurrentGraphicsSystem() )
0597 return;
0598
0599 G4Circle circle;
0600 G4VisAttributes visAttributes( CexmcTrackPointsMarkerColour );
0601 circle.SetScreenSize( CexmcSmallCircleScreenSize );
0602 circle.SetFillStyle( G4Circle::filled );
0603 circle.SetVisAttributes( visAttributes );
0604
0605 if ( tpStore->monitorTP.IsValid() )
0606 {
0607 circle.SetPosition( tpStore->monitorTP.positionWorld );
0608 visManager->Draw( circle );
0609 }
0610
0611 if ( tpStore->targetTPBeamParticle.IsValid() )
0612 {
0613 circle.SetPosition( tpStore->targetTPBeamParticle.positionWorld );
0614 visManager->Draw( circle );
0615 }
0616
0617 if ( tpStore->targetTPOutputParticle.IsValid() )
0618 {
0619 circle.SetPosition( tpStore->targetTPOutputParticle.positionWorld );
0620 visManager->Draw( circle );
0621 }
0622
0623 if ( tpStore->vetoCounterTPLeft.IsValid() )
0624 {
0625 circle.SetPosition( tpStore->vetoCounterTPLeft.positionWorld );
0626 visManager->Draw( circle );
0627 }
0628
0629 if ( tpStore->vetoCounterTPRight.IsValid() )
0630 {
0631 circle.SetPosition( tpStore->vetoCounterTPRight.positionWorld );
0632 visManager->Draw( circle );
0633 }
0634
0635 if ( tpStore->calorimeterTPLeft.IsValid() )
0636 {
0637 circle.SetPosition( tpStore->calorimeterTPLeft.positionWorld );
0638 visManager->Draw( circle );
0639 }
0640
0641 if ( tpStore->calorimeterTPRight.IsValid() )
0642 {
0643 circle.SetPosition( tpStore->calorimeterTPRight.positionWorld );
0644 visManager->Draw( circle );
0645 }
0646 }
0647
0648
0649 void CexmcEventAction::DrawReconstructionData( void )
0650 {
0651 G4VisManager * visManager( static_cast< G4VisManager * >(
0652 G4VVisManager::GetConcreteInstance() ) );
0653 if ( ! visManager || ! visManager->GetCurrentGraphicsSystem() )
0654 return;
0655
0656 G4Circle circle( reconstructor->GetTargetEPWorldPosition() );
0657 circle.SetScreenSize( CexmcSmallCircleScreenSize );
0658 circle.SetFillStyle( G4Circle::filled );
0659 G4VisAttributes visAttributes( CexmcRecTrackPointsMarkerColour );
0660 circle.SetVisAttributes( visAttributes );
0661 visManager->Draw( circle );
0662
0663 circle.SetScreenSize( CexmcBigCircleScreenSize );
0664 circle.SetPosition( reconstructor->GetCalorimeterEPLeftWorldPosition() );
0665 visManager->Draw( circle );
0666
0667 circle.SetPosition( reconstructor->GetCalorimeterEPRightWorldPosition() );
0668 visManager->Draw( circle );
0669 }
0670
0671
0672 void CexmcEventAction::UpdateRunHits(
0673 const CexmcAngularRangeList & aRangesReal,
0674 const CexmcAngularRangeList & aRangesRec,
0675 G4bool tpDigitizerHasTriggered,
0676 G4bool edDigitizerHasTriggered,
0677 G4bool edDigitizerMonitorHasTriggered,
0678 G4bool reconstructorHasFullTrigger,
0679 const CexmcAngularRange & aGap )
0680 {
0681 G4RunManager * runManager( G4RunManager::GetRunManager() );
0682 const CexmcRun * run( static_cast< const CexmcRun * >(
0683 runManager->GetCurrentRun() ) );
0684 CexmcRun * theRun( const_cast< CexmcRun * >( run ) );
0685
0686 if ( tpDigitizerHasTriggered )
0687 {
0688 for ( CexmcAngularRangeList::const_iterator k( aRangesReal.begin() );
0689 k != aRangesReal.end(); ++k )
0690 {
0691 theRun->IncrementNmbOfHitsSampledFull( k->index );
0692 if ( edDigitizerMonitorHasTriggered )
0693 theRun->IncrementNmbOfHitsSampled( k->index );
0694 if ( reconstructorHasFullTrigger )
0695 theRun->IncrementNmbOfHitsTriggeredRealRange( k->index );
0696 }
0697 if ( reconstructorHasFullTrigger )
0698 {
0699 if ( aRangesRec.empty() )
0700 {
0701 theRun->IncrementNmbOfOrphanHits( aGap.index );
0702 }
0703 else
0704 {
0705 for ( CexmcAngularRangeList::const_iterator
0706 k( aRangesRec.begin() ); k != aRangesRec.end(); ++k )
0707 {
0708 theRun->IncrementNmbOfHitsTriggeredRecRange( k->index );
0709 }
0710 }
0711 }
0712 }
0713 else
0714 {
0715 if ( edDigitizerHasTriggered )
0716 theRun->IncrementNmbOfFalseHitsTriggeredEDT();
0717 if ( reconstructorHasFullTrigger )
0718 theRun->IncrementNmbOfFalseHitsTriggeredRec();
0719 }
0720 }
0721
0722
0723 #ifdef CEXMC_USE_PERSISTENCY
0724
0725 void CexmcEventAction::SaveEvent( const G4Event * event,
0726 G4bool edDigitizerHasTriggered,
0727 const CexmcEnergyDepositStore * edStore,
0728 const CexmcTrackPointsStore * tpStore,
0729 const CexmcProductionModelData & pmData )
0730 {
0731 CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
0732 G4RunManager::GetRunManager() ) );
0733 if ( ! runManager->ProjectIsSaved() )
0734 return;
0735
0736 if ( runManager->GetEventDataVerboseLevel() == CexmcWriteNoEventData )
0737 return;
0738
0739 if ( ! edDigitizerHasTriggered && runManager->GetEventDataVerboseLevel() !=
0740 CexmcWriteEventDataOnEveryTPT )
0741 return;
0742
0743 boost::archive::binary_oarchive * archive(
0744 runManager->GetEventsArchive() );
0745 if ( archive )
0746 {
0747 CexmcEventSObject sObject = { event->GetEventID(),
0748 edDigitizerHasTriggered, edStore->monitorED,
0749 edStore->vetoCounterEDLeft, edStore->vetoCounterEDRight,
0750 edStore->calorimeterEDLeft, edStore->calorimeterEDRight,
0751 edStore->calorimeterEDLeftCollection,
0752 edStore->calorimeterEDRightCollection,
0753 tpStore->monitorTP, tpStore->targetTPBeamParticle,
0754 tpStore->targetTPOutputParticle, tpStore->targetTPNucleusParticle,
0755 tpStore->targetTPOutputParticleDecayProductParticle1,
0756 tpStore->targetTPOutputParticleDecayProductParticle2,
0757 tpStore->vetoCounterTPLeft, tpStore->vetoCounterTPRight,
0758 tpStore->calorimeterTPLeft, tpStore->calorimeterTPRight, pmData };
0759 archive->operator<<( sObject );
0760 const CexmcRun * run( static_cast< const CexmcRun * >(
0761 runManager->GetCurrentRun() ) );
0762 CexmcRun * theRun( const_cast< CexmcRun * >( run ) );
0763 theRun->IncrementNmbOfSavedEvents();
0764 }
0765 }
0766
0767
0768 void CexmcEventAction::SaveEventFast( const G4Event * event,
0769 G4bool tpDigitizerHasTriggered,
0770 G4bool edDigitizerHasTriggered,
0771 G4bool edDigitizerMonitorHasTriggered,
0772 G4double opCosThetaSCM )
0773 {
0774 CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
0775 G4RunManager::GetRunManager() ) );
0776 if ( ! runManager->ProjectIsSaved() )
0777 return;
0778
0779 if ( runManager->GetEventDataVerboseLevel() == CexmcWriteNoEventData )
0780 return;
0781
0782 boost::archive::binary_oarchive * archive(
0783 runManager->GetFastEventsArchive() );
0784 if ( archive )
0785 {
0786 if ( ! tpDigitizerHasTriggered )
0787 opCosThetaSCM = CexmcInvalidCosTheta;
0788
0789 CexmcEventFastSObject sObject = { event->GetEventID(), opCosThetaSCM,
0790 edDigitizerHasTriggered,
0791 edDigitizerMonitorHasTriggered };
0792 archive->operator<<( sObject );
0793 const CexmcRun * run( static_cast< const CexmcRun * >(
0794 runManager->GetCurrentRun() ) );
0795 CexmcRun * theRun( const_cast< CexmcRun * >( run ) );
0796 theRun->IncrementNmbOfSavedFastEvents();
0797 }
0798 }
0799
0800 #endif
0801
0802
0803 void CexmcEventAction::EndOfEventAction( const G4Event * event )
0804 {
0805 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
0806 CexmcEnergyDepositDigitizer * energyDepositDigitizer(
0807 static_cast< CexmcEnergyDepositDigitizer * >( digiManager->
0808 FindDigitizerModule( CexmcEDDigitizerName ) ) );
0809 CexmcTrackPointsDigitizer * trackPointsDigitizer(
0810 static_cast< CexmcTrackPointsDigitizer * >( digiManager->
0811 FindDigitizerModule( CexmcTPDigitizerName ) ) );
0812
0813 energyDepositDigitizer->Digitize();
0814 trackPointsDigitizer->Digitize();
0815
0816 G4bool edDigitizerMonitorHasTriggered(
0817 energyDepositDigitizer->MonitorHasTriggered() );
0818 G4bool edDigitizerHasTriggered( false );
0819
0820 CexmcEventInfo * eventInfo( static_cast< CexmcEventInfo * >(
0821 event->GetUserInformation() ) );
0822 if ( ! eventInfo || eventInfo->EdTriggerIsOk() )
0823 edDigitizerHasTriggered = energyDepositDigitizer->HasTriggered();
0824
0825 G4bool tpDigitizerHasTriggered( trackPointsDigitizer->HasTriggered() );
0826 G4bool reconstructorHasBasicTrigger( false );
0827 G4bool reconstructorHasFullTrigger( false );
0828
0829 CexmcEnergyDepositStore * edStore( MakeEnergyDepositStore(
0830 energyDepositDigitizer ) );
0831 CexmcTrackPointsStore * tpStore( MakeTrackPointsStore(
0832 trackPointsDigitizer ) );
0833
0834 try
0835 {
0836 CexmcProductionModel * productionModel(
0837 physicsManager->GetProductionModel() );
0838
0839 if ( ! productionModel )
0840 throw CexmcException( CexmcWeirdException );
0841
0842 const CexmcAngularRangeList & angularRanges(
0843 productionModel->GetAngularRanges() );
0844 const CexmcAngularRangeList & triggeredAngularRanges(
0845 productionModel->GetTriggeredAngularRanges() );
0846 const CexmcProductionModelData & pmData(
0847 productionModel->GetProductionModelData() );
0848
0849 if ( edDigitizerHasTriggered )
0850 {
0851 reconstructor->Reconstruct( edStore );
0852 reconstructorHasBasicTrigger = reconstructor->HasBasicTrigger();
0853 reconstructorHasFullTrigger = reconstructor->HasFullTrigger();
0854 }
0855
0856 CexmcAngularRangeList triggeredRecAngularRanges;
0857
0858 if ( reconstructorHasBasicTrigger )
0859 {
0860 for ( CexmcAngularRangeList::const_iterator
0861 k( angularRanges.begin() ); k != angularRanges.end(); ++k )
0862 {
0863 G4double cosTheta( reconstructor->GetProductionModelData().
0864 outputParticleSCM.cosTheta() );
0865 if ( cosTheta <= k->top && cosTheta > k->bottom )
0866 triggeredRecAngularRanges.push_back( CexmcAngularRange(
0867 k->top, k->bottom, k->index ) );
0868 }
0869 }
0870
0871 CexmcAngularRange angularGap( 0.0, 0.0, 0 );
0872 if ( triggeredRecAngularRanges.empty() )
0873 {
0874 CexmcAngularRangeList angularGaps;
0875 GetAngularGaps( angularRanges, angularGaps );
0876 for ( CexmcAngularRangeList::const_iterator
0877 k( angularGaps.begin() ); k != angularGaps.end(); ++k )
0878 {
0879 G4double cosTheta( reconstructor->GetProductionModelData().
0880 outputParticleSCM.cosTheta() );
0881 if ( cosTheta <= k->top && cosTheta > k->bottom )
0882 {
0883 angularGap = *k;
0884 break;
0885 }
0886 }
0887 }
0888
0889 UpdateRunHits( triggeredAngularRanges, triggeredRecAngularRanges,
0890 tpDigitizerHasTriggered, edDigitizerHasTriggered,
0891 edDigitizerMonitorHasTriggered,
0892 reconstructorHasFullTrigger, angularGap );
0893
0894 if ( verbose > 0 )
0895 {
0896 G4bool printMessages( verbose > 3 ||
0897 ( ( verbose == 1 ) && tpDigitizerHasTriggered ) ||
0898 ( ( verbose == 2 ) && edDigitizerHasTriggered ) ||
0899 ( ( verbose == 3 ) && ( tpDigitizerHasTriggered ||
0900 edDigitizerHasTriggered ) ) );
0901 if ( printMessages )
0902 {
0903 G4cout << "Event " << event->GetEventID() << G4endl;
0904 if ( tpDigitizerHasTriggered )
0905 {
0906 PrintTrackPoints( tpStore );
0907 PrintProductionModelData( triggeredAngularRanges, pmData );
0908 }
0909 if ( reconstructorHasBasicTrigger )
0910 PrintReconstructedData( triggeredRecAngularRanges,
0911 angularGap );
0912 if ( edDigitizerHasTriggered )
0913 PrintEnergyDeposit( edStore );
0914 }
0915 }
0916
0917 if ( verboseDraw > 0 )
0918 {
0919 G4bool drawTrajectories( verboseDraw > 3 ||
0920 ( ( verboseDraw == 1 ) && tpDigitizerHasTriggered ) ||
0921 ( ( verboseDraw == 2 ) && edDigitizerHasTriggered ) ||
0922 ( ( verboseDraw == 3 ) && ( tpDigitizerHasTriggered ||
0923 edDigitizerHasTriggered ) ) );
0924 if ( drawTrajectories )
0925 {
0926 DrawTrajectories( event );
0927 if ( tpDigitizerHasTriggered )
0928 DrawTrackPoints( tpStore );
0929 if ( reconstructorHasBasicTrigger )
0930 DrawReconstructionData();
0931 }
0932 }
0933
0934 #ifdef CEXMC_USE_PERSISTENCY
0935 if ( edDigitizerHasTriggered || tpDigitizerHasTriggered )
0936 {
0937 SaveEventFast( event, tpDigitizerHasTriggered,
0938 edDigitizerHasTriggered,
0939 edDigitizerMonitorHasTriggered,
0940 pmData.outputParticleSCM.cosTheta() );
0941 SaveEvent( event, edDigitizerHasTriggered, edStore, tpStore,
0942 pmData );
0943 }
0944 #endif
0945
0946 #ifdef CEXMC_USE_ROOT
0947
0948 if ( tpStore->targetTPOutputParticle.IsValid() )
0949 {
0950 opKinEnergy = CexmcGetKinEnergy(
0951 tpStore->targetTPOutputParticle.momentumAmp,
0952 tpStore->targetTPOutputParticle.particle->GetPDGMass() );
0953 }
0954
0955 if ( edDigitizerHasTriggered )
0956 FillEDTHistos( edStore, triggeredAngularRanges );
0957
0958
0959
0960 if ( tpDigitizerHasTriggered && edDigitizerMonitorHasTriggered )
0961 FillTPTHistos( tpStore, pmData, triggeredAngularRanges );
0962
0963 if ( reconstructorHasBasicTrigger )
0964 FillRTHistos( reconstructorHasFullTrigger, edStore, tpStore,
0965 pmData, triggeredAngularRanges );
0966 #endif
0967
0968 G4Event * theEvent( const_cast< G4Event * >( event ) );
0969 if ( eventInfo )
0970 {
0971 delete eventInfo;
0972 theEvent->SetUserInformation( NULL );
0973 }
0974 theEvent->SetUserInformation( new CexmcEventInfo(
0975 edDigitizerHasTriggered,
0976 tpDigitizerHasTriggered,
0977 reconstructorHasFullTrigger ) );
0978 }
0979 catch ( CexmcException & e )
0980 {
0981 G4cout << e.what() << G4endl;
0982 }
0983 catch ( ... )
0984 {
0985 G4cout << "Unknown exception caught" << G4endl;
0986 }
0987
0988 delete edStore;
0989 delete tpStore;
0990 }
0991
