]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0V2ShSh.cxx
Pi0 peak band adjusted to 2013 period
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALPi0V2ShSh.cxx
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"
26 #include "AliEMCALGeometry.h"
27 #include "AliESDEvent.h"
28 #include "AliESDVertex.h"
29 #include "AliESDCaloCells.h"
30 #include "AliESDCaloCluster.h"
31 #include "AliESDEvent.h"
32 #include "AliESDHeader.h"
33 #include "AliESDInputHandler.h"
34 #include "AliESDtrack.h"
35 #include "AliKFParticle.h"
36 #include "AliAODEvent.h"
37 #include "AliVCluster.h"
38 #include "AliCentrality.h"
39 #include "AliEventplane.h"
40 #include "AliOADBContainer.h"
41
42
43
44 ClassImp(AliAnalysisTaskEMCALPi0V2ShSh)
45
46 //________________________________________________________________________
47 AliAnalysisTaskEMCALPi0V2ShSh::AliAnalysisTaskEMCALPi0V2ShSh() : 
48   AliAnalysisTaskSE(), 
49   fEventPlane(0),
50   fCentralityV0M(99.),
51   fESDClusters(0),
52   fAODClusters(0),
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)
121 {
122   // Default constructor.
123   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
124 }
125
126 //________________________________________________________________________
127 AliAnalysisTaskEMCALPi0V2ShSh::AliAnalysisTaskEMCALPi0V2ShSh(const char *name) : 
128   AliAnalysisTaskSE(name), 
129   fEventPlane(0), 
130   fCentralityV0M(99.), 
131   fESDClusters(0),
132   fAODClusters(0),
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)
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 //________________________________________________________________________
212 void AliAnalysisTaskEMCALPi0V2ShSh::UserCreateOutputObjects()
213 {
214   // Create histograms, called once.
215     
216   fESDClusters = new TObjArray();
217   fAODClusters = new TObjArray();
218   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
219   fOADBContainer = new AliOADBContainer("AliEMCALgeo");
220   fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
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 //________________________________________________________________________
416 void AliAnalysisTaskEMCALPi0V2ShSh::UserExec(Option_t *) 
417 {
418   // Main loop, called for each event.
419   fESDClusters = 0;
420   fESDCells = 0;
421   fAODClusters = 0;
422   fAODCells = 0;
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) {
430     fAOD =  dynamic_cast<AliAODEvent*>(event);
431     if(!fAOD){
432       printf("ERROR: Could not retrieve the event\n");
433       return;
434     }
435   }
436
437   //Bool_t isSelected =0;      
438   //isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral));
439   //if(!isSelected) { return; }
440   
441   if(event->GetCentrality()){
442     fCentralityV0M = event->GetCentrality()->GetCentralityPercentile("V0M");
443   } else {
444     printf("ERROR: Could not retrieve the Centrality\n");
445     return;
446   }     
447   
448   fEventPlane = event->GetEventplane(); 
449   if (fEventPlane) {
450      VZEROEventPlane();
451   } else {
452      printf("ERROR: Could not retrieve the Centrality\n");
453      return;
454   }
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);
467   }
468
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();
501
502   PostData(1, fOutputList);
503 }
504
505 //________________________________________________________________________
506 void 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 //________________________________________________________________________
566 void AliAnalysisTaskEMCALPi0V2ShSh::FillClusterHists()
567 {// Fill cluster histograms.
568
569   Float_t pos[3] ; 
570
571   TObjArray *clusters = fESDClusters;
572   if (!clusters)
573     clusters = fAODClusters;
574   if (!clusters)
575     return;
576
577   const Int_t nclus = clusters->GetEntries();
578
579   if(nclus==0)
580     return;
581
582   for(Int_t iclus=0;iclus<nclus;iclus++){
583     AliVCluster *clus = (AliVCluster *) clusters->At(iclus); // retrieve cluster from esd
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 //________________________________________________________________________
647 void AliAnalysisTaskEMCALPi0V2ShSh::FillTrackHists()
648 {// Fill track histograms.
649   if(!fESD)
650     return;
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 //________________________________________________________________________
666 void 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 }