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  *      without fee, provided that the above copyright notice appears in all            *
  10  *      copies and that both the copyright notice and this permission notice            *
  11  *      appear in the supporting documentation. The authors make no claims              *
  12  *      about the suitability of this software for any purpose. It is                   *
  13  *      provided "as is" without express or implied warranty.                           *
  14  *                                                                                      *
  15  *************************************************************************************/
  16
  17 /*************************************************************************************
  18  *                                                                                      *
  19  *      Class for the Selection of Non-Heavy-Flavour-Electrons trought          *
  20  *      the invariant mass method. The selection can be done from two                   *
  21  *      different algorithms, which can be choosed calling the function         *
  22  *              "SetAlgorithm(TString Algorithm)".                                      *
  23  *                                                                                      *
  24  *              Authors: R.Bailhache, C.A.Schmidt                                       *
  25  *                                                                                      *
  26  *************************************************************************************/
  27
  28 #include "TVector2.h"
  29 #include "THnSparse.h"
  30 #include "TMath.h"
  31 #include "TLorentzVector.h"
  32 #include "TParticle.h"
  33 #include "TList.h"
  34 #include "TDatabasePDG.h"
  35
  36 #include "AliVEvent.h"
  37 #include "AliMCEvent.h"
  38 #include "AliESDEvent.h"
  39 #include "AliMCParticle.h"
  40 #include "AliAODMCParticle.h"
  41 #include "AliAODEvent.h"
  42 #include "AliAODVertex.h"
  43 #include "AliAODTrack.h"
  44 #include "AliVTrack.h"
  45 #include "AliESDtrack.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliPIDResponse.h"
  48
  49 #include "AliKFParticle.h"
  50 #include "AliKFVertex.h"
  51
  52 #include "AliHFEcuts.h"
  53 #include "AliHFEpid.h"
  54 #include "AliHFEpidQAmanager.h"
  55 #include "AliHFEtools.h"
  56
  57 #include "AliHFENonPhotonicElectron.h"
  58
  59 ClassImp(AliHFENonPhotonicElectron)
  60 //________________________________________________________________________
  61 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const char *name, const Char_t *title)
  62   :TNamed(name, title)
  63   ,fMCEvent(NULL)
  64   ,fAODArrayMCInfo(NULL)
  65   ,fHFEBackgroundCuts(0)
  66   ,fPIDBackground(0x0)
  67   ,fPIDBackgroundQA(0)
  68   ,fAlgorithmMA(kTRUE)
  69   ,fUseFilterAOD(kTRUE)
  70   ,fFilter(-1)
  71   ,fChi2OverNDFCut(3.0)
  72   ,fMaxDCA(3.0)
  73 //  ,fMaxOpeningTheta(0.02)
  74 //  ,fMaxOpeningPhi(0.1)
  75   ,fMaxOpening3D(TMath::Pi())
  76   ,fMaxInvMass(0.6)
  77   ,fSetMassConstraint(kFALSE)
  78   ,fArraytrack(NULL)
  79   ,fCounterPoolBackground(0)
  80   ,fListOutput(NULL)
  81   ,fMCSource(NULL)
  82   ,fUSign(NULL)
  83   ,fLSign(NULL)
  84 //  ,fUSignAngle(NULL)
  85 //  ,fLSignAngle(NULL)
  86 {
  87   //
  88   // Constructor
  89   //
  90   fPIDBackground = new AliHFEpid("hfePidBackground");
  91   fPIDBackgroundQA = new AliHFEpidQAmanager;
  92 }
  93
  94 //________________________________________________________________________
  95 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
  96   :TNamed()
  97   ,fMCEvent(NULL)
  98   ,fAODArrayMCInfo(NULL)
  99   ,fHFEBackgroundCuts(NULL)
 100   ,fPIDBackground(0x0)
 101   ,fPIDBackgroundQA(0)
 102   ,fAlgorithmMA(kTRUE)
 103   ,fUseFilterAOD(kTRUE)
 104   ,fFilter(-1)
 105   ,fChi2OverNDFCut(3.0)
 106   ,fMaxDCA(3.0)
 107 //  ,fMaxOpeningTheta(0.02)
 108 //  ,fMaxOpeningPhi(0.1)
 109   ,fMaxOpening3D(TMath::Pi())
 110   ,fMaxInvMass(0.6)
 111   ,fSetMassConstraint(kFALSE)
 112   ,fArraytrack(NULL)
 113   ,fCounterPoolBackground(0)
 114   ,fListOutput(NULL)
 115   ,fMCSource(NULL)
 116   ,fUSign(NULL)
 117   ,fLSign(NULL)
 118 //  ,fUSignAngle(NULL)
 119 //  ,fLSignAngle(NULL)
 120 {
 121   //
 122   // Constructor
 123   //
 124   fPIDBackground   = new AliHFEpid("hfePidBackground");
 125   fPIDBackgroundQA = new AliHFEpidQAmanager;
 126 }
 127
 128 //________________________________________________________________________
 129 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const AliHFENonPhotonicElectron &ref)
 130   :TNamed(ref)
 131   ,fMCEvent(NULL)
 132   ,fAODArrayMCInfo(NULL)
 133   ,fHFEBackgroundCuts(ref.fHFEBackgroundCuts)
 134   ,fPIDBackground(ref.fPIDBackground)
 135   ,fPIDBackgroundQA(ref.fPIDBackgroundQA)
 136   ,fAlgorithmMA(ref.fAlgorithmMA)
 137   ,fUseFilterAOD(ref.fUseFilterAOD)
 138   ,fFilter(ref.fFilter)
 139   ,fChi2OverNDFCut(ref.fChi2OverNDFCut)
 140   ,fMaxDCA(ref.fMaxDCA)
 141 //  ,fMaxOpeningTheta(ref.fMaxOpeningTheta)
 142 //  ,fMaxOpeningPhi(ref.fMaxOpeningPhi)
 143   ,fMaxOpening3D(ref.fMaxOpening3D)
 144   ,fMaxInvMass(ref.fMaxInvMass)
 145   ,fSetMassConstraint(ref.fSetMassConstraint)
 146   ,fArraytrack(NULL)
 147   ,fCounterPoolBackground(0)
 148   ,fListOutput(ref.fListOutput)
 149   ,fMCSource(ref.fMCSource)
 150   ,fUSign(ref.fUSign)
 151   ,fLSign(ref.fLSign)
 152 //  ,fUSignAngle(ref.fUSignAngle)
 153 //  ,fLSignAngle(ref.fLSignAngle)
 154 {
 155   //
 156   // Copy Constructor
 157   //
 158   ref.Copy(*this);
 159 }
 160
 161 //____________________________________________________________
 162 AliHFENonPhotonicElectron &AliHFENonPhotonicElectron::operator=(const AliHFENonPhotonicElectron &ref){
 163   //
 164   // Assignment operator
 165   //
 166   if(this == &ref) ref.Copy(*this);
 167   return *this;
 168 }
 169
 170 //_________________________________________
 171 AliHFENonPhotonicElectron::~AliHFENonPhotonicElectron()
 172 {
 173   //
 174   // Destructor
 175   //
 176   if(fArraytrack) delete fArraytrack;
 177   if(fHFEBackgroundCuts) delete fHFEBackgroundCuts;
 178   if(fPIDBackground) delete fPIDBackground;
 179 }
 180
 181 //_____________________________________________________________________________________________
 182 void AliHFENonPhotonicElectron::Init()
 183 {
 184   //
 185   // Init
 186   //
 187
 188   if(!fListOutput) fListOutput = new TList;
 189   fListOutput->SetName("HFENonPhotonicElectron");
 190   fListOutput->SetOwner();
 191
 192   if(!fHFEBackgroundCuts) fHFEBackgroundCuts = new AliESDtrackCuts();
 193
 194   // Initialize PID
 195   if(!fPIDBackground) fPIDBackground = new AliHFEpid("default pid");
 196   if(fMCEvent || fAODArrayMCInfo) fPIDBackground->SetHasMCData(kTRUE);
 197   if(!fPIDBackground->GetNumberOfPIDdetectors()) fPIDBackground->AddDetector("TPC", 0);
 198   AliInfo("PID Background QA switched on");
 199   fPIDBackgroundQA->Initialize(fPIDBackground);
 200   fListOutput->Add(fPIDBackgroundQA->MakeList("HFENP_PID_Background"));
 201   fPIDBackground->SortDetectors();
 202
 203   Int_t nBinsPt = 24;
 204   Double_t binLimPt[25] = {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, 3., 3.5, 4., 5., 6.};
 205
 206   Int_t nBinsC = 11;
 207   Double_t minC = 0.0;
 208   Double_t maxC = 11.0;
 209   Double_t binLimC[nBinsC+1];
 210   for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;
 211
 212   Int_t nBinsSource = 10;
 213   Double_t minSource = 0.;
 214   Double_t maxSource = 10.;
 215   Double_t binLimSource[nBinsSource+1];
 216   for(Int_t i=0; i<=nBinsSource; i++) binLimSource[i]=(Double_t)minSource + (maxSource-minSource)/nBinsSource*(Double_t)i ;
 217
 218   Int_t nBinsInvMass = 60;
 219   Double_t minInvMass = 0.;
 220   Double_t maxInvMass = 0.6;
 221   Double_t binLimInvMass[nBinsInvMass+1];
 222   for(Int_t i=0; i<=nBinsInvMass; i++) binLimInvMass[i]=(Double_t)minInvMass + (maxInvMass-minInvMass)/nBinsInvMass*(Double_t)i ;
 223
 224   Int_t nBinsPhi = 180;
 225   Double_t minPhi = 0.0;
 226   Double_t maxPhi = TMath::Pi();
 227   Double_t binLimPhi[nBinsPhi+1];
 228   for(Int_t i=0; i<=nBinsPhi; i++)
 229   {
 230     binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;
 231     AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
 232   }
 233
 234   Int_t nBinsAngle = 180;
 235   Double_t minAngle = 0.0;
 236   Double_t maxAngle = TMath::Pi();
 237   Double_t binLimAngle[nBinsAngle+1];
 238   for(Int_t i=0; i<=nBinsAngle; i++)
 239   {
 240     binLimAngle[i]=(Double_t)minAngle + (maxAngle-minAngle)/nBinsAngle*(Double_t)i ;
 241     AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
 242   }
 243
 244   const Int_t nDimMCSource=3;
 245   Int_t nBinMCSource[nDimMCSource] = {nBinsC,nBinsPt,nBinsSource};
 246   fMCSource = new THnSparseF("fMCSource","fMCSource",nDimMCSource,nBinMCSource);
 247   fMCSource->SetBinEdges(0,binLimC);
 248   fMCSource->SetBinEdges(1,binLimPt);
 249   fMCSource->SetBinEdges(2,binLimSource);
 250   fMCSource->Sumw2();
 251   AliDebug(2,"AliHFENonPhotonicElectron: fMCSource");
 252
 253   // ee invariant mass Unlike Sign
 254   const Int_t nDimUSign=6;
 255   Int_t nBinUSign[nDimUSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource,nBinsAngle};
 256   fUSign = new THnSparseF("fUSign","fUSign",nDimUSign,nBinUSign);
 257   fUSign->SetBinEdges(0,binLimPhi);
 258   fUSign->SetBinEdges(1,binLimC);
 259   fUSign->SetBinEdges(2,binLimPt);
 260   fUSign->SetBinEdges(3,binLimInvMass);
 261   fUSign->SetBinEdges(4,binLimSource);
 262   fUSign->SetBinEdges(5,binLimAngle);
 263   fUSign->Sumw2();
 264   AliDebug(2,"AliHFENonPhotonicElectron: fUSign");
 265
 266   // ee invariant mass Like Sign
 267   const Int_t nDimLSign=6;
 268   Int_t nBinLSign[nDimLSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource,nBinsAngle};
 269   fLSign = new THnSparseF("fLSign","fLSign",nDimLSign,nBinLSign);
 270   fLSign->SetBinEdges(0,binLimPhi);
 271   fLSign->SetBinEdges(1,binLimC);
 272   fLSign->SetBinEdges(2,binLimPt);
 273   fLSign->SetBinEdges(3,binLimInvMass);
 274   fLSign->SetBinEdges(4,binLimSource);
 275   fLSign->SetBinEdges(5,binLimAngle);
 276   fLSign->Sumw2();
 277   AliDebug(2,"AliHFENonPhotonicElectron: fLSign");
 278
 279 /*
 280   // ee angle Unlike Sign
 281   const Int_t nDimUSignAngle=3;
 282   Int_t nBinUSignAngle[nDimUSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
 283   fUSignAngle = new THnSparseF("fUSignAngle","fUSignAngle",nDimUSignAngle,nBinUSignAngle);
 284   fUSignAngle->SetBinEdges(0,binLimAngle);
 285   fUSignAngle->SetBinEdges(1,binLimC);
 286   fUSignAngle->SetBinEdges(2,binLimSource);
 287   fUSignAngle->Sumw2();
 288   AliDebug(2,"AliHFENonPhotonicElectron: fUSignAngle");
 289
 290   // ee angle Like Sign
 291   const Int_t nDimLSignAngle=3;
 292   Int_t nBinLSignAngle[nDimLSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
 293   fLSignAngle = new THnSparseF("fLSignAngle","fLSignAngle",nDimLSignAngle,nBinLSignAngle);
 294   fLSignAngle->SetBinEdges(0,binLimAngle);
 295   fLSignAngle->SetBinEdges(1,binLimC);
 296   fLSignAngle->SetBinEdges(2,binLimSource);
 297   fLSignAngle->Sumw2();
 298   AliDebug(2,"AliHFENonPhotonicElectron: fLSignAngle");
 299 */
 300
 301   fListOutput->Add(fMCSource);
 302   fListOutput->Add(fUSign);
 303   fListOutput->Add(fLSign);
 304 //  fListOutput->Add(fUSignAngle);
 305 //  fListOutput->Add(fLSignAngle);
 306
 307 }
 308
 309 //_____________________________________________________________________________________________
 310 void AliHFENonPhotonicElectron::InitRun(const AliVEvent *inputEvent,const AliPIDResponse *pidResponse)
 311 {
 312   //
 313   // Init run
 314   //
 315
 316   if(!pidResponse)
 317   {
 318     AliDebug(1, "Using default PID Response");
 319     Bool_t hasmc = kFALSE;
 320     if(fMCEvent || fAODArrayMCInfo) hasmc=kTRUE;
 321     pidResponse = AliHFEtools::GetDefaultPID(hasmc, inputEvent->IsA() == AliESDEvent::Class());
 322   }
 323
 324   if(!fPIDBackground) return;
 325   fPIDBackground->SetPIDResponse(pidResponse);
 326
 327   if(!fPIDBackground->IsInitialized())
 328   {
 329     // Initialize PID with the given run number
 330     fPIDBackground->InitializePID(inputEvent->GetRunNumber());
 331   }
 332
 333 }
 334
 335 //_____________________________________________________________________________________________
 336 Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent,Int_t binct)
 337 {
 338   //
 339   // Fill the pool of associated tracks
 340   // Return the number of associated tracks
 341   //
 342
 343   Int_t nbtracks = inputEvent->GetNumberOfTracks();
 344
 345   if( fArraytrack )
 346   {
 347     fArraytrack->~TArrayI();
 348     new(fArraytrack) TArrayI(nbtracks);
 349   }
 350   else
 351   {
 352     fArraytrack = new TArrayI(nbtracks);
 353   }
 354
 355   fCounterPoolBackground = 0;
 356
 357   for(Int_t k = 0; k < nbtracks; k++)
 358   {
 359     AliVTrack *track = (AliVTrack *) inputEvent->GetTrack(k);
 360     if(!track) continue;
 361
 362     // Track cuts
 363     Bool_t survivedbackground = kTRUE;
 364     AliAODEvent *aodeventu = dynamic_cast<AliAODEvent *>(inputEvent);
 365
 366     if(aodeventu)
 367     {
 368       /**                               **
 369        *        AOD Analysis             *
 370        **                               **/
 371
 372       AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
 373       if(aodtrack)
 374       {
 375         // filter
 376         if(fUseFilterAOD)
 377         {
 378           if(!(aodtrack->TestFilterBit(fFilter))) survivedbackground = kFALSE;
 379         }
 380
 381         // additional cuts
 382         if(survivedbackground)
 383         {
 384           AliESDtrack esdTrack(aodtrack);
 385
 386           // set the TPC cluster info
 387           esdTrack.SetTPCClusterMap(aodtrack->GetTPCClusterMap());
 388           esdTrack.SetTPCSharedMap(aodtrack->GetTPCSharedMap());
 389           esdTrack.SetTPCPointsF(aodtrack->GetTPCNclsF());
 390           AliAODVertex *vAOD = aodeventu->GetPrimaryVertex();
 391           Double_t bfield = aodeventu->GetMagneticField();
 392           Double_t pos[3],cov[6];
 393           vAOD->GetXYZ(pos);
 394           vAOD->GetCovarianceMatrix(cov);
 395           const AliESDVertex vESD(pos,cov,100.,100);
 396           esdTrack.RelateToVertex(&vESD,bfield,3.);
 397
 398           if(!fHFEBackgroundCuts->IsSelected(&esdTrack))
 399           {
 400             survivedbackground = kFALSE;
 401           }
 402         }
 403       }
 404     }
 405     else
 406     {
 407       /**                               **
 408        *        ESD Analysis             *
 409        **                               **/
 410
 411       AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
 412       if(esdtrack)
 413       {
 414         if(!fHFEBackgroundCuts->IsSelected(esdtrack)) survivedbackground = kFALSE;
 415       }
 416     }
 417
 418     // PID
 419     if(survivedbackground)
 420     {
 421       // PID track cuts
 422       AliHFEpidObject hfetrack2;
 423
 424       if(!aodeventu)    hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 425       else              hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
 426
 427       hfetrack2.SetRecTrack(track);
 428       if(binct>-1)
 429       {
 430         hfetrack2.SetCentrality((Int_t)binct);
 431         AliDebug(2,Form("centrality %d and %d",binct,hfetrack2.GetCentrality()));
 432         hfetrack2.SetPbPb();
 433       }
 434
 435       if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundQA))
 436       {
 437         fArraytrack->AddAt(k,fCounterPoolBackground);
 438         fCounterPoolBackground++;
 439         AliDebug(2,Form("fCounterPoolBackground %d, track %d",fCounterPoolBackground,k));
 440       }
 441     }
 442   } // loop tracks
 443
 444   //printf(Form("Test pool: nbtrack %d, binct %d \n fCounterPoolBackground %d \n",nbtracks,binct,fCounterPoolBackground));
 445
 446   return fCounterPoolBackground;
 447 }
 448
 449 //_____________________________________________________________________________________________
 450 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)
 451 {
 452   //
 453   // Look At Non HFE
 454   //
 455
 456   /***********************************************************************************
 457    *                                                                                    *
 458    *    iTrack1:        index of the tagged electrons in AliVEvent                      *
 459    *    track1:         tagged electron                                                 *
 460    *    vEvent:         event                                                           *
 461    *    weight:         weight in pt if not realistic                                   *
 462    *    binct:          centrality bin                                                  *
 463    *    deltaphi:       phi-phi event plane for v2                                      *
 464    *    source:         MC sources                                                      *
 465    *    indexmother:    MC index mother                                                 *
 466    *                                                                                    *
 467    *                                                                                    *
 468    *    return -1  if  nothing                                                          *
 469    *    return  2  if  opposite         charge          within the mass range           *
 470    *    return  4  if      like         charge          within the mass range           *
 471    *    return  6  if  opposite & like charge           within the mass range           *
 472    *                                                                                    *
 473    ***********************************************************************************/
 474
 475   AliAODEvent *aodeventu = dynamic_cast<AliAODEvent*>(vEvent);
 476   Int_t taggedphotonic = -1;
 477
 478   AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
 479   if(!fArraytrack) return taggedphotonic;
 480   AliDebug(2,Form("process track %d",iTrack1));
 481
 482   //Set Fill-Arrays for THnSparse
 483   Double_t valueMCSource[3]     = { -1, -1, source};                                            //Centrality    Pt              Source
 484   Double_t valueSign[6]         = { deltaphi, binct, track1->Pt(), -1, source, -1};             //DeltaPhi      Centrality      Pt      InvariantMass   Source
 485 //  Double_t valueAngle[3]      = { -1, binct, source};                                         //Angle         Centrality      Source
 486
 487   Bool_t uSignPhotonic = kFALSE;
 488   Bool_t lSignPhotonic = kFALSE;
 489   Bool_t hasdcaT1 = kFALSE;
 490   Bool_t hasdcaT2 = kFALSE;
 491
 492   Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()));
 493   Int_t pdg2 = -100;
 494   Int_t numberfound = 0;
 495   Int_t iTrack2 = 0;
 496   Int_t source2 = 0;
 497   Int_t indexmother2 = -1;
 498   Int_t fPDGtrack1 = 0;
 499   Int_t fPDGtrack2 = 0;
 500
 501   Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
 502   Double_t bfield = vEvent->GetMagneticField();
 503   Double_t xt1 = 0; //radial position track 1 at the DCA point
 504   Double_t xt2 = 0; //radial position track 2 at the DCA point
 505   Double_t p1[3] = {0,0,0};
 506   Double_t p2[3] = {0,0,0};
 507   Double_t dca12 = 0;
 508
 509   Double_t angleESD = 0;
 510   Double_t invmassESD = 0;
 511
 512   Double_t chi2OverNDF = 0;
 513   Double_t width = 0;
 514   Double_t angleAOD = 0;
 515   Double_t invmassAOD = 0;
 516
 517   Float_t fCharge1 = 0;
 518   Float_t fCharge2 = 0;
 519
 520   TLorentzVector electron1;
 521   TLorentzVector electron2;
 522   TLorentzVector mother;
 523
 524   AliVTrack *track2;
 525   AliESDtrack *esdtrack1;
 526   AliESDtrack *esdtrack2;
 527   AliKFParticle *ktrack1;
 528   AliKFParticle *ktrack2;
 529   AliKFVertex primV(*(vEvent->GetPrimaryVertex()));
 530
 531   //! FILL MCsource TODO: Use MC information!!!
 532   /**   if(fAODArrayMCInfo)     valueMCSource[3] = {fAODArrayMCInfo};
 533         if(fMCEvent)            valueMCSource[3] = {fMCEvent};          */
 534   fMCSource->Fill(&valueMCSource[0],weight);
 535
 536
 537   for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
 538   {
 539     iTrack2 = fArraytrack->At(idex);
 540     AliDebug(2,Form("track %d",iTrack2));
 541     track2 = (AliVTrack *)vEvent->GetTrack(iTrack2);
 542
 543     if(!track2)
 544     {
 545       printf("ERROR: Could not receive track %d", iTrack2);
 546       continue;
 547     }
 548
 549     if(iTrack2==iTrack1) continue;
 550     AliDebug(2,"Different");
 551
 552     fCharge1 = track1->Charge();                //Charge from track1
 553     fCharge2 = track2->Charge();                //Charge from track2
 554
 555     // Reset the MC info
 556     //valueAngle[2] = source;
 557     valueSign[4] = source;
 558
 559     // track cuts and PID already done
 560
 561     // if MC look
 562     if(fMCEvent || fAODArrayMCInfo)
 563     {
 564       source2    = FindMother(TMath::Abs(track2->GetLabel()), indexmother2);
 565       pdg2       = CheckPdg(TMath::Abs(track2->GetLabel()));
 566
 567       if(source2 >=0 )
 568       {
 569         if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0))
 570         {
 571           if(source == kElectronfromconversion)
 572           {
 573             //valueAngle[2] = kElectronfromconversionboth;
 574             valueSign[4] = kElectronfromconversionboth;
 575             numberfound++;
 576           }
 577
 578           if(source == kElectronfrompi0)
 579           {
 580             //valueAngle[2] = kElectronfrompi0both;
 581             valueSign[4] = kElectronfrompi0both;
 582           }
 583
 584           if(source == kElectronfrometa)
 585           {
 586             //valueAngle[2] = kElectronfrometaboth;
 587             valueSign[4] = kElectronfrometaboth;
 588           }
 589         }
 590       }
 591     }
 592
 593     if(fAlgorithmMA && (!aodeventu))
 594     {
 595       /**                               *
 596        *        ESD-Analysis            *
 597        **                               */
 598
 599       esdtrack1 = dynamic_cast<AliESDtrack *>(track1);                  //ESD-track1
 600       esdtrack2 = dynamic_cast<AliESDtrack *>(track2);                  //ESD-track2
 601       if((!esdtrack1) || (!esdtrack2)) continue;
 602
 603       dca12 = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);              //DCA track1-track2
 604
 605       if(dca12 > fMaxDCA) continue;                                     //! Cut on DCA
 606
 607       //Momento of the track extrapolated to DCA track-track
 608       hasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);         //Track1
 609       hasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);         //Track2
 610       if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
 611
 612       electron1.SetXYZM(esdtrack1->Px(), esdtrack1->Py(), esdtrack1->Pz(), eMass);
 613       electron2.SetXYZM(esdtrack2->Px(), esdtrack2->Py(), esdtrack2->Pz(), eMass);
 614
 615       mother      = electron1 + electron2;
 616       invmassESD  = mother.M();
 617       angleESD    = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
 618
 619       //valueAngle[0] = angleESD;
 620       valueSign[5] = angleESD;
 621
 622       //if((fCharge1*fCharge2)>0.0)     fLSignAngle->Fill(&valueAngle[0],weight);
 623       //else                            fUSignAngle->Fill(&valueAngle[0],weight);
 624
 625       if(angleESD > fMaxOpening3D) continue;                             //! Cut on Opening Angle
 626       valueSign[3] = invmassESD;
 627
 628       if((fCharge1*fCharge2)>0.0)       fLSign->Fill(&valueSign[0],weight);
 629       else                              fUSign->Fill(&valueSign[0],weight);
 630
 631       if(invmassESD > fMaxInvMass) continue;                            //! Cut on Invariant Mass
 632
 633       if((fCharge1*fCharge2)>0.0)       lSignPhotonic=kTRUE;
 634       else                              uSignPhotonic=kTRUE;
 635     }
 636     else
 637     {
 638       /**                               *
 639        *        AOD-AliKF-Analysis      *
 640        **                               */
 641
 642       fPDGtrack1 = 11;
 643       fPDGtrack2 = 11;
 644
 645       if(fCharge1>0) fPDGtrack1 = -11;
 646       if(fCharge2>0) fPDGtrack2 = -11;
 647
 648       ktrack1 = new AliKFParticle(*track1, fPDGtrack1);
 649       ktrack2 = new AliKFParticle(*track2, fPDGtrack2);
 650       AliKFParticle recoGamma(*ktrack1,*ktrack2);
 651
 652       if(recoGamma.GetNDF()<1) continue;                                //! Cut on Reconstruction
 653
 654       chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
 655       if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
 656
 657       // DCA
 658       //Double_t dca12 = ktrack1.GetDistanceFromParticle(ktrack2);
 659       //if(dca12 > fMaxDCA) continue;
 660
 661       // if set mass constraint
 662       if(fSetMassConstraint) //&& pVtx)
 663       {
 664         primV += recoGamma;
 665         primV -= *ktrack1;
 666         primV -= *ktrack2;
 667         recoGamma.SetProductionVertex(primV);
 668         recoGamma.SetMassConstraint(0,0.0001);
 669       }
 670
 671       recoGamma.GetMass(invmassAOD,width);
 672       angleAOD = ktrack1->GetAngle(*ktrack2);
 673
 674       //valueAngle[0] = angleAOD;
 675       valueSign[5] = angleAOD;
 676
 677       //if((fCharge1*fCharge2)>0.0)     fLSignAngle->Fill(&valueAngle[0],weight);
 678       //else                            fUSignAngle->Fill(&valueAngle[0],weight);
 679
 680       if(angleAOD > fMaxOpening3D) continue;                            //! Cut on Opening Angle
 681
 682       valueSign[3] = invmassAOD;
 683
 684       if((fCharge1*fCharge2)>0.0)       fLSign->Fill(&valueSign[0],weight);
 685       else                              fUSign->Fill(&valueSign[0],weight);
 686
 687       if(invmassAOD > fMaxInvMass) continue;                            //! Cut on Invariant Mass
 688
 689       if((fCharge1*fCharge2)>0.0)       lSignPhotonic=kTRUE;
 690       else                              uSignPhotonic=kTRUE;
 691     }
 692   }
 693
 694   if( uSignPhotonic &&  lSignPhotonic) taggedphotonic = 6;
 695   if(!uSignPhotonic &&  lSignPhotonic) taggedphotonic = 4;
 696   if( uSignPhotonic && !lSignPhotonic) taggedphotonic = 2;
 697
 698   return taggedphotonic;
 699 }
 700
 701 //_________________________________________________________________________
 702 Int_t AliHFENonPhotonicElectron::FindMother(Int_t tr, Int_t &indexmother){
 703   //
 704   // Find the mother if MC
 705   //
 706
 707   if(!fMCEvent && !fAODArrayMCInfo) return 0;
 708
 709   Int_t pdg = CheckPdg(tr);
 710   if(TMath::Abs(pdg)!= 11)
 711   {
 712     indexmother = -1;
 713     return kNoElectron;
 714   }
 715
 716   indexmother = IsMotherGamma(tr);
 717   if(indexmother > 0) return kElectronfromconversion;
 718   indexmother = IsMotherPi0(tr);
 719   if(indexmother > 0) return kElectronfrompi0;
 720   indexmother = IsMotherC(tr);
 721   if(indexmother > 0) return kElectronfromC;
 722   indexmother = IsMotherB(tr);
 723   if(indexmother > 0) return kElectronfromB;
 724   indexmother = IsMotherEta(tr);
 725   if(indexmother > 0) return kElectronfrometa;
 726
 727   return kElectronfromother;
 728 }
 729
 730 //________________________________________________________________________________________________
 731 Int_t AliHFENonPhotonicElectron::CheckPdg(Int_t tr) {
 732
 733   //
 734   // Return the pdg of the particle
 735   //
 736
 737   Int_t pdgcode = -1;
 738   if(tr < 0) return pdgcode;
 739
 740   if(fMCEvent)
 741   {
 742     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 743     if(!mctrack) return -1;
 744     AliMCParticle *mctrackesd = NULL;
 745     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
 746     pdgcode = mctrackesd->PdgCode();
 747   }
 748
 749   if(fAODArrayMCInfo)
 750   {
 751     if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
 752     AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 753     if(!mctrackaod) return pdgcode;
 754     pdgcode = mctrackaod->GetPdgCode();
 755   }
 756
 757   return pdgcode;
 758 }
 759
 760 //_______________________________________________________________________________________________
 761 Int_t AliHFENonPhotonicElectron::IsMotherGamma(Int_t tr) {
 762
 763   //
 764   // Return the lab of gamma mother or -1 if not gamma
 765   //
 766
 767   if(tr < 0) return -1;
 768
 769   if(fMCEvent)
 770   {
 771     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 772     if(!mctrack) return -1;
 773     AliMCParticle *mctrackesd = NULL;
 774     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
 775     TParticle *particle = 0x0;
 776     particle = mctrackesd->Particle();
 777
 778     // Take mother
 779     if(!particle) return -1;
 780     Int_t imother   = particle->GetFirstMother();
 781     if(imother < 0) return -1;
 782     AliMCParticle *mothertrack = NULL;
 783     if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
 784     TParticle * mother = mothertrack->Particle();
 785     if(!mother) return -1;
 786
 787     // Check gamma
 788     Int_t pdg = mother->GetPdgCode();
 789     if(TMath::Abs(pdg) == 22) return imother;
 790     if(TMath::Abs(pdg) == 11)
 791     {
 792       return IsMotherGamma(imother);
 793     }
 794
 795     return -1;
 796   }
 797
 798   if(fAODArrayMCInfo)
 799   {
 800     if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
 801     AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 802     if(!mctrackaod) return -1;
 803
 804     // Take mother
 805     Int_t imother = mctrackaod->GetMother();
 806     if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
 807     AliAODMCParticle *mothertrack = NULL;
 808     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
 809
 810     // Check gamma
 811     Int_t pdg = mothertrack->GetPdgCode();
 812     if(TMath::Abs(pdg) == 22) return imother;
 813     if(TMath::Abs(pdg) == 11)
 814     {
 815       return IsMotherGamma(imother);
 816     }
 817
 818     return -1;
 819   }
 820
 821   return -1;
 822 }
 823
 824 //________________________________________________________________________________________________
 825 Int_t AliHFENonPhotonicElectron::IsMotherPi0(Int_t tr) {
 826
 827   //
 828   // Return the lab of pi0 mother or -1 if not pi0
 829   //
 830
 831   if(tr < 0) return -1;
 832
 833   if(fMCEvent)
 834   {
 835     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 836     if(!mctrack) return -1;
 837     AliMCParticle *mctrackesd = NULL;
 838     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
 839     TParticle *particle = 0x0;
 840     particle = mctrackesd->Particle();
 841     // Take mother
 842     if(!particle) return -1;
 843     Int_t imother   = particle->GetFirstMother();
 844     if(imother < 0) return -1;
 845     AliMCParticle *mothertrack = NULL;
 846     if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
 847     TParticle * mother = mothertrack->Particle();
 848     if(!mother) return -1;
 849     // Check pi0
 850     Int_t pdg = mother->GetPdgCode();
 851     if(TMath::Abs(pdg) == 111) return imother;
 852     if(TMath::Abs(pdg) == 11)
 853     {
 854       return IsMotherPi0(imother);
 855     }
 856
 857     return -1;
 858   }
 859
 860   if(fAODArrayMCInfo)  {
 861
 862     if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
 863     AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 864     if(!mctrackaod) return -1;
 865
 866     // Take mother
 867     Int_t imother = mctrackaod->GetMother();
 868     if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
 869     AliAODMCParticle *mothertrack = NULL;
 870     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
 871     // Check pi0
 872     Int_t pdg = mothertrack->GetPdgCode();
 873     if(TMath::Abs(pdg) == 111) return imother;
 874     if(TMath::Abs(pdg) == 11)
 875     {
 876       return IsMotherPi0(imother);
 877     }
 878
 879     return -1;
 880   }
 881
 882   return -1;
 883 }
 884 //________________________________________________________________________________________________
 885 Int_t AliHFENonPhotonicElectron::IsMotherC(Int_t tr) {
 886
 887   //
 888   // Return the lab of signal mother or -1 if not from C
 889   //
 890
 891   if(tr < 0) return -1;
 892
 893   if(fMCEvent)
 894   {
 895     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 896     if(!mctrack) return -1;
 897     AliMCParticle *mctrackesd = NULL;
 898     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
 899     TParticle *particle = 0x0;
 900     particle = mctrackesd->Particle();
 901
 902     // Take mother
 903     if(!particle) return -1;
 904     Int_t imother   = particle->GetFirstMother();
 905     if(imother < 0) return -1;
 906     AliMCParticle *mothertrack = NULL;
 907     if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
 908     TParticle * mother = mothertrack->Particle();
 909     if(!mother) return -1;
 910
 911     // Check C
 912     Int_t pdg = mother->GetPdgCode();
 913     if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
 914     if(TMath::Abs(pdg) == 11)
 915     {
 916       return IsMotherC(imother);
 917     }
 918
 919     return -1;
 920   }
 921
 922   if(fAODArrayMCInfo)
 923   {
 924     if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
 925     AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 926     if(!mctrackaod) return -1;
 927
 928     // Take mother
 929     Int_t imother = mctrackaod->GetMother();
 930     if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
 931     AliAODMCParticle *mothertrack = NULL;
 932     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
 933
 934     // Check C
 935     Int_t pdg = mothertrack->GetPdgCode();
 936     if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
 937     if(TMath::Abs(pdg) == 11)
 938     {
 939       return IsMotherC(imother);
 940     }
 941
 942     return -1;
 943   }
 944
 945   return -1;
 946 }
 947 //_______________________________________________________________________________________________
 948 Int_t AliHFENonPhotonicElectron::IsMotherB(Int_t tr) {
 949
 950   //
 951   // Return the lab of signal mother or -1 if not B
 952   //
 953
 954   if(tr < 0) return -1;
 955
 956   if(fMCEvent)
 957   {
 958     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
 959     if(!mctrack) return -1;
 960     AliMCParticle *mctrackesd = NULL;
 961     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
 962     TParticle *particle = 0x0;
 963     particle = mctrackesd->Particle();
 964
 965     // Take mother
 966     if(!particle) return -1;
 967     Int_t imother   = particle->GetFirstMother();
 968     if(imother < 0) return -1;
 969     AliMCParticle *mothertrack = NULL;
 970     if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
 971     TParticle * mother = mothertrack->Particle();
 972     if(!mother) return -1;
 973
 974     // Check B
 975     Int_t pdg = mother->GetPdgCode();
 976     if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
 977     if(TMath::Abs(pdg) == 11)
 978     {
 979       return IsMotherB(imother);
 980     }
 981
 982     return -1;
 983   }
 984
 985   if(fAODArrayMCInfo)
 986   {
 987     if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
 988     AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
 989     if(!mctrackaod) return -1;
 990
 991     // Take mother
 992     Int_t imother = mctrackaod->GetMother();
 993     if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
 994     AliAODMCParticle *mothertrack = NULL;
 995     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
 996     // Check B
 997     Int_t pdg = mothertrack->GetPdgCode();
 998     if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
 999     if(TMath::Abs(pdg) == 11)
1000     {
1001       return IsMotherB(imother);
1002     }
1003
1004     return -1;
1005   }
1006
1007   return -1;
1008 }
1009
1010 //_______________________________________________________________________________________________
1011 Int_t AliHFENonPhotonicElectron::IsMotherEta(Int_t tr) {
1012
1013   //
1014   // Return the lab of eta mother or -1 if not eta
1015   //
1016
1017   if(tr < 0) return -1;
1018
1019   if(fMCEvent)
1020   {
1021     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
1022     if(!mctrack) return -1;
1023     AliMCParticle *mctrackesd = NULL;
1024     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
1025     TParticle *particle = 0x0;
1026     particle = mctrackesd->Particle();
1027
1028     // Take mother
1029     if(!particle) return -1;
1030     Int_t imother   = particle->GetFirstMother();
1031     if(imother < 0) return -1;
1032     AliMCParticle *mothertrack = NULL;
1033     if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
1034     TParticle * mother = mothertrack->Particle();
1035     if(!mother) return -1;
1036
1037     // Check eta
1038     Int_t pdg = mother->GetPdgCode();
1039     if(TMath::Abs(pdg) == 221) return imother;
1040     if(TMath::Abs(pdg) == 11)
1041     {
1042       return IsMotherEta(imother);
1043     }
1044
1045     return -1;
1046   }
1047
1048   if(fAODArrayMCInfo)
1049   {
1050     if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
1051     AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
1052     if(!mctrackaod) return -1;
1053
1054     // Take mother
1055     Int_t imother = mctrackaod->GetMother();
1056     if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
1057     AliAODMCParticle *mothertrack = NULL;
1058     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
1059
1060     // Check eta
1061     Int_t pdg = mothertrack->GetPdgCode();
1062     if(TMath::Abs(pdg) == 221) return imother;
1063     if(TMath::Abs(pdg) == 11)
1064     {
1065       return IsMotherEta(imother);
1066     }
1067
1068     return -1;
1069   }
1070
1071   return -1;
1072 }