]> 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
Update of the flow tasks
[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"
37
914042c2 38#include "AliAnalysisTaskFlowTPCEMCalEP.h"
1d08b6e8 39#include "TGeoGlobalMagField.h"
40#include "AliLog.h"
41#include "AliAnalysisTaskSE.h"
42#include "TRefArray.h"
43#include "TVector.h"
44#include "AliESDInputHandler.h"
45#include "AliESDpid.h"
46#include "AliESDtrackCuts.h"
47#include "AliPhysicsSelection.h"
48#include "AliESDCaloCluster.h"
49#include "AliAODCaloCluster.h"
50#include "AliEMCALRecoUtils.h"
51#include "AliEMCALGeometry.h"
52#include "AliGeomManager.h"
53#include "stdio.h"
54#include "TGeoManager.h"
55#include "iostream"
56#include "fstream"
57
58#include "AliEMCALTrack.h"
59#include "AliMagF.h"
60
61#include "AliKFParticle.h"
62#include "AliKFVertex.h"
63
3a94d6fe 64#include "AliMCEventHandler.h"
65#include "AliMCEvent.h"
66#include "AliMCParticle.h"
67#include "AliStack.h"
68
1d08b6e8 69#include "AliPID.h"
70#include "AliPIDResponse.h"
71#include "AliHFEcontainer.h"
72#include "AliHFEcuts.h"
73#include "AliHFEpid.h"
74#include "AliHFEpidBase.h"
75#include "AliHFEpidQAmanager.h"
76#include "AliHFEtools.h"
77#include "AliCFContainer.h"
78#include "AliCFManager.h"
79
80#include "AliEventplane.h"
81#include "AliCentrality.h"
82
914042c2 83ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
1d08b6e8 84//________________________________________________________________________
914042c2 85AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
1d08b6e8 86 : AliAnalysisTaskSE(name)
87 ,fESD(0)
3a94d6fe 88 ,fMC(0)
1d08b6e8 89 ,fOutputList(0)
90 ,fTrackCuts(0)
91 ,fCuts(0)
92 ,fIdentifiedAsOutInz(kFALSE)
93 ,fPassTheEventCut(kFALSE)
94 ,fRejectKinkMother(kFALSE)
3a94d6fe 95 ,fIsMC(kFALSE)
1d08b6e8 96 ,fVz(0.0)
97 ,fCFM(0)
98 ,fPID(0)
99 ,fPIDqa(0)
100 ,fOpeningAngleCut(0.1)
101 ,fInvmassCut(0.01)
102 ,fNoEvents(0)
103 ,fTrkpt(0)
104 ,fTrkEovPBef(0)
105 ,fTrkEovPAft(0)
106 ,fdEdxBef(0)
107 ,fdEdxAft(0)
108 ,fInvmassLS(0)
109 ,fInvmassULS(0)
110 ,fOpeningAngleLS(0)
111 ,fOpeningAngleULS(0)
112 ,fPhotoElecPt(0)
113 ,fSemiInclElecPt(0)
3a94d6fe 114 ,fMCphotoElecPt(0)
1d08b6e8 115 ,fTrackPtBefTrkCuts(0)
116 ,fTrackPtAftTrkCuts(0)
117 ,fTPCnsigma(0)
118 ,fCent(0)
44fdf7ac 119 ,fevPlaneV0A(0)
120 ,fevPlaneV0C(0)
121 ,fevPlaneTPC(0)
1d08b6e8 122 ,fTPCsubEPres(0)
123 ,fEPres(0)
124 ,fCorr(0)
125 ,feTPCV2(0)
126 ,feV2(0)
127 ,fphoteV2(0)
128 ,fChargPartV2(0)
6522b3b5 129 ,fGammaWeight(0)
130 ,fPi0Weight(0)
131 ,fEtaWeight(0)
1d08b6e8 132{
133 //Named constructor
134
135 fPID = new AliHFEpid("hfePid");
136 fTrackCuts = new AliESDtrackCuts();
137
138 // Define input and output slots here
139 // Input slot #0 works with a TChain
140 DefineInput(0, TChain::Class());
141 // Output slot #0 id reserved by the base class for AOD
142 // Output slot #1 writes into a TH1 container
143 // DefineOutput(1, TH1I::Class());
144 DefineOutput(1, TList::Class());
145 // DefineOutput(3, TTree::Class());
146}
147
148//________________________________________________________________________
914042c2 149AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 150 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
151 ,fESD(0)
3a94d6fe 152 ,fMC(0)
1d08b6e8 153 ,fOutputList(0)
154 ,fTrackCuts(0)
155 ,fCuts(0)
156 ,fIdentifiedAsOutInz(kFALSE)
157 ,fPassTheEventCut(kFALSE)
158 ,fRejectKinkMother(kFALSE)
3a94d6fe 159 ,fIsMC(kFALSE)
1d08b6e8 160 ,fVz(0.0)
161 ,fCFM(0)
162 ,fPID(0)
163 ,fPIDqa(0)
164 ,fOpeningAngleCut(0.1)
165 ,fInvmassCut(0.01)
166 ,fNoEvents(0)
167 ,fTrkpt(0)
168 ,fTrkEovPBef(0)
169 ,fTrkEovPAft(0)
170 ,fdEdxBef(0)
171 ,fdEdxAft(0)
172 ,fInvmassLS(0)
173 ,fInvmassULS(0)
174 ,fOpeningAngleLS(0)
175 ,fOpeningAngleULS(0)
176 ,fPhotoElecPt(0)
177 ,fSemiInclElecPt(0)
3a94d6fe 178 ,fMCphotoElecPt(0)
1d08b6e8 179 ,fTrackPtBefTrkCuts(0)
180 ,fTrackPtAftTrkCuts(0)
181 ,fTPCnsigma(0)
182 ,fCent(0)
44fdf7ac 183 ,fevPlaneV0A(0)
184 ,fevPlaneV0C(0)
185 ,fevPlaneTPC(0)
1d08b6e8 186 ,fTPCsubEPres(0)
187 ,fEPres(0)
188 ,fCorr(0)
189 ,feTPCV2(0)
190 ,feV2(0)
191 ,fphoteV2(0)
192 ,fChargPartV2(0)
6522b3b5 193 ,fGammaWeight(0)
194 ,fPi0Weight(0)
195 ,fEtaWeight(0)
1d08b6e8 196{
197 //Default constructor
198 fPID = new AliHFEpid("hfePid");
199
200 fTrackCuts = new AliESDtrackCuts();
201
202 // Constructor
203 // Define input and output slots here
204 // Input slot #0 works with a TChain
205 DefineInput(0, TChain::Class());
206 // Output slot #0 id reserved by the base class for AOD
207 // Output slot #1 writes into a TH1 container
208 // DefineOutput(1, TH1I::Class());
209 DefineOutput(1, TList::Class());
210 //DefineOutput(3, TTree::Class());
211}
212//_________________________________________
213
914042c2 214AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 215{
216 //Destructor
217
218 delete fOutputList;
219 delete fPID;
220 delete fCFM;
221 delete fPIDqa;
222 delete fTrackCuts;
223}
224//_________________________________________
225
914042c2 226void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
1d08b6e8 227{
228 //Main loop
229 //Called for each event
230
231 // create pointer to event
232 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
233 if (!fESD) {
234 printf("ERROR: fESD not available\n");
235 return;
236 }
237
238 if(!fCuts){
239 AliError("HFE cuts not available");
240 return;
241 }
242
243 if(!fPID->IsInitialized()){
244 // Initialize PID with the given run number
245 AliWarning("PID not initialised, get from Run no");
246 fPID->InitializePID(fESD->GetRunNumber());
247 }
248
3a94d6fe 249 if(fIsMC)fMC = MCEvent();
250 AliStack* stack = NULL;
251 if(fIsMC && fMC) stack = fMC->Stack();
1d08b6e8 252
253 Int_t fNOtrks = fESD->GetNumberOfTracks();
254 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
255
256 Double_t pVtxZ = -999;
257 pVtxZ = pVtx->GetZ();
258
259 if(TMath::Abs(pVtxZ)>10) return;
260 fNoEvents->Fill(0);
261
262 if(fNOtrks<2) return;
263
264 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
265 if(!pidResponse){
266 AliDebug(1, "Using default PID Response");
267 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
268 }
269
270 fPID->SetPIDResponse(pidResponse);
271
272 fCFM->SetRecEventInfo(fESD);
273
274 Float_t cent = -1.;
275 AliCentrality *centrality = fESD->GetCentrality();
276 cent = centrality->GetCentralityPercentile("V0M");
277 fCent->Fill(cent);
6522b3b5 278
1d08b6e8 279 //Event planes
280
281 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
282 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
5945f829 283 fevPlaneV0A->Fill(evPlaneV0A,cent);
1d08b6e8 284
285 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
286 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
5945f829 287 fevPlaneV0C->Fill(evPlaneV0C,cent);
1d08b6e8 288
289 AliEventplane* esdTPCep = fESD->GetEventplane();
44fdf7ac 290 TVector2 *standardQ = 0x0;
1d08b6e8 291 Double_t qx = -999., qy = -999.;
44fdf7ac 292 standardQ = esdTPCep->GetQVector();
ee89e744 293 if(!standardQ)return;
294
295 qx = standardQ->X();
296 qy = standardQ->Y();
297
1d08b6e8 298 TVector2 qVectorfortrack;
299 qVectorfortrack.Set(qx,qy);
300 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
5945f829 301 fevPlaneTPC->Fill(evPlaneTPC,cent);
44fdf7ac 302
1d08b6e8 303 TVector2 *qsub1a = esdTPCep->GetQsub1();
304 TVector2 *qsub2a = esdTPCep->GetQsub2();
305 Double_t evPlaneResTPC = -999.;
306 if(qsub1a && qsub2a)
307 {
308 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
309 }
310
311 fTPCsubEPres->Fill(evPlaneResTPC,cent);
312
313 Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
314 fEPres->Fill(evPlaneRes);
315
6522b3b5 316 // Pi0, eta and gamma weights
317
318 if(fIsMC && fMC && stack && cent>20 && cent<40){
319 Int_t nParticles = stack->GetNtrack();
320 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
321 TParticle* particle = stack->Particle(iParticle);
322 int fPDG = particle->GetPdgCode();
323 double pTMC = particle->Pt();
324 double etaMC = particle->Eta();
325 if(fabs(etaMC)>0.7)continue;
326
327 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
328 int iHijing = 1;
329 if(!MChijing)iHijing = 0;
330
331 if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
332 if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
333
334 Int_t idMother = particle->GetFirstMother();
335 if (idMother>0){
336 TParticle *mother = stack->Particle(idMother);
337 int motherPDG = mother->GetPdgCode();
338 if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
339 }
340
341 }
342 }
343
344
1d08b6e8 345 // Track loop
346 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
347 AliESDtrack* track = fESD->GetTrack(iTracks);
348 if (!track) {
349 printf("ERROR: Could not receive track %d\n", iTracks);
350 continue;
351 }
352
353 if(TMath::Abs(track->Eta())>0.7) continue;
354
355 fTrackPtBefTrkCuts->Fill(track->Pt());
b1a1875f 356
1d08b6e8 357 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
358
1d08b6e8 359 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
360 if(track->GetKinkIndex(0) != 0) continue;
361 }
362
1d08b6e8 363 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
364
1d08b6e8 365 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
366
1d08b6e8 367 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
368
369 fTrackPtAftTrkCuts->Fill(track->Pt());
370
32d22164 371 Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999., wclsE = -999., wEovP = -999.;//, m02= -999., m20= -999.;
372
1d08b6e8 373 pt = track->Pt();
44fdf7ac 374 if(pt<2) continue;
1d08b6e8 375 fTrkpt->Fill(pt);
3a94d6fe 376
1d08b6e8 377 Int_t clsId = track->GetEMCALcluster();
378 if (clsId>0){
379 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
380 if(cluster && cluster->IsEMCAL()){
381 clsE = cluster->E();
06d74a7b 382// m20 = cluster->M20();
383// m02 = cluster->M02();
1d08b6e8 384 }
385 }
b1a1875f 386
1d08b6e8 387 p = track->P();
388 phi = track->Phi();
389 dEdx = track->GetTPCsignal();
390 EovP = clsE/p;
b1a1875f 391 wEovP = wclsE/p;
1d08b6e8 392 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
393 fdEdxBef->Fill(p,dEdx);
394 fTPCnsigma->Fill(p,fTPCnSigma);
32d22164 395
396 //Remove electron candidate from the event plane
6522b3b5 397 Float_t evPlaneCorrTPC = evPlaneTPC;
32d22164 398 if(dEdx>70 && dEdx<90){
399 Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
400 Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
401 TVector2 newQVectorfortrack;
402 newQVectorfortrack.Set(qX,qY);
403 evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
404 }
405
06d74a7b 406 Bool_t fFlagPhotonicElec = kFALSE;
407 Bool_t fFlagPhotonicElecBCG = kFALSE;
5945f829 408
06d74a7b 409 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG);
410
32d22164 411 Double_t corr[10]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneCorrTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG};
1d08b6e8 412 fCorr->Fill(corr);
06d74a7b 413
32d22164 414
415 Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
6522b3b5 416 Int_t whichPart = -99;
32d22164 417 Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
418 Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
419 Bool_t MChijing;
420
421 if(fIsMC && fMC && stack){
06d74a7b 422 Int_t label = track->GetLabel();
423 if(label>0){
424 TParticle *particle = stack->Particle(label);
425 if(particle){
32d22164 426 partPDG = particle->GetPdgCode();
427 partPt = particle->Pt();
428
6522b3b5 429 if (TMath::Abs(partPDG)==11) whichPart = 0; //electron
430 if (partPDG==22) whichPart = 3; //gamma
431 if (partPDG==111) whichPart = 2; //pi0
432 if (partPDG==221) whichPart = 1; //eta
433
32d22164 434 MChijing = fMC->IsFromBGEvent(label);
06d74a7b 435
32d22164 436 int iHijing = 1;
437 if(!MChijing) iHijing = 0; // 0 if enhanced sample
438
06d74a7b 439 Int_t idMother = particle->GetFirstMother();
440 if (idMother>0){
32d22164 441 TParticle *mother = stack->Particle(idMother);
442 motherPt = mother->Pt();
443 motherPDG = mother->GetPdgCode();
444
32d22164 445 if (motherPDG==22) whichFirstMother = 3; //gamma
446 if (motherPDG==111) whichFirstMother = 2; //pi0
447 if (motherPDG==221) whichFirstMother = 1; //eta
448
449 Int_t idSecondMother = particle->GetSecondMother();
450 if (idSecondMother>0){
451 TParticle *secondMother = stack->Particle(idSecondMother);
452 secondMotherPt = secondMother->Pt();
453 secondMotherPDG = secondMother->GetPdgCode();
454
455 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
456 if (secondMotherPDG==221) whichSecondMother = 1; //eta
457
458 Int_t idThirdMother = secondMother->GetFirstMother();
459 if (idThirdMother>0){
460 TParticle *thirdMother = stack->Particle(idThirdMother);
461 thirdMotherPt = thirdMother->Pt();
462 thirdMotherPDG = thirdMother->GetPdgCode();
463
464 if (thirdMotherPDG==221) whichThirdMother = 1; //eta
465 }
466 }
467
6522b3b5 468 Double_t mc[15]={EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,whichFirstMother,whichSecondMother,whichThirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt};
469 fMCphotoElecPt->Fill(mc);
470// if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
06d74a7b 471 }
472 }
473 }
474 }
06d74a7b 475
1d08b6e8 476
477 if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
478 Int_t pidpassed = 1;
479
480 //--- track accepted
481 AliHFEpidObject hfetrack;
482 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
483 hfetrack.SetRecTrack(track);
484 hfetrack.SetPbPb();
485 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
486
44fdf7ac 487 Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 488 fChargPartV2->Fill(corrV2);
b3bd01b3 489
490 if(fTPCnSigma >= -0.5){
32d22164 491 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
b3bd01b3 492 feTPCV2->Fill(correctedV2);
493 }
494
495 if(pidpassed==0) continue;
496
32d22164 497 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 498
1d08b6e8 499 feV2->Fill(correctedV2);
500
501 fTrkEovPAft->Fill(pt,EovP);
502 fdEdxAft->Fill(p,dEdx);
503
3a94d6fe 504 if(fFlagPhotonicElec){
505 fphoteV2->Fill(correctedV2);
506 fPhotoElecPt->Fill(pt);
507 }
508
509 if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
510
1d08b6e8 511 }
512 PostData(1, fOutputList);
513}
514//_________________________________________
914042c2 515void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
1d08b6e8 516{
3a94d6fe 517 //--- Check MC
518 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
519 fIsMC = kTRUE;
520 printf("+++++ MC Data available");
521 }
1d08b6e8 522 //--------Initialize PID
3a94d6fe 523 fPID->SetHasMCData(fIsMC);
524
1d08b6e8 525 if(!fPID->GetNumberOfPIDdetectors())
526 {
527 fPID->AddDetector("TPC", 0);
528 fPID->AddDetector("EMCAL", 1);
529 }
530
531 fPID->SortDetectors();
532 fPIDqa = new AliHFEpidQAmanager();
533 fPIDqa->Initialize(fPID);
534
535 //--------Initialize correction Framework and Cuts
536 fCFM = new AliCFManager;
537 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
538 fCFM->SetNStepParticle(kNcutSteps);
539 for(Int_t istep = 0; istep < kNcutSteps; istep++)
540 fCFM->SetParticleCutsList(istep, NULL);
541
542 if(!fCuts){
543 AliWarning("Cuts not available. Default cuts will be used");
544 fCuts = new AliHFEcuts;
545 fCuts->CreateStandardCuts();
546 }
547 fCuts->Initialize(fCFM);
548
549 //---------Output Tlist
550 fOutputList = new TList();
551 fOutputList->SetOwner();
552 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
553
554 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
555 fOutputList->Add(fNoEvents);
556
557 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
558 fOutputList->Add(fTrkpt);
559
560 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
561 fOutputList->Add(fTrackPtBefTrkCuts);
562
563 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
564 fOutputList->Add(fTrackPtAftTrkCuts);
565
566 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
567 fOutputList->Add(fTPCnsigma);
568
569 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
570 fOutputList->Add(fTrkEovPBef);
571
572 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
573 fOutputList->Add(fTrkEovPAft);
574
575 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
576 fOutputList->Add(fdEdxBef);
577
578 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
579 fOutputList->Add(fdEdxAft);
580
581 fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
582 fOutputList->Add(fInvmassLS);
583
584 fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
585 fOutputList->Add(fInvmassULS);
586
587 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
588 fOutputList->Add(fOpeningAngleLS);
589
590 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
591 fOutputList->Add(fOpeningAngleULS);
592
593 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
594 fOutputList->Add(fPhotoElecPt);
595
596 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
597 fOutputList->Add(fSemiInclElecPt);
598
599 fCent = new TH1F("fCent","Centrality",100,0,100) ;
600 fOutputList->Add(fCent);
601
5945f829 602 fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 603 fOutputList->Add(fevPlaneV0A);
604
5945f829 605 fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 606 fOutputList->Add(fevPlaneV0C);
607
5945f829 608 fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 609 fOutputList->Add(fevPlaneTPC);
610
1d08b6e8 611 fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
612 fOutputList->Add(fTPCsubEPres);
613
614 Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
615 Double_t xminv1[4]={-1,-1,-1,0};
616 Double_t xmaxv1[4]={1,1,1,90};
617 fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
618 fOutputList->Add(fEPres);
619
06d74a7b 620 //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m20,m02
621 Int_t binsv2[10]={100,100,90,100,100,100,100,100,3,3};
622 Double_t xminv2[10]={0,-3.5,0,0,0,0,0,0,-1,-1};
623 Double_t xmaxv2[10]={2*TMath::Pi(),3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi(),2,2};
624 fCorr = new THnSparseD ("fCorr","Correlations",10,binsv2,xminv2,xmaxv2);
1d08b6e8 625 fOutputList->Add(fCorr);
626
627 Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
628 Double_t xminv3[5]={0,0,-1,-1,-1};
629 Double_t xmaxv3[5]={90,50,1,1,1};
630 feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
631 fOutputList->Add(feV2);
632
633 Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
634 Double_t xminv4[5]={0,0,-1,-1,-1};
635 Double_t xmaxv4[5]={90,50,1,1,1};
636 fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
637 fOutputList->Add(fphoteV2);
638
44fdf7ac 639 Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
640 Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
641 Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
642 fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
1d08b6e8 643 fOutputList->Add(fChargPartV2);
644
645 Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
646 Double_t xminv6[5]={0,0,-1,-1,-1};
647 Double_t xmaxv6[5]={90,50,1,1,1};
648 feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
649 fOutputList->Add(feTPCV2);
650
6522b3b5 651 //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
652 Int_t binsv7[15]={100,100,100,3,3,5,90,100,5,5,5,3,100,100,100};
32d22164 653 Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
6522b3b5 654 Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
32d22164 655 fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
3a94d6fe 656 fOutputList->Add(fMCphotoElecPt);
6522b3b5 657
658 fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
659 fOutputList->Add(fGammaWeight);
660
661 fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
662 fOutputList->Add(fPi0Weight);
663
664 fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
665 fOutputList->Add(fEtaWeight);
666
1d08b6e8 667 PostData(1,fOutputList);
668}
669
670//________________________________________________________________________
914042c2 671void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
1d08b6e8 672{
673 // Info("Terminate");
674 AliAnalysisTaskSE::Terminate();
675}
676
677//________________________________________________________________________
914042c2 678Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1d08b6e8 679{
680 // Check single track cuts for a given cut step
681 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
682 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
683 return kTRUE;
684}
685//_________________________________________
914042c2 686void AliAnalysisTaskFlowTPCEMCalEP::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG)
1d08b6e8 687{
688 //Identify non-heavy flavour electrons using Invariant mass method
689
690 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
691 fTrackCuts->SetRequireTPCRefit(kTRUE);
5945f829 692 fTrackCuts->SetRequireITSRefit(kTRUE);
1d08b6e8 693 fTrackCuts->SetEtaRange(-0.7,0.7);
694 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
695 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
696 fTrackCuts->SetMinNClustersTPC(100);
697
698 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
699
700 Bool_t flagPhotonicElec = kFALSE;
06d74a7b 701 Bool_t flagPhotonicElecBCG = kFALSE;
1d08b6e8 702
b1a1875f 703 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
704
5945f829 705 if(jTracks==iTracks) continue;
b1a1875f 706
1d08b6e8 707 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
708 if (!trackAsso) {
709 printf("ERROR: Could not receive track %d\n", jTracks);
710 continue;
711 }
712
5945f829 713 Double_t dEdxAsso = -999., ptAsso=-999.;
06d74a7b 714 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
715 Double_t openingAngle = -999., mass=999., width = -999;
1d08b6e8 716
717 dEdxAsso = trackAsso->GetTPCsignal();
718 ptAsso = trackAsso->Pt();
719 Int_t chargeAsso = trackAsso->Charge();
720 Int_t charge = track->Charge();
721
06d74a7b 722 if(ptAsso <0.5) continue;
1d08b6e8 723 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
06d74a7b 724 if(dEdxAsso <65 || dEdxAsso>100) continue;
1d08b6e8 725
726 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
727 if(charge>0) fPDGe1 = -11;
728 if(chargeAsso>0) fPDGe2 = -11;
729
730 if(charge == chargeAsso) fFlagLS = kTRUE;
731 if(charge != chargeAsso) fFlagULS = kTRUE;
732
733 AliKFParticle ge1(*track, fPDGe1);
734 AliKFParticle ge2(*trackAsso, fPDGe2);
735 AliKFParticle recg(ge1, ge2);
736
737 if(recg.GetNDF()<1) continue;
738 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
739 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
740
741 AliKFVertex primV(*pVtx);
742 primV += recg;
743 recg.SetProductionVertex(primV);
744
745 recg.SetMassConstraint(0,0.0001);
746
747 openingAngle = ge1.GetAngle(ge2);
748 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
749 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
750
5945f829 751 if(openingAngle > fOpeningAngleCut) continue;
06d74a7b 752
753 recg.GetMass(mass,width);
5945f829 754
1d08b6e8 755 if(fFlagLS) fInvmassLS->Fill(mass);
756 if(fFlagULS) fInvmassULS->Fill(mass);
757
5945f829 758 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
06d74a7b 759 if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
1d08b6e8 760
761 }
762 fFlagPhotonicElec = flagPhotonicElec;
06d74a7b 763 fFlagPhotonicElecBCG = flagPhotonicElecBCG;
1d08b6e8 764
765}
766//_________________________________________
914042c2 767Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 768{
769 //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
770 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
771 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
772 Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
773
774 return cos2DeltaPhi;
775}
776
777//_________________________________________
914042c2 778Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 779{
780 //Get phi-psi_EP
781 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
782 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
783
784 return dPhi;
785}