]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0V2ShSh.cxx
New version for photon isolation task
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALPi0V2ShSh.cxx
CommitLineData
fc0f5568 1// $Id: AliAnalysisTaskEMCALPi0V2ShSh.cxx$
2
3#include "AliAnalysisTaskEMCALPi0V2ShSh.h"
4
5//Root include files
6//#include <Riostream.h>
7#include <TParticle.h>
8#include <TRefArray.h>
9#include <TCanvas.h>
10#include <TChain.h>
11#include <TFile.h>
12#include <TH1F.h>
13#include <TH1D.h>
14#include <TH2F.h>
15#include <THnSparse.h>
16#include <TList.h>
17#include <TMath.h>
18#include <TVirtualFFT.h>
19
20//AliRoot include files
21#include "AliAnalysisTaskSE.h"
22#include "AliRunLoader.h"
23#include "AliAnalysisManager.h"
24#include "AliAnalysisTask.h"
25#include "AliStack.h"
c4a4ee27 26#include "AliEMCALGeometry.h"
fc0f5568 27#include "AliESDEvent.h"
28#include "AliESDVertex.h"
c4a4ee27 29#include "AliESDCaloCells.h"
fc0f5568 30#include "AliESDCaloCluster.h"
31#include "AliESDEvent.h"
32#include "AliESDHeader.h"
33#include "AliESDInputHandler.h"
34#include "AliESDtrack.h"
35#include "AliKFParticle.h"
c4a4ee27 36#include "AliAODEvent.h"
fc0f5568 37#include "AliVCluster.h"
38#include "AliCentrality.h"
39#include "AliEventplane.h"
c4a4ee27 40#include "AliOADBContainer.h"
fc0f5568 41
42
43
44ClassImp(AliAnalysisTaskEMCALPi0V2ShSh)
45
46//________________________________________________________________________
47AliAnalysisTaskEMCALPi0V2ShSh::AliAnalysisTaskEMCALPi0V2ShSh() :
48 AliAnalysisTaskSE(),
c4a4ee27 49 fEventPlane(0),
50 fCentralityV0M(99.),
51 fESDClusters(0),
52 fAODClusters(0),
c4a4ee27 53 fESDCells(0),
54 fAODCells(0),
55 fGeom(0),
56 fGeoName("EMCAL_COMPLETEV1"),
57 fOADBContainer(0),
58 fESD(0),
59 fAOD(0),
60 fOutputList(0),
61 fEPTPC(-999.),
62 fEPTPCResolution(0.),
63 fEPV0(-999.),
64 fEPV0A(-999.),
65 fEPV0C(-999.),
66 fEPV0Ar(-999.),
67 fEPV0Cr(-999.),
68 fEPV0r(-999.),
69 fEPV0A4r(-999.),
70 fEPV0A5r(-999.),
71 fEPV0A6r(-999.),
72 fEPV0A7r(-999.),
73 fEPV0C0r(-999.),
74 fEPV0C1r(-999.),
75 fEPV0C2r(-999.),
76 fEPV0C3r(-999.),
77 fHistAllcentV0(0),
78 fHistAllcentV0r(0),
79 fHistAllcentV0A(0),
80 fHistAllcentV0C(0),
81 fHistAllcentTPC(0),
82 fHistEPTPC(0),
83 fHistEPTPCResolution(0),
84 fHistEPV0(0),
85 fHistEPV0A(0),
86 fHistEPV0C(0),
87 fHistEPV0Ar(0),
88 fHistEPV0Cr(0),
89 fHistEPV0r(0),
90 fHistEPV0A4r(0),
91 fHistEPV0A7r(0),
92 fHistEPV0C0r(0),
93 fHistEPV0C3r(0),
94 fHistdifV0A_V0C0r(0),
95 fHistdifV0A_V0C3r(0),
96 fHistdifV0C0r_V0C3r(0),
97 fHistdifV0C_V0A4r(0),
98 fHistdifV0C_V0A7r(0),
99 fHistdifV0A4r_V0A7r(0),
100 fHistdifV0Ar_V0Cr(0),
101 fHistClusterEta(0),
102 fHistClusterPhi(0),
103 fHistClusterE(0),
104 fHistClusterEt(0),
105 fHistClusterN(0),
106 fHistClusterM02(0),
107 fHistClusterEN(0),
108 fHistClusterEM02(0),
109 fHistClusterPhiEta(0),
110 fHistClusterEtN(0),
111 fHistClusterEtM02(0),
112 fHistClusterdphiV0(0),
113 fHistTrackPt(0),
114 fHistTrackEta(0),
115 fHistTrackPhi(0),
116 fHistTrackPhiEta(0),
117 fClusterPbV0(0),
118 fClusterPbV0A(0),
119 fClusterPbV0C(0),
120 fClusterPbTPC(0)
fc0f5568 121{
122 // Default constructor.
c4a4ee27 123 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
fc0f5568 124}
125
126//________________________________________________________________________
127AliAnalysisTaskEMCALPi0V2ShSh::AliAnalysisTaskEMCALPi0V2ShSh(const char *name) :
128 AliAnalysisTaskSE(name),
c4a4ee27 129 fEventPlane(0),
130 fCentralityV0M(99.),
131 fESDClusters(0),
132 fAODClusters(0),
c4a4ee27 133 fESDCells(0),
134 fAODCells(0),
135 fGeom(0),
136 fGeoName("EMCAL_COMPLETEV1"),
137 fOADBContainer(0),
138 fESD(0),
139 fAOD(0),
140 fOutputList(0),
141 fEPTPC(-999.),
142 fEPTPCResolution(0.),
143 fEPV0(-999.),
144 fEPV0A(-999.),
145 fEPV0C(-999.),
146 fEPV0Ar(-999.),
147 fEPV0Cr(-999.),
148 fEPV0r(-999.),
149 fEPV0A4r(-999.),
150 fEPV0A5r(-999.),
151 fEPV0A6r(-999.),
152 fEPV0A7r(-999.),
153 fEPV0C0r(-999.),
154 fEPV0C1r(-999.),
155 fEPV0C2r(-999.),
156 fEPV0C3r(-999.),
157 fHistAllcentV0(0),
158 fHistAllcentV0r(0),
159 fHistAllcentV0A(0),
160 fHistAllcentV0C(0),
161 fHistAllcentTPC(0),
162 fHistEPTPC(0),
163 fHistEPTPCResolution(0),
164 fHistEPV0(0),
165 fHistEPV0A(0),
166 fHistEPV0C(0),
167 fHistEPV0Ar(0),
168 fHistEPV0Cr(0),
169 fHistEPV0r(0),
170 fHistEPV0A4r(0),
171 fHistEPV0A7r(0),
172 fHistEPV0C0r(0),
173 fHistEPV0C3r(0),
174 fHistdifV0A_V0C0r(0),
175 fHistdifV0A_V0C3r(0),
176 fHistdifV0C0r_V0C3r(0),
177 fHistdifV0C_V0A4r(0),
178 fHistdifV0C_V0A7r(0),
179 fHistdifV0A4r_V0A7r(0),
180 fHistdifV0Ar_V0Cr(0),
181 fHistClusterEta(0),
182 fHistClusterPhi(0),
183 fHistClusterE(0),
184 fHistClusterEt(0),
185 fHistClusterN(0),
186 fHistClusterM02(0),
187 fHistClusterEN(0),
188 fHistClusterEM02(0),
189 fHistClusterPhiEta(0),
190 fHistClusterEtN(0),
191 fHistClusterEtM02(0),
192 fHistClusterdphiV0(0),
193 fHistTrackPt(0),
194 fHistTrackEta(0),
195 fHistTrackPhi(0),
196 fHistTrackPhiEta(0),
197 fClusterPbV0(0),
198 fClusterPbV0A(0),
199 fClusterPbV0C(0),
200 fClusterPbTPC(0)
fc0f5568 201{
202 // Constructor
203 // Define input and output slots here
204 // Input slot #0 works with a TChain
205 DefineInput(0, TChain::Class());
206 // Output slot #0 id reserved by the base class for AOD
207 // Output slot #1 writes into a TH1 container
208 DefineOutput(1, TList::Class());
209}
210
211//________________________________________________________________________
212void AliAnalysisTaskEMCALPi0V2ShSh::UserCreateOutputObjects()
213{
214 // Create histograms, called once.
215
c4a4ee27 216 fESDClusters = new TObjArray();
217 fAODClusters = new TObjArray();
c4a4ee27 218 fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
219 fOADBContainer = new AliOADBContainer("AliEMCALgeo");
220 fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
fc0f5568 221
222 fOutputList = new TList();
223 fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
224
225 fHistTrackPt = new TH1F("fHistTrackPt","Track Transverse Momentum Distribution of Pb+Pb",100,0.0,30.0);
226 fHistTrackPt->GetYaxis()->SetTitle("Entries"); fHistTrackPt->GetXaxis()->SetTitle("P_{t} [GeV/c]");
227 fOutputList->Add(fHistTrackPt);
228
229 fHistTrackEta = new TH1F("fHistTrackEta","Track Pseudorapidity Distribution of Pb+Pb",100,-1.0,1.0);
230 fHistTrackEta->GetYaxis()->SetTitle("Entries"); fHistTrackEta->GetXaxis()->SetTitle("#eta");
231 fOutputList->Add(fHistTrackEta);
232
233 fHistTrackPhi = new TH1F("fHistTrackPhi","Track #phi Distribution of Pb+Pb",100,0.0,6.29);
234 fHistTrackPhi->GetYaxis()->SetTitle("Entries"); fHistTrackPhi->GetXaxis()->SetTitle("#phi [rad]");
235 fOutputList->Add(fHistTrackPhi);
236
237 fHistTrackPhiEta = new TH2F("fHistTrackPhiEta","Track Pseudorapidity vs #phi of Pb+Pb",100,-1.0,1.0,100,0.0,6.29);
238 fHistTrackPhiEta->GetXaxis()->SetTitle("#eta"); fHistTrackPhiEta->GetYaxis()->SetTitle("#phi [rad]");
239 fOutputList->Add(fHistTrackPhiEta);
240
241 fHistClusterEta = new TH1F("fHistClusterEta","Cluster Pseudorapidity Distribution of Pb+Pb",100,-1.0,1.0);
242 fHistClusterEta->GetYaxis()->SetTitle("Entries"); fHistClusterEta->GetXaxis()->SetTitle("#eta");
243 fOutputList->Add(fHistClusterEta);
244
245 fHistClusterPhi = new TH1F("fHistClusterPhi","Cluster #phi Distribution of Pb+Pb",100,0.0,6.29);
246 fHistClusterPhi->GetYaxis()->SetTitle("Entries"); fHistClusterPhi->GetXaxis()->SetTitle("#phi [rad]");
247 fOutputList->Add(fHistClusterPhi);
248
249 fHistClusterPhiEta = new TH2F("fHistClusterPhiEta","Cluster Pseudorapidity vs #phi of Pb+Pb",100,-1.0,1.0,100,0.0,6.29);
250 fHistClusterPhiEta->GetXaxis()->SetTitle("#eta"); fHistClusterPhiEta->GetYaxis()->SetTitle("#phi [rad]");
251 fOutputList->Add(fHistClusterPhiEta);
252
253 fHistClusterM02 = new TH1F("fHistClusterM02","Cluster M02 Distribution of Pb+Pb",100,0.0,3.0);
254 fHistClusterM02->GetYaxis()->SetTitle("Entries"); fHistClusterM02->GetXaxis()->SetTitle("M02");
255 fOutputList->Add(fHistClusterM02);
256
257 fHistClusterE = new TH1F("fHistClusterE","Cluster Energy Distribution of Pb+Pb",100,0.0,20.0);
258 fHistClusterE->GetYaxis()->SetTitle("Entries"); fHistClusterE->GetXaxis()->SetTitle("Energy [GeV]");
259 fOutputList->Add(fHistClusterE);
260
261 fHistClusterEt = new TH1F("fHistClusterEt","Cluster Transverse Energy Distribution of Pb+Pb",100,0.0,20.0);
262 fHistClusterEt->GetYaxis()->SetTitle("Entries"); fHistClusterEt->GetXaxis()->SetTitle("Transverse Energy [GeV]");
263 fOutputList->Add(fHistClusterEt);
264
265 fHistClusterEM02 = new TH2F("fHistClusterEM02","Cluster Energy vs M02 of Pb+Pb",100,0.0,20.0,100,0.0,3.0);
266 fHistClusterEM02->GetYaxis()->SetTitle("M02"); fHistClusterEM02->GetXaxis()->SetTitle("Energy [GeV]");
267 fOutputList->Add(fHistClusterEM02);
268
269 fHistClusterEtM02 = new TH2F("fHistClusterEtM02","Cluster Transverse Energy vs M02 of Pb+Pb",100,0.0,20.0,100,0.0,3.0);
270 fHistClusterEtM02->GetYaxis()->SetTitle("M02"); fHistClusterEtM02->GetXaxis()->SetTitle("Transverse Energy [GeV]");
271 fOutputList->Add(fHistClusterEtM02);
272
273 fHistClusterN = new TH1F("fHistClusterN","Cluster N Distribution of Pb+Pb",30,0.0,30.0);
274 fHistClusterN->GetYaxis()->SetTitle("Entries"); fHistClusterN->GetXaxis()->SetTitle("N");
275 fOutputList->Add(fHistClusterN);
276
277 fHistClusterEN = new TH2F("fHistClusterEN","N vs Cluster Energy of Pb+Pb",100,0.0,20.0,30,0.0,30.0);
278 fHistClusterEN->GetYaxis()->SetTitle("N"); fHistClusterEN->GetXaxis()->SetTitle("Energy [GeV]");
279 fOutputList->Add(fHistClusterEN);
280
281 fHistClusterEtN = new TH2F("fHistClusterEtN","N vs Cluster Transverse Energy of Pb+Pb",100,0.0,20.0,30,0.0,30.0);
282 fHistClusterEtN->GetYaxis()->SetTitle("N"); fHistClusterEtN->GetXaxis()->SetTitle("Transverse Energy [GeV]");
283 fOutputList->Add(fHistClusterEtN);
284
285 fHistClusterdphiV0 = new TH1D("fHistClusterdphiV0","Cluster dphiV0 Distribution of Pb+Pb",100,0.0,TMath::Pi());
286 fHistClusterdphiV0->GetYaxis()->SetTitle("Entries"); fHistClusterdphiV0->GetXaxis()->SetTitle("dphiV0 [rad]");
287 fOutputList->Add(fHistClusterdphiV0);
288
289 fHistEPV0 = new TH2F("fHistEPV0","V0 Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
290 fHistEPV0->GetYaxis()->SetTitle("V0 Event Plane Angle [rad]"); fHistEPV0->GetXaxis()->SetTitle("V0M Centrality");
291 fOutputList->Add(fHistEPV0);
292
293 fHistEPV0A = new TH2F("fHistEPV0A","V0A Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
294 fHistEPV0A->GetYaxis()->SetTitle("V0A Event Plane Angle [rad]"); fHistEPV0A->GetXaxis()->SetTitle("V0M Centrality");
295 fOutputList->Add(fHistEPV0A);
296
297 fHistEPV0C = new TH2F("fHistEPV0","V0C Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
298 fHistEPV0C->GetYaxis()->SetTitle("V0C Event Plane Angle [rad]"); fHistEPV0C->GetXaxis()->SetTitle("V0M Centrality");
299 fOutputList->Add(fHistEPV0C);
300
301 fHistEPV0r = new TH2F("fHistEPV0r","V0r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
302 fHistEPV0r->GetYaxis()->SetTitle("V0r Event Plane Angle [rad]"); fHistEPV0r->GetXaxis()->SetTitle("V0M Centrality");
303 fOutputList->Add(fHistEPV0r);
304
305 fHistEPV0Ar = new TH2F("fHistEPV0","V0Ar Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
306 fHistEPV0Ar->GetYaxis()->SetTitle("V0Ar Event Plane Angle [rad]"); fHistEPV0Ar->GetXaxis()->SetTitle("V0M Centrality");
307 fOutputList->Add(fHistEPV0Ar);
308
309 fHistEPV0Cr = new TH2F("fHistEPV0Cr","V0Cr Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
310 fHistEPV0Cr->GetYaxis()->SetTitle("V0Cr Event Plane Angle [rad]"); fHistEPV0Cr->GetXaxis()->SetTitle("V0M Centrality");
311 fOutputList->Add(fHistEPV0Cr);
312
313 fHistEPV0A4r = new TH2F("fHistEPV0Cr","V0A4r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
314 fHistEPV0A4r->GetYaxis()->SetTitle("V0A4r Event Plane Angle [rad]"); fHistEPV0A4r->GetXaxis()->SetTitle("V0M Centrality");
315 fOutputList->Add(fHistEPV0A4r);
316
317 fHistEPV0A7r = new TH2F("fHistEPV0Cr","V0A7r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
318 fHistEPV0A7r->GetYaxis()->SetTitle("V0A7r Event Plane Angle [rad]"); fHistEPV0A7r->GetXaxis()->SetTitle("V0M Centrality");
319 fOutputList->Add(fHistEPV0A7r);
320
321 fHistEPV0C0r = new TH2F("fHistEPV0Cr","V0C0r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
322 fHistEPV0C0r->GetYaxis()->SetTitle("V0C0r Event Plane Angle [rad]"); fHistEPV0C0r->GetXaxis()->SetTitle("V0M Centrality");
323 fOutputList->Add(fHistEPV0C0r);
324
325 fHistEPV0C3r = new TH2F("fHistEPV0Cr","V0C3r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
326 fHistEPV0C3r->GetYaxis()->SetTitle("V0C3r Event Plane Angle [rad]"); fHistEPV0C3r->GetXaxis()->SetTitle("V0M Centrality");
327 fOutputList->Add(fHistEPV0C3r);
328
329 fHistdifV0A_V0C0r = new TH2F("fHistdifV0A_V0C0r","(V0A - V0C0r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
330 fHistdifV0A_V0C0r->GetYaxis()->SetTitle("Cos[2*(V0A - V0C0r)]"); fHistdifV0A_V0C0r->GetXaxis()->SetTitle("V0M Centrality");
331 fOutputList->Add(fHistdifV0A_V0C0r);
332
333 fHistdifV0A_V0C3r = new TH2F("fHistdifV0A_V0C3r","(V0A - V0C3r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
334 fHistdifV0A_V0C3r->GetYaxis()->SetTitle("Cos[2*(V0A - V0C3r)]"); fHistdifV0A_V0C3r->GetXaxis()->SetTitle("V0M Centrality");
335 fOutputList->Add(fHistdifV0A_V0C3r);
336
337 fHistdifV0C0r_V0C3r = new TH2F("fHistdifV0C0r_V0C3r","(V0C0r - V0C3r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
338 fHistdifV0C0r_V0C3r->GetYaxis()->SetTitle("Cos[2*(V0C0r - V0C3r)]"); fHistdifV0C0r_V0C3r->GetXaxis()->SetTitle("V0M Centrality");
339 fOutputList->Add(fHistdifV0C0r_V0C3r);
340
341 fHistdifV0C_V0A4r = new TH2F("fHistdifV0C_V0A4r","(V0C - V0A4r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
342 fHistdifV0C_V0A4r->GetYaxis()->SetTitle("Cos[2*(V0C - V0C4r)]"); fHistdifV0C_V0A4r->GetXaxis()->SetTitle("V0M Centrality");
343 fOutputList->Add(fHistdifV0C_V0A4r);
344
345 fHistdifV0C_V0A7r = new TH2F("fHistdifV0C_V0A3r","(V0C - V0A7r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
346 fHistdifV0C_V0A7r->GetYaxis()->SetTitle("Cos[2*(V0C - V0A7r)]"); fHistdifV0C_V0A7r->GetXaxis()->SetTitle("V0M Centrality");
347 fOutputList->Add(fHistdifV0C_V0A7r);
348
349 fHistdifV0A4r_V0A7r = new TH2F("fHistdifV0A4r_V0A7r","(V0A4r - V0A7r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
350 fHistdifV0A4r_V0A7r->GetYaxis()->SetTitle("Cos[2*(V0A4r - V0A7r)]"); fHistdifV0A4r_V0A7r->GetXaxis()->SetTitle("V0M Centrality");
351 fOutputList->Add(fHistdifV0A4r_V0A7r);
352
353 fHistdifV0Ar_V0Cr = new TH2F("fHistdifV0Ar_V0Cr","(V0Ar - V0Cr) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
354 fHistdifV0Ar_V0Cr->GetYaxis()->SetTitle("Cos[2*(V0Ar - V0Cr)]"); fHistdifV0Ar_V0Cr->GetXaxis()->SetTitle("V0M Centrality");
355 fOutputList->Add(fHistdifV0Ar_V0Cr);
356
357 fHistAllcentV0 = new TH1F("fHistAllcentV0","V0 Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
358 fHistAllcentV0->GetXaxis()->SetTitle("V0 Event Plane Angle [rad]"); fHistAllcentV0->GetYaxis()->SetTitle("");
359 fOutputList->Add(fHistAllcentV0);
360
361 fHistAllcentV0r = new TH1F("fHistAllcentV0r","V0r Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
362 fHistAllcentV0r->GetXaxis()->SetTitle("V0r Event Plane Angle [rad]"); fHistAllcentV0r->GetYaxis()->SetTitle("");
363 fOutputList->Add(fHistAllcentV0r);
364
365 fHistAllcentV0A = new TH1F("fHistAllcentV0A","V0A Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
366 fHistAllcentV0A->GetXaxis()->SetTitle("V0A Event Plane Angle [rad]"); fHistAllcentV0A->GetYaxis()->SetTitle("");
367 fOutputList->Add(fHistAllcentV0A);
368
369 fHistAllcentV0C = new TH1F("fHistAllcentV0C","V0C Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
370 fHistAllcentV0C->GetXaxis()->SetTitle("V0C Event Plane Angle [rad]"); fHistAllcentV0C->GetYaxis()->SetTitle("");
371 fOutputList->Add(fHistAllcentV0C);
372
373 fHistAllcentTPC = new TH1F("fHistAllcentTPC","TPC Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
374 fHistAllcentTPC->GetXaxis()->SetTitle("TPC Event Plane Angle [rad]"); fHistAllcentTPC->GetYaxis()->SetTitle("");
375 fOutputList->Add(fHistAllcentTPC);
376
377 fHistEPTPC = new TH2F("fHistEPTPC","TPC Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
378 fHistEPTPC->GetYaxis()->SetTitle("Event Plane Angle [rad]"); fHistEPTPC->GetXaxis()->SetTitle("V0M Centrality");
379 fOutputList->Add(fHistEPTPC);
380
381 fHistEPTPCResolution = new TH2F("fHistEPTPCResolution","TPC Resolution vs Centrality of Pb+Pb",100,0,100,100,0.0,1.0);
382 fHistEPTPCResolution->GetYaxis()->SetTitle("TPC Resolution"); fHistEPTPCResolution->GetXaxis()->SetTitle("V0M Centrality");
383 fOutputList->Add(fHistEPTPCResolution);
384
385
386 // Et M02 V0Mcent DeltaPhi Cos[2*DeltaPhi]
387 Int_t bins[5] = { 500, 250, 100, 100, 100 }; // binning
388 Double_t min[5] = { 0.0, 0.0, 0, 0.0, -1.0}; // min x
389 Double_t max[5] = { 50.0, 2.5, 100, TMath::Pi(), 1.0}; // max x
390
391 fClusterPbV0 = new THnSparseF("fClusterPbV0","",5,bins,min,max);
392 fClusterPbV0->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0->GetAxis(1)->SetTitle("M02"); fClusterPbV0->GetAxis(2)->SetTitle("V0M Centrality");
393 fClusterPbV0->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
394 fOutputList->Add(fClusterPbV0);
395
396 fClusterPbV0A = new THnSparseF("fClusterPbV0A","",5,bins,min,max);
397 fClusterPbV0A->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0A->GetAxis(1)->SetTitle("M02"); fClusterPbV0A->GetAxis(2)->SetTitle("V0M Centrality");
398 fClusterPbV0A->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0A->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
399 fOutputList->Add(fClusterPbV0A);
400
401 fClusterPbV0C = new THnSparseF("fClusterPbV0C","",5,bins,min,max);
402 fClusterPbV0C->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0C->GetAxis(1)->SetTitle("M02"); fClusterPbV0C->GetAxis(2)->SetTitle("V0M Centrality");
403 fClusterPbV0C->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0C->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
404 fOutputList->Add(fClusterPbV0C);
405
406 fClusterPbTPC = new THnSparseF("fClusterPbTPC","",5,bins,min,max);
407 fClusterPbTPC->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbTPC->GetAxis(1)->SetTitle("M02"); fClusterPbTPC->GetAxis(2)->SetTitle("V0M Centrality");
408 fClusterPbTPC->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbTPC->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
409 fOutputList->Add(fClusterPbTPC);
410
411
412 PostData(1, fOutputList);
413}
414
415//________________________________________________________________________
416void AliAnalysisTaskEMCALPi0V2ShSh::UserExec(Option_t *)
417{
418 // Main loop, called for each event.
c4a4ee27 419 fESDClusters = 0;
420 fESDCells = 0;
421 fAODClusters = 0;
422 fAODCells = 0;
fc0f5568 423
424 // Create pointer to reconstructed event
425 AliVEvent *event = InputEvent();
426 if (!event) { Printf("ERROR: Could not retrieve event\n"); return; }
427
428 fESD = dynamic_cast<AliESDEvent*>(event);
429 if (!fESD) {
c4a4ee27 430 fAOD = dynamic_cast<AliAODEvent*>(event);
431 if(!fAOD){
432 printf("ERROR: Could not retrieve the event\n");
433 return;
434 }
fc0f5568 435 }
436
437 //Bool_t isSelected =0;
438 //isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral));
439 //if(!isSelected) { return; }
440
c4a4ee27 441 if(event->GetCentrality()){
442 fCentralityV0M = event->GetCentrality()->GetCentralityPercentile("V0M");
fc0f5568 443 } else {
444 printf("ERROR: Could not retrieve the Centrality\n");
445 return;
446 }
c4a4ee27 447
448 fEventPlane = event->GetEventplane();
fc0f5568 449 if (fEventPlane) {
450 VZEROEventPlane();
451 } else {
452 printf("ERROR: Could not retrieve the Centrality\n");
453 return;
454 }
c4a4ee27 455 Int_t runnumber = InputEvent()->GetRunNumber() ;
456
457 TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
458 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
459 if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
460 break;
461 /*if(fESD)
462 fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
463 else*/
464 // if(fVEvent->GetEMCALMatrix(mod))
465 fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
466 fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
fc0f5568 467 }
468
c4a4ee27 469 if(fESD){
470 TList *l = fESD->GetList();
471 if(1){//fDebug){
472 for(int nk=0;nk<l->GetEntries();nk++){
473 TObject *obj = (TObject*)l->At(nk);
474 TString oname = obj->GetName();
475 //if(oname.Contains("lus"))
476 printf("Object %d has a clus array named %s +++++++++\n",nk,oname.Data());
477 }
478 }
479 fESDClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
480 fESDCells = fESD->GetEMCALCells();
481 if(fDebug)
482 printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
483 if(fESDClusters->GetEntriesFast()<1){
484 printf("ERROR: There are no EMCAL clusters in this event\n");
485 return;
486 }
487 }
488 else if(fAOD){
489 fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
490 fAODCells = fAOD->GetEMCALCells();
491 if(fDebug)
492 printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
493 if(fAODClusters->GetEntriesFast()<1){
494 printf("ERROR: There are no EMCAL clusters in this event\n");
495 return;
496 }
497 }
498
499 if(fESD)
500 FillTrackHists();
fc0f5568 501
fc0f5568 502 PostData(1, fOutputList);
503}
504
505//________________________________________________________________________
506void AliAnalysisTaskEMCALPi0V2ShSh::VZEROEventPlane()
507{// Calculate the V0 Event Plane
508
509 if (fEventPlane->GetQVector()) {
510 fEPTPC = TVector2::Phi_0_2pi(fEventPlane->GetQVector()->Phi())/2.0; //if(fEPTPC>TMath::Pi()) {fEPTPC-=TMath::Pi();}
511 } else { fEPTPC = -999.; }
512
513 if (fEventPlane->GetQsub1()&&fEventPlane->GetQsub2()) {
514 fEPTPCResolution = TMath::Cos(2.0*(fEventPlane->GetQsub1()->Phi()/2.0-fEventPlane->GetQsub2()->Phi()/2.0)); }
515 else { fEPTPCResolution = -1; }
516
517 fEPV0 = TVector2::Phi_0_2pi(fEventPlane->GetEventplane("V0", fESD)); if(fEPV0>TMath::Pi()) {fEPV0-=TMath::Pi();}
518 fEPV0A = TVector2::Phi_0_2pi(fEventPlane->GetEventplane("V0A", fESD)); if(fEPV0A>TMath::Pi()) {fEPV0A-=TMath::Pi();}
519 fEPV0C = TVector2::Phi_0_2pi(fEventPlane->GetEventplane("V0C", fESD)); if(fEPV0C>TMath::Pi()) {fEPV0C-=TMath::Pi();}
520
521 Double_t qx=0, qy=0, qxr=0, qyr=0;
522 fEPV0Ar = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 4, 5, 2, qxr, qyr)); if(fEPV0Ar>TMath::Pi()) {fEPV0Ar-=TMath::Pi();}
523 fEPV0Cr = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 2, 3, 2, qx, qy)); if(fEPV0Cr>TMath::Pi()) {fEPV0Cr-=TMath::Pi();}
524 qxr += qx; qyr += qy;
525 fEPV0r = TVector2::Phi_0_2pi(TMath::ATan2(qyr,qxr))/2.0;
526
527 fEPV0A4r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 4, 2, qx, qy)); if(fEPV0A4r>TMath::Pi()) {fEPV0A4r-=TMath::Pi();}
528 fEPV0A5r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 5, 2, qx, qy)); if(fEPV0A5r>TMath::Pi()) {fEPV0A5r-=TMath::Pi();}
529 fEPV0A6r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 6, 2, qx, qy)); if(fEPV0A6r>TMath::Pi()) {fEPV0A6r-=TMath::Pi();}
530 fEPV0A7r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 7, 2, qx, qy)); if(fEPV0A7r>TMath::Pi()) {fEPV0A7r-=TMath::Pi();}
531 fEPV0C0r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 0, 2, qx, qy)); if(fEPV0C0r>TMath::Pi()) {fEPV0C0r-=TMath::Pi();}
532 fEPV0C1r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 1, 2, qx, qy)); if(fEPV0C1r>TMath::Pi()) {fEPV0C1r-=TMath::Pi();}
533 fEPV0C2r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 2, 2, qx, qy)); if(fEPV0C2r>TMath::Pi()) {fEPV0C2r-=TMath::Pi();}
534 fEPV0C3r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 3, 2, qx, qy)); if(fEPV0C3r>TMath::Pi()) {fEPV0C3r-=TMath::Pi();}
535
536 fHistEPTPC->Fill(fCentralityV0M, fEPTPC);
537 if(fEPTPCResolution!=-1) { fHistEPTPCResolution->Fill(fCentralityV0M, fEPTPCResolution); }
538 fHistEPV0->Fill(fCentralityV0M, fEPV0);
539 fHistEPV0A->Fill(fCentralityV0M, fEPV0A);
540 fHistEPV0C->Fill(fCentralityV0M, fEPV0C);
541 fHistEPV0Ar->Fill(fCentralityV0M, fEPV0Ar);
542 fHistEPV0Cr->Fill(fCentralityV0M, fEPV0Cr);
543 fHistEPV0r->Fill(fCentralityV0M, fEPV0r);
544 fHistEPV0A4r->Fill(fCentralityV0M, fEPV0A4r);
545 fHistEPV0A7r->Fill(fCentralityV0M, fEPV0A7r);
546 fHistEPV0C0r->Fill(fCentralityV0M, fEPV0C0r);
547 fHistEPV0C3r->Fill(fCentralityV0M, fEPV0C3r);
548
549 fHistAllcentV0->Fill(fEPV0);
550 fHistAllcentV0r->Fill(fEPV0r);
551 fHistAllcentV0A->Fill(fEPV0A);
552 fHistAllcentV0C->Fill(fEPV0C);
553 fHistAllcentTPC->Fill(fEPTPC);
554
555 fHistdifV0A_V0C0r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0A - fEPV0C0r)));
556 fHistdifV0A_V0C3r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0A - fEPV0C3r)));
557 fHistdifV0C0r_V0C3r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0C0r - fEPV0C3r)));
558 fHistdifV0C_V0A4r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0C - fEPV0A4r)));
559 fHistdifV0C_V0A7r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0C - fEPV0A7r)));
560 fHistdifV0A4r_V0A7r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0A4r - fEPV0A7r)));
561 fHistdifV0Ar_V0Cr->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0Ar - fEPV0Cr)));
562
563}
564
565//________________________________________________________________________
566void AliAnalysisTaskEMCALPi0V2ShSh::FillClusterHists()
567{// Fill cluster histograms.
568
569 Float_t pos[3] ;
570
c4a4ee27 571 TObjArray *clusters = fESDClusters;
572 if (!clusters)
573 clusters = fAODClusters;
574 if (!clusters)
fc0f5568 575 return;
576
c4a4ee27 577 const Int_t nclus = clusters->GetEntries();
fc0f5568 578
579 if(nclus==0)
580 return;
581
582 for(Int_t iclus=0;iclus<nclus;iclus++){
c4a4ee27 583 AliVCluster *clus = (AliVCluster *) clusters->At(iclus); // retrieve cluster from esd
fc0f5568 584 if(!clus)
585 continue;
586 if(!clus->IsEMCAL()){
587 printf("ERROR: Current Cluster is not an EMCAL Cluster\n");
588 continue;
589 }
590 clus->GetPosition(pos);
591 TVector3 vpos(pos[0],pos[1],pos[2]);
592 Double_t Transverse_Energy = ((clus->E())/ (TMath::CosH(vpos.Eta())));
593 Double_t dphiV0 = TVector2::Phi_0_2pi(vpos.Phi()-fEPV0r); if(dphiV0>TMath::Pi()) {dphiV0-=TMath::Pi();}
594 Double_t dphiV0A = TVector2::Phi_0_2pi(vpos.Phi()-fEPV0A); if(dphiV0A>TMath::Pi()) {dphiV0A-=TMath::Pi();}
595 Double_t dphiV0C = TVector2::Phi_0_2pi(vpos.Phi()-fEPV0C); if(dphiV0C>TMath::Pi()) {dphiV0C-=TMath::Pi();}
596 Double_t dphiTPC = TVector2::Phi_0_2pi(vpos.Phi()-fEPTPC); if(dphiTPC>TMath::Pi()) {dphiTPC-=TMath::Pi();}
597
598 Double_t DataV0[5];
599 DataV0[0] = Transverse_Energy;
600 DataV0[1] = clus->GetM02();
601 DataV0[2] = fCentralityV0M;
602 DataV0[3] = dphiV0;
603 DataV0[4] = TMath::Cos(2.0*(dphiV0));
604 fClusterPbV0->Fill(DataV0);
605
606 Double_t DataV0A[5];
607 DataV0A[0] = Transverse_Energy;
608 DataV0A[1] = clus->GetM02();
609 DataV0A[2] = fCentralityV0M;
610 DataV0A[3] = dphiV0A;
611 DataV0A[4] = TMath::Cos(2.0*(dphiV0A));
612 fClusterPbV0A->Fill(DataV0A);
613
614 Double_t DataV0C[5];
615 DataV0C[0] = Transverse_Energy;
616 DataV0C[1] = clus->GetM02();
617 DataV0C[2] = fCentralityV0M;
618 DataV0C[3] = dphiV0C;
619 DataV0C[4] = TMath::Cos(2.0*(dphiV0C));
620 fClusterPbV0C->Fill(DataV0C);
621
622 Double_t DataTPC[5];
623 DataTPC[0] = Transverse_Energy;
624 DataTPC[1] = clus->GetM02();
625 DataTPC[2] = fCentralityV0M;
626 DataTPC[3] = dphiTPC;
627 DataTPC[4] = TMath::Cos(2.0*(dphiTPC));
628 fClusterPbTPC->Fill(DataTPC);
629
630 fHistClusterE->Fill(clus->E());
631 fHistClusterEt->Fill(Transverse_Energy);
632 fHistClusterM02->Fill(clus->GetM02());
633 fHistClusterN->Fill(clus->GetNCells());
634 fHistClusterPhi->Fill(vpos.Phi());
635 fHistClusterEta->Fill(vpos.Eta());
636 fHistClusterPhiEta->Fill(vpos.Eta(),vpos.Phi());
637 fHistClusterEN->Fill(clus->E(),clus->GetNCells());
638 fHistClusterEM02->Fill(clus->E(),clus->GetM02());
639 fHistClusterEtN->Fill(Transverse_Energy,clus->GetNCells());
640 fHistClusterEtM02->Fill(Transverse_Energy,clus->GetM02());
641 fHistClusterdphiV0->Fill(dphiV0);
642
643 }
644}
645
646//________________________________________________________________________
647void AliAnalysisTaskEMCALPi0V2ShSh::FillTrackHists()
648{// Fill track histograms.
c4a4ee27 649 if(!fESD)
650 return;
fc0f5568 651 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
652 AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
653 if (!track){
654 printf("ERROR: Could not retreive esdtrack\n");
655 continue;
656 }
657 fHistTrackPt->Fill(track->Pt());
658 fHistTrackEta->Fill(track->Eta());
659 fHistTrackPhi->Fill(track->Phi());
660 fHistTrackPhiEta->Fill(track->Eta(),track->Phi());
661
662 }
663}
664
665//________________________________________________________________________
666void AliAnalysisTaskEMCALPi0V2ShSh::Terminate(Option_t *)
667{
668 // Draw result to screen, or perform fitting, normalizations
669 // Called once at the end of the query
670
671 /* fOutputList = dynamic_cast<TList*> (GetOutputData(1));
672 if(!fOutputList) { Printf("ERROR: could not retrieve TList fOutputList"); return; }
673
674 fHistTrackPt = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistTrackPt"));
675 if (!fHistTrackPt) { Printf("ERROR: could not retrieve fHistTrackPt"); return;}
676 fHistTrackEta = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistTrackEta"));
677 if (!fHistTrackEta) { Printf("ERROR: could not retrieve fHistTrackEta"); return;}
678 fHistTrackPhi = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistTrackPhi"));
679 if (!fHistTrackPhi) { Printf("ERROR: could not retrieve fHistTrackPhi"); return;}
680 fHistClusterEta = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterEta"));
681 if (!fHistClusterEta) { Printf("ERROR: could not retrieve fHistClusterEta"); return;}
682 fHistClusterPhi = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterPhi"));
683 if (!fHistClusterPhi) { Printf("ERROR: could not retrieve fHistClusterPhi"); return;}
684 fHistClusterE = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterE"));
685 if (!fHistClusterE) { Printf("ERROR: could not retrieve fHistClusterE"); return;}
686 fHistClusterN = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterN"));
687 if (!fHistClusterN) { Printf("ERROR: could not retrieve fHistClusterN"); return;}
688 fHistClusterM02 = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterM02"));
689 if (!fHistClusterM02) { Printf("ERROR: could not retrieve fHistClusterM02"); return;}
690
691 // Get the physics selection histograms with the selection statistics
692 //AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
693 //AliESDInputHandler *inputH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
694 //TH2F *histStat = (TH2F*)inputH->GetStatistics();
695
696
697 TCanvas *canvas1 = new TCanvas("canvas1","Track P_{T} & #eta & #phi",10,10,1020,510);
698 canvas1->Divide(3,1);
699 canvas1->cd(1)->SetLogy();
700 fHistTrackPt->DrawCopy("");
701 canvas1->cd(2);
702 fHistTrackEta->DrawCopy("");
703 canvas1->cd(3);
704 fHistTrackPhi->DrawCopy("");
705
706 TCanvas *canvas2 = new TCanvas("canvas2","Cluster #eta & #phi",10,10,1020,510);
707 canvas2->Divide(2,1);
708 canvas2->cd(1);
709 fHistClusterEta->DrawCopy("");
710 canvas2->cd(2);
711 fHistClusterPhi->DrawCopy("");
712
713 TCanvas *canvas3 = new TCanvas("canvas3","Cluster E & N & M02",10,10,1020,510);
714 canvas3->Divide(3,1);
715 canvas3->cd(1)->SetLogy();
716 fHistClusterE->DrawCopy("");
717 canvas3->cd(2)->SetLogy();
718 fHistClusterN->DrawCopy("");
719 canvas3->cd(3)->SetLogy();
720 fHistClusterM02->DrawCopy("");
721
722 TCanvas *canvas4 = new TCanvas("canvas4","Track #phi vs #eta & Cluster #phi vs #eta",10,10,1020,510);
723 canvas4->Divide(2,1);
724 canvas4->cd(1)->SetLogz();
725 fHistTrackPhiEta->DrawCopy("colorz");
726 canvas4->cd(2)->SetLogz();
727 fHistClusterPhiEta->DrawCopy("colorz");
728
729 TCanvas *canvas5 = new TCanvas("canvas5","Cluster E vs N & E vs M02",10,10,1020,510);
730 canvas5->Divide(2,1);
731 canvas5->cd(1)->SetLogz();
732 fHistClusterEN->DrawCopy("colorz");
733 canvas5->cd(2)->SetLogz();
734 fHistClusterEM02->DrawCopy("colorz");
735
736 TCanvas *canvas6 = new TCanvas("canvas6","Cluster Et & Et vs N & Et vs M02",10,10,1020,510);
737 canvas6->cd(1)->SetLogy();
738 fHistClusterEt->DrawCopy("");
739 canvas6->cd(2)->SetLogz();
740 fHistClusterEtN->DrawCopy("colorz");
741 canvas6->cd(3)->SetLogz();
742 fHistClusterEtM02->DrawCopy("colorz");
743
744 canvas1->SaveAs("lhc11h_2_Track_PT_Eta_Phi.jpg");
745 canvas2->SaveAs("lhc11h_2_Cluster_Eta_Phi.jpg");
746 canvas3->SaveAs("lhc11h_2_Cluster_E_N_M02.jpg");
747 canvas4->SaveAs("lhc11h_2_Track_PhivsEta_Cluster_PhivsEta.jpg");
748 canvas5->SaveAs("lhc11h_2_Cluster_EvsN_EvsM02.jpg");
749 canvas6->SaveAs("lhc11h_2_Cluster_Et_EtvsN_EtvsM02.jpg"); */
750}