PWGHF/hfe/AliHFEspectrum.cxx

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