PWGHF/hfe/AliHFEspectrum.cxx

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