PWG3/hfe/AliHFEelecbackground.cxx

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