BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSPhysicsUtilities.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "G4Version.hh"
20
21#include "BDSDebug.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"
29#endif
30#include "BDSPhysicsCutsAndLimits.hh"
31#include "BDSPhysicsEMDissociation.hh"
32#include "BDSPhysicsMuonSplitting.hh"
33#include "BDSPhysicsUtilities.hh"
34#include "BDSUtilities.hh"
35#include "BDSWarning.hh"
36#include "BDSEmStandardPhysicsOp4Channelling.hh" // included with bdsim
37
38#include "FTFP_BERT.hh"
39#include "globals.hh"
40#include "G4AntiNeutrinoE.hh"
41#include "G4AntiNeutrinoMu.hh"
42#include "G4AntiNeutrinoTau.hh"
43#include "G4AntiNeutron.hh"
44#include "G4AntiProton.hh"
45#include "G4Electron.hh"
46#include "G4EmParameters.hh"
47#include "G4EmStandardPhysics_option4.hh"
48#include "G4EmStandardPhysicsSS.hh"
49#include "G4DynamicParticle.hh"
50#include "G4Gamma.hh"
51#include "G4GenericBiasingPhysics.hh"
52#include "G4GenericIon.hh"
53#include "G4IonElasticPhysics.hh"
54#include "G4IonTable.hh"
55#include "G4KaonMinus.hh"
56#include "G4KaonPlus.hh"
57#include "G4KaonZero.hh"
58#include "G4KaonZeroLong.hh"
59#include "G4KaonZeroShort.hh"
60#include "G4LeptonConstructor.hh"
61#include "G4MuonMinus.hh"
62#include "G4MuonPlus.hh"
63#include "G4NeutrinoE.hh"
64#include "G4NeutrinoMu.hh"
65#include "G4NeutrinoTau.hh"
66#include "G4Neutron.hh"
67#include "G4ParticleTable.hh"
68#include "G4ParticleTableIterator.hh"
69#include "G4PionMinus.hh"
70#include "G4PionPlus.hh"
71#include "G4PionZero.hh"
72#include "G4Positron.hh"
73#include "G4ProductionCutsTable.hh"
74#include "G4Proton.hh"
75#include "G4PhysListFactory.hh"
76#include "G4ProcessManager.hh"
77#include "G4ProcessVector.hh"
78#include "G4Proton.hh"
79#include "G4String.hh"
80#include "G4UImanager.hh"
81#if G4VERSION_NUMBER > 1049
82#include "G4ParticleDefinition.hh"
83#include "G4CoupledTransportation.hh"
84#include "G4Transportation.hh"
85#include <utility>
86#endif
87
88#if G4VERSION_NUMBER > 1069
89#include "G4HadronicParameters.hh"
90#endif
91
92#include "parser/beam.h"
93#include "parser/fastlist.h"
94#include "parser/physicsbiasing.h"
95
96#include <iomanip>
97#include <map>
98#include <regex>
99#include <set>
100#include <stdexcept>
101#include <string> // for stoi
102
103G4bool BDS::IsIon(const G4ParticleDefinition* particle)
104{
105 return G4IonTable::IsIon(particle) && particle!=G4Proton::Definition();
106}
107
108G4bool BDS::IsIon(const G4DynamicParticle* particle)
109{
110 return BDS::IsIon(particle->GetDefinition()) || particle->GetTotalOccupancy()>0;
111}
112
113G4VModularPhysicsList* BDS::BuildPhysics(const G4String& physicsList, G4int verbosity)
114{
115 // this must be done BEFORE any physics processes are constructed
116 // set the upper and lower energy levels applicable for all physics processes
118
119 G4VModularPhysicsList* result = nullptr;
120
122
124 G4String physicsListNameLower = BDS::LowerCase(physicsList);
125 G4bool useGeant4Physics = BDS::StrContains(physicsListNameLower, "g4");
126 G4bool completePhysics = BDS::StrContains(physicsListNameLower, "complete");
127 if (useGeant4Physics)
128 {
129 // strip off G4_ prefix - from original as G4 factory case sensitive
130 G4String geant4PhysicsList = physicsList.substr(2);
131 G4PhysListFactory factory;
132 if (!factory.IsReferencePhysList(geant4PhysicsList))
133 {
134 G4cerr << "Unknown Geant4 physics list \"" << geant4PhysicsList << "\"" << G4endl;
135 G4cout << "Available Geant4 hadronic physics lists:" << G4endl;
136 for (const auto &name : factory.AvailablePhysLists())
137 {G4cout << "\"" << name << "\"" << G4endl;}
138 G4cout << "Available Geant4 EM physics lists:" << G4endl;
139 for (const auto &name : factory.AvailablePhysListsEM())
140 {G4cout << "\"" << name << "\"" << G4endl;}
141 throw BDSException(__METHOD_NAME__, "Unknown Geant4 physics list \"" + geant4PhysicsList + "\"");
142 }
143 else
144 {
145 result = factory.GetReferencePhysList(geant4PhysicsList);
146 if (g->G4PhysicsUseBDSIMRangeCuts())
147 {BDS::SetRangeCuts(result);}
148 if (g->MinimumKineticEnergy() > 0 || g->G4PhysicsUseBDSIMCutsAndLimits())
149 {
150 G4cout << "\nAdding cuts and limits physics process to Geant4 reference physics list" << G4endl;
151 G4cout << "This is to enforce BDSIM range cuts and the minimumKinetic energy option.\n";
152 G4cout << "This is done by default for the functionality of BDSIM tracking and should not affect the physics greatly.\n";
153 G4cout << "See the BDSIM manual about Geant4 reference physics lists for details." << G4endl;
154 result->RegisterPhysics(new BDSPhysicsCutsAndLimits(g->ParticlesToExcludeFromCutsAsSet()));
155 }
156 else if (!g->G4PhysicsUseBDSIMCutsAndLimits() && g->Circular())
157 {
158 G4String message = "g4PhysicsUseBDSIMCutsAndLimits turned off but using a circular machine - circular mechanics will be broken";
159 BDS::Warning(__METHOD_NAME__, message);
160 }
161 }
162 }
163 else if (completePhysics)
164 {// we test one by one for the exact name of very specific physics lists
165 if (BDS::StrContains(physicsListNameLower, "channelling"))
166 {
167 G4cout << "Constructing \"" << physicsListNameLower << "\" complete physics list" << G4endl;
168#if G4VERSION_NUMBER > 1039
169 G4bool useEMD = BDS::StrContains(physicsListNameLower, "emd");
170 G4bool regular = BDS::StrContains(physicsListNameLower, "regular");
171 G4bool em4 = BDS::StrContains(physicsListNameLower, "em4");
172 G4bool emss = BDS::StrContains(physicsListNameLower, "emss");
173 // we don't assign 'result' variable or proceed as that would result in the
174 // range cuts being set for a complete physics list that we wouldn't use
175 auto r = BDS::ChannellingPhysicsComplete(useEMD, regular, em4, emss);
176 r->SetVerboseLevel(verbosity);
177 return r;
178#else
179 throw BDSException(__METHOD_NAME__, "Channel physics is not supported with Geant4 versions less than 10.4");
180#endif
181 }
182 else
183 {throw BDSException(__METHOD_NAME__, "Unknown 'complete' physics list \"" + physicsList + "\"");}
184 }
185 else
186 {
187 result = new BDSModularPhysicsList(physicsList);
188 BDS::SetRangeCuts(result, verbosity); // always set our range cuts for our physics list
189 }
190
191 // force construction of the particles - does no harm and helps with
192 // usage of exotic particle beams
193 result->ConstructParticle();
194 result->SetVerboseLevel(verbosity);
195
196 // apply any user-supplied macro to adjust geant4 physics parameters
197 G4UImanager* UIManager = G4UImanager::GetUIpointer();
198 G4String physicsMacro = BDSGlobalConstants::Instance()->Geant4PhysicsMacroFileName();
199 if (!physicsMacro.empty())
200 {
201 G4bool setInExecOptions = BDSGlobalConstants::Instance()->Geant4PhysicsMacroFileNameFromExecOptions();
202 G4String physicsMacroFull = BDS::GetFullPath(physicsMacro, false, setInExecOptions);
203 G4cout << "Applying geant4 physics macro file: " << physicsMacroFull << G4endl;
204 UIManager->ApplyCommand("/control/execute " + physicsMacroFull);
205 }
206
207 G4VUserPhysicsList* resultAsUserPhysicsList = dynamic_cast<G4VUserPhysicsList*>(result);
208 if (resultAsUserPhysicsList)
209 {// have to cast as they shadow functions and aren't virtual :(
210 resultAsUserPhysicsList->DumpCutValuesTable(verbosity);
211 resultAsUserPhysicsList->SetVerboseLevel(verbosity);
212 }
213 return result;
214}
215
216G4int BDS::NBeamParametersSet(const GMAD::Beam& beamDefinition,
217 const std::set<std::string>& keys)
218{
219 G4int nSet = 0;
220 for (const auto& k : keys)
221 {nSet += beamDefinition.HasBeenSet(k) ? 1 : 0;}
222 return nSet;
223}
224
225void BDS::ConflictingParametersSet(const GMAD::Beam& beamDefinition,
226 const std::set<std::string>& keys,
227 G4int nSet,
228 G4bool warnZeroParamsSet,
229 const G4String& unitString)
230{
231 if (nSet > 1)
232 {
233 G4cerr << "Beam> More than one parameter set - there should only be one" << G4endl;
234 for (const auto& k : keys)
235 {G4cerr << std::setw(14) << std::left << k << ": " << std::setw(7) << std::right << beamDefinition.get_value(k) << " " << unitString << G4endl;}
236 throw BDSException(__METHOD_NAME__, "conflicting parameters set");
237 }
238 else if (nSet == 0 && warnZeroParamsSet)
239 {
240 G4cerr << "Beam> One of the following required to be set" << G4endl;
241 for (const auto &k : keys)
242 {G4cerr << std::setw(14) << std::left << k << ": " << std::setw(7) << std::right << beamDefinition.get_value(k) << " " << unitString << G4endl;}
243 throw BDSException(__METHOD_NAME__, "insufficient parameters set");
244 }
245}
246
248 G4double ffact,
249 BDSParticleDefinition*& designParticle,
250 BDSParticleDefinition*& beamParticle,
251 G4bool& beamDifferentFromDesignParticle)
252{
253 if (beamDefinition.particle.empty())
254 {throw BDSException("Beam> no \"particle\" specified (required).");}
255
256 // check only one of the following has been set - ie no conflicting information
257 std::set<std::string> keysDesign = {"energy", "momentum", "kineticEnergy"};
258 G4int nSetDesign = BDS::NBeamParametersSet(beamDefinition, keysDesign);
259 BDS::ConflictingParametersSet(beamDefinition, keysDesign, nSetDesign);
260 std::set<std::string> keysParticle = {"E0", "P0", "Ek0"};
261 G4int nSetParticle = BDS::NBeamParametersSet(beamDefinition, keysParticle);
262
263 designParticle = BDS::ConstructParticleDefinition(beamDefinition.particle,
264 beamDefinition.beamEnergy * CLHEP::GeV,
265 beamDefinition.beamKineticEnergy * CLHEP::GeV,
266 beamDefinition.beamMomentum * CLHEP::GeV,
267 ffact);
268
269 // ensure a default
270 std::string beamParticleName = beamDefinition.beamParticleName.empty() ? beamDefinition.particle : beamDefinition.beamParticleName;
271 beamDifferentFromDesignParticle = nSetParticle > 0 || beamParticleName != beamDefinition.particle;
272 if (nSetParticle > 0)
273 {// at least one specified so use all of the beam particle ones
274 BDS::ConflictingParametersSet(beamDefinition, keysParticle, nSetParticle, false);
275 beamParticle = BDS::ConstructParticleDefinition(beamParticleName,
276 beamDefinition.E0 * CLHEP::GeV,
277 beamDefinition.Ek0 * CLHEP::GeV,
278 beamDefinition.P0 * CLHEP::GeV,
279 ffact);
280 }
281 else if (beamDefinition.beamParticleName != beamDefinition.particle && !beamDefinition.beamParticleName.empty())
282 {// different particle name but no extra E0, Ek0 or P0 -> therefore use general beam defaults
283 beamParticle = BDS::ConstructParticleDefinition(beamParticleName,
284 beamDefinition.beamEnergy * CLHEP::GeV,
285 beamDefinition.beamKineticEnergy * CLHEP::GeV,
286 beamDefinition.beamMomentum * CLHEP::GeV,
287 ffact);
288 }
289 else
290 {
291 beamParticle = new BDSParticleDefinition(*designParticle);
292 }
293}
294
296 G4double totalEnergyIn,
297 G4double kineticEnergyIn,
298 G4double momentumIn,
299 G4double ffact)
300{
301 BDSParticleDefinition* particleDefB = nullptr; // result
302 G4String particleName = BDS::LowerCase(particleNameIn);
303
304 std::map<G4String, G4String> commonSubstitutions = { {"photon", "gamma"},
305 {"electron", "e-"},
306 {"positron", "e+"},
307 {"pion+", "pi+"},
308 {"pion-", "pi-"},
309 {"pion0", "pi0"},
310 {"antiproton", "anti_proton"},
311 {"anti-proton", "anti_proton"} };
312
313 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
314 std::regex ionParticle("(ion\\s).*");
315 if (std::regex_match(particleName, ionParticle))
316 {
317 G4GenericIon::GenericIonDefinition(); // construct general ion particle
318 auto ionDef = new BDSIonDefinition(particleName); // parse the ion definition
319
320 G4IonTable* ionTable = particleTable->GetIonTable();
322 G4int ionPDGID = G4IonTable::GetNucleusEncoding(ionDef->Z(), ionDef->A());
323 G4double mass = ionTable->GetIonMass(ionDef->Z(), ionDef->A());
324 mass += ionDef->NElectrons()*G4Electron::Definition()->GetPDGMass();
325 G4double charge = ionDef->Charge(); // correct even if overridden
326 particleDefB = new BDSParticleDefinition(particleName, mass, charge,
327 totalEnergyIn, kineticEnergyIn, momentumIn, ffact, ionDef, ionPDGID);
328 delete ionDef;
329 }
330 else
331 {
332 // swap out some common name substitutions for Geant4 ones
333 auto searchName = commonSubstitutions.find(particleName);
334 if (searchName != commonSubstitutions.end())
335 {
336 G4cout << "Substituting particle name \"" << particleName << "\" for the Geant4 name: \"" << searchName->second << "\"" << G4endl;
337 particleName = searchName->second;
338 }
339
340 BDS::ConstructBeamParticleG4(particleName); // enforce construction of some basic particles
341 G4ParticleDefinition* particleDef = nullptr;
342
343 // try and see if it's an integer and therefore PDG ID, if not search by string
344 try
345 {
346 // we try this because std::stoi can throw a std::invalid_argument or
347 // std::out_of_range exception, both of which inherit std::logic_error
348 int particleID = std::stoi(particleName);
349 // we don't use the G4ParticleTable->FindParticle(int) because it unnecessarily
350 // checks for physics readiness and throws an exception. here we just inspect
351 // the encoding dictionary ourselves. it's all typedeffed but it's std::map<G4int, G4ParticleDefinition*>
352 G4ParticleTable::G4PTblEncodingDictionary* encoding = G4ParticleTable::fEncodingDictionary;
353 auto search = encoding->find(particleID);
354 if (search != encoding->end())
355 {particleDef = search->second;}
356 else
357 {throw BDSException(__METHOD_NAME__,"PDG ID \"" + particleName + "not found in particle table");}
358 }
359 catch (const std::logic_error&) // else, usual way by string search
360 {
361 particleDef = particleTable->FindParticle(particleName);
362 if (!particleDef)
363 {particleDef = particleTable->FindParticle(particleNameIn);} // try with original case
364 }
365 if (!particleDef)
366 {
368 throw BDSException(__METHOD_NAME__, "Particle \"" + particleName + "\" not found.");
369 }
370 particleDefB = new BDSParticleDefinition(particleDef, totalEnergyIn, kineticEnergyIn, momentumIn, ffact);
371 }
372 return particleDefB;
373}
374
375void BDS::ConstructBeamParticleG4(const G4String& name)
376{
377 // note, we compare to the lower case name here, not the Geant4 one as
378 // we will already have converted to lower case when this is used
379 if (name == "proton")
380 {G4Proton::ProtonDefinition();}
381 else if (name == "anti_proton")
382 {G4AntiProton::AntiProtonDefinition();}
383 else if (name == "e-")
384 {G4Electron::ElectronDefinition();}
385 else if (name == "e+")
386 {G4Positron::PositronDefinition();}
387 else if (name == "pi-")
388 {G4PionMinus::PionMinusDefinition();}
389 else if (name == "pi+")
390 {G4PionPlus::PionPlusDefinition();}
391 else if (name == "pi0")
392 {G4PionZero::PionZeroDefinition();}
393 else if (name == "neutron")
394 {G4Neutron::NeutronDefinition();}
395 else if (name == "photon" || name == "gamma")
396 {G4Gamma::Gamma();}
397 else if (name == "mu-")
398 {G4MuonMinus::MuonMinusDefinition();}
399 else if (name == "mu+")
400 {G4MuonPlus::MuonPlusDefinition();}
401 else if (name == "kaon-")
402 {G4KaonMinus::KaonMinusDefinition();}
403 else if (name == "kaon+")
404 {G4KaonPlus::KaonPlusDefinition();}
405 else if (name == "kaon0")
406 {G4KaonZero::KaonZeroDefinition();}
407 else if (name == "kaon0l")
408 {G4KaonZeroLong::KaonZeroLongDefinition();}
409 else if (name == "kaon0s")
410 {G4KaonZeroShort::KaonZeroShortDefinition();}
411 else if (name == "nu_e")
412 {G4NeutrinoE::NeutrinoEDefinition();}
413 else if (name == "anti_nu_e")
414 {G4AntiNeutrinoE::AntiNeutrinoEDefinition();}
415 else if (name == "nu_mu")
416 {G4NeutrinoMu::NeutrinoMuDefinition();}
417 else if (name == "anti_nu_mu")
418 {G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();}
419 else if (name == "nu_tau")
420 {G4NeutrinoTau::NeutrinoTauDefinition();}
421 else if (name == "anti_nu_tau")
422 {G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();}
423 else
424 {
425 G4String msg = "Unknown common beam particle type \"" + name;
426 msg += "\" - if it doesn't work, include all \"all_particles\" in the physicsList option.";
427 BDS::Warning(__METHOD_NAME__, msg);
428 }
429}
430
432{
433 // e-, e+, v_e, v_e(bar)
434 G4Electron::ElectronDefinition();
435 G4Positron::PositronDefinition();
436 G4NeutrinoE::NeutrinoEDefinition();
437 G4AntiNeutrinoE::AntiNeutrinoEDefinition();
438
439 // p, pbar, neutron, anti-neutron
440 G4Proton::ProtonDefinition();
441 G4AntiProton::AntiProtonDefinition();
442 G4Neutron::NeutronDefinition();
443 G4AntiNeutron::AntiNeutronDefinition();
444
445 // photon
446 G4Gamma::Gamma();
447}
448
450{
451 G4LeptonConstructor::ConstructParticle();
452 G4PionPlus::PionPlusDefinition();
453 G4PionMinus::PionMinusDefinition();
454 G4KaonPlus::KaonPlusDefinition();
455 G4KaonMinus::KaonMinusDefinition();
456 G4GenericIon::GenericIonDefinition();
457}
458
459void BDS::PrintPrimaryParticleProcesses(const G4String& primaryParticleName)
460{
461 G4cout << "Register physics processes by name for the primary particle \""
462 << primaryParticleName << "\":" << G4endl;
463
464 auto particle = G4ParticleTable::GetParticleTable()->FindParticle(primaryParticleName);
465 if (!particle)
466 {// could be ion that isn't defined
467#ifdef BDSDEBUG
468 G4cout << __METHOD_NAME__ << "primary particle not defined yet - could be ion" << G4endl;
469#endif
470 return;
471 }
472 auto pl = particle->GetProcessManager()->GetProcessList();
473 for (G4int i = 0; i < (G4int)pl->length(); i++)
474 {G4cout << "\"" << (*pl)[i]->GetProcessName() << "\"" << G4endl;}
475}
476
477
478G4GenericBiasingPhysics* BDS::BuildAndAttachBiasWrapper(const GMAD::FastList<GMAD::PhysicsBiasing>& biases)
479{
480 std::set<const G4ParticleDefinition*> particlesToBias;
481 for (const auto& b : biases)
482 {
483 const G4ParticleDefinition* particle = nullptr;
484 G4String particleName = G4String(b.particle);
485 // this works for "GenericIon" too
486 particle = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
487
488 if (particle)
489 {particlesToBias.insert(particle);}
490 else
491 {throw BDSException(__METHOD_NAME__, "Unknown particle type for biasing: \"" + particleName + "\"");}
492 }
493
494 if (particlesToBias.empty()) // nothing valid to bias
495 {return nullptr;}
496
497 G4GenericBiasingPhysics* physBias = new G4GenericBiasingPhysics();
498 for (auto part : particlesToBias)
499 {
500 const G4String& particleName = part->GetParticleName();
501 G4cout << __METHOD_NAME__ << "wrapping \"" << particleName << "\" for biasing" << G4endl;
502 physBias->Bias(particleName);
503 }
504 return physBias;
505}
506
507void BDS::BuildMuonBiasing(G4VModularPhysicsList* physicsList)
508{
509 auto globals = BDSGlobalConstants::Instance();
510 G4int muonSplittingFactor = globals->MuonSplittingFactor();
511 if (muonSplittingFactor > 1)
512 {
513 G4int muonSplittingFactor2 = globals->MuonSplittingFactor2();
514 G4double muonSplittingThresholdParentEk = globals->MuonSplittingThresholdParentEk();
515 G4double muonSplittingThresholdParentEk2 = globals->MuonSplittingThresholdParentEk2();
516 G4cout << "BDSPhysicsMuonSplitting -> using muon splitting wrapper -> factor of: " << muonSplittingFactor << G4endl;
517 if (muonSplittingThresholdParentEk > 0)
518 {G4cout << "BDSPhysicsMuonSplitting -> minimum parent kinetic energy: " << muonSplittingThresholdParentEk / CLHEP::GeV << " GeV" << G4endl;}
519 if (muonSplittingFactor2 > 1)
520 {
521 G4cout << "BDSPhysicsMuonSplitting -> factor #2: " << muonSplittingFactor2 << " for muons above "
522 << muonSplittingThresholdParentEk / CLHEP::GeV << " GeV" << G4endl;
523 }
524 G4bool excludeW1P = globals->MuonSplittingExcludeWeight1Particles();
525 physicsList->RegisterPhysics(new BDSPhysicsMuonSplitting(muonSplittingFactor, muonSplittingThresholdParentEk,
526 muonSplittingFactor2, muonSplittingThresholdParentEk2,
527 excludeW1P, globals->MuonSplittingExclusionWeight()));
528 }
529}
530
532{
533 G4cout << __METHOD_NAME__ << "Defined particles: " << G4endl;
534 auto it = G4ParticleTable::GetParticleTable()->GetIterator();
535 it->reset(); // because there's only 1 iterator due to geant4 design
536 while ((*it)())
537 {G4cout << it->value()->GetParticleName() << " ";}
538 G4cout << G4endl;
539}
540
541#if G4VERSION_NUMBER > 1039
542G4VModularPhysicsList* BDS::ChannellingPhysicsComplete(G4bool useEMD,
543 G4bool regular,
544 G4bool em4,
545 G4bool emss)
546{
547 G4VModularPhysicsList* physlist = new FTFP_BERT();
548 physlist->RegisterPhysics(new G4IonElasticPhysics()); // not included by default in FTFP_BERT
549 G4GenericBiasingPhysics* biasingPhysics = new G4GenericBiasingPhysics();
550 physlist->RegisterPhysics(new BDSPhysicsChannelling());
551 if (!regular)
552 {
553 // replace the EM physics in the Geant4 provided FTFP_BERT composite physics list
554 physlist->ReplacePhysics(new BDSEmStandardPhysicsOp4Channelling());
555 }
556 else if (em4)
557 {
558 physlist->ReplacePhysics(new G4EmStandardPhysics_option4());
559 }
560 else if (emss)
561 {
562 physlist->ReplacePhysics(new G4EmStandardPhysicsSS());
563 }
564
565 // optional electromagnetic dissociation that isn't in FTFP_BERT by default
566 if (useEMD)
567 {
568 G4cout << "Adding EM Dissocation to crystal channelling physics list" << G4endl;
569 physlist->RegisterPhysics(new BDSPhysicsEMDissociation());
570 }
571
572 biasingPhysics->PhysicsBiasAllCharged();
573 physlist->RegisterPhysics(biasingPhysics);
574 if (BDSGlobalConstants::Instance()->MinimumKineticEnergy() > 0 &&
575 BDSGlobalConstants::Instance()->G4PhysicsUseBDSIMCutsAndLimits())
576 {
577 G4cout << "\nWARNING - adding cuts and limits physics process to \"COMPLETE\" physics list" << G4endl;
578 G4cout << "This is to enforce BDSIM range cuts and the minimumKinetic energy option.\n";
579 G4cout << "This can be turned off by setting option, g4PhysicsUseBDSIMCutsAndLimits=0;\n" << G4endl;
580 physlist->RegisterPhysics(new BDSPhysicsCutsAndLimits(BDSGlobalConstants::Instance()->ParticlesToExcludeFromCutsAsSet()));
581 }
582 return physlist;
583}
584#endif
585
586void BDS::SetRangeCuts(G4VModularPhysicsList* physicsList, G4int verbosity)
587{
589
590 // overwrite when explicitly set in options
591 if (globals->DefaultRangeCutsSet())
592 {
593 G4cout << __METHOD_NAME__ << "Default production range cut " << physicsList->GetDefaultCutValue() << " mm" << G4endl;
594 physicsList->SetDefaultCutValue(globals->DefaultRangeCut());
595 }
596 if (globals->ProdCutPhotonsSet())
597 {
598 G4cout << __METHOD_NAME__ << "Photon production range cut " << physicsList->GetCutValue("gamma") << " mm" << G4endl;
599 physicsList->SetCutValue(globals->ProdCutPhotons(), "gamma");
600 }
601 if (globals->ProdCutElectronsSet())
602 {
603 G4cout << __METHOD_NAME__ << "Electron production range cut " << physicsList->GetCutValue("e-") << " mm" << G4endl;
604 physicsList->SetCutValue(globals->ProdCutElectrons(),"e-");
605 }
606 if (globals->ProdCutPositronsSet())
607 {
608 G4cout << __METHOD_NAME__ << "Positron production range cut " << physicsList->GetCutValue("e+") << " mm" << G4endl;
609 physicsList->SetCutValue(globals->ProdCutPositrons(),"e+");
610 }
611 if (globals->ProdCutProtonsSet())
612 {
613 G4cout << __METHOD_NAME__ << "Proton production range cut " << physicsList->GetCutValue("proton") << " mm" << G4endl;
614 physicsList->SetCutValue(globals->ProdCutProtons(), "proton");
615 }
616
617 // inspection of the physics list doesn't work at this point and seems to always return 0 apart from the default
618 G4cout << __METHOD_NAME__ << "Range cuts from inspection of the physics list" << G4endl;
619 G4cout << __METHOD_NAME__ << "Default production range cut " << physicsList->GetDefaultCutValue() << " mm" << G4endl;
620
621#ifdef BDSDEBUG
622 G4cout << __METHOD_NAME__ << "List of all constructed particles by physics lists" << G4endl;
623 for (auto particle : *G4ParticleTable::fDictionary)
624 {G4cout << particle.second->GetParticleName() << ", ";}
625 G4cout << G4endl;
626#endif
627
628 physicsList->DumpCutValuesTable(verbosity);
629}
630
632{
634 G4double energyLimitLow = globals->PhysicsEnergyLimitLow();
635 G4double energyLimitHigh = globals->PhysicsEnergyLimitHigh();
636 G4bool setEnergyLimitLow = BDS::IsFinite(energyLimitLow);
637 G4bool setEnergyLimitHigh = BDS::IsFinite(energyLimitHigh);
638
639 if (setEnergyLimitLow || setEnergyLimitHigh)
640 {
641 auto table = G4ProductionCutsTable::GetProductionCutsTable();
642 G4double defaultEnergyLimitLow = table->GetLowEdgeEnergy();
643 G4double defaultEnergyLimitHigh = table->GetHighEdgeEnergy();
644 G4double elLow = setEnergyLimitLow ? energyLimitLow : defaultEnergyLimitLow;
645 G4double elHigh = setEnergyLimitHigh ? energyLimitHigh : defaultEnergyLimitHigh;
646 table->SetEnergyRange(elLow, elHigh);
647 if (setEnergyLimitLow)
648 {
649 G4cout << __METHOD_NAME__ << "set EM physics low energy limit: "
650 << elLow/CLHEP::MeV << " MeV" << G4endl;
651 }
652 if (setEnergyLimitHigh)
653 {
654 G4cout << __METHOD_NAME__ << "set high energy limit: "
655 << elHigh/CLHEP::TeV << " TeV" << G4endl;
656 if (elHigh > G4EmParameters::Instance()->MaxKinEnergy())
657 {
658 G4cout << __METHOD_NAME__ << "set EM physics Ek limit to " << elHigh/CLHEP::TeV << " TeV" << G4endl;
659 G4EmParameters::Instance()->SetMaxEnergy(elHigh);
660 }
661#if G4VERSION_NUMBER > 1069
662 // this was in a patch of our own before 10.7 and compensates for ion physics
663 if (elHigh > G4HadronicParameters::Instance()->GetMaxEnergy())
664 {
665 G4cout << __METHOD_NAME__ << "set hadronic physics Ek limit to " << elHigh/CLHEP::TeV << " TeV" << G4endl;
666 G4HadronicParameters::Instance()->SetMaxEnergy(elHigh);
667 }
668#endif
669 }
670 }
671}
672
673#if G4VERSION_NUMBER > 1049
675{
676 G4ParticleDefinition* particleDef = particleDefinition->ParticleDefinition();
678}
679
680void BDS::FixGeant105ThreshholdsForParticle(const G4ParticleDefinition* particleDef)
681{
682 // in the case of ions the particle definition isn't available early on so protect
683 // against this
684 if (!particleDef)
685 {return;}
686 // taken from the Geant4.10.5 field01 example
687 // used to compensate for aggressive killing in Geant4.10.5
688 G4double warningEnergy = 1.0 * CLHEP::kiloelectronvolt; // Arbitrary
689 G4double importantEnergy = 10.0 * CLHEP::kiloelectronvolt; // Arbitrary
690 G4double numberOfTrials = 1500; // Arbitrary
691 auto transportPair = BDS::FindTransportation(particleDef);
692 auto transport = transportPair.first;
693 auto coupledTransport = transportPair.second;
694
695 if (transport)
696 {
697 // Change the values of the looping particle parameters of Transportation
698 transport->SetThresholdWarningEnergy(warningEnergy);
699 transport->SetThresholdImportantEnergy(importantEnergy);
700 transport->SetThresholdTrials(numberOfTrials);
701 }
702 else if(coupledTransport)
703 {
704 // Change the values for Coupled Transport
705 coupledTransport->SetThresholdWarningEnergy(warningEnergy);
706 coupledTransport->SetThresholdImportantEnergy(importantEnergy);
707 coupledTransport->SetThresholdTrials(numberOfTrials);
708 }
709}
710
711std::pair<G4Transportation*, G4CoupledTransportation*> BDS::FindTransportation(const G4ParticleDefinition* particleDef)
712{
713 const auto* partPM = particleDef->GetProcessManager();
714
715 G4VProcess* partTransport = partPM->GetProcess("Transportation");
716 auto transport = dynamic_cast<G4Transportation*>(partTransport);
717
718 partTransport = partPM->GetProcess("CoupledTransportation");
719 auto coupledTransport = dynamic_cast<G4CoupledTransportation*>(partTransport);
720
721 return std::make_pair(transport, coupledTransport);
722}
723#endif
Copy of crystal channelling em physics from Geant4.
General exception with possible name of object and message.
Definition: BDSException.hh:35
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.
High energy muon processes.
std::string particle
beam parameters
Definition: beamBase.h:40
double beamKineticEnergy
beam parameters
Definition: beamBase.h:43
double beamEnergy
beam parameters
Definition: beamBase.h:42
double Ek0
initial beam centroid
Definition: beamBase.h:74
std::string beamParticleName
beam parameters
Definition: beamBase.h:41
double P0
initial beam centroid
Definition: beamBase.h:75
double beamMomentum
beam parameters
Definition: beamBase.h:44
double E0
initial beam centroid
Definition: beamBase.h:73
Beam class.
Definition: beam.h:44
double get_value(std::string name) const
get method (only for doubles)
Definition: beam.cc:39
bool HasBeenSet(const std::string &name) const
Whether a parameter has been set using the set_value method or not.
Definition: beam.cc:142
List with Efficient Lookup.
Definition: fastlist.h:42
void SetRangeCuts(G4VModularPhysicsList *physicsList, G4int verbosity=1)
void ConstructBeamParticleG4(const G4String &name)
void PrintDefinedParticles()
G4bool StrContains(const G4String &str, const G4String &test)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
Definition: BDSUtilities.cc:66
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 BuildMuonBiasing(G4VModularPhysicsList *physicsList)
Build muon splitting biasing and wrap the various processes in the physics list.
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.