]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/macros/PID/extractFFs.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / macros / PID / extractFFs.C
CommitLineData
e131b05f 1#include "TH1D.h"
2#include "TH2D.h"
3#include "TCanvas.h"
4#include "TFile.h"
5#include "TLegend.h"
6#include "TROOT.h"
7#include "TSystem.h"
8
9#include "AliPID.h"
10
11#include <iostream>
12
13#include "THnSparseDefinitions.h"
14
15#include "../../UserTasks/AliAnalysisTaskPID.h"
16#include "SystematicErrorUtils.h"
17
18enum axisDataProj { kPidPtProj = 0, kPidJetPtProj = 1, kPidZProj = 2, kPidXiProj = 3, kNDimsProj = 4 };
19
20
21
22//___________________________________________________________________
23void setupHist(TH1* h, TString histName, TString histTitle, TString xAxisTitle, TString yAxisTitle, Int_t color)
24{
25 if (histName != "")
26 h->SetName(histName.Data());
27 h->SetTitle(histTitle.Data());
28
29 if (xAxisTitle != "")
30 h->GetXaxis()->SetTitle(xAxisTitle.Data());
31 if (yAxisTitle != "")
32 h->GetYaxis()->SetTitle(yAxisTitle.Data());
33
34 h->SetMarkerStyle(24);
35 h->SetLineColor(color);
36 h->SetMarkerColor(color);
37
38 h->SetStats(kFALSE);
39}
40
41
42//____________________________________________________________________________________________________________________
43void normaliseYieldHist2D(TH2* hData, TH2* hNumJets, const Int_t lowerCentralityBinLimit, const Int_t upperCentralityBinLimit)
44{
45 // Normalise to 1/numJets and bin width. NOTE: jetPt binning of hData and hNumJets assumed to be the same!
46
47 if (!hData)
48 return;
49
50 for (Int_t binJetPt = 0; binJetPt <= hData->GetNbinsY() + 1; binJetPt++) {
51 // No normalisation to numJets, if no histo is provided
52 const Double_t numJets = hNumJets ? hNumJets->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, binJetPt, binJetPt) : 1.;
53 Bool_t noJets = numJets < 1e-13;
54
55 for (Int_t binObs = 0; binObs <= hData->GetNbinsX() + 1; binObs++) {
56 if (noJets) {
57 if (hData->GetBinContent(binObs, binJetPt) > 0.) {
58 printf("Error: No jets for jetPt ~ %f, but found content %f at y-coord %f!\n", hData->GetYaxis()->GetBinCenter(binJetPt),
59 hData->GetBinContent(binObs, binJetPt), hData->GetXaxis()->GetBinCenter(binObs));
60 }
61 continue;
62 }
63 const Double_t dObservable = hData->GetXaxis()->GetBinWidth(binObs);
64 const Double_t normFactor = 1. / (numJets * dObservable);
65
66 hData->SetBinContent(binObs, binJetPt, hData->GetBinContent(binObs, binJetPt) * normFactor);
67 hData->SetBinError(binObs, binJetPt, hData->GetBinError(binObs, binJetPt) * normFactor);
68 }
69 }
70}
71
72//____________________________________________________________________________________________________________________
73void normaliseYieldHist(TH1* h, Double_t numJets)
74{
75 // Normalise to 1/numJets and bin width
76
77 if (numJets <= 0) // Do not normalise
78 numJets = 1.;
79
80 for (Int_t bin = 0; bin <= h->GetNbinsX() + 1; bin++) {
81 const Double_t dObservable = h->GetBinWidth(bin);
82 const Double_t normFactor = 1. / (numJets * dObservable);
83 h->SetBinContent(bin, h->GetBinContent(bin) * normFactor);
84 h->SetBinError(bin, h->GetBinError(bin) * normFactor);
85 }
86}
87
88
89//___________________________________________________________________
90void GenerateParticleYields(THnSparse* hDataProj, AliAnalysisTaskPID* pidTask, const Double_t centrality,
91 TH2D** hIDFFtrackPt, TH2D** hIDFFz, TH2D** hIDFFxi,
92 const Bool_t setMean, const Bool_t addErrorsQuadratically, const Bool_t smearByError,
93 const Bool_t uniformSystematicError, const Bool_t takeIntoAccountSysError, Int_t nGenerations)
94{
95 if (!hDataProj || !pidTask || !hIDFFtrackPt || !hIDFFz || !hIDFFxi) {
96 printf("Cannot generate particle fractions - missing input!\n");
97 return;
98 }
99
100 if (takeIntoAccountSysError && smearByError) {
101 printf("It doesn't make sense to correlate statistical and systematic errors!\n");
102 return;
103 }
104
105 const Int_t nDimProj = hDataProj->GetNdimensions();
106 const Long64_t nBinsTHnSparseProj = hDataProj->GetNbins();
107 Double_t binContent = 0., binError2 = 0.;
108 Int_t binCoord[nDimProj];
109 Double_t binCentre[nDimProj];
110 Double_t prob[AliPID::kSPECIES];
111
112 Bool_t success = kTRUE;
113
114 // NOTE: In the following, the bin error is divided according to the fraction. The idea is to divide bin content (and error)
115 // into independent sub-samples according to the fractions. But then, adding up the errors of the sub-samples quadratically
116 // should recover the original error (think of 100 tracks, error sqrt(100), divided into 10 samples a 10 tracks => error of
117 // samples should be sqrt(10) then).
118 // This means, that one needs to add error^2 * fraction and NOT error^2 * fraction^2!!
119 // However, the errors are ignored anyway at the moment....
120
121 if (!smearByError && !takeIntoAccountSysError) {
122 nGenerations = 1;
123 // Just calculate fractions/spectra w/o any smearing
124 const Int_t smearSpeciesByError = -1;
125 const Int_t takeIntoAccountSysErrorOfSpecies = -1;
126
127 for (Long64_t bin = 0; bin < nBinsTHnSparseProj; bin++) {
128 binContent = hDataProj->GetBinContent(bin, binCoord);
129 binError2 = hDataProj->GetBinError2(bin);
130
131 // If the bin is empty, do not compute the particle fractions -> This bin contributes nothing anyway
132 if (binContent < 2. * std::numeric_limits<double>::min())
133 continue;
134
135 for (Int_t dim = 0; dim < nDimProj; dim++)
136 binCentre[dim] = hDataProj->GetAxis(dim)->GetBinCenter(binCoord[dim]);
137
138 // Since there is no smearing, the fraction will stay the same and there is no need to save the result for one iteration in a
139 // histogram (cf. case smearByError >= 0).
140 success = pidTask->GetParticleFractions(binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, prob,
141 smearSpeciesByError, takeIntoAccountSysErrorOfSpecies, uniformSystematicError);
142
143 if (success) {
144 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
145 // Track pT
146 hIDFFtrackPt[species]->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
147 hIDFFtrackPt[species]->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj])
148 + binContent * prob[species]);
149 Double_t tempError = hIDFFtrackPt[species]->GetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
150 hIDFFtrackPt[species]->SetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
151 TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
152
153 // z
154 hIDFFz[species]->SetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj],
155 hIDFFz[species]->GetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj])
156 + binContent * prob[species]);
157 tempError = hIDFFz[species]->GetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj]);
158 hIDFFz[species]->SetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj],
159 TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
160
161 // xi
162 hIDFFxi[species]->SetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
163 hIDFFxi[species]->GetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj])
164 + binContent * prob[species]);
165 tempError = hIDFFxi[species]->GetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj]);
166 hIDFFxi[species]->SetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
167 TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
168 }
169 }
170 else
171 printf("Problem obtaining fractions for: trackPt %f, jetPt %f, cent %f\n", binCentre[kPidPtProj], binCentre[kPidJetPtProj],
172 centrality);
173 }
174 }
175 else {
176 // Calculate fractions/spectra w/ smearing "nGenerations" times for every species and obtain the statistical or systematic error
177 // from the spread of the results.
178 // Note: For a given species, only the spread of the variation of exactly that species is considered, i.e. one is not
179 // looking at the change of the fraction if another species is varied (this would bias the result towards smaller errors
180 // by construction since the fraction of the considered species changes weighted with its statistical/systematic error)
181
182 TH2D* hIDFFtrackPtVaried[AliPID::kSPECIES][nGenerations];
183 TH2D* hIDFFzVaried[AliPID::kSPECIES][nGenerations];
184 TH2D* hIDFFxiVaried[AliPID::kSPECIES][nGenerations];
185
186 // Create this histo and set all entries to -1 (= not calculated yet) in every loop
187 TH2D* hPartFractionsSmeared = hDataProj->Projection(kPidJetPtProj, kPidPtProj);
188 hPartFractionsSmeared->SetName("hPartFractionsSmeared");
189
190 // Generate results with varying fractions
191 for (Int_t smearSpeciesByError = 0; smearSpeciesByError < AliPID::kSPECIES; smearSpeciesByError++) {
192 for (Int_t i = 0; i < nGenerations; i++) {
193 hIDFFtrackPtVaried[smearSpeciesByError][i] = new TH2D(*(TH2D*)hIDFFtrackPt[smearSpeciesByError]);
194 hIDFFtrackPtVaried[smearSpeciesByError][i]->SetName(Form("hIDFFtrackPtVaried_%d_%d", smearSpeciesByError, i));
195 hIDFFtrackPtVaried[smearSpeciesByError][i]->Reset();
196
197 hIDFFzVaried[smearSpeciesByError][i] = new TH2D(*(TH2D*)hIDFFz[smearSpeciesByError]);
198 hIDFFzVaried[smearSpeciesByError][i]->SetName(Form("hIDFFzVaried_%d_%d", smearSpeciesByError, i));
199 hIDFFzVaried[smearSpeciesByError][i]->Reset();
200
201 hIDFFxiVaried[smearSpeciesByError][i] = new TH2D(*(TH2D*)hIDFFxi[smearSpeciesByError]);
202 hIDFFxiVaried[smearSpeciesByError][i]->SetName(Form("hIDFFxiVaried_%d_%d", smearSpeciesByError, i));
203 hIDFFxiVaried[smearSpeciesByError][i]->Reset();
204
205 // In each iteration (moving over all bins), one want only ONE fixed particle fraction per bin in the spectra map.
206 // Otherwise, one would assing different fractions to e.g. the same trackPt bin, if only z and xi is different
207 // (but jetPt bin is still the same). This would then cause the fluctuations to cancel partially.
208 // Thus: Calculate the smeared fraction for each bin in the spectra map only once and store it in a histo.
209 // If it is requested again (i.e. histoEntry >= 0), use the value from the histo.
210
211 // Set all entries to -1 (= not calculated yet)
212 for (Int_t iX = 0; iX <= hPartFractionsSmeared->GetNbinsX(); iX++) {
213 for (Int_t iY = 0; iY <= hPartFractionsSmeared->GetNbinsY(); iY++) {
214 hPartFractionsSmeared->SetBinContent(iX, iY, -1);
215 }
216 }
217
218 for (Long64_t bin = 0; bin < nBinsTHnSparseProj; bin++) {
219 binContent = hDataProj->GetBinContent(bin, binCoord);
220 binError2 = hDataProj->GetBinError2(bin);
221
222 // If the bin is empty, do not compute the particle fractions -> This bin contributes nothing anyway
223 if (binContent < 2. * std::numeric_limits<double>::min())
224 continue;
225
226 for (Int_t iS = 0; iS < AliPID::kSPECIES; iS++)
227 prob[iS] = 0.;
228
229 for (Int_t dim = 0; dim < nDimProj; dim++)
230 binCentre[dim] = hDataProj->GetAxis(dim)->GetBinCenter(binCoord[dim]);
231
232 if (hPartFractionsSmeared->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj]) >= 0) {
233 success = kTRUE;
234 prob[smearSpeciesByError] = hPartFractionsSmeared->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
235 }
236 else {
237 const Int_t smearSpeciesByStatisticalError = smearByError ? smearSpeciesByError : -1;
238 const Int_t takeIntoAccountSysErrorOfSpecies = takeIntoAccountSysError ? smearSpeciesByError : -1;
239
240 success = pidTask->GetParticleFractions(binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, prob,
241 smearSpeciesByStatisticalError, takeIntoAccountSysErrorOfSpecies,
242 uniformSystematicError);
243
244 if (success)
245 hPartFractionsSmeared->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj], prob[smearSpeciesByError]);
246 }
247
248 // NOTE: Since only the bin content of the current species is stored in hPartFractionsSmeared,
249 // only prob[smearSpeciesByError] can be used in the following!!!!
250
251 if (success) {
252 // To make things readable...
253 TH2D* hIDFFtrackPtVariedCurr = hIDFFtrackPtVaried[smearSpeciesByError][i];
254 TH2D* hIDFFzVariedCurr = hIDFFzVaried[smearSpeciesByError][i];
255 TH2D* hIDFFxiVariedCurr = hIDFFxiVaried[smearSpeciesByError][i];
256
257 // Track pT
258 hIDFFtrackPtVariedCurr->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
259 hIDFFtrackPtVariedCurr->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj])
260 + binContent * prob[smearSpeciesByError]);
261 Double_t tempError = hIDFFtrackPtVariedCurr->GetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
262 hIDFFtrackPtVariedCurr->SetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
263 TMath::Sqrt(tempError * tempError + binError2 * prob[smearSpeciesByError]));
264
265 // z
266 hIDFFzVariedCurr->SetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj],
267 hIDFFzVariedCurr->GetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj])
268 + binContent * prob[smearSpeciesByError]);
269 tempError = hIDFFzVariedCurr->GetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj]);
270 hIDFFzVariedCurr->SetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj],
271 TMath::Sqrt(tempError * tempError + binError2 * prob[smearSpeciesByError]));
272
273 // xi
274 hIDFFxiVariedCurr->SetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
275 hIDFFxiVariedCurr->GetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj])
276 + binContent * prob[smearSpeciesByError]);
277 tempError = hIDFFxiVariedCurr->GetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj]);
278 hIDFFxiVariedCurr->SetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
279 TMath::Sqrt(tempError * tempError + binError2 * prob[smearSpeciesByError]));
280 }
281 else
282 printf("Problem obtaining fractions for: trackPt %f, jetPt %f, cent %f, smearSpeciesByError %d\n",
283 binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, smearSpeciesByError);
284 }
285 }
286 }
287
288 delete hPartFractionsSmeared;
289
290 const Int_t nHistos = nGenerations;
291
292
293 /*OLD error for each species by variation of ALL species (will bias towards smaller errors!)
294 // TODO Still does not store all the values in histogram to avoid different fractions for same map bin in same iteration.
295 // => If this is build in, one needs a histogram for EACH species!!!
296
297 // Calculate fractions/spectra w/ smearing "nGenerations" times for every species and obtain the statistical error from
298 // the spread of the results
299
300 TH2D* hIDFFtrackPtVaried[AliPID::kSPECIES][AliPID::kSPECIES * nGenerations];
301 TH2D* hIDFFzVaried[AliPID::kSPECIES][AliPID::kSPECIES * nGenerations];
302 TH2D* hIDFFxiVaried[AliPID::kSPECIES][AliPID::kSPECIES * nGenerations];
303
304 // Generate results with varying fractions
305 for (Int_t smearSpeciesByError = 0; smearSpeciesByError < AliPID::kSPECIES; smearSpeciesByError++) {
306 for (Int_t i = 0; i < nGenerations; i++) {
307 for (Int_t consideredSpecies = 0; consideredSpecies < AliPID::kSPECIES; consideredSpecies++) {
308 hIDFFtrackPtVaried[consideredSpecies][smearSpeciesByError * nGenerations + i] =
309 new TH2D(*(TH2D*)hIDFFtrackPt[consideredSpecies]);
310 hIDFFtrackPtVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->Reset();
311 hIDFFtrackPtVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->SetName(Form("hIDFFtrackPtVaried_%d_%d",
312 consideredSpecies,
313 smearSpeciesByError * nGenerations + i));
314 hIDFFzVaried[consideredSpecies][smearSpeciesByError * nGenerations + i] =
315 new TH2D(*(TH2D*)hIDFFz[consideredSpecies]);
316 hIDFFzVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->Reset();
317 hIDFFzVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->SetName(Form("hIDFFzVaried_%d_%d",
318 consideredSpecies,
319 smearSpeciesByError * nGenerations + i));
320 hIDFFxiVaried[consideredSpecies][smearSpeciesByError * nGenerations + i] =
321 new TH2D(*(TH2D*)hIDFFxi[consideredSpecies]);
322 hIDFFxiVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->Reset();
323 hIDFFxiVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->SetName(Form("hIDFFxiVaried_%d_%d",
324 consideredSpecies,
325 smearSpeciesByError * nGenerations + i));
326 }
327
328 for (Long64_t bin = 0; bin < nBinsTHnSparseProj; bin++) {
329 binContent = hDataProj->GetBinContent(bin, binCoord);
330 binError2 = hDataProj->GetBinError2(bin);
331
332 // If the bin is empty, do not compute the particle fractions -> This bin contributes nothing anyway
333 if (binContent < 2. * std::numeric_limits<double>::min())
334 continue;
335
336 for (Int_t dim = 0; dim < nDimProj; dim++)
337 binCentre[dim] = hDataProj->GetAxis(dim)->GetBinCenter(binCoord[dim]);
338
339 success = pidTask->GetParticleFractions(binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, prob,
340 smearSpeciesByError);
341
342 if (success) {
343 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
344 // To make things readable...
345 TH2D* hIDFFtrackPtVariedCurr = hIDFFtrackPtVaried[species][smearSpeciesByError * nGenerations + i];
346 TH2D* hIDFFzVariedCurr = hIDFFzVaried[species][smearSpeciesByError * nGenerations + i];
347 TH2D* hIDFFxiVariedCurr = hIDFFxiVaried[species][smearSpeciesByError * nGenerations + i];
348
349 // Track pT
350 hIDFFtrackPtVariedCurr->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
351 hIDFFtrackPtVariedCurr->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj])
352 + binContent * prob[species]);
353 Double_t tempError = hIDFFtrackPtVariedCurr->GetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
354 hIDFFtrackPtVariedCurr->SetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
355 TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
356
357 // z
358 hIDFFzVariedCurr->SetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj],
359 hIDFFzVariedCurr->GetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj])
360 + binContent * prob[species]);
361 tempError = hIDFFzVariedCurr->GetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj]);
362 hIDFFzVariedCurr->SetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj],
363 TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
364
365 // xi
366 hIDFFxiVariedCurr->SetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
367 hIDFFxiVariedCurr->GetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj])
368 + binContent * prob[species]);
369 tempError = hIDFFxiVariedCurr->GetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj]);
370 hIDFFxiVariedCurr->SetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
371 TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
372 }
373 }
374 else
375 printf("Problem obtaining fractions for: trackPt %f, jetPt %f, cent %f, smearSpeciesByError %d\n",
376 binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, smearSpeciesByError);
377 }
378 }
379 }
380
381 const Int_t nHistos = AliPID::kSPECIES * nGenerations;
382 */
383
384 // Compare results to obtain error
385 const Bool_t ignoreSigmaErrors = kTRUE;
386
387 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
388 if (!extractSystematicError(nHistos, hIDFFtrackPtVaried[species], hIDFFtrackPt[species], setMean, addErrorsQuadratically,
389 ignoreSigmaErrors))
390 printf("Failed to determine systematic error for trackPt, species %d\n", species);
391
392 if (!extractSystematicError(nHistos, hIDFFzVaried[species], hIDFFz[species], setMean, addErrorsQuadratically, ignoreSigmaErrors))
393 printf("Failed to determine systematic error for z, species %d\n", species);
394
395 if (!extractSystematicError(nHistos, hIDFFxiVaried[species], hIDFFxi[species], setMean, addErrorsQuadratically,
396 ignoreSigmaErrors))
397 printf("Failed to determine systematic error for xi, species %d\n", species);
398
399 for (Int_t i = 0; i < nHistos; i++) {
400 delete hIDFFtrackPtVaried[species][i];
401 hIDFFtrackPtVaried[species][i] = 0x0;
402
403 delete hIDFFzVaried[species][i];
404 hIDFFzVaried[species][i] = 0x0;
405
406 delete hIDFFxiVaried[species][i];
407 hIDFFxiVaried[species][i] = 0x0;
408 }
409 }
410 }
411}
412
413
414//___________________________________________________________________
415Int_t extractFFs(TString particleFractionPackagePathName, TString pathNameData, TString listName /*= ""*/,
416 Bool_t uniformSystematicError, Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
417 Double_t lowerCentrality = -2, Double_t upperCentrality = -2, Int_t rebinPt = 1, Int_t rebinZ = 1, Int_t rebinXi = 1)
418{
419 TObjArray* histList = 0x0;
420
421 if (listName == "") {
422 listName = pathNameData;
423 listName.Replace(0, listName.Last('/') + 1, "");
424 listName.ReplaceAll(".root", "");
425 }
426
427 TFile* f = TFile::Open(pathNameData.Data());
428 if (!f) {
429 std::cout << std::endl;
430 std::cout << "Failed to open file \"" << pathNameData.Data() << "\"!" << std::endl;
431 return -1;
432 }
433
434 histList = (TObjArray*)(f->Get(listName.Data()));
435 if (!histList) {
436 std::cout << std::endl;
437 std::cout << "Failed to load list \"" << listName.Data() << "\"!" << std::endl;
438 return -1;
439 }
440
441 // Extract the data histogram
442 THnSparse* hPIDdata = dynamic_cast<THnSparse*>(histList->FindObject("hPIDdataAll"));
443 if (!hPIDdata) {
444 std::cout << std::endl;
445 std::cout << "Failed to load data histo!" << std::endl;
446 return -1;
447 }
448
449 // Set proper errors, if not yet calculated
450 if (!hPIDdata->GetCalculateErrors()) {
451 std::cout << "Re-calculating errors of " << hPIDdata->GetName() << "..." << std::endl;
452 hPIDdata->Sumw2();
453 Long64_t nBinsTHnSparse = hPIDdata->GetNbins();
454 Double_t binContent = 0;
455
456 for (Long64_t bin = 0; bin < nBinsTHnSparse; bin++) {
457 binContent = hPIDdata->GetBinContent(bin);
458 hPIDdata->SetBinError(bin, TMath::Sqrt(binContent));
459 }
460 }
461
462 // If desired, restrict centrality axis
463 Int_t lowerCentralityBinLimit = -1;
464 Int_t upperCentralityBinLimit = -2; // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept
465 Bool_t restrictCentralityAxis = kFALSE;
466 Double_t actualLowerCentrality = -1.;
467 Double_t actualUpperCentrality = -1.;
468
469 if (lowerCentrality >= -1 && upperCentrality >= -1) {
470 // Add subtract a very small number to avoid problems with values right on the border between to bins
471 lowerCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(lowerCentrality + 0.001);
472 upperCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(upperCentrality - 0.001);
473
474 // Check if the values look reasonable
475 if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
476 && upperCentralityBinLimit <= hPIDdata->GetAxis(kPidCentrality)->GetNbins()) {
477 actualLowerCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinLowEdge(lowerCentralityBinLimit);
478 actualUpperCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinUpEdge(upperCentralityBinLimit);
479
480 restrictCentralityAxis = kTRUE;
481 }
482 else {
483 std::cout << std::endl;
484 std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
485 return -1;
486 }
487 }
488
489 std::cout << "centrality: ";
490 if (restrictCentralityAxis) {
491 std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
492 }
493 else {
494 std::cout << "All" << std::endl;
495 }
496
497 if (restrictCentralityAxis) {
498 hPIDdata->GetAxis(kPidCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
499 }
500
501 // If desired, restrict charge axis
502 std::cout << "Charge selection: ";
503 if (chargeMode == kAllCharged)
504 std::cout << "All charged particles" << std::endl;
505 else if (chargeMode == kNegCharge)
506 std::cout << "Negative particles only" << std::endl;
507 else if (chargeMode == kPosCharge)
508 std::cout << "Positive particles only" << std::endl;
509 else {
510 std::cout << "Unknown -> ERROR" << std::endl;
511 return -1;
512 }
513
514 const Bool_t restrictCharge = (chargeMode != kAllCharged);
515
516 const Int_t indexChargeAxisData = GetAxisByTitle(hPIDdata, "Charge (e_{0})");
517 if (indexChargeAxisData < 0 && restrictCharge) {
518 std::cout << "Error: Charge axis not found for data histogram!" << std::endl;
519 return -1;
520 }
521 Int_t lowerChargeBinLimitData = -1;
522 Int_t upperChargeBinLimitData = -2;
523 Double_t actualLowerChargeData = -999;
524 Double_t actualUpperChargeData = -999;
525
526 if (restrictCharge) {
527 // Add subtract a very small number to avoid problems with values right on the border between to bins
528 if (chargeMode == kNegCharge) {
529 lowerChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(-1. + 0.001);
530 upperChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(0. - 0.001);
531 }
532 else if (chargeMode == kPosCharge) {
533 lowerChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(0. + 0.001);
534 upperChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(1. - 0.001);
535 }
536
537 // Check if the values look reasonable
538 if (lowerChargeBinLimitData <= upperChargeBinLimitData && lowerChargeBinLimitData >= 1
539 && upperChargeBinLimitData <= hPIDdata->GetAxis(indexChargeAxisData)->GetNbins()) {
540 actualLowerChargeData = hPIDdata->GetAxis(indexChargeAxisData)->GetBinLowEdge(lowerChargeBinLimitData);
541 actualUpperChargeData = hPIDdata->GetAxis(indexChargeAxisData)->GetBinUpEdge(upperChargeBinLimitData);
542
543 std::cout << "Charge range data: " << actualLowerChargeData << " - " << actualUpperChargeData << std::endl;
544 }
545 else {
546 std::cout << std::endl;
547 std::cout << "Requested charge range out of limits or upper and lower limit are switched!" << std::endl;
548 return -1;
549 }
550
551 hPIDdata->GetAxis(indexChargeAxisData)->SetRange(lowerChargeBinLimitData, upperChargeBinLimitData);
552 }
553
554 // Just take one arbitrary selectSpecies to avoid multiple counting
555 hPIDdata->GetAxis(kPidSelectSpecies)->SetRange(1, 1);
556
557 // Get projection on dimensions relevant for FFs
558 const Int_t nDimProj = kNDimsProj;
559 Int_t dimProj[nDimProj];
560 dimProj[kPidPtProj] = kPidPt;
561 dimProj[kPidJetPtProj] = kPidJetPt;
562 dimProj[kPidZProj] = kPidZ;
563 dimProj[kPidXiProj] =kPidXi;
564
565 THnSparse* hDataProj = hPIDdata->Projection(nDimProj, dimProj, "e");
566
567 // If desired, rebin axes
568 if (rebinPt != 1 || rebinZ != 1 || rebinXi != 1) {
569 Int_t group[hDataProj->GetNdimensions()];
570
571 for (Int_t i = 0; i < hDataProj->GetNdimensions(); i++) {
572 if (i == kPidPtProj)
573 group[i] = rebinPt;
574 else if (i == kPidZProj)
575 group[i] = rebinZ;
576 else if (i == kPidXiProj)
577 group[i] = rebinXi;
578 else
579 group[i] = 1;
580 }
581
582 THnSparse* hTemp = hDataProj->Rebin(group);
583 hDataProj->SetName("temp");
584 delete hDataProj;
585
586 hDataProj = hTemp;
587 }
588
589 /* OLD - Now normalisation to nJets
590 Double_t numEvents = -1;
591 TH1* hNumEvents = dynamic_cast<TH1*>(histList->FindObject("fhEventsProcessed"));
592 if (!hNumEvents) {
593 std::cout << std::endl;
594 std::cout << "Histo with number of processed events not found! Yields will NOT be normalised to this number!" << std::endl
595 << std::endl;
596 }
597 else {
598 numEvents = hNumEvents->Integral(lowerCentralityBinLimit, upperCentralityBinLimit);
599
600 if (numEvents <= 0) {
601 numEvents = -1;
602 std::cout << std::endl;
603 std::cout << "Number of processed events < 1 in selected range! Yields will NOT be normalised to this number!"
604 << std::endl << std::endl;
605 }
606 }*/
607
608
609 TH2D* hNjetsGen = 0x0;
610 TH2D* hNjetsRec = 0x0;
611
612 TH2D* hMCgenPrimYieldPt[AliPID::kSPECIES];
613 TH2D* hMCgenPrimYieldZ[AliPID::kSPECIES];
614 TH2D* hMCgenPrimYieldXi[AliPID::kSPECIES];
615
616 hNjetsGen = (TH2D*)histList->FindObject("fh2FFJetPtGen");
617 hNjetsRec = (TH2D*)histList->FindObject("fh2FFJetPtRec");
618
619 if (!hNjetsRec) {
620 printf("Failed to load number of jets (rec) histo!\n");
621
622 // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other
623 // period is considered; also: No multiplicity information)
624 TFile* fBackward = TFile::Open("finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/AnalysisResults.root");
625
626 TString dirDataInFile = "";
627 TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0;
628
629 TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0;
630
631 TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCutsInc") : 0x0;
632
633 if (hFFJetPtRec) {
634 printf("***WARNING: For backward compatibility, using file \"finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/AnalysisResults.root\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n");
635 printf("ALSO: Using Njets for inclusive jets!!!!\n");
636
637 hNjetsRec = new TH2D("fh2FFJetPtRec", "", 1, -1, 1, hDataProj->GetAxis(kPidJetPtProj)->GetNbins(),
638 hDataProj->GetAxis(kPidJetPtProj)->GetXbins()->GetArray());
639
640 for (Int_t iJet = 1; iJet <= hNjetsRec->GetNbinsY(); iJet++) {
641 Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
642 Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
643 hNjetsRec->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin));
644 }
645 }
646
647 if (!hNjetsRec)
648 return -1;
649 }
650
651 THnSparse* hMCgeneratedYieldsPrimaries = (THnSparse*)histList->FindObject("fhMCgeneratedYieldsPrimaries");
652
653 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
654 hMCgenPrimYieldPt[i] = 0x0;
655 hMCgenPrimYieldZ[i] = 0x0;
656 hMCgenPrimYieldXi[i] = 0x0;
657 }
658
659 if (hMCgeneratedYieldsPrimaries && hMCgeneratedYieldsPrimaries->GetEntries() > 0) {
660 // Set proper errors, if not yet calculated
661 if (!hMCgeneratedYieldsPrimaries->GetCalculateErrors()) {
662 std::cout << "Re-calculating errors of " << hMCgeneratedYieldsPrimaries->GetName() << "..." << std::endl;
663
664 hMCgeneratedYieldsPrimaries->Sumw2();
665
666 Long64_t nBinsTHnSparseGenYield = hMCgeneratedYieldsPrimaries->GetNbins();
667 Double_t binContentGenYield = 0;
668 for (Long64_t bin = 0; bin < nBinsTHnSparseGenYield; bin++) {
669 binContentGenYield = hMCgeneratedYieldsPrimaries->GetBinContent(bin);
670 hMCgeneratedYieldsPrimaries->SetBinError(bin, TMath::Sqrt(binContentGenYield));
671 }
672 }
673
674 // If desired, rebin axes
675 if (rebinPt != 1 || rebinZ != 1 || rebinXi != 1) {
676 Int_t group[hMCgeneratedYieldsPrimaries->GetNdimensions()];
677
678 for (Int_t k = 0; k < hMCgeneratedYieldsPrimaries->GetNdimensions(); k++) {
679 if (k == kPidGenYieldPt)
680 group[k] = rebinPt;
681 else if (k == kPidGenYieldZ)
682 group[k] = rebinZ;
683 else if (k == kPidGenYieldXi)
684 group[k] = rebinXi;
685 else
686 group[k] = 1;
687 }
688
689 THnSparse* hTemp = hMCgeneratedYieldsPrimaries->Rebin(group);
690 hMCgeneratedYieldsPrimaries->SetName("temp");
691 delete hMCgeneratedYieldsPrimaries;
692
693 hMCgeneratedYieldsPrimaries = hTemp;
694 }
695
696
697 if (restrictCentralityAxis)
698 hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
699
700 if (restrictCharge) {
701 const Int_t indexChargeAxisGenYield = GetAxisByTitle(hMCgeneratedYieldsPrimaries, "Charge (e_{0})");
702 if (indexChargeAxisGenYield < 0) {
703 std::cout << "Error: Charge axis not found for gen yield histogram!" << std::endl;
704 return -1;
705 }
706
707 Int_t lowerChargeBinLimitGenYield = -1;
708 Int_t upperChargeBinLimitGenYield = -2;
709 Double_t actualLowerChargeGenYield = -999;
710 Double_t actualUpperChargeGenYield = -999;
711
712 // Add subtract a very small number to avoid problems with values right on the border between to bins
713 if (chargeMode == kNegCharge) {
714 lowerChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(-1. + 0.001);
715 upperChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(0. - 0.001);
716 }
717 else if (chargeMode == kPosCharge) {
718 lowerChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(0. + 0.001);
719 upperChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(1. - 0.001);
720 }
721
722 // Check if the values look reasonable
723 if (lowerChargeBinLimitGenYield <= upperChargeBinLimitGenYield && lowerChargeBinLimitGenYield >= 1
724 && upperChargeBinLimitGenYield <= hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetNbins()) {
725 actualLowerChargeGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetBinLowEdge(lowerChargeBinLimitGenYield);
726 actualUpperChargeGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetBinUpEdge(upperChargeBinLimitGenYield);
727
728 if (TMath::Abs(actualLowerChargeGenYield - actualLowerChargeData) > 1e-4 ||
729 TMath::Abs(actualUpperChargeGenYield - actualUpperChargeData) > 1e-4) {
730 std::cout << std::endl;
731 std::cout << "Error: Charge range gen yield: " << actualLowerChargeGenYield << " - " << actualUpperChargeGenYield
732 << std::endl << "differs from that of data: " << actualLowerChargeData << " - " << actualUpperChargeData
733 << std::endl;
734 return -1;
735 }
736 }
737 else {
738 std::cout << std::endl;
739 std::cout << "Requested charge range (gen yield) out of limits or upper and lower limit are switched!" << std::endl;
740 return -1;
741 }
742
743 hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->SetRange(lowerChargeBinLimitGenYield,
744 upperChargeBinLimitGenYield);
745 }
746
747 for (Int_t MCid = 0; MCid < AliPID::kSPECIES; MCid++) {
748 hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldMCpid)->SetRange(MCid + 1, MCid + 1);
749
750 hMCgenPrimYieldPt[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldJetPt, kPidGenYieldPt, "e");
751 hMCgenPrimYieldPt[MCid]->SetName(Form("hMCgenYieldsPrimPt_%s", AliPID::ParticleShortName(MCid)));
752 hMCgenPrimYieldPt[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
753 hMCgenPrimYieldPt[MCid]->GetZaxis()->SetTitle("1/N_{Jets} dN/dP_{T} (GeV/c)^{-1}");
754 hMCgenPrimYieldPt[MCid]->SetStats(kFALSE);
755
756 hMCgenPrimYieldZ[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldJetPt, kPidGenYieldZ, "e");
757 hMCgenPrimYieldZ[MCid]->SetName(Form("hMCgenYieldsPrimZ_%s", AliPID::ParticleShortName(MCid)));
758 hMCgenPrimYieldZ[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
759 hMCgenPrimYieldZ[MCid]->GetZaxis()->SetTitle("1/N_{Jets} dN/dz");
760 hMCgenPrimYieldZ[MCid]->SetStats(kFALSE);
761
762 hMCgenPrimYieldXi[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldJetPt, kPidGenYieldXi, "e");
763 hMCgenPrimYieldXi[MCid]->SetName(Form("hMCgenYieldsPrimXi_%s", AliPID::ParticleShortName(MCid)));
764 hMCgenPrimYieldXi[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
765 hMCgenPrimYieldXi[MCid]->GetZaxis()->SetTitle("1/N_{Jets} dN/d#xi");
766 hMCgenPrimYieldXi[MCid]->SetStats(kFALSE);
767
768 hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldMCpid)->SetRange(0, -1);
769 }
770 }
771
772
773 AliAnalysisTaskPID *pidTask = new AliAnalysisTaskPID("spectrumExtractorTask");
774
775 if (!pidTask->SetParticleFractionHistosFromFile(particleFractionPackagePathName, kFALSE)) {
776 printf("Failed to load particle fraction package from file \"%s\"!\n", particleFractionPackagePathName.Data());
777 return -1;
778 }
779
780 if (!pidTask->SetParticleFractionHistosFromFile(particleFractionPackagePathName, kTRUE)) {
781 printf("Failed to load particle fraction sys error package from file \"%s\"!\n", particleFractionPackagePathName.Data());
782 return -1;
783 }
784
785 printf("Loaded particle fraction package from file \"%s\"!\n", particleFractionPackagePathName.Data());
786
787 TH2D* hIDFFtrackPt[AliPID::kSPECIES];
788 TH2D* hIDFFz[AliPID::kSPECIES];
789 TH2D* hIDFFxi[AliPID::kSPECIES];
790
791 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
792 hIDFFtrackPt[i] = hDataProj->Projection(kPidJetPtProj, kPidPtProj);
793 hIDFFtrackPt[i]->Reset();
794 if (hIDFFtrackPt[i]->GetSumw2N() <= 0)
795 hIDFFtrackPt[i]->Sumw2();
796 setupHist(hIDFFtrackPt[i], Form("hIDFFtrackPt_%s", AliPID::ParticleShortName(i)), Form("%s", AliPID::ParticleShortName(i)),
797 hDataProj->GetAxis(kPidPtProj)->GetTitle(), hDataProj->GetAxis(kPidJetPtProj)->GetTitle(), kBlack);
798
799 hIDFFz[i] = hDataProj->Projection(kPidJetPtProj, kPidZProj);
800 hIDFFz[i]->Reset();
801 if (hIDFFz[i]->GetSumw2N() <= 0)
802 hIDFFz[i]->Sumw2();
803 setupHist(hIDFFz[i], Form("hIDFFz_%s", AliPID::ParticleShortName(i)), Form("%s", AliPID::ParticleShortName(i)),
804 hDataProj->GetAxis(kPidZProj)->GetTitle(), hDataProj->GetAxis(kPidJetPtProj)->GetTitle(), kBlack);
805
806 hIDFFxi[i] = hDataProj->Projection(kPidJetPtProj, kPidXiProj);
807 hIDFFxi[i]->Reset();
808 if (hIDFFxi[i]->GetSumw2N() <= 0)
809 hIDFFxi[i]->Sumw2();
810 setupHist(hIDFFxi[i], Form("hIDFFxi_%s", AliPID::ParticleShortName(i)), Form("%s", AliPID::ParticleShortName(i)),
811 hDataProj->GetAxis(kPidXiProj)->GetTitle(), hDataProj->GetAxis(kPidJetPtProj)->GetTitle(), kBlack);
812 }
813
814 const Double_t centrality = (actualLowerCentrality + actualUpperCentrality) / 2.;
815
816 // First iteration: Just take the default fractions to obtain the "mean" of the yields
817 Bool_t setMean = kTRUE;
818 Bool_t addErrorsQuadratically = kFALSE;
819 Bool_t smearByError = kFALSE;
820 Bool_t takeIntoAccountSysError = kFALSE;
821 // setMean and all the follwoing parameters are anyway irrelevant, if smearByError = kFALSE and takeIntoAccountSysError = kFALSE
822 GenerateParticleYields(hDataProj, pidTask, centrality, hIDFFtrackPt, hIDFFz, hIDFFxi, setMean, addErrorsQuadratically,
823 smearByError, uniformSystematicError, takeIntoAccountSysError, 1);
824
825
826 // Next iteration: Vary the PID map within statistical errors several times and calculate the resulting statistical error
827 // of the spectra from this. But since the mean of the fraction is some kind of "best estimate" of the truth, leave the means
828 // as they have been set during the last step.
829
830
831 // NOTE: Do NOT add the statistical errors from the last step, since the yields/fractions (rel. error is the same) already
832 // contain the stat. uncertainty of the yield of this species in the corresponding bin! I.e. one just takes some kind of weighted
833 // mean of the fractions (mean and error) (equivalent of taking the yields with corresponding error).
834 setMean = kFALSE;
835 addErrorsQuadratically = kFALSE;
836 smearByError = kTRUE;
837 takeIntoAccountSysError = kFALSE;
838 const Int_t nGenerations = 5000; //TODO set to 5000 in the end, after all testing was successful
839 GenerateParticleYields(hDataProj, pidTask, centrality, hIDFFtrackPt, hIDFFz, hIDFFxi, setMean, addErrorsQuadratically,
840 smearByError, uniformSystematicError, takeIntoAccountSysError, nGenerations);
841
842
843 // Next iteration: Vary the PID map within systematic errors several times and calculate the resulting systematic error
844 // of the spectra from this. But since the mean of the fraction is some kind of "best estimate" of the truth, leave the means
845 // as they have been set during the last step.
846 setMean = kFALSE;
847 addErrorsQuadratically = kFALSE;
848 smearByError = kFALSE;
849 takeIntoAccountSysError = kTRUE;
850 // Clone histograms with final statistical errors -> Only set systematic errors, but leave mean as it is
851 TH2D* hIDFFtrackPtSysError[AliPID::kSPECIES];
852 TH2D* hIDFFzSysError[AliPID::kSPECIES];
853 TH2D* hIDFFxiSysError[AliPID::kSPECIES];
854
855 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
856 hIDFFtrackPtSysError[i] = new TH2D(*hIDFFtrackPt[i]);
857 hIDFFtrackPtSysError[i]->SetName(Form("%s_sysError", hIDFFtrackPt[i]->GetName()));
858 hIDFFtrackPtSysError[i]->SetFillStyle(0);
859
860 hIDFFzSysError[i] = new TH2D(*hIDFFz[i]);
861 hIDFFzSysError[i]->SetName(Form("%s_sysError", hIDFFz[i]->GetName()));
862 hIDFFzSysError[i]->SetFillStyle(0);
863
864 hIDFFxiSysError[i] = new TH2D(*hIDFFxi[i]);
865 hIDFFxiSysError[i]->SetName(Form("%s_sysError", hIDFFxi[i]->GetName()));
866 hIDFFxiSysError[i]->SetFillStyle(0);
867 }
868
869 GenerateParticleYields(hDataProj, pidTask, centrality, hIDFFtrackPtSysError, hIDFFzSysError, hIDFFxiSysError, setMean,
870 addErrorsQuadratically, smearByError, uniformSystematicError, takeIntoAccountSysError, nGenerations);
871
872 delete pidTask;
873
874
875 // Normalise properly
876 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
877 normaliseYieldHist2D(hIDFFtrackPt[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
878 normaliseYieldHist2D(hIDFFz[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
879 normaliseYieldHist2D(hIDFFxi[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
880
881 normaliseYieldHist2D(hIDFFtrackPtSysError[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
882 normaliseYieldHist2D(hIDFFzSysError[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
883 normaliseYieldHist2D(hIDFFxiSysError[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
884
885 normaliseYieldHist2D(hMCgenPrimYieldPt[species], hNjetsGen, lowerCentralityBinLimit, upperCentralityBinLimit);
886 normaliseYieldHist2D(hMCgenPrimYieldZ[species], hNjetsGen, lowerCentralityBinLimit, upperCentralityBinLimit);
887 normaliseYieldHist2D(hMCgenPrimYieldXi[species], hNjetsGen, lowerCentralityBinLimit, upperCentralityBinLimit);
888 }
889
890
891
892
893 // Save results to file
894 TString chargeString = "";
895 if (chargeMode == kPosCharge)
896 chargeString = "_posCharge";
897 else if (chargeMode == kNegCharge)
898 chargeString = "_negCharge";
899
900 TString saveFileName = pathNameData;
901 saveFileName.Replace(0, pathNameData.Last('/') + 1, "");
902
903 TString savePath = pathNameData;
904 savePath.ReplaceAll(Form("/%s", saveFileName.Data()), "");
905
906 saveFileName.Prepend("output_extractedFFs_");
907 saveFileName.ReplaceAll(".root", Form("_centrality_%s%s.root",
908 restrictCentralityAxis ? Form("%.0f_%.0f.root", actualLowerCentrality, actualUpperCentrality)
909 : "all",
910 chargeString.Data()));
911
912 TString saveFilePathName = Form("%s/%s", savePath.Data(), saveFileName.Data());
913 TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE");
914
915 if (!saveFile) {
916 printf("Failed to save results to file \"%s\"!\n", saveFilePathName.Data());
917 return -1;
918 }
919
920 saveFile->cd();
921
922
923 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
924 if (hIDFFtrackPt[i])
925 hIDFFtrackPt[i]->Write();
926
927 if (hIDFFz[i])
928 hIDFFz[i]->Write();
929
930 if (hIDFFxi[i])
931 hIDFFxi[i]->Write();
932
933
934 if (hIDFFtrackPtSysError[i])
935 hIDFFtrackPtSysError[i]->Write();
936
937 if (hIDFFzSysError[i])
938 hIDFFzSysError[i]->Write();
939
940 if (hIDFFxiSysError[i])
941 hIDFFxiSysError[i]->Write();
942
943
944 if (hMCgenPrimYieldPt[i])
945 hMCgenPrimYieldPt[i]->Write();
946
947 if (hMCgenPrimYieldZ[i])
948 hMCgenPrimYieldZ[i]->Write();
949
950 if (hMCgenPrimYieldXi[i])
951 hMCgenPrimYieldXi[i]->Write();
952 }
953
954 if (hNjetsGen)
955 hNjetsGen->Write();
956
957 if (hNjetsRec)
958 hNjetsRec->Write();
959
960
961 TNamed* settings = new TNamed(
962 Form("Settings: Fraction package \"%s\", Data file \"%s\", lowerCentrality %.0f, upperCentrality %.0f, chargeMode %d, uniformSystematicError %d, rebinPt %d, rebinZ %d, rebinXi %d\n",
963 particleFractionPackagePathName.Data(), pathNameData.Data(), lowerCentrality, upperCentrality, chargeMode,
964 uniformSystematicError, rebinPt, rebinZ, rebinXi), "");
965 settings->Write();
966
967 saveFile->Close();
968
969
970
971 /*
972 Int_t t1 = hIDFFtrackPt[AliPID::kPion]->GetYaxis()->FindBin(5.1);
973 Int_t t2 = hIDFFtrackPt[AliPID::kPion]->GetYaxis()->FindBin(9.9);
974
975 printf("t1 %d, t2 %d\n", t1, t2);
976 new TCanvas();
977 TH1D* hTemp = hIDFFtrackPt[AliPID::kPion]->ProjectionX("_pfy1", t1, t2, "e");
978 hTemp->SetFillStyle(0);
979 hTemp->Draw("");
980
981 hTemp = hIDFFtrackPtSysError[AliPID::kPion]->ProjectionX("sys_pfy1", t1, t2, "e");
982 hTemp->SetFillStyle(0);
983 hTemp->Draw("E2same");
984
985 if (hMCgenPrimYieldPt[AliPID::kPion]) {
986 hTemp = hMCgenPrimYieldPt[AliPID::kPion]->ProjectionX("MC_pfy1", t1, t2, "e");
987 hTemp->SetFillStyle(0);
988 hTemp->SetMarkerStyle(23);
989 hTemp->Draw("same");
990 }
991
992
993 new TCanvas();
994 hTemp = hIDFFz[AliPID::kPion]->ProjectionX("_pfy2", t1, t2, "e");
995 hTemp->SetFillStyle(0);
996 hTemp->Draw("");
997
998 hTemp = hIDFFzSysError[AliPID::kPion]->ProjectionX("sys_pfy2", t1, t2, "e");
999 hTemp->SetFillStyle(0);
1000 hTemp->Draw("E2same");
1001
1002 if (hMCgenPrimYieldZ[AliPID::kPion]) {
1003 hTemp = hMCgenPrimYieldZ[AliPID::kPion]->ProjectionX("MC_pfy2", t1, t2, "e");
1004 hTemp->SetFillStyle(0);
1005 hTemp->SetMarkerStyle(23);
1006 hTemp->Draw("same");
1007 }
1008
1009
1010 new TCanvas();
1011 hTemp = hIDFFxi[AliPID::kPion]->ProjectionX("_pfy3", t1, t2, "e");
1012 hTemp->SetFillStyle(0);
1013 hTemp->Draw("");
1014
1015 hTemp = hIDFFxiSysError[AliPID::kPion]->ProjectionX("sys_pfy3", t1, t2, "e");
1016 hTemp->SetFillStyle(0);
1017 hTemp->Draw("E2same");
1018
1019 if (hMCgenPrimYieldXi[AliPID::kPion]) {
1020 hTemp = hMCgenPrimYieldXi[AliPID::kPion]->ProjectionX("MC_pfy3", t1, t2, "e");
1021 hTemp->SetFillStyle(0);
1022 hTemp->SetMarkerStyle(23);
1023 hTemp->Draw("same");
1024 }
1025
1026
1027
1028
1029 new TCanvas();
1030 hTemp = hIDFFtrackPt[AliPID::kProton]->ProjectionX("_pfy4", t1, t2, "e");
1031 hTemp->SetFillStyle(0);
1032 hTemp->Draw("");
1033
1034 hTemp = hIDFFtrackPtSysError[AliPID::kProton]->ProjectionX("sys_pfy4", t1, t2, "e");
1035 hTemp->SetFillStyle(0);
1036 hTemp->Draw("E2same");
1037
1038 if (hMCgenPrimYieldPt[AliPID::kProton]) {
1039 hTemp = hMCgenPrimYieldPt[AliPID::kProton]->ProjectionX("MC_pfy4", t1, t2, "e");
1040 hTemp->SetFillStyle(0);
1041 hTemp->SetMarkerStyle(23);
1042 hTemp->Draw("same");
1043 }
1044
1045
1046 new TCanvas();
1047 hTemp = hIDFFz[AliPID::kProton]->ProjectionX("_pfy5", t1, t2, "e");
1048 hTemp->SetFillStyle(0);
1049 hTemp->Draw("");
1050
1051 hTemp = hIDFFzSysError[AliPID::kProton]->ProjectionX("sys_pfy5", t1, t2, "e");
1052 hTemp->SetFillStyle(0);
1053 hTemp->Draw("E2same");
1054
1055 if (hMCgenPrimYieldZ[AliPID::kProton]) {
1056 hTemp = hMCgenPrimYieldZ[AliPID::kProton]->ProjectionX("MC_pfy5", t1, t2, "e");
1057 hTemp->SetFillStyle(0);
1058 hTemp->SetMarkerStyle(23);
1059 hTemp->Draw("same");
1060 }
1061
1062
1063 new TCanvas();
1064 hTemp = hIDFFxi[AliPID::kProton]->ProjectionX("_pfy", t1, t2, "e");
1065 hTemp->SetFillStyle(0);
1066 hTemp->Draw("");
1067
1068 hTemp = hIDFFxiSysError[AliPID::kProton]->ProjectionX("sys_pfy", t1, t2, "e");
1069 hTemp->SetFillStyle(0);
1070 hTemp->Draw("E2same");
1071
1072 if (hMCgenPrimYieldXi[AliPID::kProton]) {
1073 hTemp = hMCgenPrimYieldXi[AliPID::kProton]->ProjectionX("MC_pfy", t1, t2, "e");
1074 hTemp->SetFillStyle(0);
1075 hTemp->SetMarkerStyle(23);
1076 hTemp->Draw("same");
1077 }*/
1078
1079 return 0;
1080}