several updates for HFEv2 by TPC-EMCal EP
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowTPCEMCalEP.cxx
CommitLineData
1d08b6e8 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"
0f3a6a57 37#include "AliVEvent.h"
1d08b6e8 38
914042c2 39#include "AliAnalysisTaskFlowTPCEMCalEP.h"
1d08b6e8 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
3a94d6fe 65#include "AliMCEventHandler.h"
66#include "AliMCEvent.h"
67#include "AliMCParticle.h"
68#include "AliStack.h"
69
1d08b6e8 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
0f3a6a57 84#include "AliSelectNonHFE.h"
85
86
914042c2 87ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
1d08b6e8 88//________________________________________________________________________
914042c2 89AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
1d08b6e8 90 : AliAnalysisTaskSE(name)
91 ,fESD(0)
0f3a6a57 92 ,fAOD(0)
93 ,fVevent(0)
94 ,fpidResponse(0)
3a94d6fe 95 ,fMC(0)
1d08b6e8 96 ,fOutputList(0)
97 ,fTrackCuts(0)
98 ,fCuts(0)
0f3a6a57 99 ,fNonHFE(0)
1d08b6e8 100 ,fIdentifiedAsOutInz(kFALSE)
101 ,fPassTheEventCut(kFALSE)
102 ,fRejectKinkMother(kFALSE)
3a94d6fe 103 ,fIsMC(kFALSE)
0f3a6a57 104 ,fIsAOD(kFALSE)
1d08b6e8 105 ,fVz(0.0)
106 ,fCFM(0)
107 ,fPID(0)
108 ,fPIDqa(0)
109 ,fOpeningAngleCut(0.1)
0f3a6a57 110 ,fInvmassCut(0.05)
111 ,fChi2Cut(3.5)
112 ,fDCAcut(999)
113 ,fminCent(0)
114 ,fmaxCent(0)
115 ,fnonHFEalgorithm("KF")
1d08b6e8 116 ,fNoEvents(0)
117 ,fTrkpt(0)
118 ,fTrkEovPBef(0)
119 ,fTrkEovPAft(0)
120 ,fdEdxBef(0)
121 ,fdEdxAft(0)
1d08b6e8 122 ,fPhotoElecPt(0)
123 ,fSemiInclElecPt(0)
3a94d6fe 124 ,fMCphotoElecPt(0)
1d08b6e8 125 ,fTrackPtBefTrkCuts(0)
126 ,fTrackPtAftTrkCuts(0)
127 ,fTPCnsigma(0)
128 ,fCent(0)
44fdf7ac 129 ,fevPlaneV0A(0)
130 ,fevPlaneV0C(0)
0f3a6a57 131 ,fevPlaneV0(0)
44fdf7ac 132 ,fevPlaneTPC(0)
1d08b6e8 133 ,fTPCsubEPres(0)
134 ,fEPres(0)
135 ,fCorr(0)
136 ,feTPCV2(0)
137 ,feV2(0)
138 ,fphoteV2(0)
139 ,fChargPartV2(0)
6522b3b5 140 ,fGammaWeight(0)
141 ,fPi0Weight(0)
142 ,fEtaWeight(0)
def484dc 143 ,fD0Weight(0)
144 ,fDplusWeight(0)
145 ,fDminusWeight(0)
937b2ff2 146 ,fD0_e(0)
1c424962 147 ,fTot_pi0e(0)
148 ,fPhot_pi0e(0)
149 ,fPhotBCG_pi0e(0)
150 ,fTot_etae(0)
151 ,fPhot_etae(0)
152 ,fPhotBCG_etae(0)
0f3a6a57 153 ,fInvMass(0)
154 ,fInvMassBack(0)
155 ,fDCA(0)
156 ,fDCABack(0)
157 ,fOpAngle(0)
158 ,fOpAngleBack(0)
1d08b6e8 159{
160 //Named constructor
937b2ff2 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
1d08b6e8 169 fPID = new AliHFEpid("hfePid");
170 fTrackCuts = new AliESDtrackCuts();
0f3a6a57 171 fNonHFE = new AliSelectNonHFE();
172
173 InitParameters();
1d08b6e8 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//________________________________________________________________________
914042c2 186AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 187 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
188 ,fESD(0)
0f3a6a57 189 ,fAOD(0)
190 ,fVevent(0)
191 ,fpidResponse(0)
3a94d6fe 192 ,fMC(0)
1d08b6e8 193 ,fOutputList(0)
194 ,fTrackCuts(0)
195 ,fCuts(0)
0f3a6a57 196 ,fNonHFE(0)
1d08b6e8 197 ,fIdentifiedAsOutInz(kFALSE)
198 ,fPassTheEventCut(kFALSE)
199 ,fRejectKinkMother(kFALSE)
3a94d6fe 200 ,fIsMC(kFALSE)
0f3a6a57 201 ,fIsAOD(kFALSE)
1d08b6e8 202 ,fVz(0.0)
203 ,fCFM(0)
204 ,fPID(0)
205 ,fPIDqa(0)
206 ,fOpeningAngleCut(0.1)
0f3a6a57 207 ,fInvmassCut(0.05)
208 ,fChi2Cut(3.5)
209 ,fDCAcut(999)
210 ,fminCent(0)
211 ,fmaxCent(0)
212 ,fnonHFEalgorithm("KF")
1d08b6e8 213 ,fNoEvents(0)
214 ,fTrkpt(0)
215 ,fTrkEovPBef(0)
216 ,fTrkEovPAft(0)
217 ,fdEdxBef(0)
218 ,fdEdxAft(0)
1d08b6e8 219 ,fPhotoElecPt(0)
220 ,fSemiInclElecPt(0)
3a94d6fe 221 ,fMCphotoElecPt(0)
1d08b6e8 222 ,fTrackPtBefTrkCuts(0)
223 ,fTrackPtAftTrkCuts(0)
224 ,fTPCnsigma(0)
225 ,fCent(0)
44fdf7ac 226 ,fevPlaneV0A(0)
227 ,fevPlaneV0C(0)
0f3a6a57 228 ,fevPlaneV0(0)
44fdf7ac 229 ,fevPlaneTPC(0)
1d08b6e8 230 ,fTPCsubEPres(0)
231 ,fEPres(0)
232 ,fCorr(0)
233 ,feTPCV2(0)
234 ,feV2(0)
235 ,fphoteV2(0)
236 ,fChargPartV2(0)
6522b3b5 237 ,fGammaWeight(0)
238 ,fPi0Weight(0)
239 ,fEtaWeight(0)
def484dc 240 ,fD0Weight(0)
241 ,fDplusWeight(0)
242 ,fDminusWeight(0)
937b2ff2 243 ,fD0_e(0)
1c424962 244 ,fTot_pi0e(0)
245 ,fPhot_pi0e(0)
246 ,fPhotBCG_pi0e(0)
247 ,fTot_etae(0)
248 ,fPhot_etae(0)
249 ,fPhotBCG_etae(0)
0f3a6a57 250 ,fInvMass(0)
251 ,fInvMassBack(0)
252 ,fDCA(0)
253 ,fDCABack(0)
254 ,fOpAngle(0)
255 ,fOpAngleBack(0)
1d08b6e8 256{
0f3a6a57 257
258 //Default constructor
937b2ff2 259
0f3a6a57 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 }
937b2ff2 266
0f3a6a57 267 fPID = new AliHFEpid("hfePid");
268 fTrackCuts = new AliESDtrackCuts();
269 fNonHFE = new AliSelectNonHFE();
1d08b6e8 270
0f3a6a57 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());
1d08b6e8 283}
284//_________________________________________
285
914042c2 286AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 287{
288 //Destructor
289
290 delete fOutputList;
291 delete fPID;
292 delete fCFM;
293 delete fPIDqa;
294 delete fTrackCuts;
295}
296//_________________________________________
297
914042c2 298void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
1d08b6e8 299{
300 //Main loop
301 //Called for each event
0f3a6a57 302
1d08b6e8 303 // create pointer to event
304 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
0f3a6a57 305 if (!fESD){
1d08b6e8 306 printf("ERROR: fESD not available\n");
307 return;
308 }
0f3a6a57 309
310 fVevent = dynamic_cast<AliVEvent*>(InputEvent());
311 if(!fVevent){
312 printf("ERROR: fVEvent not available\n");
313 return;
314 }
315
1d08b6e8 316 if(!fCuts){
317 AliError("HFE cuts not available");
318 return;
319 }
0f3a6a57 320
1d08b6e8 321 if(!fPID->IsInitialized()){
322 // Initialize PID with the given run number
323 AliWarning("PID not initialised, get from Run no");
0f3a6a57 324 if(fIsAOD)fPID->InitializePID(fAOD->GetRunNumber());
325 if(!fIsAOD) fPID->InitializePID(fESD->GetRunNumber());
1d08b6e8 326 }
327
3a94d6fe 328 if(fIsMC)fMC = MCEvent();
329 AliStack* stack = NULL;
330 if(fIsMC && fMC) stack = fMC->Stack();
1d08b6e8 331
0f3a6a57 332 Int_t fNOtrks =fESD->GetNumberOfTracks();
1d08b6e8 333 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
0f3a6a57 334
1d08b6e8 335 Double_t pVtxZ = -999;
336 pVtxZ = pVtx->GetZ();
0f3a6a57 337
1d08b6e8 338 if(TMath::Abs(pVtxZ)>10) return;
339 fNoEvents->Fill(0);
0f3a6a57 340
1d08b6e8 341 if(fNOtrks<2) return;
0f3a6a57 342
343 fpidResponse = fInputHandler->GetPIDResponse();
344 if(!fpidResponse){
1d08b6e8 345 AliDebug(1, "Using default PID Response");
0f3a6a57 346 fpidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
1d08b6e8 347 }
0f3a6a57 348
349 fPID->SetPIDResponse(fpidResponse);
350
351 fCFM->SetRecEventInfo(fVevent);
352
1d08b6e8 353 Float_t cent = -1.;
354 AliCentrality *centrality = fESD->GetCentrality();
355 cent = centrality->GetCentralityPercentile("V0M");
356 fCent->Fill(cent);
0f3a6a57 357
358// cout<<"TEST: "<<fInvmassCut<< " "<< fOpeningAngleCut<<" "<<fnonHFEalgorithm<<" "<<fminCent<<" "<<fmaxCent<<" here!!!!!!!!!!!" <<endl;
535b11a0 359
1d08b6e8 360 //Event planes
0f3a6a57 361
1d08b6e8 362 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
363 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
0f3a6a57 364 fevPlaneV0A->Fill(evPlaneV0A);
365
1d08b6e8 366 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
367 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
0f3a6a57 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
1d08b6e8 374 AliEventplane* esdTPCep = fESD->GetEventplane();
44fdf7ac 375 TVector2 *standardQ = 0x0;
1d08b6e8 376 Double_t qx = -999., qy = -999.;
44fdf7ac 377 standardQ = esdTPCep->GetQVector();
ee89e744 378 if(!standardQ)return;
379
380 qx = standardQ->X();
381 qy = standardQ->Y();
0f3a6a57 382
1d08b6e8 383 TVector2 qVectorfortrack;
384 qVectorfortrack.Set(qx,qy);
385 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
0f3a6a57 386 fevPlaneTPC->Fill(evPlaneTPC);
387
388 //Event plane resolutions
389
390 // --> 2 subevent method (only for TPC EP)
391
1d08b6e8 392 TVector2 *qsub1a = esdTPCep->GetQsub1();
393 TVector2 *qsub2a = esdTPCep->GetQsub2();
394 Double_t evPlaneResTPC = -999.;
0f3a6a57 395 if(qsub1a && qsub2a){
396 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
1d08b6e8 397 }
0f3a6a57 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};
1d08b6e8 442 fEPres->Fill(evPlaneRes);
0f3a6a57 443
444 // MC
445 if(fIsMC && fMC && stack){
446 Int_t nParticles = stack->GetNtrack();
447 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
6522b3b5 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;
0f3a6a57 453
6522b3b5 454 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
455 int iHijing = 1;
456 if(!MChijing)iHijing = 0;
0f3a6a57 457
6522b3b5 458 if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
459 if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
def484dc 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
6522b3b5 464 Int_t idMother = particle->GetFirstMother();
465 if (idMother>0){
0f3a6a57 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 }
6522b3b5 472
0f3a6a57 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);
6522b3b5 478 }
6522b3b5 479
1d08b6e8 480 // Track loop
0f3a6a57 481 for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
482
483 AliVParticle* vparticle = fVevent->GetTrack(iTracks);
484 if (!vparticle){
1d08b6e8 485 printf("ERROR: Could not receive track %d\n", iTracks);
486 continue;
487 }
0f3a6a57 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
1d08b6e8 496 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
0f3a6a57 497
1d08b6e8 498 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
499 if(track->GetKinkIndex(0) != 0) continue;
500 }
0f3a6a57 501
1d08b6e8 502 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
0f3a6a57 503
1d08b6e8 504 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
0f3a6a57 505
1d08b6e8 506 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
0f3a6a57 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
1d08b6e8 513 pt = track->Pt();
0f3a6a57 514 if(pt<1.5) continue;
1d08b6e8 515 fTrkpt->Fill(pt);
0f3a6a57 516
1d08b6e8 517 Int_t clsId = track->GetEMCALcluster();
518 if (clsId>0){
519 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
520 if(cluster && cluster->IsEMCAL()){
0f3a6a57 521 clsE = cluster->E();
522 m20 = cluster->GetM20();
523 m02 = cluster->GetM02();
1d08b6e8 524 }
525 }
0f3a6a57 526
1d08b6e8 527 p = track->P();
528 phi = track->Phi();
529 dEdx = track->GetTPCsignal();
530 EovP = clsE/p;
b1a1875f 531 wEovP = wclsE/p;
1d08b6e8 532 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
533 fdEdxBef->Fill(p,dEdx);
534 fTPCnsigma->Fill(p,fTPCnSigma);
0f3a6a57 535
32d22164 536 //Remove electron candidate from the event plane
6522b3b5 537 Float_t evPlaneCorrTPC = evPlaneTPC;
32d22164 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
06d74a7b 546 Bool_t fFlagPhotonicElec = kFALSE;
547 Bool_t fFlagPhotonicElecBCG = kFALSE;
0f3a6a57 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};
1d08b6e8 568 fCorr->Fill(corr);
0f3a6a57 569
32d22164 570 Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
6522b3b5 571 Int_t whichPart = -99;
32d22164 572 Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
1c424962 573 Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
574 Double_t weight = 1.;
def484dc 575 Double_t Dweight = 1.;
32d22164 576 Bool_t MChijing;
0f3a6a57 577
1c424962 578 Bool_t pi0Decay= kFALSE;
579 Bool_t etaDecay= kFALSE;
535b11a0 580
581
582 Double_t phiD = -999.,phie = -999.,phieRec = -999.,ptD = -999.,pte = -999.,pteRec = -999.;
583
32d22164 584 if(fIsMC && fMC && stack){
06d74a7b 585 Int_t label = track->GetLabel();
586 if(label>0){
0f3a6a57 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
def484dc 601
0f3a6a57 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();
def484dc 617
0f3a6a57 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
1d08b6e8 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;
0f3a6a57 735
736 Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
737 GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 738 fChargPartV2->Fill(corrV2);
b3bd01b3 739
740 if(fTPCnSigma >= -0.5){
0f3a6a57 741 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
742 GetCos2DeltaPhi(phi,evPlaneV0C)};
743 feTPCV2->Fill(correctedV2);
b3bd01b3 744 }
0f3a6a57 745
b3bd01b3 746 if(pidpassed==0) continue;
0f3a6a57 747
748 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
749 GetCos2DeltaPhi(phi,evPlaneV0C)};
750
1d08b6e8 751 feV2->Fill(correctedV2);
1d08b6e8 752 fTrkEovPAft->Fill(pt,EovP);
753 fdEdxAft->Fill(p,dEdx);
3a94d6fe 754
0f3a6a57 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);
1d08b6e8 764}
765//_________________________________________
914042c2 766void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
1d08b6e8 767{
3a94d6fe 768 //--- Check MC
769 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
770 fIsMC = kTRUE;
771 printf("+++++ MC Data available");
772 }
1d08b6e8 773 //--------Initialize PID
3a94d6fe 774 fPID->SetHasMCData(fIsMC);
775
1d08b6e8 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
1d08b6e8 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
0f3a6a57 841 fevPlaneV0A = new TH1F("fevPlaneV0A","V0A EP",100,0,TMath::Pi());
44fdf7ac 842 fOutputList->Add(fevPlaneV0A);
843
0f3a6a57 844 fevPlaneV0C = new TH1F("fevPlaneV0C","V0C EP",100,0,TMath::Pi());
44fdf7ac 845 fOutputList->Add(fevPlaneV0C);
846
0f3a6a57 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());
44fdf7ac 851 fOutputList->Add(fevPlaneTPC);
852
0f3a6a57 853 fTPCsubEPres = new TH1F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1);
1d08b6e8 854 fOutputList->Add(fTPCsubEPres);
855
0f3a6a57 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);
1d08b6e8 860 fOutputList->Add(fEPres);
861
0f3a6a57 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);
1d08b6e8 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
44fdf7ac 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);
1d08b6e8 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
6522b3b5 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};
32d22164 895 Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
6522b3b5 896 Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
32d22164 897 fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
3a94d6fe 898 fOutputList->Add(fMCphotoElecPt);
6522b3b5 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);
1c424962 908
def484dc 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
908103bf 918 fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
937b2ff2 919 fOutputList->Add(fD0_e);
def484dc 920
937b2ff2 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 }
535b11a0 939
1c424962 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);
6522b3b5 945
1c424962 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);
535b11a0 960
0f3a6a57 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);
535b11a0 978
1d08b6e8 979 PostData(1,fOutputList);
980}
981
982//________________________________________________________________________
914042c2 983void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
1d08b6e8 984{
985 // Info("Terminate");
986 AliAnalysisTaskSE::Terminate();
987}
988
989//________________________________________________________________________
914042c2 990Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1d08b6e8 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//_________________________________________
914042c2 998Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 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}
1d08b6e8 1007//_________________________________________
914042c2 1008Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 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}
1c424962 1016//_________________________________________
0f3a6a57 1017Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT, Float_t cent) const
1c424962 1018{
937b2ff2 1019 //Get Pi0 weight
0f3a6a57 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 }
1c424962 1030 return weight;
1031}
1032//_________________________________________
0f3a6a57 1033Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT, Float_t cent) const
1c424962 1034{
937b2ff2 1035 //Get eta weight
1c424962 1036 double weight = 1.0;
1037
0f3a6a57 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));
1c424962 1040 return weight;
1041}
def484dc 1042//_________________________________________
0f3a6a57 1043Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDweight(Int_t whichD, Double_t mcDpT, Float_t cent) const
def484dc 1044{
0f3a6a57 1045 //get D weights
1046 double weight = 1.0;
937b2ff2 1047
0f3a6a57 1048 if (cent>30.0 && cent<50.0){
937b2ff2 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-
0f3a6a57 1053 }
def484dc 1054 return weight;
1055}
0f3a6a57 1056//_________________________________________
1057void 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}