PWG3/hfe/AliHFEelecbackground.cxx

   1  /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //
  16 //
  17 // First implementation of a class
  18 // to reject tagged electron coming from conversion, pi0 and eta
  19 // by calculating the e+e- invariant mass
  20 // of the tagged electron with other tracks
  21 // after looser cuts for the partner.
  22 // PostProcess should extract the background yield
  23 // If running with MC, it can be compared to the expected background yield
  24 //
  25 // Authors:
  26 //   Raphaelle Bailhache <rbailhache@ikf.uni-frankfurt.de > <R.Bailhache@gsi.de >
  27 //
  28 //
  29
  30 #include <THnSparse.h>
  31 #include <TParticle.h>
  32 #include <TFile.h>
  33 #include <TCanvas.h>
  34 #include <TAxis.h>
  35 #include <TH1D.h>
  36 #include <TH2D.h>
  37 #include <TH1F.h>
  38 #include <TH2F.h>
  39 #include <TF1.h>
  40 #include <TString.h>
  41 #include <TLegend.h>
  42 #include <TStyle.h>
  43
  44 #include <AliESDEvent.h>
  45 #include <AliAODEvent.h>
  46 #include <AliESDtrack.h>
  47 #include <AliAODTrack.h>
  48 #include "AliHFEelecbackground.h"
  49 #include <AliMCEvent.h>
  50 #include <AliStack.h>
  51 #include <AliKFParticle.h>
  52 #include "AliCFContainer.h"
  53 #include "AliHFEpid.h"
  54 #include "AliESDpid.h"
  55 #include "AliLog.h"
  56 #include "AliITSPIDResponse.h"
  57 #include "AliTPCPIDResponse.h"
  58
  59 ClassImp(AliHFEelecbackground)
  60
  61 Bool_t AliHFEelecbackground::fgUseMCPID = kFALSE;
  62 const Double_t    AliHFEelecbackground::fgkMe = 0.00051099892;
  63
  64 //___________________________________________________________________________________________
  65 AliHFEelecbackground::AliHFEelecbackground():
  66   fhtmp(0x0)
  67   ,fhtmpf(0x0)
  68   ,fhtmpp(0x0)
  69   ,fESD1(0x0)
  70   ,fAOD1(0x0)
  71   ,fMCEvent(0x0)
  72   ,fBz(0)
  73   ,fkVertex(0x0)
  74   ,fPtESD(0.0)
  75   ,fIndexTrack(0)
  76   ,fPdg(0)
  77   ,fLabMother(-1)
  78   ,fIsFrom(-1)
  79   ,fMotherGamma(-1)
  80   ,fMotherPi0(-1)
  81   ,fMotherC(-1)
  82   ,fMotherB(-1)
  83   ,fMotherEta(-1)
  84   ,fIsPartner(kFALSE)
  85   ,fIsSplittedTrack(kFALSE)
  86   ,fOpeningAngleCut(0.35)
  87   ,fInvMassCut(140.0)
  88   ,fChi2NdfCut(999999999.0)
  89   ,fUseAliKFCode(kTRUE)
  90   ,fSharedClusterCut(kFALSE)
  91   ,fRequireITSStandalone(0)
  92   ,fMinNbCls(2)
  93   ,fMinITSChi2(10.0)
  94   ,fMinNbClsSDDSPD(2)
  95   ,fPIDPartner(kFALSE)
  96   ,fPIDMethodPartner(0x0)
  97   ,fPIDMethodPartnerITS(0x0)
  98   ,fDebugLevel(0)
  99   ,fList(0x0)
 100   ,fListPostProcess(0x0)
 101 {
 102   //
 103   // Default constructor
 104   //
 105   for(Int_t k =0; k < 10; k++) {
 106     fCuts[k] = kFALSE;
 107   }
 108
 109 }
 110
 111 //_______________________________________________________________________________________________
 112 AliHFEelecbackground::AliHFEelecbackground(const AliHFEelecbackground &p):
 113   TObject(p)
 114   ,fhtmp(0x0)
 115   ,fhtmpf(0x0)
 116   ,fhtmpp(0x0)
 117   ,fESD1(0x0)
 118   ,fAOD1(0x0)
 119   ,fMCEvent(0x0)
 120   ,fBz(p.fBz)
 121   ,fkVertex(p.fkVertex)
 122   ,fPtESD(p.fPtESD)
 123   ,fIndexTrack(0)
 124   ,fPdg(0)
 125   ,fLabMother(-1)
 126   ,fIsFrom(-1)
 127   ,fMotherGamma(-1)
 128   ,fMotherPi0(-1)
 129   ,fMotherC(-1)
 130   ,fMotherB(-1)
 131   ,fMotherEta(-1)
 132   ,fIsPartner(kFALSE)
 133   ,fIsSplittedTrack(kFALSE)
 134   ,fOpeningAngleCut(0.35)
 135   ,fInvMassCut(140.0)
 136   ,fChi2NdfCut(999999999.0)
 137   ,fUseAliKFCode(kTRUE)
 138   ,fSharedClusterCut(kFALSE)
 139   ,fRequireITSStandalone(0)
 140   ,fMinNbCls(2)
 141   ,fMinITSChi2(10.0)
 142   ,fMinNbClsSDDSPD(2)
 143   ,fPIDPartner(kFALSE)
 144   ,fPIDMethodPartner(0x0)
 145   ,fPIDMethodPartnerITS(0x0)
 146   ,fDebugLevel(0)
 147   ,fList(0x0)
 148   ,fListPostProcess(0x0)
 149 {
 150   //
 151   // Copy constructor
 152   //
 153   for(Int_t k =0; k < 10; k++) {
 154     fCuts[k] = kFALSE;
 155   }
 156 }
 157
 158 //_______________________________________________________________________________________________
 159 AliHFEelecbackground&
 160 AliHFEelecbackground::operator=(const AliHFEelecbackground &)
 161 {
 162   //
 163   // Assignment operator
 164   //
 165
 166   AliInfo("Not yet implemented.");
 167   return *this;
 168 }
 169
 170 //_______________________________________________________________________________________________
 171 AliHFEelecbackground::~AliHFEelecbackground()
 172 {
 173   //
 174   // Destructor
 175   //
 176   if(fPIDMethodPartner) delete fPIDMethodPartner;
 177   if(fPIDMethodPartnerITS) delete fPIDMethodPartnerITS;
 178
 179   if(fList){
 180     fList->Clear();
 181     delete fList;
 182   }
 183
 184   if(fListPostProcess){
 185     fListPostProcess->Clear();
 186     delete fListPostProcess;
 187   }
 188 }
 189 //___________________________________________________________________________________________
 190 Bool_t AliHFEelecbackground::Load(const Char_t * filename)
 191 {
 192   //
 193   // Generic container loader
 194   //
 195
 196   if(!TFile::Open(filename)){
 197     return kFALSE;
 198   }
 199   TList *o = 0x0;
 200   if(!(o = (TList*)gFile->Get("Results"))){
 201     return kFALSE;
 202   }
 203   TList *oe = 0x0;
 204   if(!(oe = (TList*)dynamic_cast<TList *>(o->FindObject("HFEelecbackground")))){
 205     return kFALSE;
 206   }
 207   fList = (TList*)oe->Clone("HFEelecbackground");
 208   gFile->Close();
 209   return kTRUE;
 210 }
 211 //___________________________________________________________________________________________
 212 Bool_t AliHFEelecbackground::Load(TList * const outputlist)
 213 {
 214   //
 215   // Generic container loader
 216   //
 217   if(!outputlist) return kFALSE;
 218   else   fList = (TList*)outputlist->Clone("HFEelecbackground");
 219   return kTRUE;
 220 }
 221 //_______________________________________________________________________________________________
 222 void AliHFEelecbackground::Reset()
 223 {
 224   //
 225   // Reset variables
 226   //
 227   fPtESD = 0.0;
 228   fIndexTrack = -1;
 229   fPdg = -1;
 230   fLabMother = -1;
 231   fIsFrom = -1;
 232   fMotherGamma = -1;
 233   fMotherPi0 = -1;
 234   fMotherC = -1;
 235   fMotherB = -1;
 236   fMotherEta = -1;
 237   fIsPartner = kFALSE;
 238   fIsSplittedTrack = kFALSE;
 239   for(Int_t id = 0; id < 10; id++) {
 240     fCuts[id] = kFALSE;
 241   }
 242
 243 }
 244 //_______________________________________________________________________________________________
 245 void AliHFEelecbackground::CreateHistograms(TList * const qaList)
 246 {
 247   //
 248   // create histograms
 249   //
 250   if(!qaList) return;
 251
 252   fList = qaList;
 253   fList->SetName("HFEelecbackground");
 254
 255   //////////
 256   // bins
 257   /////////
 258
 259   const Int_t nBinsPt = 25;
 260   Double_t minPt = 0.01;
 261   Double_t maxPt = 10.0;
 262
 263   const Int_t nBinsPtMore = 100;
 264   Double_t minPtMore = 0.01;
 265   Double_t maxPtMore = 10.0;
 266
 267   const Int_t nBinsInv = 50;
 268   Double_t minInv = 0.0;
 269   Double_t maxInv = 0.2;
 270
 271   const Int_t nBinsOp = 50;
 272   Double_t minOp = 0.0;
 273   Double_t maxOp = 2;
 274
 275   const Int_t nBinsCh = 4;
 276   Double_t minCh = 0.0;
 277   Double_t maxCh = 4.0;
 278
 279   Double_t binLimLogPt[nBinsPt+1];
 280   Double_t binLimPt[nBinsPt+1];
 281   for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
 282   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
 283
 284   Double_t binLimLogPtMore[nBinsPtMore+1];
 285   Double_t binLimPtMore[nBinsPtMore+1];
 286   for(Int_t i=0; i<=nBinsPtMore; i++) binLimLogPtMore[i]=(Double_t)TMath::Log10(minPtMore) + (TMath::Log10(maxPtMore)-TMath::Log10(minPtMore))/nBinsPtMore*(Double_t)i ;
 287   for(Int_t i=0; i<=nBinsPtMore; i++) binLimPtMore[i]=(Double_t)TMath::Power(10,binLimLogPtMore[i]);
 288
 289   Double_t binLimInv[nBinsInv+1];
 290   for(Int_t i=0; i<=nBinsInv; i++) binLimInv[i]=(Double_t)minInv  + (maxInv-minInv)  /nBinsInv*(Double_t)i ;
 291
 292   Double_t binLimOp[nBinsOp+1];
 293   for(Int_t i=0; i<=nBinsOp; i++) binLimOp[i]=(Double_t)minOp  + (maxOp-minOp) /nBinsOp*(Double_t)i ;
 294
 295   Double_t binLimCh[nBinsCh+1];
 296   for(Int_t i=0; i<=nBinsCh; i++) binLimCh[i]=(Double_t)minCh  + (maxCh-minCh) /nBinsCh*(Double_t)i ;
 297
 298   const Int_t nvarData = 5;
 299   // pt reconstructed tagged e
 300   // pt reconstructed mother
 301   // opening angle
 302   // invariant mass
 303   // Data: charge (0-opposite sign, 1-like sign ++, 2-like sign --, 3-rotated tracks)
 304
 305   Int_t iBinData[nvarData];
 306   iBinData[0]=nBinsPt;
 307   iBinData[1]=nBinsPt;
 308   iBinData[2]=nBinsOp;
 309   iBinData[3]=nBinsInv;
 310   iBinData[4]=nBinsCh;
 311
 312   //
 313   // Opening angle and invariant mass
 314   //
 315
 316   THnSparseF *hsSparseData = new THnSparseF("OpeningangleinvmassData","",nvarData,iBinData);
 317   hsSparseData->SetBinEdges(0,&binLimPt[0]);
 318   hsSparseData->SetBinEdges(1,&binLimPt[0]);
 319   hsSparseData->SetBinEdges(2,&binLimOp[0]);
 320   hsSparseData->SetBinEdges(3,&binLimInv[0]);
 321   hsSparseData->SetBinEdges(4,&binLimCh[0]);
 322   hsSparseData->Sumw2();
 323
 324   fList->AddAt(hsSparseData,kDatai);
 325
 326   //
 327   // Radius, DCA and Chi2Ndf
 328   //
 329
 330   TH1F *dataRadiusHisto = new TH1F("DataRadius","", 200, 0.0, 200.0); // recontructed radius from the AliKF of the e+e- pair
 331   fList->AddAt(dataRadiusHisto,kDatar);
 332
 333   TH1F *dataDcaHisto = new TH1F("DataDCA","", 100, 0.0, 6.0); // dca distribution
 334   fList->AddAt(dataDcaHisto,kDatadca);
 335
 336   TH1F *dataChi2NdfHisto = new TH1F("DataChi2Ndf","", 100, 0.0, 5.0); // chi2Ndf distribution
 337   fList->AddAt(dataChi2NdfHisto,kDatachi2Ndf);
 338
 339
 340   if(HasMCData()) {
 341
 342     //
 343     // Opening angle and invariant mass with MC infos
 344     //
 345
 346     const Int_t nvarMCo = 6;
 347     // pt reconstructed tagged e
 348     // pt reconstructed mother
 349     // opening angle
 350     // invariant mass
 351     // MC: 0-FromBackground, 1-FromGamma, 2-FromPi0, 3-FromEta, 4-FromC, 5-FromB
 352     // MCSplitted: 0-not, 1-splittedOs, 2-ksplittedSs
 353
 354
 355     const Int_t nBinsMCOrigin = 6;
 356     Double_t minMCOrigin = 0.0;
 357     Double_t maxMCOrigin = 6.0;
 358
 359     Double_t binLimMCOrigin[nBinsMCOrigin+1];
 360     for(Int_t i=0; i<=nBinsMCOrigin; i++) binLimMCOrigin[i]=(Double_t)minMCOrigin  + (maxMCOrigin-minMCOrigin) /nBinsMCOrigin*(Double_t)i ;
 361
 362     const Int_t nBinsMCSplitted = 3;
 363     Double_t minMCSplitted = 0.0;
 364     Double_t maxMCSplitted = 3.0;
 365
 366     Double_t binLimMCSplitted[nBinsMCSplitted+1];
 367     for(Int_t i=0; i<=nBinsMCSplitted; i++) binLimMCSplitted[i]=(Double_t)minMCSplitted  + (maxMCSplitted-minMCSplitted) /nBinsMCSplitted*(Double_t)i ;
 368
 369     Int_t iBinMCo[nvarMCo];
 370     iBinMCo[0]=nBinsPt;
 371     iBinMCo[1]=nBinsPt;
 372     iBinMCo[2]=nBinsOp;
 373     iBinMCo[3]=nBinsInv;
 374     iBinMCo[4]=nBinsMCOrigin;
 375     iBinMCo[5]=nBinsMCSplitted;
 376
 377     THnSparseF *hsSparseMCo = new THnSparseF("OpeningangleinvmassMC","",nvarMCo,iBinMCo);
 378     hsSparseMCo->SetBinEdges(0,&binLimPt[0]);
 379     hsSparseMCo->SetBinEdges(1,&binLimPt[0]);
 380     hsSparseMCo->SetBinEdges(2,&binLimOp[0]);
 381     hsSparseMCo->SetBinEdges(3,&binLimInv[0]);
 382     hsSparseMCo->SetBinEdges(4,&binLimMCOrigin[0]);
 383     hsSparseMCo->SetBinEdges(5,&binLimMCSplitted[0]);
 384     hsSparseMCo->Sumw2();
 385
 386     fList->AddAt(hsSparseMCo,kMCo);
 387
 388     //
 389     // Radius, DCA and Chi2Ndf with MC info
 390     //
 391
 392     TH2F *mcRadiusHisto = new TH2F("MCRadius","", 200, 0.0, 200.0,6,-0.5,5.5); // recontructed radius from the AliKF of the e+e- pair
 393     fList->AddAt(mcRadiusHisto,kMCr);
 394
 395     TH2F *mcDcaHisto = new TH2F("MCDCA","", 100, 0.0, 6.0,6,-0.5,5.5); // dca distribution
 396     fList->AddAt(mcDcaHisto,kMCdca);
 397
 398     TH2F *mcChi2NdfHisto = new TH2F("MCChi2Ndf","", 100, 0.0, 5.0,6,-0.5,5.5); // chi2Ndf distribution
 399     fList->AddAt(mcChi2NdfHisto,kMCchi2Ndf);
 400
 401     //////////////////////////////////////////////////////////
 402     // if fDebugLevel 1: Rejection efficiencies of the cuts
 403     //////////////////////////////////////////////////////////
 404
 405     if(fDebugLevel > 0) {
 406
 407       if(HasMCData()) {
 408
 409         const Int_t nvarMCe = 3;
 410         // pt reconstructed tagged e
 411         // cut passed: 0-all, 1-Partner tracked, 2-Opposite-sign, 3-SingleTrackCutPart, 4-ShareCluster, 5-PID, 6-DCA, 7-chi2Ndf AliKF, 8-Openingangle, 9-Invmass
 412         // MC: 0-FromBackground, 1-FromGamma, 2-FromPi0, 3-FromEta, 4-FromC, 5-FromB
 413
 414       const Int_t nBinsMCCutPassed = 10;
 415       Double_t minMCCutPassed = -0.5;
 416       Double_t maxMCCutPassed = 9.5;
 417
 418       Double_t binLimMCCutPassed[nBinsMCCutPassed+1];
 419       for(Int_t i=0; i<=nBinsMCCutPassed; i++) binLimMCCutPassed[i]=(Double_t)minMCCutPassed  + (maxMCCutPassed-minMCCutPassed) /nBinsMCCutPassed*(Double_t)i ;
 420
 421       Int_t iBinMCe[nvarMCe];
 422       iBinMCe[0]=nBinsPt;
 423       iBinMCe[1]=nBinsMCCutPassed;
 424       iBinMCe[2]=nBinsMCOrigin;
 425
 426       THnSparseF *hsSparseMCe = new THnSparseF("CutPassedMC","",nvarMCe,iBinMCe);
 427       hsSparseMCe->SetBinEdges(0,&binLimPt[0]);
 428       hsSparseMCe->SetBinEdges(1,&binLimMCCutPassed[0]);
 429       hsSparseMCe->SetBinEdges(2,&binLimMCOrigin[0]);
 430       hsSparseMCe->Sumw2();
 431
 432       fList->AddAt(hsSparseMCe,kMCe);
 433
 434       }
 435     }
 436
 437     /////////////////////////////////////////////////////////////////
 438     // if fDebugLevel 1: PIDPartCut and ShareClusters
 439     /////////////////////////////////////////////////////////////////
 440
 441     if(fDebugLevel > 1) {
 442
 443       if(!fRequireITSStandalone){
 444
 445         //
 446         // TPC
 447         //
 448
 449         TH2F *tpcPartner0 = new TH2F("TPCPartner0","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0);
 450         fList->AddAt(tpcPartner0,kMCcutPart0);
 451         TH2F *tpcPartner1 = new TH2F("TPCPartner1","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0);
 452         fList->AddAt(tpcPartner1,kMCcutPart1);
 453
 454       }
 455       else {
 456
 457         //
 458         // ITS
 459         //
 460
 461         TH2F *itsPartner0 = new TH2F("ITSPartner0","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0);
 462         fList->AddAt(itsPartner0,kMCcutPart0);
 463         TH2F *itsPartner1 = new TH2F("ITSPartner1","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0);
 464         fList->AddAt(itsPartner1,kMCcutPart1);
 465
 466         /////////////////////////////////////////////////////
 467         // dEdx of the four layers for the track partner
 468         /////////////////////////////////////////////////////
 469         const Int_t nvarITSsignal = 5;
 470
 471         const Int_t nBinsITSsignal  = 100;
 472         Double_t minITSsignal = 0.0;
 473         Double_t maxITSsignal = 350.0;
 474
 475         Double_t binLimITSsignal[nBinsITSsignal+1];
 476         for(Int_t i=0; i<=nBinsITSsignal; i++) binLimITSsignal[i]=(Double_t)minITSsignal  + (maxITSsignal-minITSsignal) /nBinsITSsignal*(Double_t)i ;
 477
 478         Int_t iBinITSsignal[nvarITSsignal];
 479         iBinITSsignal[0]=nBinsPt;
 480         iBinITSsignal[1]=nBinsITSsignal;
 481         iBinITSsignal[2]=nBinsITSsignal;
 482         iBinITSsignal[3]=nBinsITSsignal;
 483         iBinITSsignal[4]=nBinsITSsignal;
 484
 485         THnSparseF *hsSparseITSpid = new THnSparseF("SparseITSsignal","",nvarITSsignal,iBinITSsignal);
 486         hsSparseITSpid->SetBinEdges(0,&binLimPt[0]);
 487         hsSparseITSpid->SetBinEdges(1,&binLimITSsignal[0]);
 488         hsSparseITSpid->SetBinEdges(2,&binLimITSsignal[0]);
 489         hsSparseITSpid->SetBinEdges(3,&binLimITSsignal[0]);
 490         hsSparseITSpid->SetBinEdges(4,&binLimITSsignal[0]);
 491         hsSparseITSpid->Sumw2();
 492
 493         fList->AddAt(hsSparseITSpid,kMCcutPart2);
 494
 495         ///////////////////////////////////////////////////////////////////////////////////////
 496         // dEdx of the four layers for the track partner and track to reject splitted track
 497         ///////////////////////////////////////////////////////////////////////////////////////
 498         const Int_t nvarITSsignalSplit = 5;
 499
 500         const Int_t nBinsITSSplit  = 2;
 501         Double_t minITSSplit = 0.0;
 502         Double_t maxITSSplit = 2.0;
 503
 504         Double_t binLimITSSplit[nBinsITSSplit+1];
 505         for(Int_t i=0; i<=nBinsITSSplit; i++) binLimITSSplit[i]=(Double_t)minITSSplit  + (maxITSSplit-minITSSplit) /nBinsITSSplit*(Double_t)i ;
 506
 507
 508         const Int_t nBinsITSsignalSplit  = 50;
 509         Double_t minITSsignalSplit = -25.0;
 510         Double_t maxITSsignalSplit = 25.0;
 511
 512         Double_t binLimITSsignalSplit[nBinsITSsignalSplit+1];
 513         for(Int_t i=0; i<=nBinsITSsignalSplit; i++) binLimITSsignalSplit[i]=(Double_t)minITSsignalSplit  + (maxITSsignalSplit-minITSsignalSplit) /nBinsITSsignalSplit*(Double_t)i ;
 514
 515         Int_t iBinITSsignalSplit[nvarITSsignalSplit];
 516         iBinITSsignalSplit[0]=nBinsITSSplit;
 517         for(Int_t k = 1; k < 5; k++){
 518           iBinITSsignalSplit[k]=nBinsITSsignalSplit;
 519         }
 520
 521         THnSparseF *hsSparseITSpidSplit = new THnSparseF("SparseITSsignalSplit","",nvarITSsignalSplit,iBinITSsignalSplit);
 522         hsSparseITSpidSplit->SetBinEdges(0,&binLimITSSplit[0]);
 523         for(Int_t k = 1; k < 5; k++) {
 524           hsSparseITSpidSplit->SetBinEdges(k,&binLimITSsignalSplit[0]);
 525         }
 526         hsSparseITSpidSplit->Sumw2();
 527
 528         fList->AddAt(hsSparseITSpidSplit,kMCcutPart3);
 529
 530       }
 531
 532     }
 533
 534   }
 535
 536   //qaList->Add(fList);
 537
 538 }
 539 //_______________________________________________________________________________________________
 540 void AliHFEelecbackground::PairAnalysis(AliESDtrack* const track, AliESDtrack* const trackPart)
 541 {
 542   //
 543   // calculate (tagged e-partner) dca, opening angle, invariant mass
 544   //
 545
 546   /////////////////////
 547   // pt tagged
 548   //////////////////////
 549   TVector3 v3Dtagged;
 550   Double_t pxyz[3];
 551   track->PxPyPz(&pxyz[0]);
 552   v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
 553   fPtESD = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]);
 554
 555   ////////////////////////
 556   // Take label
 557   ////////////////////////
 558   Int_t indexTrack = TMath::Abs(track->GetLabel());
 559   Int_t indexTrackPart = TMath::Abs(trackPart->GetLabel());
 560
 561   /////////////////////////
 562   // If MC data
 563   ////////////////////////
 564
 565   if(HasMCData()) {
 566
 567     // Take info track if not already done
 568     if(indexTrack!= fIndexTrack) {
 569
 570       for(Int_t id = 0; id < 10; id++) {
 571         fCuts[id] = kFALSE;
 572       }
 573
 574       fIsFrom = kElectronFromBackground;
 575
 576       fPdg = GetPdg(indexTrack);
 577       fLabMother = GetLabMother(indexTrack);
 578
 579       fMotherGamma = IsMotherGamma(indexTrack);
 580       if((fMotherGamma != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromGamma;
 581       fMotherPi0 = IsMotherPi0(indexTrack);
 582       if((fMotherPi0 != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromPi0;
 583       fMotherC = IsMotherC(indexTrack);
 584       if((fMotherC != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromC;
 585       fMotherB = IsMotherB(indexTrack);
 586       if((fMotherB != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromB;
 587       fMotherEta = IsMotherEta(indexTrack);
 588       if((fMotherEta != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromEta;
 589
 590       fIndexTrack = indexTrack;
 591
 592     }
 593
 594     // MC PID for tagged
 595     if(fgUseMCPID) {
 596       if(TMath::Abs(fPdg) != 11) return;
 597     }
 598
 599     // Look at trackPart
 600     fIsPartner = kFALSE;
 601     Int_t pdgPart = GetPdg(indexTrackPart);
 602     if(TMath::Abs(pdgPart) == 11) {
 603       Int_t labMotherPart = GetLabMother(indexTrackPart);
 604       if((labMotherPart == fLabMother) && (indexTrack != indexTrackPart) && (TMath::Abs(fPdg) == 11) && (fPdg*pdgPart < 0) && (fLabMother >=0) && (fLabMother < (((AliStack *)fMCEvent->Stack())->GetNtrack()))) fIsPartner = kTRUE;
 605       // special case of c and b
 606       Int_t motherCPart = IsMotherC(indexTrackPart);
 607       if((motherCPart != -1) && (fIsFrom == kElectronFromC) && (fPdg*pdgPart < 0)){
 608         fIsPartner = kTRUE;
 609       }
 610       Int_t motherBPart = IsMotherB(indexTrackPart);
 611       if((motherBPart != -1) && (fIsFrom == kElectronFromB) && (fPdg*pdgPart < 0)){
 612         fIsPartner = kTRUE;
 613       }
 614     }
 615
 616     // Look at splitted tracks
 617     fIsSplittedTrack = kFALSE;
 618     if(indexTrackPart == fIndexTrack) fIsSplittedTrack = kTRUE;
 619
 620   }
 621
 622   //////////////////////
 623   // Sign
 624   /////////////////////
 625   Int_t sign = -1;
 626   if((track->Charge() > 0.0) && (trackPart->Charge() > 0.0)) sign = kPp;
 627   if((track->Charge() < 0.0) && (trackPart->Charge() < 0.0)) sign = kNn;
 628   if(((track->Charge() > 0.0) && (trackPart->Charge() < 0.0)) || ((track->Charge() < 0.0) && (trackPart->Charge() > 0.0))) sign = kOs;
 629
 630   /////////////////////////
 631   // Cut effects
 632   ////////////////////////
 633   Double_t cuteffect[3];
 634
 635   if(fDebugLevel > 0) {
 636     if(HasMCData()) {
 637       cuteffect[0] = fPtESD;
 638       cuteffect[1] = 0.0;
 639       cuteffect[2] = fIsFrom;
 640       if(!fCuts[0]){
 641         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 642         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 643         fCuts[0] = kTRUE;
 644       }
 645     }
 646   }
 647
 648
 649   ///////////////////////////////
 650   // Cut effect: Partner track
 651   ///////////////////////////////
 652
 653   if(fDebugLevel > 0) {
 654     if(HasMCData() && fIsPartner) {
 655       cuteffect[1] = 1.0;
 656       if(!fCuts[1]) {
 657         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 658         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 659         fCuts[1] = kTRUE;
 660       }
 661     }
 662   }
 663
 664   ///////////////////////////////
 665   // Cut effect: Opposite sign
 666   ///////////////////////////////
 667
 668   if(fDebugLevel > 0) {
 669     if(HasMCData() && fIsPartner && (sign == kOs)) {
 670       cuteffect[1] = 2.0;
 671       if(!fCuts[2]) {
 672         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 673         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 674         fCuts[2] = kTRUE;
 675       }
 676     }
 677   }
 678
 679   ////////////////////////
 680   // Partner track cut
 681   ////////////////////////
 682   if(!SingleTrackCut(trackPart)) return;
 683
 684   if(fDebugLevel > 0) {
 685     if(HasMCData() && fIsPartner && (sign==kOs)) {
 686       cuteffect[1] = 3.0;
 687       if(!fCuts[3]) {
 688         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 689         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 690         fCuts[3] = kTRUE;
 691       }
 692     }
 693   }
 694
 695   /////////////////////////
 696   // shared clusters cut
 697   /////////////////////////
 698   if(fSharedClusterCut && ShareCluster(track,trackPart)) return;
 699
 700   if(fDebugLevel > 0) {
 701     if(HasMCData() && fIsPartner && (sign==kOs)) {
 702       cuteffect[1] = 4.0;
 703       if(!fCuts[4]) {
 704         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 705         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 706         fCuts[4] = kTRUE;
 707       }
 708     }
 709   }
 710
 711   ////////////////////////
 712   // PID Partner track
 713   ////////////////////////
 714   if(!PIDTrackCut(trackPart)) return;
 715
 716   if(fDebugLevel > 0) {
 717     if(HasMCData() && fIsPartner && (sign==kOs)) {
 718       cuteffect[1] = 5.0;
 719       if(!fCuts[5]) {
 720         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 721         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 722         fCuts[5] = kTRUE;
 723       }
 724     }
 725   }
 726
 727
 728   //////////////////////
 729   // DCA
 730   /////////////////////
 731
 732   Double_t xthis,xp;
 733   Double_t dca = track->GetDCA(trackPart,fBz,xthis,xp);
 734   if((fhtmpp = dynamic_cast<TH1F *>(fList->At(kDatadca)))) fhtmpp->Fill(dca);
 735   if(HasMCData()) {
 736     //printf("has MC data for DCA\n");
 737     //printf("IsPartner %d and isfrom %d\n",fIsPartner,fIsFrom);
 738     if(fIsFrom==kElectronFromBackground) {
 739       if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCdca)))) fhtmpf->Fill(dca,fIsFrom);
 740     }
 741     else {
 742       if(fIsPartner){
 743         if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCdca)))) fhtmpf->Fill(dca,fIsFrom);
 744       }
 745     }
 746   }
 747
 748   if(TMath::Abs(dca) > 3.0) return;
 749
 750   if(fDebugLevel > 0) {
 751     if(HasMCData() && fIsPartner && (sign==kOs)) {
 752       cuteffect[1] = 6.0;
 753       if(!fCuts[6]) {
 754         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 755         //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 756         fCuts[6] = kTRUE;
 757       }
 758     }
 759   }
 760
 761   ///////////////////////////////////
 762   // Calcul mother variables
 763   ///////////////////////////////////
 764   Double_t results[5];
 765   Double_t resultsr[5];
 766
 767
 768   if(!fUseAliKFCode) {
 769
 770     /////////////////////////////
 771     // Propagate
 772     ////////////////////////////
 773
 774     Double_t norradius = TMath::Sqrt(fkVertex->GetX()*fkVertex->GetX()+fkVertex->GetY()*fkVertex->GetY());
 775
 776     AliESDtrack *trackCopy = new AliESDtrack(*track);
 777     AliESDtrack *trackPartCopy = new AliESDtrack(*trackPart);
 778     Bool_t propagateok = kTRUE;
 779     if((!(trackPartCopy->PropagateTo(norradius,fBz))) || (!(trackCopy->PropagateTo(norradius,fBz)))) propagateok = kFALSE;
 780     if(!propagateok) {
 781       if(trackCopy) delete trackCopy;
 782       if(trackPartCopy) delete trackPartCopy;
 783       return;
 784     }
 785
 786     CalculateMotherVariable(trackCopy,trackPartCopy,&results[0]);
 787     CalculateMotherVariableR(trackCopy,trackPartCopy,&resultsr[0]);
 788
 789     if(trackCopy) delete trackCopy;
 790     if(trackPartCopy) delete trackPartCopy;
 791
 792   }
 793   else {
 794
 795     if(!CalculateMotherVariable(track,trackPart,&results[0])) return;
 796     if(fDebugLevel > 0) {
 797       if(HasMCData() && fIsPartner && (sign==kOs)) {
 798         cuteffect[1] = 7.0;
 799         if(!fCuts[7]) {
 800           if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 801           //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 802           fCuts[7] = kTRUE;
 803         }
 804       }
 805     }
 806
 807   }
 808
 809   /////////////////////////////////////
 810   // Fill
 811   /////////////////////////////////////
 812
 813   FillOutput(results, resultsr, sign);
 814
 815   if(fDebugLevel > 0) {
 816     if(HasMCData() && fIsPartner && (sign==kOs)) {
 817
 818       if(TMath::Abs(results[4]) < fOpeningAngleCut) {
 819
 820         cuteffect[1] = 8.0;
 821         if(!fCuts[8]) {
 822           if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 823           //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 824           fCuts[8] = kTRUE;
 825         }
 826         if(TMath::Abs(results[1]) < fInvMassCut) {
 827           cuteffect[1] = 9.0;
 828           if(!fCuts[9]) {
 829             if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCe)))) fhtmp->Fill(cuteffect);
 830             //if(fList->At(kMCe)) (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
 831             fCuts[9] = kTRUE;
 832           }
 833         }
 834       }
 835     }
 836   }
 837
 838
 839 }
 840 //_____________________________________________________________________________________
 841 Bool_t AliHFEelecbackground::CalculateMotherVariable(AliESDtrack* const track, AliESDtrack* const trackpart, Double_t *results)
 842 {
 843   //
 844   // variables mother and take the pt of the track
 845   //
 846   // results contain: ptmother, invmass, etamother, phimother, openingangle
 847   //
 848   // with a chi2Ndf cut for AliKF code
 849   //
 850
 851   if(!fUseAliKFCode) {
 852
 853     TVector3 v3Dtagged;
 854     TVector3 v3Dpart;
 855
 856     Double_t pxyz[3];
 857     track->PxPyPz(&pxyz[0]);
 858     v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
 859
 860     Double_t pxyzpart[3];
 861     trackpart->PxPyPz(&pxyzpart[0]);
 862     v3Dpart.SetXYZ(pxyzpart[0],pxyzpart[1],pxyzpart[2]);
 863
 864
 865     TVector3 motherrec = v3Dtagged + v3Dpart;
 866
 867     Double_t etaESDmother = motherrec.Eta();
 868     Double_t ptESDmother  = motherrec.Pt();
 869     Double_t phiESDmother = motherrec.Phi();
 870     if(phiESDmother > TMath::Pi()) phiESDmother = phiESDmother - (2*TMath::Pi());
 871
 872
 873     // openinganglepropagated
 874     Double_t openingangle = v3Dtagged.Angle(v3Dpart);
 875
 876     // invmass
 877     Double_t pESD      = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]+pxyz[2]*pxyz[2]);
 878     Double_t pESDpart  = TMath::Sqrt(pxyzpart[0]*pxyzpart[0]+pxyzpart[1]*pxyzpart[1]+pxyzpart[2]*pxyzpart[2]);
 879
 880     // e propagate
 881     Double_t eESD     = TMath::Sqrt(pESD*pESD+fgkMe*fgkMe);
 882     Double_t eESDpart = TMath::Sqrt(pESDpart*pESDpart+fgkMe*fgkMe);
 883
 884     Double_t invmass = TMath::Sqrt((eESD+eESDpart)*(eESD+eESDpart)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
 885
 886     if(!results) return kFALSE;
 887
 888     results[0] = ptESDmother;
 889     results[1] = etaESDmother;
 890     results[2] = phiESDmother;
 891     results[3] = invmass;
 892     results[4] = openingangle;
 893
 894     return kTRUE;
 895
 896   }
 897   else {
 898
 899     AliKFParticle pair;
 900     pair.Initialize();
 901
 902     // pid
 903     Int_t pid1 = -11;
 904     if(track->Charge() > 0.0) pid1 = 11;
 905     Int_t pid2 = -11;
 906     if(trackpart->Charge() > 0.0) pid2 = 11;
 907
 908
 909     // daughters
 910     AliKFParticle kf(*track,pid1);
 911     AliKFParticle kfpart(*trackpart,pid2);
 912
 913     pair.AddDaughter(kf);
 914     pair.AddDaughter(kfpart);
 915
 916     // variables
 917     Double_t openingangle = kf.GetAngle(kfpart);
 918     Double_t chi2ndf = pair.GetChi2()/pair.GetNDF();
 919     //Double_t decayLength = pair.GetDecayLength();
 920     Double_t radius = pair.GetR();
 921     //Double_t masserror = pair.GetErrMass()>0?pair.GetErrMass()/pair.GetMass():1000000;
 922     Double_t ptpair = pair.GetPt();
 923     Double_t etapair = pair.GetEta();
 924     Double_t phipair = pair.GetPhi();
 925     Double_t masspair = pair.GetMass();
 926
 927     // Put them
 928     if(!results) return kFALSE;
 929
 930     results[0] = ptpair;
 931     results[1] = etapair;
 932     results[2] = phipair;
 933     results[3] = masspair;
 934     results[4] = openingangle;
 935
 936     // chi2Ndf cut
 937     if((fhtmpp = dynamic_cast<TH1F *>(fList->At(kDatachi2Ndf)))) fhtmpp->Fill(chi2ndf);
 938     //if(fList->At(kDatachi2Ndf)) (dynamic_cast<TH1F *>(fList->At(kDatachi2Ndf)))->Fill(chi2ndf);
 939     if(HasMCData()){
 940       if(fIsFrom==kElectronFromBackground) {
 941         if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCchi2Ndf)))) fhtmpf->Fill(chi2ndf,fIsFrom);
 942       }
 943       else {
 944         if(fIsPartner){
 945           if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCchi2Ndf)))) fhtmpf->Fill(chi2ndf,fIsFrom);
 946         }
 947       }
 948     }
 949     if(chi2ndf > fChi2NdfCut) return kFALSE;
 950     else {
 951       if((fhtmpp = dynamic_cast<TH1F *>(fList->At(kDatar)))) fhtmpp->Fill(radius);
 952       //if(fList->At(kDatar)) (dynamic_cast<TH1F *>(fList->At(kDatar)))->Fill(radius);
 953       if(HasMCData()) {
 954         if(fIsFrom==kElectronFromBackground) {
 955           if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCr)))) fhtmpf->Fill(radius,fIsFrom);
 956           //if(fList->At(kMCr))) (dynamic_cast<TH2F *>(fList->At(kMCr)))->Fill(radius,fIsFrom);
 957         }
 958         else {
 959           if(fIsPartner) {
 960             if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCr)))) fhtmpf->Fill(radius,fIsFrom);
 961           }
 962         }
 963       }
 964       return kTRUE;
 965     }
 966
 967   }
 968
 969 }
 970 //_____________________________________________________________________________________
 971 void AliHFEelecbackground::CalculateMotherVariableR(AliESDtrack* const track, AliESDtrack* const trackpart, Double_t *results)
 972 {
 973   //
 974   // variables mother
 975   //
 976   // results contain: ptmother, invmass, etamother, phimother, openingangle
 977   // Implemented only for no AliKF
 978   //
 979
 980   if(!fUseAliKFCode) {
 981
 982     TVector3 v3Dtagged;
 983     TVector3 v3Dpart;
 984
 985     Double_t pxyz[3];
 986     track->PxPyPz(&pxyz[0]);
 987     v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
 988     Double_t pxyzpart[3];
 989     trackpart->PxPyPz(&pxyzpart[0]);
 990     v3Dpart.SetXYZ(pxyzpart[0],pxyzpart[1],pxyzpart[2]);
 991
 992     // rotate the partner
 993     v3Dpart.RotateZ(TMath::Pi());
 994     v3Dpart.GetXYZ(pxyzpart);
 995
 996
 997     TVector3 motherrec = v3Dtagged + v3Dpart;
 998
 999     Double_t etaESDmother = motherrec.Eta();
1000     Double_t ptESDmother  = motherrec.Pt();
1001     Double_t phiESDmother = motherrec.Phi();
1002     if(phiESDmother > TMath::Pi()) phiESDmother = phiESDmother - (2*TMath::Pi());
1003
1004
1005     // openinganglepropagated
1006     Double_t openingangle = v3Dtagged.Angle(v3Dpart);
1007     //printf("Openingangle %f\n",openingangle);
1008
1009     // invmass
1010     Double_t pESD      = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]+pxyz[2]*pxyz[2]);
1011     Double_t pESDpart  = TMath::Sqrt(pxyzpart[0]*pxyzpart[0]+pxyzpart[1]*pxyzpart[1]+pxyzpart[2]*pxyzpart[2]);
1012     // e propagate
1013     Double_t eESD     = TMath::Sqrt(pESD*pESD+fgkMe*fgkMe);
1014     Double_t eESDpart = TMath::Sqrt(pESDpart*pESDpart+fgkMe*fgkMe);
1015
1016     Double_t invmass = TMath::Sqrt((eESD+eESDpart)*(eESD+eESDpart)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
1017
1018     if(!results) return;
1019
1020     results[0] = ptESDmother;
1021     results[1] = etaESDmother;
1022     results[2] = phiESDmother;
1023     results[3] = invmass;
1024     results[4] = openingangle;
1025
1026   }
1027
1028 }
1029 //_________________________________________________________________________________
1030 void AliHFEelecbackground::FillOutput(Double_t *results, Double_t *resultsr, Int_t sign)
1031 {
1032   //
1033   // Fill the Data and MC THnSparseF
1034   //
1035
1036   if((!results) || (!resultsr)) return;
1037
1038   Double_t co[6];
1039   co[0] = fPtESD;
1040   co[1] = results[0];
1041   co[2] = TMath::Abs(results[4]);
1042   co[3] = results[3];
1043   co[4] = sign;
1044   co[5] = 0.0;
1045
1046   if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kDatai)))) fhtmp->Fill(co);
1047   //if(fList->At(kDatai))(dynamic_cast<THnSparseF *>(fList->At(kDatai)))->Fill(co);
1048
1049   if((sign==kOs) && (!fUseAliKFCode)) {
1050
1051     co[1] = resultsr[0];
1052     co[2] = TMath::Abs(resultsr[4]);
1053     co[3] = resultsr[3];
1054     co[4] = kR;
1055     co[5] = 0.0;
1056
1057     if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kDatai)))) fhtmp->Fill(co);
1058     //if(fList->At(kDatai)) (dynamic_cast<THnSparseF *>(fList->At(kDatai)))->Fill(co);
1059
1060   }
1061
1062   if(HasMCData()){
1063
1064     // Reset
1065     co[1] = results[0];
1066     co[2] = TMath::Abs(results[4]);
1067     co[3] = results[3];
1068
1069     // Origin
1070     co[4] = kElectronFromBackground;
1071     if((sign==kOs) && fIsPartner) co[4] = fIsFrom;
1072
1073     // Splitted tracks
1074     co[5] = kNotSplitted;
1075     if(fIsSplittedTrack) {
1076       if(sign==kOs){
1077         co[5] = kSplittedOs;
1078       }
1079       else {
1080         co[5] = kSplittedSs;
1081       }
1082     }
1083
1084     if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCo)))) fhtmp->Fill(co);
1085     //if(fList->At(kMCo)) (dynamic_cast<THnSparseF *>(fList->At(kMCo)))->Fill(co);
1086
1087   }
1088
1089 }
1090 //_______________________________________________________________________________________________
1091 Bool_t AliHFEelecbackground::SingleTrackCut(AliESDtrack* const trackPart) const
1092 {
1093   //
1094   // Return minimum quality for the partner
1095   //
1096
1097   if(trackPart->GetKinkIndex(0)>0) return kFALSE;
1098
1099
1100   UInt_t status = trackPart->GetStatus();
1101
1102   if(fRequireITSStandalone > 0) {
1103
1104     /////////////////////
1105     // ITS Standalone
1106     ////////////////////
1107
1108     if(fRequireITSStandalone==1) {
1109       if(((status & AliESDtrack::kITSin) == 0 || (trackPart->IsPureITSStandalone()) || ((status&AliESDtrack::kITSrefit)==0))) return kFALSE;
1110     }
1111
1112     if(fRequireITSStandalone==2) {
1113       if(!trackPart->IsPureITSStandalone() || ((status&AliESDtrack::kITSrefit)==0)) return kFALSE;
1114     }
1115
1116     // Chi2
1117     Double_t chi2 = trackPart->GetITSchi2();
1118     if(chi2 > fMinITSChi2) return kFALSE;
1119
1120     // Min Nb of clusters
1121     Int_t nbcl = trackPart->GetITSclusters(0);
1122     if(nbcl < fMinNbCls)  return kFALSE;
1123
1124     // Min Nb of points in SDD and SPD
1125     Int_t nbSDDSPD = 0;
1126     for(Int_t layer = 0; layer < 4; layer++){
1127       if(trackPart->HasPointOnITSLayer(layer)) nbSDDSPD++;
1128     }
1129     if(nbSDDSPD < fMinNbClsSDDSPD) return kFALSE;
1130
1131
1132   }
1133   else {
1134
1135     /////////
1136     // TPC
1137     /////////
1138
1139     if((status&AliESDtrack::kTPCrefit)==0) return kFALSE;
1140
1141     // Min Nb of clusters
1142     Int_t nbcl = trackPart->GetTPCclusters(0);
1143     if(nbcl < fMinNbCls) return kFALSE;
1144
1145   }
1146
1147   return kTRUE;
1148
1149 }
1150 //_______________________________________________________________________________________________
1151 Bool_t AliHFEelecbackground::PIDTrackCut(AliESDtrack* const trackPart)
1152 {
1153   //
1154   // PID for the partner using TPC or ITS
1155   //
1156
1157   if(fRequireITSStandalone > 0) {
1158
1159     /////////////////////
1160     // ITS Standalone
1161     ////////////////////
1162
1163     // signal
1164     Double_t itsSignal = trackPart->GetITSsignal();
1165     Double_t p = trackPart->P();
1166
1167     if(fDebugLevel > 1) {
1168       if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCcutPart0)))) fhtmpf->Fill(p,itsSignal);
1169       //if(fList->At(kMCcutPart0)) (dynamic_cast<TH2F *>(fList->At(kMCcutPart0)))->Fill(p,itsSignal);
1170     }
1171
1172     ///////////
1173     // PID
1174     //////////
1175     if(fPIDPartner) {
1176
1177       // Take signal trackPart
1178       Double_t dEdxSamplesPart[4];
1179       trackPart->GetITSdEdxSamples(dEdxSamplesPart);
1180
1181       // Cut at 2 sigma
1182       if(!fPIDMethodPartnerITS) fPIDMethodPartnerITS = new AliESDpid;
1183       Float_t nsigma = fPIDMethodPartnerITS->NumberOfSigmasITS(trackPart, AliPID::kElectron);
1184       if(TMath::Abs(nsigma) > 2.0) return kFALSE;
1185
1186       // fill signal
1187       if(fDebugLevel > 1) {
1188
1189         Double_t entries[5];
1190         entries[0] = p;
1191         entries[1] = dEdxSamplesPart[0];
1192         entries[2] = dEdxSamplesPart[1];
1193         entries[3] = dEdxSamplesPart[2];
1194         entries[4] = dEdxSamplesPart[3];
1195
1196         //if(fList->At(kMCcutPart1)) (dynamic_cast<TH2F *>(fList->At(kMCcutPart1)))->Fill(p,itsSignal);
1197         if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCcutPart1)))) fhtmpf->Fill(p,itsSignal);
1198         if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCcutPart2)))) fhtmp->Fill(entries);
1199         //if(fList->At(kMCcutPart2)) (dynamic_cast<THnSparseF *>(fList->At(kMCcutPart2)))->Fill(entries);
1200
1201       }
1202
1203     }
1204
1205   }
1206   else {
1207
1208     /////////
1209     // TPC
1210     /////////
1211
1212     Double_t tpcSignal = trackPart->GetTPCsignal();
1213     Double_t p = trackPart->GetInnerParam() ? trackPart->GetInnerParam()->P() : trackPart->P();
1214
1215     if(fDebugLevel > 1) {
1216       //printf("tpcSignal %f\n",tpcSignal);
1217       //if(fList->At(kMCcutPart0)) (dynamic_cast<TH2F *>(fList->At(kMCcutPart0)))->Fill(p,tpcSignal);
1218       if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCcutPart0)))) fhtmpf->Fill(p,tpcSignal);
1219     }
1220
1221     // PID
1222     if(fPIDPartner) {
1223       if(!fPIDMethodPartner) return kFALSE;
1224       AliHFEpidObject hfetrack;
1225       hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis;
1226       hfetrack.fRecTrack = trackPart;
1227       //if(HasMCData()) hfetrack.fMCtrack = mctrack;
1228       if(!fPIDMethodPartner->IsSelected(&hfetrack)) return kFALSE;
1229
1230       if(fDebugLevel > 1) {
1231         if((fhtmpf = dynamic_cast<TH2F *>(fList->At(kMCcutPart1)))) fhtmpf->Fill(p,tpcSignal);
1232         //if(fList->At(kMCcutPart1)) (dynamic_cast<TH2F *>(fList->At(kMCcutPart1)))->Fill(p,tpcSignal);
1233       }
1234     }
1235
1236   }
1237
1238   return kTRUE;
1239
1240 }
1241 //__________________________________________________________________________________________
1242 Bool_t AliHFEelecbackground::ShareCluster(AliESDtrack * const track1,AliESDtrack * const track2)
1243 {
1244   //
1245   // Look if the two tracks shared clusters in the TPC or in the ITS depending on the method
1246   //
1247   // For TPC:
1248   // hsfval: number of shared clusters
1249   // hsmval: quality of the tracks
1250   //
1251   // For ITS:
1252   // compare the dEdx in the ITS
1253   //
1254
1255   if(!fRequireITSStandalone) {
1256
1257     //////////
1258     // TPC
1259     //////////
1260
1261
1262     Int_t nh = 0;
1263     Int_t an = 0;
1264     Int_t ns = 0;
1265     Float_t hsmval = 0.0;
1266     Float_t hsfval = 0.0;
1267
1268     for(unsigned int imap=0;imap<track1->GetTPCClusterMap().GetNbits(); imap++) {
1269       if(track1->GetTPCClusterMap().TestBitNumber(imap) &&
1270          track2->GetTPCClusterMap().TestBitNumber(imap)) {
1271         // Do they share it ?
1272         if (track1->GetTPCSharedMap().TestBitNumber(imap) &&
1273             track2->GetTPCSharedMap().TestBitNumber(imap))
1274           {
1275             an++;
1276             nh+=2;
1277             ns+=2;
1278           }
1279         else {
1280           an--;
1281           nh+=2;
1282         }
1283       }
1284       else if (track1->GetTPCClusterMap().TestBitNumber(imap) ||
1285                track2->GetTPCClusterMap().TestBitNumber(imap)) {
1286         an++;
1287         nh++;
1288       }
1289     }
1290
1291     if (nh >0) {
1292       hsmval = an*1.0/nh;
1293       hsfval = ns*1.0/nh;
1294     }
1295
1296
1297     if((hsfval > 0.15) || (hsmval > -0.15)) return kTRUE; //they share cluster
1298     else return kFALSE;
1299
1300
1301   }
1302   else {
1303
1304
1305     //////////
1306     // ITS
1307     /////////
1308
1309     // Take signals
1310     Double_t dEdxSamples1[4];
1311     track1->GetITSdEdxSamples(dEdxSamples1);
1312     Double_t dEdxSamples2[4];
1313     track2->GetITSdEdxSamples(dEdxSamples2);
1314
1315     // If there are matching
1316     Int_t nbClusters = 0;
1317     Bool_t match[4] = {kTRUE,kTRUE,kTRUE,kTRUE};
1318     Double_t limit[4] = {1.5,1.5,1.5,1.5};
1319     for(Int_t layer = 0; layer < 4; layer++) {
1320       if(track1->HasPointOnITSLayer(layer+2) && track2->HasPointOnITSLayer(layer+2)) {
1321         if(TMath::Abs(dEdxSamples1[layer]-dEdxSamples2[layer])>limit[layer]) match[layer] = kFALSE;
1322         nbClusters++;
1323       }
1324     }
1325     //printf("nbClusters %d\n",nbClusters);
1326
1327     // fill signal
1328     if(fDebugLevel > 1) {
1329       Double_t entriesSplit[5];
1330       entriesSplit[0] = 0.0;
1331       if(fIsSplittedTrack) entriesSplit[0] = 1.0;
1332
1333       for(Int_t layer = 0; layer < 4; layer++) {
1334         if(track1->HasPointOnITSLayer(layer+2) && track2->HasPointOnITSLayer(layer+2)) {
1335           entriesSplit[layer+1] = dEdxSamples1[layer]-dEdxSamples2[layer];
1336         }
1337         else entriesSplit[layer+1] = -100.0;
1338       }
1339       if((fhtmp = dynamic_cast<THnSparseF *>(fList->At(kMCcutPart3)))) fhtmp->Fill(entriesSplit);
1340       //if(fList->At(kMCcutPart3)) (dynamic_cast<THnSparseF *>(fList->At(kMCcutPart3)))->Fill(entriesSplit);
1341     }
1342
1343     // Return
1344     Int_t nbClustersNotClose = 0;
1345     for(Int_t layer = 0; layer < 4; layer++) {
1346       if(!match[layer]) nbClustersNotClose++;
1347     }
1348     if((nbClusters > 1) && (nbClustersNotClose > 0.75*nbClusters)) return kFALSE;
1349     else return kTRUE;
1350
1351   }
1352
1353 }
1354 //____________________________________________________________________________________________________________
1355 void AliHFEelecbackground::SetPIDPartner() {
1356
1357   //
1358   // Init the stuff for PID on the partner track
1359   //
1360
1361   fPIDPartner = kTRUE;
1362
1363   if(fRequireITSStandalone == 0) {
1364
1365     if(!fPIDMethodPartner) {
1366       fPIDMethodPartner = new AliHFEpid();
1367       fPIDMethodPartner->InitializePID("Strategy1");     // 3 sigma cut in TPC
1368     }
1369
1370   }
1371   else {
1372
1373     if(!fPIDMethodPartnerITS) fPIDMethodPartnerITS = new AliESDpid;
1374
1375   }
1376
1377 }
1378 //______________________________________________________________________________________________
1379 void AliHFEelecbackground::SetEvent(AliESDEvent* const ESD)
1380 {
1381   //
1382   // Set the AliESD Event, the magnetic field and the primary vertex
1383   //
1384
1385   fESD1=ESD;
1386   fBz=fESD1->GetMagneticField();
1387   fkVertex = fESD1->GetPrimaryVertex();
1388
1389 }
1390 //____________________________________________________________________________________________________________
1391 Int_t AliHFEelecbackground::IsMotherGamma(Int_t tr) {
1392
1393   //
1394   // Return the lab of gamma mother or -1 if not gamma
1395   //
1396
1397   AliStack* stack = fMCEvent->Stack();
1398   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1399
1400   // Take mother
1401   TParticle * particle = stack->Particle(tr);
1402   if(!particle) return -1;
1403   Int_t imother   = particle->GetFirstMother();
1404   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
1405   TParticle * mother = stack->Particle(imother);
1406   if(!mother) return -1;
1407
1408   // Check gamma
1409   Int_t pdg = mother->GetPdgCode();
1410   if(TMath::Abs(pdg) == 22) return imother;
1411   if(TMath::Abs(pdg) == 11) {
1412     return IsMotherGamma(imother);
1413   }
1414   return -1;
1415
1416 }
1417 //
1418 //____________________________________________________________________________________________________________
1419 Int_t AliHFEelecbackground::IsMotherPi0(Int_t tr) {
1420
1421   //
1422   // Return the lab of pi0 mother or -1 if not pi0
1423   //
1424
1425   AliStack* stack = fMCEvent->Stack();
1426   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1427
1428   // Take mother
1429   TParticle * particle = stack->Particle(tr);
1430   if(!particle) return -1;
1431   Int_t imother   = particle->GetFirstMother();
1432   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
1433   TParticle * mother = stack->Particle(imother);
1434   if(!mother) return -1;
1435
1436   // Check gamma
1437   Int_t pdg = mother->GetPdgCode();
1438   if(TMath::Abs(pdg) == 111) return imother;
1439   if(TMath::Abs(pdg) == 11) {
1440     return IsMotherPi0(imother);
1441   }
1442   return -1;
1443
1444 }
1445 //____________________________________________________________________________________________________________
1446 Int_t AliHFEelecbackground::IsMotherEta(Int_t tr) {
1447
1448   //
1449   // Return the lab of pi0 mother or -1 if not pi0
1450   //
1451
1452   AliStack* stack = fMCEvent->Stack();
1453   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1454
1455   // Take mother
1456   TParticle * particle = stack->Particle(tr);
1457   if(!particle) return -1;
1458   Int_t imother   = particle->GetFirstMother();
1459   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
1460   TParticle * mother = stack->Particle(imother);
1461   if(!mother) return -1;
1462
1463   // Check gamma
1464   Int_t pdg = mother->GetPdgCode();
1465   if(TMath::Abs(pdg) == 221) return imother;
1466   if(TMath::Abs(pdg) == 11) {
1467     return IsMotherEta(imother);
1468   }
1469   return -1;
1470
1471 }
1472 //____________________________________________________________________________________________________________
1473 Int_t AliHFEelecbackground::IsMotherC(Int_t tr) {
1474
1475   //
1476   // Return the lab of signal mother or -1 if not signal
1477   //
1478
1479   AliStack* stack = fMCEvent->Stack();
1480   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1481
1482   // Take mother
1483   TParticle * particle = stack->Particle(tr);
1484   if(!particle) return -1;
1485   Int_t imother   = particle->GetFirstMother();
1486   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
1487   TParticle * mother = stack->Particle(imother);
1488   if(!mother) return -1;
1489
1490   // Check gamma
1491   Int_t pdg = mother->GetPdgCode();
1492   if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
1493   if(TMath::Abs(pdg) == 11) {
1494     return IsMotherC(imother);
1495   }
1496   return -1;
1497
1498 }
1499 //____________________________________________________________________________________________________________
1500 Int_t AliHFEelecbackground::IsMotherB(Int_t tr) {
1501
1502   //
1503   // Return the lab of signal mother or -1 if not signal
1504   //
1505
1506   AliStack* stack = fMCEvent->Stack();
1507   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1508
1509   // Take mother
1510   TParticle * particle = stack->Particle(tr);
1511   if(!particle) return -1;
1512   Int_t imother   = particle->GetFirstMother();
1513   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
1514   TParticle * mother = stack->Particle(imother);
1515   if(!mother) return -1;
1516
1517   // Check gamma
1518   Int_t pdg = mother->GetPdgCode();
1519   if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
1520   if(TMath::Abs(pdg) == 11) {
1521     return IsMotherB(imother);
1522   }
1523   return -1;
1524
1525 }
1526 //____________________________________________________________________________________________________________
1527 Int_t AliHFEelecbackground::GetPdg(Int_t tr) {
1528
1529   //
1530   // Simply pdg code
1531   //
1532
1533   AliStack* stack = fMCEvent->Stack();
1534   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1535
1536   // MC Information
1537   TParticle * particle = stack->Particle(tr);
1538   if(!particle) return -1;
1539   Int_t pdg = particle->GetPdgCode();
1540
1541   return pdg;
1542
1543 }
1544 //____________________________________________________________________________________________________________
1545 Int_t AliHFEelecbackground::GetLabMother(Int_t tr) {
1546
1547   //
1548   // Simply lab mother
1549   //
1550
1551   AliStack* stack = fMCEvent->Stack();
1552   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
1553
1554   // MC Information
1555   TParticle * particle = stack->Particle(tr);
1556   if(!particle) return -1;
1557   Int_t imother = particle->GetFirstMother();
1558
1559   return imother;
1560
1561 }
1562 //_______________________________________________________________________________________________
1563 void AliHFEelecbackground::PostProcess()
1564 {
1565   //
1566   // Post process the histos and extract the background pt spectra
1567   //
1568
1569   if(!fList) return;
1570
1571   gStyle->SetPalette(1);
1572   gStyle->SetOptStat(1111);
1573   gStyle->SetPadBorderMode(0);
1574   gStyle->SetCanvasColor(10);
1575   gStyle->SetPadLeftMargin(0.13);
1576   gStyle->SetPadRightMargin(0.13);
1577
1578   /////////////////////////
1579   // Take the THnSparseF
1580   /////////////////////////
1581   THnSparseF *hsSparseData = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassData"));
1582   THnSparseF *hsSparseMC = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassMC"));
1583   THnSparseF *hsSparseCutPassedMC = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC"));
1584
1585   /////////////////////////////////
1586   // Cuts on the opening angle
1587   ////////////////////////////////
1588   if(!hsSparseData) return;
1589   TAxis *axisOpeningAngleData = hsSparseData->GetAxis(2);
1590   Int_t binCutData = axisOpeningAngleData->FindBin(fOpeningAngleCut);
1591   hsSparseData->GetAxis(2)->SetRange(1,binCutData);
1592
1593   if(hsSparseMC) {
1594     TAxis *axisOpeningAngleMC = hsSparseMC->GetAxis(2);
1595     Int_t binCutMC = axisOpeningAngleMC->FindBin(fOpeningAngleCut);
1596     hsSparseMC->GetAxis(2)->SetRange(1,binCutMC);
1597   }
1598
1599   /////////////////////////
1600   // Prepare the histos
1601   ////////////////////////
1602
1603   TAxis *ptaxisinvmass = hsSparseData->GetAxis(3);
1604   Int_t  nbinsptinvmass = ptaxisinvmass->GetNbins();
1605
1606   TH1D **invmassosptproj = new TH1D*[nbinsptinvmass];
1607   TH1D **invmassssptproj = new TH1D*[nbinsptinvmass];
1608   TH1D **invmassrptproj = new TH1D*[nbinsptinvmass];
1609   TH1D **invmassdiffptproj = new TH1D*[nbinsptinvmass];
1610   TH1D **invmassgammaptproj = new TH1D*[nbinsptinvmass];
1611   TH1D **invmasspi0ptproj = new TH1D*[nbinsptinvmass];
1612   TH1D **invmassetaptproj = new TH1D*[nbinsptinvmass];
1613   TH1D **invmassCptproj = new TH1D*[nbinsptinvmass];
1614   TH1D **invmassBptproj = new TH1D*[nbinsptinvmass];
1615
1616   TH1D *yieldPtFound = (TH1D *) hsSparseData->Projection(0);
1617   yieldPtFound->SetName("Found yield");
1618   yieldPtFound->Reset();
1619
1620   TH1D *yieldPtSourcesMC = 0x0;
1621   TH1D *yieldPtSignalCutMC = 0x0;
1622   if(hsSparseMC) {
1623     yieldPtSourcesMC = (TH1D *) hsSparseMC->Projection(0);
1624     yieldPtSourcesMC->SetName("Found yield");
1625     yieldPtSourcesMC->Reset();
1626
1627     yieldPtSignalCutMC = (TH1D *) hsSparseMC->Projection(0);
1628     yieldPtSignalCutMC->SetName("Found yield");
1629     yieldPtSignalCutMC->Reset();
1630   }
1631
1632   ////////////
1633   // canvas
1634   ///////////
1635   Int_t nbrow = (Int_t) (nbinsptinvmass/5);
1636   TString namecanvas("InvMassSpectra");
1637   TCanvas * canvas =new TCanvas(namecanvas,namecanvas,800,800);
1638   canvas->Divide(5,nbrow+1);
1639
1640   /////////////////////////////
1641   // Loop on pt bins
1642   /////////////////////////////
1643
1644   for(Int_t k=1; k <= nbinsptinvmass; k++){
1645
1646     Double_t lowedge = ptaxisinvmass->GetBinLowEdge(k);
1647     Double_t upedge  = ptaxisinvmass->GetBinUpEdge(k);
1648
1649     // Pt bin
1650     hsSparseData->GetAxis(0)->SetRange(k,k);
1651     if(hsSparseMC) hsSparseMC->GetAxis(0)->SetRange(k,k);
1652
1653     //
1654     hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);
1655     invmassosptproj[k-1] = hsSparseData->Projection(3);
1656     hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);
1657     invmassssptproj[k-1] = hsSparseData->Projection(3);
1658     hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);
1659     invmassrptproj[k-1]  = hsSparseData->Projection(3);
1660     hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins());
1661     invmassgammaptproj[k-1] = 0x0;
1662     invmasspi0ptproj[k-1] = 0x0;
1663     invmassetaptproj[k-1] = 0x0;
1664     invmassCptproj[k-1] = 0x0;
1665     invmassBptproj[k-1] = 0x0;
1666     if(hsSparseMC) {
1667       hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);
1668       invmassgammaptproj[k-1] = hsSparseMC->Projection(3);
1669       hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);
1670       invmasspi0ptproj[k-1] = hsSparseMC->Projection(3);
1671       hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);
1672       invmassetaptproj[k-1] = hsSparseMC->Projection(3);
1673       hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);
1674       invmassCptproj[k-1] = hsSparseMC->Projection(3);
1675       hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);
1676       invmassBptproj[k-1] = hsSparseMC->Projection(3);
1677       hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins());
1678     }
1679
1680     invmassdiffptproj[k-1] = (TH1D *) invmassosptproj[k-1]->Clone();
1681     TString name("Invmassdiffptbin");
1682     name += k;
1683     invmassdiffptproj[k-1]->SetName(name);
1684     invmassdiffptproj[k-1]->Add(invmassssptproj[k-1],-1.0);
1685
1686     TString namee("p_{T} tagged from ");
1687     namee += lowedge;
1688     namee += " GeV/c to ";
1689     namee += upedge;
1690     namee += " GeV/c";
1691
1692     invmassosptproj[k-1]->SetTitle((const char*)namee);
1693     invmassssptproj[k-1]->SetTitle((const char*)namee);
1694     invmassrptproj[k-1]->SetTitle((const char*)namee);
1695     invmassdiffptproj[k-1]->SetTitle((const char*)namee);
1696     if(invmassgammaptproj[k-1]) invmassgammaptproj[k-1]->SetTitle((const char*)namee);
1697     if(invmasspi0ptproj[k-1]) invmasspi0ptproj[k-1]->SetTitle((const char*)namee);
1698     if(invmassetaptproj[k-1]) invmassetaptproj[k-1]->SetTitle((const char*)namee);
1699     if(invmassCptproj[k-1]) invmassCptproj[k-1]->SetTitle((const char*)namee);
1700     if(invmassBptproj[k-1]) invmassBptproj[k-1]->SetTitle((const char*)namee);
1701
1702     invmassosptproj[k-1]->SetStats(0);
1703     invmassssptproj[k-1]->SetStats(0);
1704     invmassrptproj[k-1]->SetStats(0);
1705     invmassdiffptproj[k-1]->SetStats(0);
1706     if(invmassgammaptproj[k-1]) invmassgammaptproj[k-1]->SetStats(0);
1707     if(invmasspi0ptproj[k-1]) invmasspi0ptproj[k-1]->SetStats(0);
1708     if(invmassetaptproj[k-1]) invmassetaptproj[k-1]->SetStats(0);
1709     if(invmassCptproj[k-1]) invmassCptproj[k-1]->SetStats(0);
1710     if(invmassBptproj[k-1]) invmassBptproj[k-1]->SetStats(0);
1711
1712     Double_t yieldf = invmassdiffptproj[k-1]->Integral();
1713     if(invmassetaptproj[k-1] && invmasspi0ptproj[k-1] && invmassgammaptproj[k-1] && invmassCptproj[k-1] && invmassBptproj[k-1]) {
1714       Double_t yieldg = invmassetaptproj[k-1]->Integral() + invmasspi0ptproj[k-1]->Integral() + invmassgammaptproj[k-1]->Integral();
1715       if(yieldPtSourcesMC) yieldPtSourcesMC->SetBinContent(k,yieldg);
1716
1717       Double_t yieldsignal = invmassCptproj[k-1]->Integral() + invmassBptproj[k-1]->Integral();
1718       if(yieldPtSignalCutMC) yieldPtSignalCutMC->SetBinContent(k,yieldsignal);
1719     }
1720
1721     yieldPtFound->SetBinContent(k,yieldf);
1722
1723     canvas->cd(k);
1724     invmassosptproj[k-1]->Draw();
1725     invmassssptproj[k-1]->Draw("same");
1726     invmassdiffptproj[k-1]->Draw("same");
1727     invmassrptproj[k-1]->Draw("same");
1728     TLegend *legiv = new TLegend(0.4,0.6,0.89,0.89);
1729     legiv->AddEntry(invmassosptproj[k-1],"Opposite signs","p");
1730     legiv->AddEntry(invmassssptproj[k-1],"Same signs","p");
1731     legiv->AddEntry(invmassdiffptproj[k-1],"(Opposite - Same) signs","p");
1732     legiv->AddEntry(invmassrptproj[k-1],"rotated","p");
1733     if(invmassgammaptproj[k-1]) legiv->AddEntry(invmassgammaptproj[k-1],"e^{+}e^{-} from #gamma","p");
1734     if(invmasspi0ptproj[k-1]) legiv->AddEntry(invmasspi0ptproj[k-1],"e^{+}e^{-} from #pi^{0}","p");
1735     if(invmassetaptproj[k-1]) legiv->AddEntry(invmassetaptproj[k-1],"e^{+}e^{-} from #eta","p");
1736     legiv->Draw("same");
1737
1738     hsSparseData->GetAxis(0)->SetRange(1,hsSparseData->GetAxis(0)->GetNbins());
1739     if(hsSparseMC) hsSparseMC->GetAxis(0)->SetRange(1,hsSparseMC->GetAxis(0)->GetNbins());
1740
1741   }
1742
1743   ////////////////////////////////////////////////////
1744   // End of plotting: do subtraction of background
1745   ///////////////////////////////////////////////////
1746
1747   yieldPtFound->SetStats(0);
1748   if(yieldPtSourcesMC) yieldPtSourcesMC->SetStats(0);
1749   if(yieldPtSignalCutMC) yieldPtSignalCutMC->SetStats(0);
1750
1751   TCanvas * canvasfin =new TCanvas("ResultsElecBackGround","ResultsElecBackGround",800,800);
1752   canvasfin->cd(1);
1753   yieldPtFound->Draw();
1754   if(yieldPtSourcesMC && yieldPtSignalCutMC) {
1755     yieldPtSourcesMC->Draw("same");
1756     yieldPtSignalCutMC->Draw("same");
1757     TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1758     lega->AddEntry(yieldPtFound,"Contributions found","l");
1759     lega->AddEntry(yieldPtSourcesMC,"Contributions of e^{+}e^{-} from #gamma, #pi^{0} and #eta","l");
1760     lega->AddEntry(yieldPtSignalCutMC,"Contributions of e^{+}e^{-} from C and B","l");
1761     lega->Draw("same");
1762   }
1763
1764   if(hsSparseCutPassedMC){
1765     hsSparseCutPassedMC->GetAxis(1)->SetRange(1,1);
1766     hsSparseCutPassedMC->GetAxis(2)->SetRange(1,4);
1767     TH1D *hsSparseCutPassedMCproj = hsSparseCutPassedMC->Projection(0);
1768
1769     TH1D *cYieldPtFound = (TH1D*)yieldPtFound->Clone("RatioEfficiency");
1770     if(hsSparseCutPassedMCproj->Integral() > 0.0) cYieldPtFound->Divide(hsSparseCutPassedMCproj);
1771
1772     TCanvas * canvasfratio =new TCanvas("RatioEfficiency","RatioEfficiency",800,800);
1773     canvasfratio->cd(1);
1774     cYieldPtFound->Draw();
1775   }
1776
1777   //////////////////////////
1778   // fListPostProcess
1779   /////////////////////////
1780
1781   if(!fListPostProcess) fListPostProcess = new TList();
1782   fListPostProcess->SetName("ListPostProcess");
1783
1784   for(Int_t k=0; k < nbinsptinvmass; k++){
1785     fListPostProcess->AddAt(invmassosptproj[k],kOos+kNOutput*k);
1786     fListPostProcess->AddAt(invmassssptproj[k],kOss+kNOutput*k);
1787     fListPostProcess->AddAt(invmassrptproj[k],kOr+kNOutput*k);
1788     fListPostProcess->AddAt(invmassdiffptproj[k],kOdiff+kNOutput*k);
1789     if(invmassgammaptproj[k]) fListPostProcess->AddAt(invmassgammaptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromGamma);
1790     if(invmasspi0ptproj[k]) fListPostProcess->AddAt(invmasspi0ptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromPi0);
1791     if(invmassetaptproj[k]) fListPostProcess->AddAt(invmassetaptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromEta);
1792     if(invmassCptproj[k]) fListPostProcess->AddAt(invmassCptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromC);
1793     if(invmassBptproj[k]) fListPostProcess->AddAt(invmassBptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromB);
1794   }
1795
1796   fListPostProcess->AddAt(yieldPtFound,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass);
1797   if(yieldPtSourcesMC) fListPostProcess->AddAt(yieldPtSourcesMC,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass+1);
1798   if(yieldPtSignalCutMC) fListPostProcess->AddAt(yieldPtSignalCutMC,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass+2);
1799
1800   // delete dynamic array
1801   delete[] invmassosptproj;
1802   delete[] invmassssptproj;
1803   delete[] invmassrptproj;
1804   delete[] invmassdiffptproj;
1805   delete[] invmassgammaptproj;
1806   delete[] invmasspi0ptproj;
1807   delete[] invmassetaptproj;
1808   delete[] invmassCptproj;
1809   delete[] invmassBptproj;
1810
1811 }
1812 //_______________________________________________________________________________________________
1813 void AliHFEelecbackground::Plot() const
1814 {
1815   //
1816   // Plot the output
1817   //
1818
1819   if(!fList) return;
1820
1821   gStyle->SetPalette(1);
1822   gStyle->SetOptStat(1111);
1823   gStyle->SetPadBorderMode(0);
1824   gStyle->SetCanvasColor(10);
1825   gStyle->SetPadLeftMargin(0.13);
1826   gStyle->SetPadRightMargin(0.13);
1827
1828
1829   /////////////////////////
1830   // Take the THnSparseF
1831   /////////////////////////
1832   THnSparseF *hsSparseData = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassData"));
1833   THnSparseF *hsSparseMC = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassMC"));
1834   if(!hsSparseData) return;
1835
1836   ////////////////////
1837   // Opening angle
1838   ////////////////////
1839
1840   // Opening angle one direction
1841   hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1);
1842   TH1D *openingangleppproj = hsSparseData->Projection(2);
1843   hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1);
1844   TH1D *openinganglennproj = hsSparseData->Projection(2);
1845   hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);
1846   TH1D *openinganglessproj = hsSparseData->Projection(2);
1847   hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);
1848   TH1D *openinganglerproj  = hsSparseData->Projection(2);
1849   hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);
1850   TH1D *openingangleosproj = hsSparseData->Projection(2);
1851   hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins());
1852
1853   TH1D *openinganglegammaproj = 0x0;
1854   TH1D *openinganglepi0proj = 0x0;
1855   TH1D *openingangleCproj = 0x0;
1856   TH1D *openingangleBproj = 0x0;
1857   TH1D *openingangleetaproj = 0x0;
1858   TH1D *openingangleSplittedTrackssproj = 0x0;
1859   TH1D *openingangleSplittedTrackosproj = 0x0;
1860   if(hsSparseMC) {
1861     hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);
1862     openinganglegammaproj = hsSparseMC->Projection(2);
1863     hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);
1864     openinganglepi0proj = hsSparseMC->Projection(2);
1865     hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);
1866     openingangleetaproj = hsSparseMC->Projection(2);
1867     hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);
1868     openingangleCproj = hsSparseMC->Projection(2);
1869     hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);
1870     openingangleBproj = hsSparseMC->Projection(2);
1871     hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins());
1872     hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1);
1873     openingangleSplittedTrackssproj = hsSparseMC->Projection(2);
1874     hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1);
1875     openingangleSplittedTrackosproj = hsSparseMC->Projection(2);
1876     hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins());
1877   }
1878
1879   // Projection pt-opening angle
1880   hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1);
1881   TH2D *openingangleppproj2D = hsSparseData->Projection(0,2);
1882   hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1);
1883   TH2D *openinganglennproj2D = hsSparseData->Projection(0,2);
1884   hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);
1885   TH2D *openinganglessproj2D = hsSparseData->Projection(0,2);
1886   hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);
1887   TH2D *openinganglerproj2D  = hsSparseData->Projection(0,2);
1888   hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);
1889   TH2D *openingangleosproj2D = hsSparseData->Projection(0,2);
1890   hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins());
1891
1892   TH2D *openinganglegammaproj2D = 0x0;
1893   TH2D *openinganglepi0proj2D = 0x0;
1894   TH2D *openingangleCproj2D = 0x0;
1895   TH2D *openingangleBproj2D = 0x0;
1896   TH2D *openingangleetaproj2D = 0x0;
1897   TH2D *openingangleSplittedTrackssproj2D = 0x0;
1898   TH2D *openingangleSplittedTrackosproj2D = 0x0;
1899   if(hsSparseMC) {
1900     hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);
1901     openinganglegammaproj2D = hsSparseMC->Projection(0,2);
1902     hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);
1903     openinganglepi0proj2D = hsSparseMC->Projection(0,2);
1904     hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);
1905     openingangleetaproj2D = hsSparseMC->Projection(0,2);
1906     hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);
1907     openingangleCproj2D = hsSparseMC->Projection(0,2);
1908     hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);
1909     openingangleBproj2D = hsSparseMC->Projection(0,2);
1910     hsSparseMC->GetAxis(4)->SetRange(1, hsSparseMC->GetAxis(4)->GetNbins());
1911     hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1);
1912     openingangleSplittedTrackssproj2D = hsSparseMC->Projection(0,2);
1913     hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1);
1914     openingangleSplittedTrackosproj2D = hsSparseMC->Projection(0,2);
1915     hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins());
1916   }
1917
1918   openingangleppproj2D->SetStats(0);
1919   openinganglennproj2D->SetStats(0);
1920   openinganglessproj2D->SetStats(0);
1921   openinganglerproj2D->SetStats(0);
1922   openingangleosproj2D->SetStats(0);
1923   if(openinganglegammaproj2D) openinganglegammaproj2D->SetStats(0);
1924   if(openinganglepi0proj2D) openinganglepi0proj2D->SetStats(0);
1925   if(openingangleCproj2D) openingangleCproj2D->SetStats(0);
1926   if(openingangleBproj2D) openingangleBproj2D->SetStats(0);
1927   if(openingangleetaproj2D) openingangleetaproj2D->SetStats(0);
1928   if(openingangleSplittedTrackssproj2D) openingangleSplittedTrackssproj2D->SetStats(0);
1929   if(openingangleSplittedTrackosproj2D) openingangleSplittedTrackosproj2D->SetStats(0);
1930
1931   openingangleppproj2D->SetTitle("openingangleppproj2D");
1932   openinganglennproj2D->SetTitle("openinganglennproj2D");
1933   openinganglessproj2D->SetTitle("openinganglessproj2D");
1934   openinganglerproj2D->SetTitle("openinganglerproj2D");
1935   openingangleosproj2D->SetTitle("openingangleosproj2D");
1936   if(openinganglegammaproj2D) openinganglegammaproj2D->SetTitle("openinganglegammaproj2D");
1937   if(openinganglepi0proj2D) openinganglepi0proj2D->SetTitle("openinganglepi0proj2D");
1938   if(openingangleCproj2D) openingangleCproj2D->SetTitle("openingangleCproj2D");
1939   if(openingangleBproj2D) openingangleBproj2D->SetTitle("openingangleBproj2D");
1940   if(openingangleetaproj2D) openingangleetaproj2D->SetTitle("openingangleetaproj2D");
1941   if(openingangleSplittedTrackssproj2D) openingangleSplittedTrackssproj2D->SetTitle("openingangleSplittedTrackssproj2D");
1942   if(openingangleSplittedTrackosproj2D) openingangleSplittedTrackosproj2D->SetTitle("openingangleSplittedTrackosproj2D");
1943
1944   openingangleppproj->SetStats(0);
1945   openinganglennproj->SetStats(0);
1946   openinganglessproj->SetStats(0);
1947   openinganglerproj->SetStats(0);
1948   openingangleosproj->SetStats(0);
1949   if(openinganglegammaproj) openinganglegammaproj->SetStats(0);
1950   if(openinganglepi0proj) openinganglepi0proj->SetStats(0);
1951   if(openingangleCproj) openingangleCproj->SetStats(0);
1952   if(openingangleBproj) openingangleBproj->SetStats(0);
1953   if(openingangleetaproj) openingangleetaproj->SetStats(0);
1954   if(openingangleSplittedTrackssproj) openingangleSplittedTrackssproj->SetStats(0);
1955   if(openingangleSplittedTrackosproj) openingangleSplittedTrackosproj->SetStats(0);
1956
1957   openingangleppproj->SetTitle("");
1958   openinganglennproj->SetTitle("");
1959   openinganglessproj->SetTitle("");
1960   openinganglerproj->SetTitle("");
1961   openingangleosproj->SetTitle("");
1962   if(openinganglegammaproj) openinganglegammaproj->SetTitle("");
1963   if(openinganglepi0proj) openinganglepi0proj->SetTitle("");
1964   if(openingangleCproj) openingangleCproj->SetTitle("");
1965   if(openingangleBproj) openingangleBproj->SetTitle("");
1966   if(openingangleetaproj) openingangleetaproj->SetTitle("");
1967   if(openingangleSplittedTrackssproj) openingangleSplittedTrackssproj->SetTitle("");
1968   if(openingangleSplittedTrackosproj) openingangleSplittedTrackosproj->SetTitle("");
1969
1970   ////////////////////////////
1971   // Invariant mass
1972   ///////////////////////////
1973
1974   // Cuts on the opening angle
1975   TAxis *axisOpeningAngleData = hsSparseData->GetAxis(2);
1976   Int_t binCutData = axisOpeningAngleData->FindBin(fOpeningAngleCut);
1977   hsSparseData->GetAxis(2)->SetRange(1,binCutData);
1978
1979   // Debug
1980   //printf("Get Bin low edge %f, Get Bin Up edge %f for hsSparseData\n",axisOpeningAngleData->GetBinLowEdge(binCutData),axisOpeningAngleData->GetBinUpEdge(binCutData));
1981
1982   // Invariant mass
1983   hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1);
1984   TH1D *invmassppproj = hsSparseData->Projection(3);
1985   hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1);
1986   TH1D *invmassnnproj = hsSparseData->Projection(3);
1987   hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);
1988   TH1D *invmassssproj = hsSparseData->Projection(3);
1989   hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);
1990   TH1D *invmassrproj  = hsSparseData->Projection(3);
1991   hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);
1992   TH1D *invmassosproj = hsSparseData->Projection(3);
1993   hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins());
1994
1995   TH1D *invmassgammaproj = 0x0;
1996   TH1D *invmasspi0proj = 0x0;
1997   TH1D *invmassCproj = 0x0;
1998   TH1D *invmassBproj = 0x0;
1999   TH1D *invmassetaproj = 0x0;
2000   TH1D *invmassSplittedTrackssproj = 0x0;
2001   TH1D *invmassSplittedTrackosproj = 0x0;
2002   if(hsSparseMC) {
2003     TAxis *axisOpeningAngleMC = hsSparseMC->GetAxis(2);
2004     Int_t binCutMC = axisOpeningAngleMC->FindBin(fOpeningAngleCut);
2005     hsSparseMC->GetAxis(2)->SetRange(1,binCutMC);
2006
2007     // Debug
2008     //printf("Get Bin low edge %f, Get Bin Up edge %f for hsSparseMC\n",axisOpeningAngleMC->GetBinLowEdge(binCutMC),axisOpeningAngleMC->GetBinUpEdge(binCutMC));
2009
2010     hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);
2011     invmassgammaproj = hsSparseMC->Projection(3);
2012     hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);
2013     invmasspi0proj = hsSparseMC->Projection(3);
2014     hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);
2015     invmassetaproj = hsSparseMC->Projection(3);
2016     hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);
2017     invmassCproj = hsSparseMC->Projection(3);
2018     hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);
2019     invmassBproj = hsSparseMC->Projection(3);
2020     hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins());
2021     hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1);
2022     invmassSplittedTrackssproj = hsSparseMC->Projection(3);
2023     hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1);
2024     invmassSplittedTrackosproj = hsSparseMC->Projection(3);
2025     hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins());
2026   }
2027
2028   invmassppproj->SetStats(0);
2029   invmassnnproj->SetStats(0);
2030   invmassssproj->SetStats(0);
2031   invmassrproj->SetStats(0);
2032   invmassosproj->SetStats(0);
2033   if(invmassgammaproj) invmassgammaproj->SetStats(0);
2034   if(invmasspi0proj) invmasspi0proj->SetStats(0);
2035   if(invmassCproj) invmassCproj->SetStats(0);
2036   if(invmassBproj) invmassBproj->SetStats(0);
2037   if(invmassetaproj) invmassetaproj->SetStats(0);
2038   if(invmassSplittedTrackssproj) invmassSplittedTrackssproj->SetStats(0);
2039   if(invmassSplittedTrackosproj) invmassSplittedTrackosproj->SetStats(0);
2040
2041   invmassppproj->SetTitle("");
2042   invmassnnproj->SetTitle("");
2043   invmassssproj->SetTitle("");
2044   invmassrproj->SetTitle("");
2045   invmassosproj->SetTitle("");
2046   if(invmassgammaproj) invmassgammaproj->SetTitle("");
2047   if(invmasspi0proj) invmasspi0proj->SetTitle("");
2048   if(invmassCproj) invmassCproj->SetTitle("");
2049   if(invmassBproj) invmassBproj->SetTitle("");
2050   if(invmassetaproj) invmassetaproj->SetTitle("");
2051   if(invmassSplittedTrackssproj) invmassSplittedTrackssproj->SetTitle("");
2052   if(invmassSplittedTrackosproj) invmassSplittedTrackosproj->SetTitle("");
2053
2054   // Projection pt-invariant mass angle
2055   hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1);
2056   TH2D *invmassppproj2D = hsSparseData->Projection(0,3);
2057   hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1);
2058   TH2D *invmassnnproj2D = hsSparseData->Projection(0,3);
2059   hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);
2060   TH2D *invmassssproj2D = hsSparseData->Projection(0,3);
2061   hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);
2062   TH2D *invmassrproj2D  = hsSparseData->Projection(0,3);
2063   hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);
2064   TH2D *invmassosproj2D = hsSparseData->Projection(0,3);
2065   hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins());
2066
2067   TH2D *invmassgammaproj2D = 0x0;
2068   TH2D *invmasspi0proj2D = 0x0;
2069   TH2D *invmassCproj2D = 0x0;
2070   TH2D *invmassBproj2D = 0x0;
2071   TH2D *invmassetaproj2D = 0x0;
2072   TH2D *invmassSplittedTrackssproj2D = 0x0;
2073   TH2D *invmassSplittedTrackosproj2D = 0x0;
2074   if(hsSparseMC) {
2075     hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);
2076     invmassgammaproj2D = hsSparseMC->Projection(0,3);
2077     hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);
2078     invmasspi0proj2D = hsSparseMC->Projection(0,3);
2079     hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);
2080     invmassetaproj2D = hsSparseMC->Projection(0,3);
2081     hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);
2082     invmassCproj2D = hsSparseMC->Projection(0,3);
2083     hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);
2084     invmassBproj2D = hsSparseMC->Projection(0,3);
2085     hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins());
2086     hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1);
2087     invmassSplittedTrackssproj2D = hsSparseMC->Projection(0,3);
2088     hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1);
2089     invmassSplittedTrackosproj2D = hsSparseMC->Projection(0,3);
2090     hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins());
2091   }
2092
2093
2094   invmassppproj2D->SetStats(0);
2095   invmassnnproj2D->SetStats(0);
2096   invmassssproj2D->SetStats(0);
2097   invmassrproj2D->SetStats(0);
2098   invmassosproj2D->SetStats(0);
2099   if(invmassgammaproj2D) invmassgammaproj2D->SetStats(0);
2100   if(invmasspi0proj2D) invmasspi0proj2D->SetStats(0);
2101   if(invmassCproj2D) invmassCproj2D->SetStats(0);
2102   if(invmassBproj2D) invmassBproj2D->SetStats(0);
2103   if(invmassetaproj2D) invmassetaproj2D->SetStats(0);
2104   if(invmassSplittedTrackssproj2D) invmassSplittedTrackssproj2D->SetStats(0);
2105   if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->SetStats(0);
2106
2107   invmassppproj2D->SetTitle("invmassppproj2D");
2108   invmassnnproj2D->SetTitle("invmassnnproj2D");
2109   invmassssproj2D->SetTitle("invmassssproj2D");
2110   invmassrproj2D->SetTitle("invmassrproj2D");
2111   invmassosproj2D->SetTitle("invmassosproj2D");
2112   if(invmassgammaproj2D) invmassgammaproj2D->SetTitle("invmassgammaproj2D");
2113   if(invmasspi0proj2D) invmasspi0proj2D->SetTitle("invmasspi0proj2D");
2114   if(invmassCproj2D) invmassCproj2D->SetTitle("invmassCproj2D");
2115   if(invmassBproj2D) invmassBproj2D->SetTitle("invmassBproj2D");
2116   if(invmassetaproj2D) invmassetaproj2D->SetTitle("invmassetaproj2D");
2117   if(invmassSplittedTrackssproj2D) invmassSplittedTrackssproj2D->SetTitle("invmassSplittedTrackssproj2D");
2118   if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->SetTitle("invmassSplittedTrackosproj2D");
2119
2120
2121   /////////////
2122   // Plot
2123   ////////////
2124
2125   // Draw histograms for opening angle
2126   TCanvas * copeningangle =new TCanvas("openingangle","Openingangle",800,800);
2127   copeningangle->cd();
2128   //openingangleppproj->Draw();
2129   //openinganglennproj->Draw("same");
2130   openinganglessproj->Draw();
2131   //openinganglerproj->Draw("same");
2132   openingangleosproj->Draw("same");
2133   if(openinganglegammaproj) openinganglegammaproj->Draw("same");
2134   if(openinganglepi0proj) openinganglepi0proj->Draw("same");
2135   //if(openingangleCproj) openingangleCproj->Draw("same");
2136   //if(openingangleBproj) openingangleBproj->Draw("same");
2137   if(openingangleetaproj) openingangleetaproj->Draw("same");
2138   if(openingangleSplittedTrackssproj) openingangleSplittedTrackssproj->Draw("same");
2139   if(openingangleSplittedTrackosproj) openingangleSplittedTrackosproj->Draw("same");
2140   TLegend *lego = new TLegend(0.4,0.6,0.89,0.89);
2141   //lego->AddEntry(openingangleppproj,"positive-positive","p");
2142   //lego->AddEntry(openinganglennproj,"negative-negative","p");
2143   lego->AddEntry(openinganglessproj,"same-sign","p");
2144   //lego->AddEntry(openinganglerproj,"rotated","p");
2145   lego->AddEntry(openingangleosproj,"positive-negative","p");
2146   if(openinganglegammaproj) lego->AddEntry(openinganglegammaproj,"e^{+}e^{-} from #gamma","p");
2147   if(openinganglepi0proj) lego->AddEntry(openinganglepi0proj,"e^{+}e^{-} from #pi^{0}","p");
2148   //if(openingangleCproj) lego->AddEntry(openingangleCproj,"e^{+}e^{-} from c","p");
2149   //if(openingangleBproj) lego->AddEntry(openingangleBproj,"e^{+}e^{-} from b","p");
2150   if(openingangleetaproj) lego->AddEntry(openingangleetaproj,"e^{+}e^{-} from #eta","p");
2151   if(openingangleSplittedTrackssproj) lego->AddEntry(openingangleSplittedTrackssproj,"Splitted tracks same sign","p");
2152   if(openingangleSplittedTrackosproj) lego->AddEntry(openingangleSplittedTrackosproj,"Splitted tracks opposite sign","p");
2153   lego->Draw("same");
2154
2155   // Draw histograms for invariant mass
2156   TCanvas * cinvmass =new TCanvas("invmass","Invmass",800,800);
2157   cinvmass->cd();
2158   //invmassppproj->Draw();
2159   //invmassnnproj->Draw("same");
2160   invmassssproj->Draw();
2161   //invmassrproj->Draw("same");
2162   invmassosproj->Draw("same");
2163   if(invmassgammaproj) invmassgammaproj->Draw("same");
2164   if(invmasspi0proj) invmasspi0proj->Draw("same");
2165   //if(invmassCproj) invmassCproj->Draw("same");
2166   //if(invmassBproj) invmassBproj->Draw("same");
2167   if(invmassetaproj) invmassetaproj->Draw("same");
2168   if(invmassSplittedTrackssproj) invmassSplittedTrackssproj->Draw("same");
2169   if(invmassSplittedTrackosproj) invmassSplittedTrackosproj->Draw("same");
2170   TLegend *legi = new TLegend(0.4,0.6,0.89,0.89);
2171   //legi->AddEntry(invmassppproj,"positive-positive","p");
2172   //legi->AddEntry(invmassnnproj,"negative-negative","p");
2173   legi->AddEntry(invmassssproj,"same-sign","p");
2174   //legi->AddEntry(invmassrproj,"rotated","p");
2175   legi->AddEntry(invmassosproj,"positive-negative","p");
2176   if(invmassgammaproj) legi->AddEntry(invmassgammaproj,"e^{+}e^{-} from #gamma","p");
2177   if(invmasspi0proj) legi->AddEntry(invmasspi0proj,"e^{+}e^{-} from #pi^{0}","p");
2178   //if(invmassCproj) legi->AddEntry(invmassCproj,"e^{+}e^{-} from c","p");
2179   //if(invmassBproj) legi->AddEntry(invmassBproj,"e^{+}e^{-} from b","p");
2180   if(invmassetaproj) legi->AddEntry(invmassetaproj,"e^{+}e^{-} from #eta","p");
2181   if(invmassSplittedTrackssproj) legi->AddEntry(invmassSplittedTrackssproj,"Splitted tracks same sign","p");
2182   if(invmassSplittedTrackosproj) legi->AddEntry(invmassSplittedTrackosproj,"Splitted tracks opposite sign","p");
2183   legi->Draw("same");
2184
2185
2186
2187   // Draw histograms for opening angle 2D
2188   TCanvas * copeningangle2D =new TCanvas("openingangle2D","Openingangle2D",800,800);
2189   copeningangle2D->Divide(6,2);
2190   copeningangle2D->cd(1);
2191   openingangleppproj2D->Draw("lego");
2192   copeningangle2D->cd(2);
2193   openinganglennproj2D->Draw("lego");
2194   copeningangle2D->cd(3);
2195   openinganglessproj2D->Draw("lego");
2196   copeningangle2D->cd(4);
2197   openinganglerproj2D->Draw("lego");
2198   copeningangle2D->cd(5);
2199   openingangleosproj2D->Draw("lego");
2200   copeningangle2D->cd(6);
2201   if(openinganglegammaproj2D) openinganglegammaproj2D->Draw("lego");
2202   copeningangle2D->cd(7);
2203   if(openinganglepi0proj2D) openinganglepi0proj2D->Draw("lego");
2204   copeningangle2D->cd(8);
2205   if(openingangleCproj2D) openingangleCproj2D->Draw("lego");
2206   copeningangle2D->cd(9);
2207   if(openingangleBproj2D) openingangleBproj2D->Draw("lego");
2208   copeningangle2D->cd(10);
2209   if(openingangleetaproj2D) openingangleetaproj2D->Draw("lego");
2210   copeningangle2D->cd(11);
2211   if(openingangleSplittedTrackssproj2D) openingangleSplittedTrackssproj2D->Draw("lego");
2212   copeningangle2D->cd(12);
2213   if(openingangleSplittedTrackosproj2D) openingangleSplittedTrackosproj2D->Draw("lego");
2214
2215   // Draw histograms for invariant mass 2D
2216   TCanvas * cinvmass2D =new TCanvas("invmass2D","Invmass2D",800,800);
2217   cinvmass2D->Divide(6,2);
2218   cinvmass2D->cd(1);
2219   invmassppproj2D->Draw("lego");
2220   cinvmass2D->cd(2);
2221   invmassnnproj2D->Draw("lego");
2222   cinvmass2D->cd(3);
2223   invmassssproj2D->Draw("lego");
2224   cinvmass2D->cd(4);
2225   invmassrproj2D->Draw("lego");
2226   cinvmass2D->cd(5);
2227   invmassosproj2D->Draw("lego");
2228   cinvmass2D->cd(6);
2229   if(invmassgammaproj2D) invmassgammaproj2D->Draw("lego");
2230   cinvmass2D->cd(7);
2231   if(invmasspi0proj2D) invmasspi0proj2D->Draw("lego");
2232   cinvmass2D->cd(8);
2233   if(invmassCproj2D) invmassCproj2D->Draw("lego");
2234   cinvmass2D->cd(9);
2235   if(invmassBproj2D) invmassBproj2D->Draw("lego");
2236   cinvmass2D->cd(10);
2237   if(invmassetaproj2D) invmassetaproj2D->Draw("lego");
2238   cinvmass2D->cd(11);
2239   if(invmassSplittedTrackssproj2D) invmassSplittedTrackssproj2D->Draw("lego");
2240   cinvmass2D->cd(12);
2241   if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->Draw("lego");
2242
2243
2244   ////////////////////////
2245   // Cut efficiencies
2246   ////////////////////////
2247
2248   THnSparseF *hsSparseMCe = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC"));
2249
2250   if(hsSparseMCe) {
2251
2252     // init histos
2253     TAxis *axissources = hsSparseMCe->GetAxis(2);
2254     Int_t  nbsources = axissources->GetNbins();
2255     TAxis *axiscuts = hsSparseMCe->GetAxis(1);
2256     Int_t  nbcuts = axiscuts->GetNbins();
2257     TH1D **histopassedcuts = new TH1D*[nbsources*nbcuts];
2258     Double_t *nbEntriesCuts = new Double_t[nbsources*nbcuts];
2259
2260     //printf("Number of cuts %d\n",nbcuts);
2261
2262     // canvas
2263     TCanvas * chsSparseMCeeff =new TCanvas("hsSparseMCeeffDebug","hsSparseMCeeffDebug",800,800);
2264     chsSparseMCeeff->Divide(3,1);
2265
2266     // histos
2267     for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
2268       hsSparseMCe->GetAxis(2)->SetRange(sourceid+1,sourceid+1);
2269       for(Int_t cut = 0; cut < nbcuts; cut++){
2270         hsSparseMCe->GetAxis(1)->SetRange(cut+1,cut+1);
2271         histopassedcuts[sourceid*nbcuts+cut] = hsSparseMCe->Projection(0);
2272         hsSparseMCe->GetAxis(1)->SetRange(1,hsSparseMCe->GetAxis(1)->GetNbins());
2273       }
2274       hsSparseMCe->GetAxis(2)->SetRange(1,hsSparseMCe->GetAxis(2)->GetNbins());
2275     }
2276
2277     // calcul efficiencies
2278     ///////////////////////
2279     // histos
2280     for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
2281       // Next is compared to the partner tracked
2282       for(Int_t cut = 2; cut < nbcuts; cut++){
2283         nbEntriesCuts[sourceid*nbcuts+cut] = histopassedcuts[sourceid*nbcuts+cut]->GetEntries();
2284         if(histopassedcuts[sourceid*nbcuts+1]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+cut]->Divide(histopassedcuts[sourceid*nbcuts+1]);
2285       }
2286       // First one is if the partner is tracked.
2287       nbEntriesCuts[sourceid*nbcuts+1] = histopassedcuts[sourceid*nbcuts+1]->GetEntries();
2288       if(histopassedcuts[sourceid*nbcuts]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+1]->Divide(histopassedcuts[sourceid*nbcuts]);
2289       // First one is input
2290       nbEntriesCuts[sourceid*nbcuts] = histopassedcuts[sourceid*nbcuts]->GetEntries();
2291     }
2292
2293     /////////////
2294     // ratios
2295     ////////////
2296     for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
2297       for(Int_t cut = 1; cut < nbcuts; cut++){
2298         if(nbEntriesCuts[sourceid*nbcuts] > 0.0) nbEntriesCuts[sourceid*nbcuts+cut] = nbEntriesCuts[sourceid*nbcuts+cut]/nbEntriesCuts[sourceid*nbcuts];
2299       }
2300     }
2301     TH1F *ratioHistoEntriesGamma = new TH1F("ratioHistoEntriesGamma","", nbcuts-1, 0.0, nbcuts-1.0);
2302     TH1F *ratioHistoEntriesPi0 = new TH1F("ratioHistoEntriesPi0","", nbcuts-1, 0.0, nbcuts-1.0);
2303     TH1F *ratioHistoEntriesC = new TH1F("ratioHistoEntriesC","", nbcuts-1, 0.0, nbcuts-1.0);
2304     for(Int_t k = 1; k < nbcuts; k++){
2305       ratioHistoEntriesGamma->SetBinContent(k,nbEntriesCuts[nbcuts+k]);
2306       ratioHistoEntriesPi0->SetBinContent(k,nbEntriesCuts[2*nbcuts+k]);
2307       ratioHistoEntriesC->SetBinContent(k,nbEntriesCuts[4*nbcuts+k]);
2308     }
2309     //
2310     TAxis *xAxisGamma = ratioHistoEntriesGamma->GetXaxis();
2311     xAxisGamma->SetBinLabel(1,"Partner tracked");
2312     xAxisGamma->SetBinLabel(2,"Opposite sign");
2313     xAxisGamma->SetBinLabel(3,"Single Track Cut");
2314     xAxisGamma->SetBinLabel(4,"Shared Clusters");
2315     xAxisGamma->SetBinLabel(5,"PID");
2316     xAxisGamma->SetBinLabel(6,"DCA");
2317     xAxisGamma->SetBinLabel(7,"Chi^{2}/Ndf");
2318     xAxisGamma->SetBinLabel(8,"Opening angle");
2319     xAxisGamma->SetBinLabel(9,"Invariant mass");
2320     //
2321     TAxis *xAxisPi0 = ratioHistoEntriesPi0->GetXaxis();
2322     xAxisPi0->SetBinLabel(1,"Partner tracked");
2323     xAxisPi0->SetBinLabel(2,"Opposite sign");
2324     xAxisPi0->SetBinLabel(3,"Single Track Cut");
2325     xAxisPi0->SetBinLabel(4,"Shared Clusters");
2326     xAxisPi0->SetBinLabel(5,"PID");
2327     xAxisPi0->SetBinLabel(6,"DCA");
2328     xAxisPi0->SetBinLabel(7,"Chi^{2}/Ndf");
2329     xAxisPi0->SetBinLabel(8,"Opening angle");
2330     xAxisPi0->SetBinLabel(9,"Invariant mass");
2331     //
2332     TAxis *xAxisC = ratioHistoEntriesC->GetXaxis();
2333     xAxisC->SetBinLabel(1,"Partner tracked");
2334     xAxisC->SetBinLabel(2,"Opposite sign");
2335     xAxisC->SetBinLabel(3,"Single Track Cut");
2336     xAxisC->SetBinLabel(4,"Shared Clusters");
2337     xAxisC->SetBinLabel(5,"PID");
2338     xAxisC->SetBinLabel(6,"DCA");
2339     xAxisC->SetBinLabel(7,"Chi^{2}/Ndf");
2340     xAxisC->SetBinLabel(8,"Opening angle");
2341     xAxisC->SetBinLabel(9,"Invariant mass");
2342     //
2343     TCanvas * cRatioHistoEntries =new TCanvas("cRatioHistoEntries","cRatioHistoEntries",800,800);
2344     cRatioHistoEntries->cd(1);
2345     ratioHistoEntriesGamma->SetStats(0);
2346     ratioHistoEntriesGamma->Draw();
2347     ratioHistoEntriesPi0->SetStats(0);
2348     ratioHistoEntriesPi0->Draw("same");
2349     ratioHistoEntriesC->SetStats(0);
2350     //ratioHistoEntriesC->Draw("same");
2351     TLegend *legEntries = new TLegend(0.4,0.6,0.89,0.89);
2352     legEntries->AddEntry(ratioHistoEntriesGamma,"#gamma","l");
2353     legEntries->AddEntry(ratioHistoEntriesPi0,"#pi^{0}","l");
2354     //legEntries->AddEntry(ratioHistoEntriesC,"c","p");
2355     legEntries->Draw("same");
2356
2357     ////////////////////
2358     // plot Debug
2359     ///////////////////
2360     Int_t source = 1;
2361     chsSparseMCeeff->cd(1);
2362     histopassedcuts[source*nbcuts+0]->SetTitle("#gamma");
2363     histopassedcuts[source*nbcuts+1]->SetTitle("#gamma");
2364     histopassedcuts[source*nbcuts+2]->SetTitle("#gamma");
2365     histopassedcuts[source*nbcuts+3]->SetTitle("#gamma");
2366     histopassedcuts[source*nbcuts+4]->SetTitle("#gamma");
2367     histopassedcuts[source*nbcuts+5]->SetTitle("#gamma");
2368     histopassedcuts[source*nbcuts+6]->SetTitle("#gamma");
2369     histopassedcuts[source*nbcuts+7]->SetTitle("#gamma");
2370     histopassedcuts[source*nbcuts+8]->SetTitle("#gamma");
2371     histopassedcuts[source*nbcuts+9]->SetTitle("#gamma");
2372     //histopassedcuts[source*nbcuts+0]->SetStats(0);
2373     histopassedcuts[source*nbcuts+1]->SetStats(0);
2374     histopassedcuts[source*nbcuts+2]->SetStats(0);
2375     histopassedcuts[source*nbcuts+3]->SetStats(0);
2376     histopassedcuts[source*nbcuts+4]->SetStats(0);
2377     histopassedcuts[source*nbcuts+5]->SetStats(0);
2378     histopassedcuts[source*nbcuts+6]->SetStats(0);
2379     histopassedcuts[source*nbcuts+7]->SetStats(0);
2380     histopassedcuts[source*nbcuts+8]->SetStats(0);
2381     histopassedcuts[source*nbcuts+9]->SetStats(0);
2382     //histopassedcuts[source*nbcuts+0]->Draw();
2383     //histopassedcuts[source*nbcuts+1]->Draw("");
2384     histopassedcuts[source*nbcuts+2]->Draw();
2385     histopassedcuts[source*nbcuts+3]->Draw("same");
2386     //histopassedcuts[source*nbcuts+4]->Draw("same");
2387     histopassedcuts[source*nbcuts+5]->Draw("same");
2388     histopassedcuts[source*nbcuts+6]->Draw("same");
2389     //histopassedcuts[source*nbcuts+7]->Draw("same");
2390     histopassedcuts[source*nbcuts+8]->Draw("same");
2391     histopassedcuts[source*nbcuts+9]->Draw("same");
2392     TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
2393     //legb->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
2394     //legb->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
2395     legb->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
2396     legb->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
2397     //legb->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
2398     legb->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
2399     legb->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
2400     //legb->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
2401     legb->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
2402     legb->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
2403     legb->Draw("same");
2404
2405     source = 2;
2406     chsSparseMCeeff->cd(2);
2407     histopassedcuts[source*nbcuts+0]->SetTitle("#pi^{0}");
2408     histopassedcuts[source*nbcuts+1]->SetTitle("#pi^{0}");
2409     histopassedcuts[source*nbcuts+2]->SetTitle("#pi^{0}");
2410     histopassedcuts[source*nbcuts+3]->SetTitle("#pi^{0}");
2411     histopassedcuts[source*nbcuts+4]->SetTitle("#pi^{0}");
2412     histopassedcuts[source*nbcuts+5]->SetTitle("#pi^{0}");
2413     histopassedcuts[source*nbcuts+6]->SetTitle("#pi^{0}");
2414     histopassedcuts[source*nbcuts+7]->SetTitle("#pi^{0}");
2415     histopassedcuts[source*nbcuts+8]->SetTitle("#pi^{0}");
2416     histopassedcuts[source*nbcuts+9]->SetTitle("#pi^{0}");
2417     //histopassedcuts[source*nbcuts+0]->SetStats(0);
2418     histopassedcuts[source*nbcuts+1]->SetStats(0);
2419     histopassedcuts[source*nbcuts+2]->SetStats(0);
2420     histopassedcuts[source*nbcuts+3]->SetStats(0);
2421     histopassedcuts[source*nbcuts+4]->SetStats(0);
2422     histopassedcuts[source*nbcuts+5]->SetStats(0);
2423     histopassedcuts[source*nbcuts+6]->SetStats(0);
2424     histopassedcuts[source*nbcuts+7]->SetStats(0);
2425     histopassedcuts[source*nbcuts+8]->SetStats(0);
2426     histopassedcuts[source*nbcuts+9]->SetStats(0);
2427     //histopassedcuts[source*nbcuts+0]->Draw();
2428     //histopassedcuts[source*nbcuts+1]->Draw();
2429     histopassedcuts[source*nbcuts+2]->Draw();
2430     histopassedcuts[source*nbcuts+3]->Draw("same");
2431     //histopassedcuts[source*nbcuts+4]->Draw("same");
2432     histopassedcuts[source*nbcuts+5]->Draw("same");
2433     histopassedcuts[source*nbcuts+6]->Draw("same");
2434     //histopassedcuts[source*nbcuts+7]->Draw("same");
2435     histopassedcuts[source*nbcuts+8]->Draw("same");
2436     histopassedcuts[source*nbcuts+9]->Draw("same");
2437     TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
2438     //legc->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
2439     //legc->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
2440     legc->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
2441     legc->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
2442     //legc->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
2443     legc->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
2444     legc->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
2445     //legc->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
2446     legc->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
2447     legc->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
2448     legc->Draw("same");
2449
2450     source = 4;
2451     chsSparseMCeeff->cd(3);
2452     histopassedcuts[source*nbcuts+0]->SetTitle("C");
2453     histopassedcuts[source*nbcuts+1]->SetTitle("C");
2454     histopassedcuts[source*nbcuts+2]->SetTitle("C");
2455     histopassedcuts[source*nbcuts+3]->SetTitle("C");
2456     histopassedcuts[source*nbcuts+4]->SetTitle("C");
2457     histopassedcuts[source*nbcuts+5]->SetTitle("C");
2458     histopassedcuts[source*nbcuts+6]->SetTitle("C");
2459     histopassedcuts[source*nbcuts+7]->SetTitle("C");
2460     histopassedcuts[source*nbcuts+8]->SetTitle("C");
2461     histopassedcuts[source*nbcuts+9]->SetTitle("C");
2462     //histopassedcuts[source*nbcuts+0]->SetStats(0);
2463     histopassedcuts[source*nbcuts+1]->SetStats(0);
2464     histopassedcuts[source*nbcuts+2]->SetStats(0);
2465     histopassedcuts[source*nbcuts+3]->SetStats(0);
2466     histopassedcuts[source*nbcuts+4]->SetStats(0);
2467     histopassedcuts[source*nbcuts+5]->SetStats(0);
2468     histopassedcuts[source*nbcuts+6]->SetStats(0);
2469     histopassedcuts[source*nbcuts+7]->SetStats(0);
2470     histopassedcuts[source*nbcuts+8]->SetStats(0);
2471     histopassedcuts[source*nbcuts+9]->SetStats(0);
2472     //histopassedcuts[source*nbcuts+0]->Draw();
2473     //histopassedcuts[source*nbcuts+1]->Draw();
2474     histopassedcuts[source*nbcuts+2]->Draw();
2475     histopassedcuts[source*nbcuts+3]->Draw("same");
2476     //histopassedcuts[source*nbcuts+4]->Draw("same");
2477     histopassedcuts[source*nbcuts+5]->Draw("same");
2478     histopassedcuts[source*nbcuts+6]->Draw("same");
2479     //histopassedcuts[source*nbcuts+7]->Draw("same");
2480     histopassedcuts[source*nbcuts+8]->Draw("same");
2481     histopassedcuts[source*nbcuts+9]->Draw("same");
2482     TLegend *lege = new TLegend(0.4,0.6,0.89,0.89);
2483     //lege->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
2484     //lege->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
2485     lege->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
2486     lege->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
2487     //lege->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
2488     lege->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
2489     lege->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
2490     //lege->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
2491     lege->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
2492     lege->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
2493     lege->Draw("same");
2494
2495     //////////////////////
2496     // Input
2497     //////////////////////
2498
2499     TCanvas * chsSparseMCein =new TCanvas("hsSparseMCeinput","hsSparseMCeinput",800,800);
2500     chsSparseMCein->cd(1);
2501     Double_t nbGamma = 0.0;
2502     source = 1;
2503     nbGamma = histopassedcuts[source*nbcuts+0]->GetEntries();
2504     histopassedcuts[source*nbcuts+0]->SetStats(0);
2505     histopassedcuts[source*nbcuts+0]->Draw();
2506     TLegend *leginput = new TLegend(0.4,0.6,0.89,0.89);
2507     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#gamma","p");
2508     Double_t nbPi0 = 0.0;
2509     source = 2;
2510     nbPi0 = histopassedcuts[source*nbcuts+0]->GetEntries();
2511     histopassedcuts[source*nbcuts+0]->SetStats(0);
2512     histopassedcuts[source*nbcuts+0]->Draw("same");
2513     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#pi^{0}","p");
2514     Double_t nbEta = 0.0;
2515     source = 3;
2516     nbEta = histopassedcuts[source*nbcuts+0]->GetEntries();
2517     histopassedcuts[source*nbcuts+0]->SetStats(0);
2518     histopassedcuts[source*nbcuts+0]->Draw("same");
2519     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#eta","p");
2520     Double_t nbC = 0.0;
2521     source = 4;
2522     nbC = histopassedcuts[source*nbcuts+0]->GetEntries();
2523     histopassedcuts[source*nbcuts+0]->SetStats(0);
2524     histopassedcuts[source*nbcuts+0]->Draw("same");
2525     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"c","p");
2526     leginput->Draw("same");
2527
2528     //printf("Gamma %f, pi^{0} %f and #eta %f, c %f\n",nbGamma,nbPi0,nbEta,nbC);
2529
2530     //////////////////////
2531     // Tracked
2532     //////////////////////
2533
2534     TCanvas * cTracked = new TCanvas("cTracked","cTracked",800,800);
2535     cTracked->cd(1);
2536     source = 1;
2537     histopassedcuts[source*nbcuts+1]->Draw();
2538     TLegend *legTracked = new TLegend(0.4,0.6,0.89,0.89);
2539     legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#gamma","p");
2540     source = 2;
2541     histopassedcuts[source*nbcuts+1]->Draw("same");
2542     legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#pi^{0}","p");
2543     legTracked->Draw("same");
2544
2545   }
2546
2547   /////////////////////////////////////
2548   // Data Radius and chi2Ndf if AliKF
2549   ////////////////////////////////////
2550
2551   TH1F *hDataRadius = dynamic_cast<TH1F *>(fList->FindObject("DataRadius"));
2552   TH1F *hDataChi2Ndf = dynamic_cast<TH1F *>(fList->FindObject("DataChi2Ndf"));
2553
2554   if(hDataRadius || hDataChi2Ndf) {
2555     TCanvas * cDataRadiusChi2Ndf =new TCanvas("CanvasDataRadiusChi2Ndf","CanvasDataRadiusChi2Ndf",800,800);
2556     cDataRadiusChi2Ndf->Divide(2,1);
2557     cDataRadiusChi2Ndf->cd(1);
2558     if(hDataRadius) hDataRadius->Draw();
2559     cDataRadiusChi2Ndf->cd(2);
2560     if(hDataChi2Ndf) hDataChi2Ndf->Draw();
2561   }
2562
2563   ///////////////////////
2564   // Data DCA
2565   //////////////////////
2566
2567   TH1F *hDataDCA = dynamic_cast<TH1F *>(fList->FindObject("DataDCA"));
2568
2569   if(hDataDCA) {
2570     TCanvas * cDataDCA =new TCanvas("CanvasDataDCA","CanvasDataDCA",800,800);
2571     cDataDCA->cd(1);
2572     hDataDCA->Draw();
2573   }
2574
2575   /////////////////////////////////////
2576   // MC Radius and chi2Ndf if AliKF
2577   ////////////////////////////////////
2578
2579   TH2F *hMCRadius = dynamic_cast<TH2F *>(fList->FindObject("MCRadius"));
2580   TH2F *hMCChi2Ndf = dynamic_cast<TH2F *>(fList->FindObject("MCChi2Ndf"));
2581
2582   if(hMCRadius || hMCChi2Ndf) {
2583     TCanvas * cMCRadiusChi2Ndf =new TCanvas("CanvasMCRadiusChi2Ndf","CanvasMCRadiusChi2Ndf",800,800);
2584     cMCRadiusChi2Ndf->Divide(2,1);
2585     cMCRadiusChi2Ndf->cd(1);
2586     //TH1D *hMCRadiusBackground = hMCRadius->ProjectionX("MCRadiusBackGround",1,1,"e");
2587     TH1D *hMCRadiusGamma = hMCRadius->ProjectionX("MCRadiusGamma",2,2,"e");
2588     TH1D *hMCRadiusPi0 = hMCRadius->ProjectionX("MCRadiusPi0",3,3,"e");
2589     TH1D *hMCRadiusEta = hMCRadius->ProjectionX("MCRadiusEta",4,4,"e");
2590     TH1D *hMCRadiusC = hMCRadius->ProjectionX("MCRadiusC",5,5,"e");
2591     TH1D *hMCRadiusB = hMCRadius->ProjectionX("MCRadiusB",6,6,"e");
2592     //hMCRadiusBackground->Draw();
2593     hMCRadiusGamma->Draw();
2594     hMCRadiusPi0->Draw("same");
2595     hMCRadiusEta->Draw("same");
2596     hMCRadiusC->Draw("same");
2597     hMCRadiusB->Draw("same");
2598     TLegend *legRadius = new TLegend(0.4,0.6,0.89,0.89);
2599     //legRadius->AddEntry(hMCRadiusBackground,"Background","p");
2600     legRadius->AddEntry(hMCRadiusGamma,"#gamma","p");
2601     legRadius->AddEntry(hMCRadiusPi0,"#pi^{0}","p");
2602     legRadius->AddEntry(hMCRadiusEta,"#eta","p");
2603     legRadius->AddEntry(hMCRadiusC,"c","p");
2604     legRadius->AddEntry(hMCRadiusB,"b","p");
2605     legRadius->Draw("same");
2606     cMCRadiusChi2Ndf->cd(2);
2607     //TH1D *hMCChi2NdfBackground = hMCChi2Ndf->ProjectionX("MCChi2NdfBackGround",1,1,"e");
2608     TH1D *hMCChi2NdfGamma = hMCChi2Ndf->ProjectionX("MCChi2NdfGamma",2,2,"e");
2609     TH1D *hMCChi2NdfPi0 = hMCChi2Ndf->ProjectionX("MCChi2NdfPi0",3,3,"e");
2610     TH1D *hMCChi2NdfEta = hMCChi2Ndf->ProjectionX("MCChi2NdfEta",4,4,"e");
2611     TH1D *hMCChi2NdfC = hMCChi2Ndf->ProjectionX("MCChi2NdfC",5,5,"e");
2612     TH1D *hMCChi2NdfB = hMCChi2Ndf->ProjectionX("MCChi2NdfB",6,6,"e");
2613     //hMCChi2NdfBackground->Draw();
2614     hMCChi2NdfGamma->Draw();
2615     hMCChi2NdfPi0->Draw("same");
2616     hMCChi2NdfEta->Draw("same");
2617     hMCChi2NdfC->Draw("same");
2618     hMCChi2NdfB->Draw("same");
2619     TLegend *legChi2Ndf = new TLegend(0.4,0.6,0.89,0.89);
2620     //legChi2Ndf->AddEntry(hMCChi2NdfBackground,"Background","p");
2621     legChi2Ndf->AddEntry(hMCChi2NdfGamma,"#gamma","p");
2622     legChi2Ndf->AddEntry(hMCChi2NdfPi0,"#pi^{0}","p");
2623     legChi2Ndf->AddEntry(hMCChi2NdfEta,"#eta","p");
2624     legChi2Ndf->AddEntry(hMCChi2NdfC,"c","p");
2625     legChi2Ndf->AddEntry(hMCChi2NdfB,"b","p");
2626     legChi2Ndf->Draw("same");
2627   }
2628
2629   ///////////////////////
2630   // MC DCA
2631   //////////////////////
2632
2633   TH2F *hMCDCA = dynamic_cast<TH2F *>(fList->FindObject("MCDCA"));
2634
2635   if(hMCDCA) {
2636     TCanvas * cMCDCA =new TCanvas("CanvasMCDCA","CanvasMCDCA",800,800);
2637     cMCDCA->cd(1);
2638     //TH1D *hMCDCABackground = hMCDCA->ProjectionX("MCDCABackGround",1,1,"e");
2639     TH1D *hMCDCAGamma = hMCDCA->ProjectionX("MCDCAGamma",2,2,"e");
2640     TH1D *hMCDCAPi0 = hMCDCA->ProjectionX("MCDCAPi0",3,3,"e");
2641     TH1D *hMCDCAEta = hMCDCA->ProjectionX("MCDCAEta",4,4,"e");
2642     TH1D *hMCDCAC = hMCDCA->ProjectionX("MCDCAC",5,5,"e");
2643     TH1D *hMCDCAB = hMCDCA->ProjectionX("MCDCAB",6,6,"e");
2644     //hMCDCABackground->Draw();
2645     hMCDCAGamma->Draw();
2646     hMCDCAPi0->Draw("same");
2647     hMCDCAEta->Draw("same");
2648     hMCDCAC->Draw("same");
2649     hMCDCAB->Draw("same");
2650     TLegend *legDCA = new TLegend(0.4,0.6,0.89,0.89);
2651     //legDCA->AddEntry(hMCDCABackground,"Background","p");
2652     legDCA->AddEntry(hMCDCAGamma,"#gamma","p");
2653     legDCA->AddEntry(hMCDCAPi0,"#pi^{0}","p");
2654     legDCA->AddEntry(hMCDCAEta,"#eta","p");
2655     legDCA->AddEntry(hMCDCAC,"c","p");
2656     legDCA->AddEntry(hMCDCAB,"b","p");
2657     legDCA->Draw("same");
2658   }
2659
2660
2661   /////////////////////////
2662   // PID Partner Signal
2663   /////////////////////////
2664   TF1 *betheBlochElectron = 0x0;
2665   TF1 *betheBlochMuon = 0x0;
2666   TF1 *betheBlochPion = 0x0;
2667   TF1 *betheBlochKaon = 0x0;
2668   TF1 *betheBlochProton = 0x0;
2669
2670   TH2F *hsignalPidPartner0 = dynamic_cast<TH2F *>(fList->FindObject("TPCPartner0"));
2671   TH2F *hsignalPidPartner1 = dynamic_cast<TH2F *>(fList->FindObject("TPCPartner1"));
2672   if((!hsignalPidPartner0) && (!hsignalPidPartner1)) {
2673     hsignalPidPartner0 = dynamic_cast<TH2F *>(fList->FindObject("ITSPartner0"));
2674     hsignalPidPartner1 = dynamic_cast<TH2F *>(fList->FindObject("ITSPartner1"));
2675
2676     betheBlochElectron = new TF1("betheBlochElectron",BetheBlochElectronITS,0.1,10.0,0);
2677     betheBlochMuon = new TF1("betheBlochMuon",BetheBlochMuonITS,0.1,10.0,0);
2678     betheBlochPion = new TF1("betheBlochPion",BetheBlochPionITS,0.1,10.0,0);
2679     betheBlochKaon = new TF1("betheBlochKaon",BetheBlochKaonITS,0.1,10.0,0);
2680     betheBlochProton = new TF1("betheBlochProton",BetheBlochProtonITS,0.1,10.0,0);
2681
2682   }
2683   else {
2684
2685     betheBlochElectron = new TF1("betheBlochElectron",BetheBlochElectronTPC,0.1,10.0,0);
2686     betheBlochMuon = new TF1("betheBlochMuon",BetheBlochMuonTPC,0.1,10.0,0);
2687     betheBlochPion = new TF1("betheBlochPion",BetheBlochPionTPC,0.1,10.0,0);
2688     betheBlochKaon = new TF1("betheBlochKaon",BetheBlochKaonTPC,0.1,10.0,0);
2689     betheBlochProton = new TF1("betheBlochProton",BetheBlochProtonTPC,0.1,10.0,0);
2690
2691   }
2692
2693
2694   if((hsignalPidPartner0) || (hsignalPidPartner1)) {
2695     TCanvas * cPidSignal =new TCanvas("cPidSignal","cPidSignal",800,800);
2696     cPidSignal->Divide(2,1);
2697     cPidSignal->cd(1);
2698     if(hsignalPidPartner0) hsignalPidPartner0->Draw("colz");
2699     if(betheBlochElectron) betheBlochElectron->Draw("same");
2700     if(betheBlochMuon) betheBlochMuon->Draw("same");
2701     if(betheBlochPion) betheBlochPion->Draw("same");
2702     if(betheBlochKaon) betheBlochKaon->Draw("same");
2703     if(betheBlochProton) betheBlochProton->Draw("same");
2704     cPidSignal->cd(2);
2705     if(hsignalPidPartner1) hsignalPidPartner1->Draw("colz");
2706     if(betheBlochElectron) betheBlochElectron->Draw("same");
2707     if(betheBlochMuon) betheBlochMuon->Draw("same");
2708     if(betheBlochPion) betheBlochPion->Draw("same");
2709     if(betheBlochKaon) betheBlochKaon->Draw("same");
2710     if(betheBlochProton) betheBlochProton->Draw("same");
2711   }
2712
2713   THnSparseF *hsSparseITSsignal = dynamic_cast<THnSparseF *>(fList->FindObject("SparseITSsignal"));
2714   if(hsSparseITSsignal) {
2715
2716
2717     TH2D *sddsdd = hsSparseITSsignal->Projection(1,2);
2718     TH2D *ssdssd = hsSparseITSsignal->Projection(3,4);
2719     TH2D *sddssda = hsSparseITSsignal->Projection(1,3);
2720     TH2D *sddssdb = hsSparseITSsignal->Projection(2,4);
2721     TH2D *sddssdc = hsSparseITSsignal->Projection(1,4);
2722     TH2D *sddssdd = hsSparseITSsignal->Projection(2,3);
2723
2724     TCanvas * cITSSignal =new TCanvas("cITSSignal","cITSSignal",800,800);
2725     cITSSignal->Divide(2,3);
2726     cITSSignal->cd(1);
2727     sddsdd->Draw("colz");
2728     cITSSignal->cd(2);
2729     ssdssd->Draw("colz");
2730     cITSSignal->cd(3);
2731     sddssda->Draw("colz");
2732     cITSSignal->cd(4);
2733     sddssdb->Draw("colz");
2734     cITSSignal->cd(5);
2735     sddssdc->Draw("colz");
2736     cITSSignal->cd(6);
2737     sddssdd->Draw("colz");
2738
2739   }
2740
2741   THnSparseF *hsSparseITSsignalSplit = dynamic_cast<THnSparseF *>(fList->FindObject("SparseITSsignalSplit"));
2742   if(hsSparseITSsignalSplit) {
2743
2744     // no splitted
2745     hsSparseITSsignalSplit->GetAxis(0)->SetRange(1,1);
2746
2747     TH1D *layerITS2 = hsSparseITSsignalSplit->Projection(1);
2748     TH1D *layerITS3 = hsSparseITSsignalSplit->Projection(2);
2749     TH1D *layerITS4 = hsSparseITSsignalSplit->Projection(3);
2750     TH1D *layerITS5 = hsSparseITSsignalSplit->Projection(4);
2751
2752     // splitted
2753     hsSparseITSsignalSplit->GetAxis(0)->SetRange(2,2);
2754
2755     TH1D *layerITS2s = hsSparseITSsignalSplit->Projection(1);
2756     TH1D *layerITS3s = hsSparseITSsignalSplit->Projection(2);
2757     TH1D *layerITS4s = hsSparseITSsignalSplit->Projection(3);
2758     TH1D *layerITS5s = hsSparseITSsignalSplit->Projection(4);
2759
2760     TCanvas * cITSSignalSplit =new TCanvas("cITSSignalSplit","cITSSignalSplit",800,800);
2761     cITSSignalSplit->Divide(2,2);
2762     cITSSignalSplit->cd(1);
2763     layerITS2->Draw();
2764     layerITS2s->Draw("same");
2765     TLegend *legITS2 = new TLegend(0.4,0.6,0.89,0.89);
2766     legITS2->AddEntry(layerITS2,"No splitted","p");
2767     legITS2->AddEntry(layerITS2s,"Splitted","p");
2768     legITS2->Draw("same");
2769     cITSSignalSplit->cd(2);
2770     layerITS3->Draw();
2771     layerITS3s->Draw("same");
2772     TLegend *legITS3 = new TLegend(0.4,0.6,0.89,0.89);
2773     legITS3->AddEntry(layerITS3,"No splitted","p");
2774     legITS3->AddEntry(layerITS3s,"Splitted","p");
2775     legITS3->Draw("same");
2776     cITSSignalSplit->cd(3);
2777     layerITS4->Draw();
2778     layerITS4s->Draw("same");
2779     TLegend *legITS4 = new TLegend(0.4,0.6,0.89,0.89);
2780     legITS4->AddEntry(layerITS4,"No splitted","p");
2781     legITS4->AddEntry(layerITS4s,"Splitted","p");
2782     legITS4->Draw("same");
2783     cITSSignalSplit->cd(4);
2784     layerITS5->Draw();
2785     layerITS5s->Draw("same");
2786     TLegend *legITS5 = new TLegend(0.4,0.6,0.89,0.89);
2787     legITS5->AddEntry(layerITS5,"No splitted","p");
2788     legITS5->AddEntry(layerITS5s,"Splitted","p");
2789     legITS5->Draw("same");
2790
2791
2792   }
2793
2794 }
2795 //_____________________________________________________________________________
2796 Double_t AliHFEelecbackground::BetheBlochElectronITS(const Double_t *x, const Double_t * /*par*/)
2797 {
2798   //
2799   // Bethe Bloch for ITS
2800   //
2801   static AliITSPIDResponse itsPidResponse;
2802   return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(0));
2803 }
2804 //_____________________________________________________________________________
2805 Double_t AliHFEelecbackground::BetheBlochMuonITS(const Double_t *x, const Double_t * /*par*/)
2806 {
2807   //
2808   // Bethe Bloch for ITS
2809   //
2810   static AliITSPIDResponse itsPidResponse;
2811   return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(1));
2812 }
2813 //_____________________________________________________________________________
2814 Double_t AliHFEelecbackground::BetheBlochPionITS(const Double_t *x, const Double_t * /*par*/)
2815 {
2816   //
2817   // Bethe Bloch for ITS
2818   //
2819   static AliITSPIDResponse itsPidResponse;
2820   return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(2));
2821 }
2822 //_____________________________________________________________________________
2823 Double_t AliHFEelecbackground::BetheBlochKaonITS(const Double_t *x, const Double_t * /*par*/)
2824 {
2825   //
2826   // Bethe Bloch for ITS
2827   //
2828   static AliITSPIDResponse itsPidResponse;
2829   return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(3));
2830 }
2831 //_____________________________________________________________________________
2832 Double_t AliHFEelecbackground::BetheBlochProtonITS(const Double_t *x, const Double_t * /*par*/)
2833 {
2834   //
2835   // Bethe Bloch for ITS
2836   //
2837   static AliITSPIDResponse itsPidResponse;
2838   return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(4));
2839 }
2840 //_____________________________________________________________________________
2841 Double_t AliHFEelecbackground::BetheBlochElectronTPC(const Double_t *x, const Double_t * /*par*/)
2842 {
2843   //
2844   // Bethe Bloch for TPC
2845   //
2846   static AliTPCPIDResponse tpcPidResponse;
2847   return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kElectron);
2848 }
2849 //_____________________________________________________________________________
2850 Double_t AliHFEelecbackground::BetheBlochMuonTPC(const Double_t *x, const Double_t * /*par*/)
2851 {
2852   //
2853   // Bethe Bloch for TPC
2854   //
2855   static AliTPCPIDResponse tpcPidResponse;
2856   return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kMuon);
2857 }
2858 //_____________________________________________________________________________
2859 Double_t AliHFEelecbackground::BetheBlochPionTPC(const Double_t *x, const Double_t * /*par*/)
2860 {
2861   //
2862   // Bethe Bloch for TPC
2863   //
2864   static AliTPCPIDResponse tpcPidResponse;
2865   return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kPion);
2866 }
2867 //_____________________________________________________________________________
2868 Double_t AliHFEelecbackground::BetheBlochKaonTPC(const Double_t *x, const Double_t * /*par*/)
2869 {
2870   //
2871   // Bethe Bloch for TPC
2872   //
2873   static AliTPCPIDResponse tpcPidResponse;
2874   return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kKaon);
2875 }
2876 //_____________________________________________________________________________
2877 Double_t AliHFEelecbackground::BetheBlochProtonTPC(const Double_t *x, const Double_t * /*par*/)
2878 {
2879   //
2880   // Bethe Bloch for TPC
2881   //
2882   static AliTPCPIDResponse tpcPidResponse;
2883   return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kProton);
2884 }
2885
2886
2887