PWGHF/hfe/AliHFENonPhotonicElectron.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
  10  *      without fee, provided that the above copyright notice appears in all            *
  11  *      copies and that both the copyright notice and this permission notice            *
  12  *      appear in the supporting documentation. The authors make no claims              *
  13  *      about the suitability of this software for any purpose. It is                   *
  14  *      provided "as is" without express or implied warranty.                           *
  15  *                                                                                      *
  16  *************************************************************************************/
  17
  18 /*************************************************************************************
  19  *                                                                                      *
  20  *      Class for the Selection of Non-Heavy-Flavour-Electrons trought          *
  21  *      the invariant mass method. The selection can be done from two                   *
  22  *      different algorithms, which can be choosed calling the function         *
  23  *              "SetAlgorithm(TString Algorithm)".                                      *
  24  *                                                                                      *
  25  *              Authors: R.Bailhache, C.A.Schmidt                                       *
  26  *                                                                                      *
  27  *************************************************************************************/
  28
  29 #include "TVector2.h"
  30 #include "THnSparse.h"
  31 #include "TMath.h"
  32 #include "TLorentzVector.h"
  33 #include "TParticle.h"
  34 #include "TList.h"
  35 #include "TDatabasePDG.h"
  36
  37 #include "AliVEvent.h"
  38 #include "AliMCEvent.h"
  39 #include "AliESDEvent.h"
  40 #include "AliMCParticle.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODVertex.h"
  44 #include "AliAODTrack.h"
  45 #include "AliVTrack.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliPIDResponse.h"
  49 #include "AliPID.h"
  50
  51 #include "AliKFParticle.h"
  52 #include "AliKFVertex.h"
  53
  54 #include "AliHFEcuts.h"
  55 #include "AliHFEpid.h"
  56 #include "AliHFEpidQAmanager.h"
  57 #include "AliHFEtools.h"
  58 #include "AliHFEmcQA.h"
  59
  60 #include "AliHFENonPhotonicElectron.h"
  61
  62 ClassImp(AliHFENonPhotonicElectron)
  63 //________________________________________________________________________
  64 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const char *name, const Char_t *title)
  65   :TNamed               (name, title)
  66   ,fIsAOD               (kFALSE)
  67   ,fMCEvent             (NULL)
  68   ,fAODArrayMCInfo      (NULL)
  69   ,fLevelBack(-1)
  70   ,fHFEBackgroundCuts   (NULL)
  71   ,fPIDBackground       (0x0)
  72   ,fPIDBackgroundQA     (0)
  73   ,fkPIDRespons         (NULL)
  74   ,fPtBinning()
  75   ,fEtaBinning()
  76   ,fAlgorithmMA         (kTRUE)
  77   ,fChi2OverNDFCut      (3.0)
  78   ,fMaxDCA              (3.0)
  79 //  ,fMaxOpeningTheta   (0.02)
  80 //  ,fMaxOpeningPhi     (0.1)
  81   ,fMaxOpening3D        (TMath::Pi())
  82   ,fMaxInvMass          (1000)
  83   ,fSetMassConstraint   (kFALSE)
  84   ,fSelectCategory1tracks(kTRUE)
  85   ,fSelectCategory2tracks(kFALSE)
  86   ,fITSmeanShift(0.)
  87   ,fITSnSigmaHigh(3.)
  88   ,fITSnSigmaLow(-3.)
  89   ,fminPt(0.1)
  90   ,fArraytrack          (NULL)
  91   ,fCounterPoolBackground       (0)
  92   ,fnumberfound                 (0)
  93   ,fListOutput          (NULL)
  94   ,fAssElectron         (NULL)
  95   ,fIncElectron         (NULL)
  96   ,fUSign               (NULL)
  97   ,fLSign               (NULL)
  98   ,fUSmatches(NULL)
  99   ,fLSmatches(NULL)
 100   ,fHnsigmaITS(NULL)
 101   ,fWeightsSource(NULL)
 102   ,fIncElectronRadius(NULL)
 103   ,fRecElectronRadius(NULL)
 104 //  ,fUSignAngle        (NULL)
 105 //  ,fLSignAngle        (NULL)
 106 {
 107   //
 108   // Constructor
 109   //
 110   fPIDBackground   = new AliHFEpid("hfePidBackground");
 111   fPIDBackgroundQA = new AliHFEpidQAmanager;
 112 }
 113
 114 //________________________________________________________________________
 115 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
 116   :TNamed               ()
 117   ,fIsAOD               (kFALSE)
 118   ,fMCEvent             (NULL)
 119   ,fAODArrayMCInfo      (NULL)
 120   ,fLevelBack(-1)
 121   ,fHFEBackgroundCuts   (NULL)
 122   ,fPIDBackground       (0x0)
 123   ,fPIDBackgroundQA     (0)
 124   ,fkPIDRespons         (NULL)
 125   ,fPtBinning()
 126   ,fEtaBinning()
 127   ,fAlgorithmMA         (kTRUE)
 128   ,fChi2OverNDFCut      (3.0)
 129   ,fMaxDCA              (3.0)
 130 //  ,fMaxOpeningTheta   (0.02)
 131 //  ,fMaxOpeningPhi     (0.1)
 132   ,fMaxOpening3D        (TMath::TwoPi())
 133   ,fMaxInvMass          (1000)
 134   ,fSetMassConstraint   (kFALSE)
 135   ,fSelectCategory1tracks(kTRUE)
 136   ,fSelectCategory2tracks(kFALSE)
 137   ,fITSmeanShift(0.)
 138   ,fITSnSigmaHigh(3.)
 139   ,fITSnSigmaLow(-3.)
 140   ,fminPt(0.1)
 141   ,fArraytrack          (NULL)
 142   ,fCounterPoolBackground       (0)
 143   ,fnumberfound                 (0)
 144   ,fListOutput          (NULL)
 145   ,fAssElectron         (NULL)
 146   ,fIncElectron         (NULL)
 147   ,fUSign               (NULL)
 148   ,fLSign               (NULL)
 149   ,fUSmatches(NULL)
 150   ,fLSmatches(NULL)
 151   ,fHnsigmaITS(NULL)
 152   ,fWeightsSource(NULL)
 153   ,fIncElectronRadius(NULL)
 154   ,fRecElectronRadius(NULL)
 155 //  ,fUSignAngle        (NULL)
 156 //  ,fLSignAngle        (NULL)
 157 {
 158   //
 159   // Constructor
 160   //
 161   fPIDBackground   = new AliHFEpid("hfePidBackground");
 162   fPIDBackgroundQA = new AliHFEpidQAmanager;
 163 }
 164
 165 //________________________________________________________________________
 166 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const AliHFENonPhotonicElectron &ref)
 167   :TNamed(ref)
 168   ,fIsAOD               (ref.fIsAOD)
 169   ,fMCEvent             (NULL)
 170   ,fAODArrayMCInfo      (NULL)
 171   ,fLevelBack           (ref.fLevelBack)
 172   ,fHFEBackgroundCuts   (ref.fHFEBackgroundCuts)
 173   ,fPIDBackground       (ref.fPIDBackground)
 174   ,fPIDBackgroundQA     (ref.fPIDBackgroundQA)
 175   ,fkPIDRespons         (ref.fkPIDRespons)
 176   ,fPtBinning(ref.fPtBinning)
 177   ,fEtaBinning(ref.fEtaBinning)
 178   ,fAlgorithmMA         (ref.fAlgorithmMA)
 179   ,fChi2OverNDFCut      (ref.fChi2OverNDFCut)
 180   ,fMaxDCA              (ref.fMaxDCA)
 181 //  ,fMaxOpeningTheta   (ref.fMaxOpeningTheta)
 182 //  ,fMaxOpeningPhi     (ref.fMaxOpeningPhi)
 183   ,fMaxOpening3D        (ref.fMaxOpening3D)
 184   ,fMaxInvMass          (ref.fMaxInvMass)
 185   ,fSetMassConstraint   (ref.fSetMassConstraint)
 186   ,fSelectCategory1tracks(ref.fSelectCategory1tracks)
 187   ,fSelectCategory2tracks(ref.fSelectCategory2tracks)
 188   ,fITSmeanShift(ref.fITSmeanShift)
 189   ,fITSnSigmaHigh(ref.fITSnSigmaHigh)
 190   ,fITSnSigmaLow(ref.fITSnSigmaLow)
 191   ,fminPt(ref.fminPt)
 192   ,fArraytrack          (NULL)
 193   ,fCounterPoolBackground       (0)
 194   ,fnumberfound                 (0)
 195   ,fListOutput          (ref.fListOutput)
 196   ,fAssElectron         (ref.fAssElectron)
 197   ,fIncElectron         (ref.fIncElectron)
 198   ,fUSign               (ref.fUSign)
 199   ,fLSign               (ref.fLSign)
 200   ,fUSmatches(ref.fUSmatches)
 201   ,fLSmatches(ref.fLSmatches)
 202   ,fHnsigmaITS(ref.fHnsigmaITS)
 203   ,fWeightsSource(ref.fWeightsSource)
 204   ,fIncElectronRadius(ref.fIncElectronRadius)
 205   ,fRecElectronRadius(ref.fRecElectronRadius)
 206 //  ,fUSignAngle        (ref.fUSignAngle)
 207 //  ,fLSignAngle        (ref.fLSignAngle)
 208 {
 209   //
 210   // Copy Constructor
 211   //
 212   ref.Copy(*this);
 213 }
 214
 215 //____________________________________________________________
 216 AliHFENonPhotonicElectron &AliHFENonPhotonicElectron::operator=(const AliHFENonPhotonicElectron &ref){
 217   //
 218   // Assignment operator
 219   //
 220   if(this == &ref) ref.Copy(*this);
 221   return *this;
 222 }
 223
 224 //_________________________________________
 225 AliHFENonPhotonicElectron::~AliHFENonPhotonicElectron()
 226 {
 227   //
 228   // Destructor
 229   //
 230   if(fArraytrack)               delete fArraytrack;
 231   //if(fHFEBackgroundCuts)      delete fHFEBackgroundCuts;
 232   if(fPIDBackground)            delete fPIDBackground;
 233   if(fPIDBackgroundQA)          delete fPIDBackgroundQA;
 234 }
 235
 236 //_____________________________________________________________________________________________
 237 void AliHFENonPhotonicElectron::Init()
 238 {
 239   //
 240   // Init
 241   //
 242
 243   //printf("Analysis Mode for AliHFENonPhotonicElectron: %s Analysis\n", fIsAOD ? "AOD" : "ESD");
 244
 245   if(!fListOutput) fListOutput = new TList;
 246   fListOutput->SetName("HFENonPhotonicElectron");
 247   fListOutput->SetOwner();
 248
 249   if(!fHFEBackgroundCuts) fHFEBackgroundCuts = new AliHFEcuts();
 250   if(fIsAOD) fHFEBackgroundCuts->SetAOD();
 251   fHFEBackgroundCuts->Initialize();
 252   if(fHFEBackgroundCuts->IsQAOn()) {
 253     fListOutput->Add(fHFEBackgroundCuts->GetQAhistograms());
 254   }
 255
 256   // Initialize PID
 257   if(!fPIDBackground) fPIDBackground = new AliHFEpid("default pid");
 258   if(fMCEvent || fAODArrayMCInfo) fPIDBackground->SetHasMCData(kTRUE); // does nothing since the fMCEvent are set afterwards at the moment
 259   if(!fPIDBackground->GetNumberOfPIDdetectors())
 260   {
 261     //fPIDBackground->AddDetector("TOF", 0);
 262     fPIDBackground->AddDetector("TPC", 0);
 263   }
 264   AliInfo("PID Background QA switched on");
 265   fPIDBackgroundQA->Initialize(fPIDBackground);
 266   fListOutput->Add(fPIDBackgroundQA->MakeList("HFENP_PID_Background"));
 267   fPIDBackground->SortDetectors();
 268
 269   const Int_t kBinsPtDefault = 35;
 270   Double_t binLimPtDefault[kBinsPtDefault+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
 271   const Int_t kBinsEtaInclusiveDefault = 8;
 272   Double_t binLimEtaInclusiveDefault[kBinsEtaInclusiveDefault+1] = {-0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8};
 273   const Int_t kBinsEtaAssociated = 15;
 274   Double_t binLimEtaAssociat[kBinsEtaAssociated+1] = {-1.5,-1.3,-1.1,-0.9,-0.7,-0.5,-0.3,-0.1,0.1,0.3,0.5,0.7,0.9,1.1,1.3,1.5};
 275
 276   if(!fPtBinning.GetSize()) fPtBinning.Set(kBinsPtDefault+1, binLimPtDefault);
 277   if(!fEtaBinning.GetSize()) fEtaBinning.Set(kBinsEtaInclusiveDefault+1, binLimEtaInclusiveDefault);
 278
 279   //Int_t nBinsP = 400;
 280   //Double_t minP = 0.0;
 281   //Double_t maxP = 20.0;
 282   //Double_t binLimP[nBinsP+1];
 283   //for(Int_t i=0; i<=nBinsP; i++) binLimP[i]=(Double_t)minP + (maxP-minP)/nBinsP*(Double_t)i ;
 284
 285   Int_t nBinsC = 11;
 286   Double_t minC = 0.0;
 287   Double_t maxC = 11.0;
 288   Double_t binLimC[nBinsC+1];
 289   for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;
 290
 291   Int_t nBinsSource = 10;
 292   Double_t minSource = 0.;
 293   Double_t maxSource = 10.;
 294   Double_t binLimSource[nBinsSource+1];
 295   for(Int_t i=0; i<=nBinsSource; i++) binLimSource[i]=(Double_t)minSource + (maxSource-minSource)/nBinsSource*(Double_t)i ;
 296
 297   Int_t nBinsInvMass = 30;
 298   Double_t minInvMass = 0.;
 299   Double_t maxInvMass = 0.3;
 300   Double_t binLimInvMass[nBinsInvMass+1];
 301   for(Int_t i=0; i<=nBinsInvMass; i++) binLimInvMass[i]=(Double_t)minInvMass + (maxInvMass-minInvMass)/nBinsInvMass*(Double_t)i ;
 302
 303   Int_t nBinsPhi = 8;
 304   Double_t minPhi = 0.0;
 305   Double_t maxPhi = TMath::Pi();
 306   Double_t binLimPhi[nBinsPhi+1];
 307   for(Int_t i=0; i<=nBinsPhi; i++)
 308   {
 309     binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;
 310     AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
 311   }
 312
 313   Int_t nBinsAngle = 72;
 314   Double_t minAngle = 0.0;
 315   Double_t maxAngle = 0.4;
 316   Double_t binLimAngle[nBinsAngle+1];
 317   for(Int_t i=0; i<=nBinsAngle; i++)
 318   {
 319     binLimAngle[i]=(Double_t)minAngle + (maxAngle-minAngle)/nBinsAngle*(Double_t)i ;
 320     AliDebug(2,Form("bin phi is %f for %d",binLimAngle[i],i));
 321   }
 322
 323   // Constrain histograms
 324   const Int_t nDimSingle=4;
 325   const Int_t nDimPair=9;
 326   Int_t nBinPair[nDimPair] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1,nBinsInvMass,nBinsSource,nBinsAngle,fPtBinning.GetSize()-1,fEtaBinning.GetSize()-1,kBinsEtaAssociated};
 327
 328   // Associated Electron
 329   Int_t nBinAssElectron[nDimSingle] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource,kBinsEtaAssociated};
 330   fAssElectron = new THnSparseF("fAssElectron","fAssElectron",nDimSingle,nBinAssElectron);
 331   fAssElectron->SetBinEdges(0,binLimC);
 332   fAssElectron->SetBinEdges(1,fPtBinning.GetArray());
 333   fAssElectron->SetBinEdges(2,binLimSource);
 334   fAssElectron->SetBinEdges(3,binLimEtaAssociat);
 335   fAssElectron->Sumw2();
 336   AliDebug(2,"AliHFENonPhotonicElectron: fAssElectron");
 337
 338   // Inclusive Electron
 339   Int_t nBinIncElectron[nDimSingle] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource,fEtaBinning.GetSize()-1};
 340   fIncElectron = new THnSparseF("fIncElectron","fIncElectron",nDimSingle,nBinIncElectron);
 341   fIncElectron->SetBinEdges(0,binLimC);
 342   fIncElectron->SetBinEdges(1,fPtBinning.GetArray());
 343   fIncElectron->SetBinEdges(2,binLimSource);
 344   fIncElectron->SetBinEdges(3,fEtaBinning.GetArray());
 345   fIncElectron->Sumw2();
 346   AliDebug(2,"AliHFENonPhotonicElectron: fIncElectron");
 347
 348   // ee invariant mass Unlike Sign
 349   fUSign = new THnSparseF("fUSign","fUSign",nDimPair,nBinPair);
 350   fUSign->SetBinEdges(0,binLimPhi);
 351   fUSign->SetBinEdges(1,binLimC);
 352   fUSign->SetBinEdges(2,fPtBinning.GetArray());
 353   fUSign->SetBinEdges(3,binLimInvMass);
 354   fUSign->SetBinEdges(4,binLimSource);
 355   fUSign->SetBinEdges(5,binLimAngle);
 356   fUSign->SetBinEdges(6,fPtBinning.GetArray());
 357   fUSign->SetBinEdges(7,fEtaBinning.GetArray());
 358   fUSign->SetBinEdges(8,binLimEtaAssociat);
 359   fUSign->Sumw2();
 360   AliDebug(2,"AliHFENonPhotonicElectron: fUSign");
 361
 362   // ee invariant mass Like Sign
 363   fLSign = new THnSparseF("fLSign","fLSign",nDimPair,nBinPair);
 364   fLSign->SetBinEdges(0,binLimPhi);
 365   fLSign->SetBinEdges(1,binLimC);
 366   fLSign->SetBinEdges(2,fPtBinning.GetArray());
 367   fLSign->SetBinEdges(3,binLimInvMass);
 368   fLSign->SetBinEdges(4,binLimSource);
 369   fLSign->SetBinEdges(5,binLimAngle);
 370   fLSign->SetBinEdges(6,fPtBinning.GetArray());
 371   fLSign->SetBinEdges(7,fEtaBinning.GetArray());
 372   fLSign->SetBinEdges(8,binLimEtaAssociat);
 373   fLSign->Sumw2();
 374   AliDebug(2,"AliHFENonPhotonicElectron: fLSign");
 375
 376   // Histograms counting the number of like sign / unlike sign matches per inclusive track
 377   const Int_t nBinsMatches = 50;
 378   Double_t binLimMatches[nBinsMatches+1];
 379   for(int ib = 0; ib <= nBinsMatches; ib++) binLimMatches[ib] = ib;
 380   const Int_t nDimMatches = 3;  // centrality, pt_inc, number of matches
 381   const Int_t nBinsMatchHist[nDimMatches] = {nBinsC, fPtBinning.GetSize()-1, nBinsMatches};
 382   fUSmatches = new THnSparseF("fUSmatches", "fUSmatches", nDimMatches, nBinsMatchHist);
 383   fUSmatches->SetBinEdges(0,binLimC);
 384   fUSmatches->SetBinEdges(1,fPtBinning.GetArray());
 385   fUSmatches->SetBinEdges(2,binLimMatches);
 386
 387   fLSmatches = new THnSparseF("fLSmatches", "fLSmatches", nDimMatches, nBinsMatchHist);
 388   fLSmatches->SetBinEdges(0,binLimC);
 389   fLSmatches->SetBinEdges(1,fPtBinning.GetArray());
 390   fLSmatches->SetBinEdges(2,binLimMatches);
 391
 392   // Histograms for radius studies
 393   Int_t nBinsradius = 50;
 394   Double_t minradius = 0.0;
 395   Double_t maxradius = 100.0;
 396   Double_t binLimradius[nBinsradius+1];
 397   for(Int_t i=0; i<=nBinsradius; i++) binLimradius[i]=(Double_t)minradius + (maxradius-minradius)/nBinsradius*(Double_t)i ;
 398   const Int_t nDimIncElectronRadius = 4;  // centrality, pt_inc, radius
 399   const Int_t nBinsIncElectronRadius[nDimIncElectronRadius] = {nBinsC, fPtBinning.GetSize()-1, nBinsradius, nBinsSource};
 400   fIncElectronRadius = new THnSparseF("fIncElectronRadius", "fIncElectronRadius", nDimIncElectronRadius, nBinsIncElectronRadius);
 401   fIncElectronRadius->SetBinEdges(0,binLimC);
 402   fIncElectronRadius->SetBinEdges(1,fPtBinning.GetArray());
 403   fIncElectronRadius->SetBinEdges(2,binLimradius);
 404   fIncElectronRadius->SetBinEdges(3,binLimSource);
 405
 406   fRecElectronRadius = new THnSparseF("fRecElectronRadius", "fRecElectronRadius", nDimIncElectronRadius, nBinsIncElectronRadius);
 407   fRecElectronRadius->SetBinEdges(0,binLimC);
 408   fRecElectronRadius->SetBinEdges(1,fPtBinning.GetArray());
 409   fRecElectronRadius->SetBinEdges(2,binLimradius);
 410   fRecElectronRadius->SetBinEdges(2,binLimSource);
 411
 412 /*
 413   // ee angle Unlike Sign
 414   const Int_t nDimUSignAngle=3;
 415   Int_t nBinUSignAngle[nDimUSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
 416   fUSignAngle = new THnSparseF("fUSignAngle","fUSignAngle",nDimUSignAngle,nBinUSignAngle);
 417   fUSignAngle->SetBinEdges(0,binLimAngle);
 418   fUSignAngle->SetBinEdges(1,binLimC);
 419   fUSignAngle->SetBinEdges(2,binLimSource);
 420   fUSignAngle->Sumw2();
 421   AliDebug(2,"AliHFENonPhotonicElectron: fUSignAngle");
 422
 423   // ee angle Like Sign
 424   const Int_t nDimLSignAngle=3;
 425   Int_t nBinLSignAngle[nDimLSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
 426   fLSignAngle = new THnSparseF("fLSignAngle","fLSignAngle",nDimLSignAngle,nBinLSignAngle);
 427   fLSignAngle->SetBinEdges(0,binLimAngle);
 428   fLSignAngle->SetBinEdges(1,binLimC);
 429   fLSignAngle->SetBinEdges(2,binLimSource);
 430   fLSignAngle->Sumw2();
 431   AliDebug(2,"AliHFENonPhotonicElectron: fLSignAngle");
 432 */
 433
 434   // control histogram for ITS PID
 435   fHnsigmaITS = new TH2F("fHnsigmaITS", "Number of sigmas in the ITS", 30, 0., 0.3, 1200, -10., 10.);
 436
 437   // control histogram for weights sources
 438   fWeightsSource = new TH2F("fWeightsSource", "Source code for weights", 11, -1.5, 9.5, 29, -1.5, 27.5);
 439
 440   fListOutput->Add(fAssElectron);
 441   fListOutput->Add(fIncElectron);
 442   fListOutput->Add(fUSign);
 443   fListOutput->Add(fLSign);
 444   fListOutput->Add(fUSmatches);
 445   fListOutput->Add(fLSmatches);
 446   fListOutput->Add(fHnsigmaITS);
 447   fListOutput->Add(fWeightsSource);
 448   fListOutput->Add(fIncElectronRadius);
 449   fListOutput->Add(fRecElectronRadius);
 450 //  fListOutput->Add(fUSignAngle);
 451 //  fListOutput->Add(fLSignAngle);
 452
 453 }
 454
 455 //_____________________________________________________________________________________________
 456 void AliHFENonPhotonicElectron::SetWithWeights(Int_t levelBack)
 457 {
 458   //
 459   // Init the HFE level
 460   //
 461   if(levelBack >= 0) fLevelBack = levelBack;
 462
 463 }
 464
 465 //_____________________________________________________________________________________________
 466 void AliHFENonPhotonicElectron::SetMCEvent(AliMCEvent *mcEvent)
 467 {
 468   //
 469   // Pass the mcEvent
 470   //
 471
 472   fMCEvent = mcEvent;
 473
 474 }
 475
 476 //_____________________________________________________________________________________________
 477 void AliHFENonPhotonicElectron::SetAODArrayMCInfo(TClonesArray *aodArrayMCInfo)
 478 {
 479   //
 480   // Pass the mcEvent info
 481   //
 482
 483   fAODArrayMCInfo = aodArrayMCInfo;
 484
 485 }
 486
 487 //_____________________________________________________________________________________________
 488 void AliHFENonPhotonicElectron::InitRun(const AliVEvent *inputEvent,const AliPIDResponse *pidResponse)
 489 {
 490   //
 491   // Init run
 492   //
 493
 494   if(!pidResponse)
 495   {
 496     AliDebug(1, "Using default PID Response");
 497     Bool_t hasmc = kFALSE;
 498     if(fMCEvent || fAODArrayMCInfo) hasmc=kTRUE;
 499     pidResponse = AliHFEtools::GetDefaultPID(hasmc, inputEvent->IsA() == AliESDEvent::Class());
 500   }
 501
 502   if(!fPIDBackground) return;
 503   fPIDBackground->SetPIDResponse(pidResponse);
 504
 505   if(!fPIDBackground->IsInitialized())
 506   {
 507     // Initialize PID with the given run number
 508     fPIDBackground->InitializePID(inputEvent->GetRunNumber());
 509   }
 510
 511 }
 512
 513 //_____________________________________________________________________________________________
 514 Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent, Int_t binct)
 515 {
 516   //
 517   // Fill the pool of associated tracks
 518   // Return the number of associated tracks
 519   //
 520
 521   fnumberfound = 0;
 522
 523   fHFEBackgroundCuts->SetRecEvent(inputEvent);
 524   Int_t nbtracks = inputEvent->GetNumberOfTracks();
 525
 526   if( fArraytrack ){
 527     fArraytrack->~TArrayI();
 528     new(fArraytrack) TArrayI(nbtracks);
 529   } else {
 530     fArraytrack = new TArrayI(nbtracks);
 531   }
 532
 533   fCounterPoolBackground = 0;
 534
 535   Bool_t isSelected(kFALSE);
 536   Bool_t isAOD = (dynamic_cast<AliAODEvent *>(inputEvent) != NULL);
 537   AliDebug(2, Form("isAOD: %s", isAOD ? "yes" : "no"));
 538   for(Int_t k = 0; k < nbtracks; k++) {
 539     AliVTrack *track = (AliVTrack *) inputEvent->GetTrack(k);
 540     if(!track) continue;
 541
 542     //
 543     isSelected = kFALSE;
 544     if(fSelectCategory1tracks && FilterCategory1Track(track, isAOD, binct)) isSelected = kTRUE;
 545     else if(fSelectCategory2tracks && FilterCategory2Track(track, isAOD)) isSelected = kTRUE;
 546
 547     if(isSelected){
 548             AliDebug(2,Form("fCounterPoolBackground %d, track %d",fCounterPoolBackground,k));
 549       fArraytrack->AddAt(k,fCounterPoolBackground);
 550       fCounterPoolBackground++;
 551     }
 552   } // loop tracks
 553
 554   //printf(Form("Associated Pool: Tracks %d, fCounterPoolBackground %d \n", nbtracks, fCounterPoolBackground));
 555
 556   return fCounterPoolBackground;
 557
 558 }
 559
 560 //_____________________________________________________________________________________________
 561 Int_t AliHFENonPhotonicElectron::CountPoolAssociated(AliVEvent *inputEvent, Int_t binct)
 562 {
 563   //
 564   // Count the pool of assiocated tracks
 565   //
 566
 567
 568   if(fnumberfound > 0) //!count only events with an inclusive electron
 569   {
 570     Double_t valueAssElectron[4] = { (Double_t)binct, -1, -1};          //Centrality    Pt      Source
 571     Int_t iTrack2 = 0;
 572     Int_t indexmother2 = -1;
 573     AliVTrack *track2 = 0x0;
 574
 575     for(Int_t ii = 0; ii < fCounterPoolBackground; ii++){
 576       iTrack2 = fArraytrack->At(ii);
 577       AliDebug(2,Form("track %d",iTrack2));
 578       track2 = (AliVTrack *)inputEvent->GetTrack(iTrack2);
 579
 580       if(!track2){
 581               //printf("ERROR: Could not receive track %d", iTrack2);
 582               continue;
 583       }
 584
 585       // if MC look
 586       if(fMCEvent || fAODArrayMCInfo) valueAssElectron[2] = FindMother(TMath::Abs(track2->GetLabel()), indexmother2) ;
 587
 588       fkPIDRespons = fPIDBackground->GetPIDResponse();
 589
 590       valueAssElectron[1] = track2->Pt() ;
 591       valueAssElectron[3] = track2->Eta() ;
 592
 593       fAssElectron->Fill( valueAssElectron) ;
 594     }
 595   //printf(Form("Associated Pool: fCounterPoolBackground %d \n", fCounterPoolBackground));
 596   }
 597   return fnumberfound;
 598 }
 599
 600 //_____________________________________________________________________________________________
 601 Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, Double_t weight, Int_t binct, Double_t deltaphi, Int_t source, Int_t indexmother,Int_t mcQAsource)
 602 {
 603   //
 604   // Look At Non HFE
 605   //
 606
 607   /***********************************************************************************
 608    *                                                                                    *
 609    *    iTrack1:        index of the tagged electrons in AliVEvent                      *
 610    *    track1:         tagged electron                                                 *
 611    *    vEvent:         event                                                           *
 612    *    weight:         weight in pt if not realistic                                   *
 613    *    binct:          centrality bin                                                  *
 614    *    deltaphi:       phi-phi event plane for v2                                      *
 615    *    source:         MC sources                                                      *
 616    *    indexmother:    MC index mother                                                 *
 617    *                                                                                    *
 618    *                                                                                    *
 619    *    return -1  if  nothing                                                          *
 620    *    return  2  if  opposite         charge          within the mass range           *
 621    *    return  4  if      like         charge          within the mass range           *
 622    *    return  6  if  opposite & like charge           within the mass range           *
 623    *                                                                                    *
 624    ***********************************************************************************/
 625
 626   //printf("weight %f and source %d\n",weight,source);
 627
 628
 629   AliAODEvent *aodeventu = dynamic_cast<AliAODEvent*>(vEvent);
 630   Int_t taggedphotonic = -1;
 631
 632   AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
 633   if(!fArraytrack) return taggedphotonic;
 634   AliDebug(2,Form("process track %d",iTrack1));
 635   AliDebug(1,Form("Inclusive source is %d\n", source));
 636
 637   fkPIDRespons = fPIDBackground->GetPIDResponse();
 638
 639   //Set Fill-Arrays for THnSparse
 640   Double_t valueIncElectron[4]  = { (Double_t)binct, track1->Pt(), (Double_t)source, track1->Eta()};    //Centrality    Pt      Source  P
 641   Double_t valueSign[9]         = { deltaphi, (Double_t)binct, track1->Pt(), -1, (Double_t)source, -1, -1, track1->Eta(), -1};                  //DeltaPhi      Centrality      Pt      InvariantMass   Source  Angle   Pt
 642   //Double_t valueAngle[3]      = { -1, binct, source};
 643   Double_t valueradius[4]       = { (Double_t)binct, track1->Pt(), 0.,(Double_t)source};
 644
 645
 646   Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()));
 647   Double_t radius = Radius(TMath::Abs(track1->GetLabel()));
 648   AliKFParticle::SetField(vEvent->GetMagneticField());
 649   AliKFVertex primV(*(vEvent->GetPrimaryVertex()));
 650   valueradius[2] = radius;
 651
 652   AliVTrack *track2(NULL);
 653   Int_t iTrack2 = 0;
 654   Int_t indexmother2 = -1;
 655   Int_t pdg2 = -100;
 656   Int_t source2 = -1;
 657   Float_t fCharge2 = 0;
 658
 659   // count number of matches with opposite/same sign track in the given mass range
 660   Int_t countsMatchLikesign(0),
 661         countsMatchUnlikesign(0);
 662
 663   // Values to fill
 664   Double_t angle(-1.);
 665   Double_t invmass(-1);
 666
 667   Float_t fCharge1 = track1->Charge();                                                  //Charge from track1
 668
 669   Bool_t kUSignPhotonic = kFALSE;
 670   Bool_t kLSignPhotonic = kFALSE;
 671
 672   //! FILL Inclusive Electron
 673   fWeightsSource->Fill(source,mcQAsource);
 674   fIncElectron->Fill(valueIncElectron,weight);
 675   fnumberfound++;
 676   fIncElectronRadius->Fill(valueradius,weight);
 677
 678   //printf(Form("Inclusive Pool: TrackNr. %d, fnumberfound %d \n", iTrack1, fnumberfound));
 679
 680   for(Int_t idex = 0; idex < fCounterPoolBackground; idex++){
 681     iTrack2 = fArraytrack->At(idex);
 682     AliDebug(2,Form("track %d",iTrack2));
 683     track2 = (AliVTrack *)vEvent->GetTrack(iTrack2);
 684
 685     if(!track2){
 686       //printf("ERROR: Could not receive track %d", iTrack2);
 687       continue;
 688     }
 689
 690     fCharge2 = track2->Charge();                //Charge from track2
 691
 692     // Reset the MC info
 693     //valueAngle[2] = source;
 694     valueSign[4] = source;
 695     valueSign[6] = track2->Pt();
 696     valueSign[8] = track2->Eta();
 697
 698     // track cuts and PID already done
 699
 700     // Checking if it is the same Track!
 701     if(iTrack2==iTrack1) continue;
 702     AliDebug(2,"Different");
 703
 704     // if MC look
 705     if(fMCEvent || fAODArrayMCInfo){
 706       AliDebug(2, "Checking for source");
 707       source2    = FindMother(TMath::Abs(track2->GetLabel()), indexmother2);
 708       AliDebug(2, Form("source is %d", source2));
 709       pdg2       = CheckPdg(TMath::Abs(track2->GetLabel()));
 710
 711       if(source == kElectronfromconversion){
 712         AliDebug(2, Form("Electron from conversion (source %d), paired with source %d", source, source2));
 713         AliDebug(2, Form("Index of the mothers: incl %d, associated %d", indexmother, indexmother2));
 714         AliDebug(2, Form("PDGs: incl %d, associated %d", pdg1, pdg2));
 715       }
 716
 717       if(source2 >=0 ){
 718               if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)){
 719           AliDebug(1, "Real pair");
 720           switch(source){
 721             case kElectronfromconversion:
 722                  valueSign[4] = kElectronfromconversionboth;
 723                  valueradius[3] = kElectronfromconversionboth;
 724                  break;
 725             case kElectronfrompi0:
 726                  valueSign[4] = kElectronfrompi0both;
 727                  valueradius[3] = kElectronfrompi0both;
 728                  break;
 729             case kElectronfrometa:
 730                  valueSign[4] = kElectronfrometaboth;
 731                  valueradius[3] = kElectronfrometaboth;
 732                  break;
 733           };
 734         }
 735       }
 736     }
 737
 738     if(fAlgorithmMA){
 739       // Use TLorentzVector
 740       if(!MakePairDCA(track1, track2, vEvent, (aodeventu != NULL), invmass, angle)) continue;
 741     } else {
 742       // Use AliKF package
 743       if(!MakePairKF(track1, track2, primV, invmass, angle)) continue;
 744     }
 745
 746     valueSign[3] = invmass;
 747     valueSign[5] = angle;
 748
 749     //if((fCharge1*fCharge2)>0.0)       fLSignAngle->Fill(&valueAngle[0],weight);
 750     //else                              fUSignAngle->Fill(&valueAngle[0],weight);
 751
 752     if(angle > fMaxOpening3D) continue;                          //! Cut on Opening Angle
 753     if(invmass > fMaxInvMass) continue;                         //! Cut on Invariant Mass
 754
 755     if((fCharge1*fCharge2)>0.0){
 756       if(invmass < 1.0) fLSign->Fill( valueSign, weight);
 757       // count like-sign background matched pairs per inclusive based on mass cut
 758       if(invmass < 0.14) countsMatchLikesign++;
 759       AliDebug(1, "Selected Like sign");
 760     } else {
 761       if(invmass < 1.0)fUSign->Fill( valueSign, weight);
 762       // count unlike-sign matched pairs per inclusive based on mass cut
 763       if(invmass < 0.14) {
 764         countsMatchUnlikesign++;
 765         fRecElectronRadius->Fill(valueradius,weight);
 766       }
 767       AliDebug(1, "Selected Unike sign");
 768     }
 769
 770     if((fCharge1*fCharge2)>0.0) kLSignPhotonic=kTRUE;
 771     else                                kUSignPhotonic=kTRUE;
 772   }
 773
 774   // Fill counted
 775   Double_t valCountsLS[3] = {(Double_t)binct, track1->Pt(), (Double_t)countsMatchLikesign},
 776     valCountsUS[3] = {(Double_t)binct, track1->Pt(), (Double_t)countsMatchUnlikesign};
 777   fUSmatches->Fill(valCountsUS);
 778   fLSmatches->Fill(valCountsLS);
 779
 780   if( kUSignPhotonic &&  kLSignPhotonic) taggedphotonic = 6;
 781   if(!kUSignPhotonic &&  kLSignPhotonic) taggedphotonic = 4;
 782   if( kUSignPhotonic && !kLSignPhotonic) taggedphotonic = 2;
 783
 784   return taggedphotonic;
 785 }
 786
 787 //_________________________________________________________________________
 788 Int_t AliHFENonPhotonicElectron::FindMother(Int_t tr, Int_t &indexmother) const {
 789   //
 790   // Find the mother if MC
 791   //
 792
 793   if(!fMCEvent && !fAODArrayMCInfo) return -1;
 794
 795   Int_t pdg = CheckPdg(tr);
 796   if(TMath::Abs(pdg)!= 11)
 797   {
 798     indexmother = -1;
 799     return kNoElectron;
 800   }
 801
 802   indexmother = IsMotherGamma(tr);
 803   if(indexmother > 0) return kElectronfromconversion;
 804   indexmother = IsMotherPi0(tr);
 805   if(indexmother > 0) return kElectronfrompi0;
 806   indexmother = IsMotherC(tr);
 807   if(indexmother > 0) return kElectronfromC;
 808   indexmother = IsMotherB(tr);
 809   if(indexmother > 0) return kElectronfromB;
 810   indexmother = IsMotherEta(tr);
 811   if(indexmother > 0) return kElectronfrometa;
 812
 813   return kElectronfromother;
 814 }
 815
 816 //________________________________________________________________________________________________
 817 Int_t AliHFENonPhotonicElectron::CheckPdg(Int_t tr) const {
 818
 819   //
 820   // Return the pdg of the particle
 821   //
 822
 823   Int_t pdgcode = -1;
 824   if(tr < 0) return pdgcode;
 825
 826   AliMCParticle *mctrackesd = NULL; AliAODMCParticle *mctrackaod = NULL;
 827   if(fMCEvent){
 828     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 829     if(mctrack){
 830       if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) pdgcode = mctrackesd->PdgCode();
 831       else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) pdgcode = mctrackaod->GetPdgCode();
 832     }
 833   } else if(fAODArrayMCInfo) {
 834     if(tr < fAODArrayMCInfo->GetEntriesFast()){
 835       mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 836       if(mctrackaod) return pdgcode = mctrackaod->GetPdgCode();
 837     }
 838   }
 839
 840   return pdgcode;
 841 }
 842
 843 //________________________________________________________________________________________________
 844 Double_t AliHFENonPhotonicElectron::Radius(Int_t tr) const {
 845
 846   //
 847   // Return the production vertex radius
 848   //
 849
 850   Double_t radius = 0.;
 851   if(tr < 0) return radius;
 852
 853   AliMCParticle *mctrackesd = NULL; AliAODMCParticle *mctrackaod = NULL;
 854   if(fMCEvent){
 855     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 856     if(mctrack){
 857       if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) radius = TMath::Sqrt(mctrackesd->Xv()*mctrackesd->Xv()+mctrackesd->Yv()*mctrackesd->Yv());
 858       else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) radius = TMath::Sqrt(mctrackaod->Xv()*mctrackaod->Xv()+mctrackaod->Yv()*mctrackaod->Yv());
 859     }
 860   } else if(fAODArrayMCInfo) {
 861     if(tr < fAODArrayMCInfo->GetEntriesFast()){
 862       mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 863       if(mctrackaod) return radius = TMath::Sqrt(mctrackaod->Xv()*mctrackaod->Xv()+mctrackaod->Yv()*mctrackaod->Yv());
 864     }
 865   }
 866
 867   return radius;
 868 }
 869
 870 //_______________________________________________________________________________________________
 871 Int_t AliHFENonPhotonicElectron::GetMotherPDG(Int_t tr, Int_t &motherIndex) const {
 872   //
 873   // Returns the mother PDG of the track (return value) and the index of the
 874   // mother track
 875   //
 876   if(tr < 0) return -1;
 877
 878   Int_t pdg(-1);
 879   AliMCParticle *mctrackesd(NULL); AliAODMCParticle *mctrackaod(NULL);
 880
 881   motherIndex = -1;
 882   if(fMCEvent) {
 883     AliDebug(2, "Using MC Event");
 884     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 885     if(mctrack){
 886       if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))){
 887         // Case ESD
 888         TParticle *particle = mctrackesd->Particle();
 889
 890         // Take mother
 891         if(particle){
 892           motherIndex   = particle->GetFirstMother();
 893           if(motherIndex >= 0){
 894             AliMCParticle *mothertrack = NULL;
 895             if((mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherIndex))))){
 896               TParticle * mother = mothertrack->Particle();
 897               pdg = mother->GetPdgCode();
 898             }
 899           }
 900         }
 901       } else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))){
 902         // Case AOD
 903         // Take mother
 904         motherIndex = mctrackaod->GetMother();
 905         if(motherIndex >= 0){
 906           AliAODMCParticle *mothertrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(motherIndex));
 907           if(mothertrack) pdg = mothertrack->GetPdgCode();
 908         }
 909       }
 910     }
 911   } else if(fAODArrayMCInfo) {
 912     AliDebug(2, "Using AOD list");
 913     if(tr < fAODArrayMCInfo->GetEntriesFast()){
 914       mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 915
 916       // Take mother
 917       if(mctrackaod){
 918         motherIndex = mctrackaod->GetMother();
 919         if(motherIndex >= 0 && motherIndex < fAODArrayMCInfo->GetEntriesFast()){
 920           AliAODMCParticle *mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(motherIndex)));
 921           if(mothertrack) pdg = mothertrack->GetPdgCode();
 922         }
 923       }
 924     }
 925   }
 926   return pdg;
 927 }
 928
 929 //_______________________________________________________________________________________________
 930 Int_t AliHFENonPhotonicElectron::IsMotherGamma(Int_t tr) const {
 931
 932   //
 933   // Return the lab of gamma mother or -1 if not gamma
 934   //
 935
 936   Int_t imother(-1), pdg(-1);
 937   pdg = GetMotherPDG(tr, imother);
 938
 939   // Check gamma
 940   if(imother >= 0){
 941     if(TMath::Abs(pdg) == 22){
 942       AliDebug(2, "Gamma Mother selected");
 943       return imother;
 944     }
 945     if(TMath::Abs(pdg) == 11){
 946       AliDebug(2, "Mother is electron - look further in hierarchy");
 947       return IsMotherGamma(imother);
 948     }
 949     AliDebug(2, "Nothing selected");
 950     return -1;
 951   }
 952   AliDebug(2, "Not mother");
 953   return -1;
 954 }
 955
 956 //________________________________________________________________________________________________
 957 Int_t AliHFENonPhotonicElectron::IsMotherPi0(Int_t tr) const {
 958
 959   //
 960   // Return the lab of pi0 mother or -1 if not pi0
 961   //
 962
 963   Int_t imother(-1), pdg(-1);
 964   pdg = GetMotherPDG(tr, imother);
 965
 966   // Check pi0
 967   if(imother >= 0){
 968     if(TMath::Abs(pdg) == 111){
 969       AliDebug(2, "Pi0 Mother selected");
 970       return imother;
 971     }
 972     if(TMath::Abs(pdg) == 11){
 973       AliDebug(2, "Mother is electron - look further in hierarchy");
 974       return IsMotherPi0(imother);
 975     }
 976     AliDebug(2, "Nothing selected");
 977     return -1;
 978   }
 979   AliDebug(2, "Not mother");
 980   return -1;
 981 }
 982 //________________________________________________________________________________________________
 983 Int_t AliHFENonPhotonicElectron::IsMotherC(Int_t tr) const {
 984
 985   //
 986   // Return the lab of signal mother or -1 if not from C
 987   //
 988
 989   Int_t imother(-1), pdg(-1);
 990   pdg = GetMotherPDG(tr, imother);
 991
 992     // Check C
 993   if(imother >= 0){
 994     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)){
 995       AliDebug(2, "Charm Mother selected");
 996       return imother;
 997     }
 998     if(TMath::Abs(pdg) == 11){
 999       AliDebug(2, "Mother is electron - look further in hierarchy");
1000       return IsMotherC(imother);
1001     }
1002     AliDebug(2, "Nothing selected");
1003     return -1;
1004   }
1005   AliDebug(2, "Not mother");
1006   return -1;
1007 }
1008
1009 //_______________________________________________________________________________________________
1010 Int_t AliHFENonPhotonicElectron::IsMotherB(Int_t tr) const {
1011
1012   //
1013   // Return the lab of signal mother or -1 if not B
1014   //
1015
1016   Int_t imother(-1), pdg(-1);
1017   pdg = GetMotherPDG(tr, imother);
1018
1019   // Check B
1020   if(imother >= 0){
1021     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)){
1022       AliDebug(2, "Bottom Mother selected");
1023       return imother;
1024     }
1025     if(TMath::Abs(pdg) == 11){
1026       return IsMotherB(imother);
1027       AliDebug(2, "Mother is electron - look further in hierarchy");
1028     }
1029     AliDebug(2, "Nothing selected");
1030     return -1;
1031   }
1032   AliDebug(2, "Not mother");
1033   return -1;
1034 }
1035
1036 //_______________________________________________________________________________________________
1037 Int_t AliHFENonPhotonicElectron::IsMotherEta(Int_t tr) const {
1038
1039   //
1040   // Return the lab of eta mother or -1 if not eta
1041   //
1042
1043   Int_t imother(-1), pdg(-1);
1044   pdg = GetMotherPDG(tr, imother);
1045
1046   // Check eta
1047   if(imother >= 0){
1048     if(TMath::Abs(pdg) == 221){
1049       AliDebug(2, "Eta mother selected");
1050       return imother;
1051     }
1052     if(TMath::Abs(pdg) == 11){
1053       AliDebug(2, "Mother is electron - look further in hierarchy");
1054       return IsMotherEta(imother);
1055     }
1056     AliDebug(2, "Nothing selected");
1057     return -1;
1058   }
1059   AliDebug(2, "Not mother");
1060   return -1;
1061 }
1062
1063 //_______________________________________________________________________________________________
1064 Bool_t AliHFENonPhotonicElectron::MakePairDCA(const AliVTrack *inclusive, const AliVTrack *associated, AliVEvent *vEvent, Bool_t isAOD, Double_t &invMass, Double_t &angle) const {
1065   //
1066   // Make Pairs of electrons using TLorentzVector
1067   //
1068   Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
1069   Double_t bfield = vEvent->GetMagneticField();
1070
1071   AliESDtrack *esdtrack1, *esdtrack2;
1072   if(isAOD){
1073     // call copy constructor for AODs
1074     esdtrack1 = new AliESDtrack(inclusive);
1075     esdtrack2 = new AliESDtrack(associated);
1076   } else {
1077     // call copy constructor for ESDs
1078     esdtrack1 = new AliESDtrack(*(static_cast<const AliESDtrack *>(inclusive)));
1079     esdtrack2 = new AliESDtrack(*(static_cast<const AliESDtrack *>(associated)));
1080   }
1081   if((!esdtrack1) || (!esdtrack2)){
1082     delete esdtrack1;
1083     delete esdtrack2;
1084     return kFALSE;
1085   }
1086
1087   Double_t xt1 = 0; //radial position track 1 at the DCA point
1088   Double_t xt2 = 0; //radial position track 2 at the DCA point
1089   Double_t dca = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);           //DCA track1-track2
1090   if(dca > fMaxDCA){
1091     // Apply DCA cut already in the function
1092     delete esdtrack1;
1093     delete esdtrack2;
1094     return kFALSE;
1095   }
1096
1097   //Momenta of the track extrapolated to DCA track-track
1098   Double_t p1[3] = {0,0,0};
1099   Double_t p2[3] = {0,0,0};
1100   Bool_t kHasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);             //Track1
1101   Bool_t kHasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);             //Track2
1102   if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation");
1103
1104   TLorentzVector electron1, electron2, mother;
1105   electron1.SetXYZM(p1[0], p1[1], p1[2], eMass);
1106   electron2.SetXYZM(p2[0], p2[1], p2[2], eMass);
1107
1108   mother      = electron1 + electron2;
1109   invMass  = mother.M();
1110   angle    = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
1111
1112   delete esdtrack1;
1113   delete esdtrack2;
1114   return kTRUE;
1115 }
1116
1117 //_______________________________________________________________________________________________
1118 Bool_t AliHFENonPhotonicElectron::MakePairKF(const AliVTrack *inclusive, const AliVTrack *associated, AliKFVertex &primV, Double_t &invMass, Double_t &angle) const {
1119   //
1120   // Make pairs of electrons using the AliKF package
1121   //
1122
1123   //printf("AOD HFE non photonic\n");
1124
1125   Int_t fPDGtrack1 = 11;
1126   if(inclusive->Charge()>0) fPDGtrack1 = -11;
1127   Int_t fPDGtrack2 = 11;
1128   if(associated->Charge()>0) fPDGtrack2 = -11;
1129
1130   AliKFParticle ktrack1(*inclusive, fPDGtrack1);
1131   AliKFParticle ktrack2(*associated, fPDGtrack2);
1132   AliKFParticle recoGamma(ktrack1,ktrack2);
1133
1134   if(recoGamma.GetNDF()<1) return kFALSE;                               //! Cut on Reconstruction
1135
1136   Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
1137   if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) return kFALSE;
1138
1139   if(fSetMassConstraint){
1140           primV += recoGamma;
1141           primV -= ktrack1;
1142           primV -= ktrack2;
1143           recoGamma.SetProductionVertex(primV);
1144           recoGamma.SetMassConstraint(0,0.0001);
1145   }
1146
1147   Double_t width(0.);
1148   recoGamma.GetMass(invMass,width);
1149   angle = ktrack1.GetAngle(ktrack2);
1150   return kTRUE;
1151 }
1152
1153 //_______________________________________________________________________________________________
1154 Bool_t AliHFENonPhotonicElectron::FilterCategory1Track(const AliVTrack * const track, Bool_t isAOD, Int_t binct){
1155   //
1156   // Selection of good associated tracks for the pool
1157   // selection is done using strong cuts
1158   // Tracking in the TPC and the ITS is a minimal requirement
1159   //
1160   Bool_t survivedbackground = kTRUE;
1161
1162   if(!fHFEBackgroundCuts->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *) track)) survivedbackground = kFALSE;
1163   AliDebug(3, Form("First cut: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
1164   if(!fHFEBackgroundCuts->CheckParticleCuts(AliHFEcuts::kStepRecPrim       + AliHFEcuts::kNcutStepsMCTrack, (TObject *) track)) survivedbackground = kFALSE;
1165   AliDebug(3, Form("Second cut: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
1166
1167   if(survivedbackground){
1168     // PID track cuts
1169     AliHFEpidObject hfetrack2;
1170
1171     if(!isAOD)  hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1172     else                hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1173
1174     hfetrack2.SetRecTrack(track);
1175     if(binct>-1){
1176      hfetrack2.SetCentrality((Int_t)binct);
1177      AliDebug(3,Form("centrality %d and %d",binct,hfetrack2.GetCentrality()));
1178      hfetrack2.SetPbPb();
1179     }
1180
1181     if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundQA)){
1182       survivedbackground = kTRUE;
1183     } else survivedbackground = kFALSE;
1184   }
1185   AliDebug(3, Form("PID: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
1186   return survivedbackground;
1187 }
1188
1189 //_______________________________________________________________________________________________
1190 Bool_t AliHFENonPhotonicElectron::FilterCategory2Track(const AliVTrack * const track, Bool_t isAOD){
1191   //
1192   // Selection of category 2 tracks: These tracks are exclusively low pt tracks below
1193   // 300 MeV/c which have at least 2 hits in the 4 outer ITS layers and are identified as
1194   // electron candidates by the ITS
1195   //
1196   if(TMath::Abs(track->Pt()) > 0.3) return kFALSE;
1197   if(TMath::Abs(track->Pt()) < fminPt) return kFALSE;
1198   Int_t nclustersITS(0), nclustersOuter(0);
1199   if(isAOD){
1200     const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
1201     if(!(aodtrack->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA) || aodtrack->TestFilterBit(AliAODTrack::kTrkITSsa))) return kFALSE;
1202     if(!aodtrack->IsOn(AliAODTrack::kITSrefit)) return kFALSE;
1203     nclustersITS = aodtrack->GetITSNcls();
1204     for(int ily = 2; ily < 5; ily++)
1205       if(aodtrack->HasPointOnITSLayer(ily)) nclustersOuter++;
1206   } else {
1207     const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
1208     if(esdtrack->GetStatus() & AliESDtrack::kITSpureSA) return kFALSE;
1209     if(esdtrack->GetStatus() & AliESDtrack::kTPCin) return kFALSE;
1210     if(!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1211     nclustersITS = esdtrack->GetITSclusters(NULL);
1212     for(int ily = 2; ily < 5; ily++)
1213       if(esdtrack->HasPointOnITSLayer(ily)) nclustersOuter++;
1214   }
1215   if(nclustersITS < 4) return kFALSE;
1216   if(nclustersOuter < 3) return kFALSE;
1217
1218   // Do ITS PID
1219   Double_t nsigmaITS = fPIDBackground->GetPIDResponse()->NumberOfSigmasITS(track, AliPID::kElectron);
1220   fHnsigmaITS->Fill(track->Pt(), nsigmaITS);
1221   if((nsigmaITS - fITSmeanShift) > fITSnSigmaHigh ) return kFALSE;
1222   if((nsigmaITS - fITSmeanShift) < fITSnSigmaLow ) return kFALSE;
1223   // if global track, we apply also TPC PID
1224   return kTRUE;
1225 }