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