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