]>
Commit | Line | Data |
---|---|---|
6d08aa8d | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | /* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */ | |
16 | ||
17 | //_________________________________________________________________________ | |
18 | // | |
19 | // Class cloned from AliAnaPhoton, main aim is shower shape studies | |
20 | // | |
21 | // | |
22 | // | |
23 | //-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS) | |
24 | ////////////////////////////////////////////////////////////////////////////// | |
25 | ||
26 | ||
27 | // --- ROOT system --- | |
28 | #include <TH3F.h> | |
29 | #include <TH2F.h> | |
30 | #include <TClonesArray.h> | |
31 | //#include <TObjString.h> | |
32 | #include <Riostream.h> | |
33 | #include "TParticle.h" | |
34 | //#include <fstream> | |
35 | ||
36 | // --- Analysis system --- | |
37 | #include "AliAnaShowerParameter.h" | |
38 | #include "AliCaloTrackReader.h" | |
39 | #include "AliStack.h" | |
40 | #include "AliCaloPID.h" | |
41 | #include "AliMCAnalysisUtils.h" | |
42 | #include "AliFiducialCut.h" | |
43 | #include "AliAODCaloCluster.h" | |
44 | #include "AliAODMCParticle.h" | |
45 | #include "AliAnaShowerParameter.h" | |
46 | #include "AliEMCALGeoUtils.h" | |
47 | #include "AliAODCaloCells.h" | |
48 | #include "AliAODEvent.h" | |
49 | ||
50 | ||
51 | ClassImp(AliAnaShowerParameter) | |
52 | ||
53 | //____________________________________________________________________________ | |
54 | AliAnaShowerParameter::AliAnaShowerParameter() : | |
55 | AliAnaPartCorrBaseClass(), fCalorimeter(""), | |
56 | fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0), | |
57 | fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0), fNCellsCut(0), | |
58 | fTimeCutMin(-1), fTimeCutMax(9999999),fEMCALGeo(0), | |
59 | fhNClusters(0),fhNCellCluster(0),fhPtCluster(0),fhPhiCluster(0),fhEtaCluster(0), | |
60 | fhLambdaCluster(0),fhDispersionCluster(0),fhELambdaCluster(0),fhELambdaCellCluster(0), | |
61 | fhNCellPhoton(0),fhLambdaPhoton(0),fhDispersionPhoton(0),fhELambdaPhoton(0),fhELambdaCellPhoton(0), | |
62 | fhNCellPi0(0),fhLambdaPi0(0),fhDispersionPi0(0),fhELambdaPi0(0),fhELambdaCellPi0(0), | |
63 | fhNCellChargedHadron(0),fhLambdaChargedHadron(0),fhDispersionChargedHadron(0),fhELambdaChargedHadron(0),fhELambdaCellChargedHadron(0), | |
64 | //MC | |
65 | fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0), | |
66 | fhEMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),fhLambdaMCPhoton(0),fhPhotTrueE(0),fhRatioEPhoton(0), | |
67 | fhEMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0),fhLambdaMCPi0(0),fhPi0TrueE(0), | |
68 | fhProductionDistance(0),fhProductionRadius(0),fhDecayAngle(0),fhRatioEPi0(0), | |
69 | fhEMCPion(0),fhPhiMCPion(0),fhEtaMCPion(0),fhLambdaMCPion(0),fhPionTrueE(0),fhRatioEPion(0), | |
70 | fhEMCProton(0),fhPhiMCProton(0),fhEtaMCProton(0),fhLambdaMCProton(0),fhProtonTrueE(0),fhRatioEProton(0), | |
71 | fhEMCAntiProton(0),fhPhiMCAntiProton(0),fhEtaMCAntiProton(0),fhLambdaMCAntiProton(0),fhAntiProtonTrueE(0),fhRatioEAntiProton(0), | |
72 | fhEMCNeutron(0),fhPhiMCNeutron(0),fhEtaMCNeutron(0),fhLambdaMCNeutron(0),fhNeutronTrueE(0),fhRatioENeutron(0), | |
73 | fhEMCEta(0),fhPhiMCEta(0),fhEtaMCEta(0),fhLambdaMCEta(0), | |
74 | fhPrimPt(0x0) | |
75 | { | |
76 | //default ctor | |
77 | ||
78 | //Initialize parameters | |
79 | InitParameters(); | |
80 | ||
81 | } | |
82 | ||
83 | //____________________________________________________________________________ | |
84 | AliAnaShowerParameter::~AliAnaShowerParameter() | |
85 | { | |
86 | //dtor | |
87 | ||
88 | } | |
89 | ||
90 | //________________________________________________________________________ | |
91 | TObjString * AliAnaShowerParameter::GetAnalysisCuts() | |
92 | { | |
93 | //Save parameters used for analysis | |
94 | TString parList ; //this will be list of parameters used for this analysis. | |
164a1d84 | 95 | const Int_t buffersize = 255; |
96 | char onePar[buffersize] ; | |
6d08aa8d | 97 | |
164a1d84 | 98 | snprintf(onePar,buffersize,"--- AliAnaShowerParameter ---\n") ; |
6d08aa8d | 99 | parList+=onePar ; |
164a1d84 | 100 | snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ; |
6d08aa8d | 101 | parList+=onePar ; |
164a1d84 | 102 | snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ; |
6d08aa8d | 103 | parList+=onePar ; |
164a1d84 | 104 | snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ; |
6d08aa8d | 105 | parList+=onePar ; |
164a1d84 | 106 | snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ; |
6d08aa8d | 107 | parList+=onePar ; |
164a1d84 | 108 | snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ; |
6d08aa8d | 109 | parList+=onePar ; |
110 | ||
111 | //Get parameters set in base class. | |
112 | parList += GetBaseParametersList() ; | |
113 | ||
114 | //Get parameters set in PID class. | |
115 | parList += GetCaloPID()->GetPIDParametersList() ; | |
116 | ||
117 | //Get parameters set in FiducialCut class (not available yet) | |
118 | //parlist += GetFidCut()->GetFidCutParametersList() | |
119 | ||
120 | return new TObjString(parList) ; | |
121 | } | |
122 | ||
123 | ||
124 | //________________________________________________________________________ | |
125 | TList * AliAnaShowerParameter::GetCreateOutputObjects() | |
126 | { | |
127 | // Create histograms to be saved in output file and | |
128 | // store them in outputContainer | |
129 | TList * outputContainer = new TList() ; | |
130 | outputContainer->SetName("PhotonHistos") ; | |
131 | ||
132 | Int_t nptbins = GetHistoPtBins(); | |
133 | Int_t nphibins = GetHistoPhiBins(); | |
134 | Int_t netabins = GetHistoEtaBins(); | |
135 | Float_t ptmax = GetHistoPtMax(); | |
136 | Float_t phimax = GetHistoPhiMax(); | |
137 | Float_t etamax = GetHistoEtaMax(); | |
138 | Float_t ptmin = GetHistoPtMin(); | |
139 | Float_t phimin = GetHistoPhiMin(); | |
140 | Float_t etamin = GetHistoEtaMin(); | |
141 | ||
142 | //General non-MC Cluster histograms | |
143 | fhNClusters = new TH1F ("hNClusters","NClusters",21,-0.5,20.5); | |
144 | fhNClusters->SetXTitle("N_{Clusters}"); | |
145 | outputContainer->Add(fhNClusters) ; | |
146 | ||
147 | fhNCellCluster = new TH2F ("hNCellCluster","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); | |
148 | fhNCellCluster->SetXTitle("p_{T}"); | |
149 | fhNCellCluster->SetYTitle("N_{Cells}"); | |
150 | outputContainer->Add(fhNCellCluster) ; | |
151 | ||
152 | fhPtCluster = new TH1F("hPtCluster","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); | |
153 | fhPtCluster->SetYTitle("N"); | |
154 | fhPtCluster->SetXTitle("p_{T #gamma}(GeV/c)"); | |
155 | outputContainer->Add(fhPtCluster) ; | |
156 | ||
157 | fhPhiCluster = new TH2F | |
158 | ("hPhiCluster","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
159 | fhPhiCluster->SetYTitle("#phi"); | |
160 | fhPhiCluster->SetXTitle("p_{T #gamma} (GeV/c)"); | |
161 | outputContainer->Add(fhPhiCluster) ; | |
162 | ||
163 | fhEtaCluster = new TH2F | |
164 | ("hEtaCluster","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
165 | fhEtaCluster->SetYTitle("#eta"); | |
166 | fhEtaCluster->SetXTitle("p_{T Cluster} (GeV/c)"); | |
167 | outputContainer->Add(fhEtaCluster) ; | |
168 | ||
169 | fhLambdaCluster = new TH3F | |
170 | ("hLambdaCluster","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
171 | fhLambdaCluster->SetZTitle("#lambda_{1}^{2}"); | |
172 | fhLambdaCluster->SetYTitle("#lambda_{0}^{2}"); | |
173 | fhLambdaCluster->SetXTitle("p_{T Cluster} (GeV/c)"); | |
174 | outputContainer->Add(fhLambdaCluster) ; | |
175 | ||
176 | fhELambdaCluster = new TH3F | |
177 | ("hELambdaCluster","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
178 | fhELambdaCluster->SetZTitle("#lambda_{1}^{2}"); | |
179 | fhELambdaCluster->SetYTitle("#lambda_{0}^{2}"); | |
180 | fhELambdaCluster->SetXTitle("E_{Cluster} (GeV)"); | |
181 | outputContainer->Add(fhELambdaCluster) ; | |
182 | ||
183 | fhDispersionCluster = new TH2F | |
184 | ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); | |
185 | fhDispersionCluster->SetYTitle("Dispersion"); | |
186 | fhDispersionCluster->SetXTitle("p_{T Cluster} (GeV/c)"); | |
187 | outputContainer->Add(fhDispersionCluster) ; | |
188 | ||
189 | fhELambdaCellCluster = new TH3F | |
190 | ("hELambdaCellCluster","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
191 | fhELambdaCellCluster->SetZTitle("N Cells"); | |
192 | fhELambdaCellCluster->SetYTitle("#lambda_{0}^{2}"); | |
193 | fhELambdaCellCluster->SetXTitle("E_{Cluster} (GeV)"); | |
194 | outputContainer->Add(fhELambdaCellCluster) ; | |
195 | ||
196 | //Identified photon histograms | |
197 | fhNCellPhoton = new TH2F ("hNCellPhoton","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); | |
198 | fhNCellPhoton->SetXTitle("p_{T}"); | |
199 | fhNCellPhoton->SetYTitle("N_{Cells}"); | |
200 | outputContainer->Add(fhNCellPhoton) ; | |
201 | ||
202 | fhLambdaPhoton = new TH3F | |
203 | ("hLambdaPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
204 | fhLambdaPhoton->SetZTitle("#lambda_{1}^{2}"); | |
205 | fhLambdaPhoton->SetYTitle("#lambda_{0}^{2}"); | |
206 | fhLambdaPhoton->SetXTitle("p_{T Photon} (GeV/c)"); | |
207 | outputContainer->Add(fhLambdaPhoton) ; | |
208 | ||
209 | fhELambdaPhoton = new TH3F | |
210 | ("hELambdaPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
211 | fhELambdaPhoton->SetZTitle("#lambda_{1}^{2}"); | |
212 | fhELambdaPhoton->SetYTitle("#lambda_{0}^{2}"); | |
213 | fhELambdaPhoton->SetXTitle("E_{Photon} (GeV)"); | |
214 | outputContainer->Add(fhELambdaPhoton) ; | |
215 | ||
216 | fhDispersionPhoton = new TH2F | |
217 | ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); | |
218 | fhDispersionPhoton->SetYTitle("Dispersion"); | |
219 | fhDispersionPhoton->SetXTitle("p_{T Photon} (GeV/c)"); | |
220 | outputContainer->Add(fhDispersionPhoton) ; | |
221 | ||
222 | fhELambdaCellPhoton = new TH3F | |
223 | ("hELambdaCellPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
224 | fhELambdaCellPhoton->SetZTitle("N Cells"); | |
225 | fhELambdaCellPhoton->SetYTitle("#lambda_{0}^{2}"); | |
226 | fhELambdaCellPhoton->SetXTitle("E_{Photon} (GeV)"); | |
227 | outputContainer->Add(fhELambdaCellPhoton) ; | |
228 | ||
229 | //Identified Pi0 histograms | |
230 | fhNCellPi0 = new TH2F ("hNCellPi0","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); | |
231 | fhNCellPi0->SetXTitle("p_{T}"); | |
232 | fhNCellPi0->SetYTitle("N_{Cells}"); | |
233 | outputContainer->Add(fhNCellPi0) ; | |
234 | ||
235 | fhLambdaPi0 = new TH3F | |
236 | ("hLambdaPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
237 | fhLambdaPi0->SetZTitle("#lambda_{1}^{2}"); | |
238 | fhLambdaPi0->SetYTitle("#lambda_{0}^{2}"); | |
239 | fhLambdaPi0->SetXTitle("p_{T Pi0} (GeV/c)"); | |
240 | outputContainer->Add(fhLambdaPi0) ; | |
241 | ||
242 | fhELambdaPi0 = new TH3F | |
243 | ("hELambdaPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
244 | fhELambdaPi0->SetZTitle("#lambda_{1}^{2}"); | |
245 | fhELambdaPi0->SetYTitle("#lambda_{0}^{2}"); | |
246 | fhELambdaPi0->SetXTitle("E_{Pi0} (GeV)"); | |
247 | outputContainer->Add(fhELambdaPi0) ; | |
248 | ||
249 | fhDispersionPi0 = new TH2F | |
250 | ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); | |
251 | fhDispersionPi0->SetYTitle("Dispersion"); | |
252 | fhDispersionPi0->SetXTitle("p_{T Pi0} (GeV/c)"); | |
253 | outputContainer->Add(fhDispersionPi0) ; | |
254 | ||
255 | fhELambdaCellPi0 = new TH3F | |
256 | ("hELambdaCellPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
257 | fhELambdaCellPi0->SetZTitle("N Cells"); | |
258 | fhELambdaCellPi0->SetYTitle("#lambda_{0}^{2}"); | |
259 | fhELambdaCellPi0->SetXTitle("E_{Pi0} (GeV)"); | |
260 | outputContainer->Add(fhELambdaCellPi0) ; | |
261 | ||
262 | //Identified charged hadron histograms | |
263 | fhNCellChargedHadron = new TH2F ("hNCellChargedHadron","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); | |
264 | fhNCellChargedHadron->SetXTitle("p_{T}"); | |
265 | fhNCellChargedHadron->SetYTitle("N_{Cells}"); | |
266 | outputContainer->Add(fhNCellChargedHadron) ; | |
267 | ||
268 | fhLambdaChargedHadron = new TH3F | |
269 | ("hLambdaChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
270 | fhLambdaChargedHadron->SetZTitle("#lambda_{1}^{2}"); | |
271 | fhLambdaChargedHadron->SetYTitle("#lambda_{0}^{2}"); | |
272 | fhLambdaChargedHadron->SetXTitle("p_{T ChargedHadron} (GeV/c)"); | |
273 | outputContainer->Add(fhLambdaChargedHadron) ; | |
274 | ||
275 | fhELambdaChargedHadron = new TH3F | |
276 | ("hELambdaChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); | |
277 | fhELambdaChargedHadron->SetZTitle("#lambda_{1}^{2}"); | |
278 | fhELambdaChargedHadron->SetYTitle("#lambda_{0}^{2}"); | |
279 | fhELambdaChargedHadron->SetXTitle("E_{ChargedHadron} (GeV)"); | |
280 | outputContainer->Add(fhELambdaChargedHadron) ; | |
281 | ||
282 | fhDispersionChargedHadron = new TH2F | |
283 | ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); | |
284 | fhDispersionChargedHadron->SetYTitle("Dispersion"); | |
285 | fhDispersionChargedHadron->SetXTitle("p_{T ChargedHadron} (GeV/c)"); | |
286 | outputContainer->Add(fhDispersionChargedHadron) ; | |
287 | ||
288 | fhELambdaCellChargedHadron = new TH3F | |
289 | ("hELambdaCellChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
290 | fhELambdaCellChargedHadron->SetZTitle("N Cells"); | |
291 | fhELambdaCellChargedHadron->SetYTitle("#lambda_{0}^{2}"); | |
292 | fhELambdaCellChargedHadron->SetXTitle("E_{ChargedHadron} (GeV)"); | |
293 | outputContainer->Add(fhELambdaCellChargedHadron) ; | |
294 | ||
295 | if(IsDataMC()){ | |
296 | fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50); | |
297 | fhDeltaE->SetXTitle("#Delta E (GeV)"); | |
298 | outputContainer->Add(fhDeltaE); | |
299 | ||
300 | fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50); | |
301 | fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)"); | |
302 | outputContainer->Add(fhDeltaPt); | |
303 | ||
304 | fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2); | |
305 | fhRatioE->SetXTitle("E_{reco}/E_{gen}"); | |
306 | outputContainer->Add(fhRatioE); | |
307 | ||
308 | fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2); | |
309 | fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}"); | |
310 | outputContainer->Add(fhRatioPt); | |
311 | ||
312 | fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); | |
313 | fh2E->SetXTitle("E_{rec} (GeV)"); | |
314 | fh2E->SetYTitle("E_{gen} (GeV)"); | |
315 | outputContainer->Add(fh2E); | |
316 | ||
317 | fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); | |
318 | fh2Pt->SetXTitle("p_{T,rec} (GeV/c)"); | |
319 | fh2Pt->SetYTitle("p_{T,gen} (GeV/c)"); | |
320 | outputContainer->Add(fh2Pt); | |
321 | ||
322 | fhEMCPhoton = new TH1F("hEMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); | |
323 | //fhEMCPhoton->SetYTitle("N"); | |
324 | fhEMCPhoton->SetXTitle("E_{#gamma}(GeV)"); | |
325 | outputContainer->Add(fhEMCPhoton) ; | |
326 | ||
327 | fhPhiMCPhoton = new TH2F | |
328 | ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
329 | fhPhiMCPhoton->SetYTitle("#phi"); | |
330 | fhPhiMCPhoton->SetXTitle("E_{ #gamma} (GeV)"); | |
331 | outputContainer->Add(fhPhiMCPhoton) ; | |
332 | ||
333 | fhEtaMCPhoton = new TH2F | |
334 | ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
335 | fhEtaMCPhoton->SetYTitle("#eta"); | |
336 | fhEtaMCPhoton->SetXTitle("E_{ #gamma} (GeV)"); | |
337 | outputContainer->Add(fhEtaMCPhoton) ; | |
338 | ||
339 | fhLambdaMCPhoton = new TH3F | |
340 | ("hLambdaMCPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
341 | fhLambdaMCPhoton->SetZTitle("N_{Cells}"); | |
342 | fhLambdaMCPhoton->SetYTitle("#lambda_{0}^{2}"); | |
343 | fhLambdaMCPhoton->SetXTitle("E_{#gamma} (GeV)"); | |
344 | outputContainer->Add(fhLambdaMCPhoton) ; | |
345 | ||
346 | fhPhotTrueE = new TH3F | |
347 | ("hPhotTrueE","#lambda_{#gamma}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); | |
348 | fhPhotTrueE->SetZTitle("#lambda_{0}^{2}"); | |
349 | fhPhotTrueE->SetYTitle("Recons. E_{#gamma}"); | |
350 | fhPhotTrueE->SetXTitle("MC Truth E_{#gamma} (GeV)"); | |
351 | outputContainer->Add(fhPhotTrueE) ; | |
352 | ||
353 | fhRatioEPhoton = new TH1F ("hRatioEPhoton","E_{Reco #gamma}/E_{MC truth #gamma}", 40,0,2); | |
354 | fhRatioEPhoton->SetXTitle("E_{reco #gamma}/E_{gen #gamma}"); | |
355 | outputContainer->Add(fhRatioEPhoton); | |
356 | ||
357 | fhEMCPi0 = new TH1F("hEMCPi0","Number of #pi^0 over calorimeter",nptbins,ptmin,ptmax); | |
358 | fhEMCPi0->SetYTitle("N"); | |
359 | fhEMCPi0->SetXTitle("E_{ #pi^0}(GeV)"); | |
360 | outputContainer->Add(fhEMCPi0) ; | |
361 | ||
362 | fhPhiMCPi0 = new TH2F | |
363 | ("hPhiMCPi0","#phi_{#pi^0}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
364 | fhPhiMCPi0->SetYTitle("#phi"); | |
365 | fhPhiMCPi0->SetXTitle("E_{ #pi^0} (GeV)"); | |
366 | outputContainer->Add(fhPhiMCPi0) ; | |
367 | ||
368 | fhEtaMCPi0 = new TH2F | |
369 | ("hEtaMCPi0","#phi_{#pi^0}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
370 | fhEtaMCPi0->SetYTitle("#eta"); | |
371 | fhEtaMCPi0->SetXTitle("E_{ #pi^0} (GeV)"); | |
372 | outputContainer->Add(fhEtaMCPi0) ; | |
373 | ||
374 | fhLambdaMCPi0 = new TH3F | |
375 | ("hLambdaMCPi0","#lambda_{#pi^0}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
376 | fhLambdaMCPi0->SetZTitle("N_{Cells}"); | |
377 | fhLambdaMCPi0->SetYTitle("#lambda_{0}^{2}"); | |
378 | fhLambdaMCPi0->SetXTitle("E_{ #pi^0} (GeV)"); | |
379 | outputContainer->Add(fhLambdaMCPi0) ; | |
380 | ||
381 | fhProductionDistance = new TH2F | |
382 | ("hProductionDistance","#lambda_{#pi^0}",nptbins,ptmin,ptmax,160,0,800); | |
383 | fhProductionDistance->SetYTitle("Distance of production of Pi0 from vertex (cm)"); | |
384 | fhProductionDistance->SetXTitle("E_{ #pi^0} (GeV)"); | |
385 | outputContainer->Add(fhProductionDistance) ; | |
386 | ||
387 | fhProductionRadius = new TH2F | |
388 | ("hProductionRadius","#lambda_{#pi^0}",nptbins,ptmin,ptmax,160,0,800); | |
389 | fhProductionRadius->SetYTitle("Radius of production of Pi0 from vertex (cm)"); | |
390 | fhProductionRadius->SetXTitle("E_{ #pi^0} (GeV)"); | |
391 | outputContainer->Add(fhProductionRadius) ; | |
392 | ||
393 | fhDecayAngle = new TH3F | |
394 | ("hDecayAngle","#lambda_{#pi^0}",nptbins,ptmin,ptmax,100,-0.5,0.5,200,0,2); | |
395 | fhDecayAngle->SetZTitle("#lambda_{0}^{2}"); | |
396 | fhDecayAngle->SetYTitle("Opening angle of the #pi^{0} decay (radians)"); | |
397 | fhDecayAngle->SetXTitle("E_{ #pi^0} (GeV)"); | |
398 | outputContainer->Add(fhDecayAngle) ; | |
399 | ||
400 | fhPi0TrueE = new TH3F | |
401 | ("hPi0TrueE","#lambda_{#pi^0}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); | |
402 | fhPi0TrueE->SetZTitle("#lambda_{0}^{2}"); | |
403 | fhPi0TrueE->SetYTitle("Recons. E_{#pi^0}"); | |
404 | fhPi0TrueE->SetXTitle("MC Truth E_{#pi^0} (GeV)"); | |
405 | outputContainer->Add(fhPi0TrueE) ; | |
406 | ||
407 | fhRatioEPi0 = new TH1F ("hRatioEPi0","E_{Reco #pi^{0}}/E_{MC truth #pi^{0}}", 40,0,2); | |
408 | fhRatioEPi0->SetXTitle("E_{reco #pi^{0}}/E_{gen #pi^{0}}"); | |
409 | outputContainer->Add(fhRatioEPi0); | |
410 | ||
411 | fhEMCPion = new TH1F("hEMCPion","Number of #pi over calorimeter",nptbins,ptmin,ptmax); | |
412 | fhEMCPion->SetYTitle("N"); | |
413 | fhEMCPion->SetXTitle("E_{ #pi}(GeV)"); | |
414 | outputContainer->Add(fhEMCPion) ; | |
415 | ||
416 | fhPhiMCPion = new TH2F | |
417 | ("hPhiMCPion","#phi_{#pi}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
418 | fhPhiMCPion->SetYTitle("#phi"); | |
419 | fhPhiMCPion->SetXTitle("E_{ #pi} (GeV)"); | |
420 | outputContainer->Add(fhPhiMCPion) ; | |
421 | ||
422 | fhEtaMCPion = new TH2F | |
423 | ("hEtaMCPion","#phi_{#pi}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
424 | fhEtaMCPion->SetYTitle("#eta"); | |
425 | fhEtaMCPion->SetXTitle("E_{ #pi} (GeV)"); | |
426 | outputContainer->Add(fhEtaMCPion) ; | |
427 | ||
428 | fhLambdaMCPion = new TH3F | |
429 | ("hLambdaMCPion","#lambda_{#pi}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
430 | fhLambdaMCPion->SetZTitle("N_{Cells}"); | |
431 | fhLambdaMCPion->SetYTitle("#lambda_{0}^{2}"); | |
432 | fhLambdaMCPion->SetXTitle("E_{ #pi} (GeV)"); | |
433 | outputContainer->Add(fhLambdaMCPion) ; | |
434 | ||
435 | fhPionTrueE = new TH3F | |
436 | ("hPionTrueE","#lambda_{#pi}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); | |
437 | fhPionTrueE->SetZTitle("#lambda_{0}^{2}"); | |
438 | fhPionTrueE->SetYTitle("Recons. E_{#pi}"); | |
439 | fhPionTrueE->SetXTitle("MC Truth E_{#pi} (GeV)"); | |
440 | outputContainer->Add(fhPionTrueE) ; | |
441 | ||
442 | fhRatioEPion = new TH1F ("hRatioEPion","E_{Reco #pi^{#pm}}/E_{MC truth #pi^{#pm}}", 40,0,2); | |
443 | fhRatioEPion->SetXTitle("E_{reco #pi^{#pm}}/E_{gen #pi^{#pm}}"); | |
444 | outputContainer->Add(fhRatioEPion); | |
445 | ||
446 | fhEMCProton = new TH1F("hEMCProton","Number of p over calorimeter",nptbins,ptmin,ptmax); | |
447 | fhEMCProton->SetYTitle("N"); | |
448 | fhEMCProton->SetXTitle("E_{ p}(GeV)"); | |
449 | outputContainer->Add(fhEMCProton) ; | |
450 | ||
451 | fhPhiMCProton = new TH2F | |
452 | ("hPhiMCProton","#phi_{p}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
453 | fhPhiMCProton->SetYTitle("#phi"); | |
454 | fhPhiMCProton->SetXTitle("E_{ p} (GeV)"); | |
455 | outputContainer->Add(fhPhiMCProton) ; | |
456 | ||
457 | fhEtaMCProton = new TH2F | |
458 | ("hEtaMCProton","#phi_{p}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
459 | fhEtaMCProton->SetYTitle("#eta"); | |
460 | fhEtaMCProton->SetXTitle("E_{ p} (GeV)"); | |
461 | outputContainer->Add(fhEtaMCProton) ; | |
462 | ||
463 | fhLambdaMCProton = new TH3F | |
464 | ("hLambdaMCProton","#lambda_{p}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
465 | fhLambdaMCProton->SetZTitle("N_{Cells}"); | |
466 | fhLambdaMCProton->SetYTitle("#lambda_{0}^{2}"); | |
467 | fhLambdaMCProton->SetXTitle("E_{ p} (GeV)"); | |
468 | outputContainer->Add(fhLambdaMCProton) ; | |
469 | ||
470 | fhProtonTrueE = new TH3F | |
471 | ("hProtonTrueE","#lambda_{p}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); | |
472 | fhProtonTrueE->SetZTitle("#lambda_{0}^{2}"); | |
473 | fhProtonTrueE->SetYTitle("Recons. E_{p}"); | |
474 | fhProtonTrueE->SetXTitle("MC Truth E_{p} (GeV)"); | |
475 | outputContainer->Add(fhProtonTrueE) ; | |
476 | ||
477 | fhRatioEProton = new TH1F ("hRatioEProton","E_{Reco p}/E_{MC truth p}", 40,0,2); | |
478 | fhRatioEProton->SetXTitle("E_{reco p}/E_{gen p}"); | |
479 | outputContainer->Add(fhRatioEProton); | |
480 | ||
481 | fhEMCAntiProton = new TH1F("hEMCAntiProton","Number of #bar{p} over calorimeter",nptbins,ptmin,ptmax); | |
482 | fhEMCAntiProton->SetYTitle("N"); | |
483 | fhEMCAntiProton->SetXTitle("E_{#bar{p}}(GeV)"); | |
484 | outputContainer->Add(fhEMCAntiProton) ; | |
485 | ||
486 | fhPhiMCAntiProton = new TH2F | |
487 | ("hPhiMCAntiProton","#phi_{#bar{p}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
488 | fhPhiMCAntiProton->SetYTitle("#phi"); | |
489 | fhPhiMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)"); | |
490 | outputContainer->Add(fhPhiMCAntiProton) ; | |
491 | ||
492 | fhEtaMCAntiProton = new TH2F | |
493 | ("hEtaMCAntiProton","#phi_{#bar{p}}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
494 | fhEtaMCAntiProton->SetYTitle("#eta"); | |
495 | fhEtaMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)"); | |
496 | outputContainer->Add(fhEtaMCAntiProton) ; | |
497 | ||
498 | fhLambdaMCAntiProton = new TH3F | |
499 | ("hLambdaMCAntiProton","#lambda_{#bar{p}}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
500 | fhLambdaMCAntiProton->SetZTitle("N_{Cells}"); | |
501 | fhLambdaMCAntiProton->SetYTitle("#lambda_{0}^{2}"); | |
502 | fhLambdaMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)"); | |
503 | outputContainer->Add(fhLambdaMCAntiProton) ; | |
504 | ||
505 | fhAntiProtonTrueE = new TH3F | |
506 | ("hAntiProtonTrueE","#lambda_{#bar{p}}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); | |
507 | fhAntiProtonTrueE->SetZTitle("#lambda_{0}^{2}"); | |
508 | fhAntiProtonTrueE->SetYTitle("Recons. E_{#bar{p}}"); | |
509 | fhAntiProtonTrueE->SetXTitle("MC Truth E_{#bar{p}} (GeV)"); | |
510 | outputContainer->Add(fhAntiProtonTrueE) ; | |
511 | ||
512 | fhRatioEAntiProton = new TH1F ("hRatioEAntiProton","E_{Reco #bar{p}}/E_{MC truth #bar{p}}", 40,0,2); | |
513 | fhRatioEAntiProton->SetXTitle("E_{reco #bar{p}}/E_{gen #bar{p}}"); | |
514 | outputContainer->Add(fhRatioEAntiProton); | |
515 | ||
516 | fhEMCNeutron = new TH1F("hEMCNeutron","Number of p over calorimeter",nptbins,ptmin,ptmax); | |
517 | fhEMCNeutron->SetYTitle("N"); | |
518 | fhEMCNeutron->SetXTitle("E_{ p}(GeV)"); | |
519 | outputContainer->Add(fhEMCNeutron) ; | |
520 | ||
521 | fhPhiMCNeutron = new TH2F | |
522 | ("hPhiMCNeutron","#phi_{p}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
523 | fhPhiMCNeutron->SetYTitle("#phi"); | |
524 | fhPhiMCNeutron->SetXTitle("E_{ p} (GeV)"); | |
525 | outputContainer->Add(fhPhiMCNeutron) ; | |
526 | ||
527 | fhEtaMCNeutron = new TH2F | |
528 | ("hEtaMCNeutron","#phi_{p}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
529 | fhEtaMCNeutron->SetYTitle("#eta"); | |
530 | fhEtaMCNeutron->SetXTitle("E_{ p} (GeV)"); | |
531 | outputContainer->Add(fhEtaMCNeutron) ; | |
532 | ||
533 | fhLambdaMCNeutron = new TH3F | |
534 | ("hLambdaMCNeutron","#lambda_{p}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
535 | fhLambdaMCNeutron->SetZTitle("N_{Cells}"); | |
536 | fhLambdaMCNeutron->SetYTitle("#lambda_{0}^{2}"); | |
537 | fhLambdaMCNeutron->SetXTitle("E_{ p} (GeV)"); | |
538 | outputContainer->Add(fhLambdaMCNeutron) ; | |
539 | ||
540 | fhNeutronTrueE = new TH3F | |
541 | ("hNeutronTrueE","#lambda_{n}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); | |
542 | fhNeutronTrueE->SetZTitle("#lambda_{0}^{2}"); | |
543 | fhNeutronTrueE->SetYTitle("Recons. E_{n}"); | |
544 | fhNeutronTrueE->SetXTitle("MC Truth E_{n} (GeV)"); | |
545 | outputContainer->Add(fhNeutronTrueE) ; | |
546 | ||
547 | fhRatioENeutron = new TH1F ("hRatioENeutron","E_{Reco n}/E_{MC truth n}", 40,0,2); | |
548 | fhRatioENeutron->SetXTitle("E_{reco n}/E_{gen n}"); | |
549 | outputContainer->Add(fhRatioENeutron); | |
550 | ||
551 | fhEMCEta = new TH1F("hEMCEta","Number of #eta over calorimeter",nptbins,ptmin,ptmax); | |
552 | fhEMCEta->SetYTitle("N"); | |
553 | fhEMCEta->SetXTitle("E_{ #eta}(GeV)"); | |
554 | outputContainer->Add(fhEMCEta) ; | |
555 | ||
556 | fhPhiMCEta = new TH2F | |
557 | ("hPhiMCEta","#phi_{#eta}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); | |
558 | fhPhiMCEta->SetYTitle("#phi"); | |
559 | fhPhiMCEta->SetXTitle("E_{ #eta} (GeV)"); | |
560 | outputContainer->Add(fhPhiMCEta) ; | |
561 | ||
562 | fhEtaMCEta = new TH2F | |
563 | ("hEtaMCEta","#phi_{#eta}",nptbins,ptmin,ptmax,netabins,etamin,etamax); | |
564 | fhEtaMCEta->SetYTitle("#eta"); | |
565 | fhEtaMCEta->SetXTitle("E_{ #eta} (GeV)"); | |
566 | outputContainer->Add(fhEtaMCEta) ; | |
567 | ||
568 | fhLambdaMCEta = new TH3F | |
569 | ("hLambdaMCEta","#lambda_{#eta}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); | |
570 | fhLambdaMCEta->SetZTitle("N_{Cells}"); | |
571 | fhLambdaMCEta->SetYTitle("#lambda_{0}^{2}"); | |
572 | fhLambdaMCEta->SetXTitle("E_{ #eta} (GeV)"); | |
573 | outputContainer->Add(fhLambdaMCEta) ; | |
574 | ||
575 | fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ; | |
576 | outputContainer->Add(fhPrimPt) ; | |
577 | ||
578 | }//Histos with MC | |
579 | ||
580 | return outputContainer ; | |
581 | ||
582 | } | |
583 | ||
584 | //____________________________________________________________________________ | |
585 | void AliAnaShowerParameter::Init() | |
586 | { | |
587 | ||
588 | //Init | |
589 | //Do some checks | |
590 | if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){ | |
591 | printf("AliAnaShowerParameter::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n"); | |
592 | abort(); | |
593 | } | |
594 | else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){ | |
595 | printf("AliAnaShowerParameter::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n"); | |
596 | abort(); | |
597 | } | |
598 | ||
599 | } | |
600 | ||
601 | //____________________________________________________________________________ | |
602 | void AliAnaShowerParameter::InitParameters() | |
603 | { | |
604 | ||
605 | //Initialize the parameters of the analysis. | |
606 | AddToHistogramsName("AnaPhoton_"); | |
607 | ||
608 | fCalorimeter = "PHOS" ; | |
609 | fMinDist = 2.; | |
610 | fMinDist2 = 4.; | |
611 | fMinDist3 = 5.; | |
612 | fMassCut = 0.03; //30 MeV | |
613 | ||
614 | fTimeCutMin = -1; | |
615 | fTimeCutMax = 9999999; | |
616 | fNCellsCut = 0; | |
617 | ||
618 | fRejectTrackMatch = kTRUE ; | |
619 | fCheckConversion = kFALSE; | |
620 | fAddConvertedPairsToAOD = kFALSE; | |
621 | ||
622 | } | |
623 | ||
624 | //__________________________________________________________________ | |
625 | void AliAnaShowerParameter::MakeAnalysisFillAOD() | |
626 | { | |
627 | //Do analysis and fill aods | |
628 | //Search for photons in fCalorimeter | |
629 | ||
630 | //Get vertex for photon momentum calculation | |
631 | ||
632 | for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) { | |
633 | if (!GetMixedEvent()) | |
634 | GetReader()->GetVertex(GetVertex(iev)); | |
635 | else | |
636 | GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); | |
637 | } | |
638 | ||
639 | //Select the Calorimeter of the photon | |
640 | TObjArray * pl = 0x0; | |
641 | if(fCalorimeter == "PHOS") | |
642 | pl = GetAODPHOS(); | |
643 | else if (fCalorimeter == "EMCAL") | |
644 | pl = GetAODEMCAL(); | |
645 | ||
164a1d84 | 646 | if(!pl){ |
647 | printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Careful cluster array NULL!!\n"); | |
648 | return; | |
649 | } | |
650 | ||
6d08aa8d | 651 | //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods |
652 | TLorentzVector mom, mom2 ; | |
653 | Int_t nCaloClusters = pl->GetEntriesFast(); | |
654 | Bool_t * indexConverted = new Bool_t[nCaloClusters]; | |
655 | for (Int_t i = 0; i < nCaloClusters; i++) | |
656 | indexConverted[i] = kFALSE; | |
657 | ||
658 | for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){ | |
659 | ||
660 | AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo)); | |
661 | Int_t evtIndex = 0 ; | |
662 | if (GetMixedEvent()) { | |
663 | evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; | |
664 | } | |
665 | //Cluster selection, not charged, with photon id and in fiducial cut | |
666 | ||
667 | //Input from second AOD? | |
668 | Int_t input = 0; | |
669 | ||
670 | //Get Momentum vector, | |
671 | if (input == 0) | |
672 | calo->GetMomentum(mom,GetVertex(evtIndex)) ;//Assume that come from vertex in straight line | |
673 | ||
674 | //Skip the cluster if it doesn't fit inside the cuts. | |
675 | if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; | |
676 | Double_t tof = calo->GetTOF()*1e9; | |
677 | if(tof < fTimeCutMin || tof > fTimeCutMax) continue; | |
678 | if(calo->GetNCells() <= fNCellsCut) continue; | |
679 | ||
680 | //Check acceptance selection | |
681 | if(IsFiducialCutOn()){ | |
682 | Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ; | |
683 | if(! in ) continue ; | |
684 | } | |
685 | ||
686 | //Create AOD for analysis | |
687 | AliAODPWG4Particle aodph = AliAODPWG4Particle(mom); | |
688 | Int_t label = calo->GetLabel(); | |
689 | aodph.SetLabel(label); | |
690 | aodph.SetInputFileIndex(input); | |
691 | ||
692 | //Set the indices of the original caloclusters | |
693 | aodph.SetCaloLabel(calo->GetID(),-1); | |
694 | aodph.SetDetector(fCalorimeter); | |
695 | if(GetDebug() > 1) | |
696 | printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta()); | |
697 | ||
698 | //Check Distance to Bad channel, set bit. | |
699 | Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel | |
700 | if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ; | |
701 | if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm) | |
702 | continue ; | |
703 | ||
704 | if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad); | |
705 | ||
706 | if (distBad > fMinDist3) aodph.SetDistToBad(2) ; | |
707 | else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; | |
708 | else aodph.SetDistToBad(0) ; | |
709 | ||
710 | //Skip matched clusters with tracks | |
711 | if(fRejectTrackMatch && IsTrackMatched(calo)) continue ; | |
712 | if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - TrackMatching cut passed \n"); | |
713 | ||
714 | //Set PID bits for later selection (AliAnaPi0 for example) | |
715 | //GetPDG already called in SetPIDBits. | |
716 | GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph); | |
717 | if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - PID Bits set \n"); | |
718 | ||
719 | //Play with the MC stack if available | |
720 | //Check origin of the candidates | |
721 | if(IsDataMC()){ | |
722 | aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex())); | |
723 | if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag()); | |
724 | } | |
725 | ||
726 | // Check if cluster comes from a conversion in the material in front of the calorimeter | |
727 | // Do invariant mass of all pairs, if mass is close to 0, then it is conversion. | |
728 | ||
729 | if(fCheckConversion && nCaloClusters > 1){ | |
730 | Bool_t bConverted = kFALSE; | |
731 | Int_t id2 = -1; | |
732 | ||
733 | //Check if set previously as converted couple, if so skip its use. | |
734 | if (indexConverted[icalo]) continue; | |
735 | ||
736 | for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) { | |
737 | //Check if set previously as converted couple, if so skip its use. | |
738 | if (indexConverted[jcalo]) continue; | |
739 | //printf("Check Conversion indeces %d and %d\n",icalo,jcalo); | |
740 | AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics | |
741 | Int_t evtIndex2 = 0 ; | |
742 | if (GetMixedEvent()) { | |
743 | evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ; | |
744 | } | |
745 | calo2->GetMomentum(mom2,GetVertex(evtIndex2)); | |
746 | //Check only certain regions | |
747 | Bool_t in2 = kTRUE; | |
748 | if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ; | |
749 | if(!in2) continue; | |
750 | ||
751 | //Get mass of pair, if small, take this pair. | |
752 | //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut); | |
753 | if((mom+mom2).M() < fMassCut){ | |
754 | bConverted = kTRUE; | |
755 | id2 = calo2->GetID(); | |
756 | indexConverted[jcalo]=kTRUE; | |
757 | break; | |
758 | } | |
759 | ||
760 | }//Mass loop | |
761 | ||
762 | if(bConverted){ | |
763 | if(fAddConvertedPairsToAOD){ | |
764 | //Create AOD of pair analysis | |
765 | TLorentzVector mpair = mom+mom2; | |
766 | AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair); | |
767 | aodpair.SetLabel(aodph.GetLabel()); | |
768 | aodpair.SetInputFileIndex(input); | |
769 | ||
770 | //printf("Index %d, Id %d\n",icalo, calo->GetID()); | |
771 | //Set the indeces of the original caloclusters | |
772 | aodpair.SetCaloLabel(calo->GetID(),id2); | |
773 | aodpair.SetDetector(fCalorimeter); | |
774 | aodpair.SetPdg(aodph.GetPdg()); | |
775 | aodpair.SetTag(aodph.GetTag()); | |
776 | ||
777 | //Add AOD with pair object to aod branch | |
778 | AddAODParticle(aodpair); | |
779 | //printf("\t \t both added pair\n"); | |
780 | } | |
781 | ||
782 | //Do not add the current calocluster | |
783 | continue; | |
784 | }//converted pair | |
785 | }//check conversion | |
786 | ||
787 | //Add AOD with photon object to aod branch | |
788 | AddAODParticle(aodph); | |
789 | ||
790 | }//loop; | |
791 | delete [] indexConverted; | |
792 | ||
793 | if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast()); | |
794 | ||
795 | } | |
796 | ||
797 | //__________________________________________________________________ | |
798 | void AliAnaShowerParameter::MakeAnalysisFillHistograms() | |
799 | { | |
800 | ||
801 | //Do analysis and fill histograms | |
802 | ||
803 | // Access MC information in stack if requested, check that it exists. | |
804 | AliStack * stack = 0x0; | |
805 | TParticle * primary = 0x0; | |
806 | TClonesArray * mcparticles0 = 0x0; | |
164a1d84 | 807 | //TClonesArray * mcparticles1 = 0x0; |
6d08aa8d | 808 | AliAODMCParticle * aodprimary = 0x0; |
809 | TObjArray * pl = 0x0; | |
810 | Int_t iNumCell=0; | |
811 | ||
812 | //Check if the stack is available when analysing MC data. | |
813 | if(IsDataMC()){ | |
814 | ||
815 | if(GetReader()->ReadStack()){ | |
816 | stack = GetMCStack() ; | |
817 | if(!stack) { | |
818 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n"); | |
819 | abort(); | |
820 | } | |
821 | ||
822 | } | |
823 | else if(GetReader()->ReadAODMCParticles()){ | |
824 | ||
825 | //Get the list of MC particles | |
826 | mcparticles0 = GetReader()->GetAODMCParticles(0); | |
827 | if(!mcparticles0 && GetDebug() > 0) { | |
828 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n"); | |
829 | } | |
830 | }// is data and MC | |
831 | } | |
832 | //Loop on stored AOD photons | |
833 | Int_t naod = GetOutputAODBranch()->GetEntriesFast(); | |
834 | fhNClusters->Fill(naod); | |
835 | if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod); | |
836 | ||
837 | for(Int_t iaod = 0; iaod < naod ; iaod++){ | |
838 | AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod)); | |
839 | ||
840 | if(ph->GetDetector() != fCalorimeter) continue; | |
841 | ||
842 | if(GetDebug() > 2) | |
843 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ; | |
844 | ||
845 | //Fill Cluster histograms | |
846 | Float_t ptcluster = ph->Pt(); | |
847 | Float_t phicluster = ph->Phi(); | |
848 | Float_t etacluster = ph->Eta(); | |
849 | Float_t ecluster = ph->E(); | |
850 | Float_t lambdaMainCluster = -1; //Can't set their values here, but need a default -1 value in case something is wrong | |
851 | Float_t lambdaSecondCluster = -1; | |
852 | Float_t dispcluster = -1; | |
853 | ||
854 | //Get the list of clusters from the AOD and loop over them to find the onces corresponding to the reconstructed 'photons'. | |
855 | if(fCalorimeter == "PHOS") | |
856 | pl = GetAODPHOS(); | |
857 | else if (fCalorimeter == "EMCAL") | |
858 | pl = GetAODEMCAL(); | |
164a1d84 | 859 | if(pl){ |
860 | //Some values are stored in AliAODCaloCluster objects only; we need to fetch them. | |
861 | for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){ | |
862 | AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo)); | |
863 | if (calo->GetLabel()==ph->GetLabel()) { //The Cluster is the right one for this particle | |
864 | lambdaMainCluster = calo->GetM02(); //lambda_0 | |
865 | lambdaSecondCluster = calo->GetM20(); //lambda_1 | |
866 | dispcluster = calo->GetDispersion(); //Dispersion | |
867 | iNumCell = calo->GetNCells(); | |
868 | if(GetDebug() > 2) | |
869 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster Lambda0: %3.2f, Lambda1: %3.2f, Dispersion: %3.2f, NCells: %d \n",lambdaMainCluster,lambdaSecondCluster,dispcluster,iNumCell) ; | |
870 | } | |
6d08aa8d | 871 | } |
872 | } | |
873 | ||
874 | //Fill the basic non-MC information about the cluster. | |
875 | fhPtCluster ->Fill(ptcluster); | |
876 | fhPhiCluster ->Fill(ptcluster,phicluster); | |
877 | fhEtaCluster ->Fill(ptcluster,etacluster); | |
878 | fhDispersionCluster->Fill(ptcluster,dispcluster); | |
879 | fhLambdaCluster ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); | |
880 | fhELambdaCluster ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); | |
881 | fhELambdaCellCluster ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
882 | fhNCellCluster->Fill(ecluster,iNumCell); | |
883 | ||
884 | //Fill the lambda distribution for identified particles | |
885 | if(ph->GetPdg() == AliCaloPID::kPhoton){ | |
886 | fhDispersionPhoton->Fill(ptcluster,dispcluster); | |
887 | fhLambdaPhoton ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); | |
888 | fhELambdaPhoton ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); | |
889 | fhELambdaCellPhoton ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
890 | fhNCellPhoton->Fill(ecluster,iNumCell); | |
891 | } | |
892 | if(ph->GetPdg() == AliCaloPID::kPi0){ | |
893 | fhDispersionPi0->Fill(ptcluster,dispcluster); | |
894 | fhLambdaPi0 ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); | |
895 | fhELambdaPi0 ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); | |
896 | fhELambdaCellPi0 ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
897 | fhNCellPi0->Fill(ecluster,iNumCell); | |
898 | } | |
899 | if(ph->GetPdg() == AliCaloPID::kChargedHadron){ | |
900 | fhDispersionChargedHadron->Fill(ptcluster,dispcluster); | |
901 | fhLambdaChargedHadron ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); | |
902 | fhELambdaChargedHadron ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); | |
903 | fhELambdaCellChargedHadron ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
904 | fhNCellChargedHadron->Fill(ecluster,iNumCell); | |
905 | } | |
906 | ||
907 | //Play with the MC data if available | |
908 | if(IsDataMC()){ | |
909 | //Get the tag from AliMCAnalysisUtils for PID | |
910 | Int_t tag =ph->GetTag(); | |
1263717b | 911 | if ( ph->GetLabel() < 0){ |
912 | if(GetDebug() > 0) | |
913 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? "); | |
914 | continue; | |
915 | } | |
6d08aa8d | 916 | |
917 | //Check if the tag matches on of the different particle types and fill the corresponding histograms | |
918 | if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)&&!(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))) //kMCPhoton; making sure that this is not a Pi0 cluster (it should not happen). | |
919 | { | |
920 | fhEMCPhoton ->Fill(ecluster); | |
921 | fhPhiMCPhoton ->Fill(ecluster,phicluster); | |
922 | fhEtaMCPhoton ->Fill(ecluster,etacluster); | |
923 | fhLambdaMCPhoton ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
6d08aa8d | 924 | Double_t ePhot; |
925 | //Get the true energy of the photon | |
1263717b | 926 | Int_t iCurrent = ph->GetLabel(); |
6d08aa8d | 927 | TParticle* pCurrent = stack->Particle(iCurrent); |
928 | ePhot = pCurrent->Energy(); | |
929 | fhPhotTrueE->Fill(ecluster,ePhot,lambdaMainCluster); | |
930 | fhRatioEPhoton->Fill(ecluster/ePhot); | |
931 | if(GetDebug() > 3) | |
932 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPhoton: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",ePhot,ecluster,lambdaMainCluster); | |
933 | }//kMCPhoton | |
934 | ||
935 | if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) //kMCPi0 : single cluster created by 2 photons from the same Pi0 decay | |
1263717b | 936 | { |
937 | ||
6d08aa8d | 938 | //Fill the basic information about the Pi0 cluster |
939 | fhEMCPi0 ->Fill(ecluster); | |
940 | fhPhiMCPi0 ->Fill(ecluster,phicluster); | |
941 | fhEtaMCPi0 ->Fill(ecluster,etacluster); | |
942 | fhLambdaMCPi0 ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
943 | ||
944 | //Check the position of the Pi0: does it come from the vertex or was it created by hadron-material collision? | |
945 | Int_t iCurrent = ph->GetLabel(); | |
946 | Double_t ePi0 = 0; | |
947 | TParticle* pCurrent = stack->Particle(iCurrent); | |
948 | while (pCurrent->GetPdgCode()!=111) | |
949 | { | |
950 | iCurrent = pCurrent->GetFirstMother(); | |
951 | pCurrent = stack->Particle(iCurrent); | |
952 | } | |
953 | Float_t fDistance = pCurrent->Rho(); | |
954 | Float_t fRadius = pCurrent->R(); | |
955 | TLorentzVector mom1 ; | |
956 | (stack->Particle((pCurrent->GetFirstDaughter())))->Momentum(mom1); | |
957 | TLorentzVector mom2 ; | |
958 | (stack->Particle((pCurrent->GetLastDaughter())))->Momentum(mom2); | |
959 | ePi0 = pCurrent->Energy(); | |
960 | Float_t fDecayAngle = mom1.Angle(mom2.Vect()); | |
961 | fhProductionDistance->Fill(ecluster,fDistance); | |
962 | fhProductionRadius->Fill(ecluster,fRadius); | |
963 | //Decay angle of the Pi0 | |
964 | fhDecayAngle->Fill(ecluster,fDecayAngle,lambdaMainCluster); | |
965 | //Compare the cluster energy and true energy | |
966 | fhPi0TrueE->Fill(ecluster,ePi0,lambdaMainCluster); | |
967 | fhRatioEPi0->Fill(ecluster/ePi0); | |
968 | if(GetDebug() > 3) | |
969 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPi0: True E: %3.2f, Reco E: %3.2f, Decay Opening Angle: %3.2f, Lambda0: %3.2f\n",ePi0,ecluster,fDecayAngle,lambdaMainCluster); | |
970 | } | |
971 | if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion))//kMCPion | |
972 | { | |
973 | fhEMCPion ->Fill(ecluster); | |
974 | fhPhiMCPion ->Fill(ecluster,phicluster); | |
975 | fhEtaMCPion ->Fill(ecluster,etacluster); | |
976 | fhLambdaMCPion ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
977 | ||
978 | Int_t iCurrent = ph->GetLabel(); | |
979 | Double_t ePion = 0; | |
980 | TParticle* pCurrent = stack->Particle(iCurrent); | |
981 | while ((TMath::Abs(pCurrent->GetPdgCode())!=211)&&(iCurrent>=0)) | |
982 | { | |
983 | if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops | |
984 | if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); | |
985 | } | |
986 | if ((TMath::Abs(pCurrent->GetPdgCode())==211)&&(iCurrent>=0)) | |
987 | { | |
988 | ePion = pCurrent->Energy(); | |
989 | fhPionTrueE->Fill(ecluster,ePion,lambdaMainCluster); | |
990 | fhRatioEPion->Fill(ecluster/ePion); | |
991 | } | |
992 | if(GetDebug() > 3) | |
993 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPion: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",ePion,ecluster,lambdaMainCluster); | |
994 | } | |
995 | if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton))//kMCProton | |
996 | { | |
997 | fhEMCProton ->Fill(ecluster); | |
998 | fhPhiMCProton ->Fill(ecluster,phicluster); | |
999 | fhEtaMCProton ->Fill(ecluster,etacluster); | |
1000 | fhLambdaMCProton ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
1001 | ||
1002 | Int_t iCurrent = ph->GetLabel(); | |
1003 | Double_t eProton = 0; | |
1004 | TParticle* pCurrent = stack->Particle(iCurrent); | |
1005 | while ((TMath::Abs(pCurrent->GetPdgCode())!=2212)&&(iCurrent>=0)) | |
1006 | { | |
1007 | if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops | |
1008 | if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); | |
1009 | } | |
1010 | if ((TMath::Abs(pCurrent->GetPdgCode())==2212)&&(iCurrent>=0)) | |
1011 | { | |
1012 | eProton = pCurrent->Energy(); | |
1013 | fhProtonTrueE->Fill(ecluster,eProton,lambdaMainCluster); | |
1014 | fhRatioEProton->Fill(ecluster/eProton); | |
1015 | } | |
1016 | if(GetDebug() > 3) | |
1017 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCProton: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",eProton,ecluster,lambdaMainCluster); | |
1018 | } | |
1019 | if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))//kMCAntiProton | |
1020 | { | |
1021 | fhEMCAntiProton ->Fill(ecluster); | |
1022 | fhPhiMCAntiProton ->Fill(ecluster,phicluster); | |
1023 | fhEtaMCAntiProton ->Fill(ecluster,etacluster); | |
1024 | fhLambdaMCAntiProton ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
1025 | ||
1026 | Int_t iCurrent = ph->GetLabel(); | |
1027 | Double_t eAntiProton = 0; | |
1028 | TParticle* pCurrent = stack->Particle(iCurrent); | |
1029 | while ((TMath::Abs(pCurrent->GetPdgCode())!=-2212)&&(iCurrent>=0)) | |
1030 | { | |
1031 | if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops | |
1032 | if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); | |
1033 | } | |
1034 | if ((TMath::Abs(pCurrent->GetPdgCode())==-2212)&&(iCurrent>=0)) | |
1035 | { | |
1036 | eAntiProton = pCurrent->Energy(); | |
1037 | fhAntiProtonTrueE->Fill(ecluster,eAntiProton,lambdaMainCluster); | |
1038 | fhRatioEAntiProton->Fill(ecluster/eAntiProton); | |
1039 | } | |
1040 | if(GetDebug() > 3) | |
1041 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCAntiProton: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",eAntiProton,ecluster,lambdaMainCluster); | |
1042 | } | |
1043 | if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCNeutron))//kMCNeutron | |
1044 | { | |
1045 | fhEMCNeutron ->Fill(ecluster); | |
1046 | fhPhiMCNeutron ->Fill(ecluster,phicluster); | |
1047 | fhEtaMCNeutron ->Fill(ecluster,etacluster); | |
1048 | fhLambdaMCNeutron ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
1049 | ||
1050 | Int_t iCurrent = ph->GetLabel(); | |
1051 | Double_t eNeutron = 0; | |
1052 | TParticle* pCurrent = stack->Particle(iCurrent); | |
1053 | while ((TMath::Abs(pCurrent->GetPdgCode())!=2112)&&(iCurrent>=0)) | |
1054 | { | |
1055 | if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops | |
1056 | if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); | |
1057 | } | |
1058 | if ((TMath::Abs(pCurrent->GetPdgCode())==2112)&&(iCurrent>=0)) | |
1059 | { | |
1060 | eNeutron = pCurrent->Energy(); | |
1061 | fhNeutronTrueE->Fill(ecluster,eNeutron,lambdaMainCluster); | |
1062 | fhRatioENeutron->Fill(ecluster/eNeutron); | |
1063 | } | |
1064 | if(GetDebug() > 3) | |
1065 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCNeutron: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",eNeutron,ecluster,lambdaMainCluster); | |
1066 | } | |
1067 | if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))//kMCEta | |
1068 | { | |
1069 | fhEMCEta ->Fill(ecluster); | |
1070 | fhPhiMCEta ->Fill(ecluster,phicluster); | |
1071 | fhEtaMCEta ->Fill(ecluster,etacluster); | |
1072 | fhLambdaMCEta ->Fill(ecluster,lambdaMainCluster,iNumCell); | |
1073 | } | |
1074 | ||
1075 | // Access MC information in stack if requested, check that it exists. | |
1076 | Int_t label =ph->GetLabel(); | |
1077 | if(label < 0) { | |
1078 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label); | |
1079 | continue; | |
1080 | } | |
1081 | ||
1082 | Float_t eprim = 0; | |
1083 | Float_t ptprim = 0; | |
1084 | if(GetReader()->ReadStack()){ | |
1085 | ||
1086 | if(label >= stack->GetNtrack()) { | |
1087 | if(GetDebug() > 2) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); | |
1088 | continue ; | |
1089 | } | |
1090 | ||
1091 | primary = stack->Particle(label); | |
1092 | if(!primary){ | |
1093 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label); | |
1094 | continue; | |
1095 | } | |
1096 | eprim = primary->Energy(); | |
1097 | ptprim = primary->Pt(); | |
1098 | } | |
1099 | else if(GetReader()->ReadAODMCParticles()){ | |
1100 | //Check which is the input | |
1101 | if(ph->GetInputFileIndex() == 0){ | |
1102 | if(!mcparticles0) continue; | |
1103 | if(label >= mcparticles0->GetEntriesFast()) { | |
1104 | if(GetDebug() > 2) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",label, mcparticles0->GetEntriesFast()); | |
1105 | continue ; | |
1106 | } | |
1107 | //Get the particle | |
1108 | aodprimary = (AliAODMCParticle*) mcparticles0->At(label); | |
1109 | ||
1110 | } | |
164a1d84 | 1111 | // else {//Second input |
1112 | // if(!mcparticles1) continue; | |
1113 | // if(label >= mcparticles1->GetEntriesFast()) { | |
1114 | // if(GetDebug() > 2) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",label, mcparticles1->GetEntriesFast()); | |
1115 | // continue ; | |
1116 | // } | |
1117 | // //Get the particle | |
1118 | // aodprimary = (AliAODMCParticle*) mcparticles1->At(label); | |
1119 | // }//second input | |
6d08aa8d | 1120 | |
1121 | if(!aodprimary){ | |
1122 | printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label); | |
1123 | continue; | |
1124 | } | |
1125 | ||
1126 | eprim = aodprimary->E(); | |
1127 | ptprim = aodprimary->Pt(); | |
1128 | } | |
1129 | //Write down MC truth information concerning all clusters | |
1130 | fh2E ->Fill(ecluster, eprim); | |
1131 | fh2Pt ->Fill(ptcluster, ptprim); | |
1132 | fhDeltaE ->Fill(eprim-ecluster); | |
1133 | fhDeltaPt->Fill(ptprim-ptcluster); | |
1134 | if(eprim > 0) fhRatioE ->Fill(ecluster/eprim); | |
1135 | if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim); | |
1136 | }//Histograms with MC | |
1137 | }// aod loop | |
1138 | } | |
1139 | ||
1140 | ||
1141 | //__________________________________________________________________ | |
1142 | void AliAnaShowerParameter::Print(const Option_t * opt) const | |
1143 | { | |
1144 | //Print some relevant parameters set for the analysis | |
1145 | ||
1146 | if(! opt) | |
1147 | return; | |
1148 | ||
1149 | printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ; | |
1150 | AliAnaPartCorrBaseClass::Print(" "); | |
1151 | ||
1152 | printf("Calorimeter = %s\n", fCalorimeter.Data()) ; | |
1153 | printf("Min Distance to Bad Channel = %2.1f\n",fMinDist); | |
1154 | printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2); | |
1155 | printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); | |
1156 | printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch); | |
1157 | printf("Check Pair Conversion = %d\n",fCheckConversion); | |
1158 | printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD); | |
1159 | printf("Conversion pair mass cut = %f\n",fMassCut); | |
1160 | printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax); | |
1161 | printf("Number of cells in cluster is > %f \n", fNCellsCut); | |
1162 | printf(" \n") ; | |
1163 | ||
1164 | } |