Coding conventions (Ionut)
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / DatabasePDG.cxx
CommitLineData
3fa37a65 1// DatabasePDG stores and handles PDG information
2// The PDG particle definitions and decay channels are read
3// in the begining from ASCII files
4// PDG definitions loaded can be selected according to their
5// mass and decay width
6
b1c2e580 7/*
8 Copyright : The FASTMC and SPHMC Collaboration
9 Author : Ionut Cristian Arsene
10 Affiliation : Oslo University, Norway & Institute for Space Sciences, Bucharest, Romania
11 e-mail : i.c.arsene@fys.uio.no
12 Date : 2007/05/30
13
14 This class is using the particle and decay lists provided by the
15 THERMINATOR (Computer Physics Communications 174 669 (2006)) and
16 SHARE (Computer Physics Communications 167 229 (2005)) collaborations.
17*/
18
19
3fa37a65 20#ifndef DATABASEPDG_H
b1c2e580 21#include "DatabasePDG.h"
22#endif
23
24#include <iostream>
25#include <cstring>
26#include <fstream>
27using std::cout;
28using std::endl;
29using std::strcmp;
30using std::strcpy;
31using std::ifstream;
32using std::ios;
33
786056a2 34DatabasePDG::DatabasePDG():
35 fNParticles(0),
36 fUseCharmParticles(kTRUE),
7b7936e9 37 fMinimumWidth(0.0),
786056a2 38 fMaximumWidth(10.),
7b7936e9 39 fMinimumMass(0.0001),
786056a2 40 fMaximumMass(10.)
41{
3fa37a65 42 // Default constructor, initialize members, set input files
786056a2 43
b1c2e580 44 strcpy(fParticleFilename, "particles.data");
45 strcpy(fDecayFilename, "tabledecay.txt");
46 for(Int_t i=0; i<kMaxParticles; i++) {
47 fParticles[i] = new ParticlePDG();
48 fStatus[i] = kFALSE;
49 }
b1c2e580 50}
51
52DatabasePDG::~DatabasePDG() {
53 for(Int_t i=0; i<kMaxParticles; i++)
54 if(fParticles[i])
55 delete fParticles[i];
56}
57
58void DatabasePDG::SetParticleFilename(Char_t *filename) {
59 strcpy(fParticleFilename, filename);
60}
61
62void DatabasePDG::SetDecayFilename(Char_t *filename) {
63 strcpy(fDecayFilename, filename);
64}
65
66Bool_t DatabasePDG::LoadData() {
67 return (LoadParticles() && LoadDecays());
68}
69
70Bool_t DatabasePDG::LoadParticles() {
3fa37a65 71 // Read particle definitions from the ascii file
72
b1c2e580 73 ifstream particleFile;
74 particleFile.open(fParticleFilename);
75 if(!particleFile) {
76 cout << "ERROR in DatabasePDG::LoadParticles() : The ASCII file containing the PDG particle list (\""
3eeebf49 77 << fParticleFilename << "\") was not found !! Exiting..." << endl;
b1c2e580 78 return kFALSE;
79 }
80
81 Char_t name[9];
82 Double_t mass, width, spin, isospin, isospinZ, q, s, aq, as, c, ac;
83 Int_t pdg;
84 Int_t goodStatusParticles = 0;
85
86 cout << "Info in DatabasePDG::LoadParticles() : Start loading particles with the following criteria:" << endl
87 << " Use particles containing charm quarks (1-yes;0-no) : " << fUseCharmParticles << endl
88 << " Mass range : (" << fMinimumMass << "; " << fMaximumMass << ")" << endl
89 << " Width range : (" << fMinimumWidth << "; " << fMaximumWidth << ")" << endl;
90
91 particleFile.exceptions(ios::failbit);
92 while(!particleFile.eof()) {
93 try {
94 particleFile >> name >> mass >> width >> spin >> isospin >> isospinZ >> q >> s >> aq >> as >> c >> ac >> pdg;
95 }
96 catch (ios::failure const &problem) {
97 cout << problem.what() << endl;
98 break;
99 }
100
101 fParticles[fNParticles]->SetName(name);
102 fParticles[fNParticles]->SetPDG(pdg);
103 fParticles[fNParticles]->SetMass(mass);
104 fParticles[fNParticles]->SetWidth(width);
105 fParticles[fNParticles]->SetSpin(spin);
106 fParticles[fNParticles]->SetIsospin(isospin);
107 fParticles[fNParticles]->SetIsospinZ(isospinZ);
108 fParticles[fNParticles]->SetLightQNumber(q);
109 fParticles[fNParticles]->SetStrangeQNumber(s);
110 fParticles[fNParticles]->SetLightAQNumber(aq);
111 fParticles[fNParticles]->SetStrangeAQNumber(as);
112 fParticles[fNParticles]->SetCharmQNumber(c);
113 fParticles[fNParticles]->SetCharmAQNumber(ac);
114
115 fStatus[fNParticles] = kTRUE;
116 // check if we want charmed particles
117 if(!fUseCharmParticles && (c>0 || ac>0)) {
118 fStatus[fNParticles] = kFALSE;
119 }
120 // check that the particle mass is inside accepted limits
121 if(!(fMinimumMass<=mass && mass<=fMaximumMass)) {
122 fStatus[fNParticles] = kFALSE;
123 }
124 // check that the particle width is inside accepted limits
125 if(!(fMinimumWidth<=width && width<=fMaximumWidth)) {
126 fStatus[fNParticles] = kFALSE;
127 }
128 if(fStatus[fNParticles]) goodStatusParticles++;
129 fNParticles++;
130 }
131 particleFile.close();
132 if(fNParticles==0) {
133 cout << "Warning in DatabasePDG::LoadParticles(): No particles were found in the file specified!!" << endl;
134 return kFALSE;
135 }
136 SortParticles();
137 cout << "Info in DatabasePDG::LoadParticles(): Particle definitions found = " << fNParticles << endl
138 << " Good status particles = " << goodStatusParticles << endl;
139 return kTRUE;
140}
141
142Bool_t DatabasePDG::LoadDecays() {
3fa37a65 143 // Read the decay channel definitions from the ascii file
144
b1c2e580 145 ifstream decayFile;
146 decayFile.open(fDecayFilename);
147 if(!decayFile) {
148 cout << "ERROR in DatabasePDG::LoadDecays() : The ASCII file containing the decays list (\""
3fa37a65 149 << fDecayFilename << "\") was not found !! Exiting..." << endl;
b1c2e580 150 return kFALSE;
151 }
152
3fa37a65 153 Int_t motherPdg, daughterPdg[3];
b1c2e580 154 Double_t branching;
155
156 decayFile.exceptions(ios::failbit);
157 while(!decayFile.eof()) {
3fa37a65 158 motherPdg = 0;
159 for(Int_t i=0; i<3; i++) daughterPdg[i] = 0;
b1c2e580 160 branching = -1.0;
161 try {
3fa37a65 162 decayFile >> motherPdg;
b1c2e580 163 for(Int_t i=0; i<3; i++)
3fa37a65 164 decayFile >> daughterPdg[i];
b1c2e580 165 decayFile >> branching;
166 }
167 catch (ios::failure const &problem) {
168 cout << problem.what() << endl;
169 break;
170 }
3fa37a65 171 if((motherPdg!=0) && (daughterPdg[0]!=0) && (branching>=0)) {
b1c2e580 172 Int_t nDaughters = 0;
173 for(Int_t i=0; i<3; i++)
3fa37a65 174 if(daughterPdg[i]!=0)
b1c2e580 175 nDaughters++;
3fa37a65 176 ParticlePDG* particle = GetPDGParticle(motherPdg);
177 DecayChannel decay(motherPdg, branching, nDaughters, daughterPdg);
b1c2e580 178 particle->AddChannel(&decay);
179 }
180 }
181 decayFile.close();
182 Int_t nDecayChannels = 0;
183 for(Int_t i=0; i<fNParticles; i++) {
184 nDecayChannels += fParticles[i]->GetNDecayChannels();
185 }
186 cout << "Info in DatabasePDG::LoadDecays(): Number of decays found in the database is " << nDecayChannels << endl;
187 return kTRUE;
188}
189
3fa37a65 190ParticlePDG* DatabasePDG::GetPDGParticleByIndex(Int_t index) const {
191 // Return a PDG particle definition based on its index in the particle list
192
b1c2e580 193 if(index<0 || index>fNParticles) {
194 cout << "Warning in DatabasePDG::GetPDGParticleByIndex(Int_t): Particle index is negative or too big !!" << endl
195 << " It must be inside this range: [0, " << fNParticles-1 << "]" << endl
196 << " Returning null pointer!!" << endl;
197 return 0x0;
198 }
199 return fParticles[index];
200}
201
3fa37a65 202Bool_t DatabasePDG::GetPDGParticleStatusByIndex(Int_t index) const {
203 // Return the status of a PDG particle definition based on its index in the particle list
204 // The status is kTRUE when a particle passed the mass, width and charm criteria
205 // and kFALSE otherwise
206
b1c2e580 207 if(index<0 || index>fNParticles) {
208 cout << "Warning in DatabasePDG::GetPDGParticleStatusByIndex(Int_t): Particle index is negative or too big !!" << endl
209 << " It must be inside this range: [0, " << fNParticles-1 << "]" << endl
210 << " Returning null pointer!!" << endl;
211 return kFALSE;
212 }
213 return fStatus[index];
214}
215
3fa37a65 216ParticlePDG* DatabasePDG::GetPDGParticle(Int_t pdg) const {
217 // Return a PDG particle definition based on the PDG code (PYTHIA convention)
218 // If more than 1 definition with the asked PDG code is found then a warning is issued
219
b1c2e580 220 Int_t nFindings = 0;
221 Int_t firstTimeIndex = 0;
222 for(Int_t i=0; i<fNParticles; i++) {
223 if(pdg == fParticles[i]->GetPDG()) {
224 if(nFindings == 0) firstTimeIndex = i;
225 nFindings++;
226 }
227 }
228 if(nFindings == 1) return fParticles[firstTimeIndex];
229 if(nFindings == 0) {
230 //cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
231 // << " was not found in the database!!" << endl;
232 return 0x0;
233 }
234 if(nFindings >= 2) {
235 cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
236 << " was found with " << nFindings << " entries in the database. Check it out !!" << endl
237 << "Returning the first instance found" << endl;
238 return fParticles[firstTimeIndex];
239 }
240 return 0x0;
241}
242
3fa37a65 243Bool_t DatabasePDG::GetPDGParticleStatus(Int_t pdg) const {
244 // Return the status of a PDG particle definition based on its PDG code (PYTHIA convention)
245 // The status is kTRUE when a particle passed the mass, width and charm criteria
246 // and kFALSE otherwise
247
b1c2e580 248 Int_t nFindings = 0;
249 Int_t firstTimeIndex = 0;
250 for(Int_t i=0; i<fNParticles; i++) {
251 if(pdg == fParticles[i]->GetPDG()) {
252 if(nFindings == 0) firstTimeIndex = i;
253 nFindings++;
254 }
255 }
256 if(nFindings == 1) return fStatus[firstTimeIndex];
257 if(nFindings == 0) {
258 //cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
259 // << " was not found in the database!!" << endl;
260 return kFALSE;
261 }
262 if(nFindings >= 2) {
263 cout << "Warning in DatabasePDG::GetPDGParticleStatus(Int_t): The particle status required for PDG = " << pdg
264 << " was found with " << nFindings << " entries in the database. Check it out !!" << endl
265 << "Returning the status of first instance found" << endl;
266 return fStatus[firstTimeIndex];
267 }
268 return kFALSE;
269}
270
3fa37a65 271ParticlePDG* DatabasePDG::GetPDGParticle(Char_t* name) const {
272 // Return a PDG particle definition based on its name
273 // If more than 1 definition with the asked PDG code is found then a warning is issued
274
b1c2e580 275 Int_t nFindings = 0;
276 Int_t firstTimeIndex = 0;
277 for(Int_t i=0; i<fNParticles; i++) {
278 if(!strcmp(name, fParticles[i]->GetName())) {
279 if(nFindings == 0) firstTimeIndex = i;
280 nFindings++;
281 }
282 }
283 if(nFindings == 1) return fParticles[firstTimeIndex];
284 if(nFindings == 0) {
285 //cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
286 // << "\" was not found in the database!!" << endl;
287 return 0x0;
288 }
289 if(nFindings >= 2) {
290 cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
291 << "\" was found with " << nFindings << " entries in the database. Check it out !!" << endl
292 << "Returning the first instance found" << endl;
293 return fParticles[firstTimeIndex];
294 }
295 return 0x0;
296}
297
3fa37a65 298Bool_t DatabasePDG::GetPDGParticleStatus(Char_t* name) const {
299 // Return the status of a PDG particle definition based on its name
300 // The status is kTRUE when a particle passed the mass, width and charm criteria
301 // and kFALSE otherwise
302
b1c2e580 303 Int_t nFindings = 0;
304 Int_t firstTimeIndex = 0;
305 for(Int_t i=0; i<fNParticles; i++) {
306 if(!strcmp(name, fParticles[i]->GetName())) {
307 if(nFindings == 0) firstTimeIndex = i;
308 nFindings++;
309 }
310 }
311 if(nFindings == 1) return fStatus[firstTimeIndex];
312 if(nFindings == 0) {
313 //cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
314 // << "\" was not found in the database!!" << endl;
315 return kFALSE;
316 }
317 if(nFindings >= 2) {
318 cout << "Warning in DatabasePDG::GetPDGParticleStatus(Char_t*): The particle status required for name \"" << name
319 << "\" was found with " << nFindings << " entries in the database. Check it out !!" << endl
320 << "Returning the first instance found" << endl;
321 return fStatus[firstTimeIndex];
322 }
323 return kFALSE;
324}
325
3fa37a65 326void DatabasePDG::DumpData(Bool_t dumpAll) const {
327 // Printout all the information in the PDG database
328
b1c2e580 329 cout << "***********************************************************************************************************" << endl;
330 cout << "Dumping all the information contained in the database..." << endl;
331 Int_t nDecays = 0;
332 Int_t nGoodStatusDecays = 0;
333 Int_t nGoodStatusParticles = 0;
334 for(Int_t currPart=0; currPart<fNParticles; currPart++) {
335 nGoodStatusParticles += (fStatus[currPart] ? 1:0);
336 nGoodStatusDecays += (fStatus[currPart] ? fParticles[currPart]->GetNDecayChannels() : 0);
337 nDecays += fParticles[currPart]->GetNDecayChannels();
338 if(!(dumpAll || (!dumpAll && fStatus[currPart]))) continue;
339 cout << "###### Particle: " << fParticles[currPart]->GetName() << " with PDG code " << fParticles[currPart]->GetPDG() << endl;
340 cout << " status = " << fStatus[currPart] << endl;
341 cout << " mass = " << fParticles[currPart]->GetMass() << " GeV" << endl;
342 cout << " width = " << fParticles[currPart]->GetWidth() << " GeV" << endl;
343 cout << " 2*spin = " << Int_t(2.*fParticles[currPart]->GetSpin()) << endl;
344 cout << " 2*isospin = " << Int_t(2.*fParticles[currPart]->GetIsospin()) << endl;
345 cout << " 2*isospin3 = " << Int_t(2.*fParticles[currPart]->GetIsospinZ()) << endl;
346 cout << " u,d quarks = " << Int_t(fParticles[currPart]->GetLightQNumber()) << endl;
347 cout << " s quarks = " << Int_t(fParticles[currPart]->GetStrangeQNumber()) << endl;
348 cout << " c quarks = " << Int_t(fParticles[currPart]->GetCharmQNumber()) << endl;
349 cout << " anti u,d quarks = " << Int_t(fParticles[currPart]->GetLightAQNumber()) << endl;
350 cout << " anti s quarks = " << Int_t(fParticles[currPart]->GetStrangeAQNumber()) << endl;
351 cout << " anti c quarks = " << Int_t(fParticles[currPart]->GetCharmAQNumber()) << endl;
352 cout << " baryon number = " << Int_t(fParticles[currPart]->GetBaryonNumber()) << endl;
353 cout << " strangeness = " << Int_t(fParticles[currPart]->GetStrangeness()) << endl;
354 cout << " charmness = " << Int_t(fParticles[currPart]->GetCharmness()) << endl;
355 cout << " electric charge = " << Int_t(fParticles[currPart]->GetElectricCharge()) << endl;
356 cout << " full branching = " << fParticles[currPart]->GetFullBranching() << endl;
357 cout << " decay modes = " << fParticles[currPart]->GetNDecayChannels() << endl;
358 for(Int_t currChannel=0; currChannel<fParticles[currPart]->GetNDecayChannels(); currChannel++) {
359 cout << " channel " << currChannel+1 << " with branching " << fParticles[currPart]->GetDecayChannel(currChannel)->GetBranching() << endl;
360 cout << " daughters PDG codes: ";
361 Double_t daughtersMass = 0.0;
362 for(Int_t currDaughter=0; currDaughter<fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters(); currDaughter++) {
363 cout << fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter) << "\t";
364 ParticlePDG *daughter = GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
365 daughtersMass += daughter->GetMass();
366 }
367 cout << endl;
368 cout << " daughters sum mass = " << daughtersMass << endl;
369 }
370 }
371 if(dumpAll) {
372 cout << "Finished dumping information for " << fNParticles << " particles with " << nDecays << " decay channels in total." << endl;
373 cout << "*************************************************************************************************************" << endl;
374 }
375 else {
376 cout << "Finished dumping information for " << nGoodStatusParticles << "(" << fNParticles << ")"
377 << " particles with " << nGoodStatusDecays << "(" << nDecays << ")" << " decay channels in total." << endl;
378 cout << "*************************************************************************************************************" << endl;
379 }
380}
381
3fa37a65 382Int_t DatabasePDG::CheckImpossibleDecays(Bool_t dump) const {
b1c2e580 383 // Check the database for impossible decays
3fa37a65 384
b1c2e580 385 Int_t nImpossibleDecays = 0;
386 for(Int_t currPart=0; currPart<fNParticles; currPart++) {
387 if(!fStatus[currPart]) continue;
388 Int_t allChannels = fParticles[currPart]->GetNDecayChannels();
389 Int_t allowedChannels = GetNAllowedChannels(fParticles[currPart], fParticles[currPart]->GetMass());
390 if(dump) {
391 cout << "Particle " << fParticles[currPart]->GetPDG() << " has " << allChannels << " decay channels specified in the database" << endl;
392 cout << " Allowed channels assuming table mass = " << allowedChannels << endl;
393 }
394 if(dump && allChannels>0 && allowedChannels == 0) {
395 cout << "**********************************************************************" << endl;
396 cout << " All channels for this particles are not allowed" << endl;
397 cout << "**********************************************************************" << endl;
398 }
399 if(dump && fParticles[currPart]->GetWidth() > 0. && allChannels == 0) {
400 cout << "**********************************************************************" << endl;
401 cout << " Particle has finite width but no decay channels specified" << endl;
402 cout << "**********************************************************************" << endl;
403 }
404 for(Int_t currChannel=0; currChannel<fParticles[currPart]->GetNDecayChannels(); currChannel++) {
405 Double_t motherMass = fParticles[currPart]->GetMass();
406 Double_t daughtersSumMass = 0.;
407 for(Int_t currDaughter=0; currDaughter<fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters(); currDaughter++) {
408 ParticlePDG *daughter = GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
409 daughtersSumMass += daughter->GetMass();
410 }
411 if(daughtersSumMass >= motherMass) {
412 nImpossibleDecays++;
413 if(dump) {
414 cout << "Imposible decay for particle " << fParticles[currPart]->GetPDG() << endl;
415 cout << " Channel: " << fParticles[currPart]->GetPDG() << " --> ";
416 for(Int_t currDaughter=0; currDaughter<fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters(); currDaughter++) {
417 ParticlePDG *daughter = GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
418 cout << daughter->GetPDG() << " ";
419 }
420 cout << endl;
421 cout << " Mother particle mass = " << motherMass << endl;
422 cout << " Daughters sum mass = " << daughtersSumMass << endl;
423 }
424 }
425 }
426 }
427 return nImpossibleDecays;
428}
429
430void DatabasePDG::SetUseCharmParticles(Bool_t flag) {
3fa37a65 431 // Switch on/off the use of charmed particles
432
b1c2e580 433 if(fNParticles>0) {
434 fUseCharmParticles = flag;
435 for(Int_t i=0; i<fNParticles; i++) {
436 if(fParticles[i]->GetCharmQNumber()>0 || fParticles[i]->GetCharmAQNumber())
437 fStatus[i] = flag;
438 }
439 SortParticles();
440 return;
441 }
442 else
443 fUseCharmParticles = flag;
444 return;
445}
446
447void DatabasePDG::SetMinimumWidth(Double_t value) {
3fa37a65 448 // Set the minimum decay width for the particle definitions to be generated in the soft fireball
b1c2e580 449 if(fNParticles>0) {
450 fMinimumWidth = value;
451 for(Int_t i=0; i<fNParticles; i++) {
452 if(fParticles[i]->GetWidth() < fMinimumWidth)
453 fStatus[i] = kFALSE;
454 }
455 SortParticles();
456 return;
457 }
458 else
459 fMinimumWidth = value;
460 return;
461}
462
463void DatabasePDG::SetMaximumWidth(Double_t value) {
3fa37a65 464 // Set the maximum decay width for the particle definitions to be generated in the soft fireball
465
b1c2e580 466 if(fNParticles>0) {
467 fMaximumWidth = value;
468 for(Int_t i=0; i<fNParticles; i++) {
469 if(fParticles[i]->GetWidth() > fMaximumWidth)
470 fStatus[i] = kFALSE;
471 }
472 SortParticles();
473 return;
474 }
475 else
476 fMaximumWidth = value;
477 return;
478}
479
480void DatabasePDG::SetWidthRange(Double_t min, Double_t max) {
3fa37a65 481 // Set the decay width range for the particle definitions to be generated in the soft fireball
482
b1c2e580 483 if(fNParticles>0) {
484 fMinimumWidth = min;
485 fMaximumWidth = max;
486 for(Int_t i=0; i<fNParticles; i++) {
487 if((fParticles[i]->GetWidth()<fMinimumWidth) || (fParticles[i]->GetWidth()>fMaximumWidth))
488 fStatus[i] = kFALSE;
489 }
490 SortParticles();
491 return;
492 }
493 else {
494 fMinimumWidth = min;
495 fMaximumWidth = max;
496 }
497 return;
498}
499
500void DatabasePDG::SetMinimumMass(Double_t value) {
3fa37a65 501 // Set the minimum mass for the particle definitions to be generated in the soft fireball
502
b1c2e580 503 if(fNParticles>0) {
504 fMinimumMass = value;
505 for(Int_t i=0; i<fNParticles; i++) {
506 if(fParticles[i]->GetMass() < fMinimumMass)
507 fStatus[i] = kFALSE;
508 }
509 SortParticles();
510 return;
511 }
512 else
513 fMinimumMass = value;
514 return;
515}
516
517void DatabasePDG::SetMaximumMass(Double_t value) {
3fa37a65 518 // Set the maximum mass for the particle definitions to be generated in the soft fireball
519
b1c2e580 520 if(fNParticles>0) {
521 fMaximumMass = value;
522 for(Int_t i=0; i<fNParticles; i++) {
523 if(fParticles[i]->GetMass() > fMaximumMass)
524 fStatus[i] = kFALSE;
525 }
526 SortParticles();
527 return;
528 }
529 else
530 fMaximumMass = value;
531 return;
532}
533
534void DatabasePDG::SetMassRange(Double_t min, Double_t max) {
3fa37a65 535 // Set the mass range for the particle definitions to be generated in the soft fireball
536
b1c2e580 537 if(fNParticles>0) {
538 fMinimumMass = min;
539 fMaximumMass = max;
540 for(Int_t i=0; i<fNParticles; i++) {
541 if((fParticles[i]->GetMass()<fMinimumMass) || (fParticles[i]->GetMass()>fMaximumMass))
542 fStatus[i] = kFALSE;
543 }
544 SortParticles();
545 return;
546 }
547 else {
548 fMinimumMass = min;
549 fMaximumMass = max;
550 }
551 return;
552}
553
554void DatabasePDG::SortParticles() {
3fa37a65 555 // Sort the particle list so that those with kTRUE status will be always on top of the list
556
b1c2e580 557 if(fNParticles<2) {
558 cout << "Warning in DatabasePDG::SortParticles() : No particles to sort. Load data first!!" << endl;
559 return;
560 }
561
562 Int_t nGoodStatus = 0;
563 for(Int_t i=0; i<fNParticles; i++)
564 if(fStatus[i]) nGoodStatus++;
565 if(nGoodStatus==fNParticles) // if all particles have good status then there is nothing to do
566 return;
567 if(nGoodStatus==0) // no good status particles, again nothing to do
568 return;
569
570 Int_t shifts = 1;
571 while(shifts) {
572 shifts = 0;
573 for(Int_t i=0; i<fNParticles-1; i++) {
574 if(!fStatus[i] && fStatus[i+1]) { // switch if false status is imediately before a true status particle
575 ParticlePDG *temporaryPointer = fParticles[i];
576 fParticles[i] = fParticles[i+1];
577 fParticles[i+1] = temporaryPointer;
578 Bool_t temporaryStatus = fStatus[i];
579 fStatus[i] = fStatus[i+1];
580 fStatus[i+1] = temporaryStatus;
581 shifts++;
582 }
583 }
584 }
585 return;
586}
587
3fa37a65 588Int_t DatabasePDG::GetNParticles(Bool_t all) const {
589 // Return the number of particle definitions in the database
590 // If all is kTRUE then return number of all particle
591 // If all is kFALSE then return the number of good status (kTRUE) particles
592
b1c2e580 593 if(all)
594 return fNParticles;
595
596 Int_t nGoodStatus = 0;
597 for(Int_t i=0; i<fNParticles; i++)
598 if(fStatus[i]) nGoodStatus++;
599 return nGoodStatus;
600}
601
602void DatabasePDG::UseThisListOfParticles(Char_t *filename, Bool_t exclusive) {
3fa37a65 603 // Read a list of PDG codes from the file "filename" and mark them with good status (kTRUE)
604 // while all the other will be marked kFALSE (only if exclusive = kTRUE)
605
b1c2e580 606 if(fNParticles<1) {
607 cout << "Error in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : You must load the data before calling this function!!" << endl;
608 return;
609 }
610
611 ifstream listFile;
612 listFile.open(filename);
613 if(!listFile) {
614 cout << "ERROR in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The ASCII file containing the PDG codes list (\""
3fa37a65 615 << filename << "\") was not found !! Exiting..." << endl;
b1c2e580 616 return;
617 }
618
619 Bool_t flaggedIndexes[kMaxParticles];
620 for(Int_t i=0; i<kMaxParticles; i++)
621 flaggedIndexes[i] = kFALSE;
622 Int_t pdg = 0;
623 listFile.exceptions(ios::failbit);
624 while(!listFile.eof()) {
625 try {
626 listFile >> pdg;
627 }
628 catch (ios::failure const &problem) {
629 cout << problem.what() << endl;
630 break;
631 }
632 Int_t found = 0;
633 for(Int_t i=0; i<fNParticles; i++) {
634 if(fParticles[i]->GetPDG()==pdg) {
635 found++;
636 flaggedIndexes[i] = kTRUE;
637 }
638 }
639 if(!found) {
640 cout << "Warning in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The particle with PDG code "
641 << pdg << " was asked but not found in the database!!" << endl;
642 }
643 if(found>1) {
644 cout << "Warning in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The particle with PDG code "
645 << pdg << " was found more than once in the database!!" << endl;
646 }
647 }
648
649 if(exclusive) {
650 for(Int_t i=0; i<kMaxParticles; i++)
651 fStatus[i] = flaggedIndexes[i];
652 }
653 else {
654 for(Int_t i=0; i<kMaxParticles; i++)
655 fStatus[i] = (fStatus[i] && flaggedIndexes[i]);
656 }
657 SortParticles();
658
659 return;
660}
661
3fa37a65 662Bool_t DatabasePDG::IsChannelAllowed(DecayChannel *channel, Double_t motherMass) const {
663 // Check if the decay channel "channel" is allowed by using the mother particle mass "motherMass"
664
b1c2e580 665 Double_t daughtersSumMass = 0.0;
666 for(Int_t i=0; i<channel->GetNDaughters(); i++)
667 daughtersSumMass += GetPDGParticle(channel->GetDaughterPDG(i))->GetMass();
7b7936e9 668 if(daughtersSumMass<=motherMass)
b1c2e580 669 return kTRUE;
670 return kFALSE;
671}
672
3fa37a65 673Int_t DatabasePDG::GetNAllowedChannels(ParticlePDG *particle, Double_t motherMass) const {
674 // Check how many decay channels are allowed for a given particle definition at a given mass
675
b1c2e580 676 Int_t nAllowedChannels = 0;
677 for(Int_t i=0; i<particle->GetNDecayChannels(); i++)
678 nAllowedChannels += (IsChannelAllowed(particle->GetDecayChannel(i), motherMass) ? 1:0);
679 return nAllowedChannels;
680}