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