coverity fix and some clean up
[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / AliParticleYield.cxx
CommitLineData
2f86df33 1#include <iostream>
5a633e19 2#include <stdio.h>
2f86df33 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"
8820ee3c 15#include "TCut.h"
2f86df33 16
17using std::endl;
18using std::left;
19using std::setw;
5e2539d9 20using std::ifstream;
21using std::ofstream;
2f86df33 22
23ClassImp(AliParticleYield)
24
25// set statics
26const char * AliParticleYield::kStatusString[] = {"Published", "Preliminary", "Final, but not published", "May change"} ;
27const char * AliParticleYield::kSystemString[] = {"pp", "p-Pb", "Pb-Pb"} ;
28Int_t AliParticleYield::fSignificantDigits = 3;
16163c99 29Float_t AliParticleYield::fEpsilon = 0.000000000000000001;
2f86df33 30
31AliParticleYield::AliParticleYield() :
32TObject(),
33fPdgCode(0),
34fPdgCode2(0),
35fPartName(""),
36fCollisionSystem(0),
37fSqrtS(0),
38fYield(0),
39fStatError(0),
40fSystError(0),
3e8c4dd8 41fNormErrorPos(0),
42fNormErrorNeg(0),
2f86df33 43fYMin(0),
44fYMax(0),
45fStatus(0),
46fMeasurementType(0),
47fCentr(""),
48fIsSum(0),
49fTag("")
50{
51 // constructor
52 AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
53}
54
55AliParticleYield::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):
56TObject(),
57fPdgCode(pdg),
58fPdgCode2(0),
59fPartName(""),
60fCollisionSystem(system),
61fSqrtS(sqrts),
62fYield(value),
63fStatError(stat),
64fSystError(syst),
3e8c4dd8 65fNormErrorPos(norm),
66fNormErrorNeg(0),
2f86df33 67fYMin(ymin),
68fYMax(ymax),
69fStatus(status),
70fMeasurementType(type),
71fCentr(centr),
72fIsSum(isSum),
73fTag(tag)
74
75{
76 // Constructor
2f86df33 77 AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
8820ee3c 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();
2f86df33 81}
82
3e8c4dd8 83AliParticleYield::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):
84TObject(),
85fPdgCode(pdg),
86fPdgCode2(0),
87fPartName(""),
88fCollisionSystem(system),
89fSqrtS(sqrts),
90fYield(value),
91fStatError(stat),
92fSystError(syst),
93fNormErrorPos(normPos),
94fNormErrorNeg(normNeg),
95fYMin(ymin),
96fYMax(ymax),
97fStatus(status),
98fMeasurementType(type),
99fCentr(centr),
100fIsSum(isSum),
101fTag(tag)
102
103{
104 // Constructor
3e8c4dd8 105 AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB
8820ee3c 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();
3e8c4dd8 109}
110
111
2f86df33 112AliParticleYield::AliParticleYield(const AliParticleYield& part) :
113TObject(),
114fPdgCode(part.fPdgCode),
115fPdgCode2(part.fPdgCode2),
116fPartName(part.fPartName),
117fCollisionSystem(part.fCollisionSystem),
118fSqrtS(part.fSqrtS),
119fYield(part.fYield),
120fStatError(part.fStatError),
121fSystError(part.fSystError),
3e8c4dd8 122fNormErrorPos(part.fNormErrorPos),
123fNormErrorNeg(part.fNormErrorNeg),
2f86df33 124fYMin(part.fYMin),
125fYMax(part.fYMax),
126fStatus(part.fStatus),
127fMeasurementType(part.fMeasurementType),
128fCentr(part.fCentr),
129fIsSum(part.fIsSum),
130fTag(part.fTag){
131 // Copy constructor
9e1ebca0 132
2f86df33 133}
134
135AliParticleYield::~AliParticleYield() {
136
137}
138
589e4fb0 139TTree * AliParticleYield::GetTreeFromArray(TClonesArray * arr) {
140 // Returns a tree from an array of Tparticles
2f86df33 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 }
589e4fb0 148 if(part) delete part;
149 return tree;
150
151
152
153}
154
155
156TTree * 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);
2f86df33 160 delete arr;
2f86df33 161 return tree;
162}
163
589e4fb0 164TClonesArray * AliParticleYield::GetEntriesMatchingSelection(TTree * tree, TCut selection) {
2f86df33 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++){
2f86df33 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
183TClonesArray * 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);
16163c99 202 if(!filestream.is_open()) {
203 Printf("Cannot open file %s\n", fileName);
204 exit(1);
205 }
2f86df33 206 TString line;
207 Int_t ipart = 0;
ff45f222 208 std::cout << "Reading " << fileName << std::endl;
209
589e4fb0 210 while (line.ReadLine(filestream) ) {
2f86df33 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
ff45f222 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 }
2f86df33 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);
3e8c4dd8 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 }
2f86df33 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
3e8c4dd8 300 AliParticleYield * part = new AliParticleYield(pdg,system,sqrts,yield,stat,syst,normPos, normNeg,ymin,ymax,status,type,centr,issum,tag);
16163c99 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
2f86df33 309 }
3e8c4dd8 310 std::cout << "<- File read" << std::endl;
2f86df33 311
312
313 return arr;
314}
315
316const 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;
90132fab 325 name ="#frac{";
2f86df33 326 name += GetLatexName(fPdgCode);
90132fab 327 name += "}{";
2f86df33 328 name += GetLatexName(fPdgCode2);
90132fab 329 name += "}";
2f86df33 330 return name.Data();
331 }
332
333 if(!pdg) pdg = fPdgCode;
334
335 switch (pdg) {
336 case 211:
101d9a81 337 if (fIsSum) {
338 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#pi^{+} + #pi^{-})/2";
339 return "(#pi^{+} + #pi^{-})";
340 }
2f86df33 341 return "#pi^{+}";
87381af7 342 break;
2f86df33 343 case -211:
344 return "#pi^{-}";
87381af7 345 break;
2f86df33 346 case 321:
101d9a81 347 if (fIsSum) {
348 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(K^{+} + K^{-})/2";
349 return "(K^{+} + K^{-})";
350 }
2f86df33 351 return "K^{+}";
87381af7 352 break;
2f86df33 353 case -321:
354 return "K^{-}";
87381af7 355 break;
2f86df33 356 case 2212:
101d9a81 357 if (fIsSum) {
358 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(p + #bar{p})/2";
359 return "(p + #bar{p})";
360 }
2f86df33 361 return "p";
87381af7 362 break;
2f86df33 363 case -2212:
364 return "#bar^{p}";
87381af7 365 break;
2f86df33 366 case 3122:
101d9a81 367 if (fIsSum) {
368 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Lambda + #bar{#Lambda})/2";
369 return "(#Lambda + #bar{#Lambda})";
370 }
87381af7 371 return "#Lambda";
372 break;
2f86df33 373 case -3122:
87381af7 374 return "#bar{#Lamnba}";
375 break;
2f86df33 376 case 3312:
101d9a81 377 if (fIsSum) {
378 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Xi^{-} + #bar{#Xi}^{+})/2";
379 return "(#Xi^{-} + #bar{#Xi}^{+})";
380 }
2f86df33 381 return "#Xi^{-}";
87381af7 382 break;
2f86df33 383 case -3312:
384 return "#bar{#Xi}^{+}";
87381af7 385 break;
2f86df33 386 case 3334:
101d9a81 387 if (fIsSum) {
388 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Omega^{-} + #bar{#Omega}^{+})/2";
389 return "(#Omega^{-} + #bar{#Omega}^{+})";
390 }
2f86df33 391 return "#Omega^{-}";
87381af7 392 break;
2f86df33 393 case -3334:
394 return "#bar{#Omega}^{+}";
87381af7 395 break;
2f86df33 396 case 310:
397 return "K^{0}_{S}";
87381af7 398 break;
2f86df33 399 case 333:
400 return "#phi";
87381af7 401 break;
2f86df33 402 case 313:
101d9a81 403 if(fIsSum) {
404 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(K* + #bar{K*})/2";
405 return "(K* + #bar{K*})";
406 }
87381af7 407 return "K*";
408 break;
409 case -313:
410 return "#bar{K*}";
411 break;
2f86df33 412 case 1000010020:
101d9a81 413 if(fIsSum) {
414 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(d + #bar{d})/2";
415 return "(d + #bar{d})";
416 }
2f86df33 417 return "d";// Deuteron
87381af7 418 break;
2f86df33 419 case -1000010020:
420 return "#bar{d}";// Deuteron
87381af7 421 break;
2f86df33 422 case 1000020030:
101d9a81 423 if(fIsSum) {
424 if (fMeasurementType & kTypeAveragePartAntiPart ) return "(^{3}He + #bar{^{3}He})/2";
425 return "(^{3}He + #bar{^{3}He})";
426 }
2f86df33 427 return "^{3}He";
87381af7 428 break;
2f86df33 429 case -1000020030:
430 return "#bar{^{3}He}";
87381af7 431 break;
2f86df33 432 case 1010010030:
101d9a81 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 }
f10825a3 437 return "{}^{3}_{#Lambda}H";
87381af7 438 break;
2f86df33 439 case -1010010030:
f10825a3 440 return "{}^{3}_{#Lambda}#bar{H}";
87381af7 441 break;
442 default:
443 AliWarning("Latex Name not know for this particle");
2f86df33 444 }
87381af7 445
2f86df33 446 return fPartName.Data();
447
448}
449
450Float_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
464void 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
5a633e19 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 }
2f86df33 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
5a633e19 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
3e8c4dd8 490 char formatA[30];// We have to rebuild the format for asymmetric uncertainties...
491 snprintf(formatA,30,"+%%%dg-%%%dg", fSignificantDigits, fSignificantDigits);
492
2f86df33 493 TIter iter(arr);
494 AliParticleYield * part = 0;
3e8c4dd8 495 TString normError ;
2f86df33 496 while ((part = (AliParticleYield*) iter.Next())){
3e8c4dd8 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 }
2f86df33 507 fileOut
5a633e19 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)
3e8c4dd8 515 << normError.Data()
5a633e19 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)
2f86df33 523 << endl;
524 }
525
526
527}
528
5a633e19 529void 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");
3e8c4dd8 538 return;
5a633e19 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
cffbaa61 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))
5a633e19 560 << endl;
561 }
562 else {
cffbaa61 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))
5a633e19 567 << endl;
568 }
569
570 }
571
572}
573
574
2f86df33 575const 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
582Double_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
601Float_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
609void 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
624void 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
639Bool_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");
8ac70929 648 isOK = kFALSE;
2f86df33 649 }
87381af7 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");
8ac70929 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;
87381af7 657 }
2f86df33 658
659 return isOK;
660}
661
87381af7 662void AliParticleYield::Print (Option_t *opt) const {
2f86df33 663 // Print
87381af7 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);
3e8c4dd8 676 }
87381af7 677 else if (sopt.Contains("justvalue")) {
678 Printf("%f +- %f +- %f ", fYield, fStatError, fSystError);
679
680 } else {
681 Printf("-------------------------------");
682 CheckTypeConsistency();
8ac70929 683 TString sumType = "";
684 if (fIsSum && (fMeasurementType & kTypeAveragePartAntiPart)) sumType = "(particle + antiparticle)/2";
685 else if (fIsSum && !(fMeasurementType & kTypeAveragePartAntiPart)) sumType = "particle + antiparticle";
87381af7 686 if(fMeasurementType & kTypeParticleRatio) {
8ac70929 687 Printf("%s [%s] (%d/%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fPdgCode2, sumType.Data(), fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
87381af7 688 }
689 else{
8ac70929 690 Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, sumType.Data(), fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
87381af7 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 }
3e8c4dd8 710 }
2f86df33 711}
712
16163c99 713Bool_t AliParticleYield::operator==(const AliParticleYield& rhs) {
714 // Check if the two particles are identical
715
716 Bool_t isEqual =
3e8c4dd8 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 ) ;
16163c99 734
735 return isEqual;
736
737}
738Bool_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
758Bool_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
780Bool_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}
3e8c4dd8 792
793
794Float_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}
90132fab 803
804AliParticleYield * 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
87d00916 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
90132fab 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 }
87381af7 841 if(!foundPart) {
dce8f600 842 Printf("ERROR<AliParticleYield::FindParticle>: Cannot find %d (System %d, sqrts = %2.2f GeV, %s, %s, Status:%d, pdg2:%d)",
87381af7 843 pdg, system, sqrts, centrality.Data(), isSum ? "part+antipart" : "", status, pdg2);
844 }
90132fab 845 return foundPart;
846
847}
848
87381af7 849void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield*part2, const char * pdgSep) {
90132fab 850 // Combines metadata from part1 and part2
87381af7 851 // pdgSep is a separator to be added in the name and pdg (e.g. + for a sum, / for a ratio)
852
90132fab 853 Int_t pdg1 = part1->GetPdgCode();
87d00916 854 Int_t pdg2 = pdg1 == part2->GetPdgCode() ? part1->GetPdgCode2() : part2->GetPdgCode();
90132fab 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;
9e1ebca0 860 Int_t type = part1->GetMeasurementType() | part2->GetMeasurementType();
90132fab 861
f10825a3 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;
90132fab 867
9e1ebca0 868
90132fab 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);
f10825a3 880 SetPartName(name);
90132fab 881
882}
883
884AliParticleYield * AliParticleYield::Add (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt){
885
9e1ebca0 886 // Computes the sum of 2 particles.
90132fab 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);
87381af7 917 part->CombineMetadata(part1, part2, "+");
f10825a3 918 part->SetIsSum(1); // CombineMetadata inherits this form part1 and part2
90132fab 919 return part;
920
921
922}
923
924void 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
934AliParticleYield * 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
9e1ebca0 948
90132fab 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();
87d00916 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
90132fab 959 Float_t norm = 0;
960 if(sopt.Contains("NQ")) {// if opt contains N, propagate the normalization error assuming it is independent
87d00916 961 norm = SumErrors(part1, part2, 2, "R")*value;
90132fab 962 }
963
964 if(correlatedError) {
8ac70929 965 std::cout << "Subtracting correlated error " << correlatedError << std::endl;
966 std::cout << "Before : " << syst << "[" << syst/value*100 <<"%]"<< std::endl;
90132fab 967 syst = TMath::Sqrt(syst/value*syst/value -correlatedError*correlatedError)*value; // FIXME: this line was never tested
8ac70929 968 std::cout << "After : " << syst << "[" << syst/value*100 <<"%]"<< std::endl;
90132fab 969 }
970
971 AliParticleYield * part = new AliParticleYield();
972 part->SetYield(value);
973 part->SetStatError(stat);
974 part->SetSystError(syst);
975 part->SetNormError(norm);
87381af7 976 part->CombineMetadata(part1, part2, "/");
977
978 part->SetMeasurementType(part->GetMeasurementType() | kTypeParticleRatio);// Set ratio bit
90132fab 979
980 return part;
981
982}
983
984Double_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
87d00916 991 // "R" it propagates it as a relative error, WARNING: it also returns a relative error!
90132fab 992 // "L" it propagates it sums the errors linearly (by default it is done in quadrature)
87d00916 993
90132fab 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
87d00916 1027 if(isRelative) return err;
90132fab 1028
1029 return err;
1030
1031}