]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGHF/hfe/AliHFEBeautySpectrum.cxx
Update HFE code
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEBeautySpectrum.cxx
1  
2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 *                                                                        *
5 * Author: The ALICE Off-line Project.                                    *
6 * Contributors are mentioned in the code where appropriate.              *
7 *                                                                        *
8 * Permission to use, copy, modify and distribute this software and its   *
9 * documentation strictly for non-commercial purposes is hereby granted   *
10 * without fee, provided that the above copyright notice appears in all   *
11 * copies and that both the copyright notice and this permission notice   *
12 * appear in the supporting documentation. The authors make no claims     *
13 * about the suitability of this software for any purpose. It is          *
14 * provided "as is" without express or implied warranty.                  *
15 **************************************************************************/
16 //
17 // Class for spectrum correction
18 // Subtraction of hadronic background, Unfolding of the data and
19 // Renormalization done here
20 // The following containers have to be set:
21 //  - Correction framework container for real data
22 //  - Correction framework container for MC (Efficiency Map)
23 //  - Correction framework container for background coming from data
24 //  - Correction framework container for background coming from MC
25 //
26 //  Author: 
27 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
28 //            Markus Fasel <M.Fasel@gsi.de>
29 //
30 #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 "AliHFEBeautySpectrum.h"
55 #include "AliHFEBeautySpectrumQA.h"
56 #include "AliHFEcuts.h"
57 #include "AliHFEcontainer.h"
58 #include "AliHFEtools.h"
59
60 ClassImp(AliHFEBeautySpectrum)
61
62 //____________________________________________________________
63 AliHFEBeautySpectrum::AliHFEBeautySpectrum(const char *name):
64 AliHFECorrectSpectrumBase(name),
65   fQA(NULL),
66   fTemporaryObjects(NULL),
67   fBackground(NULL),
68   fEfficiencyTOFPIDD(NULL),
69   fEfficiencyesdTOFPIDD(NULL),
70   fEfficiencyIPCharmD(NULL),
71   fEfficiencyIPBeautyD(NULL),
72   fEfficiencyIPBeautyesdD(NULL),
73   fEfficiencyIPConversionD(NULL),
74   fEfficiencyIPNonhfeD(NULL), 
75   fNonHFEbg(NULL),
76   fWeightCharm(NULL),
77   fInclusiveSpectrum(kFALSE),
78   fDumpToFile(kFALSE),
79   fUnSetCorrelatedErrors(kTRUE),
80   fIPanaHadronBgSubtract(kFALSE),
81   fIPanaCharmBgSubtract(kFALSE),
82   fIPanaNonHFEBgSubtract(kFALSE),
83   fIPParameterizedEff(kFALSE),
84   fNonHFEsyst(0),
85   fBeauty2ndMethod(kFALSE),
86   fIPEffCombinedSamples(kTRUE),
87   fNMCEvents(0),
88   fNMCbgEvents(0),
89   fNCentralityBinAtTheEnd(0),
90   fFillMoreCorrelationMatrix(kFALSE),
91   fHadronEffbyIPcut(NULL),
92   fEfficiencyCharmSigD(NULL),
93   fEfficiencyBeautySigD(NULL),
94   fEfficiencyBeautySigesdD(NULL),
95   fConversionEff(NULL),
96   fNonHFEEff(NULL),
97   fCharmEff(NULL),
98   fBeautyEff(NULL),
99   fConversionEffbgc(NULL),
100   fNonHFEEffbgc(NULL),      
101   fBSpectrum2ndMethod(NULL),
102   fkBeauty2ndMethodfilename(""),
103   fBeamType(0),
104   fEtaSyst(kTRUE),
105   fDebugLevel(0),
106   fWriteToFile(kFALSE),
107   fUnfoldBG(kFALSE)
108 {
109   //
110   // Default constructor
111   //
112
113   fQA = new AliHFEBeautySpectrumQA("AliHFEBeautySpectrumQA");
114
115   for(Int_t k = 0; k < 20; k++){
116     fLowBoundaryCentralityBinAtTheEnd[k] = 0;
117     fHighBoundaryCentralityBinAtTheEnd[k] = 0;
118   }
119   fNMCbgEvents = 0;  
120   fEfficiencyTOFPIDD = 0;
121   fEfficiencyesdTOFPIDD = 0;
122   fEfficiencyIPCharmD = 0;     
123   fEfficiencyIPBeautyD = 0;    
124   fEfficiencyIPBeautyesdD = 0;
125   fEfficiencyIPConversionD = 0;
126   fEfficiencyIPNonhfeD = 0;   
127   
128   fConversionEff = 0;
129   fNonHFEEff = 0;
130   fCharmEff = 0;
131   fBeautyEff = 0;
132   fEfficiencyCharmSigD = 0;
133   fEfficiencyBeautySigD = 0;
134   fEfficiencyBeautySigesdD = 0;
135   
136   
137   memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
138   memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
139 }
140 //____________________________________________________________
141 AliHFEBeautySpectrum::AliHFEBeautySpectrum(const AliHFEBeautySpectrum &ref):
142   AliHFECorrectSpectrumBase(ref),
143   fQA(ref.fQA),
144   fTemporaryObjects(ref.fTemporaryObjects),
145   fBackground(ref.fBackground),
146   fEfficiencyTOFPIDD(ref.fEfficiencyTOFPIDD),
147   fEfficiencyesdTOFPIDD(ref.fEfficiencyesdTOFPIDD),
148   fEfficiencyIPCharmD(ref.fEfficiencyIPCharmD),
149   fEfficiencyIPBeautyD(ref.fEfficiencyIPBeautyD),
150   fEfficiencyIPBeautyesdD(ref.fEfficiencyIPBeautyesdD),
151   fEfficiencyIPConversionD(ref.fEfficiencyIPConversionD),
152   fEfficiencyIPNonhfeD(ref.fEfficiencyIPNonhfeD), 
153   fNonHFEbg(ref.fNonHFEbg),
154   fWeightCharm(ref.fWeightCharm),
155   fInclusiveSpectrum(ref.fInclusiveSpectrum),
156   fDumpToFile(ref.fDumpToFile),
157   fUnSetCorrelatedErrors(ref.fUnSetCorrelatedErrors),
158   fIPanaHadronBgSubtract(ref.fIPanaHadronBgSubtract),
159   fIPanaCharmBgSubtract(ref.fIPanaCharmBgSubtract),
160   fIPanaNonHFEBgSubtract(ref.fIPanaNonHFEBgSubtract),
161   fIPParameterizedEff(ref.fIPParameterizedEff),
162   fNonHFEsyst(ref.fNonHFEsyst),
163   fBeauty2ndMethod(ref.fBeauty2ndMethod),
164   fIPEffCombinedSamples(ref.fIPEffCombinedSamples),
165   fNMCEvents(ref.fNMCEvents),
166   fNMCbgEvents(ref.fNMCbgEvents),
167   fNCentralityBinAtTheEnd(ref.fNCentralityBinAtTheEnd),
168   fFillMoreCorrelationMatrix(ref.fFillMoreCorrelationMatrix),
169   fHadronEffbyIPcut(ref.fHadronEffbyIPcut),
170   fEfficiencyCharmSigD(ref.fEfficiencyCharmSigD),
171   fEfficiencyBeautySigD(ref.fEfficiencyBeautySigD),
172   fEfficiencyBeautySigesdD(ref.fEfficiencyBeautySigesdD),
173   fConversionEff(ref.fConversionEff),
174   fNonHFEEff(ref.fNonHFEEff),
175   fCharmEff(ref.fCharmEff),
176   fBeautyEff(ref.fBeautyEff),
177   fConversionEffbgc(ref.fConversionEffbgc),
178   fNonHFEEffbgc(ref.fNonHFEEffbgc),      
179   fBSpectrum2ndMethod(ref.fBSpectrum2ndMethod),
180   fkBeauty2ndMethodfilename(ref.fkBeauty2ndMethodfilename),
181   fBeamType(ref.fBeamType),
182   fEtaSyst(ref.fEtaSyst),
183   fDebugLevel(ref.fDebugLevel),
184   fWriteToFile(ref.fWriteToFile),
185   fUnfoldBG(ref.fUnfoldBG)
186 {
187   //
188   // Copy constructor
189   //
190
191   ref.Copy(*this);  
192 }
193 //____________________________________________________________
194 AliHFEBeautySpectrum &AliHFEBeautySpectrum::operator=(const AliHFEBeautySpectrum &ref){
195  //
196   // Assignment operator
197   //
198   if(this == &ref) 
199     ref.Copy(*this);
200   return *this;
201 }
202 //____________________________________________________________
203 void AliHFEBeautySpectrum::Copy(TObject &o) const {
204   // 
205   // Copy into object o
206   //
207   AliHFEBeautySpectrum &target = dynamic_cast<AliHFEBeautySpectrum &>(o);
208   target.fQA = fQA;
209
210   target.fQA = fQA;
211   target.fTemporaryObjects = fTemporaryObjects;
212   target.fBackground = fBackground;
213   target.fWeightCharm = fWeightCharm;
214   target.fDumpToFile = fDumpToFile;
215   target.fUnSetCorrelatedErrors = fUnSetCorrelatedErrors;
216   target.fIPanaHadronBgSubtract = fIPanaHadronBgSubtract;
217   target.fIPanaCharmBgSubtract = fIPanaCharmBgSubtract;
218   target.fIPanaNonHFEBgSubtract = fIPanaNonHFEBgSubtract;
219   target.fIPParameterizedEff = fIPParameterizedEff;
220   target.fNonHFEsyst = fNonHFEsyst;
221   target.fBeauty2ndMethod = fBeauty2ndMethod;
222   target.fIPEffCombinedSamples = fIPEffCombinedSamples;
223   target.fNMCEvents = fNMCEvents;
224   target.fNCentralityBinAtTheEnd = fNCentralityBinAtTheEnd;
225   target.fFillMoreCorrelationMatrix = fFillMoreCorrelationMatrix;
226   target.fHadronEffbyIPcut = fHadronEffbyIPcut;
227   target.fConversionEffbgc = fConversionEffbgc;
228   target.fNonHFEEffbgc = fNonHFEEffbgc;      
229   target.fBSpectrum2ndMethod = fBSpectrum2ndMethod;
230   target.fkBeauty2ndMethodfilename = fkBeauty2ndMethodfilename;
231   target.fBeamType = fBeamType;
232   target.fEtaSyst = fEtaSyst;
233   target.fDebugLevel = fDebugLevel;
234   target.fWriteToFile = fWriteToFile;
235   target.fUnfoldBG = fUnfoldBG;
236 }
237 //____________________________________________________________
238 AliHFEBeautySpectrum::~AliHFEBeautySpectrum(){
239   //
240   // Destructor
241   //
242   if(fTemporaryObjects){
243     fTemporaryObjects->Clear();
244     delete fTemporaryObjects;
245   }
246   if(fQA) delete fQA;
247 }
248 //____________________________________________________________
249 Bool_t AliHFEBeautySpectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *bghfecontainer, const AliHFEcontainer */*v0hfecontainer*/,AliCFContainer */*photoniccontainerD*/){
250   //
251   // Init what we need for the correction:
252   //
253   // Raw spectrum, hadron contamination
254   // MC efficiency maps, correlation matrix
255   //
256   // This for a given dimension.
257   // If no inclusive spectrum, then take only efficiency map for beauty electron
258   // and the appropriate correlation matrix
259   //
260
261   
262   if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
263
264   Int_t kNdim = 3;
265  
266   Int_t dims[kNdim];
267   
268   switch(fNbDimensions){
269   case 1:   dims[0] = 0;
270     break;
271   case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
272     break;
273   case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
274     break;
275   default:
276     AliError("Container with this number of dimensions not foreseen (yet)");
277     return kFALSE;
278   };
279   
280   //
281   // Data container: raw spectrum + hadron contamination  
282   //
283   AliCFContainer *datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
284   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
285   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
286   AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
287   AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
288   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
289   SetContainer(datacontainerD,AliHFECorrectSpectrumBase::kDataContainer);
290   SetContainer(contaminationcontainerD,AliHFECorrectSpectrumBase::kBackgroundData);
291
292   // QA : see with MJ if it is necessary for Beauty
293   Int_t dimqa = datacontainer->GetNVar();
294   Int_t dimsqa[dimqa];
295   for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
296   AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
297   fQA->AddResultAt(datacontainerDQA,AliHFEBeautySpectrumQA::kDataProjection);
298
299   //
300   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
301   // 
302   AliCFContainer *mccontaineresd = 0x0;
303   AliCFContainer *mccontaineresdbg = 0x0;
304   AliCFContainer *mccontainermc = 0x0;
305   AliCFContainer *mccontainermcbg = 0x0;
306   AliCFContainer *nonHFEweightedContainer = 0x0;
307   AliCFContainer *convweightedContainer = 0x0;
308   AliCFContainer *nonHFEweightedContainerSig = 0x0;//mjmj
309   AliCFContainer *convweightedContainerSig = 0x0;//mjmj
310   AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
311   AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
312    
313   mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
314   mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
315   mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
316   mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
317   
318   if(fNonHFEsyst){   
319     const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
320     const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
321     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
322       for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
323         if(!(nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]))))continue;
324         convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));     
325         if(!(fConvSourceContainer[iSource][iLevel] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
326         if(!(fNonHFESourceContainer[iSource][iLevel] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
327         //            if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
328         if(fBeamType == 1)break;//no systematics yet for PbPb non-HFE backgrounds
329       }
330     }
331   }
332   // else{      
333   nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
334   convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
335   if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE; 
336   nonHFEweightedContainerSig = mchfecontainer->GetCFContainer("mesonElecs");//mjmj
337   convweightedContainerSig = mchfecontainer->GetCFContainer("conversionElecs");//mjmj
338   if((!convweightedContainerSig)||(!nonHFEweightedContainerSig)) return kFALSE;
339   //}
340   
341   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;  
342   
343   Int_t source = 1; //beauty
344   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
345   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
346   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
347   SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerMC);
348   SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
349
350   // set charm, nonHFE container to estimate BG
351   source = 0; //charm
352   mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); //mjmjmj
353   //mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
354   SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerCharmMC);
355
356   //if(!fNonHFEsyst){
357   AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
358   SetContainer(nonHFEweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESD);
359   AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
360   SetContainer(convweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESD);
361   
362   //mjmj
363   AliCFContainer *nonHFEweightedContainerSigD = GetSlicedContainer(nonHFEweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
364   SetContainer(nonHFEweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
365   if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig)){printf("merged non-HFE container not found!\n");return kFALSE;};
366   AliCFContainer *convweightedContainerSigD = GetSlicedContainer(convweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
367   SetContainer(convweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig); 
368   if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig)){printf(" merged conversion container not found!\n");return kFALSE;}; 
369   //}
370   
371   SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
372   
373   //
374   // MC container: correlation matrix
375   //
376   THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
377   //mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
378   if(!mccorrelation) return kFALSE;
379   Int_t testCentralityLow = fTestCentralityLow;
380   Int_t testCentralityHigh = fTestCentralityHigh;
381   if(fFillMoreCorrelationMatrix) {
382     testCentralityLow = fTestCentralityLow-1;
383     testCentralityHigh = fTestCentralityHigh+1;
384   }
385   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims, fChargeChoosen, testCentralityLow,testCentralityHigh);
386   if(!mccorrelationD) {
387     printf("No correlation\n");
388     return kFALSE;
389   }
390   SetCorrelation(mccorrelationD);
391
392   // QA
393   fQA->AddResultAt(mccorrelation,AliHFEBeautySpectrumQA::kCMProjection); 
394   //
395   fQA->DrawProjections();
396
397
398   return kTRUE;
399 }
400
401
402 //____________________________________________________________
403 void AliHFEBeautySpectrum::CallInputFileForBeauty2ndMethod(){
404   //
405   // get spectrum for beauty 2nd method
406   //
407   //
408     TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
409     fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
410 }
411
412 //____________________________________________________________
413 Bool_t AliHFEBeautySpectrum::Correct(Bool_t subtractcontamination, Bool_t /*subtractphotonic*/){
414   //
415   // Correct the spectrum for efficiency and unfolding for beauty analysis
416   // with both method and compare
417   //
418   
419   gStyle->SetPalette(1);
420   gStyle->SetOptStat(1111);
421   gStyle->SetPadBorderMode(0);
422   gStyle->SetCanvasColor(10);
423   gStyle->SetPadLeftMargin(0.13);
424   gStyle->SetPadRightMargin(0.13);
425
426   printf("Steps are: stepdata %d, stepMC %d, steptrue %d\n",fStepData,fStepMC,fStepTrue);
427
428   ///////////////////////////
429   // Check initialization
430   ///////////////////////////
431
432   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
433     AliInfo("You have to init before");
434     return kFALSE;
435   }
436   
437   if((fStepTrue <= 0) && (fStepMC <= 0) && (fStepData <= 0)) {
438     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
439     return kFALSE;
440   }
441  
442   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
443     
444   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
445   //////////////////////////////////
446   // Subtract hadron background
447   /////////////////////////////////
448   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
449   AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
450   TGraphErrors *gNormalizedRawSpectrum = 0x0;
451   if(subtractcontamination) {
452       if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
453       else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
454       unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
455       dataGridAfterFirstSteps = dataspectrumaftersubstraction;
456       gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
457   }
458
459   /////////////////////////////////////////////////////////////////////////////////////////
460   // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
461   /////////////////////////////////////////////////////////////////////////////////////////
462
463   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
464   
465   if(fEfficiencyFunction){
466     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
467     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
468   }
469  
470   ///////////////
471   // Unfold
472   //////////////
473   THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
474   if(!correctedspectrum){
475     AliError("No corrected spectrum\n");
476     return kFALSE;
477   }
478   
479   /////////////////////
480   // Simply correct
481   ////////////////////
482
483   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
484   
485
486   //////////
487   // Plot
488   //////////
489
490   // QA final results
491   TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
492   TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
493   fQA->AddResultAt(correctedspectrumD,AliHFEBeautySpectrumQA::kFinalResultUnfolded);
494   fQA->AddResultAt(alltogetherspectrumD,AliHFEBeautySpectrumQA::kFinalResultDirectEfficiency);
495   fQA->AddResultAt(correctedspectrum,AliHFEBeautySpectrumQA::kFinalResultUnfSparse);
496   fQA->AddResultAt(alltogetherCorrection,AliHFEBeautySpectrumQA::kFinalResultDirectEffSparse);
497   fQA->DrawResult();
498
499  
500   if(fBeamType == 0){
501     if(fNonHFEsyst){
502       CalculateNonHFEsyst();
503     }
504   }
505
506   // Dump to file if needed
507   
508   if(fDumpToFile) {
509     TFile *out;
510     out = new TFile("finalSpectrum.root","update");
511     out->cd();
512     // to do centrality dependent
513     correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
514     correctedspectrum->Write();
515     alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
516     alltogetherCorrection->Write();
517     //
518     if(unnormalizedRawSpectrum) {
519       unnormalizedRawSpectrum->SetName("beautyAfterIP");
520       unnormalizedRawSpectrum->Write();
521     }
522     
523     if(gNormalizedRawSpectrum){
524       gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
525       gNormalizedRawSpectrum->Write();
526     }
527         
528     fEfficiencyCharmSigD->SetTitle("IPEfficiencyForCharmSig");
529     fEfficiencyCharmSigD->SetName("IPEfficiencyForCharmSig");
530     fEfficiencyCharmSigD->Write();
531     fEfficiencyBeautySigD->SetTitle("IPEfficiencyForBeautySig");
532     fEfficiencyBeautySigD->SetName("IPEfficiencyForBeautySig");
533     fEfficiencyBeautySigD->Write();
534     fCharmEff->SetTitle("IPEfficiencyForCharm");
535     fCharmEff->SetName("IPEfficiencyForCharm");
536     fCharmEff->Write();
537     fBeautyEff->SetTitle("IPEfficiencyForBeauty");
538     fBeautyEff->SetName("IPEfficiencyForBeauty");
539     fBeautyEff->Write();
540     fConversionEff->SetTitle("IPEfficiencyForConversion");
541     fConversionEff->SetName("IPEfficiencyForConversion");
542     fConversionEff->Write();
543     fNonHFEEff->SetTitle("IPEfficiencyForNonHFE");
544     fNonHFEEff->SetName("IPEfficiencyForNonHFE");
545     fNonHFEEff->Write();
546     
547    
548     out->Close(); 
549     delete out;
550    }
551   return kTRUE;
552 }
553 //____________________________________________________________
554 AliCFDataGrid* AliHFEBeautySpectrum::SubtractBackground(Bool_t setBackground){
555   //
556   // Apply background subtraction for IP analysis
557   //
558   
559   Int_t nbins=1;
560   printf("subtracting background\n");
561   // Raw spectrum
562   AliCFContainer *dataContainer = GetContainer(kDataContainer);
563   if(!dataContainer){
564     AliError("Data Container not available");
565     return NULL;
566   }
567   printf("Step data: %d\n",fStepData);
568   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
569   
570   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
571   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); 
572   
573   // Background Estimate
574   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
575   if(!backgroundContainer){
576     AliError("MC background container not found");
577     return NULL;
578   }
579     
580   Int_t stepbackground = 1;
581   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
582   
583   const Int_t bgPlots = 8;
584   TH1D *incElec = 0x0;
585   TH1D *hadron = 0x0;
586   TH1D *charm = 0x0;
587   TH1D *nonHFE[bgPlots];
588   TH1D *subtracted = 0x0;
589   
590   THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();  
591   incElec = (TH1D *) sparseIncElec->Projection(0);
592   AliHFEtools::NormaliseBinWidth(incElec);
593
594   TH1D* htemp;
595   Int_t* bins=new Int_t[2];
596
597   if(fIPanaHadronBgSubtract){
598     // Hadron background
599     printf("Hadron background for IP analysis subtracted!\n");
600     htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
601     bins[0]=htemp->GetNbinsX();
602     
603     AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
604     hbgContainer->SetGrid(fHadronEffbyIPcut);
605     backgroundGrid->Multiply(hbgContainer,1);
606     // subtract hadron contamination
607     spectrumSubtracted->Add(backgroundGrid,-1.0);
608     // draw raw hadron bg spectra     
609     THnSparseF* sparseHbg = (THnSparseF *) hbgContainer->GetGrid();
610     hadron = (TH1D *) sparseHbg->Projection(0);
611     AliHFEtools::NormaliseBinWidth(hadron);
612   }
613
614   if(fIPanaCharmBgSubtract){
615     // Charm background
616     printf("Charm background for IP analysis subtracted!\n");
617     AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
618     spectrumSubtracted->Add(charmbgContainer,-1.0);
619     THnSparseF *sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();  
620     charm = (TH1D *) sparseCharmElec->Projection(0);
621     AliHFEtools::NormaliseBinWidth(charm); 
622   }
623
624   const Char_t *sourceName[bgPlots]={"Conversion","Dalitz","K0s secondary pion Dalitz","K0s secondary pion conversions","Other secondary pion Dalitz","Other secondary pion Dalitz conversions","Other Dalitz", "Other conversions"};
625   const Int_t style[2] = {21,22};
626   const Int_t color[3] = {7,12,8};
627   if(fIPanaNonHFEBgSubtract){ 
628     Int_t iPlot = 0;
629     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
630       if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
631       for(Int_t decay = 0; decay < 2; decay++){
632         AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground(decay,iSource);//conversions/Dalitz decays from different sources
633         if(iSource == 0)
634           spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
635         THnSparseF* sparseNonHFEelecs = (THnSparseF *) nonHFEbgContainer->GetGrid();
636         nonHFE[iPlot] = (TH1D *) sparseNonHFEelecs->Projection(0);
637         AliHFEtools::NormaliseBinWidth(nonHFE[iPlot]);
638         iPlot++;
639       }
640       if(!(fNonHFEsyst && (fBeamType == 1)))break;
641     } 
642   }
643
644   TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
645   TCanvas *cRaw = new TCanvas("cRaw","cRaw",1000,800);     
646   cRaw->cd();
647   gPad->SetLogx();
648   gPad->SetLogy();
649   incElec->GetXaxis()->SetRangeUser(0.4,8.);
650   incElec->Draw("p");
651   incElec->SetMarkerColor(1);
652   incElec->SetMarkerStyle(20);
653   lRaw->AddEntry(incElec,"inclusive electron spectrum");
654   
655   if(fIPanaCharmBgSubtract){
656     charm->Draw("samep");
657     charm->SetMarkerColor(3);
658     charm->SetMarkerStyle(20);
659     charm->GetXaxis()->SetRangeUser(0.4,8.);
660     lRaw->AddEntry(charm,"charm elecs");
661   }
662   
663   if(fIPanaNonHFEBgSubtract){
664     Int_t iPlot = 0;
665     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
666       if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
667       for(Int_t decay = 0; decay < 2; decay++){
668         nonHFE[iPlot]->Draw("samep");
669         if(iSource == 0){
670           nonHFE[iPlot]->SetMarkerStyle(20);
671           if(decay == 0)
672             nonHFE[iPlot]->SetMarkerColor(4);
673           else
674             nonHFE[iPlot]->SetMarkerColor(6);
675         }
676         else{     
677           nonHFE[iPlot]->SetMarkerColor(color[iSource-6]);      
678           nonHFE[iPlot]->SetMarkerStyle(style[decay]);
679         }
680         nonHFE[iPlot]->GetXaxis()->SetRangeUser(0.4,8.);
681         lRaw->AddEntry(nonHFE[iPlot],sourceName[iPlot]);
682         iPlot++;
683       }
684       if(!(fNonHFEsyst && (fBeamType == 1)))break;//only draw standard sources
685     }
686   }
687
688   THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
689   subtracted = (TH1D *) sparseSubtracted->Projection(0);
690   AliHFEtools::NormaliseBinWidth(subtracted);
691   subtracted->Draw("samep");
692   subtracted->SetMarkerStyle(24);
693   lRaw->AddEntry(subtracted,"subtracted electron spectrum");
694   lRaw->Draw("SAME");
695  
696   delete[] bins; 
697
698   TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(0);
699   AliHFEtools::NormaliseBinWidth(measuredTH1background);
700
701   if(setBackground){
702     if(fBackground) delete fBackground;
703     fBackground = backgroundGrid;
704   } else delete backgroundGrid;
705   
706   fQA->AddResultAt(subtracted,AliHFEBeautySpectrumQA::kAfterSC);
707   fQA->AddResultAt(incElec,AliHFEBeautySpectrumQA::kBeforeSC);
708   fQA->AddResultAt(measuredTH1background, AliHFEBeautySpectrumQA::kMeasBG);
709   fQA->DrawSubtractContamination();
710   
711   return spectrumSubtracted;
712 }
713 //____________________________________________________________________________________________
714 AliCFDataGrid* AliHFEBeautySpectrum::GetNonHFEBackground(Int_t decay, Int_t source){
715   //
716   // calculate non-HF electron background
717   // Arguments: decay = 0 for conversions, decay = 1 for meson decays to electrons
718   // source: gives mother either of the electron (for meson->electron) or of the gamma (conversions);
719   // source numbers as given by array in AliAnalysisTaskHFE: 
720   // {"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
721   //
722   
723   Double_t evtnorm[1] = {1.};
724   if(fNMCbgEvents) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
725   
726   Int_t stepbackground = 1;//for non-standard background, step after IP cut is taken directly
727   if(source == 0) stepbackground = 3;//for standard sources, take step before PID -> correction of PID x IP efficiencies from enhanced samples
728   AliCFContainer *backgroundContainer = 0x0;
729   if(fNonHFEsyst){
730     if(decay == 0){//conversions
731       backgroundContainer = (AliCFContainer*)fConvSourceContainer[source][0]->Clone();
732       if(source == 0)//standard conversion containers
733         for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){//add all standard sources
734           backgroundContainer->Add(fConvSourceContainer[iSource][0]);     
735         }
736     }
737     else if(decay == 1){//Dalitz decays, same procedure as above
738       backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[source][0]->Clone();
739       if(source == 0)
740         for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
741           if(iSource == 1)
742             backgroundContainer->Add(fNonHFESourceContainer[iSource][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
743           else
744             backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);   
745         }
746     }   
747   }
748   else{//careful: these containers include the non-standard electron sources (K0s, etc.). If necessary, use SetRange in species axis to fix it?
749     if(source == 0){
750       // Background Estimate
751       if(decay == 0)
752         backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
753       else if(decay == 1)
754         backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
755     } 
756   }
757   if(!backgroundContainer){
758     AliError("MC background container not found");
759     return NULL;
760   }
761   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("backgroundGrid","backgroundGrid",*backgroundContainer,stepbackground);
762  
763   Double_t rangelow = 1.;
764   Double_t rangeup = 6.;
765   if(decay == 1) rangelow = 0.9;
766
767
768   const Char_t *dmode[2]={"Conversions","Dalitz decays"};
769   TF1 *fitHagedorn = new TF1("fitHagedorn", "[0]/TMath::Power(([1]+x/[2]), [3])", rangelow, rangeup);
770   fNonHFEbg = 0x0;
771   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
772   TAxis *axis = h1->GetXaxis();
773   AliHFEtools::NormaliseBinWidth(h1);
774   if(source == 0){
775     fitHagedorn->SetParameter(0, 0.15);
776     fitHagedorn->SetParameter(1, 0.09);
777     fitHagedorn->SetParameter(2, 8.4);
778     fitHagedorn->SetParameter(3, 6.3);
779     // TCanvas *ch1conversion = new TCanvas("ch1conversion","ch1conversion",500,400);
780     // ch1conversion->cd();
781     //fHagedorn->SetLineColor(2);
782     h1->Fit(fitHagedorn,"R");
783     // h1->Draw();
784     if(!(fNonHFEbg = h1->GetFunction("fitHagedorn")))printf("electron background fit failed for %s\n",dmode[decay]);
785   }
786   
787   Int_t *nBinpp=new Int_t[1];
788   Int_t *binspp=new Int_t[1];
789   binspp[0]=kSignalPtBins;// number of pt bins
790   
791   Int_t looppt=binspp[0];
792   
793   for(Long_t iBin=1; iBin<= looppt;iBin++){       
794     Double_t iipt= h1->GetBinCenter(iBin);
795     Double_t iiwidth= axis->GetBinWidth(iBin);
796     nBinpp[0]=iBin;
797     Double_t fitFactor = backgroundGrid->GetElement(nBinpp);//if no fit available, just take bin-by-bin information
798     if(fNonHFEbg)fitFactor = fNonHFEbg->Eval(iipt)*iiwidth;
799     backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
800     backgroundGrid->SetElement(nBinpp,fitFactor*evtnorm[0]);
801   }
802   //end of workaround for statistical errors
803   if(source == 0){  
804     AliCFDataGrid *weightedBGContainer = 0x0;   
805     weightedBGContainer = new AliCFDataGrid("weightedBGContainer","weightedBGContainer",1,binspp);
806     if(decay == 0)
807       weightedBGContainer->SetGrid(GetPIDxIPEff(2));
808     else if(decay == 1)
809       weightedBGContainer->SetGrid(GetPIDxIPEff(3));
810     if(stepbackground == 3){    
811       backgroundGrid->Multiply(weightedBGContainer,1.0);
812     }  
813   }
814   delete[] nBinpp;
815   delete[] binspp;  
816   return backgroundGrid;
817 }
818 //____________________________________________________________
819 AliCFDataGrid* AliHFEBeautySpectrum::GetCharmBackground(){
820   //
821   // calculate charm background; when using centrality binning 2: centBin = 0 for 0-20%, centBin = 5 for 40-80%
822   //
823
824   Int_t nDim = 1;
825  
826   Double_t evtnorm=0;
827   if(fNMCEvents) evtnorm= double(fNEvents)/double(fNMCEvents); //mjmjmj
828   //if(fNMCbgEvents) evtnorm= double(fNEvents)/double(fNMCbgEvents);
829   printf("events for data: %d",fNEvents);
830   printf("events for MC: %d",fNMCEvents);
831   printf("normalization factor for charm: %f",evtnorm);
832   
833   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
834   if(!mcContainer){
835     AliError("MC Container not available");
836     return NULL;
837   }
838
839   if(!fCorrelation){
840     AliError("No Correlation map available");
841     return NULL;
842   }
843  
844   AliCFDataGrid *charmBackgroundGrid= 0x0;
845   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
846   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
847   Int_t* bins=new Int_t[2];
848   bins[0]=charmbgaftertofpid->GetNbinsX();
849    
850   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
851   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
852   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
853   
854   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
855   Int_t *nBinpp=new Int_t[1];
856   Int_t* binspp=new Int_t[1];
857   binspp[0]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
858   
859   Int_t looppt=binspp[0];
860   
861   for(Long_t iBin=1; iBin<= looppt;iBin++){
862     nBinpp[0]=iBin;
863     charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
864     charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
865   }
866   
867   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
868   
869   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
870   weightedCharmContainer->SetGrid(GetCharmWeights(fTestCentralityLow)); // get charm weighting factors
871   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
872   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
873   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
874   
875   // Efficiency (set efficiency to 1 for only folding) 
876   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
877   efficiencyD->CalculateEfficiency(0,0);
878   
879   // Folding
880   if(fBeamType==0)nDim = 1;
881   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
882   folding.SetMaxNumberOfIterations(1);
883   folding.Unfold();
884   
885   // Results
886   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
887   THnSparse* result=(THnSparse*)result1->Clone();
888   charmBackgroundGrid->SetGrid(result);
889   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
890
891   //Charm background evaluation plots
892   
893   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
894   cCharmBgEval->Divide(3,1);
895   
896   cCharmBgEval->cd(1);
897   
898   if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
899     
900   AliHFEtools::NormaliseBinWidth(charmbgaftertofpid);
901   charmbgaftertofpid->SetMarkerStyle(25);
902   charmbgaftertofpid->Draw("p");
903   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
904   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
905   gPad->SetLogy();
906   
907   AliHFEtools::NormaliseBinWidth(charmbgafteripcut);
908   charmbgafteripcut->SetMarkerStyle(24);
909   charmbgafteripcut->Draw("samep");
910   
911   AliHFEtools::NormaliseBinWidth(charmbgafterweight);
912   charmbgafterweight->SetMarkerStyle(24);
913   charmbgafterweight->SetMarkerColor(4);
914   charmbgafterweight->Draw("samep");
915   
916   AliHFEtools::NormaliseBinWidth(charmbgafterfolding);
917   charmbgafterfolding->SetMarkerStyle(24);
918   charmbgafterfolding->SetMarkerColor(2);
919   charmbgafterfolding->Draw("samep");
920   
921   cCharmBgEval->cd(2);
922   parametrizedcharmpidipeff->SetMarkerStyle(24);
923   parametrizedcharmpidipeff->Draw("p");
924   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
925
926   cCharmBgEval->cd(3);
927   charmweightingfc->SetMarkerStyle(24);
928   charmweightingfc->Draw("p");
929   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
930   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
931
932   cCharmBgEval->cd(1);
933   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
934   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
935   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
936   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
937   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
938   legcharmbg->Draw("same");
939
940   cCharmBgEval->cd(2);
941   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
942   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
943   legcharmbg2->Draw("same");
944
945   CorrectStatErr(charmBackgroundGrid);
946   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
947
948   delete[] bins;
949   delete[] nBinpp;
950   delete[] binspp;
951   
952   if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
953
954   return charmBackgroundGrid;
955 }
956 //____________________________________________________________
957 AliCFDataGrid *AliHFEBeautySpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){  
958   //
959   // Apply TPC pid efficiency correction from parametrisation
960   //
961
962   // Data in the right format
963   AliCFDataGrid *dataGrid = 0x0;  
964   if(bgsubpectrum) {
965     dataGrid = bgsubpectrum;
966   }
967   else {
968     
969     AliCFContainer *dataContainer = GetContainer(kDataContainer);
970     if(!dataContainer){
971       AliError("Data Container not available");
972       return NULL;
973     }
974     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
975   } 
976   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
977   result->SetName("ParametrizedEfficiencyBefore");
978   THnSparse *h = result->GetGrid();
979   Int_t nbdimensions = h->GetNdimensions();
980   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
981   AliCFContainer *dataContainer = GetContainer(kDataContainer);
982   if(!dataContainer){
983     AliError("Data Container not available");
984     return NULL;
985   }
986   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
987   dataContainerbis->Add(dataContainerbis,-1.0);
988
989   Int_t* coord = new Int_t[nbdimensions];
990   memset(coord, 0, sizeof(Int_t) * nbdimensions);
991   Double_t* points = new Double_t[nbdimensions];
992
993   ULong64_t nEntries = h->GetNbins();
994   for (ULong64_t i = 0; i < nEntries; ++i) {   
995     Double_t value = h->GetBinContent(i, coord);
996     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
997     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
998     
999     // Get the bin co-ordinates given an coord
1000     for (Int_t j = 0; j < nbdimensions; ++j)
1001       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1002
1003     if (!fEfficiencyFunction->IsInside(points))
1004          continue;
1005     TF1::RejectPoint(kFALSE);
1006
1007     // Evaulate function at points
1008     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1009     //printf("Value efficiency is %f\n",valueEfficiency);
1010
1011     if(valueEfficiency > 0.0) {
1012       h->SetBinContent(coord,value/valueEfficiency);
1013       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1014     }
1015     Double_t error = h->GetBinError(i);
1016     h->SetBinError(coord,error/valueEfficiency);
1017     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);  
1018   } 
1019
1020   delete[] coord;
1021   delete[] points;
1022
1023   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1024
1025 // QA
1026   TH1D *afterE = (TH1D *) resultt->Project(0);
1027   AliHFEtools::NormaliseBinWidth(afterE);
1028   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1029   AliHFEtools::NormaliseBinWidth(beforeE);
1030   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterPE);
1031   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforePE);
1032   fQA->AddResultAt(fEfficiencyFunction,AliHFEBeautySpectrumQA::kPEfficiency);
1033   fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kParametrized);
1034   
1035   return resultt;
1036 }
1037 //____________________________________________________________
1038 THnSparse *AliHFEBeautySpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1039   
1040   //
1041   // Unfold and eventually correct for efficiency the bgsubspectrum
1042   //
1043
1044   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1045   if(!mcContainer){
1046     AliError("MC Container not available");
1047     return NULL;
1048   }
1049
1050   if(!fCorrelation){
1051     AliError("No Correlation map available");
1052     return NULL;
1053   }
1054
1055   // Data 
1056   AliCFDataGrid *dataGrid = 0x0;  
1057   if(bgsubpectrum) {
1058     dataGrid = bgsubpectrum;
1059   }
1060   else {
1061
1062     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1063     if(!dataContainer){
1064       AliError("Data Container not available");
1065       return NULL;
1066     }
1067
1068     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1069   } 
1070   
1071   // Guessed
1072   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1073   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1074
1075   // Efficiency
1076   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1077   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1078
1079   if(!fBeauty2ndMethod)
1080   {
1081     // Consider parameterized IP cut efficiency
1082     Int_t* bins=new Int_t[1];
1083     bins[0]=kSignalPtBins;
1084     
1085     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1086     beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
1087     efficiencyD->Multiply(beffContainer,1);
1088   }
1089   
1090
1091   // Unfold 
1092   
1093   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);//check with MinJung if last arguments are correct (taken from inclusive analysis...
1094   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1095   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1096   if(fSetSmoothing) unfolding.UseSmoothing();
1097   unfolding.Unfold();
1098
1099   // Results
1100   THnSparse* result = unfolding.GetUnfolded();
1101   THnSparse* residual = unfolding.GetEstMeasured();
1102
1103  // QA
1104   TH1D *residualh = (TH1D *) residual->Projection(0);
1105   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1106   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1107   TH1D *afterE = (TH1D *) result->Projection(0);
1108
1109   AliHFEtools::NormaliseBinWidth(residualh);
1110   AliHFEtools::NormaliseBinWidth(beforeE);
1111   AliHFEtools::NormaliseBinWidth(afterE);
1112   fQA->AddResultAt(residualh,AliHFEBeautySpectrumQA::kResidualU);
1113   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterU);
1114   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeU);
1115   fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kUEfficiency);
1116   fQA->DrawUnfolding();
1117
1118   return (THnSparse *) result->Clone();
1119 }
1120
1121 //____________________________________________________________
1122 AliCFDataGrid *AliHFEBeautySpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1123   
1124   //
1125   // Apply unfolding and efficiency correction together to bgsubspectrum
1126   //
1127
1128   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1129   if(!mcContainer){
1130     AliError("MC Container not available");
1131     return NULL;
1132   }
1133
1134   // Efficiency
1135   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1136   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1137
1138
1139   if(!fBeauty2ndMethod)
1140   {
1141     // Consider parameterized IP cut efficiency
1142     Int_t* bins=new Int_t[1];
1143     bins[0]=kSignalPtBins;
1144     
1145     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1146     beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
1147     efficiencyD->Multiply(beffContainer,1);
1148   }
1149
1150   // Data in the right format
1151   AliCFDataGrid *dataGrid = 0x0;  
1152   if(bgsubpectrum) {
1153     dataGrid = bgsubpectrum;
1154   }
1155   else {
1156     
1157     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1158     if(!dataContainer){
1159       AliError("Data Container not available");
1160       return NULL;
1161     }
1162
1163     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1164   } 
1165
1166   // Correct
1167   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1168   result->ApplyEffCorrection(*efficiencyD);
1169
1170   // QA
1171   TH1D *afterE = (TH1D *) result->Project(0);
1172   AliHFEtools::NormaliseBinWidth(afterE);
1173   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1174   AliHFEtools::NormaliseBinWidth(beforeE);
1175   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1176   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterMCE);
1177   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeMCE);
1178   fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kMCEfficiency);
1179   fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kMC);
1180   
1181   return result;
1182
1183 }
1184 //____________________________________________________________
1185 void AliHFEBeautySpectrum::AddTemporaryObject(TObject *o){
1186   // 
1187   // Emulate garbage collection on class level
1188   // 
1189   if(!fTemporaryObjects) {
1190     fTemporaryObjects= new TList;
1191     fTemporaryObjects->SetOwner();
1192   }
1193   fTemporaryObjects->Add(o);
1194 }
1195
1196 //____________________________________________________________
1197 void AliHFEBeautySpectrum::ClearObject(TObject *o){
1198   //
1199   // Do a safe deletion
1200   //
1201   if(fTemporaryObjects){
1202     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
1203     delete o;
1204   } else{
1205     // just delete
1206     delete o;
1207   }
1208 }
1209 //_________________________________________________________________________
1210 THnSparse* AliHFEBeautySpectrum::GetCharmWeights(Int_t centBin){
1211  
1212   //
1213   // Measured D->e based weighting factors
1214   //
1215
1216   const Int_t nDimpp=1;
1217   Int_t nBinpp[nDimpp] = {kSignalPtBins};
1218   Double_t ptbinning1[kSignalPtBins+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
1219   Int_t looppt=nBinpp[0];
1220   
1221   fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
1222   fWeightCharm->SetBinEdges(0, ptbinning1);
1223  
1224   // //if(fBeamType == 0){// Weighting factor for pp
1225   //Double_t weightpp[kSignalPtBins]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
1226 //TPC+TOF standard cut 4800
1227   Double_t weightpp[kSignalPtBins]={0.050326, 0.045826, 0.042043, 0.039641, 0.039872, 0.041796, 0.044320, 0.046273, 0.047376, 0.047657, 0.047973, 0.048307, 0.045325, 0.044067, 0.043367, 0.040417, 0.037048, 0.033695, 0.032192, 0.029270, 0.027270, 0.024451, 0.020846, 0.019032, 0.018210, 0.017554, 0.015604, 0.015194, 0.013542, 0.013447, 0.015160, 0.015507, 0.014989, 0.012533, 0.025603}; 
1228   //}
1229   //else{
1230   //if(centBin == 0){
1231       // Weighting factor for PbPb (0-20%)
1232   Double_t weightPbPb1[kSignalPtBins]={0.641897,  0.640472,  0.615228,  0.650469,  0.737762,  0.847867,  1.009317,  1.158594,  1.307482,  1.476973,  1.551131,  1.677131,  1.785478,  1.888933,  2.017957,  2.074757,  1.926700,  1.869495,  1.546558,  1.222873,  1.160313,  0.903375,  0.799642,  0.706244,  0.705449,  0.599947,  0.719570,  0.499422,  0.703978,  0.477452,  0.325057,  0.093391,  0.096675,  0.000000,  0.000000};
1233   //}
1234   //else{
1235   // Weighting factor for PbPb (40-80%)
1236   Double_t weightPbPb2[kSignalPtBins]={0.181953,  0.173248,  0.166799,  0.182558,  0.206581,  0.236955,  0.279390,  0.329129,  0.365260,  0.423059,  0.452057,  0.482726,  0.462627,  0.537770,  0.584663,  0.579452,  0.587194,  0.499498,  0.443299,  0.398596,  0.376695,  0.322331,  0.260890,  0.374834,  0.249114,  0.310330,  0.287326,  0.243174,  0.758945,  0.138867,  0.170576,  0.107797,  0.011390,  0.000000,  0.000000};
1237   //  }
1238   //}
1239   Double_t weight[kSignalPtBins];
1240   for(Int_t ipt = 0; ipt < kSignalPtBins; ipt++){
1241     if(fBeamType == 0)weight[ipt] = weightpp[ipt];
1242     else if(centBin == 0)weight[ipt] = weightPbPb1[ipt];
1243     else weight[ipt] = weightPbPb2[ipt];
1244   }  
1245
1246   //points
1247   Double_t pt[1];
1248   Double_t contents[2];
1249   
1250   for(int i=0; i<looppt; i++){
1251     pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
1252     contents[0]=pt[0];
1253     fWeightCharm->Fill(contents,weight[i]);
1254   }
1255   
1256   
1257   Int_t nDimSparse = fWeightCharm->GetNdimensions();
1258   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1259   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1260
1261   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1262       binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
1263       nBins *= binsvar[iVar];
1264   }
1265
1266   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1267   // loop that sets 0 error in each bin
1268   for (Long_t iBin=0; iBin<nBins; iBin++) {
1269     Long_t bintmp = iBin ;
1270     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1271       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1272       bintmp /= binsvar[iVar] ;
1273     }
1274     fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
1275   }
1276   delete[] binsvar;
1277   delete[] binfill;
1278
1279   return fWeightCharm;
1280 }
1281 //____________________________________________________________________________
1282 void AliHFEBeautySpectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
1283
1284    // TOF PID efficiencies
1285   
1286    TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1287    TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,6.);
1288    TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.0);
1289
1290    if(fBeamType == 1){//not in use yet - adapt as soon as possible!
1291      fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1292      fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1293      fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1294    }
1295
1296    TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
1297    cefficiencyParamtof->cd();
1298
1299    AliCFContainer *mccontainermcD = 0x0;
1300    AliCFContainer *mccontaineresdD = 0x0;
1301    TH1D* efficiencysigTOFPIDD;
1302    TH1D* efficiencyTOFPIDD;
1303    TH1D* efficiencysigesdTOFPIDD;
1304    TH1D* efficiencyesdTOFPIDD;
1305    Int_t source = -1; //get parameterized TOF PID efficiencies
1306
1307    // signal sample
1308    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1309    AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
1310    efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1311    
1312    // mb sample for double check
1313    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1314    AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
1315    efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1316    
1317    // mb sample with reconstructed variables
1318    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1319    AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
1320    efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1321    
1322    // mb sample with reconstructed variables
1323    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1324    AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
1325    efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1326    
1327    //fill histo
1328    efficiencysigTOFPIDD = (TH1D *) efficiencymcsigParamTOFPID->Project(0);
1329    efficiencyTOFPIDD = (TH1D *) efficiencymcParamTOFPID->Project(0);
1330    efficiencysigesdTOFPIDD = (TH1D *) efficiencysigesdParamTOFPID->Project(0);
1331    efficiencyesdTOFPIDD = (TH1D *) efficiencyesdParamTOFPID->Project(0);
1332    efficiencysigTOFPIDD->SetName("efficiencysigTOFPIDD");
1333    efficiencyTOFPIDD->SetName("efficiencyTOFPIDD");
1334    efficiencysigesdTOFPIDD->SetName("efficiencysigesdTOFPIDD");
1335    efficiencyesdTOFPIDD->SetName("efficiencyesdTOFPIDD");
1336    
1337    //fit (mc enhenced sample)
1338    fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1339    efficiencysigTOFPIDD->Fit(fittofpid,"R");
1340    efficiencysigTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1341    fEfficiencyTOFPIDD = efficiencysigTOFPIDD->GetFunction("fittofpid");
1342    
1343    //fit (esd enhenced sample)
1344    efficiencysigesdTOFPIDD->Fit(fittofpid,"R");
1345    efficiencysigesdTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1346    fEfficiencyesdTOFPIDD = efficiencysigesdTOFPIDD->GetFunction("fittofpid");
1347    
1348    // draw (for PbPb, only 1st bin)
1349    //sig mc
1350    efficiencysigTOFPIDD->SetTitle("");
1351    efficiencysigTOFPIDD->SetStats(0);
1352    efficiencysigTOFPIDD->SetMarkerStyle(25);
1353    efficiencysigTOFPIDD->SetMarkerColor(2);
1354    efficiencysigTOFPIDD->SetLineColor(2);
1355    efficiencysigTOFPIDD->Draw();
1356
1357    //mb mc
1358    efficiencyTOFPIDD->SetTitle("");
1359    efficiencyTOFPIDD->SetStats(0);
1360    efficiencyTOFPIDD->SetMarkerStyle(24);
1361    efficiencyTOFPIDD->SetMarkerColor(4);
1362    efficiencyTOFPIDD->SetLineColor(4);
1363    efficiencyTOFPIDD->Draw("same");
1364
1365    //sig esd
1366    efficiencysigesdTOFPIDD->SetTitle("");
1367    efficiencysigesdTOFPIDD->SetStats(0);
1368    efficiencysigesdTOFPIDD->SetMarkerStyle(25);
1369    efficiencysigesdTOFPIDD->SetMarkerColor(3);
1370    efficiencysigesdTOFPIDD->SetLineColor(3);
1371    efficiencysigesdTOFPIDD->Draw("same");
1372
1373    //mb esd
1374    efficiencyesdTOFPIDD->SetTitle("");
1375    efficiencyesdTOFPIDD->SetStats(0);
1376    efficiencyesdTOFPIDD->SetMarkerStyle(25);
1377    efficiencyesdTOFPIDD->SetMarkerColor(1);
1378    efficiencyesdTOFPIDD->SetLineColor(1);
1379    efficiencyesdTOFPIDD->Draw("same");
1380
1381    //signal mc fit
1382    if(fEfficiencyTOFPIDD){
1383      fEfficiencyTOFPIDD->SetLineColor(2);
1384      fEfficiencyTOFPIDD->Draw("same");
1385    }
1386    //mb esd fit
1387    if(fEfficiencyesdTOFPIDD){
1388        fEfficiencyesdTOFPIDD->SetLineColor(3);
1389        fEfficiencyesdTOFPIDD->Draw("same");
1390      }
1391
1392    TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
1393    legtofeff->AddEntry(efficiencysigTOFPIDD,"TOF PID Step Efficiency","");
1394    legtofeff->AddEntry(efficiencysigTOFPIDD,"vs MC p_{t} for enhenced samples","p");
1395    legtofeff->AddEntry(efficiencyTOFPIDD,"vs MC p_{t} for mb samples","p");
1396    legtofeff->AddEntry(efficiencysigesdTOFPIDD,"vs esd p_{t} for enhenced samples","p");
1397    legtofeff->AddEntry(efficiencyesdTOFPIDD,"vs esd p_{t} for mb samples","p");
1398    legtofeff->Draw("same");
1399
1400
1401    TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
1402    cefficiencyParamIP->cd();
1403    gStyle->SetOptStat(0);
1404
1405    // IP cut efficiencies
1406    AliCFContainer *charmCombined = 0x0; 
1407    AliCFContainer *beautyCombined = 0x0;
1408    AliCFContainer *beautyCombinedesd = 0x0;
1409
1410    source = 0; //charm enhenced
1411    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1412    AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
1413    efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
1414
1415    charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");  
1416
1417    source = 1; //beauty enhenced
1418    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1419    AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
1420    efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
1421
1422    beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined"); 
1423
1424    mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1425    AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
1426    efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1427
1428    beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
1429
1430    source = 0; //charm mb
1431    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1432    AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
1433    efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
1434    
1435    charmCombined->Add(mccontainermcD); 
1436    AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined); 
1437    efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
1438
1439    source = 1; //beauty mb
1440    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1441    AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
1442    efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
1443    
1444    beautyCombined->Add(mccontainermcD);
1445    AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined); 
1446    efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
1447    
1448    mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1449    AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
1450    efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1451
1452    beautyCombinedesd->Add(mccontaineresdD);
1453    AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
1454    efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1455    
1456    source = 2; //conversion mb
1457    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1458    AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
1459    efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
1460
1461    source = 3; //non HFE except for the conversion mb
1462    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1463    AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
1464    efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1465
1466    if(fIPEffCombinedSamples){printf("combined samples taken for beauty and charm\n");
1467      fEfficiencyCharmSigD = (TH1D*)efficiencyCharmCombined->Project(0); //signal enhenced + mb 
1468      fEfficiencyBeautySigD = (TH1D*)efficiencyBeautyCombined->Project(0); //signal enhenced + mb
1469      fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautyCombinedesd->Project(0); //signal enhenced + mb
1470      }
1471      else{printf("signal enhanced samples taken for beauty and charm\n");
1472        fEfficiencyCharmSigD = (TH1D*)efficiencyCharmSig->Project(0); //signal enhenced only
1473        fEfficiencyBeautySigD = (TH1D*)efficiencyBeautySig->Project(0); //signal enhenced only
1474        fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautySigesd->Project(0); //signal enhenced only
1475      }
1476      fCharmEff = (TH1D*)efficiencyCharm->Project(0); //mb only
1477      fBeautyEff = (TH1D*)efficiencyBeauty->Project(0); //mb only
1478      fConversionEff = (TH1D*)efficiencyConv->Project(0); //mb only
1479      fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); //mb only
1480
1481    if(fBeamType==0){
1482      //AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
1483      AliCFContainer *cfcontainer = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
1484      if(!cfcontainer) return;
1485      AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(cfcontainer,1,0); //mjmj
1486      fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
1487      
1488      //AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
1489      AliCFContainer *cfcontainerr = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
1490      if(!cfcontainerr) return;
1491      AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(cfcontainerr,1,0); //mjmj
1492      fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
1493      
1494      //MHMH
1495      if(fBeamType == 0){
1496        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1497        fNonHFEEffbgc->Fit(fipfitnonhfe,"R"); 
1498        fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1499        
1500        fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1501        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);      
1502        fConversionEffbgc->Fit(fipfitnonhfe,"R");
1503        fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1504      }
1505      //MHMH
1506    }
1507    
1508    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1509    fipfit->SetLineColor(2);
1510    fEfficiencyBeautySigD->Fit(fipfit,"R");
1511    fEfficiencyBeautySigD->GetYaxis()->SetTitle("Efficiency");
1512    fEfficiencyIPBeautyD = fEfficiencyBeautySigD->GetFunction("fipfit");
1513    
1514    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1515    fipfit->SetLineColor(6);
1516    fEfficiencyBeautySigesdD->Fit(fipfit,"R");
1517    fEfficiencyBeautySigesdD->GetYaxis()->SetTitle("Efficiency");
1518    fEfficiencyIPBeautyesdD = fEfficiencyBeautySigesdD->GetFunction("fipfit");
1519    
1520    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1521    fipfit->SetLineColor(1);
1522    fEfficiencyCharmSigD->Fit(fipfit,"R");
1523    fEfficiencyCharmSigD->GetYaxis()->SetTitle("Efficiency");
1524    fEfficiencyIPCharmD = fEfficiencyCharmSigD->GetFunction("fipfit");
1525    
1526    //if(0){
1527    if(fIPParameterizedEff){
1528      // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1529      fipfitnonhfe->SetLineColor(3);
1530      if(fBeamType==1){
1531        fConversionEff->Fit(fipfitnonhfe,"R");
1532        fConversionEff->GetYaxis()->SetTitle("Efficiency");
1533        fEfficiencyIPConversionD = fConversionEff->GetFunction("fipfitnonhfe");
1534      }
1535      else{
1536        fConversionEffbgc->Fit(fipfitnonhfe,"R");
1537        fConversionEffbgc->GetYaxis()->SetTitle("Efficiency");
1538        fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1539      }       
1540      // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1541      fipfitnonhfe->SetLineColor(4);
1542      if(fBeamType==1){
1543        fNonHFEEff->Fit(fipfitnonhfe,"R");
1544        fNonHFEEff->GetYaxis()->SetTitle("Efficiency");
1545        fEfficiencyIPNonhfeD = fNonHFEEff->GetFunction("fipfitnonhfe");
1546      }
1547      else{
1548        fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1549        fNonHFEEffbgc->GetYaxis()->SetTitle("Efficiency");
1550        fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1551      }
1552    }
1553    
1554    // draw
1555    fEfficiencyCharmSigD->SetMarkerStyle(21);
1556    fEfficiencyCharmSigD->SetMarkerColor(1);
1557    fEfficiencyCharmSigD->SetLineColor(1);
1558    fEfficiencyBeautySigD->SetMarkerStyle(21);
1559    fEfficiencyBeautySigD->SetMarkerColor(2);
1560    fEfficiencyBeautySigD->SetLineColor(2);
1561    fEfficiencyBeautySigesdD->SetStats(0);
1562    fEfficiencyBeautySigesdD->SetMarkerStyle(21);
1563    fEfficiencyBeautySigesdD->SetMarkerColor(6);
1564    fEfficiencyBeautySigesdD->SetLineColor(6);
1565    fCharmEff->SetMarkerStyle(24);
1566    fCharmEff->SetMarkerColor(1);
1567    fCharmEff->SetLineColor(1);
1568    fBeautyEff->SetMarkerStyle(24);
1569    fBeautyEff->SetMarkerColor(2);
1570    fBeautyEff->SetLineColor(2);
1571    fConversionEff->SetMarkerStyle(24);
1572    fConversionEff->SetMarkerColor(3);
1573    fConversionEff->SetLineColor(3);
1574    fNonHFEEff->SetMarkerStyle(24);
1575    fNonHFEEff->SetMarkerColor(4);
1576    fNonHFEEff->SetLineColor(4);
1577
1578    fEfficiencyCharmSigD->Draw();
1579    fEfficiencyCharmSigD->GetXaxis()->SetRangeUser(0.0,7.9);
1580    fEfficiencyCharmSigD->GetYaxis()->SetRangeUser(0.0,0.5);
1581
1582    fEfficiencyBeautySigD->Draw("same");
1583    fEfficiencyBeautySigesdD->Draw("same");
1584    if(fBeamType == 1){
1585      fNonHFEEff->Draw("same");
1586      fConversionEff->Draw("same");
1587    }
1588    //fCharmEff->Draw("same");
1589    //fBeautyEff->Draw("same");
1590
1591    if(fBeamType==0){
1592      fConversionEffbgc->SetMarkerStyle(25);
1593      fConversionEffbgc->SetMarkerColor(3);
1594      fConversionEffbgc->SetLineColor(3);
1595      fNonHFEEffbgc->SetMarkerStyle(25);
1596      fNonHFEEffbgc->SetMarkerColor(4);
1597      fNonHFEEffbgc->SetLineColor(4);
1598      fConversionEffbgc->Draw("same");
1599      fNonHFEEffbgc->Draw("same");
1600    }
1601  
1602    if(fEfficiencyIPBeautyD)
1603       fEfficiencyIPBeautyD->Draw("same");
1604    if(fEfficiencyIPBeautyesdD)
1605      fEfficiencyIPBeautyesdD->Draw("same");
1606    if( fEfficiencyIPCharmD)
1607      fEfficiencyIPCharmD->Draw("same");
1608    if(fIPParameterizedEff){
1609      if(fEfficiencyIPConversionD)
1610        fEfficiencyIPConversionD->Draw("same");
1611      if(fEfficiencyIPNonhfeD)
1612        fEfficiencyIPNonhfeD->Draw("same");
1613    }
1614    TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
1615    legipeff->AddEntry(fEfficiencyBeautySigD,"IP Step Efficiency","");
1616    legipeff->AddEntry(fEfficiencyBeautySigD,"beauty e","p");
1617    legipeff->AddEntry(fEfficiencyBeautySigesdD,"beauty e(esd pt)","p");
1618    legipeff->AddEntry(fEfficiencyCharmSigD,"charm e","p");
1619    if(fBeamType == 0){
1620      legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
1621      legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
1622    }
1623    else{
1624      legipeff->AddEntry(fConversionEff,"conversion e","p");
1625      legipeff->AddEntry(fNonHFEEff,"Dalitz e","p");
1626    }
1627    legipeff->Draw("same");
1628    gPad->SetGrid();
1629    //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
1630 }
1631
1632 //____________________________________________________________________________
1633 THnSparse* AliHFEBeautySpectrum::GetBeautyIPEff(Bool_t isMCpt){
1634   //
1635   // Return beauty electron IP cut efficiency
1636   //
1637
1638   const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1639   Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
1640  
1641   Int_t nDim=1;  //dimensions of the efficiency weighting grid
1642   
1643   Int_t nBin[1] = {nPtbinning1};
1644  
1645   THnSparseF *ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
1646   
1647   ipcut->SetBinEdges(0, kPtRange);
1648   
1649   Double_t pt[1];
1650   Double_t weight;
1651   Double_t weightErr;
1652   Double_t contents[2];
1653
1654   weight = 1.0;
1655   weightErr = 1.0;
1656   
1657   Int_t looppt=nBin[0];
1658  
1659   Int_t ibin[2];
1660   
1661   for(int i=0; i<looppt; i++)
1662     {
1663       pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1664       if(isMCpt){
1665         if(fEfficiencyIPBeautyD){
1666           weight=fEfficiencyIPBeautyD->Eval(pt[0]);
1667           weightErr = 0;
1668         }
1669         else{
1670           printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1671           weight = fEfficiencyBeautySigD->GetBinContent(i+1); 
1672           weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1673         }
1674       }
1675       else{
1676         if(fEfficiencyIPBeautyesdD){
1677           weight=fEfficiencyIPBeautyesdD->Eval(pt[0]);
1678           weightErr = 0;
1679         }
1680         else{
1681           printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1682           weight = fEfficiencyBeautySigesdD->GetBinContent(i+1);
1683           weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1684         }
1685       }
1686       
1687       contents[0]=pt[0];
1688       ibin[0]=i+1;
1689       
1690       ipcut->Fill(contents,weight);
1691       ipcut->SetBinError(ibin,weightErr);
1692     }
1693   
1694   Int_t nDimSparse = ipcut->GetNdimensions();
1695   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1696   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1697
1698   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1699       binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
1700       nBins *= binsvar[iVar];
1701   }
1702
1703   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1704   // loop that sets 0 error in each bin
1705   for (Long_t iBin=0; iBin<nBins; iBin++) {
1706     Long_t bintmp = iBin ;
1707     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1708       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1709       bintmp /= binsvar[iVar] ;
1710     }
1711     //ipcut->SetBinError(binfill,0.); // put 0 everywhere
1712   }
1713
1714   delete[] binsvar;
1715   delete[] binfill;
1716
1717   return ipcut;
1718 }
1719
1720 //____________________________________________________________________________
1721 THnSparse* AliHFEBeautySpectrum::GetPIDxIPEff(Int_t source){
1722   //
1723   // Return PID x IP cut efficiency
1724   //
1725   const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1726   Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
1727   Int_t nDim=1;  //dimensions of the efficiency weighting grid
1728   
1729   Int_t nBin[1] = {nPtbinning1};
1730   
1731   THnSparseF *pideff;
1732   pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
1733   pideff->SetBinEdges(0, kPtRange);
1734   
1735   Double_t pt[1];
1736   Double_t weight;
1737   Double_t weightErr;
1738   Double_t contents[2];
1739   
1740   weight = 1.0;
1741   weightErr = 1.0;
1742   
1743   Int_t looppt=nBin[0];
1744   Int_t ibin[2];
1745       
1746   Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
1747   for(int i=0; i<looppt; i++)
1748     {
1749       pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1750       
1751       Double_t tofpideff = 1.;
1752       Double_t tofpideffesd = 1.;
1753       if(fEfficiencyTOFPIDD)
1754         tofpideff = fEfficiencyTOFPIDD->Eval(pt[0]); 
1755       else{
1756         printf("TOF PID fit failed on conversion. The result is wrong!\n");
1757       }  
1758       if(fEfficiencyesdTOFPIDD)
1759         tofpideffesd = fEfficiencyesdTOFPIDD->Eval(pt[0]);
1760       else{
1761         printf("TOF esd PID fit failed on conversion. The result is wrong!\n");
1762       }
1763       
1764       //tof pid eff x tpc pid eff x ip cut eff
1765       if(fIPParameterizedEff){
1766         if(source==0) {
1767           if(fEfficiencyIPCharmD){
1768             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1769             weightErr = 0; 
1770           }
1771           else{
1772             printf("Fit failed on charm IP cut efficiency\n");
1773             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1774             weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1); 
1775           }
1776         } 
1777         else if(source==2) {
1778           if(fEfficiencyIPConversionD){
1779             weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD->Eval(pt[0]); 
1780             weightErr = 0; 
1781           }
1782           else{
1783             printf("Fit failed on conversion IP cut efficiency\n");
1784             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1);
1785             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1786           }
1787         }
1788         else if(source==3) {
1789           if(fEfficiencyIPNonhfeD){
1790             weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD->Eval(pt[0]); 
1791             weightErr = 0; 
1792           }
1793           else{
1794             printf("Fit failed on dalitz IP cut efficiency\n");
1795             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1);
1796             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1797           }  
1798         }
1799       }
1800       else{
1801         if(source==0){ 
1802           if(fEfficiencyIPCharmD){
1803             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1804             weightErr = 0;
1805           }
1806           else{
1807             printf("Fit failed on charm IP cut efficiency\n");
1808             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1809             weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1810           }
1811         }
1812         else if(source==2){
1813           if(fBeamType==0){
1814             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
1815             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
1816           }
1817           else{
1818             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // conversion
1819             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1820           }
1821         }
1822         else if(source==3){
1823           if(fBeamType==0){
1824             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
1825             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
1826           }
1827           else{ 
1828             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // Dalitz
1829             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1830           }
1831         }
1832       }
1833       
1834       contents[0]=pt[0];
1835       ibin[0]=i+1;
1836       
1837       pideff->Fill(contents,weight);
1838       pideff->SetBinError(ibin,weightErr);
1839     }
1840   
1841   Int_t nDimSparse = pideff->GetNdimensions();
1842   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1843   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1844   
1845   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1846     binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
1847       nBins *= binsvar[iVar];
1848   }
1849
1850   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1851   // loop that sets 0 error in each bin
1852   for (Long_t iBin=0; iBin<nBins; iBin++) {
1853     Long_t bintmp = iBin ;
1854     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1855       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1856       bintmp /= binsvar[iVar] ;
1857     }
1858   }
1859
1860   delete[] binsvar;
1861   delete[] binfill;
1862
1863
1864   return pideff;
1865 }
1866
1867 //__________________________________________________________________________
1868 AliCFDataGrid *AliHFEBeautySpectrum::GetRawBspectra2ndMethod(){
1869  //
1870  // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
1871     //
1872   Int_t nDim = 1;
1873   
1874   const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
1875   Double_t kPtRange[19] = {0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 12, 16, 20};
1876   
1877   Int_t nBin[1] = {nPtbinning1};
1878   
1879   AliCFDataGrid *rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
1880       
1881   THnSparseF *brawspectra;
1882   brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
1883   
1884   brawspectra->SetBinEdges(0, kPtRange);
1885       
1886   Double_t pt[1];
1887   Double_t yields= 0.;
1888   Double_t valuesb[2];
1889             
1890   for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
1891     pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1892     
1893     yields = fBSpectrum2ndMethod->GetBinContent(i+1);
1894         
1895     valuesb[0]=pt[0];
1896     brawspectra->Fill(valuesb,yields);
1897   }
1898   
1899   
1900   
1901   Int_t nDimSparse = brawspectra->GetNdimensions();
1902   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1903   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1904     
1905   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1906     binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
1907     nBins *= binsvar[iVar];
1908   }
1909     
1910   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1911   // loop that sets 0 error in each bin
1912   for (Long_t iBin=0; iBin<nBins; iBin++) {
1913     Long_t bintmp = iBin ;
1914     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1915       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1916       bintmp /= binsvar[iVar] ;
1917     }
1918     brawspectra->SetBinError(binfill,0.); // put 0 everywhere
1919   }
1920   
1921   
1922   rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
1923   TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(0);
1924   
1925   new TCanvas;
1926   fBSpectrum2ndMethod->SetMarkerStyle(24);
1927   fBSpectrum2ndMethod->Draw("p");
1928   hRawBeautySpectra->SetMarkerStyle(25);
1929   hRawBeautySpectra->Draw("samep");
1930   
1931   delete[] binfill;
1932   delete[] binsvar; 
1933   
1934   return rawBeautyContainer;
1935 }
1936 //__________________________________________________________________________
1937 void AliHFEBeautySpectrum::CalculateNonHFEsyst(){
1938   //
1939   // Calculate non HFE sys
1940   //
1941   //
1942
1943   if(!fNonHFEsyst)
1944     return;
1945
1946   Double_t evtnorm[1] = {0.0};
1947   if(fNMCbgEvents>0) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
1948   
1949   AliCFDataGrid *convSourceGrid[kElecBgSources-3][kBgLevels];
1950   AliCFDataGrid *nonHFESourceGrid[kElecBgSources-3][kBgLevels];
1951
1952   AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
1953   AliCFDataGrid *bgNonHFEGrid[kBgLevels];
1954   AliCFDataGrid *bgConvGrid[kBgLevels];
1955
1956   Int_t stepbackground = 3;
1957   Int_t* bins=new Int_t[1];
1958   const Char_t *bgBase[2] = {"pi0","eta"};
1959  
1960   bins[0]=kSignalPtBins;//fConversionEff[centrality]->GetNbinsX();
1961    
1962   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1963   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1964
1965   for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){   
1966     for(Int_t iSource = 0; iSource < kElecBgSources-3; iSource++){
1967       convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
1968       weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1969       convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
1970       
1971       nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
1972       weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1973       nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
1974     }
1975     
1976     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
1977     for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){
1978       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
1979     }
1980     if(!fEtaSyst)
1981       bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
1982     
1983     bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone(); 
1984     for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
1985       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
1986     }
1987     if(!fEtaSyst)
1988       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
1989     
1990     bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
1991     bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
1992     if(fEtaSyst){
1993       bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
1994       bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
1995     }
1996   }
1997    
1998   //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0) 
1999   AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
2000   TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
2001   for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
2002     bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
2003     bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
2004     bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();    
2005     bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
2006
2007   //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
2008  
2009     hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
2010     hBaseErrors[iErr][0]->Scale(-1.);
2011     hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
2012     hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
2013     if(!fEtaSyst)break;
2014   }
2015    
2016   //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
2017   TH1D *hSpeciesErrors[kElecBgSources-4];
2018   for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2019     if(fEtaSyst && (iSource == 1))continue;
2020     hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
2021     TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
2022     hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
2023     hSpeciesErrors[iSource-1]->Scale(0.3);   
2024   }
2025   
2026   //Int_t firstBgSource = 0;//if eta systematics are not from scaling
2027   //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
2028   TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
2029   TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
2030   TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
2031
2032   TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
2033   TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
2034
2035   for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
2036     Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;    
2037     pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin); 
2038     pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
2039     if(fEtaSyst){
2040       etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin); 
2041       etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
2042     }
2043     else{ etaErrLow = etaErrUp = 0.;}
2044
2045     Double_t sqrsumErrs= 0;
2046     for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2047       if(fEtaSyst && (iSource == 1))continue;
2048       Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
2049       sqrsumErrs+=(scalingErr*scalingErr);
2050     }
2051     for(Int_t iErr = 0; iErr < 2; iErr++){
2052       for(Int_t iLevel = 0; iLevel < 2; iLevel++){
2053         hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
2054       }
2055       if(!fEtaSyst)break;
2056     }
2057     hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
2058     hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
2059     hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
2060
2061     hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
2062     hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));            
2063   }
2064      
2065   TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
2066   cPiErrors->cd();
2067   cPiErrors->SetLogx();
2068   cPiErrors->SetLogy();
2069   hBaseErrors[0][0]->Draw();
2070   //hBaseErrors[0][1]->SetMarkerColor(kBlack);
2071   //hBaseErrors[0][1]->SetLineColor(kBlack);
2072   //hBaseErrors[0][1]->Draw("SAME");
2073   if(fEtaSyst){
2074     hBaseErrors[1][0]->Draw("SAME");
2075     hBaseErrors[1][0]->SetMarkerColor(kBlack);
2076     hBaseErrors[1][0]->SetLineColor(kBlack);
2077   //hBaseErrors[1][1]->SetMarkerColor(13);
2078   //hBaseErrors[1][1]->SetLineColor(13);
2079   //hBaseErrors[1][1]->Draw("SAME");
2080   }
2081   //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
2082   //hOverallBinScaledErrsUp->SetLineColor(kBlue);
2083   //hOverallBinScaledErrsUp->Draw("SAME");
2084   hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
2085   hOverallBinScaledErrsLow->SetLineColor(kGreen);
2086   hOverallBinScaledErrsLow->Draw("SAME");
2087   hScalingErrors->SetLineColor(kBlue);
2088   hScalingErrors->Draw("SAME");
2089
2090   TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
2091   lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
2092   //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
2093   if(fEtaSyst){
2094   lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
2095   //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
2096   }
2097   lPiErr->AddEntry(hScalingErrors, "scaling error");
2098   lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
2099   //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
2100   lPiErr->Draw("SAME");
2101
2102   //Normalize errors
2103   TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
2104   TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
2105   hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
2106   hLowSystScaled->Scale(evtnorm[0]);
2107   TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
2108   TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
2109   //histograms to be normalized to TGraphErrors
2110   AliHFEtools::NormaliseBinWidth(hNormAllSystErrUp);
2111   AliHFEtools::NormaliseBinWidth(hNormAllSystErrLow);
2112
2113   TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
2114   cNormOvErrs->cd();
2115   cNormOvErrs->SetLogx();
2116   cNormOvErrs->SetLogy();
2117
2118   TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
2119   TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
2120   gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
2121   gOverallSystErrUp->SetMarkerColor(kBlack);
2122   gOverallSystErrUp->SetLineColor(kBlack);
2123   gOverallSystErrLow->SetMarkerColor(kRed);
2124   gOverallSystErrLow->SetLineColor(kRed);
2125   gOverallSystErrUp->Draw("AP");
2126   gOverallSystErrLow->Draw("PSAME");
2127   TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
2128   lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
2129   lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
2130   lAllSys->Draw("same");
2131
2132
2133   AliCFDataGrid *bgYieldGrid;
2134   if(fEtaSyst){
2135     bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
2136   }
2137   bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
2138
2139   TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
2140   TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
2141   hRelErrUp->Divide(hBgYield);
2142   TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
2143   hRelErrLow->Divide(hBgYield);
2144
2145   TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
2146   cRelErrs->cd();
2147   cRelErrs->SetLogx();
2148   hRelErrUp->SetTitle("Relative error of non-HFE background yield");
2149   hRelErrUp->Draw();
2150   hRelErrLow->SetLineColor(kBlack);
2151   hRelErrLow->Draw("SAME");
2152
2153   TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
2154   lRel->AddEntry(hRelErrUp, "upper");
2155   lRel->AddEntry(hRelErrLow, "lower");
2156   lRel->Draw("SAME");
2157
2158   //AliHFEtools::NormaliseBinWidth(hBgYield);
2159   //hBgYield->Scale(evtnorm[0]);
2160  
2161  
2162   //write histograms/TGraphs to file
2163   TFile *output = new TFile("systHists.root","recreate");
2164
2165   hBgYield->SetName("hBgYield");  
2166   hBgYield->Write();
2167   hRelErrUp->SetName("hRelErrUp");
2168   hRelErrUp->Write();
2169   hRelErrLow->SetName("hRelErrLow");
2170   hRelErrLow->Write();
2171   hUpSystScaled->SetName("hOverallSystErrUp");
2172   hUpSystScaled->Write();
2173   hLowSystScaled->SetName("hOverallSystErrLow");
2174   hLowSystScaled->Write();
2175   gOverallSystErrUp->SetName("gOverallSystErrUp");
2176   gOverallSystErrUp->Write();
2177   gOverallSystErrLow->SetName("gOverallSystErrLow");
2178   gOverallSystErrLow->Write(); 
2179
2180   output->Close(); 
2181   delete output;  
2182 }
2183 //____________________________________________________________
2184 void AliHFEBeautySpectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
2185   //
2186   // Unfold backgrounds to check its sanity
2187   //
2188
2189   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
2190   //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2191   if(!mcContainer){
2192     AliError("MC Container not available");
2193     return;
2194   }
2195
2196   if(!fCorrelation){
2197     AliError("No Correlation map available");
2198     return;
2199   }
2200
2201   // Data 
2202   AliCFDataGrid *dataGrid = 0x0;
2203   dataGrid = bgsubpectrum;
2204
2205   // Guessed
2206   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2207   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2208
2209   // Efficiency
2210   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2211   efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
2212
2213   // Unfold background spectra
2214   Int_t nDim=1;
2215   AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2216   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2217   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2218   if(fSetSmoothing) unfolding.UseSmoothing();
2219   unfolding.Unfold();
2220
2221   // Results
2222   THnSparse* result = unfolding.GetUnfolded();
2223   TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
2224   if(fBeamType==1)
2225   {
2226       ctest->Divide(2);
2227       ctest->cd(1);
2228       TH1D* htest1=(TH1D*)result->Projection(0);
2229       htest1->Draw();
2230       ctest->cd(2);
2231       TH1D* htest2=(TH1D*)result->Projection(1);
2232       htest2->Draw();
2233   }
2234
2235
2236   TGraphErrors* unfoldedbgspectrumD = Normalize(result);
2237   if(!unfoldedbgspectrumD) {
2238     AliError("Unfolded background spectrum doesn't exist");
2239   }
2240   else{
2241     TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
2242     if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
2243
2244     file->Close();
2245   }
2246 }