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