]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/macros/CorrectForEfficiencypptrd.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update of the hfe package
[u/mrichter/AliRoot.git] / PWG3 / hfe / macros / CorrectForEfficiencypptrd.C
CommitLineData
8c1c76e9 1void CorrectForEfficiencypptrd(const char *filedata,const char *fileMC,const char *NameList,Bool_t withtrd=kTRUE, Float_t CorrType=0, Int_t conf=2, Bool_t smoothing=kFALSE, Bool_t hadroncontaminationsubtracted=kTRUE, Bool_t UnsetCorrelatedErrors=kFALSE);\r
2void CorrectForEfficiencyBeauty(const char *filedata,const char *fileMC, const char *fileMCbg);\r
3TList *GetResults(const Char_t *testfile,const Char_t *plus="");\r
4TList *GetQA(const Char_t *testfile,const Char_t *plus="");\r
5TObject* GetSpectrum(AliCFContainer *c, Int_t step);\r
6THnSparseF* GetHadronEffbyIPcut(TList *hfeqa);\r
7void CorrectFromTheWidth(TH1D *h1);\r
8\r
9void CorrectForEfficiencypptrd(const char *filedata,const char *fileMC,const char *NameList,Bool_t withtrd, Float_t CorrType, Int_t conf, Bool_t smoothing, Bool_t hadroncontaminationsubtracted, Bool_t UnsetCorrelatedErrors) {\r
10\r
11 ///////////////////////////////////////////////////////////////////////////////////////////////////////\r
12 // Will produce finalSpectrum.root with results inside\r
13 // TGraphErrors UnfoldingCorrectedSpectrum -> unfolding procedure (doesn't work in 2D with trunk)\r
14 // TGraphErrors AlltogetherSpectrum -> corrected spectrum with other method\r
15 // TH1D RatioUnfoldingAlltogetherSpectrum \r
16 // THnSparseF UnfoldingCorrectedNotNormalizedSpectrum -> unfolding procedure not yet normalized\r
17 // AliCFDataGrid AlltogetherCorrectedNotNormalizedSpectrum -> other method not yet normalized\r
18 //\r
19 // NameList - possible extension for the full name of the data and MC TLists (e.g. "_PID2" for HFE_Results_PID2)\r
20 // CorrType - 0 use fixed TPC PID 0.5\r
21 // 1 use TPC PID efficiency from V0s \r
22 //\r
23 // Conf - TRD efficiency\r
24 // double TRDeffs[5] = {0.7,0.8,0.9,0.75,0.85};\r
25 // conf = 2, TRDeff = 0.8...\r
26 // (Only relevant for withtrd = kTRUE)\r
27 //\r
28 // f/ smoothing kFALSE\r
29 //\r
30 // g/ Flag for hadron contamination subtraction\r
31 // kTRUE means we subtract the hadron contamination\r
32 // kFALSE not\r
33 //\r
34 ///////////////////////////////////////////////////////////////////////////////////////////////////////\r
35\r
36 double TRDeffs[5] = {0.7,0.8,0.9,0.75,0.85};\r
37\r
38 gStyle->SetPalette(1);\r
39 gStyle->SetOptStat(1111);\r
40 gStyle->SetPadBorderMode(0);\r
41 gStyle->SetCanvasColor(10);\r
42 gStyle->SetPadLeftMargin(0.13);\r
43 gStyle->SetPadRightMargin(0.13);\r
44\r
45 if (CorrType>1){\r
46 printf("Wrong argument for the type of correction, it is maybe missing, or it has to be <=2!!! \n");\r
47 return;\r
48 }\r
49\r
50 /////////////////////\r
51 // Take the stuff\r
52 /////////////////////\r
53 \r
54 printf("Get results for data \n");\r
55 TList *resultsdata = GetResults(filedata,NameList);\r
56 if(!resultsdata){\r
57 printf("No output objects for data: Calculation will terminate here\n");\r
58 return;\r
59 }\r
60\r
61 ///////////////////////////////////////////////////////////\r
62 // Event container for the normalization to CINB1\r
63 // Take the number of events as function of z and number of contributors\r
64 // For those without primary vertex from tracks, take the fraction in the z range from MC\r
65 // For 10cm: LHC10f6a (88.0756%), LHC10f6 (86.7461%)\r
66 ////////////////////////////////////////////////////////////\r
67 resultsdata->Print();\r
68\r
69 AliCFContainer *containerdata = (AliCFContainer *) resultsdata->FindObject("eventContainer");\r
70 if(!containerdata) {\r
71 printf("No container \n");\r
72 return;\r
73 }\r
74 AliCFDataGrid *dataGrid = (AliCFDataGrid *) GetSpectrum(containerdata,AliHFEcuts::kEventStepZRange);\r
75 //dataGrid->PrintBinLimits();\r
76 dataGrid->PrintNBins();\r
77 \r
78 // \r
79 TH2D *spectrum_zrange = (TH2D *) dataGrid->Project(0,3);\r
80 TAxis *axis_0 = spectrum_zrange->GetXaxis();\r
81 TAxis *axis_1 = spectrum_zrange->GetYaxis();\r
82 printf("Number of bins in x %d and in y %d\n",axis_0->GetNbins(),axis_1->GetNbins());\r
83 TH1D *spectrum1Da = (TH1D *) spectrum_zrange->ProjectionX("bin0",1,1);\r
84 TH1D *spectrum1Db = (TH1D *) spectrum_zrange->ProjectionX("bin>0",2,12);\r
85 TH1D *spectrum1Dc = (TH1D *) spectrum_zrange->ProjectionX("total");\r
86 Double_t nbinbin0 = spectrum1Da->Integral();\r
87 Double_t nbinnobin0 = spectrum1Db->Integral();\r
88 Double_t nbintotal = spectrum1Dc->Integral();\r
89\r
90 //TCanvas *c1teste = new TCanvas("eventcontainertest","eventcontainertest",800,800);\r
91 //c1teste->cd(1);\r
92 //spectrum_zrange->Draw("colz");\r
93 \r
94\r
95 Float_t numberofentries = nbinnobin0+nbinbin0*0.880756;\r
96 printf("Number in bin0 %f, number out %f, total %f, normalized %f\n",nbinbin0,nbinnobin0,nbintotal,numberofentries);\r
97 \r
98 //////////////////////////////\r
99 // Take more stuff\r
100 ///////////////////////////////\r
101 \r
102 printf("Get results for MC \n");\r
103 TList *resultsmc = GetResults(fileMC,NameList);\r
104 if(!resultsmc){\r
105 printf("No output objects for mc: Calculation will terminate here\n");\r
106 return;\r
107 }\r
108\r
109 AliHFEcontainer *datahfecontainer = (AliHFEcontainer *) resultsdata->FindObject("trackContainer");\r
110 if(!datahfecontainer) {\r
111 printf("No container for data \n");\r
112 return;\r
113 }\r
114 AliCFContainer *sumcontainer = datahfecontainer->GetCFContainer("recTrackContReco");\r
115\r
116 ////////////////////////////////////////////\r
117 //Plot raw spectrum for TPC TOF scenario\r
118 ////////////////////////////////////////////\r
119 AliCFDataGrid *dataGrida = 0x0;\r
120 if(!withtrd) dataGrida = (AliCFDataGrid *) GetSpectrum(sumcontainer,AliHFEcuts::kStepHFEcutsTRD + 2);\r
121 else dataGrida = (AliCFDataGrid *) GetSpectrum(sumcontainer,AliHFEcuts::kStepHFEcutsTRD + 3);\r
122 \r
123 TH1D *spectrumpt = (TH1D *) dataGrida->Project(0);\r
124 CorrectFromTheWidth(spectrumpt);\r
125 \r
126 TH1D *total = new TH1D("total","",1,0.3,10.3);\r
127 total->SetMaximum(1.0e+07);\r
128 total->SetMinimum(1000);\r
129 total->SetXTitle("p_{T} [GeV/c]");\r
130 total->SetYTitle("dN/dp_{T}[GeV/c]^{-1}");\r
131 total->SetTitleOffset(1.5,"Y");\r
132 //total->Scale(0.9);\r
133 total->GetXaxis()->SetRangeUser(0.38,10.3);\r
134 \r
135\r
136 TCanvas *c1test = new TCanvas("rawspectrum","rawspectrum",800,800);\r
137 c1test->cd(1);\r
138 gPad->SetLeftMargin(0.13);\r
139 gPad->SetLogy();\r
140 gPad->SetTicks();\r
141 //gPad->SetGridx();\r
142 //gPad->SetGridy();\r
143 gPad->SetFillColor(10);\r
144 gPad->SetFrameFillColor(10);\r
145 gStyle->SetOptStat(0);\r
146 gStyle->SetOptFit(1111);\r
147 gStyle->SetOptTitle(0);\r
148 total->SetStats(0);\r
149 spectrumpt->SetStats(0);\r
150 total->Draw();\r
151 spectrumpt->Draw("same");\r
152 TPaveText* t1=new TPaveText(0.49,0.45,0.79,0.52,"NDC");\r
153 t1->SetFillStyle(0);\r
154 t1->SetBorderSize(0);\r
155 t1->AddText(0.,0.,"pp, #sqrt{s} = 7 TeV");\r
156 t1->SetTextColor(kRed);\r
157 //t1->SetTextSize(20);\r
158 t1->SetTextFont(42);\r
159 t1->Draw();\r
160 TPaveText* t2=new TPaveText(0.49,0.35,0.79,0.42,"NDC");\r
161 t2->SetFillStyle(0);\r
162 t2->SetBorderSize(0);\r
163 t2->AddText(0.,0.,"L = 2.6 nb^{-1}");\r
164 t2->SetTextColor(kRed);\r
165 //t1->SetTextSize(20);\r
166 t2->SetTextFont(42);\r
167 t2->Draw();\r
168 TPaveText* t3=new TPaveText(0.49,0.35,0.79,0.42,"NDC");\r
169 t3->SetFillStyle(0);\r
170 t3->SetBorderSize(0);\r
171 t3->AddText(0.,0.,"|#eta| < 0.8");\r
172 t3->SetTextColor(kRed);\r
173 //t1->SetTextSize(20);\r
174 t3->SetTextFont(42);\r
175 t3->Draw();\r
176 \r
177 /////////////////////////////\r
178 // Check number of events\r
179 /////////////////////////////\r
180\r
181\r
182 Int_t numberOfEventsafterallcuts = (Int_t) datahfecontainer->GetNumberOfEvents();\r
183 \r
184 printf("Number of events not corrected %d\n",numberOfEventsafterallcuts);\r
185 printf("Number of events corrected %f\n",numberofentries);\r
186\r
187 AliHFEcontainer *mchfecontainer = (AliHFEcontainer *) resultsmc->FindObject("trackContainer");\r
188 if(!mchfecontainer) {\r
189 printf("No mc container \n");\r
190 return;\r
191 }\r
192 \r
193 printf("Find the container V0\n"); \r
194 AliHFEcontainer *containerhfeV0 = (AliHFEcontainer *) resultsdata->FindObject("containerV0");\r
195 if(!containerhfeV0) {\r
196 printf("No hfe container \n");\r
197 return;\r
198 }\r
199 \r
200 //////////////\r
201 // Correct\r
202 /////////////\r
203\r
204 AliHFEspectrum *spectrum = new AliHFEspectrum("HFEspectrum");\r
205 spectrum->SetNbDimensions(1);\r
206 // If you want to correct positive (0) or negative (1) separately\r
207 //spectrum->SetChargeChoosen(0);\r
208 spectrum->SetSmoothing(smoothing);\r
209 spectrum->SetUnSetCorrelatedErrors(UnsetCorrelatedErrors);\r
210 spectrum->SetDebugLevel(1);\r
211 spectrum->SetNumberOfEvents((Int_t)numberofentries);\r
212 spectrum->SetDumpToFile(kTRUE);\r
213 // True step in MC (events in range +- 10 cm and no pile up)\r
214 spectrum->SetMCTruthStep(AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine);\r
215\r
216 // Step where we correct from MC (tracking + TOF)\r
217 spectrum->SetMCEffStep(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1);\r
218 printf("!!! Step we correct with MC %d\n!!!",AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1);\r
219 \r
220\r
221 // Step from data we correct for\r
222 //////////****************************//////////////////////////////\r
223 if(withtrd) spectrum->SetStepToCorrect(AliHFEcuts::kStepHFEcutsTRD + 3);\r
224 else spectrum->SetStepToCorrect(AliHFEcuts::kStepHFEcutsTRD + 2);\r
225 //////////****************************//////////////////////////////\r
226\r
227 // Steps to be corrected with V0 \r
228 if (CorrType == 1){\r
229 printf("\n This time we correct with PID efficiency from V0s !!! \n");\r
230 if(withtrd) {\r
231 spectrum->SetStepBeforeCutsV0(AliHFEcuts::kStepHFEcutsTRD + 2);\r
232 spectrum->SetStepAfterCutsV0(AliHFEcuts::kStepHFEcutsTRD + 3);\r
233 }\r
234 else {\r
235 spectrum->SetStepBeforeCutsV0(AliHFEcuts::kStepHFEcutsTRD + 1);\r
236 spectrum->SetStepAfterCutsV0(AliHFEcuts::kStepHFEcutsTRD + 2);\r
237 }\r
238 }\r
239\r
240 if (CorrType < 1){\r
241 \r
242 Double_t effTPC = 0.5,\r
243 effTRD = TRDeffs[conf-1]; \r
244 printf("effTRD = %f\n", effTRD);\r
245 // Place for a correction as an efficiency parametrized as function of pt of pt,eta,phi or p=pt\r
246 TF1 *hEfficiency = new TF1("efficiency", "[0]", 0., 20.);\r
247 if(withtrd) hEfficiency->SetParameter(0, effTPC*effTRD);\r
248 else hEfficiency->SetParameter(0, effTPC);\r
249 spectrum->SetEfficiencyFunction(hEfficiency);\r
250 } \r
251\r
252 // Give everything (data container, mc container and V0 data container)\r
253 if (CorrType == 1){\r
254 printf("\n This time we correct with PID efficiency from V0s !!! \n");\r
255 spectrum->Init(datahfecontainer,mchfecontainer,containerhfeV0);\r
256 }\r
257 else {\r
258 printf("\n This time we correct with PID efficiency from MC or fixed value !!! \n");\r
259 spectrum->Init(datahfecontainer,mchfecontainer);\r
260 }\r
261 \r
262 // kTRUE means subtract hadron background, kFALSE means do not subtract hadron background\r
263 spectrum->Correct(hadroncontaminationsubtracted);\r
264 \r
265\r
266}\r
267\r
268//_____________________________________________________________________\r
269void CorrectForEfficiencyBeauty(const char *filedata,const char *fileMC, const char *fileMCbg) {\r
270 \r
271 ///////////////////////////////////////////////////////////////////////////////////////////////////////\r
272 // Will produce finalSpectrum.root with results inside\r
273 // TGraphErrors UnfoldingCorrectedSpectrum -> unfolding procedure (doesn't work in 2D with trunk)\r
274 // TGraphErrors AlltogetherSpectrum -> corrected spectrum with other method\r
275 // TH1D RatioUnfoldingAlltogetherSpectrum \r
276 // THnSparseF UnfoldingCorrectedNotNormalizedSpectrum -> unfolding procedure not yet normalized\r
277 // AliCFDataGrid AlltogetherCorrectedNotNormalizedSpectrum -> other method not yet normalized\r
278 ///////////////////////////////////////////////////////////////////////////////////////////////////////\r
279\r
280 gStyle->SetPalette(1);\r
281 gStyle->SetOptStat(1111);\r
282 gStyle->SetPadBorderMode(0);\r
283 gStyle->SetCanvasColor(10);\r
284 gStyle->SetPadLeftMargin(0.13);\r
285 gStyle->SetPadRightMargin(0.13);\r
286\r
287 /////////////////////\r
288 // Take the stuff\r
289 /////////////////////\r
290\r
291 TList *resultsdata = GetResults(filedata,"_PID2");\r
292 //TList *resultsdata = GetResults(filedata);\r
293 if(!resultsdata){\r
294 printf("No output objects for data: Calculation will terminate here\n");\r
295 return;\r
296 }\r
297 ///////////////////////////////////////////////////////////\r
298 // Event container for the normalization to CINB1\r
299 // Take the number of events as function of z and number of contributors\r
300 // For those without primary vertex from tracks, take the fraction in the z range from MC\r
301 // For 10cm: LHC10f6a (88.0756%), LHC10f6 (86.7461%)\r
302 ////////////////////////////////////////////////////////////\r
303\r
304 AliCFContainer *containerdata = (AliCFContainer *) resultsdata->FindObject("eventContainer");\r
305 if(!containerdata) {\r
306 printf("No container \n");\r
307 return;\r
308 }\r
309 AliCFDataGrid *dataGrid = (AliCFDataGrid *) GetSpectrum(containerdata,AliHFEcuts::kEventStepZRange);\r
310 TH2D *spectrum_zrange = (TH2D *) dataGrid->Project(0,2);\r
311 TH1D *spectrum1Da = (TH1D *) spectrum_zrange->ProjectionX("bin0",1,1);\r
312 TH1D *spectrum1Db = (TH1D *) spectrum_zrange->ProjectionX("bin>0",2,12);\r
313 TH1D *spectrum1Dc = (TH1D *) spectrum_zrange->ProjectionX("total");\r
314 Double_t nbinbin0 = spectrum1Da->Integral();\r
315 Double_t nbinnobin0 = spectrum1Db->Integral();\r
316 Double_t nbintotal = spectrum1Dc->Integral();\r
317\r
318 Float_t numberofentries = nbinnobin0+nbinbin0*0.880756;\r
319 printf("Number in bin0 %f, number out %f, total %f, normalized %f\n",nbinbin0,nbinnobin0,nbintotal,numberofentries);\r
320\r
321\r
322 AliHFEcontainer *datahfecontainer = (AliHFEcontainer *) resultsdata->FindObject("trackContainer");\r
323 if(!datahfecontainer) {\r
324 printf("No container for data \n");\r
325 return;\r
326 }\r
327\r
328\r
329 //////////////////////////////\r
330 // Check MC # of events\r
331 ///////////////////////////////\r
332\r
333 TList *resultsmc = GetResults(fileMC,"_PID2");\r
334 if(!resultsmc){\r
335 printf("No output objects for mc: Calculation will terminate here\n");\r
336 return;\r
337 }\r
338\r
339 AliCFContainer *containermc = (AliCFContainer *) resultsmc->FindObject("eventContainer");\r
340 if(!containermc) {\r
341 printf("No container \n");\r
342 return;\r
343 }\r
344 AliCFDataGrid *mcGrid = (AliCFDataGrid *) GetSpectrum(containermc,AliHFEcuts::kEventStepZRange);\r
345 TH2D *spectrum_zrangemc = (TH2D *) mcGrid->Project(0,2);\r
346 TH1D *spectrum1Damc = (TH1D *) spectrum_zrangemc->ProjectionX("bin0",1,1);\r
347 TH1D *spectrum1Dbmc = (TH1D *) spectrum_zrangemc->ProjectionX("bin>0",2,12);\r
348 TH1D *spectrum1Dcmc = (TH1D *) spectrum_zrangemc->ProjectionX("total");\r
349 Double_t nbinbin0mc = spectrum1Damc->Integral();\r
350 Double_t nbinnobin0mc = spectrum1Dbmc->Integral();\r
351 Double_t nbintotalmc = spectrum1Dcmc->Integral();\r
352\r
353 Float_t numberofentriesmc = nbinnobin0mc+nbinbin0mc*0.880756;\r
354 printf("MC: Number in bin0 %f, number out %f, total %f, normalized %f\n",nbinbin0mc,nbinnobin0mc,nbintotalmc,numberofentriesmc);\r
355\r
356\r
357\r
358 /////////////////////////////\r
359 // Check number of events\r
360 /////////////////////////////\r
361\r
362 Int_t numberOfEventsafterallcuts = (Int_t) datahfecontainer->GetNumberOfEvents();\r
363\r
364 printf("Number of events not corrected %d\n",numberOfEventsafterallcuts);\r
365 printf("Number of events corrected %f\n",numberofentries);\r
366\r
367 AliHFEcontainer *mchfecontainer = (AliHFEcontainer *) resultsmc->FindObject("trackContainer");\r
368 if(!mchfecontainer) {\r
369 printf("No mc container \n");\r
370 return;\r
371 }\r
372\r
373 printf("Find the container V0\n");\r
374 AliHFEcontainer *containerhfeV0 = (AliHFEcontainer *) resultsdata->FindObject("containerV0");\r
375 if(!containerhfeV0) {\r
376 printf("No hfe container \n");\r
377 return;\r
378 }\r
379\r
380 // nonHFE backgrounds\r
381 TList *resultsmcbg = GetResults(fileMCbg,"_PID2");\r
382 if(!resultsmcbg){\r
383 printf("No output objects for mc: Calculation will terminate here\n"); return;\r
384 } \r
385 AliHFEcontainer *mcbghfecontainer = (AliHFEcontainer *) resultsmcbg->FindObject("trackContainer");\r
386 if(!mcbghfecontainer) {\r
387 printf("No mc container \n");\r
388 return;\r
389 }\r
390\r
391 AliCFContainer *containermcbg = (AliCFContainer *) resultsmcbg->FindObject("eventContainer");\r
392 if(!containermcbg) {\r
393 printf("No container \n");\r
394 return;\r
395 }\r
396\r
397 AliCFDataGrid *mcbgGrid = (AliCFDataGrid *) GetSpectrum(containermcbg,AliHFEcuts::kEventStepZRange);\r
398 TH2D *spectrum_zrangemcbg = (TH2D *) mcbgGrid->Project(0,2);\r
399 TH1D *spectrum1Damcbg = (TH1D *) spectrum_zrangemcbg->ProjectionX("bin0",1,1);\r
400 TH1D *spectrum1Dbmcbg = (TH1D *) spectrum_zrangemcbg->ProjectionX("bin>0",2,12);\r
401 TH1D *spectrum1Dcmcbg = (TH1D *) spectrum_zrangemcbg->ProjectionX("total");\r
402 Double_t nbinbin0mcbg = spectrum1Damcbg->Integral();\r
403 Double_t nbinnobin0mcbg = spectrum1Dbmcbg->Integral();\r
404 Double_t nbintotalmcbg = spectrum1Dcmcbg->Integral();\r
405\r
406 Float_t numberofentriesmcbg = nbinnobin0mcbg+nbinbin0mcbg*0.880756;\r
407 printf("MCbg: Number in bin0 %f, number out %f, total %f, normalized %f\n",nbinbin0mcbg,nbinnobin0mcbg,nbintotalmcbg,numberofentriesmcbg);\r
408\r
409\r
410 // hadron contamination after IP cuts\r
411 TList *hfeqa = GetQA(fileMC,"_PID2");\r
412 THnSparseF* hsHadronEffbyIPcut = (THnSparseF* )GetHadronEffbyIPcut(hfeqa);\r
413\r
414 //////////////\r
415 // Correct\r
416 /////////////\r
417\r
418 AliHFEspectrum *spectrum = new AliHFEspectrum("HFEspectrum");\r
419\r
420 spectrum->SetBeautyAnalysis();\r
421 // Enable background subtraction\r
422 spectrum->EnableIPanaHadronBgSubtract();\r
423 spectrum->EnableIPanaCharmBgSubtract();\r
424 spectrum->EnableIPanaConversionBgSubtract();\r
425 spectrum->EnableIPanaNonHFEBgSubtract();\r
426 //spectrum->SetNonHFEBackground2ndMethod();\r
427\r
428 spectrum->SetNbDimensions(1);\r
429 // If you want to correct positive (0) or negative (1) separately\r
430 //spectrum->SetChargeChoosen(0);\r
431 spectrum->SetDebugLevel(1);\r
432 spectrum->SetNumberOfEvents((Int_t)numberofentries);\r
433 spectrum->SetNumberOfMCEvents((Int_t)numberofentriesmc);\r
434 spectrum->SetNumberOfMC2Events((Int_t)numberofentriesmcbg);\r
435 spectrum->SetDumpToFile(smoothing);\r
436 // True step in MC (events in range +- 10 cm and no pile up)\r
437 spectrum->SetMCTruthStep(AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine);\r
438 // Step where we correct from MC (tracking + TOF)\r
439 spectrum->SetMCEffStep(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1);\r
440 // Step from data we correct for\r
441 spectrum->SetStepToCorrect(AliHFEcuts::kStepHFEcutsTRD + 3); // +3 = tracking + TOF + TPC + IPcut\r
442 // Steps to be corrected with V0 (TPC PID only)\r
443 spectrum->SetStepBeforeCutsV0(AliHFEcuts::kStepHFEcutsTRD + 1);\r
444 spectrum->SetStepAfterCutsV0(AliHFEcuts::kStepHFEcutsTRD + 2);\r
445\r
446 spectrum->SetHadronEffbyIPcut(hsHadronEffbyIPcut);\r
447 // Give everything (data container, mc container and V0 data container)\r
448 spectrum->Init(datahfecontainer,mchfecontainer,containerhfeV0,mcbghfecontainer);\r
449\r
450 // kTRUE means subtract hadron background, kFALSE means do not subtract hadron background\r
451 spectrum->CorrectBeauty(kTRUE);\r
452\r
453}\r
454\r
455//_____________________________________________________________________\r
456TList *GetResults(const Char_t *testfile,const Char_t *plus){\r
457 //\r
458 // read output\r
459 //\r
460 TFile *f = TFile::Open(testfile);\r
461 if(!f || f->IsZombie()){\r
462 printf("File not readable\n");\r
463 return 0x0;\r
464 }\r
465 TString name(plus);\r
466 printf("Name of TList %s\n",(const char*)name); \r
467 TList *l = dynamic_cast<TList *>(f->Get((const char*)name));\r
468 if(!l){\r
469 printf("Results list was not found\n");\r
470 f->Close(); delete f;\r
471 return 0x0;\r
472 } \r
473 TList *returnlist = dynamic_cast<TList *>(l->Clone());\r
474 f->Close(); delete f;\r
475 return returnlist;\r
476}\r
477\r
478//_____________________________________________________________________\r
479TList *GetQA(const Char_t *testfile,const Char_t *plus){\r
480 //\r
481 // read output\r
482 //\r
483 TFile *f = TFile::Open(testfile);\r
484 if(!f || f->IsZombie()){\r
485 printf("File not readable\n");\r
486 return 0x0;\r
487 }\r
488 TString name("HFE_QA");\r
489 name += plus;\r
490 printf("Name of TList %s\n",(const char*)name);\r
491 TList *l = dynamic_cast<TList *>(f->Get((const char*)name));\r
492 if(!l){\r
493 printf("QA list was not found\n");\r
494 f->Close(); delete f;\r
495 return 0x0;\r
496 }\r
497 TList *returnlist = dynamic_cast<TList *>(l->Clone());\r
498 f->Close(); delete f;\r
499 return returnlist;\r
500}\r
501\r
502//_____________________________________________________________________\r
503THnSparseF* GetHadronEffbyIPcut(TList *hfeqa){\r
504\r
505\r
506 // get hadron reduction factors due to the IP cuts\r
507 TList* tl=(TList*)hfeqa->FindObject("list_TaskQA");\r
508 TH1F* hbefore=(TH1F*)tl->FindObject("hadronsBeforeIPcut");\r
509 TH1F* hafter=(TH1F*)tl->FindObject("hadronsAfterIPcut");\r
510 TH1F* hreduction= (TH1F*)hafter->Clone("hreduction");\r
511 hbefore->Sumw2(); \r
512 hafter->Sumw2();\r
513 hreduction->Sumw2();\r
514 hreduction->Divide(hbefore);\r
515\r
516 Double_t* binEdges[0];\r
517 Int_t hnBin = hreduction->GetNbinsX();\r
518 Int_t nBin[1] = {hnBin};\r
519\r
520 for(int i=0; i<nBin[0]; i++){\r
521 binEdges[0][i] = hreduction->GetBinLowEdge(i+1);\r
522 }\r
523 binEdges[0][nBin[0]] = hreduction->GetBinLowEdge(nBin[0]) + hreduction->GetBinWidth(nBin[0]);\r
524\r
525 THnSparseF* hsreduction = new THnSparseF("hadroncontamin", "hadroncontamin; pt[GeV/c]", 1, nBin);\r
526 hsreduction->SetBinEdges(0, binEdges[0]);\r
527\r
528 Double_t dataE[1];\r
529 Double_t yval;\r
530 for(int i=0; i<nBin[0]; i++){\r
531 dataE[0]=hreduction->GetBinCenter(i+1);\r
532 yval=hreduction->GetBinContent(i+1);\r
533 hsreduction->Fill(dataE,yval);\r
534 }\r
535\r
536 Int_t* ibins;\r
537 ibins = new Int_t[nBin[0] + 1];\r
538 hsreduction->SetBinError(ibins,0);\r
539\r
540\r
541 return hsreduction;\r
542}\r
543\r
544//_________________________________________________________________________\r
545TObject* GetSpectrum(AliCFContainer *c, Int_t step) {\r
546 AliCFDataGrid* data = new AliCFDataGrid("data","",*c,step);\r
547 //data->SetMeasured(step);\r
548 return data;\r
549}\r
550//___________________________________________________________________________\r
551void CorrectFromTheWidth(TH1D *h1) {\r
552 //\r
553 // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}\r
554 //\r
555\r
556 TAxis *axis = h1->GetXaxis();\r
557 Int_t nbinX = h1->GetNbinsX();\r
558\r
559 for(Int_t i = 1; i <= nbinX; i++) {\r
560\r
561 Double_t width = axis->GetBinWidth(i);\r
562 Double_t content = h1->GetBinContent(i);\r
563 Double_t error = h1->GetBinError(i); \r
564 h1->SetBinContent(i,content/width); \r
565 h1->SetBinError(i,error/width);\r
566 }\r
567\r
568}\r