PWG3/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);\r
   2 void CorrectForEfficiencyBeauty(const char *filedata,const char *fileMC, const char *fileMCbg);\r
   3 TList *GetResults(const Char_t *testfile,const Char_t *plus="");\r
   4 TList *GetQA(const Char_t *testfile,const Char_t *plus="");\r
   5 TObject* GetSpectrum(AliCFContainer *c, Int_t step);\r
   6 THnSparseF* GetHadronEffbyIPcut(TList *hfeqa);\r
   7 void CorrectFromTheWidth(TH1D *h1);\r
   8 \r
   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) {\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
 269 void 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
 456 TList *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
 479 TList *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
 503 THnSparseF* 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
 545 TObject* 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
 551 void 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