Fix
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEspectrum.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15 //
16 // Class for spectrum correction
17 // Subtraction of hadronic background, Unfolding of the data and
18 // Renormalization done here
19 // The following containers have to be set:
20 //  - Correction framework container for real data
21 //  - Correction framework container for MC (Efficiency Map)
22 //  - Correction framework container for background coming from data
23 //  - Correction framework container for background coming from MC
24 //
25 //  Author: 
26 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
27 //            Markus Fasel <M.Fasel@gsi.de>
28 //
29
30 #include <TArrayD.h>
31 #include <TH1.h>
32 #include <TList.h>
33 #include <TObjArray.h>
34 #include <TROOT.h>
35 #include <TCanvas.h>
36 #include <TLegend.h>
37 #include <TStyle.h>
38 #include <TMath.h>
39 #include <TAxis.h>
40 #include <TGraphErrors.h>
41 #include <TFile.h>
42 #include <TPad.h>
43 #include <TH2D.h>
44 #include <TF1.h>
45
46 #include "AliPID.h"
47 #include "AliCFContainer.h"
48 #include "AliCFDataGrid.h"
49 #include "AliCFEffGrid.h"
50 #include "AliCFGridSparse.h"
51 #include "AliCFUnfolding.h"
52 #include "AliLog.h"
53
54 #include "AliHFEspectrum.h"
55 #include "AliHFEcuts.h"
56 #include "AliHFEcontainer.h"
57 #include "AliHFEtools.h"
58
59 ClassImp(AliHFEspectrum)
60
61 //____________________________________________________________
62 AliHFEspectrum::AliHFEspectrum(const char *name):
63   TNamed(name, ""),
64   fCFContainers(NULL),
65   fTemporaryObjects(NULL),
66   fCorrelation(NULL),
67   fBackground(NULL),
68   fEfficiencyFunction(NULL),
69   fWeightCharm(NULL),
70   fInclusiveSpectrum(kTRUE),
71   fDumpToFile(kFALSE),
72   fEtaSelected(kFALSE),
73   fUnSetCorrelatedErrors(kTRUE),
74   fSetSmoothing(kFALSE),
75   fIPanaHadronBgSubtract(kFALSE),
76   fIPanaCharmBgSubtract(kFALSE),
77   fIPanaConversionBgSubtract(kFALSE),
78   fIPanaNonHFEBgSubtract(kFALSE),
79   fNonHFEbgMethod2(kFALSE),
80   fNonHFEmode(0),
81   fNbDimensions(1),
82   fNMCEvents(0),
83   fNMCbgEvents(0),
84   fStepMC(-1),
85   fStepTrue(-1),
86   fStepData(-1),
87   fStepBeforeCutsV0(-1),
88   fStepAfterCutsV0(-1),
89   fStepGuessedUnfolding(-1),
90   fNumberOfIterations(5),
91   fChargeChoosen(-1),
92   fNCentralityBinAtTheEnd(0),
93   fHadronEffbyIPcut(NULL),
94   fConversionEff(NULL),
95   fNonHFEEff(NULL),
96   fBeamType(0),
97   fDebugLevel(0)
98 {
99   //
100   // Default constructor
101   //
102
103   for(Int_t k = 0; k < 20; k++){
104     fNEvents[k] = 0;
105     fLowBoundaryCentralityBinAtTheEnd[k] = 0;
106     fHighBoundaryCentralityBinAtTheEnd[k] = 0;
107   }
108   memset(fEtaRange, 0, sizeof(Double_t) * 2);
109 }
110
111 //____________________________________________________________
112 AliHFEspectrum::~AliHFEspectrum(){
113   //
114   // Destructor
115   //
116   if(fCFContainers) delete fCFContainers;
117   if(fTemporaryObjects){
118     fTemporaryObjects->Clear();
119     delete fTemporaryObjects;
120   }
121 }
122 //____________________________________________________________
123 Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
124   //
125   // Init what we need for the correction:
126   //
127   // Raw spectrum, hadron contamination
128   // MC efficiency maps, correlation matrix
129   // V0 efficiency if wanted
130   //
131   // This for a given dimension.
132   // If no inclusive spectrum, then take only efficiency map for beauty electron
133   // and the appropriate correlation matrix
134   //
135     Int_t kNdim = 3;
136     if(fBeamType==0) kNdim=3;
137     if(fBeamType==1) kNdim=4;
138
139     Int_t dims[kNdim];
140     // Get the requested format
141     if(fBeamType==0)
142     {
143         // pp analysis
144         switch(fNbDimensions){
145         case 1:   dims[0] = 0;
146         break;
147         case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
148         break;
149         case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
150         break;
151         default:
152             AliError("Container with this number of dimensions not foreseen (yet)");
153             return kFALSE;
154         };
155     }
156
157     if(fBeamType==1)
158     {
159         // PbPb analysis; centrality as first dimension
160       Int_t nbDimensions = fNbDimensions;
161       fNbDimensions = fNbDimensions + 1;
162         switch(nbDimensions){
163         case 1:
164           dims[0] = 5;
165           dims[1] = 0;
166           break;
167         case 2:
168           dims[0] = 5;
169           for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
170           break;
171         case 3:
172           dims[0] = 5;
173           for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
174           break;
175         default:
176             AliError("Container with this number of dimensions not foreseen (yet)");
177             return kFALSE;
178         };
179     }
180
181
182
183   // Data container: raw spectrum + hadron contamination  
184   AliCFContainer *datacontainer = 0x0; 
185   if(fInclusiveSpectrum) {
186     datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
187   }
188   else{
189     datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
190   }
191   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
192   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
193
194   AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen);
195   AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen);
196   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
197
198   SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
199   SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
200
201   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps 
202   AliCFContainer *mccontaineresd = 0x0;
203   AliCFContainer *mccontainermc = 0x0;
204   AliCFContainer *mccontainermcbg = 0x0;
205   AliCFContainer *nonHFEweightedContainer = 0x0;
206   AliCFContainer *convweightedContainer = 0x0;
207    
208   if(fInclusiveSpectrum) {
209     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
210     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
211   }
212   else {
213     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
214     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
215     mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
216     nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
217     convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
218     if((!nonHFEweightedContainer) || (!convweightedContainer)) return kFALSE;
219   }
220   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;  
221
222   Int_t source = -1;
223   if(!fInclusiveSpectrum) source = 1; //beauty
224   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen);
225   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
226   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
227   SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
228   SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
229
230   // set charm, nonHFE container to estimate BG
231   if(!fInclusiveSpectrum){
232    source = 0; //charm
233    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
234    SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
235
236    source = 2; //conversion
237    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
238    AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
239    efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
240    fConversionEff = (TH1D*)efficiencyConv->Project(0);
241
242    source = 3; //non HFE except for the conversion
243    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
244    AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
245    efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
246    fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0);
247
248    AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
249    SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
250    AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
251    SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
252   }
253  
254 // MC container: correlation matrix
255   THnSparseF *mccorrelation = 0x0;
256   if(fInclusiveSpectrum) {
257     if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
258     else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
259     else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
260     else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
261     else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
262     
263     if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
264   }
265   else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
266   //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
267   if(!mccorrelation) return kFALSE;
268   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
269   if(!mccorrelationD) {
270     printf("No correlation\n");
271     return kFALSE;
272   }
273   SetCorrelation(mccorrelationD);
274
275   // V0 container Electron, pt eta phi
276   if(v0hfecontainer) {
277     AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
278     if(!containerV0) return kFALSE;
279     AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron+1);
280     if(!containerV0Electron) return kFALSE;
281     SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
282   }
283  
284
285   if(fDebugLevel>0){
286    
287     AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
288     contaminationspectrum->SetName("contaminationspectrum");
289     TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
290     ccontaminationspectrum->Divide(3,1);
291     ccontaminationspectrum->cd(1);
292     TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
293     TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
294     TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
295     contaminationspectrum2dpteta->SetStats(0);
296     contaminationspectrum2dpteta->SetTitle("");
297     contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
298     contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
299     contaminationspectrum2dptphi->SetStats(0);
300     contaminationspectrum2dptphi->SetTitle("");
301     contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
302     contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
303     contaminationspectrum2detaphi->SetStats(0);
304     contaminationspectrum2detaphi->SetTitle("");
305     contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
306     contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
307     contaminationspectrum2dptphi->Draw("colz");
308     ccontaminationspectrum->cd(2);
309     contaminationspectrum2dpteta->Draw("colz");
310     ccontaminationspectrum->cd(3);
311     contaminationspectrum2detaphi->Draw("colz");
312
313     TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
314     ccorrelation->cd(1);
315     if(mccorrelationD) {
316       if(fBeamType==0){
317         TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
318         ptcorrelation->Draw("colz");
319       }
320       if(fBeamType==1){
321         TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
322         ptcorrelation->Draw("colz");
323       }
324     }
325   }
326
327
328   return kTRUE;
329  
330   
331 }
332
333
334 //____________________________________________________________
335 Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
336   //
337   // Correct the spectrum for efficiency and unfolding
338   // with both method and compare
339   //
340   
341   gStyle->SetPalette(1);
342   gStyle->SetOptStat(1111);
343   gStyle->SetPadBorderMode(0);
344   gStyle->SetCanvasColor(10);
345   gStyle->SetPadLeftMargin(0.13);
346   gStyle->SetPadRightMargin(0.13);
347
348   printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
349
350   ///////////////////////////
351   // Check initialization
352   ///////////////////////////
353
354   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
355     AliInfo("You have to init before");
356     return kFALSE;
357   }
358   
359   if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
360     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
361     return kFALSE;
362   }
363  
364   SetNumberOfIteration(10);
365   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
366     
367   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
368   //////////////////////////////////
369   // Subtract hadron background
370   /////////////////////////////////
371   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
372   if(subtractcontamination) {
373     dataspectrumaftersubstraction = SubtractBackground(kTRUE);
374     dataGridAfterFirstSteps = dataspectrumaftersubstraction;
375   }
376
377   ////////////////////////////////////////////////
378   // Correct for TPC efficiency from V0
379   ///////////////////////////////////////////////
380   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
381   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
382   if(dataContainerV0){printf("Got the V0 container\n");
383     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
384     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;  
385   }
386
387   //////////////////////////////////////////////////////////////////////////////
388   // Correct for efficiency parametrized (if TPC efficiency is parametrized)
389   /////////////////////////////////////////////////////////////////////////////
390   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
391   if(fEfficiencyFunction){
392     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
393     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;  
394   }
395     
396   ///////////////
397   // Unfold
398   //////////////
399   TList *listunfolded = Unfold(dataGridAfterFirstSteps);
400   if(!listunfolded){
401     printf("Unfolded failed\n");
402     return kFALSE;
403   }
404   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
405   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
406   if(!correctedspectrum){
407     AliError("No corrected spectrum\n");
408     return kFALSE;
409   }
410   if(!residualspectrum){
411     AliError("No residul spectrum\n");
412     return kFALSE;
413   }   
414   
415   /////////////////////
416   // Simply correct
417   ////////////////////
418   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
419   
420
421   //////////
422   // Plot
423   //////////
424   if(fDebugLevel > 0.0) {
425
426     Int_t ptpr = 0;
427     if(fBeamType==0) ptpr=0;
428     if(fBeamType==1) ptpr=1;
429   
430     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
431     ccorrected->Divide(2,1);
432     ccorrected->cd(1);
433     gPad->SetLogy();
434     TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
435     correctedspectrumD->SetTitle("");
436     correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
437     correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
438     correctedspectrumD->SetMarkerStyle(26);
439     correctedspectrumD->SetMarkerColor(kBlue);
440     correctedspectrumD->SetLineColor(kBlue);
441     correctedspectrumD->Draw("AP");
442     TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
443     alltogetherspectrumD->SetTitle("");
444     alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
445     alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
446     alltogetherspectrumD->SetMarkerStyle(25);
447     alltogetherspectrumD->SetMarkerColor(kBlack);
448     alltogetherspectrumD->SetLineColor(kBlack);
449     alltogetherspectrumD->Draw("P");
450     TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
451     legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
452     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
453     legcorrected->Draw("same");
454     ccorrected->cd(2);
455     TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
456     TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
457     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
458     ratiocorrected->SetName("ratiocorrected");
459     ratiocorrected->SetTitle("");
460     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
461     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
462     ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
463     ratiocorrected->SetStats(0);
464     ratiocorrected->Draw();
465
466     //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
467     //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
468     //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
469     //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
470
471     TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
472     TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
473     TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
474     TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
475
476     
477
478     if(fBeamType==1) {
479
480       TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
481       ccorrectedallspectra->Divide(2,1);
482       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
483       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
484       
485       THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
486       TAxis *cenaxisa = sparsesu->GetAxis(0);
487       THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
488       TAxis *cenaxisb = sparsed->GetAxis(0);
489       Int_t nbbin = cenaxisb->GetNbins();
490       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
491       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
492       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
493         TString titlee("corrected_centrality_bin_");
494         titlee += "[";
495         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
496         titlee += "_";
497         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
498         titlee += "[";
499         TString titleec("corrected_check_projection_bin_");
500         titleec += "[";
501         titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
502         titleec += "_";
503         titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
504         titleec += "[";
505         TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
506         titleenameunotnormalized += "[";
507         titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
508         titleenameunotnormalized += "_";
509         titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
510         titleenameunotnormalized += "[";
511         TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
512         titleenameunormalized += "[";
513         titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
514         titleenameunormalized += "_";
515         titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];      
516         titleenameunormalized += "[";
517         TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
518         titleenamednotnormalized += "[";
519         titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
520         titleenamednotnormalized += "_";
521         titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
522         titleenamednotnormalized += "[";
523         TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
524         titleenamednormalized += "[";
525         titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
526         titleenamednormalized += "_";
527         titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];      
528         titleenamednormalized += "[";
529         Int_t nbEvents = 0;
530         for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
531           printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
532           nbEvents += fNEvents[k];
533         }
534         Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
535         Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
536         printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
537         Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
538         Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
539         printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
540         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
541         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
542         TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
543         ccorrectedcheck->cd(1);
544         TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
545         TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
546         aftersuc->Draw();
547         aftersdc->Draw("same");
548         TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
549         ccorrectede->Divide(2,1);
550         ccorrectede->cd(1);
551         gPad->SetLogy();
552         TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
553         CorrectFromTheWidth(aftersu);
554         aftersu->SetName((const char*)titleenameunotnormalized);
555         unfoldingspectrac[binc] = *aftersu;
556         ccorrectede->cd(1);
557         TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
558         aftersun->SetTitle("");
559         aftersun->GetYaxis()->SetTitleOffset(1.5);
560         aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
561         aftersun->SetMarkerStyle(26);
562         aftersun->SetMarkerColor(kBlue);
563         aftersun->SetLineColor(kBlue);
564         aftersun->Draw("AP");
565         aftersun->SetName((const char*)titleenameunormalized);
566         unfoldingspectracn[binc] = *aftersun;
567         ccorrectede->cd(1);
568         TH1D *aftersd = (TH1D *) sparsed->Projection(1);
569         CorrectFromTheWidth(aftersd);
570         aftersd->SetName((const char*)titleenamednotnormalized);
571         correctedspectrac[binc] = *aftersd;
572         ccorrectede->cd(1);
573         TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
574         aftersdn->SetTitle("");
575         aftersdn->GetYaxis()->SetTitleOffset(1.5);
576         aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
577         aftersdn->SetMarkerStyle(25);
578         aftersdn->SetMarkerColor(kBlack);
579         aftersdn->SetLineColor(kBlack);
580         aftersdn->Draw("P");
581         aftersdn->SetName((const char*)titleenamednormalized);
582         correctedspectracn[binc] = *aftersdn;
583         TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
584         legcorrectedud->AddEntry(aftersun,"Corrected","p");
585         legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
586         legcorrectedud->Draw("same");
587         ccorrectedallspectra->cd(1);
588         gPad->SetLogy();
589         TH1D *aftersunn = (TH1D *) aftersun->Clone();
590         aftersunn->SetMarkerStyle(stylee[binc]);
591         aftersunn->SetMarkerColor(colorr[binc]);
592         if(binc==0) aftersunn->Draw("AP");
593         else aftersunn->Draw("P");
594         legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
595         ccorrectedallspectra->cd(2);
596         gPad->SetLogy();
597         TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
598         aftersdnn->SetMarkerStyle(stylee[binc]);
599         aftersdnn->SetMarkerColor(colorr[binc]);
600         if(binc==0) aftersdnn->Draw("AP");
601         else aftersdnn->Draw("P");
602         legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
603         ccorrectede->cd(2);
604         TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
605         TString titleee("ratiocorrected_bin_");
606         titleee += binc;
607         ratiocorrectedbinc->SetName((const char*) titleee);
608         ratiocorrectedbinc->SetTitle("");
609         ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
610         ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
611         ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
612         ratiocorrectedbinc->SetStats(0);
613         ratiocorrectedbinc->Draw();
614       }
615
616       ccorrectedallspectra->cd(1);
617       legtotal->Draw("same");
618       ccorrectedallspectra->cd(2);
619       legtotalg->Draw("same");
620       
621       cenaxisa->SetRange(0,nbbin);
622       cenaxisb->SetRange(0,nbbin);
623
624     }
625
626     // Dump to file if needed
627     if(fDumpToFile) {
628       TFile *out = new TFile("finalSpectrum.root","recreate");
629       correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
630       correctedspectrumD->Write();
631       alltogetherspectrumD->SetName("AlltogetherSpectrum");
632       alltogetherspectrumD->Write();
633       ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
634       ratiocorrected->Write();
635       correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
636       correctedspectrum->Write();
637       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
638       alltogetherCorrection->Write();
639       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
640         unfoldingspectrac[binc].Write();
641         unfoldingspectracn[binc].Write();
642         correctedspectrac[binc].Write();
643         correctedspectracn[binc].Write();
644       }
645       out->Close(); delete out;
646     }
647
648     if (unfoldingspectrac) delete[] unfoldingspectrac;
649     if (unfoldingspectracn)  delete[] unfoldingspectracn;
650     if (correctedspectrac) delete[] correctedspectrac;
651     if (correctedspectracn) delete[] correctedspectracn;
652
653   }
654
655
656   
657
658   return kTRUE;
659 }
660
661 //____________________________________________________________
662 Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
663   //
664   // Correct the spectrum for efficiency and unfolding for beauty analysis
665   // with both method and compare
666   //
667   
668   gStyle->SetPalette(1);
669   gStyle->SetOptStat(1111);
670   gStyle->SetPadBorderMode(0);
671   gStyle->SetCanvasColor(10);
672   gStyle->SetPadLeftMargin(0.13);
673   gStyle->SetPadRightMargin(0.13);
674
675   ///////////////////////////
676   // Check initialization
677   ///////////////////////////
678
679   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
680     AliInfo("You have to init before");
681     return kFALSE;
682   }
683   
684   if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
685     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
686     return kFALSE;
687   }
688  
689   SetNumberOfIteration(10);
690   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
691     
692   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
693   //////////////////////////////////
694   // Subtract hadron background
695   /////////////////////////////////
696   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
697   if(subtractcontamination) {
698     dataspectrumaftersubstraction = SubtractBackground(kTRUE);
699     dataGridAfterFirstSteps = dataspectrumaftersubstraction;
700   }
701
702   /////////////////////////////////////////////////////////////////////////////////////////
703   // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
704   /////////////////////////////////////////////////////////////////////////////////////////
705
706   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
707   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
708   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
709
710   if(fEfficiencyFunction){
711     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
712     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;  
713   }
714   else if(dataContainerV0){
715     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
716     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;  
717   }  
718  
719
720   ///////////////
721   // Unfold
722   //////////////
723   TList *listunfolded = Unfold(dataGridAfterFirstSteps);
724   if(!listunfolded){
725     printf("Unfolded failed\n");
726     return kFALSE;
727   }
728   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
729   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
730   if(!correctedspectrum){
731     AliError("No corrected spectrum\n");
732     return kFALSE;
733   }
734   if(!residualspectrum){
735     AliError("No residual spectrum\n");
736     return kFALSE;
737   }   
738   
739   /////////////////////
740   // Simply correct
741   ////////////////////
742   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
743   
744
745   //////////
746   // Plot
747   //////////
748   if(fDebugLevel > 0.0) {
749   
750     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
751     ccorrected->Divide(2,1);
752     ccorrected->cd(1);
753     gPad->SetLogy();
754     TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
755     correctedspectrumD->SetTitle("");
756     correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
757     correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
758     correctedspectrumD->SetMarkerStyle(26);
759     correctedspectrumD->SetMarkerColor(kBlue);
760     correctedspectrumD->SetLineColor(kBlue);
761     correctedspectrumD->Draw("AP");
762     TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
763     alltogetherspectrumD->SetTitle("");
764     alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
765     alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
766     alltogetherspectrumD->SetMarkerStyle(25);
767     alltogetherspectrumD->SetMarkerColor(kBlack);
768     alltogetherspectrumD->SetLineColor(kBlack);
769     alltogetherspectrumD->Draw("P");
770     TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
771     legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
772     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
773     legcorrected->Draw("same");
774     ccorrected->cd(2);
775     TH1D *correctedTH1D = correctedspectrum->Projection(0);
776     TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
777     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
778     ratiocorrected->SetName("ratiocorrected");
779     ratiocorrected->SetTitle("");
780     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
781     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
782     ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
783     ratiocorrected->SetStats(0);
784     ratiocorrected->Draw();
785
786
787     // Dump to file if needed
788
789     if(fDumpToFile) {
790       
791       TFile *out;
792       if(fNonHFEmode == 1){
793         out = new TFile("finalSpectrumLow.root","recreate");
794       }
795       else if(fNonHFEmode == 2){
796         out = new TFile("finalSpectrumUp.root","recreate");
797       }
798       else{
799         out = new TFile("finalSpectrum.root","recreate");
800       }
801       out->cd();
802       //
803       correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
804       correctedspectrumD->Write();
805       alltogetherspectrumD->SetName("AlltogetherSpectrum");
806       alltogetherspectrumD->Write();
807       ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
808       ratiocorrected->Write();
809       //
810       correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
811       correctedspectrum->Write();
812       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
813       alltogetherCorrection->Write();
814       //
815       out->Close(); 
816       delete out;
817     }
818   
819
820   }
821
822   return kTRUE;
823 }
824
825 //____________________________________________________________
826 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
827   //
828   // Apply background subtraction
829   //
830   
831   // Raw spectrum
832   AliCFContainer *dataContainer = GetContainer(kDataContainer);
833   if(!dataContainer){
834     AliError("Data Container not available");
835     return NULL;
836   }
837   printf("Step data: %d\n",fStepData);
838   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
839
840   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
841   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); 
842  
843   // Background Estimate
844   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
845   if(!backgroundContainer){
846     AliError("MC background container not found");
847     return NULL;
848   }
849  
850   Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
851   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
852
853   if(!fInclusiveSpectrum){
854     //Background subtraction for IP analysis
855     TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(0);
856     CorrectFromTheWidth(measuredTH1Draw);
857     TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
858     rawspectra->cd();
859     measuredTH1Draw->SetMarkerStyle(20);
860     measuredTH1Draw->Draw();
861     if(fIPanaHadronBgSubtract){
862       // Hadron background
863       printf("Hadron background for IP analysis subtracted!\n");
864       Int_t* bins=new Int_t[1];
865       TH1D* htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
866       bins[0]=htemp->GetNbinsX();
867       AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
868       hbgContainer->SetGrid(fHadronEffbyIPcut);
869       backgroundGrid->Multiply(hbgContainer,1);
870       // draw raw hadron bg spectra
871       TH1D *hadronbg= (TH1D *) backgroundGrid->Project(0);
872       CorrectFromTheWidth(hadronbg);
873       hadronbg->SetMarkerColor(7);
874       hadronbg->SetMarkerStyle(20);
875       rawspectra->cd();
876       hadronbg->Draw("samep");
877       // subtract hadron contamination
878       spectrumSubtracted->Add(backgroundGrid,-1.0);
879     }
880     if(fIPanaCharmBgSubtract){
881       // Charm background
882       printf("Charm background for IP analysis subtracted!\n");
883       AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
884       // draw charm bg spectra
885       TH1D *charmbg= (TH1D *) charmbgContainer->Project(0);
886       CorrectFromTheWidth(charmbg);
887       charmbg->SetMarkerColor(3);
888       charmbg->SetMarkerStyle(20);
889       rawspectra->cd();
890       charmbg->Draw("samep");
891       // subtract charm background
892       spectrumSubtracted->Add(charmbgContainer,-1.0);
893     }
894     if(fIPanaConversionBgSubtract){
895       // Conversion background
896       AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
897       // draw conversion bg spectra
898       TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(0);
899       CorrectFromTheWidth(conversionbg);
900       conversionbg->SetMarkerColor(4);
901       conversionbg->SetMarkerStyle(20);
902       rawspectra->cd();
903       conversionbg->Draw("samep");
904       // subtract conversion background
905       spectrumSubtracted->Add(conversionbgContainer,-1.0);
906       printf("Conversion background subtraction is preliminary!\n");
907     }
908     if(fIPanaNonHFEBgSubtract){
909       // NonHFE background
910       AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
911       // draw Dalitz/dielectron bg spectra
912       TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(0);
913       CorrectFromTheWidth(nonhfebg);
914       nonhfebg->SetMarkerColor(6);
915       nonhfebg->SetMarkerStyle(20);
916       rawspectra->cd();
917       nonhfebg->Draw("samep");
918       // subtract Dalitz/dielectron background
919       spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
920       printf("Non HFE background subtraction is preliminary!\n");
921     }
922     TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(0);
923     CorrectFromTheWidth(rawbgsubtracted);
924     rawbgsubtracted->SetMarkerStyle(24);
925     rawspectra->cd();
926     rawbgsubtracted->Draw("samep");
927   }
928   else{
929     // Subtract 
930     spectrumSubtracted->Add(backgroundGrid,-1.0);
931   }
932
933   if(setBackground){
934     if(fBackground) delete fBackground;
935     fBackground = backgroundGrid;
936   } else delete backgroundGrid;
937
938
939   if(fDebugLevel > 0) {
940
941     Int_t ptpr;
942     if(fBeamType==0) ptpr=0;
943     if(fBeamType==1) ptpr=1;
944     
945     TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
946     cbackgroundsubtraction->Divide(3,1);
947     cbackgroundsubtraction->cd(1);
948     gPad->SetLogy();
949     TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
950     TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
951     CorrectFromTheWidth(measuredTH1Daftersubstraction);
952     CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
953     measuredTH1Daftersubstraction->SetStats(0);
954     measuredTH1Daftersubstraction->SetTitle("");
955     measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
956     measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
957     measuredTH1Daftersubstraction->SetMarkerStyle(25);
958     measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
959     measuredTH1Daftersubstraction->SetLineColor(kBlack);
960     measuredTH1Dbeforesubstraction->SetStats(0);
961     measuredTH1Dbeforesubstraction->SetTitle("");
962     measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
963     measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
964     measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
965     measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
966     measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
967     measuredTH1Daftersubstraction->Draw();
968     measuredTH1Dbeforesubstraction->Draw("same");
969     TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
970     legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
971     legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
972     legsubstraction->Draw("same");
973     cbackgroundsubtraction->cd(2);
974     gPad->SetLogy();
975     TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
976     ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
977     ratiomeasuredcontamination->SetTitle("");
978     ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
979     ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
980     ratiomeasuredcontamination->Sumw2();
981     ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
982     ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
983     ratiomeasuredcontamination->SetStats(0);
984     ratiomeasuredcontamination->SetMarkerStyle(26);
985     ratiomeasuredcontamination->SetMarkerColor(kBlack);
986     ratiomeasuredcontamination->SetLineColor(kBlack);
987     for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
988       ratiomeasuredcontamination->SetBinError(k+1,0.0);
989     }
990     ratiomeasuredcontamination->Draw("P");
991     cbackgroundsubtraction->cd(3);
992     TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
993     CorrectFromTheWidth(measuredTH1background);
994     measuredTH1background->SetStats(0);
995     measuredTH1background->SetTitle("");
996     measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
997     measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
998     measuredTH1background->SetMarkerStyle(26);
999     measuredTH1background->SetMarkerColor(kBlack);
1000     measuredTH1background->SetLineColor(kBlack);
1001     measuredTH1background->Draw();
1002
1003     if(fBeamType==1) {
1004
1005       TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1006       cbackgrounde->Divide(2,1);
1007       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1008       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1009      
1010       THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1011       TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1012       THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1013       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1014       Int_t nbbin = cenaxisb->GetNbins();
1015       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1016       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1017       for(Int_t binc = 0; binc < nbbin; binc++){
1018         TString titlee("BackgroundSubtraction_centrality_bin_");
1019         titlee += binc;
1020         TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1021         cbackground->Divide(2,1);
1022         cbackground->cd(1);
1023         gPad->SetLogy();
1024         cenaxisa->SetRange(binc+1,binc+1);
1025         cenaxisb->SetRange(binc+1,binc+1);
1026         TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1027         TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1028         CorrectFromTheWidth(aftersubstraction);
1029         CorrectFromTheWidth(beforesubstraction);
1030         aftersubstraction->SetStats(0);
1031         aftersubstraction->SetTitle((const char*)titlee);
1032         aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1033         aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1034         aftersubstraction->SetMarkerStyle(25);
1035         aftersubstraction->SetMarkerColor(kBlack);
1036         aftersubstraction->SetLineColor(kBlack);
1037         beforesubstraction->SetStats(0);
1038         beforesubstraction->SetTitle((const char*)titlee);
1039         beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1040         beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1041         beforesubstraction->SetMarkerStyle(24);
1042         beforesubstraction->SetMarkerColor(kBlue);
1043         beforesubstraction->SetLineColor(kBlue);
1044         aftersubstraction->Draw();
1045         beforesubstraction->Draw("same");
1046         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1047         lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1048         lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1049         lega->Draw("same");
1050         cbackgrounde->cd(1);
1051         gPad->SetLogy();
1052         TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1053         aftersubtractionn->SetMarkerStyle(stylee[binc]);
1054         aftersubtractionn->SetMarkerColor(colorr[binc]);
1055         if(binc==0) aftersubtractionn->Draw();
1056         else aftersubtractionn->Draw("same");
1057         legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1058         cbackgrounde->cd(2);
1059         gPad->SetLogy();
1060         TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1061         aftersubtractionng->SetMarkerStyle(stylee[binc]);
1062         aftersubtractionng->SetMarkerColor(colorr[binc]);
1063         if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1064         if(binc==0) aftersubtractionng->Draw();
1065         else aftersubtractionng->Draw("same");
1066         legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1067         cbackground->cd(2);
1068         TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1069         ratiocontamination->SetName("ratiocontamination");
1070         ratiocontamination->SetTitle((const char*)titlee);
1071         ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1072         ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1073         ratiocontamination->Add(aftersubstraction,-1.0);
1074         ratiocontamination->Divide(beforesubstraction);
1075         Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1076         for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1077           ratiocontamination->SetBinError(nbinpt+1,0.0);
1078         }
1079         ratiocontamination->SetStats(0);
1080         ratiocontamination->SetMarkerStyle(26);
1081         ratiocontamination->SetMarkerColor(kBlack);
1082         ratiocontamination->SetLineColor(kBlack);
1083         ratiocontamination->Draw("P");
1084         
1085       }
1086      
1087       cbackgrounde->cd(1);
1088       legtotal->Draw("same");
1089       cbackgrounde->cd(2);
1090       legtotalg->Draw("same");
1091
1092       cenaxisa->SetRange(0,nbbin);
1093       cenaxisb->SetRange(0,nbbin);
1094       
1095     }
1096    
1097        
1098   }
1099   
1100   return spectrumSubtracted;
1101 }
1102
1103 //____________________________________________________________
1104 AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1105   //
1106   // calculate charm background
1107   //
1108
1109   Double_t evtnorm=0;
1110   if(fNMCbgEvents) evtnorm= double(fNEvents[0])/double(fNMCbgEvents);
1111
1112   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1113   if(!mcContainer){
1114     AliError("MC Container not available");
1115     return NULL;
1116   }
1117
1118   if(!fCorrelation){
1119     AliError("No Correlation map available");
1120     return NULL;
1121   }
1122
1123   AliCFDataGrid *charmBackgroundGrid= 0x0;
1124   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1125   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1126
1127   Int_t* bins=new Int_t[1];
1128   bins[0]=charmbgaftertofpid->GetNbinsX();
1129   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",1,bins);
1130   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1131   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
1132
1133   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
1134   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
1135
1136   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
1137   weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1138   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
1139   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1140   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
1141
1142   // Efficiency (set efficiency to 1 for only folding) 
1143   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1144   efficiencyD->CalculateEfficiency(0,0);
1145
1146   // Folding 
1147   Int_t nDim = 1;
1148   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1149   folding.SetMaxNumberOfIterations(1);
1150   folding.Unfold();
1151
1152   // Results
1153   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1154   THnSparse* result=(THnSparse*)result1->Clone();
1155   charmBackgroundGrid->SetGrid(result);
1156   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
1157
1158   //Charm background evaluation plots
1159
1160   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1161   cCharmBgEval->Divide(3,1);
1162
1163   cCharmBgEval->cd(1);
1164   charmbgaftertofpid->Scale(evtnorm);
1165   CorrectFromTheWidth(charmbgaftertofpid);
1166   charmbgaftertofpid->SetMarkerStyle(25);
1167   charmbgaftertofpid->Draw("p");
1168   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1169   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1170   gPad->SetLogy();
1171
1172   CorrectFromTheWidth(charmbgafteripcut);
1173   charmbgafteripcut->SetMarkerStyle(24);
1174   charmbgafteripcut->Draw("samep");
1175
1176   CorrectFromTheWidth(charmbgafterweight);
1177   charmbgafterweight->SetMarkerStyle(24);
1178   charmbgafterweight->SetMarkerColor(4);
1179   charmbgafterweight->Draw("samep");
1180
1181   CorrectFromTheWidth(charmbgafterfolding);
1182   charmbgafterfolding->SetMarkerStyle(24);
1183   charmbgafterfolding->SetMarkerColor(2);
1184   charmbgafterfolding->Draw("samep");
1185
1186   cCharmBgEval->cd(2);
1187   parametrizedcharmpidipeff->SetMarkerStyle(24);
1188   parametrizedcharmpidipeff->Draw("p");
1189   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1190
1191   cCharmBgEval->cd(3);
1192   charmweightingfc->SetMarkerStyle(24);
1193   charmweightingfc->Draw("p");
1194   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1195   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1196
1197   cCharmBgEval->cd(1);
1198   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1199   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1200   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1201   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1202   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1203   legcharmbg->Draw("same");
1204
1205   cCharmBgEval->cd(2);
1206   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1207   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1208   legcharmbg2->Draw("same");
1209
1210   CorrectStatErr(charmBackgroundGrid);
1211
1212   return charmBackgroundGrid;
1213 }
1214
1215 //____________________________________________________________
1216 AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1217   //
1218   // calculate conversion background
1219   //
1220
1221   Double_t evtnorm[1] = {0.0};
1222   if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1223   printf("check event!!! %lf \n",evtnorm[0]);
1224
1225   // Background Estimate
1226   AliCFContainer *backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
1227   if(!backgroundContainer){
1228     AliError("MC background container not found");
1229     return NULL;
1230   }
1231
1232   Int_t stepbackground = 3;
1233   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1234   backgroundGrid->Scale(evtnorm);
1235
1236   Int_t* bins=new Int_t[1];
1237   bins[0]=fConversionEff->GetNbinsX();
1238   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1239   weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1240   backgroundGrid->Multiply(weightedConversionContainer,1.0);
1241   
1242   CorrectStatErr(backgroundGrid);
1243   
1244   return backgroundGrid;
1245 }
1246
1247
1248 //____________________________________________________________
1249 AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1250   //
1251   // calculate non-HFE background
1252   //
1253
1254   Double_t evtnorm[1] = {0.0};
1255   if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1256   printf("check event!!! %lf \n",evtnorm[0]);
1257
1258   // Background Estimate
1259   AliCFContainer *backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1260   if(!backgroundContainer){
1261     AliError("MC background container not found");
1262     return NULL;
1263   }
1264
1265
1266   Int_t stepbackground = 3;
1267   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1268   backgroundGrid->Scale(evtnorm);
1269
1270   Int_t* bins=new Int_t[1];
1271   bins[0]=fNonHFEEff->GetNbinsX();
1272   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1273   weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1274   backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
1275
1276   CorrectStatErr(backgroundGrid);
1277
1278   return backgroundGrid;
1279 }
1280
1281 //____________________________________________________________
1282 AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1283   
1284   //
1285   // Apply TPC pid efficiency correction from parametrisation
1286   //
1287
1288   // Data in the right format
1289   AliCFDataGrid *dataGrid = 0x0;  
1290   if(bgsubpectrum) {
1291     dataGrid = bgsubpectrum;
1292   }
1293   else {
1294     
1295     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1296     if(!dataContainer){
1297       AliError("Data Container not available");
1298       return NULL;
1299     }
1300
1301     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1302   } 
1303
1304   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1305   result->SetName("ParametrizedEfficiencyBefore");
1306   THnSparse *h = result->GetGrid();
1307   Int_t nbdimensions = h->GetNdimensions();
1308   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1309
1310   AliCFContainer *dataContainer = GetContainer(kDataContainer);
1311   if(!dataContainer){
1312     AliError("Data Container not available");
1313     return NULL;
1314   }
1315   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1316   dataContainerbis->Add(dataContainerbis,-1.0);
1317
1318
1319   Int_t* coord = new Int_t[nbdimensions];
1320   memset(coord, 0, sizeof(Int_t) * nbdimensions);
1321   Double_t* points = new Double_t[nbdimensions];
1322
1323
1324   ULong64_t nEntries = h->GetNbins();
1325   for (ULong64_t i = 0; i < nEntries; ++i) {
1326     
1327     Double_t value = h->GetBinContent(i, coord);
1328     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1329     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1330     
1331     // Get the bin co-ordinates given an coord
1332     for (Int_t j = 0; j < nbdimensions; ++j)
1333       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1334
1335     if (!fEfficiencyFunction->IsInside(points))
1336          continue;
1337     TF1::RejectPoint(kFALSE);
1338
1339     // Evaulate function at points
1340     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1341     //printf("Value efficiency is %f\n",valueEfficiency);
1342
1343     if(valueEfficiency > 0.0) {
1344       h->SetBinContent(coord,value/valueEfficiency);
1345       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1346     }
1347     Double_t error = h->GetBinError(i);
1348     h->SetBinError(coord,error/valueEfficiency);
1349     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1350
1351    
1352   } 
1353
1354   delete[] coord;
1355   delete[] points;
1356
1357   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1358
1359   if(fDebugLevel > 0) {
1360     
1361     TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1362     cEfficiencyParametrized->Divide(2,1);
1363     cEfficiencyParametrized->cd(1);
1364     TH1D *afterE = (TH1D *) resultt->Project(0);
1365     TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1366     CorrectFromTheWidth(afterE);
1367     CorrectFromTheWidth(beforeE);
1368     afterE->SetStats(0);
1369     afterE->SetTitle("");
1370     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1371     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1372     afterE->SetMarkerStyle(25);
1373     afterE->SetMarkerColor(kBlack);
1374     afterE->SetLineColor(kBlack);
1375     beforeE->SetStats(0);
1376     beforeE->SetTitle("");
1377     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1378     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1379     beforeE->SetMarkerStyle(24);
1380     beforeE->SetMarkerColor(kBlue);
1381     beforeE->SetLineColor(kBlue);
1382     gPad->SetLogy();
1383     afterE->Draw();
1384     beforeE->Draw("same");
1385     TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1386     legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1387     legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1388     legefficiencyparametrized->Draw("same");
1389     cEfficiencyParametrized->cd(2);
1390     fEfficiencyFunction->Draw();
1391     //cEfficiencyParametrized->cd(3);
1392     //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1393     //ratioefficiency->Divide(afterE);
1394     //ratioefficiency->Draw();
1395
1396
1397   }
1398
1399   
1400   return resultt;
1401
1402 }
1403 //____________________________________________________________
1404 AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1405   
1406   //
1407   // Apply TPC pid efficiency correction from V0
1408   //
1409
1410   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1411   if(!v0Container){
1412     AliError("V0 Container not available");
1413     return NULL;
1414   }
1415
1416   // Efficiency
1417   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1418   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1419
1420   // Data in the right format
1421   AliCFDataGrid *dataGrid = 0x0;  
1422   if(bgsubpectrum) {
1423     dataGrid = bgsubpectrum;
1424   }
1425   else {
1426     
1427     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1428     if(!dataContainer){
1429       AliError("Data Container not available");
1430       return NULL;
1431     }
1432
1433     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1434   } 
1435
1436   // Correct
1437   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1438   result->ApplyEffCorrection(*efficiencyD);
1439
1440   if(fDebugLevel > 0) {
1441
1442     Int_t ptpr;
1443     if(fBeamType==0) ptpr=0;
1444     if(fBeamType==1) ptpr=1;
1445     
1446     TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1447     cV0Efficiency->Divide(2,1);
1448     cV0Efficiency->cd(1);
1449     TH1D *afterE = (TH1D *) result->Project(ptpr);
1450     TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
1451     afterE->SetStats(0);
1452     afterE->SetTitle("");
1453     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1454     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1455     afterE->SetMarkerStyle(25);
1456     afterE->SetMarkerColor(kBlack);
1457     afterE->SetLineColor(kBlack);
1458     beforeE->SetStats(0);
1459     beforeE->SetTitle("");
1460     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1461     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1462     beforeE->SetMarkerStyle(24);
1463     beforeE->SetMarkerColor(kBlue);
1464     beforeE->SetLineColor(kBlue);
1465     afterE->Draw();
1466     beforeE->Draw("same");
1467     TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1468     legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1469     legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1470     legV0efficiency->Draw("same");
1471     cV0Efficiency->cd(2);
1472     TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
1473     efficiencyDproj->SetTitle("");
1474     efficiencyDproj->SetStats(0);
1475     efficiencyDproj->SetMarkerStyle(25);
1476     efficiencyDproj->Draw();
1477
1478     if(fBeamType==1) {
1479
1480       TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1481       cV0Efficiencye->Divide(2,1);
1482       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1483       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1484      
1485       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1486       TAxis *cenaxisa = sparseafter->GetAxis(0);
1487       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1488       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1489       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1490       TAxis *cenaxisc = efficiencya->GetAxis(0);
1491       Int_t nbbin = cenaxisb->GetNbins();
1492       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1493       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1494       for(Int_t binc = 0; binc < nbbin; binc++){
1495         TString titlee("V0Efficiency_centrality_bin_");
1496         titlee += binc;
1497         TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1498         ccV0Efficiency->Divide(2,1);
1499         ccV0Efficiency->cd(1);
1500         gPad->SetLogy();
1501         cenaxisa->SetRange(binc+1,binc+1);
1502         cenaxisb->SetRange(binc+1,binc+1);
1503         cenaxisc->SetRange(binc+1,binc+1);
1504         TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1505         TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1506         CorrectFromTheWidth(aftere);
1507         CorrectFromTheWidth(beforee);
1508         aftere->SetStats(0);
1509         aftere->SetTitle((const char*)titlee);
1510         aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1511         aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1512         aftere->SetMarkerStyle(25);
1513         aftere->SetMarkerColor(kBlack);
1514         aftere->SetLineColor(kBlack);
1515         beforee->SetStats(0);
1516         beforee->SetTitle((const char*)titlee);
1517         beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1518         beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1519         beforee->SetMarkerStyle(24);
1520         beforee->SetMarkerColor(kBlue);
1521         beforee->SetLineColor(kBlue);
1522         aftere->Draw();
1523         beforee->Draw("same");
1524         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1525         lega->AddEntry(beforee,"Before correction","p");
1526         lega->AddEntry(aftere,"After correction","p");
1527         lega->Draw("same");
1528         cV0Efficiencye->cd(1);
1529         gPad->SetLogy();
1530         TH1D *afteree = (TH1D *) aftere->Clone();
1531         afteree->SetMarkerStyle(stylee[binc]);
1532         afteree->SetMarkerColor(colorr[binc]);
1533         if(binc==0) afteree->Draw();
1534         else afteree->Draw("same");
1535         legtotal->AddEntry(afteree,(const char*) titlee,"p");
1536         cV0Efficiencye->cd(2);
1537         gPad->SetLogy();
1538         TH1D *aftereeu = (TH1D *) aftere->Clone();
1539         aftereeu->SetMarkerStyle(stylee[binc]);
1540         aftereeu->SetMarkerColor(colorr[binc]);
1541         if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
1542         if(binc==0) aftereeu->Draw();
1543         else aftereeu->Draw("same");
1544         legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
1545         ccV0Efficiency->cd(2);
1546         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
1547         efficiencyDDproj->SetTitle("");
1548         efficiencyDDproj->SetStats(0);
1549         efficiencyDDproj->SetMarkerStyle(25);
1550         efficiencyDDproj->Draw();
1551                 
1552       }
1553      
1554       cV0Efficiencye->cd(1);
1555       legtotal->Draw("same");
1556       cV0Efficiencye->cd(2);
1557       legtotalg->Draw("same");
1558
1559       cenaxisa->SetRange(0,nbbin);
1560       cenaxisb->SetRange(0,nbbin);
1561       cenaxisc->SetRange(0,nbbin);
1562       
1563     }
1564
1565   }
1566
1567   
1568   return result;
1569
1570 }
1571 //____________________________________________________________
1572 TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1573   
1574   //
1575   // Unfold and eventually correct for efficiency the bgsubspectrum
1576   //
1577
1578   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1579   if(!mcContainer){
1580     AliError("MC Container not available");
1581     return NULL;
1582   }
1583
1584   if(!fCorrelation){
1585     AliError("No Correlation map available");
1586     return NULL;
1587   }
1588
1589   // Data 
1590   AliCFDataGrid *dataGrid = 0x0;  
1591   if(bgsubpectrum) {
1592     dataGrid = bgsubpectrum;
1593   }
1594   else {
1595
1596     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1597     if(!dataContainer){
1598       AliError("Data Container not available");
1599       return NULL;
1600     }
1601
1602     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1603   } 
1604   
1605   // Guessed
1606   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1607   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1608
1609   // Efficiency
1610   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1611   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1612   
1613   // Consider parameterized IP cut efficiency
1614   if(!fInclusiveSpectrum){
1615     Int_t* bins=new Int_t[1];
1616     bins[0]=44;
1617
1618     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1619     beffContainer->SetGrid(GetBeautyIPEff());
1620     efficiencyD->Multiply(beffContainer,1);  
1621   }
1622   
1623
1624   // Unfold 
1625   
1626   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
1627   //unfolding.SetUseCorrelatedErrors();
1628   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1629   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1630   if(fSetSmoothing) unfolding.UseSmoothing();
1631   unfolding.Unfold();
1632
1633   // Results
1634   THnSparse* result = unfolding.GetUnfolded();
1635   THnSparse* residual = unfolding.GetEstMeasured();
1636   TList *listofresults = new TList;
1637   listofresults->SetOwner();
1638   listofresults->AddAt((THnSparse*)result->Clone(),0);
1639   listofresults->AddAt((THnSparse*)residual->Clone(),1);
1640
1641   if(fDebugLevel > 0) {
1642
1643     Int_t ptpr;
1644     if(fBeamType==0) ptpr=0;
1645     if(fBeamType==1) ptpr=1;
1646     
1647     TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
1648     cresidual->Divide(2,1);
1649     cresidual->cd(1);
1650     gPad->SetLogy();
1651     TGraphErrors* residualspectrumD = Normalize(residual);
1652     if(!residualspectrumD) {
1653       AliError("Number of Events not set for the normalization");
1654       return kFALSE;
1655     }
1656     residualspectrumD->SetTitle("");
1657     residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
1658     residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1659     residualspectrumD->SetMarkerStyle(26);
1660     residualspectrumD->SetMarkerColor(kBlue);
1661     residualspectrumD->SetLineColor(kBlue);
1662     residualspectrumD->Draw("AP");
1663     AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
1664     dataGridBis->SetName("dataGridBis"); 
1665     TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
1666     measuredspectrumD->SetTitle("");  
1667     measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
1668     measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1669     measuredspectrumD->SetMarkerStyle(25);
1670     measuredspectrumD->SetMarkerColor(kBlack);
1671     measuredspectrumD->SetLineColor(kBlack);
1672     measuredspectrumD->Draw("P");
1673     TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
1674     legres->AddEntry(residualspectrumD,"Residual","p");
1675     legres->AddEntry(measuredspectrumD,"Measured","p");
1676     legres->Draw("same");
1677     cresidual->cd(2);
1678     TH1D *residualTH1D = residual->Projection(ptpr);
1679     TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
1680     TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
1681     ratioresidual->SetName("ratioresidual");
1682     ratioresidual->SetTitle("");
1683     ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
1684     ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1685     ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
1686     ratioresidual->SetStats(0);
1687     ratioresidual->Draw();
1688     
1689   }
1690
1691   return listofresults;
1692
1693 }
1694
1695 //____________________________________________________________
1696 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1697   
1698   //
1699   // Apply unfolding and efficiency correction together to bgsubspectrum
1700   //
1701
1702   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1703   if(!mcContainer){
1704     AliError("MC Container not available");
1705     return NULL;
1706   }
1707
1708   // Efficiency
1709   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1710   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1711
1712   // Consider parameterized IP cut efficiency
1713   if(!fInclusiveSpectrum){
1714     Int_t* bins=new Int_t[1];
1715     bins[0]=44;
1716   
1717     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1718     beffContainer->SetGrid(GetBeautyIPEff());
1719     efficiencyD->Multiply(beffContainer,1);
1720   }
1721
1722   // Data in the right format
1723   AliCFDataGrid *dataGrid = 0x0;  
1724   if(bgsubpectrum) {
1725     dataGrid = bgsubpectrum;
1726   }
1727   else {
1728     
1729     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1730     if(!dataContainer){
1731       AliError("Data Container not available");
1732       return NULL;
1733     }
1734
1735     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1736   } 
1737
1738   // Correct
1739   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1740   result->ApplyEffCorrection(*efficiencyD);
1741
1742   if(fDebugLevel > 0) {
1743
1744     Int_t ptpr;
1745     if(fBeamType==0) ptpr=0;
1746     if(fBeamType==1) ptpr=1;
1747     
1748     printf("Step MC: %d\n",fStepMC);
1749     printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1750     printf("Step MC true: %d\n",fStepTrue);
1751     AliCFEffGrid  *efficiencymcPID = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1752     AliCFEffGrid  *efficiencymctrackinggeo = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
1753     AliCFEffGrid  *efficiencymcall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
1754     
1755     AliCFEffGrid  *efficiencyesdall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
1756     
1757     TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
1758     cefficiency->cd(1);
1759     TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
1760     efficiencymcPIDD->SetTitle("");
1761     efficiencymcPIDD->SetStats(0);
1762     efficiencymcPIDD->SetMarkerStyle(25);
1763     efficiencymcPIDD->Draw();
1764     TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
1765     efficiencymctrackinggeoD->SetTitle("");
1766     efficiencymctrackinggeoD->SetStats(0);
1767     efficiencymctrackinggeoD->SetMarkerStyle(26);
1768     efficiencymctrackinggeoD->Draw("same");
1769     TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
1770     efficiencymcallD->SetTitle("");
1771     efficiencymcallD->SetStats(0);
1772     efficiencymcallD->SetMarkerStyle(27);
1773     efficiencymcallD->Draw("same");
1774     TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
1775     efficiencyesdallD->SetTitle("");
1776     efficiencyesdallD->SetStats(0);
1777     efficiencyesdallD->SetMarkerStyle(24);
1778     efficiencyesdallD->Draw("same");
1779     TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
1780     legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
1781     legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
1782     legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
1783     legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
1784     legeff->Draw("same");
1785
1786     if(fBeamType==1) {
1787
1788       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1789       TAxis *cenaxisa = sparseafter->GetAxis(0);
1790       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1791       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1792       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1793       TAxis *cenaxisc = efficiencya->GetAxis(0);
1794       //Int_t nbbin = cenaxisb->GetNbins();
1795       //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1796       //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1797       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
1798         TString titlee("Efficiency_centrality_bin_");
1799         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
1800         titlee += "_";
1801         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
1802         TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1803         cefficiencye->Divide(2,1);
1804         cefficiencye->cd(1);
1805         gPad->SetLogy();
1806         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1807         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1808         TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
1809         TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
1810         CorrectFromTheWidth(afterefficiency);
1811         CorrectFromTheWidth(beforeefficiency);
1812         afterefficiency->SetStats(0);
1813         afterefficiency->SetTitle((const char*)titlee);
1814         afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1815         afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1816         afterefficiency->SetMarkerStyle(25);
1817         afterefficiency->SetMarkerColor(kBlack);
1818         afterefficiency->SetLineColor(kBlack);
1819         beforeefficiency->SetStats(0);
1820         beforeefficiency->SetTitle((const char*)titlee);
1821         beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1822         beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1823         beforeefficiency->SetMarkerStyle(24);
1824         beforeefficiency->SetMarkerColor(kBlue);
1825         beforeefficiency->SetLineColor(kBlue);
1826         afterefficiency->Draw();
1827         beforeefficiency->Draw("same");
1828         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1829         lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
1830         lega->AddEntry(afterefficiency,"After efficiency correction","p");
1831         lega->Draw("same");
1832         cefficiencye->cd(2);
1833         cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
1834         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
1835         efficiencyDDproj->SetTitle("");
1836         efficiencyDDproj->SetStats(0);
1837         efficiencyDDproj->SetMarkerStyle(25);
1838         efficiencyDDproj->SetMarkerColor(2);
1839         efficiencyDDproj->Draw();
1840         cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
1841         TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
1842         efficiencyDDproja->SetTitle("");
1843         efficiencyDDproja->SetStats(0);
1844         efficiencyDDproja->SetMarkerStyle(26);
1845         efficiencyDDproja->SetMarkerColor(4);
1846         efficiencyDDproja->Draw("same");
1847         
1848       }
1849     }
1850
1851
1852   }
1853   
1854   return result;
1855
1856 }
1857
1858 //__________________________________________________________________________________
1859 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
1860   //
1861   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1862   // Give the final pt spectrum to be compared
1863   //
1864  
1865   if(fNEvents[i] > 0) {
1866
1867     Int_t ptpr = 0;
1868     if(fBeamType==0) ptpr=0;
1869     if(fBeamType==1) ptpr=1;
1870     
1871     TH1D* projection = spectrum->Projection(ptpr);
1872     CorrectFromTheWidth(projection);
1873     TGraphErrors *graphError = NormalizeTH1(projection,i);
1874     return graphError;
1875   
1876   }
1877     
1878   return 0x0;
1879   
1880
1881 }
1882 //__________________________________________________________________________________
1883 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
1884   //
1885   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1886   // Give the final pt spectrum to be compared
1887   //
1888   if(fNEvents[i] > 0) {
1889
1890     Int_t ptpr=0;
1891     if(fBeamType==0) ptpr=0;
1892     if(fBeamType==1) ptpr=1;
1893     
1894     TH1D* projection = (TH1D *) spectrum->Project(ptpr);
1895     CorrectFromTheWidth(projection);
1896     TGraphErrors *graphError = NormalizeTH1(projection,i);
1897     return graphError;
1898     
1899   }
1900
1901   return 0x0;
1902   
1903
1904 }
1905 //__________________________________________________________________________________
1906 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
1907   //
1908   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1909   // Give the final pt spectrum to be compared
1910   //
1911   Double_t chargecoefficient = 0.5;
1912   if(fChargeChoosen >= 0) chargecoefficient = 1.0;
1913
1914   Double_t etarange = fEtaSelected ? fEtaRange[1] - fEtaRange[0] : 1.6;
1915   printf("Normalizing Eta Range %f\n", etarange);
1916   if(fNEvents[i] > 0) {
1917
1918     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1919     Double_t p = 0, dp = 0; Int_t point = 1;
1920     Double_t n = 0, dN = 0;
1921     Double_t nCorr = 0, dNcorr = 0;
1922     Double_t errdN = 0, errdp = 0;
1923     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1924       point = ibin - input->GetXaxis()->GetFirst();
1925       p = input->GetXaxis()->GetBinCenter(ibin);
1926       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1927       n = input->GetBinContent(ibin);
1928       dN = input->GetBinError(ibin);
1929
1930       // New point
1931       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
1932       errdN = 1./(2. * TMath::Pi() * p);
1933       errdp = 1./(2. * TMath::Pi() * p*p) * n;
1934       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
1935       
1936       spectrumNormalized->SetPoint(point, p, nCorr);
1937       spectrumNormalized->SetPointError(point, dp, dNcorr);
1938     }
1939     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1940     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1941     spectrumNormalized->SetMarkerStyle(22);
1942     spectrumNormalized->SetMarkerColor(kBlue);
1943     spectrumNormalized->SetLineColor(kBlue);
1944     
1945     return spectrumNormalized;
1946     
1947   }
1948
1949   return 0x0;
1950   
1951
1952 }
1953 //__________________________________________________________________________________
1954 TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
1955   //
1956   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1957   // Give the final pt spectrum to be compared
1958   //
1959   Double_t chargecoefficient = 0.5;
1960   if(fChargeChoosen >= 0) chargecoefficient = 1.0;
1961
1962   Double_t etarange = fEtaSelected ? fEtaRange[1] - fEtaRange[0] : 1.6;
1963   printf("Normalizing Eta Range %f\n", etarange);
1964   if(normalization > 0) {
1965
1966     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1967     Double_t p = 0, dp = 0; Int_t point = 1;
1968     Double_t n = 0, dN = 0;
1969     Double_t nCorr = 0, dNcorr = 0;
1970     Double_t errdN = 0, errdp = 0;
1971     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1972       point = ibin - input->GetXaxis()->GetFirst();
1973       p = input->GetXaxis()->GetBinCenter(ibin);
1974       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1975       n = input->GetBinContent(ibin);
1976       dN = input->GetBinError(ibin);
1977
1978       // New point
1979       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
1980       errdN = 1./(2. * TMath::Pi() * p);
1981       errdp = 1./(2. * TMath::Pi() * p*p) * n;
1982       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
1983       
1984       spectrumNormalized->SetPoint(point, p, nCorr);
1985       spectrumNormalized->SetPointError(point, dp, dNcorr);
1986     }
1987     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1988     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
1989     spectrumNormalized->SetMarkerStyle(22);
1990     spectrumNormalized->SetMarkerColor(kBlue);
1991     spectrumNormalized->SetLineColor(kBlue);
1992     
1993     return spectrumNormalized;
1994     
1995   }
1996
1997   return 0x0;
1998   
1999
2000 }
2001 //____________________________________________________________
2002 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2003   //
2004   // Set the container for a given step to the 
2005   //
2006   if(!fCFContainers) fCFContainers = new TList;
2007   fCFContainers->AddAt(cont, type);
2008 }
2009
2010 //____________________________________________________________
2011 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2012   //
2013   // Get Correction Framework Container for given type
2014   //
2015   if(!fCFContainers) return NULL;
2016   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2017 }
2018 //____________________________________________________________
2019 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Int_t positivenegative) {
2020   //
2021   // Slice bin for a given source of electron
2022   // nDim is the number of dimension the corrections are done
2023   // dimensions are the definition of the dimensions
2024   // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2025   // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
2026   //
2027   
2028   Double_t *varMin = new Double_t[container->GetNVar()],
2029            *varMax = new Double_t[container->GetNVar()];
2030
2031   Double_t *binLimits;
2032   for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2033     
2034     binLimits = new Double_t[container->GetNBins(ivar)+1];
2035     container->GetBinLimits(ivar,binLimits);
2036     varMin[ivar] = binLimits[0];
2037     varMax[ivar] = binLimits[container->GetNBins(ivar)];
2038     // source
2039     if(ivar == 4){
2040       if((source>= 0) && (source<container->GetNBins(ivar))) {
2041         varMin[ivar] = binLimits[source];
2042         varMax[ivar] = binLimits[source];
2043       }     
2044     }
2045     // charge
2046     if(ivar == 3) {
2047       if((positivenegative>= 0) && (positivenegative<container->GetNBins(ivar))) {
2048         varMin[ivar] = binLimits[positivenegative];
2049         varMax[ivar] = binLimits[positivenegative];
2050       }
2051     }
2052     // eta
2053     if(ivar == 1){
2054       for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic))); 
2055       if(fEtaSelected){
2056         varMin[ivar] = fEtaRange[0];
2057         varMax[ivar] = fEtaRange[1];
2058       }
2059     }
2060     
2061
2062     delete[] binLimits;
2063   }
2064   
2065   AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2066   AddTemporaryObject(k);
2067   delete[] varMin; delete[] varMax;
2068
2069   return k;
2070
2071 }
2072
2073 //_________________________________________________________________________
2074 THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
2075   //
2076   // Slice correlation
2077   //
2078
2079   Int_t ndimensions = correlationmatrix->GetNdimensions();
2080   //printf("Number of dimension %d correlation map\n",ndimensions);
2081   if(ndimensions < (2*nDim)) {
2082     AliError("Problem in the dimensions");
2083     return NULL;
2084   }
2085   Int_t ndimensionsContainer = (Int_t) ndimensions/2;
2086   //printf("Number of dimension %d container\n",ndimensionsContainer);
2087
2088   Int_t *dim = new Int_t[nDim*2];
2089   for(Int_t iter=0; iter < nDim; iter++){
2090     dim[iter] = dimensions[iter];
2091     dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
2092     //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
2093   }
2094     
2095   THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
2096
2097   delete[] dim; 
2098   return k;
2099   
2100 }
2101 //___________________________________________________________________________
2102 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2103   //
2104   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2105   //
2106
2107   TAxis *axis = h1->GetXaxis();
2108   Int_t nbinX = h1->GetNbinsX();
2109
2110   for(Int_t i = 1; i <= nbinX; i++) {
2111
2112     Double_t width = axis->GetBinWidth(i);
2113     Double_t content = h1->GetBinContent(i);
2114     Double_t error = h1->GetBinError(i); 
2115     h1->SetBinContent(i,content/width); 
2116     h1->SetBinError(i,error/width);
2117   }
2118
2119 }
2120
2121 //___________________________________________________________________________
2122 void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const { 
2123   //
2124   // Correct statistical error
2125   //
2126
2127   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2128   Int_t nbinX = h1->GetNbinsX();
2129   Int_t bins[1];
2130   for(Long_t i = 1; i <= nbinX; i++) {
2131     bins[0] = i;
2132     Float_t content = h1->GetBinContent(i);
2133     if(content>0){
2134       Float_t error = TMath::Sqrt(content);
2135       backgroundGrid->SetElementError(bins, error);
2136     }
2137   }
2138 }
2139
2140 //____________________________________________________________
2141 void AliHFEspectrum::AddTemporaryObject(TObject *o){
2142   // 
2143   // Emulate garbage collection on class level
2144   // 
2145   if(!fTemporaryObjects) {
2146     fTemporaryObjects= new TList;
2147     fTemporaryObjects->SetOwner();
2148   }
2149   fTemporaryObjects->Add(o);
2150 }
2151
2152 //____________________________________________________________
2153 void AliHFEspectrum::ClearObject(TObject *o){
2154   //
2155   // Do a safe deletion
2156   //
2157   if(fTemporaryObjects){
2158     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2159     delete o;
2160   } else{
2161     // just delete
2162     delete o;
2163   }
2164 }
2165 //_________________________________________________________________________
2166 TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
2167   AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
2168   return data;
2169 }
2170 //_________________________________________________________________________
2171 TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
2172   // 
2173   // Create efficiency grid and calculate efficiency
2174   // of step to step0
2175   //
2176   TString name("eff");
2177   name += step;
2178   name+= step0;
2179   AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2180   eff->CalculateEfficiency(step,step0);
2181   return eff;
2182 }
2183
2184 //____________________________________________________________________________
2185 THnSparse* AliHFEspectrum::GetCharmWeights(){
2186   
2187   //
2188   // Measured D->e based weighting factors
2189   //
2190
2191   const Int_t nDim=1;
2192   Int_t nBin[nDim] = {44};
2193   const Double_t kPtbound[2] = {0.1, 20.};
2194
2195   Double_t* binEdges[nDim];
2196   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2197
2198   fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
2199   for(Int_t idim = 0; idim < nDim; idim++)
2200      fWeightCharm->SetBinEdges(idim, binEdges[idim]);
2201      Double_t weight[44]={
2202 1.11865,1.11444,1.13395,1.15751,1.13412,1.14446,1.15372,1.09458,1.13446,1.11707,1.1352,1.10979,1.08887,1.11164,1.10772,1.10727,1.07027,1.07016,1.04826,1.03248,1.0093,0.973534,0.932574,0.869838,0.799316,0.734015,0.643001,0.584319,0.527147,0.495661,0.470002,0.441129,0.431156,0.417242,0.39246,0.379611,0.390845,0.368857,0.34959,0.387186,0.376364,0.384437,0.349585,0.383225};
2203   //points
2204   Double_t pt[1];
2205   for(int i=0; i<nBin[0]; i++){
2206     pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2207     fWeightCharm->Fill(pt,weight[i]);
2208   }
2209   Int_t* ibins[nDim];
2210   for(Int_t ivar = 0; ivar < nDim; ivar++)
2211     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2212   fWeightCharm->SetBinError(ibins[0],0);
2213
2214   return fWeightCharm;
2215 }
2216
2217 //____________________________________________________________________________
2218 THnSparse* AliHFEspectrum::GetBeautyIPEff(){
2219   //
2220   // Return beauty electron IP cut efficiency
2221   //
2222
2223   const Int_t nDim=1;
2224   Int_t nBin[nDim] = {44};
2225   const Double_t kPtbound[2] = {0.1, 20.};
2226
2227   Double_t* binEdges[nDim];
2228   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2229
2230   THnSparseF *ipcut = new THnSparseF("beff", "b eff; pt[GeV/c]", nDim, nBin);
2231   for(Int_t idim = 0; idim < nDim; idim++)
2232      ipcut->SetBinEdges(idim, binEdges[idim]);
2233   Double_t pt[1];
2234   Double_t weight;
2235   for(int i=0; i<nBin[0]; i++){
2236     pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2237     weight=0.612*(0.385+0.111*log(pt[0])-0.00435*log(pt[0])*log(pt[0]))*tanh(5.42*pt[0]-1.22);  // for 3 sigma cut   
2238     //weight=0.786*(0.493+0.0283*log(pt[0])-0.0190*log(pt[0])*log(pt[0]))*tanh(1.84*pt[0]-0.00368);  // for 2 sigma cut   
2239     ipcut->Fill(pt,weight);
2240   }
2241   Int_t* ibins[nDim];
2242   for(Int_t ivar = 0; ivar < nDim; ivar++)
2243     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2244   ipcut->SetBinError(ibins[0],0);
2245
2246   return ipcut;
2247 }
2248
2249 //____________________________________________________________________________
2250 THnSparse* AliHFEspectrum::GetCharmEff(){
2251   //
2252   // Return charm electron IP cut efficiency
2253   //
2254
2255   const Int_t nDim=1;
2256   Int_t nBin[nDim] = {44};
2257   const Double_t kPtbound[2] = {0.1, 20.};
2258
2259   Double_t* binEdges[nDim];
2260   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2261
2262   THnSparseF *ipcut = new THnSparseF("ceff", "c eff; pt[GeV/c]", nDim, nBin);
2263   for(Int_t idim = 0; idim < nDim; idim++)
2264      ipcut->SetBinEdges(idim, binEdges[idim]);
2265   Double_t pt[1];
2266   Double_t weight;
2267   for(int i=0; i<nBin[0]; i++){
2268     pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2269     weight=0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
2270     //weight=0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
2271     ipcut->Fill(pt,weight);
2272   }
2273   Int_t* ibins[nDim];
2274   for(Int_t ivar = 0; ivar < nDim; ivar++)
2275     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2276   ipcut->SetBinError(ibins[0],0);
2277
2278   return ipcut;
2279 }
2280
2281 //____________________________________________________________________________
2282 THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
2283   //
2284   // Return PID x IP cut efficiency
2285   //
2286
2287   const Int_t nDim=1;
2288   Int_t nBin[nDim] = {44};
2289   const Double_t kPtbound[2] = {0.1, 20.};
2290
2291   Double_t* binEdges[nDim];
2292   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
2293
2294   THnSparseF *pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
2295   for(Int_t idim = 0; idim < nDim; idim++)
2296      pideff->SetBinEdges(idim, binEdges[idim]);
2297
2298   Double_t pt[1];
2299   Double_t weight = 1.0;
2300   Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
2301   for(int i=0; i<nBin[0]; i++){
2302     pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
2303     Double_t tofeff = 0.9885*(0.8551+0.0070*log(pt[0])-0.0014*log(pt[0])*log(pt[0]))*tanh(0.8884*pt[0]+0.6061); // parameterized TOF PID eff
2304
2305     if(source==0) weight = tofeff*trdtpcPidEfficiency*0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
2306     //if(source==0) weight = tofeff*trdtpcPidEfficiency*0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
2307     //if(source==2) weight = tofeff*trdtpcPidEfficiency*0.5575*(0.3594-0.3051*log(pt[0])+0.1597*log(pt[0])*log(pt[0]))*tanh(8.2436*pt[0]-0.8125);
2308     //if(source==3) weight = tofeff*trdtpcPidEfficiency*0.0937*(0.4449-0.3769*log(pt[0])+0.5031*log(pt[0])*log(pt[0]))*tanh(6.4063*pt[0]+3.1425); // multiply ip cut eff for Dalitz/dielectrons
2309     if(source==2) weight = tofeff*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // multiply ip cut eff for conversion electrons
2310     if(source==3) weight = tofeff*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // multiply ip cut eff for Dalitz/dielectrons
2311     printf("tofeff= %lf trdtpcPidEfficiency= %lf conversion= %lf nonhfe= %lf\n",tofeff,trdtpcPidEfficiency,fConversionEff->GetBinContent(i+1),fNonHFEEff->GetBinContent(i+1));
2312
2313     pideff->Fill(pt,weight);
2314   }
2315   Int_t* ibins[nDim];
2316   for(Int_t ivar = 0; ivar < nDim; ivar++)
2317     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2318   pideff->SetBinError(ibins[0],0);
2319
2320   return pideff;
2321 }