]> 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
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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++) {
1186 AliAODTrack* aodTrack = fAodEv->GetTrack(itrk); // pointer to reconstructed to track
1187 if(!aodTrack) {
1188 AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk));
1189 continue;
1190 }
1191
1192 Double_t posTrk[3] = {0,0,0};
1193 aodTrack->GetXYZ(posTrk);
1194 TVector3 vposTrk(posTrk);
1195
1196 Double_t fMass = 0.139;
1197 Double_t fStepSurface = 20.;
1198 Float_t etaproj, phiproj, pttrackproj;
1199
1200 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam());
1201 if(!trackParam) continue;
1202 AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1203
1204 double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1205 (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1206
1207 h1_dR_ClustTrk[centBin]->Fill(dR_clusttrk);
1208
1209 //uncomment this to do the track matching (2 of 3 lines, aod part)!!
1210 //if(dR_clusttrk<fdRmin_ClustTrack)
1211 //MatchesToTrack = 1;
1212
1213
1214 }//_____________________________nTracks__________________________
1215
1216 h2_cells_M02 [centBin]->Fill(aodCluster->GetNCells(),aodCluster->GetM02());
1217 h2_Ellipse [centBin]->Fill(aodCluster->GetM20(),aodCluster->GetM02());
1218 h1_Chi2 [centBin]->Fill(aodCluster->Chi2());//always -1.
1219 h1_nTrkMatch [centBin]->Fill(aodCluster->GetNTracksMatched());
1220 h1_ClusterDisp[centBin]->Fill(aodCluster->GetDispersion());
1221 h2_E_time ->Fill(aodCluster->E(),aodCluster->GetTOF());
1222
1223 // #################################################
1224 // track matching eop histograms are handled here...
1225 // #################################################
1226
1227 //uncomment this to do the track matching (3 of 3 lines, aod part)!!
1228 //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){
1229 if(isGoodAodCluster(aodCluster)){
1230 recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1231 E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ
1232 fAodEv->GetVertex(0)->GetXYZ(vertex);
1233 aodCluster->GetMomentum(Photon1,vertex);
1234 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1235 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1236 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1237 Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1238 h1_E[centBin]->Fill(E1);
1239 h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1240 h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1241 }
1242 clusterPosition.Delete();
1243 }//if(aodCluster->IsEMCAL())
1244 }//if(fAodEv)
1245
1246 }//loop over nclusters.
1247
1248 //Make same event pions...
1249 for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1250 for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){
1251 Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j];
1252 Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1253 Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1254 Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/
1255 (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt());
1256 Int_t asymCut = 0;
1257 if (pairasym<0.1) asymCut = 1;
1258 else if(pairasym<0.7) asymCut = 2;
1259 else asymCut = 3;
1260
1261 h1_M [centBin]->Fill(Parent.M());
1262 h3_MptAsymm [centBin]->Fill(Parent.M(),Parent.Pt(),asymCut);
1263 h2_dphi_deta[centBin]->Fill(deltaphi,deltaeta);
1264 }
1265 }
1266
1267 //Make mixed event...
1268 for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1269 for(unsigned int ipool=1; ipool<poolDepth; ipool++){
1270 for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){
1271 iskip = randy.Integer(Photons[0][izvtx][imult].size());
1272 if(j==iskip) continue;
1273 Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j];
1274 Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1275 Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1276 Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/
1277 (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt());
1278 Int_t asymCut = 0;
1279 if (pairasym<0.1) asymCut = 1;
1280 else if(pairasym<0.7) asymCut = 2;
1281 else asymCut = 3;
1282
1283 h1_M_mix [centBin]->Fill(Parent.M());
1284 h3_MptAsymm_mix [centBin]->Fill(Parent.M(),Parent.Pt(),asymCut);
1285 h2_dphi_deta_mix[centBin]->Fill(deltaphi,deltaeta);
1286 }
1287 }
1288 }
1289
1290 for(int ipool=poolDepth-1; ipool>0; ipool--){
1291 Photons[ipool][izvtx][imult].clear();
1292 for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++)
1293 Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]);
1294 }
1295 Photons[0][izvtx][imult].clear();
1296
1297
1298
1299 // NEW HISTO should be filled before this point, as PostData puts the
1300 // information for this iteration of the UserExec in the container
1301 PostData(1, fOutput);
1302 }
1303
1304//________________________________________________________________________
1305void AliAnalysisTaskEMCALMesonGGSDMpPb::Terminate(Option_t *) //specify what you want to have done
1306{
1307 // Called once at the end of the query.
1308
1309}
1310
1311//________________________________________________________________________
1312Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetZvtxBin(Double_t vertZ)
1313{
1314
1315 int izvtx = -1;
1316
1317 if (vertZ<-35)
1318 izvtx=0;
1319 else if(vertZ<-30)
1320 izvtx=1;
1321 else if(vertZ<-25)
1322 izvtx=2;
1323 else if(vertZ<-20)
1324 izvtx=3;
1325 else if(vertZ<-15)
1326 izvtx=4;
1327 else if(vertZ<-10)
1328 izvtx=5;
1329 else if(vertZ< -5)
1330 izvtx=6;
1331 else if(vertZ< 0)
1332 izvtx=7;
1333 else if(vertZ< 5)
1334 izvtx=8;
1335 else if(vertZ< 10)
1336 izvtx=9;
1337 else if(vertZ< 15)
1338 izvtx=10;
1339 else if(vertZ< 20)
1340 izvtx=11;
1341 else if(vertZ< 25)
1342 izvtx=12;
1343 else if(vertZ< 30)
1344 izvtx=13;
1345 else if(vertZ< 35)
1346 izvtx=14;
1347 else
1348 izvtx=15;
1349
1350 return izvtx;
1351}
1352
1353//________________________________________________________________________
1354Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetMultBin(Int_t mult){
1355
1356 int imult = -1;
1357
1358 if (mult<2)
1359 imult=0;
1360 else if(mult<25)
1361 imult=mult-2;
1362 else
1363 imult=24;
1364
1365 return imult;
1366}
1367
1368//________________________________________________________________________
1369Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::isGoodEsdCluster(AliESDCaloCluster* esdclust){
1370
1371 int pass = 1;
1372 int nMinCells = 2;
1373 double MinE = 0.4;
1374 //double MinErat = 0;
1375 //double MinEcc = 0;
1376
1377 if (!esdclust)
1378 pass = 0;
1379 if (!esdclust->IsEMCAL())
1380 pass = 0;
1381 if (esdclust->E()<MinE)
1382 pass = 0;
1383 if (esdclust->GetNCells()<nMinCells)
1384 pass = 0;
1385 //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
1386 //pass = 0;
1387 //if (esdclust->Chi2()<MinEcc) // eccentricity cut
1388 //pass = 0;//this is always -1.
1389
1390 //if(esdclust->GetM02()<0.1)
1391 // pass = 0;
1392 //if(esdclust->GetM02()>0.5)
1393 // pass = 0;
1394
1395 Float_t pos[3] = {0,0,0};
1396 esdclust->GetPosition(pos);
1397 TVector3 clusterPosition(pos);
1398 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1399 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1400 pass = 0;
1401 clusterPosition.Delete();
1402
1403 //DOING THIS BY HAND NOW...
1404 //if(!esdclust->GetNTracksMatched()==0)
1405 //pass = 0;
1406
1407 return pass;
1408}
1409
1410//________________________________________________________________________
1411Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::isGoodAodCluster(AliAODCaloCluster* aodclust){
1412
1413 int pass = 1;
1414 int nMinCells = 2;
1415 double MinE = 0.4;
1416 //double MinErat = 0;
1417 //double MinEcc = 0;
1418
1419 if (!aodclust)
1420 pass = 0;
1421 if (!aodclust->IsEMCAL())
1422 pass = 0;
1423 if (aodclust->E()<MinE)
1424 pass = 0;
1425 if (aodclust->GetNCells()<nMinCells)
1426 pass = 0;
1427 //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1428 //pass = 0;
1429 //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1430 //pass = 0;//this is always -1.
1431
1432 //if(aodclust->GetM02()<0.1)
1433 //pass = 0;
1434 //if(aodclust->GetM02()>0.5)
1435 //pass = 0;
1436
1437 Float_t pos[3] = {0,0,0};
1438 aodclust->GetPosition(pos);
1439 TVector3 clusterPosition(pos);
1440 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1441 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1442 pass = 0;
1443 clusterPosition.Delete();
1444
1445 //DOING THIS BY HAND NOW...
1446 //if(!aodclust->GetNTracksMatched()==0)
1447 //pass = 0;
1448
1449 return pass;
1450}
1451
1452//________________________________________________________________________
1453Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1454
1455 double dphi = p1.Phi() - p2.Phi();
1456
1457 if(dphi<0.5*TMath::Pi())
1458 dphi = dphi + 2.0*TMath::Pi();
1459
1460 if(dphi>1.5*TMath::Pi())
1461 dphi = dphi - 2.0*TMath::Pi();
1462
1463 return dphi;
1464}
1465
1466//________________________________________________________________________
1467Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1468
1469 double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1470
1471 return deta;
1472}
1473
1474
1475//________________________________________________________________________
1476Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1477
1478 double recalibfactor = 0.0;
1479
1480 if(iCalib==0){// no recalibration!
1481 recalibfactor = 1.0;
1482 }
1483 else if(iCalib==1){// just a scale factor:
1484 recalibfactor = 0.984;
1485 }
1486 else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1487 Double_t p[3] = {0.96968, -2.68720, -0.831607};
1488 recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1489 }
1490 else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1491 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};
1492 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1493 }
1494 else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1495 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};
1496 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1497 }
1498 else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1499 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};
1500 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) );
1501 }
1502 else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1503 Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1504 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1505 }
1506 else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1507 Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1508 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1509 }
1510 else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1511 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};
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.96968, -2.68720, -0.831607};//same SDM piece as iCalib==2
1513 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));
1514 }
1515 else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1516 Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1517 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1518 }
1519 else if(iCalib==10){// Jason played with test beam data
1520 Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1521 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1522 }
1523 else if(iCalib==11){// Jason played with test beam MC
1524 Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1525 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1526 }
1527
1528 return recalibfactor;
1529}
1530
1531
1532//________________________________________________________________________
1533Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1534{
1535 // Get maximum energy of attached cell.
1536
1537 id = -1;
1538 AliVCaloCells *fVCells=NULL;
1539 if(fEsdEv) fVCells = fEsdEv->GetEMCALCells();
1540 else if(fAodEv) fVCells = fAodEv->GetEMCALCells();
1541 if(!fVCells)
1542 return 0;
1543
1544 Double_t maxe = 0;
1545 Int_t ncells = cluster->GetNCells();
1546 for (Int_t i=0; i<ncells; i++) {
1547 Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1548 if (e>maxe) {
1549 maxe = e;
1550 id = cluster->GetCellAbsId(i);
1551 }
1552 }
1553 return maxe;
1554}
1555
1556
1557//________________________________________________________________________
1558Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){
1559
1560 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1561
1562 Int_t nPTracks= mcEvent->GetNumberOfPrimaries();
1563
1564 Int_t isPhysPrimary = 1;
1565 Int_t ismHF = 0;
1566 Int_t ismLongLivedOrK = 0;
1567
1568 if(mcP->GetMother()<0)//if it has no mother...
1569 return isPhysPrimary;
1570
1571 Int_t imTrack = mcP->GetMother();
1572 AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1573
1574 if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark...
1575 return isPhysPrimary;
1576
1577
1578 //############################################
1579 //get the PDG digits....
1580 int num = mcPm->PdgCode();
1581 int RevDigits[10] = {0};
1582 int nDigits = 0;
1583 while (num >= 1){
1584 RevDigits[nDigits++] = num%10;
1585 num = num / 10;
1586 }
1587 //##############################################
1588
1589
1590 if(RevDigits[3]>3)//Baryons
1591 ismHF = 1;
1592 else if(RevDigits[2]>3)//Mesons
1593 ismHF = 1;
1594
1595 ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1596
1597 if(!ismHF && ismLongLivedOrK)
1598 isPhysPrimary = 0;
1599 else{ // check grandmother, greatgrandmothers, etc...
1600 while(imTrack >= nPTracks){
1601
1602 if(mcPm->GetMother()<0)//if it has no mother...
1603 break;
1604
1605 if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark...
1606 return isPhysPrimary;
1607
1608 imTrack = mcPm->GetMother();
1609 mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1610
1611 //############################################
1612 //get the PDG digits....
1613 num = mcPm->PdgCode();
1614 for(int i=0; i<10; i++) RevDigits[i] = 0;
1615 nDigits = 0;
1616 while (num >= 1){
1617 RevDigits[nDigits++] = num%10;
1618 num = num / 10;
1619 }
1620 //##############################################
1621 if(RevDigits[3]>3)//Baryons
1622 ismHF = 1;
1623 else if(RevDigits[2]>3)//Mesons
1624 ismHF = 1;
1625
1626 ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1627
1628 if(!ismHF && ismLongLivedOrK)
1629 isPhysPrimary = 0;
1630
1631 }//while( >=nPTracks)
1632 }
1633
1634 return isPhysPrimary;
1635}
1636
1637
1638//________________________________________________________________________
1639Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::IsLongLivedOrK(Int_t MyPDGcode){
1640
1641 Int_t MyFlag = 0;
1642
1643 if(
1644 (TMath::Abs(MyPDGcode) == 22 ) || // Photon
1645 (TMath::Abs(MyPDGcode) == 11 ) || // Electron
1646 (TMath::Abs(MyPDGcode) == 13 ) || // Muon(-)
1647 (TMath::Abs(MyPDGcode) == 211 ) || // Pion
1648 (TMath::Abs(MyPDGcode) == 321 ) || // Kaon
1649 (TMath::Abs(MyPDGcode) == 310 ) || // K0s
1650 (TMath::Abs(MyPDGcode) == 130 ) || // K0l
1651 (TMath::Abs(MyPDGcode) == 2212) || // Proton
1652 (TMath::Abs(MyPDGcode) == 2112) || // Neutron
1653 (TMath::Abs(MyPDGcode) == 3122) || // Lambda_0
1654 (TMath::Abs(MyPDGcode) == 3112) || // Sigma Minus
1655 (TMath::Abs(MyPDGcode) == 3222) || // Sigma Plus
1656 (TMath::Abs(MyPDGcode) == 3312) || // Xsi Minus
1657 (TMath::Abs(MyPDGcode) == 3322) || // Xsi
1658 (TMath::Abs(MyPDGcode) == 3334) || // Omega
1659 (TMath::Abs(MyPDGcode) == 12 ) || // Electron Neutrino
1660 (TMath::Abs(MyPDGcode) == 14 ) || // Muon Neutrino
1661 (TMath::Abs(MyPDGcode) == 16 ) ) // Tau Neutrino
1662 MyFlag = 1;
1663
1664 return MyFlag;
1665}