Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
CustomParticleFactory.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 
6 
8 #include "G4ParticleDefinition.hh"
9 #include "CustomParticle.h"
10 #include "CustomParticleFactory.h"
11 
13 
14 #include "G4DecayTable.hh"
15 #include "G4ParticleTable.hh"
16 #include "G4PhaseSpaceDecayChannel.hh"
17 
18 #include <fstream>
19 #include <iomanip>
20 #include <iostream>
21 #include <stdexcept>
22 #include <sstream>
23 #include <iterator>
24 #include <map>
25 
27 {
28  static const std::set<G4ParticleDefinition *> particles = load();
29  return (particle!=nullptr && particles.find(particle)!=particles.end());
30 }
31 
33 {
34  // tickle the loading of the particles if it wasn't done yet
35  isCustomParticle(nullptr);
36 }
37 
38 std::set<G4ParticleDefinition *> CustomParticleFactory::load()
39 {
40  // Reading decays from file
41  std::map< std::string, std::vector<std::vector<std::string>* > > decayMap;
42  std::ifstream decayFile("decays.txt");
43  std::string line;
44  while (getline(decayFile,line)) {
45  std::istringstream is(line);
46  std::vector<std::string>* txtvec = new std::vector<std::string>(std::istream_iterator<std::string>(is),std::istream_iterator<std::string>());
47  const std::string name = (*txtvec)[0];
48  try {
49  decayMap.at(name).push_back(txtvec);
50  }
51  catch (const::std::out_of_range& e) {
52  std::vector<std::vector<std::string>* > decays;
53  decays.push_back(txtvec);
54  decayMap[name] = std::move(decays);
55  }
56  }
57  decayFile.close();
58 
59  // the existing particle table
60  G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
61 
62  std::set<G4ParticleDefinition *> particles;
63 
64  std::ifstream configFile("particles.txt");
65  G4String pType="custom";
66  G4String pSubType="";
67  G4double spectatormass = 0;
68  G4ParticleDefinition* spectator = 0;
69  bool first = true;
70  double mass;
71  int pdgCode;
72  bool stable{true};
73  G4double lifetime{-1.0};
74  std::string name;
76 
77  // This should be compatible IMO to SLHA
78  while (getline(configFile,line)) {
79  G4cout << "-------------------------------------------------" << G4endl;
80  std::string::size_type beg_idx,end_idx;
81 
82  beg_idx = line.find_first_not_of("\t #");
83  if (beg_idx > 0 && line[beg_idx-1] == '#') continue;
84  end_idx = line.find_first_of( "\t ", beg_idx);
85  if (end_idx == std::string::npos) continue;
86  char *endptr;
87  const std::string pdgCodeString = line.substr( beg_idx, end_idx - beg_idx );
88  pdgCode = strtol(pdgCodeString.c_str(), &endptr, 0);
89  if (endptr[0] != '\0') {
90  throw std::invalid_argument("Could not convert string to int: " + line.substr( beg_idx, end_idx - beg_idx ));
91  }
92 
93  beg_idx = line.find_first_not_of("\t ",end_idx);
94  end_idx = line.find_first_of( "\t #", beg_idx);
95  if (end_idx == std::string::npos) continue;
96  mass = atof(line.substr( beg_idx, end_idx - beg_idx ).c_str());
97 
98  beg_idx = line.find_first_not_of("\t# ",end_idx);
99  end_idx = line.length();
100  if (beg_idx==std::string::npos) beg_idx = end_idx;
101  name = line.substr( beg_idx, end_idx - beg_idx );
102  while (name.c_str()[0] == ' ') { name.erase(0,1); }
103  while (name[name.size()-1] == ' ') { name.erase(name.size()-1,1); }
104 
105  if (abs(pdgCode) / 1000000 == 0){
106  G4cout << "Pdg code too low " << pdgCode << " "<<abs(pdgCode) / 1000000 << std::endl;
107  continue;
108  }
109 
110  stable = !(parameters->DoDecays() == 1 || (decayMap.contains(name) && name.find("cloud")>name.size() ));
111  lifetime = stable ? -1.0 : parameters->Lifetime();
112 
113  pType="custom";
114  if (MC::isRHadron(pdgCode)) pType = "rhadron";
115  if (MC::isSlepton(pdgCode)) pType = "sLepton";
116 
117  G4cout<<"pType: "<<pType<<G4endl;
118  G4cout<<"Charge of "<<name<<" is "<<MC::charge(pdgCode)<<G4endl;
119 
120  G4ParticleDefinition* previousDefinition = theParticleTable->FindParticle(pdgCode);
121  // if the particle has somehow already been added to the G4ParticleTable remove it
122  if (previousDefinition) {
123  G4cout << "Found an existing G4ParticleDefinition for " << name << "(" << pdgCode << ") in G4ParticleTable. Assuming that we want the new definition, so removing previous definition! May cause issues elsewhere though..." << G4endl;
124  previousDefinition->DumpTable();
125  theParticleTable->Remove(previousDefinition);
126  delete previousDefinition;
127  }
128  // Arguments for constructor are as follows
129  // name mass width charge
130  // 2*spin parity C-conjugation
131  // 2*Isospin 2*Isospin3 G-parity
132  // type lepton number baryon number PDG encoding
133  // stable lifetime decay table
134  // shortlived subType anti_encoding
136  name, mass * CLHEP::GeV , 0.0*CLHEP::MeV, CLHEP::eplus * MC::charge(pdgCode),
137  MC::spin2(pdgCode), +1, 0,
138  0, 0, 0,
139  pType, 0, +1, pdgCode,
140  stable, lifetime, nullptr );
141  if (pType=="custom") {
142  spectatormass = mass;
143  spectator = particle;
144  }
145  if (first) {
146  first = false;
147  spectatormass = mass;
148  spectator = particle;
149  } else {
150  G4String cloudname = name+"cloud";
151  G4String cloudtype = pType+"cloud";
152  G4double cloudmass = mass-spectatormass;
153  const bool cloudstable = !(parameters->DoDecays() == 1 || decayMap.contains(cloudname));
154  const G4double cloudlifetime = cloudstable ? -1.0 : parameters->Lifetime();
155 
156  std::unique_ptr<G4ParticleDefinition> tmpParticle = std::unique_ptr<CustomParticle>( new CustomParticle(
157  cloudname, cloudmass * CLHEP::GeV , 0.0*CLHEP::MeV, 0 ,
158  0, +1, 0,
159  0, 0, 0,
160  cloudtype, 0, +1, 0,
161  cloudstable, cloudlifetime, nullptr ) );
162  particle->SetCloud(tmpParticle);
163  particle->SetSpectator(spectator);
164 
165  G4cout << name << " being assigned cloud " << particle->GetCloud()->GetParticleName() << " and spectator " << particle->GetSpectator()->GetParticleName() << G4endl;
166  G4cout << "Masses: "
167  << particle->GetPDGMass()/CLHEP::GeV << " Gev, "
168  << particle->GetCloud()->GetPDGMass()/CLHEP::GeV << " GeV and "
169  << particle->GetSpectator()->GetPDGMass()/CLHEP::GeV << " GeV."
170  << G4endl;
171  }
172 
173  particles.insert(particle);
174  }
175  configFile.close();
176  G4cout << "-------------------------------------------------" << G4endl;
177 
178  // Looping over custom particles to add decays
179  for (G4ParticleDefinition *part : particles) {
180  name=part->GetParticleName();
181  if ( decayMap.contains(name) ) {
182  std::vector<std::vector<std::string>* >& mydecays = decayMap.at(name);
183  // Did I get any decays?
184  if (mydecays.empty()) { continue; } // Should not happen.
185  int ndec=mydecays.size();
186  G4DecayTable* table = new G4DecayTable();
187  G4VDecayChannel** mode = new G4VDecayChannel*[ndec]{};
188  for (int i=0;i!=ndec;i++) {
189  std::vector<std::string> thisdec=*(mydecays[i]);//deliberate copy
190  if (thisdec.empty()) continue;
191  const double branch = std::stod(thisdec.back()); // Reading branching ratio
192  thisdec.pop_back(); // Removing the number from the vector
193  for (const auto & thisString:thisdec) G4cout<<thisString<<G4endl;
194  if (thisdec.size()==3) {
195  mode[i] = new G4PhaseSpaceDecayChannel(thisdec[0],branch,2,thisdec[1],thisdec[2]);
196  } else if (thisdec.size()==4) {
197  mode[i] = new G4PhaseSpaceDecayChannel(thisdec[0],branch,3,thisdec[1],thisdec[2],thisdec[3]);
198  } else {
199  G4cout<<"Decay chain invalid!!!"<<std::endl;
200  }
201  }
202  for (G4int index=0; index <ndec; index++ ) { table->Insert(mode[index]); }
203  delete [] mode;
204  part->SetDecayTable(table);
205  }
206  else { continue; }
207  }
208  return particles;
209 }
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:17
checkFileSG.line
line
Definition: checkFileSG.py:75
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
taskman.configFile
configFile
Definition: taskman.py:311
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
index
Definition: index.py:1
CustomParticleFactory::loadCustomParticles
static void loadCustomParticles()
Definition: CustomParticleFactory.cxx:32
SUSY_SimplifiedModel_PreInclude.decays
dictionary decays
Definition: SUSY_SimplifiedModel_PreInclude.py:12
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
PhysicsConfigurationHelper.h
lumiFormat.i
int i
Definition: lumiFormat.py:85
EvtGen_Fragment.decayFile
decayFile
Definition: EvtGen_Fragment.py:27
CustomParticleFactory.h
MCTruthPartClassifier::stable
@ stable
Definition: TruthClassifiers.h:148
Preparation.mode
mode
Definition: Preparation.py:107
CustomParticleFactory::isCustomParticle
static bool isCustomParticle(G4ParticleDefinition *particle)
Definition: CustomParticleFactory.cxx:26
CustomParticleFactory::load
static std::set< G4ParticleDefinition * > load()
Definition: CustomParticleFactory.cxx:38
CxxUtils::atof
double atof(std::string_view str)
Converts a string into a double / float.
Definition: Control/CxxUtils/Root/StringUtils.cxx:91
CustomParticle
Definition: CustomParticle.h:18
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
python.ext.table_printer.table
list table
Definition: table_printer.py:81
charge
double charge(const T &p)
Definition: AtlasPID.h:897
isSlepton
bool isSlepton(const T &p)
Definition: AtlasPID.h:431
spin2
int spin2(const T &p)
Definition: AtlasPID.h:986
RTTAlgmain.branch
branch
Definition: RTTAlgmain.py:61
CustomParticle.h
DeMoScan.first
bool first
Definition: DeMoScan.py:536
LArG4FSStartPointFilter.particles
list particles
Definition: LArG4FSStartPointFilter.py:84
isRHadron
bool isRHadron(const T &p)
Definition: AtlasPID.h:1012
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
python.SystemOfUnits.eplus
int eplus
Definition: SystemOfUnits.py:137
PhysicsConfigurationHelper::Instance
static const PhysicsConfigurationHelper * Instance()
Definition: PhysicsConfigurationHelper.cxx:54
HepMCHelpers.h
PhysicsConfigurationHelper
Definition: PhysicsConfigurationHelper.h:13