]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalEP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowTPCEMCalEP.cxx
... / ...
CommitLineData
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
87ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
88//________________________________________________________________________
89AliAnalysisTaskFlowTPCEMCalEP::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 ,fSetMassConstraint(kFALSE)
106 ,fVz(0.0)
107 ,fCFM(0)
108 ,fPID(0)
109 ,fPIDqa(0)
110 ,fOpeningAngleCut(1000.)
111 ,fInvmassCut(0.05)
112 ,fChi2Cut(3.5)
113 ,fDCAcut(999)
114 ,fnonHFEalgorithm("KF")
115 ,fNoEvents(0)
116 ,fTrkpt(0)
117 ,fTrkEovPBef(0)
118 ,fTrkEovPAft(0)
119 ,fdEdxBef(0)
120 ,fdEdxAft(0)
121 ,fPhotoElecPt(0)
122 ,fSemiInclElecPt(0)
123 ,fTrackPtBefTrkCuts(0)
124 ,fTrackPtAftTrkCuts(0)
125 ,fTPCnsigma(0)
126 ,fCent(0)
127 ,fTPCsubEPres(0)
128 ,fEPres(0)
129 ,fCorr(0)
130 ,fD0_e(0)
131 ,fTot_pi0e(0)
132 ,fPhot_pi0e(0)
133 ,fPhotBCG_pi0e(0)
134 ,fTot_etae(0)
135 ,fPhot_etae(0)
136 ,fPhotBCG_etae(0)
137 ,fInvMass(0)
138 ,fInvMassBack(0)
139 ,fDCA(0)
140 ,fDCABack(0)
141 ,fOpAngle(0)
142 ,fOpAngleBack(0)
143{
144 //Named constructor
145
146 for(Int_t k = 0; k < 3; k++) {
147 fevPlaneV0[k] = NULL;
148 feTPCV2[k] = NULL;
149 feV2[k] = NULL;
150 fChargPartV2[k] = NULL;
151 fMtcPartV2[k] = NULL;
152
153 fPi0Pt[k] = NULL;
154 fEtaPt[k] = NULL;
155 fInvmassLS[k] = NULL;
156 fInvmassULS[k] = NULL;
157 fOpeningAngleLS[k] = NULL;
158 fOpeningAngleULS[k] = NULL;
159 }
160
161 for(Int_t i=0; i<3; i++) {
162 for(Int_t j=0; j<8; j++) {
163 for(Int_t k=0; k<4; k++) {
164 fEoverPsig[i][j][k] = NULL;
165 fEoverPuls[i][j][k] = NULL;
166 fEoverPls[i][j][k] = NULL;
167 fEoverPbcg[i][j][k] = NULL;
168 }
169 }
170 }
171
172 for(Int_t k = 0; k < 6; k++) {
173 fDe[k]= NULL;
174 fD0e[k]= NULL;
175 fDpluse[k]= NULL;
176 fDminuse[k]= NULL;
177 }
178
179 fPID = new AliHFEpid("hfePid");
180 fTrackCuts = new AliESDtrackCuts();
181 fNonHFE = new AliSelectNonHFE();
182
183 InitParameters();
184
185 // Define input and output slots here
186 // Input slot #0 works with a TChain
187 DefineInput(0, TChain::Class());
188 // Output slot #0 id reserved by the base class for AOD
189 // Output slot #1 writes into a TH1 container
190 // DefineOutput(1, TH1I::Class());
191 DefineOutput(1, TList::Class());
192 // DefineOutput(3, TTree::Class());
193}
194
195//________________________________________________________________________
196AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
197 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
198 ,fESD(0)
199 ,fAOD(0)
200 ,fVevent(0)
201 ,fpidResponse(0)
202 ,fMC(0)
203 ,fOutputList(0)
204 ,fTrackCuts(0)
205 ,fCuts(0)
206 ,fNonHFE(0)
207 ,fIdentifiedAsOutInz(kFALSE)
208 ,fPassTheEventCut(kFALSE)
209 ,fRejectKinkMother(kFALSE)
210 ,fIsMC(kFALSE)
211 ,fIsAOD(kFALSE)
212 ,fSetMassConstraint(kFALSE)
213 ,fVz(0.0)
214 ,fCFM(0)
215 ,fPID(0)
216 ,fPIDqa(0)
217 ,fOpeningAngleCut(1000.)
218 ,fInvmassCut(0.05)
219 ,fChi2Cut(3.5)
220 ,fDCAcut(999)
221 ,fnonHFEalgorithm("KF")
222 ,fNoEvents(0)
223 ,fTrkpt(0)
224 ,fTrkEovPBef(0)
225 ,fTrkEovPAft(0)
226 ,fdEdxBef(0)
227 ,fdEdxAft(0)
228 ,fPhotoElecPt(0)
229 ,fSemiInclElecPt(0)
230 ,fTrackPtBefTrkCuts(0)
231 ,fTrackPtAftTrkCuts(0)
232 ,fTPCnsigma(0)
233 ,fCent(0)
234 ,fTPCsubEPres(0)
235 ,fEPres(0)
236 ,fCorr(0)
237 ,fD0_e(0)
238 ,fTot_pi0e(0)
239 ,fPhot_pi0e(0)
240 ,fPhotBCG_pi0e(0)
241 ,fTot_etae(0)
242 ,fPhot_etae(0)
243 ,fPhotBCG_etae(0)
244 ,fInvMass(0)
245 ,fInvMassBack(0)
246 ,fDCA(0)
247 ,fDCABack(0)
248 ,fOpAngle(0)
249 ,fOpAngleBack(0)
250{
251
252 //Default constructor
253
254 for(Int_t k = 0; k < 3; k++) {
255 fevPlaneV0[k] = NULL;
256 feTPCV2[k] = NULL;
257 feV2[k] = NULL;
258 fChargPartV2[k] = NULL;
259 fMtcPartV2[k] = NULL;
260
261 fPi0Pt[k] = NULL;
262 fEtaPt[k] = NULL;
263 fInvmassLS[k] = NULL;
264 fInvmassULS[k] = NULL;
265 fOpeningAngleLS[k] = NULL;
266 fOpeningAngleULS[k] = NULL;
267 }
268
269 for(Int_t i=0; i<3; i++) {
270 for(Int_t j=0; j<8; j++) {
271 for(Int_t k=0; k<4; k++) {
272 fEoverPsig[i][j][k] = NULL;
273 fEoverPuls[i][j][k] = NULL;
274 fEoverPls[i][j][k] = NULL;
275 fEoverPbcg[i][j][k] = NULL;
276 }
277 }
278 }
279
280 for(Int_t k = 0; k < 6; k++) {
281 fDe[k]= NULL;
282 fD0e[k]= NULL;
283 fDpluse[k]= NULL;
284 fDminuse[k]= NULL;
285 }
286
287 fPID = new AliHFEpid("hfePid");
288 fTrackCuts = new AliESDtrackCuts();
289 fNonHFE = new AliSelectNonHFE();
290
291 InitParameters();
292
293
294 // Constructor
295 // Define input and output slots here
296 // Input slot #0 works with a TChain
297 DefineInput(0, TChain::Class());
298 // Output slot #0 id reserved by the base class for AOD
299 // Output slot #1 writes into a TH1 container
300 // DefineOutput(1, TH1I::Class());
301 DefineOutput(1, TList::Class());
302 //DefineOutput(3, TTree::Class());
303}
304//_________________________________________
305
306AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
307{
308 //Destructor
309
310 delete fOutputList;
311 delete fPID;
312 delete fCFM;
313 delete fPIDqa;
314 delete fTrackCuts;
315}
316//_________________________________________
317
318void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
319{
320 //Main loop
321 //Called for each event
322
323 // create pointer to event
324 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
325 if (!fESD){
326 printf("ERROR: fESD not available\n");
327 return;
328 }
329
330 fVevent = dynamic_cast<AliVEvent*>(InputEvent());
331 if(!fVevent){
332 printf("ERROR: fVEvent not available\n");
333 return;
334 }
335
336 if(!fCuts){
337 AliError("HFE cuts not available");
338 return;
339 }
340
341 if(!fPID->IsInitialized()){
342 // Initialize PID with the given run number
343 AliWarning("PID not initialised, get from Run no");
344 if(fIsAOD)fPID->InitializePID(fAOD->GetRunNumber());
345 if(!fIsAOD) fPID->InitializePID(fESD->GetRunNumber());
346 }
347
348 Bool_t SelColl = kTRUE;
349 if(GetCollisionCandidates()==AliVEvent::kAny)
350 {
351 SelColl = kFALSE;
352 TString firedTrigger;
353 firedTrigger = fESD->GetFiredTriggerClasses();
354 if(firedTrigger.Contains("CVLN_B2-B-NOPF-ALLNOTRD") || firedTrigger.Contains("CVLN_R1-B-NOPF-ALLNOTRD") || firedTrigger.Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))SelColl=kTRUE;
355 if(!SelColl)return;
356 }
357
358 if(fIsMC)fMC = MCEvent();
359 AliStack* stack = NULL;
360 if(fIsMC && fMC) stack = fMC->Stack();
361
362 Int_t fNOtrks =fESD->GetNumberOfTracks();
363 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
364
365 Double_t pVtxZ = -999;
366 pVtxZ = pVtx->GetZ();
367
368 if(TMath::Abs(pVtxZ)>10) return;
369 fNoEvents->Fill(0);
370
371 if(fNOtrks<2) return;
372
373 fpidResponse = fInputHandler->GetPIDResponse();
374 if(!fpidResponse){
375 AliDebug(1, "Using default PID Response");
376 fpidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
377 }
378
379 fPID->SetPIDResponse(fpidResponse);
380
381 fCFM->SetRecEventInfo(fVevent);
382
383 Float_t cent = -1.;
384 Int_t iPt=8, iCent=3, iDeltaphi=4;
385 AliCentrality *centrality = fESD->GetCentrality();
386 cent = centrality->GetCentralityPercentile("V0M");
387 fCent->Fill(cent);
388
389 if (cent>=0 && cent<10) iCent=0;
390 if (cent>=10 && cent<20) iCent=1;
391 if (cent>=20 && cent<40) iCent=2;
392 if (cent<0 || cent>=40) return;
393
394 //Event planes
395
396 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
397 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
398
399 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
400 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
401
402 Double_t evPlaneV0 = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0",fESD,2));
403 if(evPlaneV0 > TMath::Pi()) evPlaneV0 = evPlaneV0 - TMath::Pi();
404 fevPlaneV0[iCent]->Fill(evPlaneV0);
405
406 AliEventplane* esdTPCep = fESD->GetEventplane();
407 TVector2 *standardQ = 0x0;
408 Double_t qx = -999., qy = -999.;
409 standardQ = esdTPCep->GetQVector();
410 if(!standardQ)return;
411
412 qx = standardQ->X();
413 qy = standardQ->Y();
414
415 TVector2 qVectorfortrack;
416 qVectorfortrack.Set(qx,qy);
417 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
418
419 //Event plane resolutions
420
421 // --> 2 subevent method (only for TPC EP)
422
423 TVector2 *qsub1a = esdTPCep->GetQsub1();
424 TVector2 *qsub2a = esdTPCep->GetQsub2();
425 Double_t evPlaneResTPC = -999.;
426 if(qsub1a && qsub2a){
427 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
428 }
429
430 fTPCsubEPres->Fill(evPlaneResTPC);
431
432 // --> 3 event method (V0, V0A, and V0C EP)
433
434 Double_t Qx2pos = 0., Qy2pos = 0., Qx2neg = 0., Qy2neg = 0., Qweight = 1;
435
436 for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
437
438 AliVParticle* vparticle = fVevent->GetTrack(iTracks);
439 if (!vparticle){
440 printf("ERROR: Could not receive track %d\n", iTracks);
441 continue;
442 }
443
444 AliESDtrack *trackEP = dynamic_cast<AliESDtrack*>(vparticle);
445
446 if (!trackEP) {
447 printf("ERROR: Could not receive track %d\n", iTracks);
448 continue;
449 }
450
451 if(TMath::Abs(trackEP->Eta())>0.8 || trackEP->Pt() < 0.15 || trackEP->Pt() > 4) continue;
452
453 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, trackEP)) continue;
454
455 if (trackEP->Pt() < 2) Qweight = trackEP->Pt()/2;
456 if (trackEP->Pt() >= 2) Qweight = 1;
457
458
459 if(trackEP->Eta()>0 && trackEP->Eta()<0.8){
460 Qx2pos += Qweight*TMath::Cos(2*trackEP->Phi());
461 Qy2pos += Qweight*TMath::Sin(2*trackEP->Phi());
462 }
463 if(trackEP->Eta()<0 && trackEP->Eta()>-0.8){
464 Qx2neg += Qweight*TMath::Cos(2*trackEP->Phi());
465 Qy2neg += Qweight*TMath::Sin(2*trackEP->Phi());
466 }
467 }//track loop only for EP
468
469 Double_t evPlaneTPCneg = TMath::ATan2(Qy2neg, Qx2neg)/2;
470 Double_t evPlaneTPCpos = TMath::ATan2(Qy2pos, Qx2pos)/2;
471
472 Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0,evPlaneTPCpos),
473 GetCos2DeltaPhi(evPlaneV0,evPlaneTPCneg),
474 GetCos2DeltaPhi(evPlaneTPCpos,evPlaneTPCneg),cent};
475 fEPres->Fill(evPlaneRes);
476
477 // Selection of primary pi0 and eta in MC to compute the weight
478 if(fIsMC && fMC && stack){
479 Int_t nParticles = stack->GetNtrack();
480 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
481 TParticle* particle = stack->Particle(iParticle);
482 int fPDG = particle->GetPdgCode();
483 double pTMC = particle->Pt();
484
485 Double_t etaMC = particle->Eta();
486 if (TMath::Abs(etaMC)>1.2)continue;
487
488 Bool_t isMotherPrimary = IsPi0EtaPrimary(particle,stack);
489 Bool_t isFromLMdecay = IsPi0EtaFromLMdecay(particle,stack);
490 Bool_t isFromHFdecay = IsPi0EtaFromHFdecay(particle,stack);
491
492 if (!isMotherPrimary || isFromHFdecay || isFromLMdecay) continue;
493
494 if(fPDG==111) fPi0Pt[iCent]->Fill(pTMC); //pi0
495 if(fPDG==221) fEtaPt[iCent]->Fill(pTMC); //eta
496 }
497 }//MC
498
499 Double_t ptRange[9] = {1.5,2,2.5,3,4,6,8,10,13};
500 Double_t deltaPhiRange[4];
501 for(Int_t j=0;j<4;j++){
502 deltaPhiRange[j] = j*(TMath::Pi()/4);
503 }
504
505 // Track loop
506 for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
507
508 AliVParticle* vparticle = fVevent->GetTrack(iTracks);
509 if (!vparticle){
510 printf("ERROR: Could not receive track %d\n", iTracks);
511 continue;
512 }
513
514 AliESDtrack *track = dynamic_cast<AliESDtrack*>(vparticle);
515
516 if(!track) continue;
517
518 if (TMath::Abs(track->Eta())>0.7) continue;
519
520 fTrackPtBefTrkCuts->Fill(track->Pt());
521
522 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
523
524 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
525 if(track->GetKinkIndex(0) != 0) continue;
526 }
527
528 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
529
530 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
531
532 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
533
534 fTrackPtAftTrkCuts->Fill(track->Pt());
535
536 Double_t clsE=-9.,p=-99.,EovP=-99.,pt=-99.,dEdx=-99.,fTPCnSigma=9.,phi=-9.,m02=-9.,m20=-9.,fEMCalnSigma=9.,dphi=9.,cosdphi=9.;
537
538 pt = track->Pt();
539 if(pt<1.5) continue;
540 fTrkpt->Fill(pt);
541 for(Int_t i=0;i<8;i++) {
542 if (pt>=ptRange[i] && pt<ptRange[i+1]){
543 iPt=i;
544 continue;
545 }
546 }
547
548 Int_t clsId = track->GetEMCALcluster();
549 if (clsId>0){
550 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
551 if(cluster && cluster->IsEMCAL()){
552 clsE = cluster->E();
553 m20 = cluster->GetM20();
554 m02 = cluster->GetM02();
555 }
556 }
557
558 p = track->P();
559 phi = track->Phi();
560 dEdx = track->GetTPCsignal();
561 EovP = clsE/p;
562 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
563 fEMCalnSigma = GetSigmaEMCal(EovP, pt, cent);
564 fdEdxBef->Fill(p,dEdx);
565 fTPCnsigma->Fill(p,fTPCnSigma);
566
567
568 dphi = GetDeltaPhi(phi,evPlaneV0);
569 cosdphi = GetCos2DeltaPhi(phi,evPlaneV0);
570 for(Int_t i=0;i<3;i++) {
571 if (dphi>=deltaPhiRange[i] && dphi<deltaPhiRange[i+1]){
572 iDeltaphi=i;
573 continue;
574 }
575 }
576
577 Bool_t fFlagPhotonicElec = kFALSE;
578 Bool_t fFlagPhotonicElecBCG = kFALSE;
579 Double_t weight = 1.;
580
581 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG,weight,iCent);
582
583 Int_t partPDG = -99;
584 Double_t partPt = -99.;
585 Bool_t MChijing;
586
587 if(fIsMC && fMC && stack){
588 if(fTPCnSigma < -1 || fTPCnSigma > 3) continue;
589 Int_t label = track->GetLabel();
590 if(label!=0){
591 TParticle *particle = stack->Particle(label);
592 if(particle){
593 partPDG = particle->GetPdgCode();
594 partPt = particle->Pt();
595 if (TMath::Abs(partPDG)!=11) continue;
596
597 MChijing = fMC->IsFromBGEvent(label);
598 int iHijing = 1;
599 if(!MChijing) iHijing = 0; // 0 if enhanced sample
600
601 Bool_t pi0Decay = IsElectronFromPi0(particle,stack,weight,cent);
602 Bool_t etaDecay = IsElectronFromEta(particle,stack,weight,cent);
603
604 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG,weight,iCent);
605
606 if (pi0Decay){
607 fTot_pi0e->Fill(partPt,weight);
608 if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
609 if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
610 }
611 if (etaDecay){
612 fTot_etae->Fill(partPt,weight);
613 if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
614 if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
615 }
616 }// end particle
617 }// end label
618 }//end MC
619
620 if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
621 Int_t pidpassed = 1;
622
623 //--- track accepted
624 AliHFEpidObject hfetrack;
625 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
626 hfetrack.SetRecTrack(track);
627 hfetrack.SetPbPb();
628 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
629
630 if (m20>0.02 && m02>0.02){
631 Double_t corr[7]={(Double_t)iCent,(Double_t)iPt,fTPCnSigma,fEMCalnSigma,m02,dphi,cosdphi};
632 fCorr->Fill(corr);
633 }
634
635 fChargPartV2[iCent]->Fill(iPt,cosdphi);
636 if (clsE>0) fMtcPartV2[iCent]->Fill(iPt,cosdphi);
637
638 if (pidpassed==0) continue;
639
640 if (fTPCnSigma>=-5 && fTPCnSigma<-3.2) fEoverPbcg[iCent][iPt][iDeltaphi]->Fill(EovP);
641 if (fTPCnSigma>=-0.5 && fTPCnSigma<3) feTPCV2[iCent]->Fill(iPt,cosdphi);
642 if (fTPCnSigma>=-0.5 && fTPCnSigma<3 && fEMCalnSigma>-1 && fEMCalnSigma<3) feV2[iCent]->Fill(iPt,cosdphi);
643
644 fTrkEovPAft->Fill(pt,EovP);
645 fdEdxAft->Fill(p,dEdx);
646
647 if(fFlagPhotonicElec){
648 fPhotoElecPt->Fill(pt);
649 }
650
651 if (!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
652
653 if (m20>0.02 && m02>0.02 && m02<0.27 && fTPCnSigma>-1 && fTPCnSigma<3){
654 fEoverPsig[iCent][iPt][iDeltaphi]->Fill(EovP);
655 if (fFlagPhotonicElec) fEoverPuls[iCent][iPt][iDeltaphi]->Fill(EovP);
656 if (fFlagPhotonicElecBCG) fEoverPls[iCent][iPt][iDeltaphi]->Fill(EovP);
657 }
658
659 }//end of track loop
660 PostData(1, fOutputList);
661}
662//_________________________________________
663void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
664{
665 //--- Check MC
666 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
667 fIsMC = kTRUE;
668 printf("+++++ MC Data available");
669 }
670 //--------Initialize PID
671 fPID->SetHasMCData(fIsMC);
672
673 if(!fPID->GetNumberOfPIDdetectors())
674 {
675 fPID->AddDetector("TPC", 0);
676 fPID->AddDetector("EMCAL", 1);
677 }
678
679 fPID->SortDetectors();
680 fPIDqa = new AliHFEpidQAmanager();
681 fPIDqa->Initialize(fPID);
682
683 //--------Initialize correction Framework and Cuts
684 fCFM = new AliCFManager;
685 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
686 fCFM->SetNStepParticle(kNcutSteps);
687 for(Int_t istep = 0; istep < kNcutSteps; istep++)
688 fCFM->SetParticleCutsList(istep, NULL);
689
690 if(!fCuts){
691 AliWarning("Cuts not available. Default cuts will be used");
692 fCuts = new AliHFEcuts;
693 fCuts->CreateStandardCuts();
694 }
695 fCuts->Initialize(fCFM);
696
697 //---------Output Tlist
698 fOutputList = new TList();
699 fOutputList->SetOwner();
700 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
701
702 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
703 fOutputList->Add(fNoEvents);
704
705 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
706 fOutputList->Add(fTrkpt);
707
708 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
709 fOutputList->Add(fTrackPtBefTrkCuts);
710
711 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
712 fOutputList->Add(fTrackPtAftTrkCuts);
713
714 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
715 fOutputList->Add(fTPCnsigma);
716
717 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
718 fOutputList->Add(fTrkEovPBef);
719
720 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
721 fOutputList->Add(fTrkEovPAft);
722
723 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
724 fOutputList->Add(fdEdxBef);
725
726 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
727 fOutputList->Add(fdEdxAft);
728
729 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
730 fOutputList->Add(fPhotoElecPt);
731
732 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
733 fOutputList->Add(fSemiInclElecPt);
734
735 fCent = new TH1F("fCent","Centrality",100,0,100) ;
736 fOutputList->Add(fCent);
737
738 fTPCsubEPres = new TH1F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1);
739 fOutputList->Add(fTPCsubEPres);
740
741 Int_t binsv1[4]={100,100,100,90}; // V0-TPCpos, V0-TPCneg, TPCpos-TPCneg, cent
742 Double_t xminv1[4]={-1,-1,-1,0};
743 Double_t xmaxv1[4]={1,1,1,90};
744 fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
745 fOutputList->Add(fEPres);
746
747 //iCent,iPt,fTPCnSigma,fEMCalnSigma,m02,dphi,cosdphi
748 Int_t binsv2[7]={3,8,100,100,100,120,100};
749 Double_t xminv2[7]={0,0,-5,-5,0,0,-1};
750 Double_t xmaxv2[7]={3,8,5,5,2,TMath::TwoPi(),1};
751 fCorr = new THnSparseD ("fCorr","Correlations",7,binsv2,xminv2,xmaxv2);
752 fOutputList->Add(fCorr);
753
754 for(Int_t i=0; i<3; i++) {
755 fevPlaneV0[i] = new TH1F(Form("fevPlaneV0%d",i),"V0 EP",100,0,TMath::Pi());
756 fOutputList->Add(fevPlaneV0[i]);
757
758 feTPCV2[i] = new TH2F(Form("feTPCV2%d",i), "", 8,0,8,100,-1,1);
759 fOutputList->Add(feTPCV2[i]);
760
761 feV2[i] = new TH2F(Form("feV2%d",i), "", 8,0,8,100,-1,1);
762 fOutputList->Add(feV2[i]);
763
764 fChargPartV2[i] = new TH2F(Form("fChargPartV2%d",i), "", 8,0,8,100,-1,1);
765 fOutputList->Add(fChargPartV2[i]);
766
767 fMtcPartV2[i] = new TH2F(Form("fMtcPartV2%d",i), "", 8,0,8,100,-1,1);
768 fOutputList->Add(fMtcPartV2[i]);
769
770 fInvmassLS[i] = new TH2F(Form("fInvmassLS%d",i), "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5,100,0,50);
771 fOutputList->Add(fInvmassLS[i]);
772
773 fInvmassULS[i] = new TH2F(Form("fInvmassULS%d",i), "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5,100,0,50);
774 fOutputList->Add(fInvmassULS[i]);
775
776 fOpeningAngleLS[i] = new TH2F(Form("fOpeningAngleLS%d",i),"Opening angle for LS pairs",100,0,1,100,0,50);
777 fOutputList->Add(fOpeningAngleLS[i]);
778
779 fOpeningAngleULS[i] = new TH2F(Form("fOpeningAngleULS%d",i),"Opening angle for ULS pairs",100,0,1,100,0,50);
780 fOutputList->Add(fOpeningAngleULS[i]);
781
782 fPi0Pt[i] = new TH1F(Form("fPi0Pt%d",i), "Pi0 weight",100,0,50);
783 fOutputList->Add(fPi0Pt[i]);
784
785 fEtaPt[i] = new TH1F(Form("fEtaPt%d",i), "Eta weight",100,0,50);
786 fOutputList->Add(fEtaPt[i]);
787 }
788
789 for(Int_t i=0; i<3; i++) {
790 for(Int_t j=0; j<8; j++) {
791 for(Int_t k=0; k<4; k++) {
792 fEoverPsig[i][j][k] = new TH1F(Form("fEoverPsig%d%d%d",i,j,k), "",100,0,3);
793 fOutputList->Add(fEoverPsig[i][j][k]);
794
795 fEoverPuls[i][j][k] = new TH1F(Form("fEoverPuls%d%d%d",i,j,k), "",100,0,3);
796 fOutputList->Add(fEoverPuls[i][j][k]);
797
798 fEoverPls[i][j][k] = new TH1F(Form("fEoverPls%d%d%d",i,j,k), "",100,0,3);
799 fOutputList->Add(fEoverPls[i][j][k]);
800
801 fEoverPbcg[i][j][k] = new TH1F(Form("fEoverPbcg%d%d%d",i,j,k), "",100,0,3);
802 fOutputList->Add(fEoverPbcg[i][j][k]);
803 }
804 }
805 }
806
807 fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
808 fOutputList->Add(fD0_e);
809
810 for(Int_t k = 0; k < 6; k++) {
811
812 TString De_name = Form("fDe%d",k);
813 TString D0e_name = Form("fD0e%d",k);
814 TString Dpluse_name = Form("fDpluse%d",k);
815 TString Dminuse_name = Form("fDminuse%d",k);
816
817 fDe[k] = new TH1F((const char*)De_name,"",100,0,50);
818 fD0e[k] = new TH1F((const char*)D0e_name,"",100,0,50);
819 fDpluse[k] = new TH1F((const char*)Dpluse_name,"",100,0,50);
820 fDminuse[k] = new TH1F((const char*)Dminuse_name,"",100,0,50);
821
822 fOutputList->Add(fDe[k]);
823 fOutputList->Add(fD0e[k]);
824 fOutputList->Add(fDpluse[k]);
825 fOutputList->Add(fDminuse[k]);
826
827 }
828
829 int nbin_v2 = 8;
830 double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
831
832 fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
833 fOutputList->Add(fTot_pi0e);
834
835 fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
836 fOutputList->Add(fPhot_pi0e);
837
838 fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
839 fOutputList->Add(fPhotBCG_pi0e);
840
841 fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
842 fOutputList->Add(fTot_etae);
843
844 fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
845 fOutputList->Add(fPhot_etae);
846
847 fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
848 fOutputList->Add(fPhotBCG_etae);
849
850 fInvMass = new TH1F("fInvMass","",200,0,0.3);
851 fOutputList->Add(fInvMass);
852
853 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
854 fOutputList->Add(fInvMassBack);
855
856 fDCA = new TH1F("fDCA","",200,0,1);
857 fOutputList->Add(fDCA);
858
859 fDCABack = new TH1F("fDCABack","",200,0,1);
860 fOutputList->Add(fDCABack);
861
862 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
863 fOutputList->Add(fOpAngle);
864
865 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
866 fOutputList->Add(fOpAngleBack);
867
868 PostData(1,fOutputList);
869}
870
871//________________________________________________________________________
872void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
873{
874 // Info("Terminate");
875 AliAnalysisTaskSE::Terminate();
876}
877
878//________________________________________________________________________
879Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
880{
881 // Check single track cuts for a given cut step
882 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
883 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
884 return kTRUE;
885}
886//_________________________________________
887Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
888{
889 //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
890 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
891 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
892 Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
893
894 return cos2DeltaPhi;
895}
896//_________________________________________
897Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
898{
899 //Get phi-psi_EP
900 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
901 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
902
903 return dPhi;
904}
905//_________________________________________
906Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT, Float_t cent) const
907{
908 //Get Pi0 weight
909 double weight = 1.0;
910 return weight;
911}
912//_________________________________________
913Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT, Float_t cent) const
914{
915 //Get eta weight
916 double weight = 1.0;
917 return weight;
918}
919//_________________________________________
920Double_t AliAnalysisTaskFlowTPCEMCalEP::GetSigmaEMCal(Double_t EoverP, Double_t pt, Float_t cent) const
921{
922 //Get sigma for EMCal PID
923 Double_t NumberOfSigmasEMCal = 99.;
924 Double_t ptRange[9] = {1.5,2,2.5,3,4,6,8,10,13};
925
926 if (cent>=0 && cent<10){
927 Double_t mean[8]={0.953184,0.957259,0.97798,0.9875,1.03409,1.06257,1.02776,1.04338};
928 Double_t sigma[8]={0.130003,0.113493,0.092966,0.0836828,0.101804,0.0893414,0.0950752,0.050427};
929 for(Int_t i=0;i<8;i++) {
930 if (pt>=ptRange[i] && pt<ptRange[i+1]){
931 NumberOfSigmasEMCal = (mean[i]-EoverP)/sigma[i];
932 continue;
933 }
934 }
935 }
936 if (cent>=10 && cent<20){
937 Double_t mean[8]={0.96905,0.952985,0.96871,0.983934,1.00047,0.988736,1.02101,1.04557};
938 Double_t sigma[8]={0.0978103,0.103215,0.0958494,0.0797962,0.0719482,0.0672677,0.0754882,0.0461192};
939 for(Int_t i=0;i<8;i++) {
940 if (pt>=ptRange[i] && pt<ptRange[i+1]){
941 NumberOfSigmasEMCal = (mean[i]-EoverP)/sigma[i];
942 continue;
943 }
944 }
945 }
946 if (cent>=20 && cent<40){
947 Double_t mean[8]={0.947362,0.951933,0.959288,0.977004,0.984502,1.02004,1.00489,0.986696};
948 Double_t sigma[8]={0.100127,0.0887731,0.0842077,0.0787335,0.0804325,0.0652376,0.0766669,0.0597849};
949 for(Int_t i=0;i<8;i++) {
950 if (pt>=ptRange[i] && pt<ptRange[i+1]){
951 NumberOfSigmasEMCal = (mean[i]-EoverP)/sigma[i];
952 continue;
953 }
954 }
955 }
956
957
958
959 return NumberOfSigmasEMCal;
960}
961//_________________________________________
962Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsPi0EtaFromHFdecay(TParticle *particle, AliStack* stack)
963{
964 // Check if the mother comes from heavy-flavour decays
965
966 Bool_t isHFdecay = kFALSE;
967 Int_t partPDG = particle->GetPdgCode();
968
969 if (TMath::Abs(partPDG)!=111 || TMath::Abs(partPDG)!=221) return isHFdecay; // particle is not pi0 or eta
970
971 Int_t idMother = particle->GetFirstMother();
972 if (idMother>0){
973 TParticle *mother = stack->Particle(idMother);
974 Int_t motherPDG = mother->GetPdgCode();
975
976 // c decays
977 if((TMath::Abs(motherPDG)==411) || (TMath::Abs(motherPDG)==421) || (TMath::Abs(motherPDG)==431) || (TMath::Abs(motherPDG)==4122) || (TMath::Abs(motherPDG)==4132) || (TMath::Abs(motherPDG)==4232) || (TMath::Abs(motherPDG)==43320)) isHFdecay = kTRUE;
978
979 // b decays
980 if((TMath::Abs(motherPDG)==511) || (TMath::Abs(motherPDG)==521) || (TMath::Abs(motherPDG)==531) || (TMath::Abs(motherPDG)==5122) || (TMath::Abs(motherPDG)==5132) || (TMath::Abs(motherPDG)==5232) || (TMath::Abs(motherPDG)==53320)) isHFdecay = kTRUE;
981 }
982
983 return isHFdecay;
984}
985//_________________________________________
986Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsPi0EtaFromLMdecay(TParticle *particle, AliStack* stack)
987{
988 // Check if the mother comes from light-meson decays
989
990 Bool_t isLMdecay = kFALSE;
991 Int_t partPDG = particle->GetPdgCode();
992
993 if (TMath::Abs(partPDG)!=111 || TMath::Abs(partPDG)!=221) return isLMdecay; // particle is not pi0 or eta
994
995 Int_t idMother = particle->GetFirstMother();
996 if (idMother>0){
997 TParticle *mother = stack->Particle(idMother);
998 Int_t motherPDG = mother->GetPdgCode();
999
1000 if(motherPDG == 111 || motherPDG == 221 || motherPDG==223 || motherPDG==333 || motherPDG==331 || (TMath::Abs(motherPDG)==113) || (TMath::Abs(motherPDG)==213) || (TMath::Abs(motherPDG)==313) || (TMath::Abs(motherPDG)==323)) isLMdecay = kTRUE;
1001 }
1002
1003 return isLMdecay;
1004}
1005//_________________________________________
1006Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsPi0EtaPrimary(TParticle *particle, AliStack* stack)
1007{
1008 // Check if the pi0 or eta are primary
1009
1010 Bool_t isprimary = kFALSE;
1011 Int_t partPDG = particle->GetPdgCode();
1012
1013 if (TMath::Abs(partPDG)!=111 || TMath::Abs(partPDG)!=221) return isprimary; // particle is not pi0 or eta
1014
1015
1016 Bool_t pi0etaprimary = particle->IsPrimary();
1017 Int_t idMother = particle->GetFirstMother();
1018
1019 if (pi0etaprimary && idMother<=0) isprimary = kTRUE;
1020
1021 return isprimary;
1022}
1023//_________________________________________
1024Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsElectronFromPi0(TParticle *particle, AliStack* stack, Double_t &weight, Float_t cent)
1025{
1026 // Check if electron comes from primary pi0 not from light-meson and heavy-flavour decays
1027
1028 Bool_t isPi0Decay = kFALSE;
1029 Int_t partPDG = particle->GetPdgCode();
1030
1031 if (TMath::Abs(partPDG)!=11) return isPi0Decay; // particle is not electron
1032
1033 Int_t idMother = particle->GetFirstMother();
1034 if (idMother>0){
1035 TParticle *mother = stack->Particle(idMother);
1036 Int_t motherPDG = mother->GetPdgCode();
1037 Double_t motherPt = mother->Pt();
1038
1039 Bool_t isMotherPi0primary = IsPi0EtaPrimary(mother,stack);
1040 Bool_t isMotherPi0fromHF = IsPi0EtaFromHFdecay(mother,stack);
1041 Bool_t isMotherPi0fromLM = IsPi0EtaFromLMdecay(mother,stack);
1042
1043 if (motherPDG==111 && (isMotherPi0primary || (!isMotherPi0fromHF && !isMotherPi0fromLM))){ // pi0 -> e
1044 isPi0Decay = kTRUE;
1045 weight = GetPi0weight(motherPt,cent);
1046 }
1047
1048 Int_t idSecondMother = particle->GetSecondMother();
1049 if (motherPDG==22 && idSecondMother>0){
1050 TParticle *secondMother = stack->Particle(idSecondMother);
1051 Int_t secondMotherPDG = secondMother->GetPdgCode();
1052 Double_t secondMotherPt = secondMother->Pt();
1053
1054 Bool_t isSecondMotherPi0primary = IsPi0EtaPrimary(secondMother,stack);
1055 Bool_t isSecondMotherPi0fromHF = IsPi0EtaFromHFdecay(secondMother,stack);
1056 Bool_t isSecondMotherPi0fromLM = IsPi0EtaFromLMdecay(secondMother,stack);
1057
1058 if (secondMotherPDG==111 && (isSecondMotherPi0primary || (!isSecondMotherPi0fromHF && !isSecondMotherPi0fromLM))){ //pi0 -> gamma -> e
1059 isPi0Decay = kTRUE;
1060 weight = GetPi0weight(secondMotherPt,cent);
1061 }
1062 }
1063 }
1064 return isPi0Decay;
1065}
1066//_________________________________________
1067Bool_t AliAnalysisTaskFlowTPCEMCalEP::IsElectronFromEta(TParticle *particle, AliStack* stack, Double_t &weight, Float_t cent)
1068{
1069 // Check if electron comes from primary eta not from light-meson and heavy-flavour decays
1070
1071 Bool_t isEtaDecay = kFALSE;
1072 Int_t partPDG = particle->GetPdgCode();
1073
1074 if (TMath::Abs(partPDG)!=11) return isEtaDecay; // particle is not electron
1075
1076 Int_t idMother = particle->GetFirstMother();
1077 if (idMother>0){
1078 TParticle *mother = stack->Particle(idMother);
1079 Int_t motherPDG = mother->GetPdgCode();
1080 Double_t motherPt = mother->Pt();
1081
1082 Bool_t isMotherEtaprimary = IsPi0EtaPrimary(mother,stack);
1083 Bool_t isMotherEtafromHF = IsPi0EtaFromHFdecay(mother,stack);
1084 Bool_t isMotherEtafromLM = IsPi0EtaFromLMdecay(mother,stack);
1085
1086 if (motherPDG==221 && (isMotherEtaprimary || (!isMotherEtafromHF && !isMotherEtafromLM))){ //primary eta -> e
1087 isEtaDecay = kTRUE;
1088 weight = GetEtaweight(motherPt,cent);
1089 }
1090
1091 Int_t idSecondMother = mother->GetFirstMother();
1092 if ((motherPDG==22 || motherPDG==111) && idSecondMother>0){
1093 TParticle *secondMother = stack->Particle(idSecondMother);
1094 Int_t secondMotherPDG = secondMother->GetPdgCode();
1095 Double_t secondMotherPt = secondMother->Pt();
1096
1097 Bool_t isSecondMotherEtaprimary = IsPi0EtaPrimary(secondMother,stack);
1098 Bool_t isSecondMotherEtafromHF = IsPi0EtaFromHFdecay(secondMother,stack);
1099 Bool_t isSecondMotherEtafromLM = IsPi0EtaFromLMdecay(secondMother,stack);
1100
1101 if (secondMotherPDG==221 && (isSecondMotherEtaprimary || (!isSecondMotherEtafromHF && !isSecondMotherEtafromLM))){ //eta -> pi0/g-> e
1102 isEtaDecay = kTRUE;
1103 weight = GetEtaweight(secondMotherPt,cent);
1104 }
1105 Int_t idThirdMother = secondMother->GetFirstMother();
1106 if (idThirdMother>0){
1107 TParticle *thirdMother = stack->Particle(idThirdMother);
1108 Int_t thirdMotherPDG = thirdMother->GetPdgCode();
1109 Double_t thirdMotherPt = thirdMother->Pt();
1110
1111 Bool_t isThirdMotherEtaprimary = IsPi0EtaPrimary(thirdMother,stack);
1112 Bool_t isThirdMotherEtafromHF = IsPi0EtaFromHFdecay(thirdMother,stack);
1113 Bool_t isThirdMotherEtafromLM = IsPi0EtaFromLMdecay(thirdMother,stack);
1114
1115 if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221 && (isThirdMotherEtaprimary || (!isThirdMotherEtafromHF && !isThirdMotherEtafromLM))){//p eta->pi0->g-> e
1116 isEtaDecay = kTRUE;
1117 weight = GetEtaweight(thirdMotherPt,cent);
1118 }
1119 }
1120 }
1121 }
1122 return isEtaDecay;
1123}
1124//_________________________________________
1125void AliAnalysisTaskFlowTPCEMCalEP::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG,Double_t weight, Int_t iCent)
1126{
1127 //Identify non-heavy flavour electrons using Invariant mass method
1128
1129 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1130 fTrackCuts->SetRequireTPCRefit(kTRUE);
1131 fTrackCuts->SetRequireITSRefit(kTRUE);
1132 fTrackCuts->SetEtaRange(-0.9,0.9);
1133 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1134 fTrackCuts->SetMaxChi2PerClusterTPC(4);
1135 fTrackCuts->SetMinNClustersTPC(80);
1136 fTrackCuts->SetMaxDCAToVertexZ(3.2);
1137 fTrackCuts->SetMaxDCAToVertexXY(2.4);
1138 fTrackCuts->SetDCAToVertex2D(kTRUE);
1139
1140 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1141
1142 Bool_t flagPhotonicElec = kFALSE;
1143 Bool_t flagPhotonicElecBCG = kFALSE;
1144
1145 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1146
1147 if(jTracks==iTracks) continue;
1148
1149 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1150 if (!trackAsso) {
1151 printf("ERROR: Could not receive track %d\n", jTracks);
1152 continue;
1153 }
1154
1155 Double_t pt=-999., ptAsso=-999., nTPCsigmaAsso=-999.;
1156 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1157 Double_t openingAngle = -999., mass=999., width = -999;
1158 Int_t chargeAsso = 0, charge = 0;
1159
1160 nTPCsigmaAsso = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1161 pt = track->Pt();
1162 ptAsso = trackAsso->Pt();
1163 chargeAsso = trackAsso->Charge();
1164 charge = track->Charge();
1165
1166 if(ptAsso <0.5) continue;
1167 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1168 if(TMath::Abs(nTPCsigmaAsso)>3) continue;
1169
1170 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1171 if(charge>0) fPDGe1 = -11;
1172 if(chargeAsso>0) fPDGe2 = -11;
1173
1174 if(charge == chargeAsso) fFlagLS = kTRUE;
1175 if(charge != chargeAsso) fFlagULS = kTRUE;
1176
1177 AliKFParticle ge1(*track, fPDGe1);
1178 AliKFParticle ge2(*trackAsso, fPDGe2);
1179 AliKFParticle recg(ge1, ge2);
1180
1181 if(recg.GetNDF()<1) continue;
1182 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1183 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1184
1185 if(fSetMassConstraint && pVtx) {
1186 AliKFVertex primV(*pVtx);
1187 primV += recg;
1188 primV -= ge1;
1189 primV -= ge2;
1190 recg.SetProductionVertex(primV);
1191 recg.SetMassConstraint(0,0.0001);
1192 }
1193
1194 openingAngle = ge1.GetAngle(ge2);
1195
1196 if(fFlagLS) fOpeningAngleLS[iCent]->Fill(openingAngle,pt);
1197 if(fFlagULS) fOpeningAngleULS[iCent]->Fill(openingAngle,pt);
1198
1199 if(openingAngle > fOpeningAngleCut) continue;
1200
1201 recg.GetMass(mass,width);
1202
1203 if(fFlagLS) fInvmassLS[iCent]->Fill(mass,pt,weight);
1204 if(fFlagULS) fInvmassULS[iCent]->Fill(mass,pt,weight);
1205
1206 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
1207 if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
1208
1209 }
1210 fFlagPhotonicElec = flagPhotonicElec;
1211 fFlagPhotonicElecBCG = flagPhotonicElecBCG;
1212
1213}
1214//_________________________________________
1215void AliAnalysisTaskFlowTPCEMCalEP::InitParameters()
1216{
1217 // Init parameters
1218
1219 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1220 fTrackCuts->SetRequireTPCRefit(kTRUE);
1221 fTrackCuts->SetRequireITSRefit(kTRUE);
1222 fTrackCuts->SetEtaRange(-0.7,0.7);
1223 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1224 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1225 fTrackCuts->SetMinNClustersTPC(100);
1226 fTrackCuts->SetPtRange(0.5,100);
1227
1228 fNonHFE->SetAODanalysis(kFALSE);
1229 fNonHFE->SetInvariantMassCut(fInvmassCut);
1230 fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
1231 fNonHFE->SetChi2OverNDFCut(fChi2Cut);
1232 fNonHFE->SetAlgorithm(fnonHFEalgorithm); //KF or DCA
1233 if (fnonHFEalgorithm == "DCA") fNonHFE->SetDCACut(fDCAcut);
1234 fNonHFE->SetTrackCuts(-3.5,3.5,fTrackCuts);
1235
1236}
1237