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
 426     if(trackEP->Eta()>0 && trackEP->Eta()<0.8){
 427       Qx2pos += Qweight*TMath::Cos(2*trackEP->Phi());
 428       Qy2pos += Qweight*TMath::Sin(2*trackEP->Phi());
 429     }
 430     if(trackEP->Eta()<0 && trackEP->Eta()>-0.8){
 431       Qx2neg += Qweight*TMath::Cos(2*trackEP->Phi());
 432       Qy2neg += Qweight*TMath::Sin(2*trackEP->Phi());
 433     }
 434   }//track loop only for EP
 435
 436   Double_t evPlaneTPCneg = TMath::ATan2(Qy2neg, Qx2neg)/2;
 437   Double_t evPlaneTPCpos = TMath::ATan2(Qy2pos, Qx2pos)/2;
 438
 439   Double_t evPlaneRes[7]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),
 440       GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0,evPlaneTPCpos),
 441       GetCos2DeltaPhi(evPlaneV0,evPlaneTPCneg),GetCos2DeltaPhi(evPlaneTPCpos,evPlaneTPCneg),cent};
 442   fEPres->Fill(evPlaneRes);
 443
 444   // MC
 445   if(fIsMC && fMC && stack){
 446     Int_t nParticles = stack->GetNtrack();
 447     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 448       TParticle* particle = stack->Particle(iParticle);
 449       int fPDG = particle->GetPdgCode();
 450       double pTMC = particle->Pt();
 451       double etaMC = particle->Eta();
 452       if(fabs(etaMC)>0.7)continue;
 453
 454       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
 455       int iHijing = 1;
 456       if(!MChijing)iHijing = 0;
 457
 458       if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
 459       if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
 460       if(fPDG==421)fD0Weight->Fill(pTMC,iHijing);//D0
 461       if(fPDG==411)fDplusWeight->Fill(pTMC,iHijing);//D+
 462       if(fPDG==-411)fDminusWeight->Fill(pTMC,iHijing);//D-
 463
 464       Int_t idMother = particle->GetFirstMother();
 465       if (idMother>0){
 466         TParticle *mother = stack->Particle(idMother);
 467         int motherPDG = mother->GetPdgCode();
 468         if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
 469       }
 470     }
 471   }
 472
 473   Double_t ptRange[8] = {2, 2.5, 3, 4, 6, 8, 10, 13};
 474   Double_t ptDmeson[7] = {2,3,4,6,8,12,16};
 475   Double_t deltaPhiRange[7];
 476   for(Int_t j=0;j<7;j++){
 477     deltaPhiRange[j] = j*(TMath::Pi()/6);
 478   }
 479
 480   // Track loop
 481   for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
 482
 483     AliVParticle* vparticle = fVevent->GetTrack(iTracks);
 484     if (!vparticle){
 485       printf("ERROR: Could not receive track %d\n", iTracks);
 486       continue;
 487     }
 488
 489     AliVTrack *vtrack = dynamic_cast<AliVTrack*>(vparticle);
 490     AliESDtrack *track = dynamic_cast<AliESDtrack*>(vparticle);
 491
 492     if (TMath::Abs(track->Eta())>0.7) continue;
 493
 494     fTrackPtBefTrkCuts->Fill(track->Pt());
 495
 496     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 497
 498     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 499       if(track->GetKinkIndex(0) != 0) continue;
 500     }
 501
 502     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 503
 504     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 505
 506     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 507
 508     fTrackPtAftTrkCuts->Fill(track->Pt());
 509
 510     Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999.,
 511       wclsE = -999., wEovP = -999., m02= -999., m20= -999.;
 512
 513     pt = track->Pt();
 514     if(pt<1.5) continue;
 515     fTrkpt->Fill(pt);
 516
 517     Int_t clsId = track->GetEMCALcluster();
 518     if (clsId>0){
 519       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
 520       if(cluster && cluster->IsEMCAL()){
 521         clsE = cluster->E();
 522         m20 = cluster->GetM20();
 523         m02 = cluster->GetM02();
 524       }
 525     }
 526
 527     p = track->P();
 528     phi = track->Phi();
 529     dEdx = track->GetTPCsignal();
 530     EovP = clsE/p;
 531     wEovP = wclsE/p;
 532     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 533     fdEdxBef->Fill(p,dEdx);
 534     fTPCnsigma->Fill(p,fTPCnSigma);
 535
 536     //Remove electron candidate from the event plane
 537     Float_t evPlaneCorrTPC = evPlaneTPC;
 538     if(dEdx>70 && dEdx<90){
 539       Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
 540       Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
 541       TVector2 newQVectorfortrack;
 542       newQVectorfortrack.Set(qX,qY);
 543       evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
 544     }
 545
 546     Bool_t fFlagPhotonicElec = kFALSE;
 547     Bool_t fFlagPhotonicElecBCG = kFALSE;
 548
 549     //Non-HFE reconstruction
 550     fNonHFE->SetPIDresponse(fpidResponse);
 551     fNonHFE->FindNonHFE(iTracks,vparticle,fVevent);
 552     Int_t *fUlsPartner = fNonHFE->GetPartnersULS(); // Pointer to the ULS partners index
 553     Int_t *fLsPartner = fNonHFE->GetPartnersLS(); // Pointer to the LS partners index
 554
 555     if (fNonHFE->IsULS()) fFlagPhotonicElec=kTRUE;
 556     if (fNonHFE->IsLS()) fFlagPhotonicElecBCG=kTRUE;
 557
 558     fNonHFE->SetHistAngleBack(fOpAngleBack);
 559     fNonHFE->SetHistAngle(fOpAngle);
 560     fNonHFE->SetHistMassBack(fInvMassBack);
 561     fNonHFE->SetHistMass(fInvMass);
 562     if (fnonHFEalgorithm == "DCA")fNonHFE->SetHistDCABack(fDCABack);
 563     if (fnonHFEalgorithm == "DCA")fNonHFE->SetHistDCA(fDCA);
 564
 565
 566     Double_t corr[11]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneV0),GetCos2DeltaPhi(phi,evPlaneV0), static_cast<Double_t>(fFlagPhotonicElec),
 567     static_cast<Double_t>(fFlagPhotonicElecBCG),m02,m20};
 568     fCorr->Fill(corr);
 569
 570     Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
 571     Int_t whichPart = -99;
 572     Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
 573     Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
 574     Double_t weight = 1.;
 575     Double_t Dweight = 1.;
 576     Bool_t MChijing;
 577
 578     Bool_t pi0Decay= kFALSE;
 579     Bool_t etaDecay= kFALSE;
 580
 581
 582     Double_t phiD = -999.,phie = -999.,phieRec = -999.,ptD = -999.,pte = -999.,pteRec = -999.;
 583
 584     if(fIsMC && fMC && stack){
 585       Int_t label = track->GetLabel();
 586       if(label>0){
 587         TParticle *particle = stack->Particle(label);
 588         if(particle){
 589           partPDG = particle->GetPdgCode();
 590           partPt = particle->Pt();
 591
 592           if (TMath::Abs(partPDG)==11) whichPart = 0; //electron
 593           if (partPDG==22) whichPart = 3; //gamma
 594           if (partPDG==111) whichPart = 2; //pi0
 595           if (partPDG==221) whichPart = 1; //eta
 596
 597           MChijing = fMC->IsFromBGEvent(label);
 598
 599           int iHijing = 1;
 600           if(!MChijing) iHijing = 0; // 0 if enhanced sample
 601
 602           Int_t idMother = particle->GetFirstMother();
 603           if (idMother>0){
 604             TParticle *mother = stack->Particle(idMother);
 605             motherPt = mother->Pt();
 606             motherPDG = mother->GetPdgCode();
 607
 608             if (motherPDG==22) whichFirstMother = 3; //gamma
 609             if (motherPDG==111) whichFirstMother = 2; //pi0
 610             if (motherPDG==221) whichFirstMother = 1; //eta
 611
 612             Int_t idSecondMother = particle->GetSecondMother();
 613             if (idSecondMother>0){
 614               TParticle *secondMother = stack->Particle(idSecondMother);
 615               secondMotherPt = secondMother->Pt();
 616               secondMotherPDG = secondMother->GetPdgCode();
 617
 618               if (secondMotherPDG==111) whichSecondMother = 2; //pi0
 619               if (secondMotherPDG==221) whichSecondMother = 1; //eta
 620
 621               Int_t idThirdMother = secondMother->GetFirstMother();
 622               if (idThirdMother>0){
 623                 TParticle *thirdMother = stack->Particle(idThirdMother);
 624                 thirdMotherPt = thirdMother->Pt();
 625                 thirdMotherPDG = thirdMother->GetPdgCode();
 626
 627                 if (thirdMotherPDG==221) whichThirdMother = 1; //eta
 628               }//third mother
 629             }//second mother
 630
 631
 632             if (TMath::Abs(partPDG)==11){
 633
 634               // D meson decay
 635               if (motherPDG==421 || TMath::Abs(motherPDG)==411){ // D
 636                 Double_t phi_D = -99., Deltaphi_De=-99.;
 637                 phi_D = mother->Phi();
 638                 Deltaphi_De = phi_D - phi;
 639
 640                 fD0_e->Fill(pt,Deltaphi_De);
 641
 642                 Dweight= GetDweight(0,motherPt,cent);
 643                 if(iHijing==1) Dweight = 1.0;
 644                 for(Int_t i=0;i<6;i++){
 645                   if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fDe[i]->Fill(pt,Dweight);
 646                 }
 647               }
 648               if (motherPDG==421){ // D0
 649                 Dweight= GetDweight(1,motherPt,cent);
 650                 if(iHijing==1) Dweight = 1.0;
 651                 for(Int_t i=0;i<6;i++){
 652                   if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fD0e[i]->Fill(pt,Dweight);
 653                 }
 654               }
 655               if (motherPDG==411){ // D+
 656                 Dweight= GetDweight(2,motherPt,cent);
 657                 if(iHijing==1) Dweight = 1.0;
 658                 for(Int_t i=0;i<6;i++){
 659                   if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fDpluse[i]->Fill(pt,Dweight);
 660                 }
 661               }
 662               if (motherPDG==-411){ //D-
 663                 Dweight= GetDweight(3,motherPt,cent);
 664                 if(iHijing==1) Dweight = 1.0;
 665                 for(Int_t i=0;i<6;i++){
 666                   if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fDminuse[i]->Fill(pt,Dweight);
 667                 }
 668               }
 669
 670               //pi0 decay
 671               if (motherPDG==111 && secondMotherPDG!=221){ //not eta -> pi0 -> e
 672                 weight = GetPi0weight(motherPt,cent);
 673                 pi0Decay = kTRUE;
 674               }
 675               if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG!=221){ //not eta -> pi0 -> gamma -> e
 676                 weight = GetPi0weight(secondMotherPt,cent);
 677                 pi0Decay = kTRUE;
 678               }
 679
 680               //eta decay
 681               if (motherPDG==221){ //eta -> e
 682                 weight = GetEtaweight(motherPt,cent);
 683                 etaDecay = kTRUE;
 684               }
 685               if (motherPDG==111 && secondMotherPDG==221){ //eta -> pi0 -> e
 686                 weight = GetEtaweight(secondMotherPt,cent);
 687                 etaDecay = kTRUE;
 688               }
 689               if (motherPDG==22 && secondMotherPDG==221){ //eta -> gamma -> e
 690                 weight = GetEtaweight(secondMotherPt,cent);
 691                 etaDecay = kTRUE;
 692               }
 693               if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221){ //eta -> pi0 -> gamma -> e
 694                 weight = GetEtaweight(thirdMotherPt,cent);
 695                 etaDecay = kTRUE;
 696               }
 697             }// end of electron
 698
 699             if (fTPCnSigma>-1 && fTPCnSigma<3 && EovP>1 && EovP<1.3 && (motherPDG==22 || motherPDG==111 || motherPDG==221)){
 700               if(iHijing==1) weight = 1.0;
 701
 702               if (pi0Decay){
 703                 fTot_pi0e->Fill(partPt,weight);
 704                 if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
 705                 if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
 706               }
 707               if (etaDecay){
 708                 fTot_etae->Fill(partPt,weight);
 709                 if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
 710                 if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
 711               }
 712             }
 713
 714             Double_t mc[15]={EovP,fTPCnSigma,partPt, static_cast<Double_t>(fFlagPhotonicElec),
 715                 static_cast<Double_t>(fFlagPhotonicElecBCG), static_cast<Double_t>(whichPart),cent,pt,
 716                   static_cast<Double_t>(whichFirstMother), static_cast<Double_t>(whichSecondMother),
 717             static_cast<Double_t>(whichThirdMother), static_cast<Double_t>(iHijing),motherPt,secondMotherPt,thirdMotherPt};
 718
 719             //fMCphotoElecPt->Fill(mc);
 720             if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
 721           }//end mother
 722         }// end particle
 723       }// end label
 724     }//end MC
 725
 726     if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
 727     Int_t pidpassed = 1;
 728
 729     //--- track accepted
 730     AliHFEpidObject hfetrack;
 731     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 732     hfetrack.SetRecTrack(track);
 733     hfetrack.SetPbPb();
 734     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 735
 736     Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
 737     GetCos2DeltaPhi(phi,evPlaneV0C)};
 738     fChargPartV2->Fill(corrV2);
 739
 740     if(fTPCnSigma >= -0.5){
 741       Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
 742       GetCos2DeltaPhi(phi,evPlaneV0C)};
 743       feTPCV2->Fill(correctedV2);
 744     }
 745
 746     if(pidpassed==0) continue;
 747
 748     Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
 749      GetCos2DeltaPhi(phi,evPlaneV0C)};
 750
 751     feV2->Fill(correctedV2);
 752     fTrkEovPAft->Fill(pt,EovP);
 753     fdEdxAft->Fill(p,dEdx);
 754
 755     if(fFlagPhotonicElec){
 756       fphoteV2->Fill(correctedV2);
 757       fPhotoElecPt->Fill(pt);
 758     }
 759
 760     if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
 761
 762   }//end of track loop
 763   PostData(1, fOutputList);
 764 }
 765 //_________________________________________
 766 void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
 767 {
 768   //--- Check MC
 769   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 770     fIsMC = kTRUE;
 771     printf("+++++ MC Data available");
 772   }
 773   //--------Initialize PID
 774   fPID->SetHasMCData(fIsMC);
 775
 776   if(!fPID->GetNumberOfPIDdetectors())
 777     {
 778       fPID->AddDetector("TPC", 0);
 779       fPID->AddDetector("EMCAL", 1);
 780     }
 781
 782   fPID->SortDetectors();
 783   fPIDqa = new AliHFEpidQAmanager();
 784   fPIDqa->Initialize(fPID);
 785
 786   //--------Initialize correction Framework and Cuts
 787   fCFM = new AliCFManager;
 788   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 789   fCFM->SetNStepParticle(kNcutSteps);
 790   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 791     fCFM->SetParticleCutsList(istep, NULL);
 792
 793   if(!fCuts){
 794     AliWarning("Cuts not available. Default cuts will be used");
 795     fCuts = new AliHFEcuts;
 796     fCuts->CreateStandardCuts();
 797   }
 798   fCuts->Initialize(fCFM);
 799
 800   //---------Output Tlist
 801   fOutputList = new TList();
 802   fOutputList->SetOwner();
 803   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 804
 805   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 806   fOutputList->Add(fNoEvents);
 807
 808   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 809   fOutputList->Add(fTrkpt);
 810
 811   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 812   fOutputList->Add(fTrackPtBefTrkCuts);
 813
 814   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 815   fOutputList->Add(fTrackPtAftTrkCuts);
 816
 817   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 818   fOutputList->Add(fTPCnsigma);
 819
 820   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 821   fOutputList->Add(fTrkEovPBef);
 822
 823   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 824   fOutputList->Add(fTrkEovPAft);
 825
 826   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
 827   fOutputList->Add(fdEdxBef);
 828
 829   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
 830   fOutputList->Add(fdEdxAft);
 831
 832   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 833   fOutputList->Add(fPhotoElecPt);
 834
 835   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
 836   fOutputList->Add(fSemiInclElecPt);
 837
 838   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 839   fOutputList->Add(fCent);
 840
 841   fevPlaneV0A = new TH1F("fevPlaneV0A","V0A EP",100,0,TMath::Pi());
 842   fOutputList->Add(fevPlaneV0A);
 843
 844   fevPlaneV0C = new TH1F("fevPlaneV0C","V0C EP",100,0,TMath::Pi());
 845   fOutputList->Add(fevPlaneV0C);
 846
 847   fevPlaneV0 = new TH1F("fevPlaneV0","V0 EP",100,0,TMath::Pi());
 848   fOutputList->Add(fevPlaneV0);
 849
 850   fevPlaneTPC = new TH1F("fevPlaneTPC","TPC EP",100,0,TMath::Pi());
 851   fOutputList->Add(fevPlaneTPC);
 852
 853   fTPCsubEPres = new TH1F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1);
 854   fOutputList->Add(fTPCsubEPres);
 855
 856   Int_t binsv1[7]={100,100,100,100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, V0-TPCpos, V0-TPCneg, TPCpos-TPCneg, cent
 857   Double_t xminv1[7]={-1,-1,-1,-1,-1,-1,0};
 858   Double_t xmaxv1[7]={1,1,1,1,1,1,90};
 859   fEPres = new THnSparseD ("fEPres","EP resolution",7,binsv1,xminv1,xmaxv1);
 860   fOutputList->Add(fEPres);
 861
 862   //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneV0),GetCos2DeltaPhi(phi,evPlaneV0), fFlagPhotonicElec, fFlagPhotonicElecBCG,m02,m20
 863   Int_t binsv2[11]={100,200,90,100,100,100,100,3,3,100,100};
 864   Double_t xminv2[11]={0,-5,0,0,0,0,-1,-1,-1,0,0};
 865   Double_t xmaxv2[11]={2*TMath::Pi(),5,90,50,3,TMath::Pi(),1,2,2,1,1};
 866   fCorr = new THnSparseD ("fCorr","Correlations",11,binsv2,xminv2,xmaxv2);
 867   fOutputList->Add(fCorr);
 868
 869   Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
 870   Double_t xminv3[5]={0,0,-1,-1,-1};
 871   Double_t xmaxv3[5]={90,50,1,1,1};
 872   feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
 873   fOutputList->Add(feV2);
 874
 875   Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 876   Double_t xminv4[5]={0,0,-1,-1,-1};
 877   Double_t xmaxv4[5]={90,50,1,1,1};
 878   fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
 879   fOutputList->Add(fphoteV2);
 880
 881   Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 882   Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
 883   Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
 884   fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
 885   fOutputList->Add(fChargPartV2);
 886
 887   Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 888   Double_t xminv6[5]={0,0,-1,-1,-1};
 889   Double_t xmaxv6[5]={90,50,1,1,1};
 890   feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
 891   fOutputList->Add(feTPCV2);
 892
 893   //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
 894   Int_t binsv7[15]={100,100,100,3,3,5,90,100,5,5,5,3,100,100,100};
 895   Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
 896   Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
 897   fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
 898   fOutputList->Add(fMCphotoElecPt);
 899
 900   fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
 901   fOutputList->Add(fGammaWeight);
 902
 903   fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
 904   fOutputList->Add(fPi0Weight);
 905
 906   fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
 907   fOutputList->Add(fEtaWeight);
 908
 909   fD0Weight = new TH2F("fD0Weight", "D0 weight",100,0,50,3,-1,2);
 910   fOutputList->Add(fD0Weight);
 911
 912   fDplusWeight = new TH2F("fDplusWeight", "D+ weight",100,0,50,3,-1,2);
 913   fOutputList->Add(fDplusWeight);
 914
 915   fDminusWeight = new TH2F("fDminusWeight", "D- weight",100,0,50,3,-1,2);
 916   fOutputList->Add(fDminusWeight);
 917
 918   fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
 919   fOutputList->Add(fD0_e);
 920
 921   for(Int_t k = 0; k < 6; k++) {
 922
 923     TString De_name = Form("fDe%d",k);
 924     TString D0e_name = Form("fD0e%d",k);
 925     TString Dpluse_name = Form("fDpluse%d",k);
 926     TString Dminuse_name = Form("fDminuse%d",k);
 927
 928     fDe[k] = new TH1F((const char*)De_name,"",100,0,50);
 929     fD0e[k] = new TH1F((const char*)D0e_name,"",100,0,50);
 930     fDpluse[k] = new TH1F((const char*)Dpluse_name,"",100,0,50);
 931     fDminuse[k] = new TH1F((const char*)Dminuse_name,"",100,0,50);
 932
 933     fOutputList->Add(fDe[k]);
 934     fOutputList->Add(fD0e[k]);
 935     fOutputList->Add(fDpluse[k]);
 936     fOutputList->Add(fDminuse[k]);
 937
 938   }
 939
 940   int nbin_v2 = 8;
 941   double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
 942
 943   fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
 944   fOutputList->Add(fTot_pi0e);
 945
 946   fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
 947   fOutputList->Add(fPhot_pi0e);
 948
 949   fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
 950   fOutputList->Add(fPhotBCG_pi0e);
 951
 952   fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
 953   fOutputList->Add(fTot_etae);
 954
 955   fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
 956   fOutputList->Add(fPhot_etae);
 957
 958   fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
 959   fOutputList->Add(fPhotBCG_etae);
 960
 961   fInvMass = new TH1F("fInvMass","",200,0,0.3);
 962   fOutputList->Add(fInvMass);
 963
 964   fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
 965   fOutputList->Add(fInvMassBack);
 966
 967   fDCA = new TH1F("fDCA","",200,0,1);
 968   fOutputList->Add(fDCA);
 969
 970   fDCABack = new TH1F("fDCABack","",200,0,1);
 971   fOutputList->Add(fDCABack);
 972
 973   fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
 974   fOutputList->Add(fOpAngle);
 975
 976   fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
 977   fOutputList->Add(fOpAngleBack);
 978
 979   PostData(1,fOutputList);
 980 }
 981
 982 //________________________________________________________________________
 983 void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
 984 {
 985   // Info("Terminate");
 986         AliAnalysisTaskSE::Terminate();
 987 }
 988
 989 //________________________________________________________________________
 990 Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 991 {
 992   // Check single track cuts for a given cut step
 993   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 994   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 995   return kTRUE;
 996 }
 997 //_________________________________________
 998 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
 999 {
1000   //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
1001   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1002   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1003   Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
1004
1005   return cos2DeltaPhi;
1006 }
1007 //_________________________________________
1008 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1009 {
1010   //Get phi-psi_EP
1011   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1012   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1013
1014   return dPhi;
1015 }
1016 //_________________________________________
1017 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT, Float_t cent) const
1018 {
1019         //Get Pi0 weight
1020   double weight = 1.0;
1021
1022   if(mcPi0pT>0.0 && mcPi0pT<5.0){
1023           if (cent>20.0 && cent<40.0) weight = (2.877091*mcPi0pT)/(TMath::Power(0.706963+mcPi0pT/3.179309,17.336628)*exp(-mcPi0pT));
1024                 if (cent>30.0 && cent<50.0) weight = (2.392024*mcPi0pT)/(TMath::Power(0.688810+mcPi0pT/3.005145,16.811845)*exp(-mcPi0pT));
1025   }
1026   else{
1027     if (cent>20.0 && cent<40.0) weight = (0.0004*mcPi0pT)/TMath::Power(-0.176181+mcPi0pT/3.989747,5.629235);
1028                 if (cent>30.0 && cent<50.0) weight = (0.000186*mcPi0pT)/TMath::Power(-0.606279+mcPi0pT/3.158642,4.365540);
1029   }
1030   return weight;
1031 }
1032 //_________________________________________
1033 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT, Float_t cent) const
1034 {
1035   //Get eta weight
1036   double weight = 1.0;
1037
1038   if (cent>20.0 && cent<40.0) weight = (0.818052*mcEtapT)/(TMath::Power(0.358651+mcEtapT/2.878631,9.494043));
1039   if (cent>30.0 && cent<50.0) weight = (0.622703*mcEtapT)/(TMath::Power(0.323045+mcEtapT/2.736407,9.180356));
1040   return weight;
1041 }
1042 //_________________________________________
1043 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDweight(Int_t whichD, Double_t mcDpT, Float_t cent) const
1044 {
1045   //get D weights
1046   double weight = 1.0;
1047
1048   if (cent>30.0 && cent<50.0){
1049     if (whichD == 0) weight = 0.271583*TMath::Landau(mcDpT,3.807103,1.536753,0); // D
1050     if (whichD == 1) weight = 0.300771*TMath::Landau(mcDpT,3.725771,1.496980,0); // D0
1051     if (whichD == 2) weight = 0.247280*TMath::Landau(mcDpT,3.746811,1.607551,0); // D+
1052     if (whichD == 3) weight = 0.249410*TMath::Landau(mcDpT,3.611508,1.632196,0); //D-
1053   }
1054   return weight;
1055 }
1056 //_________________________________________
1057 void AliAnalysisTaskFlowTPCEMCalEP::InitParameters()
1058 {
1059   // Init parameters
1060
1061   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1062   fTrackCuts->SetRequireTPCRefit(kTRUE);
1063   fTrackCuts->SetRequireITSRefit(kTRUE);
1064   fTrackCuts->SetEtaRange(-0.7,0.7);
1065   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1066   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1067   fTrackCuts->SetMinNClustersTPC(100);
1068   fTrackCuts->SetPtRange(0.5,100);
1069
1070   fNonHFE->SetAODanalysis(kFALSE);
1071   fNonHFE->SetInvariantMassCut(fInvmassCut);
1072   fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
1073   fNonHFE->SetChi2OverNDFCut(fChi2Cut);
1074   fNonHFE->SetAlgorithm(fnonHFEalgorithm); //KF or DCA
1075   if (fnonHFEalgorithm == "DCA") fNonHFE->SetDCACut(fDCAcut);
1076   fNonHFE->SetTrackCuts(-3.5,3.5,fTrackCuts);
1077
1078 }