]> 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
update
[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
32d22164 338 Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999., wclsE = -999., wEovP = -999.;//, m02= -999., m20= -999.;
339
1d08b6e8 340 pt = track->Pt();
44fdf7ac 341 if(pt<2) continue;
1d08b6e8 342 fTrkpt->Fill(pt);
3a94d6fe 343
1d08b6e8 344 Int_t clsId = track->GetEMCALcluster();
345 if (clsId>0){
346 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
347 if(cluster && cluster->IsEMCAL()){
348 clsE = cluster->E();
06d74a7b 349// m20 = cluster->M20();
350// m02 = cluster->M02();
1d08b6e8 351 }
352 }
b1a1875f 353
1d08b6e8 354 p = track->P();
355 phi = track->Phi();
356 dEdx = track->GetTPCsignal();
357 EovP = clsE/p;
b1a1875f 358 wEovP = wclsE/p;
1d08b6e8 359 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
360 fdEdxBef->Fill(p,dEdx);
361 fTPCnsigma->Fill(p,fTPCnSigma);
32d22164 362
363 //Remove electron candidate from the event plane
364 Double_t evPlaneCorrTPC = -999.;
365 if(dEdx>70 && dEdx<90){
366 Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
367 Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
368 TVector2 newQVectorfortrack;
369 newQVectorfortrack.Set(qX,qY);
370 evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
371 }
372
06d74a7b 373 Bool_t fFlagPhotonicElec = kFALSE;
374 Bool_t fFlagPhotonicElecBCG = kFALSE;
5945f829 375
06d74a7b 376 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG);
377
32d22164 378 Double_t corr[10]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneCorrTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG};
1d08b6e8 379 fCorr->Fill(corr);
06d74a7b 380
32d22164 381
382 Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
383 Int_t IsElec = 0;
384 Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
385 Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
386 Bool_t MChijing;
387
388 if(fIsMC && fMC && stack){
06d74a7b 389 Int_t label = track->GetLabel();
390 if(label>0){
391 TParticle *particle = stack->Particle(label);
392 if(particle){
32d22164 393 partPDG = particle->GetPdgCode();
394 partPt = particle->Pt();
395
06d74a7b 396 if (TMath::Abs(partPDG)==11) IsElec = 1;
32d22164 397
398 MChijing = fMC->IsFromBGEvent(label);
06d74a7b 399
32d22164 400 int iHijing = 1;
401 if(!MChijing) iHijing = 0; // 0 if enhanced sample
402
06d74a7b 403 Int_t idMother = particle->GetFirstMother();
404 if (idMother>0){
32d22164 405 TParticle *mother = stack->Particle(idMother);
406 motherPt = mother->Pt();
407 motherPDG = mother->GetPdgCode();
408
409
410
411 if (motherPDG==22) whichFirstMother = 3; //gamma
412 if (motherPDG==111) whichFirstMother = 2; //pi0
413 if (motherPDG==221) whichFirstMother = 1; //eta
414
415 Int_t idSecondMother = particle->GetSecondMother();
416 if (idSecondMother>0){
417 TParticle *secondMother = stack->Particle(idSecondMother);
418 secondMotherPt = secondMother->Pt();
419 secondMotherPDG = secondMother->GetPdgCode();
420
421 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
422 if (secondMotherPDG==221) whichSecondMother = 1; //eta
423
424 Int_t idThirdMother = secondMother->GetFirstMother();
425 if (idThirdMother>0){
426 TParticle *thirdMother = stack->Particle(idThirdMother);
427 thirdMotherPt = thirdMother->Pt();
428 thirdMotherPDG = thirdMother->GetPdgCode();
429
430 if (thirdMotherPDG==221) whichThirdMother = 1; //eta
431 }
432 }
433
434 Double_t mc[15]={EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,IsElec,cent,pt,whichFirstMother,whichSecondMother,whichThirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt};
435
436 if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// gamma, pi0, eta
06d74a7b 437 }
438 }
439 }
440 }
06d74a7b 441
1d08b6e8 442
443 if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
444 Int_t pidpassed = 1;
445
446 //--- track accepted
447 AliHFEpidObject hfetrack;
448 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
449 hfetrack.SetRecTrack(track);
450 hfetrack.SetPbPb();
451 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
452
44fdf7ac 453 Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 454 fChargPartV2->Fill(corrV2);
b3bd01b3 455
456 if(fTPCnSigma >= -0.5){
32d22164 457 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
b3bd01b3 458 feTPCV2->Fill(correctedV2);
459 }
460
461 if(pidpassed==0) continue;
462
32d22164 463 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 464
1d08b6e8 465 feV2->Fill(correctedV2);
466
467 fTrkEovPAft->Fill(pt,EovP);
468 fdEdxAft->Fill(p,dEdx);
469
3a94d6fe 470 if(fFlagPhotonicElec){
471 fphoteV2->Fill(correctedV2);
472 fPhotoElecPt->Fill(pt);
473 }
474
475 if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
476
1d08b6e8 477 }
478 PostData(1, fOutputList);
479}
480//_________________________________________
481void AliAnalysisTaskElecV2::UserCreateOutputObjects()
482{
3a94d6fe 483 //--- Check MC
484 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
485 fIsMC = kTRUE;
486 printf("+++++ MC Data available");
487 }
1d08b6e8 488 //--------Initialize PID
3a94d6fe 489 fPID->SetHasMCData(fIsMC);
490
1d08b6e8 491 if(!fPID->GetNumberOfPIDdetectors())
492 {
493 fPID->AddDetector("TPC", 0);
494 fPID->AddDetector("EMCAL", 1);
495 }
496
497 fPID->SortDetectors();
498 fPIDqa = new AliHFEpidQAmanager();
499 fPIDqa->Initialize(fPID);
500
501 //--------Initialize correction Framework and Cuts
502 fCFM = new AliCFManager;
503 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
504 fCFM->SetNStepParticle(kNcutSteps);
505 for(Int_t istep = 0; istep < kNcutSteps; istep++)
506 fCFM->SetParticleCutsList(istep, NULL);
507
508 if(!fCuts){
509 AliWarning("Cuts not available. Default cuts will be used");
510 fCuts = new AliHFEcuts;
511 fCuts->CreateStandardCuts();
512 }
513 fCuts->Initialize(fCFM);
514
515 //---------Output Tlist
516 fOutputList = new TList();
517 fOutputList->SetOwner();
518 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
519
520 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
521 fOutputList->Add(fNoEvents);
522
523 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
524 fOutputList->Add(fTrkpt);
525
526 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
527 fOutputList->Add(fTrackPtBefTrkCuts);
528
529 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
530 fOutputList->Add(fTrackPtAftTrkCuts);
531
532 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
533 fOutputList->Add(fTPCnsigma);
534
535 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
536 fOutputList->Add(fTrkEovPBef);
537
538 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
539 fOutputList->Add(fTrkEovPAft);
540
541 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
542 fOutputList->Add(fdEdxBef);
543
544 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
545 fOutputList->Add(fdEdxAft);
546
547 fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
548 fOutputList->Add(fInvmassLS);
549
550 fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
551 fOutputList->Add(fInvmassULS);
552
553 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
554 fOutputList->Add(fOpeningAngleLS);
555
556 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
557 fOutputList->Add(fOpeningAngleULS);
558
559 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
560 fOutputList->Add(fPhotoElecPt);
561
562 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
563 fOutputList->Add(fSemiInclElecPt);
564
565 fCent = new TH1F("fCent","Centrality",100,0,100) ;
566 fOutputList->Add(fCent);
567
5945f829 568 fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 569 fOutputList->Add(fevPlaneV0A);
570
5945f829 571 fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 572 fOutputList->Add(fevPlaneV0C);
573
5945f829 574 fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 575 fOutputList->Add(fevPlaneTPC);
576
1d08b6e8 577 fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
578 fOutputList->Add(fTPCsubEPres);
579
580 Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
581 Double_t xminv1[4]={-1,-1,-1,0};
582 Double_t xmaxv1[4]={1,1,1,90};
583 fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
584 fOutputList->Add(fEPres);
585
06d74a7b 586 //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m20,m02
587 Int_t binsv2[10]={100,100,90,100,100,100,100,100,3,3};
588 Double_t xminv2[10]={0,-3.5,0,0,0,0,0,0,-1,-1};
589 Double_t xmaxv2[10]={2*TMath::Pi(),3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi(),2,2};
590 fCorr = new THnSparseD ("fCorr","Correlations",10,binsv2,xminv2,xmaxv2);
1d08b6e8 591 fOutputList->Add(fCorr);
592
593 Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
594 Double_t xminv3[5]={0,0,-1,-1,-1};
595 Double_t xmaxv3[5]={90,50,1,1,1};
596 feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
597 fOutputList->Add(feV2);
598
599 Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
600 Double_t xminv4[5]={0,0,-1,-1,-1};
601 Double_t xmaxv4[5]={90,50,1,1,1};
602 fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
603 fOutputList->Add(fphoteV2);
604
44fdf7ac 605 Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
606 Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
607 Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
608 fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
1d08b6e8 609 fOutputList->Add(fChargPartV2);
610
611 Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
612 Double_t xminv6[5]={0,0,-1,-1,-1};
613 Double_t xmaxv6[5]={90,50,1,1,1};
614 feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
615 fOutputList->Add(feTPCV2);
616
32d22164 617 //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,IsElec,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
618 Int_t binsv7[15]={100,100,100,3,3,3,90,100,5,5,5,3,100,100,100};
619 Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
620 Double_t xmaxv7[15]={3,3.5,50,2,2,2,90,50,4,4,4,2,50,50,50};
621 fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
3a94d6fe 622 fOutputList->Add(fMCphotoElecPt);
623
1d08b6e8 624 PostData(1,fOutputList);
625}
626
627//________________________________________________________________________
628void AliAnalysisTaskElecV2::Terminate(Option_t *)
629{
630 // Info("Terminate");
631 AliAnalysisTaskSE::Terminate();
632}
633
634//________________________________________________________________________
635Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
636{
637 // Check single track cuts for a given cut step
638 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
639 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
640 return kTRUE;
641}
642//_________________________________________
06d74a7b 643void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG)
1d08b6e8 644{
645 //Identify non-heavy flavour electrons using Invariant mass method
646
647 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
648 fTrackCuts->SetRequireTPCRefit(kTRUE);
5945f829 649 fTrackCuts->SetRequireITSRefit(kTRUE);
1d08b6e8 650 fTrackCuts->SetEtaRange(-0.7,0.7);
651 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
652 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
653 fTrackCuts->SetMinNClustersTPC(100);
654
655 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
656
657 Bool_t flagPhotonicElec = kFALSE;
06d74a7b 658 Bool_t flagPhotonicElecBCG = kFALSE;
1d08b6e8 659
b1a1875f 660 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
661
5945f829 662 if(jTracks==iTracks) continue;
b1a1875f 663
1d08b6e8 664 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
665 if (!trackAsso) {
666 printf("ERROR: Could not receive track %d\n", jTracks);
667 continue;
668 }
669
5945f829 670 Double_t dEdxAsso = -999., ptAsso=-999.;
06d74a7b 671 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
672 Double_t openingAngle = -999., mass=999., width = -999;
1d08b6e8 673
674 dEdxAsso = trackAsso->GetTPCsignal();
675 ptAsso = trackAsso->Pt();
676 Int_t chargeAsso = trackAsso->Charge();
677 Int_t charge = track->Charge();
678
06d74a7b 679 if(ptAsso <0.5) continue;
1d08b6e8 680 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
06d74a7b 681 if(dEdxAsso <65 || dEdxAsso>100) continue;
1d08b6e8 682
683 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
684 if(charge>0) fPDGe1 = -11;
685 if(chargeAsso>0) fPDGe2 = -11;
686
687 if(charge == chargeAsso) fFlagLS = kTRUE;
688 if(charge != chargeAsso) fFlagULS = kTRUE;
689
690 AliKFParticle ge1(*track, fPDGe1);
691 AliKFParticle ge2(*trackAsso, fPDGe2);
692 AliKFParticle recg(ge1, ge2);
693
694 if(recg.GetNDF()<1) continue;
695 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
696 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
697
698 AliKFVertex primV(*pVtx);
699 primV += recg;
700 recg.SetProductionVertex(primV);
701
702 recg.SetMassConstraint(0,0.0001);
703
704 openingAngle = ge1.GetAngle(ge2);
705 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
706 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
707
5945f829 708 if(openingAngle > fOpeningAngleCut) continue;
06d74a7b 709
710 recg.GetMass(mass,width);
5945f829 711
1d08b6e8 712 if(fFlagLS) fInvmassLS->Fill(mass);
713 if(fFlagULS) fInvmassULS->Fill(mass);
714
5945f829 715 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
06d74a7b 716 if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
1d08b6e8 717
718 }
719 fFlagPhotonicElec = flagPhotonicElec;
06d74a7b 720 fFlagPhotonicElecBCG = flagPhotonicElecBCG;
1d08b6e8 721
722}
723//_________________________________________
724Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
725{
726 //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
727 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
728 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
729 Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
730
731 return cos2DeltaPhi;
732}
733
734//_________________________________________
735Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
736{
737 //Get phi-psi_EP
738 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
739 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
740
741 return dPhi;
742}
b1a1875f 743//_________________________________________
744Double_t AliAnalysisTaskElecV2::GetclusterE(Int_t iTrack, Double_t clsPhi, Double_t clsEta) const
745{
746 //Return E
747 for (Int_t jTracks = 0; jTracks < fESD->GetNumberOfTracks(); jTracks++){
748
749 if(jTracks==iTrack) continue;
750
751 AliESDtrack* wtrack = fESD->GetTrack(jTracks);
752 if (!wtrack) continue;
753
754 Double_t wclsPhi=-999., wclsEta=-999., dPhi=-999., dEta=-999., dR=-999., wclsE=-999.;
755
756 Int_t wclsId = wtrack->GetEMCALcluster();
757 if (wclsId>0){
758 AliESDCaloCluster *wcluster = fESD->GetCaloCluster(wclsId);
759 if(wcluster && wcluster->IsEMCAL()){
760 Float_t wclusterPosition[3]={0,0,0};
761 wcluster->GetPosition(wclusterPosition);
762 TVector3 clsPosVec(wclusterPosition[0],wclusterPosition[1],wclusterPosition[2]);
763 wclsPhi = clsPosVec.Phi();
764 wclsEta = clsPosVec.Eta();
765
766 dPhi = TMath::Abs(wclsPhi - clsPhi);
767 dEta = TMath::Abs(wclsEta - clsEta);
768 dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
769 if(dR>0.15){
770 wclsE = wcluster->E();
771 return wclsE;
772 }
773 }
774 }
775 }
776 return -999.;
777}