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
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
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.
21 #include "DatabasePDG.h"
34 DatabasePDG::DatabasePDG():
36 fUseCharmParticles(kTRUE),
42 // Default constructor, initialize members, set input files
44 strncpy(fParticleFilename, "particles.data", 256);
45 strncpy(fDecayFilename, "tabledecay.txt", 256);
46 for(Int_t i=0; i<kMaxParticles; i++) {
47 fParticles[i] = new ParticlePDG();
52 DatabasePDG::~DatabasePDG() {
53 for(Int_t i=0; i<kMaxParticles; i++)
58 void DatabasePDG::SetParticleFilename(Char_t *filename) {
59 strncpy(fParticleFilename, filename, 255);
62 void DatabasePDG::SetDecayFilename(Char_t *filename) {
63 strncpy(fDecayFilename, filename, 255);
66 Bool_t DatabasePDG::LoadData() {
67 return (LoadParticles() && LoadDecays());
70 Bool_t DatabasePDG::LoadParticles() {
71 // Read particle definitions from the ascii file
73 ifstream particleFile;
74 particleFile.open(fParticleFilename);
76 cout << "ERROR in DatabasePDG::LoadParticles() : The ASCII file containing the PDG particle list (\""
77 << fParticleFilename << "\") was not found !! Exiting..." << endl;
82 Double_t mass, width, spin, isospin, isospinZ, q, s, aq, as, c, ac;
84 Int_t goodStatusParticles = 0;
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;
91 particleFile.exceptions(ios::failbit);
92 while(!particleFile.eof()) {
94 particleFile >> name >> mass >> width >> spin >> isospin >> isospinZ >> q >> s >> aq >> as >> c >> ac >> pdg;
96 catch (ios::failure const &problem) {
97 cout << problem.what() << endl;
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);
115 fStatus[fNParticles] = kTRUE;
116 // check if we want charmed particles
117 if(!fUseCharmParticles && (c>0 || ac>0)) {
118 fStatus[fNParticles] = kFALSE;
120 // check that the particle mass is inside accepted limits
121 if(!(fMinimumMass<=mass && mass<=fMaximumMass)) {
122 fStatus[fNParticles] = kFALSE;
124 // check that the particle width is inside accepted limits
125 if(!(fMinimumWidth<=width && width<=fMaximumWidth)) {
126 fStatus[fNParticles] = kFALSE;
128 if(fStatus[fNParticles]) goodStatusParticles++;
131 particleFile.close();
133 cout << "Warning in DatabasePDG::LoadParticles(): No particles were found in the file specified!!" << endl;
137 cout << "Info in DatabasePDG::LoadParticles(): Particle definitions found = " << fNParticles << endl
138 << " Good status particles = " << goodStatusParticles << endl;
142 Bool_t DatabasePDG::LoadDecays() {
143 // Read the decay channel definitions from the ascii file
146 decayFile.open(fDecayFilename);
148 cout << "ERROR in DatabasePDG::LoadDecays() : The ASCII file containing the decays list (\""
149 << fDecayFilename << "\") was not found !! Exiting..." << endl;
153 Int_t motherPdg, daughterPdg[3];
156 decayFile.exceptions(ios::failbit);
157 while(!decayFile.eof()) {
159 for(Int_t i=0; i<3; i++) daughterPdg[i] = 0;
162 decayFile >> motherPdg;
163 for(Int_t i=0; i<3; i++)
164 decayFile >> daughterPdg[i];
165 decayFile >> branching;
167 catch (ios::failure const &problem) {
168 cout << problem.what() << endl;
171 if((motherPdg!=0) && (daughterPdg[0]!=0) && (branching>=0)) {
172 Int_t nDaughters = 0;
173 for(Int_t i=0; i<3; i++)
174 if(daughterPdg[i]!=0)
176 ParticlePDG* particle = GetPDGParticle(motherPdg);
177 DecayChannel decay(motherPdg, branching, nDaughters, daughterPdg);
178 particle->AddChannel(&decay);
182 Int_t nDecayChannels = 0;
183 for(Int_t i=0; i<fNParticles; i++) {
184 nDecayChannels += fParticles[i]->GetNDecayChannels();
186 cout << "Info in DatabasePDG::LoadDecays(): Number of decays found in the database is " << nDecayChannels << endl;
190 ParticlePDG* DatabasePDG::GetPDGParticleByIndex(Int_t index) const {
191 // Return a PDG particle definition based on its index in the particle list
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;
199 return fParticles[index];
202 Bool_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
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;
213 return fStatus[index];
216 ParticlePDG* 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
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;
228 if(nFindings == 1) return fParticles[firstTimeIndex];
230 //cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
231 // << " was not found in the database!!" << endl;
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];
243 Bool_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
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;
256 if(nFindings == 1) return fStatus[firstTimeIndex];
258 //cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
259 // << " was not found in the database!!" << endl;
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];
271 ParticlePDG* 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
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;
283 if(nFindings == 1) return fParticles[firstTimeIndex];
285 //cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
286 // << "\" was not found in the database!!" << endl;
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];
298 Bool_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
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;
311 if(nFindings == 1) return fStatus[firstTimeIndex];
313 //cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
314 // << "\" was not found in the database!!" << endl;
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];
326 void DatabasePDG::DumpData(Bool_t dumpAll) const {
327 // Printout all the information in the PDG database
329 cout << "***********************************************************************************************************" << endl;
330 cout << "Dumping all the information contained in the database..." << endl;
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();
368 cout << " daughters sum mass = " << daughtersMass << endl;
372 cout << "Finished dumping information for " << fNParticles << " particles with " << nDecays << " decay channels in total." << endl;
373 cout << "*************************************************************************************************************" << endl;
376 cout << "Finished dumping information for " << nGoodStatusParticles << "(" << fNParticles << ")"
377 << " particles with " << nGoodStatusDecays << "(" << nDecays << ")" << " decay channels in total." << endl;
378 cout << "*************************************************************************************************************" << endl;
382 Int_t DatabasePDG::CheckImpossibleDecays(Bool_t dump) const {
383 // Check the database for impossible decays
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());
391 cout << "Particle " << fParticles[currPart]->GetPDG() << " has " << allChannels << " decay channels specified in the database" << endl;
392 cout << " Allowed channels assuming table mass = " << allowedChannels << endl;
394 if(dump && allChannels>0 && allowedChannels == 0) {
395 cout << "**********************************************************************" << endl;
396 cout << " All channels for this particles are not allowed" << endl;
397 cout << "**********************************************************************" << endl;
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;
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();
411 if(daughtersSumMass >= motherMass) {
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() << " ";
421 cout << " Mother particle mass = " << motherMass << endl;
422 cout << " Daughters sum mass = " << daughtersSumMass << endl;
427 return nImpossibleDecays;
430 void DatabasePDG::SetUseCharmParticles(Bool_t flag) {
431 // Switch on/off the use of charmed particles
434 fUseCharmParticles = flag;
435 for(Int_t i=0; i<fNParticles; i++) {
436 if(fParticles[i]->GetCharmQNumber()>0 || fParticles[i]->GetCharmAQNumber())
443 fUseCharmParticles = flag;
447 void DatabasePDG::SetMinimumWidth(Double_t value) {
448 // Set the minimum decay width for the particle definitions to be generated in the soft fireball
450 fMinimumWidth = value;
451 for(Int_t i=0; i<fNParticles; i++) {
452 if(fParticles[i]->GetWidth() < fMinimumWidth)
459 fMinimumWidth = value;
463 void DatabasePDG::SetMaximumWidth(Double_t value) {
464 // Set the maximum decay width for the particle definitions to be generated in the soft fireball
467 fMaximumWidth = value;
468 for(Int_t i=0; i<fNParticles; i++) {
469 if(fParticles[i]->GetWidth() > fMaximumWidth)
476 fMaximumWidth = value;
480 void DatabasePDG::SetWidthRange(Double_t min, Double_t max) {
481 // Set the decay width range for the particle definitions to be generated in the soft fireball
486 for(Int_t i=0; i<fNParticles; i++) {
487 if((fParticles[i]->GetWidth()<fMinimumWidth) || (fParticles[i]->GetWidth()>fMaximumWidth))
500 void DatabasePDG::SetMinimumMass(Double_t value) {
501 // Set the minimum mass for the particle definitions to be generated in the soft fireball
504 fMinimumMass = value;
505 for(Int_t i=0; i<fNParticles; i++) {
506 if(fParticles[i]->GetMass() < fMinimumMass)
513 fMinimumMass = value;
517 void DatabasePDG::SetMaximumMass(Double_t value) {
518 // Set the maximum mass for the particle definitions to be generated in the soft fireball
521 fMaximumMass = value;
522 for(Int_t i=0; i<fNParticles; i++) {
523 if(fParticles[i]->GetMass() > fMaximumMass)
530 fMaximumMass = value;
534 void DatabasePDG::SetMassRange(Double_t min, Double_t max) {
535 // Set the mass range for the particle definitions to be generated in the soft fireball
540 for(Int_t i=0; i<fNParticles; i++) {
541 if((fParticles[i]->GetMass()<fMinimumMass) || (fParticles[i]->GetMass()>fMaximumMass))
554 void DatabasePDG::SortParticles() {
555 // Sort the particle list so that those with kTRUE status will be always on top of the list
558 cout << "Warning in DatabasePDG::SortParticles() : No particles to sort. Load data first!!" << endl;
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
567 if(nGoodStatus==0) // no good status particles, again nothing to do
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;
588 Int_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
596 Int_t nGoodStatus = 0;
597 for(Int_t i=0; i<fNParticles; i++)
598 if(fStatus[i]) nGoodStatus++;
602 void DatabasePDG::UseThisListOfParticles(Char_t *filename, Bool_t exclusive) {
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)
607 cout << "Error in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : You must load the data before calling this function!!" << endl;
612 listFile.open(filename);
614 cout << "ERROR in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The ASCII file containing the PDG codes list (\""
615 << filename << "\") was not found !! Exiting..." << endl;
619 Bool_t flaggedIndexes[kMaxParticles];
620 for(Int_t i=0; i<kMaxParticles; i++)
621 flaggedIndexes[i] = kFALSE;
623 listFile.exceptions(ios::failbit);
624 while(!listFile.eof()) {
628 catch (ios::failure const &problem) {
629 cout << problem.what() << endl;
633 for(Int_t i=0; i<fNParticles; i++) {
634 if(fParticles[i]->GetPDG()==pdg) {
636 flaggedIndexes[i] = kTRUE;
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;
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;
650 for(Int_t i=0; i<kMaxParticles; i++)
651 fStatus[i] = flaggedIndexes[i];
654 for(Int_t i=0; i<kMaxParticles; i++)
655 fStatus[i] = (fStatus[i] && flaggedIndexes[i]);
662 Bool_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"
665 Double_t daughtersSumMass = 0.0;
666 for(Int_t i=0; i<channel->GetNDaughters(); i++)
667 daughtersSumMass += GetPDGParticle(channel->GetDaughterPDG(i))->GetMass();
668 if(daughtersSumMass<=motherMass)
673 Int_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
676 Int_t nAllowedChannels = 0;
677 for(Int_t i=0; i<particle->GetNDecayChannels(); i++)
678 nAllowedChannels += (IsChannelAllowed(particle->GetDecayChannel(i), motherMass) ? 1:0);
680 return nAllowedChannels;