]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskSDMGammaMC.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 / AliAnalysisTaskSDMGammaMC.cxx
CommitLineData
3a66195e 1#include "AliAnalysisTaskSDMGammaMC.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
2942f542 66using std::cout;
67using std::endl;
3a66195e 68
69ClassImp(AliAnalysisTaskSDMGammaMC)
70
71//________________________________________________________________________
72AliAnalysisTaskSDMGammaMC::AliAnalysisTaskSDMGammaMC() :
73 AliAnalysisTaskSE(),
74 fOutput(0),
75 fMcMode(0),
76 fRecalibrator(0),
77 fPhimin(0),
78 fPhimax(0),
79 fEtamin(0),
80 fEtamax(0),
81 fTrackCuts(0),
82 fEsdEv(0),
83 fAodEv(0),
84 h1_nClusters(0),
85 h1_zvtx(0),
86 h1_trigger(0),
87 h1_E(0),
88 h1_Phi(0),
89 h2_PiMotherID(0),
90 h2_GaMotherID(0),
91 h3_gE_RecTruth(0),
92 h3_gE_RecTruth_ncellscut(0),
93 h1_Pi0TruthPt(0),
94 h1_PriPi0TruthPt(0),
95 h1_Pi0TruthPtEmcal(0),
96 h1_PriPi0TruthPtEmcal(0),
97 h2_Pi0TruthPhiEta(0),
98 h2_PriPi0TruthPhiEta(0),
99 h2_Pi0TruthPhiEtaEmcal(0),
100 h2_PriPi0TruthPhiEtaEmcal(0),
101 h1_TruthPhotonsEmcal(0),
102 h2_TruthPhotonsPhiEta(0),
103 h1_PhotonsEmcal(0),
104 h1_PhotonsNCellsCut(0),
105 h1_PhotonsTrackMatchCut(0),
106 h1_PhotonsAllCut(0),
107 h2_PhotonsPhiEtaIsEmcal(0),
108 h1_dR_RealMC(0),
109 h1_Eta(0),
110 h1_Chi2(0),
111 h1_nTrkMatch(0),
112 h1_nCells(0),
113 h1_ClusterDisp(0),
114 h2_Ellipse(0),
115 h2_EtaPt(0),
116 h2_dphi_deta(0),
117 h2_dphi_deta_mix(0),
118 h2_DispRes(0),
119 h2_cells_M02(0),
120 TriggerList(0)
121{
122 // Dummy constructor ALWAYS needed for I/O.
123}
124
125//________________________________________________________________________
126AliAnalysisTaskSDMGammaMC::AliAnalysisTaskSDMGammaMC(const char *name) :
127 AliAnalysisTaskSE(name),
128 fOutput(0),
129 fMcMode(0),
130 fRecalibrator(0),
131 fPhimin(0),
132 fPhimax(0),
133 fEtamin(0),
134 fEtamax(0),
135 fTrackCuts(0),
136 fEsdEv(0),
137 fAodEv(0),
138 h1_nClusters(0),
139 h1_zvtx(0),
140 h1_trigger(0),
141 h1_E(0),
142 h1_Phi(0),
143 h2_PiMotherID(0),
144 h2_GaMotherID(0),
145 h3_gE_RecTruth(0),
146 h3_gE_RecTruth_ncellscut(0),
147 h1_Pi0TruthPt(0),
148 h1_PriPi0TruthPt(0),
149 h1_Pi0TruthPtEmcal(0),
150 h1_PriPi0TruthPtEmcal(0),
151 h2_Pi0TruthPhiEta(0),
152 h2_PriPi0TruthPhiEta(0),
153 h2_Pi0TruthPhiEtaEmcal(0),
154 h2_PriPi0TruthPhiEtaEmcal(0),
155 h1_TruthPhotonsEmcal(0),
156 h2_TruthPhotonsPhiEta(0),
157 h1_PhotonsEmcal(0),
158 h1_PhotonsNCellsCut(0),
159 h1_PhotonsTrackMatchCut(0),
160 h1_PhotonsAllCut(0),
161 h2_PhotonsPhiEtaIsEmcal(0),
162 h1_dR_RealMC(0),
163 h1_Eta(0),
164 h1_Chi2(0),
165 h1_nTrkMatch(0),
166 h1_nCells(0),
167 h1_ClusterDisp(0),
168 h2_Ellipse(0),
169 h2_EtaPt(0),
170 h2_dphi_deta(0),
171 h2_dphi_deta_mix(0),
172 h2_DispRes(0),
173 h2_cells_M02(0),
174 TriggerList(0)
175{
176 // Constructor
177 // Define input and output slots here (never in the dummy constructor)
178 // Input slot #0 works with a TChain - it is connected to the default input container
179 // Output slot #1 writes into a TH1 container
180
181
182 DefineOutput(1, TList::Class()); // for output list
183}
184
185//________________________________________________________________________
186AliAnalysisTaskSDMGammaMC::~AliAnalysisTaskSDMGammaMC()
187{
188 // Destructor. Clean-up the output list, but not the histograms that are put inside
189 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
190 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
191 delete fOutput;
192 }
193 delete fTrackCuts;
194}
195
196//________________________________________________________________________
197void AliAnalysisTaskSDMGammaMC::UserCreateOutputObjects()
198{
199 // Create histograms
200 // Called once (on the worker node)
201
202 fOutput = new TList();
203 fOutput->SetOwner(); // IMPORTANT!
204
205 fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
206
207 cout << "__________AliAnalysisTaskSDMGammaMC: Input settings__________" << endl;
208 cout << " fMcMode: " << fMcMode << endl;
209 cout << " fRecalibrator: " << fRecalibrator << endl;
210 cout << " phi range: " << fPhimin << ", " << fPhimax << endl;
211 cout << " eta range: " << fEtamin << ", " << fEtamax << endl;
212 cout << " number of zvtx bins: " << zvtx_bins << endl;
213 cout << " number of mult bins: " << mult_bins << endl;
214 cout << " poolDepth: " << poolDepth << endl;
215 cout << endl;
216
217
218 double TotalNBins = 0.0;
219
220 // Create histograms
221 Int_t nClustersbins = 501;
222 Float_t nClusterslow = -0.5, nClustersup = 500.5;
223 h1_nClusters = new TH1F("h1_nClusters", "# of clusters", nClustersbins, nClusterslow, nClustersup);
224 h1_nClusters->GetXaxis()->SetTitle("number of clusters/evt");
225 h1_nClusters->GetYaxis()->SetTitle("counts");
226 h1_nClusters->SetMarkerStyle(kFullCircle);
227 TotalNBins+=nClustersbins;
228
229 Int_t nZvertexbins = 501;
230 Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
231 h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
232 h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
233 h1_zvtx->GetYaxis()->SetTitle("counts");
234 h1_zvtx->SetMarkerStyle(kFullCircle);
235 TotalNBins+=nZvertexbins;
236
237 h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
238 TotalNBins+=1001;
239
240 Int_t ptbins = 2000;
241 Float_t ptlow = 0.0, ptup = 20.0;
242 Int_t Ebins = 1000;
243 Float_t Elow = 0.0, Eup = 20.0;
244 h1_E = new TH1F("h1_E", "Cluster Energy in EMCal", Ebins, Elow, Eup);
245 h1_E->GetXaxis()->SetTitle("E [GeV]");
246 h1_E->GetYaxis()->SetTitle("counts");
247 h1_E->SetMarkerStyle(kFullCircle);
248 TotalNBins+=Ebins;
249
250 h1_Phi = new TH1F("h1_Phi", "phi distribution", 1000, -7, 7);
251 h1_Phi->GetXaxis()->SetTitle("#phi [rad]");
252 h1_Phi->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
253 h1_Phi->SetMarkerStyle(kFullCircle);
254 TotalNBins+=1000;
255
256 h2_PiMotherID = new TH2F("h2_PiMotherID", "Mother ID for Truth Pi0's", 100001, -0.5,100000.5, 2, 0.5,2.5);
257 h2_PiMotherID->GetXaxis()->SetTitle("#pi^{0} Mother Particle ID");
258 h2_PiMotherID->GetYaxis()->SetTitle("primary or non-primary");
259 h2_PiMotherID->GetZaxis()->SetTitle("counts");
260 h2_PiMotherID->SetMarkerStyle(kFullCircle);
261 TotalNBins+=2*100001;
262
263 h2_GaMotherID = new TH2F("h2_GaMotherID", "Mother ID for Truth #gamma's", 100001, -0.5,100000.5, 2, 0.5,2.5);
264 h2_GaMotherID->GetXaxis()->SetTitle("#gamma Mother Particle ID");
265 h2_GaMotherID->GetYaxis()->SetTitle("primary or non-primary");
266 h2_GaMotherID->GetZaxis()->SetTitle("counts");
267 h2_GaMotherID->SetMarkerStyle(kFullCircle);
268 TotalNBins+=2*100001;
269
270 h3_gE_RecTruth = new TH3F("h3_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 1000,0,4, 4,0.5,4.5);
271 h3_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
272 h3_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
273 h3_gE_RecTruth->GetZaxis()->SetTitle("category");
274 h3_gE_RecTruth->SetMarkerStyle(kFullCircle);
275 TotalNBins+=Ebins*1000*4;
276
277 h3_gE_RecTruth_ncellscut = new TH3F("h3_gE_RecTruth_ncellscut", "#gamma E_{truth}/E_{clust} vs E_{clust} (for nCells>1)", Ebins,Elow,Eup, 1000,0,4, 4,0.5,4.5);
278 h3_gE_RecTruth_ncellscut->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
279 h3_gE_RecTruth_ncellscut->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
280 h3_gE_RecTruth_ncellscut->GetZaxis()->SetTitle("category");
281 h3_gE_RecTruth_ncellscut->SetMarkerStyle(kFullCircle);
282 TotalNBins+=Ebins*1000*4;
283
284 h1_Pi0TruthPt = new TH1F("h1_Pi0TruthPt", "P_{T} distribution for Truth Pi0's", ptbins, ptlow, ptup);
285 h1_Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
286 h1_Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
287 h1_Pi0TruthPt->SetMarkerStyle(kFullCircle);
288 TotalNBins+=ptbins;
289
290 h1_PriPi0TruthPt = new TH1F("h1_PriPi0TruthPt", "P_{T} distribution for Truth Primary Pi0's", ptbins, ptlow, ptup);
291 h1_PriPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
292 h1_PriPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
293 h1_PriPi0TruthPt->SetMarkerStyle(kFullCircle);
294 TotalNBins+=ptbins;
295
296 h1_Pi0TruthPtEmcal = new TH1F("h1_Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's (hit EMCal)", ptbins, ptlow, ptup);
297 h1_Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
298 h1_Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
299 h1_Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
300 TotalNBins+=ptbins;
301
302 h1_PriPi0TruthPtEmcal = new TH1F("h1_PriPi0TruthPtEmcal", "P_{T} distribution for Truth Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
303 h1_PriPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
304 h1_PriPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
305 h1_PriPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
306 TotalNBins+=ptbins;
307
308 h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
309 h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
310 h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
311 TotalNBins+=380*200;
312
313 h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
314 h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
315 h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
316 TotalNBins+=380*200;
317
318 h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
319 h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
320 h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
321 TotalNBins+=380*150;
322
323 h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
324 h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
325 h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
326 TotalNBins+=380*150;
327
328 h1_TruthPhotonsEmcal = new TH1F("h1_TruthPhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
329 h1_TruthPhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
330 h1_TruthPhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
331 h1_TruthPhotonsEmcal->SetMarkerStyle(kFullCircle);
332 TotalNBins+=ptbins;
333
334 h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
335 "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
336 h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
337 h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
338 h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
339 TotalNBins+=380*150;
340
341 h1_PhotonsEmcal = new TH1F("h1_PhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
342 h1_PhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
343 h1_PhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
344 h1_PhotonsEmcal->SetMarkerStyle(kFullCircle);
345 TotalNBins+=ptbins;
346
347 h1_PhotonsNCellsCut = new TH1F("h1_PhotonsNCellsCut", "P_{T} distribution for #gamma's that survive NCells cut", ptbins, ptlow, ptup);
348 h1_PhotonsNCellsCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
349 h1_PhotonsNCellsCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
350 h1_PhotonsNCellsCut->SetMarkerStyle(kFullCircle);
351 TotalNBins+=ptbins;
352
353 h1_PhotonsTrackMatchCut = new TH1F("h1_PhotonsTrackMatchCut", "P_{T} distribution for #gamma's that survive TrackMatch cut", ptbins, ptlow, ptup);
354 h1_PhotonsTrackMatchCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
355 h1_PhotonsTrackMatchCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
356 h1_PhotonsTrackMatchCut->SetMarkerStyle(kFullCircle);
357 TotalNBins+=ptbins;
358
359 h1_PhotonsAllCut = new TH1F("h1_PhotonsAllCut", "P_{T} distribution for #gamma's that survive All cut", ptbins, ptlow, ptup);
360 h1_PhotonsAllCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
361 h1_PhotonsAllCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
362 h1_PhotonsAllCut->SetMarkerStyle(kFullCircle);
363 TotalNBins+=ptbins;
364
365 h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
366 "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 150,-1.5,1.5);
367 h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
368 h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
369 TotalNBins+=380*150;
370
371 h1_dR_RealMC = new TH1F("h1_dR_RealMC", "P_{T} distribution for #gamma's that survive All cut", 2000, -0.01, 10);
372 h1_dR_RealMC->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
373 h1_dR_RealMC->GetYaxis()->SetTitle("N");
374 h1_dR_RealMC->SetMarkerStyle(kFullCircle);
375 TotalNBins+=2000;
376
377 Int_t etabins = 150;
378 Float_t etalow = -1.5, etaup = 1.5;
379 h1_Eta = new TH1F("h1_Eta","#eta distribution for reconstructed",etabins, etalow, etaup);
380 h1_Eta->GetXaxis()->SetTitle("#eta");
381 h1_Eta->GetYaxis()->SetTitle("counts");
382 TotalNBins+=etabins;
383
384 Int_t chi2bins = 100;
385 Float_t chi2low = -2, chi2up = 2;
386 h1_Chi2 = new TH1F("h1_Chi2","#chi^{2} distribution for reconstructed",chi2bins, chi2low, chi2up);
387 h1_Chi2->GetXaxis()->SetTitle("#chi^{2}");
388 h1_Chi2->GetYaxis()->SetTitle("counts");
389 TotalNBins+=chi2bins;
390
391 h1_nTrkMatch = new TH1F("h1_nTrkMatch","number of matched tracks",14, -1.5, 5.5);
392 h1_nTrkMatch->GetXaxis()->SetTitle("nTracksMatched");
393 h1_nTrkMatch->GetYaxis()->SetTitle("counts");
394 TotalNBins+=14;
395
396 h1_ClusterDisp = new TH1F("h1_ClusterDisp","Dispersion of CaloCluster",1000, -1, 3);
397 h1_ClusterDisp->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
398 h1_ClusterDisp->GetYaxis()->SetTitle("counts");
399 TotalNBins+=1000;
400
401 h2_Ellipse = new TH2F("h2_Ellipse","Ellipse axis M20 vs M02",500, -0.01, 1, 500, -0.01, 1);
402 h2_Ellipse->GetXaxis()->SetTitle("cluster->GetM20()");
403 h2_Ellipse->GetYaxis()->SetTitle("cluster->GetM02()");
404 h2_Ellipse->GetZaxis()->SetTitle("counts");
405 TotalNBins+=500*500;
406
407 h2_EtaPt = new TH2F("h2_EtaPt","Cluster Energy vs ",etabins, etalow, etaup, ptbins, ptlow, ptup);
408 h2_EtaPt->GetXaxis()->SetTitle("E [GeV]");
409 h2_EtaPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
410 TotalNBins+=etabins*ptbins;
411
412 h2_dphi_deta = new TH2F("h2_dphi_deta","#Delta#phi vs #Delta#eta", 349,-1.5,5, 400,-2.0,2.0);
413 h2_dphi_deta->GetXaxis()->SetTitle("#Delta#phi");
414 h2_dphi_deta->GetYaxis()->SetTitle("#Delta#eta");
415 TotalNBins+=349*400;
416
417 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);
418 h2_dphi_deta_mix->GetXaxis()->SetTitle("#Delta#phi");
419 h2_dphi_deta_mix->GetYaxis()->SetTitle("#Delta#eta");
420 TotalNBins+=349*400;
421
422 h2_DispRes = new TH2F("h2_DispRes", "zvtx info", 500,-0.01,1, 500,-0.1,2);
423 h2_DispRes->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
424 h2_DispRes->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
425 h2_DispRes->GetZaxis()->SetTitle("counts");
426 TotalNBins+=500*500;
427
428 h2_cells_M02 = new TH2F("h2_cells_M02", "nCells vs M02", 204,-1.5,100.5, 500,-1,1.5);
429 h2_cells_M02->GetXaxis()->SetTitle("nCells");
430 h2_cells_M02->GetYaxis()->SetTitle("M02");
431 h2_cells_M02->GetZaxis()->SetTitle("counts");
432 TotalNBins+=204*500;
433
434 cout << endl << "Total number of bins in booked histograms: " << TotalNBins << endl << endl;
435
436 //TFile *f = OpenFile(1);
437 //TDirectory::TContext context(f);
438
439 fOutput->Add(h1_nClusters);
440 fOutput->Add(h1_zvtx);
441 fOutput->Add(h1_trigger);
442 fOutput->Add(h1_E);
443 fOutput->Add(h1_Phi);
444 fOutput->Add(h2_PiMotherID);
445 fOutput->Add(h2_GaMotherID);
446 fOutput->Add(h3_gE_RecTruth);
447 fOutput->Add(h3_gE_RecTruth_ncellscut);
448 fOutput->Add(h1_Pi0TruthPt);
449 fOutput->Add(h1_PriPi0TruthPt);
450 fOutput->Add(h1_Pi0TruthPtEmcal);
451 fOutput->Add(h1_PriPi0TruthPtEmcal);
452 fOutput->Add(h2_Pi0TruthPhiEta);
453 fOutput->Add(h2_PriPi0TruthPhiEta);
454 fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
455 fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
456 fOutput->Add(h1_TruthPhotonsEmcal);
457 fOutput->Add(h2_TruthPhotonsPhiEta);
458 fOutput->Add(h1_PhotonsEmcal);
459 fOutput->Add(h1_PhotonsNCellsCut);
460 fOutput->Add(h1_PhotonsTrackMatchCut);
461 fOutput->Add(h1_PhotonsAllCut);
462 fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
463 fOutput->Add(h1_dR_RealMC);
464 fOutput->Add(h1_Eta);
465 fOutput->Add(h1_Chi2);
466 fOutput->Add(h1_nTrkMatch);
467 fOutput->Add(h1_ClusterDisp);
468 fOutput->Add(h2_Ellipse);
469 fOutput->Add(h2_EtaPt);
470 fOutput->Add(h2_dphi_deta);
471 fOutput->Add(h2_dphi_deta_mix);
472 fOutput->Add(h2_DispRes);
473 fOutput->Add(h2_cells_M02);
474
475 // Post data for ALL output slots >0 here,
476 // To get at least an empty histogram
477 // 1 is the outputnumber of a certain weg of task 1
478 PostData(1, fOutput);
479}
480
481//________________________________________________________________________
482void AliAnalysisTaskSDMGammaMC::UserExec(Option_t *)
483{
484 // Main loop Called for each event
485
486 AliMCEvent *mcEvent = MCEvent();
487 Bool_t isMC = bool(mcEvent);//is this the right way to do this?
488 if (!mcEvent){
489 cout << "no MC event" << endl;
490 return;
491 }
492
493 TRandom3 randy; randy.SetSeed(0);
494
495 double recalScale = 1.0;
496
497 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
498
499 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
500 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
501 if (!aodH && !esdH) Printf("ERROR: Could not get ESD or AODInputHandler");
502
503 if(esdH) fEsdEv = esdH->GetEvent();
504 else if(aodH) fAodEv = aodH->GetEvent();
505 else{
506 AliFatal("Neither ESD nor AOD event found");
507 return;
508 }
509
510
511 // get pointer to reconstructed event
512 AliVEvent *event = InputEvent();
513 if (!event){
514 AliError("Pointer == 0, this can not happen!"); return;}
515 //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
516 //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
517 //if (!fEsdEv){
518 //AliError("Cannot get the ESD event"); return;}
519
520 Int_t iTrigger = 0;
521 if (fEsdEv) iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
522 else if (fAodEv) iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
523 //h1_trigger->Fill(iTrigger);
524
525 char saythis[500];
526 Int_t iTriggerBin = 0;
527 for(unsigned long j=0; j<TriggerList.size(); j++){
528 if(iTrigger==TriggerList[j])
529 iTriggerBin=j+1;
530 }
531 if(iTriggerBin==0){
532 TriggerList.push_back(iTrigger);
533 iTriggerBin=TriggerList.size();
534 }
535
536 h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
537 sprintf(saythis,"%d",iTrigger);
538 h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
539
540 if(fEsdEv){
541 TString trigClasses = fEsdEv->GetFiredTriggerClasses();
542 // remove "fast cluster events":
543 if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL"))
544 return;
545 }
546 else if(fAodEv){
547 TString trigClasses = fAodEv->GetFiredTriggerClasses();
548 // remove "fast cluster events":
549 if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL"))
550 return;
551 }
552
553 if (fEsdEv){
554 if(!(fEsdEv->GetPrimaryVertex()->GetStatus())) return;
555 }
556 //else if (fAodEv){
557 //if(!(fAodEv->GetPrimaryVertex()->GetStatus())) return;
558 //}
559
560 Double_t vertDisp=0.0;
561 Double_t vertZres=0.0;
562 Bool_t vertIsfromZ=0;
563 if (fEsdEv){
564 vertDisp = fEsdEv->GetPrimaryVertex()->GetDispersion();
565 vertZres = fEsdEv->GetPrimaryVertex()->GetZRes();
566 vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
567 }
568 else if (fAodEv){
569 vertDisp = 0;
570 vertZres = 0;
571 vertIsfromZ = 0;
572 }
573
574 h2_DispRes->Fill(vertDisp, vertZres);
575 // if vertex is from spd vertexZ, require more stringent cut
576 if (vertIsfromZ) {
577 if (vertDisp>0.02 || vertZres>0.25 )
578 return; // bad vertex from VertexerZ
579 }
580
581 Int_t nclusters=0;
582 if(fEsdEv){
583 //Int_t evtN = fEsdEv->GetEventNumberInFile();
584 //Int_t ntracks = fEsdEv->GetNumberOfTracks();
585 nclusters = fEsdEv->GetNumberOfCaloClusters();
586 }
587 else if(fAodEv){
588 //Int_t evtN = fAodEv->GetEventNumberInFile();
589 //Int_t ntracks = fAodEv->GetNumberOfTracks();
590 nclusters = fAodEv->GetNumberOfCaloClusters();
591 }
592
593 // EMCal cluster loop for reconstructed event
594 //numberofclusters set above!
595 Double_t vertZ=0.0;
596 if (fEsdEv) vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
597 else if (fAodEv) vertZ = fAodEv->GetPrimaryVertex()->GetZ();
598
599 h1_zvtx->Fill(vertZ);
600 //zvertex cut:
601 if(fabs(vertZ)>10.0)
602 return;
603
604 h1_nClusters->Fill(nclusters);
605
606 //cout << iskip << " " << izvtx << " " << imult << endl;
607 //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
608
609
610
611 //######################### ~~~~~~~~~~~ ##################################
612 //######################### STARTING MC ##################################
613 //######################### ~~~~~~~~~~~ ##################################
614
615 if(isMC){
616 int isPrimary = 0;
617 int isK0sDecay = 0;
618
619 if (!mcEvent){
620 cout << "no MC event" << endl;
621 return;
622 }
623
624 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
625 if (!evtVtx)
626 return;
627
628 mcEvent->PreReadAll();
629
630 Int_t nTracksMC = mcEvent->GetNumberOfTracks();
631 Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
632
633 for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
634 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
635 if (!mcP)
636 continue;
637
638 if(iTrack<nPTracksMC) isPrimary = 1;
639 else isPrimary = 0;
640
641 if(mcP->PdgCode() == 22){
642 if(isPrimary==1){
643 if(mcP->GetMother()>-1)
644 h2_GaMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 1);
645 else
646 h2_GaMotherID->Fill(0.0,1);
647 }
648 else
649 h2_GaMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 2);
650 }
651
652 // it's a pion !!
653 if(mcP->PdgCode() != 111)
654 continue;
655
656 isK0sDecay = 0;
657 if(mcP->GetMother()>-1){
658 if( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 310 ||
659 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -310 )
660 isK0sDecay = 1;
661 }
662
663 // primary particle
664 //Double_t dR_vtx = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
665 // (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
666 //if(dR_vtx <= 0.01) isPrimary = 1;
667 //else isPrimary = 0;
668
669
670 if(isPrimary==1){
671 if(mcP->GetMother()>-1)
672 h2_PiMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 1);
673 else
674 h2_PiMotherID->Fill(0.0,1);
675 }
676 else
677 h2_PiMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 2);
678
679 h1_Pi0TruthPt ->Fill(mcP->Pt());
680 h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
681
682 if(isPrimary==1){
683 h1_PriPi0TruthPt ->Fill(mcP->Pt());
684 h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
685 }
686
687 if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
688 continue;
689
690 Int_t DecayPhotonLabel[2] = {mcP->GetFirstDaughter(),
691 mcP->GetLastDaughter() };
692
693 if (DecayPhotonLabel[0]<0) continue;
694 if (DecayPhotonLabel[1]<0) DecayPhotonLabel[1]=DecayPhotonLabel[0];
695 if (DecayPhotonLabel[1]-DecayPhotonLabel[0] != 1) continue;
696
697 bool bacc = true;
698 bool binp = true;
699 bool isConv[2] = {1,1};
700 Int_t convIndices[2][2] = { {-1,-1},{-1,-1} };
701 Double_t eta_d[2] = {0.0,0.0};
702 Double_t phi_d[2] = {0.0,0.0};
703 Int_t daughter_index = -1;
704 for (Int_t iPhoton=DecayPhotonLabel[0];iPhoton<=DecayPhotonLabel[1];++iPhoton){
705 if(iPhoton==DecayPhotonLabel[0]) daughter_index=0;
706 else daughter_index=1;
707 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(iPhoton));
708 eta_d[daughter_index] = dmc->Eta();
709 phi_d[daughter_index] = dmc->Phi();
710 if(!(dmc->PdgCode()==22)) binp = false;
711 if(!(eta_d[daughter_index]>fEtamin && eta_d[daughter_index]<fEtamax &&
712 phi_d[daughter_index]>fPhimin && phi_d[daughter_index]<fPhimax )) bacc = false;
713
714 if( ((TParticle*)dmc->Particle())->GetNDaughters() != 2 ) isConv[daughter_index] = 0;
715 else{//if photon has 2 daughters.
716
717 Int_t dd1 = dmc->GetFirstDaughter();
718 Int_t dd2 = dmc->GetLastDaughter();
719 if (dd2-dd1 != 1) cout << "How can this happen???? " << endl;
720 const AliMCParticle *dd1mc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dd1));
721 const AliMCParticle *dd2mc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dd2));
722 if( dd1mc->PdgCode() != -dd2mc->PdgCode() )
723 isConv[daughter_index] = 0;
724 else if( TMath::Abs(dd1mc->PdgCode())!=11 )
725 isConv[daughter_index] = 0;
726 if(isConv[daughter_index]==1){//store the e+e- indices...
727 convIndices[daughter_index][0] = dd1;
728 convIndices[daughter_index][1] = dd2;
729 }//if this photon converted.
730 }//close else-if photon has 2 daughters.
731 }//loop over 2 decay photons (iPhoton)
732
733 if(binp && bacc){// 2 Photons hit the EMCAL!
734
735 h1_Pi0TruthPtEmcal ->Fill(mcP->Pt());
736 h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
737
738 if(isPrimary==1){
739 h1_PriPi0TruthPtEmcal ->Fill(mcP->Pt());
740 h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
741 }
742
743 Int_t PhotonClusterMatch[2][3] = { {0,-1,-1},
744 {0,-1,-1} };
745 Int_t PhotonElectronMatch[2][6] = { {0,-1,-1,-1,-1,-1},
746 {0,-1,-1,-1,-1,-1} };
747
748 for(int iCluster=0; iCluster<nclusters; iCluster++) {
749
750 AliESDCaloCluster* esdCluster=NULL;
751 AliAODCaloCluster* aodCluster=NULL;
752 if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(iCluster); // pointer to EMCal cluster
753 else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(iCluster); // pointer to EMCal cluster
754
755 Double_t clustMC_phi, clustMC_eta;
756 if(fEsdEv){
757 if(esdCluster->IsEMCAL()){
758
759 if(!isGoodEsdCluster(esdCluster))
760 continue;
761
762 Float_t pos[3] = {0,0,0};
763 esdCluster->GetPosition(pos);
764 TVector3 vpos(pos);
765 h1_Phi->Fill(vpos.Phi());
766 clustMC_phi = vpos.Phi();
767 clustMC_eta = vpos.Eta();
768
769 Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
770 (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
771 h1_dR_RealMC->Fill(dR);
772 //matches_pion_photon = 0;
773 //if(dR<=0.04) matches_pion_photon = 1;
774
775 TArrayI *TruthLabelsA = esdCluster->GetLabelsArray();
776 if(TruthLabelsA){
777 for(int itl=0; itl<TruthLabelsA->GetSize(); itl++){
778
779 for(int iPhoton=0; iPhoton<2; iPhoton++){
780 if(TruthLabelsA->At(itl)==DecayPhotonLabel[iPhoton]){
781 PhotonClusterMatch[iPhoton][0] = 1;
782 PhotonClusterMatch[iPhoton][1] = DecayPhotonLabel[iPhoton];
783 PhotonClusterMatch[iPhoton][2] = iCluster;
784 }
785 }//loop over truth labels.
786
787 AliMCParticle *elecCandidate = (AliMCParticle*)(mcEvent->GetTrack(TruthLabelsA->At(itl)));
788 if(TMath::Abs(elecCandidate->PdgCode())==11){//if we have an electron...
789 Int_t elecMother_index = elecCandidate->GetMother();
790 if(elecMother_index>1 && elecMother_index<nTracksMC){
791 AliMCParticle *elecMother = (AliMCParticle*)(mcEvent->GetTrack(elecMother_index));
792 if( TMath::Abs(elecMother->PdgCode())==22 ){//if the e's mother is a photon...
793 Int_t elecGrandMother_index = elecMother->GetMother();
794 if(elecGrandMother_index==iTrack){//if the e's gMother is THE pi0 in question...
795 AliMCParticle *elecGrandMother = (AliMCParticle*)(mcEvent->GetTrack(elecGrandMother_index));
796 if( TMath::Abs(elecGrandMother->PdgCode())!=111 ) cout << "|| This can't happen!! A pion is a pion is a pion is a pion... ||" << endl;
797
798 for(int iPhoton=0; iPhoton<2; iPhoton++){
799 //if(convIndices[iPhoton][0]==elecMother_index){
800 if(convIndices[iPhoton][0]==TruthLabelsA->At(itl) || convIndices[iPhoton][1]==TruthLabelsA->At(itl)){
801 if(PhotonElectronMatch[iPhoton][1] == DecayPhotonLabel[iPhoton]) PhotonElectronMatch[iPhoton][0] = 2;
802 else PhotonElectronMatch[iPhoton][0] = 1;
803 PhotonElectronMatch[iPhoton][1] = DecayPhotonLabel[iPhoton];
804 if(PhotonElectronMatch[iPhoton][2]==-1) PhotonElectronMatch[iPhoton][2] = iCluster;//first cluster
805 else PhotonElectronMatch[iPhoton][3] = iCluster;//second cluster
806 if(PhotonElectronMatch[iPhoton][2]==-1) PhotonElectronMatch[iPhoton][4] = elecCandidate->PdgCode();
807 else PhotonElectronMatch[iPhoton][5] = elecCandidate->PdgCode();
808 }
809 }//loop over both decay photons
810
811 }//if it's THE pi0.
812 }//if we have a photon.
813 }//if we have an electron.
814 }//if the candidate has a real mother.
815
816 }//itl (TruthLabel loop)
817 }//if(TruthLabelsA exists)
818 }//if(isEMCal)
819 }//if(esdEv)
820 }//loop over clusters.
821
822
823 for(int iPhoton=0; iPhoton<2; iPhoton++){
824 AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(DecayPhotonLabel[iPhoton]));
825 if(!truthP)
826 continue;
827 if(PhotonClusterMatch[iPhoton][0]==1){
828 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonClusterMatch[iPhoton][2]);
829 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
830 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 1);
831 if(esdCluster->GetNCells()>=2)
832 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 1);
833 }
834 else if(PhotonElectronMatch[iPhoton][0]==2 && PhotonElectronMatch[iPhoton][2] == PhotonElectronMatch[iPhoton][3]){//merged conv photon
835 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
836 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
837 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 2);
838 if(esdCluster->GetNCells()>=2)
839 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 2);
840 }
841 else if(PhotonElectronMatch[iPhoton][0]==2){//non-merged conv photon (but both hit emcal)
842 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
843 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
844 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
845 if(esdCluster->GetNCells()>=2)
846 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
847 esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][3]);
848 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
849 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
850 if(esdCluster->GetNCells()>=2)
851 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
852 }
853 else if(PhotonElectronMatch[iPhoton][0]==1){//non-merged conv photon (one missed emcal)
854 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
855 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
856 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 3);
857 if(esdCluster->GetNCells()>=2)
858 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 3);
859 }
860 }//loop over decay photons (iPhoton).
861
862 }// 2 Photons pointed at the EMCAL!
863 }//for(nTracksMC) ie. Truth Pion loop.
864
865 }//if(isMC)
866
867 //######################### ~~~~~~~~~~~~ ##################################
868 //######################### DONE WITH MC ##################################
869 //######################### ~~~~~~~~~~~~ ##################################
870
871
872
873 // NEW HISTO should be filled before this point, as PostData puts the
874 // information for this iteration of the UserExec in the container
875 PostData(1, fOutput);
876 }
877
878//________________________________________________________________________
879void AliAnalysisTaskSDMGammaMC::Terminate(Option_t *) //specify what you want to have done
880{
881 // Called once at the end of the query.
882
883}
884
885//________________________________________________________________________
886Int_t AliAnalysisTaskSDMGammaMC::GetZvtxBin(Double_t vertZ)
887{
888
889 int izvtx = -1;
890
891 if (vertZ<-35)
892 izvtx=0;
893 else if(vertZ<-30)
894 izvtx=1;
895 else if(vertZ<-25)
896 izvtx=2;
897 else if(vertZ<-20)
898 izvtx=3;
899 else if(vertZ<-15)
900 izvtx=4;
901 else if(vertZ<-10)
902 izvtx=5;
903 else if(vertZ< -5)
904 izvtx=6;
905 else if(vertZ< 0)
906 izvtx=7;
907 else if(vertZ< 5)
908 izvtx=8;
909 else if(vertZ< 10)
910 izvtx=9;
911 else if(vertZ< 15)
912 izvtx=10;
913 else if(vertZ< 20)
914 izvtx=11;
915 else if(vertZ< 25)
916 izvtx=12;
917 else if(vertZ< 30)
918 izvtx=13;
919 else if(vertZ< 35)
920 izvtx=14;
921 else
922 izvtx=15;
923
924 return izvtx;
925}
926
927//________________________________________________________________________
928Int_t AliAnalysisTaskSDMGammaMC::GetMultBin(Int_t mult){
929
930 int imult = -1;
931
932 if (mult<2)
933 imult=0;
934 else if(mult<25)
935 imult=mult-2;
936 else
937 imult=24;
938
939 return imult;
940}
941
942//________________________________________________________________________
943Int_t AliAnalysisTaskSDMGammaMC::isGoodEsdCluster(AliESDCaloCluster* esdclust){
944
945 int pass = 1;
946 int nMinCells = 1;
947 double MinE = 0.4;
948 //double MinErat = 0;
949 //double MinEcc = 0;
950
951 if (!esdclust)
952 pass = 0;
953 if (!esdclust->IsEMCAL())
954 pass = 0;//removes ~70% of clusters.
955 if (esdclust->E()<MinE)
956 pass = 0;//does nothing
957 if (esdclust->GetNCells()<nMinCells)
958 pass = 0;//does nothing
959 //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
960 //pass = 0;
961 //if (esdclust->Chi2()<MinEcc) // eccentricity cut
962 //pass = 0;//this is always -1.
963
964 /*
965 //This cuts out more than just 1 cell clusters
966 //and drains the statistics badly.
967 //haven't figured out what it does yet.
968 if(esdclust->GetM20()<0.02)
969 pass = 0;
970 */
971 //if(esdclust->GetM02()<0.1)
972 // pass = 0;
973 //if(esdclust->GetM02()>0.5)
974 // pass = 0;
975 //if(esdclust->GetNCells()<2)
976 // pass = 0;
977
978 Float_t pos[3] = {0,0,0};
979 esdclust->GetPosition(pos);
980 TVector3 clusterPosition(pos);
981 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
982 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
983 pass = 0;
984 clusterPosition.Delete();
985
986 //doing this by hand now...
987 //if(!esdclust->GetNTracksMatched()==0)
988 //pass = 0;
989
990 return pass;
991}
992
993//________________________________________________________________________
994Int_t AliAnalysisTaskSDMGammaMC::isGoodAodCluster(AliAODCaloCluster* aodclust){
995
996 int pass = 1;
997 int nMinCells = 1;
998 double MinE = 0.4;
999 //double MinErat = 0;
1000 //double MinEcc = 0;
1001
1002 if (!aodclust)
1003 pass = 0;
1004 if (!aodclust->IsEMCAL())
1005 pass = 0;//removes ~70% of clusters.
1006 if (aodclust->E()<MinE)
1007 pass = 0;//does nothing
1008 if (aodclust->GetNCells()<nMinCells)
1009 pass = 0;//does nothing
1010 //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1011 //pass = 0;
1012 //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1013 //pass = 0;//this is always -1.
1014
1015 /*
1016 //This cuts out more than just 1 cell clusters
1017 //and drains the statistics badly.
1018 //haven't figured out what it does yet.
1019 if(aodclust->GetM20()<0.02)
1020 pass = 0;
1021 if(aodclust->GetM02()<0.02)
1022 pass = 0;
1023 */
1024 //if(aodclust->GetM02()<0.1)
1025 // pass = 0;
1026 //if(aodclust->GetM02()>0.5)
1027 // pass = 0;
1028 //if(aodclust->GetNCells()<2)
1029 // pass = 0;
1030
1031 Float_t pos[3] = {0,0,0};
1032 aodclust->GetPosition(pos);
1033 TVector3 clusterPosition(pos);
1034 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1035 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1036 pass = 0;
1037 clusterPosition.Delete();
1038
1039 //if(!aodclust->GetNTracksMatched()==0)
1040 //pass = 0;
1041
1042 return pass;
1043}
1044
1045//________________________________________________________________________
1046Double_t AliAnalysisTaskSDMGammaMC::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1047
1048 double dphi = p1.Phi() - p2.Phi();
1049
1050 if(dphi<0.5*TMath::Pi())
1051 dphi = dphi + 2.0*TMath::Pi();
1052
1053 if(dphi>1.5*TMath::Pi())
1054 dphi = dphi - 2.0*TMath::Pi();
1055
1056 return dphi;
1057}
1058
1059//________________________________________________________________________
1060Double_t AliAnalysisTaskSDMGammaMC::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1061
1062 double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1063
1064 return deta;
1065}
1066
1067
1068//________________________________________________________________________
1069Double_t AliAnalysisTaskSDMGammaMC::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1070
1071 double recalibfactor = 0.0;
1072
1073 if(iCalib==0){// no recalibration!
1074 recalibfactor = 1.0;
1075 }
1076 else if(iCalib==1){// just a scale factor:
1077 recalibfactor = 0.984;
1078 }
1079 else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1080 Double_t p[3] = {0.96968, -2.68720, -0.831607};
1081 recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1082 }
1083 else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1084 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};
1085 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1086 }
1087 else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1088 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};
1089 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1090 }
1091 else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1092 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};
1093 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) );
1094 }
1095 else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1096 Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1097 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1098 }
1099 else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1100 Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1101 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1102 }
1103 else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1104 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};
1105 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));
1106 }
1107 else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1108 Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1109 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1110 }
1111 else if(iCalib==10){// Jason played with test beam data
1112 Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1113 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1114 }
1115 else if(iCalib==11){// Jason played with test beam MC
1116 Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1117 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1118 }
1119
1120 return recalibfactor;
1121}
1122