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