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