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       unnormalizedRawSpectrum->SetName("beautyAfterIP");
 828       unnormalizedRawSpectrum->Write();
 829
 830       if(gNormalizedRawSpectrum){
 831         gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
 832         gNormalizedRawSpectrum->Write();
 833       }
 834
 835       out->Close();
 836       delete out;
 837     }
 838   }
 839
 840   return kTRUE;
 841 }
 842
 843 //____________________________________________________________
 844 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
 845   //
 846   // Apply background subtraction
 847   //
 848
 849   // Raw spectrum
 850   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 851   if(!dataContainer){
 852     AliError("Data Container not available");
 853     return NULL;
 854   }
 855   printf("Step data: %d\n",fStepData);
 856   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
 857
 858   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 859   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
 860
 861   // Background Estimate
 862   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
 863   if(!backgroundContainer){
 864     AliError("MC background container not found");
 865     return NULL;
 866   }
 867
 868   Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
 869   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
 870
 871   if(!fInclusiveSpectrum){
 872     //Background subtraction for IP analysis
 873     TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(0);
 874     CorrectFromTheWidth(measuredTH1Draw);
 875     TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
 876     rawspectra->cd();
 877     rawspectra->SetLogx();
 878     rawspectra->SetLogy();
 879     TLegend *lRaw = new TLegend(0.65,0.65,0.95,0.95);
 880     measuredTH1Draw->SetMarkerStyle(20);
 881     measuredTH1Draw->Draw();
 882     lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
 883     if(fIPanaHadronBgSubtract){
 884       // Hadron background
 885       printf("Hadron background for IP analysis subtracted!\n");
 886       Int_t* bins=new Int_t[1];
 887       TH1D* htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
 888       bins[0]=htemp->GetNbinsX();
 889       AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
 890       hbgContainer->SetGrid(fHadronEffbyIPcut);
 891       backgroundGrid->Multiply(hbgContainer,1);
 892       // draw raw hadron bg spectra
 893       TH1D *hadronbg= (TH1D *) backgroundGrid->Project(0);
 894       CorrectFromTheWidth(hadronbg);
 895       hadronbg->SetMarkerColor(7);
 896       hadronbg->SetMarkerStyle(20);
 897       rawspectra->cd();
 898       hadronbg->Draw("samep");
 899       lRaw->AddEntry(hadronbg,"hadrons");
 900       // subtract hadron contamination
 901       spectrumSubtracted->Add(backgroundGrid,-1.0);
 902     }
 903     if(fIPanaCharmBgSubtract){
 904       // Charm background
 905       printf("Charm background for IP analysis subtracted!\n");
 906       AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
 907       // draw charm bg spectra
 908       TH1D *charmbg= (TH1D *) charmbgContainer->Project(0);
 909       CorrectFromTheWidth(charmbg);
 910       charmbg->SetMarkerColor(3);
 911       charmbg->SetMarkerStyle(20);
 912       rawspectra->cd();
 913       charmbg->Draw("samep");
 914       lRaw->AddEntry(charmbg,"charm elecs");
 915       // subtract charm background
 916       spectrumSubtracted->Add(charmbgContainer,-1.0);
 917     }
 918     if(fIPanaConversionBgSubtract){
 919       // Conversion background
 920       AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
 921       // draw conversion bg spectra
 922       TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(0);
 923       CorrectFromTheWidth(conversionbg);
 924       conversionbg->SetMarkerColor(4);
 925       conversionbg->SetMarkerStyle(20);
 926       rawspectra->cd();
 927       conversionbg->Draw("samep");
 928       lRaw->AddEntry(conversionbg,"conversion elecs");
 929       // subtract conversion background
 930       spectrumSubtracted->Add(conversionbgContainer,-1.0);
 931       printf("Conversion background subtraction is preliminary!\n");
 932     }
 933     if(fIPanaNonHFEBgSubtract){
 934       // NonHFE background
 935       AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
 936       // draw Dalitz/dielectron bg spectra
 937       TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(0);
 938       CorrectFromTheWidth(nonhfebg);
 939       nonhfebg->SetMarkerColor(6);
 940       nonhfebg->SetMarkerStyle(20);
 941       rawspectra->cd();
 942       nonhfebg->Draw("samep");
 943       lRaw->AddEntry(nonhfebg,"non-HF elecs");
 944       // subtract Dalitz/dielectron background
 945       spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
 946       printf("Non HFE background subtraction is preliminary!\n");
 947     }
 948     TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(0);
 949     CorrectFromTheWidth(rawbgsubtracted);
 950     rawbgsubtracted->SetMarkerStyle(24);
 951     rawspectra->cd();
 952     lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
 953     rawbgsubtracted->Draw("samep");
 954     lRaw->Draw("SAME");
 955   }
 956   else{
 957     // Subtract
 958     spectrumSubtracted->Add(backgroundGrid,-1.0);
 959   }
 960
 961   if(setBackground){
 962     if(fBackground) delete fBackground;
 963     fBackground = backgroundGrid;
 964   } else delete backgroundGrid;
 965
 966
 967   if(fDebugLevel > 0) {
 968
 969     Int_t ptpr;
 970     if(fBeamType==0) ptpr=0;
 971     if(fBeamType==1) ptpr=1;
 972
 973     TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
 974     cbackgroundsubtraction->Divide(3,1);
 975     cbackgroundsubtraction->cd(1);
 976     gPad->SetLogy();
 977     TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
 978     TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
 979     CorrectFromTheWidth(measuredTH1Daftersubstraction);
 980     CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
 981     measuredTH1Daftersubstraction->SetStats(0);
 982     measuredTH1Daftersubstraction->SetTitle("");
 983     measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 984     measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 985     measuredTH1Daftersubstraction->SetMarkerStyle(25);
 986     measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
 987     measuredTH1Daftersubstraction->SetLineColor(kBlack);
 988     measuredTH1Dbeforesubstraction->SetStats(0);
 989     measuredTH1Dbeforesubstraction->SetTitle("");
 990     measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 991     measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 992     measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
 993     measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
 994     measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
 995     measuredTH1Daftersubstraction->Draw();
 996     measuredTH1Dbeforesubstraction->Draw("same");
 997     TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
 998     legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
 999     legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
1000     legsubstraction->Draw("same");
1001     cbackgroundsubtraction->cd(2);
1002     gPad->SetLogy();
1003     TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
1004     ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
1005     ratiomeasuredcontamination->SetTitle("");
1006     ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1007     ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1008     ratiomeasuredcontamination->Sumw2();
1009     ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
1010     ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
1011     ratiomeasuredcontamination->SetStats(0);
1012     ratiomeasuredcontamination->SetMarkerStyle(26);
1013     ratiomeasuredcontamination->SetMarkerColor(kBlack);
1014     ratiomeasuredcontamination->SetLineColor(kBlack);
1015     for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
1016       ratiomeasuredcontamination->SetBinError(k+1,0.0);
1017     }
1018     ratiomeasuredcontamination->Draw("P");
1019     cbackgroundsubtraction->cd(3);
1020     TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
1021     CorrectFromTheWidth(measuredTH1background);
1022     measuredTH1background->SetStats(0);
1023     measuredTH1background->SetTitle("");
1024     measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1025     measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1026     measuredTH1background->SetMarkerStyle(26);
1027     measuredTH1background->SetMarkerColor(kBlack);
1028     measuredTH1background->SetLineColor(kBlack);
1029     measuredTH1background->Draw();
1030     if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
1031
1032     if(fBeamType==1) {
1033
1034       TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1035       cbackgrounde->Divide(2,1);
1036       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1037       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1038
1039       THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1040       TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1041       THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1042       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1043       Int_t nbbin = cenaxisb->GetNbins();
1044       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1045       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1046       for(Int_t binc = 0; binc < nbbin; binc++){
1047         TString titlee("BackgroundSubtraction_centrality_bin_");
1048         titlee += binc;
1049         TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1050         cbackground->Divide(2,1);
1051         cbackground->cd(1);
1052         gPad->SetLogy();
1053         cenaxisa->SetRange(binc+1,binc+1);
1054         cenaxisb->SetRange(binc+1,binc+1);
1055         TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1056         TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1057         CorrectFromTheWidth(aftersubstraction);
1058         CorrectFromTheWidth(beforesubstraction);
1059         aftersubstraction->SetStats(0);
1060         aftersubstraction->SetTitle((const char*)titlee);
1061         aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1062         aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1063         aftersubstraction->SetMarkerStyle(25);
1064         aftersubstraction->SetMarkerColor(kBlack);
1065         aftersubstraction->SetLineColor(kBlack);
1066         beforesubstraction->SetStats(0);
1067         beforesubstraction->SetTitle((const char*)titlee);
1068         beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1069         beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1070         beforesubstraction->SetMarkerStyle(24);
1071         beforesubstraction->SetMarkerColor(kBlue);
1072         beforesubstraction->SetLineColor(kBlue);
1073         aftersubstraction->Draw();
1074         beforesubstraction->Draw("same");
1075         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1076         lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1077         lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1078         lega->Draw("same");
1079         cbackgrounde->cd(1);
1080         gPad->SetLogy();
1081         TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1082         aftersubtractionn->SetMarkerStyle(stylee[binc]);
1083         aftersubtractionn->SetMarkerColor(colorr[binc]);
1084         if(binc==0) aftersubtractionn->Draw();
1085         else aftersubtractionn->Draw("same");
1086         legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1087         cbackgrounde->cd(2);
1088         gPad->SetLogy();
1089         TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1090         aftersubtractionng->SetMarkerStyle(stylee[binc]);
1091         aftersubtractionng->SetMarkerColor(colorr[binc]);
1092         if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1093         if(binc==0) aftersubtractionng->Draw();
1094         else aftersubtractionng->Draw("same");
1095         legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1096         cbackground->cd(2);
1097         TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1098         ratiocontamination->SetName("ratiocontamination");
1099         ratiocontamination->SetTitle((const char*)titlee);
1100         ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1101         ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1102         ratiocontamination->Add(aftersubstraction,-1.0);
1103         ratiocontamination->Divide(beforesubstraction);
1104         Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1105         for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1106           ratiocontamination->SetBinError(nbinpt+1,0.0);
1107         }
1108         ratiocontamination->SetStats(0);
1109         ratiocontamination->SetMarkerStyle(26);
1110         ratiocontamination->SetMarkerColor(kBlack);
1111         ratiocontamination->SetLineColor(kBlack);
1112         ratiocontamination->Draw("P");
1113       }
1114
1115       cbackgrounde->cd(1);
1116       legtotal->Draw("same");
1117       cbackgrounde->cd(2);
1118       legtotalg->Draw("same");
1119
1120       cenaxisa->SetRange(0,nbbin);
1121       cenaxisb->SetRange(0,nbbin);
1122       if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
1123     }
1124   }
1125
1126   return spectrumSubtracted;
1127 }
1128
1129 //____________________________________________________________
1130 AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1131   //
1132   // calculate charm background
1133   //
1134
1135   Double_t evtnorm=0;
1136   if(fNMCbgEvents) evtnorm= double(fNEvents[0])/double(fNMCbgEvents);
1137
1138   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1139   if(!mcContainer){
1140     AliError("MC Container not available");
1141     return NULL;
1142   }
1143
1144   if(!fCorrelation){
1145     AliError("No Correlation map available");
1146     return NULL;
1147   }
1148
1149   AliCFDataGrid *charmBackgroundGrid= 0x0;
1150   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1151   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1152
1153   Int_t* bins=new Int_t[1];
1154   bins[0]=charmbgaftertofpid->GetNbinsX();
1155   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",1,bins);
1156   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1157   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
1158
1159   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
1160   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
1161
1162   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
1163   weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1164   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
1165   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1166   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
1167
1168   // Efficiency (set efficiency to 1 for only folding)
1169   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1170   efficiencyD->CalculateEfficiency(0,0);
1171
1172   // Folding
1173   Int_t nDim = 1;
1174   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1175   folding.SetMaxNumberOfIterations(1);
1176   folding.Unfold();
1177
1178   // Results
1179   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1180   THnSparse* result=(THnSparse*)result1->Clone();
1181   charmBackgroundGrid->SetGrid(result);
1182   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
1183
1184   //Charm background evaluation plots
1185
1186   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1187   cCharmBgEval->Divide(3,1);
1188
1189   cCharmBgEval->cd(1);
1190   charmbgaftertofpid->Scale(evtnorm);
1191   CorrectFromTheWidth(charmbgaftertofpid);
1192   charmbgaftertofpid->SetMarkerStyle(25);
1193   charmbgaftertofpid->Draw("p");
1194   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1195   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1196   gPad->SetLogy();
1197
1198   CorrectFromTheWidth(charmbgafteripcut);
1199   charmbgafteripcut->SetMarkerStyle(24);
1200   charmbgafteripcut->Draw("samep");
1201
1202   CorrectFromTheWidth(charmbgafterweight);
1203   charmbgafterweight->SetMarkerStyle(24);
1204   charmbgafterweight->SetMarkerColor(4);
1205   charmbgafterweight->Draw("samep");
1206
1207   CorrectFromTheWidth(charmbgafterfolding);
1208   charmbgafterfolding->SetMarkerStyle(24);
1209   charmbgafterfolding->SetMarkerColor(2);
1210   charmbgafterfolding->Draw("samep");
1211
1212   cCharmBgEval->cd(2);
1213   parametrizedcharmpidipeff->SetMarkerStyle(24);
1214   parametrizedcharmpidipeff->Draw("p");
1215   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1216
1217   cCharmBgEval->cd(3);
1218   charmweightingfc->SetMarkerStyle(24);
1219   charmweightingfc->Draw("p");
1220   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1221   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1222
1223   cCharmBgEval->cd(1);
1224   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1225   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1226   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1227   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1228   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1229   legcharmbg->Draw("same");
1230
1231   cCharmBgEval->cd(2);
1232   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1233   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1234   legcharmbg2->Draw("same");
1235
1236   CorrectStatErr(charmBackgroundGrid);
1237   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
1238
1239   return charmBackgroundGrid;
1240 }
1241
1242 //____________________________________________________________
1243 AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1244   //
1245   // calculate conversion background
1246   //
1247
1248   Double_t evtnorm[1] = {0.0};
1249   if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1250   printf("check event!!! %lf \n",evtnorm[0]);
1251
1252   AliCFContainer *backgroundContainer = 0x0;
1253
1254   if(fNonHFEsyst){
1255     backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0]->Clone();
1256     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1257       backgroundContainer->Add(fConvSourceContainer[iSource][0]);
1258     }
1259   }
1260   else{
1261     // Background Estimate
1262     backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
1263   }
1264   if(!backgroundContainer){
1265     AliError("MC background container not found");
1266     return NULL;
1267   }
1268
1269   Int_t stepbackground = 3;
1270   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1271   //backgroundGrid->Scale(evtnorm);
1272   //workaround for statistical errors: "Scale" method does not work for our errors, since Sumw2 is not defined for AliCFContainer
1273   Int_t *nBin=new Int_t[1];
1274
1275   Int_t* bins=new Int_t[1];
1276   bins[0]=fConversionEff->GetNbinsX();
1277
1278   for(Long_t iBin=1; iBin<= bins[0];iBin++){
1279     nBin[0]=iBin;
1280     backgroundGrid->SetElementError(nBin, backgroundGrid->GetElementError(nBin)*evtnorm[0]);
1281     backgroundGrid->SetElement(nBin,backgroundGrid->GetElement(nBin)*evtnorm[0]);
1282   }
1283   //end of workaround for statistical errors
1284
1285   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1286   weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1287   backgroundGrid->Multiply(weightedConversionContainer,1.0);
1288
1289   return backgroundGrid;
1290 }
1291
1292
1293 //____________________________________________________________
1294 AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1295   //
1296   // calculate non-HFE background
1297   //
1298
1299   Double_t evtnorm[1] = {0.0};
1300   if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
1301   printf("check event!!! %lf \n",evtnorm[0]);
1302
1303   AliCFContainer *backgroundContainer = 0x0;
1304   if(fNonHFEsyst){
1305     backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0]->Clone();
1306     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1307       backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);
1308     }
1309   }
1310   else{
1311     // Background Estimate
1312     backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1313   }
1314   if(!backgroundContainer){
1315     AliError("MC background container not found");
1316     return NULL;
1317   }
1318
1319
1320   Int_t stepbackground = 3;
1321   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1322   //backgroundGrid->Scale(evtnorm);
1323   //workaround for statistical errors: "Scale" method does not work for our errors, since Sumw2 is not defined for AliCFContainer
1324   Int_t *nBin=new Int_t[1];
1325
1326   Int_t* bins=new Int_t[1];
1327   bins[0]=fConversionEff->GetNbinsX();
1328
1329   for(Long_t iBin=1; iBin<= bins[0];iBin++){
1330     nBin[0]=iBin;
1331     backgroundGrid->SetElementError(nBin, backgroundGrid->GetElementError(nBin)*evtnorm[0]);
1332     backgroundGrid->SetElement(nBin,backgroundGrid->GetElement(nBin)*evtnorm[0]);
1333   }
1334   //end of workaround for statistical errors
1335
1336   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1337   weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1338   backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
1339
1340   return backgroundGrid;
1341 }
1342
1343 //____________________________________________________________
1344 AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1345
1346   //
1347   // Apply TPC pid efficiency correction from parametrisation
1348   //
1349
1350   // Data in the right format
1351   AliCFDataGrid *dataGrid = 0x0;
1352   if(bgsubpectrum) {
1353     dataGrid = bgsubpectrum;
1354   }
1355   else {
1356
1357     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1358     if(!dataContainer){
1359       AliError("Data Container not available");
1360       return NULL;
1361     }
1362
1363     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1364   }
1365
1366   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1367   result->SetName("ParametrizedEfficiencyBefore");
1368   THnSparse *h = result->GetGrid();
1369   Int_t nbdimensions = h->GetNdimensions();
1370   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1371
1372   AliCFContainer *dataContainer = GetContainer(kDataContainer);
1373   if(!dataContainer){
1374     AliError("Data Container not available");
1375     return NULL;
1376   }
1377   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1378   dataContainerbis->Add(dataContainerbis,-1.0);
1379
1380
1381   Int_t* coord = new Int_t[nbdimensions];
1382   memset(coord, 0, sizeof(Int_t) * nbdimensions);
1383   Double_t* points = new Double_t[nbdimensions];
1384
1385
1386   ULong64_t nEntries = h->GetNbins();
1387   for (ULong64_t i = 0; i < nEntries; ++i) {
1388
1389     Double_t value = h->GetBinContent(i, coord);
1390     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1391     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1392
1393     // Get the bin co-ordinates given an coord
1394     for (Int_t j = 0; j < nbdimensions; ++j)
1395       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1396
1397     if (!fEfficiencyFunction->IsInside(points))
1398          continue;
1399     TF1::RejectPoint(kFALSE);
1400
1401     // Evaulate function at points
1402     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1403     //printf("Value efficiency is %f\n",valueEfficiency);
1404
1405     if(valueEfficiency > 0.0) {
1406       h->SetBinContent(coord,value/valueEfficiency);
1407       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1408     }
1409     Double_t error = h->GetBinError(i);
1410     h->SetBinError(coord,error/valueEfficiency);
1411     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1412
1413
1414   }
1415
1416   delete[] coord;
1417   delete[] points;
1418
1419   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1420
1421   if(fDebugLevel > 0) {
1422
1423     TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1424     cEfficiencyParametrized->Divide(2,1);
1425     cEfficiencyParametrized->cd(1);
1426     TH1D *afterE = (TH1D *) resultt->Project(0);
1427     TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1428     CorrectFromTheWidth(afterE);
1429     CorrectFromTheWidth(beforeE);
1430     afterE->SetStats(0);
1431     afterE->SetTitle("");
1432     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1433     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1434     afterE->SetMarkerStyle(25);
1435     afterE->SetMarkerColor(kBlack);
1436     afterE->SetLineColor(kBlack);
1437     beforeE->SetStats(0);
1438     beforeE->SetTitle("");
1439     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1440     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1441     beforeE->SetMarkerStyle(24);
1442     beforeE->SetMarkerColor(kBlue);
1443     beforeE->SetLineColor(kBlue);
1444     gPad->SetLogy();
1445     afterE->Draw();
1446     beforeE->Draw("same");
1447     TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1448     legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1449     legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1450     legefficiencyparametrized->Draw("same");
1451     cEfficiencyParametrized->cd(2);
1452     fEfficiencyFunction->Draw();
1453     //cEfficiencyParametrized->cd(3);
1454     //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1455     //ratioefficiency->Divide(afterE);
1456     //ratioefficiency->Draw();
1457
1458     if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
1459   }
1460
1461
1462   return resultt;
1463
1464 }
1465 //____________________________________________________________
1466 AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1467
1468   //
1469   // Apply TPC pid efficiency correction from V0
1470   //
1471
1472   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1473   if(!v0Container){
1474     AliError("V0 Container not available");
1475     return NULL;
1476   }
1477
1478   // Efficiency
1479   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1480   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1481
1482   // Data in the right format
1483   AliCFDataGrid *dataGrid = 0x0;
1484   if(bgsubpectrum) {
1485     dataGrid = bgsubpectrum;
1486   }
1487   else {
1488
1489     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1490     if(!dataContainer){
1491       AliError("Data Container not available");
1492       return NULL;
1493     }
1494
1495     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1496   }
1497
1498   // Correct
1499   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1500   result->ApplyEffCorrection(*efficiencyD);
1501
1502   if(fDebugLevel > 0) {
1503
1504     Int_t ptpr;
1505     if(fBeamType==0) ptpr=0;
1506     if(fBeamType==1) ptpr=1;
1507
1508     TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1509     cV0Efficiency->Divide(2,1);
1510     cV0Efficiency->cd(1);
1511     TH1D *afterE = (TH1D *) result->Project(ptpr);
1512     TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
1513     afterE->SetStats(0);
1514     afterE->SetTitle("");
1515     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1516     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1517     afterE->SetMarkerStyle(25);
1518     afterE->SetMarkerColor(kBlack);
1519     afterE->SetLineColor(kBlack);
1520     beforeE->SetStats(0);
1521     beforeE->SetTitle("");
1522     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1523     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1524     beforeE->SetMarkerStyle(24);
1525     beforeE->SetMarkerColor(kBlue);
1526     beforeE->SetLineColor(kBlue);
1527     afterE->Draw();
1528     beforeE->Draw("same");
1529     TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1530     legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1531     legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1532     legV0efficiency->Draw("same");
1533     cV0Efficiency->cd(2);
1534     TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
1535     efficiencyDproj->SetTitle("");
1536     efficiencyDproj->SetStats(0);
1537     efficiencyDproj->SetMarkerStyle(25);
1538     efficiencyDproj->Draw();
1539
1540     if(fBeamType==1) {
1541
1542       TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1543       cV0Efficiencye->Divide(2,1);
1544       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1545       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1546
1547       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1548       TAxis *cenaxisa = sparseafter->GetAxis(0);
1549       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1550       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1551       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1552       TAxis *cenaxisc = efficiencya->GetAxis(0);
1553       Int_t nbbin = cenaxisb->GetNbins();
1554       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1555       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1556       for(Int_t binc = 0; binc < nbbin; binc++){
1557         TString titlee("V0Efficiency_centrality_bin_");
1558         titlee += binc;
1559         TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1560         ccV0Efficiency->Divide(2,1);
1561         ccV0Efficiency->cd(1);
1562         gPad->SetLogy();
1563         cenaxisa->SetRange(binc+1,binc+1);
1564         cenaxisb->SetRange(binc+1,binc+1);
1565         cenaxisc->SetRange(binc+1,binc+1);
1566         TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1567         TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1568         CorrectFromTheWidth(aftere);
1569         CorrectFromTheWidth(beforee);
1570         aftere->SetStats(0);
1571         aftere->SetTitle((const char*)titlee);
1572         aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1573         aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1574         aftere->SetMarkerStyle(25);
1575         aftere->SetMarkerColor(kBlack);
1576         aftere->SetLineColor(kBlack);
1577         beforee->SetStats(0);
1578         beforee->SetTitle((const char*)titlee);
1579         beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1580         beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1581         beforee->SetMarkerStyle(24);
1582         beforee->SetMarkerColor(kBlue);
1583         beforee->SetLineColor(kBlue);
1584         aftere->Draw();
1585         beforee->Draw("same");
1586         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1587         lega->AddEntry(beforee,"Before correction","p");
1588         lega->AddEntry(aftere,"After correction","p");
1589         lega->Draw("same");
1590         cV0Efficiencye->cd(1);
1591         gPad->SetLogy();
1592         TH1D *afteree = (TH1D *) aftere->Clone();
1593         afteree->SetMarkerStyle(stylee[binc]);
1594         afteree->SetMarkerColor(colorr[binc]);
1595         if(binc==0) afteree->Draw();
1596         else afteree->Draw("same");
1597         legtotal->AddEntry(afteree,(const char*) titlee,"p");
1598         cV0Efficiencye->cd(2);
1599         gPad->SetLogy();
1600         TH1D *aftereeu = (TH1D *) aftere->Clone();
1601         aftereeu->SetMarkerStyle(stylee[binc]);
1602         aftereeu->SetMarkerColor(colorr[binc]);
1603         if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
1604         if(binc==0) aftereeu->Draw();
1605         else aftereeu->Draw("same");
1606         legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
1607         ccV0Efficiency->cd(2);
1608         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
1609         efficiencyDDproj->SetTitle("");
1610         efficiencyDDproj->SetStats(0);
1611         efficiencyDDproj->SetMarkerStyle(25);
1612         efficiencyDDproj->Draw();
1613       }
1614
1615       cV0Efficiencye->cd(1);
1616       legtotal->Draw("same");
1617       cV0Efficiencye->cd(2);
1618       legtotalg->Draw("same");
1619
1620       cenaxisa->SetRange(0,nbbin);
1621       cenaxisb->SetRange(0,nbbin);
1622       cenaxisc->SetRange(0,nbbin);
1623
1624       if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
1625     }
1626
1627   }
1628
1629
1630   return result;
1631
1632 }
1633 //____________________________________________________________
1634 TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1635
1636   //
1637   // Unfold and eventually correct for efficiency the bgsubspectrum
1638   //
1639
1640   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1641   if(!mcContainer){
1642     AliError("MC Container not available");
1643     return NULL;
1644   }
1645
1646   if(!fCorrelation){
1647     AliError("No Correlation map available");
1648     return NULL;
1649   }
1650
1651   // Data
1652   AliCFDataGrid *dataGrid = 0x0;
1653   if(bgsubpectrum) {
1654     dataGrid = bgsubpectrum;
1655   }
1656   else {
1657
1658     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1659     if(!dataContainer){
1660       AliError("Data Container not available");
1661       return NULL;
1662     }
1663
1664     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1665   }
1666
1667   // Guessed
1668   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1669   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1670
1671   // Efficiency
1672   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1673   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1674
1675   // Consider parameterized IP cut efficiency
1676   if(!fInclusiveSpectrum){
1677     Int_t* bins=new Int_t[1];
1678     bins[0]=35;
1679
1680     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1681     beffContainer->SetGrid(GetBeautyIPEff());
1682     efficiencyD->Multiply(beffContainer,1);
1683   }
1684
1685
1686   // Unfold
1687
1688   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
1689   //unfolding.SetUseCorrelatedErrors();
1690   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1691   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1692   if(fSetSmoothing) unfolding.UseSmoothing();
1693   unfolding.Unfold();
1694
1695   // Results
1696   THnSparse* result = unfolding.GetUnfolded();
1697   THnSparse* residual = unfolding.GetEstMeasured();
1698   TList *listofresults = new TList;
1699   listofresults->SetOwner();
1700   listofresults->AddAt((THnSparse*)result->Clone(),0);
1701   listofresults->AddAt((THnSparse*)residual->Clone(),1);
1702
1703   if(fDebugLevel > 0) {
1704
1705     Int_t ptpr;
1706     if(fBeamType==0) ptpr=0;
1707     if(fBeamType==1) ptpr=1;
1708
1709     TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
1710     cresidual->Divide(2,1);
1711     cresidual->cd(1);
1712     gPad->SetLogy();
1713     TGraphErrors* residualspectrumD = Normalize(residual);
1714     if(!residualspectrumD) {
1715       AliError("Number of Events not set for the normalization");
1716       return NULL;
1717     }
1718     residualspectrumD->SetTitle("");
1719     residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
1720     residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1721     residualspectrumD->SetMarkerStyle(26);
1722     residualspectrumD->SetMarkerColor(kBlue);
1723     residualspectrumD->SetLineColor(kBlue);
1724     residualspectrumD->Draw("AP");
1725     AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
1726     dataGridBis->SetName("dataGridBis");
1727     TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
1728     measuredspectrumD->SetTitle("");
1729     measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
1730     measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
1731     measuredspectrumD->SetMarkerStyle(25);
1732     measuredspectrumD->SetMarkerColor(kBlack);
1733     measuredspectrumD->SetLineColor(kBlack);
1734     measuredspectrumD->Draw("P");
1735     TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
1736     legres->AddEntry(residualspectrumD,"Residual","p");
1737     legres->AddEntry(measuredspectrumD,"Measured","p");
1738     legres->Draw("same");
1739     cresidual->cd(2);
1740     TH1D *residualTH1D = residual->Projection(ptpr);
1741     TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
1742     TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
1743     ratioresidual->SetName("ratioresidual");
1744     ratioresidual->SetTitle("");
1745     ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
1746     ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1747     ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
1748     ratioresidual->SetStats(0);
1749     ratioresidual->Draw();
1750
1751     if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
1752   }
1753
1754   return listofresults;
1755
1756 }
1757
1758 //____________________________________________________________
1759 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1760
1761   //
1762   // Apply unfolding and efficiency correction together to bgsubspectrum
1763   //
1764
1765   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1766   if(!mcContainer){
1767     AliError("MC Container not available");
1768     return NULL;
1769   }
1770
1771   // Efficiency
1772   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1773   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1774
1775   // Consider parameterized IP cut efficiency
1776   if(!fInclusiveSpectrum){
1777     Int_t* bins=new Int_t[1];
1778     bins[0]=35;
1779
1780     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1781     beffContainer->SetGrid(GetBeautyIPEff());
1782     efficiencyD->Multiply(beffContainer,1);
1783   }
1784
1785   // Data in the right format
1786   AliCFDataGrid *dataGrid = 0x0;
1787   if(bgsubpectrum) {
1788     dataGrid = bgsubpectrum;
1789   }
1790   else {
1791
1792     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1793     if(!dataContainer){
1794       AliError("Data Container not available");
1795       return NULL;
1796     }
1797
1798     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1799   }
1800
1801   // Correct
1802   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1803   result->ApplyEffCorrection(*efficiencyD);
1804
1805   if(fDebugLevel > 0) {
1806
1807     Int_t ptpr;
1808     if(fBeamType==0) ptpr=0;
1809     if(fBeamType==1) ptpr=1;
1810
1811     printf("Step MC: %d\n",fStepMC);
1812     printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1813     printf("Step MC true: %d\n",fStepTrue);
1814     AliCFEffGrid  *efficiencymcPID = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
1815     AliCFEffGrid  *efficiencymctrackinggeo = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
1816     AliCFEffGrid  *efficiencymcall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
1817
1818     AliCFEffGrid  *efficiencyesdall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
1819
1820     TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
1821     cefficiency->cd(1);
1822     TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
1823     efficiencymcPIDD->SetTitle("");
1824     efficiencymcPIDD->SetStats(0);
1825     efficiencymcPIDD->SetMarkerStyle(25);
1826     efficiencymcPIDD->Draw();
1827     TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
1828     efficiencymctrackinggeoD->SetTitle("");
1829     efficiencymctrackinggeoD->SetStats(0);
1830     efficiencymctrackinggeoD->SetMarkerStyle(26);
1831     efficiencymctrackinggeoD->Draw("same");
1832     TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
1833     efficiencymcallD->SetTitle("");
1834     efficiencymcallD->SetStats(0);
1835     efficiencymcallD->SetMarkerStyle(27);
1836     efficiencymcallD->Draw("same");
1837     TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
1838     efficiencyesdallD->SetTitle("");
1839     efficiencyesdallD->SetStats(0);
1840     efficiencyesdallD->SetMarkerStyle(24);
1841     efficiencyesdallD->Draw("same");
1842     TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
1843     legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
1844     legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
1845     legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
1846     legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
1847     legeff->Draw("same");
1848
1849     if(fBeamType==1) {
1850
1851       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1852       TAxis *cenaxisa = sparseafter->GetAxis(0);
1853       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1854       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1855       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1856       TAxis *cenaxisc = efficiencya->GetAxis(0);
1857       //Int_t nbbin = cenaxisb->GetNbins();
1858       //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1859       //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1860       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
1861         TString titlee("Efficiency_centrality_bin_");
1862         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
1863         titlee += "_";
1864         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
1865         TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1866         cefficiencye->Divide(2,1);
1867         cefficiencye->cd(1);
1868         gPad->SetLogy();
1869         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1870         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
1871         TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
1872         TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
1873         CorrectFromTheWidth(afterefficiency);
1874         CorrectFromTheWidth(beforeefficiency);
1875         afterefficiency->SetStats(0);
1876         afterefficiency->SetTitle((const char*)titlee);
1877         afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1878         afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1879         afterefficiency->SetMarkerStyle(25);
1880         afterefficiency->SetMarkerColor(kBlack);
1881         afterefficiency->SetLineColor(kBlack);
1882         beforeefficiency->SetStats(0);
1883         beforeefficiency->SetTitle((const char*)titlee);
1884         beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1885         beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1886         beforeefficiency->SetMarkerStyle(24);
1887         beforeefficiency->SetMarkerColor(kBlue);
1888         beforeefficiency->SetLineColor(kBlue);
1889         afterefficiency->Draw();
1890         beforeefficiency->Draw("same");
1891         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1892         lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
1893         lega->AddEntry(afterefficiency,"After efficiency correction","p");
1894         lega->Draw("same");
1895         cefficiencye->cd(2);
1896         cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
1897         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
1898         efficiencyDDproj->SetTitle("");
1899         efficiencyDDproj->SetStats(0);
1900         efficiencyDDproj->SetMarkerStyle(25);
1901         efficiencyDDproj->SetMarkerColor(2);
1902         efficiencyDDproj->Draw();
1903         cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
1904         TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
1905         efficiencyDDproja->SetTitle("");
1906         efficiencyDDproja->SetStats(0);
1907         efficiencyDDproja->SetMarkerStyle(26);
1908         efficiencyDDproja->SetMarkerColor(4);
1909         efficiencyDDproja->Draw("same");
1910       }
1911     }
1912
1913     if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
1914   }
1915
1916   return result;
1917
1918 }
1919
1920 //__________________________________________________________________________________
1921 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
1922   //
1923   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1924   // Give the final pt spectrum to be compared
1925   //
1926
1927   if(fNEvents[i] > 0) {
1928
1929     Int_t ptpr = 0;
1930     if(fBeamType==0) ptpr=0;
1931     if(fBeamType==1) ptpr=1;
1932
1933     TH1D* projection = spectrum->Projection(ptpr);
1934     CorrectFromTheWidth(projection);
1935     TGraphErrors *graphError = NormalizeTH1(projection,i);
1936     return graphError;
1937
1938   }
1939
1940   return 0x0;
1941
1942
1943 }
1944 //__________________________________________________________________________________
1945 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
1946   //
1947   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1948   // Give the final pt spectrum to be compared
1949   //
1950   if(fNEvents[i] > 0) {
1951
1952     Int_t ptpr=0;
1953     if(fBeamType==0) ptpr=0;
1954     if(fBeamType==1) ptpr=1;
1955
1956     TH1D* projection = (TH1D *) spectrum->Project(ptpr);
1957     CorrectFromTheWidth(projection);
1958     TGraphErrors *graphError = NormalizeTH1(projection,i);
1959     return graphError;
1960
1961   }
1962
1963   return 0x0;
1964
1965
1966 }
1967 //__________________________________________________________________________________
1968 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
1969   //
1970   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
1971   // Give the final pt spectrum to be compared
1972   //
1973   Double_t chargecoefficient = 0.5;
1974   if(fChargeChoosen != 0) chargecoefficient = 1.0;
1975
1976   Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
1977   printf("Normalizing Eta Range %f\n", etarange);
1978   if(fNEvents[i] > 0) {
1979
1980     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
1981     Double_t p = 0, dp = 0; Int_t point = 1;
1982     Double_t n = 0, dN = 0;
1983     Double_t nCorr = 0, dNcorr = 0;
1984     Double_t errdN = 0, errdp = 0;
1985     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
1986       point = ibin - input->GetXaxis()->GetFirst();
1987       p = input->GetXaxis()->GetBinCenter(ibin);
1988       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
1989       n = input->GetBinContent(ibin);
1990       dN = input->GetBinError(ibin);
1991
1992       // New point
1993       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
1994       errdN = 1./(2. * TMath::Pi() * p);
1995       errdp = 1./(2. * TMath::Pi() * p*p) * n;
1996       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
1997
1998       spectrumNormalized->SetPoint(point, p, nCorr);
1999       spectrumNormalized->SetPointError(point, dp, dNcorr);
2000     }
2001     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2002     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2003     spectrumNormalized->SetMarkerStyle(22);
2004     spectrumNormalized->SetMarkerColor(kBlue);
2005     spectrumNormalized->SetLineColor(kBlue);
2006
2007     return spectrumNormalized;
2008
2009   }
2010
2011   return 0x0;
2012
2013
2014 }
2015 //__________________________________________________________________________________
2016 TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
2017   //
2018   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2019   // Give the final pt spectrum to be compared
2020   //
2021   Double_t chargecoefficient = 0.5;
2022   if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
2023
2024   Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2025   printf("Normalizing Eta Range %f\n", etarange);
2026   if(normalization > 0) {
2027
2028     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2029     Double_t p = 0, dp = 0; Int_t point = 1;
2030     Double_t n = 0, dN = 0;
2031     Double_t nCorr = 0, dNcorr = 0;
2032     Double_t errdN = 0, errdp = 0;
2033     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2034       point = ibin - input->GetXaxis()->GetFirst();
2035       p = input->GetXaxis()->GetBinCenter(ibin);
2036       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2037       n = input->GetBinContent(ibin);
2038       dN = input->GetBinError(ibin);
2039
2040       // New point
2041       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
2042       errdN = 1./(2. * TMath::Pi() * p);
2043       errdp = 1./(2. * TMath::Pi() * p*p) * n;
2044       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2045
2046       spectrumNormalized->SetPoint(point, p, nCorr);
2047       spectrumNormalized->SetPointError(point, dp, dNcorr);
2048     }
2049     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2050     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2051     spectrumNormalized->SetMarkerStyle(22);
2052     spectrumNormalized->SetMarkerColor(kBlue);
2053     spectrumNormalized->SetLineColor(kBlue);
2054
2055     return spectrumNormalized;
2056
2057   }
2058
2059   return 0x0;
2060
2061
2062 }
2063 //____________________________________________________________
2064 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2065   //
2066   // Set the container for a given step to the
2067   //
2068   if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
2069   fCFContainers->AddAt(cont, type);
2070 }
2071
2072 //____________________________________________________________
2073 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2074   //
2075   // Get Correction Framework Container for given type
2076   //
2077   if(!fCFContainers) return NULL;
2078   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2079 }
2080 //____________________________________________________________
2081 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge) {
2082   //
2083   // Slice bin for a given source of electron
2084   // nDim is the number of dimension the corrections are done
2085   // dimensions are the definition of the dimensions
2086   // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2087   // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
2088   //
2089
2090   Double_t *varMin = new Double_t[container->GetNVar()],
2091            *varMax = new Double_t[container->GetNVar()];
2092
2093   Double_t *binLimits;
2094   for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2095
2096     binLimits = new Double_t[container->GetNBins(ivar)+1];
2097     container->GetBinLimits(ivar,binLimits);
2098     varMin[ivar] = binLimits[0];
2099     varMax[ivar] = binLimits[container->GetNBins(ivar)];
2100     // source
2101     if(ivar == 4){
2102       if((source>= 0) && (source<container->GetNBins(ivar))) {
2103               varMin[ivar] = binLimits[source];
2104               varMax[ivar] = binLimits[source];
2105       }
2106     }
2107     // charge
2108     if(ivar == 3) {
2109       if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
2110     }
2111     // eta
2112     if(ivar == 1){
2113       for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
2114         AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
2115       if(fEtaSelected){
2116         varMin[ivar] = fEtaRange[0];
2117         varMax[ivar] = fEtaRange[1];
2118       }
2119     }
2120     if(fEtaSelected){
2121       fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2122       fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2123       AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2124     }
2125
2126
2127     delete[] binLimits;
2128   }
2129
2130   AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2131   AddTemporaryObject(k);
2132   delete[] varMin; delete[] varMax;
2133
2134   return k;
2135
2136 }
2137
2138 //_________________________________________________________________________
2139 THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
2140   //
2141   // Slice correlation
2142   //
2143
2144   Int_t ndimensions = correlationmatrix->GetNdimensions();
2145   //printf("Number of dimension %d correlation map\n",ndimensions);
2146   if(ndimensions < (2*nDim)) {
2147     AliError("Problem in the dimensions");
2148     return NULL;
2149   }
2150   Int_t ndimensionsContainer = (Int_t) ndimensions/2;
2151   //printf("Number of dimension %d container\n",ndimensionsContainer);
2152
2153   Int_t *dim = new Int_t[nDim*2];
2154   for(Int_t iter=0; iter < nDim; iter++){
2155     dim[iter] = dimensions[iter];
2156     dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
2157     //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
2158   }
2159
2160   THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
2161
2162   delete[] dim;
2163   return k;
2164
2165 }
2166 //___________________________________________________________________________
2167 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2168   //
2169   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2170   //
2171
2172   TAxis *axis = h1->GetXaxis();
2173   Int_t nbinX = h1->GetNbinsX();
2174
2175   for(Int_t i = 1; i <= nbinX; i++) {
2176
2177     Double_t width = axis->GetBinWidth(i);
2178     Double_t content = h1->GetBinContent(i);
2179     Double_t error = h1->GetBinError(i);
2180     h1->SetBinContent(i,content/width);
2181     h1->SetBinError(i,error/width);
2182   }
2183
2184 }
2185
2186 //___________________________________________________________________________
2187 void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2188   //
2189   // Correct statistical error
2190   //
2191
2192   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2193   Int_t nbinX = h1->GetNbinsX();
2194   Int_t bins[1];
2195   for(Long_t i = 1; i <= nbinX; i++) {
2196     bins[0] = i;
2197     Float_t content = h1->GetBinContent(i);
2198     if(content>0){
2199       Float_t error = TMath::Sqrt(content);
2200       backgroundGrid->SetElementError(bins, error);
2201     }
2202   }
2203 }
2204
2205 //____________________________________________________________
2206 void AliHFEspectrum::AddTemporaryObject(TObject *o){
2207   //
2208   // Emulate garbage collection on class level
2209   //
2210   if(!fTemporaryObjects) {
2211     fTemporaryObjects= new TList;
2212     fTemporaryObjects->SetOwner();
2213   }
2214   fTemporaryObjects->Add(o);
2215 }
2216
2217 //____________________________________________________________
2218 void AliHFEspectrum::ClearObject(TObject *o){
2219   //
2220   // Do a safe deletion
2221   //
2222   if(fTemporaryObjects){
2223     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2224     delete o;
2225   } else{
2226     // just delete
2227     delete o;
2228   }
2229 }
2230 //_________________________________________________________________________
2231 TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
2232   AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
2233   return data;
2234 }
2235 //_________________________________________________________________________
2236 TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
2237   //
2238   // Create efficiency grid and calculate efficiency
2239   // of step to step0
2240   //
2241   TString name("eff");
2242   name += step;
2243   name+= step0;
2244   AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2245   eff->CalculateEfficiency(step,step0);
2246   return eff;
2247 }
2248
2249 //____________________________________________________________________________
2250 THnSparse* AliHFEspectrum::GetCharmWeights(){
2251
2252   //
2253   // Measured D->e based weighting factors
2254   //
2255
2256   const Int_t nDim=1;
2257   Int_t nBin[nDim] = {35};
2258
2259   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.};
2260
2261   fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
2262   for(Int_t idim = 0; idim < nDim; idim++)
2263      fWeightCharm->SetBinEdges(idim, ptbinning1);
2264      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};
2265 //{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};
2266   //points
2267   Double_t pt[1];
2268   for(int i=0; i<nBin[0]; i++){
2269     pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2270     fWeightCharm->Fill(pt,weight[i]);
2271   }
2272   Int_t* ibins[nDim];
2273   for(Int_t ivar = 0; ivar < nDim; ivar++)
2274     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2275   fWeightCharm->SetBinError(ibins[0],0);
2276
2277   return fWeightCharm;
2278 }
2279
2280 //____________________________________________________________________________
2281 THnSparse* AliHFEspectrum::GetBeautyIPEff(){
2282   //
2283   // Return beauty electron IP cut efficiency
2284   //
2285
2286   const Int_t nDim=1;
2287   Int_t nBin[nDim] = {35};
2288
2289   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.};
2290
2291   THnSparseF *ipcut = new THnSparseF("beff", "b eff; pt[GeV/c]", nDim, nBin);
2292   for(Int_t idim = 0; idim < nDim; idim++)
2293      ipcut->SetBinEdges(idim, ptbinning1);
2294   Double_t pt[1];
2295   Double_t weight;
2296   for(int i=0; i<nBin[0]; i++){
2297     pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2298     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
2299     //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
2300     ipcut->Fill(pt,weight);
2301   }
2302   Int_t* ibins[nDim];
2303   for(Int_t ivar = 0; ivar < nDim; ivar++)
2304     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2305   ipcut->SetBinError(ibins[0],0);
2306
2307   return ipcut;
2308 }
2309
2310 //____________________________________________________________________________
2311 THnSparse* AliHFEspectrum::GetCharmEff(){
2312   //
2313   // Return charm electron IP cut efficiency
2314   //
2315
2316   const Int_t nDim=1;
2317   Int_t nBin[nDim] = {35};
2318
2319   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.};
2320
2321   THnSparseF *ipcut = new THnSparseF("ceff", "c eff; pt[GeV/c]", nDim, nBin);
2322   for(Int_t idim = 0; idim < nDim; idim++)
2323      ipcut->SetBinEdges(idim, ptbinning1);
2324   Double_t pt[1];
2325   Double_t weight;
2326   for(int i=0; i<nBin[0]; i++){
2327     pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2328     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
2329     //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
2330     ipcut->Fill(pt,weight);
2331   }
2332   Int_t* ibins[nDim];
2333   for(Int_t ivar = 0; ivar < nDim; ivar++)
2334     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2335   ipcut->SetBinError(ibins[0],0);
2336
2337   return ipcut;
2338 }
2339
2340 //____________________________________________________________________________
2341 THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
2342   //
2343   // Return PID x IP cut efficiency
2344   //
2345
2346   const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
2347
2348   const Int_t nDim=1;//dimensions of the efficiency weighting grid
2349   Int_t nBin[nDim] = {nPtbinning1};
2350
2351   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
2352
2353   THnSparseF *pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
2354   for(Int_t idim = 0; idim < nDim; idim++)
2355      pideff->SetBinEdges(idim, kPtRange);
2356
2357   Double_t pt[1];
2358   Double_t weight = 1.0;
2359   Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
2360   for(int i=0; i<nBin[0]; i++){
2361     pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
2362     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
2363
2364     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
2365     //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
2366     //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);
2367     //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
2368     if(source==2) weight = tofeff*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // multiply ip cut eff for conversion electrons
2369     if(source==3) weight = tofeff*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // multiply ip cut eff for Dalitz/dielectrons
2370     //printf("tofeff= %lf trdtpcPidEfficiency= %lf conversion= %lf nonhfe= %lf\n",tofeff,trdtpcPidEfficiency,fConversionEff->GetBinContent(i+1),fNonHFEEff->GetBinContent(i+1));
2371
2372     pideff->Fill(pt,weight);
2373   }
2374   Int_t* ibins[nDim];
2375   for(Int_t ivar = 0; ivar < nDim; ivar++)
2376     ibins[ivar] = new Int_t[nBin[ivar] + 1];
2377   pideff->SetBinError(ibins[0],0);
2378
2379   return pideff;
2380 }
2381 //__________________________________________________________________________
2382 void AliHFEspectrum::CalculateNonHFEsyst(){
2383   //
2384   //Calculate systematic errors for non-HFE and conversion background,
2385   //based on correlated errors from pi0 input and uncorrelated errors
2386   //from m_t scaling
2387   //
2388   if(!fNonHFEsyst)
2389     return;
2390
2391   Double_t evtnorm[1] = {0.0};
2392   if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
2393
2394   AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
2395   AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
2396
2397   AliCFDataGrid *bgLevelGrid[kBgLevels];
2398   AliCFDataGrid *bgNonHFEGrid[kBgLevels];
2399   AliCFDataGrid *bgConvGrid[kBgLevels];
2400
2401   Int_t stepbackground = 3;
2402   Int_t* bins=new Int_t[1];
2403   bins[0]=fConversionEff->GetNbinsX();
2404   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
2405   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
2406
2407   for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
2408
2409     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
2410       convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
2411       weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
2412       convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
2413
2414       nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
2415       weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
2416       nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
2417     }
2418
2419     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
2420     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
2421       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
2422     }
2423
2424     bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
2425     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){//add other sources to get overall background from all meson decays
2426       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
2427     }
2428
2429     bgLevelGrid[iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
2430     bgLevelGrid[iLevel]->Add(bgNonHFEGrid[iLevel]);
2431   }
2432
2433   //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)
2434   AliCFDataGrid *bgErrorGrid[2];
2435   bgErrorGrid[0] = (AliCFDataGrid*)bgLevelGrid[1]->Clone();
2436   bgErrorGrid[1] = (AliCFDataGrid*)bgLevelGrid[2]->Clone();
2437   bgErrorGrid[0]->Add(bgLevelGrid[0],-1.);
2438   bgErrorGrid[1]->Add(bgLevelGrid[0],-1.);
2439
2440   //plot absolute differences between limit yields (upper/lower limit, based on pi0 errors) and best estimate
2441   TH1D* hpiErrors[2];
2442   hpiErrors[0] = (TH1D*)bgErrorGrid[0]->Project(0);
2443   hpiErrors[0]->Scale(-1.);
2444   hpiErrors[0]->SetTitle("Absolute systematic errors from non-HF meson decays and conversions");
2445   hpiErrors[1] = (TH1D*)bgErrorGrid[1]->Project(0);
2446
2447
2448
2449    //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
2450   TH1D *hSpeciesErrors[kElecBgSources-1];
2451   for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
2452     hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
2453     TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
2454     hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
2455     hSpeciesErrors[iSource-1]->Scale(0.3);
2456   }
2457
2458   TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[0]->Clone();
2459   TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[0]->Clone();
2460   TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[0]->Clone();
2461
2462   TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
2463   TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
2464
2465   for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
2466     Double_t pi0basedErrLow = hpiErrors[0]->GetBinContent(iBin);
2467     Double_t pi0basedErrUp = hpiErrors[1]->GetBinContent(iBin);
2468
2469     Double_t sqrsumErrs= 0;
2470     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
2471       Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
2472       sqrsumErrs+=(scalingErr*scalingErr);
2473     }
2474     for(Int_t iErr = 0; iErr < 2; iErr++){
2475       hpiErrors[iErr]->SetBinContent(iBin,hpiErrors[iErr]->GetBinContent(iBin)/hpiErrors[iErr]->GetBinWidth(iBin));
2476     }
2477     hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+sqrsumErrs));
2478     hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+sqrsumErrs));
2479     hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
2480
2481     hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
2482     hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
2483   }
2484
2485
2486    // /hOverallSystErrUp->GetBinWidth(iBin))
2487
2488   TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
2489   cPiErrors->cd();
2490   cPiErrors->SetLogx();
2491   cPiErrors->SetLogy();
2492   hpiErrors[0]->Draw();
2493   hpiErrors[1]->SetMarkerColor(kBlack);
2494   hpiErrors[1]->SetLineColor(kBlack);
2495   hpiErrors[1]->Draw("SAME");
2496   hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
2497   hOverallBinScaledErrsUp->SetLineColor(kBlue);
2498   hOverallBinScaledErrsUp->Draw("SAME");
2499   hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
2500   hOverallBinScaledErrsLow->SetLineColor(kGreen);
2501   hOverallBinScaledErrsLow->Draw("SAME");
2502   hScalingErrors->SetLineColor(11);
2503   hScalingErrors->Draw("SAME");
2504
2505   TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
2506   lPiErr->AddEntry(hpiErrors[0],"Lower error from pion error");
2507   lPiErr->AddEntry(hpiErrors[1],"Upper error from pion error");
2508   lPiErr->AddEntry(hScalingErrors, "scaling error");
2509   lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
2510   lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
2511   lPiErr->Draw("SAME");
2512
2513   //Normalize errors
2514   TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
2515   TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
2516   hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
2517   hLowSystScaled->Scale(evtnorm[0]);
2518   TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
2519   TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
2520   //histograms to be normalized to TGraphErrors
2521   CorrectFromTheWidth(hNormAllSystErrUp);
2522   CorrectFromTheWidth(hNormAllSystErrLow);
2523
2524   TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
2525   cNormOvErrs->cd();
2526   cNormOvErrs->SetLogx();
2527   cNormOvErrs->SetLogy();
2528
2529   TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
2530   TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
2531   gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
2532   gOverallSystErrUp->SetMarkerColor(kBlack);
2533   gOverallSystErrUp->SetLineColor(kBlack);
2534   gOverallSystErrLow->SetMarkerColor(kRed);
2535   gOverallSystErrLow->SetLineColor(kRed);
2536   gOverallSystErrUp->Draw("AP");
2537   gOverallSystErrLow->Draw("PSAME");
2538   TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
2539   lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
2540   lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
2541   lAllSys->Draw("same");
2542
2543
2544   AliCFDataGrid *bgYieldGrid;
2545   bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0]->Clone();
2546
2547   TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
2548   TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
2549   hRelErrUp->Divide(hBgYield);
2550   TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
2551   hRelErrLow->Divide(hBgYield);
2552
2553   TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
2554   cRelErrs->cd();
2555   cRelErrs->SetLogx();
2556   hRelErrUp->SetTitle("Relative error of non-HFE background yield");
2557   hRelErrUp->Draw();
2558   hRelErrLow->SetLineColor(kBlack);
2559   hRelErrLow->Draw("SAME");
2560
2561   TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
2562   lRel->AddEntry(hRelErrUp, "upper");
2563   lRel->AddEntry(hRelErrLow, "lower");
2564   lRel->Draw("SAME");
2565
2566   //CorrectFromTheWidth(hBgYield);
2567   //hBgYield->Scale(evtnorm[0]);
2568
2569
2570   //write histograms/TGraphs to file
2571   TFile *output = new TFile("systHists.root","recreate");
2572
2573   hBgYield->SetName("hBgYield");
2574   hBgYield->Write();
2575   hRelErrUp->SetName("hRelErrUp");
2576   hRelErrUp->Write();
2577   hRelErrLow->SetName("hRelErrLow");
2578   hRelErrLow->Write();
2579   hUpSystScaled->SetName("hOverallSystErrUp");
2580   hUpSystScaled->Write();
2581   hLowSystScaled->SetName("hOverallSystErrLow");
2582   hLowSystScaled->Write();
2583   gOverallSystErrUp->SetName("gOverallSystErrUp");
2584   gOverallSystErrUp->Write();
2585   gOverallSystErrLow->SetName("gOverallSystErrLow");
2586   gOverallSystErrLow->Write();
2587
2588   output->Close();
2589   delete output;
2590 }