Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CexmcRunManager.cc
Go to the documentation of this file.
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 ), areLiveHistogramsEnabled( false ),
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( 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 (1 - QGSP_BERT, 2 - QGSP_BIC_EMY): " <<
1011  sObject.basePhysicsUsed << G4endl;
1012  G4cout << " -- Production model (1 - pi0, 2 - eta): " <<
1013  sObject.productionModelType << G4endl;
1014  G4cout << " -- Geometry definition file: " << sObject.gdmlFileName <<
1015  G4endl;
1016  G4cout << " -- Angular ranges: " << sObject.angularRanges << G4endl;
1017  G4cout << " -- Eta decay modes: " << G4endl;
1018  G4Eta::Definition()->GetDecayTable()->DumpInfo();
1019  G4cout << " -- Fermi motion status (0 - disabled, 1 - enabled): " <<
1020  sObject.fermiMotionIsOn << G4endl;
1021  if ( sObject.calorimeterRegCuts.size() < 4 )
1022  throw CexmcException( CexmcWeirdException );
1023  G4cout << " -- Production cuts in calorimeter (gamma, e-, e+, p): " <<
1024  G4BestUnit( sObject.calorimeterRegCuts[ 0 ], "Length" ) << ", " <<
1025  G4BestUnit( sObject.calorimeterRegCuts[ 1 ], "Length" ) << ", " <<
1026  G4BestUnit( sObject.calorimeterRegCuts[ 2 ], "Length" ) << ", " <<
1027  G4BestUnit( sObject.calorimeterRegCuts[ 3 ], "Length" ) << G4endl;
1028  G4cout << " -- Proposed max interaction length in the target: " <<
1029  G4BestUnit( sObject.proposedMaxIL, "Length" ) << G4endl;
1030  G4cout << " -- Event count policy (0 - all, 1 - interaction, 2 - trigger)"
1031  ": " << sObject.eventCountPolicy << G4endl;
1032  G4cout << " -- Number of events (processed / effective / ordered): " <<
1033  sObject.numberOfEventsProcessed << " / " <<
1034  sObject.numberOfEventsProcessedEffective << " / " <<
1035  sObject.numberOfEventsToBeProcessed << G4endl;
1036  G4cout << " -- Incident beam particle: " << sObject.beamParticle << G4endl;
1037  G4cout << " position: " <<
1038  G4BestUnit( sObject.beamPos, "Length" ) << G4endl;
1039  G4cout << " direction: " <<
1040  G4ThreeVector( sObject.beamDir ) << G4endl;
1041  G4cout << " momentum: " <<
1042  G4BestUnit( sObject.beamMomentumAmp, "Energy" ) << G4endl;
1043  G4cout << " momentum fwhm: " << sObject.beamFwhmMomentumAmp <<
1044  G4endl;
1045  G4cout << " pos fwhm (x): " <<
1046  G4BestUnit( sObject.beamFwhmPosX, "Length" ) << G4endl;
1047  G4cout << " pos fwhm (y): " <<
1048  G4BestUnit( sObject.beamFwhmPosY, "Length" ) << G4endl;
1049  G4cout << " dir fwhm (x): " << sObject.beamFwhmDirX / deg <<
1050  " deg" << G4endl;
1051  G4cout << " dir fwhm (y): " << sObject.beamFwhmDirY / deg <<
1052  " deg" << G4endl;
1053  G4cout << " -- Monitor ED threshold: " <<
1054  G4BestUnit( sObject.monitorEDThreshold, "Energy" ) << G4endl;
1055  G4cout << " -- Veto counter (l/r) ED threshold: " <<
1056  G4BestUnit( sObject.vetoCounterEDLeftThreshold, "Energy" ) <<
1057  " / " <<
1058  G4BestUnit( sObject.vetoCounterEDRightThreshold, "Energy" ) <<
1059  G4endl;
1060  G4cout << " -- Calorimeter (l/r) ED threshold: " <<
1061  G4BestUnit( sObject.calorimeterEDLeftThreshold, "Energy" ) <<
1062  " / " <<
1063  G4BestUnit( sObject.calorimeterEDRightThreshold, "Energy" ) <<
1064  G4endl;
1065  G4cout << " -- Calorimeter trigger algorithm (0 - all, 1 - inner): " <<
1066  sObject.calorimeterTriggerAlgorithm << G4endl;
1067  G4cout << " -- Outer crystals veto algorithm "
1068  "(0 - none, 1 - max, 2 - fraction): " <<
1069  sObject.outerCrystalsVetoAlgorithm << G4endl;
1070  if ( sObject.outerCrystalsVetoAlgorithm ==
1071  CexmcFractionOfEDInOuterCrystalsVeto )
1072  {
1073  G4cout << " -- Outer crystals veto fraction: " <<
1074  sObject.outerCrystalsVetoFraction << G4endl;
1075  }
1076  G4cout << " -- Finite crystal resolution applied (0 - no, 1 - yes): " <<
1077  sObject.applyFiniteCrystalResolution << G4endl;
1078  if ( sObject.applyFiniteCrystalResolution )
1079  {
1080  G4cout << " -- Crystal resolution data: " <<
1081  sObject.crystalResolutionData;
1082  }
1083  G4cout << " -- Reconstructor settings: " << G4endl;
1084  if ( sObject.expectedMomentumAmp > 0 )
1085  {
1086  G4cout << " -- expected momentum in the target: " <<
1087  G4BestUnit( sObject.expectedMomentumAmp, "Energy" ) << G4endl;
1088  }
1089  G4cout << " -- ed collection algorithm (0 - all, 1 - adjacent): " <<
1090  sObject.edCollectionAlgorithm << G4endl;
1091  if ( sObject.edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
1092  {
1093  G4cout <<
1094  " -- inner crystal used as reference (0 - no, 1 - yes): " <<
1095  sObject.useInnerRefCrystal << G4endl;
1096  refCrystalInfoPrinted = true;
1097  }
1098  G4cout << " -- entry point definition algorithm " << G4endl;
1099  G4cout << " (0 - center of calorimeter, 1 - center of crystal with "
1100  "max ED," << G4endl;
1101  G4cout << " 2 - linear, 3 - square): " <<
1102  sObject.epDefinitionAlgorithm << G4endl;
1103  G4cout << " -- entry point depth definition algorithm "
1104  "(0 - plain, 1 - sphere): " <<
1105  sObject.epDepthDefinitionAlgorithm << G4endl;
1106  G4cout << " -- entry point depth: " <<
1107  G4BestUnit( sObject.epDepth, "Length" ) << G4endl;
1108  if ( sObject.epDefinitionAlgorithm == CexmcEntryPointByLinearEDWeights ||
1109  sObject.epDefinitionAlgorithm == CexmcEntryPointBySqrtEDWeights )
1110  {
1111  G4cout <<
1112  " -- crystal selection algorithm (0 - all, 1 - adjacent): " <<
1113  sObject.csAlgorithm << G4endl;
1114  }
1115  if ( ! refCrystalInfoPrinted &&
1116  ( sObject.epDefinitionAlgorithm ==
1117  CexmcEntryPointInTheCenterOfCrystalWithMaxED ||
1118  ( ( sObject.epDefinitionAlgorithm == CexmcEntryPointBySqrtEDWeights ||
1119  sObject.epDefinitionAlgorithm ==
1120  CexmcEntryPointByLinearEDWeights ) &&
1121  sObject.csAlgorithm == CexmcSelectAdjacentCrystals ) ) )
1122  {
1123  G4cout <<
1124  " -- inner crystal used as reference (0 - no, 1 - yes): " <<
1125  sObject.useInnerRefCrystal << G4endl;
1126  }
1127  G4cout << " -- table mass of output particle used "
1128  "(0 - no, 1 - yes): " << sObject.useTableMass << G4endl;
1129  G4cout << " -- mass cut is enabled (0 - no, 1 - yes): " <<
1130  sObject.useMassCut << G4endl;
1131  if ( sObject.useMassCut )
1132  {
1133  G4cout << " -- mass cut output particle center: " <<
1134  G4BestUnit( sObject.mCutOPCenter, "Energy" ) << G4endl;
1135  G4cout << " -- mass cut nucleus output particle center: " <<
1136  G4BestUnit( sObject.mCutNOPCenter, "Energy" ) << G4endl;
1137  G4cout << " -- mass cut output particle width of the ellipse: " <<
1138  G4BestUnit( sObject.mCutOPWidth, "Energy" ) << G4endl;
1139  G4cout << " -- mass cut nucleus output particle width of the "
1140  "ellipse: "
1141  << G4BestUnit( sObject.mCutNOPWidth, "Energy" ) << G4endl;
1142  G4cout << " -- mass cut angle of the ellipse: " <<
1143  sObject.mCutAngle / deg << " deg" << G4endl;
1144  }
1145  G4cout << " -- absorbed energy cut is enabled (0 - no, 1 - yes): " <<
1146  sObject.useAbsorbedEnergyCut << G4endl;
1147  if ( sObject.useAbsorbedEnergyCut )
1148  {
1149  G4cout << " -- absorbed energy cut left calorimeter center: " <<
1150  G4BestUnit( sObject.aeCutCLCenter, "Energy" ) << G4endl;
1151  G4cout << " -- absorbed energy cut right calorimeter center: " <<
1152  G4BestUnit( sObject.aeCutCRCenter, "Energy" ) << G4endl;
1153  G4cout << " -- absorbed energy cut left calorimeter width of the "
1154  "ellipse: " <<
1155  G4BestUnit( sObject.aeCutCLWidth, "Energy" ) << G4endl;
1156  G4cout << " -- absorbed energy cut right calorimeter width of the "
1157  "ellipse: "
1158  << G4BestUnit( sObject.aeCutCRWidth, "Energy" ) << G4endl;
1159  G4cout << " -- absorbed energy cut angle of the ellipse: " <<
1160  sObject.aeCutAngle / deg << " deg" << G4endl;
1161  }
1162  G4cout << G4endl;
1163  CexmcRunAction::PrintResults( sObject.nmbOfHitsSampled,
1164  sObject.nmbOfHitsSampledFull,
1165  sObject.nmbOfHitsTriggeredRealRange,
1166  sObject.nmbOfHitsTriggeredRecRange,
1167  sObject.nmbOfOrphanHits,
1168  sObject.angularRanges,
1169  sObject.nmbOfFalseHitsTriggeredEDT,
1170  sObject.nmbOfFalseHitsTriggeredRec );
1171  G4cout << G4endl;
1172 }
1173 
1174 
1175 void CexmcRunManager::ReadAndPrintEventsData( void ) const
1176 {
1177  if ( ! ProjectIsRead() )
1178  return;
1179 
1180  CexmcEventSObject evSObject;
1181 
1182  /* read events data */
1183  std::ifstream eventsDataFile(
1184  ( projectsDir + "/" + rProject + ".edb" ).c_str() );
1185  if ( ! eventsDataFile )
1186  throw CexmcException( CexmcReadProjectIncomplete );
1187 
1188  boost::archive::binary_iarchive evArchive( eventsDataFile );
1189 
1190  for ( int i( 0 ); i < sObject.nmbOfSavedEvents; ++i )
1191  {
1192  evArchive >> evSObject;
1193 
1194  if ( ! evSObject.edDigitizerMonitorHasTriggered )
1195  continue;
1196 
1197  CexmcEnergyDepositStore edStore( evSObject.monitorED,
1198  evSObject.vetoCounterEDLeft, evSObject.vetoCounterEDRight,
1199  evSObject.calorimeterEDLeft, evSObject.calorimeterEDRight,
1200  0, 0, 0, 0, evSObject.calorimeterEDLeftCollection,
1201  evSObject.calorimeterEDRightCollection );
1202 
1203  CexmcTrackPointsStore tpStore( evSObject.monitorTP,
1204  evSObject.targetTPBeamParticle, evSObject.targetTPOutputParticle,
1205  evSObject.targetTPNucleusParticle,
1206  evSObject.targetTPOutputParticleDecayProductParticle1,
1207  evSObject.targetTPOutputParticleDecayProductParticle2,
1208  evSObject.vetoCounterTPLeft, evSObject.vetoCounterTPRight,
1209  evSObject.calorimeterTPLeft, evSObject.calorimeterTPRight );
1210 
1211  const CexmcProductionModelData & pmData(
1212  evSObject.productionModelData );
1213 
1214  G4cout << "Event " << evSObject.eventId << G4endl;
1215  CexmcEventAction::PrintTrackPoints( &tpStore );
1216  G4cout << " --- Production model data: " << pmData;
1217  CexmcEventAction::PrintEnergyDeposit( &edStore );
1218  }
1219 }
1220 
1221 
1222 void CexmcRunManager::PrintReadData(
1223  const CexmcOutputDataTypeSet & outputData ) const
1224 {
1225  if ( ! ProjectIsRead() )
1226  return;
1227 
1228  G4bool addSpace( false );
1229 
1230  CexmcOutputDataTypeSet::const_iterator found(
1231  outputData.find( CexmcOutputGeometry ) );
1232  if ( found != outputData.end() )
1233  {
1234  G4String cmd( G4String( "cat " ) + projectsDir + "/" + rProject +
1235  gdmlFileExtension );
1236  if ( system( cmd ) != 0 )
1237  throw CexmcException( CexmcReadProjectIncomplete );
1238 
1239  if ( zipGdmlFile )
1240  {
1241  cmd = G4String( "bzip2 " ) + projectsDir + "/" + rProject +
1242  gdmlFileExtension;
1243  if ( system( cmd ) != 0 )
1244  throw CexmcException( CexmcFileCompressException );
1245  }
1246 
1247  addSpace = true;
1248  }
1249 
1250  found = outputData.find( CexmcOutputEvents );
1251  if ( found != outputData.end() )
1252  {
1253  if ( addSpace )
1254  G4cout << G4endl << G4endl;
1255 
1256  ReadAndPrintEventsData();
1257 
1258  addSpace = true;
1259  }
1260 
1261  found = outputData.find( CexmcOutputRun );
1262  if ( found != outputData.end() )
1263  {
1264  if ( addSpace )
1265  G4cout << G4endl << G4endl;
1266 
1267  G4DecayTable * etaDecayTable( G4Eta::Definition()->GetDecayTable() );
1268  for ( CexmcDecayBranchesStore::const_iterator
1269  k( sObject.etaDecayTable.GetDecayBranches().begin() );
1270  k != sObject.etaDecayTable.GetDecayBranches().end(); ++k )
1271  {
1272  etaDecayTable->GetDecayChannel( k->first )->SetBR( k->second );
1273  }
1274 
1275  PrintReadRunData();
1276  }
1277 }
1278 
1279 
1280 #ifdef CEXMC_USE_CUSTOM_FILTER
1281 
1282 void CexmcRunManager::SetCustomFilter( const G4String & cfFileName_ )
1283 {
1284  if ( customFilter )
1285  {
1286  delete customFilter;
1287  customFilter = NULL;
1288  }
1289 
1290  if ( cfFileName_.empty() )
1291  return;
1292 
1293  /* should not get here */
1294  if ( ! ProjectIsRead() )
1295  throw CexmcException( CexmcCmdIsNotAllowed );
1296 
1297  cfFileName = cfFileName_;
1298 
1299  customFilter = new CexmcCustomFilterEval( cfFileName );
1300 }
1301 
1302 #endif
1303 
1304 #endif
1305 
1306 
1307 void CexmcRunManager::RegisterScenePrimitives( void )
1308 {
1309  G4VisManager * visManager( static_cast< G4VisManager * >(
1310  G4VVisManager::GetConcreteInstance() ) );
1311  if ( ! visManager )
1312  return;
1313 
1314  G4Scene * curScene( visManager->GetCurrentScene() );
1315  if ( ! curScene )
1316  return;
1317 
1318  /* G4Scene declarations lack this kind of typedef */
1319 #if G4VERSION_NUMBER < 960
1320  typedef std::vector< G4VModel * > MList;
1321 #else
1322  typedef std::vector< G4Scene::Model > MList;
1323 #endif
1324  const MList & mList( curScene->GetRunDurationModelList() );
1325 
1326  for ( MList::const_iterator k( mList.begin() ); k != mList.end(); ++k )
1327  {
1328 #if G4VERSION_NUMBER < 960
1329  const G4String & modelDesc( ( *k )->GetGlobalDescription() );
1330 #else
1331  const G4String & modelDesc( k->fpModel->GetGlobalDescription() );
1332 #endif
1333  if ( modelDesc == CexmcScenePrimitivesDescription )
1334  return;
1335  }
1336 
1337  CexmcSetup * setup( static_cast< CexmcSetup * >( userDetector ) );
1338  if ( ! setup )
1339  throw CexmcException( CexmcWeirdException );
1340 
1341  /* BEWARE: looks like G4Scene won't delete models from its lists upon
1342  * termination! Hence destructor of the new model won't be called */
1343  curScene->AddRunDurationModel( new CexmcScenePrimitives( setup ) );
1344 }
1345 
1346 
1347 void CexmcRunManager::BeamParticleChangeHook( void )
1348 {
1349  const CexmcEventAction * eventAction(
1350  static_cast< const CexmcEventAction * >( userEventAction ) );
1351  if ( ! eventAction )
1352  throw CexmcException( CexmcWeirdException );
1353 
1354  CexmcEventAction * theEventAction( const_cast< CexmcEventAction * >(
1355  eventAction ) );
1356  theEventAction->BeamParticleChangeHook();
1357 }
1358 
1359 
1360 void CexmcRunManager::SetupConstructionHook( void )
1361 {
1362 #ifdef CEXMC_USE_PERSISTENCY
1363  /* save gdml file */
1364  G4String cmd( "" );
1365  CexmcExceptionType exceptionType( CexmcSystemException );
1366 
1367  if ( zipGdmlFile )
1368  {
1369  if ( ProjectIsRead() )
1370  {
1371  cmd = G4String( "bzip2 " ) + projectsDir + "/" + rProject +
1372  gdmlFileExtension;
1373  }
1374  else
1375  {
1376  if ( ProjectIsSaved() )
1377  cmd = G4String( "bzip2 -c " ) + gdmlFileName + " > " +
1378  projectsDir + "/" + projectId + gdmlbz2FileExtension;
1379  }
1380  exceptionType = CexmcFileCompressException;
1381  }
1382  else
1383  {
1384  if ( ! ProjectIsRead() && ProjectIsSaved() )
1385  cmd = G4String( "cp " ) + gdmlFileName + " " + projectsDir + "/" +
1386  projectId + gdmlFileExtension;
1387  /* else already saved in ReadPreinitProjectData() */
1388  }
1389 
1390  if ( ! cmd.empty() && system( cmd ) != 0 )
1391  throw CexmcException( exceptionType );
1392 #endif
1393 }
1394 
Generated on Sat May 25 2013 14:32:07 for Geant4 by   doxygen 1.8.4