]>
Commit | Line | Data |
---|---|---|
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 | 68 | using std::cout; |
69 | using std::endl; | |
70 | ||
3a66195e | 71 | ClassImp(AliAnalysisTaskEMCALMesonGGSDM) |
72 | ||
73 | //________________________________________________________________________ | |
74 | AliAnalysisTaskEMCALMesonGGSDM::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 | //________________________________________________________________________ | |
166 | AliAnalysisTaskEMCALMesonGGSDM::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 | //________________________________________________________________________ | |
263 | AliAnalysisTaskEMCALMesonGGSDM::~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 | //________________________________________________________________________ | |
274 | void 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 | //________________________________________________________________________ | |
822 | void 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 | ||
3a66195e | 1427 | for(int i=0; i<nclusters; i++) { |
1428 | ||
1429 | AliESDCaloCluster* esdCluster=NULL; | |
1430 | AliAODCaloCluster* aodCluster=NULL; | |
1431 | if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster | |
1432 | else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster | |
1433 | if(!esdCluster && !aodCluster) { | |
1434 | AliError(Form("ERROR: Could not retrieve any (ESD or AOD) Cluster %d",i)); | |
1435 | continue; | |
1436 | } | |
1437 | ||
1438 | if(fEsdEv){ | |
1439 | ||
1440 | recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); | |
1441 | ||
1442 | //uncomment this to do the track matching (1 of 3 lines, esd part)!! | |
1443 | //Bool_t MatchesToTrack = 0; | |
1444 | if(esdCluster->IsEMCAL()){ | |
1445 | ||
1446 | Float_t pos[3] = {0,0,0}; | |
1447 | Short_t maxCellID = -1; | |
1448 | Float_t celleta, cellphi; | |
1449 | esdCluster->GetPosition(pos); | |
1450 | TVector3 clusterPosition(pos); | |
1451 | h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta()); | |
1452 | GetMaxCellEnergy(esdCluster, maxCellID); | |
1453 | AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance(); | |
1454 | fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi); | |
1455 | h2_PhiEtaMaxCell->Fill(cellphi,celleta); | |
1456 | ||
1457 | // _______________Track loop for reconstructed event_____________ | |
1458 | for(Int_t itrk = 0; itrk < fEsdEv->GetNumberOfTracks(); itrk++) { | |
1459 | AliESDtrack* esdTrack = fEsdEv->GetTrack(itrk); // pointer to reconstructed to track | |
1460 | if(!esdTrack) { | |
1461 | AliError(Form("ERROR: Could not retrieve any (ESD) track %d",itrk)); | |
1462 | continue; | |
1463 | } | |
1464 | ||
1465 | Double_t posTrk[3] = {0,0,0}; | |
1466 | esdTrack->GetXYZ(posTrk); | |
1467 | TVector3 vposTrk(posTrk); | |
1468 | ||
1469 | Double_t fMass = 0.139; | |
1470 | Double_t fStepSurface = 20.; | |
1471 | Float_t etaproj, phiproj, pttrackproj; | |
1472 | ||
1473 | AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); | |
1474 | if(!trackParam) continue; | |
1475 | AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj); | |
1476 | ||
1477 | double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) + | |
1478 | (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) ); | |
1479 | ||
1480 | h1_dR_ClustTrk->Fill(dR_clusttrk); | |
1481 | ||
1482 | //uncomment this to do the track matching (2 of 3 lines, esd part)!! | |
1483 | //if(dR_clusttrk<fdRmin_ClustTrack) | |
1484 | //MatchesToTrack = 1; | |
1485 | ||
1486 | }//_____________________________nTracks__________________________ | |
1487 | ||
1488 | h2_cells_M02 ->Fill(esdCluster->GetNCells(),esdCluster->GetM02()); | |
1489 | h2_Ellipse ->Fill(esdCluster->GetM20(),esdCluster->GetM02()); | |
1490 | h1_Chi2 ->Fill(esdCluster->Chi2());//always -1. | |
1491 | h1_nTrkMatch ->Fill(esdCluster->GetNTracksMatched()); | |
1492 | h1_ClusterDisp->Fill(esdCluster->GetDispersion()); | |
1493 | h2_E_time ->Fill(esdCluster->E(),esdCluster->GetTOF()); | |
1494 | ||
1495 | TArrayI *TrackLabels = esdCluster->GetTracksMatched(); | |
1496 | if(TrackLabels){ | |
1497 | if(TrackLabels->GetSize()>0){ | |
1498 | Int_t trackindex = TrackLabels->At(0); | |
1499 | AliESDtrack* matchingT = fEsdEv->GetTrack(trackindex); // pointer to reconstructed to track | |
1500 | ||
1501 | recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); | |
1502 | h2_eop_E ->Fill(esdCluster->E()*recalScale, esdCluster->E()*recalScale/matchingT->P()); | |
1503 | h2_eop_pT->Fill(matchingT->Pt(), esdCluster->E()*recalScale/matchingT->P()); | |
1504 | } | |
1505 | } | |
1506 | ||
1507 | //uncomment this to do the track matching (2 of 3 lines, esd part)!! | |
1508 | //if(isGoodEsdCluster(esdCluster) && !MatchesToTrack){ | |
cfd87ccd | 1509 | |
3a66195e | 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); | |
cfd87ccd | 1522 | |
1523 | ||
3a66195e | 1524 | } |
1525 | clusterPosition.Delete(); | |
1526 | }//if(esdCluster->isEMCAL()) | |
1527 | }//if(fEsdEv) | |
1528 | else if(fAodEv){ | |
1529 | ||
1530 | recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator); | |
1531 | ||
1532 | //uncomment this to do the track matching (1 of 3 lines, aod part)!! | |
1533 | //Bool_t MatchesToTrack = 0; | |
1534 | if(aodCluster->IsEMCAL()){ | |
1535 | ||
1536 | Float_t pos[3] = {0,0,0}; | |
1537 | Short_t maxCellID = -1; | |
1538 | Float_t celleta, cellphi; | |
1539 | aodCluster->GetPosition(pos); | |
1540 | TVector3 clusterPosition(pos); | |
1541 | h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta()); | |
1542 | GetMaxCellEnergy(aodCluster, maxCellID); | |
1543 | AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance(); | |
1544 | fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi); | |
1545 | h2_PhiEtaMaxCell->Fill(cellphi,celleta); | |
1546 | ||
1547 | // _______________Track loop for reconstructed event_____________ | |
1548 | for(Int_t itrk = 0; itrk < fAodEv->GetNumberOfTracks(); itrk++) { | |
f15c1f69 | 1549 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fAodEv->GetTrack(itrk)); |
1550 | if(!aodTrack) AliFatal("Not a standard AOD"); // pointer to reconstructed to track | |
3a66195e | 1551 | if(!aodTrack) { |
1552 | AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk)); | |
1553 | continue; | |
1554 | } | |
1555 | ||
1556 | Double_t posTrk[3] = {0,0,0}; | |
1557 | Double_t momTrk[3] = {0,0,0}; | |
1558 | aodTrack->GetXYZ(posTrk); | |
1559 | aodTrack->GetPxPyPz(momTrk); | |
1560 | //TVector3 vposTrk(posTrk); | |
1561 | ||
1562 | //#################################################################################################### | |
1563 | // | |
1564 | // commented all this stuff just to satisfy aliroot warnings. | |
1565 | // but I may need it again if I want to do the track matching for aods. | |
1566 | /* | |
1567 | Double_t fMass = 0.139; | |
1568 | Double_t fStepSurface = 20.; | |
1569 | Float_t etaproj=0.0; | |
1570 | Float_t phiproj=0.0; | |
1571 | Float_t pttrackproj=0.0; | |
1572 | ||
1573 | Double_t cv[21] = {0.0}; | |
1574 | aodTrack->GetCovarianceXYZPxPyPz(cv); | |
1575 | AliExternalTrackParam *trackParam = new AliExternalTrackParam(posTrk,momTrk,cv,aodTrack->Charge()); | |
1576 | //AliExternalTrackParam emcalParam(*trackParam); | |
1577 | //AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam()); | |
1578 | if(!trackParam) continue; | |
1579 | ////AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj); | |
1580 | //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&emcalParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj); | |
1581 | delete trackParam; | |
1582 | ||
1583 | //Constantin's implementation... gives funny result. | |
1584 | //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(aodTrack,440.0); | |
1585 | //phiproj = aodTrack->GetTrackPhiOnEMCal(); | |
1586 | //etaproj = aodTrack->GetTrackPhiOnEMCal(); | |
1587 | ||
1588 | double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) + | |
1589 | (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) ); | |
1590 | ||
1591 | h1_dR_ClustTrk->Fill(dR_clusttrk); | |
1592 | */ | |
1593 | //#################################################################################################### | |
1594 | ||
1595 | //uncomment this to do the track matching (2 of 3 lines, aod part)!! | |
1596 | //if(dR_clusttrk<fdRmin_ClustTrack) | |
1597 | //MatchesToTrack = 1; | |
1598 | ||
1599 | ||
1600 | }//_____________________________nTracks__________________________ | |
1601 | ||
1602 | h2_cells_M02 ->Fill(aodCluster->GetNCells(),aodCluster->GetM02()); | |
1603 | h2_Ellipse ->Fill(aodCluster->GetM20(),aodCluster->GetM02()); | |
1604 | h1_Chi2 ->Fill(aodCluster->Chi2());//always -1. | |
1605 | h1_nTrkMatch ->Fill(aodCluster->GetNTracksMatched()); | |
1606 | h1_ClusterDisp->Fill(aodCluster->GetDispersion()); | |
1607 | h2_E_time ->Fill(aodCluster->E(),aodCluster->GetTOF()); | |
1608 | ||
1609 | // ################################################# | |
1610 | // track matching eop histograms are handled here... | |
1611 | // ################################################# | |
1612 | ||
1613 | //uncomment this to do the track matching (3 of 3 lines, aod part)!! | |
1614 | //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){ | |
1615 | if(isGoodAodCluster(aodCluster)){ | |
1616 | recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator); | |
1617 | E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ | |
1618 | fAodEv->GetVertex(0)->GetXYZ(vertex); | |
1619 | aodCluster->GetMomentum(Photon1,vertex); | |
1620 | Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ | |
1621 | Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ | |
1622 | Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ | |
1623 | Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) ); | |
1624 | h1_E->Fill(E1); | |
1625 | h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta()); | |
1626 | h2_PhiEtaMaxCellCut->Fill(cellphi,celleta); | |
1627 | } | |
1628 | clusterPosition.Delete(); | |
1629 | }//if(aodCluster->IsEMCAL()) | |
1630 | }//if(fAodEv) | |
1631 | ||
1632 | }//loop over nclusters. | |
cfd87ccd | 1633 | |
3a66195e | 1634 | //Make same event pions... |
1635 | for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){ | |
1636 | for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){ | |
1637 | Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j]; | |
1638 | Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]); | |
1639 | Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]); | |
1640 | Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/ | |
1641 | (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt()); | |
1642 | Int_t asymCut = 0; | |
1643 | if (pairasym<0.1) asymCut = 1; | |
1644 | else if(pairasym<0.7) asymCut = 2; | |
1645 | else asymCut = 3; | |
cfd87ccd | 1646 | |
3a66195e | 1647 | h1_M ->Fill(Parent.M()); |
1648 | h3_MptAsymm ->Fill(Parent.M(),Parent.Pt(),asymCut); | |
1649 | h2_dphi_deta->Fill(deltaphi,deltaeta); | |
1650 | } | |
1651 | } | |
1652 | ||
1653 | //Make mixed event... | |
1654 | for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){ | |
1655 | for(unsigned int ipool=1; ipool<poolDepth; ipool++){ | |
1656 | for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){ | |
1657 | iskip = randy.Integer(Photons[0][izvtx][imult].size()); | |
1658 | if(j==iskip) continue; | |
1659 | Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j]; | |
1660 | Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]); | |
1661 | Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]); | |
1662 | Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/ | |
1663 | (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt()); | |
1664 | Int_t asymCut = 0; | |
cfd87ccd | 1665 | |
3a66195e | 1666 | if (pairasym<0.1) asymCut = 1; |
1667 | else if(pairasym<0.7) asymCut = 2; | |
1668 | else asymCut = 3; | |
1669 | ||
1670 | h1_M_mix ->Fill(Parent.M()); | |
1671 | h3_MptAsymm_mix ->Fill(Parent.M(),Parent.Pt(),asymCut); | |
1672 | h2_dphi_deta_mix->Fill(deltaphi,deltaeta); | |
1673 | } | |
1674 | } | |
1675 | } | |
1676 | ||
1677 | for(int ipool=poolDepth-1; ipool>0; ipool--){ | |
1678 | Photons[ipool][izvtx][imult].clear(); | |
1679 | for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++) | |
1680 | Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]); | |
1681 | } | |
1682 | Photons[0][izvtx][imult].clear(); | |
1683 | ||
1684 | ||
1685 | ||
1686 | // NEW HISTO should be filled before this point, as PostData puts the | |
1687 | // information for this iteration of the UserExec in the container | |
1688 | PostData(1, fOutput); | |
1689 | } | |
1690 | ||
1691 | //________________________________________________________________________ | |
1692 | void AliAnalysisTaskEMCALMesonGGSDM::Terminate(Option_t *) //specify what you want to have done | |
1693 | { | |
1694 | // Called once at the end of the query. | |
1695 | ||
1696 | } | |
1697 | ||
cfd87ccd | 1698 | |
1699 | // //________________________________________________________________________ | |
1700 | // Int_t AliAnalysisTaskEMCALMesonGGSDM::GetZvtxBin(Double_t vertZ) | |
1701 | // { | |
1702 | // | |
1703 | // int izvtx = -1; | |
1704 | // | |
1705 | // if (vertZ<-35) | |
1706 | // izvtx=0; | |
1707 | // else if(vertZ<-30) | |
1708 | // izvtx=1; | |
1709 | // else if(vertZ<-25) | |
1710 | // izvtx=2; | |
1711 | // else if(vertZ<-20) | |
1712 | // izvtx=3; | |
1713 | // else if(vertZ<-15) | |
1714 | // izvtx=4; | |
1715 | // else if(vertZ<-10) | |
1716 | // izvtx=5; | |
1717 | // else if(vertZ< -5) | |
1718 | // izvtx=6; | |
1719 | // else if(vertZ< 0) | |
1720 | // izvtx=7; | |
1721 | // else if(vertZ< 5) | |
1722 | // izvtx=8; | |
1723 | // else if(vertZ< 10) | |
1724 | // izvtx=9; | |
1725 | // else if(vertZ< 15) | |
1726 | // izvtx=10; | |
1727 | // else if(vertZ< 20) | |
1728 | // izvtx=11; | |
1729 | // else if(vertZ< 25) | |
1730 | // izvtx=12; | |
1731 | // else if(vertZ< 30) | |
1732 | // izvtx=13; | |
1733 | // else if(vertZ< 35) | |
1734 | // izvtx=14; | |
1735 | // else | |
1736 | // izvtx=15; | |
1737 | // | |
1738 | // return izvtx; | |
1739 | // } | |
1740 | ||
3a66195e | 1741 | //________________________________________________________________________ |
1742 | Int_t AliAnalysisTaskEMCALMesonGGSDM::GetZvtxBin(Double_t vertZ) | |
1743 | { | |
1744 | ||
1745 | int izvtx = -1; | |
1746 | ||
cfd87ccd | 1747 | if (vertZ<-3.375) |
3a66195e | 1748 | izvtx=0; |
cfd87ccd | 1749 | else if(vertZ<-1.605) |
3a66195e | 1750 | izvtx=1; |
cfd87ccd | 1751 | else if(vertZ<-0.225) |
3a66195e | 1752 | izvtx=2; |
cfd87ccd | 1753 | else if(vertZ<1.065) |
3a66195e | 1754 | izvtx=3; |
cfd87ccd | 1755 | else if(vertZ<-2.445) |
3a66195e | 1756 | izvtx=4; |
cfd87ccd | 1757 | else if(vertZ<-4.245) |
3a66195e | 1758 | izvtx=5; |
3a66195e | 1759 | else |
cfd87ccd | 1760 | izvtx=6; |
3a66195e | 1761 | |
1762 | return izvtx; | |
1763 | } | |
1764 | ||
cfd87ccd | 1765 | |
1766 | // //________________________________________________________________________ | |
1767 | // Int_t AliAnalysisTaskEMCALMesonGGSDM::GetMultBin(Int_t mult){ | |
1768 | // | |
1769 | // int imult = -1; | |
1770 | // | |
1771 | // if (mult<2) | |
1772 | // imult=0; | |
1773 | // else if(mult<25) | |
1774 | // imult=mult-2; | |
1775 | // else | |
1776 | // imult=24; | |
1777 | // | |
1778 | // return imult; | |
1779 | // } | |
1780 | ||
3a66195e | 1781 | //________________________________________________________________________ |
1782 | Int_t AliAnalysisTaskEMCALMesonGGSDM::GetMultBin(Int_t mult){ | |
1783 | ||
1784 | int imult = -1; | |
1785 | ||
1786 | if (mult<2) | |
1787 | imult=0; | |
cfd87ccd | 1788 | else if(mult<3) |
1789 | imult=1; | |
1790 | else if(mult<4) | |
1791 | imult=2; | |
1792 | else if(mult<8) | |
1793 | imult=3; | |
1794 | else if(mult<15) | |
1795 | imult=4; | |
3a66195e | 1796 | else |
cfd87ccd | 1797 | imult=5; |
3a66195e | 1798 | |
1799 | return imult; | |
1800 | } | |
1801 | ||
cfd87ccd | 1802 | |
3a66195e | 1803 | //________________________________________________________________________ |
1804 | Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodEsdCluster(AliESDCaloCluster* esdclust){ | |
1805 | ||
1806 | int pass = 1; | |
1807 | int nMinCells = 2; | |
1808 | double MinE = 0.4; | |
1809 | //double MinErat = 0; | |
1810 | //double MinEcc = 0; | |
1811 | ||
1812 | if (!esdclust) | |
1813 | pass = 0; | |
1814 | if (!esdclust->IsEMCAL()) | |
1815 | pass = 0; | |
1816 | if (esdclust->E()<MinE) | |
1817 | pass = 0; | |
1818 | if (esdclust->GetNCells()<nMinCells) | |
1819 | pass = 0; | |
1820 | //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat) | |
1821 | //pass = 0; | |
1822 | //if (esdclust->Chi2()<MinEcc) // eccentricity cut | |
1823 | //pass = 0;//this is always -1. | |
1824 | ||
1825 | //if(esdclust->GetM02()<0.1) | |
1826 | // pass = 0; | |
1827 | //if(esdclust->GetM02()>0.5) | |
1828 | // pass = 0; | |
1829 | ||
1830 | Float_t pos[3] = {0,0,0}; | |
1831 | esdclust->GetPosition(pos); | |
1832 | TVector3 clusterPosition(pos); | |
1833 | if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax || | |
1834 | clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax ) | |
1835 | pass = 0; | |
1836 | clusterPosition.Delete(); | |
1837 | ||
1838 | //DOING THIS BY HAND NOW... | |
1839 | //if(!esdclust->GetNTracksMatched()==0) | |
1840 | //pass = 0; | |
1841 | ||
1842 | return pass; | |
1843 | } | |
1844 | ||
1845 | //________________________________________________________________________ | |
1846 | Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodAodCluster(AliAODCaloCluster* aodclust){ | |
1847 | ||
1848 | int pass = 1; | |
1849 | int nMinCells = 2; | |
1850 | double MinE = 0.4; | |
1851 | //double MinErat = 0; | |
1852 | //double MinEcc = 0; | |
1853 | ||
1854 | if (!aodclust) | |
1855 | pass = 0; | |
1856 | if (!aodclust->IsEMCAL()) | |
1857 | pass = 0; | |
1858 | if (aodclust->E()<MinE) | |
1859 | pass = 0; | |
1860 | if (aodclust->GetNCells()<nMinCells) | |
1861 | pass = 0; | |
1862 | //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat) | |
1863 | //pass = 0; | |
1864 | //if (aodclust->Chi2()<MinEcc) // eccentricity cut | |
1865 | //pass = 0;//this is always -1. | |
1866 | ||
1867 | //if(aodclust->GetM02()<0.1) | |
1868 | //pass = 0; | |
1869 | //if(aodclust->GetM02()>0.5) | |
1870 | //pass = 0; | |
1871 | ||
1872 | Float_t pos[3] = {0,0,0}; | |
1873 | aodclust->GetPosition(pos); | |
1874 | TVector3 clusterPosition(pos); | |
1875 | if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax || | |
1876 | clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax ) | |
1877 | pass = 0; | |
1878 | clusterPosition.Delete(); | |
1879 | ||
1880 | //DOING THIS BY HAND NOW... | |
1881 | //if(!aodclust->GetNTracksMatched()==0) | |
1882 | //pass = 0; | |
1883 | ||
1884 | return pass; | |
1885 | } | |
1886 | ||
1887 | //________________________________________________________________________ | |
1888 | Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){ | |
1889 | ||
1890 | double dphi = p1.Phi() - p2.Phi(); | |
1891 | ||
1892 | if(dphi<0.5*TMath::Pi()) | |
1893 | dphi = dphi + 2.0*TMath::Pi(); | |
1894 | ||
1895 | if(dphi>1.5*TMath::Pi()) | |
1896 | dphi = dphi - 2.0*TMath::Pi(); | |
1897 | ||
1898 | return dphi; | |
1899 | } | |
1900 | ||
1901 | //________________________________________________________________________ | |
1902 | Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaEta(TLorentzVector p1, TLorentzVector p2){ | |
1903 | ||
1904 | double deta = p1.PseudoRapidity() - p2.PseudoRapidity(); | |
1905 | ||
1906 | return deta; | |
1907 | } | |
1908 | ||
1909 | ||
1910 | //________________________________________________________________________ | |
1911 | Double_t AliAnalysisTaskEMCALMesonGGSDM::PrivateEnergyRecal(Double_t energy, Int_t iCalib){ | |
1912 | ||
1913 | double recalibfactor = 0.0; | |
1914 | ||
1915 | if(iCalib==0){// no recalibration! | |
1916 | recalibfactor = 1.0; | |
1917 | } | |
1918 | else if(iCalib==1){// just a scale factor: | |
1919 | recalibfactor = 0.984; | |
1920 | } | |
1921 | else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC. | |
1922 | Double_t p[3] = {0.96968, -2.68720, -0.831607}; | |
1923 | recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0); | |
1924 | } | |
1925 | else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??) | |
1926 | 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}; | |
1927 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1928 | } | |
1929 | else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??) | |
1930 | 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}; | |
1931 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1932 | } | |
1933 | else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??) | |
1934 | 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}; | |
1935 | 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) ); | |
1936 | } | |
1937 | else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT! | |
1938 | Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968}; | |
1939 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1940 | } | |
1941 | else if(iCalib==7){// kPi0MCv3 - in AliROOT! | |
1942 | Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0}; | |
1943 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1944 | } | |
1945 | else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??) | |
1946 | 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}; | |
1947 | //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 | |
1948 | 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)); | |
1949 | } | |
1950 | else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??) | |
1951 | Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286}; | |
1952 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1953 | } | |
1954 | else if(iCalib==10){// Jason played with test beam data | |
1955 | Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054}; | |
1956 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1957 | } | |
1958 | else if(iCalib==11){// Jason played with test beam MC | |
1959 | Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437}; | |
1960 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); | |
1961 | } | |
1962 | ||
1963 | return recalibfactor; | |
1964 | } | |
1965 | ||
1966 | ||
1967 | //________________________________________________________________________ | |
1968 | Double_t AliAnalysisTaskEMCALMesonGGSDM::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const | |
1969 | { | |
1970 | // Get maximum energy of attached cell. | |
1971 | ||
1972 | id = -1; | |
1973 | AliVCaloCells *fVCells=NULL; | |
1974 | if(fEsdEv) fVCells = fEsdEv->GetEMCALCells(); | |
1975 | else if(fAodEv) fVCells = fAodEv->GetEMCALCells(); | |
1976 | if(!fVCells) | |
1977 | return 0; | |
1978 | ||
1979 | Double_t maxe = 0; | |
1980 | Int_t ncells = cluster->GetNCells(); | |
1981 | for (Int_t i=0; i<ncells; i++) { | |
1982 | Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); | |
1983 | if (e>maxe) { | |
1984 | maxe = e; | |
1985 | id = cluster->GetCellAbsId(i); | |
1986 | } | |
1987 | } | |
1988 | return maxe; | |
1989 | } | |
1990 | ||
1991 | ||
1992 | //________________________________________________________________________ | |
1993 | Int_t AliAnalysisTaskEMCALMesonGGSDM::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){ | |
1994 | ||
1995 | AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack)); | |
1996 | ||
1997 | Int_t nPTracks= mcEvent->GetNumberOfPrimaries(); | |
1998 | ||
1999 | Int_t isPhysPrimary = 1; | |
2000 | Int_t ismHF = 0; | |
2001 | Int_t ismLongLivedOrK = 0; | |
2002 | ||
2003 | if(mcP->GetMother()<0)//if it has no mother... | |
2004 | return isPhysPrimary; | |
2005 | ||
2006 | Int_t imTrack = mcP->GetMother(); | |
2007 | AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack)); | |
2008 | ||
2009 | if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark... | |
2010 | return isPhysPrimary; | |
2011 | ||
2012 | ||
2013 | //############################################ | |
2014 | //get the PDG digits.... | |
2015 | int num = mcPm->PdgCode(); | |
2016 | int RevDigits[10] = {0}; | |
2017 | int nDigits = 0; | |
2018 | while (num >= 1){ | |
2019 | RevDigits[nDigits++] = num%10; | |
2020 | num = num / 10; | |
2021 | } | |
2022 | //############################################## | |
2023 | ||
2024 | ||
2025 | if(RevDigits[3]>3)//Baryons | |
2026 | ismHF = 1; | |
2027 | else if(RevDigits[2]>3)//Mesons | |
2028 | ismHF = 1; | |
2029 | ||
2030 | ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode()); | |
2031 | ||
2032 | if(!ismHF && ismLongLivedOrK) | |
2033 | isPhysPrimary = 0; | |
2034 | else{ // check grandmother, greatgrandmothers, etc... | |
2035 | while(imTrack >= nPTracks){ | |
2036 | ||
2037 | if(mcPm->GetMother()<0)//if it has no mother... | |
2038 | break; | |
2039 | ||
2040 | if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark... | |
2041 | return isPhysPrimary; | |
2042 | ||
2043 | imTrack = mcPm->GetMother(); | |
2044 | mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack)); | |
2045 | ||
2046 | //############################################ | |
2047 | //get the PDG digits.... | |
2048 | num = mcPm->PdgCode(); | |
2049 | for(int i=0; i<10; i++) RevDigits[i] = 0; | |
2050 | nDigits = 0; | |
2051 | while (num >= 1){ | |
2052 | RevDigits[nDigits++] = num%10; | |
2053 | num = num / 10; | |
2054 | } | |
2055 | //############################################## | |
2056 | if(RevDigits[3]>3)//Baryons | |
2057 | ismHF = 1; | |
2058 | else if(RevDigits[2]>3)//Mesons | |
2059 | ismHF = 1; | |
2060 | ||
2061 | ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode()); | |
2062 | ||
2063 | if(!ismHF && ismLongLivedOrK) | |
2064 | isPhysPrimary = 0; | |
2065 | ||
2066 | }//while( >=nPTracks) | |
2067 | } | |
2068 | ||
2069 | return isPhysPrimary; | |
2070 | } | |
2071 | ||
2072 | ||
2073 | //________________________________________________________________________ | |
2074 | Int_t AliAnalysisTaskEMCALMesonGGSDM::IsLongLivedOrK(Int_t MyPDGcode){ | |
2075 | ||
2076 | Int_t MyFlag = 0; | |
2077 | ||
2078 | if( | |
2079 | (TMath::Abs(MyPDGcode) == 22 ) || // Photon | |
2080 | (TMath::Abs(MyPDGcode) == 11 ) || // Electron | |
2081 | (TMath::Abs(MyPDGcode) == 13 ) || // Muon(-) | |
2082 | (TMath::Abs(MyPDGcode) == 211 ) || // Pion | |
2083 | (TMath::Abs(MyPDGcode) == 321 ) || // Kaon | |
2084 | (TMath::Abs(MyPDGcode) == 310 ) || // K0s | |
2085 | (TMath::Abs(MyPDGcode) == 130 ) || // K0l | |
2086 | (TMath::Abs(MyPDGcode) == 2212) || // Proton | |
2087 | (TMath::Abs(MyPDGcode) == 2112) || // Neutron | |
2088 | (TMath::Abs(MyPDGcode) == 3122) || // Lambda_0 | |
2089 | (TMath::Abs(MyPDGcode) == 3112) || // Sigma Minus | |
2090 | (TMath::Abs(MyPDGcode) == 3222) || // Sigma Plus | |
2091 | (TMath::Abs(MyPDGcode) == 3312) || // Xsi Minus | |
2092 | (TMath::Abs(MyPDGcode) == 3322) || // Xsi | |
2093 | (TMath::Abs(MyPDGcode) == 3334) || // Omega | |
2094 | (TMath::Abs(MyPDGcode) == 12 ) || // Electron Neutrino | |
2095 | (TMath::Abs(MyPDGcode) == 14 ) || // Muon Neutrino | |
2096 | (TMath::Abs(MyPDGcode) == 16 ) ) // Tau Neutrino | |
2097 | MyFlag = 1; | |
2098 | ||
2099 | return MyFlag; | |
2100 | } | |
2101 | ||
2102 | ||
2103 | //________________________________________________________________________ | |
2104 | Int_t AliAnalysisTaskEMCALMesonGGSDM::IsMyMCHeaderType(Int_t iTrack, char *MyType, AliMCEvent *mcEvent) const | |
2105 | { | |
2106 | ||
2107 | Int_t isMyType = 0; | |
2108 | ||
2109 | AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(mcEvent->GenEventHeader()); | |
2110 | if(!cocktail) | |
2111 | return 0; | |
2112 | ||
2113 | TList *genHeaders = cocktail->GetHeaders(); | |
2114 | ||
2115 | Int_t nGenerators = genHeaders->GetEntries(); | |
2116 | Int_t indexMyType = -1; | |
2117 | Int_t startParticle=0; | |
2118 | ||
2119 | for(Int_t igen = 0; igen < nGenerators; igen++){ | |
2120 | AliGenEventHeader* eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ; | |
2121 | TString name = eventHeader2->GetName(); | |
2122 | startParticle += eventHeader2->NProduced(); | |
2123 | //cout << name << endl; | |
2124 | if (name.Contains(MyType,TString::kIgnoreCase)){ | |
2125 | indexMyType = igen; | |
2126 | startParticle -= eventHeader2->NProduced(); | |
2127 | break; | |
2128 | } | |
2129 | } | |
2130 | ||
2131 | AliGenEventHeader *addedPi0Header = (AliGenEventHeader*)genHeaders->At(indexMyType); | |
2132 | Int_t ipi0min = startParticle; | |
2133 | Int_t ipi0max = ipi0min+addedPi0Header->NProduced()-1; | |
2134 | if(iTrack >= ipi0min && iTrack <= ipi0max) | |
2135 | isMyType = 1; | |
2136 | ||
2137 | return isMyType; | |
2138 | } | |
2139 | ||
2140 |