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