]>
Commit | Line | Data |
---|---|---|
e131b05f | 1 | #include "TH1D.h" |
2 | #include "TH2D.h" | |
3 | #include "TCanvas.h" | |
4 | #include "TFile.h" | |
5 | #include "TGraphAsymmErrors.h" | |
6 | #include "TLegend.h" | |
7 | ||
8 | #include "AliCFContainer.h" | |
9 | #include "AliCFEffGrid.h" | |
10 | #include "AliCFDataGrid.h" | |
11 | #include "AliPID.h" | |
12 | ||
13 | #include <iostream> | |
14 | ||
15 | #include "THnSparseDefinitions.h" | |
16 | ||
17 | enum type { kTrackPt = 0, kZ = 1, kXi = 2, kNtypes = 3 }; | |
18 | ||
19 | ||
20 | const Int_t nPtBinsType2 = 44; | |
21 | const Double_t binsPtType2[nPtBinsType2+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, | |
22 | 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, | |
23 | 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.4, 3.8, | |
24 | 4.5, 5.5, 6.5, 8.0, 10.0, 12.0, 14.0, 16.0, 20.0, 24.0, | |
25 | 28.0, 32.0, 36.0, 45.0, 50.0 }; | |
26 | ||
27 | ||
28 | //___________________________________________________________________ | |
29 | const Double_t* GetBins(Int_t type, Int_t& nBins) | |
30 | { | |
31 | if (type == 2) { | |
32 | nBins = nPtBinsType2; | |
33 | ||
34 | return binsPtType2; | |
35 | } | |
36 | ||
37 | return 0x0; | |
38 | } | |
39 | ||
40 | ||
41 | //___________________________________________________________________ | |
42 | Double_t trackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc) | |
43 | { | |
44 | return (1. - 0.129758 * TMath::Exp(-pTmc * 0.679612)); | |
45 | } | |
46 | ||
47 | ||
48 | //___________________________________________________________________ | |
49 | Double_t trackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc) | |
50 | { | |
51 | return TMath::Min((0.972865 + 0.0117093 * pTmc), 1.); | |
52 | } | |
53 | ||
54 | ||
55 | //___________________________________________________________________ | |
56 | Bool_t geantFlukaCorrection(AliCFContainer* data, Int_t genStepToDownscale) | |
57 | { | |
58 | // Issue: GEANT/FLUKA correction factor is for MC_pT. | |
59 | // Finally the effeciency should be DIVIDED by this correction factor. | |
60 | // To include resolution effects, it is therefore the best way to just | |
61 | // multiply the generated step with the correction factor. | |
62 | ||
63 | if (!data) { | |
64 | printf("No CFContainer for GEANT/FLUKA correction!\n"); | |
65 | return kFALSE; | |
66 | } | |
67 | ||
68 | const Int_t iPt = data->GetVar("P_{T} (GeV/c)"); | |
69 | const Int_t iMCid = data->GetVar("MC ID"); | |
70 | const Int_t iCharge = data->GetVar("Charge (e_{0})"); | |
71 | ||
72 | if (iPt < 0 || iMCid < 0 || iCharge < 0) { | |
73 | printf("Data axis for GEANT/FLUKA correction not found!\n"); | |
74 | return kFALSE; | |
75 | } | |
76 | ||
77 | if (!data->GetGrid(genStepToDownscale)) { | |
78 | printf("Step for downscaling (GEANT/FLUKA) not found!\n"); | |
79 | return kFALSE; | |
80 | } | |
81 | ||
82 | const Int_t nDim = data->GetNVar(); | |
83 | Int_t coord[nDim]; | |
84 | Double_t binCenterCoord[nDim]; | |
85 | ||
86 | Long64_t nBinsGrid = data->GetGrid(genStepToDownscale)->GetGrid()->GetNbins(); | |
87 | ||
88 | for (Long64_t iBin = 0; iBin < nBinsGrid; iBin++) { | |
89 | Double_t binContent = data->GetGrid(genStepToDownscale)->GetGrid()->GetBinContent(iBin, coord); | |
90 | Double_t binError = data->GetGrid(genStepToDownscale)->GetGrid()->GetBinError(iBin); | |
91 | ||
92 | for (Int_t iDim = 0; iDim < nDim; iDim++) | |
93 | binCenterCoord[iDim] = data->GetBinCenter(iDim, coord[iDim]); | |
94 | ||
95 | if (binCenterCoord[iCharge] < 0) { | |
96 | Double_t corrFactor = 1.; | |
97 | ||
98 | if (binCenterCoord[iMCid] - 0.5 == AliPID::kProton) | |
99 | corrFactor = trackingPtGeantFlukaCorrectionPrMinus(binCenterCoord[iPt]); | |
100 | else if (binCenterCoord[iMCid] - 0.5 == AliPID::kKaon) | |
101 | corrFactor = trackingPtGeantFlukaCorrectionKaMinus(binCenterCoord[iPt]); | |
102 | else | |
103 | continue; | |
104 | ||
105 | data->GetGrid(genStepToDownscale)->GetGrid()->SetBinContent(iBin, binContent * corrFactor); | |
106 | data->GetGrid(genStepToDownscale)->GetGrid()->SetBinError(iBin, binError * corrFactor); | |
107 | } | |
108 | } | |
109 | ||
110 | return kTRUE; | |
111 | } | |
112 | ||
113 | ||
114 | ||
115 | //___________________________________________________________________ | |
116 | void setupHist(TH1* h, TString histName, TString histTitle, TString xAxisTitle, TString yAxisTitle, Int_t color) | |
117 | { | |
118 | if (histName != "") | |
119 | h->SetName(histName.Data()); | |
120 | h->SetTitle(histTitle.Data()); | |
121 | ||
122 | if (xAxisTitle != "") | |
123 | h->GetXaxis()->SetTitle(xAxisTitle.Data()); | |
124 | if (yAxisTitle != "") | |
125 | h->GetYaxis()->SetTitle(yAxisTitle.Data()); | |
126 | ||
127 | h->SetMarkerStyle(24); | |
128 | h->SetLineColor(color); | |
129 | h->SetMarkerColor(color); | |
130 | ||
131 | h->SetStats(kFALSE); | |
132 | } | |
133 | ||
134 | ||
135 | //____________________________________________________________________________________________________________________ | |
136 | void undoJetNormalisationHist2D(TH2* hData, TH2* hNumJets, const Int_t lowerCentralityBinLimit, const Int_t upperCentralityBinLimit) | |
137 | { | |
138 | // Undo normalisation to 1/numJets. NOTE: jetPt binning of hData and hNumJets assumed to be the same! | |
139 | ||
140 | if (!hData || !hNumJets) | |
141 | return; | |
142 | ||
143 | for (Int_t binJetPt = 0; binJetPt <= hData->GetNbinsY() + 1; binJetPt++) { | |
144 | const Double_t numJets = hNumJets->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, binJetPt, binJetPt); | |
145 | Bool_t noJets = numJets < 1e-13; | |
146 | ||
147 | for (Int_t binObs = 0; binObs <= hData->GetNbinsX() + 1; binObs++) { | |
148 | if (noJets) | |
149 | continue; | |
150 | ||
151 | hData->SetBinContent(binObs, binJetPt, hData->GetBinContent(binObs, binJetPt) * numJets); | |
152 | hData->SetBinError(binObs, binJetPt, hData->GetBinError(binObs, binJetPt) * numJets); | |
153 | } | |
154 | } | |
155 | } | |
156 | ||
157 | ||
158 | //___________________________________________________________________ | |
159 | void normaliseHist(TH1* h, Double_t scale = 1.0) | |
160 | { | |
161 | if (h->GetSumw2N() <= 0) | |
162 | h->Sumw2(); | |
163 | ||
164 | h->Scale(scale); | |
165 | ||
166 | for (Int_t i = 1; i <= h->GetNbinsX(); i++) { | |
167 | Double_t normFactor = h->GetBinWidth(i); | |
168 | h->SetBinContent(i, h->GetBinContent(i) / normFactor); | |
169 | h->SetBinError(i, h->GetBinError(i) / normFactor); | |
170 | } | |
171 | } | |
172 | ||
173 | ||
174 | //___________________________________________________________________ | |
175 | void multiplyHistsDifferentBinning(TH1* h1, const TH1* h2, Double_t c1, Double_t c2, Bool_t ignoreErrorOfSecondHist = kFALSE) | |
176 | { | |
177 | Int_t nbinsx = h1->GetNbinsX(); | |
178 | Int_t nbinsy = h1->GetNbinsY(); | |
179 | Int_t nbinsz = h1->GetNbinsZ(); | |
180 | ||
181 | if (h1->GetDimension() < 2) | |
182 | nbinsy -= 1; | |
183 | if (h1->GetDimension() < 3) | |
184 | nbinsz -= 1; | |
185 | ||
186 | if (h1->GetSumw2N() == 0) | |
187 | h1->Sumw2(); | |
188 | ||
189 | Int_t bin, bin2, binx, biny, binz; | |
190 | Double_t b1,b2,w,d1,d2; | |
191 | d1 = c1*c1; | |
192 | d2 = c2*c2; | |
193 | ||
194 | // Loop over bins (including underflows/overflows) | |
195 | for (binz = 0; binz <= nbinsz + 1; binz++) { | |
196 | for (biny = 0; biny <= nbinsy + 1; biny++) { | |
197 | for (binx = 0; binx <= nbinsx + 1;binx++) { | |
198 | bin = binx + (nbinsx + 2) * (biny + (nbinsy + 2) * binz); | |
199 | bin2 = h2->FindFixBin(h1->GetXaxis()->GetBinCenter(binx), h1->GetYaxis()->GetBinCenter(biny), | |
200 | h1->GetZaxis()->GetBinCenter(binz)); | |
201 | b1 = h1->GetBinContent(bin); | |
202 | b2 = h2->GetBinContent(bin2); | |
203 | w = (c1*b1)*(c2*b2); | |
204 | h1->SetBinContent(bin, w); | |
205 | Double_t e1 = h1->GetBinError(bin); | |
206 | Double_t e2 = ignoreErrorOfSecondHist ? 0. : h2->GetBinError(bin2); | |
207 | h1->SetBinError(bin, TMath::Sqrt(d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1))); | |
208 | } | |
209 | } | |
210 | } | |
211 | } | |
212 | ||
213 | ||
214 | //___________________________________________________________________ | |
215 | void divideHistsDifferentBinning(TH1* h1, const TH1* h2, Double_t c1, Double_t c2, Bool_t ignoreErrorOfSecondHist = kFALSE) | |
216 | { | |
217 | Int_t nbinsx = h1->GetNbinsX(); | |
218 | Int_t nbinsy = h1->GetNbinsY(); | |
219 | Int_t nbinsz = h1->GetNbinsZ(); | |
220 | ||
221 | if (h1->GetDimension() < 2) | |
222 | nbinsy -= 1; | |
223 | if (h1->GetDimension() < 3) | |
224 | nbinsz -= 1; | |
225 | ||
226 | if (h1->GetSumw2N() == 0) | |
227 | h1->Sumw2(); | |
228 | ||
229 | Int_t bin, bin2, binx, biny, binz; | |
230 | Double_t b1,b2,w,d1,d2; | |
231 | d1 = c1*c1; | |
232 | d2 = c2*c2; | |
233 | ||
234 | // Loop over bins (including underflows/overflows) | |
235 | for (binz = 0; binz <= nbinsz + 1; binz++) { | |
236 | for (biny = 0; biny <= nbinsy + 1; biny++) { | |
237 | for (binx = 0; binx <= nbinsx + 1;binx++) { | |
238 | bin = binx + (nbinsx + 2) * (biny + (nbinsy + 2) * binz); | |
239 | bin2 = h2->FindFixBin(h1->GetXaxis()->GetBinCenter(binx), h1->GetYaxis()->GetBinCenter(biny), | |
240 | h1->GetZaxis()->GetBinCenter(binz)); | |
241 | b1 = h1->GetBinContent(bin); | |
242 | b2 = h2->GetBinContent(bin2); | |
243 | if (b2) | |
244 | w = (c1*b1)/(c2*b2); | |
245 | else | |
246 | w = 0; | |
247 | h1->SetBinContent(bin, w); | |
248 | ||
249 | if (!b2) { | |
250 | h1->SetBinError(bin, 0); | |
251 | continue; | |
252 | } | |
253 | ||
254 | Double_t b22 = b2*b2*d2; | |
255 | Double_t e1 = h1->GetBinError(bin); | |
256 | Double_t e2 = ignoreErrorOfSecondHist ? 0. : h2->GetBinError(bin2); | |
257 | h1->SetBinError(bin, TMath::Sqrt(d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22))); | |
258 | } | |
259 | } | |
260 | } | |
261 | } | |
262 | ||
263 | ||
264 | //___________________________________________________________________ | |
265 | // Efficiency for inclusive spectra vs. pT (or also jets, but without using the generation) | |
266 | // E.g. a 'calcEfficiency.C+("finalCuts/MC_pp/7TeV/LHC10f6a/corrected/finalisedSplines/analytical/efficiency_noClCut_preliminary/bhess_PID_Jets_Inclusive_PureGauss_efficiency.root", "finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/outputSystematicsTotal_SummedSystematicErrors__2013_10_21.root", kTRUE, kFALSE, 0, -2, -2, -1, -1, 8, 2)' -b -q | |
267 | Int_t calcEfficiency(TString pathNameEfficiency, TString pathNameData, Bool_t correctGeantFluka, Bool_t scaleStrangeness, | |
268 | Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/, | |
269 | Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/, | |
270 | Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/, | |
271 | Double_t constantCorrectionAbovePtThreshold, | |
272 | Int_t rebinEfficiency) | |
273 | { | |
274 | TString pathData = pathNameData; | |
275 | pathData.Replace(pathData.Last('/'), pathData.Length(), ""); | |
276 | ||
277 | TFile* fileEff = TFile::Open(pathNameEfficiency.Data()); | |
278 | if (!fileEff) { | |
279 | printf("Failed to open efficiency file \"%s\"\n", pathNameEfficiency.Data()); | |
280 | return -1; | |
281 | } | |
282 | ||
283 | AliCFContainer* data = (AliCFContainer*)(fileEff->Get("containerEff")); | |
284 | if (!data) { | |
285 | printf("Failed to load efficiency container!\n"); | |
286 | return -1; | |
287 | } | |
288 | ||
289 | const Int_t iMCid = data->GetVar("MC ID"); | |
290 | const Int_t iPt = data->GetVar("P_{T} (GeV/c)"); | |
291 | const Int_t iEta = data->GetVar("#eta"); | |
292 | const Int_t iCharge = data->GetVar("Charge (e_{0})"); | |
293 | const Int_t iMult = data->GetVar("Centrality Percentile"); | |
294 | const Int_t iJetPt = data->GetVar("P_{T}^{jet} (GeV/c)"); | |
295 | //const Int_t iZ = data->GetVar("z = P_{T}^{track} / P_{T}^{jet}"); | |
296 | //const Int_t iXi = data->GetVar("#xi = ln(P_{T}^{jet} / P_{T}^{track})"); | |
297 | ||
298 | TFile* fileData = TFile::Open(pathNameData.Data()); | |
299 | if (!fileData) { | |
300 | printf("Failed to open data file \"%s\"\n", pathNameData.Data()); | |
301 | return -1; | |
302 | } | |
303 | ||
304 | TH1D* hYield[AliPID::kSPECIES] = { 0x0, }; | |
305 | TH1D* hYieldSysError[AliPID::kSPECIES] = { 0x0, }; | |
306 | TGraphAsymmErrors* gYieldSysError[AliPID::kSPECIES] = { 0x0, }; | |
307 | TH1D* hMCgenPrimYield[AliPID::kSPECIES] = { 0x0, }; | |
308 | Int_t numMCgenPrimYieldHistsFound = 0; | |
309 | ||
310 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
311 | TString speciesName = AliPID::ParticleName(species); | |
312 | TString firstLetter = speciesName(0); | |
313 | firstLetter.ToUpper(); | |
314 | speciesName.Replace(0, 1, firstLetter.Data()); | |
315 | TString histName = Form("hYield%ss", speciesName.Data()); | |
316 | hYield[species] = (TH1D*)fileData->Get(histName.Data()); | |
317 | if (!hYield[species]) { | |
318 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
319 | return -1; | |
320 | } | |
321 | hYield[species]->SetFillStyle(0); | |
322 | ||
323 | TString graphName = Form("systematicErrorYields_%s", AliPID::ParticleName(species)); | |
324 | gYieldSysError[species] = (TGraphAsymmErrors*)fileData->Get(graphName.Data()); | |
325 | ||
326 | // In case of MC also retrieve the MC truth generated yields | |
327 | TString histNameMCgenYields = Form("hMCgenYieldsPrimSpecies_%s", AliPID::ParticleShortName(species)); | |
328 | hMCgenPrimYield[species] = (TH1D*)fileData->Get(histNameMCgenYields.Data()); | |
329 | if (hMCgenPrimYield[species]) | |
330 | numMCgenPrimYieldHistsFound++; | |
331 | } | |
332 | ||
333 | if (numMCgenPrimYieldHistsFound > 0 && numMCgenPrimYieldHistsFound != AliPID::kSPECIES) { | |
334 | printf("Error: Unable to retrieve all MC generated prim yield histos! Got %d.\n", numMCgenPrimYieldHistsFound); | |
335 | return -1; | |
336 | } | |
337 | ||
338 | TCanvas* cFractions = (TCanvas*)fileData->Get("cFractionsWithTotalSystematicError"); | |
339 | if (!cFractions) | |
340 | cFractions = (TCanvas*)fileData->Get("cFractions"); | |
341 | ||
342 | ||
343 | // Convert graphs with systematic errors into histograms (assume symmetric error, which is currently true) | |
344 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
345 | if (gYieldSysError[species]) { | |
346 | hYieldSysError[species] = new TH1D(*hYield[species]); | |
347 | hYieldSysError[species]->SetName(Form("%s_sysError", hYieldSysError[species]->GetName())); | |
348 | hYieldSysError[species]->SetFillStyle(0); | |
349 | ||
350 | for (Int_t binX = 1; binX <= hYieldSysError[species]->GetNbinsX(); binX++) { | |
351 | hYieldSysError[species]->SetBinContent(binX, gYieldSysError[species]->GetY()[binX - 1]); | |
352 | hYieldSysError[species]->SetBinError(binX, gYieldSysError[species]->GetErrorY(binX - 1)); | |
353 | } | |
354 | } | |
355 | } | |
356 | ||
357 | // Take pT binning from pion yield (binning for all species the same) and create a new AliCFContainer with this new binning for pT | |
358 | TH1D hDummy(*hYield[AliPID::kPion]); | |
359 | hDummy.SetName("hDummy"); | |
360 | TAxis* axis = 0x0; | |
361 | ||
362 | if (rebinEfficiency > 1) { | |
363 | Int_t nBinsNew = 0; | |
364 | const Double_t* newBins = GetBins(rebinEfficiency, nBinsNew); | |
365 | ||
366 | hDummy.SetBins(nBinsNew, newBins); | |
367 | ||
368 | axis = hDummy.GetXaxis(); | |
369 | ||
370 | //axis = hDummy.Rebin(rebinEfficiency, "", 0)->GetXaxis(); | |
371 | } | |
372 | else | |
373 | axis = hDummy.GetXaxis(); | |
374 | ||
375 | const TArrayD* binsPtCurrent = axis->GetXbins(); | |
376 | TArrayD* binsPtNew = new TArrayD(*binsPtCurrent); | |
377 | ||
378 | const Int_t nEffDims = data->GetNVar(); | |
379 | Int_t nEffBins[nEffDims]; | |
380 | ||
381 | for (Int_t iDim = 0; iDim < nEffDims; iDim++) { | |
382 | if (iDim == iPt) | |
383 | nEffBins[iDim] = axis->GetNbins(); | |
384 | else | |
385 | nEffBins[iDim] = data->GetNBins(iDim); | |
386 | } | |
387 | ||
388 | ||
389 | // Just make one large pT bin above some threshold, if desired | |
390 | if (binsPtNew->fN != 0 && constantCorrectionAbovePtThreshold > 0) { | |
391 | for (Int_t iBin = 0; iBin < nEffBins[iPt]; iBin++) { | |
392 | // Find the first bin edged really larger than the threshold. | |
393 | // If the bin edge before equals the threshold, just set the | |
394 | // current bin edge to the right end of the spectrum -> Done. | |
395 | // If the bin edge before is different, set the bin edge to the | |
396 | // threshold | |
397 | if (binsPtNew->fArray[iBin] > constantCorrectionAbovePtThreshold) { | |
398 | if (binsPtNew->fArray[iBin - 1] == constantCorrectionAbovePtThreshold) { | |
399 | binsPtNew->fArray[iBin] = binsPtNew->fArray[nEffBins[iPt]]; | |
400 | nEffBins[iPt] = iBin; | |
401 | break; | |
402 | } | |
403 | else { | |
404 | binsPtNew->fArray[iBin] = constantCorrectionAbovePtThreshold; | |
405 | } | |
406 | } | |
407 | } | |
408 | } | |
409 | ||
410 | ||
411 | AliCFContainer *dataRebinned = new AliCFContainer(Form("%s_rebinned", data->GetName()), Form("%s (rebinned)", data->GetTitle()), | |
412 | data->GetNStep(), nEffDims, nEffBins); | |
413 | ||
414 | for (Int_t iDim = 0; iDim < nEffDims; iDim++) { | |
415 | dataRebinned->SetVarTitle(iDim, data->GetVarTitle(iDim)); | |
416 | ||
417 | if (iDim == iPt) { | |
418 | if (binsPtNew->fN == 0) | |
419 | dataRebinned->SetBinLimits(iDim, axis->GetXmin(), axis->GetXmax()); | |
420 | else | |
421 | dataRebinned->SetBinLimits(iDim, binsPtNew->fArray); | |
422 | } | |
423 | else { | |
424 | dataRebinned->SetBinLimits(iDim, data->GetBinLimits(iDim)); | |
425 | } | |
426 | } | |
427 | ||
428 | for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++) | |
429 | dataRebinned->SetStepTitle(iStep, data->GetStepTitle(iStep)); | |
430 | ||
431 | Int_t coord[nEffDims]; | |
432 | Double_t binCenterCoord[nEffDims]; | |
433 | ||
434 | // Fill content from old grid into the new grid with proper binning | |
435 | for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++) { | |
436 | Long64_t nBinsGrid = data->GetGrid(iStep)->GetGrid()->GetNbins(); | |
437 | ||
438 | for (Long64_t iBin = 0; iBin < nBinsGrid; iBin++) { | |
439 | Double_t binContent = data->GetGrid(iStep)->GetGrid()->GetBinContent(iBin, coord); | |
440 | Double_t binError2 = data->GetGrid(iStep)->GetGrid()->GetBinError2(iBin); | |
441 | ||
442 | for (Int_t iDim = 0; iDim < nEffDims; iDim++) { | |
443 | binCenterCoord[iDim] = data->GetBinCenter(iDim, coord[iDim]); | |
444 | } | |
445 | ||
446 | Long64_t iBinRebinned = dataRebinned->GetGrid(iStep)->GetGrid()->GetBin(binCenterCoord); | |
447 | dataRebinned->GetGrid(iStep)->GetGrid()->AddBinContent(iBinRebinned, binContent); | |
448 | dataRebinned->GetGrid(iStep)->GetGrid()->AddBinError2(iBinRebinned, binError2); | |
449 | } | |
450 | } | |
451 | ||
452 | ||
453 | // If desired, restrict centrality axis | |
454 | Int_t lowerCentralityBinLimit = -1; | |
455 | Int_t upperCentralityBinLimit = -2; // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept | |
456 | Bool_t restrictCentralityAxis = kFALSE; | |
457 | Double_t actualLowerCentrality = -1.; | |
458 | Double_t actualUpperCentrality = -1.; | |
459 | ||
460 | if (lowerCentrality >= -1 && upperCentrality >= -1) { | |
461 | // Add subtract a very small number to avoid problems with values right on the border between to bins | |
462 | lowerCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(lowerCentrality + 0.001); | |
463 | upperCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(upperCentrality - 0.001); | |
464 | ||
465 | // Check if the values look reasonable | |
466 | if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1 | |
467 | && upperCentralityBinLimit <= dataRebinned->GetAxis(iMult, 0)->GetNbins()) { | |
468 | actualLowerCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinLowEdge(lowerCentralityBinLimit); | |
469 | actualUpperCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinUpEdge(upperCentralityBinLimit); | |
470 | ||
471 | restrictCentralityAxis = kTRUE; | |
472 | } | |
473 | else { | |
474 | std::cout << std::endl; | |
475 | std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl; | |
476 | return -1; | |
477 | } | |
478 | } | |
479 | ||
480 | std::cout << "centrality: "; | |
481 | if (restrictCentralityAxis) { | |
482 | std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl; | |
483 | dataRebinned->SetRangeUser(iMult, lowerCentralityBinLimit, upperCentralityBinLimit, kTRUE); | |
484 | } | |
485 | else { | |
486 | std::cout << "All" << std::endl; | |
487 | } | |
488 | ||
489 | ||
490 | ||
491 | // If desired, restrict jetPt axis | |
492 | Int_t lowerJetPtBinLimit = -1; | |
493 | Int_t upperJetPtBinLimit = -1; | |
494 | Bool_t restrictJetPtAxis = kFALSE; | |
495 | Double_t actualLowerJetPt = -1.; | |
496 | Double_t actualUpperJetPt = -1.; | |
497 | ||
498 | if (lowerJetPt >= 0 && upperJetPt >= 0) { | |
499 | // Add subtract a very small number to avoid problems with values right on the border between to bins | |
500 | lowerJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(lowerJetPt + 0.001); | |
501 | upperJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(upperJetPt - 0.001); | |
502 | ||
503 | // Check if the values look reasonable | |
504 | if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 && | |
505 | upperJetPtBinLimit <= dataRebinned->GetAxis(iJetPt, 0)->GetNbins()) { | |
506 | actualLowerJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinLowEdge(lowerJetPtBinLimit); | |
507 | actualUpperJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinUpEdge(upperJetPtBinLimit); | |
508 | ||
509 | restrictJetPtAxis = kTRUE; | |
510 | } | |
511 | else { | |
512 | std::cout << std::endl; | |
513 | std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl; | |
514 | return -1; | |
515 | } | |
516 | } | |
517 | ||
518 | std::cout << "jet pT: "; | |
519 | if (restrictJetPtAxis) { | |
520 | std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl; | |
521 | dataRebinned->SetRangeUser(iJetPt, lowerJetPtBinLimit, upperJetPtBinLimit, kTRUE); | |
522 | } | |
523 | else { | |
524 | std::cout << "All" << std::endl; | |
525 | } | |
526 | ||
527 | // If desired, restrict charge axis | |
528 | std::cout << "Charge selection (efficiency): "; | |
529 | if (chargeMode == kAllCharged) | |
530 | std::cout << "All charged particles" << std::endl; | |
531 | else if (chargeMode == kNegCharge) | |
532 | std::cout << "Negative particles only" << std::endl; | |
533 | else if (chargeMode == kPosCharge) | |
534 | std::cout << "Positive particles only" << std::endl; | |
535 | else { | |
536 | std::cout << "Unknown -> ERROR" << std::endl; | |
537 | return -1; | |
538 | } | |
539 | ||
540 | const Bool_t restrictCharge = (chargeMode != kAllCharged); | |
541 | ||
542 | Int_t lowerChargeBinLimit = -1; | |
543 | Int_t upperChargeBinLimit = -2; | |
544 | Double_t actualLowerCharge = -999; | |
545 | Double_t actualUpperCharge = -999; | |
546 | ||
547 | if (restrictCharge) { | |
548 | // Add subtract a very small number to avoid problems with values right on the border between to bins | |
549 | if (chargeMode == kNegCharge) { | |
550 | lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(-1. + 0.001); | |
551 | upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. - 0.001); | |
552 | } | |
553 | else if (chargeMode == kPosCharge) { | |
554 | lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. + 0.001); | |
555 | upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(1. - 0.001); | |
556 | } | |
557 | ||
558 | // Check if the values look reasonable | |
559 | if (lowerChargeBinLimit <= upperChargeBinLimit && lowerChargeBinLimit >= 1 | |
560 | && upperChargeBinLimit <= dataRebinned->GetAxis(iCharge, 0)->GetNbins()) { | |
561 | actualLowerCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinLowEdge(lowerChargeBinLimit); | |
562 | actualUpperCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinUpEdge(upperChargeBinLimit); | |
563 | ||
564 | std::cout << "Charge range (efficiency): " << actualLowerCharge << " - " << actualUpperCharge << std::endl; | |
565 | } | |
566 | else { | |
567 | std::cout << std::endl; | |
568 | std::cout << "Requested charge range (efficiency) out of limits or upper and lower limit are switched!" << std::endl; | |
569 | return -1; | |
570 | } | |
571 | ||
572 | dataRebinned->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE); | |
573 | data->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE); | |
574 | } | |
575 | ||
576 | std::cout << std::endl; | |
577 | ||
578 | ||
579 | // If jet axis is restricted (i.e. jet input), also need num jet histos for proper normalisation. | |
580 | // Load numJet histos from bhess_PID*.root file related to efficiency file. | |
581 | TH2D* hNjetsGen = 0x0; | |
582 | TH2D* hNjetsRec = 0x0; | |
583 | ||
584 | if (restrictJetPtAxis) { | |
585 | TString pathNameDataMC = pathNameEfficiency; | |
586 | pathNameDataMC.ReplaceAll("_efficiency", ""); | |
587 | ||
588 | TFile* fDataMC = TFile::Open(pathNameDataMC.Data()); | |
589 | if (!fDataMC) { | |
590 | std::cout << std::endl; | |
591 | std::cout << "Failed to open file \"" << pathNameData.Data() << "\" to obtain num of rec/gen jets!" << std::endl; | |
592 | ||
593 | return -1; | |
594 | } | |
595 | ||
596 | TString listName = pathNameDataMC; | |
597 | listName.Replace(0, listName.Last('/') + 1, ""); | |
598 | listName.ReplaceAll(".root", ""); | |
599 | ||
600 | TObjArray* histList = (TObjArray*)(fDataMC->Get(listName.Data())); | |
601 | if (!histList) { | |
602 | std::cout << std::endl; | |
603 | std::cout << "Failed to load list \"" << listName.Data() << "\" to obtain num of rec/gen jets!" << std::endl; | |
604 | return -1; | |
605 | } | |
606 | ||
607 | hNjetsGen = (TH2D*)histList->FindObject("fh2FFJetPtGen"); | |
608 | hNjetsRec = (TH2D*)histList->FindObject("fh2FFJetPtRec"); | |
609 | ||
610 | if (!hNjetsRec || !hNjetsGen) { | |
611 | std::cout << "Failed to load number of jets histos!" << std::endl; | |
612 | ||
613 | ||
614 | // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other | |
615 | // period is considered; also: No multiplicity information) | |
616 | TString pathEfficiency = pathNameEfficiency; | |
617 | pathEfficiency.Replace(pathEfficiency.Last('/'), pathEfficiency.Length(), ""); | |
618 | TString pathBackward = Form("%s/AnalysisResults.root", pathEfficiency.Data()); | |
619 | TFile* fBackward = TFile::Open(pathBackward.Data()); | |
620 | ||
621 | TString dirDataInFile = ""; | |
622 | TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0; | |
623 | ||
624 | TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0; | |
625 | ||
626 | TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCuts") : 0x0; | |
627 | TH1D* hFFJetPtGen = list ? (TH1D*)list->FindObject("fh1FFJetPtGen") : 0x0; | |
628 | ||
629 | if (hFFJetPtRec && hFFJetPtGen) { | |
630 | printf("***WARNING: For backward compatibility, using file \"%s\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n", | |
631 | pathBackward.Data()); | |
632 | ||
633 | hNjetsRec = new TH2D("fh2FFJetPtRec", "", 1, -1, 1, dataRebinned->GetAxis(iJetPt, 0)->GetNbins(), | |
634 | dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray()); | |
635 | ||
636 | for (Int_t iJet = 1; iJet <= hNjetsRec->GetNbinsY(); iJet++) { | |
637 | Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinLowEdge(iJet) + 1e-3); | |
638 | Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinUpEdge(iJet) - 1e-3); | |
639 | hNjetsRec->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin)); | |
640 | } | |
641 | ||
642 | hNjetsGen = new TH2D("fh2FFJetPtGen", "", 1, -1, 1, dataRebinned->GetAxis(iJetPt, 0)->GetNbins(), | |
643 | dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray()); | |
644 | ||
645 | for (Int_t iJet = 1; iJet <= hNjetsGen->GetNbinsY(); iJet++) { | |
646 | Int_t lowerBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinLowEdge(iJet) + 1e-3); | |
647 | Int_t upperBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinUpEdge(iJet) - 1e-3); | |
648 | hNjetsGen->SetBinContent(1, iJet, hFFJetPtGen->Integral(lowerBin, upperBin)); | |
649 | } | |
650 | } | |
651 | ||
652 | if (!hNjetsRec || ! hNjetsGen) | |
653 | return -1; | |
654 | } | |
655 | } | |
656 | ||
657 | // For normalisation to number of jets | |
658 | // NOTE: These numbers are for the efficiency only! The data will be normalised to its own number!!! | |
659 | const Double_t nJetsGen = hNjetsGen ? hNjetsGen->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit, | |
660 | upperJetPtBinLimit) : 1.; | |
661 | const Double_t nJetsRec = hNjetsRec ? hNjetsRec->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit, | |
662 | upperJetPtBinLimit) : 1.; | |
663 | ||
664 | // Secondary correction | |
665 | AliCFEffGrid* sec = new AliCFEffGrid("sec", "Secondary Contamination", *dataRebinned); | |
666 | AliCFEffGrid* secStrangeScale = new AliCFEffGrid("secStrangeScale", "Secondary Contamination with Strangeness Scaling", | |
667 | *dataRebinned); | |
668 | ||
669 | // Either one can take kStepRecWithGenCutsMeasuredObs or, what I prefer, one can take | |
670 | // kStepRecWithRecCutsMeasuredObsPrimaries => The difference is only the eta cut, which is on the rec level | |
671 | // in the latter case, i.e. one corrects for eta resolution (although the effect should be very small) | |
672 | // => TESTED: There is NO difference (only 1 bin vs. pt shows a deviation by one entry), so this cut has, | |
673 | // in principle, more or less no effect | |
674 | ||
675 | // For testing with data set w/o strangeness secStrangeScale->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObs); | |
676 | secStrangeScale->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObsStrangenessScaled); | |
677 | sec->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObs); | |
678 | ||
679 | // QA plots for secondaries | |
680 | TCanvas* cSec = new TCanvas("cSec", "Secondary Contamination", 100, 10, 1200, 800); | |
681 | cSec->Divide(2, 1); | |
682 | cSec->GetPad(1)->SetLogx(kTRUE); | |
683 | cSec->cd(1); | |
684 | TH1* hSecPt = sec->Project(iPt); | |
685 | hSecPt->SetName("hSecPt"); | |
686 | hSecPt->SetStats(kFALSE); | |
687 | hSecPt->GetXaxis()->SetRangeUser(0.15, 50.); | |
688 | hSecPt->GetXaxis()->SetMoreLogLabels(kTRUE); | |
689 | hSecPt->GetXaxis()->SetNoExponent(kTRUE); | |
690 | hSecPt->GetYaxis()->SetTitle("Primary Fraction"); | |
691 | hSecPt->Draw("E1"); | |
692 | cSec->cd(2); | |
693 | TH1* hEtaSec = sec->Project(iEta); | |
694 | hEtaSec->SetName("hEtaSec"); | |
695 | hEtaSec->SetStats(kFALSE); | |
696 | hEtaSec->GetYaxis()->SetTitle("Primary Fraction"); | |
697 | hEtaSec->Draw("E1"); | |
698 | TH2D* hSecID2Pt = (TH2D*)sec->Project(iPt, iMCid); | |
699 | hSecID2Pt->SetName("hSecID2Pt"); | |
700 | ||
701 | // Get the secondary contamination vs. pT for each species | |
702 | TH1D* hSec[AliPID::kSPECIES]; | |
703 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
704 | hSec[species] = hSecID2Pt->ProjectionX(Form("hSec_%s", AliPID::ParticleShortName(species)), species + 1, species + 1, "e"); | |
705 | hSec[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
706 | hSec[species]->SetLineColor(hYield[species]->GetLineColor()); | |
707 | hSec[species]->SetMarkerColor(hYield[species]->GetLineColor()); | |
708 | hSec[species]->SetLineWidth(2.); | |
709 | hSec[species]->GetXaxis()->SetRangeUser(0.15, 50); | |
710 | hSec[species]->GetYaxis()->SetRangeUser(0., 1.01); | |
711 | hSec[species]->GetXaxis()->SetMoreLogLabels(kTRUE); | |
712 | hSec[species]->GetXaxis()->SetNoExponent(kTRUE); | |
713 | hSec[species]->SetStats(kFALSE); | |
714 | hSec[species]->GetYaxis()->SetTitle("Primary Fraction"); | |
715 | } | |
716 | ||
717 | TCanvas* cSec2 = new TCanvas("cSec2", "Primary fraction for different species", 100, 10, 1200, 800); | |
718 | cSec2->SetGridx(1); | |
719 | cSec2->SetGridy(1); | |
720 | cSec2->SetLogx(1); | |
721 | ||
722 | hSec[0]->Draw("E1"); | |
723 | ||
724 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) { | |
725 | hSec[i]->Draw("E1 same"); | |
726 | } | |
727 | cSec2->BuildLegend()->SetFillColor(kWhite); | |
728 | ||
729 | ClearTitleFromHistoInCanvas(cSec2); | |
730 | ||
731 | // QA plots for secondaries with strangeness scaling | |
732 | TCanvas* cSecSS = new TCanvas("cSecSS", "Secondary Contamination (strangeness scaled)", 100, 10, 1200, 800); | |
733 | cSecSS->Divide(2, 1); | |
734 | cSecSS->GetPad(1)->SetLogx(kTRUE); | |
735 | cSecSS->cd(1); | |
736 | TH1* hSecSSPt = secStrangeScale->Project(iPt); | |
737 | hSecSSPt->SetName("hSecSSPt"); | |
738 | hSecSSPt->SetStats(kFALSE); | |
739 | hSecSSPt->GetXaxis()->SetRangeUser(0.15, 50.); | |
740 | hSecSSPt->GetXaxis()->SetMoreLogLabels(kTRUE); | |
741 | hSecSSPt->GetXaxis()->SetNoExponent(kTRUE); | |
742 | hSecSSPt->GetYaxis()->SetTitle("Primary Fraction"); | |
743 | hSecSSPt->Draw("E1"); | |
744 | cSecSS->cd(2); | |
745 | TH1* hEtaSecSS = secStrangeScale->Project(iEta); | |
746 | hEtaSecSS->SetName("hEtaSecSS"); | |
747 | hEtaSecSS->SetStats(kFALSE); | |
748 | hEtaSecSS->GetYaxis()->SetTitle("Primary Fraction"); | |
749 | hEtaSecSS->Draw("E1"); | |
750 | TH2D* hSecSSID2Pt = (TH2D*)secStrangeScale->Project(iPt, iMCid); | |
751 | hSecSSID2Pt->SetName("hSecSSID2Pt"); | |
752 | ||
753 | // Get the secondary contamination vs. pT for each species | |
754 | TH1D* hSecSS[AliPID::kSPECIES]; | |
755 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
756 | hSecSS[species] = hSecSSID2Pt->ProjectionX(Form("hSecSS_%s", AliPID::ParticleShortName(species)), species + 1, species + 1, "e"); | |
757 | hSecSS[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
758 | hSecSS[species]->SetLineColor(hYield[species]->GetLineColor()); | |
759 | hSecSS[species]->SetMarkerColor(hYield[species]->GetLineColor()); | |
760 | hSecSS[species]->SetLineWidth(2.); | |
761 | hSecSS[species]->GetXaxis()->SetRangeUser(0.15, 50); | |
762 | hSecSS[species]->GetYaxis()->SetRangeUser(0., 1.01); | |
763 | hSecSS[species]->GetXaxis()->SetMoreLogLabels(kTRUE); | |
764 | hSecSS[species]->GetXaxis()->SetNoExponent(kTRUE); | |
765 | hSecSS[species]->SetStats(kFALSE); | |
766 | hSecSS[species]->GetYaxis()->SetTitle("Primary Fraction"); | |
767 | } | |
768 | ||
769 | TCanvas* cSecSS2 = new TCanvas("cSecSS2", "Primary fraction for different species (strangeness scaled)", 100, 10, 1200, 800); | |
770 | cSecSS2->SetGridx(1); | |
771 | cSecSS2->SetGridy(1); | |
772 | cSecSS2->SetLogx(1); | |
773 | ||
774 | hSecSS[0]->Draw("E1"); | |
775 | ||
776 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) { | |
777 | hSecSS[i]->Draw("E1 same"); | |
778 | } | |
779 | cSecSS2->BuildLegend()->SetFillColor(kWhite); | |
780 | ||
781 | ClearTitleFromHistoInCanvas(cSecSS2); | |
782 | ||
783 | ||
784 | // Secondary correction for to-pi-ratios | |
785 | TH1D* hSecToPiRatio[AliPID::kSPECIES] = { 0x0, }; | |
786 | TH1D* hSecSSToPiRatio[AliPID::kSPECIES] = { 0x0, }; | |
787 | ||
788 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
789 | if (species == AliPID::kPion) | |
790 | continue; // Do not consider pion-to-pion ratio | |
791 | ||
792 | hSecToPiRatio[species] = new TH1D(*hSec[species]); | |
793 | hSecToPiRatio[species]->Reset(); | |
794 | hSecToPiRatio[species]->SetName(Form("hSecSSToPionRatio_%s", AliPID::ParticleShortName(species))); | |
795 | hSecToPiRatio[species]->SetTitle(Form("%s/#pi", AliPID::ParticleLatexName(species))); | |
796 | hSecToPiRatio[species]->SetLineColor(getLineColorAliPID(species)); | |
797 | hSecToPiRatio[species]->SetMarkerColor(getLineColorAliPID(species)); | |
798 | hSecToPiRatio[species]->SetMarkerStyle(24); | |
799 | hSecToPiRatio[species]->SetLineWidth(2.); | |
800 | hSecToPiRatio[species]->SetStats(kFALSE); | |
801 | hSecToPiRatio[species]->GetYaxis()->SetTitle("Primary Fraction of Ratio"); | |
802 | ||
803 | hSecSSToPiRatio[species] = new TH1D(*hSecToPiRatio[species]); | |
804 | hSecSSToPiRatio[species]->SetName(Form("hSecSSToPionRatio_%s", AliPID::ParticleShortName(species))); | |
805 | ||
806 | // Samples for different species are independent, so just divide correction factors | |
807 | hSecToPiRatio[species]->Divide(hSec[species], hSec[AliPID::kPion], 1., 1., ""); | |
808 | hSecToPiRatio[species]->GetYaxis()->SetRangeUser(0., 1.1); | |
809 | ||
810 | hSecSSToPiRatio[species]->Divide(hSecSS[species], hSecSS[AliPID::kPion], 1., 1., ""); | |
811 | hSecSSToPiRatio[species]->GetYaxis()->SetRangeUser(0., 1.1); | |
812 | } | |
813 | ||
814 | TCanvas* cSecToPiRatio = new TCanvas("cSecToPiRatio", "Primary fraction of to-#pi-ratio for different species", 100, 10, 1200, 800); | |
815 | cSecToPiRatio->SetGridx(1); | |
816 | cSecToPiRatio->SetGridy(1); | |
817 | cSecToPiRatio->SetLogx(1); | |
818 | ||
819 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
820 | if (i == AliPID::kPion) | |
821 | continue; | |
822 | ||
823 | hSecToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same")); | |
824 | } | |
825 | ||
826 | cSecToPiRatio->BuildLegend()->SetFillColor(kWhite); | |
827 | ||
828 | ClearTitleFromHistoInCanvas(cSecToPiRatio); | |
829 | ||
830 | ||
831 | ||
832 | TCanvas* cSecSSToPiRatio = new TCanvas("cSecSSToPiRatio", | |
833 | "Primary fraction of to-#pi-ratio for different species (strangeness scaled)", | |
834 | 100, 10, 1200, 800); | |
835 | cSecSSToPiRatio->SetGridx(1); | |
836 | cSecSSToPiRatio->SetGridy(1); | |
837 | cSecSSToPiRatio->SetLogx(1); | |
838 | ||
839 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
840 | if (i == AliPID::kPion) | |
841 | continue; | |
842 | ||
843 | hSecSSToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same")); | |
844 | } | |
845 | ||
846 | cSecSSToPiRatio->BuildLegend()->SetFillColor(kWhite); | |
847 | ||
848 | ClearTitleFromHistoInCanvas(cSecSSToPiRatio); | |
849 | ||
850 | ||
851 | ||
852 | ||
853 | // If desired, apply the GEANT/FLUKA correction first | |
854 | // NOTE: This will change dataRebinned! If anything else should also be calculated, it must be done before. | |
855 | // Otherwise, the GEANT/FLUKA correction factor for the efficiency will also affect it! | |
856 | const Int_t genStepEff = kStepGenWithGenCuts; | |
857 | if (correctGeantFluka) { | |
858 | printf("Applying GEANT/FLUKA correction...\n"); | |
859 | if (!geantFlukaCorrection(dataRebinned, genStepEff)) { | |
860 | printf("GEANT/FLUKA correction could not be applied!\n"); | |
861 | return kFALSE; | |
862 | } | |
863 | } | |
864 | ||
865 | // Construct the efficiency grid from the data container | |
866 | AliCFEffGrid* eff = new AliCFEffGrid("eff", "Efficiency x Acceptance x pT Resolution", *dataRebinned); | |
867 | ||
868 | // Either one can take kStepRecWithGenCutsMeasuredObs or, what I prefer, one can take | |
869 | // kStepRecWithRecCutsMeasuredObsPrimaries => The difference is only the eta cut, which is on the rec level | |
870 | // in the latter case, i.e. one corrects for eta resolution (although the effect should be very small) | |
871 | eff->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, genStepEff); | |
872 | ||
873 | // If the jet axis is restricted (i.e. jet input), scale with the corresponding number of jets. | |
874 | // Note: Since this is supposed to be a real scaling, set the error of the scale factor to zero | |
875 | // (second element of factor array) | |
876 | if (restrictJetPtAxis) { | |
877 | Double_t factor_Numerator[2] = { nJetsRec > 0 ? 1. / nJetsRec : 0., 0. }; | |
878 | Double_t factor_Denominator[2] = { nJetsGen > 0 ? 1. / nJetsGen : 0., 0. }; | |
879 | eff->GetNum()->Scale(factor_Numerator); | |
880 | eff->GetDen()->Scale(factor_Denominator); | |
881 | } | |
882 | ||
883 | //The efficiency along pt and vertex, and 2-D projection | |
884 | TCanvas* cEff = new TCanvas("cEff", "Efficiency x Acceptance x pT Resolution", 100, 10, 1200, 800); | |
885 | cEff->Divide(2, 1); | |
886 | cEff->GetPad(1)->SetLogx(kTRUE); | |
887 | cEff->cd(1); | |
888 | TH1* hEffPt = eff->Project(iPt); //the efficiency vs pt | |
889 | hEffPt->SetStats(kFALSE); | |
890 | hEffPt->GetXaxis()->SetRangeUser(0.15, 50.); | |
891 | hEffPt->GetXaxis()->SetMoreLogLabels(kTRUE); | |
892 | hEffPt->GetXaxis()->SetNoExponent(kTRUE); | |
893 | hEffPt->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution"); | |
894 | hEffPt->Draw("E1"); | |
895 | cEff->cd(2); | |
896 | TH1* hEffEta = eff->Project(iEta); //the efficiency vs eta | |
897 | hEffEta->SetStats(kFALSE); | |
898 | hEffEta->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution"); | |
899 | hEffEta->GetYaxis()->SetTitleSize(0.05); | |
900 | hEffEta->Draw("E1"); | |
901 | TH2D* hEffID2Pt = (TH2D*)eff->Project(iPt, iMCid); | |
902 | ||
903 | // Get the efficiencies vs. pT for each species | |
904 | TH1D* hEfficiency[AliPID::kSPECIES]; | |
905 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
906 | hEfficiency[species] = hEffID2Pt->ProjectionX(Form("hEfficiency_%s", AliPID::ParticleShortName(species)), species + 1, species + 1, "e"); | |
907 | hEfficiency[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
908 | hEfficiency[species]->SetLineColor(hYield[species]->GetLineColor()); | |
909 | hEfficiency[species]->SetMarkerColor(hYield[species]->GetLineColor()); | |
910 | hEfficiency[species]->SetLineWidth(2.); | |
911 | hEfficiency[species]->GetXaxis()->SetRangeUser(0.15, 50); | |
912 | hEfficiency[species]->GetYaxis()->SetRangeUser(0., 1.01); | |
913 | hEfficiency[species]->GetXaxis()->SetMoreLogLabels(kTRUE); | |
914 | hEfficiency[species]->GetXaxis()->SetNoExponent(kTRUE); | |
915 | hEfficiency[species]->SetStats(kFALSE); | |
916 | hEfficiency[species]->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution"); | |
917 | hEfficiency[species]->GetYaxis()->SetTitleSize(0.05); | |
918 | } | |
919 | ||
920 | TCanvas* cEff2 = new TCanvas("cEff2", "Efficiency x Acceptance x pT Resolution for different species", 0, 300, 900, 900); | |
921 | cEff2->SetGridx(1); | |
922 | cEff2->SetGridy(1); | |
923 | cEff2->SetLogx(1); | |
924 | ||
925 | hEfficiency[0]->Draw("E1"); | |
926 | ||
927 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) { | |
928 | hEfficiency[i]->Draw("E1 same"); | |
929 | } | |
930 | cEff2->BuildLegend()->SetFillColor(kWhite); | |
931 | ||
932 | ClearTitleFromHistoInCanvas(cEff2); | |
933 | ||
934 | ||
935 | ||
936 | ||
937 | // Efficiency correction for to-pi-ratios | |
938 | TH1D* hEfficiencyToPiRatio[AliPID::kSPECIES] = { 0x0, }; | |
939 | ||
940 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
941 | if (species == AliPID::kPion) | |
942 | continue; // Do not consider pion-to-pion ratio | |
943 | ||
944 | hEfficiencyToPiRatio[species] = new TH1D(*hEfficiency[species]); | |
945 | hEfficiencyToPiRatio[species]->Reset(); | |
946 | hEfficiencyToPiRatio[species]->SetName(Form("hEfficiencyToPionRatio_%s", AliPID::ParticleShortName(species))); | |
947 | hEfficiencyToPiRatio[species]->SetTitle(Form("%s/#pi", AliPID::ParticleLatexName(species))); | |
948 | ||
949 | // Samples for different species are independent, so just divide correction factors | |
950 | hEfficiencyToPiRatio[species]->Divide(hEfficiency[species], hEfficiency[AliPID::kPion], 1., 1., ""); | |
951 | hEfficiencyToPiRatio[species]->GetYaxis()->SetRangeUser(0., 2.0); | |
952 | } | |
953 | ||
954 | TCanvas* cEffToPiRatio = new TCanvas("cEffToPiRatio", "Efficiency x Acceptance x pT Resolution of to-#pi-ratio for different species", | |
955 | 100, 10, 1200, 800); | |
956 | cEffToPiRatio->SetGridx(1); | |
957 | cEffToPiRatio->SetGridy(1); | |
958 | cEffToPiRatio->SetLogx(1); | |
959 | ||
960 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
961 | if (i == AliPID::kPion) | |
962 | continue; | |
963 | ||
964 | hEfficiencyToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same")); | |
965 | } | |
966 | ||
967 | cEffToPiRatio->BuildLegend()->SetFillColor(kWhite); | |
968 | ||
969 | ClearTitleFromHistoInCanvas(cEffToPiRatio); | |
970 | ||
971 | ||
972 | // Correct the yields with the efficiencies and primary fractions | |
973 | TH1D* hYieldCorrected[AliPID::kSPECIES] = { 0x0, }; | |
974 | TH1D* hYieldCorrectedSysError[AliPID::kSPECIES] = { 0x0, }; | |
975 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
976 | hYieldCorrected[species] = new TH1D(*hYield[species]); | |
977 | hYieldCorrected[species]->SetName(Form("%s_corrected", hYield[species]->GetName())); | |
978 | hYieldCorrected[species]->SetTitle(Form("%s", hYield[species]->GetTitle())); | |
979 | //hYieldCorrected[species]->SetTitle(Form("%s, secondary and efficiency x acceptance x pT resolution corrected", | |
980 | // hYield[species]->GetTitle())); | |
981 | ||
982 | // Correction factor histos can have different binning than data histos (constant factor above some threshold) | |
983 | // -> Need special functions to multiply and divide such histos | |
984 | multiplyHistsDifferentBinning(hYieldCorrected[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.); | |
985 | divideHistsDifferentBinning(hYieldCorrected[species], hEfficiency[species], 1., 1.); | |
986 | ||
987 | if (hYieldSysError[species]) { | |
988 | hYieldCorrectedSysError[species] = new TH1D(*hYieldSysError[species]); | |
989 | hYieldCorrectedSysError[species]->SetName(Form("%s_corrected", hYieldSysError[species]->GetName())); | |
990 | hYieldCorrectedSysError[species]->SetTitle(Form("%s", hYieldSysError[species]->GetTitle())); | |
991 | ||
992 | multiplyHistsDifferentBinning(hYieldCorrectedSysError[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1., | |
993 | kTRUE); | |
994 | divideHistsDifferentBinning(hYieldCorrectedSysError[species], hEfficiency[species], 1., 1., kTRUE); | |
995 | } | |
996 | } | |
997 | ||
998 | // Calculate the total corrected yield. The original total yield had no error (just the number of detected tracks in a pT bin), | |
999 | // but due to the correction there is some error for the total yield. Also the error of the fractions introduces uncertainties | |
1000 | // for the yields of individual species | |
1001 | TH1D* hYieldCorrectedTotal = new TH1D(*hYieldCorrected[0]); | |
1002 | hYieldCorrectedTotal->SetLineColor(kBlack); | |
1003 | hYieldCorrectedTotal->SetMarkerColor(kBlack); | |
1004 | hYieldCorrectedTotal->SetName("hYieldCorrectedTotal"); | |
1005 | hYieldCorrectedTotal->SetTitle("Total"); | |
1006 | //hYieldCorrectedTotal->SetTitle("Total yield, secondary and efficiency x acceptance x pT resolution corrected"); | |
1007 | ||
1008 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) | |
1009 | hYieldCorrectedTotal->Add(hYieldCorrected[i], 1.); | |
1010 | ||
1011 | // Calculate the corrected fractions | |
1012 | TH1D* hFractionCorrected[AliPID::kSPECIES] = { 0x0, }; | |
1013 | TH1D* hFractionCorrectedSysError[AliPID::kSPECIES] = { 0x0, }; | |
1014 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
1015 | hFractionCorrected[species] = new TH1D(*hYield[species]); | |
1016 | TString oldName = hYield[species]->GetName(); | |
1017 | TString newName = oldName.ReplaceAll("Yield", "Fraction"); | |
1018 | TString oldTitle = hYield[species]->GetTitle(); | |
1019 | TString newTitle = oldTitle.ReplaceAll("Yield", "Fraction"); | |
1020 | hFractionCorrected[species]->SetName(Form("%s_corrected", newName.Data())); | |
1021 | hFractionCorrected[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
1022 | hFractionCorrected[species]->GetYaxis()->SetTitle("Corrected Fraction"); | |
1023 | hFractionCorrected[species]->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1024 | hFractionCorrected[species]->GetXaxis()->SetNoExponent(kTRUE); | |
1025 | ||
1026 | // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation | |
1027 | // (numerator is a "selected" subset of the denominator). It doesn't matter that the error of the histos is not just "sqrt(content)" | |
1028 | // because the error formula also works for weighted histograms (which means that the error can be more or less anything) | |
1029 | hFractionCorrected[species]->Divide(hYieldCorrected[species], hYieldCorrectedTotal, 1., 1., "B"); | |
1030 | ||
1031 | ||
1032 | // The systematic errors just need to be scaled in the same way as the fractions. | |
1033 | // So, just take the ratio to the uncorrected fraction and scale the sys. error accordingly | |
1034 | // or, in this case, just divide by the same total yield as for yield -> fractions | |
1035 | if (hYieldCorrectedSysError[species]) { | |
1036 | hFractionCorrectedSysError[species] = new TH1D(*hFractionCorrected[species]); | |
1037 | hFractionCorrectedSysError[species]->SetName(Form("%s_sysError", hFractionCorrected[species]->GetName())); | |
1038 | hFractionCorrectedSysError[species]->SetTitle(Form("%s (sys. error)", hFractionCorrected[species]->GetTitle())); | |
1039 | ||
1040 | for (Int_t binX = 1; binX <= hFractionCorrectedSysError[species]->GetNbinsX(); binX++) { | |
1041 | const Double_t corrTotalYield = hYieldCorrectedTotal->GetBinContent(binX); | |
1042 | const Double_t scaleFactor = corrTotalYield > 0 ? 1.0 / corrTotalYield : 1.; | |
1043 | hFractionCorrectedSysError[species]->SetBinError(binX, hYieldCorrectedSysError[species]->GetBinError(binX) * scaleFactor); | |
1044 | } | |
1045 | } | |
1046 | } | |
1047 | ||
1048 | // If MC is available, calculate the generated fractions | |
1049 | TH1D* hMCgenPrimYieldTotal = 0x0; | |
1050 | TH1D* hMCgenPrimFraction[AliPID::kSPECIES]; | |
1051 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) | |
1052 | hMCgenPrimFraction[i] = 0x0; | |
1053 | ||
1054 | if (numMCgenPrimYieldHistsFound > 0) { | |
1055 | hMCgenPrimYieldTotal = new TH1D(*hMCgenPrimYield[0]); | |
1056 | hMCgenPrimYieldTotal->SetLineColor(kBlack); | |
1057 | hMCgenPrimYieldTotal->SetMarkerColor(kBlack); | |
1058 | hMCgenPrimYieldTotal->SetName("hMCgenPrimYieldTotal"); | |
1059 | hMCgenPrimYieldTotal->SetTitle("Total (MC truth)"); | |
1060 | //hMCgenPrimYieldTotal->SetTitle("Total generated primary yield (MC truth)"); | |
1061 | ||
1062 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) | |
1063 | hMCgenPrimYieldTotal->Add(hMCgenPrimYield[i], 1.); | |
1064 | ||
1065 | // Calculate the MC fractions | |
1066 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
1067 | hMCgenPrimYield[species]->SetTitle(Form("%s (MC truth)", AliPID::ParticleLatexName(species))); | |
1068 | ||
1069 | hMCgenPrimFraction[species] = new TH1D(*hMCgenPrimYield[species]); | |
1070 | TString oldName = hMCgenPrimYield[species]->GetName(); | |
1071 | TString newName = oldName.ReplaceAll("Yield", "Fraction"); | |
1072 | TString oldTitle = hMCgenPrimYield[species]->GetTitle(); | |
1073 | TString newTitle = oldTitle.ReplaceAll("yield", "fraction"); | |
1074 | hMCgenPrimFraction[species]->SetName(newName.Data()); | |
1075 | hMCgenPrimFraction[species]->SetTitle(newTitle.Data()); | |
1076 | ||
1077 | // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation | |
1078 | // (numerator is a "selected" subset of the denominator). | |
1079 | hMCgenPrimFraction[species]->Divide(hMCgenPrimFraction[species], hMCgenPrimYieldTotal, 1., 1., "B"); | |
1080 | } | |
1081 | } | |
1082 | ||
1083 | TCanvas* cCorrData = new TCanvas("cCorrData", "Corrected data", 0, 300, 900, 900); | |
1084 | cCorrData->Divide(2, 1);//, 0., 0.01); | |
1085 | cCorrData->GetPad(1)->SetLogx(1); | |
1086 | cCorrData->GetPad(1)->SetLogy(1); | |
1087 | cCorrData->GetPad(2)->SetLogx(1); | |
1088 | cCorrData->GetPad(2)->SetLogy(1); | |
1089 | ||
1090 | cCorrData->GetPad(1)->SetRightMargin(0.001); | |
1091 | cCorrData->GetPad(2)->SetRightMargin(0.001); | |
1092 | ||
1093 | cCorrData->GetPad(1)->SetLeftMargin(0.2); | |
1094 | cCorrData->GetPad(2)->SetLeftMargin(0.2); | |
1095 | ||
1096 | cCorrData->cd(1); // uncorrected | |
1097 | hYield[AliPID::kPion]->GetYaxis()->SetTitleOffset(1.4); | |
1098 | hYield[AliPID::kPion]->Draw("E1"); | |
1099 | ||
1100 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1101 | if (hYieldSysError[i]) | |
1102 | hYieldSysError[i]->Draw("E2 same"); | |
1103 | ||
1104 | if (i == AliPID::kPion) continue; | |
1105 | hYield[i]->Draw("E1 same"); | |
1106 | } | |
1107 | ||
1108 | ClearTitleFromHistoInCanvas(cCorrData, 1); | |
1109 | ||
1110 | cCorrData->cd(2); // corrected | |
1111 | hYieldCorrectedTotal->GetYaxis()->SetTitleOffset(1.4); | |
1112 | hYieldCorrectedTotal->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1113 | hYieldCorrectedTotal->GetXaxis()->SetNoExponent(kTRUE); | |
1114 | hYieldCorrectedTotal->GetYaxis()->SetRangeUser(hYieldCorrected[AliPID::kMuon]->GetBinContent(hYieldCorrected[AliPID::kMuon]->FindLastBinAbove(0.)) * 0.1, | |
1115 | hYieldCorrectedTotal->GetBinContent(hYieldCorrectedTotal->GetMaximumBin()) * 10.); | |
1116 | ||
1117 | if (hMCgenPrimYieldTotal) { | |
1118 | hMCgenPrimYieldTotal->GetYaxis()->SetTitleOffset(1.4); | |
1119 | hMCgenPrimYieldTotal->GetXaxis()->SetMoreLogLabels(kTRUE); | |
1120 | hMCgenPrimYieldTotal->GetXaxis()->SetNoExponent(kTRUE); | |
1121 | hMCgenPrimYieldTotal->GetXaxis()->SetTitle(hYieldCorrectedTotal->GetXaxis()->GetTitle()); | |
1122 | hMCgenPrimYieldTotal->GetYaxis()->SetRangeUser(hYieldCorrected[AliPID::kMuon]->GetBinContent(hYieldCorrected[AliPID::kMuon]->FindLastBinAbove(0.)) * 0.1, | |
1123 | hYieldCorrectedTotal->GetBinContent(hYieldCorrectedTotal->GetMaximumBin()) * 10.); | |
1124 | } | |
1125 | ||
1126 | if (hMCgenPrimYieldTotal) { | |
1127 | hMCgenPrimYieldTotal->Draw("E1"); | |
1128 | hYieldCorrectedTotal->Draw("E1 same"); | |
1129 | } | |
1130 | else | |
1131 | hYieldCorrectedTotal->Draw("E1"); | |
1132 | ||
1133 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1134 | if (hMCgenPrimYield[i]) | |
1135 | hMCgenPrimYield[i]->Draw("E1 same"); | |
1136 | hYieldCorrected[i]->Draw("E1 same"); | |
1137 | } | |
1138 | ||
1139 | TLegend* legTemp = cCorrData->cd(2)->BuildLegend(0.25, 0.16, 0.65, 0.51); | |
1140 | legTemp->SetNColumns(2); | |
1141 | legTemp->SetFillColor(kWhite); | |
1142 | ||
1143 | // Do not include in legend | |
1144 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1145 | if (hYieldCorrectedSysError[i]) | |
1146 | hYieldCorrectedSysError[i]->Draw("E2 same"); | |
1147 | } | |
1148 | ||
1149 | ClearTitleFromHistoInCanvas(cCorrData, 2); | |
1150 | ||
1151 | TCanvas* cCorrYieldsRatio = 0x0; | |
1152 | ||
1153 | TH1D* hYieldCorrectedRatioToMC[AliPID::kSPECIES]; | |
1154 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) | |
1155 | hYieldCorrectedRatioToMC[i] = 0x0; | |
1156 | ||
1157 | TH1D* hYieldCorrectedTotalRatioToMC = 0x0; | |
1158 | ||
1159 | if (numMCgenPrimYieldHistsFound > 0) { | |
1160 | // Compare with MC truth | |
1161 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
1162 | hYieldCorrectedRatioToMC[species] = new TH1D(*hYieldCorrected[species]); | |
1163 | hYieldCorrectedRatioToMC[species]->SetName(Form("hYieldCorrectedRatioToMC_%s", AliPID::ParticleShortName(species))); | |
1164 | hYieldCorrectedRatioToMC[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
1165 | hYieldCorrectedRatioToMC[species]->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth"); | |
1166 | hYieldCorrectedRatioToMC[species]->Divide(hMCgenPrimYield[species]); | |
1167 | } | |
1168 | ||
1169 | hYieldCorrectedTotalRatioToMC = new TH1D(*hYieldCorrectedTotal); | |
1170 | hYieldCorrectedTotalRatioToMC->SetName("hYieldCorrectedTotalRatioToMC"); | |
1171 | hYieldCorrectedTotalRatioToMC->SetTitle("Total"); | |
1172 | hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth"); | |
1173 | hYieldCorrectedTotalRatioToMC->Divide(hMCgenPrimYieldTotal); | |
1174 | ||
1175 | cCorrYieldsRatio = new TCanvas("cCorrYieldsRatio", "Corrected Yields Comparison to MC", 0, 300, 900, 900); | |
1176 | cCorrYieldsRatio->SetGridx(1); | |
1177 | cCorrYieldsRatio->SetGridy(1); | |
1178 | cCorrYieldsRatio->SetLogx(1); | |
1179 | ||
1180 | hYieldCorrectedTotalRatioToMC->GetYaxis()->SetRangeUser(0.6, 1.6); | |
1181 | hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitleOffset(0.85); | |
1182 | hYieldCorrectedTotalRatioToMC->Draw("E1"); | |
1183 | ||
1184 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) | |
1185 | hYieldCorrectedRatioToMC[species]->Draw("E1 same"); | |
1186 | ||
1187 | cCorrYieldsRatio->BuildLegend()->SetFillColor(kWhite); | |
1188 | ||
1189 | ClearTitleFromHistoInCanvas(cCorrYieldsRatio); | |
1190 | } | |
1191 | ||
1192 | ||
1193 | if (cFractions) | |
1194 | cFractions->Draw(); | |
1195 | ||
1196 | TCanvas* cCorrFractions = new TCanvas("cCorrFractions", "Corrected particleFractions", 0, 300, 900, 900); | |
1197 | cCorrFractions->SetLogx(1); | |
1198 | hFractionCorrected[0]->GetYaxis()->SetRangeUser(0., 1.); | |
1199 | hFractionCorrected[0]->Draw("E1"); | |
1200 | if (hMCgenPrimFraction[0]) | |
1201 | hMCgenPrimFraction[0]->Draw("E1 same"); | |
1202 | ||
1203 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) { | |
1204 | hFractionCorrected[i]->Draw("E1 same"); | |
1205 | if (hMCgenPrimFraction[i]) | |
1206 | hMCgenPrimFraction[i]->Draw("E1 same"); | |
1207 | } | |
1208 | ||
1209 | cCorrFractions->BuildLegend()->SetFillColor(kWhite); | |
1210 | ||
1211 | // Do not include in legend!! | |
1212 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1213 | if (hFractionCorrectedSysError[i]) | |
1214 | hFractionCorrectedSysError[i]->Draw("E2 same"); | |
1215 | } | |
1216 | ||
1217 | ||
1218 | ClearTitleFromHistoInCanvas(cCorrFractions); | |
1219 | ||
1220 | // Save results to file | |
1221 | TString chargeString = ""; | |
1222 | if (chargeMode == kPosCharge) | |
1223 | chargeString = "_posCharge"; | |
1224 | else if (chargeMode == kNegCharge) | |
1225 | chargeString = "_negCharge"; | |
1226 | ||
1227 | TString saveFileName = pathNameData; | |
1228 | saveFileName.ReplaceAll(Form("%s/", pathData.Data()), ""); | |
1229 | saveFileName.Prepend("output_EfficiencyCorrection_"); | |
1230 | saveFileName.ReplaceAll(".root", Form("%s.root", chargeString.Data())); | |
1231 | ||
1232 | TString saveFilePathName = Form("%s/%s", pathData.Data(), saveFileName.Data()); | |
1233 | TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE"); | |
1234 | saveFile->cd(); | |
1235 | ||
1236 | if (cSec) | |
1237 | cSec->Write(); | |
1238 | ||
1239 | if (cSecSS) | |
1240 | cSecSS->Write(); | |
1241 | ||
1242 | if (cSec2) | |
1243 | cSec2->Write(); | |
1244 | ||
1245 | if (cSecSS2) | |
1246 | cSecSS2->Write(); | |
1247 | ||
1248 | if (cSecToPiRatio) | |
1249 | cSecToPiRatio->Write(); | |
1250 | ||
1251 | if (cSecSSToPiRatio) | |
1252 | cSecSSToPiRatio->Write(); | |
1253 | ||
1254 | if (cEff) | |
1255 | cEff->Write(); | |
1256 | ||
1257 | if (cEff2) | |
1258 | cEff2->Write(); | |
1259 | ||
1260 | if (cEffToPiRatio) | |
1261 | cEffToPiRatio->Write(); | |
1262 | ||
1263 | if (cCorrData) | |
1264 | cCorrData->Write(); | |
1265 | ||
1266 | if (cCorrYieldsRatio) | |
1267 | cCorrYieldsRatio->Write(); | |
1268 | ||
1269 | if (cFractions) | |
1270 | cFractions->Write(); | |
1271 | ||
1272 | if (cCorrFractions) | |
1273 | cCorrFractions->Write(); | |
1274 | ||
1275 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1276 | if (hSec[i]) | |
1277 | hSec[i]->Write(); | |
1278 | ||
1279 | if (hSecSS[i]) | |
1280 | hSecSS[i]->Write(); | |
1281 | } | |
1282 | ||
1283 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1284 | if (hSecToPiRatio[i]) | |
1285 | hSecToPiRatio[i]->Write(); | |
1286 | ||
1287 | if (hSecSSToPiRatio[i]) | |
1288 | hSecSSToPiRatio[i]->Write(); | |
1289 | } | |
1290 | ||
1291 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1292 | if (hEfficiency[i]) | |
1293 | hEfficiency[i]->Write(); | |
1294 | } | |
1295 | ||
1296 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1297 | if (hEfficiencyToPiRatio[i]) | |
1298 | hEfficiencyToPiRatio[i]->Write(); | |
1299 | } | |
1300 | ||
1301 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1302 | if (hYield[i]) | |
1303 | hYield[i]->Write(); | |
1304 | ||
1305 | if (hYieldSysError[i]) | |
1306 | hYieldSysError[i]->Write(); | |
1307 | } | |
1308 | ||
1309 | if (hYieldCorrectedTotal) | |
1310 | hYieldCorrectedTotal->Write(); | |
1311 | ||
1312 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1313 | if (hYieldCorrected[i]) | |
1314 | hYieldCorrected[i]->Write(); | |
1315 | ||
1316 | if (hYieldCorrectedSysError[i]) | |
1317 | hYieldCorrectedSysError[i]->Write(); | |
1318 | } | |
1319 | ||
1320 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1321 | if (hMCgenPrimYield[i]) | |
1322 | hMCgenPrimYield[i]->Write(); | |
1323 | } | |
1324 | ||
1325 | if (hMCgenPrimYieldTotal) | |
1326 | hMCgenPrimYieldTotal->Write(); | |
1327 | ||
1328 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1329 | if (hYieldCorrectedRatioToMC[i]) | |
1330 | hYieldCorrectedRatioToMC[i]->Write(); | |
1331 | } | |
1332 | ||
1333 | if (hYieldCorrectedTotalRatioToMC) | |
1334 | hYieldCorrectedTotalRatioToMC->Write(); | |
1335 | ||
1336 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1337 | if (hFractionCorrected[i]) | |
1338 | hFractionCorrected[i]->Write(); | |
1339 | ||
1340 | if (hFractionCorrectedSysError[i]) | |
1341 | hFractionCorrectedSysError[i]->Write(); | |
1342 | } | |
1343 | ||
1344 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
1345 | if (hMCgenPrimFraction[i]) | |
1346 | hMCgenPrimFraction[i]->Write(); | |
1347 | } | |
1348 | ||
1349 | TNamed* settings = new TNamed( | |
1350 | Form("Settings: Efficiency file \"%s\", Data file \"%s\", lowerCentrality %.3f, upperCentrality %.3f, lowerJetPt %.1f, upperJetPt %.1f, constantCorrectionAbovePtThreshold %.3f\n", | |
1351 | pathNameEfficiency.Data(), pathNameData.Data(), lowerCentrality, upperCentrality, lowerJetPt, | |
1352 | upperJetPt, constantCorrectionAbovePtThreshold), ""); | |
1353 | settings->Write(); | |
1354 | ||
1355 | saveFile->Close(); | |
1356 | ||
1357 | return 0; | |
1358 | } | |
1359 | ||
1360 | ||
1361 | //___________________________________________________________________ | |
1362 | // Efficiency for jet spectra (pT, z and xi) | |
1363 | // E.g. a 'calcEfficiency.C+("finalCuts/MC_pp/7TeV/LHC11a1/corrected/allJetPtMergedWeighted/bhess_PID_Jets_PureGauss_efficiency.root", "finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/output_extractedFFs_bhess_PID_Jets_centrality_all_highN_rougherXi.root", kTRUE, kFALSE, 0, -2, -2, -2, -2, 40, 80, -100, -100, 2, 1, 2)' -b -q | |
1364 | Int_t calcEfficiency(TString pathNameEfficiency, TString pathNameData, Bool_t correctGeantFluka, Bool_t scaleStrangeness, | |
1365 | Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/, | |
1366 | Double_t lowerCentralityData /*= -2*/, Double_t upperCentralityData /*= -2*/, | |
1367 | Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/, | |
1368 | Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/, | |
1369 | Double_t constantCorrectionAbovePtThreshold, | |
1370 | Double_t constantCorrectionAboveXiThreshold, | |
1371 | Int_t rebinEfficiencyPt, | |
1372 | Int_t rebinEfficiencyZ, | |
1373 | Int_t rebinEfficiencyXi) | |
1374 | { | |
1375 | TString titles[kNtypes]; | |
1376 | titles[kTrackPt] = "trackPt"; | |
1377 | titles[kZ] = "z"; | |
1378 | titles[kXi] = "xi"; | |
1379 | ||
1380 | TString pathData = pathNameData; | |
1381 | pathData.Replace(pathData.Last('/'), pathData.Length(), ""); | |
1382 | ||
1383 | TFile* fileEff = TFile::Open(pathNameEfficiency.Data()); | |
1384 | if (!fileEff) { | |
1385 | printf("Failed to open efficiency file \"%s\"\n", pathNameEfficiency.Data()); | |
1386 | return -1; | |
1387 | } | |
1388 | ||
1389 | AliCFContainer* data = (AliCFContainer*)(fileEff->Get("containerEff")); | |
1390 | if (!data) { | |
1391 | printf("Failed to load efficiency container!\n"); | |
1392 | return -1; | |
1393 | } | |
1394 | ||
1395 | const Int_t iMCid = data->GetVar("MC ID"); | |
1396 | const Int_t iPt = data->GetVar("P_{T} (GeV/c)"); | |
1397 | //const Int_t iEta = data->GetVar("#eta"); | |
1398 | const Int_t iCharge = data->GetVar("Charge (e_{0})"); | |
1399 | const Int_t iMult = data->GetVar("Centrality Percentile"); | |
1400 | const Int_t iJetPt = data->GetVar("P_{T}^{jet} (GeV/c)"); | |
1401 | const Int_t iZ = data->GetVar("z = P_{T}^{track} / P_{T}^{jet}"); | |
1402 | const Int_t iXi = data->GetVar("#xi = ln(P_{T}^{jet} / P_{T}^{track})"); | |
1403 | ||
1404 | TH2D* hNjetsGen = 0x0; | |
1405 | TH2D* hNjetsRec = 0x0; | |
1406 | ||
1407 | TH2D* hNjetsGenData = 0x0; | |
1408 | TH2D* hNjetsRecData = 0x0; | |
1409 | ||
1410 | TH1D* hYieldPt[AliPID::kSPECIES] = {0x0, }; | |
1411 | TH1D* hYieldPtSysError[AliPID::kSPECIES] = {0x0, }; | |
1412 | TH1D* hMCgenPrimYieldPt[AliPID::kSPECIES] = {0x0, }; | |
1413 | ||
1414 | TH1D* hYieldZ[AliPID::kSPECIES] = {0x0, }; | |
1415 | TH1D* hYieldZSysError[AliPID::kSPECIES] = {0x0, }; | |
1416 | TH1D* hMCgenPrimYieldZ[AliPID::kSPECIES] = {0x0, }; | |
1417 | ||
1418 | TH1D* hYieldXi[AliPID::kSPECIES] = {0x0, }; | |
1419 | TH1D* hYieldXiSysError[AliPID::kSPECIES] = {0x0, }; | |
1420 | TH1D* hMCgenPrimYieldXi[AliPID::kSPECIES] = {0x0, }; | |
1421 | ||
1422 | ||
1423 | TH2D* hTemp = 0x0; | |
1424 | ||
1425 | Int_t numMCgenPrimYieldHistsFound = 0; | |
1426 | ||
1427 | Int_t lowerJetPtBinYields = -1; | |
1428 | Int_t upperJetPtBinYields = -1; | |
1429 | ||
1430 | Double_t actualLowerJetPtYields = -10; | |
1431 | Double_t actualUpperJetPtYields = -10; | |
1432 | ||
1433 | TFile* fileData = TFile::Open(pathNameData.Data()); | |
1434 | if (!fileData) { | |
1435 | printf("Failed to open data file \"%s\"\n", pathNameData.Data()); | |
1436 | return -1; | |
1437 | } | |
1438 | ||
1439 | ||
1440 | hNjetsGenData = (TH2D*)fileData->Get("fh2FFJetPtGen"); | |
1441 | ||
1442 | hNjetsRecData = (TH2D*)fileData->Get("fh2FFJetPtRec"); | |
1443 | ||
1444 | if (!hNjetsRecData) { | |
1445 | printf("Failed to load numJet histo for data!\n"); | |
1446 | return -1; | |
1447 | } | |
1448 | ||
1449 | const Bool_t restrictCentralityData = ((lowerCentralityData >= -1) && (upperCentralityData >= -1)); | |
1450 | // Integral(lowerCentBinLimitData, uppCentBinLimitData) will not be restricted if these values are kept | |
1451 | const Int_t lowerCentralityBinLimitData = restrictCentralityData ? hNjetsRecData->GetXaxis()->FindBin(lowerCentralityData + 0.001) | |
1452 | : -1; | |
1453 | const Int_t upperCentralityBinLimitData = restrictCentralityData ? hNjetsRecData->GetXaxis()->FindBin(upperCentralityData - 0.001) | |
1454 | : -2; | |
1455 | ||
1456 | // Integral(lowerJetBinLimitData, uppJetBinLimitData) will not be restricted if these values are kept | |
1457 | const Int_t lowerJetPtBinLimitData = ((lowerJetPt >= 0) && (upperJetPt >= 0)) | |
1458 | ? hNjetsRecData->GetYaxis()->FindBin(lowerJetPt + 0.001) : -1; | |
1459 | const Int_t upperJetPtBinLimitData = ((lowerJetPt >= 0) && (upperJetPt >= 0)) | |
1460 | ? hNjetsRecData->GetYaxis()->FindBin(upperJetPt - 0.001) : -2; | |
1461 | ||
1462 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
1463 | TString histName = Form("hIDFFz_%s", AliPID::ParticleShortName(species)); | |
1464 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1465 | ||
1466 | if (hTemp && (lowerJetPtBinYields < 0 || upperJetPtBinYields < 0)) { | |
1467 | // If the jetPt range is to be restricted | |
1468 | if (lowerJetPt >= 0 && upperJetPt >= 0) { | |
1469 | lowerJetPtBinYields = hTemp->GetYaxis()->FindBin(lowerJetPt + 0.001); | |
1470 | upperJetPtBinYields = hTemp->GetYaxis()->FindBin(upperJetPt - 0.001); | |
1471 | ||
1472 | actualLowerJetPtYields = hTemp->GetYaxis()->GetBinLowEdge(lowerJetPtBinYields); | |
1473 | actualUpperJetPtYields = hTemp->GetYaxis()->GetBinUpEdge(upperJetPtBinYields); | |
1474 | } | |
1475 | else { | |
1476 | lowerJetPtBinYields = 0; | |
1477 | upperJetPtBinYields = -1; | |
1478 | } | |
1479 | } | |
1480 | ||
1481 | if ((lowerJetPtBinYields < 0 || upperJetPtBinYields < 0) && (lowerJetPt >= 0 && upperJetPt >= 0)) { | |
1482 | printf("Failed to determine jet pt bin limits!\n"); | |
1483 | return -1; | |
1484 | } | |
1485 | ||
1486 | // NOTE: The yields have already been normalised to nJets in extractFFs. However, the binning in jetPt might be different now. | |
1487 | // Thus, in order to get the correct result, the old normalisation (just the same binnnig as in the nJet histo) has to be | |
1488 | // undone and then it has to be normalised to the number of jets in the new binning again (see below). | |
1489 | ||
1490 | if (hTemp) { | |
1491 | undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1492 | hYieldZ[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields, upperJetPtBinYields, | |
1493 | "e"); | |
1494 | } | |
1495 | ||
1496 | if (!hYieldZ[species]) { | |
1497 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1498 | return -1; | |
1499 | } | |
1500 | ||
1501 | ||
1502 | ||
1503 | histName = Form("hIDFFxi_%s", AliPID::ParticleShortName(species)); | |
1504 | ||
1505 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1506 | ||
1507 | if (hTemp) { | |
1508 | undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1509 | hYieldXi[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), | |
1510 | lowerJetPtBinYields, upperJetPtBinYields, "e"); | |
1511 | } | |
1512 | ||
1513 | if (!hYieldXi[species]) { | |
1514 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1515 | return -1; | |
1516 | } | |
1517 | ||
1518 | histName = Form("hIDFFtrackPt_%s", AliPID::ParticleShortName(species)); | |
1519 | ||
1520 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1521 | ||
1522 | if (hTemp) { | |
1523 | undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1524 | hYieldPt[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), | |
1525 | lowerJetPtBinYields, upperJetPtBinYields, "e"); | |
1526 | } | |
1527 | ||
1528 | if (!hYieldPt[species]) { | |
1529 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1530 | return -1; | |
1531 | } | |
1532 | ||
1533 | ||
1534 | ||
1535 | histName = Form("hIDFFz_%s_sysError", AliPID::ParticleShortName(species)); | |
1536 | ||
1537 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1538 | ||
1539 | if (hTemp) { | |
1540 | undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1541 | hYieldZSysError[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields, | |
1542 | upperJetPtBinYields, "e"); | |
1543 | } | |
1544 | ||
1545 | if (!hYieldZSysError[species]) { | |
1546 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1547 | return -1; | |
1548 | } | |
1549 | ||
1550 | histName = Form("hIDFFxi_%s_sysError", AliPID::ParticleShortName(species)); | |
1551 | ||
1552 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1553 | ||
1554 | if (hTemp) { | |
1555 | undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1556 | hYieldXiSysError[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields, | |
1557 | upperJetPtBinYields, "e"); | |
1558 | } | |
1559 | ||
1560 | if (!hYieldXiSysError[species]) { | |
1561 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1562 | return -1; | |
1563 | } | |
1564 | ||
1565 | histName = Form("hIDFFtrackPt_%s_sysError", AliPID::ParticleShortName(species)); | |
1566 | ||
1567 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1568 | ||
1569 | if (hTemp) { | |
1570 | undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1571 | hYieldPtSysError[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields, | |
1572 | upperJetPtBinYields, "e"); | |
1573 | } | |
1574 | ||
1575 | if (!hYieldPtSysError[species]) { | |
1576 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1577 | return -1; | |
1578 | } | |
1579 | ||
1580 | // In case of MC also retrieve the MC truth generated yields | |
1581 | histName = Form("hMCgenYieldsPrimZ_%s", AliPID::ParticleShortName(species)); | |
1582 | ||
1583 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1584 | ||
1585 | if (hTemp) { | |
1586 | undoJetNormalisationHist2D(hTemp, hNjetsGenData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1587 | hMCgenPrimYieldZ[species] = hTemp->ProjectionX(histName.ReplaceAll("hMCgen", "hYieldMCgen").Data(), lowerJetPtBinYields, | |
1588 | upperJetPtBinYields, "e"); | |
1589 | } | |
1590 | ||
1591 | ||
1592 | histName = Form("hMCgenYieldsPrimXi_%s", AliPID::ParticleShortName(species)); | |
1593 | ||
1594 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1595 | ||
1596 | if (hTemp) { | |
1597 | undoJetNormalisationHist2D(hTemp, hNjetsGenData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1598 | hMCgenPrimYieldXi[species] = hTemp->ProjectionX(histName.ReplaceAll("hMCgen", "hYieldMCgen").Data(), lowerJetPtBinYields, | |
1599 | upperJetPtBinYields, "e"); | |
1600 | } | |
1601 | ||
1602 | ||
1603 | histName = Form("hMCgenYieldsPrimPt_%s", AliPID::ParticleShortName(species)); | |
1604 | ||
1605 | hTemp = (TH2D*)fileData->Get(histName.Data()); | |
1606 | ||
1607 | if (hTemp) { | |
1608 | undoJetNormalisationHist2D(hTemp, hNjetsGenData, lowerCentralityBinLimitData, upperCentralityBinLimitData); | |
1609 | hMCgenPrimYieldPt[species] = hTemp->ProjectionX(histName.ReplaceAll("hMCgen", "hYieldMCgen").Data(), lowerJetPtBinYields, | |
1610 | upperJetPtBinYields, "e"); | |
1611 | } | |
1612 | ||
1613 | if (hMCgenPrimYieldPt[species] && hMCgenPrimYieldZ[species] && hMCgenPrimYieldXi[species]) | |
1614 | numMCgenPrimYieldHistsFound++; | |
1615 | } | |
1616 | ||
1617 | if (numMCgenPrimYieldHistsFound > 0 && numMCgenPrimYieldHistsFound != AliPID::kSPECIES) { | |
1618 | printf("Error: Unable to retrieve all MC generated prim yield histos! Got %d.\n", numMCgenPrimYieldHistsFound); | |
1619 | return -1; | |
1620 | } | |
1621 | ||
1622 | ||
1623 | /*OLD: Load from AnalysisResults.root | |
1624 | ||
1625 | File* fileData = TFile::Open(pathNameData.Data()); | |
1626 | if (!fileData) { | |
1627 | printf("Failed to open data file \"%s\"\n", pathNameData.Data()); | |
1628 | return -1; | |
1629 | } | |
1630 | ||
1631 | TString dirDataInFile = ""; | |
1632 | TDirectory* dirData = (TDirectory*)fileData->Get(fileData->GetListOfKeys()->At(0)->GetName()); | |
1633 | if (!dirData) { | |
1634 | printf("No directory found inside data file \"%s\"\n", pathNameData.Data()); | |
1635 | return -1; | |
1636 | } | |
1637 | ||
1638 | TList* list = (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()); | |
1639 | if (!list) { | |
1640 | printf("No list found in directory \"%s\" inside data file \"%s\"\n", | |
1641 | dirData->GetListOfKeys()->At(0)->GetName(),pathNameData.Data()); | |
1642 | return -1; | |
1643 | } | |
1644 | ||
1645 | hNjetsGen = (TH1D*)list->FindObject("fh1FFJetPtGen"); | |
1646 | hNjetsRec = (TH1D*)list->FindObject("fh1FFJetPtRecCuts"); | |
1647 | ||
1648 | if (!hNjetsRec) { | |
1649 | printf("Failed to load number of jets (rec) histo!\n"); | |
1650 | return -1; | |
1651 | } | |
1652 | ||
1653 | ||
1654 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
1655 | TString histName = Form("fh2FFZRecCuts_%s", AliPID::ParticleName(species)); | |
1656 | hTemp = (TH2D*)list->FindObject(histName.Data()); | |
1657 | ||
1658 | if (hTemp && (lowerJetPtBinYields < 0 || upperJetPtBinYields < 0)) { | |
1659 | // If the jetPt range is to be restricted | |
1660 | if (lowerJetPt >= 0 && upperJetPt >= 0) { | |
1661 | lowerJetPtBinYields = hTemp->GetXaxis()->FindBin(lowerJetPt + 0.001); | |
1662 | upperJetPtBinYields = hTemp->GetXaxis()->FindBin(upperJetPt - 0.001); | |
1663 | ||
1664 | actualLowerJetPtYields = hTemp->GetXaxis()->GetBinLowEdge(lowerJetPtBinYields); | |
1665 | actualUpperJetPtYields = hTemp->GetXaxis()->GetBinUpEdge(upperJetPtBinYields); | |
1666 | } | |
1667 | else { | |
1668 | lowerJetPtBinYields = 0; | |
1669 | upperJetPtBinYields = -1; | |
1670 | } | |
1671 | } | |
1672 | ||
1673 | if ((lowerJetPtBinYields < 0 || upperJetPtBinYields < 0) && (lowerJetPt >= 0 && upperJetPt >= 0)) { | |
1674 | printf("Failed to determine jet pt bin limits!\n"); | |
1675 | return -1; | |
1676 | } | |
1677 | ||
1678 | if (hTemp) | |
1679 | hYieldZ[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, "e"); | |
1680 | ||
1681 | if (!hYieldZ[species]) { | |
1682 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1683 | return -1; | |
1684 | } | |
1685 | ||
1686 | ||
1687 | ||
1688 | histName = Form("fh2FFXiRecCuts_%s", AliPID::ParticleName(species)); | |
1689 | ||
1690 | hTemp = (TH2D*)list->FindObject(histName.Data()); | |
1691 | ||
1692 | if (hTemp) | |
1693 | hYieldXi[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, "e"); | |
1694 | ||
1695 | if (!hYieldXi[species]) { | |
1696 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1697 | return -1; | |
1698 | } | |
1699 | ||
1700 | histName = Form("fh2FFTrackPtRecCuts_%s", AliPID::ParticleName(species)); | |
1701 | ||
1702 | hTemp = (TH2D*)list->FindObject(histName.Data()); | |
1703 | ||
1704 | if (hTemp) | |
1705 | hYieldPt[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, "e"); | |
1706 | ||
1707 | if (!hYieldPt[species]) { | |
1708 | printf("Failed to load hist \"%s\"\n", histName.Data()); | |
1709 | return -1; | |
1710 | } | |
1711 | ||
1712 | // In case of MC also retrieve the MC truth generated yields | |
1713 | histName = Form("fh2FFZGen_%s", AliPID::ParticleName(species)); | |
1714 | ||
1715 | hTemp = (TH2D*)list->FindObject(histName.Data()); | |
1716 | if (hTemp) | |
1717 | hMCgenPrimYieldZ[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, | |
1718 | "e"); | |
1719 | ||
1720 | ||
1721 | histName = Form("fh2FFXiGen_%s", AliPID::ParticleName(species)); | |
1722 | ||
1723 | hTemp = (TH2D*)list->FindObject(histName.Data()); | |
1724 | if (hTemp) | |
1725 | hMCgenPrimYieldXi[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, | |
1726 | "e"); | |
1727 | ||
1728 | ||
1729 | histName = Form("fh2FFTrackPtGen_%s", AliPID::ParticleName(species)); | |
1730 | ||
1731 | hTemp = (TH2D*)list->FindObject(histName.Data()); | |
1732 | if (hTemp) | |
1733 | hMCgenPrimYieldPt[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, | |
1734 | "e"); | |
1735 | ||
1736 | if (hMCgenPrimYieldPt[species] && hMCgenPrimYieldZ[species] && hMCgenPrimYieldXi[species]) | |
1737 | numMCgenPrimYieldHistsFound++; | |
1738 | } | |
1739 | ||
1740 | if (numMCgenPrimYieldHistsFound > 0 && numMCgenPrimYieldHistsFound != AliPID::kSPECIES) { | |
1741 | printf("Error: Unable to retrieve all MC generated prim yield histos! Got %d.\n", numMCgenPrimYieldHistsFound); | |
1742 | return -1; | |
1743 | } | |
1744 | ||
1745 | ||
1746 | const Double_t nJetsGen = hNjetsGen ? hNjetsGen->Integral(lowerJetPtBinYields, upperJetPtBinYields) : -1; | |
1747 | const Double_t nJetsRec = hNjetsRec->Integral(lowerJetPtBinYields, upperJetPtBinYields); | |
1748 | */ | |
1749 | ||
1750 | // Take binning from pion yield (binning for all species the same) and create a new AliCFContainer with this new binning | |
1751 | TH1D hDummyPt(*hYieldPt[AliPID::kPion]); | |
1752 | hDummyPt.SetName("hDummyPt"); | |
1753 | TAxis* axisPt = 0x0; | |
1754 | ||
1755 | if (rebinEfficiencyPt > 1) { | |
1756 | Int_t nBinsNew = 0; | |
1757 | const Double_t* newBins = GetBins(rebinEfficiencyPt, nBinsNew); | |
1758 | ||
1759 | hDummyPt.SetBins(nBinsNew, newBins); | |
1760 | ||
1761 | axisPt = hDummyPt.GetXaxis(); | |
1762 | ||
1763 | //axisPt = hDummyPt.Rebin(rebinEfficiencyPt, "", 0)->GetXaxis(); | |
1764 | } | |
1765 | else | |
1766 | axisPt = hDummyPt.GetXaxis(); | |
1767 | ||
1768 | const TArrayD* binsPtCurrent = axisPt->GetXbins(); | |
1769 | TArrayD* binsPtNew = new TArrayD(*binsPtCurrent); | |
1770 | ||
1771 | TH1D hDummyZ(*hYieldZ[AliPID::kPion]); | |
1772 | hDummyZ.SetName("hDummyZ"); | |
1773 | TAxis* axisZ = 0x0; | |
1774 | ||
1775 | if (rebinEfficiencyZ > 1) | |
1776 | axisZ = hDummyZ.Rebin(rebinEfficiencyZ, "", 0)->GetXaxis(); | |
1777 | else | |
1778 | axisZ = hDummyZ.GetXaxis(); | |
1779 | ||
1780 | const TArrayD* binsZCurrent = axisZ->GetXbins(); | |
1781 | TArrayD* binsZNew = new TArrayD(*binsZCurrent); | |
1782 | ||
1783 | TH1D hDummyXi(*hYieldXi[AliPID::kPion]); | |
1784 | hDummyXi.SetName("hDummyXi"); | |
1785 | TAxis* axisXi = 0x0; | |
1786 | ||
1787 | if (rebinEfficiencyXi > 1) | |
1788 | axisXi = hDummyXi.Rebin(rebinEfficiencyXi, "", 0)->GetXaxis(); | |
1789 | else | |
1790 | axisXi = hDummyXi.GetXaxis(); | |
1791 | ||
1792 | const TArrayD* binsXiCurrent = axisXi->GetXbins(); | |
1793 | TArrayD* binsXiNew = new TArrayD(*binsXiCurrent); | |
1794 | ||
1795 | ||
1796 | const Int_t nEffDims = data->GetNVar(); | |
1797 | Int_t nEffBins[nEffDims]; | |
1798 | ||
1799 | for (Int_t iDim = 0; iDim < nEffDims; iDim++) { | |
1800 | if (iDim == iPt) | |
1801 | nEffBins[iDim] = axisPt->GetNbins(); | |
1802 | else if (iDim == iZ) | |
1803 | nEffBins[iDim] = axisZ->GetNbins(); | |
1804 | else if (iDim == iXi) | |
1805 | nEffBins[iDim] = axisXi->GetNbins(); | |
1806 | else | |
1807 | nEffBins[iDim] = data->GetNBins(iDim); | |
1808 | } | |
1809 | ||
1810 | ||
1811 | // Just make one large pT bin above some threshold, if desired | |
1812 | if (binsPtNew->fN != 0 && constantCorrectionAbovePtThreshold > 0) { | |
1813 | for (Int_t iBin = 0; iBin <= nEffBins[iPt] + 1; iBin++) { | |
1814 | // Find the first bin edged really larger than the threshold. | |
1815 | // If the bin edge before equals the threshold, just set the | |
1816 | // current bin edge to the right end of the spectrum -> Done. | |
1817 | // If the bin edge before is different, set the bin edge to the | |
1818 | // threshold | |
1819 | if (binsPtNew->fArray[iBin] > constantCorrectionAbovePtThreshold) { | |
1820 | if (binsPtNew->fArray[iBin - 1] == constantCorrectionAbovePtThreshold) { | |
1821 | binsPtNew->fArray[iBin] = binsPtNew->fArray[nEffBins[iPt]]; | |
1822 | nEffBins[iPt] = iBin; | |
1823 | break; | |
1824 | } | |
1825 | else { | |
1826 | binsPtNew->fArray[iBin] = constantCorrectionAbovePtThreshold; | |
1827 | } | |
1828 | } | |
1829 | } | |
1830 | } | |
1831 | ||
1832 | // Just make one large Xi bin above some threshold, if desired | |
1833 | if (binsXiNew->fN != 0 && constantCorrectionAboveXiThreshold > 0) { | |
1834 | for (Int_t iBin = 0; iBin <= nEffBins[iXi] + 1; iBin++) { | |
1835 | // Find the first bin edged really larger than the threshold. | |
1836 | // If the bin edge before equals the threshold, just set the | |
1837 | // current bin edge to the right end of the spectrum -> Done. | |
1838 | // If the bin edge before is different, set the bin edge to the | |
1839 | // threshold | |
1840 | if (binsXiNew->fArray[iBin] > constantCorrectionAboveXiThreshold) { | |
1841 | if (binsXiNew->fArray[iBin - 1] == constantCorrectionAboveXiThreshold) { | |
1842 | binsXiNew->fArray[iBin] = binsXiNew->fArray[nEffBins[iXi]]; | |
1843 | nEffBins[iXi] = iBin; | |
1844 | break; | |
1845 | } | |
1846 | else { | |
1847 | binsXiNew->fArray[iBin] = constantCorrectionAboveXiThreshold; | |
1848 | } | |
1849 | } | |
1850 | } | |
1851 | } | |
1852 | ||
1853 | AliCFContainer *dataRebinned = new AliCFContainer(Form("%s_rebinned", data->GetName()), Form("%s (rebinned)", data->GetTitle()), | |
1854 | data->GetNStep(), nEffDims, nEffBins); | |
1855 | ||
1856 | for (Int_t iDim = 0; iDim < nEffDims; iDim++) { | |
1857 | dataRebinned->SetVarTitle(iDim, data->GetVarTitle(iDim)); | |
1858 | ||
1859 | if (iDim == iPt) { | |
1860 | if (binsPtNew->fN == 0) | |
1861 | dataRebinned->SetBinLimits(iDim, axisPt->GetXmin(), axisPt->GetXmax()); | |
1862 | else | |
1863 | dataRebinned->SetBinLimits(iDim, binsPtNew->fArray); | |
1864 | } | |
1865 | else if (iDim == iZ) { | |
1866 | if (binsZNew->fN == 0) | |
1867 | dataRebinned->SetBinLimits(iDim, axisZ->GetXmin(), axisZ->GetXmax()); | |
1868 | else | |
1869 | dataRebinned->SetBinLimits(iDim, binsZNew->fArray); | |
1870 | } | |
1871 | else if (iDim == iXi) { | |
1872 | if (binsXiNew->fN == 0) | |
1873 | dataRebinned->SetBinLimits(iDim, axisXi->GetXmin(), axisXi->GetXmax()); | |
1874 | else | |
1875 | dataRebinned->SetBinLimits(iDim, binsXiNew->fArray); | |
1876 | } | |
1877 | else { | |
1878 | dataRebinned->SetBinLimits(iDim, data->GetBinLimits(iDim)); | |
1879 | } | |
1880 | } | |
1881 | ||
1882 | for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++) | |
1883 | dataRebinned->SetStepTitle(iStep, data->GetStepTitle(iStep)); | |
1884 | ||
1885 | Int_t coord[nEffDims]; | |
1886 | Double_t binCenterCoord[nEffDims]; | |
1887 | ||
1888 | // Fill content from old grid into the new grid with proper binning | |
1889 | for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++) { | |
1890 | Long64_t nBinsGrid = data->GetGrid(iStep)->GetGrid()->GetNbins(); | |
1891 | ||
1892 | for (Long64_t iBin = 0; iBin <= nBinsGrid + 1; iBin++) { | |
1893 | Double_t binContent = data->GetGrid(iStep)->GetGrid()->GetBinContent(iBin, coord); | |
1894 | Double_t binError2 = data->GetGrid(iStep)->GetGrid()->GetBinError2(iBin); | |
1895 | ||
1896 | for (Int_t iDim = 0; iDim < nEffDims; iDim++) { | |
1897 | binCenterCoord[iDim] = data->GetBinCenter(iDim, coord[iDim]); | |
1898 | } | |
1899 | ||
1900 | Long64_t iBinRebinned = dataRebinned->GetGrid(iStep)->GetGrid()->GetBin(binCenterCoord); | |
1901 | dataRebinned->GetGrid(iStep)->GetGrid()->AddBinContent(iBinRebinned, binContent); | |
1902 | dataRebinned->GetGrid(iStep)->GetGrid()->AddBinError2(iBinRebinned, binError2); | |
1903 | } | |
1904 | } | |
1905 | ||
1906 | ||
1907 | // If desired, restrict centrality axis | |
1908 | Int_t lowerCentralityBinLimit = -1; | |
1909 | Int_t upperCentralityBinLimit = -2; // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept | |
1910 | Bool_t restrictCentralityAxis = kFALSE; | |
1911 | Double_t actualLowerCentrality = -1.; | |
1912 | Double_t actualUpperCentrality = -1.; | |
1913 | ||
1914 | if (lowerCentrality >= -1 && upperCentrality >= -1) { | |
1915 | // Add subtract a very small number to avoid problems with values right on the border between to bins | |
1916 | lowerCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(lowerCentrality + 0.001); | |
1917 | upperCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(upperCentrality - 0.001); | |
1918 | ||
1919 | // Check if the values look reasonable | |
1920 | if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1 | |
1921 | && upperCentralityBinLimit <= dataRebinned->GetAxis(iMult, 0)->GetNbins()) { | |
1922 | actualLowerCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinLowEdge(lowerCentralityBinLimit); | |
1923 | actualUpperCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinUpEdge(upperCentralityBinLimit); | |
1924 | ||
1925 | restrictCentralityAxis = kTRUE; | |
1926 | } | |
1927 | else { | |
1928 | std::cout << std::endl; | |
1929 | std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl; | |
1930 | return -1; | |
1931 | } | |
1932 | } | |
1933 | ||
1934 | std::cout << "centrality: "; | |
1935 | if (restrictCentralityAxis) { | |
1936 | std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl; | |
1937 | dataRebinned->SetRangeUser(iMult, lowerCentralityBinLimit, upperCentralityBinLimit, kTRUE); | |
1938 | data->SetRangeUser(iMult, lowerCentralityBinLimit, upperCentralityBinLimit, kTRUE); | |
1939 | } | |
1940 | else { | |
1941 | std::cout << "All" << std::endl; | |
1942 | } | |
1943 | ||
1944 | ||
1945 | ||
1946 | // If desired, restrict jetPt axis | |
1947 | Int_t lowerJetPtBinLimit = -1; | |
1948 | Int_t upperJetPtBinLimit = -2; // Integral(lowerJetBinLimit, uppJetBinLimit) will not be restricted if these values are kept | |
1949 | Bool_t restrictJetPtAxis = kFALSE; | |
1950 | Double_t actualLowerJetPt = -1.; | |
1951 | Double_t actualUpperJetPt = -1.; | |
1952 | ||
1953 | if (lowerJetPt >= 0 && upperJetPt >= 0) { | |
1954 | // Add subtract a very small number to avoid problems with values right on the border between to bins | |
1955 | lowerJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(lowerJetPt + 0.001); | |
1956 | upperJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(upperJetPt - 0.001); | |
1957 | ||
1958 | // Check if the values look reasonable | |
1959 | if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 && | |
1960 | upperJetPtBinLimit <= dataRebinned->GetAxis(iJetPt, 0)->GetNbins()) { | |
1961 | actualLowerJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinLowEdge(lowerJetPtBinLimit); | |
1962 | actualUpperJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinUpEdge(upperJetPtBinLimit); | |
1963 | ||
1964 | restrictJetPtAxis = kTRUE; | |
1965 | } | |
1966 | else { | |
1967 | std::cout << std::endl; | |
1968 | std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl; | |
1969 | return -1; | |
1970 | } | |
1971 | } | |
1972 | ||
1973 | std::cout << "jet pT: "; | |
1974 | if (restrictJetPtAxis) { | |
1975 | std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl; | |
1976 | dataRebinned->SetRangeUser(iJetPt, lowerJetPtBinLimit, upperJetPtBinLimit, kTRUE); | |
1977 | data->SetRangeUser(iJetPt, lowerJetPtBinLimit, upperJetPtBinLimit, kTRUE); | |
1978 | } | |
1979 | else { | |
1980 | std::cout << "All" << std::endl; | |
1981 | } | |
1982 | ||
1983 | if (restrictJetPtAxis && | |
1984 | (TMath::Abs(actualLowerJetPt - actualLowerJetPtYields) > 1e-3 || TMath::Abs(actualUpperJetPt - actualUpperJetPtYields) > 1e-3)) { | |
1985 | std::cout << "ERROR: Disagreement between jetPt bounds of data and efficiency!" << std::endl; | |
1986 | return -1; | |
1987 | } | |
1988 | ||
1989 | // If desired, restrict charge axis | |
1990 | std::cout << "Charge selection (efficiency): "; | |
1991 | if (chargeMode == kAllCharged) | |
1992 | std::cout << "All charged particles" << std::endl; | |
1993 | else if (chargeMode == kNegCharge) | |
1994 | std::cout << "Negative particles only" << std::endl; | |
1995 | else if (chargeMode == kPosCharge) | |
1996 | std::cout << "Positive particles only" << std::endl; | |
1997 | else { | |
1998 | std::cout << "Unknown -> ERROR" << std::endl; | |
1999 | return -1; | |
2000 | } | |
2001 | ||
2002 | const Bool_t restrictCharge = (chargeMode != kAllCharged); | |
2003 | ||
2004 | Int_t lowerChargeBinLimit = -1; | |
2005 | Int_t upperChargeBinLimit = -2; | |
2006 | Double_t actualLowerCharge = -999; | |
2007 | Double_t actualUpperCharge = -999; | |
2008 | ||
2009 | if (restrictCharge) { | |
2010 | // Add subtract a very small number to avoid problems with values right on the border between to bins | |
2011 | if (chargeMode == kNegCharge) { | |
2012 | lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(-1. + 0.001); | |
2013 | upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. - 0.001); | |
2014 | } | |
2015 | else if (chargeMode == kPosCharge) { | |
2016 | lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. + 0.001); | |
2017 | upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(1. - 0.001); | |
2018 | } | |
2019 | ||
2020 | // Check if the values look reasonable | |
2021 | if (lowerChargeBinLimit <= upperChargeBinLimit && lowerChargeBinLimit >= 1 | |
2022 | && upperChargeBinLimit <= dataRebinned->GetAxis(iCharge, 0)->GetNbins()) { | |
2023 | actualLowerCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinLowEdge(lowerChargeBinLimit); | |
2024 | actualUpperCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinUpEdge(upperChargeBinLimit); | |
2025 | ||
2026 | std::cout << "Charge range (efficiency): " << actualLowerCharge << " - " << actualUpperCharge << std::endl; | |
2027 | } | |
2028 | else { | |
2029 | std::cout << std::endl; | |
2030 | std::cout << "Requested charge range (efficiency) out of limits or upper and lower limit are switched!" << std::endl; | |
2031 | return -1; | |
2032 | } | |
2033 | ||
2034 | dataRebinned->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE); | |
2035 | data->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE); | |
2036 | } | |
2037 | ||
2038 | std::cout << std::endl; | |
2039 | ||
2040 | ||
2041 | // Load numJet histos from bhess_PID*.root file related to efficiency file. | |
2042 | TString pathNameDataMC = pathNameEfficiency; | |
2043 | pathNameDataMC.ReplaceAll("_efficiency", ""); | |
2044 | ||
2045 | TFile* fDataMC = TFile::Open(pathNameDataMC.Data()); | |
2046 | if (!fDataMC) { | |
2047 | std::cout << std::endl; | |
2048 | std::cout << "Failed to open file \"" << pathNameData.Data() << "\" to obtain num of rec/gen jets!" << std::endl; | |
2049 | ||
2050 | return -1; | |
2051 | } | |
2052 | ||
2053 | TString listName = pathNameDataMC; | |
2054 | listName.Replace(0, listName.Last('/') + 1, ""); | |
2055 | listName.ReplaceAll(".root", ""); | |
2056 | ||
2057 | TObjArray* histList = (TObjArray*)(fDataMC->Get(listName.Data())); | |
2058 | if (!histList) { | |
2059 | std::cout << std::endl; | |
2060 | std::cout << "Failed to load list \"" << listName.Data() << "\" to obtain num of rec/gen jets!" << std::endl; | |
2061 | return -1; | |
2062 | } | |
2063 | ||
2064 | hNjetsGen = (TH2D*)histList->FindObject("fh2FFJetPtGen"); | |
2065 | hNjetsRec = (TH2D*)histList->FindObject("fh2FFJetPtRec"); | |
2066 | ||
2067 | if (!hNjetsRec || !hNjetsGen) { | |
2068 | std::cout << "Failed to load number of jets histos!" << std::endl; | |
2069 | ||
2070 | ||
2071 | // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other | |
2072 | // period is considered; also: No multiplicity information) | |
2073 | TString pathEfficiency = pathNameEfficiency; | |
2074 | pathEfficiency.Replace(pathEfficiency.Last('/'), pathEfficiency.Length(), ""); | |
2075 | TString pathBackward = Form("%s/AnalysisResults.root", pathEfficiency.Data()); | |
2076 | TFile* fBackward = TFile::Open(pathBackward.Data()); | |
2077 | ||
2078 | TString dirDataInFile = ""; | |
2079 | TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0; | |
2080 | ||
2081 | TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0; | |
2082 | ||
2083 | TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCutsInc") : 0x0; | |
2084 | TH1D* hFFJetPtGen = list ? (TH1D*)list->FindObject("fh1FFJetPtGenInc") : 0x0; | |
2085 | ||
2086 | if (hFFJetPtRec && hFFJetPtGen) { | |
2087 | printf("***WARNING: For backward compatibility, using file \"%s\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n", | |
2088 | pathBackward.Data()); | |
2089 | printf("ALSO: Using Njets for inclusive jets!!!!\n"); | |
2090 | ||
2091 | hNjetsRec = new TH2D("fh2FFJetPtRec", "", 1, -1, 1, dataRebinned->GetAxis(iJetPt, 0)->GetNbins(), | |
2092 | dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray()); | |
2093 | ||
2094 | for (Int_t iJet = 1; iJet <= hNjetsRec->GetNbinsY(); iJet++) { | |
2095 | Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinLowEdge(iJet) + 1e-3); | |
2096 | Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinUpEdge(iJet) - 1e-3); | |
2097 | hNjetsRec->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin)); | |
2098 | } | |
2099 | ||
2100 | hNjetsGen = new TH2D("fh2FFJetPtGen", "", 1, -1, 1, dataRebinned->GetAxis(iJetPt, 0)->GetNbins(), | |
2101 | dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray()); | |
2102 | ||
2103 | for (Int_t iJet = 1; iJet <= hNjetsGen->GetNbinsY(); iJet++) { | |
2104 | Int_t lowerBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinLowEdge(iJet) + 1e-3); | |
2105 | Int_t upperBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinUpEdge(iJet) - 1e-3); | |
2106 | hNjetsGen->SetBinContent(1, iJet, hFFJetPtGen->Integral(lowerBin, upperBin)); | |
2107 | } | |
2108 | } | |
2109 | ||
2110 | if (!hNjetsRec || ! hNjetsGen) | |
2111 | return -1; | |
2112 | } | |
2113 | ||
2114 | // For normalisation to number of jets | |
2115 | // NOTE: These numbers are for the efficiency only! The data will be normalised to its own number!!! | |
2116 | const Double_t nJetsGen = hNjetsGen ? hNjetsGen->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit, | |
2117 | upperJetPtBinLimit) : 1.; | |
2118 | const Double_t nJetsRec = hNjetsRec ? hNjetsRec->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit, | |
2119 | upperJetPtBinLimit) : 1.; | |
2120 | ||
2121 | ||
2122 | ||
2123 | // Secondary correction | |
2124 | AliCFEffGrid* sec = new AliCFEffGrid("sec", "Secondary Contamination", *dataRebinned); | |
2125 | sec->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObs); | |
2126 | ||
2127 | TH1 *hSecAllPt = 0x0, *hSecID2Pt = 0x0; | |
2128 | TH1D* hSecPt[AliPID::kSPECIES] = { 0x0, }; | |
2129 | TH1D* hSecToPiRatioPt[AliPID::kSPECIES] = { 0x0, }; | |
2130 | TCanvas *cSecPt = 0x0; | |
2131 | TCanvas *cSecToPiRatioPt = 0x0; | |
2132 | ||
2133 | TH1 *hSecAllZ = 0x0, *hSecID2Z = 0x0; | |
2134 | TH1D* hSecZ[AliPID::kSPECIES] = { 0x0, }; | |
2135 | TH1D* hSecToPiRatioZ[AliPID::kSPECIES] = { 0x0, }; | |
2136 | TCanvas *cSecZ = 0x0; | |
2137 | TCanvas *cSecToPiRatioZ = 0x0; | |
2138 | ||
2139 | TH1 *hSecAllXi = 0x0, *hSecID2Xi = 0x0; | |
2140 | TH1D* hSecXi[AliPID::kSPECIES] = { 0x0, }; | |
2141 | TH1D* hSecToPiRatioXi[AliPID::kSPECIES] = { 0x0, }; | |
2142 | TCanvas *cSecXi = 0x0; | |
2143 | TCanvas *cSecToPiRatioXi = 0x0; | |
2144 | ||
2145 | // Secondary correction with strangeness scaling | |
2146 | AliCFEffGrid* secStrangeScale = new AliCFEffGrid("secStrangeScale", "Secondary Contamination with strangeness scaling", | |
2147 | *dataRebinned); | |
2148 | secStrangeScale->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObsStrangenessScaled); | |
2149 | ||
2150 | TH1 *hSecSSAllPt = 0x0, *hSecSSID2Pt = 0x0; | |
2151 | TH1D* hSecSSPt[AliPID::kSPECIES] = { 0x0, }; | |
2152 | TH1D* hSecSSToPiRatioPt[AliPID::kSPECIES] = { 0x0, }; | |
2153 | TCanvas *cSecSSPt = 0x0; | |
2154 | TCanvas *cSecSSToPiRatioPt = 0x0; | |
2155 | ||
2156 | TH1 *hSecSSAllZ = 0x0, *hSecSSID2Z = 0x0; | |
2157 | TH1D* hSecSSZ[AliPID::kSPECIES] = { 0x0, }; | |
2158 | TH1D* hSecSSToPiRatioZ[AliPID::kSPECIES] = { 0x0, }; | |
2159 | TCanvas *cSecSSZ = 0x0; | |
2160 | TCanvas *cSecSSToPiRatioZ = 0x0; | |
2161 | ||
2162 | TH1 *hSecSSAllXi = 0x0, *hSecSSID2Xi = 0x0; | |
2163 | TH1D* hSecSSXi[AliPID::kSPECIES] = { 0x0, }; | |
2164 | TH1D* hSecSSToPiRatioXi[AliPID::kSPECIES] = { 0x0, }; | |
2165 | TCanvas *cSecSSXi = 0x0; | |
2166 | TCanvas *cSecSSToPiRatioXi = 0x0; | |
2167 | ||
2168 | ||
2169 | const Double_t upperTrackPt = restrictJetPtAxis ? actualUpperJetPt : 50.; | |
2170 | ||
2171 | TString strangenessString[2] = { "", "SS" }; | |
2172 | ||
2173 | // Get the secondary contamination vs. pT,z,xi for each species w/ and w/o strangeness scaling | |
2174 | ||
2175 | for (Int_t type = 0; type < kNtypes; type++) { | |
2176 | for (Int_t strangenessScaling = 0; strangenessScaling < 2; strangenessScaling++) { | |
2177 | TH1 **hSecAll, **hSecID2; | |
2178 | TH1D **hSec, **hSecToPiRatio; | |
2179 | TCanvas **cSec, **cSecToPiRatio; | |
2180 | ||
2181 | Int_t iDimProj = -1; | |
2182 | ||
2183 | AliCFEffGrid* secCurr = strangenessScaling ? secStrangeScale : sec; | |
2184 | ||
2185 | if (type == kTrackPt) { | |
2186 | hSecAll = strangenessScaling ? &hSecSSAllPt : &hSecAllPt; | |
2187 | hSecID2 = strangenessScaling ? &hSecSSID2Pt : &hSecID2Pt; | |
2188 | hSec = strangenessScaling ? &hSecSSPt[0] : &hSecPt[0]; | |
2189 | hSecToPiRatio = strangenessScaling ? &hSecSSToPiRatioPt[0] : &hSecToPiRatioPt[0]; | |
2190 | cSec = strangenessScaling ? &cSecSSPt : &cSecPt; | |
2191 | cSecToPiRatio = strangenessScaling ? &cSecSSToPiRatioPt : &cSecToPiRatioPt; | |
2192 | iDimProj = iPt; | |
2193 | } | |
2194 | else if (type == kZ) { | |
2195 | hSecAll = strangenessScaling ? &hSecSSAllZ : &hSecAllZ; | |
2196 | hSecID2 = strangenessScaling ? &hSecSSID2Z : &hSecID2Z; | |
2197 | hSec = strangenessScaling ? &hSecSSZ[0] : &hSecZ[0]; | |
2198 | hSecToPiRatio = strangenessScaling ? &hSecSSToPiRatioZ[0] : &hSecToPiRatioZ[0]; | |
2199 | cSec = strangenessScaling ? &cSecSSZ : &cSecZ; | |
2200 | cSecToPiRatio = strangenessScaling ? &cSecSSToPiRatioZ : &cSecToPiRatioZ; | |
2201 | iDimProj = iZ; | |
2202 | } | |
2203 | else if (type == kXi) { | |
2204 | hSecAll = strangenessScaling ? &hSecSSAllXi : &hSecAllXi; | |
2205 | hSecID2 = strangenessScaling ? &hSecSSID2Xi : &hSecID2Xi; | |
2206 | hSec = strangenessScaling ? &hSecSSXi[0] : &hSecXi[0]; | |
2207 | hSecToPiRatio = strangenessScaling ? &hSecSSToPiRatioXi[0] : &hSecToPiRatioXi[0]; | |
2208 | cSec = strangenessScaling ? &cSecSSXi : &cSecXi; | |
2209 | cSecToPiRatio = strangenessScaling ? &cSecSSToPiRatioXi : &cSecToPiRatioXi; | |
2210 | iDimProj = iXi; | |
2211 | } | |
2212 | else | |
2213 | continue; | |
2214 | ||
2215 | *hSecAll = secCurr->Project(iDimProj); | |
2216 | (*hSecAll)->SetName(Form("hSec%s_%s_all", strangenessString[strangenessScaling].Data(), titles[type].Data())); | |
2217 | (*hSecAll)->SetTitle("All"); | |
2218 | (*hSecAll)->SetMarkerStyle(24); | |
2219 | (*hSecAll)->SetLineWidth(2.); | |
2220 | (*hSecAll)->SetStats(kFALSE); | |
2221 | (*hSecAll)->GetYaxis()->SetTitle("Primary Fraction"); | |
2222 | (*hSecID2) = (TH2D*)secCurr->Project(iDimProj, iMCid); | |
2223 | (*hSecID2)->SetName(Form("hSecID2_%s", titles[type].Data())); | |
2224 | ||
2225 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
2226 | hSec[species] = ((TH2D*)*hSecID2)->ProjectionX(Form("hSec%s_%s_%s", strangenessString[strangenessScaling].Data(), | |
2227 | titles[type].Data(), AliPID::ParticleShortName(species)), | |
2228 | species + 1, species + 1, "e"); | |
2229 | hSec[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
2230 | hSec[species]->SetLineColor(getLineColorAliPID(species)); | |
2231 | hSec[species]->SetMarkerColor(getLineColorAliPID(species)); | |
2232 | hSec[species]->SetMarkerStyle(24); | |
2233 | hSec[species]->SetLineWidth(2.); | |
2234 | hSec[species]->GetYaxis()->SetRangeUser(0., 1.01); | |
2235 | hSec[species]->SetStats(kFALSE); | |
2236 | hSec[species]->GetYaxis()->SetTitle("Primary Fraction"); | |
2237 | } | |
2238 | ||
2239 | *cSec = new TCanvas(Form("cSec%s_%s", strangenessString[strangenessScaling].Data(), titles[type].Data()), | |
2240 | Form("Primary fraction for different species%s", strangenessScaling ? " (strangeness scaled)" : ""), | |
2241 | 100, 10, 1200, 800); | |
2242 | (*cSec)->cd(); | |
2243 | (*cSec)->SetGridx(1); | |
2244 | (*cSec)->SetGridy(1); | |
2245 | ||
2246 | if (type == kTrackPt) | |
2247 | (*hSecAll)->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2248 | ||
2249 | (*hSecAll)->Draw("E1"); | |
2250 | ||
2251 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2252 | if (type == kTrackPt) | |
2253 | hSec[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2254 | ||
2255 | hSec[i]->Draw("E1 same"); | |
2256 | } | |
2257 | (*cSec)->BuildLegend()->SetFillColor(kWhite); | |
2258 | ||
2259 | ClearTitleFromHistoInCanvas(*cSec); | |
2260 | ||
2261 | // To-pi ratios | |
2262 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
2263 | if (species == AliPID::kPion) | |
2264 | continue; // Do not consider pion-to-pion ratio | |
2265 | ||
2266 | hSecToPiRatio[species] = new TH1D(*hSec[species]); | |
2267 | hSecToPiRatio[species]->Reset(); | |
2268 | hSecToPiRatio[species]->SetName(Form("hSec%sToPionRatio_%s_%s", strangenessString[strangenessScaling].Data(), | |
2269 | titles[type].Data(), AliPID::ParticleShortName(species))); | |
2270 | hSecToPiRatio[species]->SetTitle(Form("%s/#pi", AliPID::ParticleLatexName(species))); | |
2271 | hSecToPiRatio[species]->SetLineColor(getLineColorAliPID(species)); | |
2272 | hSecToPiRatio[species]->SetMarkerColor(getLineColorAliPID(species)); | |
2273 | hSecToPiRatio[species]->SetMarkerStyle(24); | |
2274 | hSecToPiRatio[species]->SetLineWidth(2.); | |
2275 | hSecToPiRatio[species]->SetStats(kFALSE); | |
2276 | hSecToPiRatio[species]->GetYaxis()->SetTitle("Primary Fraction of Ratio"); | |
2277 | ||
2278 | // Samples for different species are independent, so just divide correction factors | |
2279 | hSecToPiRatio[species]->Divide(hSec[species], hSec[AliPID::kPion], 1., 1., ""); | |
2280 | hSecToPiRatio[species]->GetYaxis()->SetRangeUser(0., 1.1); | |
2281 | } | |
2282 | ||
2283 | *cSecToPiRatio = new TCanvas(Form("cSec%sToPiRatio_%s", strangenessString[strangenessScaling].Data(), titles[type].Data()), | |
2284 | Form("Primary fraction of to-#pi-ratio for different species%s", | |
2285 | strangenessScaling ? " (strangeness scaled)" : ""), | |
2286 | 100, 10, 1200, 800); | |
2287 | (*cSecToPiRatio)->cd(); | |
2288 | (*cSecToPiRatio)->SetGridx(1); | |
2289 | (*cSecToPiRatio)->SetGridy(1); | |
2290 | ||
2291 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2292 | if (i == AliPID::kPion) | |
2293 | continue; | |
2294 | ||
2295 | if (type == kTrackPt) | |
2296 | hSecToPiRatio[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2297 | ||
2298 | hSecToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same")); | |
2299 | } | |
2300 | (*cSecToPiRatio)->BuildLegend()->SetFillColor(kWhite); | |
2301 | ||
2302 | ClearTitleFromHistoInCanvas(*cSecToPiRatio); | |
2303 | } | |
2304 | } | |
2305 | ||
2306 | ||
2307 | // Construct the efficiency grid from the data container | |
2308 | // NOTE: Here, the scale factors 1/nJet are different. This effects the error of the ratio. | |
2309 | // Thus, the histograms must be extracted, scaled and THEN divided to arrive at the correction | |
2310 | // factor and the standard "calculateEfficiency" function CANNOT be used (multiplying the ratio | |
2311 | // of the scale factor afterwards would result in wrong errors) | |
2312 | ||
2313 | ||
2314 | ||
2315 | // If desired, apply the GEANT/FLUKA correction first | |
2316 | // NOTE: This will change dataRebinned! If anything else should also be calculated, it must be done before. | |
2317 | // Otherwise, the GEANT/FLUKA correction factor for the efficiency will also affect it! | |
2318 | const Int_t genStepEff = kStepGenWithGenCuts; | |
2319 | if (correctGeantFluka) { | |
2320 | printf("Applying GEANT/FLUKA correction...\n"); | |
2321 | if (!geantFlukaCorrection(dataRebinned, genStepEff)) { | |
2322 | printf("GEANT/FLUKA correction could not be applied!\n"); | |
2323 | return kFALSE; | |
2324 | } | |
2325 | } | |
2326 | ||
2327 | ||
2328 | // Either one can take kStepRecWithGenCutsMeasuredObs or, what I prefer, one can take | |
2329 | // kStepRecWithRecCutsMeasuredObsPrimaries => The difference is only the eta cut, which is on the rec level | |
2330 | // in the latter case, i.e. one corrects for eta resolution (although the effect should be very small) | |
2331 | const Int_t recStepEff = kStepRecWithRecCutsMeasuredObsPrimaries; | |
2332 | ||
2333 | // The efficiency vs. pT, z, xi | |
2334 | TH1* hSpecPt_RecPrimaries = 0x0; | |
2335 | TH1* hSpecPt_GenPrimaries = 0x0; | |
2336 | TH1D* hEfficiencyAllPt = 0x0; | |
2337 | ||
2338 | TH2* hSpecID2Pt_RecPrimaries = 0x0; | |
2339 | TH2* hSpecID2Pt_GenPrimaries = 0x0; | |
2340 | TH2D* hEffID2Pt = 0x0; | |
2341 | ||
2342 | TH1D* hSpecIDPt_GenPrimaries[AliPID::kSPECIES] = { 0x0, }; | |
2343 | TH1D* hSpecIDPt_RecPrimaries[AliPID::kSPECIES] = { 0x0, }; | |
2344 | TH1D* hEfficiencyPt[AliPID::kSPECIES] = { 0x0, }; | |
2345 | TH1D* hEfficiencyToPiRatioPt[AliPID::kSPECIES] = { 0x0, }; | |
2346 | ||
2347 | TCanvas *cSpecPt = 0x0, *cEffPt = 0x0, *cEffToPiRatioPt = 0x0; | |
2348 | ||
2349 | ||
2350 | TH1* hSpecZ_RecPrimaries = 0x0; | |
2351 | TH1* hSpecZ_GenPrimaries = 0x0; | |
2352 | TH1D* hEfficiencyAllZ = 0x0; | |
2353 | ||
2354 | TH2* hSpecID2Z_RecPrimaries = 0x0; | |
2355 | TH2* hSpecID2Z_GenPrimaries = 0x0; | |
2356 | TH2D* hEffID2Z = 0x0; | |
2357 | ||
2358 | TH1D* hSpecIDZ_GenPrimaries[AliPID::kSPECIES] = { 0x0, }; | |
2359 | TH1D* hSpecIDZ_RecPrimaries[AliPID::kSPECIES] = { 0x0, }; | |
2360 | TH1D* hEfficiencyZ[AliPID::kSPECIES] = { 0x0, }; | |
2361 | TH1D* hEfficiencyToPiRatioZ[AliPID::kSPECIES] = { 0x0, }; | |
2362 | ||
2363 | TCanvas *cSpecZ = 0x0, *cEffZ = 0x0, *cEffToPiRatioZ = 0x0; | |
2364 | ||
2365 | ||
2366 | TH1* hSpecXi_RecPrimaries = 0x0; | |
2367 | TH1* hSpecXi_GenPrimaries = 0x0; | |
2368 | TH1D* hEfficiencyAllXi = 0x0; | |
2369 | ||
2370 | TH2* hSpecID2Xi_RecPrimaries = 0x0; | |
2371 | TH2* hSpecID2Xi_GenPrimaries = 0x0; | |
2372 | TH2D* hEffID2Xi = 0x0; | |
2373 | ||
2374 | TH1D* hSpecIDXi_GenPrimaries[AliPID::kSPECIES] = { 0x0, }; | |
2375 | TH1D* hSpecIDXi_RecPrimaries[AliPID::kSPECIES] = { 0x0, }; | |
2376 | TH1D* hEfficiencyXi[AliPID::kSPECIES] = { 0x0, }; | |
2377 | TH1D* hEfficiencyToPiRatioXi[AliPID::kSPECIES] = { 0x0, }; | |
2378 | ||
2379 | TCanvas *cSpecXi = 0x0, *cEffXi = 0x0, *cEffToPiRatioXi = 0x0; | |
2380 | ||
2381 | TString titleYaxis[kNtypes]; | |
2382 | titleYaxis[kTrackPt] = "1/N_{Jets} dN/dP_{T} (GeV/c)^{-1}"; | |
2383 | titleYaxis[kZ] = "1/N_{Jets} dN/dz"; | |
2384 | titleYaxis[kXi] = "1/N_{Jets} dN/d#xi"; | |
2385 | ||
2386 | for (Int_t type = 0; type < kNtypes; type++) { | |
2387 | TH1** hSpec_RecPrimaries; | |
2388 | TH1** hSpec_GenPrimaries; | |
2389 | TH1D** hEfficiencyAll; | |
2390 | ||
2391 | TH2** hSpecID2_RecPrimaries; | |
2392 | TH2** hSpecID2_GenPrimaries; | |
2393 | TH2D** hEffID2; | |
2394 | ||
2395 | TH1D** hSpecID_GenPrimaries; | |
2396 | TH1D** hSpecID_RecPrimaries; | |
2397 | TH1D** hEfficiency; | |
2398 | TH1D** hEfficiencyToPiRatio; | |
2399 | ||
2400 | TCanvas **cSpec, **cEff, **cEffToPiRatio; | |
2401 | ||
2402 | Int_t iDimProj = -1; | |
2403 | ||
2404 | ||
2405 | if (type == kTrackPt) { | |
2406 | hSpec_RecPrimaries = &hSpecPt_RecPrimaries; | |
2407 | hSpec_GenPrimaries = &hSpecPt_GenPrimaries; | |
2408 | hEfficiencyAll = &hEfficiencyAllPt; | |
2409 | ||
2410 | hSpecID2_RecPrimaries = &hSpecID2Pt_RecPrimaries; | |
2411 | hSpecID2_GenPrimaries = &hSpecID2Pt_GenPrimaries; | |
2412 | hEffID2 = &hEffID2Pt; | |
2413 | ||
2414 | hSpecID_RecPrimaries = &hSpecIDPt_RecPrimaries[0]; | |
2415 | hSpecID_GenPrimaries = &hSpecIDPt_GenPrimaries[0]; | |
2416 | hEfficiency = &hEfficiencyPt[0]; | |
2417 | hEfficiencyToPiRatio = &hEfficiencyToPiRatioPt[0]; | |
2418 | ||
2419 | cSpec = &cSpecPt; | |
2420 | cEff = &cEffPt; | |
2421 | cEffToPiRatio = &cEffToPiRatioPt; | |
2422 | ||
2423 | iDimProj = iPt; | |
2424 | } | |
2425 | else if (type == kZ) { | |
2426 | hSpec_RecPrimaries = &hSpecZ_RecPrimaries; | |
2427 | hSpec_GenPrimaries = &hSpecZ_GenPrimaries; | |
2428 | hEfficiencyAll = &hEfficiencyAllZ; | |
2429 | ||
2430 | hSpecID2_RecPrimaries = &hSpecID2Z_RecPrimaries; | |
2431 | hSpecID2_GenPrimaries = &hSpecID2Z_GenPrimaries; | |
2432 | hEffID2 = &hEffID2Z; | |
2433 | ||
2434 | hSpecID_RecPrimaries = &hSpecIDZ_RecPrimaries[0]; | |
2435 | hSpecID_GenPrimaries = &hSpecIDZ_GenPrimaries[0]; | |
2436 | hEfficiency = &hEfficiencyZ[0]; | |
2437 | hEfficiencyToPiRatio = &hEfficiencyToPiRatioZ[0]; | |
2438 | ||
2439 | cSpec = &cSpecZ; | |
2440 | cEff = &cEffZ; | |
2441 | cEffToPiRatio = &cEffToPiRatioZ; | |
2442 | ||
2443 | iDimProj = iZ; | |
2444 | } | |
2445 | else if (type == kXi) { | |
2446 | hSpec_RecPrimaries = &hSpecXi_RecPrimaries; | |
2447 | hSpec_GenPrimaries = &hSpecXi_GenPrimaries; | |
2448 | hEfficiencyAll = &hEfficiencyAllXi; | |
2449 | ||
2450 | hSpecID2_RecPrimaries = &hSpecID2Xi_RecPrimaries; | |
2451 | hSpecID2_GenPrimaries = &hSpecID2Xi_GenPrimaries; | |
2452 | hEffID2 = &hEffID2Xi; | |
2453 | ||
2454 | hSpecID_RecPrimaries = &hSpecIDXi_RecPrimaries[0]; | |
2455 | hSpecID_GenPrimaries = &hSpecIDXi_GenPrimaries[0]; | |
2456 | hEfficiency = &hEfficiencyXi[0]; | |
2457 | hEfficiencyToPiRatio = &hEfficiencyToPiRatioXi[0]; | |
2458 | ||
2459 | cSpec = &cSpecXi; | |
2460 | cEff = &cEffXi; | |
2461 | cEffToPiRatio = &cEffToPiRatioXi; | |
2462 | ||
2463 | iDimProj = iXi; | |
2464 | } | |
2465 | else | |
2466 | continue; | |
2467 | ||
2468 | ||
2469 | // Divide steps: recStepEff (=kStepRecWithRecCutsMeasuredObsPrimaries) / genStepEff (=kStepGenWithGenCuts) | |
2470 | *hSpec_RecPrimaries = dataRebinned->GetGrid(recStepEff)->Project(iDimProj); | |
2471 | setupHist(*hSpec_RecPrimaries, Form("hSpec_%s_RecPrimaries", titles[type].Data()), "All, rec", "", | |
2472 | titleYaxis[type].Data(), kBlack); | |
2473 | normaliseHist(*hSpec_RecPrimaries, nJetsRec > 0 ? 1. / nJetsRec : 0.); | |
2474 | (*hSpec_RecPrimaries)->SetLineWidth(2.); | |
2475 | (*hSpec_RecPrimaries)->SetMarkerStyle(20); | |
2476 | ||
2477 | *hSpec_GenPrimaries = dataRebinned->GetGrid(genStepEff)->Project(iDimProj); | |
2478 | setupHist(*hSpec_GenPrimaries, Form("hSpec_%s_GenPrimaries", titles[type].Data()), "All, gen", "", | |
2479 | titleYaxis[type].Data(), kBlack); | |
2480 | normaliseHist(*hSpec_GenPrimaries, nJetsGen > 0 ? 1. / nJetsGen : 0.); | |
2481 | (*hSpec_GenPrimaries)->SetLineWidth(2.); | |
2482 | (*hSpec_GenPrimaries)->SetLineStyle(2); | |
2483 | ||
2484 | *hEfficiencyAll = new TH1D(*(TH1D*)(*hSpec_RecPrimaries)); | |
2485 | (*hEfficiencyAll)->SetName(Form("hEfficiencyAll_%s", titles[type].Data())); | |
2486 | (*hEfficiencyAll)->GetYaxis()->SetTitle(""); | |
2487 | (*hEfficiencyAll)->SetTitle("All"); | |
2488 | (*hEfficiencyAll)->SetLineWidth(2.0); | |
2489 | (*hEfficiencyAll)->SetStats(kFALSE); | |
2490 | (*hEfficiencyAll)->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution"); | |
2491 | (*hEfficiencyAll)->GetYaxis()->SetTitleSize(0.05); | |
2492 | (*hEfficiencyAll)->Divide(*hSpec_RecPrimaries, *hSpec_GenPrimaries, 1., 1., "B"); | |
2493 | (*hEfficiencyAll)->GetYaxis()->SetRangeUser(0., 2.0); | |
2494 | ||
2495 | *hSpecID2_RecPrimaries = (TH2*)dataRebinned->GetGrid(recStepEff)->Project(iDimProj, iMCid); | |
2496 | (*hSpecID2_RecPrimaries)->SetName(Form("hSpecID2_RecPrimaries_%s", titles[type].Data())); | |
2497 | (*hSpecID2_RecPrimaries)->Scale(nJetsRec > 0 ? 1. / nJetsRec : 0.); | |
2498 | ||
2499 | *hSpecID2_GenPrimaries = (TH2*)dataRebinned->GetGrid(genStepEff)->Project(iDimProj, iMCid); | |
2500 | (*hSpecID2_GenPrimaries)->SetName(Form("hSpecID2_GenPrimaries_%s", titles[type].Data())); | |
2501 | (*hSpecID2_GenPrimaries)->Scale(nJetsGen > 0 ? 1. / nJetsGen : 0.); | |
2502 | ||
2503 | *hEffID2 = new TH2D(*(TH2D*)(*hSpecID2_RecPrimaries)); | |
2504 | (*hEffID2)->SetName(Form("hEff_%s", titles[type].Data())); | |
2505 | (*hEffID2)->SetTitle("All"); | |
2506 | (*hEffID2)->Divide(*hSpecID2_RecPrimaries, *hSpecID2_GenPrimaries, 1., 1., "B"); | |
2507 | ||
2508 | // Get the efficiencies vs. pT,z,xi for each species | |
2509 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
2510 | hEfficiency[species] = ((TH2D*)*hEffID2)->ProjectionX(Form("hEfficiency_%s_%s", titles[type].Data(), | |
2511 | AliPID::ParticleShortName(species)), | |
2512 | species + 1, species + 1, "e"); | |
2513 | setupHist(hEfficiency[species], "", Form("%s", AliPID::ParticleLatexName(species)), "", "", | |
2514 | getLineColorAliPID(species)); | |
2515 | hEfficiency[species]->SetLineWidth(2.); | |
2516 | hEfficiency[species]->GetYaxis()->SetRangeUser(0., 2.0); | |
2517 | hEfficiency[species]->SetStats(kFALSE); | |
2518 | hEfficiency[species]->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution"); | |
2519 | hEfficiency[species]->GetYaxis()->SetTitleSize(0.05); | |
2520 | ||
2521 | hSpecID_GenPrimaries[species] = ((TH2D*)*hSpecID2_GenPrimaries)->ProjectionX(Form("hSpecID_%s_gen_%s", | |
2522 | titles[type].Data(), | |
2523 | AliPID::ParticleShortName(species)), | |
2524 | species + 1, species + 1, "e"); | |
2525 | setupHist(hSpecID_GenPrimaries[species], "", Form("%s, gen", AliPID::ParticleLatexName(species)), "", | |
2526 | titleYaxis[type].Data(), getLineColorAliPID(species)); | |
2527 | normaliseHist(hSpecID_GenPrimaries[species], 1.); | |
2528 | hSpecID_GenPrimaries[species]->SetLineWidth(2.); | |
2529 | hSpecID_GenPrimaries[species]->SetLineStyle(2.); | |
2530 | hSpecID_GenPrimaries[species]->SetStats(kFALSE); | |
2531 | ||
2532 | ||
2533 | hSpecID_RecPrimaries[species] = ((TH2D*)*hSpecID2_RecPrimaries)->ProjectionX(Form("hSpecID_%s_rec_%s", | |
2534 | titles[type].Data(), | |
2535 | AliPID::ParticleShortName(species)), | |
2536 | species + 1, species + 1, "e"); | |
2537 | setupHist(hSpecID_RecPrimaries[species], "", Form("%s, rec", AliPID::ParticleLatexName(species)), "", | |
2538 | titleYaxis[type].Data(), getLineColorAliPID(species)); | |
2539 | normaliseHist(hSpecID_RecPrimaries[species], 1.); | |
2540 | hSpecID_RecPrimaries[species]->SetLineWidth(2.); | |
2541 | hSpecID_RecPrimaries[species]->SetMarkerStyle(20); | |
2542 | hSpecID_RecPrimaries[species]->SetStats(kFALSE); | |
2543 | } | |
2544 | ||
2545 | *cSpec = new TCanvas(Form("cSpec_%s", titles[type].Data()), "Spectra for different species", 0, 300, 900, 900); | |
2546 | (*cSpec)->cd(); | |
2547 | (*cSpec)->SetGridx(1); | |
2548 | (*cSpec)->SetGridy(1); | |
2549 | (*cSpec)->SetLogy(1); | |
2550 | ||
2551 | if (type == kTrackPt) | |
2552 | (*hSpec_GenPrimaries)->GetYaxis()->SetRangeUser(1e-8, 1e2); | |
2553 | else if (type == kZ) | |
2554 | (*hSpec_GenPrimaries)->GetYaxis()->SetRangeUser(1e-5, 1e3); | |
2555 | else if (type == kXi) | |
2556 | (*hSpec_GenPrimaries)->GetYaxis()->SetRangeUser(1e-5, 1e3); | |
2557 | ||
2558 | ||
2559 | if (type == kTrackPt) | |
2560 | (*hSpec_GenPrimaries)->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2561 | (*hSpec_GenPrimaries)->Draw("E1"); | |
2562 | ||
2563 | if (type == kTrackPt) | |
2564 | (*hSpec_RecPrimaries)->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2565 | (*hSpec_RecPrimaries)->Draw("E1 same"); | |
2566 | ||
2567 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2568 | if (type == kTrackPt) { | |
2569 | hSpecID_GenPrimaries[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2570 | hSpecID_RecPrimaries[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2571 | } | |
2572 | ||
2573 | hSpecIDPt_GenPrimaries[i]->Draw("E1 same"); | |
2574 | hSpecIDPt_RecPrimaries[i]->Draw("E1 same"); | |
2575 | } | |
2576 | TLegend* leg = (*cSpec)->BuildLegend(); | |
2577 | leg->SetFillColor(kWhite); | |
2578 | leg->SetNColumns(2); | |
2579 | ||
2580 | ClearTitleFromHistoInCanvas(*cSpec); | |
2581 | ||
2582 | ||
2583 | *cEff = new TCanvas(Form("cEff_%s", titles[type].Data()), "Efficiency x Acceptance x pT Resolution for different species", | |
2584 | 0, 300, 900, 900); | |
2585 | (*cEff)->cd(); | |
2586 | (*cEff)->SetGridx(1); | |
2587 | (*cEff)->SetGridy(1); | |
2588 | ||
2589 | if (type == kTrackPt) | |
2590 | (*hEfficiencyAll)->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2591 | (*hEfficiencyAll)->Draw("E1"); | |
2592 | ||
2593 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2594 | if (type == kTrackPt) | |
2595 | hEfficiency[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2596 | hEfficiency[i]->Draw("E1 same"); | |
2597 | } | |
2598 | (*cEff)->BuildLegend()->SetFillColor(kWhite); | |
2599 | ||
2600 | ClearTitleFromHistoInCanvas(*cEff); | |
2601 | ||
2602 | ||
2603 | // To-pi ratios | |
2604 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
2605 | if (species == AliPID::kPion) | |
2606 | continue; // Do not consider pion-to-pion ratio | |
2607 | ||
2608 | hEfficiencyToPiRatio[species] = new TH1D(*hEfficiency[species]); | |
2609 | hEfficiencyToPiRatio[species]->Reset(); | |
2610 | hEfficiencyToPiRatio[species]->SetName(Form("hEfficiency_ToPionRatio_%s_%s", titles[type].Data(), | |
2611 | AliPID::ParticleShortName(species))); | |
2612 | ||
2613 | setupHist(hEfficiencyToPiRatio[species], "", Form("%s/#pi", AliPID::ParticleLatexName(species)), "", "", | |
2614 | getLineColorAliPID(species)); | |
2615 | hEfficiencyToPiRatio[species]->SetLineWidth(2.); | |
2616 | hEfficiencyToPiRatio[species]->SetStats(kFALSE); | |
2617 | hEfficiencyToPiRatio[species]->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution"); | |
2618 | hEfficiencyToPiRatio[species]->GetYaxis()->SetTitleSize(0.05); | |
2619 | ||
2620 | // Samples for different species are independent, so just divide correction factors | |
2621 | hEfficiencyToPiRatio[species]->Divide(hEfficiency[species], hEfficiency[AliPID::kPion], 1., 1., ""); | |
2622 | hEfficiencyToPiRatio[species]->GetYaxis()->SetRangeUser(0., 2.0); | |
2623 | } | |
2624 | ||
2625 | *cEffToPiRatio = new TCanvas(Form("cEffToPiRatio_%s", titles[type].Data()), | |
2626 | "Efficiency x Acceptance x pT Resolution for different to-#pi-ratios", | |
2627 | 100, 10, 1200, 800); | |
2628 | (*cEffToPiRatio)->cd(); | |
2629 | (*cEffToPiRatio)->SetGridx(1); | |
2630 | (*cEffToPiRatio)->SetGridy(1); | |
2631 | ||
2632 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2633 | if (i == AliPID::kPion) | |
2634 | continue; | |
2635 | ||
2636 | if (type == kTrackPt) | |
2637 | hEfficiencyToPiRatio[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt); | |
2638 | ||
2639 | hEfficiencyToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same")); | |
2640 | } | |
2641 | (*cEffToPiRatio)->BuildLegend()->SetFillColor(kWhite); | |
2642 | ||
2643 | ClearTitleFromHistoInCanvas(*cEffToPiRatio); | |
2644 | } | |
2645 | ||
2646 | // Save results to file | |
2647 | TString chargeString = ""; | |
2648 | if (chargeMode == kPosCharge) | |
2649 | chargeString = "_posCharge"; | |
2650 | else if (chargeMode == kNegCharge) | |
2651 | chargeString = "_negCharge"; | |
2652 | ||
2653 | TString saveFileName = pathNameData; | |
2654 | saveFileName.ReplaceAll(Form("%s/", pathData.Data()), ""); | |
2655 | saveFileName.Prepend("output_EfficiencyCorrection_"); | |
2656 | saveFileName.ReplaceAll(".root", Form("_jetPt%.1f_%.1f%s.root", actualLowerJetPt, actualUpperJetPt, chargeString.Data())); | |
2657 | ||
2658 | TString saveFilePathName = Form("%s/%s", pathData.Data(), saveFileName.Data()); | |
2659 | TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE"); | |
2660 | ||
2661 | if (!saveFile) { | |
2662 | printf("Could not save results!\n"); | |
2663 | ||
2664 | return -1; | |
2665 | } | |
2666 | ||
2667 | saveFile->cd(); | |
2668 | ||
2669 | if (cSecPt) | |
2670 | cSecPt->Write(); | |
2671 | ||
2672 | if (cSecSSPt) | |
2673 | cSecSSPt->Write(); | |
2674 | ||
2675 | if (hSecAllPt) | |
2676 | hSecAllPt->Write(); | |
2677 | ||
2678 | if (hSecSSAllPt) | |
2679 | hSecSSAllPt->Write(); | |
2680 | ||
2681 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2682 | if (hSecPt[i]) | |
2683 | hSecPt[i]->Write(); | |
2684 | ||
2685 | if (hSecSSPt[i]) | |
2686 | hSecSSPt[i]->Write(); | |
2687 | } | |
2688 | ||
2689 | if (cSecToPiRatioPt) | |
2690 | cSecToPiRatioPt->Write(); | |
2691 | ||
2692 | if (cSecSSToPiRatioPt) | |
2693 | cSecSSToPiRatioPt->Write(); | |
2694 | ||
2695 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2696 | if (hSecToPiRatioPt[i]) | |
2697 | hSecToPiRatioPt[i]->Write(); | |
2698 | ||
2699 | if (hSecSSToPiRatioPt[i]) | |
2700 | hSecSSToPiRatioPt[i]->Write(); | |
2701 | } | |
2702 | ||
2703 | ||
2704 | ||
2705 | if (cSecZ) | |
2706 | cSecZ->Write(); | |
2707 | ||
2708 | if (cSecSSZ) | |
2709 | cSecSSZ->Write(); | |
2710 | ||
2711 | if (hSecAllZ) | |
2712 | hSecAllZ->Write(); | |
2713 | ||
2714 | if (hSecSSAllZ) | |
2715 | hSecSSAllZ->Write(); | |
2716 | ||
2717 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2718 | if (hSecZ[i]) | |
2719 | hSecZ[i]->Write(); | |
2720 | ||
2721 | if (hSecSSZ[i]) | |
2722 | hSecSSZ[i]->Write(); | |
2723 | } | |
2724 | ||
2725 | if (cSecToPiRatioZ) | |
2726 | cSecToPiRatioZ->Write(); | |
2727 | ||
2728 | if (cSecSSToPiRatioZ) | |
2729 | cSecSSToPiRatioZ->Write(); | |
2730 | ||
2731 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2732 | if (hSecToPiRatioZ[i]) | |
2733 | hSecToPiRatioZ[i]->Write(); | |
2734 | ||
2735 | if (hSecSSToPiRatioZ[i]) | |
2736 | hSecSSToPiRatioZ[i]->Write(); | |
2737 | } | |
2738 | ||
2739 | ||
2740 | ||
2741 | if (cSecXi) | |
2742 | cSecXi->Write(); | |
2743 | ||
2744 | if (cSecSSXi) | |
2745 | cSecSSXi->Write(); | |
2746 | ||
2747 | if (hSecAllXi) | |
2748 | hSecAllXi->Write(); | |
2749 | ||
2750 | if (hSecSSAllXi) | |
2751 | hSecSSAllXi->Write(); | |
2752 | ||
2753 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2754 | if (hSecXi[i]) | |
2755 | hSecXi[i]->Write(); | |
2756 | ||
2757 | if (hSecSSXi[i]) | |
2758 | hSecSSXi[i]->Write(); | |
2759 | } | |
2760 | ||
2761 | if (cSecToPiRatioXi) | |
2762 | cSecToPiRatioXi->Write(); | |
2763 | ||
2764 | if (cSecSSToPiRatioXi) | |
2765 | cSecSSToPiRatioXi->Write(); | |
2766 | ||
2767 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2768 | if (hSecToPiRatioXi[i]) | |
2769 | hSecToPiRatioXi[i]->Write(); | |
2770 | ||
2771 | if (hSecSSToPiRatioXi[i]) | |
2772 | hSecSSToPiRatioXi[i]->Write(); | |
2773 | } | |
2774 | ||
2775 | ||
2776 | ||
2777 | if (cSpecPt) | |
2778 | cSpecPt->Write(); | |
2779 | ||
2780 | if (hSpecPt_RecPrimaries) | |
2781 | hSpecPt_RecPrimaries->Write(); | |
2782 | ||
2783 | if (hSpecPt_GenPrimaries) | |
2784 | hSpecPt_GenPrimaries->Write(); | |
2785 | ||
2786 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2787 | if (hSpecIDPt_GenPrimaries[i]) | |
2788 | hSpecIDPt_GenPrimaries[i]->Write(); | |
2789 | ||
2790 | if (hSpecIDPt_RecPrimaries[i]) | |
2791 | hSpecIDPt_RecPrimaries[i]->Write(); | |
2792 | } | |
2793 | ||
2794 | ||
2795 | if (cEffPt) | |
2796 | cEffPt->Write(); | |
2797 | ||
2798 | if (hEfficiencyAllPt) | |
2799 | hEfficiencyAllPt->Write(); | |
2800 | ||
2801 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2802 | if (hEfficiencyPt[i]) | |
2803 | hEfficiencyPt[i]->Write(); | |
2804 | } | |
2805 | ||
2806 | if (cEffToPiRatioPt) | |
2807 | cEffToPiRatioPt->Write(); | |
2808 | ||
2809 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2810 | if (hEfficiencyToPiRatioPt[i]) | |
2811 | hEfficiencyToPiRatioPt[i]->Write(); | |
2812 | } | |
2813 | ||
2814 | ||
2815 | ||
2816 | if (cSpecZ) | |
2817 | cSpecZ->Write(); | |
2818 | ||
2819 | if (hSpecZ_RecPrimaries) | |
2820 | hSpecZ_RecPrimaries->Write(); | |
2821 | ||
2822 | if (hSpecZ_GenPrimaries) | |
2823 | hSpecZ_GenPrimaries->Write(); | |
2824 | ||
2825 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2826 | if (hSpecIDZ_GenPrimaries[i]) | |
2827 | hSpecIDZ_GenPrimaries[i]->Write(); | |
2828 | ||
2829 | if (hSpecIDZ_RecPrimaries[i]) | |
2830 | hSpecIDZ_RecPrimaries[i]->Write(); | |
2831 | } | |
2832 | ||
2833 | ||
2834 | if (cEffZ) | |
2835 | cEffZ->Write(); | |
2836 | ||
2837 | if (hEfficiencyAllZ) | |
2838 | hEfficiencyAllZ->Write(); | |
2839 | ||
2840 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2841 | if (hEfficiencyZ[i]) | |
2842 | hEfficiencyZ[i]->Write(); | |
2843 | } | |
2844 | ||
2845 | if (cEffToPiRatioZ) | |
2846 | cEffToPiRatioZ->Write(); | |
2847 | ||
2848 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2849 | if (hEfficiencyToPiRatioZ[i]) | |
2850 | hEfficiencyToPiRatioZ[i]->Write(); | |
2851 | } | |
2852 | ||
2853 | ||
2854 | ||
2855 | if (cSpecXi) | |
2856 | cSpecXi->Write(); | |
2857 | ||
2858 | if (hSpecXi_RecPrimaries) | |
2859 | hSpecXi_RecPrimaries->Write(); | |
2860 | ||
2861 | if (hSpecXi_GenPrimaries) | |
2862 | hSpecXi_GenPrimaries->Write(); | |
2863 | ||
2864 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2865 | if (hSpecIDXi_GenPrimaries[i]) | |
2866 | hSpecIDXi_GenPrimaries[i]->Write(); | |
2867 | ||
2868 | if (hSpecIDXi_RecPrimaries[i]) | |
2869 | hSpecIDXi_RecPrimaries[i]->Write(); | |
2870 | } | |
2871 | ||
2872 | ||
2873 | if (cEffXi) | |
2874 | cEffXi->Write(); | |
2875 | ||
2876 | if (hEfficiencyAllXi) | |
2877 | hEfficiencyAllXi->Write(); | |
2878 | ||
2879 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2880 | if (hEfficiencyXi[i]) | |
2881 | hEfficiencyXi[i]->Write(); | |
2882 | } | |
2883 | ||
2884 | if (cEffToPiRatioXi) | |
2885 | cEffToPiRatioXi->Write(); | |
2886 | ||
2887 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
2888 | if (hEfficiencyToPiRatioXi[i]) | |
2889 | hEfficiencyToPiRatioXi[i]->Write(); | |
2890 | } | |
2891 | ||
2892 | ||
2893 | ||
2894 | // Correct the yields with the efficiencies and primary fractions | |
2895 | for (Int_t type = 0; type < kNtypes; type++) { | |
2896 | TH1D* hYieldCorrected[AliPID::kSPECIES]; | |
2897 | TH1D* hYieldCorrectedSysError[AliPID::kSPECIES]; | |
2898 | ||
2899 | TH1D** hYield; | |
2900 | TH1D** hYieldSysError; | |
2901 | TH1D** hMCgenPrimYield; | |
2902 | TH1D** hSecSS; | |
2903 | TH1D** hSec; | |
2904 | TH1D** hEfficiency; | |
2905 | ||
2906 | if (type == kTrackPt) { | |
2907 | hYield = &hYieldPt[0]; | |
2908 | hYieldSysError = &hYieldPtSysError[0]; | |
2909 | hMCgenPrimYield = &hMCgenPrimYieldPt[0]; | |
2910 | hSecSS = &hSecSSPt[0]; | |
2911 | hSec = &hSecPt[0]; | |
2912 | hEfficiency = &hEfficiencyPt[0]; | |
2913 | } | |
2914 | else if (type == kZ) { | |
2915 | hYield = &hYieldZ[0]; | |
2916 | hYieldSysError = &hYieldZSysError[0]; | |
2917 | hMCgenPrimYield = &hMCgenPrimYieldZ[0]; | |
2918 | hSecSS = &hSecSSZ[0]; | |
2919 | hSec = &hSecZ[0]; | |
2920 | hEfficiency = &hEfficiencyZ[0]; | |
2921 | } | |
2922 | else if (type == kXi) { | |
2923 | hYield = &hYieldXi[0]; | |
2924 | hYieldSysError = &hYieldXiSysError[0]; | |
2925 | hMCgenPrimYield = &hMCgenPrimYieldXi[0]; | |
2926 | hSecSS = &hSecSSXi[0]; | |
2927 | hSec = &hSecXi[0]; | |
2928 | hEfficiency = &hEfficiencyXi[0]; | |
2929 | } | |
2930 | else | |
2931 | continue; | |
2932 | ||
2933 | ||
2934 | // NOTE: The yields have already been normalised to nJets in extractFFs. However, the binning in jetPt might be different now. | |
2935 | // Thus, in order to get the correct result, the old normalisation (just the same binnnig as in the nJet histo) has to be | |
2936 | // undone (was done above) and then it has to be normalised to the number of jets in the new binning again (now). | |
2937 | // Also NOTE: Histos are already normalised to bin width via extractFFs! | |
2938 | const Double_t nJetsGenData = hNjetsGenData ? hNjetsGenData->Integral(lowerCentralityBinLimitData, upperCentralityBinLimitData, | |
2939 | lowerJetPtBinLimitData, upperJetPtBinLimitData) : 1.; | |
2940 | const Double_t nJetsRecData = hNjetsRecData ? hNjetsRecData->Integral(lowerCentralityBinLimitData, upperCentralityBinLimitData, | |
2941 | lowerJetPtBinLimitData, upperJetPtBinLimitData) : 1.; | |
2942 | ||
2943 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
2944 | hYield[species]->Scale((nJetsRecData > 0) ? 1. / nJetsRecData : 0.); | |
2945 | hYieldSysError[species]->Scale((nJetsRecData > 0) ? 1. / nJetsRecData : 0.); | |
2946 | ||
2947 | hYield[species]->SetStats(kFALSE); | |
2948 | hYieldSysError[species]->SetStats(kFALSE); | |
2949 | ||
2950 | if (hMCgenPrimYield[species]) { | |
2951 | hMCgenPrimYield[species]->Scale((nJetsGenData > 0) ? 1. / nJetsGenData : 0.); | |
2952 | hMCgenPrimYield[species]->SetStats(kFALSE); | |
2953 | } | |
2954 | ||
2955 | hYield[species]->GetYaxis()->SetTitle(titleYaxis[type].Data()); | |
2956 | if (hMCgenPrimYield[species]) { | |
2957 | hMCgenPrimYield[species]->GetYaxis()->SetTitle(titleYaxis[type].Data()); | |
2958 | hMCgenPrimYield[species]->SetTitle(Form("%s (MC)", AliPID::ParticleLatexName(species))); | |
2959 | hMCgenPrimYield[species]->SetMarkerStyle(24); | |
2960 | hMCgenPrimYield[species]->SetLineColor(getLineColorAliPID(species)); | |
2961 | hMCgenPrimYield[species]->SetMarkerColor(getLineColorAliPID(species)); | |
2962 | } | |
2963 | ||
2964 | hYield[species]->SetLineColor(getLineColorAliPID(species)); | |
2965 | hYield[species]->SetMarkerColor(getLineColorAliPID(species)); | |
2966 | hYield[species]->SetMarkerStyle(20); | |
2967 | ||
2968 | hYieldSysError[species]->SetLineColor(getLineColorAliPID(species)); | |
2969 | hYieldSysError[species]->SetMarkerColor(getLineColorAliPID(species)); | |
2970 | hYieldSysError[species]->SetMarkerStyle(20); | |
2971 | hYieldSysError[species]->SetFillStyle(0); | |
2972 | ||
2973 | hYieldCorrected[species] = new TH1D(*hYield[species]); | |
2974 | hYieldCorrected[species]->SetName(Form("%s_corrected", hYield[species]->GetName())); | |
2975 | hYieldCorrected[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
2976 | ||
2977 | hYieldCorrectedSysError[species] = new TH1D(*hYieldSysError[species]); | |
2978 | hYieldCorrectedSysError[species]->SetName(Form("%s_corrected", hYieldSysError[species]->GetName())); | |
2979 | hYieldCorrectedSysError[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
2980 | ||
2981 | // Correction factor histos can have different binning than data histos (constant factor above some threshold) | |
2982 | // -> Need special functions to multiply and divide such histos | |
2983 | multiplyHistsDifferentBinning(hYieldCorrected[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.); | |
2984 | divideHistsDifferentBinning(hYieldCorrected[species], hEfficiency[species], 1., 1.); | |
2985 | ||
2986 | multiplyHistsDifferentBinning(hYieldCorrectedSysError[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1., | |
2987 | kTRUE); | |
2988 | divideHistsDifferentBinning(hYieldCorrectedSysError[species], hEfficiency[species], 1., 1., kTRUE); | |
2989 | ||
2990 | //hYieldCorrected[species]->Multiply(hYieldCorrected[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1., | |
2991 | //""); | |
2992 | //hYieldCorrected[species]->Divide(hYieldCorrected[species], hEfficiency[species], 1., 1., ""); | |
2993 | } | |
2994 | ||
2995 | // Calculate the total corrected yield. The original total yield had no error (just the number of detected tracks in a pT bin), | |
2996 | // but due to the correction there is some error for the total yield. Also the error of the fractions introduces uncertainties | |
2997 | // for the yields of individual species | |
2998 | TH1D* hYieldCorrectedTotal = new TH1D(*hYieldCorrected[0]); | |
2999 | hYieldCorrectedTotal->SetLineColor(kBlack); | |
3000 | hYieldCorrectedTotal->SetMarkerColor(kBlack); | |
3001 | hYieldCorrectedTotal->SetMarkerStyle(20); | |
3002 | hYieldCorrectedTotal->SetName(Form("hYieldCorrectedTotal_%s", titles[type].Data())); | |
3003 | hYieldCorrectedTotal->SetTitle("Total"); | |
3004 | //hYieldCorrectedTotal->SetTitle("Total, secondary and efficiency x acceptance x pT resolution corrected"); | |
3005 | ||
3006 | //TODO doesn't it require the covariance matrix here to get the error of the corrected total yield? | |
3007 | // Because before the total count was taken into account just by the fit. Now the fractions and their errors | |
3008 | // influence the MC correction (and its errors), so that there is some error on the corrected total yield. | |
3009 | // However, the errors of the contributing yields are correlated! | |
3010 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) | |
3011 | hYieldCorrectedTotal->Add(hYieldCorrected[i], 1.); | |
3012 | ||
3013 | // Calculate the corrected fractions | |
3014 | TH1D* hFractionCorrected[AliPID::kSPECIES]; | |
3015 | TH1D* hFractionCorrectedSysError[AliPID::kSPECIES]; | |
3016 | ||
3017 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
3018 | hFractionCorrected[species] = new TH1D(*hYield[species]); | |
3019 | TString oldName = hYield[species]->GetName(); | |
3020 | TString newName = oldName.ReplaceAll("Yield", "Fraction"); | |
3021 | hFractionCorrected[species]->SetName(Form("%s_corrected", newName.Data())); | |
3022 | hFractionCorrected[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
3023 | hFractionCorrected[species]->GetYaxis()->SetTitle("Corrected Fraction"); | |
3024 | hFractionCorrected[species]->GetXaxis()->SetMoreLogLabels(kTRUE); | |
3025 | hFractionCorrected[species]->GetXaxis()->SetNoExponent(kTRUE); | |
3026 | hFractionCorrected[species]->SetFillStyle(0); | |
3027 | ||
3028 | ||
3029 | // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation | |
3030 | // (numerator is a "selected" subset of the denominator). It doesn't matter that the error of the histos is not just | |
3031 | // "sqrt(content)" because the error formula also works for weighted histograms (which means that the error can be more | |
3032 | // or less anything) | |
3033 | hFractionCorrected[species]->Divide(hYieldCorrected[species], hYieldCorrectedTotal, 1., 1., "B"); | |
3034 | ||
3035 | ||
3036 | ||
3037 | // The systematic errors just need to be scaled in the same way as the fractions. | |
3038 | // So, just take the ratio to the uncorrected fraction and scale the sys. error accordingly | |
3039 | // or, in this case, just divide by the same total yield as for yield -> fractions | |
3040 | hFractionCorrectedSysError[species] = new TH1D(*hFractionCorrected[species]); | |
3041 | hFractionCorrectedSysError[species]->SetName(Form("%s_sysError", hFractionCorrected[species]->GetName())); | |
3042 | hFractionCorrectedSysError[species]->SetTitle(Form("%s (sys. error)", hFractionCorrected[species]->GetTitle())); | |
3043 | ||
3044 | for (Int_t binX = 1; binX <= hFractionCorrectedSysError[species]->GetNbinsX(); binX++) { | |
3045 | const Double_t corrTotalYield = hYieldCorrectedTotal->GetBinContent(binX); | |
3046 | const Double_t scaleFactor = corrTotalYield > 0 ? 1.0 / corrTotalYield : 1.; | |
3047 | hFractionCorrectedSysError[species]->SetBinError(binX, hYieldCorrectedSysError[species]->GetBinError(binX) * scaleFactor); | |
3048 | } | |
3049 | } | |
3050 | ||
3051 | // If MC is available, calculate the generated fractions | |
3052 | TH1D* hMCgenPrimYieldTotal = 0x0; | |
3053 | TH1D* hMCgenPrimFraction[AliPID::kSPECIES]; | |
3054 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) | |
3055 | hMCgenPrimFraction[i] = 0x0; | |
3056 | ||
3057 | if (numMCgenPrimYieldHistsFound > 0) { | |
3058 | hMCgenPrimYieldTotal = new TH1D(*hMCgenPrimYield[0]); | |
3059 | hMCgenPrimYieldTotal->SetLineColor(kBlack); | |
3060 | hMCgenPrimYieldTotal->SetMarkerColor(kBlack); | |
3061 | hMCgenPrimYieldTotal->SetMarkerStyle(24); | |
3062 | hMCgenPrimYieldTotal->SetName(Form("hMCgenPrimYieldTotal%s", titles[type].Data())); | |
3063 | hMCgenPrimYieldTotal->SetTitle("Total, MC truth"); | |
3064 | //hMCgenPrimYieldTotal->SetTitle("Total generated primary yield (MC truth)"); | |
3065 | ||
3066 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) | |
3067 | hMCgenPrimYieldTotal->Add(hMCgenPrimYield[i], 1.); | |
3068 | ||
3069 | // Calculate the MC fractions | |
3070 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
3071 | hMCgenPrimYield[species]->SetLineColor(getLineColorAliPID(species)); | |
3072 | hMCgenPrimYield[species]->SetMarkerColor(getLineColorAliPID(species)); | |
3073 | hMCgenPrimYield[species]->SetMarkerStyle(24); | |
3074 | ||
3075 | hMCgenPrimYield[species]->SetTitle(Form("%s, MC truth", AliPID::ParticleLatexName(species))); | |
3076 | ||
3077 | hMCgenPrimFraction[species] = new TH1D(*hMCgenPrimYield[species]); | |
3078 | TString oldName = hMCgenPrimYield[species]->GetName(); | |
3079 | TString newName = oldName.ReplaceAll("Yield", "Fraction"); | |
3080 | hMCgenPrimFraction[species]->SetName(newName.Data()); | |
3081 | ||
3082 | // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation | |
3083 | // (numerator is a "selected" subset of the denominator). | |
3084 | hMCgenPrimFraction[species]->Divide(hMCgenPrimFraction[species], hMCgenPrimYieldTotal, 1., 1., "B"); | |
3085 | } | |
3086 | } | |
3087 | ||
3088 | TCanvas* cCorrData = new TCanvas(Form("cCorrData_%s", titles[type].Data()), "Corrected data", 0, 300, 900, 900); | |
3089 | cCorrData->Divide(2, 1, 0., 0.01); | |
3090 | cCorrData->GetPad(1)->SetLogx(1); | |
3091 | cCorrData->GetPad(1)->SetLogy(1); | |
3092 | cCorrData->GetPad(2)->SetLogx(1); | |
3093 | cCorrData->GetPad(2)->SetLogy(1); | |
3094 | ||
3095 | cCorrData->GetPad(1)->SetRightMargin(0.001); | |
3096 | cCorrData->GetPad(2)->SetRightMargin(0.001); | |
3097 | ||
3098 | cCorrData->GetPad(1)->SetLeftMargin(0.2); | |
3099 | cCorrData->GetPad(2)->SetLeftMargin(0.2); | |
3100 | ||
3101 | cCorrData->cd(1); // uncorrected | |
3102 | hYield[AliPID::kPion]->GetYaxis()->SetTitleOffset(1.4); | |
3103 | hYield[AliPID::kPion]->Draw("E1"); | |
3104 | ||
3105 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3106 | hYieldSysError[i]->Draw("E2 same"); | |
3107 | ||
3108 | if (i == AliPID::kPion) continue; | |
3109 | hYield[i]->Draw("E1 same"); | |
3110 | } | |
3111 | ||
3112 | ||
3113 | ClearTitleFromHistoInCanvas(cCorrData, 1); | |
3114 | ||
3115 | Double_t maximum = hYieldCorrectedTotal->GetBinContent(hYieldCorrectedTotal->GetMaximumBin()) * 10.; | |
3116 | Double_t minimum = maximum; | |
3117 | for (Int_t spec = 0; spec < AliPID::kSPECIES; spec++) { | |
3118 | Double_t temp = hYieldCorrected[spec]->GetBinContent(hYieldCorrected[spec]->FindLastBinAbove(0.)) * 0.1; | |
3119 | if (temp > 0) | |
3120 | minimum = temp; | |
3121 | } | |
3122 | ||
3123 | ||
3124 | cCorrData->cd(2); // corrected | |
3125 | hYieldCorrectedTotal->GetXaxis()->SetMoreLogLabels(kTRUE); | |
3126 | hYieldCorrectedTotal->GetXaxis()->SetNoExponent(kTRUE); | |
3127 | hYieldCorrectedTotal->GetYaxis()->SetRangeUser(minimum, maximum); | |
3128 | ||
3129 | if (hMCgenPrimYieldTotal) { | |
3130 | hMCgenPrimYieldTotal->GetYaxis()->SetTitleOffset(1.4); | |
3131 | hMCgenPrimYieldTotal->GetXaxis()->SetMoreLogLabels(kTRUE); | |
3132 | hMCgenPrimYieldTotal->GetXaxis()->SetNoExponent(kTRUE); | |
3133 | hMCgenPrimYieldTotal->GetXaxis()->SetTitle(hYieldCorrectedTotal->GetXaxis()->GetTitle()); | |
3134 | hMCgenPrimYieldTotal->GetYaxis()->SetRangeUser(minimum, maximum); | |
3135 | } | |
3136 | ||
3137 | if (hMCgenPrimYieldTotal) { | |
3138 | hMCgenPrimYieldTotal->Draw("E1"); | |
3139 | hYieldCorrectedTotal->Draw("E1 same"); | |
3140 | } | |
3141 | else | |
3142 | hYieldCorrectedTotal->Draw("E1"); | |
3143 | ||
3144 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3145 | if (hMCgenPrimYield[i]) | |
3146 | hMCgenPrimYield[i]->Draw("E1 same"); | |
3147 | hYieldCorrected[i]->Draw("E1 same"); | |
3148 | } | |
3149 | ||
3150 | TLegend* legTemp = cCorrData->cd(2)->BuildLegend();//0.25, 0.16, 0.65, 0.51); | |
3151 | legTemp->SetNColumns(2); | |
3152 | legTemp->SetFillColor(kWhite); | |
3153 | ||
3154 | // Do not include into legend | |
3155 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) | |
3156 | hYieldCorrectedSysError[i]->Draw("E2 same"); | |
3157 | ||
3158 | ClearTitleFromHistoInCanvas(cCorrData, 2); | |
3159 | ||
3160 | TCanvas* cCorrYieldsRatio = 0x0; | |
3161 | ||
3162 | TH1D* hYieldCorrectedRatioToMC[AliPID::kSPECIES]; | |
3163 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) | |
3164 | hYieldCorrectedRatioToMC[i] = 0x0; | |
3165 | ||
3166 | TH1D* hYieldCorrectedTotalRatioToMC = 0x0; | |
3167 | ||
3168 | if (numMCgenPrimYieldHistsFound > 0) { | |
3169 | // Compare with MC truth | |
3170 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) { | |
3171 | hYieldCorrectedRatioToMC[species] = new TH1D(*hYieldCorrected[species]); | |
3172 | hYieldCorrectedRatioToMC[species]->SetName(Form("hYieldCorrectedRatioToMC_%s_%s", titles[type].Data(), | |
3173 | AliPID::ParticleShortName(species))); | |
3174 | hYieldCorrectedRatioToMC[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species))); | |
3175 | hYieldCorrectedRatioToMC[species]->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth"); | |
3176 | hYieldCorrectedRatioToMC[species]->Divide(hMCgenPrimYield[species]); | |
3177 | } | |
3178 | ||
3179 | hYieldCorrectedTotalRatioToMC = new TH1D(*hYieldCorrectedTotal); | |
3180 | hYieldCorrectedTotalRatioToMC->SetName(Form("hYieldCorrectedTotalRatioToMC_%s", titles[type].Data())); | |
3181 | hYieldCorrectedTotalRatioToMC->SetTitle("Total yield"); | |
3182 | hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth"); | |
3183 | hYieldCorrectedTotalRatioToMC->Divide(hMCgenPrimYieldTotal); | |
3184 | ||
3185 | cCorrYieldsRatio = new TCanvas(Form("cCorrYieldsRatio_%s", titles[type].Data()), "Corrected Yields Comparison to MC", | |
3186 | 0, 300, 900, 900); | |
3187 | cCorrYieldsRatio->SetGridx(1); | |
3188 | cCorrYieldsRatio->SetGridy(1); | |
3189 | cCorrYieldsRatio->SetLogx(1); | |
3190 | ||
3191 | hYieldCorrectedTotalRatioToMC->GetYaxis()->SetRangeUser(0.6, 1.6); | |
3192 | hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitleOffset(0.85); | |
3193 | hYieldCorrectedTotalRatioToMC->Draw("E1"); | |
3194 | ||
3195 | for (Int_t species = 0; species < AliPID::kSPECIES; species++) | |
3196 | hYieldCorrectedRatioToMC[species]->Draw("E1 same"); | |
3197 | ||
3198 | cCorrYieldsRatio->BuildLegend()->SetFillColor(kWhite); | |
3199 | ||
3200 | ClearTitleFromHistoInCanvas(cCorrYieldsRatio); | |
3201 | } | |
3202 | ||
3203 | TCanvas* cCorrFractions = new TCanvas(Form("cCorrFractions_%s", titles[type].Data()), "Corrected particleFractions", | |
3204 | 0, 300, 900, 900); | |
3205 | cCorrFractions->SetLogx(1); | |
3206 | hFractionCorrected[0]->GetYaxis()->SetRangeUser(0., 1.); | |
3207 | hFractionCorrected[0]->Draw("E1"); | |
3208 | if (hMCgenPrimFraction[0]) | |
3209 | hMCgenPrimFraction[0]->Draw("E1 same"); | |
3210 | ||
3211 | for (Int_t i = 1; i < AliPID::kSPECIES; i++) { | |
3212 | hFractionCorrected[i]->Draw("E1 same"); | |
3213 | if (hMCgenPrimFraction[i]) | |
3214 | hMCgenPrimFraction[i]->Draw("E1 same"); | |
3215 | } | |
3216 | ||
3217 | cCorrFractions->BuildLegend()->SetFillColor(kWhite); | |
3218 | ||
3219 | // Do not include into legend | |
3220 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) | |
3221 | hFractionCorrectedSysError[i]->Draw("E2 same"); | |
3222 | ||
3223 | ClearTitleFromHistoInCanvas(cCorrFractions); | |
3224 | ||
3225 | // Save results | |
3226 | saveFile->cd(); | |
3227 | ||
3228 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3229 | if (hYield[i]) | |
3230 | hYield[i]->Write(); | |
3231 | } | |
3232 | ||
3233 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3234 | if (hYieldCorrected[i]) | |
3235 | hYieldCorrected[i]->Write(); | |
3236 | ||
3237 | if (hYieldCorrectedSysError[i]) | |
3238 | hYieldCorrectedSysError[i]->Write(); | |
3239 | } | |
3240 | ||
3241 | if (hYieldCorrectedTotal) | |
3242 | hYieldCorrectedTotal->Write(); | |
3243 | ||
3244 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3245 | if (hMCgenPrimYield[i]) | |
3246 | hMCgenPrimYield[i]->Write(); | |
3247 | } | |
3248 | ||
3249 | if (hMCgenPrimYieldTotal) | |
3250 | hMCgenPrimYieldTotal->Write(); | |
3251 | ||
3252 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3253 | if (hYieldCorrectedRatioToMC[i]) | |
3254 | hYieldCorrectedRatioToMC[i]->Write(); | |
3255 | } | |
3256 | ||
3257 | if (hYieldCorrectedTotalRatioToMC) | |
3258 | hYieldCorrectedTotalRatioToMC->Write(); | |
3259 | ||
3260 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3261 | if (hFractionCorrected[i]) | |
3262 | hFractionCorrected[i]->Write(); | |
3263 | ||
3264 | if (hFractionCorrectedSysError[i]) | |
3265 | hFractionCorrectedSysError[i]->Write(); | |
3266 | } | |
3267 | ||
3268 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
3269 | if (hMCgenPrimFraction[i]) | |
3270 | hMCgenPrimFraction[i]->Write(); | |
3271 | } | |
3272 | ||
3273 | if (cCorrData) | |
3274 | cCorrData->Write(); | |
3275 | ||
3276 | if (cCorrYieldsRatio) | |
3277 | cCorrYieldsRatio->Write(); | |
3278 | ||
3279 | if (cCorrFractions) | |
3280 | cCorrFractions->Write(); | |
3281 | } | |
3282 | ||
3283 | saveFile->cd(); | |
3284 | ||
3285 | TNamed* settings = new TNamed( | |
3286 | Form("Settings: Efficiency file \"%s\", numJetsForEfficiency file \"%s\", Data file \"%s\", lowerCentrality %.3f, upperCentrality %.3f, lowerJetPt %.1f, upperJetPt %.1f, constantCorrectionAbovePtThreshold %.3f\n", | |
3287 | pathNameEfficiency.Data(), pathNameDataMC.Data(), pathNameData.Data(), lowerCentrality, upperCentrality, lowerJetPt, | |
3288 | upperJetPt, constantCorrectionAbovePtThreshold), ""); | |
3289 | settings->Write(); | |
3290 | ||
3291 | saveFile->Close(); | |
3292 | ||
3293 | ||
3294 | return 0; | |
3295 | } |