PWG3/hfe/AliHFEelecbackground.cxx

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