]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSOH.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
move EMCALJetTasks from PWGGA to PWGJE
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskSOH.cxx
CommitLineData
4bbf59e5 1// $Id$
2//
3// Simulation EMCal task.
4//
5// Author: Saehanseul Oh
6
7#include "AliAnalysisTaskSOH.h"
8
9#include <TH1F.h>
10#include <TH2F.h>
9733b37f 11#include <TH3F.h>
4bbf59e5 12
13#include "AliAnalysisManager.h"
14#include "AliAnalysisTask.h"
15#include "AliEMCALRecoUtils.h"
16#include "AliEMCALTrack.h"
17#include "AliESDCaloCluster.h"
18#include "AliESDEvent.h"
19#include "AliESDInputHandler.h"
20#include "AliESDtrack.h"
21#include "AliESDtrackCuts.h"
22#include "AliExternalTrackParam.h"
23#include "AliLog.h"
24#include "AliMCEvent.h"
25
26ClassImp(AliAnalysisTaskSOH)
27
28//________________________________________________________________________
29AliAnalysisTaskSOH::AliAnalysisTaskSOH() :
30 AliAnalysisTaskSE(),
31 fESD(0),
32 fMC(0),
33 fEsdTrackCuts(0x0),
34 fHybridTrackCuts1(0x0),
35 fHybridTrackCuts2(0x0),
36 fTrackIndices(0x0),
37 fClusterIndices(0x0),
38 fClusterArray(0x0),
39 fOutputList(0x0),
40 fHEventStat(0),
9733b37f 41 fHTrkEffParGenPtEtaPhi(0),
42 fHTrkEffDetGenPtEtaPhi(0),
43 fHTrkEffDetRecPtEtaPhi(0),
44 fHTrkEffDetRecFakePtEtaPhi(0),
4bbf59e5 45 fHScaleFactor(0),
46 fHScaleFactor100HC(0),
47 fHEOverPVsPt(0x0),
48 fHEMCalResponsePion(0x0),
49 fHEMCalResponseElec(0x0),
50 fHEMCalResponseProton(0x0),
51 fHEMCalRecdPhidEta(0x0),
52 fHEMCalRecdPhidEtaP(0x0),
53 fHEMCalRecdPhidEtaM(0x0),
54 fHEMCalRecdPhidEta_Truth(0x0),
55 fHEMCalRecdPhidEtaP_Truth(0x0),
56 fHEMCalRecdPhidEtaM_Truth(0x0),
57 fHEMCalRecdPhidEtaposEta(0x0),
58 fHEMCalRecdPhidEtanegEta(0x0),
faba08b8 59 fHPhotonEdiff100HC(0x0),
a64777f6 60 fHPhotonEdiff70HC(0),
61 fHPhotonEdiff30HC(0),
faba08b8 62 fHPhotonEdiff0HC(0x0),
a64777f6 63 fHPhotonEVsClsE(0x0),
64 fHistEsub1Pch(0x0),
65 fHistEsub2Pch(0x0),
66 fHistEsub1PchRat(0x0),
67 fHistEsub2PchRat(0x0)
4bbf59e5 68{
69 // Constructor
70
4bbf59e5 71}
72
73
74//________________________________________________________________________
75AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :
76 AliAnalysisTaskSE(name),
77 fESD(0),
78 fMC(0),
79 fEsdTrackCuts(0x0),
80 fHybridTrackCuts1(0x0),
81 fHybridTrackCuts2(0x0),
82 fTrackIndices(0x0),
83 fClusterIndices(0x0),
84 fClusterArray(0x0),
85 fOutputList(0x0),
86 fHEventStat(0),
9733b37f 87 fHTrkEffParGenPtEtaPhi(0),
88 fHTrkEffDetGenPtEtaPhi(0),
89 fHTrkEffDetRecPtEtaPhi(0),
90 fHTrkEffDetRecFakePtEtaPhi(0),
4bbf59e5 91 fHScaleFactor(0),
92 fHScaleFactor100HC(0),
93 fHEOverPVsPt(0x0),
af0280a9 94 fHEMCalResponsePion(0x0),
95 fHEMCalResponseElec(0x0),
96 fHEMCalResponseProton(0x0),
4bbf59e5 97 fHEMCalRecdPhidEta(0x0),
98 fHEMCalRecdPhidEtaP(0x0),
99 fHEMCalRecdPhidEtaM(0x0),
100 fHEMCalRecdPhidEta_Truth(0x0),
101 fHEMCalRecdPhidEtaP_Truth(0x0),
102 fHEMCalRecdPhidEtaM_Truth(0x0),
103 fHEMCalRecdPhidEtaposEta(0x0),
104 fHEMCalRecdPhidEtanegEta(0x0),
faba08b8 105 fHPhotonEdiff100HC(0x0),
a64777f6 106 fHPhotonEdiff70HC(0),
107 fHPhotonEdiff30HC(0),
faba08b8 108 fHPhotonEdiff0HC(0x0),
a64777f6 109 fHPhotonEVsClsE(0x0),
110 fHistEsub1Pch(0x0),
111 fHistEsub2Pch(0x0),
112 fHistEsub1PchRat(0x0),
113 fHistEsub2PchRat(0x0)
4bbf59e5 114{
115 // Constructor
116
117 // Output slot #1 writes into a TH1 container
118 DefineOutput(1, TList::Class());
119}
120
121
122//________________________________________________________________________
123AliAnalysisTaskSOH::~AliAnalysisTaskSOH()
124{
125 // Destructor.
126
127 delete fEsdTrackCuts;
128 delete fHybridTrackCuts1;
129 delete fHybridTrackCuts2;
130 delete fTrackIndices;
131 delete fClusterIndices;
132 delete fClusterArray;
133}
134
135
136//________________________________________________________________________
137void AliAnalysisTaskSOH::UserCreateOutputObjects()
138{
139 // Create histograms, called once.
140
141 OpenFile(1);
142
143 fOutputList = new TList();
144 fOutputList->SetOwner(1);
145
146 fHEventStat = new TH1F("fHEventStat","Event statistics for analysis",8,0,8);
147 fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
148 fHEventStat->GetXaxis()->SetBinLabel(2,"cluster");
149 fHEventStat->GetXaxis()->SetBinLabel(3,"good cluster");
150 fHEventStat->GetXaxis()->SetBinLabel(4,"cls/0-truth");
151 fHEventStat->GetXaxis()->SetBinLabel(5,"cls/1-truth");
152 fHEventStat->GetXaxis()->SetBinLabel(6,"cls/2-truth");
153 fHEventStat->GetXaxis()->SetBinLabel(7,"cls/2-goodtruth");
154 fHEventStat->GetXaxis()->SetBinLabel(8,"cls/>3-truth");
155 fOutputList->Add(fHEventStat);
156
9733b37f 157 fHTrkEffParGenPtEtaPhi = new TH3F("fHTrkEffParGenPtEtaPhi","Particle level truth Phi-Eta-p_{T} distribution of primary charged pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
158 fHTrkEffParGenPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
159 fHTrkEffParGenPtEtaPhi->GetYaxis()->SetTitle("#eta");
160 fHTrkEffParGenPtEtaPhi->GetZaxis()->SetTitle("#phi");
161 fOutputList->Add(fHTrkEffParGenPtEtaPhi);
162
163 fHTrkEffDetGenPtEtaPhi = new TH3F("fHTrkEffDetGenPtEtaPhi","Detector level truth Phi-Eta-p_{T} distribution of primary charged pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
164 fHTrkEffDetGenPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
165 fHTrkEffDetGenPtEtaPhi->GetYaxis()->SetTitle("#eta");
166 fHTrkEffDetGenPtEtaPhi->GetZaxis()->SetTitle("#phi");
167 fOutputList->Add(fHTrkEffDetGenPtEtaPhi);
168
169 fHTrkEffDetRecPtEtaPhi = new TH3F("fHTrkEffDetRecPtEtaPhi","Reconstructed track Phi-Eta-p_{T} distribution of primary charged pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
170 fHTrkEffDetRecPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
171 fHTrkEffDetRecPtEtaPhi->GetYaxis()->SetTitle("#eta");
172 fHTrkEffDetRecPtEtaPhi->GetZaxis()->SetTitle("#phi");
173 fOutputList->Add(fHTrkEffDetRecPtEtaPhi);
174
175 fHTrkEffDetRecFakePtEtaPhi = new TH3F("fHTrkEffDetRecFakePtEtaPhi","Reconstructed fake track Phi-Eta-p_{T} distribution of pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
176 fHTrkEffDetRecFakePtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
177 fHTrkEffDetRecFakePtEtaPhi->GetYaxis()->SetTitle("#eta");
178 fHTrkEffDetRecFakePtEtaPhi->GetZaxis()->SetTitle("#phi");
179 fOutputList->Add(fHTrkEffDetRecFakePtEtaPhi);
de64a00e 180
4bbf59e5 181 fHScaleFactor = new TH1F("fHScaleFactor", "Scale factor distribution without hadronic correction;Scale factor",100,0,10);
182 fOutputList->Add(fHScaleFactor);
183
184 fHScaleFactor100HC = new TH1F("fHScaleFactor100HC", "Scale factor distribution with 100% hadronic correction;Scale factor",100,0,10);
185 fOutputList->Add(fHScaleFactor100HC);
186
187 fHEOverPVsPt = new TH2F("fHEOverPVsPt", "E/P vs track p_{T}; p_{T} (GeV/c); E/P", 200 , 0, 4, 200, 0, 3.2);
188 fOutputList->Add(fHEOverPVsPt);
189
190 fHEMCalResponsePion = new TH2F("fHEMCalResponsePion", "Pion E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
191 fOutputList->Add(fHEMCalResponsePion);
192
193 fHEMCalResponseElec = new TH2F("fHEMCalResponseElec", "Electron E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
194 fOutputList->Add(fHEMCalResponseElec);
195
196 fHEMCalResponseProton = new TH2F("fHEMCalResponseProton", "Proton E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
be6872c8 197 fOutputList->Add(fHEMCalResponseProton);
4bbf59e5 198
199 fHEMCalRecdPhidEta = new TH2F("fHEMCalRecdPhidEta","EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
200 fOutputList->Add(fHEMCalRecdPhidEta);
201
202 fHEMCalRecdPhidEtaP = new TH2F("fHEMCalRecdPhidEtaP","EMCAL Charge+ Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
203 fOutputList->Add(fHEMCalRecdPhidEtaP);
204
205 fHEMCalRecdPhidEtaM = new TH2F("fHEMCalRecdPhidEtaM","EMCAL Charge- Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
206 fOutputList->Add(fHEMCalRecdPhidEtaM);
207
208 fHEMCalRecdPhidEta_Truth = new TH2F("fHEMCalRecdPhidEta_Truth","EMCAL Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
209 fOutputList->Add(fHEMCalRecdPhidEta_Truth);
210
211 fHEMCalRecdPhidEtaP_Truth = new TH2F("fHEMCalRecdPhidEtaP_Truth","EMCAL Charge+ Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
212 fOutputList->Add(fHEMCalRecdPhidEtaP_Truth);
213
214 fHEMCalRecdPhidEtaM_Truth = new TH2F("fHEMCalRecdPhidEtaM_Truth","EMCAL Charge- Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
215 fOutputList->Add(fHEMCalRecdPhidEtaM_Truth);
216
217 fHEMCalRecdPhidEtaposEta = new TH2F("fHEMCalRecdPhidEtaposEta","(+eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
218 fOutputList->Add(fHEMCalRecdPhidEtaposEta);
219
220 fHEMCalRecdPhidEtanegEta = new TH2F("fHEMCalRecdPhidEtanegEta","(-eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
221 fOutputList->Add(fHEMCalRecdPhidEtanegEta);
222
a64777f6 223 fHPhotonEdiff100HC = new TH2F("fHPhotonEdiff100HC","Photon (E_{Truth}- E_{calc,100% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{calc,100% HC})/E_{Truth}",1000,0,10,600,-4.9,1.1);
4bbf59e5 224 fOutputList->Add(fHPhotonEdiff100HC);
225
a64777f6 226 fHPhotonEdiff70HC = new TH2F("fHPhotonEdiff70HC","Photon (E_{Truth}- E_{calc,70% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{calc,30% HC})/E_{Truth}",1000,0,10,600,-4.9,1.1);
227 fOutputList->Add(fHPhotonEdiff70HC);
228
229 fHPhotonEdiff30HC = new TH2F("fHPhotonEdiff30HC","Photon (E_{Truth}- E_{calc,30% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{calc,30% HC})/E_{Truth}",1000,0,10,600,-4.9,1.1);
230 fOutputList->Add(fHPhotonEdiff30HC);
231
232 fHPhotonEdiff0HC = new TH2F("fHPhotonEdiff0HC","Photon (E_{Truth}- E_{calc,0% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{cls})/E_{Truth}",1000,0,10,600,-4.9,1.1);
faba08b8 233 fOutputList->Add(fHPhotonEdiff0HC);
234
235 fHPhotonEVsClsE = new TH2F("fHPhotonEVsClsE","Cluster E vs. photon E_{Truth}; photon E_{Truth} (GeV); Cluster E (GeV)",500,0,5,500,0,5);
236 fOutputList->Add(fHPhotonEVsClsE);
237
a64777f6 238 fHistEsub1Pch =new TH2F("fHistEsub1Pch", "(subtracted E in 100% HC) vs. total track P, clusters with 1 matching track; total track P (GeV/c); E_{sub}(GeV)" , 1000, 0., 10, 1000, 0., 10.);
239 fOutputList->Add(fHistEsub1Pch);
240
241 fHistEsub2Pch =new TH2F("fHistEsub2Pch", "(subtracted E in 100% HC) vs. total track P, clusters with 2 matching tracks; total track P (GeV/c); E_{sub}(GeV)" , 1000, 0., 10, 1000, 0., 10.);
242 fOutputList->Add(fHistEsub2Pch);
243
244 fHistEsub1PchRat =new TH2F("fHistEsub1PchRat", "(subtracted E in 100% HC)/total track P vs. total track P, clusters with 1 matching track; total track P (GeV/c); E_{sub}/P_{tot}" , 1000, 0., 10, 1100, 0., 1.1);
245 fOutputList->Add(fHistEsub1PchRat);
246
247 fHistEsub2PchRat =new TH2F("fHistEsub2PchRat", "(subtracted E in 100% HC)/total track P vs. total track P, clusters with 2 matching tracks; total track P (GeV/c); E_{sub}/P_{tot}" , 1000, 0., 10, 1100, 0., 1.1);
248 fOutputList->Add(fHistEsub2PchRat);
249
4bbf59e5 250 fTrackIndices = new TArrayI();
251 fClusterIndices = new TArrayI();
252
253 fClusterArray = new TObjArray();
254 fClusterArray->SetOwner(1);
255
256 PostData(1, fOutputList);
257}
258
259//________________________________________________________________________
260void AliAnalysisTaskSOH::UserExec(Option_t *)
261{
262 // Main loop, called for each event.
263
264 // Post output data.
265 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
266 if (!fESD) {
267 printf("ERROR: fESD not available\n");
268 return;
269 }
270
271 fMC = MCEvent();
272 if (!fMC) {
273 printf("ERROR: fMC not available\n");
274 return;
275 }
276
277 fHEventStat->Fill(0.5);
278
279 if(fTrackIndices)
280 fTrackIndices->Reset();
281 if(fClusterIndices)
282 fClusterIndices->Reset();
283 if(fClusterArray)
284 fClusterArray->Delete();
285
286 ProcessTrack();
287 ProcessCluster();
288 ProcessMc();
289 ProcessScaleFactor();
290
291 PostData(1, fOutputList);
292}
293
294//________________________________________________________________________
295void AliAnalysisTaskSOH::ProcessTrack()
296{
297 // Process track.
298
299 fTrackIndices->Set(fESD->GetNumberOfTracks());
300 AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));
301
302 Int_t isMth = 0;
303 Int_t nTracks = 0;
304
305 Float_t ClsPos[3] = {-999,-999,-999};
306 Double_t emcTrkpos[3] = {-999,-999,-999};
307
308 for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)
309 {
310 AliESDtrack *esdtrack = fESD->GetTrack(itr);
311 if(!esdtrack)continue;
312 AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
313 if(!newTrack) continue;
314
315 Double_t clsE = -1;
316 Int_t clsIndex = newTrack->GetEMCALcluster();
317 if(newTrack->GetEMCALcluster()>-1)
318 {
319 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);
320 if(IsGoodCluster(cluster))
321 {
322 isMth=1;
323
324 cluster->GetPosition(ClsPos);
325 TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
326
327 AliEMCALTrack EMCTrk(*newTrack);
328 if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {
329 continue;
330 }
331 EMCTrk.GetXYZ(emcTrkpos);
332 TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);
333
334 Double_t dPhi = VClsPos.Phi() - VemcTrkPos.Phi();
335 if (dPhi < -1*TMath::Pi()) dPhi += (2*TMath::Pi());
336 else if (dPhi > TMath::Pi()) dPhi -= (2*TMath::Pi());
337
338 Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();
339
340 fHEMCalRecdPhidEta->Fill(dEta, dPhi);
341
342 if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())
343 {
344 AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());
345 if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))
346 {
347 fHEMCalRecdPhidEta_Truth->Fill(dEta, dPhi);
348 if(vParticle->Charge() > 0) fHEMCalRecdPhidEtaP_Truth->Fill(dEta, dPhi);
349 if(vParticle->Charge() < 0) fHEMCalRecdPhidEtaM_Truth->Fill(dEta, dPhi);
350 }
351 }
352
353 if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}
354 if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}
355
356 if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);
357 if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);
358
359 clsE = cluster->E();
360 if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());
361 }
362
363 Int_t ipart = newTrack->GetLabel();
364 if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
365 {
366 AliVParticle* vParticle = fMC->GetTrack(ipart);
367 if(isMth && vParticle)
368 {
369 if(TMath::Abs(vParticle->PdgCode())==211)
370 {
371 fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());
372 }
373 if(TMath::Abs(vParticle->PdgCode())==11)
374 {
375 fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());
376 }
377 if(TMath::Abs(vParticle->PdgCode())==2212)
378 {
379 fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());
380 }
381 }
382 }
383 }
de64a00e 384
385 // fake and secondary tracks
9733b37f 386 if(newTrack->GetLabel()<0 && newTrack->GetPID()==2) fHTrkEffDetRecFakePtEtaPhi->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
de64a00e 387
388 // Track Indices
4bbf59e5 389 fTrackIndices->AddAt(itr,nTracks);
390 nTracks++;
391 }
392
393 fTrackIndices->Set(nTracks);
394}
395
396//________________________________________________________________________
397void AliAnalysisTaskSOH::ProcessCluster()
398{
399 // Process cluster.
400
401 Int_t nCluster = 0;
402 TLorentzVector gamma;
403 Double_t vertex[3] = {0, 0, 0};
404 fESD->GetVertex()->GetXYZ(vertex);
405 const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
406 fClusterIndices->Set(nCaloClusters);
407
408 for(Int_t itr=0; itr<nCaloClusters; itr++)
409 {
410 fHEventStat->Fill(1.5);
411 AliESDCaloCluster *cluster = fESD->GetCaloCluster(itr);
412 if(!IsGoodCluster(cluster)) continue;
413 cluster->GetMomentum(gamma, vertex);
414 if (gamma.Pt() < 0.15) continue;
415 fHEventStat->Fill(2.5);
416
a64777f6 417 TArrayI *TrackLabels = cluster->GetTracksMatched();
418
419 if(TrackLabels->GetSize() == 1)
420 {
421 AliESDtrack *esdtrack = fESD->GetTrack(TrackLabels->GetAt(0));
422 AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
423 if(newTrack && TMath::Abs(newTrack->Eta())<0.7)
424 {
425 Double_t Esub = newTrack->P();
426 if (Esub > cluster->E()) Esub = cluster->E();
427 fHistEsub1Pch->Fill(newTrack->P(), Esub);
428 fHistEsub1PchRat->Fill(newTrack->P(), Esub/newTrack->P());
429 }
430 }
431
432 if(TrackLabels->GetSize() == 2)
433 {
434 AliESDtrack *esdtrack1 = fESD->GetTrack(TrackLabels->GetAt(0));
435 AliESDtrack *esdtrack2 = fESD->GetTrack(TrackLabels->GetAt(1));
436 AliESDtrack *newTrack1 = GetAcceptTrack(esdtrack1);
437 AliESDtrack *newTrack2 = GetAcceptTrack(esdtrack2);
438 if(newTrack1 && newTrack2 && TMath::Abs(newTrack1->Eta())<0.7 && TMath::Abs(newTrack2->Eta())<0.7)
439 {
440 Double_t Esub = newTrack1->P() + newTrack2->P();
441 if (Esub > cluster->E()) Esub = cluster->E();
442 fHistEsub2Pch->Fill(newTrack1->P() + newTrack2->P(), Esub);
443 fHistEsub2PchRat->Fill(newTrack1->P() + newTrack2->P(), Esub/(newTrack1->P() + newTrack2->P()));
444 }
445 else if(newTrack1 && !(newTrack2) && TMath::Abs(newTrack1->Eta())<0.7)
446 {
447 Double_t Esub = newTrack1->P();
448 if (Esub > cluster->E()) Esub = cluster->E();
449 fHistEsub1Pch->Fill(newTrack1->P(), Esub);
450 fHistEsub1PchRat->Fill(newTrack1->P(), Esub/newTrack1->P());
451 }
452 else if (!(newTrack1) && newTrack2 && TMath::Abs(newTrack2->Eta())<0.7)
453 {
454 Double_t Esub = newTrack2->P();
455 if (Esub > cluster->E()) Esub = cluster->E();
456 fHistEsub1Pch->Fill(newTrack2->P(), Esub);
457 fHistEsub1PchRat->Fill(newTrack2->P(), Esub/newTrack2->P());
458 }
459 else {;}
460 }
461
4bbf59e5 462 TArrayI *MCLabels = cluster->GetLabelsArray();
463
464 if(MCLabels->GetSize() == 0) fHEventStat->Fill(3.5);
465 if(MCLabels->GetSize() == 1) fHEventStat->Fill(4.5);
466 if(MCLabels->GetSize() == 2)
467 {
468 fHEventStat->Fill(5.5);
469 AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
470 AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->GetAt(1));
471 if(IsGoodMcParticle(vParticle1, MCLabels->GetAt(0)) && IsGoodMcParticle(vParticle2, MCLabels->GetAt(1)))
472 {
473 fHEventStat->Fill(6.5);
474 if((vParticle1->PdgCode()==22) && (vParticle2->PdgCode()==22)) {;}
475 else if((vParticle1->PdgCode()!=22) && (vParticle2->PdgCode()!=22)) {;}
476 else
477 {
478 fClusterIndices->AddAt(itr,nCluster);
479 nCluster++;
480 }
481 }
482 }
483 if(MCLabels->GetSize() > 2) fHEventStat->Fill(7.5);
484
485 AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cluster);
486
487 Double_t subE = 0;
488 TArrayI arrayTrackMatched(fTrackIndices->GetSize());
489 Int_t nGoodMatch = 0;
490
491 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
492 {
493 AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
494 if(itr==trk->GetEMCALcluster())
495 {
496 arrayTrackMatched[nGoodMatch] = j;
497 nGoodMatch ++;
498 subE += trk->P();
499 }
500 }
501
502 arrayTrackMatched.Set(nGoodMatch);
503 newCluster->AddTracksMatched(arrayTrackMatched);
504
505 Double_t clsE = newCluster->E();
506 Double_t newE = clsE-subE;
507 if(newE<0) newE = 0;
508 newCluster->SetDispersion(newE);
509 fClusterArray->Add(newCluster);
510 }
511
512 fClusterIndices->Set(nCluster);
513}
514//________________________________________________________________________
515void AliAnalysisTaskSOH::ProcessMc()
516{
517 // Process MC.
518
519 for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
520 {
521 AliVParticle* vParticle = fMC->GetTrack(ipart);
522 if(!IsGoodMcParticle(vParticle, ipart)) continue;
523 Int_t pdgCode = vParticle->PdgCode();
524
faba08b8 525 //tracking effciency
526 if(TMath::Abs(vParticle->Eta())<0.9)
4bbf59e5 527 {
9733b37f 528 if(TMath::Abs(pdgCode==211)) fHTrkEffParGenPtEtaPhi->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
4bbf59e5 529 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
530 {
531 AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));
532 if(esdtrack && esdtrack->GetLabel()==ipart)
533 {
faba08b8 534 if(TMath::Abs(pdgCode==211))
535 {
9733b37f 536 fHTrkEffDetGenPtEtaPhi->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
537 fHTrkEffDetRecPtEtaPhi->Fill(esdtrack->Pt(), esdtrack->Eta(), esdtrack->Phi());
faba08b8 538 }
539
540 //cluster E vs. truth photon energy
4bbf59e5 541 for(Int_t k=0; k<fClusterIndices->GetSize(); k++)
542 {
543 AliESDCaloCluster *cluster = fESD->GetCaloCluster(fClusterIndices->At(k));
544 Double_t clsE = cluster->E();
545 TArrayI *MCLabels = cluster->GetLabelsArray();
546 AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
547 AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->GetAt(1));
548
549 if(vParticle1->PdgCode()==22 && vParticle2 == vParticle)
550 {
faba08b8 551 fHPhotonEdiff0HC->Fill(vParticle1->E(), (vParticle1->E() - clsE)/vParticle1->E());
552 fHPhotonEVsClsE->Fill(vParticle1->E(), clsE);
a64777f6 553
554 if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle1->E(), 1);
555 else fHPhotonEdiff30HC->Fill(vParticle1->E(), (vParticle1->E() + 0.3*esdtrack->E() - clsE)/vParticle1->E());
556
557 if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle1->E(), 1);
558 else fHPhotonEdiff70HC->Fill(vParticle1->E(), (vParticle1->E() + 0.7*esdtrack->E() - clsE)/vParticle1->E());
559
4bbf59e5 560 if((clsE - esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle1->E(), 1);
561 else fHPhotonEdiff100HC->Fill(vParticle1->E(), (vParticle1->E() + esdtrack->E() - clsE)/vParticle1->E());
562 continue;
563 }
564 if(vParticle2->PdgCode()==22 && vParticle1 == vParticle)
565 {
faba08b8 566 fHPhotonEdiff0HC->Fill(vParticle2->E(), (vParticle2->E() - clsE)/vParticle2->E());
567 fHPhotonEVsClsE->Fill(vParticle2->E(), clsE);
a64777f6 568
569 if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle2->E(), 1);
570 else fHPhotonEdiff30HC->Fill(vParticle2->E(), (vParticle2->E() + 0.3*esdtrack->E() - clsE)/vParticle2->E());
571
572 if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle2->E(), 1);
573 else fHPhotonEdiff70HC->Fill(vParticle2->E(), (vParticle2->E() + 0.7*esdtrack->E() - clsE)/vParticle2->E());
574
4bbf59e5 575 if((clsE-esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle2->E(), 1);
576 else fHPhotonEdiff100HC->Fill(vParticle2->E(), (vParticle2->E() + esdtrack->E() - clsE)/vParticle2->E());
577 }
578 }
579 break;
580 }
581 }
582 }
583 }
584}
585
586//________________________________________________________________________
587void AliAnalysisTaskSOH::ProcessScaleFactor()
588{
589 // Scale factor.
590
591 const Double_t phiMax = 180 * TMath::DegToRad();
592 const Double_t phiMin = 80 * TMath::DegToRad();
593 const Double_t TPCArea= 2*TMath::Pi()*1.8;
594 const Double_t EMCArea = (phiMax-phiMin)*1.4;
595
596 Double_t PtEMC = 0;
597 Double_t PtTPC = 0;
598
599 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
600 {
601 AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
602 Double_t eta = trk->Eta();
603 Double_t phi = trk->Phi();
604 if(TMath::Abs(eta)<0.9) PtTPC += trk->Pt();
605 if(TMath::Abs(eta)<0.7 && phi > phiMin && phi < phiMax ) PtEMC += trk->Pt();
606 }
607
608 Double_t EtWithHadCorr = 0;
609 Double_t EtWithoutHadCorr = 0;
610 Double_t vertex[3] = {0, 0, 0};
611 fESD->GetVertex()->GetXYZ(vertex);
612 TLorentzVector gamma;
613
614 for(Int_t i=0; i<fClusterArray->GetEntriesFast(); i++)
615 {
616 AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(i);
617 cluster->GetMomentum(gamma, vertex);
618 Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
619 EtWithoutHadCorr += cluster->E() * sinTheta;
620 EtWithHadCorr += cluster->GetDispersion() * sinTheta;
621 }
622
623 if(PtTPC>0)
624 {
625 fHScaleFactor->Fill((PtEMC+EtWithoutHadCorr)/EMCArea * TPCArea/PtTPC);
626 fHScaleFactor100HC->Fill((PtEMC+EtWithHadCorr)/EMCArea * TPCArea/PtTPC);
627 }
628}
629
630//________________________________________________________________________
631AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)
632{
633 // Get accepted track.
634
635 static AliESDtrack newTrack;
636 if(fEsdTrackCuts->AcceptTrack(esdtrack))
637 {
638 esdtrack->Copy(newTrack);
639 newTrack.SetTRDQuality(0);
640 }
641 else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
642 {
643 if(esdtrack->GetConstrainedParam())
644 {
645 esdtrack->Copy(newTrack);
646 const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
647 newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
648 newTrack.SetTRDQuality(1);
649 }
650 else
651 return 0x0;
652 }
653 else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
654 {
655 if(esdtrack->GetConstrainedParam())
656 {
657 esdtrack->Copy(newTrack);
658 const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
659 newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
660 newTrack.SetTRDQuality(2);
661 }
662 else
663 return 0x0;
664 }
665 else
666 {
667 return 0x0;
668 }
669
670 return &newTrack;
671}
672
673//________________________________________________________________________
674Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)
675{
676 // Return true if good MC particle.
677
678 if(!vParticle) return kFALSE;
679 if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
680 if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
681 if(vParticle->Pt()<0.15) return kFALSE;
682 return kTRUE;
683}
684
685//________________________________________________________________________
686Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)
687{
688 // Return true if good cluster.
689
690 if(!cluster) return kFALSE;
691 if (!cluster->IsEMCAL()) return kFALSE;
692 return kTRUE;
693}
694
695//________________________________________________________________________
696void AliAnalysisTaskSOH::Terminate(Option_t *)
697{
698 // Terminate analysis.
699}