5 #include "AliParticleYield.h"
6 #include "TDatabasePDG.h"
8 #include "TClonesArray.h"
13 #include "TDirectory.h"
14 #include "TEventList.h"
23 ClassImp(AliParticleYield)
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;
31 AliParticleYield::AliParticleYield() :
52 AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
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):
60 fCollisionSystem(system),
70 fMeasurementType(type),
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();
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):
88 fCollisionSystem(system),
93 fNormErrorPos(normPos),
94 fNormErrorNeg(normNeg),
98 fMeasurementType(type),
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();
112 AliParticleYield::AliParticleYield(const AliParticleYield& part) :
114 fPdgCode(part.fPdgCode),
115 fPdgCode2(part.fPdgCode2),
116 fPartName(part.fPartName),
117 fCollisionSystem(part.fCollisionSystem),
120 fStatError(part.fStatError),
121 fSystError(part.fSystError),
122 fNormErrorPos(part.fNormErrorPos),
123 fNormErrorNeg(part.fNormErrorNeg),
126 fStatus(part.fStatus),
127 fMeasurementType(part.fMeasurementType),
135 AliParticleYield::~AliParticleYield() {
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())){
148 if(part) delete part;
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);
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".
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
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).
189 // The columns should be:
190 // PDG NAME SYSTEM SQRTS VALUE SYST STAT NORM YMIN YMAX STATUS TYPE CENTR ISSUM TAG
192 // PDG should either be an integher or the ratio of two integers (in case of particle ratios), with the following format:
194 // NO SPACES ARE ALLOWED IN NAMES AND PDG CODE, unless you use a separator which is not a whitespace
196 // A Header can be present (lines beginning with the word "PDG" are also skipped
198 const Int_t kNCols = 14; // The lines are actually 15, but the last one (TAG) can be empty, so we put 14 here.
200 TClonesArray * arr = new TClonesArray ("AliParticleYield");
201 ifstream filestream (fileName);
202 if(!filestream.is_open()) {
203 Printf("Cannot open file %s\n", fileName);
208 std::cout << "Reading " << fileName << std::endl;
210 while (line.ReadLine(filestream) ) {
211 // Strip trailing and leading whitespaces
212 line = line.Strip(TString::kLeading, ' ');
213 line = line.Strip(TString::kTrailing, ' ');
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;
224 if (line.BeginsWith("PDG")) continue;
226 // Tokenize line using custom separator
227 TObjArray * cols = line.Tokenize(separators);
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);
237 // get type first, as some operations are type-specific
238 UInt_t type = ((TObjString*)cols->At(11)) ->String().Atoi();
240 // if it's a ratio, try to get the 2 pdg codes
241 Int_t pdg =0, pdg2 = 0;
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");
249 pdg = ((TObjString*)tokens->At(0)) ->String().Atoi();
250 pdg2 = ((TObjString*)tokens->At(1)) ->String().Atoi();
254 pdg = ((TObjString*)cols->At(0)) ->String().Atoi();
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());
267 if (normString.Contains("+") && normString.Contains("-")) {
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
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
281 // symmetric error: set only normpos
282 normPos = GetError(((TObjString*)cols->At(7)) ->String(), yield);
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
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, ' ');
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);
310 std::cout << "<- File read" << std::endl;
316 const char * AliParticleYield::GetLatexName(Int_t pdg) const {
318 // Returns a TLatex compatible name for the particle
319 // if pdg == 0 uses fPdgcode;
320 // We need the pdg argument for particle ratios
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.
326 name += GetLatexName(fPdgCode);
328 name += GetLatexName(fPdgCode2);
333 if(!pdg) pdg = fPdgCode;
338 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#pi^{+} + #pi^{-})/2";
339 return "(#pi^{+} + #pi^{-})";
348 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(K^{+} + K^{-})/2";
349 return "(K^{+} + K^{-})";
358 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(p + #bar{p})/2";
359 return "(p + #bar{p})";
368 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Lambda + #bar{#Lambda})/2";
369 return "(#Lambda + #bar{#Lambda})";
374 return "#bar{#Lamnba}";
378 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Xi^{-} + #bar{#Xi}^{+})/2";
379 return "(#Xi^{-} + #bar{#Xi}^{+})";
384 return "#bar{#Xi}^{+}";
388 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Omega^{-} + #bar{#Omega}^{+})/2";
389 return "(#Omega^{-} + #bar{#Omega}^{+})";
394 return "#bar{#Omega}^{+}";
404 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(K* + #bar{K*})/2";
405 return "(K* + #bar{K*})";
414 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(d + #bar{d})/2";
415 return "(d + #bar{d})";
417 return "d";// Deuteron
420 return "#bar{d}";// Deuteron
424 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(^{3}He + #bar{^{3}He})/2";
425 return "(^{3}He + #bar{^{3}He})";
430 return "#bar{^{3}He}";
434 if (fMeasurementType & kTypeAveragePartAntiPart ) return "({}^{3}_{#Lambda}H + {}^{3}_{#Lambda}#bar{H})/2";
435 return "{}^{3}_{#Lambda}H + {}^{3}_{#Lambda}#bar{H}";
437 return "{}^{3}_{#Lambda}H";
440 return "{}^{3}_{#Lambda}#bar{H}";
443 AliWarning("Latex Name not know for this particle");
446 return fPartName.Data();
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());
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
468 // The columns should be:
469 // PDG NAME SYSTEM SQRTS VALUE STAT SYST NORM YMIN YMAX STATUS TYPE CENTR ISSUM TAG
471 Printf("<AliParticleYield::SaveAsASCIIFile> Error: no array provided");
475 Printf("<AliParticleYield::SaveAsASCIIFile> Error: no filename provided");
479 ofstream fileOut(fileName);
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;
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 .
488 snprintf(format,20,"%%%dg", fSignificantDigits);
490 char formatA[30];// We have to rebuild the format for asymmetric uncertainties...
491 snprintf(formatA,30,"+%%%dg-%%%dg", fSignificantDigits, fSignificantDigits);
494 AliParticleYield * part = 0;
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)),
505 normError = FormatCol(Form(format, RoundToSignificantFigures(part->GetNormError(),fSignificantDigits)) , colWidth , separator);
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)
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)
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
533 Printf("<AliParticleYield::WriteThermusFile> Error: no array provided");
537 Printf("<AliParticleYield::WriteThermusFile> Error: no filename provided");
541 ofstream fileOut(filename);
544 AliParticleYield * part = 0;
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);
551 // snprintf(format, 20, "%d.%d%%f", fSignificantDigits, fSignificantDigits);
552 while ((part = (AliParticleYield*) iter.Next())){
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))
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))
575 const char * AliParticleYield::FormatCol(const char * text, Int_t width, const char * sep) {
577 TString format(Form("%%-%ds %s", width, sep));
578 return Form(format.Data(), text);
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.
591 Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num));
592 Int_t power = n - (int) d;
594 Double_t magnitude = TMath::Power(10, power);
595 Long_t shifted = TMath::Nint(num*magnitude);
596 return shifted/magnitude;
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;
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
614 AliError(Form("No particle %s in TDatabasePDG", partName.Data()));
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));
620 fPdgCode = part->PdgCode();
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
629 AliError(Form("No particle with code %d in TDatabasePDG", pdgCode));
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));
635 fPartName = part->GetName();
639 Bool_t AliParticleYield::CheckTypeConsistency() const {
640 // Check for any inconsistency with the type mask. Returns true if the object is fully consistent.
643 if((fMeasurementType & kTypeOnlyTotError) && GetStatError()) {
644 AliError(Form("Error flagged as total only, but stat error is not 0 (%f)!",GetStatError()));
646 } else if (!(fMeasurementType & kTypeOnlyTotError) && (!GetStatError() || !GetSystError())) {
647 AliError("Stat or syst errors are null");
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");
654 if((fMeasurementType & kTypeAveragePartAntiPart) && !fIsSum) {
655 AliError("Average part antipart set, but fIsSum is 0! This type bit should only be set for sums.");
662 void AliParticleYield::Print (Option_t *opt) const {
664 // Available options:
666 // - justvalue (does not print normalization error)
668 if(sopt.Contains("short")) {
669 printf("[%s]: %f +- %f +- %f ", fPartName.Data(), fYield, fStatError, fSystError);
671 printf("(+%f-%f)", fNormErrorPos, fNormErrorNeg);
672 }else if(fNormErrorPos) {
673 printf("(+-%f)", fNormErrorPos);
675 printf("[0x%8.8x,%d]\n", fMeasurementType, fStatus);
677 else if (sopt.Contains("justvalue")) {
678 Printf("%f +- %f +- %f ", fYield, fStatError, fSystError);
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()) : "" );
690 Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, sumType.Data(), fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
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);
702 Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError);
705 Printf("Normalization uncertainty: +%f-%f", fNormErrorPos, fNormErrorNeg);
708 Printf("Normalization uncertainty: %f", fNormErrorPos);
713 Bool_t AliParticleYield::operator==(const AliParticleYield& rhs) {
714 // Check if the two particles are identical
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 ) ;
738 Bool_t AliParticleYield::IsTheSameMeasurement(AliParticleYield &rhs) {
740 // Check the two particles represent the same measurement (independently of the values)
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 ) ;
758 Bool_t AliParticleYield::CheckForDuplicates(TClonesArray * arr) {
760 // loop over all elements on the array and check for duplicates
762 AliParticleYield * part = 0;
763 Bool_t isDuplicate = kFALSE;
765 while ((part = (AliParticleYield*) iter.Next())) {
766 if (IsTheSameMeasurement(*part)){
767 AliWarning("Duplicated measurement found");
769 if (!((*this) == (*part))) {
772 AliFatal("The 2 particles are different!");
780 Bool_t AliParticleYield::Compare2Floats(Float_t a, Float_t b) {
781 // just a simple helper for the comparison methods
783 if(!b) return kTRUE; // They are both 0;
784 return kFALSE;// return here to avoid division by 0
786 Bool_t areEqual = (TMath::Abs((a - b)/a) < fEpsilon); // If the relative difference is < epsilon, returns true
788 Printf("Warning: %f and %f are different", a,b);
794 Float_t AliParticleYield::GetNormError() const {
795 // Returs a symmetrized error in case the normalizatione Error is asymmetric
797 AliWarning("Error is asymmetric, returining symmetrized uncertainty");
798 return (TMath::Abs(fNormErrorNeg)+TMath::Abs(fNormErrorPos))/2;
800 else return fNormErrorPos; // If the uncertainty is not asymmetric, fNormErrorPos stores it.
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
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
822 if (foundPart) { // we already found a patching particle
823 Printf("WARNING<AliParticleYield::FindParticle>: Found another particle matching the same criteria");
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");
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!");
834 } else { // First match
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);
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)
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();
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();
866 Int_t issum = part1->GetIsSum() || part2->GetIsSum() ? 1 : 0;
871 SetCollisionSystem(AliPYCSystem_t(system));
875 SetStatus(AliPYStatusCode_t(status));
876 SetMeasurementType(type);
884 AliParticleYield * AliParticleYield::Add (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt){
886 // Computes the sum of 2 particles.
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.
894 if(!part1 || !part2) {
895 Printf("WARNING<AliParticleYield::Add>: part1 or part2 is null!");
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");
908 if(correlatedError) {
909 syst = TMath::Sqrt(syst*syst -correlatedError*correlatedError*value*value); // FIXME: this line was never tested
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
924 void AliParticleYield::Scale(Float_t scale) {
925 // scales the measurement by an errorless number
927 fNormErrorNeg *= scale;
928 fNormErrorPos *= scale;
934 AliParticleYield * AliParticleYield::Divide (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt) {
935 // Computes the ratio of 2 particles.
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.
943 if(!part1 || !part2) {
944 Printf("WARNING<AliParticleYield::Divide>: part1 or part2 is null!");
951 if(part1->IsTypeRatio() || part2->IsTypeRatio()){
952 Printf("WARNING<AliParticleYield::Divide>: If computing a double ratio, some meta info may be not reliable!");
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
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;
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;
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, "/");
978 part->SetMeasurementType(part->GetMeasurementType() | kTypeParticleRatio);// Set ratio bit
984 Double_t AliParticleYield::SumErrors(AliParticleYield * part1, AliParticleYield * part2, Int_t error, Option_t * opt) {
986 // Combines 2 errors.
987 // error = 0 -> statistical error
988 // error = 1 -> systematic error
989 // error = 2 -> normalization error
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)
997 Bool_t isRelative = sopt.Contains("R");
998 Bool_t isLinear = sopt.Contains("L");
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();
1013 Printf("ERROR<AliParticleYield::SumErrors>: wrong error #:%d", error);
1017 err1 /= part1->GetYield();
1018 err2 /= part2->GetYield();
1024 err = TMath::Sqrt(err1*err1 + err2*err2);
1027 if(isRelative) return err;