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 #include "TCut.h"
  16
  17 using std::endl;
  18 using std::left;
  19 using std::setw;
  20 using std::ifstream;
  21 using std::ofstream;
  22
  23 ClassImp(AliParticleYield)
  24
  25 // set statics
  26 const char * AliParticleYield::kStatusString[] = {"Published", "Preliminary", "Final, but not published", "May change"} ;
  27 const char * AliParticleYield::kSystemString[] = {"pp", "p-Pb", "Pb-Pb"} ;
  28 Int_t AliParticleYield::fSignificantDigits = 3;
  29 Float_t AliParticleYield::fEpsilon = 0.000000000000000001;
  30
  31 AliParticleYield::AliParticleYield() :
  32 TObject(),
  33 fPdgCode(0),
  34 fPdgCode2(0),
  35 fPartName(""),
  36 fCollisionSystem(0),
  37 fSqrtS(0),
  38 fYield(0),
  39 fStatError(0),
  40 fSystError(0),
  41 fNormErrorPos(0),
  42 fNormErrorNeg(0),
  43 fYMin(0),
  44 fYMax(0),
  45 fStatus(0),
  46 fMeasurementType(0),
  47 fCentr(""),
  48 fIsSum(0),
  49 fTag("")
  50 {
  51   // constructor
  52   AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
  53 }
  54
  55 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):
  56 TObject(),
  57 fPdgCode(pdg),
  58 fPdgCode2(0),
  59 fPartName(""),
  60 fCollisionSystem(system),
  61 fSqrtS(sqrts),
  62 fYield(value),
  63 fStatError(stat),
  64 fSystError(syst),
  65 fNormErrorPos(norm),
  66 fNormErrorNeg(0),
  67 fYMin(ymin),
  68 fYMax(ymax),
  69 fStatus(status),
  70 fMeasurementType(type),
  71 fCentr(centr),
  72 fIsSum(isSum),
  73 fTag(tag)
  74
  75 {
  76   // Constructor
  77   AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
  78   TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdgCode);
  79   if(!part) AliError(Form("No particle with PDG code %d in the database", fPdgCode));
  80   else fPartName = part->GetName();
  81 }
  82
  83 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):
  84 TObject(),
  85 fPdgCode(pdg),
  86 fPdgCode2(0),
  87 fPartName(""),
  88 fCollisionSystem(system),
  89 fSqrtS(sqrts),
  90 fYield(value),
  91 fStatError(stat),
  92 fSystError(syst),
  93 fNormErrorPos(normPos),
  94 fNormErrorNeg(normNeg),
  95 fYMin(ymin),
  96 fYMax(ymax),
  97 fStatus(status),
  98 fMeasurementType(type),
  99 fCentr(centr),
 100 fIsSum(isSum),
 101 fTag(tag)
 102
 103 {
 104   // Constructor
 105   AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
 106   TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdgCode);
 107   if(!part) AliError(Form("No particle with PDG code %d in the database", fPdgCode));
 108   else fPartName = part->GetName();
 109 }
 110
 111
 112 AliParticleYield::AliParticleYield(const AliParticleYield& part) :
 113 TObject(),
 114 fPdgCode(part.fPdgCode),
 115 fPdgCode2(part.fPdgCode2),
 116 fPartName(part.fPartName),
 117 fCollisionSystem(part.fCollisionSystem),
 118 fSqrtS(part.fSqrtS),
 119 fYield(part.fYield),
 120 fStatError(part.fStatError),
 121 fSystError(part.fSystError),
 122 fNormErrorPos(part.fNormErrorPos),
 123 fNormErrorNeg(part.fNormErrorNeg),
 124 fYMin(part.fYMin),
 125 fYMax(part.fYMax),
 126 fStatus(part.fStatus),
 127 fMeasurementType(part.fMeasurementType),
 128 fCentr(part.fCentr),
 129 fIsSum(part.fIsSum),
 130 fTag(part.fTag){
 131   // Copy constructor
 132 }
 133
 134 AliParticleYield::~AliParticleYield() {
 135
 136 }
 137
 138 TTree * AliParticleYield::GetTreeFromArray(TClonesArray * arr) {
 139   // Returns a tree from an array of Tparticles
 140   AliParticleYield * part = 0;
 141   TTree * tree = new TTree ("treePart", "Particle Yields and Ratios");
 142   tree->Branch("particles", &part);
 143   TIter iterPart (arr);
 144   while ((part = (AliParticleYield*) iterPart.Next())){
 145     tree->Fill();
 146   }
 147   if(part) delete part;
 148   return tree;
 149
 150
 151
 152 }
 153
 154
 155 TTree * AliParticleYield::ReadFromASCIIFileAsTree(const char * fileName, const char * separators){
 156   // Read the table from an ASCII File and returns a tree of particles. See ReadFromASCIIFile for detailed info on the format
 157   TClonesArray * arr = ReadFromASCIIFile(fileName, separators);
 158   TTree * tree = GetTreeFromArray(arr);
 159   delete arr;
 160   return tree;
 161 }
 162
 163 TClonesArray * AliParticleYield::GetEntriesMatchingSelection(TTree * tree, TCut selection) {
 164   // 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".
 165
 166   TClonesArray * arr = new TClonesArray("AliParticleYield");
 167   AliParticleYield * part = 0;
 168   tree->SetBranchAddress("particles", &part);
 169   // 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.
 170   tree->Draw(">>particlelist", selection);// Produce selection list
 171   TEventList *elist = (TEventList*)gDirectory->Get("particlelist");
 172   Int_t npart = elist->GetN();
 173   for(Int_t ipart = 0; ipart < npart; ipart++){
 174     tree->GetEntry(elist->GetEntry(ipart));
 175     new((*arr)[ipart]) AliParticleYield(*part);// We need to clone part, because it is overwritten by the next read
 176   }
 177   elist->Delete();
 178   return arr;
 179 }
 180
 181
 182 TClonesArray * AliParticleYield::ReadFromASCIIFile(const char * fileName, const char * separators){
 183   // Read the table from an ASCII File with the format indicated
 184   // below. Returns a TClonesArray of AliParticleyields with the
 185   // content of the lines. Lines beginning by "#" are skipped.
 186   // The order of the columns is compulsory, but the separator can be set (by default whitespaces are assumed).
 187
 188   // The columns should be:
 189   // PDG   NAME  SYSTEM  SQRTS  VALUE  SYST  STAT  NORM  YMIN  YMAX  STATUS  TYPE  CENTR     ISSUM TAG
 190
 191   // PDG should either be an integher or the ratio of two integers (in case of particle ratios), with the following format:
 192   //  pdgcode/pdgcode2
 193   // NO SPACES ARE ALLOWED IN NAMES AND PDG CODE, unless you use a separator which is not a whitespace
 194
 195   // A Header can be present (lines beginning with the word "PDG" are also skipped
 196
 197   const Int_t kNCols = 14; // The lines are actually 15, but the last one (TAG) can be empty, so we put 14 here.
 198
 199   TClonesArray * arr = new TClonesArray ("AliParticleYield");
 200   ifstream filestream (fileName);
 201   if(!filestream.is_open()) {
 202     Printf("Cannot open file %s\n", fileName);
 203     exit(1);
 204   }
 205   TString line;
 206   Int_t ipart = 0;
 207   std::cout << "Reading " << fileName << std::endl;
 208
 209   while (line.ReadLine(filestream) ) {
 210     // Strip trailing and leading whitespaces
 211     line = line.Strip(TString::kLeading,  ' ');
 212     line = line.Strip(TString::kTrailing, ' ');
 213
 214     // Skip commented lines and headers
 215     if (line.BeginsWith("#")) {
 216       //print comments. It if they look like warnings, print them such that they are really visible
 217       if(line.Contains("warn", TString::kIgnoreCase)) std::cout << std::endl << "********************************************************" <<std::endl ;
 218       std::cout << " " << line.Data() << std::endl;
 219       if(line.Contains("warn", TString::kIgnoreCase)) std::cout << "********************************************************" <<std::endl << std::endl;
 220
 221       continue;
 222     }
 223     if (line.BeginsWith("PDG")) continue;
 224
 225     // Tokenize line using custom separator
 226     TObjArray * cols = line.Tokenize(separators);
 227
 228     // Check the number of columns
 229     if(cols->GetEntries() < kNCols) {
 230       Printf("Wrong number of columns in table %d vs %d expected" , cols->GetEntries(), kNCols);
 231       delete arr;
 232       return NULL;
 233     }
 234
 235     // Get Values
 236     // get type first, as some operations are type-specific
 237     UInt_t  type   = ((TObjString*)cols->At(11)) ->String().Atoi();
 238
 239     // if it's a ratio, try to get the 2 pdg codes
 240     Int_t pdg =0, pdg2 = 0;
 241
 242     if (type & kTypeParticleRatio) {
 243       TString col0 = ((TObjString*)cols->At(0))  ->String();
 244       TObjArray * tokens = col0.Tokenize("/");
 245       if(tokens->GetEntries() != 2) {
 246         Printf("ERROR: Cannot get both PDGs for ratios");
 247       } else {
 248         pdg  = ((TObjString*)tokens->At(0))  ->String().Atoi();
 249         pdg2 = ((TObjString*)tokens->At(1))  ->String().Atoi();
 250       }
 251     }
 252     else {
 253       pdg    = ((TObjString*)cols->At(0))  ->String().Atoi();
 254     }
 255     TString name   = ((TObjString*)cols->At(1))  ->String();
 256     Int_t   system = ((TObjString*)cols->At(2))  ->String().Atoi();
 257     Float_t sqrts  = ((TObjString*)cols->At(3))  ->String().Atof();
 258     Float_t yield  = ((TObjString*)cols->At(4))  ->String().Atof();
 259     // The "GetError" function can handle % errors.
 260     Float_t stat   = GetError(((TObjString*)cols->At(5))  ->String(), yield);
 261     Float_t syst   = GetError(((TObjString*)cols->At(6))  ->String(), yield);
 262     TString normString(((TObjString*)cols->At(7))->String());
 263
 264     Float_t normPos = 0;
 265     Float_t normNeg = 0;
 266     if (normString.Contains("+") && normString.Contains("-")) {
 267
 268       // If the string for the normalization uncertainty contains a + and a -, it means it is asymmetric
 269       if(normString.First("+") < normString.First("-") ) {// the + error is quoted first
 270         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
 271         normNeg = GetError(normString(normString.First("-")+1,normString.Length()), yield); // +1 -> skip sign
 272       }
 273       else {
 274         // This is the opposite case
 275         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
 276         normPos = GetError(normString(normString.First("+")+1,normString.Length()), yield); // +1 -> skip sign
 277       }
 278
 279     } else {
 280       // symmetric error: set only normpos
 281       normPos   = GetError(((TObjString*)cols->At(7))  ->String(), yield);
 282     }
 283     Float_t ymin   = ((TObjString*)cols->At(8))  ->String().Atof();
 284     Float_t ymax   = ((TObjString*)cols->At(9))  ->String().Atof();
 285     Int_t   status = ((TObjString*)cols->At(10)) ->String().Atoi();
 286     TString centr  = ((TObjString*)cols->At(12)) ->String();
 287     Int_t   issum  = ((TObjString*)cols->At(13)) ->String().Atoi();
 288     TString tag    = cols->At(14) ? ((TObjString*)cols->At(14)) ->String() : ""; // tag can be empty
 289
 290     // Cleanup strings
 291     name  = name.Strip(TString::kLeading,  ' ');
 292     name  = name.Strip(TString::kTrailing, ' ');
 293     centr = centr.Strip(TString::kLeading,  ' ');
 294     centr = centr.Strip(TString::kTrailing, ' ');
 295     tag   = tag.Strip(TString::kLeading,  ' ');
 296     tag   = tag.Strip(TString::kTrailing, ' ');
 297
 298     // add to array
 299     AliParticleYield * part = new  AliParticleYield(pdg,system,sqrts,yield,stat,syst,normPos, normNeg,ymin,ymax,status,type,centr,issum,tag);
 300     part->SetPartName(name); // Check name and PDG code consistency
 301     part->SetPdgCode2(pdg2); // Set second PDG code in case of ratios
 302     part->CheckTypeConsistency();
 303     if(!part->CheckForDuplicates(arr)) {
 304       new ((*arr)[ipart++]) AliParticleYield(*part);
 305     }
 306       //    delete part;
 307
 308   }
 309   std::cout << "<- File read" << std::endl;
 310
 311
 312   return arr;
 313 }
 314
 315 const char * AliParticleYield::GetLatexName(Int_t pdg) const {
 316
 317   // Returns a TLatex compatible name for the particle
 318   // if pdg == 0 uses fPdgcode;
 319   // We need the pdg argument for particle ratios
 320
 321   if(!pdg && fMeasurementType & kTypeParticleRatio) {
 322     // 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.
 323     static TString name;
 324     name ="#frac{";
 325     name += GetLatexName(fPdgCode);
 326     name += "}{";
 327     name += GetLatexName(fPdgCode2);
 328     name += "}";
 329     return name.Data();
 330   }
 331
 332   if(!pdg) pdg = fPdgCode;
 333
 334   switch (pdg) {
 335   case 211:
 336     if (fIsSum) return "(#pi^{+} + #pi^{-})";
 337     return "#pi^{+}";
 338     break;
 339   case -211:
 340     return "#pi^{-}";
 341     break;
 342   case 321:
 343     if (fIsSum) return "(K^{+} + K^{-})";
 344     return "K^{+}";
 345     break;
 346   case -321:
 347     return "K^{-}";
 348     break;
 349   case 2212:
 350     if (fIsSum) return "(p + #bar{p})";
 351     return "p";
 352     break;
 353   case -2212:
 354     return "#bar^{p}";
 355     break;
 356   case 3122:
 357     if (fIsSum) return "(#Lambda + #bar{#Lambda})";
 358     return "#Lambda";
 359     break;
 360   case -3122:
 361     return "#bar{#Lamnba}";
 362     break;
 363   case 3312:
 364     if (fIsSum) return "(#Xi^{-} + #bar{#Xi}^{+})";
 365     return "#Xi^{-}";
 366     break;
 367   case -3312:
 368     return "#bar{#Xi}^{+}";
 369     break;
 370   case 3334:
 371     if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
 372     return "#Omega^{-}";
 373     break;
 374   case -3334:
 375     return "#bar{#Omega}^{+}";
 376     break;
 377   case 310:
 378     return "K^{0}_{S}";
 379     break;
 380   case 333:
 381     return "#phi";
 382     break;
 383   case 313:
 384     if(fIsSum) return "(K* + #bar{K*})";
 385     return "K*";
 386     break;
 387   case -313:
 388     return "#bar{K*}";
 389     break;
 390   case 1000010020:
 391     if(fIsSum) return "(d + #bar{d})";
 392     return "d";// Deuteron
 393     break;
 394   case -1000010020:
 395     return "#bar{d}";// Deuteron
 396     break;
 397   case 1000020030:
 398     if(fIsSum) return "(^{3}He + #bar{^{3}He})";
 399     return "^{3}He";
 400     break;
 401   case -1000020030:
 402     return "#bar{^{3}He}";
 403     break;
 404   case 1010010030:
 405     if(fIsSum) return "^{3}_{#Lambda}H + #bar{^{3}_{#Lambda}H}";
 406     return "^{3}_{#Lambda}H";
 407     break;
 408   case -1010010030:
 409     return "#bar{^{3}_{#Lambda}H}";
 410     break;
 411   default:
 412     AliWarning("Latex Name not know for this particle");
 413   }
 414
 415   return fPartName.Data();
 416
 417 }
 418
 419 Float_t AliParticleYield::GetTotalError(Bool_t includeNormalization) const {
 420   // Returns the total error, including or not the normalization uncertainty
 421   // All uncertainties are supposed to be uncorrelated (e.g. summed in quadrature)
 422   // If stat and syst are stored separately, the total error is computed summing them in quadrature
 423   Float_t error = GetSystError();
 424   if (!(fMeasurementType & kTypeOnlyTotError)) error = TMath::Sqrt(error*error + GetStatError()*GetStatError());
 425   if(includeNormalization) error = TMath::Sqrt(error*error + GetNormError()*GetNormError());
 426
 427   return error;
 428
 429
 430 }
 431
 432
 433 void AliParticleYield::SaveAsASCIIFile(TClonesArray * arr, const char * fileName, const char * separator, Int_t colWidth){
 434   // Saves the array as an ASCII File with the format indicated
 435   // below.
 436
 437   // The columns should be:
 438   // PDG   NAME  SYSTEM  SQRTS  VALUE  STAT  SYST  NORM  YMIN  YMAX  STATUS  TYPE  CENTR     ISSUM TAG
 439   if(!arr) {
 440     Printf("<AliParticleYield::SaveAsASCIIFile> Error: no array provided");
 441     return;
 442   }
 443   if(!fileName) {
 444     Printf("<AliParticleYield::SaveAsASCIIFile> Error: no filename provided");
 445   }
 446
 447
 448   ofstream fileOut(fileName);
 449   //print header
 450   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;
 451
 452
 453   // This is used for float numbers in the table.
 454   // The "g" options switches between the normal or scientific notation, whathever is more appropriate.
 455   // We want to have up to fSignificantDigits digits after the .
 456   char format[20];
 457   snprintf(format,20,"%%%dg", fSignificantDigits);
 458
 459   char formatA[30];// We have to rebuild the format for asymmetric uncertainties...
 460   snprintf(formatA,30,"+%%%dg-%%%dg", fSignificantDigits, fSignificantDigits);
 461
 462   TIter iter(arr);
 463   AliParticleYield * part = 0;
 464   TString normError ;
 465   while ((part = (AliParticleYield*) iter.Next())){
 466     if(part->GetNormErrorNeg()) {
 467       normError = FormatCol(Form(formatA, // Asymmetric error format
 468                                  RoundToSignificantFigures(part->GetNormErrorPos(),fSignificantDigits),
 469                                  RoundToSignificantFigures(part->GetNormErrorNeg(),fSignificantDigits)),
 470                             colWidth,
 471                             separator);
 472     }
 473     else {
 474       normError = FormatCol(Form(format, RoundToSignificantFigures(part->GetNormError(),fSignificantDigits)) , colWidth , separator);
 475     }
 476     fileOut
 477       << FormatCol(Form("%d",part->GetPdgCode())                                                    , colWidth , separator)
 478       << FormatCol(part->GetPartName()                                                              , colWidth , separator)
 479       << FormatCol(Form("%d", part->GetCollisionSystem())                                           , colWidth , separator)
 480       << FormatCol(Form(format, part->GetSqrtS())                                                   , colWidth , separator)
 481       << FormatCol(Form(format, RoundToSignificantFigures(part->GetYield(),    fSignificantDigits)) , colWidth , separator)
 482       << FormatCol(Form(format, RoundToSignificantFigures(part->GetStatError(),fSignificantDigits)) , colWidth , separator)
 483       << FormatCol(Form(format, RoundToSignificantFigures(part->GetSystError(),fSignificantDigits)) , colWidth , separator)
 484       << normError.Data()
 485       << FormatCol(Form(format, part->GetYMin())                                                    , colWidth , separator)
 486       << FormatCol(Form(format, part->GetYMax())                                                    , colWidth , separator)
 487       << FormatCol(Form("%d",part->GetStatus()          )                                           , colWidth , separator)
 488       << FormatCol(Form("%d",part->GetMeasurementType() )                                           , colWidth , separator)
 489       << FormatCol(part->GetCentr()                                                                 , colWidth , separator)
 490       << FormatCol(Form("%d",part->GetIsSum())                                                      , colWidth , separator)
 491       << FormatCol(part->GetTag()                                                                   , colWidth , separator)
 492       << endl;
 493   }
 494
 495
 496 }
 497
 498 void AliParticleYield::WriteThermusFile(TClonesArray * arr, const char * filename, Int_t colwidth) {
 499   // Writes a txt file which can we used as input in therums fits
 500
 501   if(!arr) {
 502     Printf("<AliParticleYield::WriteThermusFile> Error: no array provided");
 503     return;
 504   }
 505   if(!filename) {
 506     Printf("<AliParticleYield::WriteThermusFile> Error: no filename provided");
 507     return;
 508   }
 509
 510   ofstream fileOut(filename);
 511
 512   TIter iter(arr);
 513   AliParticleYield * part = 0;
 514   char format[20];
 515   // This is used for float numbers in the table.
 516   // The "g" options switches between the normal or scientific notation, whathever is more appropriate.
 517   // We want to have up to fSignificantDigits digits after the .
 518   snprintf(format,20,"%%%dg", fSignificantDigits);
 519
 520   //  snprintf(format, 20, "%d.%d%%f", fSignificantDigits, fSignificantDigits);
 521   while ((part = (AliParticleYield*) iter.Next())){
 522
 523     if(part->IsTypeRatio()) {
 524       // If it's a ratio we have to write the 2 pdg codes
 525       fileOut << FormatCol(Form("%d\t%d\t",part->GetPdgCode(), part->GetPdgCode2())                              , colwidth)
 526               << part->GetTag() << "\t"
 527               << Form(format, RoundToSignificantFigures(part->GetYield()      , fSignificantDigits)) << "\t"
 528               << Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits))
 529               << endl;
 530     }
 531     else {
 532       fileOut <<Form("%d",part->GetPdgCode())                                                       << "\t"
 533               <<part->GetTag()                                                                      << "\t"
 534               <<Form(format, RoundToSignificantFigures(part->GetYield()      , fSignificantDigits)) << "\t"
 535               <<Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits))
 536               << endl;
 537     }
 538
 539   }
 540
 541 }
 542
 543
 544 const char * AliParticleYield::FormatCol(const char * text, Int_t width,  const char * sep) {
 545
 546   TString format(Form("%%-%ds %s", width, sep));
 547   return Form(format.Data(), text);
 548
 549 }
 550
 551 Double_t AliParticleYield::RoundToSignificantFigures(double num, int n) {
 552   // Rounds num to n significant digits.
 553   // Recipe from :http://stackoverflow.com/questions/202302/rounding-to-an-arbitrary-number-of-significant-digits
 554   // Basically the log is used to determine the number of leading 0s, than convert to an integer by multipliing by the expo,
 555   // round the integer and shift back.
 556   if(num == 0) {
 557     return 0;
 558   }
 559
 560   Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num));
 561   Int_t power = n - (int) d;
 562
 563   Double_t magnitude = TMath::Power(10, power);
 564   Long_t shifted = TMath::Nint(num*magnitude);
 565   return shifted/magnitude;
 566
 567 }
 568
 569
 570 Float_t AliParticleYield::GetError(TString error, Float_t yield) {
 571   // The "GetError" function can handle % errors.
 572   if(error.Contains("%")) {
 573     return yield * error.Atof()/100;
 574   }
 575   return error.Atof();
 576 }
 577
 578 void AliParticleYield::SetPdgCode(TString partName) {
 579   // Set pdg code from part name, if there was a previous name, check if it is consistent
 580   TParticlePDG * part  = TDatabasePDG::Instance()->GetParticle(partName);
 581   if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums
 582   if(!part) {
 583     AliError(Form("No particle %s in TDatabasePDG", partName.Data()));
 584     return;
 585   }
 586   if(fPdgCode != part->PdgCode() &&  !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on PDG codes is disabled case of ratios
 587     AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), partName.Data(), part->PdgCode(), fPdgCode));
 588   }
 589   fPdgCode = part->PdgCode();
 590
 591 }
 592
 593 void AliParticleYield::SetPartName(Int_t pdgCode) {
 594   // Set part name from pdg code, if there was a previous code, check if it is consistent
 595   TParticlePDG * part  = TDatabasePDG::Instance()->GetParticle(pdgCode);
 596   if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums
 597   if(!part) {
 598     AliError(Form("No particle with code %d in TDatabasePDG", pdgCode));
 599     return;
 600   }
 601   if(fPdgCode != part->PdgCode() && !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on particle names is disabled case of ratios
 602     AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), part->GetName(), part->PdgCode(), fPdgCode));
 603   }
 604   fPartName = part->GetName();
 605
 606 }
 607
 608 Bool_t AliParticleYield::CheckTypeConsistency() const {
 609   // Check for any inconsistency with the type mask. Returns true if the object is fully consistent.
 610   Bool_t isOK = kTRUE;
 611
 612   if((fMeasurementType & kTypeOnlyTotError) && GetStatError()) {
 613     AliError(Form("Error flagged as total only, but stat error is not 0 (%f)!",GetStatError()));
 614     isOK= kFALSE;
 615   } else if (!(fMeasurementType & kTypeOnlyTotError) && (!GetStatError() || !GetSystError())) {
 616     AliError("Stat or syst errors are null");
 617   }
 618   if((fMeasurementType & kTypeLinearInterpolation) && (fMeasurementType & kTypeAverageAndRefit) && (fMeasurementType & kTypeExtrPionRatio)) {
 619     AliError("Only one out of the \"Liner interpolation\",  \"Average and refit\", \"Extrapolated with constant ratio to pions\" bits can be set");
 620   }
 621
 622   return isOK;
 623 }
 624
 625 void AliParticleYield::Print (Option_t *opt) const {
 626   // Print
 627   // Available options:
 628   //  - short
 629   //  - justvalue (does not print normalization error)
 630   TString sopt(opt);
 631   if(sopt.Contains("short")) {
 632     printf("[%s]: %f +- %f +- %f ", fPartName.Data(), fYield, fStatError, fSystError);
 633     if(fNormErrorNeg) {
 634       printf("(+%f-%f)", fNormErrorPos, fNormErrorNeg);
 635     }else if(fNormErrorPos) {
 636       printf("(+-%f)", fNormErrorPos);
 637     }
 638     printf("[0x%8.8x,%d]\n", fMeasurementType, fStatus);
 639   }
 640   else if (sopt.Contains("justvalue")) {
 641     Printf("%f +- %f +- %f ", fYield, fStatError, fSystError);
 642
 643   } else {
 644     Printf("-------------------------------");
 645     CheckTypeConsistency();
 646     if(fMeasurementType & kTypeParticleRatio) {
 647       Printf("%s [%s] (%d/%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fPdgCode2, fIsSum ? "particle + antiparticle" : "", fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
 648     }
 649     else{
 650       Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fIsSum ? "particle + antiparticle" : "", fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
 651     }
 652     TString measurementType = IsTypeMeasured() ? "Measured" : "";
 653     if(fMeasurementType & kTypeLinearInterpolation) measurementType += "Interpolated";
 654     if(fMeasurementType & kTypeAverageAndRefit) measurementType     += "Averaged+Refitted";
 655     if(fMeasurementType & kTypeExtrPionRatio)   measurementType     += "Extrapolated assuming constant ratio to pions";
 656     Printf("Status: %s, %s", kStatusString[fStatus],  measurementType.Data());
 657     Printf("%s , sqrt(s) = %2.2f GeV, %2.2f < y < %2.2f %s", kSystemString[fCollisionSystem], fSqrtS, fYMin, fYMax, fCentr.Data());
 658     if(fMeasurementType & kTypeOnlyTotError) {
 659       Printf("%f +- %f (total error)", fYield, fSystError);
 660     }
 661     else {
 662       Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError);
 663     }
 664     if(fNormErrorNeg) {
 665       Printf("Normalization uncertainty: +%f-%f", fNormErrorPos, fNormErrorNeg);
 666     }
 667     else {
 668       Printf("Normalization uncertainty: %f", fNormErrorPos);
 669     }
 670   }
 671 }
 672
 673 Bool_t AliParticleYield::operator==(const AliParticleYield& rhs) {
 674   // Check if the two particles are identical
 675
 676   Bool_t isEqual =
 677     (fPdgCode         == rhs.fPdgCode              ) &&
 678     (fPdgCode2        == rhs.fPdgCode2             ) &&
 679     (fPartName        == rhs.fPartName             ) &&
 680     (fCollisionSystem == rhs.fCollisionSystem      ) &&
 681     Compare2Floats(fSqrtS,rhs.fSqrtS               ) &&
 682     Compare2Floats(fYield,rhs.fYield               ) &&
 683     Compare2Floats(fStatError,rhs.fStatError       ) &&
 684     Compare2Floats(fSystError,rhs.fSystError       ) &&
 685     Compare2Floats(fNormErrorPos,rhs.fNormErrorPos ) &&
 686     Compare2Floats(fNormErrorNeg,rhs.fNormErrorNeg ) &&
 687     Compare2Floats(fYMin,rhs.fYMin                 ) &&
 688     Compare2Floats(fYMax,rhs.fYMax                 ) &&
 689     (fStatus          == rhs.fStatus               ) &&
 690     (fMeasurementType == rhs.fMeasurementType      ) &&
 691     (fCentr           == rhs.fCentr                ) &&
 692     (fIsSum           == rhs.fIsSum                ) &&
 693     (fTag             == rhs.fTag                  ) ;
 694
 695   return isEqual;
 696
 697 }
 698 Bool_t AliParticleYield::IsTheSameMeasurement(AliParticleYield &rhs) {
 699
 700   // Check the two particles represent the same measurement (independently of the values)
 701   Bool_t isEqual =
 702     (fPdgCode         == rhs.fPdgCode         ) &&
 703     (fPdgCode2        == rhs.fPdgCode2        ) &&
 704     (fCollisionSystem == rhs.fCollisionSystem ) &&
 705     Compare2Floats(fSqrtS,rhs.fSqrtS          ) &&
 706     Compare2Floats(fYMin,rhs.fYMin            ) &&
 707     Compare2Floats(fYMax,rhs.fYMax            ) &&
 708     (fStatus          == rhs.fStatus          ) &&
 709     (fCentr           == rhs.fCentr           ) &&
 710     (fIsSum           == rhs.fIsSum           ) &&
 711     (fTag             == rhs.fTag             ) ;
 712
 713   return isEqual;
 714
 715
 716 }
 717
 718 Bool_t AliParticleYield::CheckForDuplicates(TClonesArray * arr) {
 719
 720   // loop over all elements on the array and check for duplicates
 721   TIter iter(arr);
 722   AliParticleYield * part = 0;
 723   Bool_t isDuplicate = kFALSE;
 724
 725   while ((part = (AliParticleYield*) iter.Next())) {
 726     if (IsTheSameMeasurement(*part)){
 727       AliWarning("Duplicated measurement found");
 728       isDuplicate = kTRUE;
 729       if (!((*this) == (*part))) {
 730         part->Print();
 731         Print();
 732         AliFatal("The 2 particles are different!");
 733       }
 734     }
 735   }
 736   return isDuplicate;
 737
 738 }
 739
 740 Bool_t AliParticleYield::Compare2Floats(Float_t a, Float_t b) {
 741   // just a simple helper for the comparison methods
 742   if(!a) {
 743     if(!b) return kTRUE; // They are both 0;
 744     return kFALSE;// return here to avoid division by 0
 745   }
 746   Bool_t areEqual = (TMath::Abs((a - b)/a) < fEpsilon); // If the relative difference is < epsilon, returns true
 747   if(!areEqual) {
 748     Printf("Warning: %f and %f are different", a,b);
 749   }
 750   return areEqual;
 751 }
 752
 753
 754 Float_t AliParticleYield::GetNormError() const {
 755   // Returs a symmetrized error in case the normalizatione Error is asymmetric
 756   if(fNormErrorNeg) {
 757     AliWarning("Error is asymmetric, returining symmetrized uncertainty");
 758     return (TMath::Abs(fNormErrorNeg)+TMath::Abs(fNormErrorPos))/2;
 759   }
 760   else return fNormErrorPos; // If the uncertainty is not asymmetric, fNormErrorPos stores it.
 761
 762 }
 763
 764 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){
 765   // Finds a particle in array matching the search criteria. If more than one is found, prints a warning
 766   // If status is -1, tries to return the best (lower status value)
 767   // If pdg2 is not zero, we try to match it as well (we are looking for a ratio)
 768   // The centrality is compared with TString::Contains
 769
 770   TIter iter(arr);
 771   AliParticleYield * part = 0;
 772   AliParticleYield * foundPart = 0;
 773   while ((part = dynamic_cast<AliParticleYield*>(iter.Next()))){
 774     if (part->GetPdgCode() == pdg &&                     // same pdg
 775         part->GetCollisionSystem() == system &&          // same system
 776         Compare2Floats(part->GetSqrtS(), sqrts) &&       // same energy
 777         part->GetCentr().Contains(centrality) &&         // compatible centrality
 778         (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
 779         (status < 0 || part->GetStatus() == status)  &&  // same status, if requested
 780         (isSum  < 0 || part->GetIsSum()  == isSum)       // part+antipart or not, if requested
 781         ) {
 782       if (foundPart) { // we already found a patching particle
 783         Printf("WARNING<AliParticleYield::FindParticle>: Found another particle matching the same criteria");
 784         foundPart->Print();
 785         part->Print();
 786         if (part->GetStatus() == foundPart->GetStatus()) { // Does it have the same status?! Don't know what to do!
 787           Printf("WARNING<AliParticleYield::FindParticle>: they have the same status, I cannot decide, resetting particle");
 788           foundPart = 0;
 789         }
 790         else if (part->GetStatus()< foundPart->GetStatus()) { // Is it of better quality? select it!
 791           Printf("WARNING<AliParticleYield::FindParticle>: the new one has a smaller status: selecting it!");
 792           foundPart = part;
 793         }
 794       } else { // First match
 795         foundPart = part;
 796       }
 797
 798     }
 799
 800   }
 801   if(!foundPart) {
 802     Printf("ERROR<AliParticleYield::FindParticle>: Cannot find %d (System %d, sqrts = %2.2f TeV, %s, %s, Status:%d, pdg2:%d)",
 803            pdg, system, sqrts, centrality.Data(), isSum ? "part+antipart" : "", status, pdg2);
 804   }
 805   return foundPart;
 806
 807 }
 808
 809 void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield*part2, const char * pdgSep) {
 810   // Combines metadata from part1 and part2
 811   // pdgSep is a separator to be added in the name and pdg (e.g. + for a sum, / for a ratio)
 812
 813   Int_t pdg1 = part1->GetPdgCode();
 814   Int_t pdg2 = pdg1 == part2->GetPdgCode() ? part1->GetPdgCode2() : part2->GetPdgCode();
 815   Int_t   system = part1->GetCollisionSystem() == part2->GetCollisionSystem() ? part2->GetCollisionSystem() : -1;
 816   Float_t sqrts = Compare2Floats(part1->GetSqrtS(), part2->GetSqrtS()) ? part1->GetSqrtS() : 0;
 817   Int_t ymin = part1->GetYMin() == part2->GetYMin() ? part2->GetYMin() : -1000;
 818   Int_t ymax = part1->GetYMax() == part2->GetYMax() ? part2->GetYMax() : -1000;
 819   Int_t status = part1->GetStatus() == part2->GetStatus() ? part2->GetStatus() : -1;
 820   Int_t type = part1->GetMeasurementType() == part2->GetMeasurementType() ? part2->GetMeasurementType() : -1;
 821
 822   TString centr = part1->GetCentr() == part2->GetCentr() ? part2->GetCentr() : part1->GetCentr()+"/"+part2->GetCentr();
 823   TString tag = part1->GetTag() == part2->GetTag() ? part2->GetTag() : part1->GetTag()+"/"+part2->GetTag();
 824   Int_t issum = part1->GetIsSum() == part2->GetIsSum() ? part2->GetIsSum() : -1000;
 825
 826   SetPdgCode(pdg1);
 827   SetPdgCode2(pdg2);
 828   SetCollisionSystem(AliPYCSystem_t(system));
 829   SetSqrtS(sqrts);
 830   SetYMin(ymin);
 831   SetYMax(ymax);
 832   SetStatus(AliPYStatusCode_t(status));
 833   SetMeasurementType(type);
 834   SetCentr(centr);
 835   SetTag(tag);
 836   SetIsSum(issum);
 837
 838   fPartName = TDatabasePDG::Instance()->GetParticle(fPdgCode)->GetName();
 839   if(pdg2) fPartName = fPartName + pdgSep + TDatabasePDG::Instance()->GetParticle(fPdgCode2)->GetName();
 840 }
 841
 842 AliParticleYield * AliParticleYield::Add   (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt){
 843
 844   // Computes the ratio of 2 particles.
 845   // Valid options:
 846   //  - NQ: Propagates normalization errors quadratically (by default they are propagated linearly)
 847   //  - SL: propagates STATISTICAL errors linearly
 848   //  - YQ: propagates SYSTEMATIC errors quadratically
 849   //  NB by default, statistical errors are propagated quadratically and systematic errors linearly
 850   // if "Correlated error" is non null, it is subtracted in quadrature from the result. It should be a fractional error.
 851
 852   if(!part1 || !part2) {
 853     Printf("WARNING<AliParticleYield::Add>: part1 or part2 is null!");
 854     return 0;
 855   }
 856
 857   TString sopt(opt);
 858   sopt.ToUpper();
 859
 860   Float_t value = part1->GetYield() + part2->GetYield();
 861   Float_t stat  = SumErrors(part1, part2, 0, sopt.Contains("SL") ? "L": "" ); // the option decices if it is propagated linearly pr or quadratically
 862   Float_t syst  = SumErrors(part1, part2, 1, sopt.Contains("YQ") ? "" : "L" );// the option decices if it is propagated linearly pr or quadratically
 863   Float_t norm = SumErrors(part1, part2, 2, sopt.Contains("NQ") ? "" :"L");
 864
 865
 866   if(correlatedError) {
 867     syst = TMath::Sqrt(syst*syst -correlatedError*correlatedError*value*value); // FIXME: this line was never tested
 868   }
 869
 870   AliParticleYield * part = new AliParticleYield();
 871   part->SetYield(value);
 872   part->SetStatError(stat);
 873   part->SetSystError(syst);
 874   part->SetNormError(norm);
 875   part->CombineMetadata(part1, part2, "+");
 876
 877   return part;
 878
 879
 880 }
 881
 882 void AliParticleYield::Scale(Float_t scale) {
 883   // scales the measurement by an errorless number
 884   fYield        *= scale;
 885   fNormErrorNeg *= scale;
 886   fNormErrorPos *= scale;
 887   fStatError    *= scale;
 888   fSystError    *= scale;
 889
 890 }
 891
 892 AliParticleYield * AliParticleYield::Divide (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt) {
 893   // Computes the ratio of 2 particles.
 894   // Valid options:
 895   //  - NQ: assumes normalization errors to be uncorrelated and propagates them quadratically (otherwise the normalization error on the ratio is set to 0
 896   //  - SL: propagates STATISTICAL errors linearly
 897   //  - YQ: propagates SYSTEMATIC errors quadratically
 898   //  NB by default, statistical errors are propagated quadratically and systematic errors linearly
 899   // if "Correlated error" is non null, it is subtracted in quadrature from the result.It should be a fractional error.
 900
 901   if(!part1 || !part2) {
 902     Printf("WARNING<AliParticleYield::Divide>: part1 or part2 is null!");
 903     return 0;
 904   }
 905
 906   TString sopt(opt);
 907   sopt.ToUpper();
 908   if(part1->IsTypeRatio() || part2->IsTypeRatio()){
 909     Printf("WARNING<AliParticleYield::Divide>: If computing a double ratio, some meta info may be not reliable!");
 910   }
 911
 912   Float_t value = part1->GetYield() / part2->GetYield();
 913   // Since in a ratio we propagate a relative error, we have to multiply it back for value in order to get the absolute uncertainty
 914   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
 915   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
 916   Float_t norm = 0;
 917   if(sopt.Contains("NQ")) {// if opt contains N, propagate the normalization error assuming it is independent
 918     norm  = SumErrors(part1, part2, 2, "R")*value;
 919   }
 920
 921   if(correlatedError) {
 922     syst = TMath::Sqrt(syst/value*syst/value -correlatedError*correlatedError)*value; // FIXME: this line was never tested
 923   }
 924
 925   AliParticleYield * part = new AliParticleYield();
 926   part->SetYield(value);
 927   part->SetStatError(stat);
 928   part->SetSystError(syst);
 929   part->SetNormError(norm);
 930   part->CombineMetadata(part1, part2, "/");
 931
 932   part->SetMeasurementType(part->GetMeasurementType() | kTypeParticleRatio);// Set ratio bit
 933
 934   return part;
 935
 936 }
 937
 938 Double_t AliParticleYield::SumErrors(AliParticleYield * part1, AliParticleYield * part2, Int_t error, Option_t * opt) {
 939
 940   // Combines 2 errors.
 941   //  error = 0 -> statistical error
 942   //  error = 1 -> systematic error
 943   //  error = 2 -> normalization error
 944   //  Valid options
 945   //   "R" it propagates it as a relative error, WARNING: it also returns a relative error!
 946   //   "L" it propagates it sums the errors linearly (by default it is done in quadrature)
 947
 948
 949   TString sopt(opt);
 950   sopt.ToUpper();
 951   Bool_t isRelative = sopt.Contains("R");
 952   Bool_t isLinear   = sopt.Contains("L");
 953
 954   Double_t err1 = 0;
 955   Double_t err2 = 0;
 956   Double_t err  = 0;
 957   if (error == 0) {
 958     err1 = part1->GetStatError();
 959     err2 = part2->GetStatError();
 960   } else if (error == 1) {
 961     err1 = part1->GetSystError();
 962     err2 = part2->GetSystError();
 963   } else if (error == 2) {
 964     err1 = part1->GetNormError();
 965     err2 = part2->GetNormError();
 966   } else {
 967     Printf("ERROR<AliParticleYield::SumErrors>: wrong error #:%d", error);
 968   }
 969
 970   if(isRelative) {
 971     err1  /= part1->GetYield();
 972     err2  /= part2->GetYield();
 973   }
 974
 975   if(isLinear) {
 976     err = err1+err2;
 977   } else {
 978     err = TMath::Sqrt(err1*err1 + err2*err2);
 979   }
 980
 981   if(isRelative) return err;
 982
 983   return err;
 984
 985 }