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