ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CexmcRunManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CexmcRunManager.cc
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( true ),
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