Changes to compile with Root6
[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / AliParticleYield.cxx
1 #include <iostream>
2 #include <stdio.h>
3 #include <fstream>
4 #include <iomanip>
5 #include "AliParticleYield.h"
6 #include "TDatabasePDG.h"
7 #include "AliLog.h"
8 #include "TClonesArray.h"
9 #include "TMath.h"
10 #include "AliPDG.h"
11 #include "TBranch.h"
12 #include "TTree.h"
13 #include "TDirectory.h"
14 #include "TEventList.h"
15
16 using std::endl;
17 using std::left;
18 using std::setw;
19 using std::ifstream;
20 using std::ofstream;
21
22 ClassImp(AliParticleYield)
23
24 // set statics
25 const char * AliParticleYield::kStatusString[] = {"Published", "Preliminary", "Final, but not published", "May change"} ;
26 const char * AliParticleYield::kSystemString[] = {"pp", "p-Pb", "Pb-Pb"} ;
27 Int_t AliParticleYield::fSignificantDigits = 3;
28 Float_t AliParticleYield::fEpsilon = 0.000000000000000001;
29
30 AliParticleYield::AliParticleYield() :
31 TObject(),
32 fPdgCode(0),
33 fPdgCode2(0),
34 fPartName(""),
35 fCollisionSystem(0),
36 fSqrtS(0),
37 fYield(0),
38 fStatError(0),
39 fSystError(0),
40 fNormErrorPos(0),
41 fNormErrorNeg(0),
42 fYMin(0),
43 fYMax(0),
44 fStatus(0),
45 fMeasurementType(0),
46 fCentr(""),
47 fIsSum(0),
48 fTag("")
49 {
50   // constructor
51   AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
52 }
53
54 AliParticleYield::AliParticleYield(Int_t pdg, Int_t system, Float_t sqrts, Float_t value, Float_t stat, Float_t syst, Float_t norm, Float_t ymin, Float_t ymax, Int_t status, Int_t type, TString centr, Int_t isSum, TString tag):
55 TObject(),
56 fPdgCode(pdg),
57 fPdgCode2(0),
58 fPartName(""),
59 fCollisionSystem(system),
60 fSqrtS(sqrts),
61 fYield(value),
62 fStatError(stat),
63 fSystError(syst),
64 fNormErrorPos(norm),
65 fNormErrorNeg(0),
66 fYMin(ymin),
67 fYMax(ymax),
68 fStatus(status),
69 fMeasurementType(type),
70 fCentr(centr),
71 fIsSum(isSum),
72 fTag(tag)
73
74 {
75   // Constructor
76   fPartName = TDatabasePDG::Instance()->GetParticle(fPdgCode)->GetName();
77   AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
78 }
79
80 AliParticleYield::AliParticleYield(Int_t pdg, Int_t system, Float_t sqrts, Float_t value, Float_t stat, Float_t syst, Float_t normPos, Float_t normNeg, Float_t ymin, Float_t ymax, Int_t status, Int_t type, TString centr, Int_t isSum, TString tag):
81 TObject(),
82 fPdgCode(pdg),
83 fPdgCode2(0),
84 fPartName(""),
85 fCollisionSystem(system),
86 fSqrtS(sqrts),
87 fYield(value),
88 fStatError(stat),
89 fSystError(syst),
90 fNormErrorPos(normPos),
91 fNormErrorNeg(normNeg),
92 fYMin(ymin),
93 fYMax(ymax),
94 fStatus(status),
95 fMeasurementType(type),
96 fCentr(centr),
97 fIsSum(isSum),
98 fTag(tag)
99
100 {
101   // Constructor
102   fPartName = TDatabasePDG::Instance()->GetParticle(fPdgCode)->GetName();
103   AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
104 }
105
106
107 AliParticleYield::AliParticleYield(const AliParticleYield& part) : 
108 TObject(),
109 fPdgCode(part.fPdgCode),
110 fPdgCode2(part.fPdgCode2),
111 fPartName(part.fPartName),
112 fCollisionSystem(part.fCollisionSystem),
113 fSqrtS(part.fSqrtS),
114 fYield(part.fYield),
115 fStatError(part.fStatError),
116 fSystError(part.fSystError),
117 fNormErrorPos(part.fNormErrorPos),
118 fNormErrorNeg(part.fNormErrorNeg),
119 fYMin(part.fYMin),
120 fYMax(part.fYMax),
121 fStatus(part.fStatus),
122 fMeasurementType(part.fMeasurementType),
123 fCentr(part.fCentr),
124 fIsSum(part.fIsSum),
125 fTag(part.fTag){
126   // Copy constructor
127 }
128
129 AliParticleYield::~AliParticleYield() {
130
131 }
132
133 TTree * AliParticleYield::ReadFromASCIIFileAsTree(const char * fileName, const char * separators){
134   // Read the table from an ASCII File and returns a tree of particles. See ReadFromASCIIFile for detailed info on the format
135   TClonesArray * arr = ReadFromASCIIFile(fileName, separators);
136   AliParticleYield * part = 0;
137   TTree * tree = new TTree ("treePart", "Particle Yields and Ratios");
138   tree->Branch("particles", &part);
139   TIter iterPart (arr);
140   while ((part = (AliParticleYield*) iterPart.Next())){
141     tree->Fill();
142   }
143   
144   delete arr;
145   delete part;
146   return tree;
147 }
148
149 TClonesArray * AliParticleYield::GetEntriesMatchingSelection(TTree * tree, TString selection) {
150   // Returns an array of particles from a tree created with ReadFromASCIIFileAsTree matching the selection. You can use the normal tree sintax for the selection, e.g. "fCentr == \"V0M0010\" && fStatus == 0".
151
152   TClonesArray * arr = new TClonesArray("AliParticleYield");
153   AliParticleYield * part = 0;
154   tree->SetBranchAddress("particles", &part);
155   // In order to get the array, we first create an entry list matching the selection using TTree::Draw, and them we loop over all entries in the tree.
156   tree->Draw(">>particlelist", selection);// Produce selection list
157   TEventList *elist = (TEventList*)gDirectory->Get("particlelist");
158   Int_t npart = elist->GetN();
159   for(Int_t ipart = 0; ipart < npart; ipart++){
160     tree->GetEntry(elist->GetEntry(ipart));
161     new((*arr)[ipart]) AliParticleYield(*part);// We need to clone part, because it is overwritten by the next read
162   }
163   elist->Delete();
164   return arr;
165 }
166
167
168 TClonesArray * AliParticleYield::ReadFromASCIIFile(const char * fileName, const char * separators){
169   // Read the table from an ASCII File with the format indicated
170   // below. Returns a TClonesArray of AliParticleyields with the
171   // content of the lines. Lines beginning by "#" are skipped.
172   // The order of the columns is compulsory, but the separator can be set (by default whitespaces are assumed).
173
174   // The columns should be:
175   // PDG   NAME  SYSTEM  SQRTS  VALUE  SYST  STAT  NORM  YMIN  YMAX  STATUS  TYPE  CENTR     ISSUM TAG
176
177   // PDG should either be an integher or the ratio of two integers (in case of particle ratios), with the following format:
178   //  pdgcode/pdgcode2
179   // NO SPACES ARE ALLOWED IN NAMES AND PDG CODE, unless you use a separator which is not a whitespace
180
181   // A Header can be present (lines beginning with the word "PDG" are also skipped
182
183   const Int_t kNCols = 14; // The lines are actually 15, but the last one (TAG) can be empty, so we put 14 here.
184
185   TClonesArray * arr = new TClonesArray ("AliParticleYield");
186   ifstream filestream (fileName);
187   if(!filestream.is_open()) {
188     Printf("Cannot open file %s\n", fileName);
189     exit(1);
190   }
191   TString line;
192   Int_t ipart = 0;
193   std::cout << "Reading " << fileName << std::endl;
194   
195   while (line.ReadLine(filestream) ) {
196     // Strip trailing and leading whitespaces
197     line = line.Strip(TString::kLeading,  ' ');
198     line = line.Strip(TString::kTrailing, ' ');
199
200     // Skip commented lines and headers
201     if (line.BeginsWith("#")) {
202       //print comments. It if they look like warnings, print them such that they are really visible
203       if(line.Contains("warn", TString::kIgnoreCase)) std::cout << std::endl << "********************************************************" <<std::endl ;
204       std::cout << " " << line.Data() << std::endl;      
205       if(line.Contains("warn", TString::kIgnoreCase)) std::cout << "********************************************************" <<std::endl << std::endl;
206
207       continue;
208     }
209     if (line.BeginsWith("PDG")) continue;
210
211     // Tokenize line using custom separator
212     TObjArray * cols = line.Tokenize(separators);
213
214     // Check the number of columns
215     if(cols->GetEntries() < kNCols) {
216       Printf("Wrong number of columns in table %d vs %d expected" , cols->GetEntries(), kNCols);
217       delete arr;
218       return NULL;
219     }
220
221     // Get Values
222     // get type first, as some operations are type-specific
223     UInt_t  type   = ((TObjString*)cols->At(11)) ->String().Atoi();
224
225     // if it's a ratio, try to get the 2 pdg codes
226     Int_t pdg =0, pdg2 = 0;
227     
228     if (type & kTypeParticleRatio) {      
229       TString col0 = ((TObjString*)cols->At(0))  ->String();
230       TObjArray * tokens = col0.Tokenize("/");
231       if(tokens->GetEntries() != 2) {
232         Printf("ERROR: Cannot get both PDGs for ratios");               
233       } else {
234         pdg  = ((TObjString*)tokens->At(0))  ->String().Atoi();
235         pdg2 = ((TObjString*)tokens->At(1))  ->String().Atoi();
236       }
237     }
238     else {
239       pdg    = ((TObjString*)cols->At(0))  ->String().Atoi();
240     }
241     TString name   = ((TObjString*)cols->At(1))  ->String();
242     Int_t   system = ((TObjString*)cols->At(2))  ->String().Atoi();
243     Float_t sqrts  = ((TObjString*)cols->At(3))  ->String().Atof();
244     Float_t yield  = ((TObjString*)cols->At(4))  ->String().Atof();
245     // The "GetError" function can handle % errors. 
246     Float_t stat   = GetError(((TObjString*)cols->At(5))  ->String(), yield);
247     Float_t syst   = GetError(((TObjString*)cols->At(6))  ->String(), yield);
248     TString normString(((TObjString*)cols->At(7))->String());
249
250     Float_t normPos = 0;
251     Float_t normNeg = 0;
252     if (normString.Contains("+") && normString.Contains("-")) {
253       
254       // If the string for the normalization uncertainty contains a + and a -, it means it is asymmetric
255       if(normString.First("+") < normString.First("-") ) {// the + error is quoted first
256         normPos = GetError(normString(1,normString.First("-")-1)+normString(normString.First("e"),normString.Length()), yield); // start from 1 (skip + sign). The second bit is to propagate the scientific notation to the first part of the error
257         normNeg = GetError(normString(normString.First("-")+1,normString.Length()), yield); // +1 -> skip sign
258       } 
259       else {
260         // This is the opposite case
261         normNeg = GetError(normString(1,normString.First("+")-1)+normString(normString.First("e"),normString.Length()), yield); // start from 1 (skip + sign). The second bit is to propagate the scientific notation to the first part of the error
262         normPos = GetError(normString(normString.First("+")+1,normString.Length()), yield); // +1 -> skip sign
263       }
264       
265     } else {
266       // symmetric error: set only normpos
267       normPos   = GetError(((TObjString*)cols->At(7))  ->String(), yield);
268     }
269     Float_t ymin   = ((TObjString*)cols->At(8))  ->String().Atof();
270     Float_t ymax   = ((TObjString*)cols->At(9))  ->String().Atof();
271     Int_t   status = ((TObjString*)cols->At(10)) ->String().Atoi();
272     TString centr  = ((TObjString*)cols->At(12)) ->String();
273     Int_t   issum  = ((TObjString*)cols->At(13)) ->String().Atoi();    
274     TString tag    = cols->At(14) ? ((TObjString*)cols->At(14)) ->String() : ""; // tag can be empty
275
276     // Cleanup strings
277     name  = name.Strip(TString::kLeading,  ' ');
278     name  = name.Strip(TString::kTrailing, ' ');
279     centr = centr.Strip(TString::kLeading,  ' ');
280     centr = centr.Strip(TString::kTrailing, ' ');
281     tag   = tag.Strip(TString::kLeading,  ' ');
282     tag   = tag.Strip(TString::kTrailing, ' ');
283     
284     // add to array
285     AliParticleYield * part = new  AliParticleYield(pdg,system,sqrts,yield,stat,syst,normPos, normNeg,ymin,ymax,status,type,centr,issum,tag);
286     part->SetPartName(name); // Check name and PDG code consistency   
287     part->SetPdgCode2(pdg2); // Set second PDG code in case of ratios 
288     part->CheckTypeConsistency();                                     
289     if(!part->CheckForDuplicates(arr)) {
290       new ((*arr)[ipart++]) AliParticleYield(*part); 
291     }
292       //    delete part;
293
294   }
295   std::cout << "<- File read" << std::endl;
296
297
298   return arr;
299 }
300
301 const char * AliParticleYield::GetLatexName(Int_t pdg) const {
302   
303   // Returns a TLatex compatible name for the particle
304   // if pdg == 0 uses fPdgcode;
305   // We need the pdg argument for particle ratios
306
307   if(!pdg && fMeasurementType & kTypeParticleRatio) {
308     // If it's a ratio, we try to build the ratio name. To avoid an infinite loop we have to call GetLatexname with a non-zero argument.
309     static TString name;
310     name ="";
311     name += GetLatexName(fPdgCode);
312     name += " / ";
313     name += GetLatexName(fPdgCode2);
314     return name.Data();
315   }
316
317   if(!pdg) pdg = fPdgCode;
318
319   switch (pdg) {
320   case 211:
321     if (fIsSum) return "(#pi^{+} + #pi^{-})";
322     return "#pi^{+}";
323   case -211:
324     return "#pi^{-}";
325   case 321:
326     if (fIsSum) return "(K^{+} + K^{-})";
327     return "K^{+}";
328   case -321:
329     return "K^{-}";
330   case 2212:
331     if (fIsSum) return "(p + #bar{p})";
332     return "p";
333   case -2212:
334     return "#bar^{p}";
335   case 3122:
336     if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
337     return "#Omega^{-}";
338   case -3122:
339     return "#bar{#Omega}^{+}";
340   case 3312:
341     if (fIsSum) return "(#Xi^{-} + #bar{#Xi}^{+})";
342     return "#Xi^{-}";
343   case -3312:
344     return "#bar{#Xi}^{+}";
345   case 3334:
346     if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
347     return "#Omega^{-}";
348   case -3334:
349     return "#bar{#Omega}^{+}";
350   case 310:
351     return "K^{0}_{S}";
352   case 333:
353     return "#phi";
354   case 313:
355     return "K^{*}";
356   case 1000010020:
357     if(fIsSum) return "(d + #bar{d})";
358     return "d";// Deuteron
359   case -1000010020:
360     return "#bar{d}";// Deuteron
361   case 1000020030:
362     if(fIsSum) return "(^{3}He + #bar{^{3}He})";
363     return "^{3}He";
364   case -1000020030:
365     return "#bar{^{3}He}";
366   case 1010010030:
367     if(fIsSum) return "^{3}_{#Lambda}H + #bar{^{3}_{#Lambda}H}";
368     return "^{3}_{#Lambda}H";
369   case -1010010030:
370     return "#bar{^{3}_{#Lambda}H}";    
371   }
372   AliWarning("Latex Name not know for this particle");
373   return fPartName.Data();
374
375 }
376
377 Float_t AliParticleYield::GetTotalError(Bool_t includeNormalization) const {
378   // Returns the total error, including or not the normalization uncertainty
379   // All uncertainties are supposed to be uncorrelated (e.g. summed in quadrature)
380   // If stat and syst are stored separately, the total error is computed summing them in quadrature
381   Float_t error = GetSystError();
382   if (!(fMeasurementType & kTypeOnlyTotError)) error = TMath::Sqrt(error*error + GetStatError()*GetStatError());
383   if(includeNormalization) error = TMath::Sqrt(error*error + GetNormError()*GetNormError());
384   
385   return error;
386
387
388 }
389
390
391 void AliParticleYield::SaveAsASCIIFile(TClonesArray * arr, const char * fileName, const char * separator, Int_t colWidth){
392   // Saves the array as an ASCII File with the format indicated
393   // below. 
394
395   // The columns should be:
396   // PDG   NAME  SYSTEM  SQRTS  VALUE  STAT  SYST  NORM  YMIN  YMAX  STATUS  TYPE  CENTR     ISSUM TAG
397   if(!arr) {
398     Printf("<AliParticleYield::SaveAsASCIIFile> Error: no array provided");
399     return;
400   }
401   if(!fileName) {
402     Printf("<AliParticleYield::SaveAsASCIIFile> Error: no filename provided");
403   }
404
405
406   ofstream fileOut(fileName);
407   //print header
408   fileOut << FormatCol("PDG", colWidth, separator) <<   FormatCol("NAME", colWidth, separator) <<  FormatCol("SYSTEM", colWidth, separator) <<  FormatCol("SQRTS", colWidth, separator) <<  FormatCol("VALUE", colWidth, separator) <<  FormatCol("STAT" , colWidth, separator)<<  FormatCol("SYST", colWidth, separator) <<  FormatCol("NORM", colWidth, separator) <<  FormatCol("YMIN", colWidth, separator) <<  FormatCol("YMAX", colWidth, separator) <<  FormatCol("STATUS", colWidth, separator) <<  FormatCol("TYPE", colWidth, separator) <<  FormatCol("CENTR", colWidth, separator) <<     FormatCol("ISSUM", colWidth, separator) <<  FormatCol("TAG", colWidth, separator) << endl;
409   
410
411   // This is used for float numbers in the table.
412   // The "g" options switches between the normal or scientific notation, whathever is more appropriate.
413   // We want to have up to fSignificantDigits digits after the .
414   char format[20];
415   snprintf(format,20,"%%%dg", fSignificantDigits);
416
417   char formatA[30];// We have to rebuild the format for asymmetric uncertainties...
418   snprintf(formatA,30,"+%%%dg-%%%dg", fSignificantDigits, fSignificantDigits);
419
420   TIter iter(arr);
421   AliParticleYield * part = 0;
422   TString normError ;
423   while ((part = (AliParticleYield*) iter.Next())){    
424     if(part->GetNormErrorNeg()) {
425       normError = FormatCol(Form(formatA, // Asymmetric error format  
426                                  RoundToSignificantFigures(part->GetNormErrorPos(),fSignificantDigits), 
427                                  RoundToSignificantFigures(part->GetNormErrorNeg(),fSignificantDigits)),
428                             colWidth,
429                             separator);
430     }
431     else {
432       normError = FormatCol(Form(format, RoundToSignificantFigures(part->GetNormError(),fSignificantDigits)) , colWidth , separator);
433     }
434     fileOut 
435       << FormatCol(Form("%d",part->GetPdgCode())                                                    , colWidth , separator) 
436       << FormatCol(part->GetPartName()                                                              , colWidth , separator)         
437       << FormatCol(Form("%d", part->GetCollisionSystem())                                           , colWidth , separator) 
438       << FormatCol(Form(format, part->GetSqrtS())                                                   , colWidth , separator)         
439       << FormatCol(Form(format, RoundToSignificantFigures(part->GetYield(),    fSignificantDigits)) , colWidth , separator)
440       << FormatCol(Form(format, RoundToSignificantFigures(part->GetStatError(),fSignificantDigits)) , colWidth , separator) 
441       << FormatCol(Form(format, RoundToSignificantFigures(part->GetSystError(),fSignificantDigits)) , colWidth , separator)
442       << normError.Data()       
443       << FormatCol(Form(format, part->GetYMin())                                                    , colWidth , separator) 
444       << FormatCol(Form(format, part->GetYMax())                                                    , colWidth , separator)         
445       << FormatCol(Form("%d",part->GetStatus()          )                                           , colWidth , separator) 
446       << FormatCol(Form("%d",part->GetMeasurementType() )                                           , colWidth , separator)       
447       << FormatCol(part->GetCentr()                                                                 , colWidth , separator) 
448       << FormatCol(Form("%d",part->GetIsSum())                                                      , colWidth , separator) 
449       << FormatCol(part->GetTag()                                                                   , colWidth , separator) 
450       << endl;
451   }
452
453
454 }
455
456 void AliParticleYield::WriteThermusFile(TClonesArray * arr, const char * filename, Int_t colwidth) {
457   // Writes a txt file which can we used as input in therums fits
458
459   if(!arr) {
460     Printf("<AliParticleYield::WriteThermusFile> Error: no array provided");
461     return;
462   }
463   if(!filename) {
464     Printf("<AliParticleYield::WriteThermusFile> Error: no filename provided");
465     return;
466   }
467
468   ofstream fileOut(filename);
469
470   TIter iter(arr);
471   AliParticleYield * part = 0;
472   char format[20];
473   // This is used for float numbers in the table.
474   // The "g" options switches between the normal or scientific notation, whathever is more appropriate.
475   // We want to have up to fSignificantDigits digits after the .
476   snprintf(format,20,"%%%dg", fSignificantDigits);
477   
478   //  snprintf(format, 20, "%d.%d%%f", fSignificantDigits, fSignificantDigits);
479   while ((part = (AliParticleYield*) iter.Next())){    
480     
481     if(part->IsTypeRatio()) { 
482       // If it's a ratio we have to write the 2 pdg codes
483       fileOut << FormatCol(Form("%d\t%d\t",part->GetPdgCode(), part->GetPdgCode2())                              , colwidth) 
484               << part->GetTag() << "\t"
485               << Form(format, RoundToSignificantFigures(part->GetYield()      , fSignificantDigits)) << "\t"
486               << Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits)) 
487               << endl;
488     }
489     else {
490       fileOut <<Form("%d",part->GetPdgCode())                                                       << "\t" 
491               <<part->GetTag()                                                                      << "\t"
492               <<Form(format, RoundToSignificantFigures(part->GetYield()      , fSignificantDigits)) << "\t"
493               <<Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits)) 
494               << endl;      
495     }
496   
497   }
498   
499 }
500
501
502 const char * AliParticleYield::FormatCol(const char * text, Int_t width,  const char * sep) {
503   
504   TString format(Form("%%-%ds %s", width, sep));
505   return Form(format.Data(), text);
506
507 }
508
509 Double_t AliParticleYield::RoundToSignificantFigures(double num, int n) {
510   // Rounds num to n significant digits.
511   // Recipe from :http://stackoverflow.com/questions/202302/rounding-to-an-arbitrary-number-of-significant-digits
512   // Basically the log is used to determine the number of leading 0s, than convert to an integer by multipliing by the expo, 
513   // round the integer and shift back.
514   if(num == 0) {
515     return 0;
516   }
517
518   Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num));
519   Int_t power = n - (int) d;
520
521   Double_t magnitude = TMath::Power(10, power);
522   Long_t shifted = TMath::Nint(num*magnitude);
523   return shifted/magnitude;
524
525 }
526
527
528 Float_t AliParticleYield::GetError(TString error, Float_t yield) {
529   // The "GetError" function can handle % errors. 
530   if(error.Contains("%")) {
531     return yield * error.Atof()/100;
532   }
533   return error.Atof();
534 }
535
536 void AliParticleYield::SetPdgCode(TString partName) {
537   // Set pdg code from part name, if there was a previous name, check if it is consistent
538   TParticlePDG * part  = TDatabasePDG::Instance()->GetParticle(partName);
539   if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums
540   if(!part) {
541     AliError(Form("No particle %s in TDatabasePDG", partName.Data()));
542     return;
543   }
544   if(fPdgCode != part->PdgCode() &&  !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on PDG codes is disabled case of ratios
545     AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), partName.Data(), part->PdgCode(), fPdgCode));
546   }
547   fPdgCode = part->PdgCode();
548
549 }
550
551 void AliParticleYield::SetPartName(Int_t pdgCode) {
552   // Set part name from pdg code, if there was a previous code, check if it is consistent
553   TParticlePDG * part  = TDatabasePDG::Instance()->GetParticle(pdgCode);
554   if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums
555   if(!part) {
556     AliError(Form("No particle with code %d in TDatabasePDG", pdgCode));
557     return;
558   }
559   if(fPdgCode != part->PdgCode() && !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on particle names is disabled case of ratios
560     AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), part->GetName(), part->PdgCode(), fPdgCode));
561   }
562   fPartName = part->GetName();
563
564 }
565
566 Bool_t AliParticleYield::CheckTypeConsistency() const {
567   // Check for any inconsistency with the type mask. Returns true if the object is fully consistent.
568   Bool_t isOK = kTRUE;
569
570   if((fMeasurementType & kTypeOnlyTotError) && GetStatError()) {
571     AliError(Form("Error flagged as total only, but stat error is not 0 (%f)!",GetStatError()));
572     isOK= kFALSE;
573   } else if (!(fMeasurementType & kTypeOnlyTotError) && (!GetStatError() || !GetSystError())) {
574     AliError("Stat or syst errors are null");
575   } 
576
577   return isOK;
578 }
579
580 void AliParticleYield::Print (Option_t *) const {
581   // Print
582   Printf("-------------------------------");
583   CheckTypeConsistency();
584   Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fIsSum ? "particle + antiparticle" : "", fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
585   Printf("Status: %s, %s", kStatusString[fStatus], fMeasurementType & kTypeLinearInterpolation ? "Interpolated" : "Measured");
586   Printf("%s , sqrt(s) = %2.2f GeV, %2.2f < y < %2.2f %s", kSystemString[fCollisionSystem], fSqrtS, fYMin, fYMax, fCentr.Data()); 
587   if(fMeasurementType & kTypeOnlyTotError) {
588     Printf("%f +- %f (total error)", fYield, fSystError);
589   }
590   else {
591     Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError);
592   }
593   if(fNormErrorNeg) {
594     Printf("Normalizaion uncertainty: +%f-%f", fNormErrorPos, fNormErrorNeg);    
595   }
596   else {
597     Printf("Normalizaion uncertainty: %f", fNormErrorPos);
598   }
599 }
600
601 Bool_t AliParticleYield::operator==(const AliParticleYield& rhs) {
602   // Check if the two particles are identical
603
604   Bool_t isEqual = 
605     (fPdgCode         == rhs.fPdgCode              ) &&
606     (fPdgCode2        == rhs.fPdgCode2             ) &&
607     (fPartName        == rhs.fPartName             ) &&
608     (fCollisionSystem == rhs.fCollisionSystem      ) &&
609     Compare2Floats(fSqrtS,rhs.fSqrtS               ) &&
610     Compare2Floats(fYield,rhs.fYield               ) &&
611     Compare2Floats(fStatError,rhs.fStatError       ) &&
612     Compare2Floats(fSystError,rhs.fSystError       ) &&
613     Compare2Floats(fNormErrorPos,rhs.fNormErrorPos ) &&
614     Compare2Floats(fNormErrorNeg,rhs.fNormErrorNeg ) &&
615     Compare2Floats(fYMin,rhs.fYMin                 ) &&
616     Compare2Floats(fYMax,rhs.fYMax                 ) &&
617     (fStatus          == rhs.fStatus               ) &&
618     (fMeasurementType == rhs.fMeasurementType      ) &&
619     (fCentr           == rhs.fCentr                ) &&
620     (fIsSum           == rhs.fIsSum                ) &&
621     (fTag             == rhs.fTag                  ) ;
622   
623   return isEqual;
624   
625 }
626 Bool_t AliParticleYield::IsTheSameMeasurement(AliParticleYield &rhs) {
627
628   // Check the two particles represent the same measurement (independently of the values)
629   Bool_t isEqual = 
630     (fPdgCode         == rhs.fPdgCode         ) &&
631     (fPdgCode2        == rhs.fPdgCode2        ) &&
632     (fCollisionSystem == rhs.fCollisionSystem ) &&
633     Compare2Floats(fSqrtS,rhs.fSqrtS          ) &&
634     Compare2Floats(fYMin,rhs.fYMin            ) &&
635     Compare2Floats(fYMax,rhs.fYMax            ) &&
636     (fStatus          == rhs.fStatus          ) &&
637     (fCentr           == rhs.fCentr           ) &&
638     (fIsSum           == rhs.fIsSum           ) &&
639     (fTag             == rhs.fTag             ) ;
640   
641   return isEqual;
642
643
644 }
645
646 Bool_t AliParticleYield::CheckForDuplicates(TClonesArray * arr) {
647
648   // loop over all elements on the array and check for duplicates
649   TIter iter(arr);
650   AliParticleYield * part = 0;
651   Bool_t isDuplicate = kFALSE;
652
653   while ((part = (AliParticleYield*) iter.Next())) {
654     if (IsTheSameMeasurement(*part)){
655       AliWarning("Duplicated measurement found");
656       isDuplicate = kTRUE;
657       if (!((*this) == (*part))) {
658         part->Print();
659         Print();
660         AliFatal("The 2 particles are different!");
661       }
662     }
663   }
664   return isDuplicate;
665
666
667  
668 Bool_t AliParticleYield::Compare2Floats(Float_t a, Float_t b) {
669   // just a simple helper for the comparison methods
670   if(!a) {
671     if(!b) return kTRUE; // They are both 0;
672     return kFALSE;// return here to avoid division by 0
673   }
674   Bool_t areEqual = (TMath::Abs((a - b)/a) < fEpsilon); // If the relative difference is < epsilon, returns true
675   if(!areEqual) {
676     Printf("Warning: %f and %f are different", a,b); 
677   }
678   return areEqual;
679 }
680
681
682 Float_t AliParticleYield::GetNormError() const {
683   // Returs a symmetrized error in case the normalizatione Error is asymmetric
684   if(fNormErrorNeg) {
685     AliWarning("Error is asymmetric, returining symmetrized uncertainty");
686     return (TMath::Abs(fNormErrorNeg)+TMath::Abs(fNormErrorPos))/2;
687   }
688   else return fNormErrorPos; // If the uncertainty is not asymmetric, fNormErrorPos stores it.
689
690 }