PWGHF/hfe/AliHFEspectrum.cxx

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