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