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