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 ="#frac{";
 311     name += GetLatexName(fPdgCode);
 312     name += "}{";
 313     name += GetLatexName(fPdgCode2);
 314     name += "}";
 315     return name.Data();
 316   }
 317
 318   if(!pdg) pdg = fPdgCode;
 319
 320   switch (pdg) {
 321   case 211:
 322     if (fIsSum) return "(#pi^{+} + #pi^{-})";
 323     return "#pi^{+}";
 324   case -211:
 325     return "#pi^{-}";
 326   case 321:
 327     if (fIsSum) return "(K^{+} + K^{-})";
 328     return "K^{+}";
 329   case -321:
 330     return "K^{-}";
 331   case 2212:
 332     if (fIsSum) return "(p + #bar{p})";
 333     return "p";
 334   case -2212:
 335     return "#bar^{p}";
 336   case 3122:
 337     if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
 338     return "#Omega^{-}";
 339   case -3122:
 340     return "#bar{#Omega}^{+}";
 341   case 3312:
 342     if (fIsSum) return "(#Xi^{-} + #bar{#Xi}^{+})";
 343     return "#Xi^{-}";
 344   case -3312:
 345     return "#bar{#Xi}^{+}";
 346   case 3334:
 347     if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
 348     return "#Omega^{-}";
 349   case -3334:
 350     return "#bar{#Omega}^{+}";
 351   case 310:
 352     return "K^{0}_{S}";
 353   case 333:
 354     return "#phi";
 355   case 313:
 356     return "K^{*}";
 357   case 1000010020:
 358     if(fIsSum) return "(d + #bar{d})";
 359     return "d";// Deuteron
 360   case -1000010020:
 361     return "#bar{d}";// Deuteron
 362   case 1000020030:
 363     if(fIsSum) return "(^{3}He + #bar{^{3}He})";
 364     return "^{3}He";
 365   case -1000020030:
 366     return "#bar{^{3}He}";
 367   case 1010010030:
 368     if(fIsSum) return "^{3}_{#Lambda}H + #bar{^{3}_{#Lambda}H}";
 369     return "^{3}_{#Lambda}H";
 370   case -1010010030:
 371     return "#bar{^{3}_{#Lambda}H}";
 372   }
 373   AliWarning("Latex Name not know for this particle");
 374   return fPartName.Data();
 375
 376 }
 377
 378 Float_t AliParticleYield::GetTotalError(Bool_t includeNormalization) const {
 379   // Returns the total error, including or not the normalization uncertainty
 380   // All uncertainties are supposed to be uncorrelated (e.g. summed in quadrature)
 381   // If stat and syst are stored separately, the total error is computed summing them in quadrature
 382   Float_t error = GetSystError();
 383   if (!(fMeasurementType & kTypeOnlyTotError)) error = TMath::Sqrt(error*error + GetStatError()*GetStatError());
 384   if(includeNormalization) error = TMath::Sqrt(error*error + GetNormError()*GetNormError());
 385
 386   return error;
 387
 388
 389 }
 390
 391
 392 void AliParticleYield::SaveAsASCIIFile(TClonesArray * arr, const char * fileName, const char * separator, Int_t colWidth){
 393   // Saves the array as an ASCII File with the format indicated
 394   // below.
 395
 396   // The columns should be:
 397   // PDG   NAME  SYSTEM  SQRTS  VALUE  STAT  SYST  NORM  YMIN  YMAX  STATUS  TYPE  CENTR     ISSUM TAG
 398   if(!arr) {
 399     Printf("<AliParticleYield::SaveAsASCIIFile> Error: no array provided");
 400     return;
 401   }
 402   if(!fileName) {
 403     Printf("<AliParticleYield::SaveAsASCIIFile> Error: no filename provided");
 404   }
 405
 406
 407   ofstream fileOut(fileName);
 408   //print header
 409   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;
 410
 411
 412   // This is used for float numbers in the table.
 413   // The "g" options switches between the normal or scientific notation, whathever is more appropriate.
 414   // We want to have up to fSignificantDigits digits after the .
 415   char format[20];
 416   snprintf(format,20,"%%%dg", fSignificantDigits);
 417
 418   char formatA[30];// We have to rebuild the format for asymmetric uncertainties...
 419   snprintf(formatA,30,"+%%%dg-%%%dg", fSignificantDigits, fSignificantDigits);
 420
 421   TIter iter(arr);
 422   AliParticleYield * part = 0;
 423   TString normError ;
 424   while ((part = (AliParticleYield*) iter.Next())){
 425     if(part->GetNormErrorNeg()) {
 426       normError = FormatCol(Form(formatA, // Asymmetric error format
 427                                  RoundToSignificantFigures(part->GetNormErrorPos(),fSignificantDigits),
 428                                  RoundToSignificantFigures(part->GetNormErrorNeg(),fSignificantDigits)),
 429                             colWidth,
 430                             separator);
 431     }
 432     else {
 433       normError = FormatCol(Form(format, RoundToSignificantFigures(part->GetNormError(),fSignificantDigits)) , colWidth , separator);
 434     }
 435     fileOut
 436       << FormatCol(Form("%d",part->GetPdgCode())                                                    , colWidth , separator)
 437       << FormatCol(part->GetPartName()                                                              , colWidth , separator)
 438       << FormatCol(Form("%d", part->GetCollisionSystem())                                           , colWidth , separator)
 439       << FormatCol(Form(format, part->GetSqrtS())                                                   , colWidth , separator)
 440       << FormatCol(Form(format, RoundToSignificantFigures(part->GetYield(),    fSignificantDigits)) , colWidth , separator)
 441       << FormatCol(Form(format, RoundToSignificantFigures(part->GetStatError(),fSignificantDigits)) , colWidth , separator)
 442       << FormatCol(Form(format, RoundToSignificantFigures(part->GetSystError(),fSignificantDigits)) , colWidth , separator)
 443       << normError.Data()
 444       << FormatCol(Form(format, part->GetYMin())                                                    , colWidth , separator)
 445       << FormatCol(Form(format, part->GetYMax())                                                    , colWidth , separator)
 446       << FormatCol(Form("%d",part->GetStatus()          )                                           , colWidth , separator)
 447       << FormatCol(Form("%d",part->GetMeasurementType() )                                           , colWidth , separator)
 448       << FormatCol(part->GetCentr()                                                                 , colWidth , separator)
 449       << FormatCol(Form("%d",part->GetIsSum())                                                      , colWidth , separator)
 450       << FormatCol(part->GetTag()                                                                   , colWidth , separator)
 451       << endl;
 452   }
 453
 454
 455 }
 456
 457 void AliParticleYield::WriteThermusFile(TClonesArray * arr, const char * filename, Int_t colwidth) {
 458   // Writes a txt file which can we used as input in therums fits
 459
 460   if(!arr) {
 461     Printf("<AliParticleYield::WriteThermusFile> Error: no array provided");
 462     return;
 463   }
 464   if(!filename) {
 465     Printf("<AliParticleYield::WriteThermusFile> Error: no filename provided");
 466     return;
 467   }
 468
 469   ofstream fileOut(filename);
 470
 471   TIter iter(arr);
 472   AliParticleYield * part = 0;
 473   char format[20];
 474   // This is used for float numbers in the table.
 475   // The "g" options switches between the normal or scientific notation, whathever is more appropriate.
 476   // We want to have up to fSignificantDigits digits after the .
 477   snprintf(format,20,"%%%dg", fSignificantDigits);
 478
 479   //  snprintf(format, 20, "%d.%d%%f", fSignificantDigits, fSignificantDigits);
 480   while ((part = (AliParticleYield*) iter.Next())){
 481
 482     if(part->IsTypeRatio()) {
 483       // If it's a ratio we have to write the 2 pdg codes
 484       fileOut << FormatCol(Form("%d\t%d\t",part->GetPdgCode(), part->GetPdgCode2())                              , colwidth)
 485               << part->GetTag() << "\t"
 486               << Form(format, RoundToSignificantFigures(part->GetYield()      , fSignificantDigits)) << "\t"
 487               << Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits))
 488               << endl;
 489     }
 490     else {
 491       fileOut <<Form("%d",part->GetPdgCode())                                                       << "\t"
 492               <<part->GetTag()                                                                      << "\t"
 493               <<Form(format, RoundToSignificantFigures(part->GetYield()      , fSignificantDigits)) << "\t"
 494               <<Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits))
 495               << endl;
 496     }
 497
 498   }
 499
 500 }
 501
 502
 503 const char * AliParticleYield::FormatCol(const char * text, Int_t width,  const char * sep) {
 504
 505   TString format(Form("%%-%ds %s", width, sep));
 506   return Form(format.Data(), text);
 507
 508 }
 509
 510 Double_t AliParticleYield::RoundToSignificantFigures(double num, int n) {
 511   // Rounds num to n significant digits.
 512   // Recipe from :http://stackoverflow.com/questions/202302/rounding-to-an-arbitrary-number-of-significant-digits
 513   // Basically the log is used to determine the number of leading 0s, than convert to an integer by multipliing by the expo,
 514   // round the integer and shift back.
 515   if(num == 0) {
 516     return 0;
 517   }
 518
 519   Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num));
 520   Int_t power = n - (int) d;
 521
 522   Double_t magnitude = TMath::Power(10, power);
 523   Long_t shifted = TMath::Nint(num*magnitude);
 524   return shifted/magnitude;
 525
 526 }
 527
 528
 529 Float_t AliParticleYield::GetError(TString error, Float_t yield) {
 530   // The "GetError" function can handle % errors.
 531   if(error.Contains("%")) {
 532     return yield * error.Atof()/100;
 533   }
 534   return error.Atof();
 535 }
 536
 537 void AliParticleYield::SetPdgCode(TString partName) {
 538   // Set pdg code from part name, if there was a previous name, check if it is consistent
 539   TParticlePDG * part  = TDatabasePDG::Instance()->GetParticle(partName);
 540   if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums
 541   if(!part) {
 542     AliError(Form("No particle %s in TDatabasePDG", partName.Data()));
 543     return;
 544   }
 545   if(fPdgCode != part->PdgCode() &&  !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on PDG codes is disabled case of ratios
 546     AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), partName.Data(), part->PdgCode(), fPdgCode));
 547   }
 548   fPdgCode = part->PdgCode();
 549
 550 }
 551
 552 void AliParticleYield::SetPartName(Int_t pdgCode) {
 553   // Set part name from pdg code, if there was a previous code, check if it is consistent
 554   TParticlePDG * part  = TDatabasePDG::Instance()->GetParticle(pdgCode);
 555   if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums
 556   if(!part) {
 557     AliError(Form("No particle with code %d in TDatabasePDG", pdgCode));
 558     return;
 559   }
 560   if(fPdgCode != part->PdgCode() && !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on particle names is disabled case of ratios
 561     AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), part->GetName(), part->PdgCode(), fPdgCode));
 562   }
 563   fPartName = part->GetName();
 564
 565 }
 566
 567 Bool_t AliParticleYield::CheckTypeConsistency() const {
 568   // Check for any inconsistency with the type mask. Returns true if the object is fully consistent.
 569   Bool_t isOK = kTRUE;
 570
 571   if((fMeasurementType & kTypeOnlyTotError) && GetStatError()) {
 572     AliError(Form("Error flagged as total only, but stat error is not 0 (%f)!",GetStatError()));
 573     isOK= kFALSE;
 574   } else if (!(fMeasurementType & kTypeOnlyTotError) && (!GetStatError() || !GetSystError())) {
 575     AliError("Stat or syst errors are null");
 576   }
 577
 578   return isOK;
 579 }
 580
 581 void AliParticleYield::Print (Option_t *) const {
 582   // Print
 583   Printf("-------------------------------");
 584   CheckTypeConsistency();
 585   Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fIsSum ? "particle + antiparticle" : "", fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
 586   Printf("Status: %s, %s", kStatusString[fStatus], fMeasurementType & kTypeLinearInterpolation ? "Interpolated" : "Measured");
 587   Printf("%s , sqrt(s) = %2.2f GeV, %2.2f < y < %2.2f %s", kSystemString[fCollisionSystem], fSqrtS, fYMin, fYMax, fCentr.Data());
 588   if(fMeasurementType & kTypeOnlyTotError) {
 589     Printf("%f +- %f (total error)", fYield, fSystError);
 590   }
 591   else {
 592     Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError);
 593   }
 594   if(fNormErrorNeg) {
 595     Printf("Normalization uncertainty: +%f-%f", fNormErrorPos, fNormErrorNeg);
 596   }
 597   else {
 598     Printf("Normalization uncertainty: %f", fNormErrorPos);
 599   }
 600 }
 601
 602 Bool_t AliParticleYield::operator==(const AliParticleYield& rhs) {
 603   // Check if the two particles are identical
 604
 605   Bool_t isEqual =
 606     (fPdgCode         == rhs.fPdgCode              ) &&
 607     (fPdgCode2        == rhs.fPdgCode2             ) &&
 608     (fPartName        == rhs.fPartName             ) &&
 609     (fCollisionSystem == rhs.fCollisionSystem      ) &&
 610     Compare2Floats(fSqrtS,rhs.fSqrtS               ) &&
 611     Compare2Floats(fYield,rhs.fYield               ) &&
 612     Compare2Floats(fStatError,rhs.fStatError       ) &&
 613     Compare2Floats(fSystError,rhs.fSystError       ) &&
 614     Compare2Floats(fNormErrorPos,rhs.fNormErrorPos ) &&
 615     Compare2Floats(fNormErrorNeg,rhs.fNormErrorNeg ) &&
 616     Compare2Floats(fYMin,rhs.fYMin                 ) &&
 617     Compare2Floats(fYMax,rhs.fYMax                 ) &&
 618     (fStatus          == rhs.fStatus               ) &&
 619     (fMeasurementType == rhs.fMeasurementType      ) &&
 620     (fCentr           == rhs.fCentr                ) &&
 621     (fIsSum           == rhs.fIsSum                ) &&
 622     (fTag             == rhs.fTag                  ) ;
 623
 624   return isEqual;
 625
 626 }
 627 Bool_t AliParticleYield::IsTheSameMeasurement(AliParticleYield &rhs) {
 628
 629   // Check the two particles represent the same measurement (independently of the values)
 630   Bool_t isEqual =
 631     (fPdgCode         == rhs.fPdgCode         ) &&
 632     (fPdgCode2        == rhs.fPdgCode2        ) &&
 633     (fCollisionSystem == rhs.fCollisionSystem ) &&
 634     Compare2Floats(fSqrtS,rhs.fSqrtS          ) &&
 635     Compare2Floats(fYMin,rhs.fYMin            ) &&
 636     Compare2Floats(fYMax,rhs.fYMax            ) &&
 637     (fStatus          == rhs.fStatus          ) &&
 638     (fCentr           == rhs.fCentr           ) &&
 639     (fIsSum           == rhs.fIsSum           ) &&
 640     (fTag             == rhs.fTag             ) ;
 641
 642   return isEqual;
 643
 644
 645 }
 646
 647 Bool_t AliParticleYield::CheckForDuplicates(TClonesArray * arr) {
 648
 649   // loop over all elements on the array and check for duplicates
 650   TIter iter(arr);
 651   AliParticleYield * part = 0;
 652   Bool_t isDuplicate = kFALSE;
 653
 654   while ((part = (AliParticleYield*) iter.Next())) {
 655     if (IsTheSameMeasurement(*part)){
 656       AliWarning("Duplicated measurement found");
 657       isDuplicate = kTRUE;
 658       if (!((*this) == (*part))) {
 659         part->Print();
 660         Print();
 661         AliFatal("The 2 particles are different!");
 662       }
 663     }
 664   }
 665   return isDuplicate;
 666
 667 }
 668
 669 Bool_t AliParticleYield::Compare2Floats(Float_t a, Float_t b) {
 670   // just a simple helper for the comparison methods
 671   if(!a) {
 672     if(!b) return kTRUE; // They are both 0;
 673     return kFALSE;// return here to avoid division by 0
 674   }
 675   Bool_t areEqual = (TMath::Abs((a - b)/a) < fEpsilon); // If the relative difference is < epsilon, returns true
 676   if(!areEqual) {
 677     Printf("Warning: %f and %f are different", a,b);
 678   }
 679   return areEqual;
 680 }
 681
 682
 683 Float_t AliParticleYield::GetNormError() const {
 684   // Returs a symmetrized error in case the normalizatione Error is asymmetric
 685   if(fNormErrorNeg) {
 686     AliWarning("Error is asymmetric, returining symmetrized uncertainty");
 687     return (TMath::Abs(fNormErrorNeg)+TMath::Abs(fNormErrorPos))/2;
 688   }
 689   else return fNormErrorPos; // If the uncertainty is not asymmetric, fNormErrorPos stores it.
 690
 691 }
 692
 693 AliParticleYield * AliParticleYield::FindParticle(TClonesArray * arr, Int_t pdg, Int_t system, Float_t sqrts, TString centrality, Int_t isSum, Int_t status, Int_t pdg2){
 694   // Finds a particle in array matching the search criteria. If more than one is found, prints a warning
 695   // If status is -1, tries to return the best (lower status value)
 696   // If pdg2 is not zero, we try to match it as well (we are looking for a ratio)
 697   // The centrality is compared with TString::Contains
 698
 699   TIter iter(arr);
 700   AliParticleYield * part = 0;
 701   AliParticleYield * foundPart = 0;
 702   while ((part = dynamic_cast<AliParticleYield*>(iter.Next()))){
 703     if (part->GetPdgCode() == pdg &&                     // same pdg
 704         part->GetCollisionSystem() == system &&          // same system
 705         Compare2Floats(part->GetSqrtS(), sqrts) &&       // same energy
 706         part->GetCentr().Contains(centrality) &&         // compatible centrality
 707         (part->GetPdgCode2() == pdg2) &&                 // same PDG2, if requested (we are looking for a ratio). We also need to check explicitly for pdg2=0 not to match ratios
 708         (status < 0 || part->GetStatus() == status)  &&  // same status, if requested
 709         (isSum  < 0 || part->GetIsSum()  == isSum)       // part+antipart or not, if requested
 710         ) {
 711       if (foundPart) { // we already found a patching particle
 712         Printf("WARNING<AliParticleYield::FindParticle>: Found another particle matching the same criteria");
 713         foundPart->Print();
 714         part->Print();
 715         if (part->GetStatus() == foundPart->GetStatus()) { // Does it have the same status?! Don't know what to do!
 716           Printf("WARNING<AliParticleYield::FindParticle>: they have the same status, I cannot decide, resetting particle");
 717           foundPart = 0;
 718         }
 719         else if (part->GetStatus()< foundPart->GetStatus()) { // Is it of better quality? select it!
 720           Printf("WARNING<AliParticleYield::FindParticle>: the new one has a smaller status: selecting it!");
 721           foundPart = part;
 722         }
 723       } else { // First match
 724         foundPart = part;
 725       }
 726
 727     }
 728
 729   }
 730   return foundPart;
 731
 732 }
 733
 734 void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield*part2) {
 735   // Combines metadata from part1 and part2
 736   Int_t pdg1 = part1->GetPdgCode();
 737   Int_t pdg2 = pdg1 == part2->GetPdgCode() ? part1->GetPdgCode2() : part2->GetPdgCode();
 738   Int_t   system = part1->GetCollisionSystem() == part2->GetCollisionSystem() ? part2->GetCollisionSystem() : -1;
 739   Float_t sqrts = Compare2Floats(part1->GetSqrtS(), part2->GetSqrtS()) ? part1->GetSqrtS() : 0;
 740   Int_t ymin = part1->GetYMin() == part2->GetYMin() ? part2->GetYMin() : -1000;
 741   Int_t ymax = part1->GetYMax() == part2->GetYMax() ? part2->GetYMax() : -1000;
 742   Int_t status = part1->GetStatus() == part2->GetStatus() ? part2->GetStatus() : -1;
 743   Int_t type = part1->GetMeasurementType() == part2->GetMeasurementType() ? part2->GetMeasurementType() : -1;
 744
 745   TString centr = part1->GetCentr() == part2->GetCentr() ? part2->GetCentr() : part1->GetCentr()+"/"+part2->GetCentr();
 746   TString tag = part1->GetTag() == part2->GetTag() ? part2->GetTag() : part1->GetTag()+"/"+part2->GetTag();
 747   Int_t issum = part1->GetIsSum() == part2->GetIsSum() ? part2->GetIsSum() : -1000;
 748
 749   SetPdgCode(pdg1);
 750   SetPdgCode2(pdg2);
 751   SetCollisionSystem(AliPYCSystem_t(system));
 752   SetSqrtS(sqrts);
 753   SetYMin(ymin);
 754   SetYMax(ymax);
 755   SetStatus(AliPYStatusCode_t(status));
 756   SetMeasurementType(type);
 757   SetCentr(centr);
 758   SetTag(tag);
 759   SetIsSum(issum);
 760
 761 }
 762
 763 AliParticleYield * AliParticleYield::Add   (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt){
 764
 765   // Computes the ratio of 2 particles.
 766   // Valid options:
 767   //  - NQ: Propagates normalization errors quadratically (by default they are propagated linearly)
 768   //  - SL: propagates STATISTICAL errors linearly
 769   //  - YQ: propagates SYSTEMATIC errors quadratically
 770   //  NB by default, statistical errors are propagated quadratically and systematic errors linearly
 771   // if "Correlated error" is non null, it is subtracted in quadrature from the result. It should be a fractional error.
 772
 773   if(!part1 || !part2) {
 774     Printf("WARNING<AliParticleYield::Add>: part1 or part2 is null!");
 775     return 0;
 776   }
 777
 778   TString sopt(opt);
 779   sopt.ToUpper();
 780
 781   Float_t value = part1->GetYield() + part2->GetYield();
 782   Float_t stat  = SumErrors(part1, part2, 0, sopt.Contains("SL") ? "L": "" ); // the option decices if it is propagated linearly pr or quadratically
 783   Float_t syst  = SumErrors(part1, part2, 1, sopt.Contains("YQ") ? "" : "L" );// the option decices if it is propagated linearly pr or quadratically
 784   Float_t norm = SumErrors(part1, part2, 2, sopt.Contains("NQ") ? "" :"L");
 785
 786
 787   if(correlatedError) {
 788     syst = TMath::Sqrt(syst*syst -correlatedError*correlatedError*value*value); // FIXME: this line was never tested
 789   }
 790
 791   AliParticleYield * part = new AliParticleYield();
 792   part->SetYield(value);
 793   part->SetStatError(stat);
 794   part->SetSystError(syst);
 795   part->SetNormError(norm);
 796   part->CombineMetadata(part1, part2);
 797
 798   return part;
 799
 800
 801 }
 802
 803 void AliParticleYield::Scale(Float_t scale) {
 804   // scales the measurement by an errorless number
 805   fYield        *= scale;
 806   fNormErrorNeg *= scale;
 807   fNormErrorPos *= scale;
 808   fStatError    *= scale;
 809   fSystError    *= scale;
 810
 811 }
 812
 813 AliParticleYield * AliParticleYield::Divide (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt) {
 814   // Computes the ratio of 2 particles.
 815   // Valid options:
 816   //  - NQ: assumes normalization errors to be uncorrelated and propagates them quadratically (otherwise the normalization error on the ratio is set to 0
 817   //  - SL: propagates STATISTICAL errors linearly
 818   //  - YQ: propagates SYSTEMATIC errors quadratically
 819   //  NB by default, statistical errors are propagated quadratically and systematic errors linearly
 820   // if "Correlated error" is non null, it is subtracted in quadrature from the result.It should be a fractional error.
 821
 822   if(!part1 || !part2) {
 823     Printf("WARNING<AliParticleYield::Divide>: part1 or part2 is null!");
 824     return 0;
 825   }
 826
 827   TString sopt(opt);
 828   sopt.ToUpper();
 829   if(part1->IsTypeRatio() || part2->IsTypeRatio()){
 830     Printf("WARNING<AliParticleYield::Divide>: If computing a double ratio, some meta info may be not reliable!");
 831   }
 832
 833   Float_t value = part1->GetYield() / part2->GetYield();
 834   // Since in a ratio we propagate a relative error, we have to multiply it back for value in order to get the absolute uncertainty
 835   Float_t stat  = SumErrors(part1, part2, 0, sopt.Contains("SL") ? "RL": "R" ) *value; // R means that it's a relative error, the option decices if it is propagated linearly pr or quadratically
 836   Float_t syst  = SumErrors(part1, part2, 1, sopt.Contains("YQ") ? "R" : "RL" )*value;// R means that it's a relative error, the option decices if it is propagated linearly pr or quadratically
 837   Float_t norm = 0;
 838   if(sopt.Contains("NQ")) {// if opt contains N, propagate the normalization error assuming it is independent
 839     norm  = SumErrors(part1, part2, 2, "R")*value;
 840   }
 841
 842   if(correlatedError) {
 843     syst = TMath::Sqrt(syst/value*syst/value -correlatedError*correlatedError)*value; // FIXME: this line was never tested
 844   }
 845
 846   AliParticleYield * part = new AliParticleYield();
 847   part->SetYield(value);
 848   part->SetStatError(stat);
 849   part->SetSystError(syst);
 850   part->SetNormError(norm);
 851   part->CombineMetadata(part1, part2);
 852
 853   return part;
 854
 855 }
 856
 857 Double_t AliParticleYield::SumErrors(AliParticleYield * part1, AliParticleYield * part2, Int_t error, Option_t * opt) {
 858
 859   // Combines 2 errors.
 860   //  error = 0 -> statistical error
 861   //  error = 1 -> systematic error
 862   //  error = 2 -> normalization error
 863   //  Valid options
 864   //   "R" it propagates it as a relative error, WARNING: it also returns a relative error!
 865   //   "L" it propagates it sums the errors linearly (by default it is done in quadrature)
 866
 867
 868   TString sopt(opt);
 869   sopt.ToUpper();
 870   Bool_t isRelative = sopt.Contains("R");
 871   Bool_t isLinear   = sopt.Contains("L");
 872
 873   Double_t err1 = 0;
 874   Double_t err2 = 0;
 875   Double_t err  = 0;
 876   if (error == 0) {
 877     err1 = part1->GetStatError();
 878     err2 = part2->GetStatError();
 879   } else if (error == 1) {
 880     err1 = part1->GetSystError();
 881     err2 = part2->GetSystError();
 882   } else if (error == 2) {
 883     err1 = part1->GetNormError();
 884     err2 = part2->GetNormError();
 885   } else {
 886     Printf("ERROR<AliParticleYield::SumErrors>: wrong error #:%d", error);
 887   }
 888
 889   if(isRelative) {
 890     err1  /= part1->GetYield();
 891     err2  /= part2->GetYield();
 892   }
 893
 894   if(isLinear) {
 895     err = err1+err2;
 896   } else {
 897     err = TMath::Sqrt(err1*err1 + err2*err2);
 898   }
 899
 900   if(isRelative) return err;
 901
 902   return err;
 903
 904 }