File indexing completed on 2025-01-31 09:21:54
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 <stdlib.h>
0045 #include <sys/stat.h>
0046 #include <vector>
0047 #include <fstream>
0048 #ifdef CEXMC_USE_PERSISTENCY
0049 #include <boost/archive/binary_iarchive.hpp>
0050 #include <boost/archive/binary_oarchive.hpp>
0051 #endif
0052 #include <G4Eta.hh>
0053 #include <G4DigiManager.hh>
0054 #include <G4SDManager.hh>
0055 #include <G4DecayTable.hh>
0056 #include <G4ParticleDefinition.hh>
0057 #include <G4ParticleTable.hh>
0058 #include <G4UImanager.hh>
0059 #include <G4Timer.hh>
0060 #include <G4Region.hh>
0061 #include <G4RegionStore.hh>
0062 #include <G4ProductionCuts.hh>
0063 #include <G4VisManager.hh>
0064 #include <G4Scene.hh>
0065 #include <G4VModel.hh>
0066 #include <G4Version.hh>
0067 #include "CexmcRunManager.hh"
0068 #include "CexmcRunManagerMessenger.hh"
0069 #include "CexmcRunAction.hh"
0070 #include "CexmcRun.hh"
0071 #include "CexmcPhysicsManager.hh"
0072 #include "CexmcProductionModel.hh"
0073 #include "CexmcSimpleDecayTableStore.hh"
0074 #include "CexmcPrimaryGeneratorAction.hh"
0075 #include "CexmcEnergyDepositDigitizer.hh"
0076 #include "CexmcSimpleEnergyDeposit.hh"
0077 #include "CexmcEnergyDepositInLeftRightSet.hh"
0078 #include "CexmcEnergyDepositInCalorimeter.hh"
0079 #include "CexmcTrackPoints.hh"
0080 #include "CexmcTrackPointsInLeftRightSet.hh"
0081 #include "CexmcTrackPointsInCalorimeter.hh"
0082 #include "CexmcChargeExchangeReconstructor.hh"
0083 #include "CexmcEventAction.hh"
0084 #include "CexmcParticleGun.hh"
0085 #include "CexmcEnergyDepositStore.hh"
0086 #include "CexmcTrackPointsStore.hh"
0087 #include "CexmcSetup.hh"
0088 #include "CexmcEventSObject.hh"
0089 #include "CexmcEventFastSObject.hh"
0090 #include "CexmcTrackPointInfo.hh"
0091 #include "CexmcEventInfo.hh"
0092 #include "CexmcBasicPhysicsSettings.hh"
0093 #include "CexmcSensitiveDetectorsAttributes.hh"
0094 #include "CexmcCustomFilterEval.hh"
0095 #include "CexmcScenePrimitives.hh"
0096
0097
0098 namespace
0099 {
0100 G4String gdmlFileExtension( ".gdml" );
0101 G4String gdmlbz2FileExtension( ".gdml.bz2" );
0102 }
0103
0104
0105 CexmcRunManager::CexmcRunManager( const G4String & projectId_,
0106 const G4String & rProject_,
0107 G4bool overrideExistingProject ) :
0108 basePhysicsUsed( CexmcPMFactoryInstance::GetBasePhysics() ),
0109 productionModelType( CexmcUnknownProductionModel ),
0110 gdmlFileName( "default.gdml" ), shouldGdmlFileBeValidated( false ),
0111 zipGdmlFile( false ), projectsDir( "." ), projectId( projectId_ ),
0112 rProject( rProject_ ), guiMacroName( "" ), cfFileName( "" ),
0113 eventCountPolicy( CexmcCountAllEvents ),
0114 skipInteractionsWithoutEDTonWrite( true ),
0115 evDataVerboseLevel( CexmcWriteEventDataOnEveryEDT ),
0116 rEvDataVerboseLevel( CexmcWriteNoEventData ), numberOfEventsProcessed( 0 ),
0117 numberOfEventsProcessedEffective( 0 ), curEventRead( 0 ),
0118 #ifdef CEXMC_USE_PERSISTENCY
0119 eventsArchive( NULL ), fastEventsArchive( NULL ),
0120 #ifdef CEXMC_USE_CUSTOM_FILTER
0121 customFilter( NULL ),
0122 #endif
0123 #endif
0124 physicsManager( NULL ), messenger( NULL )
0125 {
0126
0127 if ( rProject != "" && rProject == projectId )
0128 throw CexmcException( CexmcWeirdException );
0129
0130 const char * projectsDirEnv( std::getenv( "CEXMC_PROJECTS_DIR" ) );
0131
0132 if ( projectsDirEnv )
0133 projectsDir = projectsDirEnv;
0134
0135 struct stat tmp;
0136 if ( ProjectIsSaved() &&
0137 stat( ( projectsDir + "/" + projectId + ".rdb" ).c_str(), &tmp ) == 0
0138 && ! overrideExistingProject )
0139 throw CexmcException( CexmcProjectExists );
0140
0141 messenger = new CexmcRunManagerMessenger( this );
0142
0143 #ifdef CEXMC_USE_PERSISTENCY
0144 if ( ProjectIsRead() )
0145 ReadPreinitProjectData();
0146 #endif
0147 }
0148
0149
0150 CexmcRunManager::~CexmcRunManager()
0151 {
0152 #ifdef CEXMC_USE_CUSTOM_FILTER
0153 delete customFilter;
0154 #endif
0155 delete messenger;
0156 }
0157
0158
0159 #ifdef CEXMC_USE_PERSISTENCY
0160
0161 void CexmcRunManager::ReadPreinitProjectData( void )
0162 {
0163 if ( ! ProjectIsRead() )
0164 return;
0165
0166
0167 std::ifstream runDataFile( ( projectsDir + "/" + rProject + ".rdb" ).
0168 c_str() );
0169 if ( ! runDataFile )
0170 throw CexmcException( CexmcReadProjectIncomplete );
0171
0172 {
0173 boost::archive::binary_iarchive archive( runDataFile );
0174 archive >> sObject;
0175 }
0176
0177 basePhysicsUsed = sObject.basePhysicsUsed;
0178
0179 productionModelType = sObject.productionModelType;
0180
0181
0182 G4String cmd;
0183 if ( ProjectIsSaved() )
0184 {
0185 G4String fileExtension( zipGdmlFile ? gdmlbz2FileExtension :
0186 gdmlFileExtension );
0187 cmd = G4String( "cp " ) + projectsDir + "/" + rProject + fileExtension +
0188 " " + projectsDir + "/" + projectId + fileExtension;
0189 if ( system( cmd ) != 0 )
0190 throw CexmcException( CexmcReadProjectIncomplete );
0191 }
0192
0193 if ( zipGdmlFile )
0194 {
0195 cmd = G4String( "bunzip2 " ) + projectsDir + "/" + rProject +
0196 gdmlbz2FileExtension;
0197 if ( system( cmd ) != 0 )
0198 throw CexmcException( CexmcFileCompressException );
0199 }
0200
0201 gdmlFileName = projectsDir + "/" + rProject + gdmlFileExtension;
0202 }
0203
0204
0205 void CexmcRunManager::ReadProject( void )
0206 {
0207 if ( ! ProjectIsRead() )
0208 return;
0209
0210 if ( ! physicsManager )
0211 throw CexmcException( CexmcWeirdException );
0212
0213 physicsManager->GetProductionModel()->SetAngularRanges(
0214 sObject.angularRanges );
0215 G4DecayTable * etaDecayTable( G4Eta::Definition()->GetDecayTable() );
0216 for ( CexmcDecayBranchesStore::const_iterator
0217 k( sObject.etaDecayTable.GetDecayBranches().begin() );
0218 k != sObject.etaDecayTable.GetDecayBranches().end(); ++k )
0219 {
0220 etaDecayTable->GetDecayChannel( k->first )->SetBR( k->second );
0221 }
0222
0223 physicsManager->GetProductionModel()->ApplyFermiMotion(
0224 sObject.fermiMotionIsOn, false );
0225 eventCountPolicy = sObject.eventCountPolicy;
0226
0227 G4Region * region( G4RegionStore::GetInstance()->GetRegion(
0228 CexmcCalorimeterRegionName ) );
0229 if ( ! region )
0230 throw CexmcException( CexmcCalorimeterRegionNotInitialized );
0231
0232 region->GetProductionCuts()->SetProductionCuts(
0233 sObject.calorimeterRegCuts );
0234
0235 const CexmcPrimaryGeneratorAction * primaryGeneratorAction(
0236 static_cast< const CexmcPrimaryGeneratorAction * >(
0237 userPrimaryGeneratorAction ) );
0238 CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction(
0239 const_cast< CexmcPrimaryGeneratorAction * >(
0240 primaryGeneratorAction ) );
0241 CexmcParticleGun * particleGun(
0242 thePrimaryGeneratorAction->GetParticleGun() );
0243 G4ParticleDefinition * particleDefinition(
0244 G4ParticleTable::GetParticleTable()->FindParticle(
0245 sObject.beamParticle ) );
0246 if ( ! particleDefinition )
0247 throw CexmcException( CexmcWeirdException );
0248
0249 particleGun->SetParticleDefinition( particleDefinition );
0250 particleGun->SetOrigPosition( sObject.beamPos, false );
0251 particleGun->SetOrigDirection( sObject.beamDir, false );
0252 particleGun->SetOrigMomentumAmp( sObject.beamMomentumAmp, false );
0253
0254 thePrimaryGeneratorAction->SetFwhmPosX( sObject.beamFwhmPosX, false );
0255 thePrimaryGeneratorAction->SetFwhmPosY( sObject.beamFwhmPosY, false );
0256 thePrimaryGeneratorAction->SetFwhmDirX( sObject.beamFwhmDirX, false );
0257 thePrimaryGeneratorAction->SetFwhmDirY( sObject.beamFwhmDirY, false );
0258 thePrimaryGeneratorAction->SetFwhmMomentumAmp(
0259 sObject.beamFwhmMomentumAmp, false );
0260
0261 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
0262 CexmcEnergyDepositDigitizer * edDigitizer(
0263 static_cast< CexmcEnergyDepositDigitizer * >(
0264 digiManager->FindDigitizerModule( CexmcEDDigitizerName ) ) );
0265 if ( ! edDigitizer )
0266 throw CexmcException( CexmcWeirdException );
0267
0268 edDigitizer->SetMonitorThreshold( sObject.monitorEDThreshold, false );
0269 edDigitizer->SetVetoCounterLeftThreshold(
0270 sObject.vetoCounterEDLeftThreshold, false );
0271 edDigitizer->SetVetoCounterRightThreshold(
0272 sObject.vetoCounterEDRightThreshold, false );
0273 edDigitizer->SetCalorimeterLeftThreshold(
0274 sObject.calorimeterEDLeftThreshold, false );
0275 edDigitizer->SetCalorimeterRightThreshold(
0276 sObject.calorimeterEDRightThreshold, false );
0277 edDigitizer->SetCalorimeterTriggerAlgorithm(
0278 sObject.calorimeterTriggerAlgorithm, false );
0279 edDigitizer->SetOuterCrystalsVetoAlgorithm(
0280 sObject.outerCrystalsVetoAlgorithm, false );
0281 edDigitizer->SetOuterCrystalsVetoFraction(
0282 sObject.outerCrystalsVetoFraction, false );
0283 edDigitizer->ApplyFiniteCrystalResolution(
0284 sObject.applyFiniteCrystalResolution, false );
0285 edDigitizer->SetCrystalResolutionData( sObject.crystalResolutionData );
0286
0287 const CexmcEventAction * eventAction(
0288 static_cast< const CexmcEventAction * >( userEventAction ) );
0289 if ( ! eventAction )
0290 throw CexmcException( CexmcWeirdException );
0291
0292 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
0293 eventAction ) );
0294 CexmcChargeExchangeReconstructor * reconstructor(
0295 theEventAction->GetReconstructor() );
0296 if ( ! reconstructor )
0297 throw CexmcException( CexmcWeirdException );
0298
0299 reconstructor->SetCalorimeterEntryPointDefinitionAlgorithm(
0300 sObject.epDefinitionAlgorithm );
0301 reconstructor->SetCalorimeterEntryPointDepthDefinitionAlgorithm(
0302 sObject.epDepthDefinitionAlgorithm );
0303 reconstructor->SetCrystalSelectionAlgorithm( sObject.csAlgorithm );
0304 reconstructor->UseInnerRefCrystal( sObject.useInnerRefCrystal );
0305 reconstructor->SetCalorimeterEntryPointDepth( sObject.epDepth );
0306 reconstructor->UseTableMass( sObject.useTableMass );
0307 reconstructor->UseMassCut( sObject.useMassCut );
0308 reconstructor->SetMassCutOPCenter( sObject.mCutOPCenter );
0309 reconstructor->SetMassCutNOPCenter( sObject.mCutNOPCenter );
0310 reconstructor->SetMassCutOPWidth( sObject.mCutOPWidth );
0311 reconstructor->SetMassCutNOPWidth( sObject.mCutNOPWidth );
0312 reconstructor->SetMassCutEllipseAngle( sObject.mCutAngle );
0313 reconstructor->UseAbsorbedEnergyCut( sObject.useAbsorbedEnergyCut );
0314 reconstructor->SetAbsorbedEnergyCutCLCenter( sObject.aeCutCLCenter );
0315 reconstructor->SetAbsorbedEnergyCutCRCenter( sObject.aeCutCRCenter );
0316 reconstructor->SetAbsorbedEnergyCutCLWidth( sObject.aeCutCLWidth );
0317 reconstructor->SetAbsorbedEnergyCutCRWidth( sObject.aeCutCRWidth );
0318 reconstructor->SetAbsorbedEnergyCutEllipseAngle( sObject.aeCutAngle );
0319 reconstructor->SetExpectedMomentumAmp( sObject.expectedMomentumAmp );
0320 reconstructor->SetEDCollectionAlgorithm( sObject.edCollectionAlgorithm );
0321
0322 physicsManager->SetProposedMaxIL( sObject.proposedMaxIL );
0323
0324 rEvDataVerboseLevel = sObject.evDataVerboseLevel;
0325 evDataVerboseLevel = rEvDataVerboseLevel;
0326 }
0327
0328
0329 void CexmcRunManager::SaveProject( void )
0330 {
0331 if ( ! ProjectIsSaved() )
0332 return;
0333
0334
0335 if ( ! physicsManager )
0336 throw CexmcException( CexmcWeirdException );
0337
0338 CexmcSimpleDecayTableStore etaDecayTable(
0339 G4Eta::Definition()->GetDecayTable() );
0340 const CexmcPrimaryGeneratorAction * primaryGeneratorAction(
0341 static_cast< const CexmcPrimaryGeneratorAction * >(
0342 userPrimaryGeneratorAction ) );
0343 CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction(
0344 const_cast< CexmcPrimaryGeneratorAction * >(
0345 primaryGeneratorAction ) );
0346 CexmcParticleGun * particleGun(
0347 thePrimaryGeneratorAction->GetParticleGun() );
0348
0349 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
0350 CexmcEnergyDepositDigitizer * edDigitizer(
0351 static_cast< CexmcEnergyDepositDigitizer * >(
0352 digiManager->FindDigitizerModule( CexmcEDDigitizerName ) ) );
0353 if ( ! edDigitizer )
0354 throw CexmcException( CexmcWeirdException );
0355
0356 const CexmcEventAction * eventAction(
0357 static_cast< const CexmcEventAction * >( userEventAction ) );
0358 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
0359 eventAction ) );
0360 CexmcChargeExchangeReconstructor * reconstructor(
0361 theEventAction->GetReconstructor() );
0362
0363 std::vector< G4double > calorimeterRegCuts( 4 );
0364 if ( ProjectIsRead() )
0365 {
0366 calorimeterRegCuts = sObject.calorimeterRegCuts;
0367 }
0368 else
0369 {
0370 G4Region * region( G4RegionStore::GetInstance()->GetRegion(
0371 CexmcCalorimeterRegionName ) );
0372 if ( ! region )
0373 throw CexmcException( CexmcCalorimeterRegionNotInitialized );
0374
0375 calorimeterRegCuts = region->GetProductionCuts()->GetProductionCuts();
0376 }
0377
0378 CexmcNmbOfHitsInRanges nmbOfHitsSampled;
0379 CexmcNmbOfHitsInRanges nmbOfHitsSampledFull;
0380 CexmcNmbOfHitsInRanges nmbOfHitsTriggeredRealRange;
0381 CexmcNmbOfHitsInRanges nmbOfHitsTriggeredRecRange;
0382 CexmcNmbOfHitsInRanges nmbOfOrphanHits;
0383 G4int nmbOfFalseHitsTriggeredEDT( 0 );
0384 G4int nmbOfFalseHitsTriggeredRec( 0 );
0385 G4int nmbOfSavedEvents( 0 );
0386 G4int nmbOfSavedFastEvents( 0 );
0387 CexmcRun * run( static_cast< CexmcRun * >( currentRun ) );
0388
0389 if ( run )
0390 {
0391 nmbOfHitsSampled = run->GetNmbOfHitsSampled();
0392 nmbOfHitsSampledFull = run->GetNmbOfHitsSampledFull();
0393 nmbOfHitsTriggeredRealRange = run->GetNmbOfHitsTriggeredRealRange();
0394 nmbOfHitsTriggeredRecRange = run->GetNmbOfHitsTriggeredRecRange();
0395 nmbOfOrphanHits = run->GetNmbOfOrphanHits();
0396 nmbOfFalseHitsTriggeredEDT = run->GetNmbOfFalseHitsTriggeredEDT();
0397 nmbOfFalseHitsTriggeredRec = run->GetNmbOfFalseHitsTriggeredRec();
0398 nmbOfSavedEvents = run->GetNmbOfSavedEvents();
0399 nmbOfSavedFastEvents = run->GetNmbOfSavedFastEvents();
0400 }
0401
0402 CexmcRunSObject sObjectToWrite = {
0403 basePhysicsUsed, productionModelType, gdmlFileName, etaDecayTable,
0404 physicsManager->GetProductionModel()->GetAngularRanges(),
0405 physicsManager->GetProductionModel()->IsFermiMotionOn(),
0406 calorimeterRegCuts, eventCountPolicy,
0407 particleGun->GetParticleDefinition()->GetParticleName(),
0408 particleGun->GetOrigPosition(), particleGun->GetOrigDirection(),
0409 particleGun->GetOrigMomentumAmp(),
0410 primaryGeneratorAction->GetFwhmPosX(),
0411 primaryGeneratorAction->GetFwhmPosY(),
0412 primaryGeneratorAction->GetFwhmDirX(),
0413 primaryGeneratorAction->GetFwhmDirY(),
0414 primaryGeneratorAction->GetFwhmMomentumAmp(),
0415 edDigitizer->GetMonitorThreshold(),
0416 edDigitizer->GetVetoCounterLeftThreshold(),
0417 edDigitizer->GetVetoCounterRightThreshold(),
0418 edDigitizer->GetCalorimeterLeftThreshold(),
0419 edDigitizer->GetCalorimeterRightThreshold(),
0420 edDigitizer->GetCalorimeterTriggerAlgorithm(),
0421 edDigitizer->GetOuterCrystalsVetoAlgorithm(),
0422 edDigitizer->GetOuterCrystalsVetoFraction(),
0423 edDigitizer->IsFiniteCrystalResolutionApplied(),
0424 edDigitizer->GetCrystalResolutionData(),
0425 reconstructor->GetCalorimeterEntryPointDefinitionAlgorithm(),
0426 reconstructor->GetCalorimeterEntryPointDepthDefinitionAlgorithm(),
0427 reconstructor->GetCrystalSelectionAlgorithm(),
0428 reconstructor->IsInnerRefCrystalUsed(),
0429 reconstructor->GetCalorimeterEntryPointDepth(),
0430 reconstructor->IsTableMassUsed(), reconstructor->IsMassCutUsed(),
0431 reconstructor->GetMassCutOPCenter(),
0432 reconstructor->GetMassCutNOPCenter(),
0433 reconstructor->GetMassCutOPWidth(), reconstructor->GetMassCutNOPWidth(),
0434 reconstructor->GetMassCutEllipseAngle(),
0435 reconstructor->IsAbsorbedEnergyCutUsed(),
0436 reconstructor->GetAbsorbedEnergyCutCLCenter(),
0437 reconstructor->GetAbsorbedEnergyCutCRCenter(),
0438 reconstructor->GetAbsorbedEnergyCutCLWidth(),
0439 reconstructor->GetAbsorbedEnergyCutCRWidth(),
0440 reconstructor->GetAbsorbedEnergyCutEllipseAngle(),
0441 nmbOfHitsSampled, nmbOfHitsSampledFull, nmbOfHitsTriggeredRealRange,
0442 nmbOfHitsTriggeredRecRange, nmbOfOrphanHits, nmbOfFalseHitsTriggeredEDT,
0443 nmbOfFalseHitsTriggeredRec, nmbOfSavedEvents, nmbOfSavedFastEvents,
0444 numberOfEventsProcessed, numberOfEventsProcessedEffective,
0445 numberOfEventToBeProcessed, rProject, skipInteractionsWithoutEDTonWrite,
0446 cfFileName, evDataVerboseLevel, physicsManager->GetProposedMaxIL(),
0447 reconstructor->GetExpectedMomentumAmp(),
0448 reconstructor->GetEDCollectionAlgorithm(), 0 };
0449
0450 std::ofstream runDataFile( ( projectsDir + "/" + projectId + ".rdb" ).
0451 c_str() );
0452
0453 {
0454 boost::archive::binary_oarchive archive( runDataFile );
0455 archive << sObjectToWrite;
0456 }
0457 }
0458
0459 #endif
0460
0461
0462 void CexmcRunManager::DoCommonEventLoop( G4int nEvent, const G4String & cmd,
0463 G4int nSelect )
0464 {
0465 G4int iEvent( 0 );
0466 G4int iEventEffective( 0 );
0467
0468 for ( iEvent = 0; iEventEffective < nEvent; ++iEvent )
0469 {
0470 currentEvent = GenerateEvent( iEvent );
0471 eventManager->ProcessOneEvent( currentEvent );
0472 CexmcEventInfo * eventInfo( static_cast< CexmcEventInfo * >(
0473 currentEvent->GetUserInformation() ) );
0474 switch ( eventCountPolicy )
0475 {
0476 case CexmcCountAllEvents :
0477 ++iEventEffective;
0478 break;
0479 case CexmcCountEventsWithInteraction :
0480 if ( eventInfo->TpTriggerIsOk() )
0481 ++iEventEffective;
0482 break;
0483 case CexmcCountEventsWithTrigger :
0484 if ( eventInfo->EdTriggerIsOk() )
0485 ++iEventEffective;
0486 break;
0487 default :
0488 ++iEventEffective;
0489 break;
0490 }
0491 AnalyzeEvent( currentEvent );
0492 UpdateScoring();
0493 if ( iEvent < nSelect )
0494 G4UImanager::GetUIpointer()->ApplyCommand( cmd );
0495 StackPreviousEvent( currentEvent );
0496 currentEvent = 0;
0497 if ( runAborted )
0498 break;
0499 }
0500
0501 numberOfEventsProcessed = iEvent;
0502 numberOfEventsProcessedEffective = iEventEffective;
0503 }
0504
0505
0506 #ifdef CEXMC_USE_PERSISTENCY
0507
0508 void CexmcRunManager::DoReadEventLoop( G4int nEvent )
0509 {
0510 G4int iEvent( 0 );
0511 G4int iEventEffective( 0 );
0512 G4int nEventCount( 0 );
0513
0514 if ( ! ProjectIsRead() )
0515 return;
0516
0517 if ( ! physicsManager )
0518 throw CexmcException( CexmcWeirdException );
0519
0520 CexmcProductionModel * productionModel(
0521 physicsManager->GetProductionModel() );
0522 if ( ! productionModel )
0523 throw CexmcException( CexmcWeirdException );
0524
0525 CexmcSetup * setup( static_cast< CexmcSetup * >( userDetector ) );
0526 if ( ! setup )
0527 throw CexmcException( CexmcWeirdException );
0528
0529 CexmcEventSObject evSObject;
0530 CexmcEventFastSObject evFastSObject;
0531
0532
0533 std::ifstream eventsDataFile(
0534 ( projectsDir + "/" + rProject + ".edb" ).c_str() );
0535 if ( ! eventsDataFile )
0536 throw CexmcException( CexmcReadProjectIncomplete );
0537
0538 boost::archive::binary_iarchive evArchive( eventsDataFile );
0539
0540 std::ifstream eventsFastDataFile(
0541 ( projectsDir + "/" + rProject + ".fdb" ).c_str() );
0542 if ( ! eventsFastDataFile )
0543 throw CexmcException( CexmcReadProjectIncomplete );
0544
0545 boost::archive::binary_iarchive evFastArchive( eventsFastDataFile );
0546
0547 G4Event event;
0548 currentEvent = &event;
0549 G4SDManager * sdManager( G4SDManager::GetSDMpointer() );
0550 event.SetHCofThisEvent( sdManager->PrepareNewEvent() );
0551 G4HCofThisEvent * hcOfThisEvent( event.GetHCofThisEvent() );
0552
0553 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
0554
0555 G4int hcId( digiManager->GetHitsCollectionID(
0556 CexmcDetectorRoleName[ CexmcMonitorDetectorRole ] +
0557 "/" + CexmcDetectorTypeName[ CexmcEDDetector ] ) );
0558 CexmcEnergyDepositCollection * monitorED(
0559 new CexmcEnergyDepositCollection );
0560 hcOfThisEvent->AddHitsCollection( hcId, monitorED );
0561 hcId = digiManager->GetHitsCollectionID(
0562 CexmcDetectorRoleName[ CexmcVetoCounterDetectorRole ] +
0563 "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
0564 CexmcEnergyDepositCollection * vetoCounterED(
0565 new CexmcEnergyDepositCollection );
0566 hcOfThisEvent->AddHitsCollection( hcId, vetoCounterED );
0567 hcId = digiManager->GetHitsCollectionID(
0568 CexmcDetectorRoleName[ CexmcCalorimeterDetectorRole ] +
0569 "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
0570 CexmcEnergyDepositCollection * calorimeterED(
0571 new CexmcEnergyDepositCollection );
0572 hcOfThisEvent->AddHitsCollection( hcId, calorimeterED );
0573 hcId = digiManager->GetHitsCollectionID(
0574 CexmcDetectorRoleName[ CexmcMonitorDetectorRole ] +
0575 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
0576 CexmcTrackPointsCollection * monitorTP( new CexmcTrackPointsCollection );
0577 hcOfThisEvent->AddHitsCollection( hcId, monitorTP );
0578 hcId = digiManager->GetHitsCollectionID(
0579 CexmcDetectorRoleName[ CexmcVetoCounterDetectorRole ] +
0580 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
0581 CexmcTrackPointsCollection * vetoCounterTP(
0582 new CexmcTrackPointsCollection );
0583 hcOfThisEvent->AddHitsCollection( hcId, vetoCounterTP );
0584 hcId = digiManager->GetHitsCollectionID(
0585 CexmcDetectorRoleName[ CexmcCalorimeterDetectorRole ] +
0586 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
0587 CexmcTrackPointsCollection * calorimeterTP(
0588 new CexmcTrackPointsCollection );
0589 hcOfThisEvent->AddHitsCollection( hcId, calorimeterTP );
0590 hcId = digiManager->GetHitsCollectionID(
0591 CexmcDetectorRoleName[ CexmcTargetDetectorRole ] +
0592 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
0593 CexmcTrackPointsCollection * targetTP( new CexmcTrackPointsCollection );
0594 hcOfThisEvent->AddHitsCollection( hcId, targetTP );
0595
0596 #ifdef CEXMC_USE_CUSTOM_FILTER
0597 if ( customFilter )
0598 customFilter->SetAddressedData( &evFastSObject, &evSObject );
0599 #endif
0600
0601 G4int nmbOfSavedEvents( rEvDataVerboseLevel == CexmcWriteNoEventData ? 0 :
0602 sObject.nmbOfSavedFastEvents );
0603 G4bool eventDataWrittenOnEveryTPT( rEvDataVerboseLevel ==
0604 CexmcWriteEventDataOnEveryTPT );
0605
0606 for ( int i( 0 ); i < nmbOfSavedEvents; ++i )
0607 {
0608 evFastArchive >> evFastSObject;
0609
0610 if ( nEventCount < curEventRead )
0611 {
0612 if ( eventDataWrittenOnEveryTPT ||
0613 evFastSObject.edDigitizerHasTriggered )
0614 {
0615 evArchive >> evSObject;
0616 if ( evFastSObject.edDigitizerHasTriggered )
0617 ++nEventCount;
0618 }
0619 continue;
0620 }
0621
0622 ++iEvent;
0623
0624 productionModel->SetTriggeredAngularRanges(
0625 evFastSObject.opCosThetaSCM );
0626 const CexmcAngularRangeList & triggeredAngularRanges(
0627 productionModel->GetTriggeredAngularRanges() );
0628
0629 if ( ! eventDataWrittenOnEveryTPT &&
0630 ! evFastSObject.edDigitizerHasTriggered )
0631 {
0632 #ifdef CEXMC_USE_CUSTOM_FILTER
0633
0634
0635
0636
0637
0638 if ( customFilter && ! customFilter->EvalTPT() )
0639 continue;
0640 #endif
0641 SaveCurrentTPTEvent( evFastSObject, triggeredAngularRanges,
0642 ProjectIsSaved() && ! skipInteractionsWithoutEDTonWrite );
0643 continue;
0644 }
0645
0646 evArchive >> evSObject;
0647
0648 G4bool skipEDTOnThisEvent( false );
0649
0650 #ifdef CEXMC_USE_CUSTOM_FILTER
0651 if ( customFilter && ! customFilter->EvalTPT() )
0652 continue;
0653 if ( customFilter && ! customFilter->EvalEDT() )
0654 {
0655 if ( ! eventDataWrittenOnEveryTPT )
0656 {
0657 SaveCurrentTPTEvent( evFastSObject, triggeredAngularRanges,
0658 ProjectIsSaved() );
0659 continue;
0660 }
0661 skipEDTOnThisEvent = true;
0662 }
0663 #endif
0664
0665 event.SetEventID( evSObject.eventId );
0666
0667
0668
0669
0670
0671
0672 monitorED->GetMap()->operator[]( 0 ) = &evSObject.monitorED;
0673 vetoCounterED->GetMap()->operator[]( 0 ) = &evSObject.vetoCounterEDLeft;
0674 vetoCounterED->GetMap()->operator[]( 1 <<
0675 CexmcEnergyDepositInLeftRightSet::GetLeftRightBitsOffset() ) =
0676 &evSObject.vetoCounterEDRight;
0677 G4int row( 0 );
0678 G4int column( 0 );
0679 for ( CexmcEnergyDepositCalorimeterCollection::iterator
0680 k( evSObject.calorimeterEDLeftCollection.begin() );
0681 k != evSObject.calorimeterEDLeftCollection.end(); ++k )
0682 {
0683 G4int index( row <<
0684 CexmcEnergyDepositInCalorimeter::GetCopyDepth1BitsOffset() );
0685 column = 0;
0686 for ( CexmcEnergyDepositCrystalRowCollection::iterator
0687 l( k->begin() ); l != k->end(); ++l )
0688 {
0689 calorimeterED->GetMap()->operator[]( index | column ) = &*l;
0690 ++column;
0691 }
0692 ++row;
0693 }
0694 row = 0;
0695 for ( CexmcEnergyDepositCalorimeterCollection::iterator
0696 k( evSObject.calorimeterEDRightCollection.begin() );
0697 k != evSObject.calorimeterEDRightCollection.end(); ++k )
0698 {
0699 G4int index(
0700 1 << CexmcEnergyDepositInLeftRightSet::GetLeftRightBitsOffset()
0701 | row <<
0702 CexmcEnergyDepositInCalorimeter::GetCopyDepth1BitsOffset() );
0703 column = 0;
0704 for ( CexmcEnergyDepositCrystalRowCollection::iterator
0705 l( k->begin() ); l != k->end(); ++l )
0706 {
0707 calorimeterED->GetMap()->operator[]( index | column ) = &*l;
0708 ++column;
0709 }
0710 ++row;
0711 }
0712
0713 CexmcTrackPointInfo monitorTPInfo( evSObject.monitorTP );
0714 CexmcTrackPointInfo targetTPBeamParticleInfo(
0715 evSObject.targetTPBeamParticle );
0716 CexmcTrackPointInfo targetTPOutputParticleInfo(
0717 evSObject.targetTPOutputParticle );
0718 CexmcTrackPointInfo targetTPNucleusParticleInfo(
0719 evSObject.targetTPNucleusParticle );
0720 CexmcTrackPointInfo targetTPOutputParticleDecayProductParticle1Info(
0721 evSObject.targetTPOutputParticleDecayProductParticle1 );
0722 CexmcTrackPointInfo targetTPOutputParticleDecayProductParticle2Info(
0723 evSObject.targetTPOutputParticleDecayProductParticle2 );
0724 CexmcTrackPointInfo vetoCounterTPLeftInfo(
0725 evSObject.vetoCounterTPLeft );
0726 CexmcTrackPointInfo vetoCounterTPRightInfo(
0727 evSObject.vetoCounterTPRight );
0728 CexmcTrackPointInfo calorimeterTPLeftInfo(
0729 evSObject.calorimeterTPLeft );
0730 CexmcTrackPointInfo calorimeterTPRightInfo(
0731 evSObject.calorimeterTPRight );
0732
0733 if ( monitorTPInfo.IsValid() )
0734 monitorTP->GetMap()->operator[]( monitorTPInfo.trackId ) =
0735 &monitorTPInfo;
0736 if ( targetTPBeamParticleInfo.IsValid() )
0737 targetTP->GetMap()->operator[](
0738 targetTPBeamParticleInfo.trackId ) =
0739 &targetTPBeamParticleInfo;
0740 if ( targetTPOutputParticleInfo.IsValid() )
0741 targetTP->GetMap()->operator[](
0742 targetTPOutputParticleInfo.trackId ) =
0743 &targetTPOutputParticleInfo;
0744 if ( targetTPNucleusParticleInfo.IsValid() )
0745 targetTP->GetMap()->operator[](
0746 targetTPNucleusParticleInfo.trackId ) =
0747 &targetTPNucleusParticleInfo;
0748 if ( targetTPOutputParticleDecayProductParticle1Info.IsValid() )
0749 targetTP->GetMap()->operator[](
0750 targetTPOutputParticleDecayProductParticle1Info.trackId ) =
0751 &targetTPOutputParticleDecayProductParticle1Info;
0752 if ( targetTPOutputParticleDecayProductParticle2Info.IsValid() )
0753 targetTP->GetMap()->operator[](
0754 targetTPOutputParticleDecayProductParticle2Info.trackId ) =
0755 &targetTPOutputParticleDecayProductParticle2Info;
0756 if ( vetoCounterTPLeftInfo.IsValid() )
0757 vetoCounterTP->GetMap()->operator[](
0758 vetoCounterTPLeftInfo.trackId ) = &vetoCounterTPLeftInfo;
0759 if ( vetoCounterTPRightInfo.IsValid() )
0760 vetoCounterTP->GetMap()->operator[](
0761 1 << CexmcTrackPointsInLeftRightSet::GetLeftRightBitsOffset() |
0762 vetoCounterTPRightInfo.trackId ) = &vetoCounterTPRightInfo;
0763
0764 G4ThreeVector pos;
0765 if ( calorimeterTPLeftInfo.IsValid() )
0766 {
0767 pos = calorimeterTPLeftInfo.positionLocal;
0768 setup->ConvertToCrystalGeometry(
0769 calorimeterTPLeftInfo.positionLocal, row, column, pos );
0770 calorimeterTPLeftInfo.positionLocal = pos;
0771 calorimeterTP->GetMap()->operator[](
0772 row << CexmcTrackPointsInCalorimeter::
0773 GetCopyDepth1BitsOffset() |
0774 column << CexmcTrackPointsInCalorimeter::
0775 GetCopyDepth0BitsOffset() |
0776 calorimeterTPLeftInfo.trackId ) = &calorimeterTPLeftInfo;
0777 }
0778 if ( calorimeterTPRightInfo.IsValid() )
0779 {
0780 pos = calorimeterTPRightInfo.positionLocal;
0781 setup->ConvertToCrystalGeometry(
0782 calorimeterTPRightInfo.positionLocal, row, column, pos );
0783 calorimeterTPRightInfo.positionLocal = pos;
0784 calorimeterTP->GetMap()->operator[](
0785 1 << CexmcTrackPointsInLeftRightSet::GetLeftRightBitsOffset() |
0786 row << CexmcTrackPointsInCalorimeter::
0787 GetCopyDepth1BitsOffset() |
0788 column << CexmcTrackPointsInCalorimeter::
0789 GetCopyDepth0BitsOffset() |
0790 calorimeterTPRightInfo.trackId ) = &calorimeterTPRightInfo;
0791 }
0792
0793 productionModel->SetProductionModelData(
0794 evSObject.productionModelData );
0795
0796 const CexmcEventAction * eventAction(
0797 static_cast< const CexmcEventAction * >( userEventAction ) );
0798 if ( ! eventAction )
0799 {
0800
0801
0802
0803 monitorED->GetMap()->clear();
0804 vetoCounterED->GetMap()->clear();
0805 calorimeterED->GetMap()->clear();
0806 monitorTP->GetMap()->clear();
0807 targetTP->GetMap()->clear();
0808 vetoCounterTP->GetMap()->clear();
0809 calorimeterTP->GetMap()->clear();
0810 throw CexmcException( CexmcEventActionIsNotInitialized );
0811 }
0812
0813 if ( skipEDTOnThisEvent )
0814 event.SetUserInformation( new CexmcEventInfo( false, false,
0815 false ) );
0816
0817 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
0818 eventAction ) );
0819 theEventAction->EndOfEventAction( &event );
0820
0821 CexmcEventInfo * eventInfo( static_cast< CexmcEventInfo * >(
0822 event.GetUserInformation() ) );
0823
0824 if ( eventInfo->EdTriggerIsOk() )
0825 ++iEventEffective;
0826
0827 delete eventInfo;
0828 event.SetUserInformation( NULL );
0829
0830 monitorED->GetMap()->clear();
0831 vetoCounterED->GetMap()->clear();
0832 calorimeterED->GetMap()->clear();
0833 monitorTP->GetMap()->clear();
0834 targetTP->GetMap()->clear();
0835 vetoCounterTP->GetMap()->clear();
0836 calorimeterTP->GetMap()->clear();
0837
0838
0839
0840 if ( nEvent > 0 && iEventEffective == nEvent )
0841 break;
0842 }
0843
0844 curEventRead = nEventCount + iEventEffective;
0845
0846 numberOfEventsProcessed = iEvent;
0847 numberOfEventsProcessedEffective = iEventEffective;
0848
0849 #ifdef CEXMC_USE_CUSTOM_FILTER
0850 if ( customFilter )
0851 customFilter->SetAddressedData( NULL, NULL );
0852 #endif
0853 }
0854
0855
0856 void CexmcRunManager::SaveCurrentTPTEvent(
0857 const CexmcEventFastSObject & evFastSObject,
0858 const CexmcAngularRangeList & angularRanges,
0859 G4bool writeToDatabase )
0860 {
0861 CexmcRun * run( static_cast< CexmcRun * >( currentRun ) );
0862
0863 if ( ! run )
0864 return;
0865
0866 for ( CexmcAngularRangeList::const_iterator k( angularRanges.begin() );
0867 k != angularRanges.end(); ++k )
0868 {
0869 run->IncrementNmbOfHitsSampledFull( k->index );
0870 if ( evFastSObject.edDigitizerMonitorHasTriggered )
0871 run->IncrementNmbOfHitsSampled( k->index );
0872 }
0873
0874 if ( writeToDatabase )
0875 {
0876 fastEventsArchive->operator<<( evFastSObject );
0877 run->IncrementNmbOfSavedFastEvents();
0878 }
0879 }
0880
0881 #endif
0882
0883
0884
0885 void CexmcRunManager::DoEventLoop( G4int nEvent, const char * macroFile,
0886 G4int nSelect )
0887 {
0888 if ( verboseLevel > 0 )
0889 timer->Start();
0890
0891 G4String cmd;
0892 if ( macroFile != 0 )
0893 {
0894 if ( nSelect < 0 )
0895 nSelect = nEvent;
0896 cmd = "/control/execute ";
0897 cmd += macroFile;
0898 }
0899 else
0900 {
0901 nSelect = -1;
0902 }
0903
0904 numberOfEventsProcessed = 0;
0905 numberOfEventsProcessedEffective = 0;
0906
0907 #ifdef CEXMC_USE_PERSISTENCY
0908 eventsArchive = NULL;
0909 fastEventsArchive = NULL;
0910 if ( ProjectIsRead() )
0911 {
0912 if ( ProjectIsSaved() )
0913 {
0914 std::ofstream eventsDataFile(
0915 ( projectsDir + "/" + projectId + ".edb" ).c_str() );
0916 boost::archive::binary_oarchive eventsArchive_( eventsDataFile );
0917 std::ofstream fastEventsDataFile(
0918 ( projectsDir + "/" + projectId + ".fdb" ).c_str() );
0919 boost::archive::binary_oarchive fastEventsArchive_(
0920 fastEventsDataFile );
0921 eventsArchive = &eventsArchive_;
0922 fastEventsArchive = &fastEventsArchive_;
0923 DoReadEventLoop( nEvent );
0924 }
0925 else
0926 {
0927 DoReadEventLoop( nEvent );
0928 }
0929 }
0930 else
0931 {
0932 if ( ProjectIsSaved() )
0933 {
0934 std::ofstream eventsDataFile(
0935 ( projectsDir + "/" + projectId + ".edb" ).c_str() );
0936 boost::archive::binary_oarchive eventsArchive_( eventsDataFile );
0937 std::ofstream fastEventsDataFile(
0938 ( projectsDir + "/" + projectId + ".fdb" ).c_str() );
0939 boost::archive::binary_oarchive fastEventsArchive_(
0940 fastEventsDataFile );
0941 eventsArchive = &eventsArchive_;
0942 fastEventsArchive = &fastEventsArchive_;
0943 DoCommonEventLoop( nEvent, cmd, nSelect );
0944 }
0945 else
0946 {
0947 DoCommonEventLoop( nEvent, cmd, nSelect );
0948 }
0949 }
0950 eventsArchive = NULL;
0951 fastEventsArchive = NULL;
0952 #else
0953 DoCommonEventLoop( nEvent, cmd, nSelect );
0954 #endif
0955
0956 if ( verboseLevel > 0 )
0957 {
0958 timer->Stop();
0959 G4cout << "Run terminated." << G4endl;
0960 G4cout << "Run Summary" << G4endl;
0961 if ( runAborted )
0962 {
0963 G4cout << " Run Aborted after " << numberOfEventsProcessed <<
0964 " events processed." << G4endl;
0965 }
0966 else
0967 {
0968 G4cout << " Number of events processed : " <<
0969 numberOfEventsProcessed << ", effectively: " <<
0970 numberOfEventsProcessedEffective << G4endl;
0971 }
0972 G4cout << " " << *timer << G4endl;
0973 }
0974 }
0975
0976
0977 #ifdef CEXMC_USE_PERSISTENCY
0978
0979 void CexmcRunManager::PrintReadRunData( void ) const
0980 {
0981 if ( ! ProjectIsRead() )
0982 return;
0983
0984 G4bool refCrystalInfoPrinted( false );
0985
0986 G4cout << CEXMC_LINE_START << "Run data read from project '" << rProject <<
0987 "'" << G4endl;
0988 G4cout << " (archive class version " <<
0989 sObject.actualVersion << ")" << G4endl;
0990 if ( ! sObject.rProject.empty() )
0991 {
0992 G4cout << " -- Based on project '" << sObject.rProject << "'" <<
0993 G4endl;
0994 if ( ! sObject.cfFileName.empty() )
0995 G4cout << " -- Custom filter script '" << sObject.cfFileName <<
0996 "' was used" << G4endl;
0997 }
0998 G4cout << " -- Event data verbose level (0 - not saved, 1 - triggers, "
0999 "2 - interactions): " << sObject.evDataVerboseLevel << G4endl;
1000 if ( ! sObject.rProject.empty() )
1001 {
1002 if ( sObject.evDataVerboseLevel == CexmcWriteEventDataOnEveryEDT )
1003 {
1004 G4cout << " -- (fdb file contains " <<
1005 ( sObject.interactionsWithoutEDTWereSkipped ?
1006 "only interactions when an event was triggered" :
1007 "all interactions" ) << ")" << std::endl;
1008 }
1009 }
1010 G4cout << " -- Base physics used"
1011 "(1 - QGSP_BERT, 2 - QGSP_BIC_EMY, 3 - FTFP_BERT): " <<
1012 sObject.basePhysicsUsed << G4endl;
1013 G4cout << " -- Production model (1 - pi0, 2 - eta): " <<
1014 sObject.productionModelType << G4endl;
1015 G4cout << " -- Geometry definition file: " << sObject.gdmlFileName <<
1016 G4endl;
1017 G4cout << " -- Angular ranges: " << sObject.angularRanges << G4endl;
1018 G4cout << " -- Eta decay modes: " << G4endl;
1019 G4Eta::Definition()->GetDecayTable()->DumpInfo();
1020 G4cout << " -- Fermi motion status (0 - disabled, 1 - enabled): " <<
1021 sObject.fermiMotionIsOn << G4endl;
1022 if ( sObject.calorimeterRegCuts.size() < 4 )
1023 throw CexmcException( CexmcWeirdException );
1024 G4cout << " -- Production cuts in calorimeter (gamma, e-, e+, p): " <<
1025 G4BestUnit( sObject.calorimeterRegCuts[ 0 ], "Length" ) << ", " <<
1026 G4BestUnit( sObject.calorimeterRegCuts[ 1 ], "Length" ) << ", " <<
1027 G4BestUnit( sObject.calorimeterRegCuts[ 2 ], "Length" ) << ", " <<
1028 G4BestUnit( sObject.calorimeterRegCuts[ 3 ], "Length" ) << G4endl;
1029 G4cout << " -- Proposed max interaction length in the target: " <<
1030 G4BestUnit( sObject.proposedMaxIL, "Length" ) << G4endl;
1031 G4cout << " -- Event count policy (0 - all, 1 - interaction, 2 - trigger)"
1032 ": " << sObject.eventCountPolicy << G4endl;
1033 G4cout << " -- Number of events (processed / effective / ordered): " <<
1034 sObject.numberOfEventsProcessed << " / " <<
1035 sObject.numberOfEventsProcessedEffective << " / " <<
1036 sObject.numberOfEventsToBeProcessed << G4endl;
1037 G4cout << " -- Incident beam particle: " << sObject.beamParticle << G4endl;
1038 G4cout << " position: " <<
1039 G4BestUnit( sObject.beamPos, "Length" ) << G4endl;
1040 G4cout << " direction: " <<
1041 G4ThreeVector( sObject.beamDir ) << G4endl;
1042 G4cout << " momentum: " <<
1043 G4BestUnit( sObject.beamMomentumAmp, "Energy" ) << G4endl;
1044 G4cout << " momentum fwhm: " << sObject.beamFwhmMomentumAmp <<
1045 G4endl;
1046 G4cout << " pos fwhm (x): " <<
1047 G4BestUnit( sObject.beamFwhmPosX, "Length" ) << G4endl;
1048 G4cout << " pos fwhm (y): " <<
1049 G4BestUnit( sObject.beamFwhmPosY, "Length" ) << G4endl;
1050 G4cout << " dir fwhm (x): " << sObject.beamFwhmDirX / deg <<
1051 " deg" << G4endl;
1052 G4cout << " dir fwhm (y): " << sObject.beamFwhmDirY / deg <<
1053 " deg" << G4endl;
1054 G4cout << " -- Monitor ED threshold: " <<
1055 G4BestUnit( sObject.monitorEDThreshold, "Energy" ) << G4endl;
1056 G4cout << " -- Veto counter (l/r) ED threshold: " <<
1057 G4BestUnit( sObject.vetoCounterEDLeftThreshold, "Energy" ) <<
1058 " / " <<
1059 G4BestUnit( sObject.vetoCounterEDRightThreshold, "Energy" ) <<
1060 G4endl;
1061 G4cout << " -- Calorimeter (l/r) ED threshold: " <<
1062 G4BestUnit( sObject.calorimeterEDLeftThreshold, "Energy" ) <<
1063 " / " <<
1064 G4BestUnit( sObject.calorimeterEDRightThreshold, "Energy" ) <<
1065 G4endl;
1066 G4cout << " -- Calorimeter trigger algorithm (0 - all, 1 - inner): " <<
1067 sObject.calorimeterTriggerAlgorithm << G4endl;
1068 G4cout << " -- Outer crystals veto algorithm "
1069 "(0 - none, 1 - max, 2 - fraction): " <<
1070 sObject.outerCrystalsVetoAlgorithm << G4endl;
1071 if ( sObject.outerCrystalsVetoAlgorithm ==
1072 CexmcFractionOfEDInOuterCrystalsVeto )
1073 {
1074 G4cout << " -- Outer crystals veto fraction: " <<
1075 sObject.outerCrystalsVetoFraction << G4endl;
1076 }
1077 G4cout << " -- Finite crystal resolution applied (0 - no, 1 - yes): " <<
1078 sObject.applyFiniteCrystalResolution << G4endl;
1079 if ( sObject.applyFiniteCrystalResolution )
1080 {
1081 G4cout << " -- Crystal resolution data: " <<
1082 sObject.crystalResolutionData;
1083 }
1084 G4cout << " -- Reconstructor settings: " << G4endl;
1085 if ( sObject.expectedMomentumAmp > 0 )
1086 {
1087 G4cout << " -- expected momentum in the target: " <<
1088 G4BestUnit( sObject.expectedMomentumAmp, "Energy" ) << G4endl;
1089 }
1090 G4cout << " -- ed collection algorithm (0 - all, 1 - adjacent): " <<
1091 sObject.edCollectionAlgorithm << G4endl;
1092 if ( sObject.edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
1093 {
1094 G4cout <<
1095 " -- inner crystal used as reference (0 - no, 1 - yes): " <<
1096 sObject.useInnerRefCrystal << G4endl;
1097 refCrystalInfoPrinted = true;
1098 }
1099 G4cout << " -- entry point definition algorithm " << G4endl;
1100 G4cout << " (0 - center of calorimeter, 1 - center of crystal with "
1101 "max ED," << G4endl;
1102 G4cout << " 2 - linear, 3 - square): " <<
1103 sObject.epDefinitionAlgorithm << G4endl;
1104 G4cout << " -- entry point depth definition algorithm "
1105 "(0 - plain, 1 - sphere): " <<
1106 sObject.epDepthDefinitionAlgorithm << G4endl;
1107 G4cout << " -- entry point depth: " <<
1108 G4BestUnit( sObject.epDepth, "Length" ) << G4endl;
1109 if ( sObject.epDefinitionAlgorithm == CexmcEntryPointByLinearEDWeights ||
1110 sObject.epDefinitionAlgorithm == CexmcEntryPointBySqrtEDWeights )
1111 {
1112 G4cout <<
1113 " -- crystal selection algorithm (0 - all, 1 - adjacent): " <<
1114 sObject.csAlgorithm << G4endl;
1115 }
1116 if ( ! refCrystalInfoPrinted &&
1117 ( sObject.epDefinitionAlgorithm ==
1118 CexmcEntryPointInTheCenterOfCrystalWithMaxED ||
1119 ( ( sObject.epDefinitionAlgorithm == CexmcEntryPointBySqrtEDWeights ||
1120 sObject.epDefinitionAlgorithm ==
1121 CexmcEntryPointByLinearEDWeights ) &&
1122 sObject.csAlgorithm == CexmcSelectAdjacentCrystals ) ) )
1123 {
1124 G4cout <<
1125 " -- inner crystal used as reference (0 - no, 1 - yes): " <<
1126 sObject.useInnerRefCrystal << G4endl;
1127 }
1128 G4cout << " -- table mass of output particle used "
1129 "(0 - no, 1 - yes): " << sObject.useTableMass << G4endl;
1130 G4cout << " -- mass cut is enabled (0 - no, 1 - yes): " <<
1131 sObject.useMassCut << G4endl;
1132 if ( sObject.useMassCut )
1133 {
1134 G4cout << " -- mass cut output particle center: " <<
1135 G4BestUnit( sObject.mCutOPCenter, "Energy" ) << G4endl;
1136 G4cout << " -- mass cut nucleus output particle center: " <<
1137 G4BestUnit( sObject.mCutNOPCenter, "Energy" ) << G4endl;
1138 G4cout << " -- mass cut output particle width of the ellipse: " <<
1139 G4BestUnit( sObject.mCutOPWidth, "Energy" ) << G4endl;
1140 G4cout << " -- mass cut nucleus output particle width of the "
1141 "ellipse: "
1142 << G4BestUnit( sObject.mCutNOPWidth, "Energy" ) << G4endl;
1143 G4cout << " -- mass cut angle of the ellipse: " <<
1144 sObject.mCutAngle / deg << " deg" << G4endl;
1145 }
1146 G4cout << " -- absorbed energy cut is enabled (0 - no, 1 - yes): " <<
1147 sObject.useAbsorbedEnergyCut << G4endl;
1148 if ( sObject.useAbsorbedEnergyCut )
1149 {
1150 G4cout << " -- absorbed energy cut left calorimeter center: " <<
1151 G4BestUnit( sObject.aeCutCLCenter, "Energy" ) << G4endl;
1152 G4cout << " -- absorbed energy cut right calorimeter center: " <<
1153 G4BestUnit( sObject.aeCutCRCenter, "Energy" ) << G4endl;
1154 G4cout << " -- absorbed energy cut left calorimeter width of the "
1155 "ellipse: " <<
1156 G4BestUnit( sObject.aeCutCLWidth, "Energy" ) << G4endl;
1157 G4cout << " -- absorbed energy cut right calorimeter width of the "
1158 "ellipse: "
1159 << G4BestUnit( sObject.aeCutCRWidth, "Energy" ) << G4endl;
1160 G4cout << " -- absorbed energy cut angle of the ellipse: " <<
1161 sObject.aeCutAngle / deg << " deg" << G4endl;
1162 }
1163 G4cout << G4endl;
1164 CexmcRunAction::PrintResults( sObject.nmbOfHitsSampled,
1165 sObject.nmbOfHitsSampledFull,
1166 sObject.nmbOfHitsTriggeredRealRange,
1167 sObject.nmbOfHitsTriggeredRecRange,
1168 sObject.nmbOfOrphanHits,
1169 sObject.angularRanges,
1170 sObject.nmbOfFalseHitsTriggeredEDT,
1171 sObject.nmbOfFalseHitsTriggeredRec );
1172 G4cout << G4endl;
1173 }
1174
1175
1176 void CexmcRunManager::ReadAndPrintEventsData( void ) const
1177 {
1178 if ( ! ProjectIsRead() )
1179 return;
1180
1181 CexmcEventSObject evSObject;
1182
1183
1184 std::ifstream eventsDataFile(
1185 ( projectsDir + "/" + rProject + ".edb" ).c_str() );
1186 if ( ! eventsDataFile )
1187 throw CexmcException( CexmcReadProjectIncomplete );
1188
1189 boost::archive::binary_iarchive evArchive( eventsDataFile );
1190
1191 for ( int i( 0 ); i < sObject.nmbOfSavedEvents; ++i )
1192 {
1193 evArchive >> evSObject;
1194
1195 if ( ! evSObject.edDigitizerMonitorHasTriggered )
1196 continue;
1197
1198 CexmcEnergyDepositStore edStore( evSObject.monitorED,
1199 evSObject.vetoCounterEDLeft, evSObject.vetoCounterEDRight,
1200 evSObject.calorimeterEDLeft, evSObject.calorimeterEDRight,
1201 0, 0, 0, 0, evSObject.calorimeterEDLeftCollection,
1202 evSObject.calorimeterEDRightCollection );
1203
1204 CexmcTrackPointsStore tpStore( evSObject.monitorTP,
1205 evSObject.targetTPBeamParticle, evSObject.targetTPOutputParticle,
1206 evSObject.targetTPNucleusParticle,
1207 evSObject.targetTPOutputParticleDecayProductParticle1,
1208 evSObject.targetTPOutputParticleDecayProductParticle2,
1209 evSObject.vetoCounterTPLeft, evSObject.vetoCounterTPRight,
1210 evSObject.calorimeterTPLeft, evSObject.calorimeterTPRight );
1211
1212 const CexmcProductionModelData & pmData(
1213 evSObject.productionModelData );
1214
1215 G4cout << "Event " << evSObject.eventId << G4endl;
1216 CexmcEventAction::PrintTrackPoints( &tpStore );
1217 G4cout << " --- Production model data: " << pmData;
1218 CexmcEventAction::PrintEnergyDeposit( &edStore );
1219 }
1220 }
1221
1222
1223 void CexmcRunManager::PrintReadData(
1224 const CexmcOutputDataTypeSet & outputData ) const
1225 {
1226 if ( ! ProjectIsRead() )
1227 return;
1228
1229 G4bool addSpace( false );
1230
1231 CexmcOutputDataTypeSet::const_iterator found(
1232 outputData.find( CexmcOutputGeometry ) );
1233 if ( found != outputData.end() )
1234 {
1235 G4String cmd( G4String( "cat " ) + projectsDir + "/" + rProject +
1236 gdmlFileExtension );
1237 if ( system( cmd ) != 0 )
1238 throw CexmcException( CexmcReadProjectIncomplete );
1239
1240 if ( zipGdmlFile )
1241 {
1242 cmd = G4String( "bzip2 " ) + projectsDir + "/" + rProject +
1243 gdmlFileExtension;
1244 if ( system( cmd ) != 0 )
1245 throw CexmcException( CexmcFileCompressException );
1246 }
1247
1248 addSpace = true;
1249 }
1250
1251 found = outputData.find( CexmcOutputEvents );
1252 if ( found != outputData.end() )
1253 {
1254 if ( addSpace )
1255 G4cout << G4endl << G4endl;
1256
1257 ReadAndPrintEventsData();
1258
1259 addSpace = true;
1260 }
1261
1262 found = outputData.find( CexmcOutputRun );
1263 if ( found != outputData.end() )
1264 {
1265 if ( addSpace )
1266 G4cout << G4endl << G4endl;
1267
1268 G4DecayTable * etaDecayTable( G4Eta::Definition()->GetDecayTable() );
1269 for ( CexmcDecayBranchesStore::const_iterator
1270 k( sObject.etaDecayTable.GetDecayBranches().begin() );
1271 k != sObject.etaDecayTable.GetDecayBranches().end(); ++k )
1272 {
1273 etaDecayTable->GetDecayChannel( k->first )->SetBR( k->second );
1274 }
1275
1276 PrintReadRunData();
1277 }
1278 }
1279
1280
1281 #ifdef CEXMC_USE_CUSTOM_FILTER
1282
1283 void CexmcRunManager::SetCustomFilter( const G4String & cfFileName_ )
1284 {
1285 if ( customFilter )
1286 {
1287 delete customFilter;
1288 customFilter = NULL;
1289 }
1290
1291 if ( cfFileName_.empty() )
1292 return;
1293
1294
1295 if ( ! ProjectIsRead() )
1296 throw CexmcException( CexmcCmdIsNotAllowed );
1297
1298 cfFileName = cfFileName_;
1299
1300 customFilter = new CexmcCustomFilterEval( cfFileName );
1301 }
1302
1303 #endif
1304
1305 #endif
1306
1307
1308 void CexmcRunManager::RegisterScenePrimitives( void )
1309 {
1310 G4VisManager * visManager( static_cast< G4VisManager * >(
1311 G4VVisManager::GetConcreteInstance() ) );
1312 if ( ! visManager )
1313 return;
1314
1315 G4Scene * curScene( visManager->GetCurrentScene() );
1316 if ( ! curScene )
1317 return;
1318
1319
1320 #if G4VERSION_NUMBER < 960
1321 typedef std::vector< G4VModel * > MList;
1322 #else
1323 typedef std::vector< G4Scene::Model > MList;
1324 #endif
1325 const MList & mList( curScene->GetRunDurationModelList() );
1326
1327 for ( MList::const_iterator k( mList.begin() ); k != mList.end(); ++k )
1328 {
1329 #if G4VERSION_NUMBER < 960
1330 const G4String & modelDesc( ( *k )->GetGlobalDescription() );
1331 #else
1332 const G4String & modelDesc( k->fpModel->GetGlobalDescription() );
1333 #endif
1334 if ( modelDesc == CexmcScenePrimitivesDescription )
1335 return;
1336 }
1337
1338 CexmcSetup * setup( static_cast< CexmcSetup * >( userDetector ) );
1339 if ( ! setup )
1340 throw CexmcException( CexmcWeirdException );
1341
1342
1343
1344 curScene->AddRunDurationModel( new CexmcScenePrimitives( setup ) );
1345 }
1346
1347
1348 void CexmcRunManager::BeamParticleChangeHook( void )
1349 {
1350 const CexmcEventAction * eventAction(
1351 static_cast< const CexmcEventAction * >( userEventAction ) );
1352 if ( ! eventAction )
1353 throw CexmcException( CexmcWeirdException );
1354
1355 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
1356 eventAction ) );
1357 theEventAction->BeamParticleChangeHook();
1358 }
1359
1360
1361 void CexmcRunManager::SetupConstructionHook( void )
1362 {
1363 #ifdef CEXMC_USE_PERSISTENCY
1364
1365 G4String cmd( "" );
1366 CexmcExceptionType exceptionType( CexmcSystemException );
1367
1368 if ( zipGdmlFile )
1369 {
1370 if ( ProjectIsRead() )
1371 {
1372 cmd = G4String( "bzip2 " ) + projectsDir + "/" + rProject +
1373 gdmlFileExtension;
1374 }
1375 else
1376 {
1377 if ( ProjectIsSaved() )
1378 cmd = G4String( "bzip2 -c " ) + gdmlFileName + " > " +
1379 projectsDir + "/" + projectId + gdmlbz2FileExtension;
1380 }
1381 exceptionType = CexmcFileCompressException;
1382 }
1383 else
1384 {
1385 if ( ! ProjectIsRead() && ProjectIsSaved() )
1386 cmd = G4String( "cp " ) + gdmlFileName + " " + projectsDir + "/" +
1387 projectId + gdmlFileExtension;
1388
1389 }
1390
1391 if ( ! cmd.empty() && system( cmd ) != 0 )
1392 throw CexmcException( exceptionType );
1393 #endif
1394 }
1395
