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   CorrectFromTheWidth(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     CorrectFromTheWidth(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     CorrectFromTheWidth(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         CorrectFromTheWidth(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   CorrectFromTheWidth(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   CorrectFromTheWidth(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           backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);
 742         }
 743     }
 744   }
 745   else{//careful: these containers include the non-standard electron sources (K0s, etc.). If necessary, use SetRange in species axis to fix it?
 746     if(source == 0){
 747       // Background Estimate
 748       if(decay == 0)
 749         backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
 750       else if(decay == 1)
 751         backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
 752     }
 753   }
 754   if(!backgroundContainer){
 755     AliError("MC background container not found");
 756     return NULL;
 757   }
 758   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("backgroundGrid","backgroundGrid",*backgroundContainer,stepbackground);
 759
 760   Double_t rangelow = 1.;
 761   Double_t rangeup = 6.;
 762   if(decay == 1) rangelow = 0.9;
 763
 764
 765   const Char_t *dmode[2]={"Conversions","Dalitz decays"};
 766   TF1 *fitHagedorn = new TF1("fitHagedorn", "[0]/TMath::Power(([1]+x/[2]), [3])", rangelow, rangeup);
 767   fNonHFEbg = 0x0;
 768   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
 769   TAxis *axis = h1->GetXaxis();
 770   CorrectFromTheWidth(h1);
 771   if(source == 0){
 772     fitHagedorn->SetParameter(0, 0.15);
 773     fitHagedorn->SetParameter(1, 0.09);
 774     fitHagedorn->SetParameter(2, 8.4);
 775     fitHagedorn->SetParameter(3, 6.3);
 776     // TCanvas *ch1conversion = new TCanvas("ch1conversion","ch1conversion",500,400);
 777     // ch1conversion->cd();
 778     //fHagedorn->SetLineColor(2);
 779     h1->Fit(fitHagedorn,"R");
 780     // h1->Draw();
 781     if(!(fNonHFEbg = h1->GetFunction("fitHagedorn")))printf("electron background fit failed for %s\n",dmode[decay]);
 782   }
 783
 784   Int_t *nBinpp=new Int_t[1];
 785   Int_t *binspp=new Int_t[1];
 786   binspp[0]=kSignalPtBins;// number of pt bins
 787
 788   Int_t looppt=binspp[0];
 789
 790   for(Long_t iBin=1; iBin<= looppt;iBin++){
 791     Double_t iipt= h1->GetBinCenter(iBin);
 792     Double_t iiwidth= axis->GetBinWidth(iBin);
 793     nBinpp[0]=iBin;
 794     Double_t fitFactor = backgroundGrid->GetElement(nBinpp);//if no fit available, just take bin-by-bin information
 795     if(fNonHFEbg)fitFactor = fNonHFEbg->Eval(iipt)*iiwidth;
 796     backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
 797     backgroundGrid->SetElement(nBinpp,fitFactor*evtnorm[0]);
 798   }
 799   //end of workaround for statistical errors
 800   if(source == 0){
 801     AliCFDataGrid *weightedBGContainer = 0x0;
 802     weightedBGContainer = new AliCFDataGrid("weightedBGContainer","weightedBGContainer",1,binspp);
 803     if(decay == 0)
 804       weightedBGContainer->SetGrid(GetPIDxIPEff(2));
 805     else if(decay == 1)
 806       weightedBGContainer->SetGrid(GetPIDxIPEff(3));
 807     if(stepbackground == 3){
 808       backgroundGrid->Multiply(weightedBGContainer,1.0);
 809     }
 810   }
 811   delete[] nBinpp;
 812   delete[] binspp;
 813   return backgroundGrid;
 814 }
 815 //____________________________________________________________
 816 AliCFDataGrid* AliHFEBeautySpectrum::GetCharmBackground(){
 817   //
 818   // calculate charm background; when using centrality binning 2: centBin = 0 for 0-20%, centBin = 5 for 40-80%
 819   //
 820
 821   Int_t nDim = 1;
 822
 823   Double_t evtnorm=0;
 824   if(fNMCEvents) evtnorm= double(fNEvents)/double(fNMCEvents); //mjmjmj
 825   //if(fNMCbgEvents) evtnorm= double(fNEvents)/double(fNMCbgEvents);
 826   printf("events for data: %d",fNEvents);
 827   printf("events for MC: %d",fNMCEvents);
 828   printf("normalization factor for charm: %f",evtnorm);
 829
 830   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
 831   if(!mcContainer){
 832     AliError("MC Container not available");
 833     return NULL;
 834   }
 835
 836   if(!fCorrelation){
 837     AliError("No Correlation map available");
 838     return NULL;
 839   }
 840
 841   AliCFDataGrid *charmBackgroundGrid= 0x0;
 842   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
 843   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
 844   Int_t* bins=new Int_t[2];
 845   bins[0]=charmbgaftertofpid->GetNbinsX();
 846
 847   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
 848   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
 849   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
 850
 851   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
 852   Int_t *nBinpp=new Int_t[1];
 853   Int_t* binspp=new Int_t[1];
 854   binspp[0]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
 855
 856   Int_t looppt=binspp[0];
 857
 858   for(Long_t iBin=1; iBin<= looppt;iBin++){
 859     nBinpp[0]=iBin;
 860     charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
 861     charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
 862   }
 863
 864   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
 865
 866   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
 867   weightedCharmContainer->SetGrid(GetCharmWeights(fTestCentralityLow)); // get charm weighting factors
 868   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
 869   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
 870   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
 871
 872   // Efficiency (set efficiency to 1 for only folding)
 873   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 874   efficiencyD->CalculateEfficiency(0,0);
 875
 876   // Folding
 877   if(fBeamType==0)nDim = 1;
 878   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
 879   folding.SetMaxNumberOfIterations(1);
 880   folding.Unfold();
 881
 882   // Results
 883   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
 884   THnSparse* result=(THnSparse*)result1->Clone();
 885   charmBackgroundGrid->SetGrid(result);
 886   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
 887
 888   //Charm background evaluation plots
 889
 890   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
 891   cCharmBgEval->Divide(3,1);
 892
 893   cCharmBgEval->cd(1);
 894
 895   if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
 896
 897   CorrectFromTheWidth(charmbgaftertofpid);
 898   charmbgaftertofpid->SetMarkerStyle(25);
 899   charmbgaftertofpid->Draw("p");
 900   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
 901   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 902   gPad->SetLogy();
 903
 904   CorrectFromTheWidth(charmbgafteripcut);
 905   charmbgafteripcut->SetMarkerStyle(24);
 906   charmbgafteripcut->Draw("samep");
 907
 908   CorrectFromTheWidth(charmbgafterweight);
 909   charmbgafterweight->SetMarkerStyle(24);
 910   charmbgafterweight->SetMarkerColor(4);
 911   charmbgafterweight->Draw("samep");
 912
 913   CorrectFromTheWidth(charmbgafterfolding);
 914   charmbgafterfolding->SetMarkerStyle(24);
 915   charmbgafterfolding->SetMarkerColor(2);
 916   charmbgafterfolding->Draw("samep");
 917
 918   cCharmBgEval->cd(2);
 919   parametrizedcharmpidipeff->SetMarkerStyle(24);
 920   parametrizedcharmpidipeff->Draw("p");
 921   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 922
 923   cCharmBgEval->cd(3);
 924   charmweightingfc->SetMarkerStyle(24);
 925   charmweightingfc->Draw("p");
 926   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
 927   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 928
 929   cCharmBgEval->cd(1);
 930   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
 931   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
 932   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
 933   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
 934   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
 935   legcharmbg->Draw("same");
 936
 937   cCharmBgEval->cd(2);
 938   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
 939   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
 940   legcharmbg2->Draw("same");
 941
 942   CorrectStatErr(charmBackgroundGrid);
 943   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
 944
 945   delete[] bins;
 946   delete[] nBinpp;
 947   delete[] binspp;
 948
 949   if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
 950
 951   return charmBackgroundGrid;
 952 }
 953 //____________________________________________________________
 954 AliCFDataGrid *AliHFEBeautySpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
 955   //
 956   // Apply TPC pid efficiency correction from parametrisation
 957   //
 958
 959   // Data in the right format
 960   AliCFDataGrid *dataGrid = 0x0;
 961   if(bgsubpectrum) {
 962     dataGrid = bgsubpectrum;
 963   }
 964   else {
 965
 966     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 967     if(!dataContainer){
 968       AliError("Data Container not available");
 969       return NULL;
 970     }
 971     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 972   }
 973   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 974   result->SetName("ParametrizedEfficiencyBefore");
 975   THnSparse *h = result->GetGrid();
 976   Int_t nbdimensions = h->GetNdimensions();
 977   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
 978   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 979   if(!dataContainer){
 980     AliError("Data Container not available");
 981     return NULL;
 982   }
 983   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
 984   dataContainerbis->Add(dataContainerbis,-1.0);
 985
 986   Int_t* coord = new Int_t[nbdimensions];
 987   memset(coord, 0, sizeof(Int_t) * nbdimensions);
 988   Double_t* points = new Double_t[nbdimensions];
 989
 990   ULong64_t nEntries = h->GetNbins();
 991   for (ULong64_t i = 0; i < nEntries; ++i) {
 992     Double_t value = h->GetBinContent(i, coord);
 993     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
 994     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
 995
 996     // Get the bin co-ordinates given an coord
 997     for (Int_t j = 0; j < nbdimensions; ++j)
 998       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
 999
1000     if (!fEfficiencyFunction->IsInside(points))
1001          continue;
1002     TF1::RejectPoint(kFALSE);
1003
1004     // Evaulate function at points
1005     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1006     //printf("Value efficiency is %f\n",valueEfficiency);
1007
1008     if(valueEfficiency > 0.0) {
1009       h->SetBinContent(coord,value/valueEfficiency);
1010       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1011     }
1012     Double_t error = h->GetBinError(i);
1013     h->SetBinError(coord,error/valueEfficiency);
1014     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1015   }
1016
1017   delete[] coord;
1018   delete[] points;
1019
1020   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1021
1022 // QA
1023   TH1D *afterE = (TH1D *) resultt->Project(0);
1024   CorrectFromTheWidth(afterE);
1025   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1026   CorrectFromTheWidth(beforeE);
1027   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterPE);
1028   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforePE);
1029   fQA->AddResultAt(fEfficiencyFunction,AliHFEBeautySpectrumQA::kPEfficiency);
1030   fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kParametrized);
1031
1032   return resultt;
1033 }
1034 //____________________________________________________________
1035 THnSparse *AliHFEBeautySpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
1036
1037   //
1038   // Unfold and eventually correct for efficiency the bgsubspectrum
1039   //
1040
1041   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
1042   if(!mcContainer){
1043     AliError("MC Container not available");
1044     return NULL;
1045   }
1046
1047   if(!fCorrelation){
1048     AliError("No Correlation map available");
1049     return NULL;
1050   }
1051
1052   // Data
1053   AliCFDataGrid *dataGrid = 0x0;
1054   if(bgsubpectrum) {
1055     dataGrid = bgsubpectrum;
1056   }
1057   else {
1058
1059     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1060     if(!dataContainer){
1061       AliError("Data Container not available");
1062       return NULL;
1063     }
1064
1065     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1066   }
1067
1068   // Guessed
1069   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
1070   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
1071
1072   // Efficiency
1073   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1074   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1075
1076   if(!fBeauty2ndMethod)
1077   {
1078     // Consider parameterized IP cut efficiency
1079     Int_t* bins=new Int_t[1];
1080     bins[0]=kSignalPtBins;
1081
1082     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1083     beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
1084     efficiencyD->Multiply(beffContainer,1);
1085   }
1086
1087
1088   // Unfold
1089
1090   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...
1091   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
1092   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
1093   if(fSetSmoothing) unfolding.UseSmoothing();
1094   unfolding.Unfold();
1095
1096   // Results
1097   THnSparse* result = unfolding.GetUnfolded();
1098   THnSparse* residual = unfolding.GetEstMeasured();
1099
1100  // QA
1101   TH1D *residualh = (TH1D *) residual->Projection(0);
1102   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1103   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1104   TH1D *afterE = (TH1D *) result->Projection(0);
1105
1106   CorrectFromTheWidth(residualh);
1107   CorrectFromTheWidth(beforeE);
1108   CorrectFromTheWidth(afterE);
1109   fQA->AddResultAt(residualh,AliHFEBeautySpectrumQA::kResidualU);
1110   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterU);
1111   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeU);
1112   fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kUEfficiency);
1113   fQA->DrawUnfolding();
1114
1115   return (THnSparse *) result->Clone();
1116 }
1117
1118 //____________________________________________________________
1119 AliCFDataGrid *AliHFEBeautySpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
1120
1121   //
1122   // Apply unfolding and efficiency correction together to bgsubspectrum
1123   //
1124
1125   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
1126   if(!mcContainer){
1127     AliError("MC Container not available");
1128     return NULL;
1129   }
1130
1131   // Efficiency
1132   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1133   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
1134
1135
1136   if(!fBeauty2ndMethod)
1137   {
1138     // Consider parameterized IP cut efficiency
1139     Int_t* bins=new Int_t[1];
1140     bins[0]=kSignalPtBins;
1141
1142     AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
1143     beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
1144     efficiencyD->Multiply(beffContainer,1);
1145   }
1146
1147   // Data in the right format
1148   AliCFDataGrid *dataGrid = 0x0;
1149   if(bgsubpectrum) {
1150     dataGrid = bgsubpectrum;
1151   }
1152   else {
1153
1154     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1155     if(!dataContainer){
1156       AliError("Data Container not available");
1157       return NULL;
1158     }
1159
1160     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1161   }
1162
1163   // Correct
1164   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1165   result->ApplyEffCorrection(*efficiencyD);
1166
1167   // QA
1168   TH1D *afterE = (TH1D *) result->Project(0);
1169   CorrectFromTheWidth(afterE);
1170   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1171   CorrectFromTheWidth(beforeE);
1172   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
1173   fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterMCE);
1174   fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeMCE);
1175   fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kMCEfficiency);
1176   fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kMC);
1177
1178   return result;
1179
1180 }
1181 //____________________________________________________________
1182 void AliHFEBeautySpectrum::AddTemporaryObject(TObject *o){
1183   //
1184   // Emulate garbage collection on class level
1185   //
1186   if(!fTemporaryObjects) {
1187     fTemporaryObjects= new TList;
1188     fTemporaryObjects->SetOwner();
1189   }
1190   fTemporaryObjects->Add(o);
1191 }
1192
1193 //____________________________________________________________
1194 void AliHFEBeautySpectrum::ClearObject(TObject *o){
1195   //
1196   // Do a safe deletion
1197   //
1198   if(fTemporaryObjects){
1199     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
1200     delete o;
1201   } else{
1202     // just delete
1203     delete o;
1204   }
1205 }
1206 //_________________________________________________________________________
1207 THnSparse* AliHFEBeautySpectrum::GetCharmWeights(Int_t centBin){
1208
1209   //
1210   // Measured D->e based weighting factors
1211   //
1212
1213   const Int_t nDimpp=1;
1214   Int_t nBinpp[nDimpp] = {kSignalPtBins};
1215   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.};
1216   Int_t looppt=nBinpp[0];
1217
1218   fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
1219   fWeightCharm->SetBinEdges(0, ptbinning1);
1220
1221   // //if(fBeamType == 0){// Weighting factor for pp
1222   //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};
1223 //TPC+TOF standard cut 4800
1224   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};
1225   //}
1226   //else{
1227   //if(centBin == 0){
1228       // Weighting factor for PbPb (0-20%)
1229   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};
1230   //}
1231   //else{
1232   // Weighting factor for PbPb (40-80%)
1233   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};
1234   //  }
1235   //}
1236   Double_t weight[kSignalPtBins];
1237   for(Int_t ipt = 0; ipt < kSignalPtBins; ipt++){
1238     if(fBeamType == 0)weight[ipt] = weightpp[ipt];
1239     else if(centBin == 0)weight[ipt] = weightPbPb1[ipt];
1240     else weight[ipt] = weightPbPb2[ipt];
1241   }
1242
1243   //points
1244   Double_t pt[1];
1245   Double_t contents[2];
1246
1247   for(int i=0; i<looppt; i++){
1248     pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
1249     contents[0]=pt[0];
1250     fWeightCharm->Fill(contents,weight[i]);
1251   }
1252
1253
1254   Int_t nDimSparse = fWeightCharm->GetNdimensions();
1255   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1256   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1257
1258   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1259       binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
1260       nBins *= binsvar[iVar];
1261   }
1262
1263   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1264   // loop that sets 0 error in each bin
1265   for (Long_t iBin=0; iBin<nBins; iBin++) {
1266     Long_t bintmp = iBin ;
1267     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1268       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1269       bintmp /= binsvar[iVar] ;
1270     }
1271     fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
1272   }
1273   delete[] binsvar;
1274   delete[] binfill;
1275
1276   return fWeightCharm;
1277 }
1278 //____________________________________________________________________________
1279 void AliHFEBeautySpectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
1280
1281    // TOF PID efficiencies
1282
1283    TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1284    TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,6.);
1285    TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.0);
1286
1287    if(fBeamType == 1){//not in use yet - adapt as soon as possible!
1288      fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1289      fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
1290      fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1291    }
1292
1293    TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
1294    cefficiencyParamtof->cd();
1295
1296    AliCFContainer *mccontainermcD = 0x0;
1297    AliCFContainer *mccontaineresdD = 0x0;
1298    TH1D* efficiencysigTOFPIDD;
1299    TH1D* efficiencyTOFPIDD;
1300    TH1D* efficiencysigesdTOFPIDD;
1301    TH1D* efficiencyesdTOFPIDD;
1302    Int_t source = -1; //get parameterized TOF PID efficiencies
1303
1304    // signal sample
1305    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1306    AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
1307    efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1308
1309    // mb sample for double check
1310    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1311    AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
1312    efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1313
1314    // mb sample with reconstructed variables
1315    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1316    AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
1317    efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1318
1319    // mb sample with reconstructed variables
1320    if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1321    AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
1322    efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
1323
1324    //fill histo
1325    efficiencysigTOFPIDD = (TH1D *) efficiencymcsigParamTOFPID->Project(0);
1326    efficiencyTOFPIDD = (TH1D *) efficiencymcParamTOFPID->Project(0);
1327    efficiencysigesdTOFPIDD = (TH1D *) efficiencysigesdParamTOFPID->Project(0);
1328    efficiencyesdTOFPIDD = (TH1D *) efficiencyesdParamTOFPID->Project(0);
1329    efficiencysigTOFPIDD->SetName("efficiencysigTOFPIDD");
1330    efficiencyTOFPIDD->SetName("efficiencyTOFPIDD");
1331    efficiencysigesdTOFPIDD->SetName("efficiencysigesdTOFPIDD");
1332    efficiencyesdTOFPIDD->SetName("efficiencyesdTOFPIDD");
1333
1334    //fit (mc enhenced sample)
1335    fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1336    efficiencysigTOFPIDD->Fit(fittofpid,"R");
1337    efficiencysigTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1338    fEfficiencyTOFPIDD = efficiencysigTOFPIDD->GetFunction("fittofpid");
1339
1340    //fit (esd enhenced sample)
1341    efficiencysigesdTOFPIDD->Fit(fittofpid,"R");
1342    efficiencysigesdTOFPIDD->GetYaxis()->SetTitle("Efficiency");
1343    fEfficiencyesdTOFPIDD = efficiencysigesdTOFPIDD->GetFunction("fittofpid");
1344
1345    // draw (for PbPb, only 1st bin)
1346    //sig mc
1347    efficiencysigTOFPIDD->SetTitle("");
1348    efficiencysigTOFPIDD->SetStats(0);
1349    efficiencysigTOFPIDD->SetMarkerStyle(25);
1350    efficiencysigTOFPIDD->SetMarkerColor(2);
1351    efficiencysigTOFPIDD->SetLineColor(2);
1352    efficiencysigTOFPIDD->Draw();
1353
1354    //mb mc
1355    efficiencyTOFPIDD->SetTitle("");
1356    efficiencyTOFPIDD->SetStats(0);
1357    efficiencyTOFPIDD->SetMarkerStyle(24);
1358    efficiencyTOFPIDD->SetMarkerColor(4);
1359    efficiencyTOFPIDD->SetLineColor(4);
1360    efficiencyTOFPIDD->Draw("same");
1361
1362    //sig esd
1363    efficiencysigesdTOFPIDD->SetTitle("");
1364    efficiencysigesdTOFPIDD->SetStats(0);
1365    efficiencysigesdTOFPIDD->SetMarkerStyle(25);
1366    efficiencysigesdTOFPIDD->SetMarkerColor(3);
1367    efficiencysigesdTOFPIDD->SetLineColor(3);
1368    efficiencysigesdTOFPIDD->Draw("same");
1369
1370    //mb esd
1371    efficiencyesdTOFPIDD->SetTitle("");
1372    efficiencyesdTOFPIDD->SetStats(0);
1373    efficiencyesdTOFPIDD->SetMarkerStyle(25);
1374    efficiencyesdTOFPIDD->SetMarkerColor(1);
1375    efficiencyesdTOFPIDD->SetLineColor(1);
1376    efficiencyesdTOFPIDD->Draw("same");
1377
1378    //signal mc fit
1379    if(fEfficiencyTOFPIDD){
1380      fEfficiencyTOFPIDD->SetLineColor(2);
1381      fEfficiencyTOFPIDD->Draw("same");
1382    }
1383    //mb esd fit
1384    if(fEfficiencyesdTOFPIDD){
1385        fEfficiencyesdTOFPIDD->SetLineColor(3);
1386        fEfficiencyesdTOFPIDD->Draw("same");
1387      }
1388
1389    TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
1390    legtofeff->AddEntry(efficiencysigTOFPIDD,"TOF PID Step Efficiency","");
1391    legtofeff->AddEntry(efficiencysigTOFPIDD,"vs MC p_{t} for enhenced samples","p");
1392    legtofeff->AddEntry(efficiencyTOFPIDD,"vs MC p_{t} for mb samples","p");
1393    legtofeff->AddEntry(efficiencysigesdTOFPIDD,"vs esd p_{t} for enhenced samples","p");
1394    legtofeff->AddEntry(efficiencyesdTOFPIDD,"vs esd p_{t} for mb samples","p");
1395    legtofeff->Draw("same");
1396
1397
1398    TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
1399    cefficiencyParamIP->cd();
1400    gStyle->SetOptStat(0);
1401
1402    // IP cut efficiencies
1403    AliCFContainer *charmCombined = 0x0;
1404    AliCFContainer *beautyCombined = 0x0;
1405    AliCFContainer *beautyCombinedesd = 0x0;
1406
1407    source = 0; //charm enhenced
1408    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1409    AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
1410    efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1411
1412    charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
1413
1414    source = 1; //beauty enhenced
1415    mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1416    AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
1417    efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1418
1419    beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
1420
1421    mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1422    AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
1423    efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1424
1425    beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
1426
1427    source = 0; //charm mb
1428    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1429    AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
1430    efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1431
1432    charmCombined->Add(mccontainermcD);
1433    AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
1434    efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1435
1436    source = 1; //beauty mb
1437    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1438    AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
1439    efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1440
1441    beautyCombined->Add(mccontainermcD);
1442    AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
1443    efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1444
1445    mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1446    AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
1447    efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1448
1449    beautyCombinedesd->Add(mccontaineresdD);
1450    AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
1451    efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
1452
1453    source = 2; //conversion mb
1454    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1455    AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
1456    efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1457
1458    source = 3; //non HFE except for the conversion mb
1459    mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
1460    AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
1461    efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
1462
1463    if(fIPEffCombinedSamples){printf("combined samples taken for beauty and charm\n");
1464      fEfficiencyCharmSigD = (TH1D*)efficiencyCharmCombined->Project(0); //signal enhenced + mb
1465      fEfficiencyBeautySigD = (TH1D*)efficiencyBeautyCombined->Project(0); //signal enhenced + mb
1466      fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautyCombinedesd->Project(0); //signal enhenced + mb
1467      }
1468      else{printf("signal enhanced samples taken for beauty and charm\n");
1469        fEfficiencyCharmSigD = (TH1D*)efficiencyCharmSig->Project(0); //signal enhenced only
1470        fEfficiencyBeautySigD = (TH1D*)efficiencyBeautySig->Project(0); //signal enhenced only
1471        fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautySigesd->Project(0); //signal enhenced only
1472      }
1473      fCharmEff = (TH1D*)efficiencyCharm->Project(0); //mb only
1474      fBeautyEff = (TH1D*)efficiencyBeauty->Project(0); //mb only
1475      fConversionEff = (TH1D*)efficiencyConv->Project(0); //mb only
1476      fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); //mb only
1477
1478    if(fBeamType==0){
1479      //AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
1480      AliCFContainer *cfcontainer = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
1481      if(!cfcontainer) return;
1482      AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(cfcontainer,1,0); //mjmj
1483      fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
1484
1485      //AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
1486      AliCFContainer *cfcontainerr = GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig);
1487      if(!cfcontainerr) return;
1488      AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(cfcontainerr,1,0); //mjmj
1489      fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
1490
1491      //MHMH
1492      if(fBeamType == 0){
1493        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1494        fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1495        fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1496
1497        fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
1498        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1499        fConversionEffbgc->Fit(fipfitnonhfe,"R");
1500        fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1501      }
1502      //MHMH
1503    }
1504
1505    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1506    fipfit->SetLineColor(2);
1507    fEfficiencyBeautySigD->Fit(fipfit,"R");
1508    fEfficiencyBeautySigD->GetYaxis()->SetTitle("Efficiency");
1509    fEfficiencyIPBeautyD = fEfficiencyBeautySigD->GetFunction("fipfit");
1510
1511    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1512    fipfit->SetLineColor(6);
1513    fEfficiencyBeautySigesdD->Fit(fipfit,"R");
1514    fEfficiencyBeautySigesdD->GetYaxis()->SetTitle("Efficiency");
1515    fEfficiencyIPBeautyesdD = fEfficiencyBeautySigesdD->GetFunction("fipfit");
1516
1517    fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1518    fipfit->SetLineColor(1);
1519    fEfficiencyCharmSigD->Fit(fipfit,"R");
1520    fEfficiencyCharmSigD->GetYaxis()->SetTitle("Efficiency");
1521    fEfficiencyIPCharmD = fEfficiencyCharmSigD->GetFunction("fipfit");
1522
1523    //if(0){
1524    if(fIPParameterizedEff){
1525      // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1526      fipfitnonhfe->SetLineColor(3);
1527      if(fBeamType==1){
1528        fConversionEff->Fit(fipfitnonhfe,"R");
1529        fConversionEff->GetYaxis()->SetTitle("Efficiency");
1530        fEfficiencyIPConversionD = fConversionEff->GetFunction("fipfitnonhfe");
1531      }
1532      else{
1533        fConversionEffbgc->Fit(fipfitnonhfe,"R");
1534        fConversionEffbgc->GetYaxis()->SetTitle("Efficiency");
1535        fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
1536      }
1537      // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
1538      fipfitnonhfe->SetLineColor(4);
1539      if(fBeamType==1){
1540        fNonHFEEff->Fit(fipfitnonhfe,"R");
1541        fNonHFEEff->GetYaxis()->SetTitle("Efficiency");
1542        fEfficiencyIPNonhfeD = fNonHFEEff->GetFunction("fipfitnonhfe");
1543      }
1544      else{
1545        fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
1546        fNonHFEEffbgc->GetYaxis()->SetTitle("Efficiency");
1547        fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
1548      }
1549    }
1550
1551    // draw
1552    fEfficiencyCharmSigD->SetMarkerStyle(21);
1553    fEfficiencyCharmSigD->SetMarkerColor(1);
1554    fEfficiencyCharmSigD->SetLineColor(1);
1555    fEfficiencyBeautySigD->SetMarkerStyle(21);
1556    fEfficiencyBeautySigD->SetMarkerColor(2);
1557    fEfficiencyBeautySigD->SetLineColor(2);
1558    fEfficiencyBeautySigesdD->SetStats(0);
1559    fEfficiencyBeautySigesdD->SetMarkerStyle(21);
1560    fEfficiencyBeautySigesdD->SetMarkerColor(6);
1561    fEfficiencyBeautySigesdD->SetLineColor(6);
1562    fCharmEff->SetMarkerStyle(24);
1563    fCharmEff->SetMarkerColor(1);
1564    fCharmEff->SetLineColor(1);
1565    fBeautyEff->SetMarkerStyle(24);
1566    fBeautyEff->SetMarkerColor(2);
1567    fBeautyEff->SetLineColor(2);
1568    fConversionEff->SetMarkerStyle(24);
1569    fConversionEff->SetMarkerColor(3);
1570    fConversionEff->SetLineColor(3);
1571    fNonHFEEff->SetMarkerStyle(24);
1572    fNonHFEEff->SetMarkerColor(4);
1573    fNonHFEEff->SetLineColor(4);
1574
1575    fEfficiencyCharmSigD->Draw();
1576    fEfficiencyCharmSigD->GetXaxis()->SetRangeUser(0.0,7.9);
1577    fEfficiencyCharmSigD->GetYaxis()->SetRangeUser(0.0,0.5);
1578
1579    fEfficiencyBeautySigD->Draw("same");
1580    fEfficiencyBeautySigesdD->Draw("same");
1581    if(fBeamType == 1){
1582      fNonHFEEff->Draw("same");
1583      fConversionEff->Draw("same");
1584    }
1585    //fCharmEff->Draw("same");
1586    //fBeautyEff->Draw("same");
1587
1588    if(fBeamType==0){
1589      fConversionEffbgc->SetMarkerStyle(25);
1590      fConversionEffbgc->SetMarkerColor(3);
1591      fConversionEffbgc->SetLineColor(3);
1592      fNonHFEEffbgc->SetMarkerStyle(25);
1593      fNonHFEEffbgc->SetMarkerColor(4);
1594      fNonHFEEffbgc->SetLineColor(4);
1595      fConversionEffbgc->Draw("same");
1596      fNonHFEEffbgc->Draw("same");
1597    }
1598
1599    if(fEfficiencyIPBeautyD)
1600       fEfficiencyIPBeautyD->Draw("same");
1601    if(fEfficiencyIPBeautyesdD)
1602      fEfficiencyIPBeautyesdD->Draw("same");
1603    if( fEfficiencyIPCharmD)
1604      fEfficiencyIPCharmD->Draw("same");
1605    if(fIPParameterizedEff){
1606      if(fEfficiencyIPConversionD)
1607        fEfficiencyIPConversionD->Draw("same");
1608      if(fEfficiencyIPNonhfeD)
1609        fEfficiencyIPNonhfeD->Draw("same");
1610    }
1611    TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
1612    legipeff->AddEntry(fEfficiencyBeautySigD,"IP Step Efficiency","");
1613    legipeff->AddEntry(fEfficiencyBeautySigD,"beauty e","p");
1614    legipeff->AddEntry(fEfficiencyBeautySigesdD,"beauty e(esd pt)","p");
1615    legipeff->AddEntry(fEfficiencyCharmSigD,"charm e","p");
1616    if(fBeamType == 0){
1617      legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
1618      legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
1619    }
1620    else{
1621      legipeff->AddEntry(fConversionEff,"conversion e","p");
1622      legipeff->AddEntry(fNonHFEEff,"Dalitz e","p");
1623    }
1624    legipeff->Draw("same");
1625    gPad->SetGrid();
1626    //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
1627 }
1628
1629 //____________________________________________________________________________
1630 THnSparse* AliHFEBeautySpectrum::GetBeautyIPEff(Bool_t isMCpt){
1631   //
1632   // Return beauty electron IP cut efficiency
1633   //
1634
1635   const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1636   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
1637
1638   Int_t nDim=1;  //dimensions of the efficiency weighting grid
1639
1640   Int_t nBin[1] = {nPtbinning1};
1641
1642   THnSparseF *ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
1643
1644   ipcut->SetBinEdges(0, kPtRange);
1645
1646   Double_t pt[1];
1647   Double_t weight;
1648   Double_t weightErr;
1649   Double_t contents[2];
1650
1651   weight = 1.0;
1652   weightErr = 1.0;
1653
1654   Int_t looppt=nBin[0];
1655
1656   Int_t ibin[2];
1657
1658   for(int i=0; i<looppt; i++)
1659     {
1660       pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1661       if(isMCpt){
1662         if(fEfficiencyIPBeautyD){
1663           weight=fEfficiencyIPBeautyD->Eval(pt[0]);
1664           weightErr = 0;
1665         }
1666         else{
1667           printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1668           weight = fEfficiencyBeautySigD->GetBinContent(i+1);
1669           weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1670         }
1671       }
1672       else{
1673         if(fEfficiencyIPBeautyesdD){
1674           weight=fEfficiencyIPBeautyesdD->Eval(pt[0]);
1675           weightErr = 0;
1676         }
1677         else{
1678           printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
1679           weight = fEfficiencyBeautySigesdD->GetBinContent(i+1);
1680           weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
1681         }
1682       }
1683
1684       contents[0]=pt[0];
1685       ibin[0]=i+1;
1686
1687       ipcut->Fill(contents,weight);
1688       ipcut->SetBinError(ibin,weightErr);
1689     }
1690
1691   Int_t nDimSparse = ipcut->GetNdimensions();
1692   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1693   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1694
1695   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1696       binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
1697       nBins *= binsvar[iVar];
1698   }
1699
1700   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1701   // loop that sets 0 error in each bin
1702   for (Long_t iBin=0; iBin<nBins; iBin++) {
1703     Long_t bintmp = iBin ;
1704     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1705       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1706       bintmp /= binsvar[iVar] ;
1707     }
1708     //ipcut->SetBinError(binfill,0.); // put 0 everywhere
1709   }
1710
1711   delete[] binsvar;
1712   delete[] binfill;
1713
1714   return ipcut;
1715 }
1716
1717 //____________________________________________________________________________
1718 THnSparse* AliHFEBeautySpectrum::GetPIDxIPEff(Int_t source){
1719   //
1720   // Return PID x IP cut efficiency
1721   //
1722   const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
1723   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
1724   Int_t nDim=1;  //dimensions of the efficiency weighting grid
1725
1726   Int_t nBin[1] = {nPtbinning1};
1727
1728   THnSparseF *pideff;
1729   pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
1730   pideff->SetBinEdges(0, kPtRange);
1731
1732   Double_t pt[1];
1733   Double_t weight;
1734   Double_t weightErr;
1735   Double_t contents[2];
1736
1737   weight = 1.0;
1738   weightErr = 1.0;
1739
1740   Int_t looppt=nBin[0];
1741   Int_t ibin[2];
1742
1743   Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
1744   for(int i=0; i<looppt; i++)
1745     {
1746       pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1747
1748       Double_t tofpideff = 1.;
1749       Double_t tofpideffesd = 1.;
1750       if(fEfficiencyTOFPIDD)
1751         tofpideff = fEfficiencyTOFPIDD->Eval(pt[0]);
1752       else{
1753         printf("TOF PID fit failed on conversion. The result is wrong!\n");
1754       }
1755       if(fEfficiencyesdTOFPIDD)
1756         tofpideffesd = fEfficiencyesdTOFPIDD->Eval(pt[0]);
1757       else{
1758         printf("TOF esd PID fit failed on conversion. The result is wrong!\n");
1759       }
1760
1761       //tof pid eff x tpc pid eff x ip cut eff
1762       if(fIPParameterizedEff){
1763         if(source==0) {
1764           if(fEfficiencyIPCharmD){
1765             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1766             weightErr = 0;
1767           }
1768           else{
1769             printf("Fit failed on charm IP cut efficiency\n");
1770             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1771             weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1772           }
1773         }
1774         else if(source==2) {
1775           if(fEfficiencyIPConversionD){
1776             weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD->Eval(pt[0]);
1777             weightErr = 0;
1778           }
1779           else{
1780             printf("Fit failed on conversion IP cut efficiency\n");
1781             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1);
1782             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1783           }
1784         }
1785         else if(source==3) {
1786           if(fEfficiencyIPNonhfeD){
1787             weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD->Eval(pt[0]);
1788             weightErr = 0;
1789           }
1790           else{
1791             printf("Fit failed on dalitz IP cut efficiency\n");
1792             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1);
1793             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1794           }
1795         }
1796       }
1797       else{
1798         if(source==0){
1799           if(fEfficiencyIPCharmD){
1800             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
1801             weightErr = 0;
1802           }
1803           else{
1804             printf("Fit failed on charm IP cut efficiency\n");
1805             weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
1806             weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
1807           }
1808         }
1809         else if(source==2){
1810           if(fBeamType==0){
1811             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
1812             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
1813           }
1814           else{
1815             weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // conversion
1816             weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
1817           }
1818         }
1819         else if(source==3){
1820           if(fBeamType==0){
1821             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
1822             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
1823           }
1824           else{
1825             weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // Dalitz
1826             weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
1827           }
1828         }
1829       }
1830
1831       contents[0]=pt[0];
1832       ibin[0]=i+1;
1833
1834       pideff->Fill(contents,weight);
1835       pideff->SetBinError(ibin,weightErr);
1836     }
1837
1838   Int_t nDimSparse = pideff->GetNdimensions();
1839   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1840   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1841
1842   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1843     binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
1844       nBins *= binsvar[iVar];
1845   }
1846
1847   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1848   // loop that sets 0 error in each bin
1849   for (Long_t iBin=0; iBin<nBins; iBin++) {
1850     Long_t bintmp = iBin ;
1851     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1852       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1853       bintmp /= binsvar[iVar] ;
1854     }
1855   }
1856
1857   delete[] binsvar;
1858   delete[] binfill;
1859
1860
1861   return pideff;
1862 }
1863
1864 //__________________________________________________________________________
1865 AliCFDataGrid *AliHFEBeautySpectrum::GetRawBspectra2ndMethod(){
1866  //
1867  // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
1868     //
1869   Int_t nDim = 1;
1870
1871   const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
1872   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};
1873
1874   Int_t nBin[1] = {nPtbinning1};
1875
1876   AliCFDataGrid *rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
1877
1878   THnSparseF *brawspectra;
1879   brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
1880
1881   brawspectra->SetBinEdges(0, kPtRange);
1882
1883   Double_t pt[1];
1884   Double_t yields= 0.;
1885   Double_t valuesb[2];
1886
1887   for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
1888     pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
1889
1890     yields = fBSpectrum2ndMethod->GetBinContent(i+1);
1891
1892     valuesb[0]=pt[0];
1893     brawspectra->Fill(valuesb,yields);
1894   }
1895
1896
1897
1898   Int_t nDimSparse = brawspectra->GetNdimensions();
1899   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
1900   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
1901
1902   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1903     binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
1904     nBins *= binsvar[iVar];
1905   }
1906
1907   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
1908   // loop that sets 0 error in each bin
1909   for (Long_t iBin=0; iBin<nBins; iBin++) {
1910     Long_t bintmp = iBin ;
1911     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
1912       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
1913       bintmp /= binsvar[iVar] ;
1914     }
1915     brawspectra->SetBinError(binfill,0.); // put 0 everywhere
1916   }
1917
1918
1919   rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
1920   TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(0);
1921
1922   new TCanvas;
1923   fBSpectrum2ndMethod->SetMarkerStyle(24);
1924   fBSpectrum2ndMethod->Draw("p");
1925   hRawBeautySpectra->SetMarkerStyle(25);
1926   hRawBeautySpectra->Draw("samep");
1927
1928   delete[] binfill;
1929   delete[] binsvar;
1930
1931   return rawBeautyContainer;
1932 }
1933 //__________________________________________________________________________
1934 void AliHFEBeautySpectrum::CalculateNonHFEsyst(){
1935   //
1936   // Calculate non HFE sys
1937   //
1938   //
1939
1940   if(!fNonHFEsyst)
1941     return;
1942
1943   Double_t evtnorm[1] = {0.0};
1944   if(fNMCbgEvents>0) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
1945
1946   AliCFDataGrid *convSourceGrid[kElecBgSources-3][kBgLevels];
1947   AliCFDataGrid *nonHFESourceGrid[kElecBgSources-3][kBgLevels];
1948
1949   AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
1950   AliCFDataGrid *bgNonHFEGrid[kBgLevels];
1951   AliCFDataGrid *bgConvGrid[kBgLevels];
1952
1953   Int_t stepbackground = 3;
1954   Int_t* bins=new Int_t[1];
1955   const Char_t *bgBase[2] = {"pi0","eta"};
1956
1957   bins[0]=kSignalPtBins;//fConversionEff[centrality]->GetNbinsX();
1958
1959   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
1960   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
1961
1962   for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1963     for(Int_t iSource = 0; iSource < kElecBgSources-3; iSource++){
1964       convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
1965       weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1966       convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
1967
1968       nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
1969       weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1970       nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
1971     }
1972
1973     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
1974     for(Int_t iSource = 2; iSource < kElecBgSources-3; iSource++){
1975       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
1976     }
1977     if(!fEtaSyst)
1978       bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
1979
1980     bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
1981     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)
1982       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
1983     }
1984     if(!fEtaSyst)
1985       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
1986
1987     bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
1988     bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
1989     if(fEtaSyst){
1990       bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
1991       bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
1992     }
1993   }
1994
1995   //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)
1996   AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
1997   TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
1998   for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
1999     bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
2000     bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
2001     bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
2002     bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
2003
2004   //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
2005
2006     hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
2007     hBaseErrors[iErr][0]->Scale(-1.);
2008     hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
2009     hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
2010     if(!fEtaSyst)break;
2011   }
2012
2013   //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
2014   TH1D *hSpeciesErrors[kElecBgSources-4];
2015   for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2016     if(fEtaSyst && (iSource == 1))continue;
2017     hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
2018     TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
2019     hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
2020     hSpeciesErrors[iSource-1]->Scale(0.3);
2021   }
2022
2023   //Int_t firstBgSource = 0;//if eta systematics are not from scaling
2024   //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
2025   TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
2026   TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
2027   TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
2028
2029   TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
2030   TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
2031
2032   for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
2033     Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
2034     pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
2035     pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
2036     if(fEtaSyst){
2037       etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
2038       etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
2039     }
2040     else{ etaErrLow = etaErrUp = 0.;}
2041
2042     Double_t sqrsumErrs= 0;
2043     for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
2044       if(fEtaSyst && (iSource == 1))continue;
2045       Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
2046       sqrsumErrs+=(scalingErr*scalingErr);
2047     }
2048     for(Int_t iErr = 0; iErr < 2; iErr++){
2049       for(Int_t iLevel = 0; iLevel < 2; iLevel++){
2050         hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
2051       }
2052       if(!fEtaSyst)break;
2053     }
2054     hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
2055     hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
2056     hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
2057
2058     hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
2059     hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
2060   }
2061
2062   TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
2063   cPiErrors->cd();
2064   cPiErrors->SetLogx();
2065   cPiErrors->SetLogy();
2066   hBaseErrors[0][0]->Draw();
2067   //hBaseErrors[0][1]->SetMarkerColor(kBlack);
2068   //hBaseErrors[0][1]->SetLineColor(kBlack);
2069   //hBaseErrors[0][1]->Draw("SAME");
2070   if(fEtaSyst){
2071     hBaseErrors[1][0]->Draw("SAME");
2072     hBaseErrors[1][0]->SetMarkerColor(kBlack);
2073     hBaseErrors[1][0]->SetLineColor(kBlack);
2074   //hBaseErrors[1][1]->SetMarkerColor(13);
2075   //hBaseErrors[1][1]->SetLineColor(13);
2076   //hBaseErrors[1][1]->Draw("SAME");
2077   }
2078   //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
2079   //hOverallBinScaledErrsUp->SetLineColor(kBlue);
2080   //hOverallBinScaledErrsUp->Draw("SAME");
2081   hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
2082   hOverallBinScaledErrsLow->SetLineColor(kGreen);
2083   hOverallBinScaledErrsLow->Draw("SAME");
2084   hScalingErrors->SetLineColor(kBlue);
2085   hScalingErrors->Draw("SAME");
2086
2087   TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
2088   lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
2089   //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
2090   if(fEtaSyst){
2091   lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
2092   //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
2093   }
2094   lPiErr->AddEntry(hScalingErrors, "scaling error");
2095   lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
2096   //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
2097   lPiErr->Draw("SAME");
2098
2099   //Normalize errors
2100   TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
2101   TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
2102   hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
2103   hLowSystScaled->Scale(evtnorm[0]);
2104   TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
2105   TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
2106   //histograms to be normalized to TGraphErrors
2107   CorrectFromTheWidth(hNormAllSystErrUp);
2108   CorrectFromTheWidth(hNormAllSystErrLow);
2109
2110   TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
2111   cNormOvErrs->cd();
2112   cNormOvErrs->SetLogx();
2113   cNormOvErrs->SetLogy();
2114
2115   TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
2116   TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
2117   gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
2118   gOverallSystErrUp->SetMarkerColor(kBlack);
2119   gOverallSystErrUp->SetLineColor(kBlack);
2120   gOverallSystErrLow->SetMarkerColor(kRed);
2121   gOverallSystErrLow->SetLineColor(kRed);
2122   gOverallSystErrUp->Draw("AP");
2123   gOverallSystErrLow->Draw("PSAME");
2124   TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
2125   lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
2126   lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
2127   lAllSys->Draw("same");
2128
2129
2130   AliCFDataGrid *bgYieldGrid;
2131   if(fEtaSyst){
2132     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.
2133   }
2134   bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
2135
2136   TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
2137   TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
2138   hRelErrUp->Divide(hBgYield);
2139   TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
2140   hRelErrLow->Divide(hBgYield);
2141
2142   TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
2143   cRelErrs->cd();
2144   cRelErrs->SetLogx();
2145   hRelErrUp->SetTitle("Relative error of non-HFE background yield");
2146   hRelErrUp->Draw();
2147   hRelErrLow->SetLineColor(kBlack);
2148   hRelErrLow->Draw("SAME");
2149
2150   TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
2151   lRel->AddEntry(hRelErrUp, "upper");
2152   lRel->AddEntry(hRelErrLow, "lower");
2153   lRel->Draw("SAME");
2154
2155   //CorrectFromTheWidth(hBgYield);
2156   //hBgYield->Scale(evtnorm[0]);
2157
2158
2159   //write histograms/TGraphs to file
2160   TFile *output = new TFile("systHists.root","recreate");
2161
2162   hBgYield->SetName("hBgYield");
2163   hBgYield->Write();
2164   hRelErrUp->SetName("hRelErrUp");
2165   hRelErrUp->Write();
2166   hRelErrLow->SetName("hRelErrLow");
2167   hRelErrLow->Write();
2168   hUpSystScaled->SetName("hOverallSystErrUp");
2169   hUpSystScaled->Write();
2170   hLowSystScaled->SetName("hOverallSystErrLow");
2171   hLowSystScaled->Write();
2172   gOverallSystErrUp->SetName("gOverallSystErrUp");
2173   gOverallSystErrUp->Write();
2174   gOverallSystErrLow->SetName("gOverallSystErrLow");
2175   gOverallSystErrLow->Write();
2176
2177   output->Close();
2178   delete output;
2179 }
2180 //____________________________________________________________
2181 void AliHFEBeautySpectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
2182   //
2183   // Unfold backgrounds to check its sanity
2184   //
2185
2186   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
2187   //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2188   if(!mcContainer){
2189     AliError("MC Container not available");
2190     return;
2191   }
2192
2193   if(!fCorrelation){
2194     AliError("No Correlation map available");
2195     return;
2196   }
2197
2198   // Data
2199   AliCFDataGrid *dataGrid = 0x0;
2200   dataGrid = bgsubpectrum;
2201
2202   // Guessed
2203   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2204   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2205
2206   // Efficiency
2207   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2208   efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
2209
2210   // Unfold background spectra
2211   Int_t nDim=1;
2212   AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2213   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2214   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2215   if(fSetSmoothing) unfolding.UseSmoothing();
2216   unfolding.Unfold();
2217
2218   // Results
2219   THnSparse* result = unfolding.GetUnfolded();
2220   TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
2221   if(fBeamType==1)
2222   {
2223       ctest->Divide(2);
2224       ctest->cd(1);
2225       TH1D* htest1=(TH1D*)result->Projection(0);
2226       htest1->Draw();
2227       ctest->cd(2);
2228       TH1D* htest2=(TH1D*)result->Projection(1);
2229       htest2->Draw();
2230   }
2231
2232
2233   TGraphErrors* unfoldedbgspectrumD = Normalize(result);
2234   if(!unfoldedbgspectrumD) {
2235     AliError("Unfolded background spectrum doesn't exist");
2236   }
2237   else{
2238     TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
2239     if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
2240
2241     file->Close();
2242   }
2243 }