]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalEP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[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)
d6d30be7 105 ,fSetMassConstraint(kFALSE)
1d08b6e8 106 ,fVz(0.0)
107 ,fCFM(0)
108 ,fPID(0)
109 ,fPIDqa(0)
d6d30be7 110 ,fOpeningAngleCut(1000.)
0f3a6a57 111 ,fInvmassCut(0.05)
112 ,fChi2Cut(3.5)
113 ,fDCAcut(999)
0f3a6a57 114 ,fnonHFEalgorithm("KF")
1d08b6e8 115 ,fNoEvents(0)
116 ,fTrkpt(0)
117 ,fTrkEovPBef(0)
118 ,fTrkEovPAft(0)
119 ,fdEdxBef(0)
120 ,fdEdxAft(0)
1d08b6e8 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)
937b2ff2 130 ,fD0_e(0)
1c424962 131 ,fTot_pi0e(0)
132 ,fPhot_pi0e(0)
133 ,fPhotBCG_pi0e(0)
134 ,fTot_etae(0)
135 ,fPhot_etae(0)
136 ,fPhotBCG_etae(0)
0f3a6a57 137 ,fInvMass(0)
138 ,fInvMassBack(0)
139 ,fDCA(0)
140 ,fDCABack(0)
141 ,fOpAngle(0)
142 ,fOpAngleBack(0)
1d08b6e8 143{
144 //Named constructor
937b2ff2 145
d6d30be7 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
937b2ff2 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
1d08b6e8 179 fPID = new AliHFEpid("hfePid");
180 fTrackCuts = new AliESDtrackCuts();
0f3a6a57 181 fNonHFE = new AliSelectNonHFE();
182
183 InitParameters();
1d08b6e8 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//________________________________________________________________________
914042c2 196AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 197 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
198 ,fESD(0)
0f3a6a57 199 ,fAOD(0)
200 ,fVevent(0)
201 ,fpidResponse(0)
3a94d6fe 202 ,fMC(0)
1d08b6e8 203 ,fOutputList(0)
204 ,fTrackCuts(0)
205 ,fCuts(0)
0f3a6a57 206 ,fNonHFE(0)
1d08b6e8 207 ,fIdentifiedAsOutInz(kFALSE)
208 ,fPassTheEventCut(kFALSE)
209 ,fRejectKinkMother(kFALSE)
3a94d6fe 210 ,fIsMC(kFALSE)
0f3a6a57 211 ,fIsAOD(kFALSE)
d6d30be7 212 ,fSetMassConstraint(kFALSE)
1d08b6e8 213 ,fVz(0.0)
214 ,fCFM(0)
215 ,fPID(0)
216 ,fPIDqa(0)
d6d30be7 217 ,fOpeningAngleCut(1000.)
0f3a6a57 218 ,fInvmassCut(0.05)
219 ,fChi2Cut(3.5)
220 ,fDCAcut(999)
0f3a6a57 221 ,fnonHFEalgorithm("KF")
1d08b6e8 222 ,fNoEvents(0)
223 ,fTrkpt(0)
224 ,fTrkEovPBef(0)
225 ,fTrkEovPAft(0)
226 ,fdEdxBef(0)
227 ,fdEdxAft(0)
1d08b6e8 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)
937b2ff2 237 ,fD0_e(0)
1c424962 238 ,fTot_pi0e(0)
239 ,fPhot_pi0e(0)
240 ,fPhotBCG_pi0e(0)
241 ,fTot_etae(0)
242 ,fPhot_etae(0)
243 ,fPhotBCG_etae(0)
0f3a6a57 244 ,fInvMass(0)
245 ,fInvMassBack(0)
246 ,fDCA(0)
247 ,fDCABack(0)
248 ,fOpAngle(0)
249 ,fOpAngleBack(0)
1d08b6e8 250{
0f3a6a57 251
252 //Default constructor
937b2ff2 253
d6d30be7 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
0f3a6a57 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 }
937b2ff2 286
0f3a6a57 287 fPID = new AliHFEpid("hfePid");
288 fTrackCuts = new AliESDtrackCuts();
289 fNonHFE = new AliSelectNonHFE();
1d08b6e8 290
0f3a6a57 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());
1d08b6e8 303}
304//_________________________________________
305
914042c2 306AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 307{
308 //Destructor
309
310 delete fOutputList;
311 delete fPID;
312 delete fCFM;
313 delete fPIDqa;
314 delete fTrackCuts;
315}
316//_________________________________________
317
914042c2 318void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
1d08b6e8 319{
320 //Main loop
321 //Called for each event
0f3a6a57 322
1d08b6e8 323 // create pointer to event
324 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
0f3a6a57 325 if (!fESD){
1d08b6e8 326 printf("ERROR: fESD not available\n");
327 return;
328 }
0f3a6a57 329
330 fVevent = dynamic_cast<AliVEvent*>(InputEvent());
331 if(!fVevent){
332 printf("ERROR: fVEvent not available\n");
333 return;
334 }
335
1d08b6e8 336 if(!fCuts){
337 AliError("HFE cuts not available");
338 return;
339 }
0f3a6a57 340
1d08b6e8 341 if(!fPID->IsInitialized()){
342 // Initialize PID with the given run number
343 AliWarning("PID not initialised, get from Run no");
0f3a6a57 344 if(fIsAOD)fPID->InitializePID(fAOD->GetRunNumber());
345 if(!fIsAOD) fPID->InitializePID(fESD->GetRunNumber());
1d08b6e8 346 }
347
5528bcf6 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
3a94d6fe 358 if(fIsMC)fMC = MCEvent();
359 AliStack* stack = NULL;
360 if(fIsMC && fMC) stack = fMC->Stack();
1d08b6e8 361
0f3a6a57 362 Int_t fNOtrks =fESD->GetNumberOfTracks();
1d08b6e8 363 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
0f3a6a57 364
1d08b6e8 365 Double_t pVtxZ = -999;
366 pVtxZ = pVtx->GetZ();
0f3a6a57 367
1d08b6e8 368 if(TMath::Abs(pVtxZ)>10) return;
369 fNoEvents->Fill(0);
0f3a6a57 370
1d08b6e8 371 if(fNOtrks<2) return;
0f3a6a57 372
373 fpidResponse = fInputHandler->GetPIDResponse();
374 if(!fpidResponse){
1d08b6e8 375 AliDebug(1, "Using default PID Response");
0f3a6a57 376 fpidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
1d08b6e8 377 }
0f3a6a57 378
379 fPID->SetPIDResponse(fpidResponse);
380
381 fCFM->SetRecEventInfo(fVevent);
382
1d08b6e8 383 Float_t cent = -1.;
d6d30be7 384 Int_t iPt=8, iCent=3, iDeltaphi=4;
1d08b6e8 385 AliCentrality *centrality = fESD->GetCentrality();
386 cent = centrality->GetCentralityPercentile("V0M");
387 fCent->Fill(cent);
0f3a6a57 388
d6d30be7 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
1d08b6e8 394 //Event planes
0f3a6a57 395
1d08b6e8 396 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
397 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
0f3a6a57 398
1d08b6e8 399 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
400 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
0f3a6a57 401
402 Double_t evPlaneV0 = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0",fESD,2));
403 if(evPlaneV0 > TMath::Pi()) evPlaneV0 = evPlaneV0 - TMath::Pi();
d6d30be7 404 fevPlaneV0[iCent]->Fill(evPlaneV0);
0f3a6a57 405
1d08b6e8 406 AliEventplane* esdTPCep = fESD->GetEventplane();
44fdf7ac 407 TVector2 *standardQ = 0x0;
1d08b6e8 408 Double_t qx = -999., qy = -999.;
44fdf7ac 409 standardQ = esdTPCep->GetQVector();
ee89e744 410 if(!standardQ)return;
411
412 qx = standardQ->X();
413 qy = standardQ->Y();
0f3a6a57 414
1d08b6e8 415 TVector2 qVectorfortrack;
416 qVectorfortrack.Set(qx,qy);
417 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
0f3a6a57 418
419 //Event plane resolutions
420
421 // --> 2 subevent method (only for TPC EP)
422
1d08b6e8 423 TVector2 *qsub1a = esdTPCep->GetQsub1();
424 TVector2 *qsub2a = esdTPCep->GetQsub2();
425 Double_t evPlaneResTPC = -999.;
0f3a6a57 426 if(qsub1a && qsub2a){
427 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
1d08b6e8 428 }
0f3a6a57 429
430 fTPCsubEPres->Fill(evPlaneResTPC);
431
432 // --> 3 event method (V0, V0A, and V0C EP)
433
5528bcf6 434 Double_t Qx2pos = 0., Qy2pos = 0., Qx2neg = 0., Qy2neg = 0., Qweight = 1;
0f3a6a57 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;
c4de3211 456 if (trackEP->Pt() >= 2) Qweight = 1;
457
0f3a6a57 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
5528bcf6 472 Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0,evPlaneTPCpos),
473 GetCos2DeltaPhi(evPlaneV0,evPlaneTPCneg),
474 GetCos2DeltaPhi(evPlaneTPCpos,evPlaneTPCneg),cent};
1d08b6e8 475 fEPres->Fill(evPlaneRes);
0f3a6a57 476
d6d30be7 477 // Selection of primary pi0 and eta in MC to compute the weight
0f3a6a57 478 if(fIsMC && fMC && stack){
479 Int_t nParticles = stack->GetNtrack();
480 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
6522b3b5 481 TParticle* particle = stack->Particle(iParticle);
482 int fPDG = particle->GetPdgCode();
483 double pTMC = particle->Pt();
d6d30be7 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
0f3a6a57 496 }
d6d30be7 497 }//MC
6522b3b5 498
d6d30be7 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);
6522b3b5 503 }
6522b3b5 504
1d08b6e8 505 // Track loop
0f3a6a57 506 for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
507
508 AliVParticle* vparticle = fVevent->GetTrack(iTracks);
509 if (!vparticle){
1d08b6e8 510 printf("ERROR: Could not receive track %d\n", iTracks);
511 continue;
512 }
0f3a6a57 513
0f3a6a57 514 AliESDtrack *track = dynamic_cast<AliESDtrack*>(vparticle);
515
e5ad2d91 516 if(!track) continue;
517
d6d30be7 518 if (TMath::Abs(track->Eta())>0.7) continue;
0f3a6a57 519
520 fTrackPtBefTrkCuts->Fill(track->Pt());
521
1d08b6e8 522 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
0f3a6a57 523
1d08b6e8 524 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
525 if(track->GetKinkIndex(0) != 0) continue;
526 }
0f3a6a57 527
1d08b6e8 528 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
0f3a6a57 529
1d08b6e8 530 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
0f3a6a57 531
1d08b6e8 532 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
0f3a6a57 533
534 fTrackPtAftTrkCuts->Fill(track->Pt());
535
d6d30be7 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.;
0f3a6a57 537
1d08b6e8 538 pt = track->Pt();
0f3a6a57 539 if(pt<1.5) continue;
1d08b6e8 540 fTrkpt->Fill(pt);
d6d30be7 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 }
0f3a6a57 547
1d08b6e8 548 Int_t clsId = track->GetEMCALcluster();
549 if (clsId>0){
550 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
551 if(cluster && cluster->IsEMCAL()){
0f3a6a57 552 clsE = cluster->E();
553 m20 = cluster->GetM20();
554 m02 = cluster->GetM02();
1d08b6e8 555 }
556 }
0f3a6a57 557
1d08b6e8 558 p = track->P();
d6d30be7 559 phi = track->Phi();
1d08b6e8 560 dEdx = track->GetTPCsignal();
561 EovP = clsE/p;
562 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
d6d30be7 563 fEMCalnSigma = GetSigmaEMCal(EovP, pt, cent);
1d08b6e8 564 fdEdxBef->Fill(p,dEdx);
565 fTPCnsigma->Fill(p,fTPCnSigma);
d6d30be7 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 }
32d22164 575 }
576
06d74a7b 577 Bool_t fFlagPhotonicElec = kFALSE;
578 Bool_t fFlagPhotonicElecBCG = kFALSE;
1c424962 579 Double_t weight = 1.;
535b11a0 580
d6d30be7 581 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG,weight,iCent);
582
583 Int_t partPDG = -99;
584 Double_t partPt = -99.;
585 Bool_t MChijing;
535b11a0 586
32d22164 587 if(fIsMC && fMC && stack){
d6d30be7 588 if(fTPCnSigma < -1 || fTPCnSigma > 3) continue;
06d74a7b 589 Int_t label = track->GetLabel();
d6d30be7 590 if(label!=0){
591 TParticle *particle = stack->Particle(label);
0f3a6a57 592 if(particle){
593 partPDG = particle->GetPdgCode();
594 partPt = particle->Pt();
d6d30be7 595 if (TMath::Abs(partPDG)!=11) continue;
0f3a6a57 596
597 MChijing = fMC->IsFromBGEvent(label);
0f3a6a57 598 int iHijing = 1;
599 if(!MChijing) iHijing = 0; // 0 if enhanced sample
def484dc 600
d6d30be7 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 }
0f3a6a57 616 }// end particle
617 }// end label
618 }//end MC
619
1d08b6e8 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;
b3bd01b3 629
d6d30be7 630 if (m20>0.02 && m02>0.02){
7279a41e 631 Double_t corr[7]={(Double_t)iCent,(Double_t)iPt,fTPCnSigma,fEMCalnSigma,m02,dphi,cosdphi};
d6d30be7 632 fCorr->Fill(corr);
b3bd01b3 633 }
d6d30be7 634
635 fChargPartV2[iCent]->Fill(iPt,cosdphi);
636 if (clsE>0) fMtcPartV2[iCent]->Fill(iPt,cosdphi);
0f3a6a57 637
d6d30be7 638 if (pidpassed==0) continue;
0f3a6a57 639
d6d30be7 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);
0f3a6a57 643
1d08b6e8 644 fTrkEovPAft->Fill(pt,EovP);
645 fdEdxAft->Fill(p,dEdx);
3a94d6fe 646
0f3a6a57 647 if(fFlagPhotonicElec){
0f3a6a57 648 fPhotoElecPt->Fill(pt);
649 }
650
d6d30be7 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 }
0f3a6a57 658
659 }//end of track loop
660 PostData(1, fOutputList);
1d08b6e8 661}
662//_________________________________________
914042c2 663void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
1d08b6e8 664{
3a94d6fe 665 //--- Check MC
666 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
667 fIsMC = kTRUE;
668 printf("+++++ MC Data available");
669 }
1d08b6e8 670 //--------Initialize PID
3a94d6fe 671 fPID->SetHasMCData(fIsMC);
672
1d08b6e8 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
1d08b6e8 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
0f3a6a57 738 fTPCsubEPres = new TH1F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1);
1d08b6e8 739 fOutputList->Add(fTPCsubEPres);
740
5528bcf6 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);
1d08b6e8 745 fOutputList->Add(fEPres);
746
d6d30be7 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);
1d08b6e8 752 fOutputList->Add(fCorr);
1d08b6e8 753
d6d30be7 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]);
1c424962 794
d6d30be7 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]);
def484dc 797
d6d30be7 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]);
def484dc 800
d6d30be7 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 }
def484dc 806
908103bf 807 fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
937b2ff2 808 fOutputList->Add(fD0_e);
def484dc 809
937b2ff2 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 }
535b11a0 828
1c424962 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);
6522b3b5 834
1c424962 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);
535b11a0 849
0f3a6a57 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);
535b11a0 867
1d08b6e8 868 PostData(1,fOutputList);
869}
870
871//________________________________________________________________________
914042c2 872void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
1d08b6e8 873{
874 // Info("Terminate");
875 AliAnalysisTaskSE::Terminate();
876}
877
878//________________________________________________________________________
914042c2 879Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1d08b6e8 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//_________________________________________
914042c2 887Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 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}
1d08b6e8 896//_________________________________________
914042c2 897Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 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}
1c424962 905//_________________________________________
0f3a6a57 906Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT, Float_t cent) const
1c424962 907{
d6d30be7 908 //Get Pi0 weight
0f3a6a57 909 double weight = 1.0;
1c424962 910 return weight;
911}
912//_________________________________________
0f3a6a57 913Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT, Float_t cent) const
1c424962 914{
937b2ff2 915 //Get eta weight
1c424962 916 double weight = 1.0;
1c424962 917 return weight;
918}
def484dc 919//_________________________________________
d6d30be7 920Double_t AliAnalysisTaskFlowTPCEMCalEP::GetSigmaEMCal(Double_t EoverP, Double_t pt, Float_t cent) const
def484dc 921{
d6d30be7 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;
937b2ff2 1020
d6d30be7 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 }
0f3a6a57 1063 }
d6d30be7 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
def484dc 1213}
0f3a6a57 1214//_________________________________________
1215void AliAnalysisTaskFlowTPCEMCalEP::InitParameters()
1216{
1217 // Init parameters
1218
1219 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1220 fTrackCuts->SetRequireTPCRefit(kTRUE);
1221 fTrackCuts->SetRequireITSRefit(kTRUE);
d6d30be7 1222 fTrackCuts->SetEtaRange(-0.7,0.7);
0f3a6a57 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}
ece431b3 1237