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   ////////////////////////
2248   // Cut efficiencies
2249   ////////////////////////
2250
2251   THnSparseF *hsSparseMCe = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC"));
2252
2253   if(hsSparseMCe) {
2254
2255     // init histos
2256     TAxis *axissources = hsSparseMCe->GetAxis(2);
2257     Int_t  nbsources = axissources->GetNbins();
2258     TAxis *axiscuts = hsSparseMCe->GetAxis(1);
2259     Int_t  nbcuts = axiscuts->GetNbins();
2260     TH1D **histopassedcuts = new TH1D*[nbsources*nbcuts];
2261     Double_t *nbEntriesCuts = new Double_t[nbsources*nbcuts];
2262
2263     //printf("Number of cuts %d\n",nbcuts);
2264
2265     // canvas
2266     TCanvas * chsSparseMCeeff =new TCanvas("hsSparseMCeeffDebug","hsSparseMCeeffDebug",800,800);
2267     chsSparseMCeeff->Divide(3,1);
2268
2269     // histos
2270     for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
2271       hsSparseMCe->GetAxis(2)->SetRange(sourceid+1,sourceid+1);
2272       for(Int_t cut = 0; cut < nbcuts; cut++){
2273         hsSparseMCe->GetAxis(1)->SetRange(cut+1,cut+1);
2274         histopassedcuts[sourceid*nbcuts+cut] = hsSparseMCe->Projection(0);
2275         hsSparseMCe->GetAxis(1)->SetRange(1,hsSparseMCe->GetAxis(1)->GetNbins());
2276       }
2277       hsSparseMCe->GetAxis(2)->SetRange(1,hsSparseMCe->GetAxis(2)->GetNbins());
2278     }
2279
2280     // calcul efficiencies
2281     ///////////////////////
2282     // histos
2283     for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
2284       // Next is compared to the partner tracked
2285       for(Int_t cut = 2; cut < nbcuts; cut++){
2286         nbEntriesCuts[sourceid*nbcuts+cut] = histopassedcuts[sourceid*nbcuts+cut]->GetEntries();
2287         if(histopassedcuts[sourceid*nbcuts+1]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+cut]->Divide(histopassedcuts[sourceid*nbcuts+1]);
2288       }
2289       // First one is if the partner is tracked.
2290       nbEntriesCuts[sourceid*nbcuts+1] = histopassedcuts[sourceid*nbcuts+1]->GetEntries();
2291       if(histopassedcuts[sourceid*nbcuts]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+1]->Divide(histopassedcuts[sourceid*nbcuts]);
2292       // First one is input
2293       nbEntriesCuts[sourceid*nbcuts] = histopassedcuts[sourceid*nbcuts]->GetEntries();
2294     }
2295
2296     /////////////
2297     // ratios
2298     ////////////
2299     for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
2300       for(Int_t cut = 1; cut < nbcuts; cut++){
2301         if(nbEntriesCuts[sourceid*nbcuts] > 0.0) nbEntriesCuts[sourceid*nbcuts+cut] = nbEntriesCuts[sourceid*nbcuts+cut]/nbEntriesCuts[sourceid*nbcuts];
2302       }
2303     }
2304     TH1F *ratioHistoEntriesGamma = new TH1F("ratioHistoEntriesGamma","", nbcuts-1, 0.0, nbcuts-1.0);
2305     TH1F *ratioHistoEntriesPi0 = new TH1F("ratioHistoEntriesPi0","", nbcuts-1, 0.0, nbcuts-1.0);
2306     TH1F *ratioHistoEntriesC = new TH1F("ratioHistoEntriesC","", nbcuts-1, 0.0, nbcuts-1.0);
2307     for(Int_t k = 1; k < nbcuts; k++){
2308       ratioHistoEntriesGamma->SetBinContent(k,nbEntriesCuts[nbcuts+k]);
2309       ratioHistoEntriesPi0->SetBinContent(k,nbEntriesCuts[2*nbcuts+k]);
2310       ratioHistoEntriesC->SetBinContent(k,nbEntriesCuts[4*nbcuts+k]);
2311     }
2312     //
2313     TAxis *xAxisGamma = ratioHistoEntriesGamma->GetXaxis();
2314     xAxisGamma->SetBinLabel(1,"Partner tracked");
2315     xAxisGamma->SetBinLabel(2,"Opposite sign");
2316     xAxisGamma->SetBinLabel(3,"Single Track Cut");
2317     xAxisGamma->SetBinLabel(4,"Shared Clusters");
2318     xAxisGamma->SetBinLabel(5,"PID");
2319     xAxisGamma->SetBinLabel(6,"DCA");
2320     xAxisGamma->SetBinLabel(7,"Chi^{2}/Ndf");
2321     xAxisGamma->SetBinLabel(8,"Opening angle");
2322     xAxisGamma->SetBinLabel(9,"Invariant mass");
2323     //
2324     TAxis *xAxisPi0 = ratioHistoEntriesPi0->GetXaxis();
2325     xAxisPi0->SetBinLabel(1,"Partner tracked");
2326     xAxisPi0->SetBinLabel(2,"Opposite sign");
2327     xAxisPi0->SetBinLabel(3,"Single Track Cut");
2328     xAxisPi0->SetBinLabel(4,"Shared Clusters");
2329     xAxisPi0->SetBinLabel(5,"PID");
2330     xAxisPi0->SetBinLabel(6,"DCA");
2331     xAxisPi0->SetBinLabel(7,"Chi^{2}/Ndf");
2332     xAxisPi0->SetBinLabel(8,"Opening angle");
2333     xAxisPi0->SetBinLabel(9,"Invariant mass");
2334     //
2335     TAxis *xAxisC = ratioHistoEntriesC->GetXaxis();
2336     xAxisC->SetBinLabel(1,"Partner tracked");
2337     xAxisC->SetBinLabel(2,"Opposite sign");
2338     xAxisC->SetBinLabel(3,"Single Track Cut");
2339     xAxisC->SetBinLabel(4,"Shared Clusters");
2340     xAxisC->SetBinLabel(5,"PID");
2341     xAxisC->SetBinLabel(6,"DCA");
2342     xAxisC->SetBinLabel(7,"Chi^{2}/Ndf");
2343     xAxisC->SetBinLabel(8,"Opening angle");
2344     xAxisC->SetBinLabel(9,"Invariant mass");
2345     //
2346     TCanvas * cRatioHistoEntries =new TCanvas("cRatioHistoEntries","cRatioHistoEntries",800,800);
2347     cRatioHistoEntries->cd(1);
2348     ratioHistoEntriesGamma->SetStats(0);
2349     ratioHistoEntriesGamma->Draw();
2350     ratioHistoEntriesPi0->SetStats(0);
2351     ratioHistoEntriesPi0->Draw("same");
2352     ratioHistoEntriesC->SetStats(0);
2353     //ratioHistoEntriesC->Draw("same");
2354     TLegend *legEntries = new TLegend(0.4,0.6,0.89,0.89);
2355     legEntries->AddEntry(ratioHistoEntriesGamma,"#gamma","l");
2356     legEntries->AddEntry(ratioHistoEntriesPi0,"#pi^{0}","l");
2357     //legEntries->AddEntry(ratioHistoEntriesC,"c","p");
2358     legEntries->Draw("same");
2359
2360     ////////////////////
2361     // plot Debug
2362     ///////////////////
2363     Int_t source = 1;
2364     chsSparseMCeeff->cd(1);
2365     histopassedcuts[source*nbcuts+0]->SetTitle("#gamma");
2366     histopassedcuts[source*nbcuts+1]->SetTitle("#gamma");
2367     histopassedcuts[source*nbcuts+2]->SetTitle("#gamma");
2368     histopassedcuts[source*nbcuts+3]->SetTitle("#gamma");
2369     histopassedcuts[source*nbcuts+4]->SetTitle("#gamma");
2370     histopassedcuts[source*nbcuts+5]->SetTitle("#gamma");
2371     histopassedcuts[source*nbcuts+6]->SetTitle("#gamma");
2372     histopassedcuts[source*nbcuts+7]->SetTitle("#gamma");
2373     histopassedcuts[source*nbcuts+8]->SetTitle("#gamma");
2374     histopassedcuts[source*nbcuts+9]->SetTitle("#gamma");
2375     //histopassedcuts[source*nbcuts+0]->SetStats(0);
2376     histopassedcuts[source*nbcuts+1]->SetStats(0);
2377     histopassedcuts[source*nbcuts+2]->SetStats(0);
2378     histopassedcuts[source*nbcuts+3]->SetStats(0);
2379     histopassedcuts[source*nbcuts+4]->SetStats(0);
2380     histopassedcuts[source*nbcuts+5]->SetStats(0);
2381     histopassedcuts[source*nbcuts+6]->SetStats(0);
2382     histopassedcuts[source*nbcuts+7]->SetStats(0);
2383     histopassedcuts[source*nbcuts+8]->SetStats(0);
2384     histopassedcuts[source*nbcuts+9]->SetStats(0);
2385     //histopassedcuts[source*nbcuts+0]->Draw();
2386     //histopassedcuts[source*nbcuts+1]->Draw("");
2387     histopassedcuts[source*nbcuts+2]->Draw();
2388     histopassedcuts[source*nbcuts+3]->Draw("same");
2389     //histopassedcuts[source*nbcuts+4]->Draw("same");
2390     histopassedcuts[source*nbcuts+5]->Draw("same");
2391     histopassedcuts[source*nbcuts+6]->Draw("same");
2392     //histopassedcuts[source*nbcuts+7]->Draw("same");
2393     histopassedcuts[source*nbcuts+8]->Draw("same");
2394     histopassedcuts[source*nbcuts+9]->Draw("same");
2395     TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
2396     //legb->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
2397     //legb->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
2398     legb->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
2399     legb->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
2400     //legb->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
2401     legb->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
2402     legb->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
2403     //legb->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
2404     legb->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
2405     legb->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
2406     legb->Draw("same");
2407
2408     source = 2;
2409     chsSparseMCeeff->cd(2);
2410     histopassedcuts[source*nbcuts+0]->SetTitle("#pi^{0}");
2411     histopassedcuts[source*nbcuts+1]->SetTitle("#pi^{0}");
2412     histopassedcuts[source*nbcuts+2]->SetTitle("#pi^{0}");
2413     histopassedcuts[source*nbcuts+3]->SetTitle("#pi^{0}");
2414     histopassedcuts[source*nbcuts+4]->SetTitle("#pi^{0}");
2415     histopassedcuts[source*nbcuts+5]->SetTitle("#pi^{0}");
2416     histopassedcuts[source*nbcuts+6]->SetTitle("#pi^{0}");
2417     histopassedcuts[source*nbcuts+7]->SetTitle("#pi^{0}");
2418     histopassedcuts[source*nbcuts+8]->SetTitle("#pi^{0}");
2419     histopassedcuts[source*nbcuts+9]->SetTitle("#pi^{0}");
2420     //histopassedcuts[source*nbcuts+0]->SetStats(0);
2421     histopassedcuts[source*nbcuts+1]->SetStats(0);
2422     histopassedcuts[source*nbcuts+2]->SetStats(0);
2423     histopassedcuts[source*nbcuts+3]->SetStats(0);
2424     histopassedcuts[source*nbcuts+4]->SetStats(0);
2425     histopassedcuts[source*nbcuts+5]->SetStats(0);
2426     histopassedcuts[source*nbcuts+6]->SetStats(0);
2427     histopassedcuts[source*nbcuts+7]->SetStats(0);
2428     histopassedcuts[source*nbcuts+8]->SetStats(0);
2429     histopassedcuts[source*nbcuts+9]->SetStats(0);
2430     //histopassedcuts[source*nbcuts+0]->Draw();
2431     //histopassedcuts[source*nbcuts+1]->Draw();
2432     histopassedcuts[source*nbcuts+2]->Draw();
2433     histopassedcuts[source*nbcuts+3]->Draw("same");
2434     //histopassedcuts[source*nbcuts+4]->Draw("same");
2435     histopassedcuts[source*nbcuts+5]->Draw("same");
2436     histopassedcuts[source*nbcuts+6]->Draw("same");
2437     //histopassedcuts[source*nbcuts+7]->Draw("same");
2438     histopassedcuts[source*nbcuts+8]->Draw("same");
2439     histopassedcuts[source*nbcuts+9]->Draw("same");
2440     TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
2441     //legc->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
2442     //legc->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
2443     legc->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
2444     legc->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
2445     //legc->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
2446     legc->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
2447     legc->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
2448     //legc->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
2449     legc->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
2450     legc->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
2451     legc->Draw("same");
2452
2453     source = 4;
2454     chsSparseMCeeff->cd(3);
2455     histopassedcuts[source*nbcuts+0]->SetTitle("C");
2456     histopassedcuts[source*nbcuts+1]->SetTitle("C");
2457     histopassedcuts[source*nbcuts+2]->SetTitle("C");
2458     histopassedcuts[source*nbcuts+3]->SetTitle("C");
2459     histopassedcuts[source*nbcuts+4]->SetTitle("C");
2460     histopassedcuts[source*nbcuts+5]->SetTitle("C");
2461     histopassedcuts[source*nbcuts+6]->SetTitle("C");
2462     histopassedcuts[source*nbcuts+7]->SetTitle("C");
2463     histopassedcuts[source*nbcuts+8]->SetTitle("C");
2464     histopassedcuts[source*nbcuts+9]->SetTitle("C");
2465     //histopassedcuts[source*nbcuts+0]->SetStats(0);
2466     histopassedcuts[source*nbcuts+1]->SetStats(0);
2467     histopassedcuts[source*nbcuts+2]->SetStats(0);
2468     histopassedcuts[source*nbcuts+3]->SetStats(0);
2469     histopassedcuts[source*nbcuts+4]->SetStats(0);
2470     histopassedcuts[source*nbcuts+5]->SetStats(0);
2471     histopassedcuts[source*nbcuts+6]->SetStats(0);
2472     histopassedcuts[source*nbcuts+7]->SetStats(0);
2473     histopassedcuts[source*nbcuts+8]->SetStats(0);
2474     histopassedcuts[source*nbcuts+9]->SetStats(0);
2475     //histopassedcuts[source*nbcuts+0]->Draw();
2476     //histopassedcuts[source*nbcuts+1]->Draw();
2477     histopassedcuts[source*nbcuts+2]->Draw();
2478     histopassedcuts[source*nbcuts+3]->Draw("same");
2479     //histopassedcuts[source*nbcuts+4]->Draw("same");
2480     histopassedcuts[source*nbcuts+5]->Draw("same");
2481     histopassedcuts[source*nbcuts+6]->Draw("same");
2482     //histopassedcuts[source*nbcuts+7]->Draw("same");
2483     histopassedcuts[source*nbcuts+8]->Draw("same");
2484     histopassedcuts[source*nbcuts+9]->Draw("same");
2485     TLegend *lege = new TLegend(0.4,0.6,0.89,0.89);
2486     //lege->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
2487     //lege->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
2488     lege->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
2489     lege->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
2490     //lege->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
2491     lege->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
2492     lege->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
2493     //lege->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
2494     lege->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
2495     lege->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
2496     lege->Draw("same");
2497
2498     //////////////////////
2499     // Input
2500     //////////////////////
2501
2502     TCanvas * chsSparseMCein =new TCanvas("hsSparseMCeinput","hsSparseMCeinput",800,800);
2503     chsSparseMCein->cd(1);
2504     Double_t nbGamma = 0.0;
2505     source = 1;
2506     nbGamma = histopassedcuts[source*nbcuts+0]->GetEntries();
2507     histopassedcuts[source*nbcuts+0]->SetStats(0);
2508     histopassedcuts[source*nbcuts+0]->Draw();
2509     TLegend *leginput = new TLegend(0.4,0.6,0.89,0.89);
2510     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#gamma","p");
2511     Double_t nbPi0 = 0.0;
2512     source = 2;
2513     nbPi0 = histopassedcuts[source*nbcuts+0]->GetEntries();
2514     histopassedcuts[source*nbcuts+0]->SetStats(0);
2515     histopassedcuts[source*nbcuts+0]->Draw("same");
2516     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#pi^{0}","p");
2517     Double_t nbEta = 0.0;
2518     source = 3;
2519     nbEta = histopassedcuts[source*nbcuts+0]->GetEntries();
2520     histopassedcuts[source*nbcuts+0]->SetStats(0);
2521     histopassedcuts[source*nbcuts+0]->Draw("same");
2522     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#eta","p");
2523     Double_t nbC = 0.0;
2524     source = 4;
2525     nbC = histopassedcuts[source*nbcuts+0]->GetEntries();
2526     histopassedcuts[source*nbcuts+0]->SetStats(0);
2527     histopassedcuts[source*nbcuts+0]->Draw("same");
2528     leginput->AddEntry(histopassedcuts[source*nbcuts+0],"c","p");
2529     leginput->Draw("same");
2530
2531     //printf("Gamma %f, pi^{0} %f and #eta %f, c %f\n",nbGamma,nbPi0,nbEta,nbC);
2532
2533     //////////////////////
2534     // Tracked
2535     //////////////////////
2536
2537     TCanvas * cTracked = new TCanvas("cTracked","cTracked",800,800);
2538     cTracked->cd(1);
2539     source = 1;
2540     histopassedcuts[source*nbcuts+1]->Draw();
2541     TLegend *legTracked = new TLegend(0.4,0.6,0.89,0.89);
2542     legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#gamma","p");
2543     source = 2;
2544     histopassedcuts[source*nbcuts+1]->Draw("same");
2545     legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#pi^{0}","p");
2546     legTracked->Draw("same");
2547
2548   }
2549
2550   /////////////////////////////////////
2551   // Data Radius and chi2Ndf if AliKF
2552   ////////////////////////////////////
2553
2554   TH1F *hDataRadius = dynamic_cast<TH1F *>(fList->FindObject("DataRadius"));
2555   TH1F *hDataChi2Ndf = dynamic_cast<TH1F *>(fList->FindObject("DataChi2Ndf"));
2556
2557   if(hDataRadius || hDataChi2Ndf) {
2558     TCanvas * cDataRadiusChi2Ndf =new TCanvas("CanvasDataRadiusChi2Ndf","CanvasDataRadiusChi2Ndf",800,800);
2559     cDataRadiusChi2Ndf->Divide(2,1);
2560     cDataRadiusChi2Ndf->cd(1);
2561     if(hDataRadius) hDataRadius->Draw();
2562     cDataRadiusChi2Ndf->cd(2);
2563     if(hDataChi2Ndf) hDataChi2Ndf->Draw();
2564   }
2565
2566   ///////////////////////
2567   // Data DCA
2568   //////////////////////
2569
2570   TH1F *hDataDCA = dynamic_cast<TH1F *>(fList->FindObject("DataDCA"));
2571
2572   if(hDataDCA) {
2573     TCanvas * cDataDCA =new TCanvas("CanvasDataDCA","CanvasDataDCA",800,800);
2574     cDataDCA->cd(1);
2575     hDataDCA->Draw();
2576   }
2577
2578   /////////////////////////////////////
2579   // MC Radius and chi2Ndf if AliKF
2580   ////////////////////////////////////
2581
2582   TH2F *hMCRadius = dynamic_cast<TH2F *>(fList->FindObject("MCRadius"));
2583   TH2F *hMCChi2Ndf = dynamic_cast<TH2F *>(fList->FindObject("MCChi2Ndf"));
2584
2585   if(hMCRadius || hMCChi2Ndf) {
2586     TCanvas * cMCRadiusChi2Ndf =new TCanvas("CanvasMCRadiusChi2Ndf","CanvasMCRadiusChi2Ndf",800,800);
2587     cMCRadiusChi2Ndf->Divide(2,1);
2588     cMCRadiusChi2Ndf->cd(1);
2589     //TH1D *hMCRadiusBackground = hMCRadius->ProjectionX("MCRadiusBackGround",1,1,"e");
2590     TH1D *hMCRadiusGamma = hMCRadius->ProjectionX("MCRadiusGamma",2,2,"e");
2591     TH1D *hMCRadiusPi0 = hMCRadius->ProjectionX("MCRadiusPi0",3,3,"e");
2592     TH1D *hMCRadiusEta = hMCRadius->ProjectionX("MCRadiusEta",4,4,"e");
2593     TH1D *hMCRadiusC = hMCRadius->ProjectionX("MCRadiusC",5,5,"e");
2594     TH1D *hMCRadiusB = hMCRadius->ProjectionX("MCRadiusB",6,6,"e");
2595     //hMCRadiusBackground->Draw();
2596     hMCRadiusGamma->Draw();
2597     hMCRadiusPi0->Draw("same");
2598     hMCRadiusEta->Draw("same");
2599     hMCRadiusC->Draw("same");
2600     hMCRadiusB->Draw("same");
2601     TLegend *legRadius = new TLegend(0.4,0.6,0.89,0.89);
2602     //legRadius->AddEntry(hMCRadiusBackground,"Background","p");
2603     legRadius->AddEntry(hMCRadiusGamma,"#gamma","p");
2604     legRadius->AddEntry(hMCRadiusPi0,"#pi^{0}","p");
2605     legRadius->AddEntry(hMCRadiusEta,"#eta","p");
2606     legRadius->AddEntry(hMCRadiusC,"c","p");
2607     legRadius->AddEntry(hMCRadiusB,"b","p");
2608     legRadius->Draw("same");
2609     cMCRadiusChi2Ndf->cd(2);
2610     //TH1D *hMCChi2NdfBackground = hMCChi2Ndf->ProjectionX("MCChi2NdfBackGround",1,1,"e");
2611     TH1D *hMCChi2NdfGamma = hMCChi2Ndf->ProjectionX("MCChi2NdfGamma",2,2,"e");
2612     TH1D *hMCChi2NdfPi0 = hMCChi2Ndf->ProjectionX("MCChi2NdfPi0",3,3,"e");
2613     TH1D *hMCChi2NdfEta = hMCChi2Ndf->ProjectionX("MCChi2NdfEta",4,4,"e");
2614     TH1D *hMCChi2NdfC = hMCChi2Ndf->ProjectionX("MCChi2NdfC",5,5,"e");
2615     TH1D *hMCChi2NdfB = hMCChi2Ndf->ProjectionX("MCChi2NdfB",6,6,"e");
2616     //hMCChi2NdfBackground->Draw();
2617     hMCChi2NdfGamma->Draw();
2618     hMCChi2NdfPi0->Draw("same");
2619     hMCChi2NdfEta->Draw("same");
2620     hMCChi2NdfC->Draw("same");
2621     hMCChi2NdfB->Draw("same");
2622     TLegend *legChi2Ndf = new TLegend(0.4,0.6,0.89,0.89);
2623     //legChi2Ndf->AddEntry(hMCChi2NdfBackground,"Background","p");
2624     legChi2Ndf->AddEntry(hMCChi2NdfGamma,"#gamma","p");
2625     legChi2Ndf->AddEntry(hMCChi2NdfPi0,"#pi^{0}","p");
2626     legChi2Ndf->AddEntry(hMCChi2NdfEta,"#eta","p");
2627     legChi2Ndf->AddEntry(hMCChi2NdfC,"c","p");
2628     legChi2Ndf->AddEntry(hMCChi2NdfB,"b","p");
2629     legChi2Ndf->Draw("same");
2630   }
2631
2632   ///////////////////////
2633   // MC DCA
2634   //////////////////////
2635
2636   TH2F *hMCDCA = dynamic_cast<TH2F *>(fList->FindObject("MCDCA"));
2637
2638   if(hMCDCA) {
2639     TCanvas * cMCDCA =new TCanvas("CanvasMCDCA","CanvasMCDCA",800,800);
2640     cMCDCA->cd(1);
2641     //TH1D *hMCDCABackground = hMCDCA->ProjectionX("MCDCABackGround",1,1,"e");
2642     TH1D *hMCDCAGamma = hMCDCA->ProjectionX("MCDCAGamma",2,2,"e");
2643     TH1D *hMCDCAPi0 = hMCDCA->ProjectionX("MCDCAPi0",3,3,"e");
2644     TH1D *hMCDCAEta = hMCDCA->ProjectionX("MCDCAEta",4,4,"e");
2645     TH1D *hMCDCAC = hMCDCA->ProjectionX("MCDCAC",5,5,"e");
2646     TH1D *hMCDCAB = hMCDCA->ProjectionX("MCDCAB",6,6,"e");
2647     //hMCDCABackground->Draw();
2648     hMCDCAGamma->Draw();
2649     hMCDCAPi0->Draw("same");
2650     hMCDCAEta->Draw("same");
2651     hMCDCAC->Draw("same");
2652     hMCDCAB->Draw("same");
2653     TLegend *legDCA = new TLegend(0.4,0.6,0.89,0.89);
2654     //legDCA->AddEntry(hMCDCABackground,"Background","p");
2655     legDCA->AddEntry(hMCDCAGamma,"#gamma","p");
2656     legDCA->AddEntry(hMCDCAPi0,"#pi^{0}","p");
2657     legDCA->AddEntry(hMCDCAEta,"#eta","p");
2658     legDCA->AddEntry(hMCDCAC,"c","p");
2659     legDCA->AddEntry(hMCDCAB,"b","p");
2660     legDCA->Draw("same");
2661   }
2662
2663
2664   /////////////////////////
2665   // PID Partner Signal
2666   /////////////////////////
2667   TF1 *betheBlochElectron = 0x0;
2668   TF1 *betheBlochMuon = 0x0;
2669   TF1 *betheBlochPion = 0x0;
2670   TF1 *betheBlochKaon = 0x0;
2671   TF1 *betheBlochProton = 0x0;
2672
2673   TH2F *hsignalPidPartner0 = dynamic_cast<TH2F *>(fList->FindObject("TPCPartner0"));
2674   TH2F *hsignalPidPartner1 = dynamic_cast<TH2F *>(fList->FindObject("TPCPartner1"));
2675   if((!hsignalPidPartner0) && (!hsignalPidPartner1)) {
2676     hsignalPidPartner0 = dynamic_cast<TH2F *>(fList->FindObject("ITSPartner0"));
2677     hsignalPidPartner1 = dynamic_cast<TH2F *>(fList->FindObject("ITSPartner1"));
2678
2679     betheBlochElectron = new TF1("betheBlochElectron",BetheBlochElectronITS,0.1,10.0,0);
2680     betheBlochMuon = new TF1("betheBlochMuon",BetheBlochMuonITS,0.1,10.0,0);
2681     betheBlochPion = new TF1("betheBlochPion",BetheBlochPionITS,0.1,10.0,0);
2682     betheBlochKaon = new TF1("betheBlochKaon",BetheBlochKaonITS,0.1,10.0,0);
2683     betheBlochProton = new TF1("betheBlochProton",BetheBlochProtonITS,0.1,10.0,0);
2684
2685   }
2686   else {
2687
2688     betheBlochElectron = new TF1("betheBlochElectron",BetheBlochElectronTPC,0.1,10.0,0);
2689     betheBlochMuon = new TF1("betheBlochMuon",BetheBlochMuonTPC,0.1,10.0,0);
2690     betheBlochPion = new TF1("betheBlochPion",BetheBlochPionTPC,0.1,10.0,0);
2691     betheBlochKaon = new TF1("betheBlochKaon",BetheBlochKaonTPC,0.1,10.0,0);
2692     betheBlochProton = new TF1("betheBlochProton",BetheBlochProtonTPC,0.1,10.0,0);
2693
2694   }
2695
2696
2697   if((hsignalPidPartner0) || (hsignalPidPartner1)) {
2698     TCanvas * cPidSignal =new TCanvas("cPidSignal","cPidSignal",800,800);
2699     cPidSignal->Divide(2,1);
2700     cPidSignal->cd(1);
2701     if(hsignalPidPartner0) hsignalPidPartner0->Draw("colz");
2702     if(betheBlochElectron) betheBlochElectron->Draw("same");
2703     if(betheBlochMuon) betheBlochMuon->Draw("same");
2704     if(betheBlochPion) betheBlochPion->Draw("same");
2705     if(betheBlochKaon) betheBlochKaon->Draw("same");
2706     if(betheBlochProton) betheBlochProton->Draw("same");
2707     cPidSignal->cd(2);
2708     if(hsignalPidPartner1) hsignalPidPartner1->Draw("colz");
2709     if(betheBlochElectron) betheBlochElectron->Draw("same");
2710     if(betheBlochMuon) betheBlochMuon->Draw("same");
2711     if(betheBlochPion) betheBlochPion->Draw("same");
2712     if(betheBlochKaon) betheBlochKaon->Draw("same");
2713     if(betheBlochProton) betheBlochProton->Draw("same");
2714   }
2715
2716   THnSparseF *hsSparseITSsignal = dynamic_cast<THnSparseF *>(fList->FindObject("SparseITSsignal"));
2717   if(hsSparseITSsignal) {
2718
2719
2720     TH2D *sddsdd = hsSparseITSsignal->Projection(1,2);
2721     TH2D *ssdssd = hsSparseITSsignal->Projection(3,4);
2722     TH2D *sddssda = hsSparseITSsignal->Projection(1,3);
2723     TH2D *sddssdb = hsSparseITSsignal->Projection(2,4);
2724     TH2D *sddssdc = hsSparseITSsignal->Projection(1,4);
2725     TH2D *sddssdd = hsSparseITSsignal->Projection(2,3);
2726
2727     TCanvas * cITSSignal =new TCanvas("cITSSignal","cITSSignal",800,800);
2728     cITSSignal->Divide(2,3);
2729     cITSSignal->cd(1);
2730     sddsdd->Draw("colz");
2731     cITSSignal->cd(2);
2732     ssdssd->Draw("colz");
2733     cITSSignal->cd(3);
2734     sddssda->Draw("colz");
2735     cITSSignal->cd(4);
2736     sddssdb->Draw("colz");
2737     cITSSignal->cd(5);
2738     sddssdc->Draw("colz");
2739     cITSSignal->cd(6);
2740     sddssdd->Draw("colz");
2741
2742   }
2743
2744   THnSparseF *hsSparseITSsignalSplit = dynamic_cast<THnSparseF *>(fList->FindObject("SparseITSsignalSplit"));
2745   if(hsSparseITSsignalSplit) {
2746
2747     // no splitted
2748     hsSparseITSsignalSplit->GetAxis(0)->SetRange(1,1);
2749
2750     TH1D *layerITS2 = hsSparseITSsignalSplit->Projection(1);
2751     TH1D *layerITS3 = hsSparseITSsignalSplit->Projection(2);
2752     TH1D *layerITS4 = hsSparseITSsignalSplit->Projection(3);
2753     TH1D *layerITS5 = hsSparseITSsignalSplit->Projection(4);
2754
2755     // splitted
2756     hsSparseITSsignalSplit->GetAxis(0)->SetRange(2,2);
2757
2758     TH1D *layerITS2s = hsSparseITSsignalSplit->Projection(1);
2759     TH1D *layerITS3s = hsSparseITSsignalSplit->Projection(2);
2760     TH1D *layerITS4s = hsSparseITSsignalSplit->Projection(3);
2761     TH1D *layerITS5s = hsSparseITSsignalSplit->Projection(4);
2762
2763     TCanvas * cITSSignalSplit =new TCanvas("cITSSignalSplit","cITSSignalSplit",800,800);
2764     cITSSignalSplit->Divide(2,2);
2765     cITSSignalSplit->cd(1);
2766     layerITS2->Draw();
2767     layerITS2s->Draw("same");
2768     TLegend *legITS2 = new TLegend(0.4,0.6,0.89,0.89);
2769     legITS2->AddEntry(layerITS2,"No splitted","p");
2770     legITS2->AddEntry(layerITS2s,"Splitted","p");
2771     legITS2->Draw("same");
2772     cITSSignalSplit->cd(2);
2773     layerITS3->Draw();
2774     layerITS3s->Draw("same");
2775     TLegend *legITS3 = new TLegend(0.4,0.6,0.89,0.89);
2776     legITS3->AddEntry(layerITS3,"No splitted","p");
2777     legITS3->AddEntry(layerITS3s,"Splitted","p");
2778     legITS3->Draw("same");
2779     cITSSignalSplit->cd(3);
2780     layerITS4->Draw();
2781     layerITS4s->Draw("same");
2782     TLegend *legITS4 = new TLegend(0.4,0.6,0.89,0.89);
2783     legITS4->AddEntry(layerITS4,"No splitted","p");
2784     legITS4->AddEntry(layerITS4s,"Splitted","p");
2785     legITS4->Draw("same");
2786     cITSSignalSplit->cd(4);
2787     layerITS5->Draw();
2788     layerITS5s->Draw("same");
2789     TLegend *legITS5 = new TLegend(0.4,0.6,0.89,0.89);
2790     legITS5->AddEntry(layerITS5,"No splitted","p");
2791     legITS5->AddEntry(layerITS5s,"Splitted","p");
2792     legITS5->Draw("same");
2793
2794
2795   }
2796
2797 }
2798 //_____________________________________________________________________________
2799 Double_t AliHFEelecbackground::BetheBlochElectronITS(const Double_t *x, const Double_t * /*par*/)
2800 {
2801   //
2802   // Bethe Bloch for ITS
2803   //
2804   static AliITSPIDResponse itsPidResponse;
2805   return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(0));
2806 }
2807 //_____________________________________________________________________________
2808 Double_t AliHFEelecbackground::BetheBlochMuonITS(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(1));
2815 }
2816 //_____________________________________________________________________________
2817 Double_t AliHFEelecbackground::BetheBlochPionITS(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(2));
2824 }
2825 //_____________________________________________________________________________
2826 Double_t AliHFEelecbackground::BetheBlochKaonITS(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(3));
2833 }
2834 //_____________________________________________________________________________
2835 Double_t AliHFEelecbackground::BetheBlochProtonITS(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(4));
2842 }
2843 //_____________________________________________________________________________
2844 Double_t AliHFEelecbackground::BetheBlochElectronTPC(const Double_t *x, const Double_t * /*par*/)
2845 {
2846   //
2847   // Bethe Bloch for TPC
2848   //
2849   static AliTPCPIDResponse tpcPidResponse;
2850   return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kElectron);
2851 }
2852 //_____________________________________________________________________________
2853 Double_t AliHFEelecbackground::BetheBlochMuonTPC(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::kMuon);
2860 }
2861 //_____________________________________________________________________________
2862 Double_t AliHFEelecbackground::BetheBlochPionTPC(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::kPion);
2869 }
2870 //_____________________________________________________________________________
2871 Double_t AliHFEelecbackground::BetheBlochKaonTPC(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::kKaon);
2878 }
2879 //_____________________________________________________________________________
2880 Double_t AliHFEelecbackground::BetheBlochProtonTPC(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::kProton);
2887 }
2888
2889
2890