| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 G4Stri
106 const G4Stri
107 G4bool over
108 basePhysicsUsed( CexmcPMFactoryInstance::G
109 productionModelType( CexmcUnknownProductio
110 gdmlFileName( "default.gdml" ), shouldGdml
111 zipGdmlFile( false ), projectsDir( "." ),
112 rProject( rProject_ ), guiMacroName( "" ),
113 eventCountPolicy( CexmcCountAllEvents ),
114 skipInteractionsWithoutEDTonWrite( true ),
115 evDataVerboseLevel( CexmcWriteEventDataOnE
116 rEvDataVerboseLevel( CexmcWriteNoEventData
117 numberOfEventsProcessedEffective( 0 ), cur
118 #ifdef CEXMC_USE_PERSISTENCY
119 eventsArchive( NULL ), fastEventsArchive(
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 cr
127 if ( rProject != "" && rProject == project
128 throw CexmcException( CexmcWeirdExcept
129
130 const char * projectsDirEnv( std::getenv(
131
132 if ( projectsDirEnv )
133 projectsDir = projectsDirEnv;
134
135 struct stat tmp;
136 if ( ProjectIsSaved() &&
137 stat( ( projectsDir + "/" + projectId
138 && ! overrideExistingProject )
139 throw CexmcException( CexmcProjectExis
140
141 messenger = new CexmcRunManagerMessenger(
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(
162 {
163 if ( ! ProjectIsRead() )
164 return;
165
166 /* read run data */
167 std::ifstream runDataFile( ( projectsDir
168 c_str() );
169 if ( ! runDataFile )
170 throw CexmcException( CexmcReadProject
171
172 {
173 boost::archive::binary_iarchive archi
174 archive >> sObject;
175 }
176
177 basePhysicsUsed = sObject.basePhysicsUsed;
178
179 productionModelType = sObject.productionMo
180
181 /* read gdml file */
182 G4String cmd;
183 if ( ProjectIsSaved() )
184 {
185 G4String fileExtension( zipGdmlFile ?
186
187 cmd = G4String( "cp " ) + projectsDir
188 " " + projectsDir
189 if ( system( cmd ) != 0 )
190 throw CexmcException( CexmcReadPro
191 }
192
193 if ( zipGdmlFile )
194 {
195 cmd = G4String( "bunzip2 " ) + project
196
197 if ( system( cmd ) != 0 )
198 throw CexmcException( CexmcFileCom
199 }
200
201 gdmlFileName = projectsDir + "/" + rProjec
202 }
203
204
205 void CexmcRunManager::ReadProject( void )
206 {
207 if ( ! ProjectIsRead() )
208 return;
209
210 if ( ! physicsManager )
211 throw CexmcException( CexmcWeirdExcept
212
213 physicsManager->GetProductionModel()->SetA
214
215 G4DecayTable * etaDecayTable( G4Eta::Defi
216 for ( CexmcDecayBranchesStore::const_itera
217 k( sObject.etaDecayTable.GetDecayB
218 k != sObject.etaDecayTable.GetDeca
219 {
220 etaDecayTable->GetDecayChannel( k->fir
221 }
222
223 physicsManager->GetProductionModel()->Appl
224 sO
225 eventCountPolicy = sObject.eventCountPolic
226
227 G4Region * region( G4RegionStore::GetInst
228
229 if ( ! region )
230 throw CexmcException( CexmcCalorimeter
231
232 region->GetProductionCuts()->SetProduction
233
234
235 const CexmcPrimaryGeneratorAction * prima
236 static_cast< const
237
238 CexmcPrimaryGeneratorAction * thePr
239 const_cast< CexmcP
240 pr
241 CexmcParticleGun * particleGun(
242 thePrimaryGene
243 G4ParticleDefinition * particleDefinition
244 G4ParticleTable::GetPartic
245
246 if ( ! particleDefinition )
247 throw CexmcException( CexmcWeirdExcept
248
249 particleGun->SetParticleDefinition( partic
250 particleGun->SetOrigPosition( sObject.beam
251 particleGun->SetOrigDirection( sObject.bea
252 particleGun->SetOrigMomentumAmp( sObject.b
253
254 thePrimaryGeneratorAction->SetFwhmPosX( sO
255 thePrimaryGeneratorAction->SetFwhmPosY( sO
256 thePrimaryGeneratorAction->SetFwhmDirX( sO
257 thePrimaryGeneratorAction->SetFwhmDirY( sO
258 thePrimaryGeneratorAction->SetFwhmMomentum
259 sObjec
260
261 G4DigiManager * digiManager
262 CexmcEnergyDepositDigitizer * edDigitizer
263 static_cast< CexmcEnergyDepositDig
264 digiManager->FindDigitizerModu
265 if ( ! edDigitizer )
266 throw CexmcException( CexmcWeirdExcept
267
268 edDigitizer->SetMonitorThreshold( sObject.
269 edDigitizer->SetVetoCounterLeftThreshold(
270 sObject.vetoCo
271 edDigitizer->SetVetoCounterRightThreshold(
272 sObject.vetoCo
273 edDigitizer->SetCalorimeterLeftThreshold(
274 sObject.calori
275 edDigitizer->SetCalorimeterRightThreshold(
276 sObject.calori
277 edDigitizer->SetCalorimeterTriggerAlgorith
278 sObject.calori
279 edDigitizer->SetOuterCrystalsVetoAlgorithm
280 sObject.outerC
281 edDigitizer->SetOuterCrystalsVetoFraction(
282 sObject.outerC
283 edDigitizer->ApplyFiniteCrystalResolution(
284 sObject.applyF
285 edDigitizer->SetCrystalResolutionData( sOb
286
287 const CexmcEventAction * eventAction(
288 static_cast< const CexmcEventA
289 if ( ! eventAction )
290 throw CexmcException( CexmcWeirdExcept
291
292 CexmcEventAction * theEventAction(
293
294 CexmcChargeExchangeReconstructor * recons
295 theEve
296 if ( ! reconstructor )
297 throw CexmcException( CexmcWeirdExcept
298
299 reconstructor->SetCalorimeterEntryPointDef
300 sObjec
301 reconstructor->SetCalorimeterEntryPointDep
302 sObjec
303 reconstructor->SetCrystalSelectionAlgorith
304 reconstructor->UseInnerRefCrystal( sObject
305 reconstructor->SetCalorimeterEntryPointDep
306 reconstructor->UseTableMass( sObject.useTa
307 reconstructor->UseMassCut( sObject.useMass
308 reconstructor->SetMassCutOPCenter( sObject
309 reconstructor->SetMassCutNOPCenter( sObjec
310 reconstructor->SetMassCutOPWidth( sObject.
311 reconstructor->SetMassCutNOPWidth( sObject
312 reconstructor->SetMassCutEllipseAngle( sOb
313 reconstructor->UseAbsorbedEnergyCut( sObje
314 reconstructor->SetAbsorbedEnergyCutCLCente
315 reconstructor->SetAbsorbedEnergyCutCRCente
316 reconstructor->SetAbsorbedEnergyCutCLWidth
317 reconstructor->SetAbsorbedEnergyCutCRWidth
318 reconstructor->SetAbsorbedEnergyCutEllipse
319 reconstructor->SetExpectedMomentumAmp( sOb
320 reconstructor->SetEDCollectionAlgorithm( s
321
322 physicsManager->SetProposedMaxIL( sObject.
323
324 rEvDataVerboseLevel = sObject.evDataVerbos
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( CexmcWeirdExcept
337
338 CexmcSimpleDecayTableStore etaDecayTable(
339 G4Eta::Def
340 const CexmcPrimaryGeneratorAction * prima
341 static_cast< const
342
343 CexmcPrimaryGeneratorAction * thePrimaryG
344 const_cast< CexmcP
345 pr
346 CexmcParticleGun * particleGun(
347 thePrimaryGenerato
348
349 G4DigiManager * digiManager
350 CexmcEnergyDepositDigitizer * edDigitizer
351 static_cast< CexmcEnergyDepositDig
352 digiManager->FindDigitizerModu
353 if ( ! edDigitizer )
354 throw CexmcException( CexmcWeirdExcept
355
356 const CexmcEventAction * eventAction(
357 static_cast< const CexmcEventA
358 CexmcEventAction * theEventAction(
359
360 CexmcChargeExchangeReconstructor * recons
361 theEve
362
363 std::vector< G4double > calorimeterRegCut
364 if ( ProjectIsRead() )
365 {
366 calorimeterRegCuts = sObject.calorimet
367 }
368 else
369 {
370 G4Region * region( G4RegionStore::Get
371
372 if ( ! region )
373 throw CexmcException( CexmcCalorim
374
375 calorimeterRegCuts = region->GetProduc
376 }
377
378 CexmcNmbOfHitsInRanges nmbOfHitsSampled;
379 CexmcNmbOfHitsInRanges nmbOfHitsSampledFu
380 CexmcNmbOfHitsInRanges nmbOfHitsTriggered
381 CexmcNmbOfHitsInRanges nmbOfHitsTriggered
382 CexmcNmbOfHitsInRanges nmbOfOrphanHits;
383 G4int nmbOfFalseHitsTrig
384 G4int nmbOfFalseHitsTrig
385 G4int nmbOfSavedEvents(
386 G4int nmbOfSavedFastEven
387 CexmcRun * run( static_cast<
388
389 if ( run )
390 {
391 nmbOfHitsSampled = run->GetNmbOfHitsSa
392 nmbOfHitsSampledFull = run->GetNmbOfHi
393 nmbOfHitsTriggeredRealRange = run->Get
394 nmbOfHitsTriggeredRecRange = run->GetN
395 nmbOfOrphanHits = run->GetNmbOfOrphanH
396 nmbOfFalseHitsTriggeredEDT = run->GetN
397 nmbOfFalseHitsTriggeredRec = run->GetN
398 nmbOfSavedEvents = run->GetNmbOfSavedE
399 nmbOfSavedFastEvents = run->GetNmbOfSa
400 }
401
402 CexmcRunSObject sObjectToWrite = {
403 basePhysicsUsed, productionModelType,
404 physicsManager->GetProductionModel()->
405 physicsManager->GetProductionModel()->
406 calorimeterRegCuts, eventCountPolicy,
407 particleGun->GetParticleDefinition()->
408 particleGun->GetOrigPosition(), partic
409 particleGun->GetOrigMomentumAmp(),
410 primaryGeneratorAction->GetFwhmPosX(),
411 primaryGeneratorAction->GetFwhmPosY(),
412 primaryGeneratorAction->GetFwhmDirX(),
413 primaryGeneratorAction->GetFwhmDirY(),
414 primaryGeneratorAction->GetFwhmMomentu
415 edDigitizer->GetMonitorThreshold(),
416 edDigitizer->GetVetoCounterLeftThresho
417 edDigitizer->GetVetoCounterRightThresh
418 edDigitizer->GetCalorimeterLeftThresho
419 edDigitizer->GetCalorimeterRightThresh
420 edDigitizer->GetCalorimeterTriggerAlgo
421 edDigitizer->GetOuterCrystalsVetoAlgor
422 edDigitizer->GetOuterCrystalsVetoFract
423 edDigitizer->IsFiniteCrystalResolution
424 edDigitizer->GetCrystalResolutionData(
425 reconstructor->GetCalorimeterEntryPoin
426 reconstructor->GetCalorimeterEntryPoin
427 reconstructor->GetCrystalSelectionAlgo
428 reconstructor->IsInnerRefCrystalUsed()
429 reconstructor->GetCalorimeterEntryPoin
430 reconstructor->IsTableMassUsed(), reco
431 reconstructor->GetMassCutOPCenter(),
432 reconstructor->GetMassCutNOPCenter(),
433 reconstructor->GetMassCutOPWidth(), re
434 reconstructor->GetMassCutEllipseAngle(
435 reconstructor->IsAbsorbedEnergyCutUsed
436 reconstructor->GetAbsorbedEnergyCutCLC
437 reconstructor->GetAbsorbedEnergyCutCRC
438 reconstructor->GetAbsorbedEnergyCutCLW
439 reconstructor->GetAbsorbedEnergyCutCRW
440 reconstructor->GetAbsorbedEnergyCutEll
441 nmbOfHitsSampled, nmbOfHitsSampledFull
442 nmbOfHitsTriggeredRecRange, nmbOfOrpha
443 nmbOfFalseHitsTriggeredRec, nmbOfSaved
444 numberOfEventsProcessed, numberOfEvent
445 numberOfEventToBeProcessed, rProject,
446 cfFileName, evDataVerboseLevel, physic
447 reconstructor->GetExpectedMomentumAmp(
448 reconstructor->GetEDCollectionAlgorith
449
450 std::ofstream runDataFile( ( projectsDir
451 c_str(
452
453 {
454 boost::archive::binary_oarchive archi
455 archive << sObjectToWrite;
456 }
457 }
458
459 #endif
460
461
462 void CexmcRunManager::DoCommonEventLoop( G4in
463 G4in
464 {
465 G4int iEvent( 0 );
466 G4int iEventEffective( 0 );
467
468 for ( iEvent = 0; iEventEffective < nEvent
469 {
470 currentEvent = GenerateEvent( iEvent )
471 eventManager->ProcessOneEvent( current
472 CexmcEventInfo * eventInfo( static_ca
473 curren
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()->Apply
495 StackPreviousEvent( currentEvent );
496 currentEvent = 0;
497 if ( runAborted )
498 break;
499 }
500
501 numberOfEventsProcessed = iEvent;
502 numberOfEventsProcessedEffective = iEventE
503 }
504
505
506 #ifdef CEXMC_USE_PERSISTENCY
507
508 void CexmcRunManager::DoReadEventLoop( G4int
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( CexmcWeirdExcept
519
520 CexmcProductionModel * productionModel(
521 physic
522 if ( ! productionModel )
523 throw CexmcException( CexmcWeirdExcept
524
525 CexmcSetup * setup( static_cast< CexmcSet
526 if ( ! setup )
527 throw CexmcException( CexmcWeirdExcept
528
529 CexmcEventSObject evSObject;
530 CexmcEventFastSObject evFastSObject;
531
532 /* read events data */
533 std::ifstream eventsDataFile(
534 ( projectsDir + "/" +
535 if ( ! eventsDataFile )
536 throw CexmcException( CexmcReadProject
537
538 boost::archive::binary_iarchive evArchive
539
540 std::ifstream eventsFastDataFile(
541 ( projectsDir + "/" +
542 if ( ! eventsFastDataFile )
543 throw CexmcException( CexmcReadProject
544
545 boost::archive::binary_iarchive evFastArc
546
547 G4Event event;
548 currentEvent = &event;
549 G4SDManager * sdManager( G4SDManager:
550 event.SetHCofThisEvent( sdManager->Prepare
551 G4HCofThisEvent * hcOfThisEvent( event.Ge
552
553 G4DigiManager * digiManager( G4DigiManage
554
555 G4int hcId( digiManager->GetHitsCollectio
556 CexmcDetectorRoleName[
557 "/" + CexmcDetectorTyp
558 CexmcEnergyDepositCollection * monitorED(
559 ne
560 hcOfThisEvent->AddHitsCollection( hcId, mo
561 hcId = digiManager->GetHitsCollectionID(
562 CexmcDetectorRoleName[
563 "/" + CexmcDetectorTyp
564 CexmcEnergyDepositCollection * vetoCounte
565 ne
566 hcOfThisEvent->AddHitsCollection( hcId, ve
567 hcId = digiManager->GetHitsCollectionID(
568 CexmcDetectorRoleName[
569 "/" + CexmcDetectorTyp
570 CexmcEnergyDepositCollection * calorimete
571 ne
572 hcOfThisEvent->AddHitsCollection( hcId, ca
573 hcId = digiManager->GetHitsCollectionID(
574 CexmcDetectorRoleName[
575 "/" + CexmcDetectorTyp
576 CexmcTrackPointsCollection * monitorTP( n
577 hcOfThisEvent->AddHitsCollection( hcId, mo
578 hcId = digiManager->GetHitsCollectionID(
579 CexmcDetectorRoleName[
580 "/" + CexmcDetectorTyp
581 CexmcTrackPointsCollection * vetoCounterT
582 ne
583 hcOfThisEvent->AddHitsCollection( hcId, ve
584 hcId = digiManager->GetHitsCollectionID(
585 CexmcDetectorRoleName[
586 "/" + CexmcDetectorTyp
587 CexmcTrackPointsCollection * calorimeterT
588 ne
589 hcOfThisEvent->AddHitsCollection( hcId, ca
590 hcId = digiManager->GetHitsCollectionID(
591 CexmcDetectorRoleName[
592 "/" + CexmcDetectorTyp
593 CexmcTrackPointsCollection * targetTP( ne
594 hcOfThisEvent->AddHitsCollection( hcId, ta
595
596 #ifdef CEXMC_USE_CUSTOM_FILTER
597 if ( customFilter )
598 customFilter->SetAddressedData( &evFas
599 #endif
600
601 G4int nmbOfSavedEvents( rEvDataVerboseLe
602 sObject.nmbOfSave
603 G4bool eventDataWrittenOnEveryTPT( rEvDat
604 CexmcW
605
606 for ( int i( 0 ); i < nmbOfSavedEvents; +
607 {
608 evFastArchive >> evFastSObject;
609
610 if ( nEventCount < curEventRead )
611 {
612 if ( eventDataWrittenOnEveryTPT ||
613 evFastSObject.edDigitizerHasT
614 {
615 evArchive >> evSObject;
616 if ( evFastSObject.edDigitizer
617 ++nEventCount;
618 }
619 continue;
620 }
621
622 ++iEvent;
623
624 productionModel->SetTriggeredAngularRa
625
626 const CexmcAngularRangeList & trigger
627 productionMode
628
629 if ( ! eventDataWrittenOnEveryTPT &&
630 ! evFastSObject.edDigitizerHasTri
631 {
632 #ifdef CEXMC_USE_CUSTOM_FILTER
633 /* user must be aware that using t
634 * scripts for poor event data set
635 * errors! This is because most of
636 * events with EDT trigger. There
637 * was written on every TPT event.
638 if ( customFilter && ! customFilte
639 continue;
640 #endif
641 SaveCurrentTPTEvent( evFastSObject
642 ProjectIsSaved() && ! ski
643 continue;
644 }
645
646 evArchive >> evSObject;
647
648 G4bool skipEDTOnThisEvent( false );
649
650 #ifdef CEXMC_USE_CUSTOM_FILTER
651 if ( customFilter && ! customFilter->E
652 continue;
653 if ( customFilter && ! customFilter->E
654 {
655 if ( ! eventDataWrittenOnEveryTPT
656 {
657 SaveCurrentTPTEvent( evFastSOb
658 ProjectIs
659 continue;
660 }
661 skipEDTOnThisEvent = true;
662 }
663 #endif
664
665 event.SetEventID( evSObject.eventId );
666
667 /* AAA: beware! If something throws an
668 * then there would be segmentation vi
669 * would try to delete local variable
670 * See also BBB */
671
672 monitorED->GetMap()->operator[]( 0 ) =
673 vetoCounterED->GetMap()->operator[]( 0
674 vetoCounterED->GetMap()->operator[]( 1
675 CexmcEnergyDepositInLeftRightS
676
677 G4int row( 0 );
678 G4int column( 0 );
679 for ( CexmcEnergyDepositCalorimeterCol
680 k( evSObject.calorimeterEDLeft
681 k != evSObject.calorimeter
682 {
683 G4int index( row <<
684 CexmcEnergyDepositInCalorimete
685 column = 0;
686 for ( CexmcEnergyDepositCrystalRow
687 l( k->begin() ); l != k->e
688 {
689 calorimeterED->GetMap()->opera
690 ++column;
691 }
692 ++row;
693 }
694 row = 0;
695 for ( CexmcEnergyDepositCalorimeterCol
696 k( evSObject.calorimeterEDRigh
697 k != evSObject.calorimeter
698 {
699 G4int index(
700 1 << CexmcEnergyDepositInLeftR
701 | row <<
702 CexmcEnergyDepositInCalorimete
703 column = 0;
704 for ( CexmcEnergyDepositCrystalRow
705 l( k->begin() ); l != k->e
706 {
707 calorimeterED->GetMap()->opera
708 ++column;
709 }
710 ++row;
711 }
712
713 CexmcTrackPointInfo monitorTPInfo( ev
714 CexmcTrackPointInfo targetTPBeamParti
715 evSObj
716 CexmcTrackPointInfo targetTPOutputPar
717 evSObj
718 CexmcTrackPointInfo targetTPNucleusPa
719 evSObj
720 CexmcTrackPointInfo targetTPOutputPar
721 evSObject.targetTPOutp
722 CexmcTrackPointInfo targetTPOutputPar
723 evSObject.targetTPOutp
724 CexmcTrackPointInfo vetoCounterTPLeft
725
726 CexmcTrackPointInfo vetoCounterTPRigh
727
728 CexmcTrackPointInfo calorimeterTPLeft
729
730 CexmcTrackPointInfo calorimeterTPRigh
731
732
733 if ( monitorTPInfo.IsValid() )
734 monitorTP->GetMap()->operator[]( m
735
736 if ( targetTPBeamParticleInfo.IsValid(
737 targetTP->GetMap()->operator[](
738 targetTPBeamParticleInfo.t
739
740 if ( targetTPOutputParticleInfo.IsVali
741 targetTP->GetMap()->operator[](
742 targetTPOutputParticleInfo
743
744 if ( targetTPNucleusParticleInfo.IsVal
745 targetTP->GetMap()->operator[](
746 targetTPNucleusParticleInf
747
748 if ( targetTPOutputParticleDecayProduc
749 targetTP->GetMap()->operator[](
750 targetTPOutputParticleDeca
751 &targetTPOutputPar
752 if ( targetTPOutputParticleDecayProduc
753 targetTP->GetMap()->operator[](
754 targetTPOutputParticleDeca
755 &targetTPOutputPar
756 if ( vetoCounterTPLeftInfo.IsValid() )
757 vetoCounterTP->GetMap()->operator[
758 vetoCounterTPLeftInfo.trac
759 if ( vetoCounterTPRightInfo.IsValid()
760 vetoCounterTP->GetMap()->operator[
761 1 << CexmcTrackPointsInLeftRig
762 vetoCounterTPRightInfo.tra
763
764 G4ThreeVector pos;
765 if ( calorimeterTPLeftInfo.IsValid() )
766 {
767 pos = calorimeterTPLeftInfo.positi
768 setup->ConvertToCrystalGeometry(
769 calorimeterTPLeftInfo.posi
770 calorimeterTPLeftInfo.positionLoca
771 calorimeterTP->GetMap()->operator[
772 row << CexmcTrackPointsInCalor
773
774 column << CexmcTrackPointsInCa
775
776 calorimeterTPLeftInfo.trackId
777 }
778 if ( calorimeterTPRightInfo.IsValid()
779 {
780 pos = calorimeterTPRightInfo.posit
781 setup->ConvertToCrystalGeometry(
782 calorimeterTPRightInfo.pos
783 calorimeterTPRightInfo.positionLoc
784 calorimeterTP->GetMap()->operator[
785 1 << CexmcTrackPointsInLeftRig
786 row << CexmcTrackPointsInCalor
787
788 column << CexmcTrackPointsInCa
789
790 calorimeterTPRightInfo.trackId
791 }
792
793 productionModel->SetProductionModelDat
794 ev
795
796 const CexmcEventAction * eventAction(
797 static_cast< const CexmcEventA
798 if ( ! eventAction )
799 {
800 /* BBB: all hits collections must
801 * anything from here, otherwise ~
802 * local variable evSObject's fiel
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( CexmcEventAc
811 }
812
813 if ( skipEDTOnThisEvent )
814 event.SetUserInformation( new Cexm
815
816
817 CexmcEventAction * theEventAction( co
818
819 theEventAction->EndOfEventAction( &eve
820
821 CexmcEventInfo * eventInfo( static_ca
822
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 ==
841 break;
842 }
843
844 curEventRead = nEventCount + iEventEffecti
845
846 numberOfEventsProcessed = iEvent;
847 numberOfEventsProcessedEffective = iEventE
848
849 #ifdef CEXMC_USE_CUSTOM_FILTER
850 if ( customFilter )
851 customFilter->SetAddressedData( NULL,
852 #endif
853 }
854
855
856 void CexmcRunManager::SaveCurrentTPTEvent(
857 const CexmcEve
858 const CexmcAng
859 G4bool writeT
860 {
861 CexmcRun * run( static_cast< CexmcRun * >
862
863 if ( ! run )
864 return;
865
866 for ( CexmcAngularRangeList::const_iterato
867
868 {
869 run->IncrementNmbOfHitsSampledFull( k-
870 if ( evFastSObject.edDigitizerMonitorH
871 run->IncrementNmbOfHitsSampled( k-
872 }
873
874 if ( writeToDatabase )
875 {
876 fastEventsArchive->operator<<( evFastS
877 run->IncrementNmbOfSavedFastEvents();
878 }
879 }
880
881 #endif
882
883
884 /* mostly adopted from G4RunManager::DoEventLo
885 void CexmcRunManager::DoEventLoop( G4int nEv
886 G4int nSe
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 + "/" +
916 boost::archive::binary_oarchive e
917 std::ofstream fastEventsDataFile
918 ( projectsDir + "/" +
919 boost::archive::binary_oarchive f
920
921 eventsArchive = &eventsArchive_;
922 fastEventsArchive = &fastEventsArc
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 + "/" +
936 boost::archive::binary_oarchive e
937 std::ofstream fastEventsDataFile
938 ( projectsDir + "/" +
939 boost::archive::binary_oarchive f
940
941 eventsArchive = &eventsArchive_;
942 fastEventsArchive = &fastEventsArc
943 DoCommonEventLoop( nEvent, cmd, nS
944 }
945 else
946 {
947 DoCommonEventLoop( nEvent, cmd, nS
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 " <
964 " events processed." <<
965 }
966 else
967 {
968 G4cout << " Number of events proc
969 numberOfEventsProcessed
970 numberOfEventsProcessedE
971 }
972 G4cout << " " << *timer << G4endl;
973 }
974 }
975
976
977 #ifdef CEXMC_USE_PERSISTENCY
978
979 void CexmcRunManager::PrintReadRunData( void
980 {
981 if ( ! ProjectIsRead() )
982 return;
983
984 G4bool refCrystalInfoPrinted( false );
985
986 G4cout << CEXMC_LINE_START << "Run data re
987 "'" << G4endl;
988 G4cout << " (archive class v
989 sObject.actualVersion << ")" <<
990 if ( ! sObject.rProject.empty() )
991 {
992 G4cout << " -- Based on project '" <<
993 G4endl;
994 if ( ! sObject.cfFileName.empty() )
995 G4cout << " -- Custom filter scri
996 "' was used" << G4endl;
997 }
998 G4cout << " -- Event data verbose level (
999 "2 - interactions): " << sObject
1000 if ( ! sObject.rProject.empty() )
1001 {
1002 if ( sObject.evDataVerboseLevel == Ce
1003 {
1004 G4cout << " -- (fdb file contain
1005 ( sObject.interactionsW
1006 "only interactions wh
1007 "all interactions" )
1008 }
1009 }
1010 G4cout << " -- Base physics used"
1011 "(1 - QGSP_BERT, 2 - QGSP_BIC_E
1012 sObject.basePhysicsUsed << G4en
1013 G4cout << " -- Production model (1 - pi0
1014 sObject.productionModelType <<
1015 G4cout << " -- Geometry definition file:
1016 G4endl;
1017 G4cout << " -- Angular ranges: " << sObj
1018 G4cout << " -- Eta decay modes: " << G4e
1019 G4Eta::Definition()->GetDecayTable()->Dum
1020 G4cout << " -- Fermi motion status (0 -
1021 sObject.fermiMotionIsOn << G4en
1022 if ( sObject.calorimeterRegCuts.size() <
1023 throw CexmcException( CexmcWeirdExcep
1024 G4cout << " -- Production cuts in calori
1025 G4BestUnit( sObject.calorimeter
1026 G4BestUnit( sObject.calorimeter
1027 G4BestUnit( sObject.calorimeter
1028 G4BestUnit( sObject.calorimeter
1029 G4cout << " -- Proposed max interaction
1030 G4BestUnit( sObject.proposedMax
1031 G4cout << " -- Event count policy (0 - a
1032 ": " << sObject.eventCountPolic
1033 G4cout << " -- Number of events (process
1034 sObject.numberOfEventsProcessed
1035 sObject.numberOfEventsProcessed
1036 sObject.numberOfEventsToBeProce
1037 G4cout << " -- Incident beam particle: "
1038 G4cout << " position: "
1039 G4BestUnit( sObject.beamPos, "L
1040 G4cout << " direction: "
1041 G4ThreeVector( sObject.beamDir
1042 G4cout << " momentum: "
1043 G4BestUnit( sObject.beamMomentu
1044 G4cout << " momentum fwhm: "
1045 G4endl;
1046 G4cout << " pos fwhm (x): "
1047 G4BestUnit( sObject.beamFwhmPos
1048 G4cout << " pos fwhm (y): "
1049 G4BestUnit( sObject.beamFwhmPos
1050 G4cout << " dir fwhm (x): "
1051 " deg" << G4endl;
1052 G4cout << " dir fwhm (y): "
1053 " deg" << G4endl;
1054 G4cout << " -- Monitor ED threshold: " <
1055 G4BestUnit( sObject.monitorEDTh
1056 G4cout << " -- Veto counter (l/r) ED thr
1057 G4BestUnit( sObject.vetoCounter
1058 " / " <<
1059 G4BestUnit( sObject.vetoCounter
1060 G4endl;
1061 G4cout << " -- Calorimeter (l/r) ED thre
1062 G4BestUnit( sObject.calorimeter
1063 " / " <<
1064 G4BestUnit( sObject.calorimeter
1065 G4endl;
1066 G4cout << " -- Calorimeter trigger algor
1067 sObject.calorimeterTriggerAlgor
1068 G4cout << " -- Outer crystals veto algor
1069 "(0 - none, 1 - max, 2 - fracti
1070 sObject.outerCrystalsVetoAlgori
1071 if ( sObject.outerCrystalsVetoAlgorithm =
1072 CexmcFractionOfEDInOuterCrystalsVeto
1073 {
1074 G4cout << " -- Outer crystals veto f
1075 sObject.outerCrystalsVetoFr
1076 }
1077 G4cout << " -- Finite crystal resolution
1078 sObject.applyFiniteCrystalResol
1079 if ( sObject.applyFiniteCrystalResolution
1080 {
1081 G4cout << " -- Crystal resolution da
1082 sObject.crystalResolutionDa
1083 }
1084 G4cout << " -- Reconstructor settings: "
1085 if ( sObject.expectedMomentumAmp > 0 )
1086 {
1087 G4cout << " -- expected momentum
1088 G4BestUnit( sObject.expecte
1089 }
1090 G4cout << " -- ed collection algorith
1091 sObject.edCollectionAlgorithm <
1092 if ( sObject.edCollectionAlgorithm == Cex
1093 {
1094 G4cout <<
1095 " -- inner crystal used as re
1096 sObject.useInnerRefCrystal << G4e
1097 refCrystalInfoPrinted = true;
1098 }
1099 G4cout << " -- entry point definition
1100 G4cout << " (0 - center of calorim
1101 "max ED," << G4endl;
1102 G4cout << " 2 - linear, 3 - squar
1103 sObject.epDefinitionAlgorithm <
1104 G4cout << " -- entry point depth defi
1105 "(0 - plain, 1 - sphere
1106 sObject.epDepthDefi
1107 G4cout << " -- entry point depth: " <
1108 G4BestUnit( sObject.epDepth, "L
1109 if ( sObject.epDefinitionAlgorithm == Cex
1110 sObject.epDefinitionAlgorithm == Cex
1111 {
1112 G4cout <<
1113 " -- crystal selection algori
1114 sObject.csAlgorithm << G4endl;
1115 }
1116 if ( ! refCrystalInfoPrinted &&
1117 ( sObject.epDefinitionAlgorithm ==
1118 CexmcEntryPoi
1119 ( ( sObject.epDefinitionAlgorithm ==
1120 sObject.epDefinitionAlgorithm ==
1121 Cexmc
1122 sObject.csAlgorithm == CexmcSe
1123 {
1124 G4cout <<
1125 " -- inner crystal used as re
1126 sObject.useInnerRefCrystal << G4e
1127 }
1128 G4cout << " -- table mass of output p
1129 "(0 - no, 1 - yes): " <
1130 G4cout << " -- mass cut is enabled (0
1131 sObject.useMassCut << G4endl;
1132 if ( sObject.useMassCut )
1133 {
1134 G4cout << " -- mass cut output pa
1135 G4BestUnit( sObject.mCutOPC
1136 G4cout << " -- mass cut nucleus o
1137 G4BestUnit( sObject.mCutNOP
1138 G4cout << " -- mass cut output pa
1139 G4BestUnit( sObject.mCutOPW
1140 G4cout << " -- mass cut nucleus o
1141 "ellipse: "
1142 << G4BestUnit( sObject.mCutNOP
1143 G4cout << " -- mass cut angle of
1144 sObject.mCutAngle / deg <<
1145 }
1146 G4cout << " -- absorbed energy cut is
1147 sObject.useAbsorbedEnergyCut <<
1148 if ( sObject.useAbsorbedEnergyCut )
1149 {
1150 G4cout << " -- absorbed energy cu
1151 G4BestUnit( sObject.aeCutCL
1152 G4cout << " -- absorbed energy cu
1153 G4BestUnit( sObject.aeCutCR
1154 G4cout << " -- absorbed energy cu
1155 "ellipse: " <<
1156 G4BestUnit( sObject.aeCutCL
1157 G4cout << " -- absorbed energy cu
1158 "ellipse: "
1159 << G4BestUnit( sObject.aeCutCR
1160 G4cout << " -- absorbed energy cu
1161 sObject.aeCutAngle / deg <<
1162 }
1163 G4cout << G4endl;
1164 CexmcRunAction::PrintResults( sObject.nmb
1165 sObject.nmb
1166 sObject.nmb
1167 sObject.nmb
1168 sObject.nmb
1169 sObject.ang
1170 sObject.nmb
1171 sObject.nmb
1172 G4cout << G4endl;
1173 }
1174
1175
1176 void CexmcRunManager::ReadAndPrintEventsData
1177 {
1178 if ( ! ProjectIsRead() )
1179 return;
1180
1181 CexmcEventSObject evSObject;
1182
1183 /* read events data */
1184 std::ifstream eventsDataFile(
1185 ( projectsDir + "/" +
1186 if ( ! eventsDataFile )
1187 throw CexmcException( CexmcReadProjec
1188
1189 boost::archive::binary_iarchive evArchiv
1190
1191 for ( int i( 0 ); i < sObject.nmbOfSaved
1192 {
1193 evArchive >> evSObject;
1194
1195 if ( ! evSObject.edDigitizerMonitorHa
1196 continue;
1197
1198 CexmcEnergyDepositStore edStore( evS
1199 evSObject.vetoCounterEDLeft, evSO
1200 evSObject.calorimeterEDLeft, evSO
1201 0, 0, 0, 0, evSObject.calorimeter
1202 evSObject.calorimeterEDRightColle
1203
1204 CexmcTrackPointsStore tpStore( evS
1205 evSObject.targetTPBeamParticle, e
1206 evSObject.targetTPNucleusParticle
1207 evSObject.targetTPOutputParticleD
1208 evSObject.targetTPOutputParticleD
1209 evSObject.vetoCounterTPLeft, evSO
1210 evSObject.calorimeterTPLeft, evSO
1211
1212 const CexmcProductionModelData & pmD
1213 e
1214
1215 G4cout << "Event " << evSObject.event
1216 CexmcEventAction::PrintTrackPoints( &
1217 G4cout << " --- Production model data
1218 CexmcEventAction::PrintEnergyDeposit(
1219 }
1220 }
1221
1222
1223 void CexmcRunManager::PrintReadData(
1224 const CexmcOutput
1225 {
1226 if ( ! ProjectIsRead() )
1227 return;
1228
1229 G4bool addSpace( false );
1230
1231 CexmcOutputDataTypeSet::const_iterator f
1232 outputDat
1233 if ( found != outputData.end() )
1234 {
1235 G4String cmd( G4String( "cat " ) + p
1236 gdmlFileExtension );
1237 if ( system( cmd ) != 0 )
1238 throw CexmcException( CexmcReadPr
1239
1240 if ( zipGdmlFile )
1241 {
1242 cmd = G4String( "bzip2 " ) + proj
1243
1244 if ( system( cmd ) != 0 )
1245 throw CexmcException( CexmcFi
1246 }
1247
1248 addSpace = true;
1249 }
1250
1251 found = outputData.find( CexmcOutputEvent
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:
1269 for ( CexmcDecayBranchesStore::const_
1270 k( sObject.etaDecayTable.GetD
1271 k != sObject.etaDecayTable.Ge
1272 {
1273 etaDecayTable->GetDecayChannel( k
1274 }
1275
1276 PrintReadRunData();
1277 }
1278 }
1279
1280
1281 #ifdef CEXMC_USE_CUSTOM_FILTER
1282
1283 void CexmcRunManager::SetCustomFilter( const
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( CexmcCmdIsNotAl
1297
1298 cfFileName = cfFileName_;
1299
1300 customFilter = new CexmcCustomFilterEval(
1301 }
1302
1303 #endif
1304
1305 #endif
1306
1307
1308 void CexmcRunManager::RegisterScenePrimitive
1309 {
1310 G4VisManager * visManager( static_cast<
1311 G4VVisMan
1312 if ( ! visManager )
1313 return;
1314
1315 G4Scene * curScene( visManager->Get
1316 if ( ! curScene )
1317 return;
1318
1319 /* G4Scene declarations lack this kind of
1320 #if G4VERSION_NUMBER < 960
1321 typedef std::vector< G4VModel * > ML
1322 #else
1323 typedef std::vector< G4Scene::Model > ML
1324 #endif
1325 const MList & mList( curScene->GetRunDur
1326
1327 for ( MList::const_iterator k( mList.beg
1328 {
1329 #if G4VERSION_NUMBER < 960
1330 const G4String & modelDesc( ( *k )->
1331 #else
1332 const G4String & modelDesc( k->fpMod
1333 #endif
1334 if ( modelDesc == CexmcScenePrimitive
1335 return;
1336 }
1337
1338 CexmcSetup * setup( static_cast< CexmcSe
1339 if ( ! setup )
1340 throw CexmcException( CexmcWeirdExcep
1341
1342 /* BEWARE: looks like G4Scene won't delet
1343 * termination! Hence destructor of the n
1344 curScene->AddRunDurationModel( new CexmcS
1345 }
1346
1347
1348 void CexmcRunManager::BeamParticleChangeHook
1349 {
1350 const CexmcEventAction * eventAction(
1351 static_cast< const CexmcEvent
1352 if ( ! eventAction )
1353 throw CexmcException( CexmcWeirdExcep
1354
1355 CexmcEventAction * theEventAction(
1356
1357 theEventAction->BeamParticleChangeHook();
1358 }
1359
1360
1361 void CexmcRunManager::SetupConstructionHook(
1362 {
1363 #ifdef CEXMC_USE_PERSISTENCY
1364 /* save gdml file */
1365 G4String cmd( "" );
1366 CexmcExceptionType exceptionType( CexmcS
1367
1368 if ( zipGdmlFile )
1369 {
1370 if ( ProjectIsRead() )
1371 {
1372 cmd = G4String( "bzip2 " ) + proj
1373
1374 }
1375 else
1376 {
1377 if ( ProjectIsSaved() )
1378 cmd = G4String( "bzip2 -c " )
1379 projectsDir + "/" + p
1380 }
1381 exceptionType = CexmcFileCompressExce
1382 }
1383 else
1384 {
1385 if ( ! ProjectIsRead() && ProjectIsSa
1386 cmd = G4String( "cp " ) + gdmlFil
1387
1388 /* else already saved in ReadPreinitP
1389 }
1390
1391 if ( ! cmd.empty() && system( cmd ) != 0
1392 throw CexmcException( exceptionType )
1393 #endif
1394 }
1395
1396