10 #include "TParticle.h"
15 #include "THashList.h"
17 #include "AliAnalysisManager.h"
18 #include "AliMCEventHandler.h"
19 #include "AliMCEvent.h"
21 #include "AliAnalysisTaskSE.h"
22 #include "AliPHOSHijingEfficiency.h"
23 #include "AliCaloPhoton.h"
24 #include "AliPHOSGeometry.h"
25 #include "AliPHOSEsdCluster.h"
26 #include "AliPHOSCalibData.h"
27 #include "AliESDEvent.h"
28 #include "AliESDCaloCells.h"
29 #include "AliESDVertex.h"
30 #include "AliESDtrackCuts.h"
33 #include "AliCDBManager.h"
34 #include "AliCentrality.h"
35 #include "AliESDtrackCuts.h"
36 #include "AliEventplane.h"
38 #include "AliOADBContainer.h"
40 // Calculates pi0 and photon efficiencies using Hijing events
41 // Also calculates pi0 contamination: non-vertex pi0s
42 // And photon contamination:
43 // - off-vertex photons
44 // - misidentified photons
45 // Authors: Dmitri Peressounko
48 ClassImp(AliPHOSHijingEfficiency)
50 //________________________________________________________________________
51 Double_t rnlin(Double_t x)
53 //This is standard PHOS non-linearity used in reconstruction
54 //and re-calibration by 2.2%
55 return 1.022*(0.0241+1.0504*x+0.000249*x*x) ;
59 //________________________________________________________________________
60 AliPHOSHijingEfficiency::AliPHOSHijingEfficiency(const char *name)
61 : AliAnalysisTaskSE(name),
63 fOutputContainer(0x0),
66 fRP(0.),fRPV0A(0.),fRPV0C(0.),fHaveTPCRP(0),
74 for(Int_t i=0;i<1;i++){
75 for(Int_t j=0;j<10;j++)
76 for(Int_t k=0;k<11;k++)
77 fPHOSEvents[i][j][k]=0 ;
80 // Output slots #0 write into a TH1 container
81 DefineOutput(1,TList::Class());
83 // Set bad channel map
84 for(Int_t i=0; i<6; i++){
85 fPHOSBadMap[i]=new TH2I(Form("PHOS_BadMap_mod%d",i),"Bad Modules map",64,0.,64.,56,0.,56.) ;
91 //________________________________________________________________________
92 void AliPHOSHijingEfficiency::UserCreateOutputObjects()
98 const Int_t nRuns=200 ;
101 if(fOutputContainer != NULL){
102 delete fOutputContainer;
104 fOutputContainer = new THashList();
105 fOutputContainer->SetOwner(kTRUE);
106 PostData(1, fOutputContainer);
108 //========QA histograms=======
111 fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 10,0.,10.,nRuns,0.,float(nRuns))) ;
112 fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 10,0.,10.)) ;
114 //vertex distribution
115 fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
118 fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
119 fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
120 fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
121 fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
122 fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,20.,100,0.,100.)) ;
123 fOutputContainer->Add(new TH2F("hCluEvsCluM","ClusterMult vs E",200,0.,20.,100,0.,20.)) ;
125 fOutputContainer->Add(new TH3F("hCPVr","CPV radius",100,0.,20.,100,0.,20.,10,0.,100.));
128 fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
129 fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
130 fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
132 fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
133 fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
134 fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
136 fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
137 fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
138 fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
140 fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
141 fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
142 fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
144 //Single photon and pi0 spectrum
145 Int_t nPtPhot = 300 ;
146 Double_t ptPhotMax = 30 ;
151 //PHOS calibration QA
152 fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
153 fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
154 fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
155 fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
156 fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
157 fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
160 for(Int_t cent=0; cent<6; cent++){
162 fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad2_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
163 fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad4_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
164 fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad6_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
166 fOutputContainer->Add(new TH1F(Form("hPhotAll_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
167 fOutputContainer->Add(new TH1F(Form("hPhotAllcore_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
168 fOutputContainer->Add(new TH1F(Form("hPhotAllwou_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
169 fOutputContainer->Add(new TH1F(Form("hPhotDisp_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
170 fOutputContainer->Add(new TH1F(Form("hPhotDisp2_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
171 fOutputContainer->Add(new TH1F(Form("hPhotDispcore_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
172 fOutputContainer->Add(new TH1F(Form("hPhotDisp2core_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
173 fOutputContainer->Add(new TH1F(Form("hPhotDispwou_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
174 fOutputContainer->Add(new TH1F(Form("hPhotCPV_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
175 fOutputContainer->Add(new TH1F(Form("hPhotCPVcore_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
176 fOutputContainer->Add(new TH1F(Form("hPhotCPV2_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
177 fOutputContainer->Add(new TH1F(Form("hPhotBoth_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
178 fOutputContainer->Add(new TH1F(Form("hPhotBothcore_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
179 fOutputContainer->Add(new TH1F(Form("hPhotBoth2_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
180 fOutputContainer->Add(new TH1F(Form("hPhotBoth2core_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
183 fOutputContainer->Add(new TH2F(Form("hPi0All_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
184 fOutputContainer->Add(new TH2F(Form("hPi0Allcore_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
185 fOutputContainer->Add(new TH2F(Form("hPi0Allwou_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
186 fOutputContainer->Add(new TH2F(Form("hPi0Disp_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
187 fOutputContainer->Add(new TH2F(Form("hPi0Disp2_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
188 fOutputContainer->Add(new TH2F(Form("hPi0Dispcore_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
189 fOutputContainer->Add(new TH2F(Form("hPi0Disp2core_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
190 fOutputContainer->Add(new TH2F(Form("hPi0Dispwou_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
191 fOutputContainer->Add(new TH2F(Form("hPi0CPV_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
192 fOutputContainer->Add(new TH2F(Form("hPi0CPVcore_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
193 fOutputContainer->Add(new TH2F(Form("hPi0CPV2_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
194 fOutputContainer->Add(new TH2F(Form("hPi0Both_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
195 fOutputContainer->Add(new TH2F(Form("hPi0Bothcore_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
196 fOutputContainer->Add(new TH2F(Form("hPi0Both2_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
197 fOutputContainer->Add(new TH2F(Form("hPi0Both2core_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
200 snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
201 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
202 snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
203 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
204 snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
205 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
206 snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
207 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
208 snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
209 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
210 snprintf(key,55,"hPi0Both2_a07_cen%d",cent) ;
211 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
213 snprintf(key,55,"hSingleAll_cen%d",cent) ;
214 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
215 snprintf(key,55,"hSingleAllcore_cen%d",cent) ;
216 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
217 snprintf(key,55,"hSingleAllwou_cen%d",cent) ;
218 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
219 snprintf(key,55,"hSingleDisp_cen%d",cent) ;
220 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
221 snprintf(key,55,"hSingleDisp2_cen%d",cent) ;
222 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
223 snprintf(key,55,"hSingleDispcore_cen%d",cent) ;
224 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
225 snprintf(key,55,"hSingleDisp2core_cen%d",cent) ;
226 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
227 snprintf(key,55,"hSingleDispwou_cen%d",cent) ;
228 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
229 snprintf(key,55,"hSingleCPV_cen%d",cent) ;
230 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
231 snprintf(key,55,"hSingleCPVcore_cen%d",cent) ;
232 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
233 snprintf(key,55,"hSingleCPV2_cen%d",cent) ;
234 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
235 snprintf(key,55,"hSingleBoth_cen%d",cent) ;
236 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
237 snprintf(key,55,"hSingleBothcore_cen%d",cent) ;
238 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
239 snprintf(key,55,"hSingleBoth2_cen%d",cent) ;
240 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
241 snprintf(key,55,"hSingleBoth2core_cen%d",cent) ;
242 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
245 snprintf(key,55,"hMiPi0All_cen%d",cent) ;
246 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
247 snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ;
248 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
249 snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ;
250 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
251 snprintf(key,55,"hMiPi0Disp_cen%d",cent) ;
252 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
253 snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ;
254 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
255 snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ;
256 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
257 snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ;
258 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
259 snprintf(key,55,"hMiPi0Disp2core_cen%d",cent) ;
260 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
261 snprintf(key,55,"hMiPi0CPV_cen%d",cent) ;
262 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
263 snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ;
264 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
265 snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ;
266 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
267 snprintf(key,55,"hMiPi0Both_cen%d",cent) ;
268 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
269 snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ;
270 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
271 snprintf(key,55,"hMiPi0Both2_cen%d",cent) ;
272 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
273 snprintf(key,55,"hMiPi0Both2core_cen%d",cent) ;
274 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
276 snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
277 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
278 snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
279 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
280 snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
281 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
282 snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
283 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
284 snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
285 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
286 snprintf(key,55,"hMiPi0Both2_a07_cen%d",cent) ;
287 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
289 snprintf(key,55,"hMiSingleAll_cen%d",cent) ;
290 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
291 snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ;
292 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
293 snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ;
294 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
295 snprintf(key,55,"hMiSingleDisp_cen%d",cent) ;
296 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
297 snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ;
298 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
299 snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ;
300 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
301 snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ;
302 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
303 snprintf(key,55,"hMiSingleDisp2core_cen%d",cent) ;
304 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
305 snprintf(key,55,"hMiSingleCPV_cen%d",cent) ;
306 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
307 snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ;
308 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
309 snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ;
310 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
311 snprintf(key,55,"hMiSingleBoth_cen%d",cent) ;
312 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
313 snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ;
314 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
315 snprintf(key,55,"hMiSingleBoth2_cen%d",cent) ;
316 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
317 snprintf(key,55,"hMiSingleBoth2core_cen%d",cent) ;
318 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
323 for(Int_t cent=0; cent<6; cent++){
324 snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
325 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
326 snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
327 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
328 snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
329 fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
330 snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
331 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
332 snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
333 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
334 snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
335 fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
336 snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
337 fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
338 snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
339 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
340 snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
341 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
342 snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
343 fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
344 snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
345 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
346 snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
347 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
349 fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
350 fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
351 fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
357 fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
358 fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
359 fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
360 fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
361 fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
362 fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
363 // fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
364 // fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
366 for(Int_t cen=0; cen<6; cen++){
367 fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
368 fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
369 fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
370 fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
371 fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
372 fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
373 fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
374 fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
375 fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
376 fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
377 fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
378 fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
379 fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
381 //pairs from common parents
382 fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
383 fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
384 fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
385 fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
386 fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
387 fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
388 fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
390 //common parent - pi0
391 fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
392 fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
393 fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
394 fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
395 fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
396 fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
397 fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
398 fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
399 fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
404 //Photon contaminations
405 fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
406 fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
407 fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
408 fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
409 fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.));
411 const Int_t nTypes=24 ;
412 char partTypes[nTypes][55] ;
413 snprintf(partTypes[0],55,"hGammaNoPrim") ; //
414 snprintf(partTypes[1],55,"hGammaPhot") ; //
415 snprintf(partTypes[2],55,"hGammaEl") ; //
416 snprintf(partTypes[3],55,"hGammaPi0") ; //
417 snprintf(partTypes[4],55,"hGammaEta") ; //
418 snprintf(partTypes[5],55,"hhGammaOmega") ; //
419 snprintf(partTypes[6],55,"hGammaPipm") ; //
420 snprintf(partTypes[7],55,"hGammaP") ; //
421 snprintf(partTypes[8],55,"hGammaPbar") ; //
422 snprintf(partTypes[9],55,"hGammaN") ; //
423 snprintf(partTypes[10],55,"hGammaNbar") ; //
424 snprintf(partTypes[11],55,"hGammaK0S") ; //
425 snprintf(partTypes[12],55,"hGammaK0L") ; //
426 snprintf(partTypes[13],55,"hGammaKpm") ; //
427 snprintf(partTypes[14],55,"hGammaKstar") ; //
428 snprintf(partTypes[15],55,"hGammaDelta") ; //
429 snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
430 snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
431 snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
432 snprintf(partTypes[19],55,"hGammaPipmEl") ; //
433 snprintf(partTypes[20],55,"hGammaPipmOther") ; //
434 snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
435 snprintf(partTypes[22],55,"hGammaPipmp") ; //
436 snprintf(partTypes[23],55,"hGammaPipmn") ; //
438 const Int_t nPID=12 ;
440 snprintf(cPID[0],25,"All") ;
441 snprintf(cPID[1],25,"Allcore") ;
442 snprintf(cPID[2],25,"CPV") ;
443 snprintf(cPID[3],25,"CPVcore") ;
444 snprintf(cPID[4],25,"CPV2") ;
445 snprintf(cPID[5],25,"CPV2core") ;
446 snprintf(cPID[6],25,"Disp") ;
447 snprintf(cPID[7],25,"Dispcore") ;
448 snprintf(cPID[8],25,"Disp2") ;
449 snprintf(cPID[9],25,"Disp2core") ;
450 snprintf(cPID[10],25,"Both") ;
451 snprintf(cPID[11],25,"Bothcore") ;
453 for(Int_t itype=0; itype<nTypes; itype++){
454 for(Int_t iPID=0; iPID<nPID; iPID++){
455 for(Int_t cen=0; cen<5; cen++){
456 fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
463 PostData(1, fOutputContainer);
465 printf("Init done \n") ;
469 //________________________________________________________________________
470 void AliPHOSHijingEfficiency::UserExec(Option_t *)
472 // Main loop, called for each event
475 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
476 if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
477 fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
481 AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
483 Printf("ERROR: Could not retrieve event");
484 PostData(1, fOutputContainer);
488 fRunNumber=ConvertRunNumber(event->GetRunNumber()) ;
489 FillHistogram("hSelEvents",0.5,fRunNumber-0.5) ;
490 FillHistogram("hTotSelEvents",0.5) ;
492 // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
493 if(fEventCounter == 0) {
494 // Initialize the PHOS geometry
495 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
497 //We have to apply re-calibration for pass1 LCH10h
498 // Initialize decalibration factors in the form of the OCDB object
499 const Float_t wDecalib = 0.055;
501 fPHOSCalibData = new AliPHOSCalibData();
502 fPHOSCalibData->SetName("PhosCalibData");
504 for(Int_t module=1; module<=5; module++) {
505 for(Int_t column=1; column<=56; column++) {
506 for(Int_t row=1; row<=64; row++) {
507 Float_t adcChannelEmc = gRandom->Gaus(1.,wDecalib);
508 fPHOSCalibData->SetADCchannelEmc(module,column,row,adcChannelEmc);
512 for(Int_t mod=0; mod<5; mod++) {
513 if(!event->GetPHOSMatrix(mod)) continue;
514 fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
515 Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
522 // Checks if we have a primary vertex
523 // Get primary vertices form ESD
524 const AliESDVertex *esdVertex5 = event->GetPrimaryVertex();
526 Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0)
527 Double_t vtx5[3] ={0.,0.,0.};
529 vtx5[0] = esdVertex5->GetX();
530 vtx5[1] = esdVertex5->GetY();
531 vtx5[2] = esdVertex5->GetZ();
534 FillHistogram("hZvertex",esdVertex5->GetZ(),fRunNumber-0.5);
535 if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
536 PostData(1, fOutputContainer);
539 FillHistogram("hSelEvents",1.5,fRunNumber-0.5) ;
540 FillHistogram("hTotSelEvents",1.5) ;
543 if(event->IsPileupFromSPD()){
544 PostData(1, fOutputContainer);
548 FillHistogram("hSelEvents",2.5,fRunNumber-0.5) ;
549 FillHistogram("hTotSelEvents",2.5) ;
554 Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
558 //No dependence on zVtx observed, save memory
561 AliCentrality *centrality = event->GetCentrality();
562 fCentrality=centrality->GetCentralityPercentile("V0M");
564 if( fCentrality <= 0. || fCentrality>80. ){
565 PostData(1, fOutputContainer);
569 FillHistogram("hSelEvents",3.5,fRunNumber-0.5) ;
570 FillHistogram("hTotSelEvents",3.5) ;
574 else if(fCentrality<10.)
576 else if(fCentrality<20.)
578 else if(fCentrality<40.)
580 else if(fCentrality<60.)
582 else if(fCentrality<80.)
585 PostData(1, fOutputContainer);
591 FillHistogram("hSelEvents",4.5,fRunNumber-0.5) ;
592 FillHistogram("hTotSelEvents",4.5) ;
593 //All event selections done
594 FillHistogram("hCentrality",fCentrality,fRunNumber-0.5) ;
595 //Reaction plane is defined in the range (0;pi)
597 Double_t averageRP = 0.;//fRPV0A+fRPV0C+fRP ;
598 Int_t irp=Int_t(10.*(averageRP)/TMath::Pi());
601 if(!fPHOSEvents[zvtx][fCenBin][irp])
602 fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
603 TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
608 fPHOSEvent->Clear() ;
610 fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
613 const Double_t logWeight=4.5 ;
615 TVector3 vertex(vtx5);
617 Int_t multClust = event->GetNumberOfCaloClusters();
618 Int_t inPHOS=0,inMod1=0, inMod2=0, inMod3=0 ;
619 Double_t avrgEm1=0,avrgEm2=0,avrgEm3=0; //average cluster energy
621 AliESDCaloCells * cells = event->GetPHOSCells() ;
622 FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ;
623 FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ;
627 for (Int_t i=0; i<multClust; i++) {
628 AliESDCaloCluster *clu = event->GetCaloCluster(i);
629 if ( !clu->IsPHOS() || clu->E()<0.3) continue;
632 clu->GetPosition(position);
633 TVector3 global(position) ;
635 fPHOSGeo->GlobalPos2RelId(global,relId) ;
636 Int_t mod = relId[0] ;
637 Int_t cellX = relId[2];
638 Int_t cellZ = relId[3] ;
639 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
643 FillHistogram("hCluEvsClu",clu->E(),clu->GetNCells()) ;
644 FillHistogram("hCluEvsCluM",clu->E(),clu->GetM02()) ;
646 //Apply re-Calibreation
647 AliPHOSEsdCluster cluPHOS1(*clu);
648 cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
649 Reclusterize(&cluPHOS1) ;
650 cluPHOS1.EvalAll(logWeight,vertex); // recalculate the cluster parameters
651 cluPHOS1.SetE(rnlin(cluPHOS1.E()));// Users's nonlinearity
652 if(cluPHOS1.E()<0.3) continue;
653 if(clu->GetNCells()<3) continue ;
654 if(clu->GetM02()<0.2) continue ;
658 avrgEm1+=cluPHOS1.E() ;
662 avrgEm2+=cluPHOS1.E() ;
666 avrgEm3+=cluPHOS1.E() ;
671 cluPHOS1.GetMomentum(pv1 ,vtx0);
673 Double_t ecore=CoreEnergy(&cluPHOS1) ;
675 FillHistogram(Form("hCluLowM%d",mod),cellX,cellZ,1.);
676 if(cluPHOS1.E()>1.5){
677 FillHistogram(Form("hCluHighM%d",mod),cellX,cellZ,1.);
680 if(inPHOS>=fPHOSEvent->GetSize()){
681 fPHOSEvent->Expand(inPHOS+50) ;
683 new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
684 AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
686 pv1*= ecore/pv1.E() ;
688 ph->SetNCells(clu->GetNCells());
689 Double_t m02=0.,m20=0.;
690 EvalLambdas(&cluPHOS1,0,m02, m20);
691 ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
692 EvalLambdas(&cluPHOS1,1,m02, m20);
693 ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
695 Double_t distBC=clu->GetDistanceToBadChannel();
697 FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ;
699 FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ;
701 FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ;
705 FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
708 Double_t dx=clu->GetTrackDx() ;
709 Double_t dz=clu->GetTrackDz() ;
711 Bool_t cpvBit=kTRUE ; //No track matched by default
712 Bool_t cpvBit2=kTRUE ; //More Strict criterion
713 TArrayI * itracks = clu->GetTracksMatched() ;
714 if(itracks->GetSize()>0){
715 Int_t iTr = itracks->At(0);
716 if(iTr>=0 && iTr<event->GetNumberOfTracks()){
717 AliESDtrack *track = event->GetTrack(iTr) ;
718 Double_t pt = track->Pt() ;
719 Short_t charge = track->Charge() ;
720 Double_t r=TestCPV(dx, dz, pt,charge) ;
723 FillHistogram("hCPVr",r,pv1.E(),fCentrality);
726 ph->SetCPVBit(cpvBit) ;
727 ph->SetCPV2Bit(cpvBit2) ;
729 FillHistogram(Form("hCluVetoM%d",mod),cellX,cellZ,1.);
731 ph->SetEMCx(float(cellX)) ;
732 ph->SetEMCz(float(cellZ)) ;
733 ph->SetLambdas(cluPHOS1.GetM20(),cluPHOS1.GetM02()) ;
734 ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolded
736 Int_t primary=FindPrimary(clu,sure) ; //номер праймари частицы в стеке
737 ph->SetPrimary(primary) ;
738 ph->SetWeight(PrimaryWeight(primary)) ;
742 FillHistogram("hCenPHOS",fCentrality,inPHOS) ;
744 PostData(1, fOutputContainer);
750 for (Int_t i1=0; i1<inPHOS; i1++) {
751 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
752 Double_t w=ph1->GetWeight() ;
753 Double_t pt = ph1->Pt() ;
754 Double_t ptv =ph1->GetMomV2()->Pt() ;
756 FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ;
757 FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptv,w) ;
758 if(ph1->IsntUnfolded()){
759 FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ;
762 FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ;
763 FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptv,w) ;
766 FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ;
769 FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ;
771 FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptv,w) ;
772 if(ph1->IsntUnfolded()){
773 FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ;
776 FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ;
777 FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptv,w) ;
780 if(ph1->IsDisp2OK()){
781 FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ;
782 FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptv,w) ;
784 FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ;
785 FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptv,w) ;
791 const Double_t alphaCut=0.7 ;
793 for (Int_t i1=0; i1<inPHOS-1; i1++) {
794 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
795 Double_t w1=ph1->GetWeight() ;
796 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
797 AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
798 Double_t w2=ph2->GetWeight() ;
800 if(FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary())!=-1)
803 TLorentzVector p12 = *ph1 + *ph2;
804 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
806 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
807 Double_t m = p12.M() ;
808 Double_t pt = p12.Pt() ;
809 Double_t mv = p12.M() ;
810 Double_t ptv = p12.Pt() ;
812 FillHistogram(Form("hPi0All_cen%d",fCenBin),m,pt,w) ;
813 FillHistogram(Form("hPi0Allcore_cen%d",fCenBin),mv,ptv,w) ;
814 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
815 FillHistogram(Form("hPi0Allwou_cen%d",fCenBin),m,pt,w) ;
819 Double_t pt1 = ph1->Pt() ;
820 Double_t ptv1 = ph1->GetMomV2()->Pt() ;
821 Double_t pt2 = ph2->Pt() ;
822 Double_t ptv2 = ph2->GetMomV2()->Pt() ;
823 FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt1,w) ;
824 FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt2,w) ;
825 FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv1,w) ;
826 FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv2,w) ;
827 if(ph1->IsntUnfolded())
828 FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt1,w) ;
829 if(ph2->IsntUnfolded())
830 FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt2,w) ;
832 FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt1,w) ;
833 FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv1,w) ;
836 FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt2,w) ;
837 FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv2,w) ;
840 FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt1,w) ;
843 FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt2,w) ;
846 FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt1,w) ;
847 if(ph1->IsntUnfolded()){
848 FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt1,w) ;
850 FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv1,w) ;
853 FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt2,w) ;
854 if(ph1->IsntUnfolded()){
855 FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt2,w) ;
857 FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv2,w) ;
859 if(ph1->IsDisp2OK()){
860 FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt1,w) ;
861 FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv1,w) ;
863 if(ph2->IsDisp2OK()){
864 FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt2,w) ;
865 FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv2,w) ;
867 if(ph1->IsDispOK() && ph1->IsCPVOK()){
868 FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt1,w) ;
869 FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv1,w) ;
871 if(ph2->IsDispOK() && ph2->IsCPVOK()){
872 FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt2,w) ;
873 FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv2,w) ;
875 if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
876 FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt1,w) ;
877 FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv1,w) ;
879 if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
880 FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt2,w) ;
881 FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv2,w) ;
886 FillHistogram(Form("hPi0All_a07_cen%d",fCenBin),m,pt,w) ;
889 if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 CPV
890 FillHistogram(Form("hPi0CPV_cen%d",fCenBin),m,pt,w) ;
891 FillHistogram(Form("hPi0CPVcore_cen%d",fCenBin),mv,ptv,w) ;
894 FillHistogram(Form("hPi0CPV_a07_cen%d",fCenBin),m,pt,w) ;
897 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
898 FillHistogram(Form("hPi0CPV2_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
900 FillHistogram(Form("hPi0CPV2_a07_cen%d",fCenBin),m,pt,w) ;
903 if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
904 FillHistogram(Form("hPi0Disp2_cen%d",fCenBin),m,pt,w) ;
905 FillHistogram(Form("hPi0Disp2core_cen%d",fCenBin),mv,ptv,w) ;
906 if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
907 FillHistogram(Form("hPi0Both2_cen%d",fCenBin),m,pt,w) ;
908 FillHistogram(Form("hPi0Both2core_cen%d",fCenBin),mv,ptv,w) ;
910 FillHistogram(Form("hPi0Both2_a07_cen%d",fCenBin),m,pt,w) ;
916 if(ph1->IsDispOK() && ph2->IsDispOK()){ //pi0 Disp
917 FillHistogram(Form("hPi0Disp_cen%d",fCenBin),m,pt,w) ;
918 FillHistogram(Form("hPi0Dispcore_cen%d",fCenBin),mv,ptv,w) ;
920 FillHistogram(Form("hPi0Disp_a07_cen%d",fCenBin),m,pt,w) ;
923 if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
924 FillHistogram(Form("hPi0Both_cen%d",fCenBin),m,pt,w) ;
925 FillHistogram(Form("hPi0Bothcore_cen%d",fCenBin),mv,ptv,w) ;
928 FillHistogram(Form("hPi0Both_a07_cen%d",fCenBin),m,pt,w) ;
936 for (Int_t i1=0; i1<inPHOS; i1++) {
937 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
938 Double_t w1=ph1->GetWeight() ;
939 for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
940 TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
941 for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
942 AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
943 Double_t w2=ph2->GetWeight() ;
944 TLorentzVector p12 = *ph1 + *ph2;
945 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
947 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
948 Double_t m = p12.M() ;
949 Double_t pt = p12.Pt() ;
950 Double_t mv = p12.M() ;
951 Double_t ptv = p12.Pt() ;
954 FillHistogram(Form("hMiPi0All_cen%d",fCenBin),m,pt,w) ;
955 FillHistogram(Form("hMiPi0Allcore_cen%d",fCenBin),mv,ptv,w) ;
956 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
957 FillHistogram(Form("hMiPi0Allwou_cen%d",fCenBin),m,pt,w) ;
960 Double_t pt1 = ph1->Pt() ;
961 Double_t ptv1 = ph1->GetMomV2()->Pt() ;
962 Double_t pt2 = ph2->Pt() ;
963 Double_t ptv2 = ph2->GetMomV2()->Pt() ;
964 FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt1,w) ;
965 FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt2,w) ;
966 FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv1,w) ;
967 FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv2,w) ;
968 if(ph1->IsntUnfolded())
969 FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt1,w) ;
970 if(ph2->IsntUnfolded())
971 FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt2,w) ;
973 FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt1,w) ;
974 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv1,w) ;
977 FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt2,w) ;
978 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv2,w) ;
981 FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt1,w) ;
984 FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt2,w) ;
987 FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt1,w) ;
988 if(ph1->IsntUnfolded()){
989 FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt1,w) ;
991 FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv1,w) ;
994 FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt2,w) ;
995 if(ph2->IsntUnfolded()){
996 FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt2,w) ;
998 FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv2,w) ;
1000 if(ph1->IsDisp2OK()){
1001 FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt1,w) ;
1002 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv1,w) ;
1004 if(ph2->IsDisp2OK()){
1005 FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt2,w) ;
1006 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv2,w) ;
1008 if(ph1->IsDispOK() && ph1->IsCPVOK()){
1009 FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt1,w) ;
1010 FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv1,w) ;
1012 if(ph2->IsDispOK() && ph2->IsCPVOK()){
1013 FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt2,w) ;
1014 FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv2,w) ;
1016 if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
1017 FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt1,w) ;
1018 FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv1,w) ;
1020 if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
1021 FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt2,w) ;
1022 FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv2,w) ;
1028 FillHistogram(Form("hMiPi0All_a07_cen%d",fCenBin),m,pt,w) ;
1030 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1031 FillHistogram(Form("hMiPi0CPV_cen%d",fCenBin),m,pt,w) ;
1032 FillHistogram(Form("hMiPi0CPVcore_cen%d",fCenBin),mv,ptv,w) ;
1035 FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCenBin),m,pt,w) ;
1038 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1039 FillHistogram(Form("hMiPi0CPV2_cen%d",fCenBin),m,pt,w) ;
1042 FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCenBin),m,pt,w) ;
1045 if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1046 FillHistogram(Form("hMiPi0Disp2_cen%d",fCenBin),m,pt,w) ;
1047 FillHistogram(Form("hMiPi0Disp2core_cen%d",fCenBin),mv,ptv,w) ;
1048 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1049 FillHistogram(Form("hMiPi0Both2_cen%d",fCenBin),m,pt,w) ;
1050 FillHistogram(Form("hMiPi0Both2core_cen%d",fCenBin),mv,ptv,w) ;
1053 if(ph1->IsDispOK() && ph2->IsDispOK()){
1055 FillHistogram(Form("hMiPi0Disp_cen%d",fCenBin),m,pt,w) ;
1056 FillHistogram(Form("hMiPi0Dispcore_cen%d",fCenBin),mv,ptv,w) ;
1057 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1058 FillHistogram(Form("hMiPi0Dispwou_cen%d",fCenBin),m,pt,w) ;
1062 FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCenBin),m,pt,w) ;
1064 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1065 FillHistogram(Form("hMiPi0Both_cen%d",fCenBin),m,pt,w) ;
1066 FillHistogram(Form("hMiPi0Bothcore_cen%d",fCenBin),mv,ptv,w) ;
1069 FillHistogram(Form("hMiPi0Both_a07_cen%d",fCenBin),m,pt,w) ;
1081 //Now we either add current events to stack or remove
1082 //If no photons in current event - no need to add it to mixed
1083 const Int_t kMixEvents[6]={5,5,5,10,10,30} ;
1084 if(fPHOSEvent->GetEntriesFast()>0){
1085 prevPHOS->AddFirst(fPHOSEvent) ;
1087 if(prevPHOS->GetSize()>kMixEvents[fCenBin]){//Remove redundant events
1088 TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1089 prevPHOS->RemoveLast() ;
1093 // Post output data.
1094 PostData(1, fOutputContainer);
1097 //________________________________________________________________________
1098 void AliPHOSHijingEfficiency::FillSecondaries(){
1101 const Double_t rcut = 1. ;
1102 //Fill spectra of primary particles
1103 //with proper weight
1104 for(Int_t i=0; i<fStack->GetNtrack(); i++){
1105 TParticle * p = fStack->Particle(i) ;
1108 if(TMath::Abs(p->Y())>0.5)
1110 Double_t w = PrimaryParticleWeight(p) ;
1111 Int_t primPdgCode=p->GetPdgCode() ;
1112 switch(primPdgCode){
1113 case 22: FillHistogram(Form("hPrimPhot_cen%d",fCenBin),p->Pt(),w);
1117 FillHistogram(Form("hPrimEl_cen%d",fCenBin),p->Pt(),w);
1120 FillHistogram(Form("hPrimPi0_cen%d",fCenBin),p->Pt(),w);
1123 FillHistogram(Form("hPrimEta_cen%d",fCenBin),p->Pt(),w);
1127 FillHistogram(Form("hPrimPipm_cen%d",fCenBin),p->Pt(),w);
1130 FillHistogram(Form("hPrimP_cen%d",fCenBin),p->Pt(),w);
1133 FillHistogram(Form("hPrimPbar_cen%d",fCenBin),p->Pt(),w);
1136 FillHistogram(Form("hPrimN_cen%d",fCenBin),p->Pt(),w);
1139 FillHistogram(Form("hPrimNbar_cen%d",fCenBin),p->Pt(),w);
1142 FillHistogram(Form("hPrimK0S_cen%d",fCenBin),p->Pt(),w);
1145 FillHistogram(Form("hPrimK0L_cen%d",fCenBin),p->Pt(),w);
1149 FillHistogram(Form("hPrimKpm_cen%d",fCenBin),p->Pt(),w);
1152 FillHistogram(Form("hPrimOther_cen%d",fCenBin),p->Pt(),w);
1156 //Origins of secondary pi0s
1157 for(Int_t i=0; i<fStack->GetNtrack(); i++){
1158 TParticle * p = fStack->Particle(i) ;
1159 if(p->GetPdgCode()!=111)
1161 FillHistogram("Vertex",p->Pt(),p->R());
1164 Double_t phi=p->Phi() ;
1165 while(phi<0.)phi+=TMath::TwoPi() ;
1166 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1167 FillHistogram("hSecondPi0RphiZ",p->R(),phi,p->Vz()) ;
1168 Double_t w = PrimaryParticleWeight(p) ;
1169 FillHistogram("hSecondPi0RE",p->R(),p->Pt(),w) ;
1174 Int_t inPHOS=fPHOSEvent->GetEntries() ;
1175 for (Int_t i1=0; i1<inPHOS-1; i1++) {
1176 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
1177 Double_t w1=ph1->GetWeight() ;
1178 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1179 AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
1180 TLorentzVector p12 = *ph1 + *ph2;
1181 Double_t w2=ph2->GetWeight() ;
1182 Double_t w = TMath::Sqrt(w1*w2) ;
1183 FillHistogram(Form("hParentAll_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1184 Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ;
1186 TParticle * particle = (TParticle *)fStack->Particle(prim);
1187 FillHistogram("hMass_R",p12.M(),p12.Pt(),TMath::Sqrt(particle->R()*particle->R()+particle->Vz()*particle->Vz())) ;
1190 Int_t pdgCode=particle->GetPdgCode() ;
1191 if(pdgCode!=111){ //common parent - not pi0
1193 FillHistogram(Form("hParentGamma_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1195 if(pdgCode==11 || pdgCode==-11){
1196 FillHistogram(Form("hParentEl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1199 if(InPi0mass(p12.M() ,p12.Pt())){
1200 printf("Common parent: %d \n",pdgCode) ;
1202 FillHistogram(Form("hParentOther_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1205 }//Common parent not pi0
1206 else{ //common parent - pi0
1207 FillHistogram(Form("hParentPi0_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1208 FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),particle->R(),w) ;
1209 FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Vz(),w) ;
1210 if(particle->R()<rcut && TMath::Abs(particle->Vz())<10.){
1211 FillHistogram(Form("hParentDirPi0_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1214 //Common particle pi0, created off-vertex
1215 Int_t primPi0=particle->GetFirstMother();
1217 FillHistogram(Form("hParentPi0NoPrim_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1220 Int_t primPdgCode=fStack->Particle(primPi0)->GetPdgCode() ;
1221 switch(primPdgCode){
1222 case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //eta
1224 case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //omega
1227 case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1230 case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1232 case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0S
1234 case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0L
1238 FillHistogram(Form("hParentPi0pn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn
1242 FillHistogram(Form("hParentPi0antipn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn
1245 FillHistogram(Form("hParentPi0Other_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1248 }//common parent - pi0
1249 }//there is common primary
1250 }//seond photon loop
1251 }//first photon loop
1254 //Now look at photon contaiminations
1255 for (Int_t i1=0; i1<inPHOS-1; i1++) {
1256 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
1257 Int_t iprim = ph1->GetPrimary() ;
1259 FillAllHistograms(Form("hGammaNoPrim_cen%d",fCenBin),ph1) ; //
1261 //Find primary at vertex
1262 TParticle * primPHOS = fStack->Particle(iprim) ;
1263 Int_t iprimV=primPHOS->GetFirstMother();
1264 TParticle * primVtx = primPHOS ;
1265 while((iprimV>-1) && primVtx->R()>rcut){
1266 primVtx = fStack->Particle(iprimV) ;
1267 iprimV=primVtx->GetFirstMother();
1271 Int_t primPdgCode=primVtx->GetPdgCode() ;
1272 switch(primPdgCode){
1273 case 22: FillAllHistograms("hGammaPhot",ph1);
1277 FillAllHistograms("hGammaEl",ph1);
1280 FillAllHistograms("hGammaPi0",ph1);
1283 FillAllHistograms("hGammaEta",ph1);
1285 case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
1289 FillAllHistograms("hGammaPipm",ph1);
1290 //Find particle entered PHOS
1291 if(primVtx == primPHOS)
1292 FillAllHistograms("hGammaPipmDirect",ph1);
1294 Int_t primPdgPHOS=primPHOS->GetPdgCode() ;
1295 if(primPdgPHOS==22){
1296 FillAllHistograms("hGammaPipmGamma",ph1);
1297 FillHistogram("hPipmGammaConvR",ph1->Pt(),primPHOS->R());
1298 FillHistogram("hPipmGammaConvRZ",primPHOS->Vz(),primPHOS->R());
1301 if(TMath::Abs(primPdgPHOS)==11){
1302 FillAllHistograms("hGammaPipmEl",ph1);
1303 FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R());
1306 if(TMath::Abs(primPdgPHOS)==2212){
1307 FillAllHistograms("hGammaPipmp",ph1);
1308 FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1311 if(TMath::Abs(primPdgPHOS)==2112){
1312 FillAllHistograms("hGammaPipmn",ph1);
1313 FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1316 FillAllHistograms("hGammaPipmOther",ph1);
1317 FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R());
1321 FillAllHistograms("hGammaP",ph1);
1324 FillAllHistograms("hGammaPbar",ph1);
1327 FillAllHistograms("hGammaN",ph1);
1330 FillAllHistograms("hGammaNbar",ph1) ; // pn
1333 FillAllHistograms("hGammaK0S",ph1) ; // pn
1336 FillAllHistograms("hGammaK0L",ph1) ; // pn
1340 FillAllHistograms("hGammaKpm",ph1) ; // pn
1345 case 313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
1348 case -2224 : //Deltas
1349 case 2224 : //Deltas
1350 case -2214 : //Deltas
1351 case 2214 : //Deltas
1352 case -2114 : //Deltas
1353 case 2114 : //Deltas
1354 case -1114 : //Deltas
1355 case 1114 : //Deltas
1356 FillAllHistograms("hGammaDelta",ph1) ; // pn
1359 if(primVtx->GetPDG()->Charge())
1360 FillAllHistograms("hGammaOtherCharged",ph1) ; //
1362 FillAllHistograms("hGammaOtherNeutral",ph1) ; //
1369 //________________________________________________________________________
1370 void AliPHOSHijingEfficiency::Terminate(Option_t *)
1372 // Draw result to the screen
1373 // Called once at the end of the query
1377 //________________________________________________________________________
1378 Bool_t AliPHOSHijingEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
1380 //Check if this channel belogs to the good ones
1382 if(strcmp(det,"PHOS")==0){
1384 AliError(Form("No bad map for PHOS module %d ",mod)) ;
1387 if(!fPHOSBadMap[mod]){
1388 AliError(Form("No Bad map for PHOS module %d",mod)) ;
1391 if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
1397 AliError(Form("Can not find bad channels for detector %s ",det)) ;
1401 //_____________________________________________________________________________
1402 void AliPHOSHijingEfficiency::FillAllHistograms(const char * key,AliCaloPhoton * ph)const{
1403 //Fill All PID histograms
1405 Double_t w=ph->GetWeight() ;
1406 Double_t pt = ph->Pt() ;
1407 Double_t ptC=ph->GetMomV2()->Pt() ;
1408 FillHistogram(Form("%s_All_cen%d",key,fCenBin),pt,w) ;
1409 FillHistogram(Form("%s_Allcore_cen%d",key,fCenBin),ptC,w) ;
1411 FillHistogram(Form("%s_CPV_cen%d",key,fCenBin),pt,w) ;
1412 FillHistogram(Form("%s_CPVcore_cen%d",key,fCenBin),ptC,w) ;
1415 FillHistogram(Form("%s_CPV2_cen%d",key,fCenBin),pt,w) ;
1416 FillHistogram(Form("%s_CPV2core_cen%d",key,fCenBin),ptC,w) ;
1419 FillHistogram(Form("%s_Disp_cen%d",key,fCenBin),pt,w) ;
1420 FillHistogram(Form("%s_Dispcore_cen%d",key,fCenBin),ptC,w) ;
1421 if(ph->IsDisp2OK()){
1422 FillHistogram(Form("%s_Disp2_cen%d",key,fCenBin),pt,w) ;
1423 FillHistogram(Form("%s_Disp2core_cen%d",key,fCenBin),ptC,w) ;
1426 FillHistogram(Form("%s_Both_cen%d",key,fCenBin),pt,w) ;
1427 FillHistogram(Form("%s_Bothcore_cen%d",key,fCenBin),ptC,w) ;
1431 //_____________________________________________________________________________
1432 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x)const{
1434 TObject * tmp = fOutputContainer->FindObject(key) ;
1436 AliInfo(Form("can not find histogram <%s> ",key)) ;
1439 if(tmp->IsA() == TClass::GetClass("TH1I")){
1440 ((TH1I*)tmp)->Fill(x) ;
1443 if(tmp->IsA() == TClass::GetClass("TH1F")){
1444 ((TH1F*)tmp)->Fill(x) ;
1447 if(tmp->IsA() == TClass::GetClass("TH1D")){
1448 ((TH1D*)tmp)->Fill(x) ;
1451 AliInfo(Form("can not find 1D histogram <%s> ",key)) ;
1453 //_____________________________________________________________________________
1454 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
1456 TObject * tmp = fOutputContainer->FindObject(key) ;
1458 AliInfo(Form("can not find histogram <%s> ",key)) ;
1461 if(tmp->IsA() == TClass::GetClass("TH1F")){
1462 ((TH1F*)tmp)->Fill(x,y) ;
1465 if(tmp->IsA() == TClass::GetClass("TH2F")){
1466 ((TH2F*)tmp)->Fill(x,y) ;
1469 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1472 //_____________________________________________________________________________
1473 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1474 //Fills 1D histograms with key
1475 TObject * tmp = fOutputContainer->FindObject(key) ;
1477 AliInfo(Form("can not find histogram <%s> ",key)) ;
1480 if(tmp->IsA() == TClass::GetClass("TH2F")){
1481 ((TH2F*)tmp)->Fill(x,y,z) ;
1484 if(tmp->IsA() == TClass::GetClass("TH3F")){
1485 ((TH3F*)tmp)->Fill(x,y,z) ;
1489 //_____________________________________________________________________________
1490 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{
1491 //Fills 1D histograms with key
1492 TObject * tmp = fOutputContainer->FindObject(key) ;
1494 AliInfo(Form("can not find histogram <%s> ",key)) ;
1497 if(tmp->IsA() == TClass::GetClass("TH3F")){
1498 ((TH3F*)tmp)->Fill(x,y,z,w) ;
1503 //___________________________________________________________________________
1504 Int_t AliPHOSHijingEfficiency::ConvertRunNumber(Int_t run){
1507 case 139517 : return 137;
1508 case 139514 : return 136;
1509 case 139513 : return 135;
1510 case 139511 : return 134;
1511 case 139510 : return 133;
1512 case 139507 : return 132;
1513 case 139505 : return 131;
1514 case 139504 : return 130;
1515 case 139503 : return 129;
1516 case 139470 : return 128;
1517 case 139467 : return 127;
1518 case 139466 : return 126;
1519 case 139465 : return 125;
1520 case 139440 : return 124;
1521 case 139439 : return 123;
1522 case 139438 : return 122;
1523 case 139437 : return 121;
1524 case 139360 : return 120;
1525 case 139329 : return 119;
1526 case 139328 : return 118;
1527 case 139314 : return 117;
1528 case 139311 : return 116;
1529 case 139310 : return 115;
1530 case 139309 : return 114;
1531 case 139308 : return 113;
1532 case 139173 : return 112;
1533 case 139172 : return 111;
1534 case 139110 : return 110;
1535 case 139107 : return 109;
1536 case 139105 : return 108;
1537 case 139104 : return 107;
1538 case 139042 : return 106;
1539 case 139038 : return 105;
1540 case 139037 : return 104;
1541 case 139036 : return 103;
1542 case 139029 : return 102;
1543 case 139028 : return 101;
1544 case 138983 : return 100;
1545 case 138982 : return 99;
1546 case 138980 : return 98;
1547 case 138979 : return 97;
1548 case 138978 : return 96;
1549 case 138977 : return 95;
1550 case 138976 : return 94;
1551 case 138973 : return 93;
1552 case 138972 : return 92;
1553 case 138965 : return 91;
1554 case 138924 : return 90;
1555 case 138872 : return 89;
1556 case 138871 : return 88;
1557 case 138870 : return 87;
1558 case 138837 : return 86;
1559 case 138830 : return 85;
1560 case 138828 : return 84;
1561 case 138826 : return 83;
1562 case 138796 : return 82;
1563 case 138795 : return 81;
1564 case 138742 : return 80;
1565 case 138732 : return 79;
1566 case 138730 : return 78;
1567 case 138666 : return 77;
1568 case 138662 : return 76;
1569 case 138653 : return 75;
1570 case 138652 : return 74;
1571 case 138638 : return 73;
1572 case 138624 : return 72;
1573 case 138621 : return 71;
1574 case 138583 : return 70;
1575 case 138582 : return 69;
1576 case 138579 : return 68;
1577 case 138578 : return 67;
1578 case 138534 : return 66;
1579 case 138469 : return 65;
1580 case 138442 : return 64;
1581 case 138439 : return 63;
1582 case 138438 : return 62;
1583 case 138396 : return 61;
1584 case 138364 : return 60;
1585 case 138359 : return 59;
1586 case 138275 : return 58;
1587 case 138225 : return 57;
1588 case 138201 : return 56;
1589 case 138200 : return 55;
1590 case 138197 : return 54;
1591 case 138192 : return 53;
1592 case 138190 : return 52;
1593 case 138154 : return 51;
1594 case 138153 : return 50;
1595 case 138151 : return 49;
1596 case 138150 : return 48;
1597 case 138126 : return 47;
1598 case 138125 : return 46;
1599 case 137848 : return 45;
1600 case 137847 : return 44;
1601 case 137844 : return 43;
1602 case 137843 : return 42;
1603 case 137752 : return 41;
1604 case 137751 : return 40;
1605 case 137748 : return 39;
1606 case 137724 : return 38;
1607 case 137722 : return 37;
1608 case 137718 : return 36;
1609 case 137704 : return 35;
1610 case 137693 : return 34;
1611 case 137692 : return 33;
1612 case 137691 : return 32;
1613 case 137689 : return 31;
1614 case 137686 : return 30;
1615 case 137685 : return 29;
1616 case 137639 : return 28;
1617 case 137638 : return 27;
1618 case 137608 : return 26;
1619 case 137595 : return 25;
1620 case 137549 : return 24;
1621 case 137546 : return 23;
1622 case 137544 : return 22;
1623 case 137541 : return 21;
1624 case 137539 : return 20;
1625 case 137531 : return 19;
1626 case 137530 : return 18;
1627 case 137443 : return 17;
1628 case 137441 : return 16;
1629 case 137440 : return 15;
1630 case 137439 : return 14;
1631 case 137434 : return 13;
1632 case 137432 : return 12;
1633 case 137431 : return 11;
1634 case 137430 : return 10;
1635 case 137366 : return 9;
1636 case 137243 : return 8;
1637 case 137236 : return 7;
1638 case 137235 : return 6;
1639 case 137232 : return 5;
1640 case 137231 : return 4;
1641 case 137165 : return 3;
1642 case 137162 : return 2;
1643 case 137161 : return 1;
1644 default : return 199;
1648 //_____________________________________________________________________________
1649 Bool_t AliPHOSHijingEfficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1652 Double_t l1Mean = 1.170014 -0.059465/(1.+0.019343*pt+0.147544*pt*pt) ;
1653 Double_t l2Mean = 1.626270 + 0.761554*exp(-1.213839*pt)-0.020027*pt ;
1654 Double_t l1Sigma = 0.133409 + 0.261307*exp(-0.636874*pt)-0.002849*pt ;
1655 Double_t l2Sigma = 0.289698 + 0.459400*exp(-1.214242*pt)-0.012578*pt ;
1656 Double_t c=-0.124103 ;
1658 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1659 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1660 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1661 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1662 Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1664 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1665 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1666 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1667 return (R2<2.5*2.5) ;
1670 //_____________________________________________________________________________
1671 Bool_t AliPHOSHijingEfficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1674 Double_t l1Mean = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
1675 Double_t l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
1676 Double_t l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
1677 Double_t l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
1678 Double_t c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
1681 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1682 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1683 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1684 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1685 Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1687 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1688 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1689 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1690 return (R2<2.5*2.5) ;
1693 //____________________________________________________________________________
1694 Double_t AliPHOSHijingEfficiency::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1695 //Parameterization of LHC10h period
1700 Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1701 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1702 Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1703 AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
1704 if(!event) //not possible, but required by Coverity
1706 Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1708 if(mf<0.){ //field --
1711 meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1713 meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1718 meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1720 meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1723 Double_t rz=(dz-meanZ)/sz ;
1724 Double_t rx=(dx-meanX)/sx ;
1725 return TMath::Sqrt(rx*rx+rz*rz) ;
1727 //____________________________________________________________________________
1728 Bool_t AliPHOSHijingEfficiency::TestTOF(Double_t t, Double_t e){
1730 Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09
1731 +2.24611e-09*2.24611e-09/e
1732 +5.65054e-09*5.65054e-09/e/e) ;
1733 sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part
1734 Double_t mean=1.1e-9 ;
1735 if(TMath::Abs(t-mean)<2.*sigma)
1738 if(TMath::Abs(t-mean+100.e-9)<2.*sigma)
1743 //____________________________________________________________________________
1744 Double_t AliPHOSHijingEfficiency::CoreEnergy(AliESDCaloCluster * clu){
1745 //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1747 //Can not use already calculated coordinates?
1748 //They have incidence correction...
1749 const Double_t distanceCut =3.5 ;
1750 const Double_t logWeight=4.5 ;
1752 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1753 // Calculates the center of gravity in the local PHOS-module coordinates
1755 Double_t xc[100]={0} ;
1756 Double_t zc[100]={0} ;
1759 Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1760 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1764 fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1765 fPHOSGeo->RelPosInModule(relid, xi, zi);
1768 if (clu->E()>0 && elist[iDigit]>0) {
1769 Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1770 x += xc[iDigit] * w ;
1771 z += zc[iDigit] * w ;
1780 for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1781 Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1782 if(distance < distanceCut)
1783 coreE += elist[iDigit] ;
1785 //Apply non-linearity correction
1786 return rnlin(coreE) ;
1788 //____________________________________________________________________________
1789 Bool_t AliPHOSHijingEfficiency::AreNeibors(Int_t id1,Int_t id2){
1792 fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
1795 fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
1797 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
1798 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
1799 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
1801 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
1802 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
1806 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
1807 return 0; // Difference in row numbers is too large to look further
1815 //____________________________________________________________________________
1816 void AliPHOSHijingEfficiency::Reclusterize(AliESDCaloCluster * clu){
1817 //Re-clusterize to make continues cluster
1819 const Int_t oldMulDigit=clu->GetNCells() ;
1820 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1821 UShort_t * dlist = clu->GetCellsAbsId();
1823 Int_t index[oldMulDigit] ;
1824 Bool_t used[oldMulDigit] ;
1825 for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
1829 for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
1830 if(eMax<elist[iDigit]){
1836 if(inClu==0){ //empty cluster
1839 used[index[0]]=kTRUE ; //mark as used
1840 for(Int_t i=0; i<inClu; i++){
1841 for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
1842 if(used[iDigit]) //already used
1844 if(AreNeibors(dlist[index[i]],dlist[iDigit])){
1845 index[inClu]= iDigit ;
1847 used[iDigit]=kTRUE ;
1852 if(inClu==oldMulDigit) //no need to modify
1855 clu->SetNCells(inClu);
1857 UShort_t tmpD[oldMulDigit] ;
1858 Double_t tmpE[oldMulDigit] ;
1859 for(Int_t i=0; i<oldMulDigit; i++){
1863 //change order of digits in list so that
1864 //first inClu cells were true ones
1865 for(Int_t i=0; i<inClu; i++){
1866 dlist[i]=tmpD[index[i]] ;
1867 elist[i]=tmpE[index[i]] ;
1873 //____________________________________________________________________________
1874 void AliPHOSHijingEfficiency::EvalLambdas(AliESDCaloCluster * clu, Int_t iR,Double_t &m02, Double_t &m20){
1875 //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
1884 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1885 // Calculates the center of gravity in the local PHOS-module coordinates
1887 Double_t xc[100]={0} ;
1888 Double_t zc[100]={0} ;
1891 Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1892 const Double_t logWeight=4.5 ;
1893 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1897 fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1898 fPHOSGeo->RelPosInModule(relid, xi, zi);
1901 if (clu->E()>0 && elist[iDigit]>0) {
1902 Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1903 x += xc[iDigit] * w ;
1904 z += zc[iDigit] * w ;
1917 Double_t xCut = 0. ;
1918 Double_t zCut = 0. ;
1919 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1920 if (clu->E()>0 && elist[iDigit]>0.) {
1921 Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1922 Double_t xi= xc[iDigit] ;
1923 Double_t zi= zc[iDigit] ;
1924 if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
1927 dxx += w * xi * xi ;
1928 dzz += w * zi * zi ;
1929 dxz += w * xi * zi ;
1941 dxx -= xCut * xCut ;
1942 dzz -= zCut * zCut ;
1943 dxz -= xCut * zCut ;
1945 m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
1946 m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
1953 //___________________________________________________________________________
1954 void AliPHOSHijingEfficiency::ProcessMC(){
1955 //fill histograms for efficiensy etc. calculation
1956 const Double_t rcut = 1. ; //cut for primary particles
1957 //---------First pi0/eta-----------------------------
1961 if(!fStack) return ;
1962 for(Int_t i=0;i<fStack->GetNtrack();i++){
1963 TParticle* particle = fStack->Particle(i);
1964 if(particle->GetPdgCode() == 111)
1965 snprintf(partName,10,"pi0") ;
1967 if(particle->GetPdgCode() == 221)
1968 snprintf(partName,10,"eta") ;
1970 if(particle->GetPdgCode() == 22)
1971 snprintf(partName,10,"gamma") ;
1976 Double_t r=particle->R() ;
1977 Double_t pt = particle->Pt() ;
1978 //Distribution over vertex
1979 FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
1984 //Total number of pi0 with creation radius <1 cm
1985 Double_t w = PrimaryParticleWeight(particle) ;
1986 snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
1987 FillHistogram(hkey,pt,w) ;
1988 if(TMath::Abs(particle->Y())<0.12){
1989 snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
1990 FillHistogram(hkey,pt,w) ;
1993 snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
1994 FillHistogram(hkey,particle->Y(),w) ;
1996 Double_t phi=particle->Phi() ;
1997 while(phi<0.)phi+=TMath::TwoPi() ;
1998 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1999 snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
2000 FillHistogram(hkey,phi,w) ;
2003 //___________________________________________________________________________
2004 Int_t AliPHOSHijingEfficiency::FindPrimary(AliESDCaloCluster*clu, Bool_t&sure){
2005 //Finds primary and estimates if it unique one?
2006 //First check can it be photon/electron
2007 const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
2008 Int_t n=clu->GetNLabels() ;
2009 for(Int_t i=0; i<n; i++){
2010 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
2011 Int_t pdg = p->GetPdgCode() ;
2012 if(pdg==22 || pdg==11 || pdg == -11){
2013 if(p->Energy()>emFraction*clu->E()){
2015 return clu->GetLabelAt(i);
2020 Double_t* Ekin= new Double_t[n] ;
2021 for(Int_t i=0; i<n; i++){
2022 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
2023 Ekin[i]=p->P() ; // estimate of kinetic energy
2024 if(p->GetPdgCode()==-2212 || p->GetPdgCode()==-2112){
2025 Ekin[i]+=1.8 ; //due to annihilation
2029 Double_t eMax=0.,eSubMax=0. ;
2030 for(Int_t i=0; i<n; i++){
2037 if(eSubMax>0.8*eMax)//not obvious primary
2042 return clu->GetLabelAt(iMax) ;
2044 //___________________________________________________________________________
2045 Double_t AliPHOSHijingEfficiency::PrimaryWeight(Int_t primary){
2046 //Check who is the primary and introduce weight to correct primary spectrum
2048 if(primary<0 || primary>=fStack->GetNtrack())
2050 //trace primaries up to IP
2051 TParticle* particle = fStack->Particle(primary);
2052 Double_t r=particle->R() ;
2053 Int_t mother = particle->GetFirstMother() ;
2055 if(r<1. && particle->GetPdgCode()==111)
2057 particle = fStack->Particle(mother);
2058 mother = particle->GetFirstMother() ;
2062 return TMath::Max(0.,PrimaryParticleWeight(particle)) ;
2064 //________________________________________________________________________
2065 Double_t AliPHOSHijingEfficiency::PrimaryParticleWeight(TParticle * particle){
2067 Int_t pdg = particle->GetPdgCode() ;
2069 if(pdg == 111 || TMath::Abs(pdg)==211){
2073 if(TMath::Abs(pdg)<1000){ //Kaon-like
2080 Double_t x = particle->Pt() ;
2082 if(fCenBin==0) //0-5
2083 return (1.662990+1.140890*x-0.192088*x*x)/(1.-0.806630*x+0.304771*x*x)+0.141690*x ;
2084 if(fCenBin==1) //5-10
2085 return (1.474351+0.791492*x-0.066369*x*x)/(1.-0.839338*x+0.317312*x*x)+0.093289*x ;
2086 if(fCenBin==2) //10-20
2087 return (1.174728+0.959681*x-0.137695*x*x)/(1.-0.788873*x+0.299538*x*x)+0.128759*x ;
2088 if(fCenBin==3) //20-40
2089 return (0.927335+0.475349*x+0.004364*x*x)/(1.-0.817966*x+0.309787*x*x)+0.086899*x ;
2090 if(fCenBin==4) //40-60
2091 return (0.676878+0.190680*x+0.077031*x*x)/(1.-0.790623*x+0.305183*x*x)+0.064510*x ;
2092 if(fCenBin==5) //60-80
2093 return (0.684726-0.606262*x+0.409963*x*x)/(1.-1.080061*x+0.456933*x*x)+0.005151*x ;
2096 if(fCenBin==0) //0-5
2097 return (-0.417131+2.253936*x-0.337731*x*x)/(1.-0.909892*x+0.316820*x*x)+0.157312*x ;
2098 if(fCenBin==1) //5-10
2099 return (-0.352275+1.844466*x-0.248598*x*x)/(1.-0.897048*x+0.316462*x*x)+0.132461*x ;
2100 if(fCenBin==2) //10-20
2101 return (-0.475481+1.975108*x-0.336013*x*x)/(1.-0.801028*x+0.276705*x*x)+0.188164*x ;
2102 if(fCenBin==3) //20-40
2103 return (-0.198954+1.068789*x-0.103540*x*x)/(1.-0.848354*x+0.299209*x*x)+0.112939*x ;
2104 if(fCenBin==4) //40-60
2105 return (-0.111052+0.664041*x-0.019717*x*x)/(1.-0.804916*x+0.300779*x*x)+0.095784*x ;
2106 if(fCenBin==5) //0-5
2107 return (0.202788-0.439832*x+0.564585*x*x)/(1.-1.254029*x+0.679444*x*x)+0.016235*x ;
2110 if(fCenBin==0) //0-5
2111 return (-1.312732+2.743568*x-0.375775*x*x)/(1.-0.717533*x+0.164694*x*x)+0.164445*x ;
2112 if(fCenBin==1) //5-10
2113 return (-1.229425+2.585889*x-0.330164*x*x)/(1.-0.715892*x+0.167386*x*x)+0.133085*x ;
2114 if(fCenBin==2) //10-20
2115 return (-1.135677+2.397489*x-0.320355*x*x)/(1.-0.709312*x+0.164350*x*x)+0.146095*x ;
2116 if(fCenBin==3) //20-40
2117 return (-0.889993+1.928263*x-0.220785*x*x)/(1.-0.715991*x+0.174729*x*x)+0.095098*x ;
2118 if(fCenBin==4) //40-60
2119 return (-0.539237+1.329118*x-0.115439*x*x)/(1.-0.722906*x+0.186832*x*x)+0.059267*x ;
2120 if(fCenBin==5) //60-80
2121 return (-0.518126+1.327628*x-0.130881*x*x)/(1.-0.665649*x+0.184300*x*x)+0.081701*x ;
2125 //________________________________________________________________________
2126 Int_t AliPHOSHijingEfficiency::FindCommonParent(Int_t iPart, Int_t jPart){
2127 //check if there is a common parent for particles i and j
2128 // -1: no common parent or wrong iPart/jPart
2130 if(iPart==-1 || iPart>=fStack->GetNtrack() ||
2131 jPart==-1 || jPart>=fStack->GetNtrack()) return -1;
2139 iprim2=((TParticle *)fStack->Particle(iprim2))->GetFirstMother();
2141 iprim1=((TParticle *)fStack->Particle(iprim1))->GetFirstMother();
2145 //________________________________________________________________________
2146 Bool_t AliPHOSHijingEfficiency::HaveParent(Int_t iPart, Int_t pdgParent){
2147 //check if there is a common parent for particles i and j
2148 // -1: no common parent or wrong iPart/jPart
2150 if(iPart==-1 || iPart>=fStack->GetNtrack()) return -1;
2154 TParticle * tmp = fStack->Particle(iprim1) ;
2155 if(tmp->GetPdgCode()==pdgParent)
2157 iprim1=tmp->GetFirstMother();
2161 //________________________________________________________________________
2162 Bool_t AliPHOSHijingEfficiency::InPi0mass(Double_t m, Double_t /*pt*/){
2164 return TMath::Abs(m-0.135)<0.007*2.5 ;