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
  38 #include "AliAnalysisTaskFlowTPCEMCalEP.h"
  39 #include "TGeoGlobalMagField.h"
  40 #include "AliLog.h"
  41 #include "AliAnalysisTaskSE.h"
  42 #include "TRefArray.h"
  43 #include "TVector.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliESDpid.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliPhysicsSelection.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "AliAODCaloCluster.h"
  50 #include "AliEMCALRecoUtils.h"
  51 #include "AliEMCALGeometry.h"
  52 #include "AliGeomManager.h"
  53 #include "stdio.h"
  54 #include "TGeoManager.h"
  55 #include "iostream"
  56 #include "fstream"
  57
  58 #include "AliEMCALTrack.h"
  59 #include "AliMagF.h"
  60
  61 #include "AliKFParticle.h"
  62 #include "AliKFVertex.h"
  63
  64 #include "AliMCEventHandler.h"
  65 #include "AliMCEvent.h"
  66 #include "AliMCParticle.h"
  67 #include "AliStack.h"
  68
  69 #include "AliPID.h"
  70 #include "AliPIDResponse.h"
  71 #include "AliHFEcontainer.h"
  72 #include "AliHFEcuts.h"
  73 #include "AliHFEpid.h"
  74 #include "AliHFEpidBase.h"
  75 #include "AliHFEpidQAmanager.h"
  76 #include "AliHFEtools.h"
  77 #include "AliCFContainer.h"
  78 #include "AliCFManager.h"
  79
  80 #include "AliEventplane.h"
  81 #include "AliCentrality.h"
  82
  83 ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
  84 //________________________________________________________________________
  85 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
  86   : AliAnalysisTaskSE(name)
  87   ,fESD(0)
  88   ,fMC(0)
  89   ,fOutputList(0)
  90   ,fTrackCuts(0)
  91   ,fCuts(0)
  92   ,fIdentifiedAsOutInz(kFALSE)
  93   ,fPassTheEventCut(kFALSE)
  94   ,fRejectKinkMother(kFALSE)
  95   ,fIsMC(kFALSE)
  96   ,fVz(0.0)
  97   ,fCFM(0)
  98   ,fPID(0)
  99   ,fPIDqa(0)
 100   ,fOpeningAngleCut(0.1)
 101   ,fInvmassCut(0.01)
 102   ,fNoEvents(0)
 103   ,fTrkpt(0)
 104   ,fTrkEovPBef(0)
 105   ,fTrkEovPAft(0)
 106   ,fdEdxBef(0)
 107   ,fdEdxAft(0)
 108   ,fInvmassLS(0)
 109   ,fInvmassULS(0)
 110   ,fOpeningAngleLS(0)
 111   ,fOpeningAngleULS(0)
 112   ,fPhotoElecPt(0)
 113   ,fSemiInclElecPt(0)
 114   ,fMCphotoElecPt(0)
 115   ,fTrackPtBefTrkCuts(0)
 116   ,fTrackPtAftTrkCuts(0)
 117   ,fTPCnsigma(0)
 118   ,fCent(0)
 119   ,fevPlaneV0A(0)
 120   ,fevPlaneV0C(0)
 121   ,fevPlaneTPC(0)
 122   ,fTPCsubEPres(0)
 123   ,fEPres(0)
 124   ,fCorr(0)
 125   ,feTPCV2(0)
 126   ,feV2(0)
 127   ,fphoteV2(0)
 128   ,fChargPartV2(0)
 129   ,fGammaWeight(0)
 130   ,fPi0Weight(0)
 131   ,fEtaWeight(0)
 132   ,fD0Weight(0)
 133   ,fDplusWeight(0)
 134   ,fDminusWeight(0)
 135   ,fD0_e(0)
 136   ,fTot_pi0e(0)
 137   ,fPhot_pi0e(0)
 138   ,fPhotBCG_pi0e(0)
 139   ,fTot_etae(0)
 140   ,fPhot_etae(0)
 141   ,fPhotBCG_etae(0)
 142   ,fCocktail(0)
 143 {
 144   //Named constructor
 145
 146   for(Int_t k = 0; k < 6; k++) {
 147     fDe[k]= NULL;
 148     fD0e[k]= NULL;
 149     fDpluse[k]= NULL;
 150     fDminuse[k]= NULL;
 151   }
 152
 153   fPID = new AliHFEpid("hfePid");
 154   fTrackCuts = new AliESDtrackCuts();
 155
 156   // Define input and output slots here
 157   // Input slot #0 works with a TChain
 158   DefineInput(0, TChain::Class());
 159   // Output slot #0 id reserved by the base class for AOD
 160   // Output slot #1 writes into a TH1 container
 161   // DefineOutput(1, TH1I::Class());
 162   DefineOutput(1, TList::Class());
 163   //  DefineOutput(3, TTree::Class());
 164 }
 165
 166 //________________________________________________________________________
 167 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
 168   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
 169   ,fESD(0)
 170   ,fMC(0)
 171   ,fOutputList(0)
 172   ,fTrackCuts(0)
 173   ,fCuts(0)
 174   ,fIdentifiedAsOutInz(kFALSE)
 175   ,fPassTheEventCut(kFALSE)
 176   ,fRejectKinkMother(kFALSE)
 177   ,fIsMC(kFALSE)
 178   ,fVz(0.0)
 179   ,fCFM(0)
 180   ,fPID(0)
 181   ,fPIDqa(0)
 182   ,fOpeningAngleCut(0.1)
 183   ,fInvmassCut(0.01)
 184   ,fNoEvents(0)
 185   ,fTrkpt(0)
 186   ,fTrkEovPBef(0)
 187   ,fTrkEovPAft(0)
 188   ,fdEdxBef(0)
 189   ,fdEdxAft(0)
 190   ,fInvmassLS(0)
 191   ,fInvmassULS(0)
 192   ,fOpeningAngleLS(0)
 193   ,fOpeningAngleULS(0)
 194   ,fPhotoElecPt(0)
 195   ,fSemiInclElecPt(0)
 196   ,fMCphotoElecPt(0)
 197   ,fTrackPtBefTrkCuts(0)
 198   ,fTrackPtAftTrkCuts(0)
 199   ,fTPCnsigma(0)
 200   ,fCent(0)
 201   ,fevPlaneV0A(0)
 202   ,fevPlaneV0C(0)
 203   ,fevPlaneTPC(0)
 204   ,fTPCsubEPres(0)
 205   ,fEPres(0)
 206   ,fCorr(0)
 207   ,feTPCV2(0)
 208   ,feV2(0)
 209   ,fphoteV2(0)
 210   ,fChargPartV2(0)
 211   ,fGammaWeight(0)
 212   ,fPi0Weight(0)
 213   ,fEtaWeight(0)
 214   ,fD0Weight(0)
 215   ,fDplusWeight(0)
 216   ,fDminusWeight(0)
 217   ,fD0_e(0)
 218   ,fTot_pi0e(0)
 219   ,fPhot_pi0e(0)
 220   ,fPhotBCG_pi0e(0)
 221   ,fTot_etae(0)
 222   ,fPhot_etae(0)
 223   ,fPhotBCG_etae(0)
 224   ,fCocktail(0)
 225 {
 226         //Default constructor
 227
 228         for(Int_t k = 0; k < 6; k++) {
 229           fDe[k]= NULL;
 230           fD0e[k]= NULL;
 231           fDpluse[k]= NULL;
 232           fDminuse[k]= NULL;
 233         }
 234
 235         fPID = new AliHFEpid("hfePid");
 236
 237         fTrackCuts = new AliESDtrackCuts();
 238
 239         // Constructor
 240         // Define input and output slots here
 241         // Input slot #0 works with a TChain
 242         DefineInput(0, TChain::Class());
 243         // Output slot #0 id reserved by the base class for AOD
 244         // Output slot #1 writes into a TH1 container
 245         // DefineOutput(1, TH1I::Class());
 246         DefineOutput(1, TList::Class());
 247         //DefineOutput(3, TTree::Class());
 248 }
 249 //_________________________________________
 250
 251 AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
 252 {
 253   //Destructor
 254
 255   delete fOutputList;
 256   delete fPID;
 257   delete fCFM;
 258   delete fPIDqa;
 259   delete fTrackCuts;
 260 }
 261 //_________________________________________
 262
 263 void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
 264 {
 265   //Main loop
 266   //Called for each event
 267
 268   // create pointer to event
 269   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 270   if (!fESD) {
 271     printf("ERROR: fESD not available\n");
 272     return;
 273   }
 274
 275   if(!fCuts){
 276     AliError("HFE cuts not available");
 277     return;
 278   }
 279
 280   if(!fPID->IsInitialized()){
 281     // Initialize PID with the given run number
 282     AliWarning("PID not initialised, get from Run no");
 283     fPID->InitializePID(fESD->GetRunNumber());
 284   }
 285
 286   if(fIsMC)fMC = MCEvent();
 287   AliStack* stack = NULL;
 288   if(fIsMC && fMC) stack = fMC->Stack();
 289
 290   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 291   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 292
 293   Double_t pVtxZ = -999;
 294   pVtxZ = pVtx->GetZ();
 295
 296   if(TMath::Abs(pVtxZ)>10) return;
 297   fNoEvents->Fill(0);
 298
 299   if(fNOtrks<2) return;
 300
 301   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 302   if(!pidResponse){
 303     AliDebug(1, "Using default PID Response");
 304     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 305   }
 306
 307   fPID->SetPIDResponse(pidResponse);
 308
 309   fCFM->SetRecEventInfo(fESD);
 310
 311   Float_t cent = -1.;
 312   AliCentrality *centrality = fESD->GetCentrality();
 313   cent = centrality->GetCentralityPercentile("V0M");
 314   fCent->Fill(cent);
 315
 316   int kCent = 3;
 317   if(cent>=20 && cent<=40) kCent = 0;
 318   if(cent>=30 && cent<=50) kCent = 1;
 319
 320   //Event planes
 321
 322   Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
 323   if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
 324   fevPlaneV0A->Fill(evPlaneV0A,cent);
 325
 326   Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
 327   if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
 328   fevPlaneV0C->Fill(evPlaneV0C,cent);
 329
 330   AliEventplane* esdTPCep = fESD->GetEventplane();
 331   TVector2 *standardQ = 0x0;
 332   Double_t qx = -999., qy = -999.;
 333   standardQ = esdTPCep->GetQVector();
 334   if(!standardQ)return;
 335
 336   qx = standardQ->X();
 337   qy = standardQ->Y();
 338
 339   TVector2 qVectorfortrack;
 340   qVectorfortrack.Set(qx,qy);
 341   Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
 342   fevPlaneTPC->Fill(evPlaneTPC,cent);
 343
 344   TVector2 *qsub1a = esdTPCep->GetQsub1();
 345   TVector2 *qsub2a = esdTPCep->GetQsub2();
 346   Double_t evPlaneResTPC = -999.;
 347   if(qsub1a && qsub2a)
 348   {
 349           evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
 350   }
 351
 352   fTPCsubEPres->Fill(evPlaneResTPC,cent);
 353
 354   Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
 355   fEPres->Fill(evPlaneRes);
 356
 357   if(fIsMC && fMC && stack && cent>20 && cent<40){
 358    Int_t nParticles = stack->GetNtrack();
 359    for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 360       TParticle* particle = stack->Particle(iParticle);
 361       int fPDG = particle->GetPdgCode();
 362       double pTMC = particle->Pt();
 363       double etaMC = particle->Eta();
 364       if(fabs(etaMC)>0.7)continue;
 365
 366       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
 367       int iHijing = 1;
 368       if(!MChijing)iHijing = 0;
 369
 370       if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
 371       if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
 372       if(fPDG==421)fD0Weight->Fill(pTMC,iHijing);//D0
 373       if(fPDG==411)fDplusWeight->Fill(pTMC,iHijing);//D+
 374       if(fPDG==-411)fDminusWeight->Fill(pTMC,iHijing);//D-
 375
 376
 377       Int_t idMother = particle->GetFirstMother();
 378       if (idMother>0){
 379         TParticle *mother = stack->Particle(idMother);
 380         int motherPDG = mother->GetPdgCode();
 381         if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
 382       }
 383
 384    }
 385   }
 386
 387
 388   // Track loop
 389   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 390     AliESDtrack* track = fESD->GetTrack(iTracks);
 391     if (!track) {
 392       printf("ERROR: Could not receive track %d\n", iTracks);
 393       continue;
 394     }
 395
 396     if(TMath::Abs(track->Eta())>0.7) continue;
 397
 398     fTrackPtBefTrkCuts->Fill(track->Pt());
 399
 400     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 401
 402     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 403       if(track->GetKinkIndex(0) != 0) continue;
 404     }
 405
 406     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 407
 408     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 409
 410     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 411
 412     fTrackPtAftTrkCuts->Fill(track->Pt());
 413
 414     Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999., wclsE = -999., wEovP = -999., m02= -999., m20= -999.;
 415
 416     pt = track->Pt();
 417     if(pt<2) continue;
 418     fTrkpt->Fill(pt);
 419
 420     Int_t clsId = track->GetEMCALcluster();
 421     if (clsId>0){
 422       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
 423       if(cluster && cluster->IsEMCAL()){
 424         clsE = cluster->E();
 425         m20 = cluster->GetM20();
 426         m02 = cluster->GetM02();
 427       }
 428     }
 429
 430     p = track->P();
 431     phi = track->Phi();
 432     dEdx = track->GetTPCsignal();
 433     EovP = clsE/p;
 434     wEovP = wclsE/p;
 435     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 436     fdEdxBef->Fill(p,dEdx);
 437     fTPCnsigma->Fill(p,fTPCnSigma);
 438
 439     //Remove electron candidate from the event plane
 440     Float_t evPlaneCorrTPC = evPlaneTPC;
 441     if(dEdx>70 && dEdx<90){
 442       Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
 443       Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
 444       TVector2 newQVectorfortrack;
 445       newQVectorfortrack.Set(qX,qY);
 446       evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
 447     }
 448
 449     Bool_t fFlagPhotonicElec = kFALSE;
 450     Bool_t fFlagPhotonicElecBCG = kFALSE;
 451
 452     SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG);
 453
 454     Double_t corr[12]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneCorrTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C), static_cast<Double_t>(fFlagPhotonicElec), static_cast<Double_t>(fFlagPhotonicElecBCG),m02,m20};
 455     fCorr->Fill(corr);
 456
 457
 458     Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
 459     Int_t whichPart = -99;
 460     Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
 461     Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
 462     Double_t weight = 1.;
 463     Double_t Dweight = 1.;
 464     Bool_t MChijing;
 465
 466     Bool_t pi0Decay= kFALSE;
 467     Bool_t etaDecay= kFALSE;
 468
 469
 470     Double_t phiD = -999.,phie = -999.,phieRec = -999.,ptD = -999.,pte = -999.,pteRec = -999.;
 471
 472     if(fIsMC && fMC && stack){
 473       Int_t label = track->GetLabel();
 474       if(label>0){
 475         TParticle *particle = stack->Particle(label);
 476         if(particle){
 477             partPDG = particle->GetPdgCode();
 478             partPt = particle->Pt();
 479
 480             if (TMath::Abs(partPDG)==11) whichPart = 0; //electron
 481             if (partPDG==22) whichPart = 3; //gamma
 482             if (partPDG==111) whichPart = 2; //pi0
 483             if (partPDG==221) whichPart = 1; //eta
 484
 485             MChijing = fMC->IsFromBGEvent(label);
 486
 487             int iHijing = 1;
 488             if(!MChijing) iHijing = 0; // 0 if enhanced sample
 489
 490             Int_t idMother = particle->GetFirstMother();
 491             if (idMother>0){
 492               TParticle *mother = stack->Particle(idMother);
 493               motherPt = mother->Pt();
 494               motherPDG = mother->GetPdgCode();
 495
 496              //*************** for cocktail ********************************
 497              if (TMath::Abs(partPDG)==11 && (motherPDG==421 || TMath::Abs(motherPDG)==411)){
 498                 pteRec = track->Pt();
 499                 phieRec = track->Phi();
 500                 phie = particle->Phi();
 501                 pte = particle->Pt();
 502                 phiD = mother->Phi();
 503                 ptD = mother->Pt();
 504
 505                 Double_t cocktail[9]={phiD,phie,phieRec,ptD,pte,pteRec,evPlaneV0A, static_cast<Double_t>(iHijing), static_cast<Double_t>(kCent)};
 506                 fCocktail->Fill(cocktail);
 507              }
 508              //*************************************************************
 509               if (motherPDG==22) whichFirstMother = 3; //gamma
 510               if (motherPDG==111) whichFirstMother = 2; //pi0
 511               if (motherPDG==221) whichFirstMother = 1; //eta
 512
 513               Int_t idSecondMother = particle->GetSecondMother();
 514               if (idSecondMother>0){
 515                 TParticle *secondMother = stack->Particle(idSecondMother);
 516                 secondMotherPt = secondMother->Pt();
 517                 secondMotherPDG = secondMother->GetPdgCode();
 518
 519                 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
 520                 if (secondMotherPDG==221) whichSecondMother = 1; //eta
 521
 522                 Int_t idThirdMother = secondMother->GetFirstMother();
 523                 if (idThirdMother>0){
 524                   TParticle *thirdMother = stack->Particle(idThirdMother);
 525                   thirdMotherPt = thirdMother->Pt();
 526                   thirdMotherPDG = thirdMother->GetPdgCode();
 527
 528                   if (thirdMotherPDG==221) whichThirdMother = 1; //eta
 529                 }
 530               }
 531
 532
 533               if (cent>30 && cent<50 && TMath::Abs(partPDG)==11){
 534
 535                 if (motherPDG==421 || TMath::Abs(motherPDG)==411){ // D
 536                   Double_t phi_D = -99., Deltaphi_De=-99.;
 537                   phi_D = mother->Phi();
 538                   Deltaphi_De = phi_D - phi;
 539
 540                   fD0_e->Fill(pt,Deltaphi_De);
 541
 542                   Dweight  = GetDweight(0,motherPt);
 543                   if(iHijing==1) Dweight = 1.0;
 544
 545                   if (motherPt>=2 && motherPt<3) fDe[0]->Fill(pt,Dweight);
 546                   if (motherPt>=3 && motherPt<4) fDe[1]->Fill(pt,Dweight);
 547                   if (motherPt>=4 && motherPt<6) fDe[2]->Fill(pt,Dweight);
 548                   if (motherPt>=6 && motherPt<8) fDe[3]->Fill(pt,Dweight);
 549                   if (motherPt>=8 && motherPt<12) fDe[4]->Fill(pt,Dweight);
 550                   if (motherPt>=12 && motherPt<16) fDe[5]->Fill(pt,Dweight);
 551                 }
 552                 if (motherPDG==421){ // D0
 553
 554                   Dweight  = GetDweight(1,motherPt);
 555                   if(iHijing==1) Dweight = 1.0;
 556
 557                   if (motherPt>=2 && motherPt<3) fD0e[0]->Fill(pt,Dweight);
 558                   if (motherPt>=3 && motherPt<4) fD0e[1]->Fill(pt,Dweight);
 559                   if (motherPt>=4 && motherPt<6) fD0e[2]->Fill(pt,Dweight);
 560                   if (motherPt>=6 && motherPt<8) fD0e[3]->Fill(pt,Dweight);
 561                   if (motherPt>=8 && motherPt<12) fD0e[4]->Fill(pt,Dweight);
 562                   if (motherPt>=12 && motherPt<16) fD0e[5]->Fill(pt,Dweight);
 563                 }
 564                 if (motherPDG==411){ // D+
 565                   Dweight  = GetDweight(2,motherPt);
 566                   if(iHijing==1) Dweight = 1.0;
 567
 568                   if (motherPt>=2 && motherPt<3) fDpluse[0]->Fill(pt,Dweight);
 569                   if (motherPt>=3 && motherPt<4) fDpluse[1]->Fill(pt,Dweight);
 570                   if (motherPt>=4 && motherPt<6) fDpluse[2]->Fill(pt,Dweight);
 571                   if (motherPt>=6 && motherPt<8) fDpluse[3]->Fill(pt,Dweight);
 572                   if (motherPt>=8 && motherPt<12) fDpluse[4]->Fill(pt,Dweight);
 573                   if (motherPt>=12 && motherPt<16) fDpluse[5]->Fill(pt,Dweight);
 574                 }
 575                 if (motherPDG==-411){ //D-
 576                   Dweight  = GetDweight(3,motherPt);
 577                   if(iHijing==1) Dweight = 1.0;
 578
 579                   if (motherPt>=2 && motherPt<3) fDminuse[0]->Fill(pt,Dweight);
 580                   if (motherPt>=3 && motherPt<4) fDminuse[1]->Fill(pt,Dweight);
 581                   if (motherPt>=4 && motherPt<6) fDminuse[2]->Fill(pt,Dweight);
 582                   if (motherPt>=6 && motherPt<8) fDminuse[3]->Fill(pt,Dweight);
 583                   if (motherPt>=8 && motherPt<12) fDminuse[4]->Fill(pt,Dweight);
 584                   if (motherPt>=12 && motherPt<16) fDminuse[5]->Fill(pt,Dweight);
 585                 }
 586               }
 587
 588               //pi0 decay
 589               if (TMath::Abs(partPDG)==11 && motherPDG==111 && secondMotherPDG!=221){ //not eta -> pi0 -> e
 590                 weight = GetPi0weight(motherPt);
 591                 pi0Decay = kTRUE;
 592               }
 593               if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG!=221){ //not eta -> pi0 -> gamma -> e
 594                 weight = GetPi0weight(secondMotherPt);
 595                 pi0Decay = kTRUE;
 596               }
 597
 598               //eta decay
 599               if (TMath::Abs(partPDG)==11 && motherPDG==221){ //eta -> e
 600                 weight = GetEtaweight(motherPt);
 601                 etaDecay = kTRUE;
 602               }
 603               if (TMath::Abs(partPDG)==11 && motherPDG==111 && secondMotherPDG==221){ //eta -> pi0 -> e
 604                 weight = GetEtaweight(secondMotherPt);
 605                 etaDecay = kTRUE;
 606               }
 607               if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==221){ //eta -> gamma -> e
 608                 weight = GetEtaweight(secondMotherPt);
 609                 etaDecay = kTRUE;
 610               }
 611               if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221){ //eta -> pi0 -> gamma -> e
 612                 weight = GetEtaweight(thirdMotherPt);
 613                 etaDecay = kTRUE;
 614               }
 615
 616               if (cent>30 && cent<50 && fTPCnSigma>-1 && fTPCnSigma<3 && EovP>1 && EovP<1.3 && (motherPDG==22 || motherPDG==111 || motherPDG==221)){
 617                   if(iHijing==1) weight = 1.0;
 618
 619                   if (pi0Decay){
 620                       fTot_pi0e->Fill(partPt,weight);
 621                       if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
 622                       if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
 623                   }
 624
 625                   if (etaDecay){
 626                       fTot_etae->Fill(partPt,weight);
 627                       if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
 628                       if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
 629                   }
 630
 631               }
 632
 633
 634               Double_t mc[15]={EovP,fTPCnSigma,partPt, static_cast<Double_t>(fFlagPhotonicElec), static_cast<Double_t>(fFlagPhotonicElecBCG), static_cast<Double_t>(whichPart),cent,pt, static_cast<Double_t>(whichFirstMother), static_cast<Double_t>(whichSecondMother), static_cast<Double_t>(whichThirdMother), static_cast<Double_t>(iHijing),motherPt,secondMotherPt,thirdMotherPt};
 635
 636 //            fMCphotoElecPt->Fill(mc);
 637               if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
 638           }
 639         }
 640       }
 641     }
 642
 643
 644     if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
 645     Int_t pidpassed = 1;
 646
 647     //--- track accepted
 648     AliHFEpidObject hfetrack;
 649     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 650     hfetrack.SetRecTrack(track);
 651     hfetrack.SetPbPb();
 652     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 653
 654     Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 655     fChargPartV2->Fill(corrV2);
 656
 657     if(fTPCnSigma >= -0.5){
 658         Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 659         feTPCV2->Fill(correctedV2);
 660     }
 661
 662     if(pidpassed==0) continue;
 663
 664     Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 665
 666     feV2->Fill(correctedV2);
 667
 668     fTrkEovPAft->Fill(pt,EovP);
 669     fdEdxAft->Fill(p,dEdx);
 670
 671   if(fFlagPhotonicElec){
 672     fphoteV2->Fill(correctedV2);
 673     fPhotoElecPt->Fill(pt);
 674   }
 675
 676   if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
 677
 678  }
 679  PostData(1, fOutputList);
 680 }
 681 //_________________________________________
 682 void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
 683 {
 684   //--- Check MC
 685   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 686     fIsMC = kTRUE;
 687     printf("+++++ MC Data available");
 688   }
 689   //--------Initialize PID
 690   fPID->SetHasMCData(fIsMC);
 691
 692   if(!fPID->GetNumberOfPIDdetectors())
 693     {
 694       fPID->AddDetector("TPC", 0);
 695       fPID->AddDetector("EMCAL", 1);
 696     }
 697
 698   fPID->SortDetectors();
 699   fPIDqa = new AliHFEpidQAmanager();
 700   fPIDqa->Initialize(fPID);
 701
 702   //--------Initialize correction Framework and Cuts
 703   fCFM = new AliCFManager;
 704   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 705   fCFM->SetNStepParticle(kNcutSteps);
 706   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 707     fCFM->SetParticleCutsList(istep, NULL);
 708
 709   if(!fCuts){
 710     AliWarning("Cuts not available. Default cuts will be used");
 711     fCuts = new AliHFEcuts;
 712     fCuts->CreateStandardCuts();
 713   }
 714   fCuts->Initialize(fCFM);
 715
 716   //---------Output Tlist
 717   fOutputList = new TList();
 718   fOutputList->SetOwner();
 719   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 720
 721   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 722   fOutputList->Add(fNoEvents);
 723
 724   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 725   fOutputList->Add(fTrkpt);
 726
 727   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 728   fOutputList->Add(fTrackPtBefTrkCuts);
 729
 730   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 731   fOutputList->Add(fTrackPtAftTrkCuts);
 732
 733   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 734   fOutputList->Add(fTPCnsigma);
 735
 736   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 737   fOutputList->Add(fTrkEovPBef);
 738
 739   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 740   fOutputList->Add(fTrkEovPAft);
 741
 742   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
 743   fOutputList->Add(fdEdxBef);
 744
 745   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
 746   fOutputList->Add(fdEdxAft);
 747
 748   fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
 749   fOutputList->Add(fInvmassLS);
 750
 751   fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
 752   fOutputList->Add(fInvmassULS);
 753
 754   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 755   fOutputList->Add(fOpeningAngleLS);
 756
 757   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 758   fOutputList->Add(fOpeningAngleULS);
 759
 760   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 761   fOutputList->Add(fPhotoElecPt);
 762
 763   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
 764   fOutputList->Add(fSemiInclElecPt);
 765
 766   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 767   fOutputList->Add(fCent);
 768
 769   fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
 770   fOutputList->Add(fevPlaneV0A);
 771
 772   fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
 773   fOutputList->Add(fevPlaneV0C);
 774
 775   fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
 776   fOutputList->Add(fevPlaneTPC);
 777
 778   fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
 779   fOutputList->Add(fTPCsubEPres);
 780
 781   Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
 782   Double_t xminv1[4]={-1,-1,-1,0};
 783   Double_t xmaxv1[4]={1,1,1,90};
 784   fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
 785   fOutputList->Add(fEPres);
 786
 787   //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m20,m02
 788   Int_t binsv2[12]={100,100,90,100,100,100,100,100,3,3,100,100};
 789   Double_t xminv2[12]={0,-3.5,0,0,0,0,0,0,-1,-1,0,0};
 790   Double_t xmaxv2[12]={2*TMath::Pi(),3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi(),2,2,1,1};
 791   fCorr = new THnSparseD ("fCorr","Correlations",12,binsv2,xminv2,xmaxv2);
 792   fOutputList->Add(fCorr);
 793
 794   Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
 795   Double_t xminv3[5]={0,0,-1,-1,-1};
 796   Double_t xmaxv3[5]={90,50,1,1,1};
 797   feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
 798   fOutputList->Add(feV2);
 799
 800   Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 801   Double_t xminv4[5]={0,0,-1,-1,-1};
 802   Double_t xmaxv4[5]={90,50,1,1,1};
 803   fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
 804   fOutputList->Add(fphoteV2);
 805
 806   Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 807   Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
 808   Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
 809   fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
 810   fOutputList->Add(fChargPartV2);
 811
 812   Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 813   Double_t xminv6[5]={0,0,-1,-1,-1};
 814   Double_t xmaxv6[5]={90,50,1,1,1};
 815   feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
 816   fOutputList->Add(feTPCV2);
 817
 818   //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
 819   Int_t binsv7[15]={100,100,100,3,3,5,90,100,5,5,5,3,100,100,100};
 820   Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
 821   Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
 822   fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
 823   fOutputList->Add(fMCphotoElecPt);
 824
 825   fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
 826   fOutputList->Add(fGammaWeight);
 827
 828   fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
 829   fOutputList->Add(fPi0Weight);
 830
 831   fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
 832   fOutputList->Add(fEtaWeight);
 833
 834   fD0Weight = new TH2F("fD0Weight", "D0 weight",100,0,50,3,-1,2);
 835   fOutputList->Add(fD0Weight);
 836
 837   fDplusWeight = new TH2F("fDplusWeight", "D+ weight",100,0,50,3,-1,2);
 838   fOutputList->Add(fDplusWeight);
 839
 840   fDminusWeight = new TH2F("fDminusWeight", "D- weight",100,0,50,3,-1,2);
 841   fOutputList->Add(fDminusWeight);
 842
 843   fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
 844   fOutputList->Add(fD0_e);
 845
 846   for(Int_t k = 0; k < 6; k++) {
 847
 848     TString De_name = Form("fDe%d",k);
 849     TString D0e_name = Form("fD0e%d",k);
 850     TString Dpluse_name = Form("fDpluse%d",k);
 851     TString Dminuse_name = Form("fDminuse%d",k);
 852
 853     fDe[k] = new TH1F((const char*)De_name,"",100,0,50);
 854     fD0e[k] = new TH1F((const char*)D0e_name,"",100,0,50);
 855     fDpluse[k] = new TH1F((const char*)Dpluse_name,"",100,0,50);
 856     fDminuse[k] = new TH1F((const char*)Dminuse_name,"",100,0,50);
 857
 858     fOutputList->Add(fDe[k]);
 859     fOutputList->Add(fD0e[k]);
 860     fOutputList->Add(fDpluse[k]);
 861     fOutputList->Add(fDminuse[k]);
 862
 863   }
 864
 865   int nbin_v2 = 8;
 866   double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
 867
 868   fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
 869   fOutputList->Add(fTot_pi0e);
 870
 871   fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
 872   fOutputList->Add(fPhot_pi0e);
 873
 874   fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
 875   fOutputList->Add(fPhotBCG_pi0e);
 876
 877   fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
 878   fOutputList->Add(fTot_etae);
 879
 880   fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
 881   fOutputList->Add(fPhot_etae);
 882
 883   fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
 884   fOutputList->Add(fPhotBCG_etae);
 885
 886   //phiD,phie,phieRec,ptD,pte,pteRec,evPlaneV0,iHijing,kCent
 887   Int_t binsv8[9]={100,100,100,100,100,100,100,4,4};
 888   Double_t xminv8[9]={0,0,0,0,0,0,0,-2,-2};
 889   Double_t xmaxv8[9]={2*TMath::Pi(),2*TMath::Pi(),2*TMath::Pi(),50,50,50,TMath::Pi(),2,2};
 890   fCocktail = new THnSparseD ("fCocktail","Correlations",9,binsv8,xminv8,xmaxv8);
 891   fOutputList->Add(fCocktail);
 892
 893
 894   PostData(1,fOutputList);
 895 }
 896
 897 //________________________________________________________________________
 898 void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
 899 {
 900   // Info("Terminate");
 901         AliAnalysisTaskSE::Terminate();
 902 }
 903
 904 //________________________________________________________________________
 905 Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 906 {
 907   // Check single track cuts for a given cut step
 908   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 909   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 910   return kTRUE;
 911 }
 912 //_________________________________________
 913 void AliAnalysisTaskFlowTPCEMCalEP::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG)
 914 {
 915   //Identify non-heavy flavour electrons using Invariant mass method
 916
 917   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 918   fTrackCuts->SetRequireTPCRefit(kTRUE);
 919   fTrackCuts->SetRequireITSRefit(kTRUE);
 920   fTrackCuts->SetEtaRange(-0.7,0.7);
 921   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 922   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 923   fTrackCuts->SetMinNClustersTPC(100);
 924
 925   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 926
 927   Bool_t flagPhotonicElec = kFALSE;
 928   Bool_t flagPhotonicElecBCG = kFALSE;
 929
 930   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 931
 932     if(jTracks==iTracks) continue;
 933
 934     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
 935     if (!trackAsso) {
 936       printf("ERROR: Could not receive track %d\n", jTracks);
 937       continue;
 938     }
 939
 940     Double_t dEdxAsso = -999., ptAsso=-999.;
 941     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 942     Double_t openingAngle = -999., mass=999., width = -999;
 943
 944     dEdxAsso = trackAsso->GetTPCsignal();
 945     ptAsso = trackAsso->Pt();
 946     Int_t chargeAsso = trackAsso->Charge();
 947     Int_t charge = track->Charge();
 948
 949     if(ptAsso <0.5) continue;
 950     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
 951     if(dEdxAsso <65 || dEdxAsso>100) continue;
 952
 953     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 954     if(charge>0) fPDGe1 = -11;
 955     if(chargeAsso>0) fPDGe2 = -11;
 956
 957     if(charge == chargeAsso) fFlagLS = kTRUE;
 958     if(charge != chargeAsso) fFlagULS = kTRUE;
 959
 960     AliKFParticle ge1(*track, fPDGe1);
 961     AliKFParticle ge2(*trackAsso, fPDGe2);
 962     AliKFParticle recg(ge1, ge2);
 963
 964     if(recg.GetNDF()<1) continue;
 965     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 966     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 967
 968     AliKFVertex primV(*pVtx);
 969     primV += recg;
 970     recg.SetProductionVertex(primV);
 971
 972     recg.SetMassConstraint(0,0.0001);
 973
 974     openingAngle = ge1.GetAngle(ge2);
 975     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 976     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 977
 978     if(openingAngle > fOpeningAngleCut) continue;
 979
 980     recg.GetMass(mass,width);
 981
 982     if(fFlagLS) fInvmassLS->Fill(mass);
 983     if(fFlagULS) fInvmassULS->Fill(mass);
 984
 985     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
 986     if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
 987
 988   }
 989   fFlagPhotonicElec = flagPhotonicElec;
 990   fFlagPhotonicElecBCG = flagPhotonicElecBCG;
 991
 992 }
 993 //_________________________________________
 994 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
 995 {
 996   //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
 997   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 998   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 999   Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
1000
1001   return cos2DeltaPhi;
1002 }
1003 //_________________________________________
1004 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1005 {
1006   //Get phi-psi_EP
1007   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1008   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1009
1010   return dPhi;
1011 }
1012 //_________________________________________
1013 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT) const
1014 {
1015         //Get Pi0 weight
1016         double weight = 1.0;
1017
1018         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1019         {
1020
1021 //              weight = (2.877091*mcPi0pT)/(TMath::Power(0.706963+mcPi0pT/3.179309,17.336628)*exp(-mcPi0pT));
1022                 weight = (2.392024*mcPi0pT)/(TMath::Power(0.688810+mcPi0pT/3.005145,16.811845)*exp(-mcPi0pT));
1023         }
1024         else
1025         {
1026 //              weight = (0.0004*mcPi0pT)/TMath::Power(-0.176181+mcPi0pT/3.989747,5.629235);
1027                 weight = (0.000186*mcPi0pT)/TMath::Power(-0.606279+mcPi0pT/3.158642,4.365540);
1028         }
1029   return weight;
1030 }
1031 //_________________________________________
1032 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT) const
1033 {
1034   //Get eta weight
1035   double weight = 1.0;
1036
1037 //   weight = (0.818052*mcEtapT)/(TMath::Power(0.358651+mcEtapT/2.878631,9.494043));
1038   weight = (0.622703*mcEtapT)/(TMath::Power(0.323045+mcEtapT/2.736407,9.180356));
1039
1040   return weight;
1041 }
1042 //_________________________________________
1043 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDweight(Int_t whichD, Double_t mcDpT) const
1044 {
1045     //get D weights
1046     double weight = 1.0;
1047
1048     if (whichD == 0) weight = 0.271583*TMath::Landau(mcDpT,3.807103,1.536753,0); // D
1049     if (whichD == 1) weight = 0.300771*TMath::Landau(mcDpT,3.725771,1.496980,0); // D0
1050     if (whichD == 2) weight = 0.247280*TMath::Landau(mcDpT,3.746811,1.607551,0); // D+
1051     if (whichD == 3) weight = 0.249410*TMath::Landau(mcDpT,3.611508,1.632196,0); //D-
1052
1053   return weight;
1054 }