ATLAS Offline Software
RHadronsPhysicsTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // class header
6 #include "RHadronsPhysicsTool.h"
7 // package headers
9 #include "FullModelHadronicProcess.hh"
10 #include "RHadronPythiaDecayer.h"
11 // Geant4 headers
12 #include "globals.hh"
13 #include "G4ParticleTable.hh"
14 #include "G4VProcess.hh"
15 #include "G4Transportation.hh"
16 #include "G4hMultipleScattering.hh"
17 #include "G4hIonisation.hh"
18 #include "G4ProcessManager.hh"
19 #include "G4Decay.hh"
20 #include "G4BaryonConstructor.hh"
21 
22 // STL headers
23 #include <memory>
24 #include <string>
25 
26 
27 //-----------------------------------------------------------------------------
28 // Implementation file for class : RHadronsPhysicsTool
29 //
30 // 2015-05-14 Edoardo Farina
31 //-----------------------------------------------------------------------------
32 
33 //=============================================================================
34 // Standard constructor, initializes variables
35 //=============================================================================
37  const std::string& nam,const IInterface* parent )
38  : base_class ( type, nam , parent )
39 {
40  m_physicsOptionType = G4AtlasPhysicsOption::Type::BSMPhysics;
41 }
42 
43 //=============================================================================
44 // Destructor
45 //=============================================================================
46 
48 {
49 }
50 
51 //=============================================================================
52 // Initialize
53 //=============================================================================
55 {
56  ATH_MSG_DEBUG("RHadronsPhysicsTool::initialize()");
57  return StatusCode::SUCCESS;
58 }
59 
60 auto RHadronsPhysicsTool::GetPhysicsOption() -> UPPhysicsConstructor {
61  return std::make_unique<RHadronsPhysicsTool::PhysicsConstructor>(
62  name(), this->msgLevel(), m_standardpdgidtodecay.value());
63 }
64 
66  ATH_MSG_DEBUG("RHadronsPhysicsTool::ConstructParticle() - start");
68  ATH_MSG_DEBUG("RHadronsPhysicsTool::ConstructParticle() - end");
69 }
71  ATH_MSG_DEBUG("RHadronsPhysicsTool::ConstructProcess() - start");
72  G4Decay* pythiaDecayProcess = new G4Decay();
73  pythiaDecayProcess->SetExtDecayer( new RHadronPythiaDecayer("RHadronPythiaDecayer") );
74  G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
75  particleIterator->reset();
76 
77  //First deal with the standard particles that G4 doesn't know about...
78  //G4Etac::Definition();
79  G4BaryonConstructor::ConstructParticle();
80  ATH_MSG_DEBUG("RHadronsPhysicsTool::ConstructProcess() - m_standardpdgidtodecay = " << m_standardpdgidtodecay);
81  G4ProcessManager *templateProcessMgr = G4ParticleTable::GetParticleTable()->FindParticle(4122)->GetProcessManager();
82  for (const int pid : m_standardpdgidtodecay) {
83  ATH_MSG_VERBOSE ( "Adding decay for "<<pid );
84  G4ParticleDefinition *particle = G4ParticleTable::GetParticleTable()->FindParticle( pid );
85  if (particle) {
86  ATH_MSG_VERBOSE ( particle->GetParticleName()<<" is standard for Pythia, lifetime is "<<particle->GetPDGLifeTime() );
87  G4ProcessManager *processMgr = particle->GetProcessManager();
88  if (!processMgr) {
89  ATH_MSG_VERBOSE ( "No process manager found for " << particle->GetParticleName() << " (" << pid << "). Copying process manager from 4122 (one we know works) to this particle" );
90  particle->SetProcessManager(new G4ProcessManager(*templateProcessMgr));
91  processMgr = particle->GetProcessManager();
92  }
93  // Look for native G4Decay process
94  G4ProcessVector *fullProcessList = processMgr->GetProcessList(); // NB G4ProcessVector does not support range-based for loops in G4 10.6....
95  std::vector< G4VProcess * > existingDecayProcesses; existingDecayProcesses.reserve(2);
96  for (unsigned int i=0; i<fullProcessList->size(); ++i) {
97  G4VProcess* process = (*fullProcessList)[i];
98  if (process->GetProcessType() == fDecay) {
99  ATH_MSG_VERBOSE ( "Found a pre-existing decay process for " <<particle->GetParticleName() << " (" << pid << "). Will remove in favour of using RHadronPythiaDecayer." );
100  existingDecayProcesses.push_back(process);
101  }
102  }
103  // Actually remove the existing decay processes
104  for (G4VProcess* process : existingDecayProcesses) {
105  processMgr->RemoveProcess(process);
106  }
107  // Cross-check
108  for (unsigned int i=0; i<fullProcessList->size(); ++i) {
109  G4VProcess* process = (*fullProcessList)[i];
110  if (process->GetProcessType() == fDecay) {
111  ATH_MSG_WARNING ( "There is another decay process for " <<particle->GetParticleName() << " (" << pid << ") already defined!" );
112  processMgr->DumpInfo();
113  }
114  }
115  ATH_MSG_VERBOSE ( "Adding decay process for " <<particle->GetParticleName() << " (" << pid << ") using RHadronPythiaDecayer." );
116  processMgr->AddProcess(pythiaDecayProcess);
117  processMgr->SetProcessOrdering(pythiaDecayProcess, idxPostStep); processMgr->SetProcessOrdering(pythiaDecayProcess, idxAtRest);
118  //processMgr->DumpInfo();
119  } else {
120  ATH_MSG_WARNING ( "Particle with pdgid "<<pid<<" has no definition in G4?" );
121  }
122  } // Loop over all particles that we need to define
123 
124  // Now add RHadrons... Keep a vector of those we've already dealt with
125  std::vector<int> handled;
126  // Use the G4 particle iterator
127  while ((*particleIterator)()) {
128  G4ParticleDefinition *particle = particleIterator->value();
130  const int pid = particle->GetPDGEncoding();
131  if (find(handled.begin(), handled.end(), pid) == handled.end()) {
132  handled.push_back(pid);
133  ATH_MSG_VERBOSE ( particle->GetParticleName() << " (" << particle->GetPDGEncoding() << ") " << " is a Custom Particle. Attempting to add a decay process." );
134  G4ProcessManager *processMgr = particle->GetProcessManager();
135  // Look for native G4Decay process
136  G4ProcessVector *fullProcessList = processMgr->GetProcessList(); // NB G4ProcessVector does not support range-based for loops in G4 10.6....
137  std::vector< G4VProcess * > existingDecayProcesses; existingDecayProcesses.reserve(2);
138  for (unsigned int i=0; i<fullProcessList->size(); ++i) {
139  G4VProcess* process = (*fullProcessList)[i];
140  if (process->GetProcessType() == fDecay) {
141  ATH_MSG_VERBOSE ( "Found a pre-existing decay process for " <<particle->GetParticleName() << " (" << pid << "). Will remove in favour of using RHadronPythiaDecayer." );
142  existingDecayProcesses.push_back(process);
143  }
144  }
145  // Actually remove the existing decay process(es)
146  for (G4VProcess* process : existingDecayProcesses) {
147  processMgr->RemoveProcess(process);
148  }
149  // Cross-check
150  for (unsigned int i=0; i<fullProcessList->size(); ++i) {
151  G4VProcess* process = (*fullProcessList)[i];
152  if (process->GetProcessType() == fDecay) {
153  ATH_MSG_WARNING ( "There is another decay process for " <<particle->GetParticleName() << " (" << pid << ") already defined!" );
154  processMgr->DumpInfo();
155  }
156  }
157  if (particle->GetParticleType()=="rhadron" ) {
158  processMgr->AddDiscreteProcess(new FullModelHadronicProcess());
159  if (pythiaDecayProcess->IsApplicable(*particle)) {
160  ATH_MSG_VERBOSE ( "Adding decay..." );
161  processMgr->AddProcess(pythiaDecayProcess);
162  // set ordering for PostStepDoIt and AtRestDoIt
163  processMgr->SetProcessOrdering(pythiaDecayProcess, idxPostStep);
164  processMgr->SetProcessOrdering(pythiaDecayProcess, idxAtRest);
165  } else {
166  ATH_MSG_WARNING ( "No decay allowed for " << particle->GetParticleName() );
167  if (!particle->GetPDGStable() && particle->GetPDGLifeTime()<0.1*CLHEP::ns) {
168  ATH_MSG_WARNING ( "Gonna decay it anyway!!!" );
169  processMgr->AddProcess(pythiaDecayProcess);
170  // set ordering for PostStepDoIt and AtRestDoIt
171  processMgr->SetProcessOrdering(pythiaDecayProcess, idxPostStep);
172  processMgr->SetProcessOrdering(pythiaDecayProcess, idxAtRest);
173  }
174  }
175  }
176  if (particle->GetPDGCharge()/CLHEP::eplus != 0) {
177  processMgr->AddProcess(new G4hMultipleScattering, -1, 1, 1);
178  processMgr->AddProcess(new G4hIonisation, -1, 2, 2);
179  ATH_MSG_VERBOSE ( "Processes for charged particle added." );
180  } else {
181  ATH_MSG_VERBOSE ( "It is neutral!!" );
182  }
183  processMgr->DumpInfo();
184  } else {
185  ATH_MSG_VERBOSE ( "Skipping already handled particle: "<<particle->GetParticleName() );
186  } // If it has not been handled yet
187  } // If this is one of our custom particles (RHadrons)
188  } // End of the particle iterator
189  ATH_MSG_DEBUG("RHadronsPhysicsTool::ConstructProcess() - list of handled RHadrons = " << handled);
190  ATH_MSG_DEBUG("RHadronsPhysicsTool::ConstructProcess() - end");
191 }
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:79
RHadronPythiaDecayer.h
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
CustomParticleFactory::loadCustomParticles
static void loadCustomParticles()
Definition: CustomParticleFactory.cxx:32
RHadronsPhysicsTool::initialize
virtual StatusCode initialize()
Initialize method.
Definition: RHadronsPhysicsTool.cxx:54
RHadronsPhysicsTool::PhysicsConstructor::ConstructParticle
virtual void ConstructParticle() override
Definition: RHadronsPhysicsTool.cxx:65
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
RHadronsPhysicsTool::RHadronsPhysicsTool
RHadronsPhysicsTool(const std::string &type, const std::string &name, const IInterface *parent)
Standard constructor.
Definition: RHadronsPhysicsTool.cxx:36
SUSY_SimplifiedModel_PostInclude.process
string process
Definition: SUSY_SimplifiedModel_PostInclude.py:43
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
python.SystemOfUnits.eplus
float eplus
Definition: SystemOfUnits.py:155
lumiFormat.i
int i
Definition: lumiFormat.py:85
RHadronsPhysicsTool.h
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
RHadronPythiaDecayer
Definition: RHadronPythiaDecayer.h:18
CustomParticleFactory.h
ParticleGun_EoverP_Config.pid
pid
Definition: ParticleGun_EoverP_Config.py:62
test_pyathena.parent
parent
Definition: test_pyathena.py:15
CustomParticleFactory::isCustomParticle
static bool isCustomParticle(G4ParticleDefinition *particle)
Definition: CustomParticleFactory.cxx:26
RHadronsPhysicsTool::PhysicsConstructor::ConstructProcess
virtual void ConstructProcess() override
Definition: RHadronsPhysicsTool.cxx:70
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
RHadronsPhysicsTool::~RHadronsPhysicsTool
virtual ~RHadronsPhysicsTool()
Destructor.
Definition: RHadronsPhysicsTool.cxx:47
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
RHadronsPhysicsTool::GetPhysicsOption
virtual UPPhysicsConstructor GetPhysicsOption()
Implements.
Definition: RHadronsPhysicsTool.cxx:60
G4AtlasPhysicsOption::BSMPhysics
@ BSMPhysics
Definition: IPhysicsOptionTool.h:23
RHadronsPhysicsTool::m_standardpdgidtodecay
IntegerArrayProperty m_standardpdgidtodecay
Definition: RHadronsPhysicsTool.h:89
python.SystemOfUnits.ns
float ns
Definition: SystemOfUnits.py:146