LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
G4Helper.cxx
Go to the documentation of this file.
1 
9 #include "nug4/G4Base/G4Helper.h"
12 
14 
15 #include "Geant4/G4UImanager.hh"
16 #include "Geant4/G4VUserDetectorConstruction.hh"
17 #include "Geant4/G4VUserPrimaryGeneratorAction.hh"
18 #include "Geant4/G4VUserPhysicsList.hh"
19 #include "Geant4/G4UserLimits.hh"
20 #include "Geant4/G4UserRunAction.hh"
21 #include "Geant4/G4UserEventAction.hh"
22 #include "Geant4/G4UserTrackingAction.hh"
23 #include "Geant4/G4UserSteppingAction.hh"
24 #include "Geant4/G4VisExecutive.hh"
25 #include "Geant4/G4StepLimiterPhysics.hh"
26 #include "Geant4/G4LogicalVolumeStore.hh"
27 
28 #include <boost/algorithm/string.hpp>
29 
30 #include "Geant4/QGSP_BERT.hh"
31 
32 #include "Geant4/G4VModularPhysicsList.hh"
33 
34 #include "Geant4/G4PhysListFactoryAlt.hh"
35 #include "Geant4/G4PhysicsConstructorRegistry.hh"
36 
37 #include <Rtypes.h>
38 
39 #include <iostream>
40 #include <cstring>
41 #include <sys/stat.h>
42 
44 
45 namespace g4b{
46 
47  static G4VisExecutive* vm_ = 0;
48 
49  //------------------------------------------------
50  // Constructor
52  : fCheckOverlaps (false )
53  , fValidateGDMLSchema(false )
54  , fUseStepLimits (false )
55  , fUIManager (nullptr)
56  , fConvertMCTruth (nullptr)
57  , fDetector (nullptr)
58  {
59  fParallelWorlds.clear();
60  }
61 
62  //------------------------------------------------
63  // Constructor
64  G4Helper::G4Helper(std::string const& g4macropath,
65  std::string const& g4physicslist,
66  std::string const& gdmlFile)
67  : fG4MacroPath (g4macropath )
68  , fG4PhysListName (g4physicslist)
69  , fGDMLFile (gdmlFile )
70  , fCheckOverlaps (false )
71  , fValidateGDMLSchema(true )
72  , fUseStepLimits (false )
73  , fUIManager (nullptr )
74  , fConvertMCTruth (nullptr )
75  , fDetector (nullptr )
76  {
77  // Geant4 run manager. Nothing happens in Geant4 until this object
78  // is created.
79  fRunManager = new G4RunManager;
80 
81  // Get the pointer to the User Interface manager
82  fUIManager = G4UImanager::GetUIpointer();
83 
84  fParallelWorlds.clear();
85  }
86 
87  //------------------------------------------------
88  // Destructor
90  {
91  if( vm_ ) delete vm_;
92 
93  if ( fRunManager != 0 ){
94  // In SetUserAction(), we set all the G4 user-action classes to be the
95  // same action: G4Base::UserActionManager This is convenient, but
96  // it creates a problem here: First the G4RunManager deletes the
97  // G4UserRunAction, then it tries to delete the
98  // G4UserEventAction... but that's the same object, which has
99  // already been deleted. Crash.
100 
101  // To keep this from happening, handle the UserActionManager
102  // clean-up manually, then tell the G4RunManager that all those
103  // classes no longer exist.
104 
106  bool wasStacking = uaManager->DoesAnyActionProvideStacking();
107  uaManager->Close();
108 
109  // Each one of these G4RunManager::SetUserAction methods calls a
110  // different method, based on the type of the argument. We want
111  // to use "0" (a null pointer), but we have to cast that "0" to a
112  // particular type in order for the right SetUserAction method to
113  // be called.
114 
115  fRunManager->SetUserAction( static_cast<G4UserRunAction*>(0) );
116  fRunManager->SetUserAction( static_cast<G4UserEventAction*>(0) );
117  fRunManager->SetUserAction( static_cast<G4UserTrackingAction*>(0) );
118  fRunManager->SetUserAction( static_cast<G4UserSteppingAction*>(0) );
119  if ( wasStacking ) {
120  fRunManager->SetUserAction( static_cast<G4UserStackingAction*>(0) );
121  }
122 
123  delete fRunManager;
124  }
125  else{
126  MF_LOG_ERROR("G4Helper")
127  << "G4Helper never initialized; probably because there were no input primary events";
128  }
129 
130  for(size_t i = 0; i < fParallelWorlds.size(); ++i){
131  if(fParallelWorlds[i]) delete fParallelWorlds[i];
132  }
133  fParallelWorlds.clear();
134 
135  }
136 
137  //------------------------------------------------
138  void G4Helper::SetPhysicsList(std::string physicsString)
139  {
158 
159  G4VUserPhysicsList* physics = 0;
160  std::string bywhom = "User";
161  std::string factoryname = "g4alt::G4PhysListFactory";
162  bool list_known_ctors = true;
163 
164  // physics list name is the first part, anything afterwards is
165  // extra physics processes to be added to the base list
166  // ie. "QGSP_BERT ; myspace::MonopolePhysics ; MyOtherSpecialPhysics "
167  std::vector< std::string > pstrings;
168  // don't use ":" as a separator because it's used in namespaces
169  boost::algorithm::split( pstrings, physicsString, boost::is_any_of(";"),
170  boost::token_compress_on );
171  // trim lead/trail space
172  for (unsigned int j=0; j < pstrings.size(); ++j )
173  boost::algorithm::trim(pstrings[j]);
174 
175  if ( pstrings.size() < 1 ) pstrings.push_back(""); // non-empty
176  std::string phListName = pstrings[0];
177 
178  //for (unsigned int j=0; j < pstrings.size(); ++j )
179  // std::cout << "G4Helper pstrings[" << j << "] = \""
180  // << pstrings[j] << "\"" << std::endl;
181 
182  // user extensible physics list factory
183  g4alt::G4PhysListFactory factory;
184  factoryname = "g4alt::G4PhysListFactory";
185 
186  if ( factory.IsReferencePhysList(phListName) ) {
187  bywhom = factoryname;
188  physics = factory.GetReferencePhysList(phListName);
189  }
190  else {
191  // in the case of non-default name
192  if ( phListName != "" ) {
193  std::cerr << std::endl << factoryname
194  << " failed to find ReferencePhysList \""
195  << phListName << "\"" << std::endl;
196  factory.PrintAvailablePhysLists();
197 /* #else
198  std::vector<G4String> list = factory.AvailablePhysLists();
199  MF_LOG_VERBATIM("G4Helper")
200  << "For reference: PhysicsLists in G4PhysListFactory are: ";
201  for (size_t indx=0; indx < list.size(); ++indx ) {
202  MF_LOG_VERBATIM("G4Helper")
203  << " [" << std::setw(2) << indx << "] "
204  << "\"" << list[indx] << "\"";
205  }
206 #endif */
207  }
208 
209  if ( ! physics ) {
210  MF_LOG_ERROR("G4Helper")
211  << bywhom << " could not construct \""
212  << phListName
213  << "\","
214  << std::endl
215  << "FAIL hard if we don't succeed"; //"fall back to using QGSP_BERT";
216 
217  physics = 0; //new QGSP_BERT; // FAIL if not succeed
218  phListName = "<none>"; // "QGSP_BERT";
219 
220  } else {
221  MF_LOG_VERBATIM("G4Helper")
222  << bywhom
223  << " constructed G4VUserPhysicsList \""
224  << phListName
225  << "\"";
226  }
227 
228  }
229 
230  // Extend the physics list with additional physics processes
231  // Already used pstrings[0] entry for physics list name.
232  // The rest should be semi-colon separated list of:
233  // physicsProcessName ( optional UI command , more UI commands )
234  for (unsigned int k=1; k < pstrings.size(); ++k ) {
235  std::string physProcAddition = pstrings[k];
236 
237  // break off UI commands from process name
238  std::vector< std::string > physProcParts;
239  boost::algorithm::split( physProcParts, physProcAddition,
240  boost::is_any_of("(,)"),
241  boost::token_compress_on );
242  // trim lead/trail spaces
243  for (unsigned int j=0; j < physProcParts.size(); ++j )
244  boost::algorithm::trim(physProcParts[j]);
245 
246  // element 0 is the physics process name
247  std::string physProcName = physProcParts[0];
248  if ( physProcName == "" ) continue; // not real, user has trailing ";"
249 
250  G4PhysicsConstructorRegistry* pcRegistry =
251  G4PhysicsConstructorRegistry::Instance();
252 
253  if ( ! pcRegistry->IsKnownPhysicsConstructor(physProcName) ) {
254 
255  MF_LOG_VERBATIM("G4Helper")
256  << "G4PhysicsProcessFactorySingleton could not "
257  << "construct a \""
258  << physProcName
259  << "\"";
260 
261  if ( ! list_known_ctors ) continue;
262  // make sure we only do this once
263  list_known_ctors = false;
264  std::vector<G4String> list = pcRegistry->AvailablePhysicsConstructors();
265  MF_LOG_VERBATIM("G4Helper")
266  << "For reference: PhysicsProcesses in "
267  << "G4PhysicsProcessFactorySingleton are: ";
268 
269  if ( list.empty() )
270  MF_LOG_VERBATIM("G4Helper")
271  << " ... no registered processes";
272  else {
273  for (size_t indx=0; indx < list.size(); ++indx ) {
274  MF_LOG_VERBATIM("G4Helper")
275  << " [" << std::setw(2) << indx << "] "
276  << "\"" << list[indx] << "\"";
277  }
278  }
279  continue;
280  }
281 
282  MF_LOG_VERBATIM("G4Helper")
283  << "Adding \""
284  << physProcName
285  << "\" physics process to \""
286  << phListName
287  << "\"";
288 
289  // construct physics process, add it to the base physics list
290  G4VPhysicsConstructor* pctor =
291  pcRegistry->GetPhysicsConstructor(physProcName);
292 
293 
294  G4VModularPhysicsList* mpl = dynamic_cast<G4VModularPhysicsList*>(physics);
295  if ( ! pctor ) MF_LOG_VERBATIM("G4Helper") << " ... failed with null pointer";
296  else if ( ! mpl ) MF_LOG_VERBATIM("G4Helper") << " ... failed, physics list wasn't a G4VModularPhysicsList";
297  else mpl->RegisterPhysics(pctor);
298 
299  // Handle associated UI commands
300  // One must do it here for cases where values need to be set *before*
301  // one calls SetUserInitialization(physics)
302 
303  for ( unsigned int i=1; i < physProcParts.size(); ++i ) {
304  if ( physProcParts[i] == "" ) continue;
305  MF_LOG_VERBATIM("G4Helper")
306  << physProcParts[i];
307 
308  fUIManager->ApplyCommand(physProcParts[i]);
309  }
310 
311  }
312 
313  // User should have called SetVolumeStepLimit before calling this
314  // method, otherwise fUseStepLimits is always false.
315  if(fUseStepLimits){
316  auto mpl = dynamic_cast<G4VModularPhysicsList*>(physics);
317  if(mpl) mpl->RegisterPhysics(new G4StepLimiterPhysics());
318  else
319  MF_LOG_WARNING("G4Helper")
320  << "Step limits requested, but unable to register G4StepLimiterPhysics"
321  << "\n NO STEP LIMITS WILL BE APPLIED";
322  }
323 
324  // pass off (possibly augmented) physics list to run manager
325  // which calls G4RunManagerKernel->SetPhysics() on it
326  // which itself call ConstructParticle() for the list
327  fRunManager->SetUserInitialization(physics);
328  }
329 
330  //------------------------------------------------
331  void G4Helper::SetParallelWorlds(std::vector<G4VUserParallelWorld*> pworlds)
332  {
333  for(auto const& pw : pworlds){
334  MF_LOG_DEBUG("G4Helper") << pw->GetName();
335  fParallelWorlds.push_back(pw);
336  }
337 
338  return;
339  }
340 
341  //------------------------------------------------
342  void G4Helper::ConstructDetector(std::string const& gdmlFile)
343  {
344  // Build the Geant4 detector description.
345  bool checkOverlaps = fCheckOverlaps;
346  bool validateGDMLSchema = fValidateGDMLSchema;
347  fDetector = new DetectorConstruction(gdmlFile,
348  checkOverlaps,
349  validateGDMLSchema);
350 
351  return;
352  }
353 
354  //------------------------------------------------
355  void G4Helper::SetVolumeStepLimit(std::string const& volumeName,
356  double maxStepSize)
357  {
358  // get the logical volume for the desired volume name
359  G4LogicalVolume* logVol = G4LogicalVolumeStore::GetInstance()->GetVolume(volumeName);
360 
361  if(logVol){
362  // the logical volume takes ownership of the G4UserLimits pointer
363  G4UserLimits *stepLimit = new G4UserLimits(maxStepSize);
364  logVol->SetUserLimits(stepLimit);
365 
366  fUseStepLimits = true;
367  }
368  else{
369  MF_LOG_WARNING("G4Helper")
370  << "Unable to find volume "
371  << volumeName
372  << " and set step size limit";
373  }
374 
375  return;
376  }
377 
378 
379  //------------------------------------------------
382  {
384 
385  for(auto pWorld : fParallelWorlds)
386  fDetector->RegisterParallelWorld(pWorld);
387 
388  // define the physics list to use
390 
391  // Pass the detector geometry on to Geant4.
392  fRunManager->SetUserInitialization(fDetector);
393 
394  // Tell the Geant4 run manager how to generate events. The
395  // ConvertMCTruthToG4 class will "generate" events by
396  // converting MCTruth objects from the input into G4Events.
398  fRunManager->SetUserAction(fConvertMCTruth);
399 
400  return;
401  }
402 
403 
404  //------------------------------------------------
407  {
408  // Geant4 comes with "user hooks" that allows users to perform
409  // special tasks at the beginning and end of runs, events, tracks,
410  // steps. By using the UserActionManager, we've separated each
411  // set of user tasks into their own class; e.g., there can be one
412  // class for processing particles, one class for histograms, etc.
413 
414  // Use the UserActionManager to handle all the Geant4 user hooks.
416 
417  // Tell the run manager about our user-action classes. We convert
418  // the UserActionManager into different types so Geant4's run and
419  // event managers will initialize them properly.
420  G4UserRunAction* runAction = (G4UserRunAction* ) uaManager;
421  G4UserEventAction* eventAction = (G4UserEventAction* ) uaManager;
422  G4UserTrackingAction* trackingAction = (G4UserTrackingAction*) uaManager;
423  G4UserSteppingAction* steppingAction = (G4UserSteppingAction*) uaManager;
424  fRunManager->SetUserAction( runAction );
425  fRunManager->SetUserAction( eventAction );
426  fRunManager->SetUserAction( trackingAction );
427  fRunManager->SetUserAction( steppingAction );
428 
429  if ( uaManager->DoesAnyActionProvideStacking() ) {
430  G4UserStackingAction* stackingAction = (G4UserStackingAction*) uaManager;
431  fRunManager->SetUserAction( stackingAction );
432  }
433 
436  fRunManager->Initialize();
437 
438  if(!vm_) vm_ = new G4VisExecutive();
439  vm_->Initialize();
440 
441  // Tell the manager to execute the contents of the Geant4 macro
442  // file.
443  if ( ! fG4MacroPath.empty() ) {
444  G4String command = "/control/execute " + G4String(fG4MacroPath);
445  fUIManager->ApplyCommand(command);
446  }
447 
448 
449  return;
450  }
451 
452  //------------------------------------------------
454  {
455  return this->G4Run( primary.get() );
456  }
457 
458  //------------------------------------------------
459  bool G4Helper::G4Run(const simb::MCTruth* primary)
460  {
461  // Get the event converter ready.
463 
464  // Pass the MCTruth to our event generator.
465  fConvertMCTruth->Append( primary );
466 
467  // Start the simulation for this event. Note: The following
468  // statement increments the G4RunManager's run number. Because of
469  // this, it's important for events to use the run/event number
470  // from the EventDataModel Header, not G4's internal numbers.
471  fUIManager->ApplyCommand("/run/beamOn 1");
472 
473  return true;
474  }
475 
476  //------------------------------------------------
477  bool G4Helper::G4Run(std::vector<const simb::MCTruth*> &primaries)
478  {
479  // Get the event converter ready.
481 
482  // Pass all the MCTruths to our event generator.
483  for(auto primary : primaries)
484  fConvertMCTruth->Append( primary );
485 
486  // Start the simulation for this event. Note: The following
487  // statement increments the G4RunManager's run number. Because of
488  // this, it's important for events to use the run/event number
489  // from the EventDataModel Header, not G4's internal numbers.
490  fUIManager->ApplyCommand("/run/beamOn 1");
491 
492  return true;
493  }
494 
495 } // namespace
Build Geant4 geometry from GDML.
virtual ~G4Helper()
Definition: G4Helper.cxx:89
virtual bool DoesAnyActionProvideStacking()
void Reset()
Get ready to convert a new set of MCTruth objects.
bool fValidateGDMLSchema
Have G4GDML validate geometry schema?
Definition: G4Helper.h:114
bool G4Run(std::vector< const simb::MCTruth * > &primaries)
Definition: G4Helper.cxx:477
void SetUserAction()
Initialization for the Geant4 Monte Carlo.
Definition: G4Helper.cxx:406
void ConstructDetector(std::string const &gdmlFile)
Definition: G4Helper.cxx:342
#define MF_LOG_ERROR(category)
void SetParallelWorlds(std::vector< G4VUserParallelWorld * > pworlds)
Definition: G4Helper.cxx:331
std::vector< G4VUserParallelWorld * > fParallelWorlds
list of parallel worlds
Definition: G4Helper.h:122
Use Geant4 to run the LArSoft detector simulation.
void InitPhysics()
Initialization for the Geant4 Monte Carlo.
Definition: G4Helper.cxx:381
void SetPhysicsList(std::string physicsList)
Definition: G4Helper.cxx:138
static G4VisExecutive * vm_
Definition: G4Helper.cxx:47
std::string fGDMLFile
Name of the gdml file containing the detector Geometry.
Definition: G4Helper.h:112
bool fCheckOverlaps
Have G4GDML check for overlaps?
Definition: G4Helper.h:113
basic interface to Geant4 for ART-based software
void SetVolumeStepLimit(std::string const &volumeName, double maxStepSize)
Definition: G4Helper.cxx:355
std::string fG4MacroPath
to be executed before main MC processing.
Definition: G4Helper.h:109
DetectorConstruction * fDetector
DetectorConstruction object.
Definition: G4Helper.h:121
G4Helper()
Standard constructor and destructor for an FMWK module.
Definition: G4Helper.cxx:51
static UserActionManager * Instance()
#define MF_LOG_VERBATIM(category)
G4RunManager * fRunManager
Geant4&#39;s run manager.
Definition: G4Helper.h:117
std::string fG4PhysListName
Name of physics list to use.
Definition: G4Helper.h:111
#define MF_LOG_DEBUG(id)
void Append(art::Ptr< simb::MCTruth > &mct)
Event generator information.
Definition: MCTruth.h:32
bool fUseStepLimits
Set in SetVolumeStepLimit.
Definition: G4Helper.h:115
#define MF_LOG_WARNING(category)
T const * get() const
Definition: Ptr.h:138
ConvertMCTruthToG4 * fConvertMCTruth
Geant4 event generator.
Definition: G4Helper.h:119
G4UImanager * fUIManager
Geant4&#39;s user-interface manager.
Definition: G4Helper.h:118