Belle II Software development
EvtGenDatabasePDG.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <framework/particledb/EvtGenDatabasePDG.h>
10#include <framework/particledb/EvtGenParticlePDG.h>
11#include <framework/logging/Logger.h>
12#include <framework/gearbox/Unit.h>
13#include <framework/gearbox/Const.h>
14#include <framework/utilities/FileSystem.h>
15#include <fstream>
16#include <iomanip>
17#include <cstring>
18#include <string>
19#include <THashList.h>
20#include <TExMap.h>
21#include <TInterpreter.h>
22
23using namespace Belle2;
24
26{
27 static bool instanceCreated = false;
28 if (!instanceCreated) {
29 //hope we are the first to create a TDatabasePDG instance (setting the instance pointer in the process)
30 auto instance = new EvtGenDatabasePDG();
31 //ok, now load the data
32 instance->ReadEvtGenTable();
33 instanceCreated = true;
34 if (instance != TDatabasePDG::Instance())
35 B2FATAL("EvtGenDatabasePDG created too late!");
36 }
37 //this cast should be safe as we stop execution above if we couldn't create the first instance ourselves
38 return static_cast<EvtGenDatabasePDG*>(TDatabasePDG::Instance());
39}
40
41EvtGenParticlePDG* EvtGenDatabasePDG::AddParticle(const char* name, const char* title, Double_t mass, Bool_t stable, Double_t width,
42 Double_t charge, const char* ParticleClass, Int_t PDGcode, Int_t Anti, Int_t TrackingCode,
43 Double_t Lifetime, Double_t Spin, Double_t maxWidth, Int_t pythiaID)
44{
45 //
46 // Particle definition normal constructor. If the particle is set to be
47 // stable, the decay width parameter does have no meaning and can be set to
48 // any value. The parameters granularity, LowerCutOff and HighCutOff are
49 // used for the construction of the mean free path look up tables. The
50 // granularity will be the number of logwise energy points for which the
51 // mean free path will be calculated.
52 //
53
54 TParticlePDG* old = GetParticle(PDGcode);
55
56 if (old) {
57 B2ERROR("EvtGenDatabasePDG::AddParticle: particle with PDGcode=" << PDGcode << " already defined");
58 return nullptr;
59 }
60
61 if (std::strpbrk(name, " \t\n\v\f\r")) {
62 B2ERROR("EvtGenDatabasePDG::AddParticle: invalid particle name '" << name << "'. Names may not contain Whitespace");
63 return nullptr;
64 }
65
66 auto* p = new EvtGenParticlePDG(name, title, mass, stable, width, charge, ParticleClass,
67 PDGcode, Anti, TrackingCode, Lifetime, Spin, maxWidth, pythiaID);
68
69 fParticleList->Add(p);
70 if (fPdgMap)
71 fPdgMap->Add((Long_t)PDGcode, (Long_t)p);
72
73 TParticleClassPDG* pclass = GetParticleClass(ParticleClass);
74
75 if (!pclass) {
76 pclass = new TParticleClassPDG(ParticleClass);
77 fListOfClasses->Add(pclass);
78 }
79
80 pclass->AddParticle(p);
81
82 return p;
83}
84
85void EvtGenDatabasePDG::ReadEvtGenTable(const char* filename)
86{
87 //ok, now load the data
88 std::string defaultFilename = FileSystem::findFile("data/framework/particledb/evt.pdl");
89 if (defaultFilename.empty())
90 B2FATAL("Cannot find the default particle database file (evt.pdl).");
91 if (!filename) {
92 filename = defaultFilename.c_str();
93 } else {
94 if (filename != defaultFilename) {
95 B2RESULT("Loading non-standard evt.pdl file '" << filename << "'");
96 }
97 }
98
99 // do we have lists already?
100 if (fParticleList == nullptr) {
101 //if not create them
102 fParticleList = new THashList;
103 fListOfClasses = new TObjArray;
104 } else {
105 //otherwise clear them
106 fParticleList->Delete();
107 fListOfClasses->Delete();
108 //and also reset the map from pdgcode to particle
109 delete fPdgMap;
110 fPdgMap = nullptr;
111 }
112 std::ifstream indec(filename);
113
114 char cmnd[100], xxxx[100], pname[100];
115 int stdhepid; //aka PDG code
116 double mass, pwidth, pmaxwidth;
117 int chg3, spin2;
118 double ctau;
119 int lundkc;
120
121 if (!indec) {
122 B2ERROR("EvtGenDatabasePDG::ReadEvtGenTable: Could not particle data file '" << filename << "'");
123 return;
124 }
125
126 std::map<int, TParticlePDG*> pdgToPartMap;
127 do {
128
129 char ch, ch1;
130
131 do {
132
133 indec.get(ch);
134 if (ch == '\n') indec.get(ch);
135 if (ch != '*') {
136 indec.putback(ch);
137 } else {
138 while (indec.get(ch1), ch1 != '\n');
139 }
140 } while (ch == '*');
141
142 indec >> cmnd;
143
144 if (strcmp(cmnd, "end")) {
145
146 if (!strcmp(cmnd, "add")) {
147
148 indec >> xxxx >> xxxx >> pname >> stdhepid;
149 indec >> mass >> pwidth >> pmaxwidth >> chg3 >> spin2 >> ctau >> lundkc;
150
151 const double c_mm_per_s = Const::speedOfLight / (Unit::mm / Unit::s);
152 TParticlePDG* part = AddParticle(pname, pname, mass, false, pwidth, chg3, "Unknown", stdhepid, 0, 0,
153 ctau / c_mm_per_s, spin2 / 2.0, pmaxwidth, lundkc);
154 pdgToPartMap[stdhepid] = part;
155 if (!part) {
156 B2FATAL("EvtGenDatabasePDG::ReadPDGTable: Problem creating particle '" << pname << "'");
157 }
158 }
159 }
160 } while (strcmp(cmnd, "end"));
161
162 //fix up particle <-> antiparticle settings
163 for (auto entry : pdgToPartMap) {
164 int pdg = entry.first;
165 TParticlePDG* part = entry.second;
166
167 const auto& antiPartIter = pdgToPartMap.find(-pdg);
168 if (antiPartIter != pdgToPartMap.end()) {
169 //set anti-particle
170 part->SetAntiParticle(antiPartIter->second);
171 } else {
172 //we are our own antiparticle
173 part->SetAntiParticle(part);
174 }
175 }
176}
177
179{
180 const double c_mm_per_s = Const::speedOfLight / (Unit::mm / Unit::s);
181 if (fParticleList) {
182 for (TObject* object : *fParticleList) {
183 auto* part = dynamic_cast<EvtGenParticlePDG*>(object);
184 if (!part) {
185 B2FATAL("EvtGenDatabasePDG::WriteEvtgenTable: Particle "
186 << object->GetName() <<
187 " does not inherit from EventGenParticlePDG");
188 }
189 out << "add p Particle " << part->GetName() << " " << part->PdgCode()
190 << " " << std::scientific << std::setprecision(7) << part->Mass()
191 << " " << std::scientific << std::setprecision(7) << part->Width()
192 << " " << std::scientific << std::setprecision(7) << part->MaxWidth()
193 << " " << (int)(part->Charge())
194 << " " << (int)(part->Spin() * 2)
195 << " " << std::scientific << std::setprecision(7) << part->Lifetime() * c_mm_per_s
196 << " " << part->PythiaID()
197 << std::endl;
198 }
199 }
200 out << "end" << std::endl << std::flush;
201}
202
203void EvtGenDatabasePDG::WriteEvtGenTable(const char* filename)
204{
205 std::ofstream out(filename);
206 if (!out) {
207 B2FATAL("Cannot write evtgen pdl to '" << filename << "'");
208 }
209 WriteEvtGenTable(out);
210}
211
213{
214 if (!fParticleList)
215 return true;
216 bool testPassed = true;
217 for (TObject* object : *fParticleList) {
218 EvtGenParticlePDG* particle = dynamic_cast<EvtGenParticlePDG*>(object);
219 if (!particle) {
220 B2FATAL("EvtGenDatabasePDG::testParticleData: Particle"
221 << object->GetName() <<
222 " does not inherit from EventGenParticlePDG.");
223 }
224 /*
225 * Check width and mass range. The mass range in EvtGen is
226 * [M - max_Dm, M + 15 * Gamma]. Thus, max_Dm values are required to
227 * be less than or equal to 15 * Gamma, except for zero width. Zero width
228 * and non-zero max_Dm is a valid case for some special particles like
229 * quark pairs.
230 */
231 if ((particle->MaxWidth() > 15.000000000001 * particle->Width()) &&
232 (particle->Width() != 0)) {
233 B2ERROR("Maximum Delta_M is greater than 15 widths for "
234 << particle->GetName()
235 << LogVar("Width", particle->Width())
236 << LogVar("max_Dm", particle->MaxWidth()));
237 testPassed = false;
238 }
239 }
240 return testPassed;
241}
static const double speedOfLight
[cm/ns]
Definition Const.h:695
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
void WriteEvtGenTable(std::ostream &out)
Write current database as EvtGen table to a stream.
bool testParticleData()
Test particle data.
void ReadEvtGenTable(const char *filename=nullptr)
Read EvtGen table.
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
Helper class for setting additional TParticlePDG members and storing the ones it doesn't have yet.
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...
static const double mm
[millimeters]
Definition Unit.h:70
static const double s
[second]
Definition Unit.h:95
Class to store variables with their name which were sent to the logging service.
Abstract base class for different kinds of events.