| 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: CexmcRunManager.cc
30 *
31 * Description: run manager
32 *
33 * Version: 1.0
34 * Created: 03.11.2009 20:27:46
35 * Revision: none
36 * Compiler: gcc
37 *
38 * Author: Alexey Radkov (),
39 * Company: PNPI
40 *
41 * ============================================================================
42 */
43
44 #include <stdlib.h>
45 #include <sys/stat.h>
46 #include <vector>
47 #include <fstream>
48 #ifdef CEXMC_USE_PERSISTENCY
49 #include <boost/archive/binary_iarchive.hpp>
50 #include <boost/archive/binary_oarchive.hpp>
51 #endif
52 #include <G4Eta.hh>
53 #include <G4DigiManager.hh>
54 #include <G4SDManager.hh>
55 #include <G4DecayTable.hh>
56 #include <G4ParticleDefinition.hh>
57 #include <G4ParticleTable.hh>
58 #include <G4UImanager.hh>
59 #include <G4Timer.hh>
60 #include <G4Region.hh>
61 #include <G4RegionStore.hh>
62 #include <G4ProductionCuts.hh>
63 #include <G4VisManager.hh>
64 #include <G4Scene.hh>
65 #include <G4VModel.hh>
66 #include <G4Version.hh>
67 #include "CexmcRunManager.hh"
68 #include "CexmcRunManagerMessenger.hh"
69 #include "CexmcRunAction.hh"
70 #include "CexmcRun.hh"
71 #include "CexmcPhysicsManager.hh"
72 #include "CexmcProductionModel.hh"
73 #include "CexmcSimpleDecayTableStore.hh"
74 #include "CexmcPrimaryGeneratorAction.hh"
75 #include "CexmcEnergyDepositDigitizer.hh"
76 #include "CexmcSimpleEnergyDeposit.hh"
77 #include "CexmcEnergyDepositInLeftRightSet.hh"
78 #include "CexmcEnergyDepositInCalorimeter.hh"
79 #include "CexmcTrackPoints.hh"
80 #include "CexmcTrackPointsInLeftRightSet.hh"
81 #include "CexmcTrackPointsInCalorimeter.hh"
82 #include "CexmcChargeExchangeReconstructor.hh"
83 #include "CexmcEventAction.hh"
84 #include "CexmcParticleGun.hh"
85 #include "CexmcEnergyDepositStore.hh"
86 #include "CexmcTrackPointsStore.hh"
87 #include "CexmcSetup.hh"
88 #include "CexmcEventSObject.hh"
89 #include "CexmcEventFastSObject.hh"
90 #include "CexmcTrackPointInfo.hh"
91 #include "CexmcEventInfo.hh"
92 #include "CexmcBasicPhysicsSettings.hh"
93 #include "CexmcSensitiveDetectorsAttributes.hh"
94 #include "CexmcCustomFilterEval.hh"
95 #include "CexmcScenePrimitives.hh"
96
97
98 namespace
99 {
100 G4String gdmlFileExtension( ".gdml" );
101 G4String gdmlbz2FileExtension( ".gdml.bz2" );
102 }
103
104
105 CexmcRunManager::CexmcRunManager( const G4String & projectId_,
106 const G4String & rProject_,
107 G4bool overrideExistingProject ) :
108 basePhysicsUsed( CexmcPMFactoryInstance::GetBasePhysics() ),
109 productionModelType( CexmcUnknownProductionModel ),
110 gdmlFileName( "default.gdml" ), shouldGdmlFileBeValidated( false ),
111 zipGdmlFile( false ), projectsDir( "." ), projectId( projectId_ ),
112 rProject( rProject_ ), guiMacroName( "" ), cfFileName( "" ),
113 eventCountPolicy( CexmcCountAllEvents ),
114 skipInteractionsWithoutEDTonWrite( true ),
115 evDataVerboseLevel( CexmcWriteEventDataOnEveryEDT ),
116 rEvDataVerboseLevel( CexmcWriteNoEventData ), numberOfEventsProcessed( 0 ),
117 numberOfEventsProcessedEffective( 0 ), curEventRead( 0 ),
118 #ifdef CEXMC_USE_PERSISTENCY
119 eventsArchive( NULL ), fastEventsArchive( NULL ),
120 #ifdef CEXMC_USE_CUSTOM_FILTER
121 customFilter( NULL ),
122 #endif
123 #endif
124 physicsManager( NULL ), messenger( NULL )
125 {
126 /* this exception must be caught before creating the object! */
127 if ( rProject != "" && rProject == projectId )
128 throw CexmcException( CexmcWeirdException );
129
130 const char * projectsDirEnv( std::getenv( "CEXMC_PROJECTS_DIR" ) );
131
132 if ( projectsDirEnv )
133 projectsDir = projectsDirEnv;
134
135 struct stat tmp;
136 if ( ProjectIsSaved() &&
137 stat( ( projectsDir + "/" + projectId + ".rdb" ).c_str(), &tmp ) == 0
138 && ! overrideExistingProject )
139 throw CexmcException( CexmcProjectExists );
140
141 messenger = new CexmcRunManagerMessenger( this );
142
143 #ifdef CEXMC_USE_PERSISTENCY
144 if ( ProjectIsRead() )
145 ReadPreinitProjectData();
146 #endif
147 }
148
149
150 CexmcRunManager::~CexmcRunManager()
151 {
152 #ifdef CEXMC_USE_CUSTOM_FILTER
153 delete customFilter;
154 #endif
155 delete messenger;
156 }
157
158
159 #ifdef CEXMC_USE_PERSISTENCY
160
161 void CexmcRunManager::ReadPreinitProjectData( void )
162 {
163 if ( ! ProjectIsRead() )
164 return;
165
166 /* read run data */
167 std::ifstream runDataFile( ( projectsDir + "/" + rProject + ".rdb" ).
168 c_str() );
169 if ( ! runDataFile )
170 throw CexmcException( CexmcReadProjectIncomplete );
171
172 {
173 boost::archive::binary_iarchive archive( runDataFile );
174 archive >> sObject;
175 }
176
177 basePhysicsUsed = sObject.basePhysicsUsed;
178
179 productionModelType = sObject.productionModelType;
180
181 /* read gdml file */
182 G4String cmd;
183 if ( ProjectIsSaved() )
184 {
185 G4String fileExtension( zipGdmlFile ? gdmlbz2FileExtension :
186 gdmlFileExtension );
187 cmd = G4String( "cp " ) + projectsDir + "/" + rProject + fileExtension +
188 " " + projectsDir + "/" + projectId + fileExtension;
189 if ( system( cmd ) != 0 )
190 throw CexmcException( CexmcReadProjectIncomplete );
191 }
192
193 if ( zipGdmlFile )
194 {
195 cmd = G4String( "bunzip2 " ) + projectsDir + "/" + rProject +
196 gdmlbz2FileExtension;
197 if ( system( cmd ) != 0 )
198 throw CexmcException( CexmcFileCompressException );
199 }
200
201 gdmlFileName = projectsDir + "/" + rProject + gdmlFileExtension;
202 }
203
204
205 void CexmcRunManager::ReadProject( void )
206 {
207 if ( ! ProjectIsRead() )
208 return;
209
210 if ( ! physicsManager )
211 throw CexmcException( CexmcWeirdException );
212
213 physicsManager->GetProductionModel()->SetAngularRanges(
214 sObject.angularRanges );
215 G4DecayTable * etaDecayTable( G4Eta::Definition()->GetDecayTable() );
216 for ( CexmcDecayBranchesStore::const_iterator
217 k( sObject.etaDecayTable.GetDecayBranches().begin() );
218 k != sObject.etaDecayTable.GetDecayBranches().end(); ++k )
219 {
220 etaDecayTable->GetDecayChannel( k->first )->SetBR( k->second );
221 }
222
223 physicsManager->GetProductionModel()->ApplyFermiMotion(
224 sObject.fermiMotionIsOn, false );
225 eventCountPolicy = sObject.eventCountPolicy;
226
227 G4Region * region( G4RegionStore::GetInstance()->GetRegion(
228 CexmcCalorimeterRegionName ) );
229 if ( ! region )
230 throw CexmcException( CexmcCalorimeterRegionNotInitialized );
231
232 region->GetProductionCuts()->SetProductionCuts(
233 sObject.calorimeterRegCuts );
234
235 const CexmcPrimaryGeneratorAction * primaryGeneratorAction(
236 static_cast< const CexmcPrimaryGeneratorAction * >(
237 userPrimaryGeneratorAction ) );
238 CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction(
239 const_cast< CexmcPrimaryGeneratorAction * >(
240 primaryGeneratorAction ) );
241 CexmcParticleGun * particleGun(
242 thePrimaryGeneratorAction->GetParticleGun() );
243 G4ParticleDefinition * particleDefinition(
244 G4ParticleTable::GetParticleTable()->FindParticle(
245 sObject.beamParticle ) );
246 if ( ! particleDefinition )
247 throw CexmcException( CexmcWeirdException );
248
249 particleGun->SetParticleDefinition( particleDefinition );
250 particleGun->SetOrigPosition( sObject.beamPos, false );
251 particleGun->SetOrigDirection( sObject.beamDir, false );
252 particleGun->SetOrigMomentumAmp( sObject.beamMomentumAmp, false );
253
254 thePrimaryGeneratorAction->SetFwhmPosX( sObject.beamFwhmPosX, false );
255 thePrimaryGeneratorAction->SetFwhmPosY( sObject.beamFwhmPosY, false );
256 thePrimaryGeneratorAction->SetFwhmDirX( sObject.beamFwhmDirX, false );
257 thePrimaryGeneratorAction->SetFwhmDirY( sObject.beamFwhmDirY, false );
258 thePrimaryGeneratorAction->SetFwhmMomentumAmp(
259 sObject.beamFwhmMomentumAmp, false );
260
261 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
262 CexmcEnergyDepositDigitizer * edDigitizer(
263 static_cast< CexmcEnergyDepositDigitizer * >(
264 digiManager->FindDigitizerModule( CexmcEDDigitizerName ) ) );
265 if ( ! edDigitizer )
266 throw CexmcException( CexmcWeirdException );
267
268 edDigitizer->SetMonitorThreshold( sObject.monitorEDThreshold, false );
269 edDigitizer->SetVetoCounterLeftThreshold(
270 sObject.vetoCounterEDLeftThreshold, false );
271 edDigitizer->SetVetoCounterRightThreshold(
272 sObject.vetoCounterEDRightThreshold, false );
273 edDigitizer->SetCalorimeterLeftThreshold(
274 sObject.calorimeterEDLeftThreshold, false );
275 edDigitizer->SetCalorimeterRightThreshold(
276 sObject.calorimeterEDRightThreshold, false );
277 edDigitizer->SetCalorimeterTriggerAlgorithm(
278 sObject.calorimeterTriggerAlgorithm, false );
279 edDigitizer->SetOuterCrystalsVetoAlgorithm(
280 sObject.outerCrystalsVetoAlgorithm, false );
281 edDigitizer->SetOuterCrystalsVetoFraction(
282 sObject.outerCrystalsVetoFraction, false );
283 edDigitizer->ApplyFiniteCrystalResolution(
284 sObject.applyFiniteCrystalResolution, false );
285 edDigitizer->SetCrystalResolutionData( sObject.crystalResolutionData );
286
287 const CexmcEventAction * eventAction(
288 static_cast< const CexmcEventAction * >( userEventAction ) );
289 if ( ! eventAction )
290 throw CexmcException( CexmcWeirdException );
291
292 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
293 eventAction ) );
294 CexmcChargeExchangeReconstructor * reconstructor(
295 theEventAction->GetReconstructor() );
296 if ( ! reconstructor )
297 throw CexmcException( CexmcWeirdException );
298
299 reconstructor->SetCalorimeterEntryPointDefinitionAlgorithm(
300 sObject.epDefinitionAlgorithm );
301 reconstructor->SetCalorimeterEntryPointDepthDefinitionAlgorithm(
302 sObject.epDepthDefinitionAlgorithm );
303 reconstructor->SetCrystalSelectionAlgorithm( sObject.csAlgorithm );
304 reconstructor->UseInnerRefCrystal( sObject.useInnerRefCrystal );
305 reconstructor->SetCalorimeterEntryPointDepth( sObject.epDepth );
306 reconstructor->UseTableMass( sObject.useTableMass );
307 reconstructor->UseMassCut( sObject.useMassCut );
308 reconstructor->SetMassCutOPCenter( sObject.mCutOPCenter );
309 reconstructor->SetMassCutNOPCenter( sObject.mCutNOPCenter );
310 reconstructor->SetMassCutOPWidth( sObject.mCutOPWidth );
311 reconstructor->SetMassCutNOPWidth( sObject.mCutNOPWidth );
312 reconstructor->SetMassCutEllipseAngle( sObject.mCutAngle );
313 reconstructor->UseAbsorbedEnergyCut( sObject.useAbsorbedEnergyCut );
314 reconstructor->SetAbsorbedEnergyCutCLCenter( sObject.aeCutCLCenter );
315 reconstructor->SetAbsorbedEnergyCutCRCenter( sObject.aeCutCRCenter );
316 reconstructor->SetAbsorbedEnergyCutCLWidth( sObject.aeCutCLWidth );
317 reconstructor->SetAbsorbedEnergyCutCRWidth( sObject.aeCutCRWidth );
318 reconstructor->SetAbsorbedEnergyCutEllipseAngle( sObject.aeCutAngle );
319 reconstructor->SetExpectedMomentumAmp( sObject.expectedMomentumAmp );
320 reconstructor->SetEDCollectionAlgorithm( sObject.edCollectionAlgorithm );
321
322 physicsManager->SetProposedMaxIL( sObject.proposedMaxIL );
323
324 rEvDataVerboseLevel = sObject.evDataVerboseLevel;
325 evDataVerboseLevel = rEvDataVerboseLevel;
326 }
327
328
329 void CexmcRunManager::SaveProject( void )
330 {
331 if ( ! ProjectIsSaved() )
332 return;
333
334 /* save run data */
335 if ( ! physicsManager )
336 throw CexmcException( CexmcWeirdException );
337
338 CexmcSimpleDecayTableStore etaDecayTable(
339 G4Eta::Definition()->GetDecayTable() );
340 const CexmcPrimaryGeneratorAction * primaryGeneratorAction(
341 static_cast< const CexmcPrimaryGeneratorAction * >(
342 userPrimaryGeneratorAction ) );
343 CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction(
344 const_cast< CexmcPrimaryGeneratorAction * >(
345 primaryGeneratorAction ) );
346 CexmcParticleGun * particleGun(
347 thePrimaryGeneratorAction->GetParticleGun() );
348
349 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
350 CexmcEnergyDepositDigitizer * edDigitizer(
351 static_cast< CexmcEnergyDepositDigitizer * >(
352 digiManager->FindDigitizerModule( CexmcEDDigitizerName ) ) );
353 if ( ! edDigitizer )
354 throw CexmcException( CexmcWeirdException );
355
356 const CexmcEventAction * eventAction(
357 static_cast< const CexmcEventAction * >( userEventAction ) );
358 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
359 eventAction ) );
360 CexmcChargeExchangeReconstructor * reconstructor(
361 theEventAction->GetReconstructor() );
362
363 std::vector< G4double > calorimeterRegCuts( 4 );
364 if ( ProjectIsRead() )
365 {
366 calorimeterRegCuts = sObject.calorimeterRegCuts;
367 }
368 else
369 {
370 G4Region * region( G4RegionStore::GetInstance()->GetRegion(
371 CexmcCalorimeterRegionName ) );
372 if ( ! region )
373 throw CexmcException( CexmcCalorimeterRegionNotInitialized );
374
375 calorimeterRegCuts = region->GetProductionCuts()->GetProductionCuts();
376 }
377
378 CexmcNmbOfHitsInRanges nmbOfHitsSampled;
379 CexmcNmbOfHitsInRanges nmbOfHitsSampledFull;
380 CexmcNmbOfHitsInRanges nmbOfHitsTriggeredRealRange;
381 CexmcNmbOfHitsInRanges nmbOfHitsTriggeredRecRange;
382 CexmcNmbOfHitsInRanges nmbOfOrphanHits;
383 G4int nmbOfFalseHitsTriggeredEDT( 0 );
384 G4int nmbOfFalseHitsTriggeredRec( 0 );
385 G4int nmbOfSavedEvents( 0 );
386 G4int nmbOfSavedFastEvents( 0 );
387 CexmcRun * run( static_cast< CexmcRun * >( currentRun ) );
388
389 if ( run )
390 {
391 nmbOfHitsSampled = run->GetNmbOfHitsSampled();
392 nmbOfHitsSampledFull = run->GetNmbOfHitsSampledFull();
393 nmbOfHitsTriggeredRealRange = run->GetNmbOfHitsTriggeredRealRange();
394 nmbOfHitsTriggeredRecRange = run->GetNmbOfHitsTriggeredRecRange();
395 nmbOfOrphanHits = run->GetNmbOfOrphanHits();
396 nmbOfFalseHitsTriggeredEDT = run->GetNmbOfFalseHitsTriggeredEDT();
397 nmbOfFalseHitsTriggeredRec = run->GetNmbOfFalseHitsTriggeredRec();
398 nmbOfSavedEvents = run->GetNmbOfSavedEvents();
399 nmbOfSavedFastEvents = run->GetNmbOfSavedFastEvents();
400 }
401
402 CexmcRunSObject sObjectToWrite = {
403 basePhysicsUsed, productionModelType, gdmlFileName, etaDecayTable,
404 physicsManager->GetProductionModel()->GetAngularRanges(),
405 physicsManager->GetProductionModel()->IsFermiMotionOn(),
406 calorimeterRegCuts, eventCountPolicy,
407 particleGun->GetParticleDefinition()->GetParticleName(),
408 particleGun->GetOrigPosition(), particleGun->GetOrigDirection(),
409 particleGun->GetOrigMomentumAmp(),
410 primaryGeneratorAction->GetFwhmPosX(),
411 primaryGeneratorAction->GetFwhmPosY(),
412 primaryGeneratorAction->GetFwhmDirX(),
413 primaryGeneratorAction->GetFwhmDirY(),
414 primaryGeneratorAction->GetFwhmMomentumAmp(),
415 edDigitizer->GetMonitorThreshold(),
416 edDigitizer->GetVetoCounterLeftThreshold(),
417 edDigitizer->GetVetoCounterRightThreshold(),
418 edDigitizer->GetCalorimeterLeftThreshold(),
419 edDigitizer->GetCalorimeterRightThreshold(),
420 edDigitizer->GetCalorimeterTriggerAlgorithm(),
421 edDigitizer->GetOuterCrystalsVetoAlgorithm(),
422 edDigitizer->GetOuterCrystalsVetoFraction(),
423 edDigitizer->IsFiniteCrystalResolutionApplied(),
424 edDigitizer->GetCrystalResolutionData(),
425 reconstructor->GetCalorimeterEntryPointDefinitionAlgorithm(),
426 reconstructor->GetCalorimeterEntryPointDepthDefinitionAlgorithm(),
427 reconstructor->GetCrystalSelectionAlgorithm(),
428 reconstructor->IsInnerRefCrystalUsed(),
429 reconstructor->GetCalorimeterEntryPointDepth(),
430 reconstructor->IsTableMassUsed(), reconstructor->IsMassCutUsed(),
431 reconstructor->GetMassCutOPCenter(),
432 reconstructor->GetMassCutNOPCenter(),
433 reconstructor->GetMassCutOPWidth(), reconstructor->GetMassCutNOPWidth(),
434 reconstructor->GetMassCutEllipseAngle(),
435 reconstructor->IsAbsorbedEnergyCutUsed(),
436 reconstructor->GetAbsorbedEnergyCutCLCenter(),
437 reconstructor->GetAbsorbedEnergyCutCRCenter(),
438 reconstructor->GetAbsorbedEnergyCutCLWidth(),
439 reconstructor->GetAbsorbedEnergyCutCRWidth(),
440 reconstructor->GetAbsorbedEnergyCutEllipseAngle(),
441 nmbOfHitsSampled, nmbOfHitsSampledFull, nmbOfHitsTriggeredRealRange,
442 nmbOfHitsTriggeredRecRange, nmbOfOrphanHits, nmbOfFalseHitsTriggeredEDT,
443 nmbOfFalseHitsTriggeredRec, nmbOfSavedEvents, nmbOfSavedFastEvents,
444 numberOfEventsProcessed, numberOfEventsProcessedEffective,
445 numberOfEventToBeProcessed, rProject, skipInteractionsWithoutEDTonWrite,
446 cfFileName, evDataVerboseLevel, physicsManager->GetProposedMaxIL(),
447 reconstructor->GetExpectedMomentumAmp(),
448 reconstructor->GetEDCollectionAlgorithm(), 0 };
449
450 std::ofstream runDataFile( ( projectsDir + "/" + projectId + ".rdb" ).
451 c_str() );
452
453 {
454 boost::archive::binary_oarchive archive( runDataFile );
455 archive << sObjectToWrite;
456 }
457 }
458
459 #endif
460
461
462 void CexmcRunManager::DoCommonEventLoop( G4int nEvent, const G4String & cmd,
463 G4int nSelect )
464 {
465 G4int iEvent( 0 );
466 G4int iEventEffective( 0 );
467
468 for ( iEvent = 0; iEventEffective < nEvent; ++iEvent )
469 {
470 currentEvent = GenerateEvent( iEvent );
471 eventManager->ProcessOneEvent( currentEvent );
472 CexmcEventInfo * eventInfo( static_cast< CexmcEventInfo * >(
473 currentEvent->GetUserInformation() ) );
474 switch ( eventCountPolicy )
475 {
476 case CexmcCountAllEvents :
477 ++iEventEffective;
478 break;
479 case CexmcCountEventsWithInteraction :
480 if ( eventInfo->TpTriggerIsOk() )
481 ++iEventEffective;
482 break;
483 case CexmcCountEventsWithTrigger :
484 if ( eventInfo->EdTriggerIsOk() )
485 ++iEventEffective;
486 break;
487 default :
488 ++iEventEffective;
489 break;
490 }
491 AnalyzeEvent( currentEvent );
492 UpdateScoring();
493 if ( iEvent < nSelect )
494 G4UImanager::GetUIpointer()->ApplyCommand( cmd );
495 StackPreviousEvent( currentEvent );
496 currentEvent = 0;
497 if ( runAborted )
498 break;
499 }
500
501 numberOfEventsProcessed = iEvent;
502 numberOfEventsProcessedEffective = iEventEffective;
503 }
504
505
506 #ifdef CEXMC_USE_PERSISTENCY
507
508 void CexmcRunManager::DoReadEventLoop( G4int nEvent )
509 {
510 G4int iEvent( 0 );
511 G4int iEventEffective( 0 );
512 G4int nEventCount( 0 );
513
514 if ( ! ProjectIsRead() )
515 return;
516
517 if ( ! physicsManager )
518 throw CexmcException( CexmcWeirdException );
519
520 CexmcProductionModel * productionModel(
521 physicsManager->GetProductionModel() );
522 if ( ! productionModel )
523 throw CexmcException( CexmcWeirdException );
524
525 CexmcSetup * setup( static_cast< CexmcSetup * >( userDetector ) );
526 if ( ! setup )
527 throw CexmcException( CexmcWeirdException );
528
529 CexmcEventSObject evSObject;
530 CexmcEventFastSObject evFastSObject;
531
532 /* read events data */
533 std::ifstream eventsDataFile(
534 ( projectsDir + "/" + rProject + ".edb" ).c_str() );
535 if ( ! eventsDataFile )
536 throw CexmcException( CexmcReadProjectIncomplete );
537
538 boost::archive::binary_iarchive evArchive( eventsDataFile );
539
540 std::ifstream eventsFastDataFile(
541 ( projectsDir + "/" + rProject + ".fdb" ).c_str() );
542 if ( ! eventsFastDataFile )
543 throw CexmcException( CexmcReadProjectIncomplete );
544
545 boost::archive::binary_iarchive evFastArchive( eventsFastDataFile );
546
547 G4Event event;
548 currentEvent = &event;
549 G4SDManager * sdManager( G4SDManager::GetSDMpointer() );
550 event.SetHCofThisEvent( sdManager->PrepareNewEvent() );
551 G4HCofThisEvent * hcOfThisEvent( event.GetHCofThisEvent() );
552
553 G4DigiManager * digiManager( G4DigiManager::GetDMpointer() );
554
555 G4int hcId( digiManager->GetHitsCollectionID(
556 CexmcDetectorRoleName[ CexmcMonitorDetectorRole ] +
557 "/" + CexmcDetectorTypeName[ CexmcEDDetector ] ) );
558 CexmcEnergyDepositCollection * monitorED(
559 new CexmcEnergyDepositCollection );
560 hcOfThisEvent->AddHitsCollection( hcId, monitorED );
561 hcId = digiManager->GetHitsCollectionID(
562 CexmcDetectorRoleName[ CexmcVetoCounterDetectorRole ] +
563 "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
564 CexmcEnergyDepositCollection * vetoCounterED(
565 new CexmcEnergyDepositCollection );
566 hcOfThisEvent->AddHitsCollection( hcId, vetoCounterED );
567 hcId = digiManager->GetHitsCollectionID(
568 CexmcDetectorRoleName[ CexmcCalorimeterDetectorRole ] +
569 "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
570 CexmcEnergyDepositCollection * calorimeterED(
571 new CexmcEnergyDepositCollection );
572 hcOfThisEvent->AddHitsCollection( hcId, calorimeterED );
573 hcId = digiManager->GetHitsCollectionID(
574 CexmcDetectorRoleName[ CexmcMonitorDetectorRole ] +
575 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
576 CexmcTrackPointsCollection * monitorTP( new CexmcTrackPointsCollection );
577 hcOfThisEvent->AddHitsCollection( hcId, monitorTP );
578 hcId = digiManager->GetHitsCollectionID(
579 CexmcDetectorRoleName[ CexmcVetoCounterDetectorRole ] +
580 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
581 CexmcTrackPointsCollection * vetoCounterTP(
582 new CexmcTrackPointsCollection );
583 hcOfThisEvent->AddHitsCollection( hcId, vetoCounterTP );
584 hcId = digiManager->GetHitsCollectionID(
585 CexmcDetectorRoleName[ CexmcCalorimeterDetectorRole ] +
586 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
587 CexmcTrackPointsCollection * calorimeterTP(
588 new CexmcTrackPointsCollection );
589 hcOfThisEvent->AddHitsCollection( hcId, calorimeterTP );
590 hcId = digiManager->GetHitsCollectionID(
591 CexmcDetectorRoleName[ CexmcTargetDetectorRole ] +
592 "/" + CexmcDetectorTypeName[ CexmcTPDetector ] );
593 CexmcTrackPointsCollection * targetTP( new CexmcTrackPointsCollection );
594 hcOfThisEvent->AddHitsCollection( hcId, targetTP );
595
596 #ifdef CEXMC_USE_CUSTOM_FILTER
597 if ( customFilter )
598 customFilter->SetAddressedData( &evFastSObject, &evSObject );
599 #endif
600
601 G4int nmbOfSavedEvents( rEvDataVerboseLevel == CexmcWriteNoEventData ? 0 :
602 sObject.nmbOfSavedFastEvents );
603 G4bool eventDataWrittenOnEveryTPT( rEvDataVerboseLevel ==
604 CexmcWriteEventDataOnEveryTPT );
605
606 for ( int i( 0 ); i < nmbOfSavedEvents; ++i )
607 {
608 evFastArchive >> evFastSObject;
609
610 if ( nEventCount < curEventRead )
611 {
612 if ( eventDataWrittenOnEveryTPT ||
613 evFastSObject.edDigitizerHasTriggered )
614 {
615 evArchive >> evSObject;
616 if ( evFastSObject.edDigitizerHasTriggered )
617 ++nEventCount;
618 }
619 continue;
620 }
621
622 ++iEvent;
623
624 productionModel->SetTriggeredAngularRanges(
625 evFastSObject.opCosThetaSCM );
626 const CexmcAngularRangeList & triggeredAngularRanges(
627 productionModel->GetTriggeredAngularRanges() );
628
629 if ( ! eventDataWrittenOnEveryTPT &&
630 ! evFastSObject.edDigitizerHasTriggered )
631 {
632 #ifdef CEXMC_USE_CUSTOM_FILTER
633 /* user must be aware that using tpt commands in custom filter
634 * scripts for poor event data sets can easily lead to logical
635 * errors! This is because most of tpt data is only available for
636 * events with EDT trigger. There is no such problem if event data
637 * was written on every TPT event. */
638 if ( customFilter && ! customFilter->EvalTPT() )
639 continue;
640 #endif
641 SaveCurrentTPTEvent( evFastSObject, triggeredAngularRanges,
642 ProjectIsSaved() && ! skipInteractionsWithoutEDTonWrite );
643 continue;
644 }
645
646 evArchive >> evSObject;
647
648 G4bool skipEDTOnThisEvent( false );
649
650 #ifdef CEXMC_USE_CUSTOM_FILTER
651 if ( customFilter && ! customFilter->EvalTPT() )
652 continue;
653 if ( customFilter && ! customFilter->EvalEDT() )
654 {
655 if ( ! eventDataWrittenOnEveryTPT )
656 {
657 SaveCurrentTPTEvent( evFastSObject, triggeredAngularRanges,
658 ProjectIsSaved() );
659 continue;
660 }
661 skipEDTOnThisEvent = true;
662 }
663 #endif
664
665 event.SetEventID( evSObject.eventId );
666
667 /* AAA: beware! If something throws an exception between AAA and CCC
668 * then there would be segmentation violation in ~THitsMap() as it
669 * would try to delete local variable evSObject's fields.
670 * See also BBB */
671
672 monitorED->GetMap()->operator[]( 0 ) = &evSObject.monitorED;
673 vetoCounterED->GetMap()->operator[]( 0 ) = &evSObject.vetoCounterEDLeft;
674 vetoCounterED->GetMap()->operator[]( 1 <<
675 CexmcEnergyDepositInLeftRightSet::GetLeftRightBitsOffset() ) =
676 &evSObject.vetoCounterEDRight;
677 G4int row( 0 );
678 G4int column( 0 );
679 for ( CexmcEnergyDepositCalorimeterCollection::iterator
680 k( evSObject.calorimeterEDLeftCollection.begin() );
681 k != evSObject.calorimeterEDLeftCollection.end(); ++k )
682 {
683 G4int index( row <<
684 CexmcEnergyDepositInCalorimeter::GetCopyDepth1BitsOffset() );
685 column = 0;
686 for ( CexmcEnergyDepositCrystalRowCollection::iterator
687 l( k->begin() ); l != k->end(); ++l )
688 {
689 calorimeterED->GetMap()->operator[]( index | column ) = &*l;
690 ++column;
691 }
692 ++row;
693 }
694 row = 0;
695 for ( CexmcEnergyDepositCalorimeterCollection::iterator
696 k( evSObject.calorimeterEDRightCollection.begin() );
697 k != evSObject.calorimeterEDRightCollection.end(); ++k )
698 {
699 G4int index(
700 1 << CexmcEnergyDepositInLeftRightSet::GetLeftRightBitsOffset()
701 | row <<
702 CexmcEnergyDepositInCalorimeter::GetCopyDepth1BitsOffset() );
703 column = 0;
704 for ( CexmcEnergyDepositCrystalRowCollection::iterator
705 l( k->begin() ); l != k->end(); ++l )
706 {
707 calorimeterED->GetMap()->operator[]( index | column ) = &*l;
708 ++column;
709 }
710 ++row;
711 }
712
713 CexmcTrackPointInfo monitorTPInfo( evSObject.monitorTP );
714 CexmcTrackPointInfo targetTPBeamParticleInfo(
715 evSObject.targetTPBeamParticle );
716 CexmcTrackPointInfo targetTPOutputParticleInfo(
717 evSObject.targetTPOutputParticle );
718 CexmcTrackPointInfo targetTPNucleusParticleInfo(
719 evSObject.targetTPNucleusParticle );
720 CexmcTrackPointInfo targetTPOutputParticleDecayProductParticle1Info(
721 evSObject.targetTPOutputParticleDecayProductParticle1 );
722 CexmcTrackPointInfo targetTPOutputParticleDecayProductParticle2Info(
723 evSObject.targetTPOutputParticleDecayProductParticle2 );
724 CexmcTrackPointInfo vetoCounterTPLeftInfo(
725 evSObject.vetoCounterTPLeft );
726 CexmcTrackPointInfo vetoCounterTPRightInfo(
727 evSObject.vetoCounterTPRight );
728 CexmcTrackPointInfo calorimeterTPLeftInfo(
729 evSObject.calorimeterTPLeft );
730 CexmcTrackPointInfo calorimeterTPRightInfo(
731 evSObject.calorimeterTPRight );
732
733 if ( monitorTPInfo.IsValid() )
734 monitorTP->GetMap()->operator[]( monitorTPInfo.trackId ) =
735 &monitorTPInfo;
736 if ( targetTPBeamParticleInfo.IsValid() )
737 targetTP->GetMap()->operator[](
738 targetTPBeamParticleInfo.trackId ) =
739 &targetTPBeamParticleInfo;
740 if ( targetTPOutputParticleInfo.IsValid() )
741 targetTP->GetMap()->operator[](
742 targetTPOutputParticleInfo.trackId ) =
743 &targetTPOutputParticleInfo;
744 if ( targetTPNucleusParticleInfo.IsValid() )
745 targetTP->GetMap()->operator[](
746 targetTPNucleusParticleInfo.trackId ) =
747 &targetTPNucleusParticleInfo;
748 if ( targetTPOutputParticleDecayProductParticle1Info.IsValid() )
749 targetTP->GetMap()->operator[](
750 targetTPOutputParticleDecayProductParticle1Info.trackId ) =
751 &targetTPOutputParticleDecayProductParticle1Info;
752 if ( targetTPOutputParticleDecayProductParticle2Info.IsValid() )
753 targetTP->GetMap()->operator[](
754 targetTPOutputParticleDecayProductParticle2Info.trackId ) =
755 &targetTPOutputParticleDecayProductParticle2Info;
756 if ( vetoCounterTPLeftInfo.IsValid() )
757 vetoCounterTP->GetMap()->operator[](
758 vetoCounterTPLeftInfo.trackId ) = &vetoCounterTPLeftInfo;
759 if ( vetoCounterTPRightInfo.IsValid() )
760 vetoCounterTP->GetMap()->operator[](
761 1 << CexmcTrackPointsInLeftRightSet::GetLeftRightBitsOffset() |
762 vetoCounterTPRightInfo.trackId ) = &vetoCounterTPRightInfo;
763
764 G4ThreeVector pos;
765 if ( calorimeterTPLeftInfo.IsValid() )
766 {
767 pos = calorimeterTPLeftInfo.positionLocal;
768 setup->ConvertToCrystalGeometry(
769 calorimeterTPLeftInfo.positionLocal, row, column, pos );
770 calorimeterTPLeftInfo.positionLocal = pos;
771 calorimeterTP->GetMap()->operator[](
772 row << CexmcTrackPointsInCalorimeter::
773 GetCopyDepth1BitsOffset() |
774 column << CexmcTrackPointsInCalorimeter::
775 GetCopyDepth0BitsOffset() |
776 calorimeterTPLeftInfo.trackId ) = &calorimeterTPLeftInfo;
777 }
778 if ( calorimeterTPRightInfo.IsValid() )
779 {
780 pos = calorimeterTPRightInfo.positionLocal;
781 setup->ConvertToCrystalGeometry(
782 calorimeterTPRightInfo.positionLocal, row, column, pos );
783 calorimeterTPRightInfo.positionLocal = pos;
784 calorimeterTP->GetMap()->operator[](
785 1 << CexmcTrackPointsInLeftRightSet::GetLeftRightBitsOffset() |
786 row << CexmcTrackPointsInCalorimeter::
787 GetCopyDepth1BitsOffset() |
788 column << CexmcTrackPointsInCalorimeter::
789 GetCopyDepth0BitsOffset() |
790 calorimeterTPRightInfo.trackId ) = &calorimeterTPRightInfo;
791 }
792
793 productionModel->SetProductionModelData(
794 evSObject.productionModelData );
795
796 const CexmcEventAction * eventAction(
797 static_cast< const CexmcEventAction * >( userEventAction ) );
798 if ( ! eventAction )
799 {
800 /* BBB: all hits collections must be cleared before throwing
801 * anything from here, otherwise ~THitsMap() will try to delete
802 * local variable evSObject's fields like monitorED etc. */
803 monitorED->GetMap()->clear();
804 vetoCounterED->GetMap()->clear();
805 calorimeterED->GetMap()->clear();
806 monitorTP->GetMap()->clear();
807 targetTP->GetMap()->clear();
808 vetoCounterTP->GetMap()->clear();
809 calorimeterTP->GetMap()->clear();
810 throw CexmcException( CexmcEventActionIsNotInitialized );
811 }
812
813 if ( skipEDTOnThisEvent )
814 event.SetUserInformation( new CexmcEventInfo( false, false,
815 false ) );
816
817 CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
818 eventAction ) );
819 theEventAction->EndOfEventAction( &event );
820
821 CexmcEventInfo * eventInfo( static_cast< CexmcEventInfo * >(
822 event.GetUserInformation() ) );
823
824 if ( eventInfo->EdTriggerIsOk() )
825 ++iEventEffective;
826
827 delete eventInfo;
828 event.SetUserInformation( NULL );
829
830 monitorED->GetMap()->clear();
831 vetoCounterED->GetMap()->clear();
832 calorimeterED->GetMap()->clear();
833 monitorTP->GetMap()->clear();
834 targetTP->GetMap()->clear();
835 vetoCounterTP->GetMap()->clear();
836 calorimeterTP->GetMap()->clear();
837
838 /* CCC: see AAA */
839
840 if ( nEvent > 0 && iEventEffective == nEvent )
841 break;
842 }
843
844 curEventRead = nEventCount + iEventEffective;
845
846 numberOfEventsProcessed = iEvent;
847 numberOfEventsProcessedEffective = iEventEffective;
848
849 #ifdef CEXMC_USE_CUSTOM_FILTER
850 if ( customFilter )
851 customFilter->SetAddressedData( NULL, NULL );
852 #endif
853 }
854
855
856 void CexmcRunManager::SaveCurrentTPTEvent(
857 const CexmcEventFastSObject & evFastSObject,
858 const CexmcAngularRangeList & angularRanges,
859 G4bool writeToDatabase )
860 {
861 CexmcRun * run( static_cast< CexmcRun * >( currentRun ) );
862
863 if ( ! run )
864 return;
865
866 for ( CexmcAngularRangeList::const_iterator k( angularRanges.begin() );
867 k != angularRanges.end(); ++k )
868 {
869 run->IncrementNmbOfHitsSampledFull( k->index );
870 if ( evFastSObject.edDigitizerMonitorHasTriggered )
871 run->IncrementNmbOfHitsSampled( k->index );
872 }
873
874 if ( writeToDatabase )
875 {
876 fastEventsArchive->operator<<( evFastSObject );
877 run->IncrementNmbOfSavedFastEvents();
878 }
879 }
880
881 #endif
882
883
884 /* mostly adopted from G4RunManager::DoEventLoop() */
885 void CexmcRunManager::DoEventLoop( G4int nEvent, const char * macroFile,
886 G4int nSelect )
887 {
888 if ( verboseLevel > 0 )
889 timer->Start();
890
891 G4String cmd;
892 if ( macroFile != 0 )
893 {
894 if ( nSelect < 0 )
895 nSelect = nEvent;
896 cmd = "/control/execute ";
897 cmd += macroFile;
898 }
899 else
900 {
901 nSelect = -1;
902 }
903
904 numberOfEventsProcessed = 0;
905 numberOfEventsProcessedEffective = 0;
906
907 #ifdef CEXMC_USE_PERSISTENCY
908 eventsArchive = NULL;
909 fastEventsArchive = NULL;
910 if ( ProjectIsRead() )
911 {
912 if ( ProjectIsSaved() )
913 {
914 std::ofstream eventsDataFile(
915 ( projectsDir + "/" + projectId + ".edb" ).c_str() );
916 boost::archive::binary_oarchive eventsArchive_( eventsDataFile );
917 std::ofstream fastEventsDataFile(
918 ( projectsDir + "/" + projectId + ".fdb" ).c_str() );
919 boost::archive::binary_oarchive fastEventsArchive_(
920 fastEventsDataFile );
921 eventsArchive = &eventsArchive_;
922 fastEventsArchive = &fastEventsArchive_;
923 DoReadEventLoop( nEvent );
924 }
925 else
926 {
927 DoReadEventLoop( nEvent );
928 }
929 }
930 else
931 {
932 if ( ProjectIsSaved() )
933 {
934 std::ofstream eventsDataFile(
935 ( projectsDir + "/" + projectId + ".edb" ).c_str() );
936 boost::archive::binary_oarchive eventsArchive_( eventsDataFile );
937 std::ofstream fastEventsDataFile(
938 ( projectsDir + "/" + projectId + ".fdb" ).c_str() );
939 boost::archive::binary_oarchive fastEventsArchive_(
940 fastEventsDataFile );
941 eventsArchive = &eventsArchive_;
942 fastEventsArchive = &fastEventsArchive_;
943 DoCommonEventLoop( nEvent, cmd, nSelect );
944 }
945 else
946 {
947 DoCommonEventLoop( nEvent, cmd, nSelect );
948 }
949 }
950 eventsArchive = NULL;
951 fastEventsArchive = NULL;
952 #else
953 DoCommonEventLoop( nEvent, cmd, nSelect );
954 #endif
955
956 if ( verboseLevel > 0 )
957 {
958 timer->Stop();
959 G4cout << "Run terminated." << G4endl;
960 G4cout << "Run Summary" << G4endl;
961 if ( runAborted )
962 {
963 G4cout << " Run Aborted after " << numberOfEventsProcessed <<
964 " events processed." << G4endl;
965 }
966 else
967 {
968 G4cout << " Number of events processed : " <<
969 numberOfEventsProcessed << ", effectively: " <<
970 numberOfEventsProcessedEffective << G4endl;
971 }
972 G4cout << " " << *timer << G4endl;
973 }
974 }
975
976
977 #ifdef CEXMC_USE_PERSISTENCY
978
979 void CexmcRunManager::PrintReadRunData( void ) const
980 {
981 if ( ! ProjectIsRead() )
982 return;
983
984 G4bool refCrystalInfoPrinted( false );
985
986 G4cout << CEXMC_LINE_START << "Run data read from project '" << rProject <<
987 "'" << G4endl;
988 G4cout << " (archive class version " <<
989 sObject.actualVersion << ")" << G4endl;
990 if ( ! sObject.rProject.empty() )
991 {
992 G4cout << " -- Based on project '" << sObject.rProject << "'" <<
993 G4endl;
994 if ( ! sObject.cfFileName.empty() )
995 G4cout << " -- Custom filter script '" << sObject.cfFileName <<
996 "' was used" << G4endl;
997 }
998 G4cout << " -- Event data verbose level (0 - not saved, 1 - triggers, "
999 "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 /* read events data */
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 /* should not get here */
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 /* G4Scene declarations lack this kind of typedef */
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 /* BEWARE: looks like G4Scene won't delete models from its lists upon
1343 * termination! Hence destructor of the new model won't be called */
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 /* save gdml file */
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 /* else already saved in ReadPreinitProjectData() */
1389 }
1390
1391 if ( ! cmd.empty() && system( cmd ) != 0 )
1392 throw CexmcException( exceptionType );
1393 #endif
1394 }
1395
1396