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