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