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