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