PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalEP.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 // Class for heavy-flavour electron v2 with EMCal triggered events
  17 // Author: Denise Godoy
  18
  19
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH2F.h"
  23 #include "TMath.h"
  24 #include "TCanvas.h"
  25 #include "THnSparse.h"
  26 #include "TLorentzVector.h"
  27 #include "TString.h"
  28 #include "TFile.h"
  29
  30 #include "AliAnalysisTask.h"
  31 #include "AliAnalysisManager.h"
  32
  33 #include "AliESDEvent.h"
  34 #include "AliESDHandler.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAODHandler.h"
  37 #include "AliVEvent.h"
  38
  39 #include "AliAnalysisTaskFlowTPCEMCalEP.h"
  40 #include "TGeoGlobalMagField.h"
  41 #include "AliLog.h"
  42 #include "AliAnalysisTaskSE.h"
  43 #include "TRefArray.h"
  44 #include "TVector.h"
  45 #include "AliESDInputHandler.h"
  46 #include "AliESDpid.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliPhysicsSelection.h"
  49 #include "AliESDCaloCluster.h"
  50 #include "AliAODCaloCluster.h"
  51 #include "AliEMCALRecoUtils.h"
  52 #include "AliEMCALGeometry.h"
  53 #include "AliGeomManager.h"
  54 #include "stdio.h"
  55 #include "TGeoManager.h"
  56 #include "iostream"
  57 #include "fstream"
  58
  59 #include "AliEMCALTrack.h"
  60 #include "AliMagF.h"
  61
  62 #include "AliKFParticle.h"
  63 #include "AliKFVertex.h"
  64
  65 #include "AliMCEventHandler.h"
  66 #include "AliMCEvent.h"
  67 #include "AliMCParticle.h"
  68 #include "AliStack.h"
  69
  70 #include "AliPID.h"
  71 #include "AliPIDResponse.h"
  72 #include "AliHFEcontainer.h"
  73 #include "AliHFEcuts.h"
  74 #include "AliHFEpid.h"
  75 #include "AliHFEpidBase.h"
  76 #include "AliHFEpidQAmanager.h"
  77 #include "AliHFEtools.h"
  78 #include "AliCFContainer.h"
  79 #include "AliCFManager.h"
  80
  81 #include "AliEventplane.h"
  82 #include "AliCentrality.h"
  83
  84 #include "AliSelectNonHFE.h"
  85
  86
  87 ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
  88 //________________________________________________________________________
  89 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
  90   : AliAnalysisTaskSE(name)
  91   ,fESD(0)
  92   ,fAOD(0)
  93   ,fVevent(0)
  94   ,fpidResponse(0)
  95   ,fMC(0)
  96   ,fOutputList(0)
  97   ,fTrackCuts(0)
  98   ,fCuts(0)
  99   ,fNonHFE(0)
 100   ,fIdentifiedAsOutInz(kFALSE)
 101   ,fPassTheEventCut(kFALSE)
 102   ,fRejectKinkMother(kFALSE)
 103   ,fIsMC(kFALSE)
 104   ,fIsAOD(kFALSE)
 105   ,fSetMassConstraint(kFALSE)
 106   ,fVz(0.0)
 107   ,fCFM(0)
 108   ,fPID(0)
 109   ,fPIDqa(0)
 110   ,fOpeningAngleCut(1000.)
 111   ,fInvmassCut(0.05)
 112   ,fChi2Cut(3.5)
 113   ,fDCAcut(999)
 114   ,fnonHFEalgorithm("KF")
 115   ,fNoEvents(0)
 116   ,fTrkpt(0)
 117   ,fTrkEovPBef(0)
 118   ,fTrkEovPAft(0)
 119   ,fdEdxBef(0)
 120   ,fdEdxAft(0)
 121   ,fPhotoElecPt(0)
 122   ,fSemiInclElecPt(0)
 123   ,fTrackPtBefTrkCuts(0)
 124   ,fTrackPtAftTrkCuts(0)
 125   ,fTPCnsigma(0)
 126   ,fCent(0)
 127   ,fTPCsubEPres(0)
 128   ,fEPres(0)
 129   ,fCorr(0)
 130   ,fD0_e(0)
 131   ,fTot_pi0e(0)
 132   ,fPhot_pi0e(0)
 133   ,fPhotBCG_pi0e(0)
 134   ,fTot_etae(0)
 135   ,fPhot_etae(0)
 136   ,fPhotBCG_etae(0)
 137   ,fInvMass(0)
 138   ,fInvMassBack(0)
 139   ,fDCA(0)
 140   ,fDCABack(0)
 141   ,fOpAngle(0)
 142   ,fOpAngleBack(0)
 143 {
 144   //Named constructor
 145
 146   for(Int_t k = 0; k < 3; k++) {
 147     fevPlaneV0[k] = NULL;
 148     feTPCV2[k] = NULL;
 149     feV2[k] = NULL;
 150     fChargPartV2[k] = NULL;
 151     fMtcPartV2[k] = NULL;
 152
 153     fPi0Pt[k] = NULL;
 154     fEtaPt[k] = NULL;
 155     fInvmassLS[k] = NULL;
 156     fInvmassULS[k] = NULL;
 157     fOpeningAngleLS[k] = NULL;
 158     fOpeningAngleULS[k] = NULL;
 159   }
 160
 161   for(Int_t i=0; i<3; i++) {
 162     for(Int_t j=0; j<8; j++) {
 163       for(Int_t k=0; k<4; k++) {
 164         fEoverPsig[i][j][k] = NULL;
 165         fEoverPuls[i][j][k] = NULL;
 166         fEoverPls[i][j][k] = NULL;
 167         fEoverPbcg[i][j][k] = NULL;
 168       }
 169     }
 170   }
 171
 172   for(Int_t k = 0; k < 6; k++) {
 173     fDe[k]= NULL;
 174     fD0e[k]= NULL;
 175     fDpluse[k]= NULL;
 176     fDminuse[k]= NULL;
 177   }
 178
 179   fPID = new AliHFEpid("hfePid");
 180   fTrackCuts = new AliESDtrackCuts();
 181   fNonHFE = new AliSelectNonHFE();
 182
 183   InitParameters();
 184
 185   // Define input and output slots here
 186   // Input slot #0 works with a TChain
 187   DefineInput(0, TChain::Class());
 188   // Output slot #0 id reserved by the base class for AOD
 189   // Output slot #1 writes into a TH1 container
 190   // DefineOutput(1, TH1I::Class());
 191   DefineOutput(1, TList::Class());
 192   //  DefineOutput(3, TTree::Class());
 193 }
 194
 195 //________________________________________________________________________
 196 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
 197   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
 198   ,fESD(0)
 199   ,fAOD(0)
 200   ,fVevent(0)
 201   ,fpidResponse(0)
 202   ,fMC(0)
 203   ,fOutputList(0)
 204   ,fTrackCuts(0)
 205   ,fCuts(0)
 206   ,fNonHFE(0)
 207   ,fIdentifiedAsOutInz(kFALSE)
 208   ,fPassTheEventCut(kFALSE)
 209   ,fRejectKinkMother(kFALSE)
 210   ,fIsMC(kFALSE)
 211   ,fIsAOD(kFALSE)
 212   ,fSetMassConstraint(kFALSE)
 213   ,fVz(0.0)
 214   ,fCFM(0)
 215   ,fPID(0)
 216   ,fPIDqa(0)
 217   ,fOpeningAngleCut(1000.)
 218   ,fInvmassCut(0.05)
 219   ,fChi2Cut(3.5)
 220   ,fDCAcut(999)
 221   ,fnonHFEalgorithm("KF")
 222   ,fNoEvents(0)
 223   ,fTrkpt(0)
 224   ,fTrkEovPBef(0)
 225   ,fTrkEovPAft(0)
 226   ,fdEdxBef(0)
 227   ,fdEdxAft(0)
 228   ,fPhotoElecPt(0)
 229   ,fSemiInclElecPt(0)
 230   ,fTrackPtBefTrkCuts(0)
 231   ,fTrackPtAftTrkCuts(0)
 232   ,fTPCnsigma(0)
 233   ,fCent(0)
 234   ,fTPCsubEPres(0)
 235   ,fEPres(0)
 236   ,fCorr(0)
 237   ,fD0_e(0)
 238   ,fTot_pi0e(0)
 239   ,fPhot_pi0e(0)
 240   ,fPhotBCG_pi0e(0)
 241   ,fTot_etae(0)
 242   ,fPhot_etae(0)
 243   ,fPhotBCG_etae(0)
 244   ,fInvMass(0)
 245   ,fInvMassBack(0)
 246   ,fDCA(0)
 247   ,fDCABack(0)
 248   ,fOpAngle(0)
 249   ,fOpAngleBack(0)
 250 {
 251
 252   //Default constructor
 253
 254   for(Int_t k = 0; k < 3; k++) {
 255     fevPlaneV0[k] = NULL;
 256     feTPCV2[k] = NULL;
 257     feV2[k] = NULL;
 258     fChargPartV2[k] = NULL;
 259     fMtcPartV2[k] = NULL;
 260
 261     fPi0Pt[k] = NULL;
 262     fEtaPt[k] = NULL;
 263     fInvmassLS[k] = NULL;
 264     fInvmassULS[k] = NULL;
 265     fOpeningAngleLS[k] = NULL;
 266     fOpeningAngleULS[k] = NULL;
 267   }
 268
 269   for(Int_t i=0; i<3; i++) {
 270     for(Int_t j=0; j<8; j++) {
 271       for(Int_t k=0; k<4; k++) {
 272         fEoverPsig[i][j][k] = NULL;
 273         fEoverPuls[i][j][k] = NULL;
 274         fEoverPls[i][j][k] = NULL;
 275         fEoverPbcg[i][j][k] = NULL;
 276       }
 277     }
 278   }
 279
 280   for(Int_t k = 0; k < 6; k++) {
 281     fDe[k]= NULL;
 282     fD0e[k]= NULL;
 283     fDpluse[k]= NULL;
 284     fDminuse[k]= NULL;
 285   }
 286
 287   fPID = new AliHFEpid("hfePid");
 288   fTrackCuts = new AliESDtrackCuts();
 289   fNonHFE = new AliSelectNonHFE();
 290
 291   InitParameters();
 292
 293
 294   // Constructor
 295   // Define input and output slots here
 296   // Input slot #0 works with a TChain
 297   DefineInput(0, TChain::Class());
 298   // Output slot #0 id reserved by the base class for AOD
 299   // Output slot #1 writes into a TH1 container
 300   // DefineOutput(1, TH1I::Class());
 301   DefineOutput(1, TList::Class());
 302   //DefineOutput(3, TTree::Class());
 303 }
 304 //_________________________________________
 305
 306 AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
 307 {
 308   //Destructor
 309
 310   delete fOutputList;
 311   delete fPID;
 312   delete fCFM;
 313   delete fPIDqa;
 314   delete fTrackCuts;
 315 }
 316 //_________________________________________
 317
 318 void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
 319 {
 320   //Main loop
 321   //Called for each event
 322
 323   // create pointer to event
 324   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 325   if (!fESD){
 326     printf("ERROR: fESD not available\n");
 327     return;
 328   }
 329
 330   fVevent = dynamic_cast<AliVEvent*>(InputEvent());
 331   if(!fVevent){
 332     printf("ERROR: fVEvent not available\n");
 333     return;
 334   }
 335
 336   if(!fCuts){
 337     AliError("HFE cuts not available");
 338     return;
 339   }
 340
 341   if(!fPID->IsInitialized()){
 342     // Initialize PID with the given run number
 343     AliWarning("PID not initialised, get from Run no");
 344     if(fIsAOD)fPID->InitializePID(fAOD->GetRunNumber());
 345     if(!fIsAOD) fPID->InitializePID(fESD->GetRunNumber());
 346   }
 347
 348   Bool_t SelColl = kTRUE;
 349   if(GetCollisionCandidates()==AliVEvent::kAny)
 350   {
 351      SelColl = kFALSE;
 352      TString firedTrigger;
 353      firedTrigger = fESD->GetFiredTriggerClasses();
 354      if(firedTrigger.Contains("CVLN_B2-B-NOPF-ALLNOTRD") || firedTrigger.Contains("CVLN_R1-B-NOPF-ALLNOTRD") || firedTrigger.Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))SelColl=kTRUE;
 355      if(!SelColl)return;
 356   }
 357
 358   if(fIsMC)fMC = MCEvent();
 359   AliStack* stack = NULL;
 360   if(fIsMC && fMC) stack = fMC->Stack();
 361
 362   Int_t fNOtrks =fESD->GetNumberOfTracks();
 363   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 364
 365   Double_t pVtxZ = -999;
 366   pVtxZ = pVtx->GetZ();
 367
 368   if(TMath::Abs(pVtxZ)>10) return;
 369   fNoEvents->Fill(0);
 370
 371   if(fNOtrks<2) return;
 372
 373   fpidResponse = fInputHandler->GetPIDResponse();
 374   if(!fpidResponse){
 375     AliDebug(1, "Using default PID Response");
 376     fpidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 377   }
 378
 379   fPID->SetPIDResponse(fpidResponse);
 380
 381   fCFM->SetRecEventInfo(fVevent);
 382
 383   Float_t cent = -1.;
 384   Int_t iPt=8, iCent=3, iDeltaphi=4;
 385   AliCentrality *centrality = fESD->GetCentrality();
 386   cent = centrality->GetCentralityPercentile("V0M");
 387   fCent->Fill(cent);
 388
 389   if (cent>=0  && cent<10) iCent=0;
 390   if (cent>=10 && cent<20) iCent=1;
 391   if (cent>=20 && cent<40) iCent=2;
 392   if (cent<0 || cent>=40) return;
 393
 394   //Event planes
 395
 396   Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
 397   if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
 398
 399   Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
 400   if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
 401
 402   Double_t evPlaneV0 = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0",fESD,2));
 403   if(evPlaneV0 > TMath::Pi()) evPlaneV0 = evPlaneV0 - TMath::Pi();
 404   fevPlaneV0[iCent]->Fill(evPlaneV0);
 405
 406   AliEventplane* esdTPCep = fESD->GetEventplane();
 407   TVector2 *standardQ = 0x0;
 408   Double_t qx = -999., qy = -999.;
 409   standardQ = esdTPCep->GetQVector();
 410   if(!standardQ)return;
 411
 412   qx = standardQ->X();
 413   qy = standardQ->Y();
 414
 415   TVector2 qVectorfortrack;
 416   qVectorfortrack.Set(qx,qy);
 417   Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
 418
 419   //Event plane resolutions
 420
 421   // --> 2 subevent method (only for TPC EP)
 422
 423   TVector2 *qsub1a = esdTPCep->GetQsub1();
 424   TVector2 *qsub2a = esdTPCep->GetQsub2();
 425   Double_t evPlaneResTPC = -999.;
 426   if(qsub1a && qsub2a){
 427     evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
 428   }
 429
 430   fTPCsubEPres->Fill(evPlaneResTPC);
 431
 432   // --> 3 event method (V0, V0A, and V0C EP)
 433
 434   Double_t Qx2pos = 0., Qy2pos = 0., Qx2neg = 0., Qy2neg = 0., Qweight = 1;
 435
 436   for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
 437
 438     AliVParticle* vparticle = fVevent->GetTrack(iTracks);
 439     if (!vparticle){
 440       printf("ERROR: Could not receive track %d\n", iTracks);
 441       continue;
 442     }
 443
 444     AliESDtrack *trackEP = dynamic_cast<AliESDtrack*>(vparticle);
 445
 446     if (!trackEP) {
 447       printf("ERROR: Could not receive track %d\n", iTracks);
 448      continue;
 449     }
 450
 451     if(TMath::Abs(trackEP->Eta())>0.8 || trackEP->Pt() < 0.15 || trackEP->Pt() > 4) continue;
 452
 453     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, trackEP)) continue;
 454
 455     if (trackEP->Pt() < 2) Qweight = trackEP->Pt()/2;
 456     if (trackEP->Pt() >= 2) Qweight = 1;
 457
 458
 459     if(trackEP->Eta()>0 && trackEP->Eta()<0.8){
 460       Qx2pos += Qweight*TMath::Cos(2*trackEP->Phi());
 461       Qy2pos += Qweight*TMath::Sin(2*trackEP->Phi());
 462     }
 463     if(trackEP->Eta()<0 && trackEP->Eta()>-0.8){
 464       Qx2neg += Qweight*TMath::Cos(2*trackEP->Phi());
 465       Qy2neg += Qweight*TMath::Sin(2*trackEP->Phi());
 466     }
 467   }//track loop only for EP
 468
 469   Double_t evPlaneTPCneg = TMath::ATan2(Qy2neg, Qx2neg)/2;
 470   Double_t evPlaneTPCpos = TMath::ATan2(Qy2pos, Qx2pos)/2;
 471
 472   Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0,evPlaneTPCpos),
 473                           GetCos2DeltaPhi(evPlaneV0,evPlaneTPCneg),
 474                           GetCos2DeltaPhi(evPlaneTPCpos,evPlaneTPCneg),cent};
 475   fEPres->Fill(evPlaneRes);
 476
 477   // Selection of primary pi0 and eta in MC to compute the weight
 478   if(fIsMC && fMC && stack){
 479     Int_t nParticles = stack->GetNtrack();
 480     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 481       TParticle* particle = stack->Particle(iParticle);
 482       int fPDG = particle->GetPdgCode();
 483       double pTMC = particle->Pt();
 484
 485       Double_t etaMC = particle->Eta();
 486       if (TMath::Abs(etaMC)>1.2)continue;
 487
 488       Bool_t isMotherPrimary = IsPi0EtaPrimary(particle,stack);
 489       Bool_t isFromLMdecay = IsPi0EtaFromLMdecay(particle,stack);
 490       Bool_t isFromHFdecay = IsPi0EtaFromHFdecay(particle,stack);
 491
 492       if (!isMotherPrimary || isFromHFdecay || isFromLMdecay) continue;
 493
 494       if(fPDG==111) fPi0Pt[iCent]->Fill(pTMC); //pi0
 495       if(fPDG==221) fEtaPt[iCent]->Fill(pTMC); //eta
 496     }
 497   }//MC
 498
 499   Double_t ptRange[9] = {1.5,2,2.5,3,4,6,8,10,13};
 500   Double_t deltaPhiRange[4];
 501   for(Int_t j=0;j<4;j++){
 502     deltaPhiRange[j] = j*(TMath::Pi()/4);
 503   }
 504
 505   // Track loop
 506   for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
 507
 508     AliVParticle* vparticle = fVevent->GetTrack(iTracks);
 509     if (!vparticle){
 510       printf("ERROR: Could not receive track %d\n", iTracks);
 511       continue;
 512     }
 513
 514     AliESDtrack *track = dynamic_cast<AliESDtrack*>(vparticle);
 515
 516     if(!track) continue;
 517
 518     if (TMath::Abs(track->Eta())>0.7) continue;
 519
 520     fTrackPtBefTrkCuts->Fill(track->Pt());
 521
 522     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 523
 524     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 525       if(track->GetKinkIndex(0) != 0) continue;
 526     }
 527
 528     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 529
 530     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 531
 532     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 533
 534     fTrackPtAftTrkCuts->Fill(track->Pt());
 535
 536     Double_t clsE=-9.,p=-99.,EovP=-99.,pt=-99.,dEdx=-99.,fTPCnSigma=9.,phi=-9.,m02=-9.,m20=-9.,fEMCalnSigma=9.,dphi=9.,cosdphi=9.;
 537
 538     pt = track->Pt();
 539     if(pt<1.5) continue;
 540     fTrkpt->Fill(pt);
 541     for(Int_t i=0;i<8;i++) {
 542       if (pt>=ptRange[i] && pt<ptRange[i+1]){
 543         iPt=i;
 544         continue;
 545       }
 546     }
 547
 548     Int_t clsId = track->GetEMCALcluster();
 549     if (clsId>0){
 550       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
 551       if(cluster && cluster->IsEMCAL()){
 552         clsE = cluster->E();
 553         m20 = cluster->GetM20();
 554         m02 = cluster->GetM02();
 555       }
 556     }
 557
 558     p = track->P();
 559     phi = track->Phi();
 560     dEdx = track->GetTPCsignal();
 561     EovP = clsE/p;
 562     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 563     fEMCalnSigma = GetSigmaEMCal(EovP, pt, cent);
 564     fdEdxBef->Fill(p,dEdx);
 565     fTPCnsigma->Fill(p,fTPCnSigma);
 566
 567
 568     dphi = GetDeltaPhi(phi,evPlaneV0);
 569     cosdphi = GetCos2DeltaPhi(phi,evPlaneV0);
 570     for(Int_t i=0;i<3;i++) {
 571       if (dphi>=deltaPhiRange[i] && dphi<deltaPhiRange[i+1]){
 572         iDeltaphi=i;
 573         continue;
 574       }
 575     }
 576
 577     Bool_t fFlagPhotonicElec = kFALSE;
 578     Bool_t fFlagPhotonicElecBCG = kFALSE;
 579     Double_t weight = 1.;
 580
 581     SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG,weight,iCent);
 582
 583     Int_t partPDG = -99;
 584     Double_t partPt = -99.;
 585     Bool_t MChijing;
 586
 587     if(fIsMC && fMC && stack){
 588       if(fTPCnSigma < -1 || fTPCnSigma > 3) continue;
 589       Int_t label = track->GetLabel();
 590       if(label!=0){
 591         TParticle *particle = stack->Particle(label);
 592         if(particle){
 593           partPDG = particle->GetPdgCode();
 594           partPt = particle->Pt();
 595           if (TMath::Abs(partPDG)!=11) continue;
 596
 597           MChijing = fMC->IsFromBGEvent(label);
 598           int iHijing = 1;
 599           if(!MChijing) iHijing = 0; // 0 if enhanced sample
 600
 601           Bool_t pi0Decay = IsElectronFromPi0(particle,stack,weight,cent);
 602           Bool_t etaDecay = IsElectronFromEta(particle,stack,weight,cent);
 603
 604           SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG,weight,iCent);
 605
 606           if (pi0Decay){
 607             fTot_pi0e->Fill(partPt,weight);
 608             if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
 609             if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
 610           }
 611           if (etaDecay){
 612             fTot_etae->Fill(partPt,weight);
 613             if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
 614             if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
 615           }
 616         }// end particle
 617       }// end label
 618     }//end MC
 619
 620     if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
 621     Int_t pidpassed = 1;
 622
 623     //--- track accepted
 624     AliHFEpidObject hfetrack;
 625     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 626     hfetrack.SetRecTrack(track);
 627     hfetrack.SetPbPb();
 628     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 629
 630     if (m20>0.02 && m02>0.02){
 631       Double_t corr[7]={(Double_t)iCent,(Double_t)iPt,fTPCnSigma,fEMCalnSigma,m02,dphi,cosdphi};
 632       fCorr->Fill(corr);
 633     }
 634
 635     fChargPartV2[iCent]->Fill(iPt,cosdphi);
 636     if (clsE>0) fMtcPartV2[iCent]->Fill(iPt,cosdphi);
 637
 638     if (pidpassed==0) continue;
 639
 640     if (fTPCnSigma>=-5 && fTPCnSigma<-3.2) fEoverPbcg[iCent][iPt][iDeltaphi]->Fill(EovP);
 641     if (fTPCnSigma>=-0.5 && fTPCnSigma<3) feTPCV2[iCent]->Fill(iPt,cosdphi);
 642     if (fTPCnSigma>=-0.5 && fTPCnSigma<3 && fEMCalnSigma>-1 && fEMCalnSigma<3) feV2[iCent]->Fill(iPt,cosdphi);
 643
 644     fTrkEovPAft->Fill(pt,EovP);
 645     fdEdxAft->Fill(p,dEdx);
 646
 647     if(fFlagPhotonicElec){
 648       fPhotoElecPt->Fill(pt);
 649     }
 650
 651     if (!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
 652
 653     if (m20>0.02 && m02>0.02 && m02<0.27 && fTPCnSigma>-1 && fTPCnSigma<3){
 654       fEoverPsig[iCent][iPt][iDeltaphi]->Fill(EovP);
 655       if (fFlagPhotonicElec) fEoverPuls[iCent][iPt][iDeltaphi]->Fill(EovP);
 656       if (fFlagPhotonicElecBCG) fEoverPls[iCent][iPt][iDeltaphi]->Fill(EovP);
 657     }
 658
 659   }//end of track loop
 660   PostData(1, fOutputList);
 661 }
 662 //_________________________________________
 663 void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
 664 {
 665   //--- Check MC
 666   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 667     fIsMC = kTRUE;
 668     printf("+++++ MC Data available");
 669   }
 670   //--------Initialize PID
 671   fPID->SetHasMCData(fIsMC);
 672
 673   if(!fPID->GetNumberOfPIDdetectors())
 674     {
 675       fPID->AddDetector("TPC", 0);
 676       fPID->AddDetector("EMCAL", 1);
 677     }
 678
 679   fPID->SortDetectors();
 680   fPIDqa = new AliHFEpidQAmanager();
 681   fPIDqa->Initialize(fPID);
 682
 683   //--------Initialize correction Framework and Cuts
 684   fCFM = new AliCFManager;
 685   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 686   fCFM->SetNStepParticle(kNcutSteps);
 687   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 688     fCFM->SetParticleCutsList(istep, NULL);
 689
 690   if(!fCuts){
 691     AliWarning("Cuts not available. Default cuts will be used");
 692     fCuts = new AliHFEcuts;
 693     fCuts->CreateStandardCuts();
 694   }
 695   fCuts->Initialize(fCFM);
 696
 697   //---------Output Tlist
 698   fOutputList = new TList();
 699   fOutputList->SetOwner();
 700   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 701
 702   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 703   fOutputList->Add(fNoEvents);
 704
 705   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 706   fOutputList->Add(fTrkpt);
 707
 708   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 709   fOutputList->Add(fTrackPtBefTrkCuts);
 710
 711   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 712   fOutputList->Add(fTrackPtAftTrkCuts);
 713
 714   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 715   fOutputList->Add(fTPCnsigma);
 716
 717   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 718   fOutputList->Add(fTrkEovPBef);
 719
 720   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 721   fOutputList->Add(fTrkEovPAft);
 722
 723   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
 724   fOutputList->Add(fdEdxBef);
 725
 726   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
 727   fOutputList->Add(fdEdxAft);
 728
 729   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 730   fOutputList->Add(fPhotoElecPt);
 731
 732   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
 733   fOutputList->Add(fSemiInclElecPt);
 734
 735   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 736   fOutputList->Add(fCent);
 737
 738   fTPCsubEPres = new TH1F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1);
 739   fOutputList->Add(fTPCsubEPres);
 740
 741   Int_t binsv1[4]={100,100,100,90}; // V0-TPCpos, V0-TPCneg, TPCpos-TPCneg, cent
 742   Double_t xminv1[4]={-1,-1,-1,0};
 743   Double_t xmaxv1[4]={1,1,1,90};
 744   fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
 745   fOutputList->Add(fEPres);
 746
 747   //iCent,iPt,fTPCnSigma,fEMCalnSigma,m02,dphi,cosdphi
 748   Int_t binsv2[7]={3,8,100,100,100,120,100};
 749   Double_t xminv2[7]={0,0,-5,-5,0,0,-1};
 750   Double_t xmaxv2[7]={3,8,5,5,2,TMath::TwoPi(),1};
 751   fCorr = new THnSparseD ("fCorr","Correlations",7,binsv2,xminv2,xmaxv2);
 752   fOutputList->Add(fCorr);
 753
 754   for(Int_t i=0; i<3; i++) {
 755     fevPlaneV0[i] = new TH1F(Form("fevPlaneV0%d",i),"V0 EP",100,0,TMath::Pi());
 756     fOutputList->Add(fevPlaneV0[i]);
 757
 758     feTPCV2[i] = new TH2F(Form("feTPCV2%d",i), "", 8,0,8,100,-1,1);
 759     fOutputList->Add(feTPCV2[i]);
 760
 761     feV2[i] = new TH2F(Form("feV2%d",i), "", 8,0,8,100,-1,1);
 762     fOutputList->Add(feV2[i]);
 763
 764     fChargPartV2[i] = new TH2F(Form("fChargPartV2%d",i), "", 8,0,8,100,-1,1);
 765     fOutputList->Add(fChargPartV2[i]);
 766
 767     fMtcPartV2[i] = new TH2F(Form("fMtcPartV2%d",i), "", 8,0,8,100,-1,1);
 768     fOutputList->Add(fMtcPartV2[i]);
 769
 770     fInvmassLS[i] = new TH2F(Form("fInvmassLS%d",i), "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5,100,0,50);
 771     fOutputList->Add(fInvmassLS[i]);
 772
 773     fInvmassULS[i] = new TH2F(Form("fInvmassULS%d",i), "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5,100,0,50);
 774     fOutputList->Add(fInvmassULS[i]);
 775
 776     fOpeningAngleLS[i] = new TH2F(Form("fOpeningAngleLS%d",i),"Opening angle for LS pairs",100,0,1,100,0,50);
 777     fOutputList->Add(fOpeningAngleLS[i]);
 778
 779     fOpeningAngleULS[i] = new TH2F(Form("fOpeningAngleULS%d",i),"Opening angle for ULS pairs",100,0,1,100,0,50);
 780     fOutputList->Add(fOpeningAngleULS[i]);
 781
 782     fPi0Pt[i] = new TH1F(Form("fPi0Pt%d",i), "Pi0 weight",100,0,50);
 783     fOutputList->Add(fPi0Pt[i]);
 784
 785     fEtaPt[i] = new TH1F(Form("fEtaPt%d",i), "Eta weight",100,0,50);
 786     fOutputList->Add(fEtaPt[i]);
 787   }
 788
 789   for(Int_t i=0; i<3; i++) {
 790     for(Int_t j=0; j<8; j++) {
 791       for(Int_t k=0; k<4; k++) {
 792         fEoverPsig[i][j][k] = new TH1F(Form("fEoverPsig%d%d%d",i,j,k), "",100,0,3);
 793         fOutputList->Add(fEoverPsig[i][j][k]);
 794
 795         fEoverPuls[i][j][k] = new TH1F(Form("fEoverPuls%d%d%d",i,j,k), "",100,0,3);
 796         fOutputList->Add(fEoverPuls[i][j][k]);
 797
 798         fEoverPls[i][j][k] = new TH1F(Form("fEoverPls%d%d%d",i,j,k), "",100,0,3);
 799         fOutputList->Add(fEoverPls[i][j][k]);
 800
 801         fEoverPbcg[i][j][k] = new TH1F(Form("fEoverPbcg%d%d%d",i,j,k), "",100,0,3);
 802         fOutputList->Add(fEoverPbcg[i][j][k]);
 803       }
 804     }
 805   }
 806
 807   fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
 808   fOutputList->Add(fD0_e);
 809
 810   for(Int_t k = 0; k < 6; k++) {
 811
 812     TString De_name = Form("fDe%d",k);
 813     TString D0e_name = Form("fD0e%d",k);
 814     TString Dpluse_name = Form("fDpluse%d",k);
 815     TString Dminuse_name = Form("fDminuse%d",k);
 816
 817     fDe[k] = new TH1F((const char*)De_name,"",100,0,50);
 818     fD0e[k] = new TH1F((const char*)D0e_name,"",100,0,50);
 819     fDpluse[k] = new TH1F((const char*)Dpluse_name,"",100,0,50);
 820     fDminuse[k] = new TH1F((const char*)Dminuse_name,"",100,0,50);
 821
 822     fOutputList->Add(fDe[k]);
 823     fOutputList->Add(fD0e[k]);
 824     fOutputList->Add(fDpluse[k]);
 825     fOutputList->Add(fDminuse[k]);
 826
 827   }
 828
 829   int nbin_v2 = 8;
 830   double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
 831
 832   fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
 833   fOutputList->Add(fTot_pi0e);
 834
 835   fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
 836   fOutputList->Add(fPhot_pi0e);
 837
 838   fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
 839   fOutputList->Add(fPhotBCG_pi0e);
 840
 841   fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
 842   fOutputList->Add(fTot_etae);
 843
 844   fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
 845   fOutputList->Add(fPhot_etae);
 846
 847   fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
 848   fOutputList->Add(fPhotBCG_etae);
 849
 850   fInvMass = new TH1F("fInvMass","",200,0,0.3);
 851   fOutputList->Add(fInvMass);
 852
 853   fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
 854   fOutputList->Add(fInvMassBack);
 855
 856   fDCA = new TH1F("fDCA","",200,0,1);
 857   fOutputList->Add(fDCA);
 858
 859   fDCABack = new TH1F("fDCABack","",200,0,1);
 860   fOutputList->Add(fDCABack);
 861
 862   fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
 863   fOutputList->Add(fOpAngle);
 864
 865   fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
 866   fOutputList->Add(fOpAngleBack);
 867
 868   PostData(1,fOutputList);
 869 }
 870
 871 //________________________________________________________________________
 872 void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
 873 {
 874   // Info("Terminate");
 875         AliAnalysisTaskSE::Terminate();
 876 }
 877
 878 //________________________________________________________________________
 879 Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 880 {
 881   // Check single track cuts for a given cut step
 882   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 883   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 884   return kTRUE;
 885 }
 886 //_________________________________________
 887 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
 888 {
 889   //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
 890   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 891   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 892   Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
 893
 894   return cos2DeltaPhi;
 895 }
 896 //_________________________________________
 897 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
 898 {
 899   //Get phi-psi_EP
 900   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 901   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 902
 903   return dPhi;
 904 }
 905 //_________________________________________
 906 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT, Float_t cent) const
 907 {
 908   //Get Pi0 weight
 909   double weight = 1.0;
 910   return weight;
 911 }
 912 //_________________________________________
 913 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT, Float_t cent) const
 914 {
 915   //Get eta weight
 916   double weight = 1.0;
 917   return weight;
 918 }
 919 //_________________________________________
 920 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetSigmaEMCal(Double_t EoverP, Double_t pt, Float_t cent) const
 921 {
 922   //Get sigma for EMCal PID
 923   Double_t NumberOfSigmasEMCal = 99.;
 924   Double_t ptRange[9] = {1.5,2,2.5,3,4,6,8,10,13};
 925
 926   if (cent>=0  && cent<10){
 927     Double_t mean[8]={0.953184,0.957259,0.97798,0.9875,1.03409,1.06257,1.02776,1.04338};
 928     Double_t sigma[8]={0.130003,0.113493,0.092966,0.0836828,0.101804,0.0893414,0.0950752,0.050427};
 929     for(Int_t i=0;i<8;i++) {
 930       if (pt>=ptRange[i] && pt<ptRange[i+1]){
 931         NumberOfSigmasEMCal = (mean[i]-EoverP)/sigma[i];
 932         continue;
 933       }
 934     }
 935   }
 936   if (cent>=10 && cent<20){
 937     Double_t mean[8]={0.96905,0.952985,0.96871,0.983934,1.00047,0.988736,1.02101,1.04557};
 938     Double_t sigma[8]={0.0978103,0.103215,0.0958494,0.0797962,0.0719482,0.0672677,0.0754882,0.0461192};
 939     for(Int_t i=0;i<8;i++) {
 940       if (pt>=ptRange[i] && pt<ptRange[i+1]){
 941         NumberOfSigmasEMCal = (mean[i]-EoverP)/sigma[i];
 942         continue;
 943       }
 944     }
 945   }
 946   if (cent>=20 && cent<40){
 947     Double_t mean[8]={0.947362,0.951933,0.959288,0.977004,0.984502,1.02004,1.00489,0.986696};
 948     Double_t sigma[8]={0.100127,0.0887731,0.0842077,0.0787335,0.0804325,0.0652376,0.0766669,0.0597849};
 949     for(Int_t i=0;i<8;i++) {
 950       if (pt>=ptRange[i] && pt<ptRange[i+1]){
 951         NumberOfSigmasEMCal = (mean[i]-EoverP)/sigma[i];
 952         continue;
 953       }
 954     }
 955   }
 956
 957
 958
 959   return NumberOfSigmasEMCal;
 960 }
 961 //_________________________________________
 962 Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsPi0EtaFromHFdecay(TParticle *particle, AliStack* stack)
 963 {
 964   // Check if the mother comes from heavy-flavour decays
 965
 966   Bool_t isHFdecay = kFALSE;
 967   Int_t partPDG = particle->GetPdgCode();
 968
 969   if (TMath::Abs(partPDG)!=111 || TMath::Abs(partPDG)!=221) return isHFdecay; // particle is not pi0 or eta
 970
 971   Int_t idMother = particle->GetFirstMother();
 972   if (idMother>0){
 973     TParticle *mother = stack->Particle(idMother);
 974     Int_t motherPDG = mother->GetPdgCode();
 975
 976     // c decays
 977     if((TMath::Abs(motherPDG)==411) || (TMath::Abs(motherPDG)==421) || (TMath::Abs(motherPDG)==431) || (TMath::Abs(motherPDG)==4122) || (TMath::Abs(motherPDG)==4132) || (TMath::Abs(motherPDG)==4232) || (TMath::Abs(motherPDG)==43320)) isHFdecay = kTRUE;
 978
 979     // b decays
 980     if((TMath::Abs(motherPDG)==511) || (TMath::Abs(motherPDG)==521) || (TMath::Abs(motherPDG)==531) || (TMath::Abs(motherPDG)==5122) || (TMath::Abs(motherPDG)==5132) || (TMath::Abs(motherPDG)==5232) || (TMath::Abs(motherPDG)==53320)) isHFdecay = kTRUE;
 981   }
 982
 983   return isHFdecay;
 984 }
 985 //_________________________________________
 986 Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsPi0EtaFromLMdecay(TParticle *particle, AliStack* stack)
 987 {
 988   // Check if the mother comes from light-meson decays
 989
 990   Bool_t isLMdecay = kFALSE;
 991   Int_t partPDG = particle->GetPdgCode();
 992
 993   if (TMath::Abs(partPDG)!=111 || TMath::Abs(partPDG)!=221) return isLMdecay; // particle is not pi0 or eta
 994
 995   Int_t idMother = particle->GetFirstMother();
 996   if (idMother>0){
 997     TParticle *mother = stack->Particle(idMother);
 998     Int_t motherPDG = mother->GetPdgCode();
 999
1000     if(motherPDG == 111 || motherPDG == 221 || motherPDG==223 || motherPDG==333 || motherPDG==331 || (TMath::Abs(motherPDG)==113) || (TMath::Abs(motherPDG)==213) || (TMath::Abs(motherPDG)==313) || (TMath::Abs(motherPDG)==323)) isLMdecay = kTRUE;
1001   }
1002
1003   return isLMdecay;
1004 }
1005 //_________________________________________
1006 Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsPi0EtaPrimary(TParticle *particle, AliStack* stack)
1007 {
1008   // Check if the pi0 or eta are primary
1009
1010   Bool_t isprimary = kFALSE;
1011   Int_t partPDG = particle->GetPdgCode();
1012
1013   if (TMath::Abs(partPDG)!=111 || TMath::Abs(partPDG)!=221) return isprimary; // particle is not pi0 or eta
1014
1015
1016   Bool_t pi0etaprimary = particle->IsPrimary();
1017   Int_t idMother = particle->GetFirstMother();
1018
1019   if (pi0etaprimary && idMother<=0) isprimary = kTRUE;
1020
1021   return isprimary;
1022 }
1023 //_________________________________________
1024 Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsElectronFromPi0(TParticle *particle, AliStack* stack, Double_t &weight, Float_t cent)
1025 {
1026   // Check if electron comes from primary pi0 not from light-meson and heavy-flavour decays
1027
1028   Bool_t isPi0Decay = kFALSE;
1029   Int_t partPDG = particle->GetPdgCode();
1030
1031   if (TMath::Abs(partPDG)!=11) return isPi0Decay; // particle is not electron
1032
1033   Int_t idMother = particle->GetFirstMother();
1034   if (idMother>0){
1035     TParticle *mother = stack->Particle(idMother);
1036     Int_t motherPDG = mother->GetPdgCode();
1037     Double_t motherPt = mother->Pt();
1038
1039     Bool_t isMotherPi0primary = IsPi0EtaPrimary(mother,stack);
1040     Bool_t isMotherPi0fromHF = IsPi0EtaFromHFdecay(mother,stack);
1041     Bool_t isMotherPi0fromLM = IsPi0EtaFromLMdecay(mother,stack);
1042
1043     if (motherPDG==111 && (isMotherPi0primary || (!isMotherPi0fromHF && !isMotherPi0fromLM))){ // pi0 -> e
1044       isPi0Decay = kTRUE;
1045       weight = GetPi0weight(motherPt,cent);
1046     }
1047
1048     Int_t idSecondMother = particle->GetSecondMother();
1049     if (motherPDG==22 && idSecondMother>0){
1050       TParticle *secondMother = stack->Particle(idSecondMother);
1051       Int_t secondMotherPDG = secondMother->GetPdgCode();
1052       Double_t secondMotherPt = secondMother->Pt();
1053
1054       Bool_t isSecondMotherPi0primary = IsPi0EtaPrimary(secondMother,stack);
1055       Bool_t isSecondMotherPi0fromHF = IsPi0EtaFromHFdecay(secondMother,stack);
1056       Bool_t isSecondMotherPi0fromLM = IsPi0EtaFromLMdecay(secondMother,stack);
1057
1058       if (secondMotherPDG==111 && (isSecondMotherPi0primary || (!isSecondMotherPi0fromHF && !isSecondMotherPi0fromLM))){ //pi0 -> gamma -> e
1059         isPi0Decay = kTRUE;
1060         weight = GetPi0weight(secondMotherPt,cent);
1061       }
1062     }
1063   }
1064   return isPi0Decay;
1065 }
1066 //_________________________________________
1067 Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsElectronFromEta(TParticle *particle, AliStack* stack, Double_t &weight, Float_t cent)
1068 {
1069   // Check if electron comes from primary eta not from light-meson and heavy-flavour decays
1070
1071   Bool_t isEtaDecay = kFALSE;
1072   Int_t partPDG = particle->GetPdgCode();
1073
1074   if (TMath::Abs(partPDG)!=11) return isEtaDecay; // particle is not electron
1075
1076   Int_t idMother = particle->GetFirstMother();
1077   if (idMother>0){
1078     TParticle *mother = stack->Particle(idMother);
1079     Int_t motherPDG = mother->GetPdgCode();
1080     Double_t motherPt = mother->Pt();
1081
1082     Bool_t isMotherEtaprimary = IsPi0EtaPrimary(mother,stack);
1083     Bool_t isMotherEtafromHF = IsPi0EtaFromHFdecay(mother,stack);
1084     Bool_t isMotherEtafromLM = IsPi0EtaFromLMdecay(mother,stack);
1085
1086     if (motherPDG==221  && (isMotherEtaprimary || (!isMotherEtafromHF && !isMotherEtafromLM))){ //primary eta -> e
1087       isEtaDecay = kTRUE;
1088       weight = GetEtaweight(motherPt,cent);
1089     }
1090
1091     Int_t idSecondMother = mother->GetFirstMother();
1092     if ((motherPDG==22 || motherPDG==111) && idSecondMother>0){
1093       TParticle *secondMother = stack->Particle(idSecondMother);
1094       Int_t secondMotherPDG = secondMother->GetPdgCode();
1095       Double_t secondMotherPt = secondMother->Pt();
1096
1097       Bool_t isSecondMotherEtaprimary = IsPi0EtaPrimary(secondMother,stack);
1098       Bool_t isSecondMotherEtafromHF = IsPi0EtaFromHFdecay(secondMother,stack);
1099       Bool_t isSecondMotherEtafromLM = IsPi0EtaFromLMdecay(secondMother,stack);
1100
1101       if (secondMotherPDG==221  && (isSecondMotherEtaprimary || (!isSecondMotherEtafromHF && !isSecondMotherEtafromLM))){ //eta -> pi0/g-> e
1102         isEtaDecay = kTRUE;
1103         weight = GetEtaweight(secondMotherPt,cent);
1104       }
1105       Int_t idThirdMother = secondMother->GetFirstMother();
1106       if (idThirdMother>0){
1107         TParticle *thirdMother = stack->Particle(idThirdMother);
1108         Int_t thirdMotherPDG = thirdMother->GetPdgCode();
1109         Double_t thirdMotherPt = thirdMother->Pt();
1110
1111         Bool_t isThirdMotherEtaprimary = IsPi0EtaPrimary(thirdMother,stack);
1112         Bool_t isThirdMotherEtafromHF = IsPi0EtaFromHFdecay(thirdMother,stack);
1113         Bool_t isThirdMotherEtafromLM = IsPi0EtaFromLMdecay(thirdMother,stack);
1114
1115         if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221 && (isThirdMotherEtaprimary || (!isThirdMotherEtafromHF && !isThirdMotherEtafromLM))){//p eta->pi0->g-> e
1116           isEtaDecay = kTRUE;
1117           weight = GetEtaweight(thirdMotherPt,cent);
1118         }
1119       }
1120     }
1121   }
1122   return isEtaDecay;
1123 }
1124 //_________________________________________
1125 void AliAnalysisTaskFlowTPCEMCalEP::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG,Double_t weight, Int_t iCent)
1126 {
1127   //Identify non-heavy flavour electrons using Invariant mass method
1128
1129   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1130   fTrackCuts->SetRequireTPCRefit(kTRUE);
1131   fTrackCuts->SetRequireITSRefit(kTRUE);
1132   fTrackCuts->SetEtaRange(-0.9,0.9);
1133   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1134   fTrackCuts->SetMaxChi2PerClusterTPC(4);
1135   fTrackCuts->SetMinNClustersTPC(80);
1136   fTrackCuts->SetMaxDCAToVertexZ(3.2);
1137   fTrackCuts->SetMaxDCAToVertexXY(2.4);
1138   fTrackCuts->SetDCAToVertex2D(kTRUE);
1139
1140   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1141
1142   Bool_t flagPhotonicElec = kFALSE;
1143   Bool_t flagPhotonicElecBCG = kFALSE;
1144
1145   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1146
1147     if(jTracks==iTracks) continue;
1148
1149     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1150     if (!trackAsso) {
1151       printf("ERROR: Could not receive track %d\n", jTracks);
1152       continue;
1153     }
1154
1155     Double_t pt=-999., ptAsso=-999., nTPCsigmaAsso=-999.;
1156     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1157     Double_t openingAngle = -999., mass=999., width = -999;
1158     Int_t chargeAsso = 0, charge = 0;
1159
1160     nTPCsigmaAsso = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1161     pt = track->Pt();
1162     ptAsso = trackAsso->Pt();
1163     chargeAsso = trackAsso->Charge();
1164     charge = track->Charge();
1165
1166     if(ptAsso <0.5) continue;
1167     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1168     if(TMath::Abs(nTPCsigmaAsso)>3) continue;
1169
1170     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1171     if(charge>0) fPDGe1 = -11;
1172     if(chargeAsso>0) fPDGe2 = -11;
1173
1174     if(charge == chargeAsso) fFlagLS = kTRUE;
1175     if(charge != chargeAsso) fFlagULS = kTRUE;
1176
1177     AliKFParticle ge1(*track, fPDGe1);
1178     AliKFParticle ge2(*trackAsso, fPDGe2);
1179     AliKFParticle recg(ge1, ge2);
1180
1181     if(recg.GetNDF()<1) continue;
1182     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1183     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1184
1185     if(fSetMassConstraint && pVtx) {
1186       AliKFVertex primV(*pVtx);
1187       primV += recg;
1188       primV -= ge1;
1189       primV -= ge2;
1190       recg.SetProductionVertex(primV);
1191       recg.SetMassConstraint(0,0.0001);
1192     }
1193
1194     openingAngle = ge1.GetAngle(ge2);
1195
1196     if(fFlagLS) fOpeningAngleLS[iCent]->Fill(openingAngle,pt);
1197     if(fFlagULS) fOpeningAngleULS[iCent]->Fill(openingAngle,pt);
1198
1199     if(openingAngle > fOpeningAngleCut) continue;
1200
1201     recg.GetMass(mass,width);
1202
1203     if(fFlagLS) fInvmassLS[iCent]->Fill(mass,pt,weight);
1204     if(fFlagULS) fInvmassULS[iCent]->Fill(mass,pt,weight);
1205
1206     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
1207     if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
1208
1209   }
1210   fFlagPhotonicElec = flagPhotonicElec;
1211   fFlagPhotonicElecBCG = flagPhotonicElecBCG;
1212
1213 }
1214 //_________________________________________
1215 void AliAnalysisTaskFlowTPCEMCalEP::InitParameters()
1216 {
1217   // Init parameters
1218
1219   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1220   fTrackCuts->SetRequireTPCRefit(kTRUE);
1221   fTrackCuts->SetRequireITSRefit(kTRUE);
1222   fTrackCuts->SetEtaRange(-0.7,0.7);
1223   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1224   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1225   fTrackCuts->SetMinNClustersTPC(100);
1226   fTrackCuts->SetPtRange(0.5,100);
1227
1228   fNonHFE->SetAODanalysis(kFALSE);
1229   fNonHFE->SetInvariantMassCut(fInvmassCut);
1230   fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
1231   fNonHFE->SetChi2OverNDFCut(fChi2Cut);
1232   fNonHFE->SetAlgorithm(fnonHFEalgorithm); //KF or DCA
1233   if (fnonHFEalgorithm == "DCA") fNonHFE->SetDCACut(fDCAcut);
1234   fNonHFE->SetTrackCuts(-3.5,3.5,fTrackCuts);
1235
1236 }
1237