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