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