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