]>
Commit | Line | Data |
---|---|---|
2f86df33 | 1 | #include <iostream> |
5a633e19 | 2 | #include <stdio.h> |
2f86df33 | 3 | #include <fstream> |
4 | #include <iomanip> | |
5 | #include "AliParticleYield.h" | |
6 | #include "TDatabasePDG.h" | |
7 | #include "AliLog.h" | |
8 | #include "TClonesArray.h" | |
9 | #include "TMath.h" | |
10 | #include "AliPDG.h" | |
11 | #include "TBranch.h" | |
12 | #include "TTree.h" | |
13 | #include "TDirectory.h" | |
14 | #include "TEventList.h" | |
15 | ||
16 | using std::endl; | |
17 | using std::left; | |
18 | using std::setw; | |
19 | ||
20 | ClassImp(AliParticleYield) | |
21 | ||
22 | // set statics | |
23 | const char * AliParticleYield::kStatusString[] = {"Published", "Preliminary", "Final, but not published", "May change"} ; | |
24 | const char * AliParticleYield::kSystemString[] = {"pp", "p-Pb", "Pb-Pb"} ; | |
25 | Int_t AliParticleYield::fSignificantDigits = 3; | |
16163c99 | 26 | Float_t AliParticleYield::fEpsilon = 0.000000000000000001; |
2f86df33 | 27 | |
28 | AliParticleYield::AliParticleYield() : | |
29 | TObject(), | |
30 | fPdgCode(0), | |
31 | fPdgCode2(0), | |
32 | fPartName(""), | |
33 | fCollisionSystem(0), | |
34 | fSqrtS(0), | |
35 | fYield(0), | |
36 | fStatError(0), | |
37 | fSystError(0), | |
38 | fNormError(0), | |
39 | fYMin(0), | |
40 | fYMax(0), | |
41 | fStatus(0), | |
42 | fMeasurementType(0), | |
43 | fCentr(""), | |
44 | fIsSum(0), | |
45 | fTag("") | |
46 | { | |
47 | // constructor | |
48 | AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB | |
49 | } | |
50 | ||
51 | 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): | |
52 | TObject(), | |
53 | fPdgCode(pdg), | |
54 | fPdgCode2(0), | |
55 | fPartName(""), | |
56 | fCollisionSystem(system), | |
57 | fSqrtS(sqrts), | |
58 | fYield(value), | |
59 | fStatError(stat), | |
60 | fSystError(syst), | |
61 | fNormError(norm), | |
62 | fYMin(ymin), | |
63 | fYMax(ymax), | |
64 | fStatus(status), | |
65 | fMeasurementType(type), | |
66 | fCentr(centr), | |
67 | fIsSum(isSum), | |
68 | fTag(tag) | |
69 | ||
70 | { | |
71 | // Constructor | |
72 | fPartName = TDatabasePDG::Instance()->GetParticle(fPdgCode)->GetName(); | |
73 | AliPDG::AddParticlesToPdgDataBase(); // Make sure that ALICE-defined particles were added to the PDG DB | |
74 | } | |
75 | ||
76 | AliParticleYield::AliParticleYield(const AliParticleYield& part) : | |
77 | TObject(), | |
78 | fPdgCode(part.fPdgCode), | |
79 | fPdgCode2(part.fPdgCode2), | |
80 | fPartName(part.fPartName), | |
81 | fCollisionSystem(part.fCollisionSystem), | |
82 | fSqrtS(part.fSqrtS), | |
83 | fYield(part.fYield), | |
84 | fStatError(part.fStatError), | |
85 | fSystError(part.fSystError), | |
86 | fNormError(part.fNormError), | |
87 | fYMin(part.fYMin), | |
88 | fYMax(part.fYMax), | |
89 | fStatus(part.fStatus), | |
90 | fMeasurementType(part.fMeasurementType), | |
91 | fCentr(part.fCentr), | |
92 | fIsSum(part.fIsSum), | |
93 | fTag(part.fTag){ | |
94 | // Copy constructor | |
95 | } | |
96 | ||
97 | AliParticleYield::~AliParticleYield() { | |
98 | ||
99 | } | |
100 | ||
101 | TTree * AliParticleYield::ReadFromASCIIFileAsTree(const char * fileName, const char * separators){ | |
102 | // Read the table from an ASCII File and returns a tree of particles. See ReadFromASCIIFile for detailed info on the format | |
103 | TClonesArray * arr = ReadFromASCIIFile(fileName, separators); | |
104 | AliParticleYield * part = 0; | |
105 | TTree * tree = new TTree ("treePart", "Particle Yields and Ratios"); | |
106 | tree->Branch("particles", &part); | |
107 | TIter iterPart (arr); | |
108 | while ((part = (AliParticleYield*) iterPart.Next())){ | |
109 | tree->Fill(); | |
110 | } | |
111 | ||
112 | delete arr; | |
113 | delete part; | |
114 | return tree; | |
115 | } | |
116 | ||
117 | TClonesArray * AliParticleYield::GetEntriesMatchingSelection(TTree * tree, TString selection) { | |
118 | // 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". | |
119 | ||
120 | TClonesArray * arr = new TClonesArray("AliParticleYield"); | |
121 | AliParticleYield * part = 0; | |
122 | tree->SetBranchAddress("particles", &part); | |
123 | // 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. | |
124 | tree->Draw(">>particlelist", selection);// Produce selection list | |
125 | TEventList *elist = (TEventList*)gDirectory->Get("particlelist"); | |
126 | Int_t npart = elist->GetN(); | |
127 | for(Int_t ipart = 0; ipart < npart; ipart++){ | |
16163c99 | 128 | // std::cout << ipart << " " << elist->GetEntry(ipart) << std::endl; |
2f86df33 | 129 | tree->GetEntry(elist->GetEntry(ipart)); |
130 | new((*arr)[ipart]) AliParticleYield(*part);// We need to clone part, because it is overwritten by the next read | |
131 | } | |
132 | elist->Delete(); | |
133 | return arr; | |
134 | } | |
135 | ||
136 | ||
137 | TClonesArray * AliParticleYield::ReadFromASCIIFile(const char * fileName, const char * separators){ | |
138 | // Read the table from an ASCII File with the format indicated | |
139 | // below. Returns a TClonesArray of AliParticleyields with the | |
140 | // content of the lines. Lines beginning by "#" are skipped. | |
141 | // The order of the columns is compulsory, but the separator can be set (by default whitespaces are assumed). | |
142 | ||
143 | // The columns should be: | |
144 | // PDG NAME SYSTEM SQRTS VALUE SYST STAT NORM YMIN YMAX STATUS TYPE CENTR ISSUM TAG | |
145 | ||
146 | // PDG should either be an integher or the ratio of two integers (in case of particle ratios), with the following format: | |
147 | // pdgcode/pdgcode2 | |
148 | // NO SPACES ARE ALLOWED IN NAMES AND PDG CODE, unless you use a separator which is not a whitespace | |
149 | ||
150 | // A Header can be present (lines beginning with the word "PDG" are also skipped | |
151 | ||
152 | const Int_t kNCols = 14; // The lines are actually 15, but the last one (TAG) can be empty, so we put 14 here. | |
153 | ||
154 | TClonesArray * arr = new TClonesArray ("AliParticleYield"); | |
155 | ifstream filestream (fileName); | |
16163c99 | 156 | if(!filestream.is_open()) { |
157 | Printf("Cannot open file %s\n", fileName); | |
158 | exit(1); | |
159 | } | |
2f86df33 | 160 | TString line; |
161 | Int_t ipart = 0; | |
162 | while (line.ReadLine(filestream) ) { | |
163 | // Strip trailing and leading whitespaces | |
164 | line = line.Strip(TString::kLeading, ' '); | |
165 | line = line.Strip(TString::kTrailing, ' '); | |
166 | ||
167 | // Skip commented lines and headers | |
168 | if (line.BeginsWith("#")) continue; | |
169 | if (line.BeginsWith("PDG")) continue; | |
170 | ||
171 | // Tokenize line using custom separator | |
172 | TObjArray * cols = line.Tokenize(separators); | |
173 | ||
174 | // Check the number of columns | |
175 | if(cols->GetEntries() < kNCols) { | |
176 | Printf("Wrong number of columns in table %d vs %d expected" , cols->GetEntries(), kNCols); | |
177 | delete arr; | |
178 | return NULL; | |
179 | } | |
180 | ||
181 | // Get Values | |
182 | // get type first, as some operations are type-specific | |
183 | UInt_t type = ((TObjString*)cols->At(11)) ->String().Atoi(); | |
184 | ||
185 | // if it's a ratio, try to get the 2 pdg codes | |
186 | Int_t pdg =0, pdg2 = 0; | |
187 | ||
188 | if (type & kTypeParticleRatio) { | |
189 | TString col0 = ((TObjString*)cols->At(0)) ->String(); | |
190 | TObjArray * tokens = col0.Tokenize("/"); | |
191 | if(tokens->GetEntries() != 2) { | |
192 | Printf("ERROR: Cannot get both PDGs for ratios"); | |
193 | } else { | |
194 | pdg = ((TObjString*)tokens->At(0)) ->String().Atoi(); | |
195 | pdg2 = ((TObjString*)tokens->At(1)) ->String().Atoi(); | |
196 | } | |
197 | } | |
198 | else { | |
199 | pdg = ((TObjString*)cols->At(0)) ->String().Atoi(); | |
200 | } | |
201 | TString name = ((TObjString*)cols->At(1)) ->String(); | |
202 | Int_t system = ((TObjString*)cols->At(2)) ->String().Atoi(); | |
203 | Float_t sqrts = ((TObjString*)cols->At(3)) ->String().Atof(); | |
204 | Float_t yield = ((TObjString*)cols->At(4)) ->String().Atof(); | |
205 | // The "GetError" function can handle % errors. | |
206 | Float_t stat = GetError(((TObjString*)cols->At(5)) ->String(), yield); | |
207 | Float_t syst = GetError(((TObjString*)cols->At(6)) ->String(), yield); | |
208 | Float_t norm = GetError(((TObjString*)cols->At(7)) ->String(), yield); | |
209 | Float_t ymin = ((TObjString*)cols->At(8)) ->String().Atof(); | |
210 | Float_t ymax = ((TObjString*)cols->At(9)) ->String().Atof(); | |
211 | Int_t status = ((TObjString*)cols->At(10)) ->String().Atoi(); | |
212 | TString centr = ((TObjString*)cols->At(12)) ->String(); | |
213 | Int_t issum = ((TObjString*)cols->At(13)) ->String().Atoi(); | |
214 | TString tag = cols->At(14) ? ((TObjString*)cols->At(14)) ->String() : ""; // tag can be empty | |
215 | ||
216 | // Cleanup strings | |
217 | name = name.Strip(TString::kLeading, ' '); | |
218 | name = name.Strip(TString::kTrailing, ' '); | |
219 | centr = centr.Strip(TString::kLeading, ' '); | |
220 | centr = centr.Strip(TString::kTrailing, ' '); | |
221 | tag = tag.Strip(TString::kLeading, ' '); | |
222 | tag = tag.Strip(TString::kTrailing, ' '); | |
223 | ||
224 | // add to array | |
16163c99 | 225 | AliParticleYield * part = new AliParticleYield(pdg,system,sqrts,yield,stat,syst,norm,ymin,ymax,status,type,centr,issum,tag); |
226 | part->SetPartName(name); // Check name and PDG code consistency | |
227 | part->SetPdgCode2(pdg2); // Set second PDG code in case of ratios | |
228 | part->CheckTypeConsistency(); | |
229 | if(!part->CheckForDuplicates(arr)) { | |
230 | new ((*arr)[ipart++]) AliParticleYield(*part); | |
231 | } | |
232 | // delete part; | |
233 | ||
2f86df33 | 234 | } |
235 | ||
236 | ||
237 | return arr; | |
238 | } | |
239 | ||
240 | const char * AliParticleYield::GetLatexName(Int_t pdg) const { | |
241 | ||
242 | // Returns a TLatex compatible name for the particle | |
243 | // if pdg == 0 uses fPdgcode; | |
244 | // We need the pdg argument for particle ratios | |
245 | ||
246 | if(!pdg && fMeasurementType & kTypeParticleRatio) { | |
247 | // 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. | |
248 | static TString name; | |
249 | name =""; | |
250 | name += GetLatexName(fPdgCode); | |
251 | name += " / "; | |
252 | name += GetLatexName(fPdgCode2); | |
253 | return name.Data(); | |
254 | } | |
255 | ||
256 | if(!pdg) pdg = fPdgCode; | |
257 | ||
258 | switch (pdg) { | |
259 | case 211: | |
260 | if (fIsSum) return "(#pi^{+} + #pi^{-})"; | |
261 | return "#pi^{+}"; | |
262 | case -211: | |
263 | return "#pi^{-}"; | |
264 | case 321: | |
265 | if (fIsSum) return "(K^{+} + K^{-})"; | |
266 | return "K^{+}"; | |
267 | case -321: | |
268 | return "K^{-}"; | |
269 | case 2212: | |
270 | if (fIsSum) return "(p + #bar{p})"; | |
271 | return "p"; | |
272 | case -2212: | |
273 | return "#bar^{p}"; | |
274 | case 3122: | |
275 | if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})"; | |
276 | return "#Omega^{-}"; | |
277 | case -3122: | |
278 | return "#bar{#Omega}^{+}"; | |
279 | case 3312: | |
280 | if (fIsSum) return "(#Xi^{-} + #bar{#Xi}^{+})"; | |
281 | return "#Xi^{-}"; | |
282 | case -3312: | |
283 | return "#bar{#Xi}^{+}"; | |
284 | case 3334: | |
285 | if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})"; | |
286 | return "#Omega^{-}"; | |
287 | case -3334: | |
288 | return "#bar{#Omega}^{+}"; | |
289 | case 310: | |
290 | return "K^{0}_{S}"; | |
291 | case 333: | |
292 | return "#phi"; | |
293 | case 313: | |
294 | return "K^{*}"; | |
295 | case 1000010020: | |
296 | if(fIsSum) return "(d + #bar{d})"; | |
297 | return "d";// Deuteron | |
298 | case -1000010020: | |
299 | return "#bar{d}";// Deuteron | |
300 | case 1000020030: | |
301 | if(fIsSum) return "(^{3}He + #bar{^{3}He})"; | |
302 | return "^{3}He"; | |
303 | case -1000020030: | |
304 | return "#bar{^{3}He}"; | |
305 | case 1010010030: | |
306 | if(fIsSum) return "^{3}_{#Lambda}He + #bar{^{3}_{#Lambda}He}"; | |
307 | return "^{3}_{#Lambda}He"; | |
308 | case -1010010030: | |
309 | return "#bar{^{3}_{#Lambda}He}"; | |
310 | } | |
311 | AliWarning("Latex Name not know for this particle"); | |
312 | return fPartName.Data(); | |
313 | ||
314 | } | |
315 | ||
316 | Float_t AliParticleYield::GetTotalError(Bool_t includeNormalization) const { | |
317 | // Returns the total error, including or not the normalization uncertainty | |
318 | // All uncertainties are supposed to be uncorrelated (e.g. summed in quadrature) | |
319 | // If stat and syst are stored separately, the total error is computed summing them in quadrature | |
320 | Float_t error = GetSystError(); | |
321 | if (!(fMeasurementType & kTypeOnlyTotError)) error = TMath::Sqrt(error*error + GetStatError()*GetStatError()); | |
322 | if(includeNormalization) error = TMath::Sqrt(error*error + GetNormError()*GetNormError()); | |
323 | ||
324 | return error; | |
325 | ||
326 | ||
327 | } | |
328 | ||
329 | ||
330 | void AliParticleYield::SaveAsASCIIFile(TClonesArray * arr, const char * fileName, const char * separator, Int_t colWidth){ | |
331 | // Saves the array as an ASCII File with the format indicated | |
332 | // below. | |
333 | ||
334 | // The columns should be: | |
335 | // PDG NAME SYSTEM SQRTS VALUE STAT SYST NORM YMIN YMAX STATUS TYPE CENTR ISSUM TAG | |
5a633e19 | 336 | if(!arr) { |
337 | Printf("<AliParticleYield::SaveAsASCIIFile> Error: no array provided"); | |
338 | return; | |
339 | } | |
340 | if(!fileName) { | |
341 | Printf("<AliParticleYield::SaveAsASCIIFile> Error: no filename provided"); | |
342 | } | |
2f86df33 | 343 | |
344 | ||
345 | ofstream fileOut(fileName); | |
346 | //print header | |
347 | 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; | |
348 | ||
5a633e19 | 349 | |
350 | // This is used for float numbers in the table. | |
351 | // The "g" options switches between the normal or scientific notation, whathever is more appropriate. | |
352 | // We want to have up to fSignificantDigits digits after the . | |
353 | char format[20]; | |
354 | snprintf(format,20,"%%%dg", fSignificantDigits); | |
355 | ||
2f86df33 | 356 | TIter iter(arr); |
357 | AliParticleYield * part = 0; | |
358 | while ((part = (AliParticleYield*) iter.Next())){ | |
359 | fileOut | |
5a633e19 | 360 | << FormatCol(Form("%d",part->GetPdgCode()) , colWidth , separator) |
361 | << FormatCol(part->GetPartName() , colWidth , separator) | |
362 | << FormatCol(Form("%d", part->GetCollisionSystem()) , colWidth , separator) | |
363 | << FormatCol(Form(format, part->GetSqrtS()) , colWidth , separator) | |
364 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetYield(), fSignificantDigits)) , colWidth , separator) | |
365 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetStatError(),fSignificantDigits)) , colWidth , separator) | |
366 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetSystError(),fSignificantDigits)) , colWidth , separator) | |
367 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetNormError(),fSignificantDigits)) , colWidth , separator) | |
368 | << FormatCol(Form(format, part->GetYMin()) , colWidth , separator) | |
369 | << FormatCol(Form(format, part->GetYMax()) , colWidth , separator) | |
370 | << FormatCol(Form("%d",part->GetStatus() ) , colWidth , separator) | |
371 | << FormatCol(Form("%d",part->GetMeasurementType() ) , colWidth , separator) | |
372 | << FormatCol(part->GetCentr() , colWidth , separator) | |
373 | << FormatCol(Form("%d",part->GetIsSum()) , colWidth , separator) | |
374 | << FormatCol(part->GetTag() , colWidth , separator) | |
2f86df33 | 375 | << endl; |
376 | } | |
377 | ||
378 | ||
379 | } | |
380 | ||
5a633e19 | 381 | void AliParticleYield::WriteThermusFile(TClonesArray * arr, const char * filename, Int_t colwidth) { |
382 | // Writes a txt file which can we used as input in therums fits | |
383 | ||
384 | if(!arr) { | |
385 | Printf("<AliParticleYield::WriteThermusFile> Error: no array provided"); | |
386 | return; | |
387 | } | |
388 | if(!filename) { | |
389 | Printf("<AliParticleYield::WriteThermusFile> Error: no filename provided"); | |
390 | } | |
391 | ||
392 | ofstream fileOut(filename); | |
393 | ||
394 | TIter iter(arr); | |
395 | AliParticleYield * part = 0; | |
396 | char format[20]; | |
397 | // This is used for float numbers in the table. | |
398 | // The "g" options switches between the normal or scientific notation, whathever is more appropriate. | |
399 | // We want to have up to fSignificantDigits digits after the . | |
400 | snprintf(format,20,"%%%dg", fSignificantDigits); | |
401 | ||
402 | // snprintf(format, 20, "%d.%d%%f", fSignificantDigits, fSignificantDigits); | |
403 | while ((part = (AliParticleYield*) iter.Next())){ | |
404 | ||
405 | if(part->IsTypeRatio()) { | |
406 | // If it's a ratio we have to write the 2 pdg codes | |
407 | fileOut << FormatCol(Form("%d %d ",part->GetPdgCode(), part->GetPdgCode2()) , colwidth) | |
408 | << FormatCol(part->GetTag() , colwidth) | |
409 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetYield() , fSignificantDigits)) , colwidth) | |
410 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits)) , colwidth) | |
411 | << endl; | |
412 | } | |
413 | else { | |
414 | fileOut << FormatCol(Form("%d",part->GetPdgCode()) , colwidth) | |
415 | << FormatCol(part->GetTag() , colwidth) | |
416 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetYield() , fSignificantDigits)) , colwidth) | |
417 | << FormatCol(Form(format, RoundToSignificantFigures(part->GetTotalError() , fSignificantDigits)) , colwidth) | |
418 | << endl; | |
419 | } | |
420 | ||
421 | } | |
422 | ||
423 | } | |
424 | ||
425 | ||
2f86df33 | 426 | const char * AliParticleYield::FormatCol(const char * text, Int_t width, const char * sep) { |
427 | ||
428 | TString format(Form("%%-%ds %s", width, sep)); | |
429 | return Form(format.Data(), text); | |
430 | ||
431 | } | |
432 | ||
433 | Double_t AliParticleYield::RoundToSignificantFigures(double num, int n) { | |
434 | // Rounds num to n significant digits. | |
435 | // Recipe from :http://stackoverflow.com/questions/202302/rounding-to-an-arbitrary-number-of-significant-digits | |
436 | // Basically the log is used to determine the number of leading 0s, than convert to an integer by multipliing by the expo, | |
437 | // round the integer and shift back. | |
438 | if(num == 0) { | |
439 | return 0; | |
440 | } | |
441 | ||
442 | Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num)); | |
443 | Int_t power = n - (int) d; | |
444 | ||
445 | Double_t magnitude = TMath::Power(10, power); | |
446 | Long_t shifted = TMath::Nint(num*magnitude); | |
447 | return shifted/magnitude; | |
448 | ||
449 | } | |
450 | ||
451 | ||
452 | Float_t AliParticleYield::GetError(TString error, Float_t yield) { | |
453 | // The "GetError" function can handle % errors. | |
454 | if(error.Contains("%")) { | |
455 | return yield * error.Atof()/100; | |
456 | } | |
457 | return error.Atof(); | |
458 | } | |
459 | ||
460 | void AliParticleYield::SetPdgCode(TString partName) { | |
461 | // Set pdg code from part name, if there was a previous name, check if it is consistent | |
462 | TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(partName); | |
463 | if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums | |
464 | if(!part) { | |
465 | AliError(Form("No particle %s in TDatabasePDG", partName.Data())); | |
466 | return; | |
467 | } | |
468 | if(fPdgCode != part->PdgCode() && !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on PDG codes is disabled case of ratios | |
469 | AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), partName.Data(), part->PdgCode(), fPdgCode)); | |
470 | } | |
471 | fPdgCode = part->PdgCode(); | |
472 | ||
473 | } | |
474 | ||
475 | void AliParticleYield::SetPartName(Int_t pdgCode) { | |
476 | // Set part name from pdg code, if there was a previous code, check if it is consistent | |
477 | TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(pdgCode); | |
478 | if(IsTypeRatio() || fIsSum) return; // Name check does not make sense for ratios and sums | |
479 | if(!part) { | |
480 | AliError(Form("No particle with code %d in TDatabasePDG", pdgCode)); | |
481 | return; | |
482 | } | |
483 | if(fPdgCode != part->PdgCode() && !(fMeasurementType&kTypeParticleRatio)) { // The consistency check on particle names is disabled case of ratios | |
484 | AliError(Form("Name and pdg code are not consistent! fPartName: %s partName %s, pdgcode %d fPdgCode %d", fPartName.Data(), part->GetName(), part->PdgCode(), fPdgCode)); | |
485 | } | |
486 | fPartName = part->GetName(); | |
487 | ||
488 | } | |
489 | ||
490 | Bool_t AliParticleYield::CheckTypeConsistency() const { | |
491 | // Check for any inconsistency with the type mask. Returns true if the object is fully consistent. | |
492 | Bool_t isOK = kTRUE; | |
493 | ||
494 | if((fMeasurementType & kTypeOnlyTotError) && GetStatError()) { | |
495 | AliError(Form("Error flagged as total only, but stat error is not 0 (%f)!",GetStatError())); | |
496 | isOK= kFALSE; | |
497 | } else if (!(fMeasurementType & kTypeOnlyTotError) && (!GetStatError() || !GetSystError())) { | |
498 | AliError("Stat or syst errors are null"); | |
499 | } | |
500 | ||
501 | return isOK; | |
502 | } | |
503 | ||
504 | void AliParticleYield::Print (Option_t *) const { | |
505 | ||
506 | Printf("-------------------------------"); | |
507 | CheckTypeConsistency(); | |
508 | Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fIsSum ? "particle + antiparticle" : "", fTag.Length() ? Form("[%s]", fTag.Data()) : "" ); | |
509 | Printf("Status: %s, %s", kStatusString[fStatus], fMeasurementType & kTypeLinearInterpolation ? "Interpolated" : "Measured"); | |
510 | Printf("%s , sqrt(s) = %2.2f GeV, %2.2f < y < %2.2f %s", kSystemString[fCollisionSystem], fSqrtS, fYMin, fYMax, fCentr.Data()); | |
511 | if(fMeasurementType & kTypeOnlyTotError) { | |
512 | Printf("%f +- %f (total error)", fYield, fSystError); | |
513 | } | |
514 | else { | |
515 | Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError); | |
516 | } | |
517 | Printf("Normalizaion uncertainty: %f", fNormError); | |
518 | } | |
519 | ||
16163c99 | 520 | Bool_t AliParticleYield::operator==(const AliParticleYield& rhs) { |
521 | // Check if the two particles are identical | |
522 | ||
523 | Bool_t isEqual = | |
524 | (fPdgCode == rhs.fPdgCode ) && | |
525 | (fPdgCode2 == rhs.fPdgCode2 ) && | |
526 | (fPartName == rhs.fPartName ) && | |
527 | (fCollisionSystem == rhs.fCollisionSystem ) && | |
528 | Compare2Floats(fSqrtS,rhs.fSqrtS ) && | |
529 | Compare2Floats(fYield,rhs.fYield ) && | |
530 | Compare2Floats(fStatError,rhs.fStatError ) && | |
531 | Compare2Floats(fSystError,rhs.fSystError ) && | |
532 | Compare2Floats(fNormError,rhs.fNormError ) && | |
533 | Compare2Floats(fYMin,rhs.fYMin ) && | |
534 | Compare2Floats(fYMax,rhs.fYMax ) && | |
535 | (fStatus == rhs.fStatus ) && | |
536 | (fMeasurementType == rhs.fMeasurementType ) && | |
537 | (fCentr == rhs.fCentr ) && | |
538 | (fIsSum == rhs.fIsSum ) && | |
539 | (fTag == rhs.fTag ) ; | |
540 | ||
541 | return isEqual; | |
542 | ||
543 | } | |
544 | Bool_t AliParticleYield::IsTheSameMeasurement(AliParticleYield &rhs) { | |
545 | ||
546 | // Check the two particles represent the same measurement (independently of the values) | |
547 | Bool_t isEqual = | |
548 | (fPdgCode == rhs.fPdgCode ) && | |
549 | (fPdgCode2 == rhs.fPdgCode2 ) && | |
550 | (fCollisionSystem == rhs.fCollisionSystem ) && | |
551 | Compare2Floats(fSqrtS,rhs.fSqrtS ) && | |
552 | Compare2Floats(fYMin,rhs.fYMin ) && | |
553 | Compare2Floats(fYMax,rhs.fYMax ) && | |
554 | (fStatus == rhs.fStatus ) && | |
555 | (fCentr == rhs.fCentr ) && | |
556 | (fIsSum == rhs.fIsSum ) && | |
557 | (fTag == rhs.fTag ) ; | |
558 | ||
559 | return isEqual; | |
560 | ||
561 | ||
562 | } | |
563 | ||
564 | Bool_t AliParticleYield::CheckForDuplicates(TClonesArray * arr) { | |
565 | ||
566 | // loop over all elements on the array and check for duplicates | |
567 | TIter iter(arr); | |
568 | AliParticleYield * part = 0; | |
569 | Bool_t isDuplicate = kFALSE; | |
570 | ||
571 | while ((part = (AliParticleYield*) iter.Next())) { | |
572 | if (IsTheSameMeasurement(*part)){ | |
573 | AliWarning("Duplicated measurement found"); | |
574 | isDuplicate = kTRUE; | |
575 | if (!((*this) == (*part))) { | |
576 | part->Print(); | |
577 | Print(); | |
578 | AliFatal("The 2 particles are different!"); | |
579 | } | |
580 | } | |
581 | } | |
582 | return isDuplicate; | |
583 | ||
584 | } | |
585 | ||
586 | Bool_t AliParticleYield::Compare2Floats(Float_t a, Float_t b) { | |
587 | // just a simple helper for the comparison methods | |
588 | if(!a) { | |
589 | if(!b) return kTRUE; // They are both 0; | |
590 | return kFALSE;// return here to avoid division by 0 | |
591 | } | |
592 | Bool_t areEqual = (TMath::Abs((a - b)/a) < fEpsilon); // If the relative difference is < epsilon, returns true | |
593 | if(!areEqual) { | |
594 | Printf("Warning: %f and %f are different", a,b); | |
595 | } | |
596 | return areEqual; | |
597 | } |