CbmRoot
ThermalParticleSystem.cxx
Go to the documentation of this file.
1 #include "ThermalParticleSystem.h"
2 #include <algorithm>
3 #include <fstream>
4 #include <iostream>
5 #include <queue>
6 #include <sstream>
7 
8 using namespace std;
9 
10 namespace ThermalParticleSystemNameSpace {
11  bool cmpParticleMass(const ThermalParticle& a, const ThermalParticle& b) {
12  return (a.fMass < b.fMass);
13  }
14 } // namespace ThermalParticleSystemNameSpace
15 
16 ThermalParticleSystem::ThermalParticleSystem(std::string InputFile)
17  : fParticles(), PDGtoID(), IDtoPDG(), NumberOfParticles(0) {
18  // NumberOfParticles = 0;
19  fParticles.resize(0);
20  PDGtoID.clear();
21  IDtoPDG.resize(0);
22  if (InputFile.size() < 4) return;
23  if (InputFile.substr(InputFile.size() - 4) == ".txt")
24  LoadTHERMUSDatabase(InputFile);
25  else
26  LoadDatabase(InputFile);
27 }
28 
29 
30 ThermalParticleSystem::~ThermalParticleSystem(void) {}
31 
32 void ThermalParticleSystem::FillResonanceDecays() {
33  for (unsigned int i = 0; i < fParticles.size(); ++i) {
34  fParticles[i].ResonanceBR.resize(0);
35  }
36  for (int i = fParticles.size() - 1; i >= 0; i--)
37  if (!fParticles[i].fStable) { GoResonance(i, i, 1.); }
38 }
39 
40 void ThermalParticleSystem::GoResonance(int ind, int startind, double BR) {
41  if (ind != startind && fParticles[ind].ResonanceBR.size() > 0
42  && fParticles[ind]
43  .ResonanceBR[fParticles[ind].ResonanceBR.size() - 1]
44  .second
45  == startind) {
46  fParticles[ind].ResonanceBR[fParticles[ind].ResonanceBR.size() - 1].first +=
47  BR;
48  } else if (ind != startind)
49  fParticles[ind].ResonanceBR.push_back(make_pair(BR, startind));
50  if (!fParticles[ind].fStable) {
51  for (unsigned int i = 0; i < fParticles[ind].fDecays.size(); ++i)
52  for (unsigned int j = 0; j < fParticles[ind].fDecays[i].fDaughters.size();
53  ++j) {
54  GoResonance(PDGtoID[fParticles[ind].fDecays[i].fDaughters[j]],
55  startind,
56  BR * fParticles[ind].fDecays[i].fBratio);
57  }
58  }
59 }
60 
61 void ThermalParticleSystem::LoadTHERMUSDatabase(std::string InputFile) {
62  NumberOfParticles = 0;
63  ifstream fin(InputFile.c_str());
64  fParticles.resize(0);
65  PDGtoID.clear();
66  IDtoPDG.resize(0);
67  string decayprefix = "";
68  for (int i = InputFile.size() - 1; i >= 0; --i)
69  if (InputFile[i] == '\\' || InputFile[i] == '/') {
70  decayprefix = InputFile.substr(0, i + 1);
71  break;
72  }
73  if (fin.is_open()) {
74  int stable, pdgid, spin, stat, str, bary, chg;
75  double mass, width, threshold, abss;
76  string name, decayname;
77  while (fin >> stable) {
78  fin >> name >> pdgid >> spin >> stat >> mass >> str >> bary >> chg >> abss
79  >> width >> threshold;
80  if (fin.peek() != '\n') fin >> decayname;
81  fParticles.push_back(ThermalParticle(stable,
82  name,
83  pdgid,
84  spin,
85  stat,
86  mass,
87  str,
88  bary,
89  chg,
90  abss,
91  width,
92  threshold,
93  decayname));
94  NumberOfParticles++;
95  IDtoPDG.push_back(pdgid);
96  PDGtoID[pdgid] = IDtoPDG.size() - 1;
97  fParticles[fParticles.size() - 1].ReadDecays(decayprefix + name
98  + "_decay.txt");
99  if (!(bary == 0 && chg == 0 && str == 0) || name.substr(0, 3) == "K0S"
100  || name.substr(0, 3) == "K0L") {
101  if (bary == 0 && name[name.size() - 1] == '+')
102  name[name.size() - 1] = '-';
103  else if (bary == 0 && name[name.size() - 1] == '-')
104  name[name.size() - 1] = '+';
105  else
106  name += string("bar");
107  fParticles.push_back(ThermalParticle(stable,
108  name,
109  -pdgid,
110  spin,
111  stat,
112  mass,
113  -str,
114  -bary,
115  -chg,
116  abss,
117  width,
118  threshold,
119  decayname));
120  IDtoPDG.push_back(-pdgid);
121  PDGtoID[-pdgid] = IDtoPDG.size() - 1;
122  }
123  }
124  }
125  sort(fParticles.begin(),
126  fParticles.end(),
127  ThermalParticleSystemNameSpace::cmpParticleMass);
128  PDGtoID.clear();
129  IDtoPDG.resize(0);
130  for (unsigned int i = 0; i < fParticles.size(); ++i) {
131  IDtoPDG.push_back(fParticles[i].fPDGID);
132  PDGtoID[fParticles[i].fPDGID] = IDtoPDG.size() - 1;
133  }
134  for (unsigned int i = 0; i < fParticles.size(); ++i) {
135  if (fParticles[i].fPDGID < 0)
136  fParticles[i].SetDecays(GetDecaysFromAntiParticle(
137  fParticles[PDGtoID[-fParticles[i].fPDGID]].fDecays));
138  }
139  fin.close();
140  FillResonanceDecays();
141 }
142 
143 std::vector<std::string>&
144 split(const std::string& s, char delim, std::vector<std::string>& elems) {
145  std::stringstream ss(s);
146  std::string item;
147  while (std::getline(ss, item, delim)) {
148  elems.push_back(item);
149  }
150  return elems;
151 }
152 
153 
154 std::vector<std::string> split(const std::string& s, char delim) {
155  std::vector<std::string> elems;
156  split(s, delim, elems);
157  return elems;
158 }
159 
160 void ThermalParticleSystem::LoadDatabase(std::string InputFile) {
161  vector<vector<ParticleDecay>> decays(0);
162  vector<int> pdgids(0);
163  map<int, int> decaymap;
164  decaymap.clear();
165  NumberOfParticles = 0;
166  fParticles.resize(0);
167  PDGtoID.clear();
168  IDtoPDG.resize(0);
169  string decayprefix = "";
170 
171  for (int i = InputFile.size() - 1; i >= 0; --i)
172  if (InputFile[i] == '\\' || InputFile[i] == '/') {
173  decayprefix = InputFile.substr(0, i + 1);
174  break;
175  }
176 
177  ifstream fin((decayprefix + "decays.dat").c_str());
178  if (fin.is_open()) {
179  int decaypartnumber = 0;
180  fin >> decaypartnumber;
181  decays.reserve(decaypartnumber);
182  pdgids.reserve(decaypartnumber);
183 
184  for (unsigned int i = 0; i < static_cast<unsigned int>(decaypartnumber);
185  ++i) {
186  int pdgid, decaysnumber, tmpid, daughters;
187  double bratio;
188  fin >> pdgid >> decaysnumber;
189  decaymap[pdgid] = i;
190  decays.push_back(vector<ParticleDecay>(0));
191  pdgids.push_back(pdgid);
192  for (unsigned int j = 0; j < static_cast<unsigned int>(decaysnumber);
193  ++j) {
194  ParticleDecay decay;
195  fin >> bratio;
196  decay.fBratio = bratio / 100.;
197  fin >> daughters;
198  decay.fDaughters.reserve(daughters);
199  for (unsigned int k = 0; k < static_cast<unsigned int>(daughters);
200  ++k) {
201  fin >> tmpid;
202  decay.fDaughters.push_back(tmpid);
203  }
204  decays[i].push_back(decay);
205  }
206  }
207  fin.close();
208  }
209 
210  // fin.open(InputFile);
211  ifstream fin2(InputFile.c_str());
212  if (fin2.is_open()) {
213  string tmp;
214  char tmpc[500];
215  fin2.getline(tmpc, 500);
216  tmp = string(tmpc);
217  while (1) {
218  vector<string> fields = split(tmp, ';');
219  if (fields.size() < 14) break;
220  int stable, pdgid, spin, stat, str, bary, chg, charm;
221  double mass, width, threshold, abss, absc, radius = 0.5;
222  string name, decayname = "";
223  stable = atoi(fields[0].c_str());
224  name = fields[1];
225  pdgid = atoi(fields[2].c_str());
226  spin = atoi(fields[3].c_str());
227  stat = atoi(fields[4].c_str());
228  mass = atof(fields[5].c_str());
229  str = atoi(fields[6].c_str());
230  bary = atoi(fields[7].c_str());
231  chg = atoi(fields[8].c_str());
232  charm = atoi(fields[9].c_str());
233  abss = atof(fields[10].c_str());
234  absc = atof(fields[11].c_str());
235  width = atof(fields[12].c_str());
236  threshold = atof(fields[13].c_str());
237  if (fields.size() >= 15) radius = atof(fields[14].c_str());
238  if (fields.size() == 16) decayname = fields[15];
239 
240  fParticles.push_back(ThermalParticle(stable,
241  name,
242  pdgid,
243  spin,
244  stat,
245  mass,
246  str,
247  bary,
248  chg,
249  abss,
250  width,
251  threshold,
252  decayname,
253  charm,
254  absc,
255  radius));
256  NumberOfParticles++;
257  IDtoPDG.push_back(pdgid);
258  PDGtoID[pdgid] = IDtoPDG.size() - 1;
259 
260  if (decaymap.count(pdgid) != 0)
261  fParticles[fParticles.size() - 1].fDecays = decays[decaymap[pdgid]];
262 
263  if (!(bary == 0 && chg == 0 && str == 0 && charm == 0)
264  || name.substr(0, 3) == "K0S" || name.substr(0, 3) == "K0L") {
265  if (bary == 0 && name[name.size() - 1] == '+')
266  name[name.size() - 1] = '-';
267  else if (bary == 0 && name[name.size() - 1] == '-')
268  name[name.size() - 1] = '+';
269  else
270  name += string("bar");
271  fParticles.push_back(ThermalParticle(stable,
272  name,
273  -pdgid,
274  spin,
275  stat,
276  mass,
277  -str,
278  -bary,
279  -chg,
280  abss,
281  width,
282  threshold,
283  decayname,
284  -charm,
285  absc,
286  radius));
287  IDtoPDG.push_back(-pdgid);
288  PDGtoID[-pdgid] = IDtoPDG.size() - 1;
289  }
290 
291  fin2.getline(tmpc, 500);
292  tmp = string(tmpc);
293  }
294 
295  sort(fParticles.begin(),
296  fParticles.end(),
297  ThermalParticleSystemNameSpace::cmpParticleMass);
298  PDGtoID.clear();
299  IDtoPDG.resize(0);
300  for (unsigned int i = 0; i < fParticles.size(); ++i) {
301  IDtoPDG.push_back(fParticles[i].fPDGID);
302  PDGtoID[fParticles[i].fPDGID] = IDtoPDG.size() - 1;
303  }
304  for (unsigned int i = 0; i < fParticles.size(); ++i) {
305  if (fParticles[i].fPDGID < 0)
306  fParticles[i].SetDecays(GetDecaysFromAntiParticle(
307  fParticles[PDGtoID[-fParticles[i].fPDGID]].fDecays));
308  }
309  for (unsigned int i = 0; i < fParticles.size(); ++i)
310  fParticles[i].fDecaysOrig = fParticles[i].fDecays;
311  fin2.close();
312  FillResonanceDecays();
313  }
314 }
315 
316 std::string ThermalParticleSystem::GetNameFromPDG(int pdgid) {
317  if (PDGtoID.count(pdgid) == 0)
318  return string("???");
319  else
320  return fParticles[PDGtoID[pdgid]].fName;
321 }
ThermalParticleSystem::fParticles
std::vector< ThermalParticle > fParticles
Definition: ThermalParticleSystem.h:9
ThermalParticle::fMass
double fMass
Definition: ThermalParticle.h:23
split
std::vector< std::string > & split(const std::string &s, char delim, std::vector< std::string > &elems)
Definition: ThermalParticleSystem.cxx:144
ThermalParticleSystem::GoResonance
void GoResonance(int ind, int startPDG, double BR)
Definition: ThermalParticleSystem.cxx:40
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
ThermalParticleSystem::FillResonanceDecays
void FillResonanceDecays()
Definition: ThermalParticleSystem.cxx:32
ThermalParticleSystem.h
ThermalParticleSystem::PDGtoID
std::map< int, int > PDGtoID
Definition: ThermalParticleSystem.h:10
ParticleDecay::fBratio
double fBratio
Definition: ThermalParticle.h:8
ThermalParticleSystem::IDtoPDG
std::vector< int > IDtoPDG
Definition: ThermalParticleSystem.h:11
ThermalParticleSystemNameSpace
Definition: ThermalParticleSystem.cxx:10
ThermalParticleSystem::GetDecaysFromAntiParticle
std::vector< ParticleDecay > GetDecaysFromAntiParticle(const std::vector< ParticleDecay > &Decays)
Definition: ThermalParticleSystem.h:15
ThermalParticleSystem::NumberOfParticles
int NumberOfParticles
Definition: ThermalParticleSystem.h:27
ParticleDecay
Definition: ThermalParticle.h:7
ThermalParticleSystem::ThermalParticleSystem
ThermalParticleSystem(std::string InputFile="")
Definition: ThermalParticleSystem.cxx:16
ThermalParticle
Definition: ThermalParticle.h:15
ThermalParticleSystem::LoadTHERMUSDatabase
void LoadTHERMUSDatabase(std::string InputFile="")
Definition: ThermalParticleSystem.cxx:61
ParticleDecay::fDaughters
std::vector< int > fDaughters
Definition: ThermalParticle.h:9
ThermalParticleSystem::~ThermalParticleSystem
~ThermalParticleSystem(void)
Definition: ThermalParticleSystem.cxx:30
ThermalParticleSystemNameSpace::cmpParticleMass
bool cmpParticleMass(const ThermalParticle &a, const ThermalParticle &b)
Definition: ThermalParticleSystem.cxx:11
ThermalParticleSystem::LoadDatabase
void LoadDatabase(std::string InputFile="")
Definition: ThermalParticleSystem.cxx:160
ThermalParticleSystem::GetNameFromPDG
std::string GetNameFromPDG(int pdgid)
Definition: ThermalParticleSystem.cxx:316