Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALMesonGGSDM.cxx
CommitLineData
3a66195e 1#include "AliAnalysisTaskEMCALMesonGGSDM.h"
2
3#include <vector>
4#include <Riostream.h>
5#include <TChain.h>
6#include <TTree.h>
7#include <TF1.h>
8#include <TH1F.h>
9#include <TH2F.h>
10#include <TH3F.h>
11#include <TH1D.h>
12#include <TH2D.h>
13#include <TH3D.h>
14#include <TCanvas.h>
15#include <TList.h>
16#include <TFile.h>
17#include <TLorentzVector.h>
18#include <TNtuple.h>
19#include <TRandom3.h>
20
21#include "AliAnalysisTaskSE.h"
22#include "AliAnalysisManager.h"
23#include "AliStack.h"
24#include "AliESDtrackCuts.h"
25#include "AliESDEvent.h"
26#include "AliESDInputHandler.h"
27#include "AliAODEvent.h"
28#include "AliMCEvent.h"
29#include "AliEMCALGeometry.h"
30#include "AliInputEventHandler.h"
31#include "AliESDInputHandler.h"
32#include "AliAODInputHandler.h"
33
34#include "AliEMCALRecoUtils.h"
35#include "AliExternalTrackParam.h"
36
37// ROOT includes
38#include <TGeoManager.h>
39#include <TGeoMatrix.h>
40#include <TGeoBBox.h>
41#include <TH2F.h>
42#include <TArrayI.h>
43#include <TArrayF.h>
44#include <TObjArray.h>
45
46// STEER includes
47#include "AliVCluster.h"
48#include "AliVCaloCells.h"
49#include "AliLog.h"
50#include "AliPID.h"
51#include "AliESDEvent.h"
52#include "AliAODEvent.h"
53#include "AliESDtrack.h"
54#include "AliAODTrack.h"
55#include "AliExternalTrackParam.h"
56#include "AliESDfriendTrack.h"
57#include "AliTrackerBase.h"
58
59// EMCAL includes
60#include "AliEMCALRecoUtils.h"
61#include "AliEMCALGeometry.h"
62#include "AliTrackerBase.h"
63#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
64#include "AliEMCALPIDUtils.h"
65
66#include "AliGenCocktailEventHeader.h"
67
2942f542 68using std::cout;
69using std::endl;
70
3a66195e 71ClassImp(AliAnalysisTaskEMCALMesonGGSDM)
72
73//________________________________________________________________________
74AliAnalysisTaskEMCALMesonGGSDM::AliAnalysisTaskEMCALMesonGGSDM() :
75 AliAnalysisTaskSE(),
76 fOutput(0),
77 fMcMode(0),
78 fMyMCType(0),
79 fRecalibrator(0),
80 fdRmin_ClustTrack(0),
81 fPhimin(0),
82 fPhimax(0),
83 fEtamin(0),
84 fEtamax(0),
85 fTrackCuts(0),
86 fEsdEv(0),
87 fAodEv(0),
88 h1_nClusters(0),
89 h1_zvtx(0),
90 h1_trigger(0),
91 h1_M(0),
92 h1_M_mix(0),
93 h1_E(0),
94 h2_PhiEtaCluster(0),
95 h2_PhiEtaClusterCut(0),
96 h2_PhiEtaMaxCell(0),
97 h2_PhiEtaMaxCellCut(0),
98 h1_dR_ClustTrk(0),
99 h2_gE_RecTruth(0),
100 h2_eop_E(0),
101 h2_eop_pT(0),
102 h2_E_time(0),
103 h1_Pi0TruthPt(0),
104 h1_K0Pi0TruthPt(0),
105 h1_PriPi0TruthPt(0),
106 h1_PhysPi0TruthPt(0),
107 h1_Pi0TruthPtEmcal(0),
108 h1_K0Pi0TruthPtEmcal(0),
109 h1_PriPi0TruthPtEmcal(0),
110 h1_PhysPi0TruthPtEmcal(0),
111 h1_Pi0TruthPtPhi2piEta065(0),
112 h1_K0Pi0TruthPtPhi2piEta065(0),
113 h1_PriPi0TruthPtPhi2piEta065(0),
114 h1_PhysPi0TruthPtPhi2piEta065(0),
115 h1_Pi0TruthPtPhi2piEta1(0),
116 h1_K0Pi0TruthPtPhi2piEta1(0),
117 h1_PriPi0TruthPtPhi2piEta1(0),
118 h1_PhysPi0TruthPtPhi2piEta1(0),
119 h2_Pi0TruthPhiEta(0),
120 h2_PriPi0TruthPhiEta(0),
121 h2_Pi0TruthPhiEtaEmcal(0),
122 h2_PriPi0TruthPhiEtaEmcal(0),
123 h1_TruthPhotonsEmcal(0),
124 h2_TruthPhotonsPhiEta(0),
125 h1_PhotonsEmcal(0),
126 h1_PhotonsNCellsCut(0),
127 h1_PhotonsTrackMatchCut(0),
128 h1_PhotonsAllCut(0),
129 h2_PhotonsPhiEtaIsEmcal(0),
130 h1_dR_RealMC(0),
131 h2_Mpt_Pri(0),
132 h2_Mpt_Sec(0),
133 h3_MptR_Sec(0),
134 h3_MptR_K0s(0),
135 h3_MptR_Mat(0),
136 h2_PtR_MatM(0),
137 h2_Mpt_Pri_conv(0),
138 h2_Mpt_Sec_conv(0),
139 h3_MptR_Sec_conv(0),
140 h3_MptR_K0s_conv(0),
141 h3_MptR_Mat_conv(0),
142 h1_eConversionR(0),
143 h1_PriPi0Mother(0),
144 h1_SecPi0Mother(0),
145 h1_Chi2(0),
146 h1_nTrkMatch(0),
147 h1_nCells(0),
148 h1_ClusterDisp(0),
149 h2_Ellipse(0),
150 h2_EtaPt(0),
151//h2_Mpt(0),
152 h3_MptAsymm(0),
153//h2_Mpt_mix(0),
154 h3_MptAsymm_mix(0),
155 h2_dphi_deta(0),
156 h2_dphi_deta_mix(0),
157 h2_DispRes(0),
158 h2_cells_M02(0),
159 TriggerList(0),
160 fHelperClass(0)
161{
162 // Dummy constructor ALWAYS needed for I/O.
163}
164
165//________________________________________________________________________
166AliAnalysisTaskEMCALMesonGGSDM::AliAnalysisTaskEMCALMesonGGSDM(const char *name) :
167 AliAnalysisTaskSE(name),
168 fOutput(0),
169 fMcMode(0),
170 fMyMCType(0),
171 fRecalibrator(0),
172 fdRmin_ClustTrack(0),
173 fPhimin(0),
174 fPhimax(0),
175 fEtamin(0),
176 fEtamax(0),
177 fTrackCuts(0),
178 fEsdEv(0),
179 fAodEv(0),
180 h1_nClusters(0),
181 h1_zvtx(0),
182 h1_trigger(0),
183 h1_M(0),
184 h1_M_mix(0),
185 h1_E(0),
186 h2_PhiEtaCluster(0),
187 h2_PhiEtaClusterCut(0),
188 h2_PhiEtaMaxCell(0),
189 h2_PhiEtaMaxCellCut(0),
190 h1_dR_ClustTrk(0),
191 h2_gE_RecTruth(0),
192 h2_eop_E(0),
193 h2_eop_pT(0),
194 h2_E_time(0),
195 h1_Pi0TruthPt(0),
196 h1_K0Pi0TruthPt(0),
197 h1_PriPi0TruthPt(0),
198 h1_PhysPi0TruthPt(0),
199 h1_Pi0TruthPtEmcal(0),
200 h1_K0Pi0TruthPtEmcal(0),
201 h1_PriPi0TruthPtEmcal(0),
202 h1_PhysPi0TruthPtEmcal(0),
203 h1_Pi0TruthPtPhi2piEta065(0),
204 h1_K0Pi0TruthPtPhi2piEta065(0),
205 h1_PriPi0TruthPtPhi2piEta065(0),
206 h1_PhysPi0TruthPtPhi2piEta065(0),
207 h1_Pi0TruthPtPhi2piEta1(0),
208 h1_K0Pi0TruthPtPhi2piEta1(0),
209 h1_PriPi0TruthPtPhi2piEta1(0),
210 h1_PhysPi0TruthPtPhi2piEta1(0),
211 h2_Pi0TruthPhiEta(0),
212 h2_PriPi0TruthPhiEta(0),
213 h2_Pi0TruthPhiEtaEmcal(0),
214 h2_PriPi0TruthPhiEtaEmcal(0),
215 h1_TruthPhotonsEmcal(0),
216 h2_TruthPhotonsPhiEta(0),
217 h1_PhotonsEmcal(0),
218 h1_PhotonsNCellsCut(0),
219 h1_PhotonsTrackMatchCut(0),
220 h1_PhotonsAllCut(0),
221 h2_PhotonsPhiEtaIsEmcal(0),
222 h1_dR_RealMC(0),
223 h2_Mpt_Pri(0),
224 h2_Mpt_Sec(0),
225 h3_MptR_Sec(0),
226 h3_MptR_K0s(0),
227 h3_MptR_Mat(0),
228 h2_PtR_MatM(0),
229 h2_Mpt_Pri_conv(0),
230 h2_Mpt_Sec_conv(0),
231 h3_MptR_Sec_conv(0),
232 h3_MptR_K0s_conv(0),
233 h3_MptR_Mat_conv(0),
234 h1_eConversionR(0),
235 h1_PriPi0Mother(0),
236 h1_SecPi0Mother(0),
237 h1_Chi2(0),
238 h1_nTrkMatch(0),
239 h1_nCells(0),
240 h1_ClusterDisp(0),
241 h2_Ellipse(0),
242 h2_EtaPt(0),
243//h2_Mpt(0),
244 h3_MptAsymm(0),
245//h2_Mpt_mix(0),
246 h3_MptAsymm_mix(0),
247 h2_dphi_deta(0),
248 h2_dphi_deta_mix(0),
249 h2_DispRes(0),
250 h2_cells_M02(0),
251 TriggerList(0),
252 fHelperClass(0)
253{
254 // Constructor
255 // Define input and output slots here (never in the dummy constructor)
256 // Input slot #0 works with a TChain - it is connected to the default input container
257 // Output slot #1 writes into a TH1 container
258
259 DefineOutput(1, TList::Class()); // for output list
260}
261
262//________________________________________________________________________
263AliAnalysisTaskEMCALMesonGGSDM::~AliAnalysisTaskEMCALMesonGGSDM()
264{
265 // Destructor. Clean-up the output list, but not the histograms that are put inside
266 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
267 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
268 delete fOutput;
269 }
270 delete fTrackCuts;
271}
272
273//________________________________________________________________________
274void AliAnalysisTaskEMCALMesonGGSDM::UserCreateOutputObjects()
275{
276 // Create histograms
277 // Called once (on the worker node)
278
279 fOutput = new TList();
280 fOutput->SetOwner(); // IMPORTANT!
281
282 fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
283
284 cout << "__________AliAnalysisTaskEMCALMesonGGSDM: Input settings__________" << endl;
285 cout << " fMcMode: " << fMcMode << endl;
286 cout << " fMyMCType: " << fMyMCType << endl;
287 cout << " fRecalibrator: " << fRecalibrator << endl;
288 cout << " dRmin_ClustTrack: " << fdRmin_ClustTrack << endl;
289 cout << " phi range: " << fPhimin << ", " << fPhimax << endl;
290 cout << " eta range: " << fEtamin << ", " << fEtamax << endl;
291 cout << " number of zvtx bins: " << zvtx_bins << endl;
292 cout << " number of mult bins: " << mult_bins << endl;
293 cout << " poolDepth: " << poolDepth << endl;
294 cout << endl;
295
296
297 //AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
298 //AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
299 //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
300 //
301 //fPIDCombined = new AliPIDCombined();
302 //fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);
303 //fPIDCombined->SetDetectorMask(AliPIDResponse::kDetEMCAL);
304 //fPIDCombined->SetEnablePriors(kFALSE);
305
306 double TotalNBins = 0.0;
307
308 // Create histograms
309 Int_t nClustersbins = 501;
310 Float_t nClusterslow = -0.5, nClustersup = 500.5;
311 h1_nClusters = new TH1F("h1_nClusters", "# of clusters", nClustersbins, nClusterslow, nClustersup);
312 h1_nClusters->GetXaxis()->SetTitle("number of clusters/evt");
313 h1_nClusters->GetYaxis()->SetTitle("counts");
314 h1_nClusters->SetMarkerStyle(kFullCircle);
315 TotalNBins+=nClustersbins;
316
317 Int_t nZvertexbins = 501;
318 Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
319 h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
320 h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
321 h1_zvtx->GetYaxis()->SetTitle("counts");
322 h1_zvtx->SetMarkerStyle(kFullCircle);
323 TotalNBins+=nZvertexbins;
324
325 h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
326 TotalNBins+=1001;
327
328 Int_t Mbins = 3000;
329 Float_t Mlow = 0.0, Mup = 3.0;
330 h1_M = new TH1F("h1_M", "Invariant Mass", Mbins, Mlow, Mup);
331 h1_M->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
332 h1_M->GetYaxis()->SetTitle("counts");
333 h1_M->SetMarkerStyle(kFullCircle);
334 TotalNBins+=Mbins;
335
336 h1_M_mix = new TH1F("h1_M_mix", "Invariant Mass (mixed events)", Mbins, Mlow, Mup);
337 h1_M_mix->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
338 h1_M_mix->GetYaxis()->SetTitle("counts");
339 h1_M_mix->SetMarkerStyle(kFullCircle);
340 TotalNBins+=Mbins;
341
342 Int_t ptbins = 2000;
343 Float_t ptlow = 0.0, ptup = 20.0;
344 Int_t Ebins = 1000;
345 Float_t Elow = 0.0, Eup = 20.0;
346 h1_E = new TH1F("h1_E", "Cluster Energy in EMCal", Ebins, Elow, Eup);
347 h1_E->GetXaxis()->SetTitle("E [GeV]");
348 h1_E->GetYaxis()->SetTitle("counts");
349 h1_E->SetMarkerStyle(kFullCircle);
350 TotalNBins+=Ebins;
351
352 h2_PhiEtaCluster = new TH2F("h2_PhiEtaCluster", "cluster phi vs eta", 400,1.362,3.178, 300,-0.728,0.728);
353 h2_PhiEtaCluster->GetXaxis()->SetTitle("#phi [rad]");
354 h2_PhiEtaCluster->GetYaxis()->SetTitle("#eta");
355 h2_PhiEtaCluster->GetZaxis()->SetTitle("hits");
356 h2_PhiEtaCluster->SetMarkerStyle(kFullCircle);
357 TotalNBins+=400*300;
358
359 h2_PhiEtaClusterCut = new TH2F("h2_PhiEtaClusterCut", "cluster phi vs eta (after cuts)", 400,1.362,3.178, 300,-0.728,0.728);
360 h2_PhiEtaClusterCut->GetXaxis()->SetTitle("#phi [rad]");
361 h2_PhiEtaClusterCut->GetYaxis()->SetTitle("#eta");
362 h2_PhiEtaClusterCut->GetZaxis()->SetTitle("hits");
363 h2_PhiEtaClusterCut->SetMarkerStyle(kFullCircle);
364 TotalNBins+=400*300;
365
366// eta binning
367 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};
368
369 // phi binning
370 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};
371
372 h2_PhiEtaMaxCell = new TH2F("h2_PhiEtaMaxCell", "maxcell phi vs eta", 124,PhiBins, 96,EtaBins);
373 h2_PhiEtaMaxCell->GetXaxis()->SetTitle("#phi [rad]");
374 h2_PhiEtaMaxCell->GetYaxis()->SetTitle("#eta");
375 h2_PhiEtaMaxCell->GetZaxis()->SetTitle("hits");
376 h2_PhiEtaMaxCell->SetMarkerStyle(kFullCircle);
377 TotalNBins+=96*124;
378
379 h2_PhiEtaMaxCellCut = new TH2F("h2_PhiEtaMaxCellCut", "maxcell phi vs eta (after cuts)", 124,PhiBins, 96,EtaBins);
380 h2_PhiEtaMaxCellCut->GetXaxis()->SetTitle("#phi [rad]");
381 h2_PhiEtaMaxCellCut->GetYaxis()->SetTitle("#eta");
382 h2_PhiEtaMaxCellCut->GetZaxis()->SetTitle("hits");
383 h2_PhiEtaMaxCellCut->SetMarkerStyle(kFullCircle);
384 TotalNBins+=96*124;
385
386 h1_dR_ClustTrk = new TH1F("h1_dR_ClustTrk", "Cluster-Track matching", 5000, -0.01, 5);
387 h1_dR_ClustTrk->GetXaxis()->SetTitle("dR [sqrt(d#phi^{2}+d#eta^{2})]");
388 h1_dR_ClustTrk->GetYaxis()->SetTitle("N");
389 h1_dR_ClustTrk->SetMarkerStyle(kFullCircle);
390 TotalNBins+=5000;
391
392 h2_gE_RecTruth = new TH2F("h2_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 500,0,2);
393 h2_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
394 h2_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
395 h2_gE_RecTruth->GetZaxis()->SetTitle("counts");
396 h2_gE_RecTruth->SetMarkerStyle(kFullCircle);
397 TotalNBins+=Ebins*500;
398
399 h2_eop_E = new TH2F("h2_eop_E","E/p vs E (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
400 h2_eop_E->GetXaxis()->SetTitle("cluster Energy [GeV]");
401 h2_eop_E->GetYaxis()->SetTitle("E/p");
402 TotalNBins+=Ebins*1200;
403
404 h2_eop_pT = new TH2F("h2_eop_pT","E/p vs p_{T} (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
405 h2_eop_pT->GetXaxis()->SetTitle("cluster Energy [GeV]");
406 h2_eop_pT->GetYaxis()->SetTitle("E/p");
407 TotalNBins+=Ebins*1200;
408
409 h2_E_time = new TH2F("h2_E_time","cluster energy vs time", Ebins, Elow, Eup, 1000,-1e-6,1e-6);
410 h2_E_time->GetXaxis()->SetTitle("cluster Energy [GeV]");
411 h2_E_time->GetYaxis()->SetTitle("time [s]");
412 TotalNBins+=Ebins*1000;
413
414 h1_Pi0TruthPt = new TH1F("h1_Pi0TruthPt", "P_{T} distribution for Truth Pi0's", ptbins, ptlow, ptup);
415 h1_Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
416 h1_Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
417 h1_Pi0TruthPt->SetMarkerStyle(kFullCircle);
418 TotalNBins+=ptbins;
419
420 h1_K0Pi0TruthPt = new TH1F("h1_K0Pi0TruthPt", "P_{T} distribution for Truth Pi0's from K^{0}_{s} decays", ptbins, ptlow, ptup);
421 h1_K0Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
422 h1_K0Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
423 h1_K0Pi0TruthPt->SetMarkerStyle(kFullCircle);
424 TotalNBins+=ptbins;
425
426 h1_PriPi0TruthPt = new TH1F("h1_PriPi0TruthPt", "P_{T} distribution for Truth Primary Pi0's", ptbins, ptlow, ptup);
427 h1_PriPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
428 h1_PriPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
429 h1_PriPi0TruthPt->SetMarkerStyle(kFullCircle);
430 TotalNBins+=ptbins;
431
432 h1_PhysPi0TruthPt = new TH1F("h1_PhysPi0TruthPt", "P_{T} distribution for Truth Physical Primary Pi0's", ptbins, ptlow, ptup);
433 h1_PhysPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
434 h1_PhysPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
435 h1_PhysPi0TruthPt->SetMarkerStyle(kFullCircle);
436 TotalNBins+=ptbins;
437
438 h1_Pi0TruthPtEmcal = new TH1F("h1_Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's (hit EMCal)", ptbins, ptlow, ptup);
439 h1_Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
440 h1_Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
441 h1_Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
442 TotalNBins+=ptbins;
443
444 h1_K0Pi0TruthPtEmcal = new TH1F("h1_K0Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's from K^{0}_{s} decays (hit EMCal)", ptbins, ptlow, ptup);
445 h1_K0Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
446 h1_K0Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
447 h1_K0Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
448 TotalNBins+=ptbins;
449
450 h1_PriPi0TruthPtEmcal = new TH1F("h1_PriPi0TruthPtEmcal", "P_{T} distribution for Truth Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
451 h1_PriPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
452 h1_PriPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
453 h1_PriPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
454 TotalNBins+=ptbins;
455
456 h1_PhysPi0TruthPtEmcal = new TH1F("h1_PhysPi0TruthPtEmcal", "P_{T} distribution for Truth Physical Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
457 h1_PhysPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
458 h1_PhysPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
459 h1_PhysPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
460 TotalNBins+=ptbins;
461
462 h1_Pi0TruthPtPhi2piEta065 = new TH1F("h1_Pi0TruthPtPhi2piEta065",
463 "P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
464 h1_Pi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
465 h1_Pi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
466 h1_Pi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
467 TotalNBins+=ptbins;
468
469 h1_K0Pi0TruthPtPhi2piEta065 = new TH1F("h1_K0Pi0TruthPtPhi2piEta065",
470 "P_{T} for Truth Pi0's (K^{0}_{s} decays) [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
471 h1_K0Pi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
472 h1_K0Pi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
473 h1_K0Pi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
474 TotalNBins+=ptbins;
475
476 h1_PriPi0TruthPtPhi2piEta065 = new TH1F("h1_PriPi0TruthPtPhi2piEta065",
477 "P_{T} for Primary Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
478 h1_PriPi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
479 h1_PriPi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
480 h1_PriPi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
481 TotalNBins+=ptbins;
482
483 h1_PhysPi0TruthPtPhi2piEta065 = new TH1F("h1_PhysPi0TruthPtPhi2piEta065",
484 "P_{T} for Truth Pi0's (not from material) [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
485 h1_PhysPi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
486 h1_PhysPi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
487 h1_PhysPi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
488 TotalNBins+=ptbins;
489
490 h1_Pi0TruthPtPhi2piEta1 = new TH1F("h1_Pi0TruthPtPhi2piEta1",
491 "P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
492 h1_Pi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
493 h1_Pi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
494 h1_Pi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
495 TotalNBins+=ptbins;
496
497 h1_K0Pi0TruthPtPhi2piEta1 = new TH1F("h1_K0Pi0TruthPtPhi2piEta1",
498 "P_{T} for Truth Pi0's (k^{0}_{s} decays) [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
499 h1_K0Pi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
500 h1_K0Pi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
501 h1_K0Pi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
502 TotalNBins+=ptbins;
503
504 h1_PriPi0TruthPtPhi2piEta1 = new TH1F("h1_PriPi0TruthPtPhi2piEta1",
505 "P_{T} for Primary Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
506 h1_PriPi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
507 h1_PriPi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
508 h1_PriPi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
509 TotalNBins+=ptbins;
510
511 h1_PhysPi0TruthPtPhi2piEta1 = new TH1F("h1_PhysPi0TruthPtPhi2piEta1",
512 "P_{T} for Truth Pi0's (not from material) [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
513 h1_PhysPi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
514 h1_PhysPi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
515 h1_PhysPi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
516 TotalNBins+=ptbins;
517
518 h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
519 h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
520 h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
521 TotalNBins+=380*200;
522
523 h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
524 h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
525 h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
526 TotalNBins+=380*200;
527
528 h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
529 h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
530 h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
531 TotalNBins+=380*150;
532
533 h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
534 h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
535 h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
536 TotalNBins+=380*150;
537
538 h1_TruthPhotonsEmcal = new TH1F("h1_TruthPhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
539 h1_TruthPhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
540 h1_TruthPhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
541 h1_TruthPhotonsEmcal->SetMarkerStyle(kFullCircle);
542 TotalNBins+=ptbins;
543
544 h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
545 "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
546 h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
547 h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
548 h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
549 TotalNBins+=380*150;
550
551 h1_PhotonsEmcal = new TH1F("h1_PhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
552 h1_PhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
553 h1_PhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
554 h1_PhotonsEmcal->SetMarkerStyle(kFullCircle);
555 TotalNBins+=ptbins;
556
557 h1_PhotonsNCellsCut = new TH1F("h1_PhotonsNCellsCut", "P_{T} distribution for #gamma's that survive NCells cut", ptbins, ptlow, ptup);
558 h1_PhotonsNCellsCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
559 h1_PhotonsNCellsCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
560 h1_PhotonsNCellsCut->SetMarkerStyle(kFullCircle);
561 TotalNBins+=ptbins;
562
563 h1_PhotonsTrackMatchCut = new TH1F("h1_PhotonsTrackMatchCut", "P_{T} distribution for #gamma's that survive TrackMatch cut", ptbins, ptlow, ptup);
564 h1_PhotonsTrackMatchCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
565 h1_PhotonsTrackMatchCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
566 h1_PhotonsTrackMatchCut->SetMarkerStyle(kFullCircle);
567 TotalNBins+=ptbins;
568
569 h1_PhotonsAllCut = new TH1F("h1_PhotonsAllCut", "P_{T} distribution for #gamma's that survive All cut", ptbins, ptlow, ptup);
570 h1_PhotonsAllCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
571 h1_PhotonsAllCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
572 h1_PhotonsAllCut->SetMarkerStyle(kFullCircle);
573 TotalNBins+=ptbins;
574
575 h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
576 "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 100,-1.0,1.0);
577 h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
578 h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
579 TotalNBins+=380*100;
580
581 h1_dR_RealMC = new TH1F("h1_dR_RealMC", "P_{T} distribution for #gamma's that survive All cut", 2000, -0.01, 10);
582 h1_dR_RealMC->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
583 h1_dR_RealMC->GetYaxis()->SetTitle("N");
584 h1_dR_RealMC->SetMarkerStyle(kFullCircle);
585 TotalNBins+=2000;
586
587 h2_Mpt_Pri = new TH2F("h2_Mpt_Pri", "mass vs pT for primary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
588 h2_Mpt_Pri->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
589 h2_Mpt_Pri->GetYaxis()->SetTitle("p_{T} [GeV/c]");
590 h2_Mpt_Pri->SetMarkerStyle(kFullCircle);
591 TotalNBins+=Mbins*ptbins;
592
593 h2_Mpt_Sec = new TH2F("h2_Mpt_Sec", "mass vs pT for secondary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
594 h2_Mpt_Sec->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
595 h2_Mpt_Sec->GetYaxis()->SetTitle("p_{T} [GeV/c]");
596 h2_Mpt_Sec->SetMarkerStyle(kFullCircle);
597 TotalNBins+=Mbins*ptbins;
598
599 h3_MptR_Sec = new TH3F("h3_MptR_Sec", "mass vs pT vs production radius for secondary pions", 500,0,0.5, 100,0,20, 300,0,600);
600 h3_MptR_Sec->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
601 h3_MptR_Sec->GetYaxis()->SetTitle("p_{T} [GeV/c]");
602 h3_MptR_Sec->GetZaxis()->SetTitle("production radius [cm]");
603 h3_MptR_Sec->SetMarkerStyle(kFullCircle);
604 TotalNBins+=500*100*300;
605
606 h3_MptR_K0s = new TH3F("h3_MptR_K0s", "mass vs pT vs production radius for K0s pions", 500,0,0.5, 100,0,20, 300,0,600);
607 h3_MptR_K0s->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
608 h3_MptR_K0s->GetYaxis()->SetTitle("p_{T} [GeV/c]");
609 h3_MptR_K0s->GetZaxis()->SetTitle("production radius [cm]");
610 h3_MptR_K0s->SetMarkerStyle(kFullCircle);
611 TotalNBins+=500*100*300;
612
613 h3_MptR_Mat = new TH3F("h3_MptR_Mat", "mass vs pT vs production radius for material pions", 500,0,0.5, 100,0,20, 300,0,600);
614 h3_MptR_Mat->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
615 h3_MptR_Mat->GetYaxis()->SetTitle("p_{T} [GeV/c]");
616 h3_MptR_Mat->GetZaxis()->SetTitle("production radius [cm]");
617 h3_MptR_Mat->SetMarkerStyle(kFullCircle);
618 TotalNBins+=500*100*300;
619
620 h2_PtR_MatM = new TH2F("h2_PtR_MatM", "pT vs production radius for merged material pions (pi mass assumed)", 100,0,20, 300,0,600);
621 h2_PtR_MatM->GetXaxis()->SetTitle("p_{T} [GeV/c]");
622 h2_PtR_MatM->GetYaxis()->SetTitle("production radius [cm]");
623 h2_PtR_MatM->SetMarkerStyle(kFullCircle);
624 TotalNBins+=100*300;
625
626 h2_Mpt_Pri_conv = new TH2F("h2_Mpt_Pri_conv", "mass vs pT for primary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
627 h2_Mpt_Pri_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
628 h2_Mpt_Pri_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
629 h2_Mpt_Pri_conv->SetMarkerStyle(kFullCircle);
630 TotalNBins+=Mbins*ptbins;
631
632 h2_Mpt_Sec_conv = new TH2F("h2_Mpt_Sec_conv", "mass vs pT for secondary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
633 h2_Mpt_Sec_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
634 h2_Mpt_Sec_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
635 h2_Mpt_Sec_conv->SetMarkerStyle(kFullCircle);
636 TotalNBins+=Mbins*ptbins;
637
638 h3_MptR_Sec_conv = new TH3F("h3_MptR_Sec_conv", "mass vs pT vs production radius for secondary pions", 500,0,0.5, 100,0,20, 300,0,600);
639 h3_MptR_Sec_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
640 h3_MptR_Sec_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
641 h3_MptR_Sec_conv->GetZaxis()->SetTitle("production radius [cm]");
642 h3_MptR_Sec_conv->SetMarkerStyle(kFullCircle);
643 TotalNBins+=500*100*300;
644
645 h3_MptR_K0s_conv = new TH3F("h3_MptR_K0s_conv", "mass vs pT vs production radius for K0s pions", 500,0,0.5, 100,0,20, 300,0,600);
646 h3_MptR_K0s_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
647 h3_MptR_K0s_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
648 h3_MptR_K0s_conv->GetZaxis()->SetTitle("production radius [cm]");
649 h3_MptR_K0s_conv->SetMarkerStyle(kFullCircle);
650 TotalNBins+=500*100*300;
651
652 h3_MptR_Mat_conv = new TH3F("h3_MptR_Mat_conv", "mass vs pT vs production radius for material pions", 500,0,0.5, 100,0,20, 300,0,600);
653 h3_MptR_Mat_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
654 h3_MptR_Mat_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
655 h3_MptR_Mat_conv->GetZaxis()->SetTitle("production radius [cm]");
656 h3_MptR_Mat_conv->SetMarkerStyle(kFullCircle);
657 TotalNBins+=500*100*300;
658
659 h1_eConversionR = new TH1F("h1_eConversionR", "conversion point (radius)", 600,0,600);
660 h1_eConversionR->GetXaxis()->SetTitle("production radius [cm]");
661 h1_eConversionR->SetMarkerStyle(kFullCircle);
662 TotalNBins+=600;
663
664 h1_PriPi0Mother = new TH1F("h1_PriPi0Mother", "primary pi0 mother ID", 12001,-6000.5,6000.5);
665 h1_PriPi0Mother->GetXaxis()->SetTitle("#pi^{0} mother ID");
666 h1_PriPi0Mother->SetMarkerStyle(kFullCircle);
667 TotalNBins+=12001;
668
669 h1_SecPi0Mother = new TH1F("h1_SecPi0Mother", "secondray pi0 mother ID", 12001,-6000.5,6000.5);
670 h1_SecPi0Mother->GetXaxis()->SetTitle("#pi^{0} mother ID");
671 h1_SecPi0Mother->SetMarkerStyle(kFullCircle);
672 TotalNBins+=12001;
673
674 Int_t chi2bins = 100;
675 Float_t chi2low = -2, chi2up = 2;
676 h1_Chi2 = new TH1F("h1_Chi2","#chi^{2} distribution for reconstructed",chi2bins, chi2low, chi2up);
677 h1_Chi2->GetXaxis()->SetTitle("#chi^{2}");
678 h1_Chi2->GetYaxis()->SetTitle("counts");
679 TotalNBins+=chi2bins;
680
681 h1_nTrkMatch = new TH1F("h1_nTrkMatch","number of matched tracks",14, -1.5, 5.5);
682 h1_nTrkMatch->GetXaxis()->SetTitle("nTracksMatched");
683 h1_nTrkMatch->GetYaxis()->SetTitle("counts");
684 TotalNBins+=14;
685
686 h1_ClusterDisp = new TH1F("h1_ClusterDisp","Dispersion of CaloCluster",1000, -1, 3);
687 h1_ClusterDisp->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
688 h1_ClusterDisp->GetYaxis()->SetTitle("counts");
689 TotalNBins+=1000;
690
691 h2_Ellipse = new TH2F("h2_Ellipse","Ellipse axis M20 vs M02",500, -0.01, 1, 500, -0.01, 1);
692 h2_Ellipse->GetXaxis()->SetTitle("cluster->GetM20()");
693 h2_Ellipse->GetYaxis()->SetTitle("cluster->GetM02()");
694 h2_Ellipse->GetZaxis()->SetTitle("counts");
695 TotalNBins+=500*500;
696
697 Int_t etabins = 150;
698 Float_t etalow = -1.5, etaup = 1.5;
699 h2_EtaPt = new TH2F("h2_EtaPt","Cluster Energy vs ",etabins, etalow, etaup, ptbins, ptlow, ptup);
700 h2_EtaPt->GetXaxis()->SetTitle("E [GeV]");
701 h2_EtaPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
702 TotalNBins+=etabins*ptbins;
703
704 h3_MptAsymm = new TH3F("h3_MptAsymm","mass vs p_{T} vs Asymm cut",Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
705 h3_MptAsymm->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
706 h3_MptAsymm->GetYaxis()->SetTitle("p_{T} [GeV/c]");
707 h3_MptAsymm->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
708 TotalNBins+=Mbins*ptbins*3.0;
709
710 h3_MptAsymm_mix = new TH3F("h3_MptAsymm_mix","mass vs p_{T} vs Asymm cut (mixed events)",Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
711 h3_MptAsymm_mix->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
712 h3_MptAsymm_mix->GetYaxis()->SetTitle("p_{T} [GeV/c]");
713 h3_MptAsymm_mix->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
714 TotalNBins+=Mbins*ptbins*3.0;
715
716 h2_dphi_deta = new TH2F("h2_dphi_deta","#Delta#phi vs #Delta#eta", 349,-1.5,5, 400,-2.0,2.0);
717 h2_dphi_deta->GetXaxis()->SetTitle("#Delta#phi");
718 h2_dphi_deta->GetYaxis()->SetTitle("#Delta#eta");
719 TotalNBins+=349*400;
720
721 h2_dphi_deta_mix = new TH2F("h2_dphi_deta_mix","#Delta#phi vs #Delta#eta (mixed events)", 349,-1.5,5, 400,-2.0,2.0);
722 h2_dphi_deta_mix->GetXaxis()->SetTitle("#Delta#phi");
723 h2_dphi_deta_mix->GetYaxis()->SetTitle("#Delta#eta");
724 TotalNBins+=349*400;
725
726 h2_DispRes = new TH2F("h2_DispRes", "zvtx info", 500,-0.01,1, 500,-0.1,2);
727 h2_DispRes->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
728 h2_DispRes->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
729 h2_DispRes->GetZaxis()->SetTitle("counts");
730 TotalNBins+=500*500;
731
732 h2_cells_M02 = new TH2F("h2_cells_M02", "nCells vs M02", 204,-1.5,100.5, 500,-1,1.5);
733 h2_cells_M02->GetXaxis()->SetTitle("nCells");
734 h2_cells_M02->GetYaxis()->SetTitle("M02");
735 h2_cells_M02->GetZaxis()->SetTitle("counts");
736 TotalNBins+=204*500;
737
738 cout << endl << "Total number of bins in booked histograms: " << TotalNBins << endl << endl;
739
740 // Initialize helper class (for vertex selection & pile up correction)
741 fHelperClass = new AliAnalysisUtils();
742
743 //TFile *f = OpenFile(1);
744 //TDirectory::TContext context(f);
745
746 fOutput->Add(h1_nClusters);
747 fOutput->Add(h1_zvtx);
748 fOutput->Add(h1_trigger);
749 fOutput->Add(h1_M);
750 fOutput->Add(h1_M_mix);
751 fOutput->Add(h1_E);
752 fOutput->Add(h2_PhiEtaCluster);
753 fOutput->Add(h2_PhiEtaClusterCut);
754 fOutput->Add(h2_PhiEtaMaxCell);
755 fOutput->Add(h2_PhiEtaMaxCellCut);
756 fOutput->Add(h1_dR_ClustTrk);
757 fOutput->Add(h2_gE_RecTruth);
758 fOutput->Add(h2_eop_E);
759 fOutput->Add(h2_eop_pT);
760 fOutput->Add(h2_E_time);
761 fOutput->Add(h1_Pi0TruthPt);
762 fOutput->Add(h1_K0Pi0TruthPt);
763 fOutput->Add(h1_PriPi0TruthPt);
764 fOutput->Add(h1_PhysPi0TruthPt);
765 fOutput->Add(h1_Pi0TruthPtEmcal);
766 fOutput->Add(h1_K0Pi0TruthPtEmcal);
767 fOutput->Add(h1_PriPi0TruthPtEmcal);
768 fOutput->Add(h1_PhysPi0TruthPtEmcal);
769 fOutput->Add(h1_Pi0TruthPtPhi2piEta065);
770 fOutput->Add(h1_K0Pi0TruthPtPhi2piEta065);
771 fOutput->Add(h1_PriPi0TruthPtPhi2piEta065);
772 fOutput->Add(h1_PhysPi0TruthPtPhi2piEta065);
773 fOutput->Add(h1_Pi0TruthPtPhi2piEta1);
774 fOutput->Add(h1_K0Pi0TruthPtPhi2piEta1);
775 fOutput->Add(h1_PriPi0TruthPtPhi2piEta1);
776 fOutput->Add(h1_PhysPi0TruthPtPhi2piEta1);
777 fOutput->Add(h2_Pi0TruthPhiEta);
778 fOutput->Add(h2_PriPi0TruthPhiEta);
779 fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
780 fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
781 fOutput->Add(h1_TruthPhotonsEmcal);
782 fOutput->Add(h2_TruthPhotonsPhiEta);
783 fOutput->Add(h1_PhotonsEmcal);
784 fOutput->Add(h1_PhotonsNCellsCut);
785 fOutput->Add(h1_PhotonsTrackMatchCut);
786 fOutput->Add(h1_PhotonsAllCut);
787 fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
788 fOutput->Add(h1_dR_RealMC);
789 fOutput->Add(h2_Mpt_Pri);
790 fOutput->Add(h2_Mpt_Sec);
791 fOutput->Add(h3_MptR_Sec);
792 fOutput->Add(h3_MptR_K0s);
793 fOutput->Add(h3_MptR_Mat);
794 fOutput->Add(h2_PtR_MatM);
795 fOutput->Add(h2_Mpt_Pri_conv);
796 fOutput->Add(h2_Mpt_Sec_conv);
797 fOutput->Add(h3_MptR_Sec_conv);
798 fOutput->Add(h3_MptR_K0s_conv);
799 fOutput->Add(h3_MptR_Mat_conv);
800 fOutput->Add(h1_eConversionR);
801 fOutput->Add(h1_PriPi0Mother);
802 fOutput->Add(h1_SecPi0Mother);
803 fOutput->Add(h1_Chi2);
804 fOutput->Add(h1_nTrkMatch);
805 fOutput->Add(h1_ClusterDisp);
806 fOutput->Add(h2_Ellipse);
807 fOutput->Add(h2_EtaPt);
808 fOutput->Add(h3_MptAsymm);
809 fOutput->Add(h3_MptAsymm_mix);
810 fOutput->Add(h2_dphi_deta);
811 fOutput->Add(h2_dphi_deta_mix);
812 fOutput->Add(h2_DispRes);
813 fOutput->Add(h2_cells_M02);
814
815 // Post data for ALL output slots >0 here,
816 // To get at least an empty histogram
817 // 1 is the outputnumber of a certain weg of task 1
818 PostData(1, fOutput);
819}
820
821//________________________________________________________________________
822void AliAnalysisTaskEMCALMesonGGSDM::UserExec(Option_t *)
823{
824 // Main loop Called for each event
825
826 AliMCEvent *mcEvent = MCEvent();
827 Bool_t isMC = bool(mcEvent);//is this the right way to do this?
828
829 TRandom3 randy; randy.SetSeed(0);
830 unsigned int iskip = -1;
831 TLorentzVector ParentMix;
832
833 double recalScale = 1.0;
834
835 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
836
837 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
838 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
839 if (!aodH && !esdH) Printf("ERROR: Could not get ESD or AODInputHandler");
840
841 if(esdH) fEsdEv = esdH->GetEvent();
842 else if(aodH) fAodEv = aodH->GetEvent();
843 else{
844 AliFatal("Neither ESD nor AOD event found");
845 return;
846 }
847
848 // get pointer to reconstructed event
849 AliVEvent *event = InputEvent();
850 if (!event){
851 AliError("Pointer == 0, this can not happen!"); return;}
852 //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
853 //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
854 //if (!fEsdEv){
855 //AliError("Cannot get the ESD event"); return;}
856
857 /*
858 fHelperClass->SetCutOnZVertexSPD(kFALSE);//does the zvtx have to match the spd vertex?
859 fHelperClass->SetMaxVtxZ(1.0e6);//i set this myself later..
860 // simply makes sure that there is at least 1 contributer to the zvtx determination.
861 // this should only remove the *extra* events at zvtx==0.
862 if(!fHelperClass->IsVertexSelected2013pA(event))
863 return;
864 */
865
866 Int_t iTrigger = 0;
867 if (fEsdEv) iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
868 else if (fAodEv) iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
869
870 char saythis[500];
871 Int_t iTriggerBin = 0;
872 for(unsigned long j=0; j<TriggerList.size(); j++){
873 if(iTrigger==TriggerList[j])
874 iTriggerBin=j+1;
875 }
876 if(iTriggerBin==0){
877 TriggerList.push_back(iTrigger);
878 iTriggerBin=TriggerList.size();
879 }
880
881 h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
882 sprintf(saythis,"%d",iTrigger);
883 h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
884
885 if(fEsdEv){
886 TString trigClasses = fEsdEv->GetFiredTriggerClasses();
887 // remove "fast cluster events":
888 if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL"))
889 return;
890 }
891 else if(fAodEv){
892 TString trigClasses = fAodEv->GetFiredTriggerClasses();
893 // remove "fast cluster events":
894 if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL"))
895 return;
896 }
897
898 if (fEsdEv){
899 if(!(fEsdEv->GetPrimaryVertex()->GetStatus())) return;
900 }
901 //else if (fAodEv){
902 //if(!(fAodEv->GetPrimaryVertex()->GetStatus())) return;
903 //}
904
905 Double_t vertDisp=0.0;
906 Double_t vertZres=0.0;
907 Bool_t vertIsfromZ=0;
908 if (fEsdEv){
909 vertDisp = fEsdEv->GetPrimaryVertex()->GetDispersion();
910 vertZres = fEsdEv->GetPrimaryVertex()->GetZRes();
911 vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
912 }
913 else if (fAodEv){
914 vertDisp = 0;
915 vertZres = 0;
916 vertIsfromZ = 0;
917 }
918
919 h2_DispRes->Fill(vertDisp, vertZres);
920 // if vertex is from spd vertexZ, require more stringent cut
921 if (vertIsfromZ) {
922 if (vertDisp>0.02 || vertZres>0.25 )
923 return; // bad vertex from VertexerZ
924 }
925
926
927 Int_t nclusters=0;
928 if(fEsdEv){
929 //Int_t evtN = fEsdEv->GetEventNumberInFile();
930 //Int_t ntracks = fEsdEv->GetNumberOfTracks();
931 nclusters = fEsdEv->GetNumberOfCaloClusters();
932 }
933 else if(fAodEv){
934 //Int_t evtN = fAodEv->GetEventNumberInFile();
935 //Int_t ntracks = fAodEv->GetNumberOfTracks();
936 nclusters = fAodEv->GetNumberOfCaloClusters();
937 }
938
939 // EMCal cluster loop for reconstructed event
940 //numberofclusters set above!
941 TLorentzVector Photon1, Photon2, Parent;
942 Double_t vertex[3];
943 Double_t E1=0.0;
944 Double_t E2=0.0;
945 Double_t vertZ=0.0;
946 if (fEsdEv) vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
947 else if (fAodEv) vertZ = fAodEv->GetPrimaryVertex()->GetZ();
948
949 h1_zvtx->Fill(vertZ);
950 //zvertex cut:
951 if(fabs(vertZ)>10.0)
952 return;
953
954 h1_nClusters->Fill(nclusters);
955
956 int izvtx = GetZvtxBin(vertZ);
957 int imult = GetMultBin(nclusters);
958
959 //cout << iskip << " " << izvtx << " " << imult << endl;
960 //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
961
962
963
964 //######################### ~~~~~~~~~~~ ##################################
965 //######################### STARTING MC ##################################
966 //######################### ~~~~~~~~~~~ ##################################
967
968 if(isMC){
969 int isPrimary = 0;
970 int isMaterialSec = 0;
971 int isK0sDecay = 0;
972
973 if (!mcEvent){
974 cout << "no MC event" << endl;
975 return;
976 }
977
978 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
979 if (!evtVtx)
980 return;
981
982 mcEvent->PreReadAll();
983
984 Int_t nTracksMC = mcEvent->GetNumberOfTracks();
985 Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
986 Int_t MyEMCPion = -1;
987 for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
988 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
989 if (!mcP)
990 continue;
991
992 if(iTrack<nPTracksMC) isPrimary = 1;
993 else isPrimary = 0;
994
995 isK0sDecay = 0;
996 if(mcP->GetMother()>-1){
997 if( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 310 ||
998 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -310 )
999 isK0sDecay = 1;
1000 }
1001
1002 // it's a pion !!
1003 if(mcP->PdgCode() != 111)
1004 continue;
1005
1006 MyEMCPion = 0;
1007 if(strcmp(fMyMCType,"ANY")==0)
1008 MyEMCPion = 1;
1009 else if(IsMyMCHeaderType(iTrack, fMyMCType, mcEvent)){
1010 //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << " i: " << iTrack << " pdg: " << mcP->PdgCode() << " pT: " << mcP->Pt() << endl;
1011 //cout << "iTrack: " << iTrack << " nPrimaryMC: " << nPTracksMC << " nTracksMC: " << nTracksMC << endl;
1012 MyEMCPion = 1;
1013 }
1014
1015 //if(MyEMCPion)
1016 //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << " i: " << iTrack << " pdg: " << mcP->PdgCode() << " pT: " << mcP->Pt() << endl;
1017 //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << " i: " << iTrack << " pdg: " << static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack+1))->PdgCode() << " pT: " << mcP->Pt() << endl;
1018
1019 if(MyEMCPion!=1 && isPrimary==1)
1020 continue;
1021
1022
1023 if(isPrimary==1 && mcP->GetMother()>-1)
1024 h1_PriPi0Mother->Fill( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() );
1025 else if(isPrimary==0 && mcP->GetMother()>-1)
1026 h1_SecPi0Mother->Fill( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() );
1027
1028 Int_t daughter[2] = {-1,-1};
1029 daughter[0] = mcP->GetFirstDaughter();
1030 daughter[1] = mcP->GetLastDaughter();
1031
1032 if (daughter[0]<0) continue;
1033 if (daughter[1]<0) daughter[1]=daughter[0];
1034 if (daughter[1]-daughter[0] != 1) continue;
1035
1036 Int_t eIndexofConvertedPhoton[2] = {-1,-1};
1037
1038 bool bacc = true;
1039 bool binp = true;
1040 Double_t eta_d[2] = {0.0,0.0};
1041 Double_t phi_d[2] = {0.0,0.0};
1042 for (Int_t daughter_index=0; daughter_index<2; daughter_index++){
1043 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(daughter[daughter_index]));
1044 eta_d[daughter_index] = dmc->Eta();
1045 phi_d[daughter_index] = dmc->Phi();
1046 if(!(dmc->PdgCode()==22)) binp = false;
1047 if(!(dmc->PdgCode()==22 &&
1048 eta_d[daughter_index]>fEtamin && eta_d[daughter_index]<fEtamax &&
1049 phi_d[daughter_index]>fPhimin && phi_d[daughter_index]<fPhimax)) bacc = false;
1050
1051 if(dmc->GetFirstDaughter()>0 && dmc->GetLastDaughter()>0) {
1052 // get the photons's daughters...
1053 const AliMCParticle *dmcd1 = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dmc->GetFirstDaughter()));
1054 const AliMCParticle *dmcd2 = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dmc->GetLastDaughter()));
1055 Double_t productionR1 = TMath::Sqrt(dmcd1->Xv()*dmcd1->Xv() + dmcd1->Yv()*dmcd1->Yv());
1056 if(bacc) h1_eConversionR->Fill(productionR1);
1057 // check if this is a conversion...
1058 if( (dmcd1->PdgCode()== -1.0*dmcd2->PdgCode()) &&
1059 (dmcd1->PdgCode()==11 || dmcd1->PdgCode()==-11) &&
1060 productionR1<440.0){
1061 //find the conv e with highest energy, assign it to be that photon decay product.
1062 if( dmcd1->E() > dmcd2->E() )
1063 eIndexofConvertedPhoton[daughter_index] = dmc->GetFirstDaughter();
1064 else
1065 eIndexofConvertedPhoton[daughter_index] = dmc->GetLastDaughter();
1066 }
1067 }
1068 }
1069
1070 if(binp!=true)
1071 continue;
1072
1073 // primary particle
1074 //Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1075 // (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1076 //if(dR <= 0.01) isPrimary = 1;
1077 //else isPrimary = 0;
1078
1079 isMaterialSec = 0;
1080 if(isPrimary!=1){
1081 if(((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 2212 || //proton
1082 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -2212 || //anti-proton
1083 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 2112 || //neutron
1084 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -2112 || //anti-neutron
1085 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 321 || //K+
1086 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -321 || //K-
1087 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 211 || //pi+
1088 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -211 //pi-
1089 )
1090 isMaterialSec = 1;
1091 }
1092
1093 h1_Pi0TruthPt ->Fill(mcP->Pt());
1094 if(isK0sDecay) h1_K0Pi0TruthPt->Fill(mcP->Pt());
1095 h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
1096
1097 if(isPrimary==1){
1098 h1_PriPi0TruthPt ->Fill(mcP->Pt());
1099 h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
1100 }
1101 if(isPrimary!=1 && isMaterialSec!=1) h1_PhysPi0TruthPt->Fill(mcP->Pt());
1102
1103 if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
1104 continue;
1105
1106 h1_Pi0TruthPtPhi2piEta1 ->Fill(mcP->Pt());
1107 if(isPrimary==1)
1108 h1_PriPi0TruthPtPhi2piEta1 ->Fill(mcP->Pt());
1109 if(isK0sDecay)
1110 h1_K0Pi0TruthPtPhi2piEta1 ->Fill(mcP->Pt());
1111 if(isPrimary!=1 && isMaterialSec!=1)
1112 h1_PhysPi0TruthPtPhi2piEta1 ->Fill(mcP->Pt());
1113
1114 if(mcP->Eta()>fEtamin && mcP->Eta()<fEtamax){
1115 h1_Pi0TruthPtPhi2piEta065 ->Fill(mcP->Pt());
1116 if(isPrimary==1)
1117 h1_PriPi0TruthPtPhi2piEta065 ->Fill(mcP->Pt());
1118 if(isK0sDecay)
1119 h1_K0Pi0TruthPtPhi2piEta065 ->Fill(mcP->Pt());
1120 if(isPrimary!=1 && isMaterialSec!=1)
1121 h1_PhysPi0TruthPtPhi2piEta065 ->Fill(mcP->Pt());
1122 }
1123
1124
1125 if(binp && bacc){// 2 Photons hit the EMCAL!
1126
1127 Int_t Nfoundphotons = 0;
1128 Int_t iFoundphotons[10] = {0,0,0,0,0,0,0,0,0,0};
1129 Int_t Nfoundelectrons = 0;
1130 Int_t iFoundelectrons[10] = {0,0,0,0,0,0,0,0,0,0};
1131 for (Int_t daughter_index=0; daughter_index<2; daughter_index++){//both truth photons. (also includes conversions..)
1132
1133 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(daughter[daughter_index]));
1134
1135 h1_TruthPhotonsEmcal->Fill(dmc->Pt());
1136 h2_TruthPhotonsPhiEta->Fill(dmc->Phi(),dmc->Eta());
1137
1138 for(int i=0; i<nclusters; i++) {
1139 Bool_t matches_pion_photon = 0;
1140
1141 AliESDCaloCluster* esdCluster=NULL;
1142 AliAODCaloCluster* aodCluster=NULL;
1143 if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
1144 else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
1145
1146 Double_t clustMC_phi, clustMC_eta;
1147
1148 if(fEsdEv){
1149
1150 if(esdCluster->IsEMCAL()){
1151
1152 Float_t pos[3] = {0,0,0};
1153 esdCluster->GetPosition(pos);
1154 TVector3 vpos(pos);
1155 //h1_Phi->Fill(vpos.Phi());
1156 clustMC_phi = vpos.Phi();
1157 clustMC_eta = vpos.Eta();
1158
1159 Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
1160 (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
1161 h1_dR_RealMC->Fill(dR);
1162 matches_pion_photon = 0;
1163 //if(dR<=0.04) matches_pion_photon = 1;
1164
1165
1166 TArrayI *TruthLabelsA = esdCluster->GetLabelsArray();
1167 if(TruthLabelsA){
1168 Int_t trackindex = TruthLabelsA->At(0);
1169 if( trackindex==daughter[daughter_index] ){
1170 matches_pion_photon = 1;
1171 iFoundphotons[Nfoundphotons] = i;
1172 Nfoundphotons++;
1173 }
1174 else if( trackindex==eIndexofConvertedPhoton[daughter_index] ){
1175 iFoundelectrons[Nfoundelectrons] = i;
1176 Nfoundelectrons++;
1177 }
1178 AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(trackindex));
1179
1180 if(matches_pion_photon){
1181
1182 h1_PhotonsEmcal->Fill(esdCluster->E());
1183 h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
1184 if(esdCluster->GetNCells()>=2)
1185 h1_PhotonsNCellsCut->Fill(esdCluster->E());
1186 if(esdCluster->GetNTracksMatched()==0)
1187 h1_PhotonsTrackMatchCut->Fill(esdCluster->E());
1188 if(esdCluster->GetNCells()>=2 && esdCluster->GetNTracksMatched()==0)
1189 h1_PhotonsAllCut->Fill(esdCluster->E());
1190
1191 if(esdCluster->GetNCells()>=2){
1192 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1193 h2_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()));
1194 }
1195 }//if(matches_pion_photon)
1196
1197 }//if Truthlabels exists
1198 vpos.Delete();
1199 }//if(IsEMCAL())
1200
1201 }//if(fEsdEv)
1202 else if(fAodEv){
1203
1204 if(aodCluster->IsEMCAL()){
1205
1206 Float_t pos[3] = {0,0,0};
1207 aodCluster->GetPosition(pos);
1208 TVector3 vpos(pos);
1209 //h1_Phi->Fill(vpos.Phi());
1210 clustMC_phi = vpos.Phi();
1211 clustMC_eta = vpos.Eta();
1212
1213 Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
1214 (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
1215 h1_dR_RealMC->Fill(dR);
1216 matches_pion_photon = 0;
1217 if(dR<=0.04) matches_pion_photon = 1;
1218
1219 //TArrayI *TruthLabelsA = aodCluster->GetLabelsArray();
1220 //if(TruthLabelsA){
1221 // Int_t trackindex = TruthLabelsA->At(0);
1222 // if( trackindex==daughter[daughter_index] )
1223 // matches_pion_photon = 1;
1224 // AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(trackindex));
1225
1226 if(matches_pion_photon){
1227
1228 h1_PhotonsEmcal->Fill(aodCluster->E());
1229 h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
1230 if(aodCluster->GetNCells()>=2)
1231 h1_PhotonsNCellsCut->Fill(aodCluster->E());
1232 if(aodCluster->GetNTracksMatched()==0)
1233 h1_PhotonsTrackMatchCut->Fill(aodCluster->E());
1234 if(aodCluster->GetNCells()>=2 && aodCluster->GetNTracksMatched()==0)
1235 h1_PhotonsAllCut->Fill(aodCluster->E());
1236
1237 //if(aodCluster->GetNCells()>=2){
1238 //recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1239 //h2_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()));
1240 //}
1241 }//if(matches_pion_photon)
1242
1243 //}//if Truthlabels exists
1244 vpos.Delete();
1245 }//if(IsEMCAL())
1246
1247 }//if(fAodEv)
1248
1249 }//loop over nclusters.
1250
1251 }//both truth photons.
1252
1253 if(Nfoundphotons>1){
1254 AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster
1255 AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundphotons[1]); // pointer to EMCal cluster
1256
1257 if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){
1258
1259 recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1260 E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1261 fEsdEv->GetVertex()->GetXYZ(vertex);
1262 esdCluster1->GetMomentum(Photon1,vertex);
1263 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1264 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1265 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1266
1267 recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator);
1268 E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ
1269 fEsdEv->GetVertex()->GetXYZ(vertex);
1270 esdCluster2->GetMomentum(Photon2,vertex);
1271 Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ
1272 Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ
1273 Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ
1274
1275 Parent = TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2);
1276
1277 //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() );
1278 Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1279 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1280 if(isPrimary==1){
1281 //cout << "Primary production vertex: " << productionR << endl;
1282 h2_Mpt_Pri->Fill(Parent.M(),Parent.Pt());
1283 }
1284 else{
1285 //cout << "Secondary production vertex: " << productionR << endl;
1286 h2_Mpt_Sec ->Fill(Parent.M(),Parent.Pt());
1287 h3_MptR_Sec->Fill(Parent.M(),Parent.Pt(),productionR);
1288 if(isK0sDecay)
1289 h3_MptR_K0s->Fill(Parent.M(),Parent.Pt(),productionR);
1290 if(isMaterialSec)
1291 h3_MptR_Mat->Fill(Parent.M(),Parent.Pt(),productionR);
1292 }
1293 }//both good clusters
1294 }//found 2 photons.
1295 else if(Nfoundphotons==1){
1296 int mergedPion = 0;
1297 AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster
1298
1299 TArrayI *TruthLabelsA = esdCluster1->GetLabelsArray();
1300 if(TruthLabelsA){
1301 if(TruthLabelsA->GetSize()>1){
1302 Int_t trackindex[2];
1303 trackindex[0] = TruthLabelsA->At(0);
1304 trackindex[1] = TruthLabelsA->At(1);
1305 if( (trackindex[0]==daughter[0] && trackindex[1]==daughter[1]) ||
1306 (trackindex[0]==daughter[1] && trackindex[1]==daughter[0]) ){
1307 mergedPion = 1;
1308 }
1309 if(mergedPion==1){
1310 recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1311 E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1312 Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1313 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1314 h2_PtR_MatM->Fill(E1,productionR);
1315 }//if merged pion.
1316 }//truthlabel.size > 1
1317 }//if truthlabels
1318 }// Nfoundphotons==1
1319
1320 if(Nfoundphotons==1 && Nfoundelectrons==1){
1321 AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster
1322 AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundelectrons[0]); // pointer to EMCal cluster
1323
1324 if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){
1325
1326 recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1327 E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1328 fEsdEv->GetVertex()->GetXYZ(vertex);
1329 esdCluster1->GetMomentum(Photon1,vertex);
1330 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1331 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1332 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1333
1334 recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator);
1335 E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ
1336 fEsdEv->GetVertex()->GetXYZ(vertex);
1337 esdCluster2->GetMomentum(Photon2,vertex);
1338 Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ
1339 Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ
1340 Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ
1341
1342 Parent = TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2);
1343
1344 //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() );
1345 Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1346 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1347 if(isPrimary==1){
1348 //cout << "Primary production vertex: " << productionR << endl;
1349 h2_Mpt_Pri_conv->Fill(Parent.M(),Parent.Pt());
1350 }
1351 else{
1352 //cout << "Secondary production vertex: " << productionR << endl;
1353 h2_Mpt_Sec_conv ->Fill(Parent.M(),Parent.Pt());
1354 h3_MptR_Sec_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1355 if(isK0sDecay)
1356 h3_MptR_K0s_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1357 if(isMaterialSec)
1358 h3_MptR_Mat_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1359 }
1360 }//both good clusters
1361 }// Nfoundphotons==1 && Nfoundelectrons==1
1362 else if(Nfoundelectrons==2){
1363 AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundelectrons[0]); // pointer to EMCal cluster
1364 AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundelectrons[1]); // pointer to EMCal cluster
1365
1366 if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){
1367
1368 recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1369 E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1370 fEsdEv->GetVertex()->GetXYZ(vertex);
1371 esdCluster1->GetMomentum(Photon1,vertex);
1372 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1373 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1374 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1375
1376 recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator);
1377 E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ
1378 fEsdEv->GetVertex()->GetXYZ(vertex);
1379 esdCluster2->GetMomentum(Photon2,vertex);
1380 Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ
1381 Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ
1382 Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ
1383
1384 Parent = TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2);
1385
1386 //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() );
1387 Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1388 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1389 if(isPrimary==1){
1390 //cout << "Primary production vertex: " << productionR << endl;
1391 h2_Mpt_Pri_conv->Fill(Parent.M(),Parent.Pt());
1392 }
1393 else{
1394 //cout << "Secondary production vertex: " << productionR << endl;
1395 h2_Mpt_Sec_conv ->Fill(Parent.M(),Parent.Pt());
1396 h3_MptR_Sec_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1397 if(isK0sDecay)
1398 h3_MptR_K0s_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1399 if(isMaterialSec)
1400 h3_MptR_Mat_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1401 }
1402 }//both good clusters
1403 }// Nfoundelectrons==2
1404
1405 h1_Pi0TruthPtEmcal ->Fill(mcP->Pt());
1406 if(isK0sDecay)
1407 h1_K0Pi0TruthPtEmcal ->Fill(mcP->Pt());
1408 h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1409
1410 if(isPrimary==1){
1411 h1_PriPi0TruthPtEmcal ->Fill(mcP->Pt());
1412 h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1413 }
1414 if(isPrimary!=1 && isMaterialSec!=1) h1_PhysPi0TruthPtEmcal->Fill(mcP->Pt());
1415
1416
1417 }// 2 Photons hit the EMCAL!
1418
1419 }//for(nTracksMC)
1420
1421 }//if(isMC)
1422
1423 //######################### ~~~~~~~~~~~~ ##################################
1424 //######################### DONE WITH MC ##################################
1425 //######################### ~~~~~~~~~~~~ ##################################
1426
1427
1428 for(int i=0; i<nclusters; i++) {
1429
1430 AliESDCaloCluster* esdCluster=NULL;
1431 AliAODCaloCluster* aodCluster=NULL;
1432 if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
1433 else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
1434 if(!esdCluster && !aodCluster) {
1435 AliError(Form("ERROR: Could not retrieve any (ESD or AOD) Cluster %d",i));
1436 continue;
1437 }
1438
1439 if(fEsdEv){
1440
1441 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1442
1443 //uncomment this to do the track matching (1 of 3 lines, esd part)!!
1444 //Bool_t MatchesToTrack = 0;
1445 if(esdCluster->IsEMCAL()){
1446
1447 Float_t pos[3] = {0,0,0};
1448 Short_t maxCellID = -1;
1449 Float_t celleta, cellphi;
1450 esdCluster->GetPosition(pos);
1451 TVector3 clusterPosition(pos);
1452 h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1453 GetMaxCellEnergy(esdCluster, maxCellID);
1454 AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1455 fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1456 h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1457
1458 // _______________Track loop for reconstructed event_____________
1459 for(Int_t itrk = 0; itrk < fEsdEv->GetNumberOfTracks(); itrk++) {
1460 AliESDtrack* esdTrack = fEsdEv->GetTrack(itrk); // pointer to reconstructed to track
1461 if(!esdTrack) {
1462 AliError(Form("ERROR: Could not retrieve any (ESD) track %d",itrk));
1463 continue;
1464 }
1465
1466 Double_t posTrk[3] = {0,0,0};
1467 esdTrack->GetXYZ(posTrk);
1468 TVector3 vposTrk(posTrk);
1469
1470 Double_t fMass = 0.139;
1471 Double_t fStepSurface = 20.;
1472 Float_t etaproj, phiproj, pttrackproj;
1473
1474 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1475 if(!trackParam) continue;
1476 AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1477
1478 double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1479 (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1480
1481 h1_dR_ClustTrk->Fill(dR_clusttrk);
1482
1483 //uncomment this to do the track matching (2 of 3 lines, esd part)!!
1484 //if(dR_clusttrk<fdRmin_ClustTrack)
1485 //MatchesToTrack = 1;
1486
1487 }//_____________________________nTracks__________________________
1488
1489 h2_cells_M02 ->Fill(esdCluster->GetNCells(),esdCluster->GetM02());
1490 h2_Ellipse ->Fill(esdCluster->GetM20(),esdCluster->GetM02());
1491 h1_Chi2 ->Fill(esdCluster->Chi2());//always -1.
1492 h1_nTrkMatch ->Fill(esdCluster->GetNTracksMatched());
1493 h1_ClusterDisp->Fill(esdCluster->GetDispersion());
1494 h2_E_time ->Fill(esdCluster->E(),esdCluster->GetTOF());
1495
1496 TArrayI *TrackLabels = esdCluster->GetTracksMatched();
1497 if(TrackLabels){
1498 if(TrackLabels->GetSize()>0){
1499 Int_t trackindex = TrackLabels->At(0);
1500 AliESDtrack* matchingT = fEsdEv->GetTrack(trackindex); // pointer to reconstructed to track
1501
1502 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1503 h2_eop_E ->Fill(esdCluster->E()*recalScale, esdCluster->E()*recalScale/matchingT->P());
1504 h2_eop_pT->Fill(matchingT->Pt(), esdCluster->E()*recalScale/matchingT->P());
1505 }
1506 }
1507
1508 //uncomment this to do the track matching (2 of 3 lines, esd part)!!
1509 //if(isGoodEsdCluster(esdCluster) && !MatchesToTrack){
1510 if(isGoodEsdCluster(esdCluster)){
1511 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1512 E1 = esdCluster->E()*recalScale;// TOTAL HACK - JJ
1513 fEsdEv->GetVertex()->GetXYZ(vertex);
1514 esdCluster->GetMomentum(Photon1,vertex);
1515 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1516 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1517 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1518 Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1519 h1_E->Fill(E1);
1520 h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1521 h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1522 }
1523 clusterPosition.Delete();
1524 }//if(esdCluster->isEMCAL())
1525 }//if(fEsdEv)
1526 else if(fAodEv){
1527
1528 recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1529
1530 //uncomment this to do the track matching (1 of 3 lines, aod part)!!
1531 //Bool_t MatchesToTrack = 0;
1532 if(aodCluster->IsEMCAL()){
1533
1534 Float_t pos[3] = {0,0,0};
1535 Short_t maxCellID = -1;
1536 Float_t celleta, cellphi;
1537 aodCluster->GetPosition(pos);
1538 TVector3 clusterPosition(pos);
1539 h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1540 GetMaxCellEnergy(aodCluster, maxCellID);
1541 AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1542 fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1543 h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1544
1545 // _______________Track loop for reconstructed event_____________
1546 for(Int_t itrk = 0; itrk < fAodEv->GetNumberOfTracks(); itrk++) {
1547 AliAODTrack* aodTrack = fAodEv->GetTrack(itrk); // pointer to reconstructed to track
1548 if(!aodTrack) {
1549 AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk));
1550 continue;
1551 }
1552
1553 Double_t posTrk[3] = {0,0,0};
1554 Double_t momTrk[3] = {0,0,0};
1555 aodTrack->GetXYZ(posTrk);
1556 aodTrack->GetPxPyPz(momTrk);
1557 //TVector3 vposTrk(posTrk);
1558
1559 //####################################################################################################
1560 //
1561 // commented all this stuff just to satisfy aliroot warnings.
1562 // but I may need it again if I want to do the track matching for aods.
1563 /*
1564 Double_t fMass = 0.139;
1565 Double_t fStepSurface = 20.;
1566 Float_t etaproj=0.0;
1567 Float_t phiproj=0.0;
1568 Float_t pttrackproj=0.0;
1569
1570 Double_t cv[21] = {0.0};
1571 aodTrack->GetCovarianceXYZPxPyPz(cv);
1572 AliExternalTrackParam *trackParam = new AliExternalTrackParam(posTrk,momTrk,cv,aodTrack->Charge());
1573 //AliExternalTrackParam emcalParam(*trackParam);
1574 //AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam());
1575 if(!trackParam) continue;
1576 ////AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1577 //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&emcalParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1578 delete trackParam;
1579
1580 //Constantin's implementation... gives funny result.
1581 //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(aodTrack,440.0);
1582 //phiproj = aodTrack->GetTrackPhiOnEMCal();
1583 //etaproj = aodTrack->GetTrackPhiOnEMCal();
1584
1585 double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1586 (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1587
1588 h1_dR_ClustTrk->Fill(dR_clusttrk);
1589 */
1590 //####################################################################################################
1591
1592 //uncomment this to do the track matching (2 of 3 lines, aod part)!!
1593 //if(dR_clusttrk<fdRmin_ClustTrack)
1594 //MatchesToTrack = 1;
1595
1596
1597 }//_____________________________nTracks__________________________
1598
1599 h2_cells_M02 ->Fill(aodCluster->GetNCells(),aodCluster->GetM02());
1600 h2_Ellipse ->Fill(aodCluster->GetM20(),aodCluster->GetM02());
1601 h1_Chi2 ->Fill(aodCluster->Chi2());//always -1.
1602 h1_nTrkMatch ->Fill(aodCluster->GetNTracksMatched());
1603 h1_ClusterDisp->Fill(aodCluster->GetDispersion());
1604 h2_E_time ->Fill(aodCluster->E(),aodCluster->GetTOF());
1605
1606 // #################################################
1607 // track matching eop histograms are handled here...
1608 // #################################################
1609
1610 //uncomment this to do the track matching (3 of 3 lines, aod part)!!
1611 //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){
1612 if(isGoodAodCluster(aodCluster)){
1613 recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1614 E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ
1615 fAodEv->GetVertex(0)->GetXYZ(vertex);
1616 aodCluster->GetMomentum(Photon1,vertex);
1617 Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1618 Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1619 Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1620 Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1621 h1_E->Fill(E1);
1622 h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1623 h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1624 }
1625 clusterPosition.Delete();
1626 }//if(aodCluster->IsEMCAL())
1627 }//if(fAodEv)
1628
1629 }//loop over nclusters.
1630
1631 //Make same event pions...
1632 for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1633 for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){
1634 Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j];
1635 Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1636 Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1637 Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/
1638 (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt());
1639 Int_t asymCut = 0;
1640 if (pairasym<0.1) asymCut = 1;
1641 else if(pairasym<0.7) asymCut = 2;
1642 else asymCut = 3;
1643
1644 h1_M ->Fill(Parent.M());
1645 h3_MptAsymm ->Fill(Parent.M(),Parent.Pt(),asymCut);
1646 h2_dphi_deta->Fill(deltaphi,deltaeta);
1647 }
1648 }
1649
1650 //Make mixed event...
1651 for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1652 for(unsigned int ipool=1; ipool<poolDepth; ipool++){
1653 for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){
1654 iskip = randy.Integer(Photons[0][izvtx][imult].size());
1655 if(j==iskip) continue;
1656 Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j];
1657 Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1658 Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1659 Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/
1660 (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt());
1661 Int_t asymCut = 0;
1662 if (pairasym<0.1) asymCut = 1;
1663 else if(pairasym<0.7) asymCut = 2;
1664 else asymCut = 3;
1665
1666 h1_M_mix ->Fill(Parent.M());
1667 h3_MptAsymm_mix ->Fill(Parent.M(),Parent.Pt(),asymCut);
1668 h2_dphi_deta_mix->Fill(deltaphi,deltaeta);
1669 }
1670 }
1671 }
1672
1673 for(int ipool=poolDepth-1; ipool>0; ipool--){
1674 Photons[ipool][izvtx][imult].clear();
1675 for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++)
1676 Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]);
1677 }
1678 Photons[0][izvtx][imult].clear();
1679
1680
1681
1682 // NEW HISTO should be filled before this point, as PostData puts the
1683 // information for this iteration of the UserExec in the container
1684 PostData(1, fOutput);
1685 }
1686
1687//________________________________________________________________________
1688void AliAnalysisTaskEMCALMesonGGSDM::Terminate(Option_t *) //specify what you want to have done
1689{
1690 // Called once at the end of the query.
1691
1692}
1693
1694//________________________________________________________________________
1695Int_t AliAnalysisTaskEMCALMesonGGSDM::GetZvtxBin(Double_t vertZ)
1696{
1697
1698 int izvtx = -1;
1699
1700 if (vertZ<-35)
1701 izvtx=0;
1702 else if(vertZ<-30)
1703 izvtx=1;
1704 else if(vertZ<-25)
1705 izvtx=2;
1706 else if(vertZ<-20)
1707 izvtx=3;
1708 else if(vertZ<-15)
1709 izvtx=4;
1710 else if(vertZ<-10)
1711 izvtx=5;
1712 else if(vertZ< -5)
1713 izvtx=6;
1714 else if(vertZ< 0)
1715 izvtx=7;
1716 else if(vertZ< 5)
1717 izvtx=8;
1718 else if(vertZ< 10)
1719 izvtx=9;
1720 else if(vertZ< 15)
1721 izvtx=10;
1722 else if(vertZ< 20)
1723 izvtx=11;
1724 else if(vertZ< 25)
1725 izvtx=12;
1726 else if(vertZ< 30)
1727 izvtx=13;
1728 else if(vertZ< 35)
1729 izvtx=14;
1730 else
1731 izvtx=15;
1732
1733 return izvtx;
1734}
1735
1736//________________________________________________________________________
1737Int_t AliAnalysisTaskEMCALMesonGGSDM::GetMultBin(Int_t mult){
1738
1739 int imult = -1;
1740
1741 if (mult<2)
1742 imult=0;
1743 else if(mult<25)
1744 imult=mult-2;
1745 else
1746 imult=24;
1747
1748 return imult;
1749}
1750
1751//________________________________________________________________________
1752Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodEsdCluster(AliESDCaloCluster* esdclust){
1753
1754 int pass = 1;
1755 int nMinCells = 2;
1756 double MinE = 0.4;
1757 //double MinErat = 0;
1758 //double MinEcc = 0;
1759
1760 if (!esdclust)
1761 pass = 0;
1762 if (!esdclust->IsEMCAL())
1763 pass = 0;
1764 if (esdclust->E()<MinE)
1765 pass = 0;
1766 if (esdclust->GetNCells()<nMinCells)
1767 pass = 0;
1768 //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
1769 //pass = 0;
1770 //if (esdclust->Chi2()<MinEcc) // eccentricity cut
1771 //pass = 0;//this is always -1.
1772
1773 //if(esdclust->GetM02()<0.1)
1774 // pass = 0;
1775 //if(esdclust->GetM02()>0.5)
1776 // pass = 0;
1777
1778 Float_t pos[3] = {0,0,0};
1779 esdclust->GetPosition(pos);
1780 TVector3 clusterPosition(pos);
1781 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1782 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1783 pass = 0;
1784 clusterPosition.Delete();
1785
1786 //DOING THIS BY HAND NOW...
1787 //if(!esdclust->GetNTracksMatched()==0)
1788 //pass = 0;
1789
1790 return pass;
1791}
1792
1793//________________________________________________________________________
1794Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodAodCluster(AliAODCaloCluster* aodclust){
1795
1796 int pass = 1;
1797 int nMinCells = 2;
1798 double MinE = 0.4;
1799 //double MinErat = 0;
1800 //double MinEcc = 0;
1801
1802 if (!aodclust)
1803 pass = 0;
1804 if (!aodclust->IsEMCAL())
1805 pass = 0;
1806 if (aodclust->E()<MinE)
1807 pass = 0;
1808 if (aodclust->GetNCells()<nMinCells)
1809 pass = 0;
1810 //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1811 //pass = 0;
1812 //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1813 //pass = 0;//this is always -1.
1814
1815 //if(aodclust->GetM02()<0.1)
1816 //pass = 0;
1817 //if(aodclust->GetM02()>0.5)
1818 //pass = 0;
1819
1820 Float_t pos[3] = {0,0,0};
1821 aodclust->GetPosition(pos);
1822 TVector3 clusterPosition(pos);
1823 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1824 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1825 pass = 0;
1826 clusterPosition.Delete();
1827
1828 //DOING THIS BY HAND NOW...
1829 //if(!aodclust->GetNTracksMatched()==0)
1830 //pass = 0;
1831
1832 return pass;
1833}
1834
1835//________________________________________________________________________
1836Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1837
1838 double dphi = p1.Phi() - p2.Phi();
1839
1840 if(dphi<0.5*TMath::Pi())
1841 dphi = dphi + 2.0*TMath::Pi();
1842
1843 if(dphi>1.5*TMath::Pi())
1844 dphi = dphi - 2.0*TMath::Pi();
1845
1846 return dphi;
1847}
1848
1849//________________________________________________________________________
1850Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1851
1852 double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1853
1854 return deta;
1855}
1856
1857
1858//________________________________________________________________________
1859Double_t AliAnalysisTaskEMCALMesonGGSDM::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1860
1861 double recalibfactor = 0.0;
1862
1863 if(iCalib==0){// no recalibration!
1864 recalibfactor = 1.0;
1865 }
1866 else if(iCalib==1){// just a scale factor:
1867 recalibfactor = 0.984;
1868 }
1869 else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1870 Double_t p[3] = {0.96968, -2.68720, -0.831607};
1871 recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1872 }
1873 else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1874 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};
1875 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1876 }
1877 else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1878 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};
1879 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1880 }
1881 else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1882 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};
1883 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) );
1884 }
1885 else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1886 Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1887 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1888 }
1889 else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1890 Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1891 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1892 }
1893 else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1894 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};
1895 //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
1896 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));
1897 }
1898 else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1899 Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1900 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1901 }
1902 else if(iCalib==10){// Jason played with test beam data
1903 Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1904 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1905 }
1906 else if(iCalib==11){// Jason played with test beam MC
1907 Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1908 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1909 }
1910
1911 return recalibfactor;
1912}
1913
1914
1915//________________________________________________________________________
1916Double_t AliAnalysisTaskEMCALMesonGGSDM::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1917{
1918 // Get maximum energy of attached cell.
1919
1920 id = -1;
1921 AliVCaloCells *fVCells=NULL;
1922 if(fEsdEv) fVCells = fEsdEv->GetEMCALCells();
1923 else if(fAodEv) fVCells = fAodEv->GetEMCALCells();
1924 if(!fVCells)
1925 return 0;
1926
1927 Double_t maxe = 0;
1928 Int_t ncells = cluster->GetNCells();
1929 for (Int_t i=0; i<ncells; i++) {
1930 Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1931 if (e>maxe) {
1932 maxe = e;
1933 id = cluster->GetCellAbsId(i);
1934 }
1935 }
1936 return maxe;
1937}
1938
1939
1940//________________________________________________________________________
1941Int_t AliAnalysisTaskEMCALMesonGGSDM::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){
1942
1943 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1944
1945 Int_t nPTracks= mcEvent->GetNumberOfPrimaries();
1946
1947 Int_t isPhysPrimary = 1;
1948 Int_t ismHF = 0;
1949 Int_t ismLongLivedOrK = 0;
1950
1951 if(mcP->GetMother()<0)//if it has no mother...
1952 return isPhysPrimary;
1953
1954 Int_t imTrack = mcP->GetMother();
1955 AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1956
1957 if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark...
1958 return isPhysPrimary;
1959
1960
1961 //############################################
1962 //get the PDG digits....
1963 int num = mcPm->PdgCode();
1964 int RevDigits[10] = {0};
1965 int nDigits = 0;
1966 while (num >= 1){
1967 RevDigits[nDigits++] = num%10;
1968 num = num / 10;
1969 }
1970 //##############################################
1971
1972
1973 if(RevDigits[3]>3)//Baryons
1974 ismHF = 1;
1975 else if(RevDigits[2]>3)//Mesons
1976 ismHF = 1;
1977
1978 ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1979
1980 if(!ismHF && ismLongLivedOrK)
1981 isPhysPrimary = 0;
1982 else{ // check grandmother, greatgrandmothers, etc...
1983 while(imTrack >= nPTracks){
1984
1985 if(mcPm->GetMother()<0)//if it has no mother...
1986 break;
1987
1988 if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark...
1989 return isPhysPrimary;
1990
1991 imTrack = mcPm->GetMother();
1992 mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1993
1994 //############################################
1995 //get the PDG digits....
1996 num = mcPm->PdgCode();
1997 for(int i=0; i<10; i++) RevDigits[i] = 0;
1998 nDigits = 0;
1999 while (num >= 1){
2000 RevDigits[nDigits++] = num%10;
2001 num = num / 10;
2002 }
2003 //##############################################
2004 if(RevDigits[3]>3)//Baryons
2005 ismHF = 1;
2006 else if(RevDigits[2]>3)//Mesons
2007 ismHF = 1;
2008
2009 ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
2010
2011 if(!ismHF && ismLongLivedOrK)
2012 isPhysPrimary = 0;
2013
2014 }//while( >=nPTracks)
2015 }
2016
2017 return isPhysPrimary;
2018}
2019
2020
2021//________________________________________________________________________
2022Int_t AliAnalysisTaskEMCALMesonGGSDM::IsLongLivedOrK(Int_t MyPDGcode){
2023
2024 Int_t MyFlag = 0;
2025
2026 if(
2027 (TMath::Abs(MyPDGcode) == 22 ) || // Photon
2028 (TMath::Abs(MyPDGcode) == 11 ) || // Electron
2029 (TMath::Abs(MyPDGcode) == 13 ) || // Muon(-)
2030 (TMath::Abs(MyPDGcode) == 211 ) || // Pion
2031 (TMath::Abs(MyPDGcode) == 321 ) || // Kaon
2032 (TMath::Abs(MyPDGcode) == 310 ) || // K0s
2033 (TMath::Abs(MyPDGcode) == 130 ) || // K0l
2034 (TMath::Abs(MyPDGcode) == 2212) || // Proton
2035 (TMath::Abs(MyPDGcode) == 2112) || // Neutron
2036 (TMath::Abs(MyPDGcode) == 3122) || // Lambda_0
2037 (TMath::Abs(MyPDGcode) == 3112) || // Sigma Minus
2038 (TMath::Abs(MyPDGcode) == 3222) || // Sigma Plus
2039 (TMath::Abs(MyPDGcode) == 3312) || // Xsi Minus
2040 (TMath::Abs(MyPDGcode) == 3322) || // Xsi
2041 (TMath::Abs(MyPDGcode) == 3334) || // Omega
2042 (TMath::Abs(MyPDGcode) == 12 ) || // Electron Neutrino
2043 (TMath::Abs(MyPDGcode) == 14 ) || // Muon Neutrino
2044 (TMath::Abs(MyPDGcode) == 16 ) ) // Tau Neutrino
2045 MyFlag = 1;
2046
2047 return MyFlag;
2048}
2049
2050
2051//________________________________________________________________________
2052Int_t AliAnalysisTaskEMCALMesonGGSDM::IsMyMCHeaderType(Int_t iTrack, char *MyType, AliMCEvent *mcEvent) const
2053{
2054
2055 Int_t isMyType = 0;
2056
2057 AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(mcEvent->GenEventHeader());
2058 if(!cocktail)
2059 return 0;
2060
2061 TList *genHeaders = cocktail->GetHeaders();
2062
2063 Int_t nGenerators = genHeaders->GetEntries();
2064 Int_t indexMyType = -1;
2065 Int_t startParticle=0;
2066
2067 for(Int_t igen = 0; igen < nGenerators; igen++){
2068 AliGenEventHeader* eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
2069 TString name = eventHeader2->GetName();
2070 startParticle += eventHeader2->NProduced();
2071 //cout << name << endl;
2072 if (name.Contains(MyType,TString::kIgnoreCase)){
2073 indexMyType = igen;
2074 startParticle -= eventHeader2->NProduced();
2075 break;
2076 }
2077 }
2078
2079 AliGenEventHeader *addedPi0Header = (AliGenEventHeader*)genHeaders->At(indexMyType);
2080 Int_t ipi0min = startParticle;
2081 Int_t ipi0max = ipi0min+addedPi0Header->NProduced()-1;
2082 if(iTrack >= ipi0min && iTrack <= ipi0max)
2083 isMyType = 1;
2084
2085 return isMyType;
2086}
2087
2088