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