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