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