BDSIMClass.cc

/* 
Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway, 
University of London 2001 - 2022.
 
This file is part of BDSIM.
 
BDSIM is free software: you can redistribute it and/or modify 
it under the terms of the GNU General Public License as published 
by the Free Software Foundation version 3 of the License.
 
BDSIM is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
 
You should have received a copy of the GNU General Public License
along with BDSIM.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "BDSIMClass.hh"
 
#include "BDSExecOptions.hh"     // executable command line options 
#include "BDSGlobalConstants.hh" //  global parameters
 
#include <algorithm>
#include <csignal>
#include <cstdlib>
#include <cstdio>
 
#include "G4EventManager.hh" // Geant4 includes
#include "G4GenericBiasingPhysics.hh"
#include "G4GeometryManager.hh"
#include "G4GeometryTolerance.hh"
#include "G4PhysicsListHelper.hh"
#include "G4ParallelWorldPhysics.hh"
#include "G4ParticleDefinition.hh"
#include "G4SteppingManager.hh"
#include "G4TrackingManager.hh"
#include "G4Version.hh"
#include "G4VModularPhysicsList.hh"
 
#include "CLHEP/Units/SystemOfUnits.h"
 
#include "BDSAcceleratorModel.hh"
#include "BDSAperturePointsLoader.hh"
#include "BDSBeamPipeFactory.hh"
#include "BDSBunch.hh"
#include "BDSBunchFactory.hh"
#include "BDSCavityFactory.hh"
#include "BDSColours.hh"
#include "BDSComponentFactoryUser.hh"
#include "BDSDebug.hh"
#include "BDSDetectorConstruction.hh"   
#include "BDSEventAction.hh"
#include "BDSException.hh"
#include "BDSFieldFactory.hh"
#include "BDSFieldLoader.hh"
#include "BDSGeometryFactory.hh"
#include "BDSGeometryFactorySQL.hh"
#include "BDSGeometryWriter.hh"
#include "BDSIonDefinition.hh"
#include "BDSMaterials.hh"
#include "BDSOutput.hh" 
#include "BDSOutputFactory.hh"
#include "BDSParallelWorldUtilities.hh"
#include "BDSParser.hh" // Parser
#include "BDSParticleDefinition.hh"
#include "BDSPhysicsMuonSplitting.hh"
#include "BDSPhysicsUtilities.hh"
#include "BDSPrimaryGeneratorAction.hh"
#include "BDSRandom.hh" // for random number generator from CLHEP
#include "BDSRunAction.hh"
#include "BDSRunManager.hh"
#include "BDSSamplerRegistry.hh"
#include "BDSSDManager.hh"
#include "BDSSteppingAction.hh"
#include "BDSStackingAction.hh"
#include "BDSTemporaryFiles.hh"
#include "BDSTrackingAction.hh"
#include "BDSUtilities.hh"
#include "BDSVisManager.hh"
#include "BDSWarning.hh"
 
BDSIM::BDSIM():
  ignoreSIGINT(false),
  usualPrintOut(true),
  initialised(false),
  initialisationResult(1),
  argcCache(0),
  argvCache(nullptr),
  parser(nullptr),
  bdsOutput(nullptr),
  bdsBunch(nullptr),
  runManager(nullptr),
  userComponentFactory(nullptr),
  userPhysicsList(nullptr),
  realWorld(nullptr)
{;}
 
BDSIM::BDSIM(int argc, char** argv, bool usualPrintOutIn):
  ignoreSIGINT(false),
  usualPrintOut(usualPrintOutIn),
  initialised(false),
  initialisationResult(1),
  argcCache(argc),
  argvCache(argv),
  parser(nullptr),
  bdsOutput(nullptr),
  bdsBunch(nullptr),
  runManager(nullptr),
  userComponentFactory(nullptr),
  userPhysicsList(nullptr),
  realWorld(nullptr)
{
  initialisationResult = Initialise();
}
 
int BDSIM::Initialise(int argc, char** argv, bool usualPrintOutIn)
{
  argcCache = argc;
  argvCache = argv;
  usualPrintOut = usualPrintOutIn;
  initialisationResult = Initialise();
  return initialisationResult;
}
 
int BDSIM::Initialise()
{
  const BDSExecOptions* execOptions = new BDSExecOptions(argcCache,argvCache);
  if (usualPrintOut)
    {execOptions->Print();}
  ignoreSIGINT = execOptions->IgnoreSIGINT(); // different sig catching for cmake
 
  execOptions->PrintCopyright();  
#ifdef BDSDEBUG
  G4cout << __METHOD_NAME__ << "DEBUG mode is on." << G4endl;
#endif
 
  parser = BDSParser::Instance(execOptions->InputFileName());
  parser->AmalgamateOptions(execOptions->Options());
  parser->AmalgamateBeam(execOptions->Beam(), execOptions->Options().recreate);
  parser->CheckOptions();
  
  BDSMaterials::Instance()->PrepareRequiredMaterials(execOptions->Options().verbose);
 
  delete execOptions; 
 
  BDSGlobalConstants* globals = BDSGlobalConstants::Instance();
 
  BDSRandom::CreateRandomNumberGenerator(globals->RandomEngine());
  BDSRandom::SetSeed(); // set the seed from options
 
  bdsOutput = BDSOutputFactory::CreateOutput(globals->OutputFormat(),
                                             globals->OutputFileName());
 
  if (!BDS::Geant4EnvironmentIsSet())
    {
      G4cerr << "No Geant4 environmental variables found - please source geant4.sh environment" << G4endl;
      G4cout << "A common fault is the wrong Geant4 environment as compared to the one BDSIM was compiled with." << G4endl;
      return 1;
    }
 
  runManager = new BDSRunManager();
 
  realWorld = new BDSDetectorConstruction(userComponentFactory);
  
  auto parallelWorldsRequiringPhysics = BDS::ConstructAndRegisterParallelWorlds(realWorld,
                                                                                realWorld->BuildSamplerWorld(),
                                                                                realWorld->BuildPlacementFieldsWorld());
  runManager->SetUserInitialization(realWorld);  
 
#ifdef BDSDEBUG
  G4cout << __METHOD_NAME__ << "> Constructing physics processes" << G4endl;
#endif
  G4String physicsListName = parser->GetOptions().physicsList;
 
#if G4VERSION_NUMBER > 1049
  // from 10.5 onwards they have a looping particle killer that warnings and kills particles
  // deemed to be looping that are <100 MeV. This is unrelated to the primary energy so troublesome.
  // set to the 'low' limits here ~10keV. This must be done before any physics is created as the
  // parameters are copied into the transportation physics process for each particle and it's very
  // hard to sift through and fix afterwards
  G4PhysicsListHelper::GetPhysicsListHelper()->UseLowLooperThresholds();
#endif
  // sampler physics process for parallel world tracking must be instantiated BEFORE
  // regular physics.
  // Note, we purposively don't create a parallel world process for the curvilinear
  // world as we don't need the track information from it - unreliable that way. We
  // query the geometry directly using our BDSAuxiliaryNavigator class.
  auto parallelWorldPhysics = BDS::ConstructParallelWorldPhysics(parallelWorldsRequiringPhysics);
  G4int physicsVerbosity = globals->PhysicsVerbosity();
  G4VModularPhysicsList* physList;
  if (userPhysicsList)
    {
      G4cout << "Using externally registered user defined physics list" << G4endl;
      physList = userPhysicsList;
    }
  else
    {physList = BDS::BuildPhysics(physicsListName, physicsVerbosity);}
 
  // create geometry sampler and register importance sampling biasing. Has to be here
  // before physicsList is "initialised" in run manager.
  if (BDSGlobalConstants::Instance()->UseImportanceSampling())
    {BDS::RegisterImportanceBiasing(parallelWorldsRequiringPhysics,physList);}
 
  // Construction of the physics lists defines the necessary particles and therefore
  // we can calculate the beam rigidity for the particle the beam is designed w.r.t. This
  // must happen before the geometry is constructed (which is called by
  // runManager->Initialize()).
  BDSParticleDefinition* designParticle = nullptr;
  BDSParticleDefinition* beamParticle = nullptr;
  G4bool beamDifferentFromDesignParticle = false;
  BDS::ConstructDesignAndBeamParticle(BDSParser::Instance()->GetBeam(),
                                      globals->FFact(),
                                      designParticle,
                                      beamParticle,
                                      beamDifferentFromDesignParticle);
  G4double minEK = globals->MinimumKineticEnergy();
  if (beamParticle->KineticEnergy() < minEK && BDS::IsFinite(minEK))
    {throw BDSException("option, minimumKineticEnergy is higher than kinetic energy of the beam - all primary particles wil be killed!");}
  if (usualPrintOut)
    {
      G4cout << "Design particle properties: " << G4endl << *designParticle;
      if (beamDifferentFromDesignParticle)
        {G4cout << "Beam particle properties: " << G4endl << *beamParticle;}
    }
  // update rigidity where needed
  realWorld->SetDesignParticle(designParticle);
  BDSFieldFactory::SetDesignParticle(designParticle);
  BDSGeometryFactorySQL::SetDefaultRigidity(designParticle->BRho()); // used for sql field loading
  
  // Muon splitting - optional - should be done *after* biasing to work with it
  G4int muonSplittingFactor = globals->MuonSplittingFactor();
  if (muonSplittingFactor > 1)
    {
      G4int muonSplittingFactor2 = globals->MuonSplittingFactor2();
      G4double muonSplittingThresholdParentEk = globals->MuonSplittingThresholdParentEk();
      G4double muonSplittingThresholdParentEk2 = globals->MuonSplittingThresholdParentEk2();
      G4cout << "BDSPhysicsMuonSplitting -> using muon splitting wrapper -> factor of: " << muonSplittingFactor << G4endl;
      if (muonSplittingThresholdParentEk > 0)
        {G4cout << "BDSPhysicsMuonSplitting -> minimum parent kinetic energy: " << muonSplittingThresholdParentEk / CLHEP::GeV << " GeV" << G4endl;}
      if (muonSplittingFactor2 > 1)
        {
          G4cout << "BDSPhysicsMuonSplitting -> factor #2: " << muonSplittingFactor2 << " for muons above "
                 << muonSplittingThresholdParentEk / CLHEP::GeV << " GeV" << G4endl;
        }
      G4bool excludeW1P = globals->MuonSplittingExcludeWeight1Particles();
      physList->RegisterPhysics(new BDSPhysicsMuonSplitting(muonSplittingFactor,  muonSplittingThresholdParentEk,
                                                            muonSplittingFactor2, muonSplittingThresholdParentEk2,
                                                            excludeW1P, globals->MuonSplittingExclusionWeight()));
    }
  
  BDS::RegisterSamplerPhysics(parallelWorldPhysics, physList);
  auto biasPhysics = BDS::BuildAndAttachBiasWrapper(parser->GetBiasing());
  if (biasPhysics)//could be nullptr and can't be passed to geant4 like this
    {physList->RegisterPhysics(biasPhysics);}
  runManager->SetUserInitialization(physList);
 
  bdsBunch = BDSBunchFactory::CreateBunch(beamParticle,
                                          parser->GetBeam(),
                                          globals->BeamlineTransform(),
                                          globals->BeamlineS(),
                                          globals->GeneratePrimariesOnly());
  G4cout << "Bunch distribution: " << bdsBunch->Name() << G4endl;
  delete beamParticle;
  if (bdsBunch->ExpectChangingParticleType())
    {BDS::ConstructExtendedParticleSet();}
  
  if (globals->GeneratePrimariesOnly())
    {      
      // output creation is duplicated below but with this if loop, we exit so ok.
      bdsOutput->NewFile();
      const G4int nToGenerate = globals->NGenerate();
      const G4int printModulo = globals->PrintModuloEvents();
      bdsBunch->BeginOfRunAction(nToGenerate);
      auto flagsCache(G4cout.flags());
      for (G4int i = 0; i < nToGenerate; i++)
        {
          if (i%printModulo == 0)
            {G4cout << "\r Primary> " << std::fixed << i << " of " << nToGenerate << G4endl;}
          BDSParticleCoordsFullGlobal coords = bdsBunch->GetNextParticleValid();
          // always pull particle definition in case it's updated
          const BDSParticleDefinition* pDef = bdsBunch->ParticleDefinition();
          bdsOutput->FillEventPrimaryOnly(coords, pDef);
        }
      G4cout.flags(flagsCache); // restore cout flags
      // Write options now the file is open
      const GMAD::OptionsBase* ob = BDSParser::Instance()->GetOptionsBase();
      bdsOutput->FillOptions(ob);
      
      // Write beam
      const GMAD::BeamBase* bb = BDSParser::Instance()->GetBeamBase();
      bdsOutput->FillBeam(bb);
 
      bdsOutput->CloseFile();
      return 0;
    }
  
  G4GeometryTolerance* theGeometryTolerance = G4GeometryTolerance::GetInstance();
  if (usualPrintOut)
    {
      G4cout << __METHOD_NAME__ << "Geometry Tolerances: "     << G4endl;
      G4cout << __METHOD_NAME__ << std::setw(12) << "Surface: " << std::setw(7) << theGeometryTolerance->GetSurfaceTolerance() << " mm"   << G4endl;
      G4cout << __METHOD_NAME__ << std::setw(12) << "Angular: " << std::setw(7) << theGeometryTolerance->GetAngularTolerance() << " rad"  << G4endl;
      G4cout << __METHOD_NAME__ << std::setw(12) << "Radial: "  << std::setw(7) << theGeometryTolerance->GetRadialTolerance()  << " mm"   << G4endl;
    }
#ifdef BDSDEBUG 
  G4cout << __METHOD_NAME__ << "Registering user action - Event Action" << G4endl;
#endif
  BDSEventAction* eventAction = new BDSEventAction(bdsOutput);
  runManager->SetUserAction(eventAction);
 
#ifdef BDSDEBUG 
  G4cout << __METHOD_NAME__ << "Registering user action - Run Action"<<G4endl;
#endif
  runManager->SetUserAction(new BDSRunAction(bdsOutput,
                                             bdsBunch,
                                             bdsBunch->ParticleDefinition()->IsAnIon(),
                                             eventAction,
                                             globals->StoreTrajectorySamplerID()));
 
#ifdef BDSDEBUG 
  G4cout << __METHOD_NAME__ << "Registering user action - Stepping Action"<<G4endl;
#endif
  // Only add steppingaction if it is actually used, so do check here (for performance reasons)
  G4int verboseSteppingEventStart = globals->VerboseSteppingEventStart();
  G4int verboseSteppingEventStop  = BDS::VerboseEventStop(verboseSteppingEventStart,
                                                          globals->VerboseSteppingEventContinueFor());
  if (globals->VerboseSteppingBDSIM())
    {
      runManager->SetUserAction(new BDSSteppingAction(true,
                                                      verboseSteppingEventStart,
                                                      verboseSteppingEventStop));
    }
  
#ifdef BDSDEBUG 
  G4cout << __METHOD_NAME__ << "Registering user action - Tracking Action"<<G4endl;
#endif
  runManager->SetUserAction(new BDSTrackingAction(globals->Batch(),
                                                  globals->StoreTrajectory(),
                                                  globals->StoreTrajectoryOptions(),
                                                  eventAction,
                                                  verboseSteppingEventStart,
                                                  verboseSteppingEventStop,
                                                  globals->VerboseSteppingPrimaryOnly(),
                                                  globals->VerboseSteppingLevel()));
 
#ifdef BDSDEBUG 
  G4cout << __METHOD_NAME__ << "Registering user action - Stacking Action"<<G4endl;
#endif
  runManager->SetUserAction(new BDSStackingAction(globals));
 
#ifdef BDSDEBUG 
  G4cout << __METHOD_NAME__ << "Registering user action - Primary Generator"<<G4endl;
#endif
  auto primaryGeneratorAction = new BDSPrimaryGeneratorAction(bdsBunch, parser->GetBeam());
  runManager->SetUserAction(primaryGeneratorAction);
  BDSFieldFactory::SetPrimaryGeneratorAction(primaryGeneratorAction);
 
  runManager->Initialize();
 
  if (globals->UseImportanceSampling())
    {BDS::AddIStore(parallelWorldsRequiringPhysics);}
 
  realWorld->BuildPhysicsBias();
 
  if (usualPrintOut && globals->PhysicsVerbose())
    {
      BDS::PrintPrimaryParticleProcesses(bdsBunch->ParticleDefinition()->Name());
      BDS::PrintDefinedParticles();
    }
 
  runManager->SetVerboseLevel(std::min(globals->VerboseRunLevel(), globals->PhysicsVerbosity()));
  G4EventManager::GetEventManager()->SetVerboseLevel(globals->VerboseEventLevel());
  G4EventManager::GetEventManager()->GetTrackingManager()->SetVerboseLevel(globals->VerboseTrackingLevel());
  
  G4bool bCloseGeometry = G4GeometryManager::GetInstance()->CloseGeometry();
  if (!bCloseGeometry)
    { 
      G4cerr << __METHOD_NAME__ << "error - geometry not closed." << G4endl;
      return 1;
    }
 
  if (globals->ExportGeometry())
    {
      BDSGeometryWriter geometrywriter;
      geometrywriter.ExportGeometry(globals->ExportType(),
                                    globals->ExportFileName());
    }
 
  initialised = true;
  return 0;
}
 
void BDSIM::BeamOn(int nGenerate)
{
  if (initialisationResult > 1 || !initialised)
    {return;} // a mode where we don't do anything
 
  G4cout.precision(10);
  struct sigaction act;
  act.sa_handler = &BDS::HandleAborts;
  sigemptyset(&act.sa_mask);
  act.sa_flags = 0;
  if (!ignoreSIGINT)
    {sigaction(SIGINT,  &act, nullptr);}
  sigaction(SIGABRT, &act, nullptr);
  sigaction(SIGTERM, &act, nullptr);
  sigaction(SIGSEGV, &act, nullptr);
  
  try
    {
      if (!BDSGlobalConstants::Instance()->Batch())   // Interactive mode
        {
          BDSVisManager visManager = BDSVisManager(BDSGlobalConstants::Instance()->VisMacroFileName(),
                                                   BDSGlobalConstants::Instance()->Geant4MacroFileName(),
                                                   realWorld);
          visManager.StartSession(argcCache, argvCache);
        }
      else
        {// batch mode
          if (nGenerate < 0)
            {runManager->BeamOn(BDSGlobalConstants::Instance()->NGenerate());}
          else
            {runManager->BeamOn(nGenerate);}
        }
    }
  catch (const BDSException& exception)
    {
      // don't do this for now in case it's dangerous and we try tracking with open geometry
      //G4GeometryManager::GetInstance()->OpenGeometry();
      throw exception;
    } 
}
 
BDSIM::~BDSIM()
{
  G4GeometryManager::GetInstance()->OpenGeometry();
    
#ifdef BDSDEBUG
  G4cout << __METHOD_NAME__ << "deleting..." << G4endl;
#endif
  delete bdsOutput;
 
  try
    {
      // order important here because of singletons relying on each other
          delete BDSSDManager::Instance();
      delete BDSBeamPipeFactory::Instance();
      delete BDSCavityFactory::Instance();
      delete BDSGeometryFactory::Instance();
      delete BDSAcceleratorModel::Instance();
      delete BDSTemporaryFiles::Instance();
      delete BDSFieldFactory::Instance(); // this uses BDSGlobalConstants which uses BDSMaterials
      delete BDSGlobalConstants::Instance();
      delete BDSMaterials::Instance();
      
      // instances not used in this file, but no other good location for deletion
      if (initialisationResult < 2)
        {
          delete BDSColours::Instance();
          delete BDSFieldLoader::Instance();
          delete BDSSamplerRegistry::Instance();
    BDSAperturePointsCache::Instance()->ClearCachedFiles();
        }
    }
  catch (...)
    {;} // ignore any exception as this is a destructor
      
  delete runManager;
  delete bdsBunch;
  delete parser;
 
  if (usualPrintOut)
    {G4cout << __METHOD_NAME__ << "End of Run. Thank you for using BDSIM!" << G4endl;}
}
 
void BDSIM::RegisterUserComponent(const G4String& componentTypeName,
                                  BDSComponentConstructor* componentConstructor)
{
  if (initialised)
    {BDS::Warning(__METHOD_NAME__, "BDSIM kernel already initialised - this component will not be available");}
  if (!userComponentFactory)
    {userComponentFactory = new BDSComponentFactoryUser();}
 
  userComponentFactory->RegisterComponent(componentTypeName,
                                          componentConstructor);
}