19#include "G4Version.hh"
22#include "BDSException.hh"
23#include "BDSGlobalConstants.hh"
24#include "BDSModularPhysicsList.hh"
25#include "BDSIonDefinition.hh"
26#include "BDSParticleDefinition.hh"
27#if G4VERSION_NUMBER > 1039
28#include "BDSPhysicsChannelling.hh"
30#include "BDSPhysicsCutsAndLimits.hh"
31#include "BDSPhysicsEMDissociation.hh"
32#include "BDSPhysicsUtilities.hh"
33#include "BDSUtilities.hh"
34#include "BDSWarning.hh"
35#include "BDSEmStandardPhysicsOp4Channelling.hh"
37#include "FTFP_BERT.hh"
39#include "G4AntiNeutrinoE.hh"
40#include "G4AntiNeutrinoMu.hh"
41#include "G4AntiNeutrinoTau.hh"
42#include "G4AntiNeutron.hh"
43#include "G4AntiProton.hh"
44#include "G4Electron.hh"
45#include "G4EmParameters.hh"
46#include "G4EmStandardPhysics_option4.hh"
47#include "G4EmStandardPhysicsSS.hh"
48#include "G4DynamicParticle.hh"
50#include "G4GenericBiasingPhysics.hh"
51#include "G4GenericIon.hh"
52#include "G4IonElasticPhysics.hh"
53#include "G4IonTable.hh"
54#include "G4KaonMinus.hh"
55#include "G4KaonPlus.hh"
56#include "G4MuonMinus.hh"
57#include "G4KaonZeroLong.hh"
58#include "G4LeptonConstructor.hh"
59#include "G4MuonPlus.hh"
60#include "G4NeutrinoE.hh"
61#include "G4NeutrinoMu.hh"
62#include "G4NeutrinoTau.hh"
63#include "G4Neutron.hh"
64#include "G4ParticleTable.hh"
65#include "G4ParticleTableIterator.hh"
66#include "G4PionMinus.hh"
67#include "G4PionPlus.hh"
68#include "G4PionZero.hh"
69#include "G4Positron.hh"
70#include "G4ProductionCutsTable.hh"
72#include "G4PhysListFactory.hh"
73#include "G4ProcessManager.hh"
74#include "G4ProcessVector.hh"
77#include "G4UImanager.hh"
78#if G4VERSION_NUMBER > 1049
79#include "G4ParticleDefinition.hh"
80#include "G4CoupledTransportation.hh"
81#include "G4Transportation.hh"
85#if G4VERSION_NUMBER > 1069
86#include "G4HadronicParameters.hh"
89#include "parser/beam.h"
90#include "parser/fastlist.h"
91#include "parser/physicsbiasing.h"
102 return G4IonTable::IsIon(particle) && particle!=G4Proton::Definition();
107 return BDS::IsIon(particle->GetDefinition()) || particle->GetTotalOccupancy()>0;
116 G4VModularPhysicsList* result =
nullptr;
123 G4bool completePhysics =
BDS::StrContains(physicsListNameLower,
"complete");
124 if (useGeant4Physics)
127 G4String geant4PhysicsList = physicsList.substr(2);
128 G4PhysListFactory factory;
129 if (!factory.IsReferencePhysList(geant4PhysicsList))
131 G4cerr <<
"Unknown Geant4 physics list \"" << geant4PhysicsList <<
"\"" << G4endl;
132 G4cout <<
"Available Geant4 hadronic physics lists:" << G4endl;
133 for (
const auto &name : factory.AvailablePhysLists())
134 {G4cout <<
"\"" << name <<
"\"" << G4endl;}
135 G4cout <<
"Available Geant4 EM physics lists:" << G4endl;
136 for (
const auto &name : factory.AvailablePhysListsEM())
137 {G4cout <<
"\"" << name <<
"\"" << G4endl;}
138 throw BDSException(__METHOD_NAME__,
"Unknown Geant4 physics list \"" + geant4PhysicsList +
"\"");
142 result = factory.GetReferencePhysList(geant4PhysicsList);
143 if (g->G4PhysicsUseBDSIMRangeCuts())
145 if (g->MinimumKineticEnergy() > 0 || g->G4PhysicsUseBDSIMCutsAndLimits())
147 G4cout <<
"\nAdding cuts and limits physics process to Geant4 reference physics list" << G4endl;
148 G4cout <<
"This is to enforce BDSIM range cuts and the minimumKinetic energy option.\n";
149 G4cout <<
"This is done by default for the functionality of BDSIM tracking and should not affect the physics greatly.\n";
150 G4cout <<
"See the BDSIM manual about Geant4 reference physics lists for details." << G4endl;
153 else if (!g->G4PhysicsUseBDSIMCutsAndLimits() && g->Circular())
155 G4String message =
"g4PhysicsUseBDSIMCutsAndLimits turned off but using a circular machine - circular mechanics will be broken";
156 BDS::Warning(__METHOD_NAME__, message);
160 else if (completePhysics)
164 G4cout <<
"Constructing \"" << physicsListNameLower <<
"\" complete physics list" << G4endl;
165#if G4VERSION_NUMBER > 1039
173 r->SetVerboseLevel(verbosity);
176 throw BDSException(__METHOD_NAME__,
"Channel physics is not supported with Geant4 versions less than 10.4");
180 {
throw BDSException(__METHOD_NAME__,
"Unknown 'complete' physics list \"" + physicsList +
"\"");}
190 result->ConstructParticle();
191 result->SetVerboseLevel(verbosity);
194 G4UImanager* UIManager = G4UImanager::GetUIpointer();
196 if (!physicsMacro.empty())
199 G4String physicsMacroFull =
BDS::GetFullPath(physicsMacro,
false, setInExecOptions);
200 G4cout <<
"Applying geant4 physics macro file: " << physicsMacroFull << G4endl;
201 UIManager->ApplyCommand(
"/control/execute " + physicsMacroFull);
204 G4VUserPhysicsList* resultAsUserPhysicsList =
dynamic_cast<G4VUserPhysicsList*
>(result);
205 if (resultAsUserPhysicsList)
207 resultAsUserPhysicsList->DumpCutValuesTable(verbosity);
208 resultAsUserPhysicsList->SetVerboseLevel(verbosity);
214 const std::set<std::string>& keys)
217 for (
const auto& k : keys)
218 {nSet += beamDefinition.
HasBeenSet(k) ? 1 : 0;}
223 const std::set<std::string>& keys,
225 G4bool warnZeroParamsSet,
226 const G4String& unitString)
230 G4cerr <<
"Beam> More than one parameter set - there should only be one" << G4endl;
231 for (
const auto& k : keys)
232 {G4cerr << std::setw(14) << std::left << k <<
": " << std::setw(7) << std::right << beamDefinition.
get_value(k) <<
" " << unitString << G4endl;}
233 throw BDSException(__METHOD_NAME__,
"conflicting parameters set");
235 else if (nSet == 0 && warnZeroParamsSet)
237 G4cerr <<
"Beam> One of the following required to be set" << G4endl;
238 for (
const auto &k : keys)
239 {G4cerr << std::setw(14) << std::left << k <<
": " << std::setw(7) << std::right << beamDefinition.
get_value(k) <<
" " << unitString << G4endl;}
240 throw BDSException(__METHOD_NAME__,
"insufficient parameters set");
248 G4bool& beamDifferentFromDesignParticle)
250 if (beamDefinition.
particle.empty())
251 {
throw BDSException(
"Beam> no \"particle\" specified (required).");}
254 std::set<std::string> keysDesign = {
"energy",
"momentum",
"kineticEnergy"};
257 std::set<std::string> keysParticle = {
"E0",
"P0",
"Ek0"};
268 beamDifferentFromDesignParticle = nSetParticle > 0 || beamParticleName != beamDefinition.
particle;
269 if (nSetParticle > 0)
273 beamDefinition.
E0 * CLHEP::GeV,
274 beamDefinition.
Ek0 * CLHEP::GeV,
275 beamDefinition.
P0 * CLHEP::GeV,
293 G4double totalEnergyIn,
294 G4double kineticEnergyIn,
301 std::map<G4String, G4String> commonSubstitutions = { {
"photon",
"gamma"},
308 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
309 std::regex ionParticle(
"(ion\\s).*");
310 if (std::regex_match(particleName, ionParticle))
312 G4GenericIon::GenericIonDefinition();
315 G4IonTable* ionTable = particleTable->GetIonTable();
317 G4int ionPDGID = G4IonTable::GetNucleusEncoding(ionDef->Z(), ionDef->A());
318 G4double mass = ionTable->GetIonMass(ionDef->Z(), ionDef->A());
319 mass += ionDef->NElectrons()*G4Electron::Definition()->GetPDGMass();
320 G4double charge = ionDef->Charge();
322 totalEnergyIn, kineticEnergyIn, momentumIn, ffact, ionDef, ionPDGID);
328 auto searchName = commonSubstitutions.find(particleName);
329 if (searchName != commonSubstitutions.end())
331 G4cout <<
"Substituting particle name \"" << particleName <<
"\" for the Geant4 name: \"" << searchName->second <<
"\"" << G4endl;
332 particleName = searchName->second;
336 G4ParticleDefinition* particleDef =
nullptr;
343 int particleID = std::stoi(particleName);
347 G4ParticleTable::G4PTblEncodingDictionary* encoding = G4ParticleTable::fEncodingDictionary;
348 auto search = encoding->find(particleID);
349 if (search != encoding->end())
350 {particleDef = search->second;}
352 {
throw BDSException(__METHOD_NAME__,
"PDG ID \"" + particleName +
"not found in particle table");}
354 catch (
const std::logic_error&)
355 {particleDef = particleTable->FindParticle(particleName);}
359 throw BDSException(__METHOD_NAME__,
"Particle \"" + particleName +
"\" not found.");
361 particleDefB =
new BDSParticleDefinition(particleDef, totalEnergyIn, kineticEnergyIn, momentumIn, ffact);
368 if (name ==
"proton")
369 {G4Proton::ProtonDefinition();}
370 else if (name ==
"antiproton")
371 {G4AntiProton::AntiProtonDefinition();}
372 else if (name ==
"e-")
373 {G4Electron::ElectronDefinition();}
374 else if (name ==
"e+")
375 {G4Positron::PositronDefinition();}
376 else if (name ==
"pi-")
377 {G4PionMinus::PionMinusDefinition();}
378 else if (name ==
"pi+")
379 {G4PionPlus::PionPlusDefinition();}
380 else if (name ==
"pi0")
381 {G4PionZero::PionZeroDefinition();}
382 else if (name ==
"neutron")
383 {G4Neutron::NeutronDefinition();}
384 else if (name ==
"photon" || name ==
"gamma")
386 else if (name ==
"mu-")
387 {G4MuonMinus::MuonMinusDefinition();}
388 else if (name ==
"mu+")
389 {G4MuonPlus::MuonPlusDefinition();}
390 else if (name ==
"kaon-")
391 {G4KaonMinus::KaonMinusDefinition();}
392 else if (name ==
"kaon+")
393 {G4KaonPlus::KaonPlusDefinition();}
394 else if (name ==
"kaon0L")
395 {G4KaonZeroLong::KaonZeroLongDefinition();}
396 else if (name ==
"nu_e")
397 {G4NeutrinoE::NeutrinoEDefinition();}
398 else if (name ==
"anti_nu_e")
399 {G4AntiNeutrinoE::AntiNeutrinoEDefinition();}
400 else if (name ==
"nu_mu")
401 {G4NeutrinoMu::NeutrinoMuDefinition();}
402 else if (name ==
"anti_nu_mu")
403 {G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();}
404 else if (name ==
"nu_tau")
405 {G4NeutrinoTau::NeutrinoTauDefinition();}
406 else if (name ==
"anti_nu_tau")
407 {G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();}
410 G4String msg =
"Unknown common particle type \"" + name;
411 msg +=
"\" - if it doesn't work, include all \"all_particles\" in the physicsList option.";
412 BDS::Warning(__METHOD_NAME__, msg);
419 G4Electron::ElectronDefinition();
420 G4Positron::PositronDefinition();
421 G4NeutrinoE::NeutrinoEDefinition();
422 G4AntiNeutrinoE::AntiNeutrinoEDefinition();
425 G4Proton::ProtonDefinition();
426 G4AntiProton::AntiProtonDefinition();
427 G4Neutron::NeutronDefinition();
428 G4AntiNeutron::AntiNeutronDefinition();
436 G4LeptonConstructor::ConstructParticle();
437 G4PionPlus::PionPlusDefinition();
438 G4PionMinus::PionMinusDefinition();
439 G4KaonPlus::KaonPlusDefinition();
440 G4KaonMinus::KaonMinusDefinition();
441 G4GenericIon::GenericIonDefinition();
446 G4cout <<
"Register physics processes by name for the primary particle \""
447 << primaryParticleName <<
"\":" << G4endl;
449 auto particle = G4ParticleTable::GetParticleTable()->FindParticle(primaryParticleName);
453 G4cout << __METHOD_NAME__ <<
"primary particle not defined yet - could be ion" << G4endl;
457 auto pl = particle->GetProcessManager()->GetProcessList();
458 for (G4int i = 0; i < (G4int)pl->length(); i++)
459 {G4cout <<
"\"" << (*pl)[i]->GetProcessName() <<
"\"" << G4endl;}
465 std::set<const G4ParticleDefinition*> particlesToBias;
466 for (
const auto& b : biases)
468 const G4ParticleDefinition* particle =
nullptr;
469 G4String particleName = G4String(b.particle);
471 particle = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
474 {particlesToBias.insert(particle);}
476 {
throw BDSException(__METHOD_NAME__,
"Unknown particle type for biasing: \"" + particleName +
"\"");}
479 if (particlesToBias.empty())
482 G4GenericBiasingPhysics* physBias =
new G4GenericBiasingPhysics();
483 for (
auto part : particlesToBias)
485 const G4String& particleName = part->GetParticleName();
486 G4cout << __METHOD_NAME__ <<
"wrapping \"" << particleName <<
"\" for biasing" << G4endl;
487 physBias->Bias(particleName);
494 G4cout << __METHOD_NAME__ <<
"Defined particles: " << G4endl;
495 auto it = G4ParticleTable::GetParticleTable()->GetIterator();
498 {G4cout << it->value()->GetParticleName() <<
" ";}
502#if G4VERSION_NUMBER > 1039
508 G4VModularPhysicsList* physlist =
new FTFP_BERT();
509 physlist->RegisterPhysics(
new G4IonElasticPhysics());
510 G4GenericBiasingPhysics* biasingPhysics =
new G4GenericBiasingPhysics();
519 physlist->ReplacePhysics(
new G4EmStandardPhysics_option4());
523 physlist->ReplacePhysics(
new G4EmStandardPhysicsSS());
529 G4cout <<
"Adding EM Dissocation to crystal channelling physics list" << G4endl;
533 biasingPhysics->PhysicsBiasAllCharged();
534 physlist->RegisterPhysics(biasingPhysics);
538 G4cout <<
"\nWARNING - adding cuts and limits physics process to \"COMPLETE\" physics list" << G4endl;
539 G4cout <<
"This is to enforce BDSIM range cuts and the minimumKinetic energy option.\n";
540 G4cout <<
"This can be turned off by setting option, g4PhysicsUseBDSIMCutsAndLimits=0;\n" << G4endl;
552 if (globals->DefaultRangeCutsSet())
554 G4cout << __METHOD_NAME__ <<
"Default production range cut " << physicsList->GetDefaultCutValue() <<
" mm" << G4endl;
555 physicsList->SetDefaultCutValue(globals->DefaultRangeCut());
557 if (globals->ProdCutPhotonsSet())
559 G4cout << __METHOD_NAME__ <<
"Photon production range cut " << physicsList->GetCutValue(
"gamma") <<
" mm" << G4endl;
560 physicsList->SetCutValue(globals->ProdCutPhotons(),
"gamma");
562 if (globals->ProdCutElectronsSet())
564 G4cout << __METHOD_NAME__ <<
"Electron production range cut " << physicsList->GetCutValue(
"e-") <<
" mm" << G4endl;
565 physicsList->SetCutValue(globals->ProdCutElectrons(),
"e-");
567 if (globals->ProdCutPositronsSet())
569 G4cout << __METHOD_NAME__ <<
"Positron production range cut " << physicsList->GetCutValue(
"e+") <<
" mm" << G4endl;
570 physicsList->SetCutValue(globals->ProdCutPositrons(),
"e+");
572 if (globals->ProdCutProtonsSet())
574 G4cout << __METHOD_NAME__ <<
"Proton production range cut " << physicsList->GetCutValue(
"proton") <<
" mm" << G4endl;
575 physicsList->SetCutValue(globals->ProdCutProtons(),
"proton");
579 G4cout << __METHOD_NAME__ <<
"Range cuts from inspection of the physics list" << G4endl;
580 G4cout << __METHOD_NAME__ <<
"Default production range cut " << physicsList->GetDefaultCutValue() <<
" mm" << G4endl;
583 G4cout << __METHOD_NAME__ <<
"List of all constructed particles by physics lists" << G4endl;
584 for (
auto particle : *G4ParticleTable::fDictionary)
585 {G4cout << particle.second->GetParticleName() <<
", ";}
589 physicsList->DumpCutValuesTable(verbosity);
595 G4double energyLimitLow = globals->PhysicsEnergyLimitLow();
596 G4double energyLimitHigh = globals->PhysicsEnergyLimitHigh();
600 if (setEnergyLimitLow || setEnergyLimitHigh)
602 auto table = G4ProductionCutsTable::GetProductionCutsTable();
603 G4double defaultEnergyLimitLow = table->GetLowEdgeEnergy();
604 G4double defaultEnergyLimitHigh = table->GetHighEdgeEnergy();
605 G4double elLow = setEnergyLimitLow ? energyLimitLow : defaultEnergyLimitLow;
606 G4double elHigh = setEnergyLimitHigh ? energyLimitHigh : defaultEnergyLimitHigh;
607 table->SetEnergyRange(elLow, elHigh);
608 if (setEnergyLimitLow)
610 G4cout << __METHOD_NAME__ <<
"set EM physics low energy limit: "
611 << elLow/CLHEP::MeV <<
" MeV" << G4endl;
613 if (setEnergyLimitHigh)
615 G4cout << __METHOD_NAME__ <<
"set high energy limit: "
616 << elHigh/CLHEP::TeV <<
" TeV" << G4endl;
617 if (elHigh > G4EmParameters::Instance()->MaxKinEnergy())
619 G4cout << __METHOD_NAME__ <<
"set EM physics Ek limit to " << elHigh/CLHEP::TeV <<
" TeV" << G4endl;
620 G4EmParameters::Instance()->SetMaxEnergy(elHigh);
622#if G4VERSION_NUMBER > 1069
624 if (elHigh > G4HadronicParameters::Instance()->GetMaxEnergy())
626 G4cout << __METHOD_NAME__ <<
"set hadronic physics Ek limit to " << elHigh/CLHEP::TeV <<
" TeV" << G4endl;
627 G4HadronicParameters::Instance()->SetMaxEnergy(elHigh);
634#if G4VERSION_NUMBER > 1049
649 G4double warningEnergy = 1.0 * CLHEP::kiloelectronvolt;
650 G4double importantEnergy = 10.0 * CLHEP::kiloelectronvolt;
651 G4double numberOfTrials = 1500;
653 auto transport = transportPair.first;
654 auto coupledTransport = transportPair.second;
659 transport->SetThresholdWarningEnergy(warningEnergy);
660 transport->SetThresholdImportantEnergy(importantEnergy);
661 transport->SetThresholdTrials(numberOfTrials);
663 else if(coupledTransport)
666 coupledTransport->SetThresholdWarningEnergy(warningEnergy);
667 coupledTransport->SetThresholdImportantEnergy(importantEnergy);
668 coupledTransport->SetThresholdTrials(numberOfTrials);
674 const auto* partPM = particleDef->GetProcessManager();
676 G4VProcess* partTransport = partPM->GetProcess(
"Transportation");
677 auto transport =
dynamic_cast<G4Transportation*
>(partTransport);
679 partTransport = partPM->GetProcess(
"CoupledTransportation");
680 auto coupledTransport =
dynamic_cast<G4CoupledTransportation*
>(partTransport);
682 return std::make_pair(transport, coupledTransport);
Copy of crystal channelling em physics from Geant4.
General exception with possible name of object and message.
A class that holds global options and constants.
static BDSGlobalConstants * Instance()
Access method.
Class to parse an ion particle definition.
Modular physics list based on Geant4 constructors.
Wrapper for particle definition.
G4ParticleDefinition * ParticleDefinition() const
Accessor.
Channelling physics process.
Physics processes required for user tracking limits.
Electromagnetic dissociation for high energy ions.
std::string particle
beam parameters
double beamKineticEnergy
beam parameters
double beamEnergy
beam parameters
double Ek0
initial beam centroid
std::string beamParticleName
beam parameters
double P0
initial beam centroid
double beamMomentum
beam parameters
double E0
initial beam centroid
double get_value(std::string name) const
get method (only for doubles)
bool HasBeenSet(const std::string &name) const
Whether a parameter has been set using the set_value method or not.
void SetRangeCuts(G4VModularPhysicsList *physicsList, G4int verbosity=1)
void ConstructBeamParticleG4(const G4String &name)
Ensure required beam particle has been constructed for Geant4 purposes.
void PrintDefinedParticles()
G4bool StrContains(const G4String &str, const G4String &test)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
std::pair< G4Transportation *, G4CoupledTransportation * > FindTransportation(const G4ParticleDefinition *particleDef)
Taken from Geant4 field01 example. Get the two possible transportation processes.
G4String GetFullPath(G4String filename, bool excludeNameFromPath=false, bool useCWDForPrefix=false)
BDSParticleDefinition * ConstructParticleDefinition(const G4String &particleNameIn, G4double totalEnergyIn, G4double kineticEnergyIn, G4double momentumIn, G4double ffact=1)
G4String LowerCase(const G4String &str)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
void ConflictingParametersSet(const GMAD::Beam &beamDefinition, const std::set< std::string > &keys, G4int nSet, G4bool warnZeroParamsSet=true, const G4String &unitString="")
Throw an exception if too few or too many parameters are set for the supplied keys.
void PrintPrimaryParticleProcesses(const G4String &primaryParticleName)
void CheckAndSetEnergyValidityRange()
void FixGeant105ThreshholdsForBeamParticle(const BDSParticleDefinition *particleDefinition)
Apply FixGeant105ThreshholdsForParticle to the beam particle definition.
void FixGeant105ThreshholdsForParticle(const G4ParticleDefinition *particleDefinition)
void ConstructExtendedParticleSet()
void ConstructDesignAndBeamParticle(const GMAD::Beam &beamDefinition, G4double ffact, BDSParticleDefinition *&designParticle, BDSParticleDefinition *&beamParticle, G4bool &beamDifferentFromDesignParticle)
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())
G4int NBeamParametersSet(const GMAD::Beam &beamDefinition, const std::set< std::string > &keys)
Count how many out of the set of keys in a beam instance are set.
G4VModularPhysicsList * BuildPhysics(const G4String &physicsList, G4int verbosity=1)
G4bool IsIon(const G4ParticleDefinition *particle)
Whether a particle is an ion. A proton is counted NOT as an ion.
void ConstructMinimumParticleSet()
G4VModularPhysicsList * ChannellingPhysicsComplete(G4bool useEMD=false, G4bool regular=false, G4bool em4=false, G4bool emss=false)
Build the physics required for channelling to work correctly.