Update hfe package
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEspectrum.cxx
1
2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 *                                                                        *
5 * Author: The ALICE Off-line Project.                                    *
6 * Contributors are mentioned in the code where appropriate.              *
7 *                                                                        *
8 * Permission to use, copy, modify and distribute this software and its   *
9 * documentation strictly for non-commercial purposes is hereby granted   *
10 * without fee, provided that the above copyright notice appears in all   *
11 * copies and that both the copyright notice and this permission notice   *
12 * appear in the supporting documentation. The authors make no claims     *
13 * about the suitability of this software for any purpose. It is          *
14 * provided "as is" without express or implied warranty.                  *
15 **************************************************************************/
16 //
17 // Class for spectrum correction
18 // Subtraction of hadronic background, Unfolding of the data and
19 // Renormalization done here
20 // The following containers have to be set:
21 //  - Correction framework container for real data
22 //  - Correction framework container for MC (Efficiency Map)
23 //  - Correction framework container for background coming from data
24 //  - Correction framework container for background coming from MC
25 //
26 //  Author: 
27 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
28 //            Markus Fasel <M.Fasel@gsi.de>
29 //
30
31 #include <TArrayD.h>
32 #include <TH1.h>
33 #include <TList.h>
34 #include <TObjArray.h>
35 #include <TROOT.h>
36 #include <TCanvas.h>
37 #include <TLegend.h>
38 #include <TStyle.h>
39 #include <TMath.h>
40 #include <TAxis.h>
41 #include <TGraphErrors.h>
42 #include <TFile.h>
43 #include <TPad.h>
44 #include <TH2D.h>
45 #include <TF1.h>
46
47 #include "AliPID.h"
48 #include "AliCFContainer.h"
49 #include "AliCFDataGrid.h"
50 #include "AliCFEffGrid.h"
51 #include "AliCFGridSparse.h"
52 #include "AliCFUnfolding.h"
53 #include "AliLog.h"
54
55 #include "AliHFEspectrum.h"
56 #include "AliHFEcuts.h"
57 #include "AliHFEcontainer.h"
58 #include "AliHFEtools.h"
59
60 ClassImp(AliHFEspectrum)
61
62 //____________________________________________________________
63 AliHFEspectrum::AliHFEspectrum(const char *name):
64   TNamed(name, ""),
65   fCFContainers(new TObjArray(kDataContainerV0+1)),
66   fTemporaryObjects(NULL),
67   fCorrelation(NULL),
68   fBackground(NULL),
69   fEfficiencyFunction(NULL),
70   fWeightCharm(NULL),
71   fInclusiveSpectrum(kTRUE),
72   fDumpToFile(kFALSE),
73   fEtaSelected(kFALSE),
74   fUnSetCorrelatedErrors(kTRUE),
75   fSetSmoothing(kFALSE),
76   fIPanaHadronBgSubtract(kFALSE),
77   fIPanaCharmBgSubtract(kFALSE),
78   fIPanaConversionBgSubtract(kFALSE),
79   fIPanaNonHFEBgSubtract(kFALSE),
80   fIPParameterizedEff(kFALSE),
81   fNonHFEsyst(0),
82   fBeauty2ndMethod(kFALSE),
83   fIPEffCombinedSamples(kTRUE),
84   fNbDimensions(1),
85   fNMCEvents(0),
86   fStepMC(-1),
87   fStepTrue(-1),
88   fStepData(-1),
89   fStepBeforeCutsV0(-1),
90   fStepAfterCutsV0(-1),
91   fStepGuessedUnfolding(-1),
92   fNumberOfIterations(5),
93   fChargeChoosen(kAllCharge),
94   fNCentralityBinAtTheEnd(0),
95   fTestCentralityLow(-1),
96   fTestCentralityHigh(-1),
97   fFillMoreCorrelationMatrix(kFALSE),
98   fHadronEffbyIPcut(NULL),
99   fConversionEffbgc(NULL),
100   fNonHFEEffbgc(NULL),      
101   fBSpectrum2ndMethod(NULL),
102   fkBeauty2ndMethodfilename(""),
103   fBeamType(0),
104   fEtaSyst(kTRUE),
105   fDebugLevel(0),
106   fWriteToFile(kFALSE),
107   fUnfoldBG(kFALSE)
108 {
109   //
110   // Default constructor
111   //
112
113   for(Int_t k = 0; k < 20; k++){
114       fNEvents[k] = 0;
115       fNMCbgEvents[k] = 0;
116       fLowBoundaryCentralityBinAtTheEnd[k] = 0;
117       fHighBoundaryCentralityBinAtTheEnd[k] = 0;
118       if(k<kCentrality)
119       {
120           fEfficiencyTOFPIDD[k] = 0;
121           fEfficiencyesdTOFPIDD[k] = 0;
122           fEfficiencyIPCharmD[k] = 0;     
123           fEfficiencyIPBeautyD[k] = 0;    
124           fEfficiencyIPBeautyesdD[k] = 0;
125           fEfficiencyIPConversionD[k] = 0;
126           fEfficiencyIPNonhfeD[k] = 0;   
127
128           fConversionEff[k] = 0;
129           fNonHFEEff[k] = 0;
130           fCharmEff[k] = 0;
131           fBeautyEff[k] = 0;
132           fEfficiencyCharmSigD[k] = 0;
133           fEfficiencyBeautySigD[k] = 0;
134           fEfficiencyBeautySigesdD[k] = 0;
135       }
136   }
137   memset(fEtaRange, 0, sizeof(Double_t) * 2);
138   memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
139   memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
140   memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
141 }
142
143 //____________________________________________________________
144 AliHFEspectrum::~AliHFEspectrum(){
145   //
146   // Destructor
147   //
148   if(fCFContainers) delete fCFContainers;
149   if(fTemporaryObjects){
150     fTemporaryObjects->Clear();
151     delete fTemporaryObjects;
152   }
153 }
154 //____________________________________________________________
155 Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
156   //
157   // Init what we need for the correction:
158   //
159   // Raw spectrum, hadron contamination
160   // MC efficiency maps, correlation matrix
161   // V0 efficiency if wanted
162   //
163   // This for a given dimension.
164   // If no inclusive spectrum, then take only efficiency map for beauty electron
165   // and the appropriate correlation matrix
166   //
167
168   
169   if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
170
171   Int_t kNdim = 3;
172   Int_t kNcentr =1;
173   Int_t ptpr =0;
174   if(fBeamType==0) kNdim=3;
175   if(fBeamType==1)
176   {
177       kNdim=4;
178       kNcentr=11;
179       ptpr=1;
180   }
181
182   Int_t dims[kNdim];
183   // Get the requested format
184   if(fBeamType==0){
185     // pp analysis
186     switch(fNbDimensions){
187       case 1:   dims[0] = 0;
188                 break;
189       case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
190                 break;
191       case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
192                 break;
193       default:
194                 AliError("Container with this number of dimensions not foreseen (yet)");
195                 return kFALSE;
196     };
197   }
198
199   if(fBeamType==1){
200     // PbPb analysis; centrality as first dimension
201     Int_t nbDimensions = fNbDimensions;
202     fNbDimensions = fNbDimensions + 1;
203     switch(nbDimensions){
204       case 1: dims[0] = 5;
205               dims[1] = 0;
206               break;
207       case 2: dims[0] = 5;
208               for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
209               break;
210       case 3: dims[0] = 5;
211               for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
212               break;
213       default:
214               AliError("Container with this number of dimensions not foreseen (yet)");
215               return kFALSE;
216     };
217   }
218
219   // Data container: raw spectrum + hadron contamination  
220   AliCFContainer *datacontainer = 0x0; 
221   if(fInclusiveSpectrum) {
222     datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
223   }
224   else{
225
226       datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
227   }
228   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
229   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
230
231   AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
232   AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
233   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
234
235   SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
236   SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
237
238   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps 
239   AliCFContainer *mccontaineresd = 0x0;
240   AliCFContainer *mccontaineresdbg = 0x0;
241   AliCFContainer *mccontainermc = 0x0;
242   AliCFContainer *mccontainermcbg = 0x0;
243   AliCFContainer *nonHFEweightedContainer = 0x0;
244   AliCFContainer *convweightedContainer = 0x0;
245   AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
246   AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
247    
248   if(fInclusiveSpectrum) {
249     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
250     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
251   }
252   else {
253     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
254     mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
255     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
256     mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
257
258     if(fNonHFEsyst){   
259       const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
260       const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
261       for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
262         for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
263           nonHFEtempContainer =  bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
264           convtempContainer =  bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
265           for(Int_t icentr=0;icentr<kNcentr;icentr++)
266           {
267               if(fBeamType==0)
268               {
269                   fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
270                   fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
271               }
272               if(fBeamType==1)
273               {
274               
275                 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
276                 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
277               }
278 //            if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
279           }
280           if(fBeamType == 1)break;
281         }
282       }
283     }
284     // else{      
285       nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
286       convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
287       if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;  
288       //}
289   }
290   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;  
291   
292   Int_t source = -1;
293   if(!fInclusiveSpectrum) source = 1; //beauty
294   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
295   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
296   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
297   SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
298   SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
299
300   // set charm, nonHFE container to estimate BG
301   if(!fInclusiveSpectrum){
302    source = 0; //charm
303    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
304    SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
305
306    //if(!fNonHFEsyst){
307      AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
308      SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
309      AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
310      SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
311      //}
312
313      SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
314
315   }
316   // MC container: correlation matrix
317   THnSparseF *mccorrelation = 0x0;
318   if(fInclusiveSpectrum) {
319     if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
320     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
321     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
322     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
323     else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
324     
325     if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
326   }
327   else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
328   //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
329   if(!mccorrelation) return kFALSE;
330   Int_t testCentralityLow = fTestCentralityLow;
331   Int_t testCentralityHigh = fTestCentralityHigh;
332   if(fFillMoreCorrelationMatrix) {
333     testCentralityLow = fTestCentralityLow-1;
334     testCentralityHigh = fTestCentralityHigh+1;
335   }
336   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh);
337   if(!mccorrelationD) {
338     printf("No correlation\n");
339     return kFALSE;
340   }
341   SetCorrelation(mccorrelationD);
342
343   // V0 container Electron, pt eta phi
344   if(v0hfecontainer) {
345     AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
346     if(!containerV0) return kFALSE;
347     AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
348     if(!containerV0Electron) return kFALSE;
349     SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
350   }
351  
352
353   if(fDebugLevel>0){
354    
355     AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
356     contaminationspectrum->SetName("contaminationspectrum");
357     TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
358     ccontaminationspectrum->Divide(3,1);
359     ccontaminationspectrum->cd(1);
360     TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
361     TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
362     TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
363     contaminationspectrum2dpteta->SetStats(0);
364     contaminationspectrum2dpteta->SetTitle("");
365     contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
366     contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
367     contaminationspectrum2dptphi->SetStats(0);
368     contaminationspectrum2dptphi->SetTitle("");
369     contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
370     contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
371     contaminationspectrum2detaphi->SetStats(0);
372     contaminationspectrum2detaphi->SetTitle("");
373     contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
374     contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
375     contaminationspectrum2dptphi->Draw("colz");
376     ccontaminationspectrum->cd(2);
377     contaminationspectrum2dpteta->Draw("colz");
378     ccontaminationspectrum->cd(3);
379     contaminationspectrum2detaphi->Draw("colz");
380
381     TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
382     ccorrelation->cd(1);
383     if(mccorrelationD) {
384       if(fBeamType==0){
385         TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
386         ptcorrelation->Draw("colz");
387       }
388       if(fBeamType==1){
389         TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
390         ptcorrelation->Draw("colz");
391       }
392     }
393     if(fWriteToFile) ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
394   }
395
396   TFile *file = TFile::Open("tests.root","recreate");
397   datacontainerD->Write("data");
398   mccontainermcD->Write("mcefficiency");
399   mccorrelationD->Write("correlationmatrix");
400   file->Close();
401
402   return kTRUE;
403 }
404
405
406 //____________________________________________________________
407 void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
408   //
409   // get spectrum for beauty 2nd method
410   //
411   //
412     TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
413     fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
414 }
415
416
417 //____________________________________________________________
418 Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
419   //
420   // Correct the spectrum for efficiency and unfolding
421   // with both method and compare
422   //
423   
424   gStyle->SetPalette(1);
425   gStyle->SetOptStat(1111);
426   gStyle->SetPadBorderMode(0);
427   gStyle->SetCanvasColor(10);
428   gStyle->SetPadLeftMargin(0.13);
429   gStyle->SetPadRightMargin(0.13);
430
431   printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
432
433   ///////////////////////////
434   // Check initialization
435   ///////////////////////////
436
437   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
438     AliInfo("You have to init before");
439     return kFALSE;
440   }
441   
442   if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
443     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
444     return kFALSE;
445   }
446  
447   SetNumberOfIteration(10);
448   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
449     
450   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
451   //////////////////////////////////
452   // Subtract hadron background
453   /////////////////////////////////
454   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
455   if(subtractcontamination) {
456     dataspectrumaftersubstraction = SubtractBackground(kTRUE);
457     dataGridAfterFirstSteps = dataspectrumaftersubstraction;
458   }
459
460   ////////////////////////////////////////////////
461   // Correct for TPC efficiency from V0
462   ///////////////////////////////////////////////
463   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
464   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
465   if(dataContainerV0){printf("Got the V0 container\n");
466     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
467     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;  
468   }
469
470   //////////////////////////////////////////////////////////////////////////////
471   // Correct for efficiency parametrized (if TPC efficiency is parametrized)
472   /////////////////////////////////////////////////////////////////////////////
473   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
474   if(fEfficiencyFunction){
475     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
476     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;  
477   }
478     
479   ///////////////
480   // Unfold
481   //////////////
482   TList *listunfolded = Unfold(dataGridAfterFirstSteps);
483   if(!listunfolded){
484     printf("Unfolded failed\n");
485     return kFALSE;
486   }
487   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
488   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
489   if(!correctedspectrum){
490     AliError("No corrected spectrum\n");
491     return kFALSE;
492   }
493   if(!residualspectrum){
494     AliError("No residul spectrum\n");
495     return kFALSE;
496   }   
497   
498   /////////////////////
499   // Simply correct
500   ////////////////////
501   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
502   
503
504   //////////
505   // Plot
506   //////////
507   if(fDebugLevel > 0.0) {
508
509     Int_t ptpr = 0;
510     if(fBeamType==0) ptpr=0;
511     if(fBeamType==1) ptpr=1;
512   
513     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
514     ccorrected->Divide(2,1);
515     ccorrected->cd(1);
516     gPad->SetLogy();
517     TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
518     correctedspectrumD->SetTitle("");
519     correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
520     correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
521     correctedspectrumD->SetMarkerStyle(26);
522     correctedspectrumD->SetMarkerColor(kBlue);
523     correctedspectrumD->SetLineColor(kBlue);
524     correctedspectrumD->Draw("AP");
525     TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
526     alltogetherspectrumD->SetTitle("");
527     alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
528     alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
529     alltogetherspectrumD->SetMarkerStyle(25);
530     alltogetherspectrumD->SetMarkerColor(kBlack);
531     alltogetherspectrumD->SetLineColor(kBlack);
532     alltogetherspectrumD->Draw("P");
533     TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
534     legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
535     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
536     legcorrected->Draw("same");
537     ccorrected->cd(2);
538     TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
539     TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
540     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
541     ratiocorrected->SetName("ratiocorrected");
542     ratiocorrected->SetTitle("");
543     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
544     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
545     ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
546     ratiocorrected->SetStats(0);
547     ratiocorrected->Draw();
548     if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps");
549
550     //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
551     //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
552     //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
553     //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
554
555     TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
556     TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
557     TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
558     TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
559
560     if(fBeamType==1) {
561
562       TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
563       ccorrectedallspectra->Divide(2,1);
564       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
565       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
566       
567       THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
568       TAxis *cenaxisa = sparsesu->GetAxis(0);
569       THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
570       TAxis *cenaxisb = sparsed->GetAxis(0);
571       Int_t nbbin = cenaxisb->GetNbins();
572       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
573       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
574       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
575         TString titlee("corrected_centrality_bin_");
576         titlee += "[";
577         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
578         titlee += "_";
579         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
580         titlee += "[";
581         TString titleec("corrected_check_projection_bin_");
582         titleec += "[";
583         titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
584         titleec += "_";
585         titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
586         titleec += "[";
587         TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
588         titleenameunotnormalized += "[";
589         titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
590         titleenameunotnormalized += "_";
591         titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
592         titleenameunotnormalized += "[";
593         TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
594         titleenameunormalized += "[";
595         titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
596         titleenameunormalized += "_";
597         titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];      
598         titleenameunormalized += "[";
599         TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
600         titleenamednotnormalized += "[";
601         titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
602         titleenamednotnormalized += "_";
603         titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
604         titleenamednotnormalized += "[";
605         TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
606         titleenamednormalized += "[";
607         titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
608         titleenamednormalized += "_";
609         titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];      
610         titleenamednormalized += "[";
611         Int_t nbEvents = 0;
612         for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
613           printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
614           nbEvents += fNEvents[k];
615         }
616         Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
617         Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
618         printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
619         Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
620         Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
621         printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
622         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
623         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
624         TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
625         ccorrectedcheck->cd(1);
626         TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
627         TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
628         aftersuc->Draw();
629         aftersdc->Draw("same");
630         TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
631         ccorrectede->Divide(2,1);
632         ccorrectede->cd(1);
633         gPad->SetLogy();
634         TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
635         CorrectFromTheWidth(aftersu);
636         aftersu->SetName((const char*)titleenameunotnormalized);
637         unfoldingspectrac[binc] = *aftersu;
638         ccorrectede->cd(1);
639         TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
640         aftersun->SetTitle("");
641         aftersun->GetYaxis()->SetTitleOffset(1.5);
642         aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
643         aftersun->SetMarkerStyle(26);
644         aftersun->SetMarkerColor(kBlue);
645         aftersun->SetLineColor(kBlue);
646         aftersun->Draw("AP");
647         aftersun->SetName((const char*)titleenameunormalized);
648         unfoldingspectracn[binc] = *aftersun;
649         ccorrectede->cd(1);
650         TH1D *aftersd = (TH1D *) sparsed->Projection(1);
651         CorrectFromTheWidth(aftersd);
652         aftersd->SetName((const char*)titleenamednotnormalized);
653         correctedspectrac[binc] = *aftersd;
654         ccorrectede->cd(1);
655         TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
656         aftersdn->SetTitle("");
657         aftersdn->GetYaxis()->SetTitleOffset(1.5);
658         aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
659         aftersdn->SetMarkerStyle(25);
660         aftersdn->SetMarkerColor(kBlack);
661         aftersdn->SetLineColor(kBlack);
662         aftersdn->Draw("P");
663         aftersdn->SetName((const char*)titleenamednormalized);
664         correctedspectracn[binc] = *aftersdn;
665         TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
666         legcorrectedud->AddEntry(aftersun,"Corrected","p");
667         legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
668         legcorrectedud->Draw("same");
669         ccorrectedallspectra->cd(1);
670         gPad->SetLogy();
671         TH1D *aftersunn = (TH1D *) aftersun->Clone();
672         aftersunn->SetMarkerStyle(stylee[binc]);
673         aftersunn->SetMarkerColor(colorr[binc]);
674         if(binc==0) aftersunn->Draw("AP");
675         else aftersunn->Draw("P");
676         legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
677         ccorrectedallspectra->cd(2);
678         gPad->SetLogy();
679         TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
680         aftersdnn->SetMarkerStyle(stylee[binc]);
681         aftersdnn->SetMarkerColor(colorr[binc]);
682         if(binc==0) aftersdnn->Draw("AP");
683         else aftersdnn->Draw("P");
684         legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
685         ccorrectede->cd(2);
686         TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
687         TString titleee("ratiocorrected_bin_");
688         titleee += binc;
689         ratiocorrectedbinc->SetName((const char*) titleee);
690         ratiocorrectedbinc->SetTitle("");
691         ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
692         ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
693         ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
694         ratiocorrectedbinc->SetStats(0);
695         ratiocorrectedbinc->Draw();
696       }
697
698       ccorrectedallspectra->cd(1);
699       legtotal->Draw("same");
700       ccorrectedallspectra->cd(2);
701       legtotalg->Draw("same");
702       
703       cenaxisa->SetRange(0,nbbin);
704       cenaxisb->SetRange(0,nbbin);
705       if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps");
706     }
707
708     // Dump to file if needed
709     if(fDumpToFile) {
710       TFile *out = new TFile("finalSpectrum.root","recreate");
711       correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
712       correctedspectrumD->Write();
713       alltogetherspectrumD->SetName("AlltogetherSpectrum");
714       alltogetherspectrumD->Write();
715       ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
716       ratiocorrected->Write();
717       correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
718       correctedspectrum->Write();
719       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
720       alltogetherCorrection->Write();
721       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
722               unfoldingspectrac[binc].Write();
723               unfoldingspectracn[binc].Write();
724               correctedspectrac[binc].Write();
725               correctedspectracn[binc].Write();
726       }
727       out->Close(); delete out;
728     }
729
730     if (unfoldingspectrac) delete[] unfoldingspectrac;
731     if (unfoldingspectracn)  delete[] unfoldingspectracn;
732     if (correctedspectrac) delete[] correctedspectrac;
733     if (correctedspectracn) delete[] correctedspectracn;
734
735   }
736
737   return kTRUE;
738 }
739
740 //____________________________________________________________
741 Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
742   //
743   // Correct the spectrum for efficiency and unfolding for beauty analysis
744   // with both method and compare
745   //
746   
747   gStyle->SetPalette(1);
748   gStyle->SetOptStat(1111);
749   gStyle->SetPadBorderMode(0);
750   gStyle->SetCanvasColor(10);
751   gStyle->SetPadLeftMargin(0.13);
752   gStyle->SetPadRightMargin(0.13);
753
754   ///////////////////////////
755   // Check initialization
756   ///////////////////////////
757
758   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
759     AliInfo("You have to init before");
760     return kFALSE;
761   }
762   
763   if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
764     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
765     return kFALSE;
766   }
767  
768   SetNumberOfIteration(10);
769   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
770     
771   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
772   //////////////////////////////////
773   // Subtract hadron background
774   /////////////////////////////////
775   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
776   AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
777   TGraphErrors *gNormalizedRawSpectrum = 0x0;
778   if(subtractcontamination) {
779       if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
780       else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
781       unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
782       dataGridAfterFirstSteps = dataspectrumaftersubstraction;
783       gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
784   }
785
786   printf("after normalize getting IP \n");
787
788   /////////////////////////////////////////////////////////////////////////////////////////
789   // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
790   /////////////////////////////////////////////////////////////////////////////////////////
791
792   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
793   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
794   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
795
796   if(fEfficiencyFunction){
797     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
798     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
799   }
800   else if(dataContainerV0){
801     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
802     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;  
803   }  
804  
805
806
807   ///////////////
808   // Unfold
809   //////////////
810   TList *listunfolded = Unfold(dataGridAfterFirstSteps);
811   if(!listunfolded){
812     printf("Unfolded failed\n");
813     return kFALSE;
814   }
815   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
816   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
817   if(!correctedspectrum){
818     AliError("No corrected spectrum\n");
819     return kFALSE;
820   }
821   if(!residualspectrum){
822     AliError("No residual spectrum\n");
823     return kFALSE;
824   }   
825   
826   /////////////////////
827   // Simply correct
828   ////////////////////
829
830   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
831   
832
833   //////////
834   // Plot
835   //////////
836
837   if(fDebugLevel > 0.0) {
838
839     Int_t ptpr = 0;
840     if(fBeamType==0) ptpr=0;
841     if(fBeamType==1) ptpr=1;
842   
843     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
844     ccorrected->Divide(2,1);
845     ccorrected->cd(1);
846     gPad->SetLogy();
847     TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
848     correctedspectrumD->SetTitle("");
849     correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
850     correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
851     correctedspectrumD->SetMarkerStyle(26);
852     correctedspectrumD->SetMarkerColor(kBlue);
853     correctedspectrumD->SetLineColor(kBlue);
854     correctedspectrumD->Draw("AP");
855     TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
856     alltogetherspectrumD->SetTitle("");
857     alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
858     alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
859     alltogetherspectrumD->SetMarkerStyle(25);
860     alltogetherspectrumD->SetMarkerColor(kBlack);
861     alltogetherspectrumD->SetLineColor(kBlack);
862     alltogetherspectrumD->Draw("P");
863     TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
864     legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
865     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
866     legcorrected->Draw("same");
867     ccorrected->cd(2);
868     TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
869     TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
870     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
871     ratiocorrected->SetName("ratiocorrected");
872     ratiocorrected->SetTitle("");
873     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
874     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
875     ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
876     ratiocorrected->SetStats(0);
877     ratiocorrected->Draw();
878     if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
879
880     if(fBeamType == 0){
881         if(fNonHFEsyst){
882             CalculateNonHFEsyst(0);
883         }
884     }
885
886     // Dump to file if needed
887
888     if(fDumpToFile) {
889         // to do centrality dependent
890
891       TFile *out;
892       out = new TFile("finalSpectrum.root","recreate");
893       out->cd();
894       //
895       correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
896       correctedspectrumD->Write();
897       alltogetherspectrumD->SetName("AlltogetherSpectrum");
898       alltogetherspectrumD->Write();
899       ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
900       ratiocorrected->Write();
901       //
902       correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
903       correctedspectrum->Write();
904       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
905       alltogetherCorrection->Write();
906       //
907       if(unnormalizedRawSpectrum) {
908       unnormalizedRawSpectrum->SetName("beautyAfterIP");
909       unnormalizedRawSpectrum->Write();
910       }
911
912       if(gNormalizedRawSpectrum){
913         gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
914         gNormalizedRawSpectrum->Write();
915       }
916
917       if(fBeamType==0) {
918           Int_t countpp=0;
919           fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
920           fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
921           fEfficiencyCharmSigD[countpp]->Write();
922           fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
923           fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
924           fEfficiencyBeautySigD[countpp]->Write();
925           fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
926           fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
927           fCharmEff[countpp]->Write();
928           fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
929           fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
930           fBeautyEff[countpp]->Write();
931           fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
932           fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
933           fConversionEff[countpp]->Write();
934           fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
935           fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
936           fNonHFEEff[countpp]->Write();
937       }
938
939       if(fBeamType==1) {
940
941           TGraphErrors* correctedspectrumDc[kCentrality];
942           TGraphErrors* alltogetherspectrumDc[kCentrality];
943           for(Int_t i=0;i<kCentrality-2;i++)
944           {
945               correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
946               correctedspectrumDc[i] = Normalize(correctedspectrum,i);
947               if(correctedspectrumDc[i]){
948                 correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
949                 correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
950                 correctedspectrumDc[i]->Write();
951               }
952               alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
953               alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
954               if(alltogetherspectrumDc[i]){
955                 alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
956                 alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
957                 alltogetherspectrumDc[i]->Write();
958               }
959             
960               TH1D *centrcrosscheck = correctedspectrum->Projection(0);
961               centrcrosscheck->SetTitle(Form("centrality_%i",i));
962               centrcrosscheck->SetName(Form("centrality_%i",i));
963               centrcrosscheck->Write();
964
965               TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
966               TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
967
968               TH1D *centrcrosscheck2 =  (TH1D *) alltogetherCorrection->Project(0);
969               centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
970               centrcrosscheck2->SetName(Form("centrality2_%i",i));
971               centrcrosscheck2->Write();
972
973               TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
974               ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
975               ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
976               ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
977               ratiocorrectedc->Write();
978
979               fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
980               fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
981               fEfficiencyCharmSigD[i]->Write();
982               fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
983               fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
984               fEfficiencyBeautySigD[i]->Write();
985               fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
986               fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
987               fCharmEff[i]->Write();
988               fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
989               fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
990               fBeautyEff[i]->Write();
991               fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
992               fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
993               fConversionEff[i]->Write();
994               fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
995               fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
996               fNonHFEEff[i]->Write();
997           }
998
999       }
1000
1001       out->Close(); 
1002       delete out;
1003     }
1004   }
1005
1006   return kTRUE;
1007 }
1008
1009 //____________________________________________________________
1010 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
1011   //
1012   // Apply background subtraction
1013   //
1014
1015     Int_t ptpr = 0;
1016     Int_t nbins=1;
1017     if(fBeamType==0)
1018     {
1019         ptpr=0;
1020         nbins=1;
1021     }
1022     if(fBeamType==1)
1023     {
1024         ptpr=1;
1025         nbins=2;
1026     }
1027
1028   // Raw spectrum
1029   AliCFContainer *dataContainer = GetContainer(kDataContainer);
1030   if(!dataContainer){
1031     AliError("Data Container not available");
1032     return NULL;
1033   }
1034   printf("Step data: %d\n",fStepData);
1035   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
1036
1037   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
1038   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); 
1039  
1040
1041   // Background Estimate
1042   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
1043   if(!backgroundContainer){
1044     AliError("MC background container not found");
1045     return NULL;
1046   }
1047  
1048   Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
1049   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
1050
1051   if(!fInclusiveSpectrum){
1052     //Background subtraction for IP analysis
1053
1054     TH1D *incElecCent[kCentrality-1];
1055     TH1D *charmCent[kCentrality-1];
1056     TH1D *convCent[kCentrality-1];
1057     TH1D *nonHFECent[kCentrality-1]; 
1058     TH1D *subtractedCent[kCentrality-1];
1059     TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
1060     CorrectFromTheWidth(measuredTH1Draw);
1061     if(fBeamType==1){
1062       THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1063       for(Int_t icent =  1; icent < kCentrality-1; icent++){
1064         sparseIncElec->GetAxis(0)->SetRange(icent,icent);
1065         incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr);
1066         CorrectFromTheWidth(incElecCent[icent-1]);
1067       }
1068     }
1069     TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400);
1070     rawspectra->cd();
1071     rawspectra->SetLogy();
1072     gStyle->SetOptStat(0);
1073     TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
1074     measuredTH1Draw->SetMarkerStyle(20);
1075     measuredTH1Draw->Draw();
1076     measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9);
1077     lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
1078     TH1D* htemp;
1079     Int_t* bins=new Int_t[2];
1080     if(fIPanaHadronBgSubtract){
1081       // Hadron background
1082         printf("Hadron background for IP analysis subtracted!\n");
1083         if(fBeamType==0)
1084         {
1085
1086             htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
1087             bins[0]=htemp->GetNbinsX();
1088         }
1089         if(fBeamType==1)
1090         {
1091             htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
1092             bins[0]=htemp->GetNbinsX();
1093             htemp  = (TH1D *) fHadronEffbyIPcut->Projection(1);
1094             bins[1]=htemp->GetNbinsX();
1095         }
1096       AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
1097       hbgContainer->SetGrid(fHadronEffbyIPcut);
1098       backgroundGrid->Multiply(hbgContainer,1);
1099       // draw raw hadron bg spectra
1100       TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
1101       CorrectFromTheWidth(hadronbg);
1102       hadronbg->SetMarkerColor(7);
1103       hadronbg->SetMarkerStyle(20);
1104       rawspectra->cd();
1105       hadronbg->Draw("samep");
1106       lRaw->AddEntry(hadronbg,"hadrons");
1107       // subtract hadron contamination
1108       spectrumSubtracted->Add(backgroundGrid,-1.0);
1109     }
1110     if(fIPanaCharmBgSubtract){
1111       // Charm background
1112       printf("Charm background for IP analysis subtracted!\n");
1113       AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
1114       // draw charm bg spectra
1115       TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
1116       CorrectFromTheWidth(charmbg);
1117       charmbg->SetMarkerColor(3);
1118       charmbg->SetMarkerStyle(20);
1119       rawspectra->cd();
1120       charmbg->Draw("samep");
1121       lRaw->AddEntry(charmbg,"charm elecs");
1122       // subtract charm background
1123       spectrumSubtracted->Add(charmbgContainer,-1.0);
1124       if(fBeamType==1){
1125         THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
1126         for(Int_t icent =  1; icent < kCentrality-1; icent++){ 
1127           sparseCharmElec->GetAxis(0)->SetRange(icent,icent);
1128           charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr);
1129           CorrectFromTheWidth(charmCent[icent-1]);
1130         }
1131       }
1132     }
1133     if(fIPanaConversionBgSubtract){
1134       // Conversion background
1135       AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground(); 
1136       // draw conversion bg spectra
1137       TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
1138       CorrectFromTheWidth(conversionbg);
1139       conversionbg->SetMarkerColor(4);
1140       conversionbg->SetMarkerStyle(20);
1141       rawspectra->cd();
1142       conversionbg->Draw("samep");
1143       lRaw->AddEntry(conversionbg,"conversion elecs");
1144       // subtract conversion background
1145       spectrumSubtracted->Add(conversionbgContainer,-1.0);
1146       if(fBeamType==1){
1147         THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid();
1148         for(Int_t icent =  1; icent < kCentrality-1; icent++){
1149           sparseconvElec->GetAxis(0)->SetRange(icent,icent);
1150           convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr);
1151           CorrectFromTheWidth(convCent[icent-1]);
1152         }
1153       }
1154     }
1155     if(fIPanaNonHFEBgSubtract){
1156       // NonHFE background
1157       AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
1158       // draw Dalitz/dielectron bg spectra
1159       TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
1160       CorrectFromTheWidth(nonhfebg);
1161       nonhfebg->SetMarkerColor(6);
1162       nonhfebg->SetMarkerStyle(20);
1163       rawspectra->cd();
1164       nonhfebg->Draw("samep");
1165       lRaw->AddEntry(nonhfebg,"non-HF elecs");
1166       // subtract Dalitz/dielectron background
1167       spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
1168       if(fBeamType==1){
1169         THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid();
1170         for(Int_t icent =  1; icent < kCentrality-1; icent++){
1171           sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent);
1172           nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr);
1173           CorrectFromTheWidth(nonHFECent[icent-1]);
1174         }
1175       }
1176     }
1177
1178     TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
1179     CorrectFromTheWidth(rawbgsubtracted);
1180     rawbgsubtracted->SetMarkerStyle(24);
1181     rawspectra->cd();
1182     lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
1183     rawbgsubtracted->Draw("samep");
1184     lRaw->Draw("SAME");
1185     gPad->SetGrid();
1186     //rawspectra->SaveAs("rawspectra.eps");
1187     
1188     if(fBeamType==1){
1189       THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1190       for(Int_t icent =  1; icent < kCentrality-1; icent++){
1191         sparseSubtracted->GetAxis(0)->SetRange(icent,icent);
1192         subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr);
1193         CorrectFromTheWidth(subtractedCent[icent-1]);
1194       }
1195     
1196       TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85);
1197       TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800);
1198       centRaw->Divide(3,3);
1199       for(Int_t icent = 1; icent < kCentrality-1; icent++){
1200         centRaw->cd(icent);
1201         gPad->SetLogx();
1202         gPad->SetLogy();
1203         incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.);
1204         incElecCent[icent-1]->Draw("p");
1205         incElecCent[icent-1]->SetMarkerColor(1);
1206         incElecCent[icent-1]->SetMarkerStyle(20);
1207         charmCent[icent-1]->Draw("samep");
1208         charmCent[icent-1]->SetMarkerColor(3);
1209         charmCent[icent-1]->SetMarkerStyle(20);
1210         convCent[icent-1]->Draw("samep");
1211         convCent[icent-1]->SetMarkerColor(4);
1212         convCent[icent-1]->SetMarkerStyle(20);
1213         nonHFECent[icent-1]->Draw("samep");
1214         nonHFECent[icent-1]->SetMarkerColor(6);
1215         nonHFECent[icent-1]->SetMarkerStyle(20);
1216         subtractedCent[icent-1]->Draw("samep");
1217         subtractedCent[icent-1]->SetMarkerStyle(24);
1218         if(icent == 1){
1219           lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum");
1220           lCentRaw->AddEntry(charmCent[0],"charm elecs");
1221           lCentRaw->AddEntry(convCent[0],"conversion elecs");
1222           lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs");
1223           lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum");
1224           lCentRaw->Draw("SAME");
1225         }
1226       }
1227     }
1228
1229     delete[] bins; 
1230
1231   }
1232   else{
1233     // Subtract 
1234     spectrumSubtracted->Add(backgroundGrid,-1.0);
1235   }
1236
1237   if(setBackground){
1238     if(fBackground) delete fBackground;
1239     fBackground = backgroundGrid;
1240   } else delete backgroundGrid;
1241
1242
1243   if(fDebugLevel > 0) {
1244
1245     Int_t ptprd;
1246     if(fBeamType==0) ptprd=0;
1247     if(fBeamType==1) ptprd=1;
1248     
1249     TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
1250     cbackgroundsubtraction->Divide(3,1);
1251     cbackgroundsubtraction->cd(1);
1252     gPad->SetLogy();
1253     TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
1254     TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
1255     CorrectFromTheWidth(measuredTH1Daftersubstraction);
1256     CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
1257     measuredTH1Daftersubstraction->SetStats(0);
1258     measuredTH1Daftersubstraction->SetTitle("");
1259     measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1260     measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1261     measuredTH1Daftersubstraction->SetMarkerStyle(25);
1262     measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
1263     measuredTH1Daftersubstraction->SetLineColor(kBlack);
1264     measuredTH1Dbeforesubstraction->SetStats(0);
1265     measuredTH1Dbeforesubstraction->SetTitle("");
1266     measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1267     measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1268     measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
1269     measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
1270     measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
1271     measuredTH1Daftersubstraction->Draw();
1272     measuredTH1Dbeforesubstraction->Draw("same");
1273     TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
1274     legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
1275     legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
1276     legsubstraction->Draw("same");
1277     cbackgroundsubtraction->cd(2);
1278     gPad->SetLogy();
1279     TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
1280     ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
1281     ratiomeasuredcontamination->SetTitle("");
1282     ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1283     ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1284     ratiomeasuredcontamination->Sumw2();
1285     ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
1286     ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
1287     ratiomeasuredcontamination->SetStats(0);
1288     ratiomeasuredcontamination->SetMarkerStyle(26);
1289     ratiomeasuredcontamination->SetMarkerColor(kBlack);
1290     ratiomeasuredcontamination->SetLineColor(kBlack);
1291     for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
1292       ratiomeasuredcontamination->SetBinError(k+1,0.0);
1293     }
1294     ratiomeasuredcontamination->Draw("P");
1295     cbackgroundsubtraction->cd(3);
1296     TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
1297     CorrectFromTheWidth(measuredTH1background);
1298     measuredTH1background->SetStats(0);
1299     measuredTH1background->SetTitle("");
1300     measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1301     measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1302     measuredTH1background->SetMarkerStyle(26);
1303     measuredTH1background->SetMarkerColor(kBlack);
1304     measuredTH1background->SetLineColor(kBlack);
1305     measuredTH1background->Draw();
1306     if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
1307
1308     if(fBeamType==1) {
1309
1310       TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1311       cbackgrounde->Divide(2,1);
1312       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1313       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1314      
1315       THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1316       TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1317       THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1318       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1319       Int_t nbbin = cenaxisb->GetNbins();
1320       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1321       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1322       for(Int_t binc = 0; binc < nbbin; binc++){
1323         TString titlee("BackgroundSubtraction_centrality_bin_");
1324         titlee += binc;
1325         TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1326         cbackground->Divide(2,1);
1327         cbackground->cd(1);
1328         gPad->SetLogy();
1329         cenaxisa->SetRange(binc+1,binc+1);
1330         cenaxisb->SetRange(binc+1,binc+1);
1331         TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1332         TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1333         CorrectFromTheWidth(aftersubstraction);
1334         CorrectFromTheWidth(beforesubstraction);
1335         aftersubstraction->SetStats(0);
1336         aftersubstraction->SetTitle((const char*)titlee);
1337         aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1338         aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1339         aftersubstraction->SetMarkerStyle(25);
1340         aftersubstraction->SetMarkerColor(kBlack);
1341         aftersubstraction->SetLineColor(kBlack);
1342         beforesubstraction->SetStats(0);
1343         beforesubstraction->SetTitle((const char*)titlee);
1344         beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1345         beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1346         beforesubstraction->SetMarkerStyle(24);
1347         beforesubstraction->SetMarkerColor(kBlue);
1348         beforesubstraction->SetLineColor(kBlue);
1349         aftersubstraction->Draw();
1350         beforesubstraction->Draw("same");
1351         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1352         lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1353         lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1354         lega->Draw("same");
1355         cbackgrounde->cd(1);
1356         gPad->SetLogy();
1357         TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1358         aftersubtractionn->SetMarkerStyle(stylee[binc]);
1359         aftersubtractionn->SetMarkerColor(colorr[binc]);
1360         if(binc==0) aftersubtractionn->Draw();
1361         else aftersubtractionn->Draw("same");
1362         legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1363         cbackgrounde->cd(2);
1364         gPad->SetLogy();
1365         TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1366         aftersubtractionng->SetMarkerStyle(stylee[binc]);
1367         aftersubtractionng->SetMarkerColor(colorr[binc]);
1368         if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1369         if(binc==0) aftersubtractionng->Draw();
1370         else aftersubtractionng->Draw("same");
1371         legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1372         cbackground->cd(2);
1373         TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1374         ratiocontamination->SetName("ratiocontamination");
1375         ratiocontamination->SetTitle((const char*)titlee);
1376         ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1377         ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1378         ratiocontamination->Add(aftersubstraction,-1.0);
1379         ratiocontamination->Divide(beforesubstraction);
1380         Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1381         for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1382           ratiocontamination->SetBinError(nbinpt+1,0.0);
1383         }
1384         ratiocontamination->SetStats(0);
1385         ratiocontamination->SetMarkerStyle(26);
1386         ratiocontamination->SetMarkerColor(kBlack);
1387         ratiocontamination->SetLineColor(kBlack);
1388         ratiocontamination->Draw("P");
1389       }
1390      
1391       cbackgrounde->cd(1);
1392       legtotal->Draw("same");
1393       cbackgrounde->cd(2);
1394       legtotalg->Draw("same");
1395
1396       cenaxisa->SetRange(0,nbbin);
1397       cenaxisb->SetRange(0,nbbin);
1398       if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
1399     }
1400   }
1401   
1402   return spectrumSubtracted;
1403 }
1404
1405 //____________________________________________________________
1406 AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1407   //
1408   // calculate charm background
1409   //
1410     Int_t ptpr = 0;
1411     Int_t nDim = 1;
1412     if(fBeamType==0)
1413     {
1414         ptpr=0;
1415     }
1416     if(fBeamType==1)
1417     {
1418         ptpr=1;
1419         nDim=2;
1420     }
1421
1422   Double_t evtnorm=0;
1423   if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
1424
1425   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1426   if(!mcContainer){
1427     AliError("MC Container not available");
1428     return NULL;
1429   }
1430
1431   if(!fCorrelation){
1432     AliError("No Correlation map available");
1433     return NULL;
1434   }
1435
1436   AliCFDataGrid *charmBackgroundGrid= 0x0;
1437   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1438
1439   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1440   Int_t* bins=new Int_t[2];
1441   bins[0]=charmbgaftertofpid->GetNbinsX();
1442
1443   if(fBeamType==1)
1444   {
1445       bins[0]=12;
1446       charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
1447       bins[1]=charmbgaftertofpid->GetNbinsX();
1448
1449   }
1450
1451   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
1452   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1453   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
1454
1455   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
1456
1457   Int_t *nBinpp=new Int_t[1];
1458   Int_t* binspp=new Int_t[1];
1459   binspp[0]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
1460
1461   Int_t *nBinPbPb=new Int_t[2];
1462   Int_t* binsPbPb=new Int_t[2];
1463   binsPbPb[1]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
1464   binsPbPb[0]=12;
1465
1466   Int_t looppt=binspp[0];
1467   if(fBeamType==1) looppt=binsPbPb[1];
1468
1469   for(Long_t iBin=1; iBin<= looppt;iBin++){
1470       if(fBeamType==0)
1471       {
1472           nBinpp[0]=iBin;
1473           charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
1474           charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
1475       }
1476       if(fBeamType==1)
1477       {
1478           // loop over centrality
1479           for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1480               nBinPbPb[0]=iiBin;
1481               nBinPbPb[1]=iBin;
1482               Double_t evtnormPbPb=0;
1483               if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1484               charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1485               charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1486           }
1487       }
1488   }
1489
1490   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
1491
1492   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
1493   weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1494   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
1495   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1496   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
1497
1498   // Efficiency (set efficiency to 1 for only folding) 
1499   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1500   efficiencyD->CalculateEfficiency(0,0);
1501
1502   // Folding
1503   if(fBeamType==0)nDim = 1;
1504   if(fBeamType==1)nDim = 2;
1505   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1506   folding.SetMaxNumberOfIterations(1);
1507   folding.Unfold();
1508
1509   // Results
1510   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1511   THnSparse* result=(THnSparse*)result1->Clone();
1512   charmBackgroundGrid->SetGrid(result);
1513   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
1514
1515   //Charm background evaluation plots
1516
1517   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1518   cCharmBgEval->Divide(3,1);
1519
1520   cCharmBgEval->cd(1);
1521
1522   if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
1523   if(fBeamType==1)
1524   {
1525       Double_t evtnormPbPb=0;
1526       for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
1527       {
1528           if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
1529       }
1530       charmbgaftertofpid->Scale(evtnormPbPb);
1531   }
1532
1533   CorrectFromTheWidth(charmbgaftertofpid);
1534   charmbgaftertofpid->SetMarkerStyle(25);
1535   charmbgaftertofpid->Draw("p");
1536   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1537   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1538   gPad->SetLogy();
1539
1540   CorrectFromTheWidth(charmbgafteripcut);
1541   charmbgafteripcut->SetMarkerStyle(24);
1542   charmbgafteripcut->Draw("samep");
1543
1544   CorrectFromTheWidth(charmbgafterweight);
1545   charmbgafterweight->SetMarkerStyle(24);
1546   charmbgafterweight->SetMarkerColor(4);
1547   charmbgafterweight->Draw("samep");
1548
1549   CorrectFromTheWidth(charmbgafterfolding);
1550   charmbgafterfolding->SetMarkerStyle(24);
1551   charmbgafterfolding->SetMarkerColor(2);
1552   charmbgafterfolding->Draw("samep");
1553
1554   cCharmBgEval->cd(2);
1555   parametrizedcharmpidipeff->SetMarkerStyle(24);
1556   parametrizedcharmpidipeff->Draw("p");
1557   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1558
1559   cCharmBgEval->cd(3);
1560   charmweightingfc->SetMarkerStyle(24);
1561   charmweightingfc->Draw("p");
1562   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1563   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1564
1565   cCharmBgEval->cd(1);
1566   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1567   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1568   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1569   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1570   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1571   legcharmbg->Draw("same");
1572
1573   cCharmBgEval->cd(2);
1574   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1575   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1576   legcharmbg2->Draw("same");
1577
1578   CorrectStatErr(charmBackgroundGrid);
1579   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
1580
1581   delete[] bins;
1582   delete[] nBinpp;
1583   delete[] binspp;
1584   delete[] nBinPbPb;
1585   delete[] binsPbPb;
1586
1587   if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
1588
1589   return charmBackgroundGrid;
1590 }
1591
1592 //____________________________________________________________
1593 AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1594   //
1595   // calculate conversion background
1596   //
1597   
1598   Double_t evtnorm[1] = {0.0};
1599   if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1600   printf("check event!!! %lf \n",evtnorm[0]);
1601   
1602   AliCFContainer *backgroundContainer = 0x0;
1603   
1604   if(fNonHFEsyst){
1605     backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
1606     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1607       backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
1608     }  
1609   }
1610   else{    
1611     // Background Estimate
1612     backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);    
1613   } 
1614   if(!backgroundContainer){
1615     AliError("MC background container not found");
1616     return NULL;
1617   }
1618   
1619   Int_t stepbackground = 3;
1620
1621   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1622   Int_t *nBinpp=new Int_t[1];
1623   Int_t* binspp=new Int_t[1];
1624   binspp[0]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1625
1626   Int_t *nBinPbPb=new Int_t[2];
1627   Int_t* binsPbPb=new Int_t[2];
1628   binsPbPb[1]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1629   binsPbPb[0]=12;
1630
1631   Int_t looppt=binspp[0];
1632   if(fBeamType==1) looppt=binsPbPb[1];
1633
1634   for(Long_t iBin=1; iBin<= looppt;iBin++){
1635       if(fBeamType==0)
1636       {
1637           nBinpp[0]=iBin;
1638           backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1639           backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1640       }
1641       if(fBeamType==1)
1642       {
1643           // loop over centrality
1644           for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1645               nBinPbPb[0]=iiBin;
1646               nBinPbPb[1]=iBin;
1647               Double_t evtnormPbPb=0;
1648               if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1649               backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1650               backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1651           }
1652       }
1653   }
1654   //end of workaround for statistical errors
1655
1656   AliCFDataGrid *weightedConversionContainer;
1657   if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
1658   else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
1659   weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1660   backgroundGrid->Multiply(weightedConversionContainer,1.0);
1661
1662   delete[] nBinpp;
1663   delete[] binspp;
1664   delete[] nBinPbPb;
1665   delete[] binsPbPb; 
1666
1667   return backgroundGrid;
1668 }
1669
1670
1671 //____________________________________________________________
1672 AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1673   //
1674   // calculate non-HFE background
1675   //
1676
1677   Double_t evtnorm[1] = {0.0};
1678   if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1679   printf("check event!!! %lf \n",evtnorm[0]);
1680   
1681   AliCFContainer *backgroundContainer = 0x0;
1682   if(fNonHFEsyst){
1683     backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
1684     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1685       backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
1686     } 
1687   }
1688   else{    
1689     // Background Estimate 
1690     backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1691   }
1692   if(!backgroundContainer){
1693     AliError("MC background container not found");
1694     return NULL;
1695   }
1696   
1697   
1698   Int_t stepbackground = 3;
1699
1700   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1701   Int_t *nBinpp=new Int_t[1];
1702   Int_t* binspp=new Int_t[1];
1703   binspp[0]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1704
1705   Int_t *nBinPbPb=new Int_t[2];
1706   Int_t* binsPbPb=new Int_t[2];
1707   binsPbPb[1]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1708   binsPbPb[0]=12;
1709
1710   Int_t looppt=binspp[0];
1711   if(fBeamType==1) looppt=binsPbPb[1];
1712
1713
1714   for(Long_t iBin=1; iBin<= looppt;iBin++){
1715       if(fBeamType==0)
1716       {
1717           nBinpp[0]=iBin;
1718           backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1719           backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1720       }
1721       if(fBeamType==1)
1722       {
1723           for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
1724               nBinPbPb[0]=iiBin;
1725               nBinPbPb[1]=iBin;
1726               Double_t evtnormPbPb=0;
1727               if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1728               backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1729               backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1730           }
1731       }
1732   }
1733   //end of workaround for statistical errors
1734   AliCFDataGrid *weightedNonHFEContainer;
1735   if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
1736   else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
1737   weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1738   backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
1739
1740   delete[] nBinpp;
1741   delete[] binspp;
1742   delete[] nBinPbPb;
1743   delete[] binsPbPb; 
1744
1745   return backgroundGrid;
1746 }
1747
1748 //____________________________________________________________
1749 AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1750   
1751   //
1752   // Apply TPC pid efficiency correction from parametrisation
1753   //
1754
1755   // Data in the right format
1756   AliCFDataGrid *dataGrid = 0x0;  
1757   if(bgsubpectrum) {
1758     dataGrid = bgsubpectrum;
1759   }
1760   else {
1761     
1762     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1763     if(!dataContainer){
1764       AliError("Data Container not available");
1765       return NULL;
1766     }
1767     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1768   } 
1769   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1770   result->SetName("ParametrizedEfficiencyBefore");
1771   THnSparse *h = result->GetGrid();
1772   Int_t nbdimensions = h->GetNdimensions();
1773   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1774   AliCFContainer *dataContainer = GetContainer(kDataContainer);
1775   if(!dataContainer){
1776     AliError("Data Container not available");
1777     return NULL;
1778   }
1779   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1780   dataContainerbis->Add(dataContainerbis,-1.0);
1781
1782
1783   Int_t* coord = new Int_t[nbdimensions];
1784   memset(coord, 0, sizeof(Int_t) * nbdimensions);
1785   Double_t* points = new Double_t[nbdimensions];
1786
1787   ULong64_t nEntries = h->GetNbins();
1788   for (ULong64_t i = 0; i < nEntries; ++i) {
1789     
1790     Double_t value = h->GetBinContent(i, coord);
1791     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1792     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1793     
1794     // Get the bin co-ordinates given an coord
1795     for (Int_t j = 0; j < nbdimensions; ++j)
1796       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1797
1798     if (!fEfficiencyFunction->IsInside(points))
1799          continue;
1800     TF1::RejectPoint(kFALSE);
1801
1802     // Evaulate function at points
1803     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1804     //printf("Value efficiency is %f\n",valueEfficiency);
1805
1806     if(valueEfficiency > 0.0) {
1807       h->SetBinContent(coord,value/valueEfficiency);
1808       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1809     }
1810     Double_t error = h->GetBinError(i);
1811     h->SetBinError(coord,error/valueEfficiency);
1812     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1813
1814    
1815   } 
1816
1817   delete[] coord;
1818   delete[] points;
1819
1820   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1821
1822   if(fDebugLevel > 0) {
1823
1824     TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1825     cEfficiencyParametrized->Divide(2,1);
1826     cEfficiencyParametrized->cd(1);
1827     TH1D *afterE = (TH1D *) resultt->Project(0);
1828     TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1829     CorrectFromTheWidth(afterE);
1830     CorrectFromTheWidth(beforeE);
1831     afterE->SetStats(0);
1832     afterE->SetTitle("");
1833     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1834     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1835     afterE->SetMarkerStyle(25);
1836     afterE->SetMarkerColor(kBlack);
1837     afterE->SetLineColor(kBlack);
1838     beforeE->SetStats(0);
1839     beforeE->SetTitle("");
1840     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1841     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1842     beforeE->SetMarkerStyle(24);
1843     beforeE->SetMarkerColor(kBlue);
1844     beforeE->SetLineColor(kBlue);
1845     gPad->SetLogy();
1846     afterE->Draw();
1847     beforeE->Draw("same");
1848     TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1849     legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1850     legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1851     legefficiencyparametrized->Draw("same");
1852     cEfficiencyParametrized->cd(2);
1853     fEfficiencyFunction->Draw();
1854     //cEfficiencyParametrized->cd(3);
1855     //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1856     //ratioefficiency->Divide(afterE);
1857     //ratioefficiency->Draw();
1858
1859     if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
1860
1861   }
1862
1863   return resultt;
1864
1865 }
1866 //____________________________________________________________
1867 AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1868   
1869   //
1870   // Apply TPC pid efficiency correction from V0
1871   //
1872
1873   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1874   if(!v0Container){
1875     AliError("V0 Container not available");
1876     return NULL;
1877   }
1878
1879   // Efficiency
1880   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1881   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1882
1883   // Data in the right format
1884   AliCFDataGrid *dataGrid = 0x0;  
1885   if(bgsubpectrum) {
1886     dataGrid = bgsubpectrum;
1887   }
1888   else {
1889     
1890     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1891     if(!dataContainer){
1892       AliError("Data Container not available");
1893       return NULL;
1894     }
1895
1896     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1897   } 
1898
1899   // Correct
1900   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1901   result->ApplyEffCorrection(*efficiencyD);
1902
1903   if(fDebugLevel > 0) {
1904
1905     Int_t ptpr;
1906     if(fBeamType==0) ptpr=0;
1907     if(fBeamType==1) ptpr=1;
1908     
1909     TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1910     cV0Efficiency->Divide(2,1);
1911     cV0Efficiency->cd(1);
1912     TH1D *afterE = (TH1D *) result->Project(ptpr);
1913     TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
1914     afterE->SetStats(0);
1915     afterE->SetTitle("");
1916     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1917     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1918     afterE->SetMarkerStyle(25);
1919     afterE->SetMarkerColor(kBlack);
1920     afterE->SetLineColor(kBlack);
1921     beforeE->SetStats(0);
1922     beforeE->SetTitle("");
1923     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1924     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1925     beforeE->SetMarkerStyle(24);
1926     beforeE->SetMarkerColor(kBlue);
1927     beforeE->SetLineColor(kBlue);
1928     afterE->Draw();
1929     beforeE->Draw("same");
1930     TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1931     legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1932     legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1933     legV0efficiency->Draw("same");
1934     cV0Efficiency->cd(2);
1935     TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
1936     efficiencyDproj->SetTitle("");
1937     efficiencyDproj->SetStats(0);
1938     efficiencyDproj->SetMarkerStyle(25);
1939     efficiencyDproj->Draw();
1940
1941     if(fBeamType==1) {
1942
1943       TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1944       cV0Efficiencye->Divide(2,1);
1945       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1946       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1947      
1948       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1949       TAxis *cenaxisa = sparseafter->GetAxis(0);
1950       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1951       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1952       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1953       TAxis *cenaxisc = efficiencya->GetAxis(0);
1954       Int_t nbbin = cenaxisb->GetNbins();
1955       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1956       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1957       for(Int_t binc = 0; binc < nbbin; binc++){
1958         TString titlee("V0Efficiency_centrality_bin_");
1959         titlee += binc;
1960         TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1961         ccV0Efficiency->Divide(2,1);
1962         ccV0Efficiency->cd(1);
1963         gPad->SetLogy();
1964         cenaxisa->SetRange(binc+1,binc+1);
1965         cenaxisb->SetRange(binc+1,binc+1);
1966         cenaxisc->SetRange(binc+1,binc+1);
1967         TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1968         TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1969         CorrectFromTheWidth(aftere);
1970         CorrectFromTheWidth(beforee);
1971         aftere->SetStats(0);
1972         aftere->SetTitle((const char*)titlee);
1973         aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1974         aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1975         aftere->SetMarkerStyle(25);
1976         aftere->SetMarkerColor(kBlack);
1977         aftere->SetLineColor(kBlack);
1978         beforee->SetStats(0);
1979         beforee->SetTitle((const char*)titlee);
1980         beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1981         beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1982         beforee->SetMarkerStyle(24);
1983         beforee->SetMarkerColor(kBlue);
1984         beforee->SetLineColor(kBlue);
1985         aftere->Draw();
1986         beforee->Draw("same");
1987         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1988         lega->AddEntry(beforee,"Before correction","p");
1989         lega->AddEntry(aftere,"After correction","p");
1990         lega->Draw("same");
1991         cV0Efficiencye->cd(1);
1992         gPad->SetLogy();
1993         TH1D *afteree = (TH1D *) aftere->Clone();
1994         afteree->SetMarkerStyle(stylee[binc]);
1995         afteree->SetMarkerColor(colorr[binc]);
1996         if(binc==0) afteree->Draw();
1997         else afteree->Draw("same");
1998         legtotal->AddEntry(afteree,(const char*) titlee,"p");
1999         cV0Efficiencye->cd(2);
2000         gPad->SetLogy();
2001         TH1D *aftereeu = (TH1D *) aftere->Clone();
2002         aftereeu->SetMarkerStyle(stylee[binc]);
2003         aftereeu->SetMarkerColor(colorr[binc]);
2004         if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
2005         if(binc==0) aftereeu->Draw();
2006         else aftereeu->Draw("same");
2007         legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
2008         ccV0Efficiency->cd(2);
2009         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
2010         efficiencyDDproj->SetTitle("");
2011         efficiencyDDproj->SetStats(0);
2012         efficiencyDDproj->SetMarkerStyle(25);
2013         efficiencyDDproj->Draw();
2014       }
2015      
2016       cV0Efficiencye->cd(1);
2017       legtotal->Draw("same");
2018       cV0Efficiencye->cd(2);
2019       legtotalg->Draw("same");
2020
2021       cenaxisa->SetRange(0,nbbin);
2022       cenaxisb->SetRange(0,nbbin);
2023       cenaxisc->SetRange(0,nbbin);
2024
2025       if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
2026     }
2027
2028   }
2029
2030   
2031   return result;
2032
2033 }
2034 //____________________________________________________________
2035 TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
2036   
2037   //
2038   // Unfold and eventually correct for efficiency the bgsubspectrum
2039   //
2040
2041   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2042   if(!mcContainer){
2043     AliError("MC Container not available");
2044     return NULL;
2045   }
2046
2047   if(!fCorrelation){
2048     AliError("No Correlation map available");
2049     return NULL;
2050   }
2051
2052   // Data 
2053   AliCFDataGrid *dataGrid = 0x0;  
2054   if(bgsubpectrum) {
2055     dataGrid = bgsubpectrum;
2056   }
2057   else {
2058
2059     AliCFContainer *dataContainer = GetContainer(kDataContainer);
2060     if(!dataContainer){
2061       AliError("Data Container not available");
2062       return NULL;
2063     }
2064
2065     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2066   } 
2067   
2068   // Guessed
2069   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2070   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2071
2072   // Efficiency
2073   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2074   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2075
2076   if(!fBeauty2ndMethod)
2077   {
2078       // Consider parameterized IP cut efficiency
2079       if(!fInclusiveSpectrum){
2080           Int_t* bins=new Int_t[1];
2081           bins[0]=35;
2082
2083           AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2084           beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
2085           efficiencyD->Multiply(beffContainer,1);
2086       }
2087   }
2088   
2089
2090   // Unfold 
2091   
2092   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2093   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2094   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2095   if(fSetSmoothing) unfolding.UseSmoothing();
2096   unfolding.Unfold();
2097
2098   // Results
2099   THnSparse* result = unfolding.GetUnfolded();
2100   THnSparse* residual = unfolding.GetEstMeasured();
2101   TList *listofresults = new TList;
2102   listofresults->SetOwner();
2103   listofresults->AddAt((THnSparse*)result->Clone(),0);
2104   listofresults->AddAt((THnSparse*)residual->Clone(),1);
2105
2106   if(fDebugLevel > 0) {
2107
2108     Int_t ptpr;
2109     if(fBeamType==0) ptpr=0;
2110     if(fBeamType==1) ptpr=1;
2111     
2112     TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
2113     cresidual->Divide(2,1);
2114     cresidual->cd(1);
2115     gPad->SetLogy();
2116     TGraphErrors* residualspectrumD = Normalize(residual);
2117     if(!residualspectrumD) {
2118       AliError("Number of Events not set for the normalization");
2119       return NULL;
2120     }
2121     residualspectrumD->SetTitle("");
2122     residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
2123     residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2124     residualspectrumD->SetMarkerStyle(26);
2125     residualspectrumD->SetMarkerColor(kBlue);
2126     residualspectrumD->SetLineColor(kBlue);
2127     residualspectrumD->Draw("AP");
2128     AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
2129     dataGridBis->SetName("dataGridBis"); 
2130     TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
2131     measuredspectrumD->SetTitle("");  
2132     measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
2133     measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2134     measuredspectrumD->SetMarkerStyle(25);
2135     measuredspectrumD->SetMarkerColor(kBlack);
2136     measuredspectrumD->SetLineColor(kBlack);
2137     measuredspectrumD->Draw("P");
2138     TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
2139     legres->AddEntry(residualspectrumD,"Residual","p");
2140     legres->AddEntry(measuredspectrumD,"Measured","p");
2141     legres->Draw("same");
2142     cresidual->cd(2);
2143     TH1D *residualTH1D = residual->Projection(ptpr);
2144     TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
2145     TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
2146     ratioresidual->SetName("ratioresidual");
2147     ratioresidual->SetTitle("");
2148     ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
2149     ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2150     ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
2151     ratioresidual->SetStats(0);
2152     ratioresidual->Draw();
2153
2154     if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
2155   }
2156
2157   return listofresults;
2158
2159 }
2160
2161 //____________________________________________________________
2162 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
2163   
2164   //
2165   // Apply unfolding and efficiency correction together to bgsubspectrum
2166   //
2167
2168   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
2169   if(!mcContainer){
2170     AliError("MC Container not available");
2171     return NULL;
2172   }
2173
2174   // Efficiency
2175   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2176   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2177
2178
2179   if(!fBeauty2ndMethod)
2180   {
2181       // Consider parameterized IP cut efficiency
2182       if(!fInclusiveSpectrum){
2183           Int_t* bins=new Int_t[1];
2184           bins[0]=35;
2185
2186           AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2187           beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
2188           efficiencyD->Multiply(beffContainer,1);
2189       }
2190   }
2191
2192   // Data in the right format
2193   AliCFDataGrid *dataGrid = 0x0;  
2194   if(bgsubpectrum) {
2195     dataGrid = bgsubpectrum;
2196   }
2197   else {
2198     
2199     AliCFContainer *dataContainer = GetContainer(kDataContainer);
2200     if(!dataContainer){
2201       AliError("Data Container not available");
2202       return NULL;
2203     }
2204
2205     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2206   } 
2207
2208   // Correct
2209   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
2210   result->ApplyEffCorrection(*efficiencyD);
2211
2212   if(fDebugLevel > 0) {
2213
2214     Int_t ptpr;
2215     if(fBeamType==0) ptpr=0;
2216     if(fBeamType==1) ptpr=1;
2217     
2218     printf("Step MC: %d\n",fStepMC);
2219     printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2220     printf("Step MC true: %d\n",fStepTrue);
2221     AliCFEffGrid  *efficiencymcPID = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2222     AliCFEffGrid  *efficiencymctrackinggeo = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
2223     AliCFEffGrid  *efficiencymcall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
2224     
2225     AliCFEffGrid  *efficiencyesdall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
2226     
2227     TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
2228     cefficiency->cd(1);
2229     TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
2230     efficiencymcPIDD->SetTitle("");
2231     efficiencymcPIDD->SetStats(0);
2232     efficiencymcPIDD->SetMarkerStyle(25);
2233     efficiencymcPIDD->Draw();
2234     TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
2235     efficiencymctrackinggeoD->SetTitle("");
2236     efficiencymctrackinggeoD->SetStats(0);
2237     efficiencymctrackinggeoD->SetMarkerStyle(26);
2238     efficiencymctrackinggeoD->Draw("same");
2239     TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
2240     efficiencymcallD->SetTitle("");
2241     efficiencymcallD->SetStats(0);
2242     efficiencymcallD->SetMarkerStyle(27);
2243     efficiencymcallD->Draw("same");
2244     TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
2245     efficiencyesdallD->SetTitle("");
2246     efficiencyesdallD->SetStats(0);
2247     efficiencyesdallD->SetMarkerStyle(24);
2248     efficiencyesdallD->Draw("same");
2249     TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
2250     legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
2251     legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
2252     legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
2253     legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
2254     legeff->Draw("same");
2255
2256     if(fBeamType==1) {
2257
2258       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
2259       TAxis *cenaxisa = sparseafter->GetAxis(0);
2260       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
2261       TAxis *cenaxisb = sparsebefore->GetAxis(0);
2262       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
2263       TAxis *cenaxisc = efficiencya->GetAxis(0);
2264       //Int_t nbbin = cenaxisb->GetNbins();
2265       //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
2266       //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
2267       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
2268         TString titlee("Efficiency_centrality_bin_");
2269         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
2270         titlee += "_";
2271         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
2272         TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
2273         cefficiencye->Divide(2,1);
2274         cefficiencye->cd(1);
2275         gPad->SetLogy();
2276         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2277         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2278         TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
2279         TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
2280         CorrectFromTheWidth(afterefficiency);
2281         CorrectFromTheWidth(beforeefficiency);
2282         afterefficiency->SetStats(0);
2283         afterefficiency->SetTitle((const char*)titlee);
2284         afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2285         afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2286         afterefficiency->SetMarkerStyle(25);
2287         afterefficiency->SetMarkerColor(kBlack);
2288         afterefficiency->SetLineColor(kBlack);
2289         beforeefficiency->SetStats(0);
2290         beforeefficiency->SetTitle((const char*)titlee);
2291         beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2292         beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2293         beforeefficiency->SetMarkerStyle(24);
2294         beforeefficiency->SetMarkerColor(kBlue);
2295         beforeefficiency->SetLineColor(kBlue);
2296         afterefficiency->Draw();
2297         beforeefficiency->Draw("same");
2298         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2299         lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
2300         lega->AddEntry(afterefficiency,"After efficiency correction","p");
2301         lega->Draw("same");
2302         cefficiencye->cd(2);
2303         cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
2304         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
2305         efficiencyDDproj->SetTitle("");
2306         efficiencyDDproj->SetStats(0);
2307         efficiencyDDproj->SetMarkerStyle(25);
2308         efficiencyDDproj->SetMarkerColor(2);
2309         efficiencyDDproj->Draw();
2310         cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
2311         TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
2312         efficiencyDDproja->SetTitle("");
2313         efficiencyDDproja->SetStats(0);
2314         efficiencyDDproja->SetMarkerStyle(26);
2315         efficiencyDDproja->SetMarkerColor(4);
2316         efficiencyDDproja->Draw("same");
2317       }
2318     }
2319
2320     if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
2321   }
2322   
2323   return result;
2324
2325 }
2326
2327 //__________________________________________________________________________________
2328 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
2329   //
2330   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2331   // Give the final pt spectrum to be compared
2332   //
2333
2334   if(fNEvents[i] > 0) {
2335
2336     Int_t ptpr = 0;
2337     if(fBeamType==0) ptpr=0;
2338     if(fBeamType==1) ptpr=1;
2339     
2340     TH1D* projection = spectrum->Projection(ptpr);
2341     CorrectFromTheWidth(projection);
2342     TGraphErrors *graphError = NormalizeTH1(projection,i);
2343     return graphError;
2344   
2345   }
2346     
2347   return 0x0;
2348   
2349
2350 }
2351 //__________________________________________________________________________________
2352 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
2353   //
2354   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2355   // Give the final pt spectrum to be compared
2356   //
2357   if(fNEvents[i] > 0) {
2358
2359     Int_t ptpr=0;
2360     if(fBeamType==0) ptpr=0;
2361     if(fBeamType==1) ptpr=1;
2362     
2363     TH1D* projection = (TH1D *) spectrum->Project(ptpr);
2364     CorrectFromTheWidth(projection);
2365     TGraphErrors *graphError = NormalizeTH1(projection,i);
2366
2367     return graphError;
2368     
2369   }
2370
2371   return 0x0;
2372   
2373
2374 }
2375 //__________________________________________________________________________________
2376 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
2377   //
2378   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2379   // Give the final pt spectrum to be compared
2380   //
2381   Double_t chargecoefficient = 0.5;
2382   if(fChargeChoosen != 0) chargecoefficient = 1.0;
2383
2384   Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2385   printf("Normalizing Eta Range %f\n", etarange);
2386   if(fNEvents[i] > 0) {
2387
2388     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2389     Double_t p = 0, dp = 0; Int_t point = 1;
2390     Double_t n = 0, dN = 0;
2391     Double_t nCorr = 0, dNcorr = 0;
2392     Double_t errdN = 0, errdp = 0;
2393     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2394       point = ibin - input->GetXaxis()->GetFirst();
2395       p = input->GetXaxis()->GetBinCenter(ibin);
2396       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2397       n = input->GetBinContent(ibin);
2398       dN = input->GetBinError(ibin);
2399       // New point
2400       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
2401       errdN = 1./(2. * TMath::Pi() * p);
2402       errdp = 1./(2. * TMath::Pi() * p*p) * n;
2403       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2404       
2405       spectrumNormalized->SetPoint(point, p, nCorr);
2406       spectrumNormalized->SetPointError(point, dp, dNcorr);
2407     }
2408     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2409     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2410     spectrumNormalized->SetMarkerStyle(22);
2411     spectrumNormalized->SetMarkerColor(kBlue);
2412     spectrumNormalized->SetLineColor(kBlue);
2413
2414     return spectrumNormalized;
2415     
2416   }
2417
2418   return 0x0;
2419   
2420
2421 }
2422 //__________________________________________________________________________________
2423 TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
2424   //
2425   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2426   // Give the final pt spectrum to be compared
2427   //
2428   Double_t chargecoefficient = 0.5;
2429   if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
2430
2431   Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2432   printf("Normalizing Eta Range %f\n", etarange);
2433   if(normalization > 0) {
2434
2435     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2436     Double_t p = 0, dp = 0; Int_t point = 1;
2437     Double_t n = 0, dN = 0;
2438     Double_t nCorr = 0, dNcorr = 0;
2439     Double_t errdN = 0, errdp = 0;
2440     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2441       point = ibin - input->GetXaxis()->GetFirst();
2442       p = input->GetXaxis()->GetBinCenter(ibin);
2443       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2444       n = input->GetBinContent(ibin);
2445       dN = input->GetBinError(ibin);
2446       // New point
2447       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
2448       errdN = 1./(2. * TMath::Pi() * p);
2449       errdp = 1./(2. * TMath::Pi() * p*p) * n;
2450       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2451       
2452       spectrumNormalized->SetPoint(point, p, nCorr);
2453       spectrumNormalized->SetPointError(point, dp, dNcorr);
2454     }
2455     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2456     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2457     spectrumNormalized->SetMarkerStyle(22);
2458     spectrumNormalized->SetMarkerColor(kBlue);
2459     spectrumNormalized->SetLineColor(kBlue);
2460     
2461     return spectrumNormalized;
2462     
2463   }
2464
2465   return 0x0;
2466   
2467
2468 }
2469 //____________________________________________________________
2470 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2471   //
2472   // Set the container for a given step to the 
2473   //
2474   if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
2475   fCFContainers->AddAt(cont, type);
2476 }
2477
2478 //____________________________________________________________
2479 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2480   //
2481   // Get Correction Framework Container for given type
2482   //
2483   if(!fCFContainers) return NULL;
2484   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2485 }
2486 //____________________________________________________________
2487 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
2488   //
2489   // Slice bin for a given source of electron
2490   // nDim is the number of dimension the corrections are done
2491   // dimensions are the definition of the dimensions
2492   // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2493   // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
2494   // centrality (-1 means we do not cut on centrality)
2495   //
2496   
2497   Double_t *varMin = new Double_t[container->GetNVar()],
2498            *varMax = new Double_t[container->GetNVar()];
2499
2500   Double_t *binLimits;
2501   for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2502     
2503     binLimits = new Double_t[container->GetNBins(ivar)+1];
2504     container->GetBinLimits(ivar,binLimits);
2505     varMin[ivar] = binLimits[0];
2506     varMax[ivar] = binLimits[container->GetNBins(ivar)];
2507     // source
2508     if(ivar == 4){
2509       if((source>= 0) && (source<container->GetNBins(ivar))) {
2510               varMin[ivar] = binLimits[source];
2511               varMax[ivar] = binLimits[source];
2512       }     
2513     }
2514     // charge
2515     if(ivar == 3) {
2516       if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
2517     }
2518     // eta
2519     if(ivar == 1){
2520       for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) 
2521         AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic))); 
2522       if(fEtaSelected){
2523         varMin[ivar] = fEtaRange[0];
2524         varMax[ivar] = fEtaRange[1];
2525       }
2526     }
2527     if(fEtaSelected){
2528       fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2529       fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2530       AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2531     }
2532     if(ivar == 5){
2533         if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
2534             varMin[ivar] = binLimits[centralitylow];
2535             varMax[ivar] = binLimits[centralityhigh];
2536
2537             TAxis *axistest = container->GetAxis(5,0);
2538             printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
2539             printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
2540             Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
2541             Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
2542             printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
2543         
2544         }
2545
2546     }
2547
2548
2549     delete[] binLimits;
2550   }
2551   
2552   AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2553   AddTemporaryObject(k);
2554   delete[] varMin; delete[] varMax;
2555
2556   return k;
2557
2558 }
2559
2560 //_________________________________________________________________________
2561 THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
2562   //
2563   // Slice correlation
2564   //
2565
2566   Int_t ndimensions = correlationmatrix->GetNdimensions();
2567   //printf("Number of dimension %d correlation map\n",ndimensions);
2568   if(ndimensions < (2*nDim)) {
2569     AliError("Problem in the dimensions");
2570     return NULL;
2571   }
2572   
2573   // Cut in centrality is centrality > -1
2574   if((centralitylow >=0) && (centralityhigh >=0)) {
2575
2576     TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
2577     TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
2578
2579     Int_t bins0 = axiscentrality0->GetNbins();
2580     Int_t bins1 = axiscentrality1->GetNbins();
2581     
2582     printf("Number of centrality bins: %d and %d\n",bins0,bins1);
2583     if(bins0 != bins1) {
2584       AliError("Problem in the dimensions");
2585       return NULL;
2586     }
2587     
2588     if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
2589       axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
2590       axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
2591
2592       Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
2593       Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
2594       Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
2595       Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
2596       printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
2597       printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
2598
2599       TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
2600       ctestcorrelation->cd(1);
2601       TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
2602       projection->Draw("colz");
2603
2604     }
2605     
2606   }
2607
2608
2609   Int_t ndimensionsContainer = (Int_t) ndimensions/2;
2610   //printf("Number of dimension %d container\n",ndimensionsContainer);
2611
2612   Int_t *dim = new Int_t[nDim*2];
2613   for(Int_t iter=0; iter < nDim; iter++){
2614     dim[iter] = dimensions[iter];
2615     dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
2616     //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
2617   }
2618     
2619   THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
2620
2621   delete[] dim; 
2622   return k;
2623   
2624 }
2625 //___________________________________________________________________________
2626 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2627   //
2628   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2629   //
2630
2631   TAxis *axis = h1->GetXaxis();
2632   Int_t nbinX = h1->GetNbinsX();
2633
2634   for(Int_t i = 1; i <= nbinX; i++) {
2635
2636     Double_t width = axis->GetBinWidth(i);
2637     Double_t content = h1->GetBinContent(i);
2638     Double_t error = h1->GetBinError(i); 
2639     h1->SetBinContent(i,content/width); 
2640     h1->SetBinError(i,error/width);
2641   }
2642
2643 }
2644
2645 //___________________________________________________________________________
2646 void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const { 
2647   //
2648   // Correct statistical error
2649   //
2650
2651   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2652   Int_t nbinX = h1->GetNbinsX();
2653   Int_t bins[1];
2654   for(Long_t i = 1; i <= nbinX; i++) {
2655     bins[0] = i;
2656     Float_t content = h1->GetBinContent(i);
2657     if(content>0){
2658       Float_t error = TMath::Sqrt(content);
2659       backgroundGrid->SetElementError(bins, error);
2660     }
2661   }
2662 }
2663
2664 //____________________________________________________________
2665 void AliHFEspectrum::AddTemporaryObject(TObject *o){
2666   // 
2667   // Emulate garbage collection on class level
2668   // 
2669   if(!fTemporaryObjects) {
2670     fTemporaryObjects= new TList;
2671     fTemporaryObjects->SetOwner();
2672   }
2673   fTemporaryObjects->Add(o);
2674 }
2675
2676 //____________________________________________________________
2677 void AliHFEspectrum::ClearObject(TObject *o){
2678   //
2679   // Do a safe deletion
2680   //
2681   if(fTemporaryObjects){
2682     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2683     delete o;
2684   } else{
2685     // just delete
2686     delete o;
2687   }
2688 }
2689 //_________________________________________________________________________
2690 TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
2691   AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
2692   return data;
2693 }
2694 //_________________________________________________________________________
2695 TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
2696   // 
2697   // Create efficiency grid and calculate efficiency
2698   // of step to step0
2699   //
2700   TString name("eff");
2701   name += step;
2702   name+= step0;
2703   AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2704   eff->CalculateEfficiency(step,step0);
2705   return eff;
2706 }
2707
2708 //____________________________________________________________________________
2709 THnSparse* AliHFEspectrum::GetCharmWeights(){
2710   
2711   //
2712   // Measured D->e based weighting factors
2713   //
2714
2715   const Int_t nDimpp=1;
2716   Int_t nBinpp[nDimpp] = {35};
2717   Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
2718   const Int_t nDimPbPb=2;
2719   Int_t nBinPbPb[nDimPbPb] = {11,35};
2720   Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2721   Int_t loopcentr=1;
2722   Int_t looppt=nBinpp[0];
2723   if(fBeamType==0)
2724   {
2725       fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
2726       fWeightCharm->SetBinEdges(0, ptbinning1);
2727   }
2728   if(fBeamType==1)
2729   {
2730       fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
2731       fWeightCharm->SetBinEdges(1, ptbinning1);
2732       fWeightCharm->SetBinEdges(0, kCentralityRange);
2733       loopcentr=nBinPbPb[0];
2734   }
2735
2736   // Weighting factor for pp
2737   Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
2738   
2739   // Weighting factor for PbPb (0-20%)
2740   //Double_t weight[35]={0.641897,  0.640472,  0.615228,  0.650469,  0.737762,  0.847867,  1.009317,  1.158594,  1.307482,  1.476973,  1.551131,  1.677131,  1.785478,  1.888933,  2.017957,  2.074757,  1.926700,  1.869495,  1.546558,  1.222873,  1.160313,  0.903375,  0.799642,  0.706244,  0.705449,  0.599947,  0.719570,  0.499422,  0.703978,  0.477452,  0.325057,  0.093391,  0.096675,  0.000000,  0.000000};
2741
2742   // Weighting factor for PbPb (40-80%)
2743   //Double_t weight[35]={0.181953,  0.173248,  0.166799,  0.182558,  0.206581,  0.236955,  0.279390,  0.329129,  0.365260,  0.423059,  0.452057,  0.482726,  0.462627,  0.537770,  0.584663,  0.579452,  0.587194,  0.499498,  0.443299,  0.398596,  0.376695,  0.322331,  0.260890,  0.374834,  0.249114,  0.310330,  0.287326,  0.243174,  0.758945,  0.138867,  0.170576,  0.107797,  0.011390,  0.000000,  0.000000};
2744
2745   //points
2746   Double_t pt[1];
2747   Double_t contents[2];
2748
2749   for(int icentr=0; icentr<loopcentr; icentr++)
2750   {
2751       for(int i=0; i<looppt; i++){
2752           pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2753           if(fBeamType==1)
2754           {
2755               contents[0]=icentr;
2756               contents[1]=pt[0];
2757           }
2758           if(fBeamType==0)
2759           {
2760               contents[0]=pt[0];
2761           }
2762
2763           fWeightCharm->Fill(contents,weight[i]);
2764       }
2765   }
2766
2767   Int_t nDimSparse = fWeightCharm->GetNdimensions();
2768   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2769   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2770
2771   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2772       binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
2773       nBins *= binsvar[iVar];
2774   }
2775
2776   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2777   // loop that sets 0 error in each bin
2778   for (Long_t iBin=0; iBin<nBins; iBin++) {
2779     Long_t bintmp = iBin ;
2780     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2781       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2782       bintmp /= binsvar[iVar] ;
2783     }
2784     fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
2785   }
2786
2787   delete[] binsvar;
2788   delete[] binfill;
2789
2790   return fWeightCharm;
2791 }
2792
2793 //____________________________________________________________________________
2794 void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
2795
2796    // TOF PID efficiencies
2797    Int_t ptpr;
2798    if(fBeamType==0) ptpr=0;
2799    if(fBeamType==1) ptpr=1;
2800
2801    Int_t loopcentr=1;
2802    const Int_t nCentralitybinning=11; //number of centrality bins
2803    if(fBeamType==1)
2804    {
2805      loopcentr=nCentralitybinning;
2806    }
2807
2808    TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2809    TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2810    TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
2811
2812    TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
2813    cefficiencyParamtof->cd();
2814
2815    AliCFContainer *mccontainermcD = 0x0;
2816    AliCFContainer *mccontaineresdD = 0x0;
2817    TH1D* efficiencysigTOFPIDD[nCentralitybinning];
2818    TH1D* efficiencyTOFPIDD[nCentralitybinning];
2819    TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
2820    TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
2821    Int_t source = -1; //get parameterized TOF PID efficiencies
2822
2823    for(int icentr=0; icentr<loopcentr; icentr++) {
2824       // signal sample
2825       if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2826       else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2827       AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
2828       efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2829
2830       // mb sample for double check
2831       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2832       else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2833       AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
2834       efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2835
2836       // mb sample with reconstructed variables
2837       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2838       else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2839       AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
2840       efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2841
2842       // mb sample with reconstructed variables
2843       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2844       else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2845       AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
2846       efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2847
2848       //fill histo
2849       efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
2850       efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
2851       efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
2852       efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
2853       efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
2854       efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
2855       efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
2856       efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
2857
2858       //fit (mc enhenced sample)
2859       fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2860       efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
2861       efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2862       fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
2863
2864       //fit (esd enhenced sample)
2865       efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
2866       efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2867       fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
2868
2869    }
2870
2871    // draw (for PbPb, only 1st bin)
2872    //sig mc
2873    efficiencysigTOFPIDD[0]->SetTitle("");
2874    efficiencysigTOFPIDD[0]->SetStats(0);
2875    efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
2876    efficiencysigTOFPIDD[0]->SetMarkerColor(2);
2877    efficiencysigTOFPIDD[0]->SetLineColor(2);
2878    efficiencysigTOFPIDD[0]->Draw();
2879
2880    //mb mc
2881    efficiencyTOFPIDD[0]->SetTitle("");
2882    efficiencyTOFPIDD[0]->SetStats(0);
2883    efficiencyTOFPIDD[0]->SetMarkerStyle(24);
2884    efficiencyTOFPIDD[0]->SetMarkerColor(4);
2885    efficiencyTOFPIDD[0]->SetLineColor(4);
2886    efficiencyTOFPIDD[0]->Draw("same");
2887
2888    //sig esd
2889    efficiencysigesdTOFPIDD[0]->SetTitle("");
2890    efficiencysigesdTOFPIDD[0]->SetStats(0);
2891    efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
2892    efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
2893    efficiencysigesdTOFPIDD[0]->SetLineColor(3);
2894    efficiencysigesdTOFPIDD[0]->Draw("same");
2895
2896    //mb esd
2897    efficiencyesdTOFPIDD[0]->SetTitle("");
2898    efficiencyesdTOFPIDD[0]->SetStats(0);
2899    efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
2900    efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
2901    efficiencyesdTOFPIDD[0]->SetLineColor(1);
2902    efficiencyesdTOFPIDD[0]->Draw("same");
2903
2904    //signal mc fit
2905    if(fEfficiencyTOFPIDD[0]){
2906      fEfficiencyTOFPIDD[0]->SetLineColor(2);
2907      fEfficiencyTOFPIDD[0]->Draw("same");
2908    }
2909    //mb esd fit
2910    if(fEfficiencyesdTOFPIDD[0]){
2911        fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
2912        fEfficiencyesdTOFPIDD[0]->Draw("same");
2913      }
2914
2915    TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
2916    legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
2917    legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
2918    legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
2919    legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
2920    legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
2921    legtofeff->Draw("same");
2922
2923
2924    TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
2925    cefficiencyParamIP->cd();
2926    gStyle->SetOptStat(0);
2927
2928    // IP cut efficiencies
2929    for(int icentr=0; icentr<loopcentr; icentr++)  {
2930
2931      AliCFContainer *charmCombined = 0x0; 
2932      AliCFContainer *beautyCombined = 0x0;
2933      AliCFContainer *beautyCombinedesd = 0x0;
2934
2935      printf("centrality printing %i \n",icentr);
2936
2937      source = 0; //charm enhenced
2938      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2939      else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2940      AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
2941      efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
2942
2943      charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");  
2944
2945      source = 1; //beauty enhenced
2946      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2947      else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2948      AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
2949      efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
2950
2951      beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined"); 
2952
2953      if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2954      else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2955      AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
2956      efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2957
2958      beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
2959
2960      source = 0; //charm mb
2961      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2962      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2963      AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
2964      efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
2965
2966      charmCombined->Add(mccontainermcD); 
2967      AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined); 
2968      efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
2969
2970      source = 1; //beauty mb
2971      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2972      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2973      AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
2974      efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
2975
2976      beautyCombined->Add(mccontainermcD);
2977      AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined); 
2978      efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
2979
2980      if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2981      else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2982      AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
2983      efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2984
2985      beautyCombinedesd->Add(mccontaineresdD);
2986      AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
2987      efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2988
2989      source = 2; //conversion mb
2990      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2991      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2992      AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
2993      efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
2994
2995      source = 3; //non HFE except for the conversion mb
2996      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2997      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2998      AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
2999      efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3000
3001      if(fIPEffCombinedSamples){
3002        fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb 
3003        fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
3004        fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
3005      }
3006      else{
3007        fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
3008        fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
3009        fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
3010      }
3011      fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
3012      fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
3013      fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
3014      fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
3015
3016    }
3017
3018    if(fBeamType==0){
3019      AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
3020      fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
3021
3022      AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
3023      fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
3024    }
3025
3026    for(int icentr=0; icentr<loopcentr; icentr++)  {
3027      fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3028      fipfit->SetLineColor(2);
3029      fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
3030      fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3031      if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
3032      else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
3033
3034      fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3035      fipfit->SetLineColor(6);
3036      fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
3037      fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
3038      if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
3039      else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
3040
3041      fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3042      fipfit->SetLineColor(1);
3043      fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
3044      fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3045      fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
3046      
3047      if(fIPParameterizedEff){
3048        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3049        fipfitnonhfe->SetLineColor(3);
3050        fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
3051        fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3052        fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
3053
3054        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3055        fipfitnonhfe->SetLineColor(4);
3056        fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
3057        fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3058        fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
3059      }
3060    }
3061
3062    // draw (for PbPb, only 1st bin)
3063    fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
3064    fEfficiencyCharmSigD[0]->SetMarkerColor(1);
3065    fEfficiencyCharmSigD[0]->SetLineColor(1);
3066    fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
3067    fEfficiencyBeautySigD[0]->SetMarkerColor(2);
3068    fEfficiencyBeautySigD[0]->SetLineColor(2);
3069    fEfficiencyBeautySigesdD[0]->SetStats(0);
3070    fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
3071    fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
3072    fEfficiencyBeautySigesdD[0]->SetLineColor(6);
3073    fCharmEff[0]->SetMarkerStyle(24);
3074    fCharmEff[0]->SetMarkerColor(1);
3075    fCharmEff[0]->SetLineColor(1);
3076    fBeautyEff[0]->SetMarkerStyle(24);
3077    fBeautyEff[0]->SetMarkerColor(2);
3078    fBeautyEff[0]->SetLineColor(2);
3079    fConversionEff[0]->SetMarkerStyle(24);
3080    fConversionEff[0]->SetMarkerColor(3);
3081    fConversionEff[0]->SetLineColor(3);
3082    fNonHFEEff[0]->SetMarkerStyle(24);
3083    fNonHFEEff[0]->SetMarkerColor(4);
3084    fNonHFEEff[0]->SetLineColor(4);
3085
3086    fEfficiencyCharmSigD[0]->Draw();
3087    fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
3088    fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
3089
3090    fEfficiencyBeautySigD[0]->Draw("same");
3091    fEfficiencyBeautySigesdD[0]->Draw("same");
3092    //fCharmEff[0]->Draw("same");
3093    //fBeautyEff[0]->Draw("same");
3094
3095    if(fBeamType==0){
3096      fConversionEffbgc->SetMarkerStyle(25);
3097      fConversionEffbgc->SetMarkerColor(3);
3098      fConversionEffbgc->SetLineColor(3);
3099      fNonHFEEffbgc->SetMarkerStyle(25);
3100      fNonHFEEffbgc->SetMarkerColor(4);
3101      fNonHFEEffbgc->SetLineColor(4);
3102      fConversionEffbgc->Draw("same");
3103      fNonHFEEffbgc->Draw("same");
3104    }
3105    else{
3106      fConversionEff[0]->Draw("same");
3107      fNonHFEEff[0]->Draw("same");
3108    }
3109    if(fEfficiencyIPBeautyD[0])
3110       fEfficiencyIPBeautyD[0]->Draw("same");
3111    if(fEfficiencyIPBeautyesdD[0])
3112      fEfficiencyIPBeautyesdD[0]->Draw("same");
3113    if( fEfficiencyIPCharmD[0])
3114      fEfficiencyIPCharmD[0]->Draw("same");
3115    if(fIPParameterizedEff){
3116      fEfficiencyIPConversionD[0]->Draw("same");
3117      fEfficiencyIPNonhfeD[0]->Draw("same");
3118    }
3119    TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
3120    legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
3121    legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
3122    legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
3123    legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
3124    legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
3125    legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
3126    //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
3127    //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
3128    legipeff->Draw("same");
3129    gPad->SetGrid();
3130    //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
3131 }
3132
3133 //____________________________________________________________________________
3134 THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
3135   //
3136   // Return beauty electron IP cut efficiency
3137   //
3138
3139   const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3140   const Int_t nCentralitybinning=11;//number of centrality bins
3141   Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3142   Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3143   Int_t ptpr = 0;
3144   Int_t nDim=1;  //dimensions of the efficiency weighting grid
3145   if(fBeamType==1)
3146   {
3147     ptpr=1;
3148     nDim=2; //dimensions of the efficiency weighting grid
3149   }
3150   Int_t nBin[1] = {nPtbinning1};
3151   Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3152
3153
3154   THnSparseF *ipcut;
3155   if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
3156   else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3157  
3158   for(Int_t idim = 0; idim < nDim; idim++)
3159   {
3160     if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
3161     if(nDim==2)
3162       {
3163         ipcut->SetBinEdges(0, kCentralityRange);
3164         ipcut->SetBinEdges(1, kPtRange);
3165       }
3166   }
3167   Double_t pt[1];
3168   Double_t weight[nCentralitybinning];
3169   Double_t weightErr[nCentralitybinning];
3170   Double_t contents[2];
3171
3172   for(Int_t a=0;a<11;a++)
3173   {
3174       weight[a] = 1.0;
3175       weightErr[a] = 1.0;
3176   }
3177
3178
3179   Int_t looppt=nBin[0];
3180   Int_t loopcentr=1;
3181   Int_t ibin[2];
3182   if(fBeamType==1)
3183   {
3184       loopcentr=nBinPbPb[0];
3185   }
3186
3187
3188   for(int icentr=0; icentr<loopcentr; icentr++)
3189   {
3190       for(int i=0; i<looppt; i++)
3191       {
3192           pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3193           if(isMCpt){
3194             if(fEfficiencyIPBeautyD[icentr]){
3195               weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
3196               weightErr[icentr] = 0;
3197             }
3198             else{
3199               printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3200               weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1); 
3201               weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3202             }
3203           }
3204           else{
3205             if(fEfficiencyIPBeautyesdD[icentr]){
3206               weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
3207               weightErr[icentr] = 0;
3208             }
3209             else{
3210               printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3211               weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
3212               weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3213             }
3214           }
3215
3216           if(fBeamType==1){
3217               contents[0]=icentr;
3218               contents[1]=pt[0];
3219               ibin[0]=icentr;
3220               ibin[1]=i+1;
3221           }
3222           if(fBeamType==0){
3223               contents[0]=pt[0];
3224               ibin[0]=i+1;
3225           }
3226           ipcut->Fill(contents,weight[icentr]);
3227           ipcut->SetBinError(ibin,weightErr[icentr]);
3228       }
3229   } 
3230
3231   Int_t nDimSparse = ipcut->GetNdimensions();
3232   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3233   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3234
3235   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3236       binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
3237       nBins *= binsvar[iVar];
3238   }
3239
3240   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3241   // loop that sets 0 error in each bin
3242   for (Long_t iBin=0; iBin<nBins; iBin++) {
3243     Long_t bintmp = iBin ;
3244     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3245       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3246       bintmp /= binsvar[iVar] ;
3247     }
3248     //ipcut->SetBinError(binfill,0.); // put 0 everywhere
3249   }
3250
3251   delete[] binsvar;
3252   delete[] binfill;
3253
3254   return ipcut;
3255 }
3256
3257 //____________________________________________________________________________
3258 THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
3259   //
3260   // Return PID x IP cut efficiency
3261   //
3262     const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3263     const Int_t nCentralitybinning=11;//number of centrality bins
3264     Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3265     Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3266     Int_t ptpr = 0;
3267     Int_t nDim=1;  //dimensions of the efficiency weighting grid
3268     if(fBeamType==1)
3269     {
3270         ptpr=1;
3271         nDim=2; //dimensions of the efficiency weighting grid
3272     }
3273     Int_t nBin[1] = {nPtbinning1};
3274     Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3275
3276     THnSparseF *pideff;
3277     if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
3278     else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3279     for(Int_t idim = 0; idim < nDim; idim++)
3280     {
3281
3282         if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
3283         if(nDim==2)
3284         {
3285             pideff->SetBinEdges(0, kCentralityRange);
3286             pideff->SetBinEdges(1, kPtRange);
3287         }
3288     }
3289
3290   Double_t pt[1];
3291   Double_t weight[nCentralitybinning];
3292   Double_t weightErr[nCentralitybinning];
3293   Double_t contents[2];
3294
3295   for(Int_t a=0;a<nCentralitybinning;a++)
3296   {
3297       weight[a] = 1.0;
3298       weightErr[a] = 1.0;
3299   }
3300
3301   Int_t looppt=nBin[0];
3302   Int_t loopcentr=1;
3303   Int_t ibin[2];
3304   if(fBeamType==1)
3305   {
3306       loopcentr=nBinPbPb[0];
3307   }
3308
3309   for(int icentr=0; icentr<loopcentr; icentr++)
3310   {
3311       Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
3312       for(int i=0; i<looppt; i++)
3313       {
3314           pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3315
3316           Double_t tofpideff = 0.;
3317           Double_t tofpideffesd = 0.;
3318           if(fEfficiencyTOFPIDD[icentr])
3319             tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]); 
3320           else{
3321             printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3322           }  
3323           if(fEfficiencyesdTOFPIDD[icentr])
3324             tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
3325           else{
3326             printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3327           }
3328
3329           //tof pid eff x tpc pid eff x ip cut eff
3330           if(fIPParameterizedEff){
3331             if(source==0) {
3332               if(fEfficiencyIPCharmD[icentr]){
3333                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3334                 weightErr[icentr] = 0; 
3335               }
3336               else{
3337                 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3338                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3339                 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1); 
3340               }
3341             } 
3342             else if(source==2) {
3343               if(fEfficiencyIPConversionD[icentr]){
3344                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]); 
3345                 weightErr[icentr] = 0; 
3346               }
3347               else{
3348                 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
3349                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
3350                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3351               }
3352             }
3353             else if(source==3) {
3354               if(fEfficiencyIPNonhfeD[icentr]){
3355                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]); 
3356                 weightErr[icentr] = 0; 
3357               }
3358               else{
3359                 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
3360                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
3361                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3362               }  
3363             }
3364           }
3365           else{
3366             if(source==0){ 
3367               if(fEfficiencyIPCharmD[icentr]){
3368                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3369                 weightErr[icentr] = 0;
3370               }
3371               else{
3372                 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3373                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3374                 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3375               }
3376             }
3377             else if(source==2){
3378               if(fBeamType==0){
3379                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
3380                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
3381               }
3382               else{
3383                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
3384                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3385               }
3386             }
3387             else if(source==3){
3388               if(fBeamType==0){
3389                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
3390                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
3391               }
3392               else{ 
3393                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
3394                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3395               }
3396             }
3397           }
3398
3399           if(fBeamType==1){
3400               contents[0]=icentr;
3401               contents[1]=pt[0];
3402               ibin[0]=icentr;
3403               ibin[1]=i+1;
3404           }
3405           if(fBeamType==0){
3406               contents[0]=pt[0];
3407               ibin[0]=i+1;
3408           }
3409
3410           pideff->Fill(contents,weight[icentr]);
3411           pideff->SetBinError(ibin,weightErr[icentr]);
3412       }
3413   }
3414
3415   Int_t nDimSparse = pideff->GetNdimensions();
3416   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3417   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3418
3419   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3420       binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
3421       nBins *= binsvar[iVar];