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