PWGHF/hfe/AliHFEspectrum.cxx

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