]> 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
1) Expend pt range for photons and pi0 from 30 to 40 GeV/c
[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>
048f6758 12#include <THnSparse.h>
4bbf59e5 13
14#include "AliAnalysisManager.h"
15#include "AliAnalysisTask.h"
16#include "AliEMCALRecoUtils.h"
17#include "AliEMCALTrack.h"
18#include "AliESDCaloCluster.h"
19#include "AliESDEvent.h"
20#include "AliESDInputHandler.h"
21#include "AliESDtrack.h"
22#include "AliESDtrackCuts.h"
23#include "AliExternalTrackParam.h"
24#include "AliLog.h"
25#include "AliMCEvent.h"
26
27ClassImp(AliAnalysisTaskSOH)
28
29//________________________________________________________________________
30AliAnalysisTaskSOH::AliAnalysisTaskSOH() :
31 AliAnalysisTaskSE(),
32 fESD(0),
33 fMC(0),
34 fEsdTrackCuts(0x0),
35 fHybridTrackCuts1(0x0),
36 fHybridTrackCuts2(0x0),
37 fTrackIndices(0x0),
38 fClusterIndices(0x0),
39 fClusterArray(0x0),
40 fOutputList(0x0),
41 fHEventStat(0),
9733b37f 42 fHTrkEffParGenPtEtaPhi(0),
43 fHTrkEffDetGenPtEtaPhi(0),
44 fHTrkEffDetRecPtEtaPhi(0),
45 fHTrkEffDetRecFakePtEtaPhi(0),
fb7b5505 46 fHTrkEffParGenPt_all(0),
47 fHTrkEffDetRecPt_all(0),
48 fHTrkEffParGenPt_Dch(0),
49 fHTrkEffDetRecPt_Dch(0),
50 fHTrkEffParGenPt_Dn(0),
51 fHTrkEffDetRecPt_Dn(0),
52 fHTrkEffParGenPt_Ds(0),
53 fHTrkEffDetRecPt_Ds(0),
54 fHTrkEffParGenPt_cL(0),
55 fHTrkEffDetRecPt_cL(0),
56 fHTrkEffParGenPt_JPsi(0),
57 fHTrkEffDetRecPt_JPsi(0),
4bbf59e5 58 fHScaleFactor(0),
59 fHScaleFactor100HC(0),
60 fHEOverPVsPt(0x0),
61 fHEMCalResponsePion(0x0),
62 fHEMCalResponseElec(0x0),
63 fHEMCalResponseProton(0x0),
64 fHEMCalRecdPhidEta(0x0),
65 fHEMCalRecdPhidEtaP(0x0),
66 fHEMCalRecdPhidEtaM(0x0),
67 fHEMCalRecdPhidEta_Truth(0x0),
68 fHEMCalRecdPhidEtaP_Truth(0x0),
69 fHEMCalRecdPhidEtaM_Truth(0x0),
70 fHEMCalRecdPhidEtaposEta(0x0),
71 fHEMCalRecdPhidEtanegEta(0x0),
faba08b8 72 fHPhotonEdiff100HC(0x0),
a64777f6 73 fHPhotonEdiff70HC(0),
74 fHPhotonEdiff30HC(0),
faba08b8 75 fHPhotonEdiff0HC(0x0),
a64777f6 76 fHPhotonEVsClsE(0x0),
77 fHistEsub1Pch(0x0),
78 fHistEsub2Pch(0x0),
79 fHistEsub1PchRat(0x0),
048f6758 80 fHistEsub2PchRat(0x0),
81 fHClsEoverMcE_All(0x0),
82 fHClsEoverMcE_Photon(0x0),
83 fHClsEoverMcE_Elec(0x0),
84 fHClsEoverMcE_Pion(0x0)
4bbf59e5 85{
86 // Constructor
87
4bbf59e5 88}
89
90
91//________________________________________________________________________
92AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :
93 AliAnalysisTaskSE(name),
94 fESD(0),
95 fMC(0),
96 fEsdTrackCuts(0x0),
97 fHybridTrackCuts1(0x0),
98 fHybridTrackCuts2(0x0),
99 fTrackIndices(0x0),
100 fClusterIndices(0x0),
101 fClusterArray(0x0),
102 fOutputList(0x0),
103 fHEventStat(0),
9733b37f 104 fHTrkEffParGenPtEtaPhi(0),
105 fHTrkEffDetGenPtEtaPhi(0),
106 fHTrkEffDetRecPtEtaPhi(0),
107 fHTrkEffDetRecFakePtEtaPhi(0),
fb7b5505 108 fHTrkEffParGenPt_all(0),
109 fHTrkEffDetRecPt_all(0),
110 fHTrkEffParGenPt_Dch(0),
111 fHTrkEffDetRecPt_Dch(0),
112 fHTrkEffParGenPt_Dn(0),
113 fHTrkEffDetRecPt_Dn(0),
114 fHTrkEffParGenPt_Ds(0),
115 fHTrkEffDetRecPt_Ds(0),
116 fHTrkEffParGenPt_cL(0),
117 fHTrkEffDetRecPt_cL(0),
118 fHTrkEffParGenPt_JPsi(0),
119 fHTrkEffDetRecPt_JPsi(0),
4bbf59e5 120 fHScaleFactor(0),
121 fHScaleFactor100HC(0),
122 fHEOverPVsPt(0x0),
af0280a9 123 fHEMCalResponsePion(0x0),
124 fHEMCalResponseElec(0x0),
125 fHEMCalResponseProton(0x0),
4bbf59e5 126 fHEMCalRecdPhidEta(0x0),
127 fHEMCalRecdPhidEtaP(0x0),
128 fHEMCalRecdPhidEtaM(0x0),
129 fHEMCalRecdPhidEta_Truth(0x0),
130 fHEMCalRecdPhidEtaP_Truth(0x0),
131 fHEMCalRecdPhidEtaM_Truth(0x0),
132 fHEMCalRecdPhidEtaposEta(0x0),
133 fHEMCalRecdPhidEtanegEta(0x0),
faba08b8 134 fHPhotonEdiff100HC(0x0),
a64777f6 135 fHPhotonEdiff70HC(0),
136 fHPhotonEdiff30HC(0),
faba08b8 137 fHPhotonEdiff0HC(0x0),
a64777f6 138 fHPhotonEVsClsE(0x0),
139 fHistEsub1Pch(0x0),
140 fHistEsub2Pch(0x0),
141 fHistEsub1PchRat(0x0),
048f6758 142 fHistEsub2PchRat(0x0),
143 fHClsEoverMcE_All(0x0),
144 fHClsEoverMcE_Photon(0x0),
145 fHClsEoverMcE_Elec(0x0),
146 fHClsEoverMcE_Pion(0x0)
4bbf59e5 147{
148 // Constructor
149
150 // Output slot #1 writes into a TH1 container
151 DefineOutput(1, TList::Class());
152}
153
154
155//________________________________________________________________________
156AliAnalysisTaskSOH::~AliAnalysisTaskSOH()
157{
158 // Destructor.
159
160 delete fEsdTrackCuts;
161 delete fHybridTrackCuts1;
162 delete fHybridTrackCuts2;
163 delete fTrackIndices;
164 delete fClusterIndices;
165 delete fClusterArray;
166}
167
168
169//________________________________________________________________________
170void AliAnalysisTaskSOH::UserCreateOutputObjects()
171{
172 // Create histograms, called once.
173
174 OpenFile(1);
175
176 fOutputList = new TList();
177 fOutputList->SetOwner(1);
178
179 fHEventStat = new TH1F("fHEventStat","Event statistics for analysis",8,0,8);
180 fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
181 fHEventStat->GetXaxis()->SetBinLabel(2,"cluster");
182 fHEventStat->GetXaxis()->SetBinLabel(3,"good cluster");
183 fHEventStat->GetXaxis()->SetBinLabel(4,"cls/0-truth");
184 fHEventStat->GetXaxis()->SetBinLabel(5,"cls/1-truth");
185 fHEventStat->GetXaxis()->SetBinLabel(6,"cls/2-truth");
186 fHEventStat->GetXaxis()->SetBinLabel(7,"cls/2-goodtruth");
187 fHEventStat->GetXaxis()->SetBinLabel(8,"cls/>3-truth");
188 fOutputList->Add(fHEventStat);
189
9733b37f 190 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);
191 fHTrkEffParGenPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
192 fHTrkEffParGenPtEtaPhi->GetYaxis()->SetTitle("#eta");
193 fHTrkEffParGenPtEtaPhi->GetZaxis()->SetTitle("#phi");
194 fOutputList->Add(fHTrkEffParGenPtEtaPhi);
195
196 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);
197 fHTrkEffDetGenPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
198 fHTrkEffDetGenPtEtaPhi->GetYaxis()->SetTitle("#eta");
199 fHTrkEffDetGenPtEtaPhi->GetZaxis()->SetTitle("#phi");
200 fOutputList->Add(fHTrkEffDetGenPtEtaPhi);
201
202 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);
203 fHTrkEffDetRecPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
204 fHTrkEffDetRecPtEtaPhi->GetYaxis()->SetTitle("#eta");
205 fHTrkEffDetRecPtEtaPhi->GetZaxis()->SetTitle("#phi");
206 fOutputList->Add(fHTrkEffDetRecPtEtaPhi);
207
208 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);
209 fHTrkEffDetRecFakePtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
210 fHTrkEffDetRecFakePtEtaPhi->GetYaxis()->SetTitle("#eta");
211 fHTrkEffDetRecFakePtEtaPhi->GetZaxis()->SetTitle("#phi");
212 fOutputList->Add(fHTrkEffDetRecFakePtEtaPhi);
de64a00e 213
fb7b5505 214 fHTrkEffParGenPt_all = new TH1F("fHTrkEffParGenPt_all", "Truth p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
215 fOutputList->Add(fHTrkEffParGenPt_all);
216
217 fHTrkEffDetRecPt_all = new TH1F("fHTrkEffDetRecPt_all", "Reconstructed track p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
218 fOutputList->Add(fHTrkEffDetRecPt_all);
219
220 fHTrkEffParGenPt_Dch = new TH1F("fHTrkEffParGenPt_Dch", "Charged D meson truth p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
221 fOutputList->Add(fHTrkEffParGenPt_Dch);
222
223 fHTrkEffDetRecPt_Dch = new TH1F("fHTrkEffDetRecPt_Dch", "Charged D meson reconstructed track p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
224 fOutputList->Add(fHTrkEffDetRecPt_Dch);
225
226 fHTrkEffParGenPt_Dn = new TH1F("fHTrkEffParGenPt_Dn", "Neutral D meson truth p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
227 fOutputList->Add(fHTrkEffParGenPt_Dn);
228
229 fHTrkEffDetRecPt_Dn = new TH1F("fHTrkEffDetRecPt_Dn", "Neutral D meson reconstructed track p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
230 fOutputList->Add(fHTrkEffDetRecPt_Dn);
231
232 fHTrkEffParGenPt_Ds = new TH1F("fHTrkEffParGenPt_Ds", "D_{s} truth p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
233 fOutputList->Add(fHTrkEffParGenPt_Ds);
234
235 fHTrkEffDetRecPt_Ds = new TH1F("fHTrkEffDetRecPt_Ds", "D_{s} reconstructed track p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
236 fOutputList->Add(fHTrkEffDetRecPt_Ds);
237
238 fHTrkEffParGenPt_cL = new TH1F("fHTrkEffParGenPt_cL", "#Lambda_{c} truth p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
239 fOutputList->Add(fHTrkEffParGenPt_cL);
240
241 fHTrkEffDetRecPt_cL = new TH1F("fHTrkEffDetRecPt_cL", "#Lambda_{c} reconstructed track p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
242 fOutputList->Add(fHTrkEffDetRecPt_cL);
243
244 fHTrkEffParGenPt_JPsi = new TH1F("fHTrkEffParGenPt_JPsi", "J/#Psi truth p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
245 fOutputList->Add(fHTrkEffParGenPt_JPsi);
246
247 fHTrkEffDetRecPt_JPsi = new TH1F("fHTrkEffDetRecPt_JPsi", "J/#Psi Reconstructed track p_{T} distribution;p_{T}^{rec} (GeV/c)",500,0.0,50);
248 fOutputList->Add(fHTrkEffDetRecPt_JPsi);
249
4bbf59e5 250 fHScaleFactor = new TH1F("fHScaleFactor", "Scale factor distribution without hadronic correction;Scale factor",100,0,10);
251 fOutputList->Add(fHScaleFactor);
252
253 fHScaleFactor100HC = new TH1F("fHScaleFactor100HC", "Scale factor distribution with 100% hadronic correction;Scale factor",100,0,10);
254 fOutputList->Add(fHScaleFactor100HC);
255
256 fHEOverPVsPt = new TH2F("fHEOverPVsPt", "E/P vs track p_{T}; p_{T} (GeV/c); E/P", 200 , 0, 4, 200, 0, 3.2);
257 fOutputList->Add(fHEOverPVsPt);
258
259 fHEMCalResponsePion = new TH2F("fHEMCalResponsePion", "Pion E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
260 fOutputList->Add(fHEMCalResponsePion);
261
262 fHEMCalResponseElec = new TH2F("fHEMCalResponseElec", "Electron E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
263 fOutputList->Add(fHEMCalResponseElec);
264
265 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 266 fOutputList->Add(fHEMCalResponseProton);
4bbf59e5 267
268 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);
269 fOutputList->Add(fHEMCalRecdPhidEta);
270
271 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);
272 fOutputList->Add(fHEMCalRecdPhidEtaP);
273
274 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);
275 fOutputList->Add(fHEMCalRecdPhidEtaM);
276
277 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);
278 fOutputList->Add(fHEMCalRecdPhidEta_Truth);
279
280 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);
281 fOutputList->Add(fHEMCalRecdPhidEtaP_Truth);
282
283 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);
284 fOutputList->Add(fHEMCalRecdPhidEtaM_Truth);
285
286 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);
287 fOutputList->Add(fHEMCalRecdPhidEtaposEta);
288
289 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);
290 fOutputList->Add(fHEMCalRecdPhidEtanegEta);
291
a64777f6 292 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 293 fOutputList->Add(fHPhotonEdiff100HC);
294
a64777f6 295 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);
296 fOutputList->Add(fHPhotonEdiff70HC);
297
298 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);
299 fOutputList->Add(fHPhotonEdiff30HC);
300
301 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 302 fOutputList->Add(fHPhotonEdiff0HC);
303
304 fHPhotonEVsClsE = new TH2F("fHPhotonEVsClsE","Cluster E vs. photon E_{Truth}; photon E_{Truth} (GeV); Cluster E (GeV)",500,0,5,500,0,5);
305 fOutputList->Add(fHPhotonEVsClsE);
306
a64777f6 307 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.);
308 fOutputList->Add(fHistEsub1Pch);
309
310 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.);
311 fOutputList->Add(fHistEsub2Pch);
312
313 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);
314 fOutputList->Add(fHistEsub1PchRat);
315
316 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);
317 fOutputList->Add(fHistEsub2PchRat);
318
048f6758 319 Int_t bins[4] = {150, 150, 100, 200};
320 Double_t xmin[4] = {1.3, -0.8, 0, 0};
321 Double_t xmax[4] = {3.2, 0.8, 10, 2};
322
323 fHClsEoverMcE_All = new THnSparseF("fHClsEoverMcE_All", "Cluster E/MC E, clusters with 1 matching particle; #phi; #eta; E (GeV); ClsE/McE", 4, bins, xmin, xmax);
324 fOutputList->Add(fHClsEoverMcE_All);
325
326 fHClsEoverMcE_Photon = new THnSparseF("fHClsEoverMcE_Photon", "Cluster E/MC E, clusters with 1 matching photon; #phi; #eta; E (GeV); ClsE/McE", 4, bins, xmin, xmax);
327 fOutputList->Add(fHClsEoverMcE_Photon);
328
329 fHClsEoverMcE_Elec = new THnSparseF("fHClsEoverMcE_Elec", "Cluster E/MC E, clusters with 1 matching electron; #phi; #eta; E (GeV); ClsE/McE", 4, bins, xmin, xmax);
330 fOutputList->Add(fHClsEoverMcE_Elec);
331
332 fHClsEoverMcE_Pion = new THnSparseF("fHClsEoverMcE_Pion", "Cluster E/MC E, clusters with 1 matching pion; #phi; #eta; E (GeV); ClsE/McE",4, bins, xmin, xmax);
333 fOutputList->Add(fHClsEoverMcE_Pion);
334
4bbf59e5 335 fTrackIndices = new TArrayI();
336 fClusterIndices = new TArrayI();
337
338 fClusterArray = new TObjArray();
339 fClusterArray->SetOwner(1);
340
341 PostData(1, fOutputList);
342}
343
344//________________________________________________________________________
345void AliAnalysisTaskSOH::UserExec(Option_t *)
346{
347 // Main loop, called for each event.
348
349 // Post output data.
350 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
351 if (!fESD) {
352 printf("ERROR: fESD not available\n");
353 return;
354 }
355
356 fMC = MCEvent();
357 if (!fMC) {
358 printf("ERROR: fMC not available\n");
359 return;
360 }
361
362 fHEventStat->Fill(0.5);
363
364 if(fTrackIndices)
365 fTrackIndices->Reset();
366 if(fClusterIndices)
367 fClusterIndices->Reset();
368 if(fClusterArray)
369 fClusterArray->Delete();
370
371 ProcessTrack();
372 ProcessCluster();
373 ProcessMc();
374 ProcessScaleFactor();
375
376 PostData(1, fOutputList);
377}
378
379//________________________________________________________________________
380void AliAnalysisTaskSOH::ProcessTrack()
381{
382 // Process track.
383
384 fTrackIndices->Set(fESD->GetNumberOfTracks());
385 AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));
386
387 Int_t isMth = 0;
388 Int_t nTracks = 0;
389
390 Float_t ClsPos[3] = {-999,-999,-999};
391 Double_t emcTrkpos[3] = {-999,-999,-999};
392
393 for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)
394 {
395 AliESDtrack *esdtrack = fESD->GetTrack(itr);
396 if(!esdtrack)continue;
397 AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
398 if(!newTrack) continue;
399
400 Double_t clsE = -1;
401 Int_t clsIndex = newTrack->GetEMCALcluster();
402 if(newTrack->GetEMCALcluster()>-1)
403 {
404 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);
405 if(IsGoodCluster(cluster))
406 {
407 isMth=1;
408
409 cluster->GetPosition(ClsPos);
410 TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
411
412 AliEMCALTrack EMCTrk(*newTrack);
413 if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {
414 continue;
415 }
416 EMCTrk.GetXYZ(emcTrkpos);
417 TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);
418
419 Double_t dPhi = VClsPos.Phi() - VemcTrkPos.Phi();
420 if (dPhi < -1*TMath::Pi()) dPhi += (2*TMath::Pi());
421 else if (dPhi > TMath::Pi()) dPhi -= (2*TMath::Pi());
422
423 Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();
424
425 fHEMCalRecdPhidEta->Fill(dEta, dPhi);
426
427 if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())
428 {
429 AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());
430 if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))
431 {
432 fHEMCalRecdPhidEta_Truth->Fill(dEta, dPhi);
433 if(vParticle->Charge() > 0) fHEMCalRecdPhidEtaP_Truth->Fill(dEta, dPhi);
434 if(vParticle->Charge() < 0) fHEMCalRecdPhidEtaM_Truth->Fill(dEta, dPhi);
435 }
436 }
437
438 if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}
439 if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}
440
441 if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);
442 if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);
443
444 clsE = cluster->E();
445 if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());
446 }
447
448 Int_t ipart = newTrack->GetLabel();
449 if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
450 {
451 AliVParticle* vParticle = fMC->GetTrack(ipart);
452 if(isMth && vParticle)
453 {
454 if(TMath::Abs(vParticle->PdgCode())==211)
455 {
456 fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());
457 }
458 if(TMath::Abs(vParticle->PdgCode())==11)
459 {
460 fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());
461 }
462 if(TMath::Abs(vParticle->PdgCode())==2212)
463 {
464 fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());
465 }
466 }
467 }
468 }
fb7b5505 469 if(newTrack->GetLabel()<0 && newTrack->GetPID()==2) fHTrkEffDetRecFakePtEtaPhi->Fill(esdtrack->Pt(), esdtrack->Eta(), esdtrack->Phi());
de64a00e 470
de64a00e 471 // Track Indices
4bbf59e5 472 fTrackIndices->AddAt(itr,nTracks);
473 nTracks++;
474 }
475
476 fTrackIndices->Set(nTracks);
477}
478
479//________________________________________________________________________
480void AliAnalysisTaskSOH::ProcessCluster()
481{
482 // Process cluster.
483
484 Int_t nCluster = 0;
485 TLorentzVector gamma;
486 Double_t vertex[3] = {0, 0, 0};
487 fESD->GetVertex()->GetXYZ(vertex);
488 const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
489 fClusterIndices->Set(nCaloClusters);
048f6758 490 Float_t ClsPos[3] = {-999,-999,-999};
4bbf59e5 491
492 for(Int_t itr=0; itr<nCaloClusters; itr++)
493 {
494 fHEventStat->Fill(1.5);
495 AliESDCaloCluster *cluster = fESD->GetCaloCluster(itr);
496 if(!IsGoodCluster(cluster)) continue;
497 cluster->GetMomentum(gamma, vertex);
498 if (gamma.Pt() < 0.15) continue;
499 fHEventStat->Fill(2.5);
500
048f6758 501 cluster->GetPosition(ClsPos);
502 TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
503
a64777f6 504 TArrayI *TrackLabels = cluster->GetTracksMatched();
505
506 if(TrackLabels->GetSize() == 1)
507 {
508 AliESDtrack *esdtrack = fESD->GetTrack(TrackLabels->GetAt(0));
509 AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
510 if(newTrack && TMath::Abs(newTrack->Eta())<0.7)
511 {
512 Double_t Esub = newTrack->P();
513 if (Esub > cluster->E()) Esub = cluster->E();
514 fHistEsub1Pch->Fill(newTrack->P(), Esub);
515 fHistEsub1PchRat->Fill(newTrack->P(), Esub/newTrack->P());
516 }
517 }
518
519 if(TrackLabels->GetSize() == 2)
520 {
521 AliESDtrack *esdtrack1 = fESD->GetTrack(TrackLabels->GetAt(0));
522 AliESDtrack *esdtrack2 = fESD->GetTrack(TrackLabels->GetAt(1));
523 AliESDtrack *newTrack1 = GetAcceptTrack(esdtrack1);
524 AliESDtrack *newTrack2 = GetAcceptTrack(esdtrack2);
525 if(newTrack1 && newTrack2 && TMath::Abs(newTrack1->Eta())<0.7 && TMath::Abs(newTrack2->Eta())<0.7)
526 {
527 Double_t Esub = newTrack1->P() + newTrack2->P();
528 if (Esub > cluster->E()) Esub = cluster->E();
529 fHistEsub2Pch->Fill(newTrack1->P() + newTrack2->P(), Esub);
530 fHistEsub2PchRat->Fill(newTrack1->P() + newTrack2->P(), Esub/(newTrack1->P() + newTrack2->P()));
531 }
532 else if(newTrack1 && !(newTrack2) && TMath::Abs(newTrack1->Eta())<0.7)
533 {
534 Double_t Esub = newTrack1->P();
535 if (Esub > cluster->E()) Esub = cluster->E();
536 fHistEsub1Pch->Fill(newTrack1->P(), Esub);
537 fHistEsub1PchRat->Fill(newTrack1->P(), Esub/newTrack1->P());
538 }
539 else if (!(newTrack1) && newTrack2 && TMath::Abs(newTrack2->Eta())<0.7)
540 {
541 Double_t Esub = newTrack2->P();
542 if (Esub > cluster->E()) Esub = cluster->E();
543 fHistEsub1Pch->Fill(newTrack2->P(), Esub);
544 fHistEsub1PchRat->Fill(newTrack2->P(), Esub/newTrack2->P());
545 }
546 else {;}
547 }
548
4bbf59e5 549 TArrayI *MCLabels = cluster->GetLabelsArray();
550
551 if(MCLabels->GetSize() == 0) fHEventStat->Fill(3.5);
048f6758 552 if(MCLabels->GetSize() == 1)
553 {
554 fHEventStat->Fill(4.5);
555 AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
556 if(IsGoodMcParticle(vParticle1, MCLabels->GetAt(0)))
557 {
558 Double_t Entries[4] = {VClsPos.Phi(), VClsPos.Eta(), vParticle1->E(), cluster->E()/vParticle1->E()};
559 fHClsEoverMcE_All->Fill(Entries);
560 if(vParticle1->PdgCode() == 22)
561 {
562 fHClsEoverMcE_Photon->Fill(Entries);
563 }
564 if(TMath::Abs(vParticle1->PdgCode()) == 11)
565 {
566 fHClsEoverMcE_Elec->Fill(Entries);
567 }
568 if(TMath::Abs(vParticle1->PdgCode()) == 211)
569 {
570 fHClsEoverMcE_Pion->Fill(Entries);
571 }
572 }
573 }
4bbf59e5 574 if(MCLabels->GetSize() == 2)
575 {
576 fHEventStat->Fill(5.5);
577 AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
578 AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->GetAt(1));
579 if(IsGoodMcParticle(vParticle1, MCLabels->GetAt(0)) && IsGoodMcParticle(vParticle2, MCLabels->GetAt(1)))
580 {
581 fHEventStat->Fill(6.5);
582 if((vParticle1->PdgCode()==22) && (vParticle2->PdgCode()==22)) {;}
583 else if((vParticle1->PdgCode()!=22) && (vParticle2->PdgCode()!=22)) {;}
584 else
585 {
586 fClusterIndices->AddAt(itr,nCluster);
587 nCluster++;
588 }
589 }
590 }
591 if(MCLabels->GetSize() > 2) fHEventStat->Fill(7.5);
592
593 AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cluster);
594
595 Double_t subE = 0;
596 TArrayI arrayTrackMatched(fTrackIndices->GetSize());
597 Int_t nGoodMatch = 0;
598
599 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
600 {
601 AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
602 if(itr==trk->GetEMCALcluster())
603 {
604 arrayTrackMatched[nGoodMatch] = j;
605 nGoodMatch ++;
606 subE += trk->P();
607 }
608 }
609
610 arrayTrackMatched.Set(nGoodMatch);
611 newCluster->AddTracksMatched(arrayTrackMatched);
612
613 Double_t clsE = newCluster->E();
614 Double_t newE = clsE-subE;
615 if(newE<0) newE = 0;
616 newCluster->SetDispersion(newE);
617 fClusterArray->Add(newCluster);
618 }
619
620 fClusterIndices->Set(nCluster);
621}
622//________________________________________________________________________
623void AliAnalysisTaskSOH::ProcessMc()
624{
625 // Process MC.
626
048f6758 627 for(Int_t i=0; i<fTrackIndices->GetSize(); i++)
628 {
629 AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(i));
630 Int_t index = esdtrack->GetLabel();
631 if(index < 0 || index > fESD->GetNumberOfTracks()) continue;
632 AliVParticle* vParticle = fMC->GetTrack(index);
633 if(!IsGoodMcParticle(vParticle, index) && esdtrack->GetPID() == 2) fHTrkEffDetRecFakePtEtaPhi->Fill(esdtrack->Pt(), esdtrack->Eta(), esdtrack->Phi());
634 }
635
fb7b5505 636 //tracking effciency
4bbf59e5 637 for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
638 {
639 AliVParticle* vParticle = fMC->GetTrack(ipart);
4bbf59e5 640 Int_t pdgCode = vParticle->PdgCode();
fb7b5505 641 AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(ipart);
642 if(!IsGoodMcParticle(vParticle, ipart)) continue;
643
644 if(McParticle->GetMother() > 0 && TMath::Abs(vParticle->Eta())<0.9 && TMath::Abs(pdgCode)==211)
645 {
646 AliMCParticle *McMother = (AliMCParticle*)fMC->GetTrack(McParticle->GetMother());
647 if(TMath::Abs(McMother->PdgCode()) == 411) fHTrkEffParGenPt_Dch->Fill(vParticle->Pt());
648 else if(TMath::Abs(McMother->PdgCode()) == 421) fHTrkEffParGenPt_Dn->Fill(vParticle->Pt());
649 else if(TMath::Abs(McMother->PdgCode()) == 431) fHTrkEffParGenPt_Ds->Fill(vParticle->Pt());
650 else if(TMath::Abs(McMother->PdgCode()) == 4122) fHTrkEffParGenPt_cL->Fill(vParticle->Pt());
651 else if(TMath::Abs(McMother->PdgCode()) == 443) fHTrkEffParGenPt_JPsi->Fill(vParticle->Pt());
652 else {;}
653 }
654
faba08b8 655 if(TMath::Abs(vParticle->Eta())<0.9)
4bbf59e5 656 {
fb7b5505 657 fHTrkEffParGenPt_all->Fill(vParticle->Pt());
9733b37f 658 if(TMath::Abs(pdgCode==211)) fHTrkEffParGenPtEtaPhi->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
fb7b5505 659
4bbf59e5 660 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
661 {
662 AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));
663 if(esdtrack && esdtrack->GetLabel()==ipart)
664 {
fb7b5505 665 fHTrkEffDetRecPt_all->Fill(esdtrack->Pt());
faba08b8 666 if(TMath::Abs(pdgCode==211))
667 {
9733b37f 668 fHTrkEffDetGenPtEtaPhi->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
669 fHTrkEffDetRecPtEtaPhi->Fill(esdtrack->Pt(), esdtrack->Eta(), esdtrack->Phi());
faba08b8 670 }
fb7b5505 671 if(McParticle->GetMother() > 0 && TMath::Abs(pdgCode) == 211)
672 {
673 AliMCParticle *McMother = (AliMCParticle*)fMC->GetTrack(McParticle->GetMother());
674 if(TMath::Abs(McMother->PdgCode()) == 411) fHTrkEffDetRecPt_Dch->Fill(esdtrack->Pt());
675 else if(TMath::Abs(McMother->PdgCode()) == 421) fHTrkEffDetRecPt_Dn->Fill(esdtrack->Pt());
676 else if(TMath::Abs(McMother->PdgCode()) == 431) fHTrkEffDetRecPt_Ds->Fill(esdtrack->Pt());
677 else if(TMath::Abs(McMother->PdgCode()) == 4122) fHTrkEffDetRecPt_cL->Fill(esdtrack->Pt());
678 else if(TMath::Abs(McMother->PdgCode()) == 443) fHTrkEffDetRecPt_JPsi->Fill(esdtrack->Pt());
679 else {;}
680 }
faba08b8 681 //cluster E vs. truth photon energy
4bbf59e5 682 for(Int_t k=0; k<fClusterIndices->GetSize(); k++)
683 {
684 AliESDCaloCluster *cluster = fESD->GetCaloCluster(fClusterIndices->At(k));
685 Double_t clsE = cluster->E();
686 TArrayI *MCLabels = cluster->GetLabelsArray();
687 AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
688 AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->GetAt(1));
689
690 if(vParticle1->PdgCode()==22 && vParticle2 == vParticle)
691 {
faba08b8 692 fHPhotonEdiff0HC->Fill(vParticle1->E(), (vParticle1->E() - clsE)/vParticle1->E());
693 fHPhotonEVsClsE->Fill(vParticle1->E(), clsE);
a64777f6 694
695 if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle1->E(), 1);
696 else fHPhotonEdiff30HC->Fill(vParticle1->E(), (vParticle1->E() + 0.3*esdtrack->E() - clsE)/vParticle1->E());
697
698 if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle1->E(), 1);
699 else fHPhotonEdiff70HC->Fill(vParticle1->E(), (vParticle1->E() + 0.7*esdtrack->E() - clsE)/vParticle1->E());
700
4bbf59e5 701 if((clsE - esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle1->E(), 1);
702 else fHPhotonEdiff100HC->Fill(vParticle1->E(), (vParticle1->E() + esdtrack->E() - clsE)/vParticle1->E());
703 continue;
704 }
705 if(vParticle2->PdgCode()==22 && vParticle1 == vParticle)
706 {
faba08b8 707 fHPhotonEdiff0HC->Fill(vParticle2->E(), (vParticle2->E() - clsE)/vParticle2->E());
708 fHPhotonEVsClsE->Fill(vParticle2->E(), clsE);
a64777f6 709
710 if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle2->E(), 1);
711 else fHPhotonEdiff30HC->Fill(vParticle2->E(), (vParticle2->E() + 0.3*esdtrack->E() - clsE)/vParticle2->E());
712
713 if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle2->E(), 1);
714 else fHPhotonEdiff70HC->Fill(vParticle2->E(), (vParticle2->E() + 0.7*esdtrack->E() - clsE)/vParticle2->E());
715
4bbf59e5 716 if((clsE-esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle2->E(), 1);
717 else fHPhotonEdiff100HC->Fill(vParticle2->E(), (vParticle2->E() + esdtrack->E() - clsE)/vParticle2->E());
718 }
719 }
720 break;
721 }
722 }
723 }
724 }
725}
726
727//________________________________________________________________________
728void AliAnalysisTaskSOH::ProcessScaleFactor()
729{
730 // Scale factor.
731
732 const Double_t phiMax = 180 * TMath::DegToRad();
733 const Double_t phiMin = 80 * TMath::DegToRad();
734 const Double_t TPCArea= 2*TMath::Pi()*1.8;
735 const Double_t EMCArea = (phiMax-phiMin)*1.4;
736
737 Double_t PtEMC = 0;
738 Double_t PtTPC = 0;
739
740 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
741 {
742 AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
743 Double_t eta = trk->Eta();
744 Double_t phi = trk->Phi();
745 if(TMath::Abs(eta)<0.9) PtTPC += trk->Pt();
746 if(TMath::Abs(eta)<0.7 && phi > phiMin && phi < phiMax ) PtEMC += trk->Pt();
747 }
748
749 Double_t EtWithHadCorr = 0;
750 Double_t EtWithoutHadCorr = 0;
751 Double_t vertex[3] = {0, 0, 0};
752 fESD->GetVertex()->GetXYZ(vertex);
753 TLorentzVector gamma;
754
755 for(Int_t i=0; i<fClusterArray->GetEntriesFast(); i++)
756 {
757 AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(i);
758 cluster->GetMomentum(gamma, vertex);
759 Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
760 EtWithoutHadCorr += cluster->E() * sinTheta;
761 EtWithHadCorr += cluster->GetDispersion() * sinTheta;
762 }
763
764 if(PtTPC>0)
765 {
766 fHScaleFactor->Fill((PtEMC+EtWithoutHadCorr)/EMCArea * TPCArea/PtTPC);
767 fHScaleFactor100HC->Fill((PtEMC+EtWithHadCorr)/EMCArea * TPCArea/PtTPC);
768 }
769}
770
771//________________________________________________________________________
772AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)
773{
774 // Get accepted track.
775
776 static AliESDtrack newTrack;
777 if(fEsdTrackCuts->AcceptTrack(esdtrack))
778 {
779 esdtrack->Copy(newTrack);
780 newTrack.SetTRDQuality(0);
781 }
782 else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
783 {
784 if(esdtrack->GetConstrainedParam())
785 {
786 esdtrack->Copy(newTrack);
787 const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
788 newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
789 newTrack.SetTRDQuality(1);
790 }
791 else
792 return 0x0;
793 }
794 else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
795 {
796 if(esdtrack->GetConstrainedParam())
797 {
798 esdtrack->Copy(newTrack);
799 const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
800 newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
801 newTrack.SetTRDQuality(2);
802 }
803 else
804 return 0x0;
805 }
806 else
807 {
808 return 0x0;
809 }
810
811 return &newTrack;
812}
813
814//________________________________________________________________________
815Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)
816{
817 // Return true if good MC particle.
818
819 if(!vParticle) return kFALSE;
820 if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
821 if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
822 if(vParticle->Pt()<0.15) return kFALSE;
823 return kTRUE;
824}
825
826//________________________________________________________________________
827Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)
828{
829 // Return true if good cluster.
830
831 if(!cluster) return kFALSE;
832 if (!cluster->IsEMCAL()) return kFALSE;
833 return kTRUE;
834}
835
836//________________________________________________________________________
837void AliAnalysisTaskSOH::Terminate(Option_t *)
838{
839 // Terminate analysis.
840}