Belle II Software  release-05-02-19
EvtGenDatabasePDG.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2015 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <framework/particledb/EvtGenDatabasePDG.h>
12 #include <framework/particledb/EvtGenParticlePDG.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Unit.h>
15 #include <framework/gearbox/Const.h>
16 #include <framework/utilities/FileSystem.h>
17 #include <fstream>
18 #include <iomanip>
19 #include <cstring>
20 #include <string>
21 #include <THashList.h>
22 #include <TExMap.h>
23 #include <TInterpreter.h>
24 
25 using namespace Belle2;
26 
28 {
29  static bool instanceCreated = false;
30  if (!instanceCreated) {
31  // problem with 6.14: ROOT cannot find includes this early. Fix promised for 6.16.
32  // (https://root-forum.cern.ch/t/problem-with-loadclassinfo-for-class-instantiated-very-early/30831)
33  gInterpreter->Declare("#include <framework/particledb/EvtGenParticlePDG.h>");
34  //hope we are the first to create a TDatabasePDG instance (setting the instance pointer in the process)
35  auto instance = new EvtGenDatabasePDG();
36  //ok, now load the data
37  instance->ReadEvtGenTable();
38  instanceCreated = true;
39  if (instance != TDatabasePDG::Instance())
40  B2FATAL("EvtGenDatabasePDG created too late!");
41  }
42  //this cast should be safe as we stop execution above if we couldn't create the first instance ourselves
43  return static_cast<EvtGenDatabasePDG*>(TDatabasePDG::Instance());
44 }
45 
46 EvtGenParticlePDG* EvtGenDatabasePDG::AddParticle(const char* name, const char* title, Double_t mass, Bool_t stable, Double_t width,
47  Double_t charge, const char* ParticleClass, Int_t PDGcode, Int_t Anti, Int_t TrackingCode,
48  Double_t Lifetime, Double_t Spin, Double_t maxWidth, Int_t pythiaID)
49 {
50  //
51  // Particle definition normal constructor. If the particle is set to be
52  // stable, the decay width parameter does have no meaning and can be set to
53  // any value. The parameters granularity, LowerCutOff and HighCutOff are
54  // used for the construction of the mean free path look up tables. The
55  // granularity will be the number of logwise energy points for which the
56  // mean free path will be calculated.
57  //
58 
59  TParticlePDG* old = GetParticle(PDGcode);
60 
61  if (old) {
62  B2ERROR("EvtGenDatabasePDG::AddParticle: particle with PDGcode=" << PDGcode << " already defined");
63  return nullptr;
64  }
65 
66  if (std::strpbrk(name, " \t\n\v\f\r")) {
67  B2ERROR("EvtGenDatabasePDG::AddParticle: invalid particle name '" << name << "'. Names may not contain Whitespace");
68  return nullptr;
69  }
70 
71  auto* p = new EvtGenParticlePDG(name, title, mass, stable, width, charge, ParticleClass,
72  PDGcode, Anti, TrackingCode, Lifetime, Spin, maxWidth, pythiaID);
73 
74  fParticleList->Add(p);
75  if (fPdgMap)
76  fPdgMap->Add((Long_t)PDGcode, (Long_t)p);
77 
78  TParticleClassPDG* pclass = GetParticleClass(ParticleClass);
79 
80  if (!pclass) {
81  pclass = new TParticleClassPDG(ParticleClass);
82  fListOfClasses->Add(pclass);
83  }
84 
85  pclass->AddParticle(p);
86 
87  return p;
88 }
89 
90 void EvtGenDatabasePDG::ReadEvtGenTable(const char* filename)
91 {
92  //ok, now load the data
93  std::string defaultFilename = FileSystem::findFile("data/framework/particledb/evt.pdl");
94  if (defaultFilename.empty())
95  B2FATAL("Cannot find the default particle database file (evt.pdl).");
96  if (!filename) {
97  filename = defaultFilename.c_str();
98  } else {
99  if (filename != defaultFilename) {
100  B2RESULT("Loading non-standard evt.pdl file '" << filename << "'");
101  }
102  }
103 
104  // do we have lists already?
105  if (fParticleList == nullptr) {
106  //if not create them
107  fParticleList = new THashList;
108  fListOfClasses = new TObjArray;
109  } else {
110  //otherwise clear them
111  fParticleList->Delete();
112  fListOfClasses->Delete();
113  //and also reset the map from pdgcode to particle
114  delete fPdgMap;
115  fPdgMap = nullptr;
116  }
117  std::ifstream indec(filename);
118 
119  char cmnd[100], xxxx[100], pname[100];
120  int stdhepid; //aka PDG code
121  double mass, pwidth, pmaxwidth;
122  int chg3, spin2;
123  double ctau;
124  int lundkc;
125 
126  if (!indec) {
127  B2ERROR("EvtGenDatabasePDG::ReadEvtGenTable: Could not particle data file '" << filename << "'");
128  return;
129  }
130 
131  std::map<int, TParticlePDG*> pdgToPartMap;
132  do {
133 
134  char ch, ch1;
135 
136  do {
137 
138  indec.get(ch);
139  if (ch == '\n') indec.get(ch);
140  if (ch != '*') {
141  indec.putback(ch);
142  } else {
143  while (indec.get(ch1), ch1 != '\n');
144  }
145  } while (ch == '*');
146 
147  indec >> cmnd;
148 
149  if (strcmp(cmnd, "end")) {
150 
151  if (!strcmp(cmnd, "add")) {
152 
153  indec >> xxxx >> xxxx >> pname >> stdhepid;
154  indec >> mass >> pwidth >> pmaxwidth >> chg3 >> spin2 >> ctau >> lundkc;
155 
156  const double c_mm_per_s = Const::speedOfLight / (Unit::mm / Unit::s);
157  TParticlePDG* part = AddParticle(pname, pname, mass, false, pwidth, chg3, "Unknown", stdhepid, 0, 0,
158  ctau / c_mm_per_s, spin2 / 2.0, pmaxwidth, lundkc);
159  pdgToPartMap[stdhepid] = part;
160  if (!part) {
161  B2FATAL("EvtGenDatabasePDG::ReadPDGTable: Problem creating particle '" << pname << "'");
162  }
163  }
164  }
165  } while (strcmp(cmnd, "end"));
166 
167  //fix up particle <-> antiparticle settings
168  for (auto entry : pdgToPartMap) {
169  int pdg = entry.first;
170  TParticlePDG* part = entry.second;
171 
172  const auto& antiPartIter = pdgToPartMap.find(-pdg);
173  if (antiPartIter != pdgToPartMap.end()) {
174  //set anti-particle
175  part->SetAntiParticle(antiPartIter->second);
176  } else {
177  //we are our own antiparticle
178  part->SetAntiParticle(part);
179  }
180  }
181 }
182 
183 void EvtGenDatabasePDG::WriteEvtGenTable(std::ostream& out)
184 {
185  const double c_mm_per_s = Const::speedOfLight / (Unit::mm / Unit::s);
186  if (fParticleList) {
187  for (TObject* obj : *fParticleList) {
188  auto* part = dynamic_cast<EvtGenParticlePDG*>(obj);
189  if (!part) {
190  B2FATAL("EvtGenDatabasePDG::WriteEvtgenTable: Particle does not inherit from EventGenParticlePDG");
191  }
192  out << "add p Particle " << part->GetName() << " " << part->PdgCode()
193  << " " << std::scientific << std::setprecision(7) << part->Mass()
194  << " " << std::scientific << std::setprecision(7) << part->Width()
195  << " " << std::scientific << std::setprecision(7) << part->MaxWidth()
196  << " " << (int)(part->Charge())
197  << " " << (int)(part->Spin() * 2)
198  << " " << std::scientific << std::setprecision(7) << part->Lifetime() * c_mm_per_s
199  << " " << part->PythiaID()
200  << std::endl;
201  }
202  }
203  out << "end" << std::endl << std::flush;
204 }
205 
206 void EvtGenDatabasePDG::WriteEvtGenTable(const char* filename)
207 {
208  std::ofstream out(filename);
209  if (!out) {
210  B2FATAL("Cannot write evtgen pdl to '" << filename << "'");
211  }
212  WriteEvtGenTable(out);
213 }
Belle2::EvtGenDatabasePDG::ReadEvtGenTable
void ReadEvtGenTable(const char *filename=nullptr)
Read EvtGen table.
Definition: EvtGenDatabasePDG.cc:90
Belle2::Unit::s
static const double s
[second]
Definition: Unit.h:105
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::EvtGenDatabasePDG::WriteEvtGenTable
void WriteEvtGenTable(std::ostream &out)
Write current database as EvtGen table to a stream.
Definition: EvtGenDatabasePDG.cc:183
Belle2::EvtGenDatabasePDG
Replacement for TDatabasePDG that is filled from EvtGen's evt.pdl.
Definition: EvtGenDatabasePDG.h:31
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::EvtGenDatabasePDG::EvtGenDatabasePDG
EvtGenDatabasePDG()
singleton.
Definition: EvtGenDatabasePDG.h:54
Belle2::EvtGenDatabasePDG::AddParticle
TParticlePDG * AddParticle(const char *name, const char *title, Double_t mass, Bool_t stable, Double_t width, Double_t charge, const char *ParticleClass, Int_t PDGcode, Int_t Anti, Int_t TrackingCode) override
override old AddParticle
Definition: EvtGenDatabasePDG.h:42
Belle2::FileSystem::findFile
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
Definition: FileSystem.cc:147
Belle2::EvtGenParticlePDG
Helper class for setting additional TParticlePDG members and storing the ones it doesn't have yet.
Definition: EvtGenParticlePDG.h:32
Belle2::EvtGenDatabasePDG::Instance
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
Definition: EvtGenDatabasePDG.cc:27