]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_PbPb_MC/AliPHOSHijingEfficiency.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Improved primary selection
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb_MC / AliPHOSHijingEfficiency.cxx
CommitLineData
04c118da 1#include "TChain.h"
2#include "TTree.h"
3#include "TObjArray.h"
4#include "TF1.h"
5#include "TFile.h"
6#include "TH1F.h"
7#include "TH2F.h"
8#include "TH2I.h"
9#include "TH3F.h"
10#include "TParticle.h"
11#include "TCanvas.h"
12#include "TStyle.h"
13#include "TRandom.h"
14#include "TROOT.h"
15#include "THashList.h"
16
17#include "AliAnalysisManager.h"
18#include "AliMCEventHandler.h"
19#include "AliMCEvent.h"
20#include "AliStack.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"
31#include "AliLog.h"
32#include "AliPID.h"
33#include "AliCDBManager.h"
34#include "AliCentrality.h"
35#include "AliESDtrackCuts.h"
36#include "AliEventplane.h"
37#include "TProfile.h"
38#include "AliOADBContainer.h"
39
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
46// Date : 05.02.2013
47
48ClassImp(AliPHOSHijingEfficiency)
49
50//________________________________________________________________________
51Double_t rnlin(Double_t x)
52{
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) ;
56
57}
58
59//________________________________________________________________________
60AliPHOSHijingEfficiency::AliPHOSHijingEfficiency(const char *name)
61: AliAnalysisTaskSE(name),
62 fStack(0x0),
63 fOutputContainer(0x0),
64 fPHOSEvent(0x0),
65 fPHOSCalibData(0x0),
66 fRP(0.),fRPV0A(0.),fRPV0C(0.),fHaveTPCRP(0),
67 fRunNumber(0),
68 fCentrality(0.),
69 fCenBin(0),
70 fPHOSGeo(0),
71 fEventCounter(0)
72{
73 // Constructor
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 ;
78 }
79
80 // Output slots #0 write into a TH1 container
81 DefineOutput(1,TList::Class());
82
83 // Set bad channel map
84 char key[55] ;
85 for(Int_t i=0; i<6; i++){
86 sprintf(key,"PHOS_BadMap_mod%d",i) ;
87 fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
88 }
89
90
91}
92
93//________________________________________________________________________
94void AliPHOSHijingEfficiency::UserCreateOutputObjects()
95{
96
97
98 // Create histograms
99 // Called once
100 const Int_t nRuns=200 ;
101
102 // ESD histograms
103 if(fOutputContainer != NULL){
104 delete fOutputContainer;
105 }
106 fOutputContainer = new THashList();
107 fOutputContainer->SetOwner(kTRUE);
108 PostData(1, fOutputContainer);
109
110 //========QA histograms=======
111
112 //Event selection
113 fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 10,0.,10.,nRuns,0.,float(nRuns))) ;
114 fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 10,0.,10.)) ;
115
116 //vertex distribution
117 fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
118
119 //Centrality
120 fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
121 fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
122 fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
123 fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
124 fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,20.,100,0.,100.)) ;
125 fOutputContainer->Add(new TH2F("hCluEvsCluM","ClusterMult vs E",200,0.,20.,100,0.,20.)) ;
126
127 fOutputContainer->Add(new TH3F("hCPVr","CPV radius",100,0.,20.,100,0.,20.,10,0.,100.));
128
129 //Bad Map
130 fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
131 fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
132 fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
133
134 fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
135 fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
136 fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
137
138 fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
139 fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
140 fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
141
142 fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
143 fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
144 fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
145
146 //Single photon and pi0 spectrum
147 Int_t nPtPhot = 300 ;
148 Double_t ptPhotMax = 30 ;
149 Int_t nM = 500;
150 Double_t mMin = 0.0;
151 Double_t mMax = 1.0;
152
153 //PHOS calibration QA
154 fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
155 fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
156 fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
157 fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
158 fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
159 fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
160
161 char key[55] ;
162 for(Int_t cent=0; cent<6; cent++){
163
164 fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad2_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
165 fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad4_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
166 fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad6_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
167
168 fOutputContainer->Add(new TH1F(Form("hPhotAll_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
169 fOutputContainer->Add(new TH1F(Form("hPhotAllcore_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
170 fOutputContainer->Add(new TH1F(Form("hPhotAllwou_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
171 fOutputContainer->Add(new TH1F(Form("hPhotDisp_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
172 fOutputContainer->Add(new TH1F(Form("hPhotDisp2_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
173 fOutputContainer->Add(new TH1F(Form("hPhotDispcore_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
174 fOutputContainer->Add(new TH1F(Form("hPhotDisp2core_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
175 fOutputContainer->Add(new TH1F(Form("hPhotDispwou_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
176 fOutputContainer->Add(new TH1F(Form("hPhotCPV_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
177 fOutputContainer->Add(new TH1F(Form("hPhotCPVcore_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
178 fOutputContainer->Add(new TH1F(Form("hPhotCPV2_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
179 fOutputContainer->Add(new TH1F(Form("hPhotBoth_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
180 fOutputContainer->Add(new TH1F(Form("hPhotBothcore_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
181 fOutputContainer->Add(new TH1F(Form("hPhotBoth2_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
182 fOutputContainer->Add(new TH1F(Form("hPhotBoth2core_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
183
184
185 fOutputContainer->Add(new TH2F(Form("hPi0All_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
186 fOutputContainer->Add(new TH2F(Form("hPi0Allcore_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
187 fOutputContainer->Add(new TH2F(Form("hPi0Allwou_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
188 fOutputContainer->Add(new TH2F(Form("hPi0Disp_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
189 fOutputContainer->Add(new TH2F(Form("hPi0Disp2_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
190 fOutputContainer->Add(new TH2F(Form("hPi0Dispcore_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
191 fOutputContainer->Add(new TH2F(Form("hPi0Disp2core_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
192 fOutputContainer->Add(new TH2F(Form("hPi0Dispwou_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
193 fOutputContainer->Add(new TH2F(Form("hPi0CPV_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
194 fOutputContainer->Add(new TH2F(Form("hPi0CPVcore_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
195 fOutputContainer->Add(new TH2F(Form("hPi0CPV2_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
196 fOutputContainer->Add(new TH2F(Form("hPi0Both_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
197 fOutputContainer->Add(new TH2F(Form("hPi0Bothcore_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
198 fOutputContainer->Add(new TH2F(Form("hPi0Both2_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
199 fOutputContainer->Add(new TH2F(Form("hPi0Both2core_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
200
201
202 snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
203 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
204 snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
205 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
206 snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
207 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
208 snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
209 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
210 snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
211 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
212 snprintf(key,55,"hPi0Both2_a07_cen%d",cent) ;
213 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
214
215 snprintf(key,55,"hSingleAll_cen%d",cent) ;
216 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
217 snprintf(key,55,"hSingleAllcore_cen%d",cent) ;
218 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
219 snprintf(key,55,"hSingleAllwou_cen%d",cent) ;
220 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
221 snprintf(key,55,"hSingleDisp_cen%d",cent) ;
222 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
223 snprintf(key,55,"hSingleDisp2_cen%d",cent) ;
224 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
225 snprintf(key,55,"hSingleDispcore_cen%d",cent) ;
226 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
227 snprintf(key,55,"hSingleDisp2core_cen%d",cent) ;
228 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
229 snprintf(key,55,"hSingleDispwou_cen%d",cent) ;
230 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
231 snprintf(key,55,"hSingleCPV_cen%d",cent) ;
232 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
233 snprintf(key,55,"hSingleCPVcore_cen%d",cent) ;
234 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
235 snprintf(key,55,"hSingleCPV2_cen%d",cent) ;
236 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
237 snprintf(key,55,"hSingleBoth_cen%d",cent) ;
238 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
239 snprintf(key,55,"hSingleBothcore_cen%d",cent) ;
240 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
241 snprintf(key,55,"hSingleBoth2_cen%d",cent) ;
242 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
243 snprintf(key,55,"hSingleBoth2core_cen%d",cent) ;
244 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
245
246
247 snprintf(key,55,"hMiPi0All_cen%d",cent) ;
248 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
249 snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ;
250 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
251 snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ;
252 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
253 snprintf(key,55,"hMiPi0Disp_cen%d",cent) ;
254 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
255 snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ;
256 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
257 snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ;
258 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
259 snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ;
260 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
261 snprintf(key,55,"hMiPi0Disp2core_cen%d",cent) ;
262 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
263 snprintf(key,55,"hMiPi0CPV_cen%d",cent) ;
264 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
265 snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ;
266 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
267 snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ;
268 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
269 snprintf(key,55,"hMiPi0Both_cen%d",cent) ;
270 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
271 snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ;
272 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
273 snprintf(key,55,"hMiPi0Both2_cen%d",cent) ;
274 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
275 snprintf(key,55,"hMiPi0Both2core_cen%d",cent) ;
276 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
277
278 snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
279 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
280 snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
281 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
282 snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
283 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
284 snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
285 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
286 snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
287 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
288 snprintf(key,55,"hMiPi0Both2_a07_cen%d",cent) ;
289 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
290
291 snprintf(key,55,"hMiSingleAll_cen%d",cent) ;
292 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
293 snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ;
294 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
295 snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ;
296 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
297 snprintf(key,55,"hMiSingleDisp_cen%d",cent) ;
298 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
299 snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ;
300 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
301 snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ;
302 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
303 snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ;
304 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
305 snprintf(key,55,"hMiSingleDisp2core_cen%d",cent) ;
306 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
307 snprintf(key,55,"hMiSingleCPV_cen%d",cent) ;
308 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
309 snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ;
310 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
311 snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ;
312 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
313 snprintf(key,55,"hMiSingleBoth_cen%d",cent) ;
314 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
315 snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ;
316 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
317 snprintf(key,55,"hMiSingleBoth2_cen%d",cent) ;
318 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
319 snprintf(key,55,"hMiSingleBoth2core_cen%d",cent) ;
320 fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
321 }
322
323
324 //MC
325 for(Int_t cent=0; cent<6; cent++){
326 snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
327 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
328 snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
329 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
330 snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
331 fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
332 snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
333 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
334 snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
335 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
336 snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
337 fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
338 snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
339 fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
340 snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
341 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
342 snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
343 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
344 snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
345 fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
346 snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
347 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
348 snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
349 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
350 }
351 fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
352 fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
353 fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
354
355
356 Int_t nPt = 200;
357 Double_t ptMin = 0;
358 Double_t ptMax = 20;
359 fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
360 fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
361 fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
362 fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
363 fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
364 fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
365// fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
366// fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
367
368 for(Int_t cen=0; cen<6; cen++){
369 fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
370 fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
371 fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
372 fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
373 fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
374 fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
375 fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
376 fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
377 fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
378 fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
379 fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
380 fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
381 fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
382
383 //pairs from common parents
384 fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
385 fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
386 fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
387 fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
388 fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
389 fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
390 fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
391
392 //common parent - pi0
393 fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
394 fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
395 fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
396 fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
397 fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
398 fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
399 fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
400 fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
401 fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
402
403 }
404
405
406 //Photon contaminations
407 fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
408 fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
409 fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
410 fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
411 fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.));
412
413 const Int_t nTypes=24 ;
414 char partTypes[nTypes][55] ;
415 snprintf(partTypes[0],55,"hGammaNoPrim") ; //
416 snprintf(partTypes[1],55,"hGammaPhot") ; //
417 snprintf(partTypes[2],55,"hGammaEl") ; //
418 snprintf(partTypes[3],55,"hGammaPi0") ; //
419 snprintf(partTypes[4],55,"hGammaEta") ; //
420 snprintf(partTypes[5],55,"hhGammaOmega") ; //
421 snprintf(partTypes[6],55,"hGammaPipm") ; //
422 snprintf(partTypes[7],55,"hGammaP") ; //
423 snprintf(partTypes[8],55,"hGammaPbar") ; //
424 snprintf(partTypes[9],55,"hGammaN") ; //
425 snprintf(partTypes[10],55,"hGammaNbar") ; //
426 snprintf(partTypes[11],55,"hGammaK0S") ; //
427 snprintf(partTypes[12],55,"hGammaK0L") ; //
428 snprintf(partTypes[13],55,"hGammaKpm") ; //
429 snprintf(partTypes[14],55,"hGammaKstar") ; //
430 snprintf(partTypes[15],55,"hGammaDelta") ; //
431 snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
432 snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
433 snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
434 snprintf(partTypes[19],55,"hGammaPipmEl") ; //
435 snprintf(partTypes[20],55,"hGammaPipmOther") ; //
436 snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
437 snprintf(partTypes[22],55,"hGammaPipmp") ; //
438 snprintf(partTypes[23],55,"hGammaPipmn") ; //
439
440 const Int_t nPID=12 ;
441 char cPID[25][12] ;
442 snprintf(cPID[0],25,"All") ;
443 snprintf(cPID[1],25,"Allcore") ;
444 snprintf(cPID[2],25,"CPV") ;
445 snprintf(cPID[3],25,"CPVcore") ;
446 snprintf(cPID[4],25,"CPV2") ;
447 snprintf(cPID[5],25,"CPV2core") ;
448 snprintf(cPID[6],25,"Disp") ;
449 snprintf(cPID[7],25,"Dispcore") ;
450 snprintf(cPID[8],25,"Disp2") ;
451 snprintf(cPID[9],25,"Disp2core") ;
452 snprintf(cPID[10],25,"Both") ;
453 snprintf(cPID[11],25,"Bothcore") ;
454
455 for(Int_t itype=0; itype<nTypes; itype++){
456 for(Int_t iPID=0; iPID<nPID; iPID++){
457 for(Int_t cen=0; cen<5; cen++){
458 fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
459 }
460 }
461 }
462
463
464
465 PostData(1, fOutputContainer);
466
467 printf("Init done \n") ;
468
469}
470
471//________________________________________________________________________
472void AliPHOSHijingEfficiency::UserExec(Option_t *)
473{
474 // Main loop, called for each event
475 // Analyze ESD/AOD
476 fStack=0 ;
477 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
478 if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
479 fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
480 }
481
482
483 AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
484 if (!event) {
485 Printf("ERROR: Could not retrieve event");
486 PostData(1, fOutputContainer);
487 return;
488 }
489
490 fRunNumber=ConvertRunNumber(event->GetRunNumber()) ;
491 FillHistogram("hSelEvents",0.5,fRunNumber-0.5) ;
492 FillHistogram("hTotSelEvents",0.5) ;
493
494 // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
495 char skey[55] ;
496 if(fEventCounter == 0) {
497 // Initialize the PHOS geometry
498 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
499
500 //We have to apply re-calibration for pass1 LCH10h
501 // Initialize decalibration factors in the form of the OCDB object
502 const Float_t wDecalib = 0.055;
503
504 fPHOSCalibData = new AliPHOSCalibData();
505 fPHOSCalibData->SetName("PhosCalibData");
506
507 for(Int_t module=1; module<=5; module++) {
508 for(Int_t column=1; column<=56; column++) {
509 for(Int_t row=1; row<=64; row++) {
510 Float_t adcChannelEmc = gRandom->Gaus(1.,wDecalib);
511 fPHOSCalibData->SetADCchannelEmc(module,column,row,adcChannelEmc);
512 }
513 }
514 }
515 for(Int_t mod=0; mod<5; mod++) {
516 if(!event->GetPHOSMatrix(mod)) continue;
517 fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
518 Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
519 }
520
521 fEventCounter++ ;
522 }
523
524
525 // Checks if we have a primary vertex
526 // Get primary vertices form ESD
527 const AliESDVertex *esdVertex5 = event->GetPrimaryVertex();
528
529 Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0)
530 Double_t vtx5[3] ={0.,0.,0.};
531
532 vtx5[0] = esdVertex5->GetX();
533 vtx5[1] = esdVertex5->GetY();
534 vtx5[2] = esdVertex5->GetZ();
535
536
537 FillHistogram("hZvertex",esdVertex5->GetZ(),fRunNumber-0.5);
538 if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
539 PostData(1, fOutputContainer);
540 return;
541 }
542 FillHistogram("hSelEvents",1.5,fRunNumber-0.5) ;
543 FillHistogram("hTotSelEvents",1.5) ;
544
545/*
546 if(event->IsPileupFromSPD()){
547 PostData(1, fOutputContainer);
548 return;
549 }
550*/
551 FillHistogram("hSelEvents",2.5,fRunNumber-0.5) ;
552 FillHistogram("hTotSelEvents",2.5) ;
553
554
555 /*
556 //Vtx class z-bin
557 Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
558 if(zvtx<0)zvtx=0 ;
559 if(zvtx>9)zvtx=9 ;
560 */
561 //No dependence on zVtx observed, save memory
562 Int_t zvtx=0 ;
563
564 AliCentrality *centrality = event->GetCentrality();
565 fCentrality=centrality->GetCentralityPercentile("V0M");
566
567 if( fCentrality <= 0. || fCentrality>80. ){
568 PostData(1, fOutputContainer);
569 return;
570 }
571
572 FillHistogram("hSelEvents",3.5,fRunNumber-0.5) ;
573 FillHistogram("hTotSelEvents",3.5) ;
574
575 if(fCentrality<5.)
576 fCenBin=0 ;
577 else if(fCentrality<10.)
578 fCenBin=1 ;
579 else if(fCentrality<20.)
580 fCenBin=2 ;
581 else if(fCentrality<40.)
582 fCenBin=3 ;
583 else if(fCentrality<60.)
584 fCenBin=4 ;
585 else if(fCentrality<80.)
586 fCenBin=5 ;
587 else{
588 PostData(1, fOutputContainer);
589 return;
590 }
591
592 fRP=0. ;
593
594 FillHistogram("hSelEvents",4.5,fRunNumber-0.5) ;
595 FillHistogram("hTotSelEvents",4.5) ;
596 //All event selections done
597 FillHistogram("hCentrality",fCentrality,fRunNumber-0.5) ;
598 //Reaction plane is defined in the range (0;pi)
599 //We have 10 bins
600 Double_t averageRP = 0.;//fRPV0A+fRPV0C+fRP ;
601 Int_t irp=Int_t(10.*(averageRP)/TMath::Pi());
602 if(irp>9)irp=9 ;
603
604 if(!fPHOSEvents[zvtx][fCenBin][irp])
605 fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
606 TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
607
608 ProcessMC() ;
609
610 if(fPHOSEvent)
611 fPHOSEvent->Clear() ;
612 else
613 fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
614
615 //For re-calibration
616 const Double_t logWeight=4.5 ;
617
618 TVector3 vertex(vtx5);
619
620 Int_t multClust = event->GetNumberOfCaloClusters();
621 Int_t inPHOS=0,inMod1=0, inMod2=0, inMod3=0 ;
622 Double_t avrgEm1=0,avrgEm2=0,avrgEm3=0; //average cluster energy
623
624 AliESDCaloCells * cells = event->GetPHOSCells() ;
625 FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ;
626 FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ;
627
628
629 TVector3 localPos ;
630 for (Int_t i=0; i<multClust; i++) {
631 AliESDCaloCluster *clu = event->GetCaloCluster(i);
632 if ( !clu->IsPHOS() || clu->E()<0.3) continue;
633
634 Float_t position[3];
635 clu->GetPosition(position);
636 TVector3 global(position) ;
637 Int_t relId[4] ;
638 fPHOSGeo->GlobalPos2RelId(global,relId) ;
639 Int_t mod = relId[0] ;
640 Int_t cellX = relId[2];
641 Int_t cellZ = relId[3] ;
642 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
643 continue ;
644
645
646 FillHistogram("hCluEvsClu",clu->E(),clu->GetNCells()) ;
647 FillHistogram("hCluEvsCluM",clu->E(),clu->GetM02()) ;
648
649 //Apply re-Calibreation
650 AliPHOSEsdCluster cluPHOS1(*clu);
651 cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
652 Reclusterize(&cluPHOS1) ;
653 cluPHOS1.EvalAll(logWeight,vertex); // recalculate the cluster parameters
654 cluPHOS1.SetE(rnlin(cluPHOS1.E()));// Users's nonlinearity
655 if(cluPHOS1.E()<0.3) continue;
656 if(clu->GetNCells()<3) continue ;
657 if(clu->GetM02()<0.2) continue ;
658
659
660 if(mod==1){
661 avrgEm1+=cluPHOS1.E() ;
662 inMod1++ ;
663 }
664 if(mod==2){
665 avrgEm2+=cluPHOS1.E() ;
666 inMod2++ ;
667 }
668 if(mod==3){
669 avrgEm3+=cluPHOS1.E() ;
670 inMod3++ ;
671 }
672
673 TLorentzVector pv1 ;
674 cluPHOS1.GetMomentum(pv1 ,vtx0);
675
676 Double_t ecore=CoreEnergy(&cluPHOS1) ;
677
678 sprintf(skey,"hCluLowM%d",mod) ;
679 FillHistogram(skey,cellX,cellZ,1.);
680 if(cluPHOS1.E()>1.5){
681 sprintf(skey,"hCluHighM%d",mod) ;
682 FillHistogram(skey,cellX,cellZ,1.);
683 }
684
685 if(inPHOS>=fPHOSEvent->GetSize()){
686 fPHOSEvent->Expand(inPHOS+50) ;
687 }
688 new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
689 AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
690 ph->SetModule(mod) ;
691 pv1*= ecore/pv1.E() ;
692 ph->SetMomV2(&pv1) ;
693 ph->SetNCells(clu->GetNCells());
694 Double_t m02=0.,m20=0.;
695 EvalLambdas(&cluPHOS1,0,m02, m20);
696 ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
697 EvalLambdas(&cluPHOS1,1,m02, m20);
698 ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
699
700 Double_t distBC=clu->GetDistanceToBadChannel();
701 if(distBC>2.)
702 FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ;
703 if(distBC>4.)
704 FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ;
705 if(distBC>6.)
706 FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ;
707
708 if(ph->IsDispOK()){
709 sprintf(skey,"hCluDispM%d",mod) ;
710 FillHistogram(skey,cellX,cellZ,1.);
711 }
712 //Track matching
713 Double_t dx=clu->GetTrackDx() ;
714 Double_t dz=clu->GetTrackDz() ;
715
716 Bool_t cpvBit=kTRUE ; //No track matched by default
717 Bool_t cpvBit2=kTRUE ; //More Strict criterion
718 TArrayI * itracks = clu->GetTracksMatched() ;
719 if(itracks->GetSize()>0){
720 Int_t iTr = itracks->At(0);
721 if(iTr>=0 && iTr<event->GetNumberOfTracks()){
722 AliESDtrack *track = event->GetTrack(iTr) ;
723 Double_t pt = track->Pt() ;
724 Short_t charge = track->Charge() ;
725 Double_t r=TestCPV(dx, dz, pt,charge) ;
726 cpvBit=(r>2.) ;
727 cpvBit2=(r>4.) ;
728 FillHistogram("hCPVr",r,pv1.E(),fCentrality);
729 }
730 }
731 ph->SetCPVBit(cpvBit) ;
732 ph->SetCPV2Bit(cpvBit2) ;
733 if(cpvBit){
734 sprintf(skey,"hCluVetoM%d",mod) ;
735 FillHistogram(skey,cellX,cellZ,1.);
736 }
737 ph->SetEMCx(float(cellX)) ;
738 ph->SetEMCz(float(cellZ)) ;
739 ph->SetLambdas(cluPHOS1.GetM20(),cluPHOS1.GetM02()) ;
740 ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolded
741 Bool_t sure= kTRUE;
742 Int_t primary=FindPrimary(clu,sure) ; //номер праймари частицы в стеке
743 ph->SetPrimary(primary) ;
744 ph->SetWeight(PrimaryWeight(primary)) ;
745 inPHOS++ ;
746 }
747
748 FillHistogram("hCenPHOS",fCentrality,inPHOS) ;
749 if(inPHOS==0){
750 PostData(1, fOutputContainer);
751 fEventCounter++;
752 return ;
753 }
754
755
756 for (Int_t i1=0; i1<inPHOS; i1++) {
757 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
758 Double_t w=ph1->GetWeight() ;
759 Double_t pt = ph1->Pt() ;
760 Double_t ptv =ph1->GetMomV2()->Pt() ;
761
762 FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ;
763 FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptv,w) ;
764 if(ph1->IsntUnfolded()){
765 FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ;
766 }
767 if(ph1->IsCPVOK()){
768 FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ;
769 FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptv,w) ;
770 }
771 if(ph1->IsCPV2OK()){
772 FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ;
773 }
774 if(ph1->IsDispOK()){
775 FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ;
776
777 FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptv,w) ;
778 if(ph1->IsntUnfolded()){
779 FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ;
780 }
781 if(ph1->IsCPVOK()){
782 FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ;
783 FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptv,w) ;
784 }
785 }
786 if(ph1->IsDisp2OK()){
787 FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ;
788 FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptv,w) ;
789 if(ph1->IsCPVOK()){
790 FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ;
791 FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptv,w) ;
792 }
793 }
794 }
795
796
797 const Double_t alphaCut=0.7 ;
798 //pi0
799 for (Int_t i1=0; i1<inPHOS-1; i1++) {
800 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
801 Double_t w1=ph1->GetWeight() ;
802 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
803 AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
804 Double_t w2=ph2->GetWeight() ;
805 Double_t w = TMath::Sqrt(w1*w2) ;
806 TLorentzVector p12 = *ph1 + *ph2;
807 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
808
809 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
810 Double_t m = p12.M() ;
811 Double_t pt = p12.Pt() ;
812 Double_t mv = p12.M() ;
813 Double_t ptv = p12.Pt() ;
814
815 FillHistogram(Form("hPi0All_cen%d",fCenBin),pt,m,w) ;
816 FillHistogram(Form("hPi0Allcore_cen%d",fCenBin),ptv,mv,w) ;
817 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
818 FillHistogram(Form("hPi0Allwou_cen%d",fCenBin),pt,m,w) ;
819 }
820
821
822 Double_t pt1 = ph1->Pt() ;
823 Double_t ptv1 = ph1->GetMomV2()->Pt() ;
824 Double_t pt2 = ph2->Pt() ;
825 Double_t ptv2 = ph2->GetMomV2()->Pt() ;
826 FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt1,w1) ;
827 FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt2,w2) ;
828 FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv1,w1) ;
829 FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv2,w2) ;
830 if(ph1->IsntUnfolded())
831 FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt1,w1) ;
832 if(ph2->IsntUnfolded())
833 FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt2,w2) ;
834 if(ph1->IsCPVOK()){
835 FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt1,w1) ;
836 FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv1,w1) ;
837 }
838 if(ph2->IsCPVOK()){
839 FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt2,w2) ;
840 FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv2,w2) ;
841 }
842 if(ph1->IsCPV2OK()){
843 FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt1,w1) ;
844 }
845 if(ph2->IsCPV2OK()){
846 FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt2,w2) ;
847 }
848 if(ph1->IsDispOK()){
849 FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt1,w1) ;
850 if(ph1->IsntUnfolded()){
851 FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt1,w1) ;
852 }
853 FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv1,w1) ;
854 }
855 if(ph2->IsDispOK()){
856 FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt2,w2) ;
857 if(ph1->IsntUnfolded()){
858 FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt2,w2) ;
859 }
860 FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv2,w2) ;
861 }
862 if(ph1->IsDisp2OK()){
863 FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt1,w1) ;
864 FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv1,w1) ;
865 }
866 if(ph2->IsDisp2OK()){
867 FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt2,w2) ;
868 FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv2,w2) ;
869 }
870 if(ph1->IsDispOK() && ph1->IsCPVOK()){
871 FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt1,w1) ;
872 FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv1,w1) ;
873 }
874 if(ph2->IsDispOK() && ph2->IsCPVOK()){
875 FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt2,w2) ;
876 FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv2,w2) ;
877 }
878 if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
879 FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt1,w1) ;
880 FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv1,w1) ;
881 }
882 if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
883 FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt2,w2) ;
884 FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv2,w2) ;
885 }
886
887
888 if(a<alphaCut){
889 FillHistogram(Form("hPi0All_a07_cen%d",fCenBin),pt,m,w) ;
890 }
891
892 if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 CPV
893 FillHistogram(Form("hPi0CPV_cen%d",fCenBin),pt,m,w) ;
894 FillHistogram(Form("hPi0CPVcore_cen%d",fCenBin),ptv,mv,w) ;
895
896 if(a<alphaCut){
897 FillHistogram(Form("hPi0CPV_a07_cen%d",fCenBin),pt,m,w) ;
898 }
899 }
900 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
901 FillHistogram(Form("hPi0CPV2_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
902 if(a<alphaCut){
903 FillHistogram(Form("hPi0CPV2_a07_cen%d",fCenBin),pt,m,w) ;
904 }
905 }
906 if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
907 FillHistogram(Form("hPi0Disp2_cen%d",fCenBin),pt,m,w) ;
908 FillHistogram(Form("hPi0Disp2core_cen%d",fCenBin),ptv,mv,w) ;
909 if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
910 FillHistogram(Form("hPi0Both2_cen%d",fCenBin),pt,m,w) ;
911 FillHistogram(Form("hPi0Both2core_cen%d",fCenBin),ptv,mv,w) ;
912 if(a<alphaCut){
913 FillHistogram(Form("hPi0Both2_a07_cen%d",fCenBin),pt,m,w) ;
914 }
915 }
916 }
917
918
919 if(ph1->IsDispOK() && ph2->IsDispOK()){ //pi0 Disp
920 FillHistogram(Form("hPi0Disp_cen%d",fCenBin),pt,m,w) ;
921 FillHistogram(Form("hPi0Dispcore_cen%d",fCenBin),ptv,mv,w) ;
922 if(a<alphaCut){
923 FillHistogram(Form("hPi0Disp_a07_cen%d",fCenBin),pt,m,w) ;
924 }
925
926 if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
927 FillHistogram(Form("hPi0Both_cen%d",fCenBin),pt,m,w) ;
928 FillHistogram(Form("hPi0Bothcore_cen%d",fCenBin),ptv,mv,w) ;
929
930 if(a<alphaCut){
931 FillHistogram(Form("hPi0Both_a07_cen%d",fCenBin),pt,m,w) ;
932 }
933 }
934 }
935 } // end of loop i2
936 } // end of loop i1
937
938 //now mixed
939 for (Int_t i1=0; i1<inPHOS; i1++) {
940 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
941 for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
942 TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
943 for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
944 AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
945 TLorentzVector p12 = *ph1 + *ph2;
946 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
947
948 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
949 Double_t m = p12.M() ;
950 Double_t pt = p12.Pt() ;
951 Double_t mv = p12.M() ;
952 Double_t ptv = p12.Pt() ;
953
954 FillHistogram(Form("hMiPi0All_cen%d",fCenBin),pt,m) ;
955 FillHistogram(Form("hMiPi0Allcore_cen%d",fCenBin),ptv,mv) ;
956 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
957 FillHistogram(Form("hMiPi0Allwou_cen%d",fCenBin),pt,m) ;
958 }
959
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) ;
965 FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt2) ;
966 FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv1) ;
967 FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv2) ;
968 if(ph1->IsntUnfolded())
969 FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt1) ;
970 if(ph2->IsntUnfolded())
971 FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt2) ;
972 if(ph1->IsCPVOK()){
973 FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt1) ;
974 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv1) ;
975 }
976 if(ph2->IsCPVOK()){
977 FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt2) ;
978 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv2) ;
979 }
980 if(ph1->IsCPV2OK()){
981 FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt1) ;
982 }
983 if(ph2->IsCPV2OK()){
984 FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt2) ;
985 }
986 if(ph1->IsDispOK()){
987 FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt1) ;
988 if(ph1->IsntUnfolded()){
989 FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt1) ;
990 }
991 FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv1) ;
992 }
993 if(ph2->IsDispOK()){
994 FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt2) ;
995 if(ph1->IsntUnfolded()){
996 FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt2) ;
997 }
998 FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv2) ;
999 }
1000 if(ph1->IsDisp2OK()){
1001 FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt1) ;
1002 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv1) ;
1003 }
1004 if(ph2->IsDisp2OK()){
1005 FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt2) ;
1006 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv2) ;
1007 }
1008 if(ph1->IsDispOK() && ph1->IsCPVOK()){
1009 FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt1) ;
1010 FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv1) ;
1011 }
1012 if(ph2->IsDispOK() && ph2->IsCPVOK()){
1013 FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt2) ;
1014 FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv2) ;
1015 }
1016 if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
1017 FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt1) ;
1018 FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv1) ;
1019 }
1020 if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
1021 FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt2) ;
1022 FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv2) ;
1023 }
1024
1025
1026
1027 if(a<alphaCut){
1028 FillHistogram(Form("hMiPi0All_a07_cen%d",fCenBin),pt,m) ;
1029 }
1030 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1031 FillHistogram(Form("hMiPi0CPV_cen%d",fCenBin),p12.M(), p12.Pt()) ;
1032 FillHistogram(Form("hMiPi0CPVcore_cen%d",fCenBin),ptv,mv) ;
1033
1034 if(a<alphaCut){
1035 FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCenBin),pt,m) ;
1036 }
1037 }
1038 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1039 FillHistogram(Form("hMiPi0CPV2_cen%d",fCenBin),pt,m) ;
1040
1041 if(a<alphaCut){
1042 FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCenBin),pt,m) ;
1043 }
1044 }
1045 if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1046 FillHistogram(Form("hMiPi0Disp2_cen%d",fCenBin),pt,m) ;
1047 FillHistogram(Form("hMiPi0Disp2core_cen%d",fCenBin),ptv,mv) ;
1048 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1049 FillHistogram(Form("hMiPi0Both2_cen%d",fCenBin),pt,m) ;
1050 FillHistogram(Form("hMiPi0Both2core_cen%d",fCenBin),pv12.M(),pv12.Pt()) ;
1051 }
1052 }
1053 if(ph1->IsDispOK() && ph2->IsDispOK()){
1054
1055 FillHistogram(Form("hMiPi0Disp_cen%d",fCenBin),pt,m) ;
1056 FillHistogram(Form("hMiPi0Dispcore_cen%d",fCenBin),pv12.M(),pv12.Pt()) ;
1057 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1058 FillHistogram(Form("hMiPi0Dispwou_cen%d",fCenBin),p12.M(),p12.Pt()) ;
1059 }
1060
1061 if(a<alphaCut){
1062 FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCenBin),pt,m) ;
1063 }
1064 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1065 FillHistogram(Form("hMiPi0Both_cen%d",fCenBin),pt,m) ;
1066 FillHistogram(Form("hMiPi0Bothcore_cen%d",fCenBin),pv12.M(),pv12.Pt()) ;
1067
1068 if(a<alphaCut){
1069 FillHistogram(Form("hMiPi0Both_a07_cen%d",fCenBin),pt,m) ;
1070 }
1071 }
1072 }
1073 } // end of loop i2
1074 }
1075 } // end of loop i1
1076
1077
1078 FillSecondaries() ;
1079
1080
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) ;
1086 fPHOSEvent=0;
1087 if(prevPHOS->GetSize()>kMixEvents[fCenBin]){//Remove redundant events
1088 TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1089 prevPHOS->RemoveLast() ;
1090 delete tmp ;
1091 }
1092 }
1093 // Post output data.
1094 PostData(1, fOutputContainer);
1095 fEventCounter++;
1096}
1097//________________________________________________________________________
1098void AliPHOSHijingEfficiency::FillSecondaries(){
1099 //Sort secondaires
1100
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) ;
1106 if(p->R()>rcut)
1107 continue ;
1108 if(TMath::Abs(p->Y())>0.5)
1109 continue ;
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);
1114 break ;
1115 case 11:
1116 case -11:
1117 FillHistogram(Form("hPrimEl_cen%d",fCenBin),p->Pt(),w);
1118 break ;
1119 case 111:
1120 FillHistogram(Form("hPrimPi0_cen%d",fCenBin),p->Pt(),w);
1121 break ;
1122 case 221:
1123 FillHistogram(Form("hPrimEta_cen%d",fCenBin),p->Pt(),w);
1124 break ;
1125 case 211:
1126 case -211:
1127 FillHistogram(Form("hPrimPipm_cen%d",fCenBin),p->Pt(),w);
1128 break ;
1129 case 2212: //p
1130 FillHistogram(Form("hPrimP_cen%d",fCenBin),p->Pt(),w);
1131 break ;
1132 case -2212: //pbar
1133 FillHistogram(Form("hPrimPbar_cen%d",fCenBin),p->Pt(),w);
1134 break ;
1135 case 2112: //n
1136 FillHistogram(Form("hPrimN_cen%d",fCenBin),p->Pt(),w);
1137 break ;
1138 case -2112: //nbar
1139 FillHistogram(Form("hPrimNbar_cen%d",fCenBin),p->Pt(),w);
1140 break ;
1141 case 310: //nbar
1142 FillHistogram(Form("hPrimK0S_cen%d",fCenBin),p->Pt(),w);
1143 break ;
1144 case 130: //nbar
1145 FillHistogram(Form("hPrimK0L_cen%d",fCenBin),p->Pt(),w);
1146 break ;
1147 case 321: //K+
1148 case -321: //K-
1149 FillHistogram(Form("hPrimKpm_cen%d",fCenBin),p->Pt(),w);
1150 break ;
1151 default: //other
1152 FillHistogram(Form("hPrimOther_cen%d",fCenBin),p->Pt(),w);
1153 }
1154 }
1155
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)
1160 continue ;
1161 FillHistogram("Vertex",p->Pt(),p->R());
1162 if(p->R()<rcut)
1163 continue ;
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) ;
1170 }
1171
1172 TLorentzVector p1;
1173
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()) ;
1185 if(prim>-1){
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())) ;
1188
1189
1190 Int_t pdgCode=particle->GetPdgCode() ;
1191 if(pdgCode!=111){ //common parent - not pi0
1192 if(pdgCode==22)
1193 FillHistogram(Form("hParentGamma_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1194 else{
1195 if(pdgCode==11 || pdgCode==-11){
1196 FillHistogram(Form("hParentEl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1197 }
1198 else{
1199 if(InPi0mass(p12.M() ,p12.Pt())){
1200 printf("Common parent: %d \n",pdgCode) ;
1201 }
1202 FillHistogram(Form("hParentOther_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1203 }
1204 }//Not photons
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) ;
1212 continue ;
1213 }
1214 //Common particle pi0, created off-vertex
1215 Int_t primPi0=particle->GetFirstMother();
1216 if(primPi0==-1){
1217 FillHistogram(Form("hParentPi0NoPrim_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1218 }
1219 else{
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
1223 break ;
1224 case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //omega
1225 break ;
1226 case 211: //pi+-
1227 case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1228 break ;
1229 case 321: //K+-
1230 case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1231 break ;
1232 case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0S
1233 break ;
1234 case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0L
1235 break ;
1236 case 2212: //p
1237 case 2112: //n
1238 FillHistogram(Form("hParentPi0pn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn
1239 break ;
1240 case -2212: //pbar
1241 case -2112: //nbar
1242 FillHistogram(Form("hParentPi0antipn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn
1243 break ;
1244 default: //other
1245 FillHistogram(Form("hParentPi0Other_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1246 }//switch
1247 }//pi0 with primary
1248 }//common parent - pi0
1249 }//there is common primary
1250 }//seond photon loop
1251 }//first photon loop
1252
1253
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() ;
1258 if(iprim<0)
1259 FillAllHistograms(Form("hGammaNoPrim_cen%d",fCenBin),ph1) ; //
1260 else{
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();
1268 }
1269
1270 //photon
1271 Int_t primPdgCode=primVtx->GetPdgCode() ;
1272 switch(primPdgCode){
1273 case 22: FillAllHistograms("hGammaPhot",ph1);
1274 break ;
1275 case 11:
1276 case -11:
1277 FillAllHistograms("hGammaEl",ph1);
1278 break ;
1279 case 111:
1280 FillAllHistograms("hGammaPi0",ph1);
1281 break ;
1282 case 221:
1283 FillAllHistograms("hGammaEta",ph1);
1284 break ;
1285 case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
1286 break ;
1287 case 211:
1288 case -211:
1289 FillAllHistograms("hGammaPipm",ph1);
1290 //Find particle entered PHOS
1291 if(primVtx == primPHOS)
1292 FillAllHistograms("hGammaPipmDirect",ph1);
1293 else{
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());
1299 break ;
1300 }
1301 if(TMath::Abs(primPdgPHOS)==11){
1302 FillAllHistograms("hGammaPipmEl",ph1);
1303 FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R());
1304 break ;
1305 }
1306 if(TMath::Abs(primPdgPHOS)==2212){
1307 FillAllHistograms("hGammaPipmp",ph1);
1308 FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1309 break ;
1310 }
1311 if(TMath::Abs(primPdgPHOS)==2112){
1312 FillAllHistograms("hGammaPipmn",ph1);
1313 FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1314 break ;
1315 }
1316 FillAllHistograms("hGammaPipmOther",ph1);
1317 FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R());
1318 }
1319 break ;
1320 case 2212: //p
1321 FillAllHistograms("hGammaP",ph1);
1322 break ;
1323 case -2212: //pbar
1324 FillAllHistograms("hGammaPbar",ph1);
1325 break ;
1326 case 2112: //n
1327 FillAllHistograms("hGammaN",ph1);
1328 break ;
1329 case -2112: //nbar
1330 FillAllHistograms("hGammaNbar",ph1) ; // pn
1331 break ;
1332 case 310: //nbar
1333 FillAllHistograms("hGammaK0S",ph1) ; // pn
1334 break ;
1335 case 130: //nbar
1336 FillAllHistograms("hGammaK0L",ph1) ; // pn
1337 break ;
1338 case 321: //K+
1339 case -321: //K-
1340 FillAllHistograms("hGammaKpm",ph1) ; // pn
1341 break ;
1342 case -323:
1343 case 323:
1344 case -313:
1345 case 313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
1346 break ;
1347
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
1357 break ;
1358 default: //other
1359 if(primVtx->GetPDG()->Charge())
1360 FillAllHistograms("hGammaOtherCharged",ph1) ; //
1361 else
1362 FillAllHistograms("hGammaOtherNeutral",ph1) ; //
1363 }
1364 }
1365
1366 }//single photons
1367
1368}
1369//________________________________________________________________________
1370void AliPHOSHijingEfficiency::Terminate(Option_t *)
1371{
1372 // Draw result to the screen
1373 // Called once at the end of the query
1374
1375}
1376
1377//________________________________________________________________________
1378Bool_t AliPHOSHijingEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
1379{
1380 //Check if this channel belogs to the good ones
1381
1382 if(strcmp(det,"PHOS")==0){
1383 if(mod>5 || mod<1){
1384 AliError(Form("No bad map for PHOS module %d ",mod)) ;
1385 return kTRUE ;
1386 }
1387 if(!fPHOSBadMap[mod]){
1388 AliError(Form("No Bad map for PHOS module %d",mod)) ;
1389 return kTRUE ;
1390 }
1391 if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
1392 return kFALSE ;
1393 else
1394 return kTRUE ;
1395 }
1396 else{
1397 AliError(Form("Can not find bad channels for detector %s ",det)) ;
1398 }
1399 return kTRUE ;
1400}
1401//_____________________________________________________________________________
1402void AliPHOSHijingEfficiency::FillAllHistograms(const char * key,AliCaloPhoton * ph)const{
1403 //Fill All PID histograms
1404
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) ;
1410 if(ph->IsCPVOK()){
1411 FillHistogram(Form("%s_CPV_cen%d",key,fCenBin),pt,w) ;
1412 FillHistogram(Form("%s_CPVcore_cen%d",key,fCenBin),ptC,w) ;
1413 }
1414 if(ph->IsCPV2OK()){
1415 FillHistogram(Form("%s_CPV2_cen%d",key,fCenBin),pt,w) ;
1416 FillHistogram(Form("%s_CPV2core_cen%d",key,fCenBin),ptC,w) ;
1417 }
1418 if(ph->IsDispOK()){
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) ;
1424 }
1425 if(ph->IsCPVOK()){
1426 FillHistogram(Form("%s_Both_cen%d",key,fCenBin),pt,w) ;
1427 FillHistogram(Form("%s_Bothcore_cen%d",key,fCenBin),ptC,w) ;
1428 }
1429 }
1430}
1431//_____________________________________________________________________________
1432void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x)const{
1433 //FillHistogram
1434 TObject * tmp = fOutputContainer->FindObject(key) ;
1435 if(!tmp){
1436 AliInfo(Form("can not find histogram <%s> ",key)) ;
1437 return ;
1438 }
1439 if(tmp->IsA() == TClass::GetClass("TH1I")){
1440 ((TH1I*)tmp)->Fill(x) ;
1441 return ;
1442 }
1443 if(tmp->IsA() == TClass::GetClass("TH1F")){
1444 ((TH1F*)tmp)->Fill(x) ;
1445 return ;
1446 }
1447 if(tmp->IsA() == TClass::GetClass("TH1D")){
1448 ((TH1D*)tmp)->Fill(x) ;
1449 return ;
1450 }
1451 AliInfo(Form("can not find 1D histogram <%s> ",key)) ;
1452}
1453//_____________________________________________________________________________
1454void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
1455 //FillHistogram
1456 TObject * tmp = fOutputContainer->FindObject(key) ;
1457 if(!tmp){
1458 AliInfo(Form("can not find histogram <%s> ",key)) ;
1459 return ;
1460 }
1461 if(tmp->IsA() == TClass::GetClass("TH1F")){
1462 ((TH1F*)tmp)->Fill(x,y) ;
1463 return ;
1464 }
1465 if(tmp->IsA() == TClass::GetClass("TH2F")){
1466 ((TH2F*)tmp)->Fill(x,y) ;
1467 return ;
1468 }
1469 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1470}
1471
1472//_____________________________________________________________________________
1473void 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) ;
1476 if(!tmp){
1477 AliInfo(Form("can not find histogram <%s> ",key)) ;
1478 return ;
1479 }
1480 if(tmp->IsA() == TClass::GetClass("TH2F")){
1481 ((TH2F*)tmp)->Fill(x,y,z) ;
1482 return ;
1483 }
1484 if(tmp->IsA() == TClass::GetClass("TH3F")){
1485 ((TH3F*)tmp)->Fill(x,y,z) ;
1486 return ;
1487 }
1488}
1489//_____________________________________________________________________________
1490void 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) ;
1493 if(!tmp){
1494 AliInfo(Form("can not find histogram <%s> ",key)) ;
1495 return ;
1496 }
1497 if(tmp->IsA() == TClass::GetClass("TH3F")){
1498 ((TH3F*)tmp)->Fill(x,y,z,w) ;
1499 return ;
1500 }
1501}
1502
1503//___________________________________________________________________________
1504Int_t AliPHOSHijingEfficiency::ConvertRunNumber(Int_t run){
1505
1506 switch(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;
1645 }
1646
1647}
1648//_____________________________________________________________________________
1649Bool_t AliPHOSHijingEfficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1650
1651 //For R=3.5
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 ;
1657/*
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) ;
1663*/
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) ;
1668
1669}
1670//_____________________________________________________________________________
1671Bool_t AliPHOSHijingEfficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1672
1673//For R=4.5
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) ;
1679
1680/*
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) ;
1686*/
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) ;
1691
1692}
1693//____________________________________________________________________________
1694Double_t AliPHOSHijingEfficiency::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1695 //Parameterization of LHC10h period
1696 //_true if neutral_
1697
1698 Double_t meanX=0;
1699 Double_t meanZ=0.;
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 Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1705
1706 if(mf<0.){ //field --
1707 meanZ = -0.468318 ;
1708 if(charge>0)
1709 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)) ;
1710 else
1711 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)) ;
1712 }
1713 else{ //Field ++
1714 meanZ= -0.468318;
1715 if(charge>0)
1716 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)) ;
1717 else
1718 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)) ;
1719 }
1720
1721 Double_t rz=(dz-meanZ)/sz ;
1722 Double_t rx=(dx-meanX)/sx ;
1723 return TMath::Sqrt(rx*rx+rz*rz) ;
1724}
1725//____________________________________________________________________________
1726Bool_t AliPHOSHijingEfficiency::TestTOF(Double_t t, Double_t e){
1727
1728 Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09
1729 +2.24611e-09*2.24611e-09/e
1730 +5.65054e-09*5.65054e-09/e/e) ;
1731 sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part
1732 Double_t mean=1.1e-9 ;
1733 if(TMath::Abs(t-mean)<2.*sigma)
1734 return kTRUE ;
1735 else
1736 if(TMath::Abs(t-mean+100.e-9)<2.*sigma)
1737 return kTRUE ;
1738
1739 return kFALSE ;
1740}
1741//____________________________________________________________________________
1742Double_t AliPHOSHijingEfficiency::CoreEnergy(AliESDCaloCluster * clu){
1743 //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1744
1745 //Can not use already calculated coordinates?
1746 //They have incidence correction...
1747 const Double_t distanceCut =3.5 ;
1748 const Double_t logWeight=4.5 ;
1749
1750 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1751// Calculates the center of gravity in the local PHOS-module coordinates
1752 Float_t wtot = 0;
1753 Double_t xc[100]={0} ;
1754 Double_t zc[100]={0} ;
1755 Double_t x = 0 ;
1756 Double_t z = 0 ;
1757 Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1758 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1759 Int_t relid[4] ;
1760 Float_t xi ;
1761 Float_t zi ;
1762 fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1763 fPHOSGeo->RelPosInModule(relid, xi, zi);
1764 xc[iDigit]=xi ;
1765 zc[iDigit]=zi ;
1766 if (clu->E()>0 && elist[iDigit]>0) {
1767 Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1768 x += xc[iDigit] * w ;
1769 z += zc[iDigit] * w ;
1770 wtot += w ;
1771 }
1772 }
1773 if (wtot>0) {
1774 x /= wtot ;
1775 z /= wtot ;
1776 }
1777 Double_t coreE=0. ;
1778 for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1779 Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1780 if(distance < distanceCut)
1781 coreE += elist[iDigit] ;
1782 }
1783 //Apply non-linearity correction
1784 return rnlin(coreE) ;
1785}
1786//____________________________________________________________________________
1787Bool_t AliPHOSHijingEfficiency::AreNeibors(Int_t id1,Int_t id2){
1788
1789 Int_t relid1[4] ;
1790 fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
1791
1792 Int_t relid2[4] ;
1793 fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
1794
1795 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
1796 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
1797 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
1798
1799 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
1800 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
1801 return 1 ;
1802 }
1803 else {
1804 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
1805 return 0; // Difference in row numbers is too large to look further
1806 }
1807 return 0 ;
1808
1809 }
1810
1811 return 0 ;
1812}
1813//____________________________________________________________________________
1814void AliPHOSHijingEfficiency::Reclusterize(AliESDCaloCluster * clu){
1815 //Re-clusterize to make continues cluster
1816
1817 const Int_t oldMulDigit=clu->GetNCells() ;
1818 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1819 UShort_t * dlist = clu->GetCellsAbsId();
1820
1821 Int_t index[oldMulDigit] ;
1822 Bool_t used[oldMulDigit] ;
1823 for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
1824 Int_t inClu=0 ;
1825 Double_t eMax=0. ;
1826 //find maximum
1827 for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
1828 if(eMax<elist[iDigit]){
1829 eMax=elist[iDigit];
1830 index[0]=iDigit ;
1831 inClu=1 ;
1832 }
1833 }
1834 if(inClu==0){ //empty cluster
1835 return ;
1836 }
1837 used[index[0]]=kTRUE ; //mark as used
1838 for(Int_t i=0; i<inClu; i++){
1839 for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
1840 if(used[iDigit]) //already used
1841 continue ;
1842 if(AreNeibors(dlist[index[i]],dlist[iDigit])){
1843 index[inClu]= iDigit ;
1844 inClu++ ;
1845 used[iDigit]=kTRUE ;
1846 }
1847 }
1848 }
1849
1850 if(inClu==oldMulDigit) //no need to modify
1851 return ;
1852
1853 clu->SetNCells(inClu);
1854 //copy
1855 UShort_t tmpD[oldMulDigit] ;
1856 Double_t tmpE[oldMulDigit] ;
1857 for(Int_t i=0; i<oldMulDigit; i++){
1858 tmpD[i]=dlist[i] ;
1859 tmpE[i]=elist[i] ;
1860 }
1861 //change order of digits in list so that
1862 //first inClu cells were true ones
1863 for(Int_t i=0; i<inClu; i++){
1864 dlist[i]=tmpD[index[i]] ;
1865 elist[i]=tmpE[index[i]] ;
1866 }
1867
1868
1869}
1870
1871//____________________________________________________________________________
1872void AliPHOSHijingEfficiency::EvalLambdas(AliESDCaloCluster * clu, Int_t iR,Double_t &m02, Double_t &m20){
1873 //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
1874
1875 Double_t rCut=0. ;
1876 if(iR==0)
1877 rCut=3.5 ;
1878 else
1879 rCut=4.5 ;
1880
1881
1882 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1883// Calculates the center of gravity in the local PHOS-module coordinates
1884 Float_t wtot = 0;
1885 Double_t xc[100]={0} ;
1886 Double_t zc[100]={0} ;
1887 Double_t x = 0 ;
1888 Double_t z = 0 ;
1889 Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1890 const Double_t logWeight=4.5 ;
1891 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1892 Int_t relid[4] ;
1893 Float_t xi ;
1894 Float_t zi ;
1895 fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1896 fPHOSGeo->RelPosInModule(relid, xi, zi);
1897 xc[iDigit]=xi ;
1898 zc[iDigit]=zi ;
1899 if (clu->E()>0 && elist[iDigit]>0) {
1900 Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1901 x += xc[iDigit] * w ;
1902 z += zc[iDigit] * w ;
1903 wtot += w ;
1904 }
1905 }
1906 if (wtot>0) {
1907 x /= wtot ;
1908 z /= wtot ;
1909 }
1910
1911 wtot = 0. ;
1912 Double_t dxx = 0.;
1913 Double_t dzz = 0.;
1914 Double_t dxz = 0.;
1915 Double_t xCut = 0. ;
1916 Double_t zCut = 0. ;
1917 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1918 if (clu->E()>0 && elist[iDigit]>0.) {
1919 Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1920 Double_t xi= xc[iDigit] ;
1921 Double_t zi= zc[iDigit] ;
1922 if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
1923 xCut += w * xi ;
1924 zCut += w * zi ;
1925 dxx += w * xi * xi ;
1926 dzz += w * zi * zi ;
1927 dxz += w * xi * zi ;
1928 wtot += w ;
1929 }
1930 }
1931
1932 }
1933 if (wtot>0) {
1934 xCut/= wtot ;
1935 zCut/= wtot ;
1936 dxx /= wtot ;
1937 dzz /= wtot ;
1938 dxz /= wtot ;
1939 dxx -= xCut * xCut ;
1940 dzz -= zCut * zCut ;
1941 dxz -= xCut * zCut ;
1942
1943 m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
1944 m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
1945 }
1946 else {
1947 m20=m02=0.;
1948 }
1949
1950}
1951//___________________________________________________________________________
1952void AliPHOSHijingEfficiency::ProcessMC(){
1953 //fill histograms for efficiensy etc. calculation
1954 const Double_t rcut = 1. ; //cut for primary particles
1955 //---------First pi0/eta-----------------------------
1956 char partName[10] ;
1957 char hkey[55] ;
1958
1959 if(!fStack) return ;
1960 for(Int_t i=0;i<fStack->GetNtrack();i++){
1961 TParticle* particle = fStack->Particle(i);
1962 if(particle->GetPdgCode() == 111)
1963 snprintf(partName,10,"pi0") ;
1964 else
1965 if(particle->GetPdgCode() == 221)
1966 snprintf(partName,10,"eta") ;
1967 else
1968 if(particle->GetPdgCode() == 22)
1969 snprintf(partName,10,"gamma") ;
1970 else
1971 continue ;
1972
1973 //Primary particle
1974 Double_t r=particle->R() ;
1975 Double_t pt = particle->Pt() ;
1976 //Distribution over vertex
1977 FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
1978
1979 if(r >rcut)
1980 continue ;
1981
1982 //Total number of pi0 with creation radius <1 cm
1983 Double_t w = PrimaryParticleWeight(particle) ;
1984 snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
1985 FillHistogram(hkey,pt,w) ;
1986 if(TMath::Abs(particle->Y())<0.12){
1987 snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
1988 FillHistogram(hkey,pt,w) ;
1989 }
1990
1991 snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
1992 FillHistogram(hkey,particle->Y(),w) ;
1993
1994 Double_t phi=particle->Phi() ;
1995 while(phi<0.)phi+=TMath::TwoPi() ;
1996 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1997 snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
1998 FillHistogram(hkey,phi,w) ;
1999 }
2000}
2001//___________________________________________________________________________
2002Int_t AliPHOSHijingEfficiency::FindPrimary(AliESDCaloCluster*clu, Bool_t&sure){
2003 //Finds primary and estimates if it unique one?
3c77bbdf 2004 //First check can it be photon/electron
2005 const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
04c118da 2006 Int_t n=clu->GetNLabels() ;
3c77bbdf 2007 for(Int_t i=0; i<n; i++){
2008 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
2009 Int_t pdg = p->GetPdgCode() ;
2010 if(pdg==22 || pdg==11 || pdg == -11){
2011 if(p->Energy()>emFraction*clu->E()){
2012 sure=kTRUE ;
2013 return clu->GetLabelAt(i);
2014 }
2015 }
2016 }
2017
04c118da 2018 Double_t* Ekin= new Double_t[n] ;
2019 for(Int_t i=0; i<n; i++){
2020 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
2021 Ekin[i]=p->P() ; // estimate of kinetic energy
2022 if(p->GetPdgCode()==-2212 || p->GetPdgCode()==-2112){
2023 Ekin[i]+=1.8 ; //due to annihilation
2024 }
2025 }
2026 Int_t iMax=0;
2027 Double_t eMax=0.,eSubMax=0. ;
2028 for(Int_t i=0; i<n; i++){
2029 if(Ekin[i]>eMax){
2030 eSubMax=eMax;
2031 eMax=Ekin[i];
2032 iMax=i;
2033 }
2034 }
2035 if(eSubMax>0.8*eMax)//not obvious primary
2036 sure=kFALSE;
2037 else
2038 sure=kTRUE;
2039 delete Ekin;
2040 return clu->GetLabelAt(iMax) ;
2041}
2042//___________________________________________________________________________
2043Double_t AliPHOSHijingEfficiency::PrimaryWeight(Int_t primary){
2044 //Check who is the primary and introduce weight to correct primary spectrum
2045
2046 if(primary<0 || primary>=fStack->GetNtrack())
2047 return 1 ;
2048 //trace primaries up to IP
2049 TParticle* particle = fStack->Particle(primary);
2050 Double_t r=particle->R() ;
2051 Int_t mother = particle->GetFirstMother() ;
2052 while(mother>-1){
2053 if(r<1. && particle->GetPdgCode()==111)
2054 break ;
2055 particle = fStack->Particle(mother);
2056 mother = particle->GetFirstMother() ;
2057 r=particle->R() ;
2058 }
2059
2060 return TMath::Max(0.,PrimaryParticleWeight(particle)) ;
2061}
2062//________________________________________________________________________
2063Double_t AliPHOSHijingEfficiency::PrimaryParticleWeight(TParticle * particle){
2064
2065 Int_t pdg = particle->GetPdgCode() ;
2066 Int_t type=0 ;
2067 if(pdg == 111 || TMath::Abs(pdg)==211){
2068 type =1 ;
2069 }
2070 else{
2071 if(TMath::Abs(pdg)<1000){ //Kaon-like
2072 type =2 ;
2073 }
2074 else
2075 type = 3; //baryons
2076 }
2077
2078 Double_t x = particle->Pt() ;
2079 if(type==1){
2080 if(fCenBin==0) //0-5
2081 return (1.662990+1.140890*x-0.192088*x*x)/(1.-0.806630*x+0.304771*x*x)+0.141690*x ;
2082 if(fCenBin==1) //5-10
2083 return (1.474351+0.791492*x-0.066369*x*x)/(1.-0.839338*x+0.317312*x*x)+0.093289*x ;
2084 if(fCenBin==2) //10-20
2085 return (1.174728+0.959681*x-0.137695*x*x)/(1.-0.788873*x+0.299538*x*x)+0.128759*x ;
2086 if(fCenBin==3) //20-40
2087 return (0.927335+0.475349*x+0.004364*x*x)/(1.-0.817966*x+0.309787*x*x)+0.086899*x ;
2088 if(fCenBin==4) //40-60
2089 return (0.676878+0.190680*x+0.077031*x*x)/(1.-0.790623*x+0.305183*x*x)+0.064510*x ;
2090 if(fCenBin==5) //60-80
2091 return (0.684726-0.606262*x+0.409963*x*x)/(1.-1.080061*x+0.456933*x*x)+0.005151*x ;
2092 }
2093 if(type==2){
2094 if(fCenBin==0) //0-5
2095 return (-0.417131+2.253936*x-0.337731*x*x)/(1.-0.909892*x+0.316820*x*x)+0.157312*x ;
2096 if(fCenBin==1) //5-10
2097 return (-0.352275+1.844466*x-0.248598*x*x)/(1.-0.897048*x+0.316462*x*x)+0.132461*x ;
2098 if(fCenBin==2) //10-20
2099 return (-0.475481+1.975108*x-0.336013*x*x)/(1.-0.801028*x+0.276705*x*x)+0.188164*x ;
2100 if(fCenBin==3) //20-40
2101 return (-0.198954+1.068789*x-0.103540*x*x)/(1.-0.848354*x+0.299209*x*x)+0.112939*x ;
2102 if(fCenBin==4) //40-60
2103 return (-0.111052+0.664041*x-0.019717*x*x)/(1.-0.804916*x+0.300779*x*x)+0.095784*x ;
2104 if(fCenBin==5) //0-5
2105 return (0.202788-0.439832*x+0.564585*x*x)/(1.-1.254029*x+0.679444*x*x)+0.016235*x ;
2106 }
2107 if(type==3){
2108 if(fCenBin==0) //0-5
2109 return (-1.312732+2.743568*x-0.375775*x*x)/(1.-0.717533*x+0.164694*x*x)+0.164445*x ;
2110 if(fCenBin==1) //5-10
2111 return (-1.229425+2.585889*x-0.330164*x*x)/(1.-0.715892*x+0.167386*x*x)+0.133085*x ;
2112 if(fCenBin==2) //10-20
2113 return (-1.135677+2.397489*x-0.320355*x*x)/(1.-0.709312*x+0.164350*x*x)+0.146095*x ;
2114 if(fCenBin==3) //20-40
2115 return (-0.889993+1.928263*x-0.220785*x*x)/(1.-0.715991*x+0.174729*x*x)+0.095098*x ;
2116 if(fCenBin==4) //40-60
2117 return (-0.539237+1.329118*x-0.115439*x*x)/(1.-0.722906*x+0.186832*x*x)+0.059267*x ;
2118 if(fCenBin==5) //60-80
2119 return (-0.518126+1.327628*x-0.130881*x*x)/(1.-0.665649*x+0.184300*x*x)+0.081701*x ;
2120 }
2121 return 1. ;
2122}
2123//________________________________________________________________________
2124Int_t AliPHOSHijingEfficiency::FindCommonParent(Int_t iPart, Int_t jPart){
2125 //check if there is a common parent for particles i and j
2126 // -1: no common parent or wrong iPart/jPart
2127
2128 if(iPart==-1 || iPart>=fStack->GetNtrack() ||
2129 jPart==-1 || jPart>=fStack->GetNtrack()) return -1;
2130
2131 Int_t iprim1=iPart;
2132 while(iprim1>-1){
2133 Int_t iprim2=jPart;
2134 while(iprim2>-1){
2135 if(iprim1==iprim2)
2136 return iprim1 ;
2137 iprim2=((TParticle *)fStack->Particle(iprim2))->GetFirstMother();
2138 }
2139 iprim1=((TParticle *)fStack->Particle(iprim1))->GetFirstMother();
2140 }
2141 return -1;
2142}
2143//________________________________________________________________________
2144Bool_t AliPHOSHijingEfficiency::HaveParent(Int_t iPart, Int_t pdgParent){
2145 //check if there is a common parent for particles i and j
2146 // -1: no common parent or wrong iPart/jPart
2147
2148 if(iPart==-1 || iPart>=fStack->GetNtrack()) return -1;
2149
2150 Int_t iprim1=iPart;
2151 while(iprim1>-1){
2152 TParticle * tmp = fStack->Particle(iprim1) ;
2153 if(tmp->GetPdgCode()==pdgParent)
2154 return kTRUE ;
2155 iprim1=tmp->GetFirstMother();
2156 }
2157 return kFALSE;
2158}
2159//________________________________________________________________________
2160Bool_t AliPHOSHijingEfficiency::InPi0mass(Double_t m, Double_t /*pt*/){
2161
2162 return TMath::Abs(m-0.135)<0.007*2.5 ;
2163}