PWGHF/hfe/AliHFEelecbackground.cxx

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