]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskEMCALMesonGGSDMpPb.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing left-over "return".
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALMesonGGSDMpPb.cxx
CommitLineData
3a66195e 1#include "AliAnalysisTaskEMCALMesonGGSDMpPb.h"
2
3// ROOT includes
4#include <vector>
5#include <Riostream.h>
6#include <TChain.h>
7#include <TTree.h>
8#include <TF1.h>
9#include <TH1F.h>
10#include <TH2F.h>
11#include <TH3F.h>
12#include <TH1D.h>
13#include <TH2D.h>
14#include <TH3D.h>
15#include <TCanvas.h>
16#include <TList.h>
17#include <TFile.h>
18#include <TLorentzVector.h>
19#include <TNtuple.h>
20#include <TRandom3.h>
21#include <TGeoManager.h>
22#include <TGeoMatrix.h>
23#include <TGeoBBox.h>
24#include <TArrayI.h>
25#include <TArrayF.h>
26#include <TObjArray.h>
27
28// STEER? includes
29#include "AliAnalysisTaskSE.h"
30#include "AliAnalysisManager.h"
31#include "AliVCluster.h"
32#include "AliVCaloCells.h"
33#include "AliLog.h"
34#include "AliPID.h"
35#include "AliStack.h"
36#include "AliESDtrack.h"
37#include "AliESDtrackCuts.h"
38#include "AliESDEvent.h"
39#include "AliAODEvent.h"
40#include "AliMCEvent.h"
41#include "AliInputEventHandler.h"
42#include "AliESDInputHandler.h"
43#include "AliAODInputHandler.h"
44#include "AliAODTrack.h"
45#include "AliExternalTrackParam.h"
46#include "AliESDfriendTrack.h"
47#include "AliTrackerBase.h"
48
49// EMCAL includes
50#include "AliEMCALRecoUtils.h"
51#include "AliEMCALGeometry.h"
52#include "AliTrackerBase.h"
53#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
54#include "AliEMCALPIDUtils.h"
55#include "AliExternalTrackParam.h"
56
57#include "AliCentrality.h"
58
2942f542 59using std::cout;
60using std::endl;
3a66195e 61
62
63ClassImp(AliAnalysisTaskEMCALMesonGGSDMpPb)
64
65//________________________________________________________________________
66AliAnalysisTaskEMCALMesonGGSDMpPb::AliAnalysisTaskEMCALMesonGGSDMpPb() :
67 AliAnalysisTaskSE(),
68 fOutput(0),
69 fMcMode(0),
70 fRecalibrator(0),
71 fdRmin_ClustTrack(0),
72 fPhimin(0),
73 fPhimax(0),
74 fEtamin(0),
75 fEtamax(0),
76 fTrackCuts(0),
77 fEsdEv(0),
78 fAodEv(0),
79 h1_zvtx(0),
80 h1_trigger(0),
81 h1_centrality(0),
82 h2_PhiEtaCluster(0),
83 h2_PhiEtaClusterCut(0),
84 h2_PhiEtaMaxCell(0),
85 h2_PhiEtaMaxCellCut(0),
86 h2_gE_RecTruth(0),
87 h2_eop_E(0),
88 h2_eop_pT(0),
89 h2_E_time(0),
90 h2_Pi0TruthPhiEta(0),
91 h2_PriPi0TruthPhiEta(0),
92 h2_Pi0TruthPhiEtaEmcal(0),
93 h2_PriPi0TruthPhiEtaEmcal(0),
94 h2_Pi0TruthPhiEta_Phi2piEta065(0),
95 h2_Pi0TruthPhiEta_Phi2piEta1(0),
96 h2_TruthPhotonsPhiEta(0),
97 h2_PhotonsPhiEtaIsEmcal(0),
98 TriggerList(0),
99 fHelperClass(0)
100{
101 // Dummy constructor ALWAYS needed for I/O.
102 for(int i=0; i<cent_bins; i++){
103 h1_nClusters[i] = 0;
104 h1_M[i] = 0;
105 h1_M_mix[i] = 0;
106 h1_E[i] = 0;
107 h1_dR_ClustTrk[i] = 0;
108 h1_Pi0TruthPt[i] = 0;
109 h1_PriPi0TruthPt[i] = 0;
110 h1_Pi0TruthPtEmcal[i] = 0;
111 h1_PriPi0TruthPtEmcal[i] = 0;
112 h1_Pi0TruthPtPhi2piEta065[i] = 0;
113 h1_Pi0TruthPtPhi2piEta1[i] = 0;
114 h1_TruthPhotonsEmcal[i] = 0;
115 h1_PhotonsEmcal[i] = 0;
116 h1_PhotonsNCellsCut[i] = 0;
117 h1_PhotonsTrackMatchCut[i] = 0;
118 h1_PhotonsAllCut[i] = 0;
119 h1_dR_RealMC[i] = 0;
120 h1_Chi2[i] = 0;
121 h1_nTrkMatch[i] = 0;
122 h1_nCells[i] = 0;
123 h1_ClusterDisp[i] = 0;
124 h2_Ellipse[i] = 0;
125 h2_EtaPt[i] = 0;
126 h3_MptAsymm[i] = 0;
127 h3_MptAsymm_mix[i] = 0;
128 h2_dphi_deta[i] = 0;
129 h2_dphi_deta_mix[i] = 0;
130 h2_DispRes[i] = 0;
131 h2_cells_M02[i] = 0;
132 }
133}
134
135//________________________________________________________________________
136AliAnalysisTaskEMCALMesonGGSDMpPb::AliAnalysisTaskEMCALMesonGGSDMpPb(const char *name) :
137 AliAnalysisTaskSE(name),
138 fOutput(0),
139 fMcMode(0),
140 fRecalibrator(0),
141 fdRmin_ClustTrack(0),
142 fPhimin(0),
143 fPhimax(0),
144 fEtamin(0),
145 fEtamax(0),
146 fTrackCuts(0),
147 fEsdEv(0),
148 fAodEv(0),
149 h1_zvtx(0),
150 h1_trigger(0),
151 h1_centrality(0),
152 h2_PhiEtaCluster(0),
153 h2_PhiEtaClusterCut(0),
154 h2_PhiEtaMaxCell(0),
155 h2_PhiEtaMaxCellCut(0),
156 h2_gE_RecTruth(0),
157 h2_eop_E(0),
158 h2_eop_pT(0),
159 h2_E_time(0),
160 h2_Pi0TruthPhiEta(0),
161 h2_PriPi0TruthPhiEta(0),
162 h2_Pi0TruthPhiEtaEmcal(0),
163 h2_PriPi0TruthPhiEtaEmcal(0),
164 h2_Pi0TruthPhiEta_Phi2piEta065(0),
165 h2_Pi0TruthPhiEta_Phi2piEta1(0),
166 h2_TruthPhotonsPhiEta(0),
167 h2_PhotonsPhiEtaIsEmcal(0),
168 TriggerList(0),
169 fHelperClass(0)
170{
171 // Constructor
172 // Define input and output slots here (never in the dummy constructor)
173 // Input slot #0 works with a TChain - it is connected to the default input container
174 // Output slot #1 writes into a TH1 container
175 DefineOutput(1, TList::Class()); // for output list
176 for(int i=0; i<cent_bins; i++){
177 h1_nClusters[i] = 0;
178 h1_M[i] = 0;
179 h1_M_mix[i] = 0;
180 h1_E[i] = 0;
181 h1_dR_ClustTrk[i] = 0;
182 h1_Pi0TruthPt[i] = 0;
183 h1_PriPi0TruthPt[i] = 0;
184 h1_Pi0TruthPtEmcal[i] = 0;
185 h1_PriPi0TruthPtEmcal[i] = 0;
186 h1_Pi0TruthPtPhi2piEta065[i] = 0;
187 h1_Pi0TruthPtPhi2piEta1[i] = 0;
188 h1_TruthPhotonsEmcal[i] = 0;
189 h1_PhotonsEmcal[i] = 0;
190 h1_PhotonsNCellsCut[i] = 0;
191 h1_PhotonsTrackMatchCut[i] = 0;
192 h1_PhotonsAllCut[i] = 0;
193 h1_dR_RealMC[i] = 0;
194 h1_Chi2[i] = 0;
195 h1_nTrkMatch[i] = 0;
196 h1_nCells[i] = 0;
197 h1_ClusterDisp[i] = 0;
198 h2_Ellipse[i] = 0;
199 h2_EtaPt[i] = 0;
200 h3_MptAsymm[i] = 0;
201 h3_MptAsymm_mix[i] = 0;
202 h2_dphi_deta[i] = 0;
203 h2_dphi_deta_mix[i] = 0;
204 h2_DispRes[i] = 0;
205 h2_cells_M02[i] = 0;
206 }
207}
208
209//________________________________________________________________________
210AliAnalysisTaskEMCALMesonGGSDMpPb::~AliAnalysisTaskEMCALMesonGGSDMpPb()
211{
212 // Destructor. Clean-up the output list, but not the histograms that are put inside
213 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
214 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
215 delete fOutput;
216 }
217 delete fTrackCuts;
218}
219
220//________________________________________________________________________
221void AliAnalysisTaskEMCALMesonGGSDMpPb::UserCreateOutputObjects()
222{
223 // Create histograms
224 // Called once (on the worker node)
225
226 fOutput = new TList();
227 fOutput->SetOwner(); // IMPORTANT!
228
229 fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
230
231 cout << "__________AliAnalysisTaskEMCALMesonGGSDMpPb: Input settings__________" << endl;
232 cout << " fMcMode: " << fMcMode << endl;
233 cout << " fRecalibrator: " << fRecalibrator << endl;
234 cout << " dRmin_ClustTrack: " << fdRmin_ClustTrack << endl;
235 cout << " phi range: " << fPhimin << ", " << fPhimax << endl;
236 cout << " eta range: " << fEtamin << ", " << fEtamax << endl;
237 cout << " number of zvtx bins: " << zvtx_bins << endl;
238 cout << " number of mult bins: " << mult_bins << endl;
239 cout << " poolDepth: " << poolDepth << endl;
240 cout << endl;
241
242 char saythis1[500];
243 char saythis2[500];
244
245 double TotalNBins = 0.0;
246
247 // Create histograms
248 Int_t nClustersbins = 501;
249 Float_t nClusterslow = -0.5, nClustersup = 500.5;
250 for(int i=0; i<cent_bins; i++){
251 sprintf(saythis1,"h1_nClusters_%d",i);
252 sprintf(saythis2,"# of clusters");
253 h1_nClusters[i] = new TH1F(saythis1, saythis2, nClustersbins, nClusterslow, nClustersup);
254 h1_nClusters[i]->GetXaxis()->SetTitle("number of clusters/evt");
255 h1_nClusters[i]->GetYaxis()->SetTitle("counts");
256 h1_nClusters[i]->SetMarkerStyle(kFullCircle);
257 TotalNBins+=nClustersbins;
258 }
259
260 Int_t nZvertexbins = 501;
261 Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
262 h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
263 h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
264 h1_zvtx->GetYaxis()->SetTitle("counts");
265 h1_zvtx->SetMarkerStyle(kFullCircle);
266 TotalNBins+=nZvertexbins;
267
268 h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
269 TotalNBins+=1001;
270
271 h1_centrality = new TH1F("h1_centrality", "centrality", 1001,-0.1,100.1);
272 TotalNBins+=1001;
273
274 Int_t Mbins = 3000;
275 Float_t Mlow = 0.0, Mup = 3.0;
276 for(int i=0; i<cent_bins; i++){
277 sprintf(saythis1,"h1_M_%d",i);
278 sprintf(saythis2,"Invariant Mass");
279 h1_M[i] = new TH1F(saythis1, saythis2, Mbins, Mlow, Mup);
280 h1_M[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
281 h1_M[i]->GetYaxis()->SetTitle("counts");
282 h1_M[i]->SetMarkerStyle(kFullCircle);
283 TotalNBins+=Mbins;
284 }
285
286 for(int i=0; i<cent_bins; i++){
287 sprintf(saythis1,"h1_M_mix_%d",i);
288 sprintf(saythis2,"Invariant Mass (mixed events)");
289 h1_M_mix[i] = new TH1F(saythis1, saythis2, Mbins, Mlow, Mup);
290 h1_M_mix[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
291 h1_M_mix[i]->GetYaxis()->SetTitle("counts");
292 h1_M_mix[i]->SetMarkerStyle(kFullCircle);
293 TotalNBins+=Mbins;
294 }
295
296 Int_t ptbins = 2000;
297 Float_t ptlow = 0.0, ptup = 20.0;
298 Int_t Ebins = 500;
299 Float_t Elow = 0.0, Eup = 20.0;
300 for(int i=0; i<cent_bins; i++){
301 sprintf(saythis1,"h1_E_%d",i);
302 sprintf(saythis2,"Cluster Energy in EMCal");
303 h1_E[i] = new TH1F(saythis1, saythis2, Ebins, Elow, Eup);
304 h1_E[i]->GetXaxis()->SetTitle("E [GeV]");
305 h1_E[i]->GetYaxis()->SetTitle("counts");
306 h1_E[i]->SetMarkerStyle(kFullCircle);
307 TotalNBins+=Ebins;
308 }
309
310 h2_PhiEtaCluster = new TH2F("h2_PhiEtaCluster", "cluster phi vs eta", 400,1.362,3.178, 300,-0.728,0.728);
311 h2_PhiEtaCluster->GetXaxis()->SetTitle("#phi [rad]");
312 h2_PhiEtaCluster->GetYaxis()->SetTitle("#eta");
313 h2_PhiEtaCluster->GetZaxis()->SetTitle("hits");
314 h2_PhiEtaCluster->SetMarkerStyle(kFullCircle);
315 TotalNBins+=400*300;
316
317 h2_PhiEtaClusterCut = new TH2F("h2_PhiEtaClusterCut", "cluster phi vs eta (after cuts)", 400,1.362,3.178, 300,-0.728,0.728);
318 h2_PhiEtaClusterCut->GetXaxis()->SetTitle("#phi [rad]");
319 h2_PhiEtaClusterCut->GetYaxis()->SetTitle("#eta");
320 h2_PhiEtaClusterCut->GetZaxis()->SetTitle("hits");
321 h2_PhiEtaClusterCut->SetMarkerStyle(kFullCircle);
322 TotalNBins+=400*300;
323
324// eta binning
325 Double_t EtaBins[97] = {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465};
326
327 // phi binning
328 Double_t PhiBins[125] = {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732, 1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128};
329
330 h2_PhiEtaMaxCell = new TH2F("h2_PhiEtaMaxCell", "maxcell phi vs eta", 124,PhiBins, 96,EtaBins);
331 h2_PhiEtaMaxCell->GetXaxis()->SetTitle("#phi [rad]");
332 h2_PhiEtaMaxCell->GetYaxis()->SetTitle("#eta");
333 h2_PhiEtaMaxCell->GetZaxis()->SetTitle("hits");
334 h2_PhiEtaMaxCell->SetMarkerStyle(kFullCircle);
335 TotalNBins+=96*124;
336
337 h2_PhiEtaMaxCellCut = new TH2F("h2_PhiEtaMaxCellCut", "maxcell phi vs eta (after cuts)", 124,PhiBins, 96,EtaBins);
338 h2_PhiEtaMaxCellCut->GetXaxis()->SetTitle("#phi [rad]");
339 h2_PhiEtaMaxCellCut->GetYaxis()->SetTitle("#eta");
340 h2_PhiEtaMaxCellCut->GetZaxis()->SetTitle("hits");
341 h2_PhiEtaMaxCellCut->SetMarkerStyle(kFullCircle);
342 TotalNBins+=96*124;
343
344 for(int i=0; i<cent_bins; i++){
345 sprintf(saythis1,"h1_dR_ClustTrk_%d",i);
346 sprintf(saythis2,"Cluster-Track matching");
347 h1_dR_ClustTrk[i] = new TH1F(saythis1, saythis2, 5000, -0.01, 5);
348 h1_dR_ClustTrk[i]->GetXaxis()->SetTitle("dR [sqrt(d#phi^{2}+d#eta^{2})]");
349 h1_dR_ClustTrk[i]->GetYaxis()->SetTitle("N");
350 h1_dR_ClustTrk[i]->SetMarkerStyle(kFullCircle);
351 TotalNBins+=5000;
352 }
353
354 h2_gE_RecTruth = new TH2F("h2_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 500,0,2);
355 h2_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
356 h2_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
357 h2_gE_RecTruth->GetZaxis()->SetTitle("counts");
358 h2_gE_RecTruth->SetMarkerStyle(kFullCircle);
359 TotalNBins+=Ebins*500;
360
361 h2_eop_E = new TH2F("h2_eop_E","E/p vs E (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
362 h2_eop_E->GetXaxis()->SetTitle("cluster Energy [GeV]");
363 h2_eop_E->GetYaxis()->SetTitle("E/p");
364 TotalNBins+=Ebins*1200;
365
366 h2_eop_pT = new TH2F("h2_eop_pT","E/p vs p_{T} (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
367 h2_eop_pT->GetXaxis()->SetTitle("cluster Energy [GeV]");
368 h2_eop_pT->GetYaxis()->SetTitle("E/p");
369 TotalNBins+=Ebins*1200;
370
371 h2_E_time = new TH2F("h2_E_time","cluster energy vs time", Ebins, Elow, Eup, 1000,-1e-6,1e-6);
372 h2_E_time->GetXaxis()->SetTitle("cluster Energy [GeV]");
373 h2_E_time->GetYaxis()->SetTitle("time [s]");
374 TotalNBins+=Ebins*1000;
375
376 for(int i=0; i<cent_bins; i++){
377 sprintf(saythis1,"h1_Pi0TruthPt_%d",i);
378 sprintf(saythis2,"P_{T} distribution for Truth Pi0's");
379 h1_Pi0TruthPt[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
380 h1_Pi0TruthPt[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
381 h1_Pi0TruthPt[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
382 h1_Pi0TruthPt[i]->SetMarkerStyle(kFullCircle);
383 TotalNBins+=ptbins;
384 }
385
386 for(int i=0; i<cent_bins; i++){
387 sprintf(saythis1,"h1_PriPi0TruthPt_%d",i);
388 sprintf(saythis2,"P_{T} distribution for Truth Primary Pi0's");
389 h1_PriPi0TruthPt[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
390 h1_PriPi0TruthPt[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
391 h1_PriPi0TruthPt[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
392 h1_PriPi0TruthPt[i]->SetMarkerStyle(kFullCircle);
393 TotalNBins+=ptbins;
394 }
395
396 for(int i=0; i<cent_bins; i++){
397 sprintf(saythis1,"h1_Pi0TruthPtEmcal_%d",i);
398 sprintf(saythis2,"P_{T} distribution for Truth Pi0's (hit EMCal)");
399 h1_Pi0TruthPtEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
400 h1_Pi0TruthPtEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
401 h1_Pi0TruthPtEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
402 h1_Pi0TruthPtEmcal[i]->SetMarkerStyle(kFullCircle);
403 TotalNBins+=ptbins;
404 }
405
406 for(int i=0; i<cent_bins; i++){
407 sprintf(saythis1,"h1_PriPi0TruthPtEmcal_%d",i);
408 sprintf(saythis2,"P_{T} distribution for Truth Primary Pi0's (hit EMCal)");
409 h1_PriPi0TruthPtEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
410 h1_PriPi0TruthPtEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
411 h1_PriPi0TruthPtEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
412 h1_PriPi0TruthPtEmcal[i]->SetMarkerStyle(kFullCircle);
413 TotalNBins+=ptbins;
414 }
415
416 for(int i=0; i<cent_bins; i++){
417 sprintf(saythis1,"h1_Pi0TruthPtPhi2piEta065_%d",i);
418 sprintf(saythis2,"P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]");
419 h1_Pi0TruthPtPhi2piEta065[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
420 h1_Pi0TruthPtPhi2piEta065[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
421 h1_Pi0TruthPtPhi2piEta065[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
422 h1_Pi0TruthPtPhi2piEta065[i]->SetMarkerStyle(kFullCircle);
423 TotalNBins+=ptbins;
424 }
425
426 for(int i=0; i<cent_bins; i++){
427 sprintf(saythis1,"h1_Pi0TruthPtPhi2piEta1_%d",i);
428 sprintf(saythis2,"P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]");
429 h1_Pi0TruthPtPhi2piEta1[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
430 h1_Pi0TruthPtPhi2piEta1[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
431 h1_Pi0TruthPtPhi2piEta1[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
432 h1_Pi0TruthPtPhi2piEta1[i]->SetMarkerStyle(kFullCircle);
433 TotalNBins+=ptbins;
434 }
435
436 h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
437 h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
438 h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
439 TotalNBins+=380*200;
440
441 h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
442 h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
443 h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
444 TotalNBins+=380*200;
445
446 h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-1.5,1.5);
447 h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
448 h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
449 TotalNBins+=380*150;
450
451 h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
452 h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
453 h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
454 TotalNBins+=380*150;
455
456 h2_Pi0TruthPhiEta_Phi2piEta065 = new TH2F("h2_Pi0TruthPhiEta_Phi2piEta065",
457 "Pi0Truth Phi vs Eta [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", 380,-0.02,6.30, 150,-5,5);
458 h2_Pi0TruthPhiEta_Phi2piEta065->GetXaxis()->SetTitle("#phi [rad]");
459 h2_Pi0TruthPhiEta_Phi2piEta065->GetYaxis()->SetTitle("#eta ");
460 TotalNBins+=380*150;
461
462 h2_Pi0TruthPhiEta_Phi2piEta1 = new TH2F("h2_Pi0TruthPhiEta_Phi2piEta1",
463 "Pi0Truth Phi vs Eta [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", 380,-0.02,6.30, 150,-1.5,1.5);
464 h2_Pi0TruthPhiEta_Phi2piEta1->GetXaxis()->SetTitle("#phi [rad]");
465 h2_Pi0TruthPhiEta_Phi2piEta1->GetYaxis()->SetTitle("#eta ");
466 TotalNBins+=380*150;
467
468 for(int i=0; i<cent_bins; i++){
469 sprintf(saythis1,"h1_TruthPhotonsEmcal_%d",i);
470 sprintf(saythis2,"P_{T} distribution for photons (in EMCal)");
471 h1_TruthPhotonsEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
472 h1_TruthPhotonsEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
473 h1_TruthPhotonsEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
474 h1_TruthPhotonsEmcal[i]->SetMarkerStyle(kFullCircle);
475 TotalNBins+=ptbins;
476 }
477
478 h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
479 "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
480 h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
481 h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
482 h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
483 TotalNBins+=380*150;
484
485 for(int i=0; i<cent_bins; i++){
486 sprintf(saythis1,"h1_PhotonsEmcal_%d",i);
487 sprintf(saythis2,"P_{T} distribution for photons (in EMCal)");
488 h1_PhotonsEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
489 h1_PhotonsEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
490 h1_PhotonsEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
491 h1_PhotonsEmcal[i]->SetMarkerStyle(kFullCircle);
492 TotalNBins+=ptbins;
493 }
494
495 for(int i=0; i<cent_bins; i++){
496 sprintf(saythis1,"h1_PhotonsNCellsCut_%d",i);
497 sprintf(saythis2,"P_{T} distribution for #gamma's that survive NCells cut");
498 h1_PhotonsNCellsCut[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
499 h1_PhotonsNCellsCut[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
500 h1_PhotonsNCellsCut[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
501 h1_PhotonsNCellsCut[i]->SetMarkerStyle(kFullCircle);
502 TotalNBins+=ptbins;
503 }
504
505 for(int i=0; i<cent_bins; i++){
506 sprintf(saythis1,"h1_PhotonsTrackMatchCut_%d",i);
507 sprintf(saythis2,"P_{T} distribution for #gamma's that survive TrackMatch cut");
508 h1_PhotonsTrackMatchCut[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
509 h1_PhotonsTrackMatchCut[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
510 h1_PhotonsTrackMatchCut[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
511 h1_PhotonsTrackMatchCut[i]->SetMarkerStyle(kFullCircle);
512 TotalNBins+=ptbins;
513 }
514
515 for(int i=0; i<cent_bins; i++){
516 sprintf(saythis1,"h1_PhotonsAllCut_%d",i);
517 sprintf(saythis2,"P_{T} distribution for #gamma's that survive All cut");
518 h1_PhotonsAllCut[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
519 h1_PhotonsAllCut[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
520 h1_PhotonsAllCut[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
521 h1_PhotonsAllCut[i]->SetMarkerStyle(kFullCircle);
522 TotalNBins+=ptbins;
523 }
524
525 h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
526 "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 100,-1.0,1.0);
527 h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
528 h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
529 TotalNBins+=380*100;
530
531 for(int i=0; i<cent_bins; i++){
532 sprintf(saythis1,"h1_dR_RealMC_%d",i);
533 sprintf(saythis2,"P_{T} distribution for #gamma's that survive All cut");
534 h1_dR_RealMC[i] = new TH1F(saythis1, saythis2, 2000, -0.01, 10);
535 h1_dR_RealMC[i]->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
536 h1_dR_RealMC[i]->GetYaxis()->SetTitle("N");
537 h1_dR_RealMC[i]->SetMarkerStyle(kFullCircle);
538 TotalNBins+=2000;
539 }
540
541 Int_t chi2bins = 100;
542 Float_t chi2low = -2, chi2up = 2;
543 for(int i=0; i<cent_bins; i++){
544 sprintf(saythis1,"h1_Chi2_%d",i);
545 sprintf(saythis2,"#chi^{2} distribution for reconstructed");
546 h1_Chi2[i] = new TH1F(saythis1,saythis2,chi2bins, chi2low, chi2up);
547 h1_Chi2[i]->GetXaxis()->SetTitle("#chi^{2}");
548 h1_Chi2[i]->GetYaxis()->SetTitle("counts");
549 TotalNBins+=chi2bins;
550 }
551
552 for(int i=0; i<cent_bins; i++){
553 sprintf(saythis1,"h1_nTrkMatch_%d",i);
554 sprintf(saythis2,"number of matched tracks");
555 h1_nTrkMatch[i] = new TH1F(saythis1,saythis2,14, -1.5, 5.5);
556 h1_nTrkMatch[i]->GetXaxis()->SetTitle("nTracksMatched");
557 h1_nTrkMatch[i]->GetYaxis()->SetTitle("counts");
558 TotalNBins+=14;
559 }
560
561 for(int i=0; i<cent_bins; i++){
562 sprintf(saythis1,"h1_ClusterDisp_%d",i);
563 sprintf(saythis2,"Dispersion of CaloCluster");
564 h1_ClusterDisp[i] = new TH1F(saythis1,saythis2,1000, -1, 3);
565 h1_ClusterDisp[i]->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
566 h1_ClusterDisp[i]->GetYaxis()->SetTitle("counts");
567 TotalNBins+=1000;
568 }
569
570 for(int i=0; i<cent_bins; i++){
571 sprintf(saythis1,"h2_Ellipse_%d",i);
572 sprintf(saythis2,"Ellipse axis M20 vs M02");
573 h2_Ellipse[i] = new TH2F(saythis1,saythis2,500, -0.01, 1, 500, -0.01, 1);
574 h2_Ellipse[i]->GetXaxis()->SetTitle("cluster->GetM20()");
575 h2_Ellipse[i]->GetYaxis()->SetTitle("cluster->GetM02()");
576 h2_Ellipse[i]->GetZaxis()->SetTitle("counts");
577 TotalNBins+=500*500;
578 }
579
580 Int_t etabins = 150;
581 Float_t etalow = -1.5, etaup = 1.5;
582 for(int i=0; i<cent_bins; i++){
583 sprintf(saythis1,"h2_EtaPt_%d",i);
584 sprintf(saythis2,"Cluster Energy vs ");
585 h2_EtaPt[i] = new TH2F(saythis1,saythis2,etabins, etalow, etaup, ptbins, ptlow, ptup);
586 h2_EtaPt[i]->GetXaxis()->SetTitle("E [GeV]");
587 h2_EtaPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
588 TotalNBins+=etabins*ptbins;
589 }
590
591 for(int i=0; i<cent_bins; i++){
592 sprintf(saythis1,"h3_MptAsymm_%d",i);
593 sprintf(saythis2,"mass vs p_{T} vs Asymm cut");
594 h3_MptAsymm[i] = new TH3F(saythis1,saythis2,Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
595 h3_MptAsymm[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
596 h3_MptAsymm[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
597 h3_MptAsymm[i]->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
598 TotalNBins+=Mbins*ptbins*3.0;
599 }
600
601 for(int i=0; i<cent_bins; i++){
602 sprintf(saythis1,"h3_MptAsymm_mix_%d",i);
603 sprintf(saythis2,"mass vs p_{T} vs Asymm cut (mixed events)");
604 h3_MptAsymm_mix[i] = new TH3F(saythis1,saythis2,Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
605 h3_MptAsymm_mix[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
606 h3_MptAsymm_mix[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
607 h3_MptAsymm_mix[i]->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
608 TotalNBins+=Mbins*ptbins*3.0;
609 }
610
611 for(int i=0; i<cent_bins; i++){
612 sprintf(saythis1,"h2_dphi_deta_%d",i);
613 sprintf(saythis2,"#Delta#phi vs #Delta#eta");
614 h2_dphi_deta[i] = new TH2F(saythis1,saythis2, 349,-1.5,5, 400,-2.0,2.0);
615 h2_dphi_deta[i]->GetXaxis()->SetTitle("#Delta#phi");
616 h2_dphi_deta[i]->GetYaxis()->SetTitle("#Delta#eta");
617 TotalNBins+=349*400;
618 }
619
620 for(int i=0; i<cent_bins; i++){
621 sprintf(saythis1,"h2_dphi_deta_mix_%d",i);
622 sprintf(saythis2,"#Delta#phi vs #Delta#eta (mixed events)");
623 h2_dphi_deta_mix[i] = new TH2F(saythis1,saythis2, 349,-1.5,5, 400,-2.0,2.0);
624 h2_dphi_deta_mix[i]->GetXaxis()->SetTitle("#Delta#phi");
625 h2_dphi_deta_mix[i]->GetYaxis()->SetTitle("#Delta#eta");
626 TotalNBins+=349*400;
627 }
628
629 for(int i=0; i<cent_bins; i++){
630 sprintf(saythis1,"h2_DispRes_%d",i);
631 sprintf(saythis2,"zvtx info");
632 h2_DispRes[i] = new TH2F(saythis1, saythis2, 500,-0.01,1, 500,-0.1,2);
633 h2_DispRes[i]->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
634 h2_DispRes[i]->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
635 h2_DispRes[i]->GetZaxis()->SetTitle("counts");
636 TotalNBins+=500*500;
637 }
638
639 for(int i=0; i<cent_bins; i++){
640 sprintf(saythis1,"h2_cells_M02_%d",i);
641 sprintf(saythis2,"nCells vs M02");
642 h2_cells_M02[i] = new TH2F(saythis1, saythis2, 204,-1.5,100.5, 500,-1,1.5);
643 h2_cells_M02[i]->GetXaxis()->SetTitle("nCells");
644 h2_cells_M02[i]->GetYaxis()->SetTitle("M02");
645 h2_cells_M02[i]->GetZaxis()->SetTitle("counts");
646 TotalNBins+=204*500;
647 }
648
649 cout << endl << "Total number of bins in booked histograms: " << TotalNBins << endl << endl;
650
651 // Initialize helper class (for vertex selection & pile up correction)
652 fHelperClass = new AliAnalysisUtils();
653
654 //TFile *f = OpenFile(1);
655 //TDirectory::TContext context(f);
656
657 fOutput->Add(h1_zvtx);
658 fOutput->Add(h1_trigger);
659 fOutput->Add(h1_centrality);
660 fOutput->Add(h2_PhiEtaCluster);
661 fOutput->Add(h2_PhiEtaClusterCut);
662 fOutput->Add(h2_PhiEtaMaxCell);
663 fOutput->Add(h2_PhiEtaMaxCellCut);
664 fOutput->Add(h2_gE_RecTruth);
665 fOutput->Add(h2_eop_E);
666 fOutput->Add(h2_eop_pT);
667 fOutput->Add(h2_E_time);
668
669 for(int i=0; i<cent_bins; i++){
670 fOutput->Add(h1_nClusters[i]);
671 fOutput->Add(h1_M[i]);
672 fOutput->Add(h1_M_mix[i]);
673 fOutput->Add(h1_E[i]);
674 fOutput->Add(h1_dR_ClustTrk[i]);
675 fOutput->Add(h1_Pi0TruthPt[i]);
676 fOutput->Add(h1_PriPi0TruthPt[i]);
677 fOutput->Add(h1_Pi0TruthPtEmcal[i]);
678 fOutput->Add(h1_PriPi0TruthPtEmcal[i]);
679 fOutput->Add(h1_Pi0TruthPtPhi2piEta065[i]);
680 fOutput->Add(h1_Pi0TruthPtPhi2piEta1[i]);
681 fOutput->Add(h1_TruthPhotonsEmcal[i]);
682 fOutput->Add(h1_PhotonsEmcal[i]);
683 fOutput->Add(h1_PhotonsNCellsCut[i]);
684 fOutput->Add(h1_PhotonsTrackMatchCut[i]);
685 fOutput->Add(h1_PhotonsAllCut[i]);
686 fOutput->Add(h1_dR_RealMC[i]);
687 fOutput->Add(h1_Chi2[i]);
688 fOutput->Add(h1_nTrkMatch[i]);
689 fOutput->Add(h1_ClusterDisp[i]);
690 fOutput->Add(h2_Ellipse[i]);
691 fOutput->Add(h2_EtaPt[i]);
692 fOutput->Add(h3_MptAsymm[i]);
693 fOutput->Add(h3_MptAsymm_mix[i]);
694 fOutput->Add(h2_dphi_deta[i]);
695 fOutput->Add(h2_dphi_deta_mix[i]);
696 fOutput->Add(h2_DispRes[i]);
697 fOutput->Add(h2_cells_M02[i]);
698 }
699 fOutput->Add(h2_Pi0TruthPhiEta);
700 fOutput->Add(h2_PriPi0TruthPhiEta);
701 fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
702 fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
703 fOutput->Add(h2_Pi0TruthPhiEta_Phi2piEta065);
704 fOutput->Add(h2_Pi0TruthPhiEta_Phi2piEta1);
705 fOutput->Add(h2_TruthPhotonsPhiEta);
706 fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
707
708 // Post data for ALL output slots >0 here,
709 // To get at least an empty histogram
710 // 1 is the outputnumber of a certain weg of task 1
711 PostData(1, fOutput);
712}
713
714//________________________________________________________________________
715void AliAnalysisTaskEMCALMesonGGSDMpPb::UserExec(Option_t *)
716{
717 // Main loop Called for each event
718
719 AliMCEvent *mcEvent = MCEvent();
720 Bool_t isMC = bool(mcEvent);//is this the right way to do this?
721
722 TRandom3 randy; randy.SetSeed(0);
723 unsigned int iskip = -1;
724 TLorentzVector ParentMix;
725
726 double recalScale = 1.0;
727
728 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
729
730 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
731 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
732 if (!aodH && !esdH) Printf("ERROR: Could not get ESD or AODInputHandler");
733
734 if(esdH) fEsdEv = esdH->GetEvent();
735 else if(aodH) fAodEv = aodH->GetEvent();
736 else{
737 AliFatal("Neither ESD nor AOD event found");
738 return;
739 }
740
741 // get pointer to reconstructed event
742 AliVEvent *event = InputEvent();
743 if (!event){
744 AliError("Pointer == 0, this can not happen!"); return;}
745 //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
746 //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
747 //if (!fEsdEv){
748 //AliError("Cannot get the ESD event"); return;}
749
750 fHelperClass->SetCutOnZVertexSPD(kFALSE);//does the zvtx have to match the spd vertex?
751 fHelperClass->SetMaxVtxZ(1.0e6);//i set this myself later..
752 // simply makes sure that there is at least 1 contributer to the zvtx determination.
753 // this should only remove the *extra* events at zvtx==0.
754 if(!fHelperClass->IsVertexSelected2013pA(event))
755 return;
756
757 Int_t iTrigger = 0;
758 if (fEsdEv) iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
759 else if (fAodEv) iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
760
761 char saythis[500];
762 Int_t iTriggerBin = 0;
763 for(unsigned long j=0; j<TriggerList.size(); j++){
764 if(iTrigger==TriggerList[j])
765 iTriggerBin=j+1;
766 }
767 if(iTriggerBin==0){
768 TriggerList.push_back(iTrigger);
769 iTriggerBin=TriggerList.size();
770 }
771
772 h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
773 sprintf(saythis,"%d",iTrigger);
774 h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
775
776 Double_t centralityVZERO=0.0;
777 Int_t centBin = 0;
778 AliCentrality *aliCent=NULL;
779
780 Int_t nclusters=0;
781 if(fEsdEv){
782 //Int_t evtN = fEsdEv->GetEventNumberInFile();
783 //Int_t ntracks = fEsdEv->GetNumberOfTracks();
784 nclusters = fEsdEv->GetNumberOfCaloClusters();
785 aliCent = fEsdEv->GetCentrality();
786 }
787 else if(fAodEv){
788 //Int_t evtN = fAodEv->GetEventNumberInFile();
789 //Int_t ntracks = fAodEv->GetNumberOfTracks();
790 nclusters = fAodEv->GetNumberOfCaloClusters();
791 aliCent = fAodEv->GetCentrality();
792 }
793 //centBin = aliCent->GetCentralityClass10("V0M");
794 //centralityVZERO = aliCent->GetCentralityPercentile("V0M");
795 //centBin = aliCent->GetCentralityClass10("V0C");
796 //centralityVZERO = aliCent->GetCentralityPercentile("V0C");
797 centBin = aliCent->GetCentralityClass10("V0A");
798 centralityVZERO = aliCent->GetCentralityPercentile("V0A");
799
800 if (centralityVZERO<20.0)
801 centBin = 0;
802 else if(centralityVZERO<40.0)
803 centBin = 1;
804 else if(centralityVZERO<60.0)
805 centBin = 2;
806 else
807 centBin = 3;
808
809 //cout << "Centrality: " << centBin << " " << centralityVZERO << endl;
810
811 if (fEsdEv){
812 if(!(fEsdEv->GetPrimaryVertex()->GetStatus())) return;
813 }
814 //else if (fAodEv){
815 //if(!(fAodEv->GetPrimaryVertex()->GetStatus())) return;
816 //}
817
818 Double_t vertDisp=0.0;
819 Double_t vertZres=0.0;
820 Bool_t vertIsfromZ=0;
821 if (fEsdEv){
822 vertDisp = fEsdEv->GetPrimaryVertex()->GetDispersion();
823 vertZres = fEsdEv->GetPrimaryVertex()->GetZRes();
824 vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
825 }
826 else if (fAodEv){
827 vertDisp = 0;
828 vertZres = 0;
829 vertIsfromZ = 0;
830 }
831
832 h2_DispRes[centBin]->Fill(vertDisp, vertZres);
833 // if vertex is from spd vertexZ, require more stringent cut
834 if (vertIsfromZ) {
835 if (vertDisp>0.02 || vertZres>0.25 )
836 return; // bad vertex from VertexerZ
837 }
838
839 // EMCal cluster loop for reconstructed event
840 //numberofclusters set above!
841 TLorentzVector Photon1, Photon2, Parent;
842 Double_t vertex[3];
843 Double_t E1=0.0;
844 Double_t vertZ=0.0;
845 if (fEsdEv) vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
846 else if (fAodEv) vertZ = fAodEv->GetPrimaryVertex()->GetZ();
847
848 h1_zvtx->Fill(vertZ);
849 //zvertex cut:
850 if(fabs(vertZ)>10.0)
851 return;
852
853 h1_nClusters[centBin]->Fill(nclusters);
854 h1_centrality->Fill(centralityVZERO);
855
856 int izvtx = GetZvtxBin(vertZ);
857 int imult = GetMultBin(nclusters);
858
859 //cout << iskip << " " << izvtx << " " << imult << endl;
860 //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
861
862
863
864 //######################### ~~~~~~~~~~~ ##################################
865 //######################### STARTING MC ##################################
866 //######################### ~~~~~~~~~~~ ##################################
867
868 if(isMC){
869 int isPrimary = 0;
870
871 if (!mcEvent){
872 cout << "no MC event" << endl;
873 return;
874 }
875
876 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
877 if (!evtVtx)
878 return;
879
880 mcEvent->PreReadAll();
881
882 Int_t nTracksMC = mcEvent->GetNumberOfTracks();
883 Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
884 //cout << "We have " << nPTracksMC << " primaries of " << nTracksMC << " total tracks." << endl;
885
886 for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
887 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
888 if (!mcP)
889 continue;
890
891
892 // it's a pion !!
893 if(mcP->PdgCode() != 111)
894 continue;
895
896 /*
897 // primary particle
898 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
899 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
900 if(dR <= 0.01) isPrimary = 1;
901 else isPrimary = 0;
902 */
903
904 if(iTrack<nPTracksMC) isPrimary = 1;
905 else isPrimary = 0;
906
907 h1_Pi0TruthPt [centBin]->Fill(mcP->Pt());
908 h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
909
910 if(isPrimary==1){
911 h1_PriPi0TruthPt [centBin]->Fill(mcP->Pt());
912 h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
913 }
914
915 if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
916 continue;
917
918 h1_Pi0TruthPtPhi2piEta1 [centBin]->Fill(mcP->Pt());
919 h2_Pi0TruthPhiEta_Phi2piEta1->Fill(mcP->Phi(),mcP->Eta());
920
921 if(mcP->Eta()>fEtamin && mcP->Eta()<fEtamax){
922 h1_Pi0TruthPtPhi2piEta065 [centBin]->Fill(mcP->Pt());
923 h2_Pi0TruthPhiEta_Phi2piEta065->Fill(mcP->Phi(),mcP->Eta());
924 }
925
926
927 Int_t d1 = mcP->GetFirstDaughter();
928 Int_t d2 = mcP->GetLastDaughter();
929
930 if (d1<0) continue;
931 if (d2<0) d2=d1;
932 if (d2-d1 != 1) continue;
933
934 bool bacc = true;
935 bool binp = true;
936 for (Int_t i=d1;i<=d2;++i){
937 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(i));
938 Double_t eta_d = dmc->Eta();
939 Double_t phi_d = dmc->Phi();
940 if(!(dmc->PdgCode()==22)){
941 binp = false;
942 }
943 if(!(dmc->PdgCode()==22 && eta_d>fEtamin && eta_d<fEtamax && phi_d>fPhimin && phi_d<fPhimax)){
944 bacc = false;
945 }
946 }
947
948 if(binp && bacc){// 2 Photons hit the EMCAL!
949
950 for (Int_t j=d1;j<=d2;++j){//both truth photons.
951
952 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(j));
953 Double_t eta_d = dmc->Eta();
954 Double_t phi_d = dmc->Phi();
955
956 if( dmc->PdgCode()==22 &&
957 dmc->Eta()>fEtamin && dmc->Eta()<fEtamax &&
958 dmc->Phi()>fPhimin && dmc->Phi()<fPhimax ){
959 h1_TruthPhotonsEmcal[centBin]->Fill(dmc->Pt());
960 h2_TruthPhotonsPhiEta->Fill(dmc->Phi(),dmc->Eta());
961 }
962
963 for(int i=0; i<nclusters; i++) {
964
965 Bool_t matches_pion_photon = 0;
966
967 AliESDCaloCluster* esdCluster=NULL;
968 AliAODCaloCluster* aodCluster=NULL;
969 if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
970 else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
971
972 Double_t clustMC_phi, clustMC_eta;
973
974 if(fEsdEv){
975
976 if(esdCluster->IsEMCAL()){
977
978 Float_t pos[3] = {0,0,0};
979 esdCluster->GetPosition(pos);
980 TVector3 vpos(pos);
981 //h1_Phi->Fill(vpos.Phi());
982 clustMC_phi = vpos.Phi();
983 clustMC_eta = vpos.Eta();
984
985 Double_t dR = TMath::Sqrt((eta_d-clustMC_eta)*(eta_d-clustMC_eta) +
986 (phi_d-clustMC_phi)*(phi_d-clustMC_phi));
987 h1_dR_RealMC[centBin]->Fill(dR);
988 if(dR<=0.04) matches_pion_photon = 1;
989
990 vpos.Delete();
991 }
992 if(matches_pion_photon){
993 if(esdCluster->IsEMCAL()){
994 h1_PhotonsEmcal[centBin]->Fill(esdCluster->E());
995 h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
996 }
997 if(esdCluster->IsEMCAL() && esdCluster->GetNCells()>=2)
998 h1_PhotonsNCellsCut[centBin]->Fill(esdCluster->E());
999 if(esdCluster->IsEMCAL() && esdCluster->GetNTracksMatched()==0)
1000 h1_PhotonsTrackMatchCut[centBin]->Fill(esdCluster->E());
1001 if(esdCluster->IsEMCAL() && esdCluster->GetNCells()>=2 && esdCluster->GetNTracksMatched()==0)
1002 h1_PhotonsAllCut[centBin]->Fill(esdCluster->E());
1003 }//if(matches_pion_photon)
1004
1005 }//if(fEsdEv)
1006 else if(fAodEv){
1007
1008 if(aodCluster->IsEMCAL()){
1009
1010 Float_t pos[3] = {0,0,0};
1011 aodCluster->GetPosition(pos);
1012 TVector3 vpos(pos);
1013 //h1_Phi->Fill(vpos.Phi());
1014 clustMC_phi = vpos.Phi();
1015 clustMC_eta = vpos.Eta();
1016
1017 Double_t dR = TMath::Sqrt((eta_d-clustMC_eta)*(eta_d-clustMC_eta) +
1018 (phi_d-clustMC_phi)*(phi_d-clustMC_phi));
1019 h1_dR_RealMC[centBin]->Fill(dR);
1020 if(dR<=0.04) matches_pion_photon = 1;
1021
1022 vpos.Delete();
1023 }
1024 if(matches_pion_photon){
1025 if(aodCluster->IsEMCAL()){
1026 h1_PhotonsEmcal[centBin]->Fill(aodCluster->E());
1027 h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
1028 }
1029 if(aodCluster->IsEMCAL() && aodCluster->GetNCells()>=2)
1030 h1_PhotonsNCellsCut[centBin]->Fill(aodCluster->E());
1031 if(aodCluster->IsEMCAL() && aodCluster->GetNTracksMatched()==0)
1032 h1_PhotonsTrackMatchCut[centBin]->Fill(aodCluster->E());
1033 if(aodCluster->IsEMCAL() && aodCluster->GetNCells()>=2 && aodCluster->GetNTracksMatched()==0)
1034 h1_PhotonsAllCut[centBin]->Fill(aodCluster->E());
1035
1036 }//if(matches_pion_photon)
1037
1038 }//if(fAodEv)
1039
1040 }//loop over nclusters.
1041
1042 }//both truth photons.
1043
1044 }// 2 Photons hit the EMCAL!
1045
1046
1047 if(binp && bacc){// 2 Photons hit the EMCAL!
1048 h1_Pi0TruthPtEmcal [centBin]->Fill(mcP->Pt());
1049 h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1050
1051 if(isPrimary==1){
1052 h1_PriPi0TruthPtEmcal [centBin]->Fill(mcP->Pt());
1053 h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1054 }
1055
1056 }//2 photons hit the EMCAL!
1057
1058 }//for(nTracksMC)
1059
1060 }//if(isMC)
1061
1062 //######################### ~~~~~~~~~~~~ ##################################
1063 //######################### DONE WITH MC ##################################
1064 //######################### ~~~~~~~~~~~~ ##################################
1065
1066
1067 for(int i=0; i<nclusters; i++) {
1068
1069 AliESDCaloCluster* esdCluster=NULL;
1070 AliAODCaloCluster* aodCluster=NULL;
1071 if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
1072 else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
1073 if(!esdCluster && !aodCluster) {
1074 AliError(Form("ERROR: Could not retrieve any (ESD or AOD) Cluster %d",i));
1075 continue;
1076 }
1077
1078 if(fEsdEv){
1079
1080 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1081
1082 //uncomment this to do the track matching (1 of 3 lines, esd part)!!
1083 //Bool_t MatchesToTrack = 0;
1084 if(esdCluster->IsEMCAL()){
1085
1086 Float_t pos[3] = {0,0,0};
1087 Short_t maxCellID = -1;
1088 Float_t celleta, cellphi;
1089 esdCluster->GetPosition(pos);
1090 TVector3 clusterPosition(pos);
1091 h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1092 GetMaxCellEnergy(esdCluster, maxCellID);
1093 AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1094 fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1095 h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1096
1097 // _______________Track loop for reconstructed event_____________
1098 for(Int_t itrk = 0; itrk < fEsdEv->GetNumberOfTracks(); itrk++) {
1099 AliESDtrack* esdTrack = fEsdEv->GetTrack(itrk); // pointer to reconstructed to track
1100 if(!esdTrack) {
1101 AliError(Form("ERROR: Could not retrieve any (ESD) track %d",itrk));
1102 continue;
1103 }
1104
1105 Double_t posTrk[3] = {0,0,0};
1106 esdTrack->GetXYZ(posTrk);
1107 TVector3 vposTrk(posTrk);
1108
1109 Double_t fMass = 0.139;
1110 Double_t fStepSurface = 20.;
1111 Float_t etaproj, phiproj, pttrackproj;
1112
1113 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1114 if(!trackParam) continue;
1115 AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1116
1117 double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1118 (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1119
1120 h1_dR_ClustTrk[centBin]->Fill(dR_clusttrk);
1121
1122 //uncomment this to do the track matching (2 of 3 lines)!!
1123 //if(dR_clusttrk<fdRmin_ClustTrack)
1124 //MatchesToTrack = 1;
1125
1126 }//_____________________________nTracks__________________________
1127
1128 h2_cells_M02 [centBin]->Fill(esdCluster->GetNCells(),esdCluster->GetM02());
1129 h2_Ellipse [centBin]->Fill(esdCluster->GetM20(),esdCluster->GetM02());
1130 h1_Chi2 [centBin]->Fill(esdCluster->Chi2());//always -1.
1131 h1_nTrkMatch [centBin]->Fill(esdCluster->GetNTracksMatched());
1132 h1_ClusterDisp[centBin]->Fill(esdCluster->GetDispersion());
1133 h2_E_time ->Fill(esdCluster->E(),esdCluster->GetTOF());
1134
1135 TArrayI *TrackLabels = esdCluster->GetTracksMatched();
1136 if(TrackLabels){
1137 if(TrackLabels->GetSize()>0){
1138 Int_t trackindex = TrackLabels->At(0);
1139 AliESDtrack* matchingT = fEsdEv->GetTrack(trackindex); // pointer to reconstructed to track
1140
1141 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1142 h2_eop_E ->Fill(esdCluster->E()*recalScale, esdCluster->E()*recalScale/matchingT->P());
1143 h2_eop_pT->Fill(matchingT->Pt(), esdCluster->E()*recalScale/matchingT->P());
1144 }
1145 }
1146
1147 //uncomment this to do the track matching (3 of 3 lines)!!
1148 //if(isGoodEsdCluster(esdCluster) && !MatchesToTrack){
1149 if(isGoodEsdCluster(esdCluster)){
1150 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1151 E1 = esdCluster->E()*recalScale;// TOTAL HACK - JJ
1152 fEsdEv->GetVertex()->GetXYZ(vertex);
1153 esdCluster->GetMomentum(Photon1,vertex);
1154 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1155 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1156 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1157 Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1158 h1_E[centBin]->Fill(E1);
1159 h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1160 h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1161 }
1162 clusterPosition.Delete();
1163 }//if(esdCluster->isEMCAL())
1164 }//if(fEsdEv)
1165 else if(fAodEv){
1166
1167 recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1168
1169 //uncomment this to do the track matching (1 of 3 lines, aod part)!!
1170 //Bool_t MatchesToTrack = 0;
1171 if(aodCluster->IsEMCAL()){
1172
1173 Float_t pos[3] = {0,0,0};
1174 Short_t maxCellID = -1;
1175 Float_t celleta, cellphi;
1176 aodCluster->GetPosition(pos);
1177 TVector3 clusterPosition(pos);
1178 h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1179 GetMaxCellEnergy(aodCluster, maxCellID);
1180 AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1181 fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1182 h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1183
1184 // _______________Track loop for reconstructed event_____________
1185 for(Int_t itrk = 0; itrk < fAodEv->GetNumberOfTracks(); itrk++) {
f15c1f69 1186 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fAodEv->GetTrack(itrk));
1187 if(!aodTrack) AliFatal("Not a standard AOD"); // pointer to reconstructed to track
3a66195e 1188 if(!aodTrack) {
1189 AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk));
1190 continue;
1191 }
1192
1193 Double_t posTrk[3] = {0,0,0};
1194 aodTrack->GetXYZ(posTrk);
1195 TVector3 vposTrk(posTrk);
1196
1197 Double_t fMass = 0.139;
1198 Double_t fStepSurface = 20.;
1199 Float_t etaproj, phiproj, pttrackproj;
1200
1201 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam());
1202 if(!trackParam) continue;
1203 AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1204
1205 double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1206 (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1207
1208 h1_dR_ClustTrk[centBin]->Fill(dR_clusttrk);
1209
1210 //uncomment this to do the track matching (2 of 3 lines, aod part)!!
1211 //if(dR_clusttrk<fdRmin_ClustTrack)
1212 //MatchesToTrack = 1;
1213
1214
1215 }//_____________________________nTracks__________________________
1216
1217 h2_cells_M02 [centBin]->Fill(aodCluster->GetNCells(),aodCluster->GetM02());
1218 h2_Ellipse [centBin]->Fill(aodCluster->GetM20(),aodCluster->GetM02());
1219 h1_Chi2 [centBin]->Fill(aodCluster->Chi2());//always -1.
1220 h1_nTrkMatch [centBin]->Fill(aodCluster->GetNTracksMatched());
1221 h1_ClusterDisp[centBin]->Fill(aodCluster->GetDispersion());
1222 h2_E_time ->Fill(aodCluster->E(),aodCluster->GetTOF());
1223
1224 // #################################################
1225 // track matching eop histograms are handled here...
1226 // #################################################
1227
1228 //uncomment this to do the track matching (3 of 3 lines, aod part)!!
1229 //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){
1230 if(isGoodAodCluster(aodCluster)){
1231 recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1232 E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ
1233 fAodEv->GetVertex(0)->GetXYZ(vertex);
1234 aodCluster->GetMomentum(Photon1,vertex);
1235 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1236 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1237 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1238 Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1239 h1_E[centBin]->Fill(E1);
1240 h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1241 h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1242 }
1243 clusterPosition.Delete();
1244 }//if(aodCluster->IsEMCAL())
1245 }//if(fAodEv)
1246
1247 }//loop over nclusters.
1248
1249 //Make same event pions...
1250 for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1251 for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){
1252 Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j];
1253 Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1254 Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1255 Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/
1256 (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt());
1257 Int_t asymCut = 0;
1258 if (pairasym<0.1) asymCut = 1;
1259 else if(pairasym<0.7) asymCut = 2;
1260 else asymCut = 3;
1261
1262 h1_M [centBin]->Fill(Parent.M());
1263 h3_MptAsymm [centBin]->Fill(Parent.M(),Parent.Pt(),asymCut);
1264 h2_dphi_deta[centBin]->Fill(deltaphi,deltaeta);
1265 }
1266 }
1267
1268 //Make mixed event...
1269 for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1270 for(unsigned int ipool=1; ipool<poolDepth; ipool++){
1271 for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){
1272 iskip = randy.Integer(Photons[0][izvtx][imult].size());
1273 if(j==iskip) continue;
1274 Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j];
1275 Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1276 Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1277 Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/
1278 (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt());
1279 Int_t asymCut = 0;
1280 if (pairasym<0.1) asymCut = 1;
1281 else if(pairasym<0.7) asymCut = 2;
1282 else asymCut = 3;
1283
1284 h1_M_mix [centBin]->Fill(Parent.M());
1285 h3_MptAsymm_mix [centBin]->Fill(Parent.M(),Parent.Pt(),asymCut);
1286 h2_dphi_deta_mix[centBin]->Fill(deltaphi,deltaeta);
1287 }
1288 }
1289 }
1290
1291 for(int ipool=poolDepth-1; ipool>0; ipool--){
1292 Photons[ipool][izvtx][imult].clear();
1293 for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++)
1294 Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]);
1295 }
1296 Photons[0][izvtx][imult].clear();
1297
1298
1299
1300 // NEW HISTO should be filled before this point, as PostData puts the
1301 // information for this iteration of the UserExec in the container
1302 PostData(1, fOutput);
1303 }
1304
1305//________________________________________________________________________
1306void AliAnalysisTaskEMCALMesonGGSDMpPb::Terminate(Option_t *) //specify what you want to have done
1307{
1308 // Called once at the end of the query.
1309
1310}
1311
1312//________________________________________________________________________
1313Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetZvtxBin(Double_t vertZ)
1314{
1315
1316 int izvtx = -1;
1317
1318 if (vertZ<-35)
1319 izvtx=0;
1320 else if(vertZ<-30)
1321 izvtx=1;
1322 else if(vertZ<-25)
1323 izvtx=2;
1324 else if(vertZ<-20)
1325 izvtx=3;
1326 else if(vertZ<-15)
1327 izvtx=4;
1328 else if(vertZ<-10)
1329 izvtx=5;
1330 else if(vertZ< -5)
1331 izvtx=6;
1332 else if(vertZ< 0)
1333 izvtx=7;
1334 else if(vertZ< 5)
1335 izvtx=8;
1336 else if(vertZ< 10)
1337 izvtx=9;
1338 else if(vertZ< 15)
1339 izvtx=10;
1340 else if(vertZ< 20)
1341 izvtx=11;
1342 else if(vertZ< 25)
1343 izvtx=12;
1344 else if(vertZ< 30)
1345 izvtx=13;
1346 else if(vertZ< 35)
1347 izvtx=14;
1348 else
1349 izvtx=15;
1350
1351 return izvtx;
1352}
1353
1354//________________________________________________________________________
1355Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetMultBin(Int_t mult){
1356
1357 int imult = -1;
1358
1359 if (mult<2)
1360 imult=0;
1361 else if(mult<25)
1362 imult=mult-2;
1363 else
1364 imult=24;
1365
1366 return imult;
1367}
1368
1369//________________________________________________________________________
1370Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::isGoodEsdCluster(AliESDCaloCluster* esdclust){
1371
1372 int pass = 1;
1373 int nMinCells = 2;
1374 double MinE = 0.4;
1375 //double MinErat = 0;
1376 //double MinEcc = 0;
1377
1378 if (!esdclust)
1379 pass = 0;
1380 if (!esdclust->IsEMCAL())
1381 pass = 0;
1382 if (esdclust->E()<MinE)
1383 pass = 0;
1384 if (esdclust->GetNCells()<nMinCells)
1385 pass = 0;
1386 //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
1387 //pass = 0;
1388 //if (esdclust->Chi2()<MinEcc) // eccentricity cut
1389 //pass = 0;//this is always -1.
1390
1391 //if(esdclust->GetM02()<0.1)
1392 // pass = 0;
1393 //if(esdclust->GetM02()>0.5)
1394 // pass = 0;
1395
1396 Float_t pos[3] = {0,0,0};
1397 esdclust->GetPosition(pos);
1398 TVector3 clusterPosition(pos);
1399 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1400 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1401 pass = 0;
1402 clusterPosition.Delete();
1403
1404 //DOING THIS BY HAND NOW...
1405 //if(!esdclust->GetNTracksMatched()==0)
1406 //pass = 0;
1407
1408 return pass;
1409}
1410
1411//________________________________________________________________________
1412Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::isGoodAodCluster(AliAODCaloCluster* aodclust){
1413
1414 int pass = 1;
1415 int nMinCells = 2;
1416 double MinE = 0.4;
1417 //double MinErat = 0;
1418 //double MinEcc = 0;
1419
1420 if (!aodclust)
1421 pass = 0;
1422 if (!aodclust->IsEMCAL())
1423 pass = 0;
1424 if (aodclust->E()<MinE)
1425 pass = 0;
1426 if (aodclust->GetNCells()<nMinCells)
1427 pass = 0;
1428 //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1429 //pass = 0;
1430 //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1431 //pass = 0;//this is always -1.
1432
1433 //if(aodclust->GetM02()<0.1)
1434 //pass = 0;
1435 //if(aodclust->GetM02()>0.5)
1436 //pass = 0;
1437
1438 Float_t pos[3] = {0,0,0};
1439 aodclust->GetPosition(pos);
1440 TVector3 clusterPosition(pos);
1441 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1442 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1443 pass = 0;
1444 clusterPosition.Delete();
1445
1446 //DOING THIS BY HAND NOW...
1447 //if(!aodclust->GetNTracksMatched()==0)
1448 //pass = 0;
1449
1450 return pass;
1451}
1452
1453//________________________________________________________________________
1454Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1455
1456 double dphi = p1.Phi() - p2.Phi();
1457
1458 if(dphi<0.5*TMath::Pi())
1459 dphi = dphi + 2.0*TMath::Pi();
1460
1461 if(dphi>1.5*TMath::Pi())
1462 dphi = dphi - 2.0*TMath::Pi();
1463
1464 return dphi;
1465}
1466
1467//________________________________________________________________________
1468Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1469
1470 double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1471
1472 return deta;
1473}
1474
1475
1476//________________________________________________________________________
1477Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1478
1479 double recalibfactor = 0.0;
1480
1481 if(iCalib==0){// no recalibration!
1482 recalibfactor = 1.0;
1483 }
1484 else if(iCalib==1){// just a scale factor:
1485 recalibfactor = 0.984;
1486 }
1487 else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1488 Double_t p[3] = {0.96968, -2.68720, -0.831607};
1489 recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1490 }
1491 else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1492 Double_t p[7] = {1.00000e+00, 3.04925e-02, 4.69043e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.00046e+00};
1493 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1494 }
1495 else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1496 Double_t p[7] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01};
1497 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1498 }
1499 else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1500 Double_t p[10] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01, 0.96968, -2.68720, -0.831607};
1501 recalibfactor = ( (p[6]/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) ) / ( p[7] + exp(p[8] + p[9]*energy/2.0) );
1502 }
1503 else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1504 Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1505 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1506 }
1507 else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1508 Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1509 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1510 }
1511 else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1512 Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.964, -3.132, -0.435};
1513 //Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.96968, -2.68720, -0.831607};//same SDM piece as iCalib==2
1514 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) * (p[7] + exp(p[8]+p[9]*energy*2.0));
1515 }
1516 else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1517 Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1518 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1519 }
1520 else if(iCalib==10){// Jason played with test beam data
1521 Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1522 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1523 }
1524 else if(iCalib==11){// Jason played with test beam MC
1525 Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1526 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1527 }
1528
1529 return recalibfactor;
1530}
1531
1532
1533//________________________________________________________________________
1534Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1535{
1536 // Get maximum energy of attached cell.
1537
1538 id = -1;
1539 AliVCaloCells *fVCells=NULL;
1540 if(fEsdEv) fVCells = fEsdEv->GetEMCALCells();
1541 else if(fAodEv) fVCells = fAodEv->GetEMCALCells();
1542 if(!fVCells)
1543 return 0;
1544
1545 Double_t maxe = 0;
1546 Int_t ncells = cluster->GetNCells();
1547 for (Int_t i=0; i<ncells; i++) {
1548 Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1549 if (e>maxe) {
1550 maxe = e;
1551 id = cluster->GetCellAbsId(i);
1552 }
1553 }
1554 return maxe;
1555}
1556
1557
1558//________________________________________________________________________
1559Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){
1560
1561 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1562
1563 Int_t nPTracks= mcEvent->GetNumberOfPrimaries();
1564
1565 Int_t isPhysPrimary = 1;
1566 Int_t ismHF = 0;
1567 Int_t ismLongLivedOrK = 0;
1568
1569 if(mcP->GetMother()<0)//if it has no mother...
1570 return isPhysPrimary;
1571
1572 Int_t imTrack = mcP->GetMother();
1573 AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1574
1575 if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark...
1576 return isPhysPrimary;
1577
1578
1579 //############################################
1580 //get the PDG digits....
1581 int num = mcPm->PdgCode();
1582 int RevDigits[10] = {0};
1583 int nDigits = 0;
1584 while (num >= 1){
1585 RevDigits[nDigits++] = num%10;
1586 num = num / 10;
1587 }
1588 //##############################################
1589
1590
1591 if(RevDigits[3]>3)//Baryons
1592 ismHF = 1;
1593 else if(RevDigits[2]>3)//Mesons
1594 ismHF = 1;
1595
1596 ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1597
1598 if(!ismHF && ismLongLivedOrK)
1599 isPhysPrimary = 0;
1600 else{ // check grandmother, greatgrandmothers, etc...
1601 while(imTrack >= nPTracks){
1602
1603 if(mcPm->GetMother()<0)//if it has no mother...
1604 break;
1605
1606 if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark...
1607 return isPhysPrimary;
1608
1609 imTrack = mcPm->GetMother();
1610 mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1611
1612 //############################################
1613 //get the PDG digits....
1614 num = mcPm->PdgCode();
1615 for(int i=0; i<10; i++) RevDigits[i] = 0;
1616 nDigits = 0;
1617 while (num >= 1){
1618 RevDigits[nDigits++] = num%10;
1619 num = num / 10;
1620 }
1621 //##############################################
1622 if(RevDigits[3]>3)//Baryons
1623 ismHF = 1;
1624 else if(RevDigits[2]>3)//Mesons
1625 ismHF = 1;
1626
1627 ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1628
1629 if(!ismHF && ismLongLivedOrK)
1630 isPhysPrimary = 0;
1631
1632 }//while( >=nPTracks)
1633 }
1634
1635 return isPhysPrimary;
1636}
1637
1638
1639//________________________________________________________________________
1640Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::IsLongLivedOrK(Int_t MyPDGcode){
1641
1642 Int_t MyFlag = 0;
1643
1644 if(
1645 (TMath::Abs(MyPDGcode) == 22 ) || // Photon
1646 (TMath::Abs(MyPDGcode) == 11 ) || // Electron
1647 (TMath::Abs(MyPDGcode) == 13 ) || // Muon(-)
1648 (TMath::Abs(MyPDGcode) == 211 ) || // Pion
1649 (TMath::Abs(MyPDGcode) == 321 ) || // Kaon
1650 (TMath::Abs(MyPDGcode) == 310 ) || // K0s
1651 (TMath::Abs(MyPDGcode) == 130 ) || // K0l
1652 (TMath::Abs(MyPDGcode) == 2212) || // Proton
1653 (TMath::Abs(MyPDGcode) == 2112) || // Neutron
1654 (TMath::Abs(MyPDGcode) == 3122) || // Lambda_0
1655 (TMath::Abs(MyPDGcode) == 3112) || // Sigma Minus
1656 (TMath::Abs(MyPDGcode) == 3222) || // Sigma Plus
1657 (TMath::Abs(MyPDGcode) == 3312) || // Xsi Minus
1658 (TMath::Abs(MyPDGcode) == 3322) || // Xsi
1659 (TMath::Abs(MyPDGcode) == 3334) || // Omega
1660 (TMath::Abs(MyPDGcode) == 12 ) || // Electron Neutrino
1661 (TMath::Abs(MyPDGcode) == 14 ) || // Muon Neutrino
1662 (TMath::Abs(MyPDGcode) == 16 ) ) // Tau Neutrino
1663 MyFlag = 1;
1664
1665 return MyFlag;
1666}