PWGHF/hfe/AliHFEBeautySpectrum.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 #include <TArrayD.h>
  31 #include <TH1.h>
  32 #include <TList.h>
  33 #include <TObjArray.h>
  34 #include <TROOT.h>
  35 #include <TCanvas.h>
  36 #include <TLegend.h>
  37 #include <TStyle.h>
  38 #include <TMath.h>
  39 #include <TAxis.h>
  40 #include <TGraphErrors.h>
  41 #include <TFile.h>
  42 #include <TPad.h>
  43 #include <TH2D.h>
  44 #include <TF1.h>
  45
  46 #include "AliPID.h"
  47 #include "AliCFContainer.h"
  48 #include "AliCFDataGrid.h"
  49 #include "AliCFEffGrid.h"
  50 #include "AliCFGridSparse.h"
  51 #include "AliCFUnfolding.h"
  52 #include "AliLog.h"
  53
  54 #include "AliHFEBeautySpectrum.h"
  55 #include "AliHFEBeautySpectrumQA.h"
  56 #include "AliHFEcuts.h"
  57 #include "AliHFEcontainer.h"
  58 #include "AliHFEtools.h"
  59
  60 ClassImp(AliHFEBeautySpectrum)
  61
  62 //____________________________________________________________
  63 AliHFEBeautySpectrum::AliHFEBeautySpectrum(const char *name):
  64 AliHFECorrectSpectrumBase(name),
  65   fQA(NULL),
  66   fTemporaryObjects(NULL),
  67   fBackground(NULL),
  68   fEfficiencyTOFPIDD(NULL),
  69   fEfficiencyesdTOFPIDD(NULL),
  70   fEfficiencyIPCharmD(NULL),
  71   fEfficiencyIPBeautyD(NULL),
  72   fEfficiencyIPBeautyesdD(NULL),
  73   fEfficiencyIPConversionD(NULL),
  74   fEfficiencyIPNonhfeD(NULL),
  75   fNonHFEbg(NULL),
  76   fWeightCharm(NULL),
  77   fInclusiveSpectrum(kFALSE),
  78   fDumpToFile(kFALSE),
  79   fUnSetCorrelatedErrors(kTRUE),
  80   fIPanaHadronBgSubtract(kFALSE),
  81   fIPanaCharmBgSubtract(kFALSE),
  82   fIPanaNonHFEBgSubtract(kFALSE),
  83   fIPParameterizedEff(kFALSE),
  84   fNonHFEsyst(0),
  85   fBeauty2ndMethod(kFALSE),
  86   fIPEffCombinedSamples(kTRUE),
  87   fNMCEvents(0),
  88   fNMCbgEvents(0),
  89   fNCentralityBinAtTheEnd(0),
  90   fFillMoreCorrelationMatrix(kFALSE),
  91   fHadronEffbyIPcut(NULL),
  92   fEfficiencyCharmSigD(NULL),
  93   fEfficiencyBeautySigD(NULL),
  94   fEfficiencyBeautySigesdD(NULL),
  95   fConversionEff(NULL),
  96   fNonHFEEff(NULL),
  97   fCharmEff(NULL),
  98   fBeautyEff(NULL),
  99   fConversionEffbgc(NULL),
 100   fNonHFEEffbgc(NULL),
 101   fBSpectrum2ndMethod(NULL),
 102   fkBeauty2ndMethodfilename(""),
 103   fBeamType(0),
 104   fEtaSyst(kTRUE),
 105   fDebugLevel(0),
 106   fWriteToFile(kFALSE),
 107   fUnfoldBG(kFALSE)
 108 {
 109   //
 110   // Default constructor
 111   //
 112
 113   fQA = new AliHFEBeautySpectrumQA("AliHFEBeautySpectrumQA");
 114
 115   for(Int_t k = 0; k < 20; k++){
 116     fLowBoundaryCentralityBinAtTheEnd[k] = 0;
 117     fHighBoundaryCentralityBinAtTheEnd[k] = 0;
 118   }
 119   fNMCbgEvents = 0;
 120   fEfficiencyTOFPIDD = 0;
 121   fEfficiencyesdTOFPIDD = 0;
 122   fEfficiencyIPCharmD = 0;
 123   fEfficiencyIPBeautyD = 0;
 124   fEfficiencyIPBeautyesdD = 0;
 125   fEfficiencyIPConversionD = 0;
 126   fEfficiencyIPNonhfeD = 0;
 127
 128   fConversionEff = 0;
 129   fNonHFEEff = 0;
 130   fCharmEff = 0;
 131   fBeautyEff = 0;
 132   fEfficiencyCharmSigD = 0;
 133   fEfficiencyBeautySigD = 0;
 134   fEfficiencyBeautySigesdD = 0;
 135
 136
 137   memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
 138   memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
 139 }
 140 //____________________________________________________________
 141 AliHFEBeautySpectrum::AliHFEBeautySpectrum(const AliHFEBeautySpectrum &ref):
 142   AliHFECorrectSpectrumBase(ref),
 143   fQA(ref.fQA),
 144   fTemporaryObjects(ref.fTemporaryObjects),
 145   fBackground(ref.fBackground),
 146   fEfficiencyTOFPIDD(ref.fEfficiencyTOFPIDD),
 147   fEfficiencyesdTOFPIDD(ref.fEfficiencyesdTOFPIDD),
 148   fEfficiencyIPCharmD(ref.fEfficiencyIPCharmD),
 149   fEfficiencyIPBeautyD(ref.fEfficiencyIPBeautyD),
 150   fEfficiencyIPBeautyesdD(ref.fEfficiencyIPBeautyesdD),
 151   fEfficiencyIPConversionD(ref.fEfficiencyIPConversionD),
 152   fEfficiencyIPNonhfeD(ref.fEfficiencyIPNonhfeD),
 153   fNonHFEbg(ref.fNonHFEbg),
 154   fWeightCharm(ref.fWeightCharm),
 155   fInclusiveSpectrum(ref.fInclusiveSpectrum),
 156   fDumpToFile(ref.fDumpToFile),
 157   fUnSetCorrelatedErrors(ref.fUnSetCorrelatedErrors),
 158   fIPanaHadronBgSubtract(ref.fIPanaHadronBgSubtract),
 159   fIPanaCharmBgSubtract(ref.fIPanaCharmBgSubtract),
 160   fIPanaNonHFEBgSubtract(ref.fIPanaNonHFEBgSubtract),
 161   fIPParameterizedEff(ref.fIPParameterizedEff),
 162   fNonHFEsyst(ref.fNonHFEsyst),
 163   fBeauty2ndMethod(ref.fBeauty2ndMethod),
 164   fIPEffCombinedSamples(ref.fIPEffCombinedSamples),
 165   fNMCEvents(ref.fNMCEvents),
 166   fNMCbgEvents(ref.fNMCbgEvents),
 167   fNCentralityBinAtTheEnd(ref.fNCentralityBinAtTheEnd),
 168   fFillMoreCorrelationMatrix(ref.fFillMoreCorrelationMatrix),
 169   fHadronEffbyIPcut(ref.fHadronEffbyIPcut),
 170   fEfficiencyCharmSigD(ref.fEfficiencyCharmSigD),
 171   fEfficiencyBeautySigD(ref.fEfficiencyBeautySigD),
 172   fEfficiencyBeautySigesdD(ref.fEfficiencyBeautySigesdD),
 173   fConversionEff(ref.fConversionEff),
 174   fNonHFEEff(ref.fNonHFEEff),
 175   fCharmEff(ref.fCharmEff),
 176   fBeautyEff(ref.fBeautyEff),
 177   fConversionEffbgc(ref.fConversionEffbgc),
 178   fNonHFEEffbgc(ref.fNonHFEEffbgc),
 179   fBSpectrum2ndMethod(ref.fBSpectrum2ndMethod),
 180   fkBeauty2ndMethodfilename(ref.fkBeauty2ndMethodfilename),
 181   fBeamType(ref.fBeamType),
 182   fEtaSyst(ref.fEtaSyst),
 183   fDebugLevel(ref.fDebugLevel),
 184   fWriteToFile(ref.fWriteToFile),
 185   fUnfoldBG(ref.fUnfoldBG)
 186 {
 187   //
 188   // Copy constructor
 189   //
 190
 191   ref.Copy(*this);
 192 }
 193 //____________________________________________________________
 194 AliHFEBeautySpectrum &AliHFEBeautySpectrum::operator=(const AliHFEBeautySpectrum &ref){
 195  //
 196   // Assignment operator
 197   //
 198   if(this == &ref)
 199     ref.Copy(*this);
 200   return *this;
 201 }
 202 //____________________________________________________________
 203 void AliHFEBeautySpectrum::Copy(TObject &o) const {
 204   //
 205   // Copy into object o
 206   //
 207   AliHFEBeautySpectrum &target = dynamic_cast<AliHFEBeautySpectrum &>(o);
 208   target.fQA = fQA;
 209
 210   target.fQA = fQA;
 211   target.fTemporaryObjects = fTemporaryObjects;
 212   target.fBackground = fBackground;
 213   target.fWeightCharm = fWeightCharm;
 214   target.fDumpToFile = fDumpToFile;
 215   target.fUnSetCorrelatedErrors = fUnSetCorrelatedErrors;
 216   target.fIPanaHadronBgSubtract = fIPanaHadronBgSubtract;
 217   target.fIPanaCharmBgSubtract = fIPanaCharmBgSubtract;
 218   target.fIPanaNonHFEBgSubtract = fIPanaNonHFEBgSubtract;
 219   target.fIPParameterizedEff = fIPParameterizedEff;
 220   target.fNonHFEsyst = fNonHFEsyst;
 221   target.fBeauty2ndMethod = fBeauty2ndMethod;
 222   target.fIPEffCombinedSamples = fIPEffCombinedSamples;
 223   target.fNMCEvents = fNMCEvents;
 224   target.fNCentralityBinAtTheEnd = fNCentralityBinAtTheEnd;
 225   target.fFillMoreCorrelationMatrix = fFillMoreCorrelationMatrix;
 226   target.fHadronEffbyIPcut = fHadronEffbyIPcut;
 227   target.fConversionEffbgc = fConversionEffbgc;
 228   target.fNonHFEEffbgc = fNonHFEEffbgc;
 229   target.fBSpectrum2ndMethod = fBSpectrum2ndMethod;
 230   target.fkBeauty2ndMethodfilename = fkBeauty2ndMethodfilename;
 231   target.fBeamType = fBeamType;
 232   target.fEtaSyst = fEtaSyst;
 233   target.fDebugLevel = fDebugLevel;
 234   target.fWriteToFile = fWriteToFile;
 235   target.fUnfoldBG = fUnfoldBG;
 236 }
 237 //____________________________________________________________
 238 AliHFEBeautySpectrum::~AliHFEBeautySpectrum(){
 239   //
 240   // Destructor
 241   //
 242   if(fTemporaryObjects){
 243     fTemporaryObjects->Clear();
 244     delete fTemporaryObjects;
 245   }
 246   if(fQA) delete fQA;
 247 }
 248 //____________________________________________________________
 249 Bool_t AliHFEBeautySpectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *bghfecontainer, const AliHFEcontainer */*v0hfecontainer*/,AliCFContainer */*photoniccontainerD*/){
 250   //
 251   // Init what we need for the correction:
 252   //
 253   // Raw spectrum, hadron contamination
 254   // MC efficiency maps, correlation matrix
 255   //
 256   // This for a given dimension.
 257   // If no inclusive spectrum, then take only efficiency map for beauty electron
 258   // and the appropriate correlation matrix
 259   //
 260
 261
 262   if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
 263
 264   Int_t kNdim = 3;
 265
 266   Int_t dims[kNdim];
 267
 268   switch(fNbDimensions){
 269   case 1:   dims[0] = 0;
 270     break;
 271   case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
 272     break;
 273   case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
 274     break;
 275   default:
 276     AliError("Container with this number of dimensions not foreseen (yet)");
 277     return kFALSE;
 278   };
 279
 280   //
 281   // Data container: raw spectrum + hadron contamination
 282   //
 283   AliCFContainer *datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
 284   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
 285   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
 286   AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 287   AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 288   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
 289   SetContainer(datacontainerD,AliHFECorrectSpectrumBase::kDataContainer);
 290   SetContainer(contaminationcontainerD,AliHFECorrectSpectrumBase::kBackgroundData);
 291
 292   // QA : see with MJ if it is necessary for Beauty
 293   Int_t dimqa = datacontainer->GetNVar();
 294   Int_t dimsqa[dimqa];
 295   for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
 296   AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 297   fQA->AddResultAt(datacontainerDQA,AliHFEBeautySpectrumQA::kDataProjection);
 298
 299   //
 300   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
 301   //
 302   AliCFContainer *mccontaineresd = 0x0;
 303   AliCFContainer *mccontaineresdbg = 0x0;
 304   AliCFContainer *mccontainermc = 0x0;
 305   AliCFContainer *mccontainermcbg = 0x0;
 306   AliCFContainer *nonHFEweightedContainer = 0x0;
 307   AliCFContainer *convweightedContainer = 0x0;
 308   AliCFContainer *nonHFEweightedContainerSig = 0x0;//mjmj
 309   AliCFContainer *convweightedContainerSig = 0x0;//mjmj
 310   AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
 311   AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
 312
 313   mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
 314   mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
 315   mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
 316   mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
 317
 318   if(fNonHFEsyst){
 319     const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
 320     const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
 321     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
 322       for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
 323         if(!(nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]))))continue;
 324         convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
 325         if(!(fConvSourceContainer[iSource][iLevel] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
 326         if(!(fNonHFESourceContainer[iSource][iLevel] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
 327         //            if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
 328         if(fBeamType == 1)break;//no systematics yet for PbPb non-HFE backgrounds
 329       }
 330     }
 331   }
 332   // else{
 333   nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
 334   convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
 335   if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
 336   nonHFEweightedContainerSig = mchfecontainer->GetCFContainer("mesonElecs");//mjmj
 337   convweightedContainerSig = mchfecontainer->GetCFContainer("conversionElecs");//mjmj
 338   if((!convweightedContainerSig)||(!nonHFEweightedContainerSig)) return kFALSE;
 339   //}
 340
 341   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
 342
 343   Int_t source = 1; //beauty
 344   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 345   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 346   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
 347   SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerMC);
 348   SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
 349
 350   // set charm, nonHFE container to estimate BG
 351   source = 0; //charm
 352   mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); //mjmjmj
 353   //mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 354   SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerCharmMC);
 355
 356   //if(!fNonHFEsyst){
 357   AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 358   SetContainer(nonHFEweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESD);
 359   AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 360   SetContainer(convweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESD);
 361
 362   //mjmj
 363   AliCFContainer *nonHFEweightedContainerSigD = GetSlicedContainer(nonHFEweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 364   SetContainer(nonHFEweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
 365   if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig)){printf("merged non-HFE container not found!\n");return kFALSE;};
 366   AliCFContainer *convweightedContainerSigD = GetSlicedContainer(convweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 367   SetContainer(convweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
 368   if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig)){printf(" merged conversion container not found!\n");return kFALSE;};
 369   //}
 370
 371   SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
 372
 373   //
 374   // MC container: correlation matrix
 375   //
 376   THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
 377   //mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
 378   if(!mccorrelation) return kFALSE;
 379   Int_t testCentralityLow = fTestCentralityLow;
 380   Int_t testCentralityHigh = fTestCentralityHigh;
 381   if(fFillMoreCorrelationMatrix) {
 382     testCentralityLow = fTestCentralityLow-1;
 383     testCentralityHigh = fTestCentralityHigh+1;
 384   }
 385   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims, fChargeChoosen, testCentralityLow,testCentralityHigh);
 386   if(!mccorrelationD) {
 387     printf("No correlation\n");
 388     return kFALSE;
 389   }
 390   SetCorrelation(mccorrelationD);
 391
 392   // QA
 393   fQA->AddResultAt(mccorrelation,AliHFEBeautySpectrumQA::kCMProjection);
 394   //
 395   fQA->DrawProjections();
 396
 397
 398   return kTRUE;
 399 }
 400
 401
 402 //____________________________________________________________
 403 void AliHFEBeautySpectrum::CallInputFileForBeauty2ndMethod(){
 404   //
 405   // get spectrum for beauty 2nd method
 406   //
 407   //
 408     TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
 409     fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
 410 }
 411
 412 //____________________________________________________________
 413 Bool_t AliHFEBeautySpectrum::Correct(Bool_t subtractcontamination, Bool_t /*subtractphotonic*/){
 414   //
 415   // Correct the spectrum for efficiency and unfolding for beauty analysis
 416   // with both method and compare
 417   //
 418
 419   gStyle->SetPalette(1);
 420   gStyle->SetOptStat(1111);
 421   gStyle->SetPadBorderMode(0);
 422   gStyle->SetCanvasColor(10);
 423   gStyle->SetPadLeftMargin(0.13);
 424   gStyle->SetPadRightMargin(0.13);
 425
 426   printf("Steps are: stepdata %d, stepMC %d, steptrue %d\n",fStepData,fStepMC,fStepTrue);
 427
 428   ///////////////////////////
 429   // Check initialization
 430   ///////////////////////////
 431
 432   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
 433     AliInfo("You have to init before");
 434     return kFALSE;
 435   }
 436
 437   if((fStepTrue <= 0) && (fStepMC <= 0) && (fStepData <= 0)) {
 438     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
 439     return kFALSE;
 440   }
 441
 442   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
 443
 444   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
 445   //////////////////////////////////
 446   // Subtract hadron background
 447   /////////////////////////////////
 448   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
 449   AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
 450   TGraphErrors *gNormalizedRawSpectrum = 0x0;
 451   if(subtractcontamination) {
 452       if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
 453       else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
 454       unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
 455       dataGridAfterFirstSteps = dataspectrumaftersubstraction;
 456       gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
 457   }
 458
 459   /////////////////////////////////////////////////////////////////////////////////////////
 460   // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
 461   /////////////////////////////////////////////////////////////////////////////////////////
 462
 463   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
 464
 465   if(fEfficiencyFunction){
 466     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
 467     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
 468   }
 469
 470   ///////////////
 471   // Unfold
 472   //////////////
 473   THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
 474   if(!correctedspectrum){
 475     AliError("No corrected spectrum\n");
 476     return kFALSE;
 477   }
 478
 479   /////////////////////
 480   // Simply correct
 481   ////////////////////
 482
 483   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
 484
 485
 486   //////////
 487   // Plot
 488   //////////
 489
 490   // QA final results
 491   TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
 492   TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
 493   fQA->AddResultAt(correctedspectrumD,AliHFEBeautySpectrumQA::kFinalResultUnfolded);
 494   fQA->AddResultAt(alltogetherspectrumD,AliHFEBeautySpectrumQA::kFinalResultDirectEfficiency);
 495   fQA->AddResultAt(correctedspectrum,AliHFEBeautySpectrumQA::kFinalResultUnfSparse);
 496   fQA->AddResultAt(alltogetherCorrection,AliHFEBeautySpectrumQA::kFinalResultDirectEffSparse);
 497   fQA->DrawResult();
 498
 499
 500   if(fBeamType == 0){
 501     if(fNonHFEsyst){
 502       CalculateNonHFEsyst();
 503     }
 504   }
 505
 506   // Dump to file if needed
 507
 508   if(fDumpToFile) {
 509     TFile *out;
 510     out = new TFile("finalSpectrum.root","update");
 511     out->cd();
 512     // to do centrality dependent
 513     correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
 514     correctedspectrum->Write();
 515     alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
 516     alltogetherCorrection->Write();
 517     //
 518     if(unnormalizedRawSpectrum) {
 519       unnormalizedRawSpectrum->SetName("beautyAfterIP");
 520       unnormalizedRawSpectrum->Write();
 521     }
 522
 523     if(gNormalizedRawSpectrum){
 524       gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
 525       gNormalizedRawSpectrum->Write();
 526     }
 527
 528     fEfficiencyCharmSigD->SetTitle("IPEfficiencyForCharmSig");
 529     fEfficiencyCharmSigD->SetName("IPEfficiencyForCharmSig");
 530     fEfficiencyCharmSigD->Write();
 531     fEfficiencyBeautySigD->SetTitle("IPEfficiencyForBeautySig");
 532     fEfficiencyBeautySigD->SetName("IPEfficiencyForBeautySig");
 533     fEfficiencyBeautySigD->Write();
 534     fCharmEff->SetTitle("IPEfficiencyForCharm");
 535     fCharmEff->SetName("IPEfficiencyForCharm");
 536     fCharmEff->Write();
 537     fBeautyEff->SetTitle("IPEfficiencyForBeauty");
 538     fBeautyEff->SetName("IPEfficiencyForBeauty");
 539     fBeautyEff->Write();
 540     fConversionEff->SetTitle("IPEfficiencyForConversion");
 541     fConversionEff->SetName("IPEfficiencyForConversion");
 542     fConversionEff->Write();
 543     fNonHFEEff->SetTitle("IPEfficiencyForNonHFE");
 544     fNonHFEEff->SetName("IPEfficiencyForNonHFE");
 545     fNonHFEEff->Write();
 546
 547
 548     out->Close();
 549     delete out;
 550    }
 551   return kTRUE;
 552 }
 553 //____________________________________________________________
 554 AliCFDataGrid* AliHFEBeautySpectrum::SubtractBackground(Bool_t setBackground){
 555   //
 556   // Apply background subtraction for IP analysis
 557   //
 558
 559   Int_t nbins=1;
 560   printf("subtracting background\n");
 561   // Raw spectrum
 562   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 563   if(!dataContainer){
 564     AliError("Data Container not available");
 565     return NULL;
 566   }
 567   printf("Step data: %d\n",fStepData);
 568   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
 569
 570   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 571   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
 572
 573   // Background Estimate
 574   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
 575   if(!backgroundContainer){
 576     AliError("MC background container not found");
 577     return NULL;
 578   }
 579
 580   Int_t stepbackground = 1;
 581   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
 582
 583   const Int_t bgPlots = 8;
 584   TH1D *incElec = 0x0;
 585   TH1D *hadron = 0x0;
 586   TH1D *charm = 0x0;
 587   TH1D *nonHFE[bgPlots];
 588   TH1D *subtracted = 0x0;
 589
 590   THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
 591   incElec = (TH1D *) sparseIncElec->Projection(0);
 592   AliHFEtools::NormaliseBinWidth(incElec);
 593
 594   TH1D* htemp;
 595   Int_t* bins=new Int_t[2];
 596
 597   if(fIPanaHadronBgSubtract){
 598     // Hadron background
 599     printf("Hadron background for IP analysis subtracted!\n");
 600     htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
 601     bins[0]=htemp->GetNbinsX();
 602
 603     AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
 604     hbgContainer->SetGrid(fHadronEffbyIPcut);
 605     backgroundGrid->Multiply(hbgContainer,1);
 606     // subtract hadron contamination
 607     spectrumSubtracted->Add(backgroundGrid,-1.0);
 608     // draw raw hadron bg spectra
 609     THnSparseF* sparseHbg = (THnSparseF *) hbgContainer->GetGrid();
 610     hadron = (TH1D *) sparseHbg->Projection(0);
 611     AliHFEtools::NormaliseBinWidth(hadron);
 612   }
 613
 614   if(fIPanaCharmBgSubtract){
 615     // Charm background
 616     printf("Charm background for IP analysis subtracted!\n");
 617     AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
 618     spectrumSubtracted->Add(charmbgContainer,-1.0);
 619     THnSparseF *sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
 620     charm = (TH1D *) sparseCharmElec->Projection(0);
 621     AliHFEtools::NormaliseBinWidth(charm);
 622   }
 623
 624   const Char_t *sourceName[bgPlots]={"Conversion","Dalitz","K0s secondary pion Dalitz","K0s secondary pion conversions","Other secondary pion Dalitz","Other secondary pion Dalitz conversions","Other Dalitz", "Other conversions"};
 625   const Int_t style[2] = {21,22};
 626   const Int_t color[3] = {7,12,8};
 627   if(fIPanaNonHFEBgSubtract){
 628     Int_t iPlot = 0;
 629     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
 630       if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
 631       for(Int_t decay = 0; decay < 2; decay++){
 632         AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground(decay,iSource);//conversions/Dalitz decays from different sources
 633         if(iSource == 0)
 634           spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
 635         THnSparseF* sparseNonHFEelecs = (THnSparseF *) nonHFEbgContainer->GetGrid();
 636         nonHFE[iPlot] = (TH1D *) sparseNonHFEelecs->Projection(0);
 637         AliHFEtools::NormaliseBinWidth(nonHFE[iPlot]);
 638         iPlot++;
 639       }
 640       if(!(fNonHFEsyst && (fBeamType == 1)))break;
 641     }
 642   }
 643
 644   TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
 645   TCanvas *cRaw = new TCanvas("cRaw","cRaw",1000,800);
 646   cRaw->cd();
 647   gPad->SetLogx();
 648   gPad->SetLogy();
 649   incElec->GetXaxis()->SetRangeUser(0.4,8.);
 650   incElec->Draw("p");
 651   incElec->SetMarkerColor(1);
 652   incElec->SetMarkerStyle(20);
 653   lRaw->AddEntry(incElec,"inclusive electron spectrum");
 654
 655   if(fIPanaCharmBgSubtract){
 656     charm->Draw("samep");
 657     charm->SetMarkerColor(3);
 658     charm->SetMarkerStyle(20);
 659     charm->GetXaxis()->SetRangeUser(0.4,8.);
 660     lRaw->AddEntry(charm,"charm elecs");
 661   }
 662
 663   if(fIPanaNonHFEBgSubtract){
 664     Int_t iPlot = 0;
 665     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
 666       if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
 667       for(Int_t decay = 0; decay < 2; decay++){
 668         nonHFE[iPlot]->Draw("samep");
 669         if(iSource == 0){
 670           nonHFE[iPlot]->SetMarkerStyle(20);
 671           if(decay == 0)
 672             nonHFE[iPlot]->SetMarkerColor(4);
 673           else
 674             nonHFE[iPlot]->SetMarkerColor(6);
 675         }
 676         else{
 677           nonHFE[iPlot]->SetMarkerColor(color[iSource-6]);
 678           nonHFE[iPlot]->SetMarkerStyle(style[decay]);
 679         }
 680         nonHFE[iPlot]->GetXaxis()->SetRangeUser(0.4,8.);
 681         lRaw->AddEntry(nonHFE[iPlot],sourceName[iPlot]);
 682         iPlot++;
 683       }
 684       if(!(fNonHFEsyst && (fBeamType == 1)))break;//only draw standard sources
 685     }
 686   }
 687
 688   THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
 689   subtracted = (TH1D *) sparseSubtracted->Projection(0);
 690   AliHFEtools::NormaliseBinWidth(subtracted);
 691   subtracted->Draw("samep");
 692   subtracted->SetMarkerStyle(24);
 693   lRaw->AddEntry(subtracted,"subtracted electron spectrum");
 694   lRaw->Draw("SAME");
 695
 696   delete[] bins;
 697
 698   TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(0);
 699   AliHFEtools::NormaliseBinWidth(measuredTH1background);
 700
 701   if(setBackground){
 702     if(fBackground) delete fBackground;
 703     fBackground = backgroundGrid;
 704   } else delete backgroundGrid;
 705
 706   fQA->AddResultAt(subtracted,AliHFEBeautySpectrumQA::kAfterSC);
 707   fQA->AddResultAt(incElec,AliHFEBeautySpectrumQA::kBeforeSC);
 708   fQA->AddResultAt(measuredTH1background, AliHFEBeautySpectrumQA::kMeasBG);
 709   fQA->DrawSubtractContamination();
 710
 711   return spectrumSubtracted;
 712 }
 713 //____________________________________________________________________________________________
 714 AliCFDataGrid* AliHFEBeautySpectrum::GetNonHFEBackground(Int_t decay, Int_t source){
 715   //
 716   // calculate non-HF electron background
 717   // Arguments: decay = 0 for conversions, decay = 1 for meson decays to electrons
 718   // source: gives mother either of the electron (for meson->electron) or of the gamma (conversions);
 719   // source numbers as given by array in AliAnalysisTaskHFE:
 720   // {"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
 721   //
 722
 723   Double_t evtnorm[1] = {1.};
 724   if(fNMCbgEvents) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
 725
 726   Int_t stepbackground = 1;//for non-standard background, step after IP cut is taken directly
 727   if(source == 0) stepbackground = 3;//for standard sources, take step before PID -> correction of PID x IP efficiencies from enhanced samples
 728   AliCFContainer *backgroundContainer = 0x0;
 729   if(fNonHFEsyst){
 730     if(decay == 0){//conversions
 731       backgroundContainer = (AliCFContainer*)fConvSourceContainer[source][0]->Clone();
 732       if(source == 0)//standard conversion containers
 733         for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){//add all standard sources
 734           backgroundContainer->Add(fConvSourceContainer[iSource][0]);
 735         }
 736     }
 737     else if(decay == 1){//Dalitz decays, same procedure as above
 738       backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[source][0]->Clone();
 739       if(source == 0)
 740         for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
 741           if(iSource == 1)
 742             backgroundContainer->Add(fNonHFESourceContainer[iSource][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
 743           else
 744             backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);
 745         }
 746     }
 747   }
 748   else{//careful: these containers include the non-standard electron sources (K0s, etc.). If necessary, use SetRange in species axis to fix it?
 749     if(source == 0){
 750       // Background Estimate
 751       if(decay == 0)
 752         backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
 753       else if(decay == 1)
 754         backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
 755     }
 756   }
 757   if(!backgroundContainer){
 758     AliError("MC background container not found");
 759     return NULL;
 760   }
 761   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("backgroundGrid","backgroundGrid",*backgroundContainer,stepbackground);
 762
 763   Double_t rangelow = 1.;
 764   Double_t rangeup = 6.;
 765   if(decay == 1) rangelow = 0.9;
 766
 767
 768   const Char_t *dmode[2]={"Conversions","Dalitz decays"};
 769   TF1 *fitHagedorn = new TF1("fitHagedorn", "[0]/TMath::Power(([1]+x/[2]), [3])", rangelow, rangeup);
 770   fNonHFEbg = 0x0;
 771   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
 772   TAxis *axis = h1->GetXaxis();
 773   AliHFEtools::NormaliseBinWidth(h1);
 774   if(source == 0){
 775     fitHagedorn->SetParameter(0, 0.15);
 776     fitHagedorn->SetParameter(1, 0.09);
 777     fitHagedorn->SetParameter(2, 8.4);
 778     fitHagedorn->SetParameter(3, 6.3);
 779     // TCanvas *ch1conversion = new TCanvas("ch1conversion","ch1conversion",500,400);
 780     // ch1conversion->cd();
 781     //fHagedorn->SetLineColor(2);
 782     h1->Fit(fitHagedorn,"R");
 783     // h1->Draw();
 784     if(!(fNonHFEbg = h1->GetFunction("fitHagedorn")))printf("electron background fit failed for %s\n",dmode[decay]);
 785   }
 786
 787   Int_t *nBinpp=new Int_t[1];
 788   Int_t *binspp=new Int_t[1];
 789   binspp[0]=kSignalPtBins;// number of pt bins
 790
 791   Int_t looppt=binspp[0];
 792
 793   for(Long_t iBin=1; iBin<= looppt;iBin++){
 794     Double_t iipt= h1->GetBinCenter(iBin);
 795     Double_t iiwidth= axis->GetBinWidth(iBin);
 796     nBinpp[0]=iBin;
 797     Double_t fitFactor = backgroundGrid->GetElement(nBinpp);//if no fit available, just take bin-by-bin information
 798     if(fNonHFEbg)fitFactor = fNonHFEbg->Eval(iipt)*iiwidth;
 799     backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
 800     backgroundGrid->SetElement(nBinpp,fitFactor*evtnorm[0]);
 801   }
 802   //end of workaround for statistical errors
 803   if(source == 0){
 804     AliCFDataGrid *weightedBGContainer = 0x0;
 805     weightedBGContainer = new AliCFDataGrid("weightedBGContainer","weightedBGContainer",1,binspp);
 806     if(decay == 0)
 807       weightedBGContainer->SetGrid(GetPIDxIPEff(2));
 808     else if(decay == 1)
 809       weightedBGContainer->SetGrid(GetPIDxIPEff(3));
 810     if(stepbackground == 3){
 811       backgroundGrid->Multiply(weightedBGContainer,1.0);
 812     }
 813   }
 814   delete[] nBinpp;
 815   delete[] binspp;
 816   return backgroundGrid;
 817 }
 818 //____________________________________________________________
 819 AliCFDataGrid* AliHFEBeautySpectrum::GetCharmBackground(){
 820   //
 821   // calculate charm background; when using centrality binning 2: centBin = 0 for 0-20%, centBin = 5 for 40-80%
 822   //
 823
 824   Int_t nDim = 1;
 825
 826   Double_t evtnorm=0;
 827   if(fNMCEvents) evtnorm= double(fNEvents)/double(fNMCEvents); //mjmjmj
 828   //if(fNMCbgEvents) evtnorm= double(fNEvents)/double(fNMCbgEvents);
 829   printf("events for data: %d",fNEvents);
 830   printf("events for MC: %d",fNMCEvents);
 831   printf("normalization factor for charm: %f",evtnorm);
 832
 833   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
 834   if(!mcContainer){
 835     AliError("MC Container not available");
 836     return NULL;
 837   }
 838
 839   if(!fCorrelation){
 840     AliError("No Correlation map available");
 841     return NULL;
 842   }
 843
 844   AliCFDataGrid *charmBackgroundGrid= 0x0;
 845   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
 846   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
 847   Int_t* bins=new Int_t[2];
 848   bins[0]=charmbgaftertofpid->GetNbinsX();
 849
 850   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
 851   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
 852   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
 853
 854   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
 855   Int_t *nBinpp=new Int_t[1];
 856   Int_t* binspp=new Int_t[1];
 857   binspp[0]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
 858
 859   Int_t looppt=binspp[0];
 860
 861   for(Long_t iBin=1; iBin<= looppt;iBin++){
 862     nBinpp[0]=iBin;
 863     charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
 864     charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
 865   }
 866
 867   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
 868
 869   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
 870   weightedCharmContainer->SetGrid(GetCharmWeights(fTestCentralityLow)); // get charm weighting factors
 871   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
 872   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
 873   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
 874
 875   // Efficiency (set efficiency to 1 for only folding)
 876   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 877   efficiencyD->CalculateEfficiency(0,0);
 878
 879   // Folding
 880   if(fBeamType==0)nDim = 1;
 881   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
 882   folding.SetMaxNumberOfIterations(1);
 883   folding.Unfold();
 884
 885   // Results
 886   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
 887   THnSparse* result=(THnSparse*)result1->Clone();
 888   charmBackgroundGrid->SetGrid(result);
 889   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
 890
 891   //Charm background evaluation plots
 892
 893   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
 894   cCharmBgEval->Divide(3,1);
 895
 896   cCharmBgEval->cd(1);
 897
 898   if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
 899
 900   AliHFEtools::NormaliseBinWidth(charmbgaftertofpid);
 901   charmbgaftertofpid->SetMarkerStyle(25);
 902   charmbgaftertofpid->Draw("p");
 903   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
 904   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 905   gPad->SetLogy();
 906
 907   AliHFEtools::NormaliseBinWidth(charmbgafteripcut);
 908   charmbgafteripcut->SetMarkerStyle(24);
 909   charmbgafteripcut->Draw("samep");
 910
 911   AliHFEtools::NormaliseBinWidth(charmbgafterweight);
 912   charmbgafterweight->SetMarkerStyle(24);
 913   charmbgafterweight->SetMarkerColor(4);
 914   charmbgafterweight->Draw("samep");
 915
 916   AliHFEtools::NormaliseBinWidth(charmbgafterfolding);
 917   charmbgafterfolding->SetMarkerStyle(24);
 918   charmbgafterfolding->SetMarkerColor(2);
 919   charmbgafterfolding->Draw("samep");
 920
 921   cCharmBgEval->cd(2);
 922   parametrizedcharmpidipeff->SetMarkerStyle(24);
 923   parametrizedcharmpidipeff->Draw("p");
 924   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 925
 926   cCharmBgEval->cd(3);
 927   charmweightingfc->SetMarkerStyle(24);
 928   charmweightingfc->Draw("p");
 929   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
 930   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 931
 932   cCharmBgEval->cd(1);
 933   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
 934   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
 935   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
 936   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
 937   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
 938   legcharmbg->Draw("same");
 939
 940   cCharmBgEval->cd(2);
 941   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
 942   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
 943   legcharmbg2->Draw("same");
 944
 945   CorrectStatErr(charmBackgroundGrid);
 946   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
 947
 948   delete[] bins;
 949   delete[] nBinpp;
 950   delete[] binspp;
 951
 952   if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
 953
 954   return charmBackgroundGrid;
 955 }
 956 //____________________________________________________________
 957 AliCFDataGrid *AliHFEBeautySpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
 958   //
 959   // Apply TPC pid efficiency correction from parametrisation
 960   //
 961
 962   // Data in the right format
 963   AliCFDataGrid *dataGrid = 0x0;
 964   if(bgsubpectrum) {
 965     dataGrid = bgsubpectrum;
 966   }
 967   else {
 968
 969     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 970     if(!dataContainer){
 971       AliError("Data Container not available");
 972       return NULL;
 973     }
 974     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 975   }
 976   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 977   result->SetName("ParametrizedEfficiencyBefore");
 978   THnSparse *h = result->GetGrid();
 979   Int_t nbdimensions = h->GetNdimensions();
 980   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
 981   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 982   if(!dataContainer){
 983     AliError("Data Container not available");
 984     return NULL;
 985   }
 986   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
 987   dataContainerbis->Add(dataContainerbis,-1.0);
 988
 989   Int_t* coord = new Int_t[nbdimensions];
 990   memset(coord, 0, sizeof(Int_t) * nbdimensions);
 991   Double_t* points = new Double_t[nbdimensions];
 992
 993   ULong64_t nEntries = h->GetNbins();
 994   for (ULong64_t i = 0; i < nEntries; ++i) {
 995     Double_t value = h->GetBinContent(i, coord);
 996     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
 997     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
 998
 999     // Get the bin co-ordinates given an coord
1000     for (Int_t j = 0; j < nbdimensions; ++j)
1001       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1002
1003     if (!fEfficiencyFunction->IsInside(points))
1004          continue;
1005     TF1::RejectPoint(kFALSE);
1006
1007     // Evaulate function at points
1008     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1009     //printf("Value efficiency is %f\n",valueEfficiency);
1010
1011     if(valueEfficiency > 0.0) {
1012       h->SetBinContent(coord,value/valueEfficiency);
1013       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1014     }
1015     Double_t error = h->GetBinError(i);
1016     h->SetBinError(coord,error/valueEfficiency);
1017     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1018   }
1019
1020   delete[] coord;
1021   delete[] points;
1022
1023   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1024
1025 // QA
1026   TH1D *afterE = (TH1D *) resultt->Project(0);
1027   AliHFEtools::NormaliseBinWidth(afterE);
1028   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1029   AliHFEtools::NormaliseBinWidth(beforeE);
1030   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterPE);
1031   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforePE);
1032   fQA->AddResultAt(fEfficiencyFunction,AliHFEBeautySpectrumQA::kPEfficiency);
1033   fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kParametrized);
1034
1035   return resultt;
1036 }
1037 //____________________________________________________________
1038 THnSparse *AliHFEBeautySpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1039
1040   //
1041   // Unfold and eventually correct for efficiency the bgsubspectrum
1042   //
1043
1044   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1045   if(!mcContainer){
1046     AliError("MC Container not available");
1047     return NULL;
1048   }
1049
1050   if(!fCorrelation){
1051     AliError("No Correlation map available");
1052     return NULL;
1053   }
1054
1055   // Data
1056   AliCFDataGrid *dataGrid = 0x0;
1057   if(bgsubpectrum) {
1058     dataGrid = bgsubpectrum;
1059   }
1060   else {
1061
1062     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1063     if(!dataContainer){
1064       AliError("Data Container not available");
1065       return NULL;
1066     }
1067
1068     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1069   }
1070
1071   // Guessed
1072   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1073   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1074
1075   // Efficiency
1076   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1077   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1078
1079   if(!fBeauty2ndMethod)
1080   {
1081     // Consider parameterized IP cut efficiency
1082     Int_t* bins=new Int_t[1];
1083     bins[0]=kSignalPtBins;
1084
1085     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1086     beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
1087     efficiencyD->Multiply(beffContainer,1);
1088   }
1089
1090
1091   // Unfold
1092
1093   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);//check with MinJung if last arguments are correct (taken from inclusive analysis...
1094   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1095   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1096   if(fSetSmoothing) unfolding.UseSmoothing();
1097   unfolding.Unfold();
1098
1099   // Results
1100   THnSparse* result = unfolding.GetUnfolded();
1101   THnSparse* residual = unfolding.GetEstMeasured();
1102
1103  // QA
1104   TH1D *residualh = (TH1D *) residual->Projection(0);
1105   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1106   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1107   TH1D *afterE = (TH1D *) result->Projection(0);
1108
1109   AliHFEtools::NormaliseBinWidth(residualh);
1110   AliHFEtools::NormaliseBinWidth(beforeE);
1111   AliHFEtools::NormaliseBinWidth(afterE);
1112   fQA->AddResultAt(residualh,AliHFEBeautySpectrumQA::kResidualU);
1113   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterU);
1114   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeU);
1115   fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kUEfficiency);
1116   fQA->DrawUnfolding();
1117
1118   return (THnSparse *) result->Clone();
1119 }
1120
1121 //____________________________________________________________
1122 AliCFDataGrid *AliHFEBeautySpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1123
1124   //
1125   // Apply unfolding and efficiency correction together to bgsubspectrum
1126   //
1127
1128   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1129   if(!mcContainer){
1130     AliError("MC Container not available");
1131     return NULL;
1132   }
1133
1134   // Efficiency
1135   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1136   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1137
1138
1139   if(!fBeauty2ndMethod)
1140   {
1141     // Consider parameterized IP cut efficiency
1142     Int_t* bins=new Int_t[1];
1143     bins[0]=kSignalPtBins;
1144
1145     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1146     beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
1147     efficiencyD->Multiply(beffContainer,1);
1148   }
1149
1150   // Data in the right format
1151   AliCFDataGrid *dataGrid = 0x0;
1152   if(bgsubpectrum) {
1153     dataGrid = bgsubpectrum;
1154   }
1155   else {
1156
1157     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1158     if(!dataContainer){
1159       AliError("Data Container not available");
1160       return NULL;
1161     }
1162
1163     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1164   }
1165
1166   // Correct
1167   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1168   result->ApplyEffCorrection(*efficiencyD);
1169
1170   // QA
1171   TH1D *afterE = (TH1D *) result->Project(0);
1172   AliHFEtools::NormaliseBinWidth(afterE);
1173   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1174   AliHFEtools::NormaliseBinWidth(beforeE);
1175   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1176   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterMCE);
1177   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeMCE);
1178   fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kMCEfficiency);
1179   fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kMC);
1180
1181   return result;
1182
1183 }
1184 //____________________________________________________________
1185 void AliHFEBeautySpectrum::AddTemporaryObject(TObject *o){
1186   //
1187   // Emulate garbage collection on class level
1188   //
1189   if(!fTemporaryObjects) {
1190     fTemporaryObjects= new TList;
1191     fTemporaryObjects->SetOwner();
1192   }
1193   fTemporaryObjects->Add(o);
1194 }
1195
1196 //____________________________________________________________
1197 void AliHFEBeautySpectrum::ClearObject(TObject *o){
1198   //
1199   // Do a safe deletion
1200   //
1201   if(fTemporaryObjects){
1202     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
1203     delete o;
1204   } else{
1205     // just delete
1206     delete o;
1207   }
1208 }
1209 //_________________________________________________________________________
1210 THnSparse* AliHFEBeautySpectrum::GetCharmWeights(Int_t centBin){
1211
1212   //
1213   // Measured D->e based weighting factors
1214   //
1215
1216   const Int_t nDimpp=1;
1217   Int_t nBinpp[nDimpp] = {kSignalPtBins};
1218   Double_t ptbinning1[kSignalPtBins+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.};
1219   Int_t looppt=nBinpp[0];
1220
1221   fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
1222   fWeightCharm->SetBinEdges(0, ptbinning1);
1223
1224   // //if(fBeamType == 0){// Weighting factor for pp
1225   //Double_t weightpp[kSignalPtBins]={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};
1226 //TPC+TOF standard cut 4800
1227   Double_t weightpp[kSignalPtBins]={0.050326, 0.045826, 0.042043, 0.039641, 0.039872, 0.041796, 0.044320, 0.046273, 0.047376, 0.047657, 0.047973, 0.048307, 0.045325, 0.044067, 0.043367, 0.040417, 0.037048, 0.033695, 0.032192, 0.029270, 0.027270, 0.024451, 0.020846, 0.019032, 0.018210, 0.017554, 0.015604, 0.015194, 0.013542, 0.013447, 0.015160, 0.015507, 0.014989, 0.012533, 0.025603};
1228   //}
1229   //else{
1230   //if(centBin == 0){
1231       // Weighting factor for PbPb (0-20%)
1232   Double_t weightPbPb1[kSignalPtBins]={0.641897,  0.640472,  0.615228,  0.650469,  0.737762,  0.847867,  1.009317,  1.158594,  1.307482,  1.476973,  1.551131,  1.677131,  1.785478,  1.888933,  2.017957,  2.074757,  1.926700,  1.869495,  1.546558,  1.222873,  1.160313,  0.903375,  0.799642,  0.706244,  0.705449,  0.599947,  0.719570,  0.499422,  0.703978,  0.477452,  0.325057,  0.093391,  0.096675,  0.000000,  0.000000};
1233   //}
1234   //else{
1235   // Weighting factor for PbPb (40-80%)
1236   Double_t weightPbPb2[kSignalPtBins]={0.181953,  0.173248,  0.166799,  0.182558,  0.206581,  0.236955,  0.279390,  0.329129,  0.365260,  0.423059,  0.452057,  0.482726,  0.462627,  0.537770,  0.584663,  0.579452,  0.587194,  0.499498,  0.443299,  0.398596,  0.376695,  0.322331,  0.260890,  0.374834,  0.249114,  0.310330,  0.287326,  0.243174,  0.758945,  0.138867,  0.170576,  0.107797,  0.011390,  0.000000,  0.000000};
1237   //  }
1238   //}
1239   Double_t weight[kSignalPtBins];
1240   for(Int_t ipt = 0; ipt < kSignalPtBins; ipt++){
1241     if(fBeamType == 0)weight[ipt] = weightpp[ipt];
1242     else if(centBin == 0)weight[ipt] = weightPbPb1[ipt];
1243     else weight[ipt] = weightPbPb2[ipt];
1244   }
1245
1246   //points
1247   Double_t pt[1];
1248   Double_t contents[2];
1249
1250   for(int i=0; i<looppt; i++){
1251     pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
1252     contents[0]=pt[0];
1253     fWeightCharm->Fill(contents,weight[i]);
1254   }
1255
1256
1257   Int_t nDimSparse = fWeightCharm->GetNdimensions();
1258   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1259   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1260
1261   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1262       binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
1263       nBins *= binsvar[iVar];
1264   }
1265
1266   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1267   // loop that sets 0 error in each bin
1268   for (Long_t iBin=0; iBin<nBins; iBin++) {
1269     Long_t bintmp = iBin ;
1270     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1271       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1272       bintmp /= binsvar[iVar] ;
1273     }
1274     fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
1275   }
1276   delete[] binsvar;
1277   delete[] binfill;
1278
1279   return fWeightCharm;
1280 }
1281 //____________________________________________________________________________
1282 void AliHFEBeautySpectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
1283
1284    // TOF PID efficiencies
1285
1286    TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1287    TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,6.);
1288    TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.0);
1289
1290    if(fBeamType == 1){//not in use yet - adapt as soon as possible!
1291      fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1292      fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1293      fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1294    }
1295
1296    TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
1297    cefficiencyParamtof->cd();
1298
1299    AliCFContainer *mccontainermcD = 0x0;
1300    AliCFContainer *mccontaineresdD = 0x0;
1301    TH1D* efficiencysigTOFPIDD;
1302    TH1D* efficiencyTOFPIDD;
1303    TH1D* efficiencysigesdTOFPIDD;
1304    TH1D* efficiencyesdTOFPIDD;
1305    Int_t source = -1; //get parameterized TOF PID efficiencies
1306
1307    // signal sample
1308    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1309    AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
1310    efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1311
1312    // mb sample for double check
1313    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1314    AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
1315    efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1316
1317    // mb sample with reconstructed variables
1318    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1319    AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
1320    efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1321
1322    // mb sample with reconstructed variables
1323    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1324    AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
1325    efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1326
1327    //fill histo
1328    efficiencysigTOFPIDD = (TH1D *) efficiencymcsigParamTOFPID->Project(0);
1329    efficiencyTOFPIDD = (TH1D *) efficiencymcParamTOFPID->Project(0);
1330    efficiencysigesdTOFPIDD = (TH1D *) efficiencysigesdParamTOFPID->Project(0);
1331    efficiencyesdTOFPIDD = (TH1D *) efficiencyesdParamTOFPID->Project(0);
1332    efficiencysigTOFPIDD->SetName("efficiencysigTOFPIDD");
1333    efficiencyTOFPIDD->SetName("efficiencyTOFPIDD");
1334    efficiencysigesdTOFPIDD->SetName("efficiencysigesdTOFPIDD");
1335    efficiencyesdTOFPIDD->SetName("efficiencyesdTOFPIDD");
1336
1337    //fit (mc enhenced sample)
1338    fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1339    efficiencysigTOFPIDD->Fit(fittofpid,"R");
1340    efficiencysigTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1341    fEfficiencyTOFPIDD = efficiencysigTOFPIDD->GetFunction("fittofpid");
1342
1343    //fit (esd enhenced sample)
1344    efficiencysigesdTOFPIDD->Fit(fittofpid,"R");
1345    efficiencysigesdTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1346    fEfficiencyesdTOFPIDD = efficiencysigesdTOFPIDD->GetFunction("fittofpid");
1347
1348    // draw (for PbPb, only 1st bin)
1349    //sig mc
1350    efficiencysigTOFPIDD->SetTitle("");
1351    efficiencysigTOFPIDD->SetStats(0);
1352    efficiencysigTOFPIDD->SetMarkerStyle(25);
1353    efficiencysigTOFPIDD->SetMarkerColor(2);
1354    efficiencysigTOFPIDD->SetLineColor(2);
1355    efficiencysigTOFPIDD->Draw();
1356
1357    //mb mc
1358    efficiencyTOFPIDD->SetTitle("");
1359    efficiencyTOFPIDD->SetStats(0);
1360    efficiencyTOFPIDD->SetMarkerStyle(24);
1361    efficiencyTOFPIDD->SetMarkerColor(4);
1362    efficiencyTOFPIDD->SetLineColor(4);
1363    efficiencyTOFPIDD->Draw("same");
1364
1365    //sig esd
1366    efficiencysigesdTOFPIDD->SetTitle("");
1367    efficiencysigesdTOFPIDD->SetStats(0);
1368    efficiencysigesdTOFPIDD->SetMarkerStyle(25);
1369    efficiencysigesdTOFPIDD->SetMarkerColor(3);
1370    efficiencysigesdTOFPIDD->SetLineColor(3);
1371    efficiencysigesdTOFPIDD->Draw("same");
1372
1373    //mb esd
1374    efficiencyesdTOFPIDD->SetTitle("");
1375    efficiencyesdTOFPIDD->SetStats(0);
1376    efficiencyesdTOFPIDD->SetMarkerStyle(25);
1377    efficiencyesdTOFPIDD->SetMarkerColor(1);
1378    efficiencyesdTOFPIDD->SetLineColor(1);
1379    efficiencyesdTOFPIDD->Draw("same");
1380
1381    //signal mc fit
1382    if(fEfficiencyTOFPIDD){
1383      fEfficiencyTOFPIDD->SetLineColor(2);
1384      fEfficiencyTOFPIDD->Draw("same");
1385    }
1386    //mb esd fit
1387    if(fEfficiencyesdTOFPIDD){
1388        fEfficiencyesdTOFPIDD->SetLineColor(3);
1389        fEfficiencyesdTOFPIDD->Draw("same");
1390      }
1391
1392    TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
1393    legtofeff->AddEntry(efficiencysigTOFPIDD,"TOF PID Step Efficiency","");
1394    legtofeff->AddEntry(efficiencysigTOFPIDD,"vs MC p_{t} for enhenced samples","p");
1395    legtofeff->AddEntry(efficiencyTOFPIDD,"vs MC p_{t} for mb samples","p");
1396    legtofeff->AddEntry(efficiencysigesdTOFPIDD,"vs esd p_{t} for enhenced samples","p");
1397    legtofeff->AddEntry(efficiencyesdTOFPIDD,"vs esd p_{t} for mb samples","p");
1398    legtofeff->Draw("same");
1399
1400
1401    TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
1402    cefficiencyParamIP->cd();
1403    gStyle->SetOptStat(0);
1404
1405    // IP cut efficiencies
1406    AliCFContainer *charmCombined = 0x0;
1407    AliCFContainer *beautyCombined = 0x0;
1408    AliCFContainer *beautyCombinedesd = 0x0;
1409
1410    source = 0; //charm enhenced
1411    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1412    AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
1413    efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1414
1415    charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
1416
1417    source = 1; //beauty enhenced
1418    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1419    AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
1420    efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1421
1422    beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
1423
1424    mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1425    AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
1426    efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1427
1428    beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
1429
1430    source = 0; //charm mb
1431    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1432    AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
1433    efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1434
1435    charmCombined->Add(mccontainermcD);
1436    AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
1437    efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1438
1439    source = 1; //beauty mb
1440    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1441    AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
1442    efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1443
1444    beautyCombined->Add(mccontainermcD);
1445    AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
1446    efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1447
1448    mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1449    AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
1450    efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1451
1452    beautyCombinedesd->Add(mccontaineresdD);
1453    AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
1454    efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1455
1456    source = 2; //conversion mb
1457    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1458    AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
1459    efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1460
1461    source = 3; //non HFE except for the conversion mb
1462    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1463    AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
1464    efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1465
1466    if(fIPEffCombinedSamples){printf("combined samples taken for beauty and charm\n");
1467      fEfficiencyCharmSigD = (TH1D*)efficiencyCharmCombined->Project(0); //signal enhenced + mb
1468      fEfficiencyBeautySigD = (TH1D*)efficiencyBeautyCombined->Project(0); //signal enhenced + mb
1469      fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautyCombinedesd->Project(0); //signal enhenced + mb
1470      }
1471      else{printf("signal enhanced samples taken for beauty and charm\n");
1472        fEfficiencyCharmSigD = (TH1D*)efficiencyCharmSig->Project(0); //signal enhenced only
1473        fEfficiencyBeautySigD = (TH1D*)efficiencyBeautySig->Project(0); //signal enhenced only
1474        fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautySigesd->Project(0); //signal enhenced only
1475      }
1476      fCharmEff = (TH1D*)efficiencyCharm->Project(0); //mb only
1477      fBeautyEff = (TH1D*)efficiencyBeauty->Project(0); //mb only
1478      fConversionEff = (TH1D*)efficiencyConv->Project(0); //mb only
1479      fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); //mb only
1480
1481    if(fBeamType==0){
1482      //AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
1483      AliCFContainer *cfcontainer = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
1484      if(!cfcontainer) return;
1485      AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(cfcontainer,1,0); //mjmj
1486      fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
1487
1488      //AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
1489      AliCFContainer *cfcontainerr = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
1490      if(!cfcontainerr) return;
1491      AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(cfcontainerr,1,0); //mjmj
1492      fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
1493
1494      //MHMH
1495      if(fBeamType == 0){
1496        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1497        fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1498        fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1499
1500        fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1501        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1502        fConversionEffbgc->Fit(fipfitnonhfe,"R");
1503        fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1504      }
1505      //MHMH
1506    }
1507
1508    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1509    fipfit->SetLineColor(2);
1510    fEfficiencyBeautySigD->Fit(fipfit,"R");
1511    fEfficiencyBeautySigD->GetYaxis()->SetTitle("Efficiency");
1512    fEfficiencyIPBeautyD = fEfficiencyBeautySigD->GetFunction("fipfit");
1513
1514    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1515    fipfit->SetLineColor(6);
1516    fEfficiencyBeautySigesdD->Fit(fipfit,"R");
1517    fEfficiencyBeautySigesdD->GetYaxis()->SetTitle("Efficiency");
1518    fEfficiencyIPBeautyesdD = fEfficiencyBeautySigesdD->GetFunction("fipfit");
1519
1520    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1521    fipfit->SetLineColor(1);
1522    fEfficiencyCharmSigD->Fit(fipfit,"R");
1523    fEfficiencyCharmSigD->GetYaxis()->SetTitle("Efficiency");
1524    fEfficiencyIPCharmD = fEfficiencyCharmSigD->GetFunction("fipfit");
1525
1526    //if(0){
1527    if(fIPParameterizedEff){
1528      // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1529      fipfitnonhfe->SetLineColor(3);
1530      if(fBeamType==1){
1531        fConversionEff->Fit(fipfitnonhfe,"R");
1532        fConversionEff->GetYaxis()->SetTitle("Efficiency");
1533        fEfficiencyIPConversionD = fConversionEff->GetFunction("fipfitnonhfe");
1534      }
1535      else{
1536        fConversionEffbgc->Fit(fipfitnonhfe,"R");
1537        fConversionEffbgc->GetYaxis()->SetTitle("Efficiency");
1538        fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1539      }
1540      // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1541      fipfitnonhfe->SetLineColor(4);
1542      if(fBeamType==1){
1543        fNonHFEEff->Fit(fipfitnonhfe,"R");
1544        fNonHFEEff->GetYaxis()->SetTitle("Efficiency");
1545        fEfficiencyIPNonhfeD = fNonHFEEff->GetFunction("fipfitnonhfe");
1546      }
1547      else{
1548        fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1549        fNonHFEEffbgc->GetYaxis()->SetTitle("Efficiency");
1550        fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1551      }
1552    }
1553
1554    // draw
1555    fEfficiencyCharmSigD->SetMarkerStyle(21);
1556    fEfficiencyCharmSigD->SetMarkerColor(1);
1557    fEfficiencyCharmSigD->SetLineColor(1);
1558    fEfficiencyBeautySigD->SetMarkerStyle(21);
1559    fEfficiencyBeautySigD->SetMarkerColor(2);
1560    fEfficiencyBeautySigD->SetLineColor(2);
1561    fEfficiencyBeautySigesdD->SetStats(0);
1562    fEfficiencyBeautySigesdD->SetMarkerStyle(21);
1563    fEfficiencyBeautySigesdD->SetMarkerColor(6);
1564    fEfficiencyBeautySigesdD->SetLineColor(6);
1565    fCharmEff->SetMarkerStyle(24);
1566    fCharmEff->SetMarkerColor(1);
1567    fCharmEff->SetLineColor(1);
1568    fBeautyEff->SetMarkerStyle(24);
1569    fBeautyEff->SetMarkerColor(2);
1570    fBeautyEff->SetLineColor(2);
1571    fConversionEff->SetMarkerStyle(24);
1572    fConversionEff->SetMarkerColor(3);
1573    fConversionEff->SetLineColor(3);
1574    fNonHFEEff->SetMarkerStyle(24);
1575    fNonHFEEff->SetMarkerColor(4);
1576    fNonHFEEff->SetLineColor(4);
1577
1578    fEfficiencyCharmSigD->Draw();
1579    fEfficiencyCharmSigD->GetXaxis()->SetRangeUser(0.0,7.9);
1580    fEfficiencyCharmSigD->GetYaxis()->SetRangeUser(0.0,0.5);
1581
1582    fEfficiencyBeautySigD->Draw("same");
1583    fEfficiencyBeautySigesdD->Draw("same");
1584    if(fBeamType == 1){
1585      fNonHFEEff->Draw("same");
1586      fConversionEff->Draw("same");
1587    }
1588    //fCharmEff->Draw("same");
1589    //fBeautyEff->Draw("same");
1590
1591    if(fBeamType==0){
1592      fConversionEffbgc->SetMarkerStyle(25);
1593      fConversionEffbgc->SetMarkerColor(3);
1594      fConversionEffbgc->SetLineColor(3);
1595      fNonHFEEffbgc->SetMarkerStyle(25);
1596      fNonHFEEffbgc->SetMarkerColor(4);
1597      fNonHFEEffbgc->SetLineColor(4);
1598      fConversionEffbgc->Draw("same");
1599      fNonHFEEffbgc->Draw("same");
1600    }
1601
1602    if(fEfficiencyIPBeautyD)
1603       fEfficiencyIPBeautyD->Draw("same");
1604    if(fEfficiencyIPBeautyesdD)
1605      fEfficiencyIPBeautyesdD->Draw("same");
1606    if( fEfficiencyIPCharmD)
1607      fEfficiencyIPCharmD->Draw("same");
1608    if(fIPParameterizedEff){
1609      if(fEfficiencyIPConversionD)
1610        fEfficiencyIPConversionD->Draw("same");
1611      if(fEfficiencyIPNonhfeD)
1612        fEfficiencyIPNonhfeD->Draw("same");
1613    }
1614    TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
1615    legipeff->AddEntry(fEfficiencyBeautySigD,"IP Step Efficiency","");
1616    legipeff->AddEntry(fEfficiencyBeautySigD,"beauty e","p");
1617    legipeff->AddEntry(fEfficiencyBeautySigesdD,"beauty e(esd pt)","p");
1618    legipeff->AddEntry(fEfficiencyCharmSigD,"charm e","p");
1619    if(fBeamType == 0){
1620      legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
1621      legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
1622    }
1623    else{
1624      legipeff->AddEntry(fConversionEff,"conversion e","p");
1625      legipeff->AddEntry(fNonHFEEff,"Dalitz e","p");
1626    }
1627    legipeff->Draw("same");
1628    gPad->SetGrid();
1629    //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
1630 }
1631
1632 //____________________________________________________________________________
1633 THnSparse* AliHFEBeautySpectrum::GetBeautyIPEff(Bool_t isMCpt){
1634   //
1635   // Return beauty electron IP cut efficiency
1636   //
1637
1638   const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1639   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
1640
1641   Int_t nDim=1;  //dimensions of the efficiency weighting grid
1642
1643   Int_t nBin[1] = {nPtbinning1};
1644
1645   THnSparseF *ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
1646
1647   ipcut->SetBinEdges(0, kPtRange);
1648
1649   Double_t pt[1];
1650   Double_t weight;
1651   Double_t weightErr;
1652   Double_t contents[2];
1653
1654   weight = 1.0;
1655   weightErr = 1.0;
1656
1657   Int_t looppt=nBin[0];
1658
1659   Int_t ibin[2];
1660
1661   for(int i=0; i<looppt; i++)
1662     {
1663       pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1664       if(isMCpt){
1665         if(fEfficiencyIPBeautyD){
1666           weight=fEfficiencyIPBeautyD->Eval(pt[0]);
1667           weightErr = 0;
1668         }
1669         else{
1670           printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1671           weight = fEfficiencyBeautySigD->GetBinContent(i+1);
1672           weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1673         }
1674       }
1675       else{
1676         if(fEfficiencyIPBeautyesdD){
1677           weight=fEfficiencyIPBeautyesdD->Eval(pt[0]);
1678           weightErr = 0;
1679         }
1680         else{
1681           printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1682           weight = fEfficiencyBeautySigesdD->GetBinContent(i+1);
1683           weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1684         }
1685       }
1686
1687       contents[0]=pt[0];
1688       ibin[0]=i+1;
1689
1690       ipcut->Fill(contents,weight);
1691       ipcut->SetBinError(ibin,weightErr);
1692     }
1693
1694   Int_t nDimSparse = ipcut->GetNdimensions();
1695   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1696   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1697
1698   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1699       binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
1700       nBins *= binsvar[iVar];
1701   }
1702
1703   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1704   // loop that sets 0 error in each bin
1705   for (Long_t iBin=0; iBin<nBins; iBin++) {
1706     Long_t bintmp = iBin ;
1707     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1708       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1709       bintmp /= binsvar[iVar] ;
1710     }
1711     //ipcut->SetBinError(binfill,0.); // put 0 everywhere
1712   }
1713
1714   delete[] binsvar;
1715   delete[] binfill;
1716
1717   return ipcut;
1718 }
1719
1720 //____________________________________________________________________________
1721 THnSparse* AliHFEBeautySpectrum::GetPIDxIPEff(Int_t source){
1722   //
1723   // Return PID x IP cut efficiency
1724   //
1725   const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1726   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
1727   Int_t nDim=1;  //dimensions of the efficiency weighting grid
1728
1729   Int_t nBin[1] = {nPtbinning1};
1730
1731   THnSparseF *pideff;
1732   pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
1733   pideff->SetBinEdges(0, kPtRange);
1734
1735   Double_t pt[1];
1736   Double_t weight;
1737   Double_t weightErr;
1738   Double_t contents[2];
1739
1740   weight = 1.0;
1741   weightErr = 1.0;
1742
1743   Int_t looppt=nBin[0];
1744   Int_t ibin[2];
1745
1746   Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
1747   for(int i=0; i<looppt; i++)
1748     {
1749       pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1750
1751       Double_t tofpideff = 1.;
1752       Double_t tofpideffesd = 1.;
1753       if(fEfficiencyTOFPIDD)
1754         tofpideff = fEfficiencyTOFPIDD->Eval(pt[0]);
1755       else{
1756         printf("TOF PID fit failed on conversion. The result is wrong!\n");
1757       }
1758       if(fEfficiencyesdTOFPIDD)
1759         tofpideffesd = fEfficiencyesdTOFPIDD->Eval(pt[0]);
1760       else{
1761         printf("TOF esd PID fit failed on conversion. The result is wrong!\n");
1762       }
1763
1764       //tof pid eff x tpc pid eff x ip cut eff
1765       if(fIPParameterizedEff){
1766         if(source==0) {
1767           if(fEfficiencyIPCharmD){
1768             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1769             weightErr = 0;
1770           }
1771           else{
1772             printf("Fit failed on charm IP cut efficiency\n");
1773             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1774             weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1775           }
1776         }
1777         else if(source==2) {
1778           if(fEfficiencyIPConversionD){
1779             weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD->Eval(pt[0]);
1780             weightErr = 0;
1781           }
1782           else{
1783             printf("Fit failed on conversion IP cut efficiency\n");
1784             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1);
1785             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1786           }
1787         }
1788         else if(source==3) {
1789           if(fEfficiencyIPNonhfeD){
1790             weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD->Eval(pt[0]);
1791             weightErr = 0;
1792           }
1793           else{
1794             printf("Fit failed on dalitz IP cut efficiency\n");
1795             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1);
1796             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1797           }
1798         }
1799       }
1800       else{
1801         if(source==0){
1802           if(fEfficiencyIPCharmD){
1803             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1804             weightErr = 0;
1805           }
1806           else{
1807             printf("Fit failed on charm IP cut efficiency\n");
1808             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1809             weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1810           }
1811         }
1812         else if(source==2){
1813           if(fBeamType==0){
1814             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
1815             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
1816           }
1817           else{
1818             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // conversion
1819             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1820           }
1821         }
1822         else if(source==3){
1823           if(fBeamType==0){
1824             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
1825             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
1826           }
1827           else{
1828             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // Dalitz
1829             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1830           }
1831         }
1832       }
1833
1834       contents[0]=pt[0];
1835       ibin[0]=i+1;
1836
1837       pideff->Fill(contents,weight);
1838       pideff->SetBinError(ibin,weightErr);
1839     }
1840
1841   Int_t nDimSparse = pideff->GetNdimensions();
1842   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1843   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1844
1845   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1846     binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
1847       nBins *= binsvar[iVar];
1848   }
1849
1850   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1851   // loop that sets 0 error in each bin
1852   for (Long_t iBin=0; iBin<nBins; iBin++) {
1853     Long_t bintmp = iBin ;
1854     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1855       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1856       bintmp /= binsvar[iVar] ;
1857     }
1858   }
1859
1860   delete[] binsvar;
1861   delete[] binfill;
1862
1863
1864   return pideff;
1865 }
1866
1867 //__________________________________________________________________________
1868 AliCFDataGrid *AliHFEBeautySpectrum::GetRawBspectra2ndMethod(){
1869  //
1870  // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
1871     //
1872   Int_t nDim = 1;
1873
1874   const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
1875   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};
1876
1877   Int_t nBin[1] = {nPtbinning1};
1878
1879   AliCFDataGrid *rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
1880
1881   THnSparseF *brawspectra;
1882   brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
1883
1884   brawspectra->SetBinEdges(0, kPtRange);
1885
1886   Double_t pt[1];
1887   Double_t yields= 0.;
1888   Double_t valuesb[2];
1889
1890   for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
1891     pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1892
1893     yields = fBSpectrum2ndMethod->GetBinContent(i+1);
1894
1895     valuesb[0]=pt[0];
1896     brawspectra->Fill(valuesb,yields);
1897   }
1898
1899
1900
1901   Int_t nDimSparse = brawspectra->GetNdimensions();
1902   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1903   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1904
1905   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1906     binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
1907     nBins *= binsvar[iVar];
1908   }
1909
1910   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1911   // loop that sets 0 error in each bin
1912   for (Long_t iBin=0; iBin<nBins; iBin++) {
1913     Long_t bintmp = iBin ;
1914     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1915       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1916       bintmp /= binsvar[iVar] ;
1917     }
1918     brawspectra->SetBinError(binfill,0.); // put 0 everywhere
1919   }
1920
1921
1922   rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
1923   TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(0);
1924
1925   new TCanvas;
1926   fBSpectrum2ndMethod->SetMarkerStyle(24);
1927   fBSpectrum2ndMethod->Draw("p");
1928   hRawBeautySpectra->SetMarkerStyle(25);
1929   hRawBeautySpectra->Draw("samep");
1930
1931   delete[] binfill;
1932   delete[] binsvar;
1933
1934   return rawBeautyContainer;
1935 }
1936 //__________________________________________________________________________
1937 void AliHFEBeautySpectrum::CalculateNonHFEsyst(){
1938   //
1939   // Calculate non HFE sys
1940   //
1941   //
1942
1943   if(!fNonHFEsyst)
1944     return;
1945
1946   Double_t evtnorm[1] = {0.0};
1947   if(fNMCbgEvents>0) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
1948
1949   AliCFDataGrid *convSourceGrid[kElecBgSources-3][kBgLevels];
1950   AliCFDataGrid *nonHFESourceGrid[kElecBgSources-3][kBgLevels];
1951
1952   AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
1953   AliCFDataGrid *bgNonHFEGrid[kBgLevels];
1954   AliCFDataGrid *bgConvGrid[kBgLevels];
1955
1956   Int_t stepbackground = 3;
1957   Int_t* bins=new Int_t[1];
1958   const Char_t *bgBase[2] = {"pi0","eta"};
1959
1960   bins[0]=kSignalPtBins;//fConversionEff[centrality]->GetNbinsX();
1961
1962   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1963   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1964
1965   for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1966     for(Int_t iSource = 0; iSource < kElecBgSources-3; iSource++){
1967       convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
1968       weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1969       convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
1970
1971       nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
1972       weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1973       nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
1974     }
1975
1976     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
1977     for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){
1978       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
1979     }
1980     if(!fEtaSyst)
1981       bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
1982
1983     bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
1984     for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
1985       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
1986     }
1987     if(!fEtaSyst)
1988       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
1989
1990     bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
1991     bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
1992     if(fEtaSyst){
1993       bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
1994       bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
1995     }
1996   }
1997
1998   //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0)
1999   AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
2000   TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
2001   for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
2002     bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
2003     bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
2004     bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
2005     bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
2006
2007   //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
2008
2009     hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
2010     hBaseErrors[iErr][0]->Scale(-1.);
2011     hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
2012     hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
2013     if(!fEtaSyst)break;
2014   }
2015
2016   //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
2017   TH1D *hSpeciesErrors[kElecBgSources-4];
2018   for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2019     if(fEtaSyst && (iSource == 1))continue;
2020     hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
2021     TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
2022     hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
2023     hSpeciesErrors[iSource-1]->Scale(0.3);
2024   }
2025
2026   //Int_t firstBgSource = 0;//if eta systematics are not from scaling
2027   //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
2028   TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
2029   TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
2030   TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
2031
2032   TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
2033   TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
2034
2035   for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
2036     Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
2037     pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
2038     pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
2039     if(fEtaSyst){
2040       etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
2041       etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
2042     }
2043     else{ etaErrLow = etaErrUp = 0.;}
2044
2045     Double_t sqrsumErrs= 0;
2046     for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2047       if(fEtaSyst && (iSource == 1))continue;
2048       Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
2049       sqrsumErrs+=(scalingErr*scalingErr);
2050     }
2051     for(Int_t iErr = 0; iErr < 2; iErr++){
2052       for(Int_t iLevel = 0; iLevel < 2; iLevel++){
2053         hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
2054       }
2055       if(!fEtaSyst)break;
2056     }
2057     hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
2058     hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
2059     hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
2060
2061     hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
2062     hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
2063   }
2064
2065   TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
2066   cPiErrors->cd();
2067   cPiErrors->SetLogx();
2068   cPiErrors->SetLogy();
2069   hBaseErrors[0][0]->Draw();
2070   //hBaseErrors[0][1]->SetMarkerColor(kBlack);
2071   //hBaseErrors[0][1]->SetLineColor(kBlack);
2072   //hBaseErrors[0][1]->Draw("SAME");
2073   if(fEtaSyst){
2074     hBaseErrors[1][0]->Draw("SAME");
2075     hBaseErrors[1][0]->SetMarkerColor(kBlack);
2076     hBaseErrors[1][0]->SetLineColor(kBlack);
2077   //hBaseErrors[1][1]->SetMarkerColor(13);
2078   //hBaseErrors[1][1]->SetLineColor(13);
2079   //hBaseErrors[1][1]->Draw("SAME");
2080   }
2081   //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
2082   //hOverallBinScaledErrsUp->SetLineColor(kBlue);
2083   //hOverallBinScaledErrsUp->Draw("SAME");
2084   hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
2085   hOverallBinScaledErrsLow->SetLineColor(kGreen);
2086   hOverallBinScaledErrsLow->Draw("SAME");
2087   hScalingErrors->SetLineColor(kBlue);
2088   hScalingErrors->Draw("SAME");
2089
2090   TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
2091   lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
2092   //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
2093   if(fEtaSyst){
2094   lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
2095   //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
2096   }
2097   lPiErr->AddEntry(hScalingErrors, "scaling error");
2098   lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
2099   //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
2100   lPiErr->Draw("SAME");
2101
2102   //Normalize errors
2103   TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
2104   TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
2105   hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
2106   hLowSystScaled->Scale(evtnorm[0]);
2107   TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
2108   TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
2109   //histograms to be normalized to TGraphErrors
2110   AliHFEtools::NormaliseBinWidth(hNormAllSystErrUp);
2111   AliHFEtools::NormaliseBinWidth(hNormAllSystErrLow);
2112
2113   TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
2114   cNormOvErrs->cd();
2115   cNormOvErrs->SetLogx();
2116   cNormOvErrs->SetLogy();
2117
2118   TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
2119   TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
2120   gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
2121   gOverallSystErrUp->SetMarkerColor(kBlack);
2122   gOverallSystErrUp->SetLineColor(kBlack);
2123   gOverallSystErrLow->SetMarkerColor(kRed);
2124   gOverallSystErrLow->SetLineColor(kRed);
2125   gOverallSystErrUp->Draw("AP");
2126   gOverallSystErrLow->Draw("PSAME");
2127   TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
2128   lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
2129   lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
2130   lAllSys->Draw("same");
2131
2132
2133   AliCFDataGrid *bgYieldGrid;
2134   if(fEtaSyst){
2135     bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
2136   }
2137   bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
2138
2139   TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
2140   TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
2141   hRelErrUp->Divide(hBgYield);
2142   TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
2143   hRelErrLow->Divide(hBgYield);
2144
2145   TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
2146   cRelErrs->cd();
2147   cRelErrs->SetLogx();
2148   hRelErrUp->SetTitle("Relative error of non-HFE background yield");
2149   hRelErrUp->Draw();
2150   hRelErrLow->SetLineColor(kBlack);
2151   hRelErrLow->Draw("SAME");
2152
2153   TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
2154   lRel->AddEntry(hRelErrUp, "upper");
2155   lRel->AddEntry(hRelErrLow, "lower");
2156   lRel->Draw("SAME");
2157
2158   //AliHFEtools::NormaliseBinWidth(hBgYield);
2159   //hBgYield->Scale(evtnorm[0]);
2160
2161
2162   //write histograms/TGraphs to file
2163   TFile *output = new TFile("systHists.root","recreate");
2164
2165   hBgYield->SetName("hBgYield");
2166   hBgYield->Write();
2167   hRelErrUp->SetName("hRelErrUp");
2168   hRelErrUp->Write();
2169   hRelErrLow->SetName("hRelErrLow");
2170   hRelErrLow->Write();
2171   hUpSystScaled->SetName("hOverallSystErrUp");
2172   hUpSystScaled->Write();
2173   hLowSystScaled->SetName("hOverallSystErrLow");
2174   hLowSystScaled->Write();
2175   gOverallSystErrUp->SetName("gOverallSystErrUp");
2176   gOverallSystErrUp->Write();
2177   gOverallSystErrLow->SetName("gOverallSystErrLow");
2178   gOverallSystErrLow->Write();
2179
2180   output->Close();
2181   delete output;
2182 }
2183 //____________________________________________________________
2184 void AliHFEBeautySpectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
2185   //
2186   // Unfold backgrounds to check its sanity
2187   //
2188
2189   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
2190   //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2191   if(!mcContainer){
2192     AliError("MC Container not available");
2193     return;
2194   }
2195
2196   if(!fCorrelation){
2197     AliError("No Correlation map available");
2198     return;
2199   }
2200
2201   // Data
2202   AliCFDataGrid *dataGrid = 0x0;
2203   dataGrid = bgsubpectrum;
2204
2205   // Guessed
2206   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2207   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2208
2209   // Efficiency
2210   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2211   efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
2212
2213   // Unfold background spectra
2214   Int_t nDim=1;
2215   AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2216   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2217   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2218   if(fSetSmoothing) unfolding.UseSmoothing();
2219   unfolding.Unfold();
2220
2221   // Results
2222   THnSparse* result = unfolding.GetUnfolded();
2223   TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
2224   if(fBeamType==1)
2225   {
2226       ctest->Divide(2);
2227       ctest->cd(1);
2228       TH1D* htest1=(TH1D*)result->Projection(0);
2229       htest1->Draw();
2230       ctest->cd(2);
2231       TH1D* htest2=(TH1D*)result->Projection(1);
2232       htest2->Draw();
2233   }
2234
2235
2236   TGraphErrors* unfoldedbgspectrumD = Normalize(result);
2237   if(!unfoldedbgspectrumD) {
2238     AliError("Unfolded background spectrum doesn't exist");
2239   }
2240   else{
2241     TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
2242     if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
2243
2244     file->Close();
2245   }
2246 }