]>
Commit | Line | Data |
---|---|---|
58c00d93 | 1 | // $Id$ |
2 | ||
3 | //_________________________________________________________________________ | |
4 | // Class for easier handling of track-cluster electron PID. | |
5 | // | |
6 | // Author: Tomas Aronsson (Yale) | |
7 | /////////////////////////////////////////////////////////////////////////// | |
8 | ||
9 | #include <Riostream.h> | |
10 | #include <TLorentzVector.h> | |
11 | #include <TObjArray.h> | |
12 | #include <AliEMCALGeometry.h> | |
13 | #include <AliESD.h> | |
14 | #include <AliESD.h> | |
15 | #include <AliESDCaloCells.h> | |
16 | #include <AliESDCaloCluster.h> | |
17 | #include <AliESDtrack.h> | |
18 | #include <AliESDtrack.h> | |
19 | #include "AliEMCALClusterParams.h" | |
c2aad3ae | 20 | using std::endl; |
21 | using std::cout; | |
22 | using std::cerr; | |
58c00d93 | 23 | |
24 | ClassImp(AliEMCALClusterParams) | |
25 | ||
26 | //____________________________________________________________________________ | |
27 | AliEMCALClusterParams::AliEMCALClusterParams(AliESDtrack *trackin, | |
28 | AliESDCaloCluster *clusin, | |
29 | AliEMCALGeometry *geometryin, | |
30 | AliESDCaloCells *cellsin) : | |
31 | TObject(), | |
32 | fTrack(trackin), | |
33 | fCluster(clusin), | |
34 | fGeom(geometryin), | |
35 | fCells(cellsin) | |
36 | { | |
37 | // Constructor. | |
38 | } | |
39 | ||
40 | //____________________________________________________________________________ | |
41 | Double_t AliEMCALClusterParams::GetPe() const | |
42 | { | |
43 | Double_t pe=fTrack->Pt()/fCluster->E(); | |
44 | return pe; | |
45 | } | |
46 | ||
47 | //____________________________________________________________________________ | |
48 | Int_t AliEMCALClusterParams::IsElectron() const | |
49 | { | |
50 | Double_t ep=fCluster->E()/fTrack->Pt(); | |
51 | // if (fCluster->GetNCells() < 2 ) return 0; | |
52 | // if (fCluster->GetNCells() > 35 ) return 0; | |
53 | // if (fCluster->E() < 0 ) return 0; | |
54 | // if (fCluster->GetDispersion() > 1.08 ) return 0; | |
55 | // if (fCluster->GetM20() > 0.42 ) return 0; | |
56 | // if (fCluster->GetM02() > 0.4 ) return 0; | |
57 | // if (fCluster->GetM20() < 0 ) return 0; | |
58 | // if (fCluster->GetM02() < 0.06 ) return 0; | |
59 | Int_t istpce=0; | |
60 | Double_t dEdx=0; | |
61 | dEdx=fTrack->GetTPCsignal(); | |
62 | //if(dEdx0>75.&&dEdx0<95&&fRunnumber>141794&&fRunnumber<146861) istpce0=1; | |
63 | //if(dEdx0>56.&&dEdx0<73&&fRunnumber>151564&&fRunnumber<155385) istpce0=1; | |
64 | //if(dEdx0 > 70.&&fRunnumber<141795) istpce0=1; | |
65 | if(dEdx>70.0&&dEdx<95.0) istpce=1; | |
66 | ||
67 | if (ep>0.8&&ep<1.2&&istpce) return 1; | |
68 | else return 0; | |
69 | } | |
70 | ||
71 | //____________________________________________________________________________ | |
72 | void AliEMCALClusterParams::LoopThroughCells() const | |
73 | { | |
74 | Int_t nclusfCells=fCluster->GetNCells(); | |
75 | UShort_t *celllist; | |
76 | celllist=fCluster->GetCellsAbsId(); | |
77 | for (Int_t i=0;i<nclusfCells;i++) { | |
78 | Int_t iSupMod = -1; | |
79 | Int_t iTower = -1; | |
80 | Int_t iIphi = -1; | |
81 | Int_t iIeta = -1; | |
82 | Int_t iphi = -1; | |
83 | Int_t ieta = -1; | |
84 | ||
85 | printf("Cell %d id: %d. ",i,celllist[i]); | |
86 | ||
87 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
88 | printf("iphi,eta: %d %d, ",iIphi,iIeta); | |
89 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
90 | printf("smod %d tower %d ieta %d iphi %d \n",iSupMod,iTower,iphi,ieta); | |
91 | } | |
92 | } | |
93 | ||
94 | //____________________________________________________________________________ | |
95 | void AliEMCALClusterParams::PrintClusterParameters() const | |
96 | { | |
97 | printf("N fCells: %d Energy: %f Dispersion %f M02: %f M20: %f p/E: %f \n", | |
98 | fCluster->GetNCells(),fCluster->E(),fCluster->GetDispersion(), | |
99 | fCluster->GetM02(),fCluster->GetM20(),fTrack->Pt()/fCluster->E()); | |
100 | } | |
101 | ||
102 | //===================================== UNWEIGHTED PARAMETERS================================== | |
103 | ||
104 | //____________________________________________________________________________ | |
105 | void AliEMCALClusterParams::GetCentroid(Double_t &xback, Double_t &yback, Double_t &rback) const | |
106 | { | |
107 | Int_t nclusfCells=fCluster->GetNCells(); | |
108 | UShort_t *celllist; | |
109 | celllist=fCluster->GetCellsAbsId(); | |
110 | ||
111 | //Calculates mean x,y and r of the fCluster | |
112 | Double_t etai=0, phii=0; | |
113 | Double_t xsum=0; | |
114 | Double_t ysum=0; | |
115 | Double_t rsum=0; | |
116 | Double_t esum=0; | |
117 | ||
118 | for (Int_t i=0;i<nclusfCells;i++) { | |
119 | Int_t iSupMod = -1; | |
120 | Int_t iTower = -1; | |
121 | Int_t iIphi = -1; | |
122 | Int_t iIeta = -1; | |
123 | Int_t iphi = -1; | |
124 | Int_t ieta = -1; | |
125 | ||
126 | //printf("Cell %d id: %d. ",i,celllist[i]); | |
127 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
128 | // printf("iphi,eta: %d %d, ",iIphi,iIeta); | |
129 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
130 | //printf("smod %d tower %d ieta %d iphi %d \n",iSupMod,iTower,iphi,ieta); | |
131 | ||
132 | etai=(Double_t)ieta+1.; | |
133 | phii=(Double_t)iphi+1.; | |
134 | xsum+=etai*fCells->GetCellAmplitude(celllist[i]); | |
135 | ysum+=phii*fCells->GetCellAmplitude(celllist[i]); | |
136 | esum+=fCells->GetCellAmplitude(celllist[i]); | |
137 | rsum+=sqrt(etai*etai+phii*phii)*fCells->GetCellAmplitude(celllist[i]); | |
138 | } | |
139 | yback=ysum/esum; | |
140 | xback=xsum/esum; | |
141 | rback=rsum/esum; | |
142 | } | |
143 | ||
144 | //____________________________________________________________________________ | |
145 | Double_t AliEMCALClusterParams::GetR(Double_t x, Double_t y) const | |
146 | { | |
147 | // Takes fCluster, and cell position (x,y) and returns distance of cell from centroid | |
148 | ||
149 | Double_t rmean=-99; | |
150 | Double_t xmean=-99; | |
151 | Double_t ymean=-99; | |
152 | ||
153 | GetCentroid(xmean,ymean,rmean); | |
154 | return sqrt(pow((x-xmean),2)+pow((y-ymean),2)); | |
155 | } | |
156 | ||
157 | //_____________________________________________________________________________ | |
158 | Double_t AliEMCALClusterParams::GetRfactor() const | |
159 | { | |
160 | Int_t nclusfCells=fCluster->GetNCells(); | |
161 | UShort_t *celllist; | |
162 | celllist=fCluster->GetCellsAbsId(); | |
163 | ||
164 | //Calculates mean x,y and r of the fCluster | |
165 | Double_t etai=0, phii=0; | |
166 | Double_t rsum=0; | |
167 | Double_t esum=0; | |
168 | ||
169 | for (Int_t i=0;i<nclusfCells;i++) { | |
170 | Int_t iSupMod = -1; | |
171 | Int_t iTower = -1; | |
172 | Int_t iIphi = -1; | |
173 | Int_t iIeta = -1; | |
174 | Int_t iphi = -1; | |
175 | Int_t ieta = -1; | |
176 | ||
177 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
178 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
179 | ||
180 | etai=(Double_t)ieta+1.; | |
181 | phii=(Double_t)iphi+1.; | |
182 | esum+=fCells->GetCellAmplitude(celllist[i]); | |
183 | rsum+=GetR(etai,phii)*fCells->GetCellAmplitude(celllist[i]);//GetR! | |
184 | } | |
185 | ||
186 | rsum=rsum/esum; | |
187 | return rsum; | |
188 | } | |
189 | ||
190 | //____________________________________________________________________________ | |
191 | Double_t AliEMCALClusterParams::ElectronFraction(Double_t r, Double_t tce) const | |
192 | { | |
193 | // In determination of the 'K-factor', we measure the 'dispersion' of the [energy fraction of the cell | |
194 | // as a function of cell distance from the centroid] from the the distribution for electrons | |
195 | // the functional form for electrons is saved here. | |
196 | // tce is TOTAL FCLUSTER ENERGY, r is distance of cell to centroid. | |
197 | ||
198 | Double_t parm[20][11] = {{0.968421,-0.952347,0.196982,-495.028,495.028,0.512768,3.21555,-134.368,950.132,0.258585,-0.0819375}, | |
199 | {0.964309,-0.955084,0.212321,-401.273,401.273,0.600581,-2.44701,-43.9126,1009.65,0.19388,-0.0527486}, | |
200 | {0.949579,-0.938399,0.169353,-398.208,398.208,0.801709,-16.1048,399.842,-2319.64,0.16054,-0.044954}, | |
201 | {0.937618,-0.926665,0.173489,-291.715,291.715,0.723635,-5.65824,222.248,-1525.94,0.125206,-0.0108058}, | |
202 | {0.922977,-0.91741,0.174686,-269.155,269.155,0.777305,-0.656328,86.8073,-733.079,0.123508,-0.0159208}, | |
203 | {0.88858,-0.86354,0.102552,-280.217,280.217,0.746299,7.06471,-93.7562,323.688,0.128045,-0.0306688}, | |
204 | {0.859556,-0.83609,0.108273,-275.632,275.632,0.752012,4.78987,-66.0216,262.093,0.102209,-0.0201003}, | |
205 | {0.789254,-0.749649,0.0960812,-331.001,331.001,0.703363,6.70229,-129.452,725.871,0.0861534,-0.0127965}, | |
206 | {0.813629,-0.784791,0.0847636,-287.469,287.469,0.7138,7.15787,-142.766,820.634,0.0625165,-0.00910334}, | |
207 | {0.828972,-0.805683,0.0812484,-522.356,522.356,0.813623,1.77465,-34.6065,147.756,0.0467835,-0.00425002}, | |
208 | {0.900904,-0.890754,0.048338,-1180.84,1180.84,0.870175,0.0385033,-6.99334,19.6102,0.0232503,-0.00446416}, | |
209 | {0.881545,-0.870425,0.0382378,-750,750,0.857701,0.181933,-7.72458,12.6329,0.0198669,-0.00436687}, | |
210 | {0.893767,-0.878626,0.0465724,-750,750,0.867552,0.0868854,-9.17702,27.6996,0.013101,-0.00261109}, | |
211 | {0.88893,-0.882229,0.0896708,-1166.94,828.101,0.864913,-0.245158,-1.5693,-10.5061,0.00658265,-0.000176185}, | |
212 | {0.875546,-0.867819,0.0433809,-750,750,0.853546,0.103699,-4.53405,0.67608,0.00667051,-0.00110355}, | |
213 | {0.879879,-0.874221,0.0510928,-750,750,0.846377,0.246414,-4.28576,-22.6443,0.00537529,-0.000852894}, | |
214 | {0.889361,-0.884282,0.0475369,-750,750,0.845565,0.419865,-7.83763,-2.71567,0.00348465,-0.000407448}, | |
215 | {0.893711,-0.890189,0.0649599,-847.351,847.351,0.844507,0.379177,-8.07506,7.69156,0.0005,-2.89084e-13}, | |
216 | {0.893706,-0.891587,0.0409521,-750,750,0.841632,0.269466,-4.04288,-14.8299,0.00211674,-0.000173327}, | |
217 | {0.942361,-0.941082,0.0703767,-750,750,0.83202,0.829046,-11.3403,29.0905,0.00025,5.67537e-05}}; | |
218 | ||
219 | Double_t par[11]={0}; | |
220 | if (tce>0&&tce<0.6) { | |
221 | for (Int_t l=0;l<11;l++) | |
222 | par[l]=parm[0][l]; | |
223 | } else if (tce>0.6&&tce<0.8) { | |
224 | for (Int_t l=0;l<11;l++) | |
225 | par[l]=parm[1][l]; | |
226 | } else if (tce>0.8&&tce<1.12) { | |
227 | for (Int_t l=0;l<11;l++) | |
228 | par[l]=parm[2][l]; | |
229 | } else if (tce>1.12&&tce<1.37) { | |
230 | for (Int_t l=0;l<11;l++) | |
231 | par[l]=parm[3][l]; | |
232 | } else if (tce>1.37&&tce<1.75) { | |
233 | for (Int_t l=0;l<11;l++) | |
234 | par[l]=parm[4][l]; | |
235 | } else if (tce>1.75&&tce<2.5) { | |
236 | for (Int_t l=0;l<11;l++) | |
237 | par[l]=parm[5][l]; | |
238 | } else if (tce>2.5&&tce<3.5) { | |
239 | for (Int_t l=0;l<11;l++) | |
240 | par[l]=parm[6][l]; | |
241 | } else if (tce>3.5&&tce<4.5) { | |
242 | for (Int_t l=0;l<11;l++) | |
243 | par[l]=parm[7][l]; | |
244 | } else if (tce>4.5&&tce<5.5) { | |
245 | for (Int_t l=0;l<11;l++) | |
246 | par[l]=parm[8][l]; | |
247 | } else if (tce>5.5&&tce<8) { | |
248 | for (Int_t l=0;l<11;l++) | |
249 | par[l]=parm[9][l]; | |
250 | } else if (tce<8&&tce>12) { | |
251 | for (Int_t l=0;l<11;l++) | |
252 | par[l]=parm[10][l]; | |
253 | } else if (tce>12&&tce<20) { | |
254 | for (Int_t l=0;l<11;l++) | |
255 | par[l]=parm[11][l]; | |
256 | } else if (tce>20&&tce<40) { | |
257 | for (Int_t l=0;l<11;l++) | |
258 | par[l]=parm[12][l]; | |
259 | } else if (tce>40&&tce<63) { | |
260 | for (Int_t l=0;l<11;l++) | |
261 | par[l]=parm[13][l]; | |
262 | } else if (tce>63&&tce<88) { | |
263 | for (Int_t l=0;l<11;l++) | |
264 | par[l]=parm[14][l]; | |
265 | } else if (tce>88&&tce<113) { | |
266 | for (Int_t l=0;l<11;l++) | |
267 | par[l]=parm[15][l]; | |
268 | } else if (tce>113&&tce<138) { | |
269 | for (Int_t l=0;l<11;l++) | |
270 | par[l]=parm[16][l]; | |
271 | } else if (tce>138&&tce<163) { | |
272 | for (Int_t l=0;l<11;l++) | |
273 | par[l]=parm[17][l]; | |
274 | } else if (tce>163&&tce<188) { | |
275 | for (Int_t l=0;l<11;l++) | |
276 | par[l]=parm[18][l]; | |
277 | } else if (tce>188) { | |
278 | for (Int_t l=0;l<11;l++) | |
279 | par[l]=parm[19][l]; | |
280 | } | |
281 | ||
282 | Double_t fr=0.; | |
283 | if (r<0.1) { | |
284 | fr =par[5]+par[6]*r+par[7]*r*r+par[8]*r*r*r; | |
285 | } else if (r>0.1&&r<0.92) { | |
286 | fr = par[0]+par[1]*r+par[2]*exp(par[3]-par[4]*r); | |
287 | } else if (r>0.92&&r<6) { | |
288 | fr = par[9] + par[10]*r; | |
289 | } | |
290 | if (r>6) { | |
291 | cout<<"Crappy f(r_i) value, more than 6!"<<endl; | |
292 | } | |
293 | ||
294 | return fr; //returns fraction of energy at r,E: f(r,E), correct? | |
295 | } | |
296 | ||
297 | //____________________________________________________________________________ | |
298 | Double_t AliEMCALClusterParams::GetKfactor() const | |
299 | { | |
300 | // Determines the k-factor, uses the funtion ElectronFraction(r,totalenergyoffCluster) | |
301 | ||
302 | Double_t totalFClusterEnergy = fCluster->E(); | |
303 | //Calculates mean x,y and r of the fCluster | |
304 | Double_t etai=0, phii=0; | |
305 | ||
306 | Double_t kfactor=0; | |
307 | ||
308 | Int_t nclusfCells=fCluster->GetNCells(); | |
309 | UShort_t *celllist; | |
310 | celllist=fCluster->GetCellsAbsId(); | |
311 | ||
312 | for (Int_t i=0;i<nclusfCells;i++) { | |
313 | Int_t iSupMod = -1; | |
314 | Int_t iTower = -1; | |
315 | Int_t iIphi = -1; | |
316 | Int_t iIeta = -1; | |
317 | Int_t iphi = -1; | |
318 | Int_t ieta = -1; | |
319 | ||
320 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
321 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
322 | ||
323 | etai=(Double_t)ieta+1.; | |
324 | phii=(Double_t)iphi+1.; | |
325 | ||
326 | kfactor+=GetR(etai,phii)*pow(fCells->GetCellAmplitude(celllist[i])/(totalFClusterEnergy) - | |
327 | ElectronFraction(GetR(etai,phii),totalFClusterEnergy),2); | |
328 | } | |
329 | ||
330 | return kfactor; | |
331 | } | |
332 | ||
333 | //____________________________________________________________________________ | |
334 | Double_t AliEMCALClusterParams::GetDispersionX() const | |
335 | { | |
336 | Double_t rmean=-99; | |
337 | Double_t xmean=-99; | |
338 | Double_t ymean=-99; | |
339 | Double_t esum=0; | |
340 | Double_t dsum=0; | |
341 | ||
342 | GetCentroid(xmean,ymean,rmean); | |
343 | ||
344 | Int_t nclusfCells=fCluster->GetNCells(); | |
345 | UShort_t *celllist; | |
346 | celllist=fCluster->GetCellsAbsId(); | |
347 | ||
348 | for (Int_t i=0;i<nclusfCells;i++) { | |
349 | Int_t iSupMod = -1; | |
350 | Int_t iTower = -1; | |
351 | Int_t iIphi = -1; | |
352 | Int_t iIeta = -1; | |
353 | Int_t iphi = -1; | |
354 | Int_t ieta = -1; | |
355 | ||
356 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
357 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
358 | ||
359 | Double_t etai=(Double_t)ieta+1.; | |
360 | ||
361 | dsum+=(pow(etai-xmean,2))*fCells->GetCellAmplitude(celllist[i]); | |
362 | esum+=fCells->GetCellAmplitude(celllist[i]); | |
363 | } | |
364 | ||
365 | return sqrt(dsum/esum); | |
366 | } | |
367 | ||
368 | //____________________________________________________________________________ | |
369 | Double_t AliEMCALClusterParams::GetDispersionY() const | |
370 | { | |
371 | Double_t rmean=-99; | |
372 | Double_t xmean=-99; | |
373 | Double_t ymean=-99; | |
374 | Double_t esum=0; | |
375 | Double_t dsum=0; | |
376 | ||
377 | GetCentroid(xmean,ymean,rmean); | |
378 | ||
379 | Int_t nclusfCells=fCluster->GetNCells(); | |
380 | UShort_t *celllist; | |
381 | celllist=fCluster->GetCellsAbsId(); | |
382 | ||
383 | for (Int_t i=0;i<nclusfCells;i++) { | |
384 | Int_t iSupMod = -1; | |
385 | Int_t iTower = -1; | |
386 | Int_t iIphi = -1; | |
387 | Int_t iIeta = -1; | |
388 | Int_t iphi = -1; | |
389 | Int_t ieta = -1; | |
390 | ||
391 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
392 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
393 | ||
394 | Double_t phii=(Double_t)iphi+1.; | |
395 | ||
396 | dsum+=(pow(phii-ymean,2))*fCells->GetCellAmplitude(celllist[i]); | |
397 | esum+=fCells->GetCellAmplitude(celllist[i]); | |
398 | } | |
399 | ||
400 | return sqrt(dsum/esum); | |
401 | } | |
402 | ||
403 | //____________________________________________________________________________ | |
404 | Double_t AliEMCALClusterParams::GetDispersionMax() const | |
405 | { | |
406 | Double_t dispX = GetDispersionX(); | |
407 | Double_t dispY = GetDispersionY(); | |
408 | if (dispY > dispX) { | |
409 | return dispY; | |
410 | } | |
411 | else { | |
412 | return dispX; | |
413 | } | |
414 | } | |
415 | ||
416 | //____________________________________________________________________________ | |
417 | void AliEMCALClusterParams::GetEllipseParameters(Double_t ¶m1, Double_t ¶m2) const | |
418 | { | |
419 | Double_t sumxx=0; | |
420 | Double_t sumyy=0; | |
421 | Double_t sumx=0; | |
422 | Double_t sumy=0; | |
423 | Double_t sumxy=0; | |
424 | Double_t esum=0; | |
425 | ||
426 | Int_t nclusfCells=fCluster->GetNCells(); | |
427 | UShort_t *celllist; | |
428 | celllist=fCluster->GetCellsAbsId(); | |
429 | ||
430 | for (Int_t i=0;i<nclusfCells;i++) { | |
431 | Int_t iSupMod = -1; | |
432 | Int_t iTower = -1; | |
433 | Int_t iIphi = -1; | |
434 | Int_t iIeta = -1; | |
435 | Int_t iphi = -1; | |
436 | Int_t ieta = -1; | |
437 | ||
438 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
439 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
440 | ||
441 | Double_t phii=(Double_t)iphi+1.; | |
442 | Double_t etai=(Double_t)ieta+1.; | |
443 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
444 | ||
445 | sumxx += amp*pow(etai,2); | |
446 | sumyy += amp*pow(phii,2); | |
447 | sumx += amp*etai; | |
448 | sumy += amp*phii; | |
449 | sumxy += amp*etai*phii; | |
450 | esum+=amp; | |
451 | } | |
452 | sumxx=sumxx/esum; | |
453 | sumyy = sumyy/esum; | |
454 | sumx = sumx/esum; | |
455 | sumy = sumy/esum; | |
456 | sumxx=sumxx-sumx*sumx; | |
457 | sumyy = sumyy-sumy*sumy; | |
458 | sumxy = sumxy/esum - sumx*sumy; | |
459 | ||
460 | param1 = 0.5*(sumxx+sumyy) + sqrt(0.25*pow((sumxx-sumyy),2) + sumxy*sumxy); | |
461 | param2 = 0.5*(sumxx+sumyy) - sqrt(0.25*pow((sumxx-sumyy),2) + sumxy*sumxy); | |
462 | } | |
463 | ||
464 | //____________________________________________________________________________ | |
465 | Double_t AliEMCALClusterParams::GetDispersion() const | |
466 | { | |
467 | Double_t rmean=-99; | |
468 | Double_t xmean=-99; | |
469 | Double_t ymean=-99; | |
470 | Double_t esum=0; | |
471 | Double_t dsum=0; | |
472 | ||
473 | GetCentroid(xmean,ymean,rmean); | |
474 | ||
475 | Int_t nclusfCells=fCluster->GetNCells(); | |
476 | UShort_t *celllist; | |
477 | celllist=fCluster->GetCellsAbsId(); | |
478 | ||
479 | ||
480 | for (Int_t i=0;i<nclusfCells;i++) { | |
481 | Int_t iSupMod = -1; | |
482 | Int_t iTower = -1; | |
483 | Int_t iIphi = -1; | |
484 | Int_t iIeta = -1; | |
485 | Int_t iphi = -1; | |
486 | Int_t ieta = -1; | |
487 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
488 | ||
489 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
490 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
491 | ||
492 | Double_t phii=(Double_t)iphi+1.; | |
493 | Double_t etai=(Double_t)ieta+1.; | |
494 | ||
495 | dsum+=(pow(phii-ymean,2)+pow(etai-xmean,2))*amp; | |
496 | esum+=amp; | |
497 | } | |
498 | ||
499 | return sqrt(dsum/esum); | |
500 | } | |
501 | ||
502 | //============================LOG WEIGHTED PARAMETERS======================================================== | |
503 | ||
504 | //____________________________________________________________________________ | |
505 | void AliEMCALClusterParams::GetWeightedCentroid(Double_t &xback, Double_t &yback, Double_t &rback) const | |
506 | { | |
507 | Float_t logWeight = 4.5; | |
508 | Double_t totalFClusterEnergy=fCluster->E(); | |
509 | ||
510 | //Calculates mean x,y and r of the fCluster | |
511 | Double_t etai=0, phii=0; | |
512 | Double_t xsum=0; | |
513 | Double_t ysum=0; | |
514 | Double_t rsum=0; | |
515 | Double_t esum=0; | |
516 | ||
517 | Int_t nclusfCells=fCluster->GetNCells(); | |
518 | UShort_t *celllist; | |
519 | celllist=fCluster->GetCellsAbsId(); | |
520 | ||
521 | for (Int_t i=0;i<nclusfCells;i++) { | |
522 | Int_t iSupMod = -1; | |
523 | Int_t iTower = -1; | |
524 | Int_t iIphi = -1; | |
525 | Int_t iIeta = -1; | |
526 | Int_t iphi = -1; | |
527 | Int_t ieta = -1; | |
528 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
529 | ||
530 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
531 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
532 | ||
533 | etai=(Double_t)ieta+1.; | |
534 | phii=(Double_t)iphi+1.; | |
535 | Double_t w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy )); | |
536 | xsum+=etai*w; | |
537 | ysum+=phii*w; | |
538 | esum+=w; | |
539 | rsum+=sqrt(etai*etai+phii*phii)*w; | |
540 | ||
541 | } | |
542 | yback=ysum/esum; | |
543 | xback=xsum/esum; | |
544 | rback=rsum/esum; | |
545 | } | |
546 | ||
547 | //____________________________________________________________________________ | |
548 | Double_t AliEMCALClusterParams::GetWeightedR(Double_t x, Double_t y) const | |
549 | { | |
550 | //takes fCluster, and cell position (x,y) and returns distance of cell from | |
551 | //centroid | |
552 | Double_t rmean=-99; | |
553 | Double_t xmean=-99; | |
554 | Double_t ymean=-99; | |
555 | ||
556 | GetWeightedCentroid(xmean,ymean,rmean); | |
557 | return sqrt(pow((x-xmean),2)+pow((y-ymean),2)); | |
558 | } | |
559 | ||
560 | //_____________________________________________________________________________ | |
561 | Double_t AliEMCALClusterParams::GetWeightedRfactor() const | |
562 | { | |
563 | Double_t rsum=0; | |
564 | Double_t esum=0; | |
565 | Double_t etai=0, phii=0; | |
566 | ||
567 | Float_t logWeight = 4.5; | |
568 | ||
569 | Double_t totalFClusterEnergy=fCluster->E(); | |
570 | Int_t nclusfCells=fCluster->GetNCells(); | |
571 | UShort_t *celllist; | |
572 | celllist=fCluster->GetCellsAbsId(); | |
573 | ||
574 | for (Int_t i=0;i<nclusfCells;i++) { | |
575 | Int_t iSupMod = -1; | |
576 | Int_t iTower = -1; | |
577 | Int_t iIphi = -1; | |
578 | Int_t iIeta = -1; | |
579 | Int_t iphi = -1; | |
580 | Int_t ieta = -1; | |
581 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
582 | ||
583 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
584 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
585 | ||
586 | etai=(Double_t)ieta+1.; | |
587 | phii=(Double_t)iphi+1.; | |
588 | Double_t w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy )); | |
589 | esum+=w; | |
590 | rsum+=GetWeightedR(etai,phii)*w; | |
591 | } | |
592 | rsum=rsum/esum; | |
593 | ||
594 | return rsum; | |
595 | } | |
596 | ||
597 | //____________________________________________________________________________ | |
598 | Double_t AliEMCALClusterParams::ElectronfractionWeighted(Double_t r, Double_t tce) const | |
599 | { | |
600 | // In determination of the 'K-factor', we measure the 'dispersion' of the [energy fraction of the cell | |
601 | // as a function of cell distance from the centroid] from the the distribution for electrons | |
602 | // the functional form for electrons is saved here. | |
603 | ||
604 | Double_t parm[20][10] = {{9.65749,-25.2283,146.575,-207.63,-3.53453,86.296,-176.854,110.51,8.21828,0.097827}, | |
605 | {9.69895,-19.8042,112.42,-156.762,16.8874,-27.4961,28.4271,-9.31945,8.18897,-0.213499}, | |
606 | {9.58816,-15.0126,92.8682,-137.142,19.3093,-42.9118,59.9102,-29.9932,8.10929,-0.118502}, | |
607 | {9.66433,-17.3556,105.064,-154.824,19.3254,-43.845,63.0948,-32.7479,7.3532,0.437176}, | |
608 | {9.63564,-12.3028,76.0095,-113.329,16.337,-29.9111,42.2855,-22.8323,7.56083,0.14301}, | |
609 | {9.57174,-8.92634,58.7273,-91.3475,9.81528,2.85217,-10.6939,4.76399,7.93151,-0.320945}, | |
610 | {9.49951,-3.10053,18.9051,-26.8457,11.5665,-6.0263,3.44864,-2.37582,7.36371,-0.057224}, | |
611 | {9.43673,-2.63179,13.7039,-19.2528,9.83544,0.91246,-6.30083,2.46559,7.36203,-0.146645}, | |
612 | {9.5378,-2.35688,12.6874,-19.3792,7.64049,12.0476,-24.2138,11.5242,7.2368,-0.337072}, | |
613 | {9.52255,-0.989157,7.61241,-14.1598,6.28615,18.4146,-33.4616,15.6595,6.66628,-0.0979268}, | |
614 | {9.92079,-1.74431,9.58206,-16.8904,20.9936,-35.2392,29.6651,-9.03178,5.74954,-0.082485}, | |
615 | {9.88709,-1.11991,6.64989,-13.221,15.0524,-15.8677,9.13398,-1.97541,5.5779,-0.146283}, | |
616 | {9.91925,-1.71077,9.34349,-16.919,15.913,-16.3563,8.1726,-1.38982,5.04338,-0.120024}, | |
617 | {9.90981,-2.53254,13.2516,-21.0921,15.1033,-14.4442,6.87168,-1.18726,4.25787,0.0881321}, | |
618 | {9.7964,-0.0431799,2.45754,-9.00884,15.2986,-15.7366,8.37136,-1.66767,3.83555,0.0409475}, | |
619 | {9.82242,-0.718889,5.19571,-11.6806,15.4565,-16.0309,8.56511,-1.70815,3.46559,0.0840781}, | |
620 | {9.88062,-1.21927,7.07515,-13.7956,14.2168,-13.5604,6.90633,-1.3678,3.06548,0.162603}, | |
621 | {9.98456,-2.58626,12.6175,-20.0417,12.9291,-10.2085,4.02983,-0.606252,3.34551,0.040163}, | |
622 | {9.89892,-1.45526,8.17534,-14.9205,12.5899,-9.74255,4.08869,-0.718114,2.51632,0.245529}, | |
623 | {11.2705,-5.29437,0.810794,-0.00331736,10.2653,-0.725442,-3.84838,1.23406,3.85982,-0.188577}}; | |
624 | ||
625 | Double_t par[10]={0}; | |
626 | if (tce>0&&tce<0.6) { | |
627 | for (Int_t l=0;l<10;l++) | |
628 | par[l]=parm[0][l]; | |
629 | } else if (tce>0.6&&tce<0.8) { | |
630 | for (Int_t l=0;l<10;l++) | |
631 | par[l]=parm[1][l]; | |
632 | } else if (tce>0.8&&tce<1.12) { | |
633 | for (Int_t l=0;l<10;l++) | |
634 | par[l]=parm[2][l]; | |
635 | } else if (tce>1.12&&tce<1.37) { | |
636 | for (Int_t l=0;l<10;l++) | |
637 | par[l]=parm[3][l]; | |
638 | } else if (tce>1.37&&tce<1.75) { | |
639 | for (Int_t l=0;l<10;l++) | |
640 | par[l]=parm[4][l]; | |
641 | } else if (tce>1.75&&tce<2.5) { | |
642 | for (Int_t l=0;l<10;l++) | |
643 | par[l]=parm[5][l]; | |
644 | } else if (tce>2.5&&tce<3.5) { | |
645 | for (Int_t l=0;l<10;l++) | |
646 | par[l]=parm[6][l]; | |
647 | } else if (tce>3.5&&tce<4.5) { | |
648 | for (Int_t l=0;l<10;l++) | |
649 | par[l]=parm[7][l]; | |
650 | } else if (tce>4.5&&tce<5.5) { | |
651 | for (Int_t l=0;l<10;l++) | |
652 | par[l]=parm[8][l]; | |
653 | } else if (tce>5.5&&tce<8) { | |
654 | for (Int_t l=0;l<10;l++) | |
655 | par[l]=parm[9][l]; | |
656 | } else if (tce<8&&tce>12) { | |
657 | for (Int_t l=0;l<10;l++) | |
658 | par[l]=parm[10][l]; | |
659 | } else if (tce>12&&tce<20) { | |
660 | for (Int_t l=0;l<10;l++) | |
661 | par[l]=parm[8][l]; | |
662 | } else if (tce>20&&tce<40) { | |
663 | for (Int_t l=0;l<10;l++) | |
664 | par[l]=parm[12][l]; | |
665 | } else if (tce>40&&tce<63) { | |
666 | for (Int_t l=0;l<10;l++) | |
667 | par[l]=parm[13][l]; | |
668 | } else if (tce>63&&tce<88) { | |
669 | for (Int_t l=0;l<10;l++) | |
670 | par[l]=parm[14][l]; | |
671 | } else if (tce>88&&tce<113) { | |
672 | for (Int_t l=0;l<10;l++) | |
673 | par[l]=parm[15][l]; | |
674 | } else if (tce>113&&tce<138) { | |
675 | for (Int_t l=0;l<10;l++) | |
676 | par[l]=parm[16][l]; | |
677 | } else if (tce>138&&tce<163) { | |
678 | for (Int_t l=0;l<10;l++) | |
679 | par[l]=parm[17][l]; | |
680 | } else if (tce>163&&tce<188) { | |
681 | for (Int_t l=0;l<10;l++) | |
682 | par[l]=parm[18][l]; | |
683 | } else if (tce>188) { | |
684 | for (Int_t l=0;l<10;l++) | |
685 | par[l]=parm[19][l]; | |
686 | } | |
687 | ||
688 | Double_t fr=0.; | |
689 | if (r<0.5) { | |
690 | fr = par[4]+par[5]*r+par[6]*r*r+par[7]*r*r*r; | |
691 | } else if (r>0.5&&r<1.0) { | |
692 | fr = par[0]+par[1]*r+par[2]*r*r+par[3]*r*r*r; | |
693 | } else if (r>1.0&&r<6) { | |
694 | fr = par[8] + par[9]*r; | |
695 | } | |
696 | if (r>6) { | |
697 | cout<<"Crappy f(r_i) value, more than 6!"<<endl; | |
698 | } | |
699 | ||
700 | //cout<<"Weighter fr is: "<<fr<<" and tce,r: "<<tce<<" "<<r<<endl; | |
701 | return fr; | |
702 | } | |
703 | ||
704 | //____________________________________________________________________________ | |
705 | Double_t AliEMCALClusterParams::GetWeightedKfactor() const | |
706 | { | |
707 | Double_t logWeight = 4.5; | |
708 | //determines the k-factor | |
709 | Double_t kfactor=0; | |
710 | ||
711 | Double_t totalFClusterEnergy=fCluster->E(); | |
712 | Int_t nclusfCells=fCluster->GetNCells(); | |
713 | UShort_t *celllist; | |
714 | celllist=fCluster->GetCellsAbsId(); | |
715 | ||
716 | for (Int_t i=0;i<nclusfCells;i++) { | |
717 | Int_t iSupMod = -1; | |
718 | Int_t iTower = -1; | |
719 | Int_t iIphi = -1; | |
720 | Int_t iIeta = -1; | |
721 | Int_t iphi = -1; | |
722 | Int_t ieta = -1; | |
723 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
724 | ||
725 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
726 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
727 | ||
728 | Double_t etai=(Double_t)ieta+1.; | |
729 | Double_t phii=(Double_t)iphi+1.; | |
730 | ||
731 | Double_t w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy )); | |
732 | kfactor+=GetWeightedR(etai,phii)*pow(w - ElectronfractionWeighted(GetWeightedR(etai,phii),totalFClusterEnergy),2); | |
733 | } | |
734 | ||
735 | return kfactor; | |
736 | } | |
737 | ||
738 | //____________________________________________________________________________ | |
739 | Double_t AliEMCALClusterParams::GetWeightedDispersionX() const | |
740 | { | |
741 | Float_t logWeight = 4.5; | |
742 | // Calculates the dispersion of the shower at the origin of the RecPoint | |
743 | // in cell units - Nov 16,2006 | |
744 | ||
745 | Double_t d = 0., wtot = 0., w = 0.; | |
746 | Int_t nstat=0; | |
747 | ||
748 | // Calculates the dispersion in cell units | |
749 | Double_t etai=0, etaMean=0.0; | |
750 | ||
751 | // Calculate mean values | |
752 | Double_t totalFClusterEnergy=fCluster->E(); | |
753 | Int_t nclusfCells=fCluster->GetNCells(); | |
754 | UShort_t *celllist; | |
755 | celllist=fCluster->GetCellsAbsId(); | |
756 | ||
757 | Int_t ncell=0;//cell counter | |
758 | ||
759 | for (Int_t i=0;i<nclusfCells;i++) { | |
760 | Int_t iSupMod = -1; | |
761 | Int_t iTower = -1; | |
762 | Int_t iIphi = -1; | |
763 | Int_t iIeta = -1; | |
764 | Int_t iphi = -1; | |
765 | Int_t ieta = -1; | |
766 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
767 | ||
768 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
769 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
770 | ||
771 | ncell++; | |
772 | ||
773 | etai=(Double_t)ieta+1.; | |
774 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy )); //Calc weight | |
775 | if (w>0.0) { | |
776 | etaMean += etai*w; | |
777 | wtot += w; | |
778 | } | |
779 | } | |
780 | ||
781 | if (wtot>0) | |
782 | etaMean = etaMean/wtot; | |
783 | ||
784 | // Calculate dispersion | |
785 | // Loop over fCells in the newly created fCluster | |
786 | Int_t ncell1=0;//cell counter | |
787 | nstat=0; | |
788 | ||
789 | for (Int_t i=0;i<nclusfCells;i++) { | |
790 | Int_t iSupMod = -1; | |
791 | Int_t iTower = -1; | |
792 | Int_t iIphi = -1; | |
793 | Int_t iIeta = -1; | |
794 | Int_t iphi = -1; | |
795 | Int_t ieta = -1; | |
796 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
797 | ||
798 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
799 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
800 | ||
801 | ncell1++; | |
802 | etai=(Double_t)ieta+1.; | |
803 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy)); | |
804 | ||
805 | if (w>0.0) { | |
806 | nstat++; | |
807 | d += w*((etai-etaMean)*(etai-etaMean)); //Add squares | |
808 | } | |
809 | } | |
810 | ||
811 | if ( wtot > 0 && nstat>1) d /= wtot; | |
812 | else d = 0.; | |
813 | ||
814 | return TMath::Sqrt(d); | |
815 | } | |
816 | ||
817 | //____________________________________________________________________________ | |
818 | Double_t AliEMCALClusterParams::GetWeightedDispersionY() const | |
819 | { | |
820 | Float_t logWeight = 4.5; | |
821 | // Calculates the dispersion of the shower at the origin of the RecPoint | |
822 | // in cell units - Nov 16,2006 | |
823 | ||
824 | Double_t d = 0., wtot = 0., w = 0.; | |
825 | Int_t nstat=0; | |
826 | ||
827 | // Calculates the dispersion in cell units | |
828 | Double_t phii=0, phiMean=0.0; | |
829 | ||
830 | // Calculate mean values | |
831 | Double_t totalFClusterEnergy=fCluster->E(); | |
832 | Int_t nclusfCells=fCluster->GetNCells(); | |
833 | UShort_t *celllist; | |
834 | celllist=fCluster->GetCellsAbsId(); | |
835 | ||
836 | Int_t ncell=0;//cell counter | |
837 | ||
838 | for (Int_t i=0;i<nclusfCells;i++) { | |
839 | Int_t iSupMod = -1; | |
840 | Int_t iTower = -1; | |
841 | Int_t iIphi = -1; | |
842 | Int_t iIeta = -1; | |
843 | Int_t iphi = -1; | |
844 | Int_t ieta = -1; | |
845 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
846 | ||
847 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
848 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
849 | ||
850 | ncell++; | |
851 | ||
852 | phii=(Double_t)iphi+1.; | |
853 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy )); //Calc weight | |
854 | if (w>0.0) { | |
855 | phiMean += phii*w; | |
856 | wtot += w; | |
857 | } | |
858 | } | |
859 | ||
860 | if (wtot>0) | |
861 | phiMean = phiMean/wtot; | |
862 | ||
863 | // Calculate dispersion | |
864 | // Loop over fCells in the newly created fCluster | |
865 | Int_t ncell1=0;//cell counter | |
866 | nstat=0; | |
867 | for (Int_t i=0;i<nclusfCells;i++) { | |
868 | Int_t iSupMod = -1; | |
869 | Int_t iTower = -1; | |
870 | Int_t iIphi = -1; | |
871 | Int_t iIeta = -1; | |
872 | Int_t iphi = -1; | |
873 | Int_t ieta = -1; | |
874 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
875 | ||
876 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
877 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
878 | ||
879 | ncell1++; | |
880 | ||
881 | phii=(Double_t)iphi+1.; | |
882 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy)); | |
883 | ||
884 | if (w>0.0) { | |
885 | nstat++; | |
886 | d += w*((phii-phiMean)*(phii-phiMean)); //Add squares | |
887 | } | |
888 | } | |
889 | ||
890 | if ( wtot > 0 && nstat>1) d /= wtot; | |
891 | else d = 0.; | |
892 | ||
893 | return TMath::Sqrt(d); | |
894 | } | |
895 | ||
896 | //____________________________________________________________________________ | |
897 | Double_t AliEMCALClusterParams::GetWeightedDispersionMax() const | |
898 | { | |
899 | Double_t dispX = GetWeightedDispersionX(); | |
900 | Double_t dispY = GetWeightedDispersionY(); | |
901 | if (dispY > dispX) { | |
902 | return dispY; | |
903 | } | |
904 | else { | |
905 | return dispX; | |
906 | } | |
907 | } | |
908 | ||
909 | //____________________________________________________________________________ | |
910 | void AliEMCALClusterParams::GetWeightedEllipseParameters(Double_t ¶m1, Double_t ¶m2) const | |
911 | { | |
912 | Float_t logWeight=4.5; | |
913 | ||
914 | Double_t wtot = 0.; | |
915 | Double_t x = 0.; | |
916 | Double_t z = 0.; | |
917 | Double_t dxx = 0.; | |
918 | Double_t dzz = 0.; | |
919 | Double_t dxz = 0.; | |
920 | ||
921 | Double_t etai =0, phii=0, w=0; | |
922 | ||
923 | Int_t ncell=0;//cell counter | |
924 | ||
925 | Double_t totalFClusterEnergy=fCluster->E(); | |
926 | Int_t nclusfCells=fCluster->GetNCells(); | |
927 | UShort_t *celllist; | |
928 | celllist=fCluster->GetCellsAbsId(); | |
929 | ||
930 | for (Int_t i=0;i<nclusfCells;i++) { | |
931 | Int_t iSupMod = -1; | |
932 | Int_t iTower = -1; | |
933 | Int_t iIphi = -1; | |
934 | Int_t iIeta = -1; | |
935 | Int_t iphi = -1; | |
936 | Int_t ieta = -1; | |
937 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
938 | ||
939 | ncell++; | |
940 | etai = phii = 0.; | |
941 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
942 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
943 | ||
944 | etai=(Double_t)ieta; | |
945 | phii=(Double_t)iphi; | |
946 | ||
947 | //Weight! | |
948 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy ) ); //Energy of tower/total clus E, in GeV | |
949 | dxx += w * etai * etai; | |
950 | x += w * etai; | |
951 | dzz += w * phii * phii; | |
952 | z += w * phii; | |
953 | dxz += w * etai * phii; | |
954 | wtot += w; | |
955 | } | |
956 | ||
957 | if ( wtot > 0 ) { | |
958 | dxx /= wtot; | |
959 | x /= wtot; | |
960 | dxx -= x * x; | |
961 | dzz /= wtot; | |
962 | z /= wtot; | |
963 | dzz -= z * z; | |
964 | dxz /= wtot; | |
965 | dxz -= x * z; | |
966 | param1 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; | |
967 | param2 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; | |
968 | } else { | |
969 | param1= 0.; | |
970 | param2= 0.; | |
971 | } | |
972 | } | |
973 | ||
974 | //____________________________________________________________________________ | |
975 | Double_t AliEMCALClusterParams::GetWeightedDispersion(Double_t &dispersionback) const | |
976 | { | |
977 | Float_t logWeight = 4.5; | |
978 | // Calculates the dispersion of the shower at the origin of the RecPoint | |
979 | // in cell units - Nov 16,2006 | |
980 | ||
981 | Double_t d = 0., wtot = 0., w = 0.; | |
982 | Int_t nstat=0; | |
983 | ||
984 | ||
985 | // Calculates the dispersion in cell units | |
986 | Double_t etai, phii, etaMean=0.0, phiMean=0.0; | |
987 | ||
988 | // Calculate mean values | |
989 | Int_t ncell=0;//cell counter | |
990 | ||
991 | Double_t totalFClusterEnergy=fCluster->E(); | |
992 | Int_t nclusfCells=fCluster->GetNCells(); | |
993 | UShort_t *celllist; | |
994 | celllist=fCluster->GetCellsAbsId(); | |
995 | ||
996 | for (Int_t i=0;i<nclusfCells;i++) { | |
997 | Int_t iSupMod = -1; | |
998 | Int_t iTower = -1; | |
999 | Int_t iIphi = -1; | |
1000 | Int_t iIeta = -1; | |
1001 | Int_t iphi = -1; | |
1002 | Int_t ieta = -1; | |
1003 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
1004 | ||
1005 | ncell++; | |
1006 | etai = phii = 0.; | |
1007 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
1008 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
1009 | ||
1010 | etai=(Double_t)ieta+1.; | |
1011 | phii=(Double_t)iphi+1.; | |
1012 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy )); //Calc weight | |
1013 | if (w>0.0) { | |
1014 | phiMean += phii*w; | |
1015 | etaMean += etai*w; | |
1016 | wtot += w; | |
1017 | } | |
1018 | //digit->Delete(); | |
1019 | } | |
1020 | ||
1021 | if (wtot>0.0) { | |
1022 | phiMean = phiMean/wtot; | |
1023 | etaMean = etaMean/wtot; | |
1024 | } | |
1025 | ||
1026 | // Calculate dispersion | |
1027 | Int_t ncell1=0;//cell counter | |
1028 | nstat=0; | |
1029 | ||
1030 | for (Int_t i=0;i<nclusfCells;i++) { | |
1031 | Int_t iSupMod = -1; | |
1032 | Int_t iTower = -1; | |
1033 | Int_t iIphi = -1; | |
1034 | Int_t iIeta = -1; | |
1035 | Int_t iphi = -1; | |
1036 | Int_t ieta = -1; | |
1037 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
1038 | ||
1039 | ncell1++; | |
1040 | etai = phii = 0.; | |
1041 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
1042 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
1043 | ||
1044 | etai=(Double_t)ieta+1.; | |
1045 | phii=(Double_t)iphi+1.; | |
1046 | w = TMath::Max(0.,logWeight+TMath::Log(amp/totalFClusterEnergy)); | |
1047 | ||
1048 | if (w>0.0) { | |
1049 | nstat++; | |
1050 | d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean)); //Add squares | |
1051 | } | |
1052 | //delete digit; | |
1053 | } | |
1054 | ||
1055 | //delete digit; | |
1056 | ||
1057 | if ( wtot > 0 && nstat>1) d /= wtot; | |
1058 | else d = 0.; | |
1059 | ||
1060 | //printf("AliEMCALRecPoint::EvalDispersion() : Dispersion %f \n",d); | |
1061 | dispersionback=d; | |
1062 | return TMath::Sqrt(d); | |
1063 | ||
1064 | } | |
1065 | ||
1066 | //____________________________________________________________________________ | |
1067 | void AliEMCALClusterParams::RecalculateClusterShowerShapeParameters(Double_t &m02, Double_t &m20, Double_t &dispersion) const | |
1068 | { | |
1069 | // Calculates new center of gravity in the local EMCAL-module coordinates | |
1070 | // and tranfers into global ALICE coordinates | |
1071 | // Calculates Dispersion and main axis | |
1072 | ||
1073 | Double_t fW0=4.5; | |
1074 | ||
1075 | Int_t nstat = 0; | |
1076 | Float_t wtot = 0.; | |
1077 | Double_t eCell = 0.; | |
1078 | ||
1079 | Int_t iSupMod = -1; | |
1080 | Int_t iTower = -1; | |
1081 | Int_t iIphi = -1; | |
1082 | Int_t iIeta = -1; | |
1083 | Int_t iphi = -1; | |
1084 | Int_t ieta = -1; | |
1085 | Double_t etai = -1.; | |
1086 | Double_t phii = -1.; | |
1087 | ||
1088 | Double_t w = 0.; | |
1089 | Double_t d = 0.; | |
1090 | Double_t dxx = 0.; | |
1091 | Double_t dzz = 0.; | |
1092 | Double_t dxz = 0.; | |
1093 | Double_t xmean = 0.; | |
1094 | Double_t zmean = 0.; | |
1095 | ||
1096 | Double_t totalFClusterEnergy=fCluster->E(); | |
1097 | Int_t nclusfCells=fCluster->GetNCells(); | |
1098 | UShort_t *celllist; | |
1099 | celllist=fCluster->GetCellsAbsId(); | |
1100 | ||
1101 | //Loop on fCells | |
1102 | for (Int_t i=0;i<nclusfCells;i++) { | |
1103 | //Get from the absid the supermodule, tower and eta/phi numbers | |
1104 | ||
1105 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
1106 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
1107 | ||
1108 | Double_t amp=fCells->GetCellAmplitude(celllist[i]); | |
1109 | eCell =amp; | |
1110 | ||
1111 | if (totalFClusterEnergy > 0 && eCell > 0) { | |
1112 | w = TMath::Max( 0., fW0 + TMath::Log( eCell / totalFClusterEnergy )); | |
1113 | etai=(Double_t)ieta; | |
1114 | phii=(Double_t)iphi; | |
1115 | if (w > 0.0) { | |
1116 | wtot += w; | |
1117 | nstat++; | |
1118 | //Shower shape | |
1119 | dxx += w * etai * etai; | |
1120 | xmean+= w * etai; | |
1121 | dzz += w * phii * phii; | |
1122 | zmean+= w * phii; | |
1123 | dxz += w * etai * phii; | |
1124 | } | |
1125 | } | |
1126 | else | |
1127 | AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, fCluster->E())); | |
1128 | }//cell loop | |
1129 | ||
1130 | //Normalize to the weight | |
1131 | if (wtot > 0) { | |
1132 | xmean /= wtot; | |
1133 | zmean /= wtot; | |
1134 | } | |
1135 | else | |
1136 | AliError(Form("Wrong weight %f\n", wtot)); | |
1137 | ||
1138 | //Calculate dispersion | |
1139 | for (Int_t i=0;i<nclusfCells;i++) { | |
1140 | //Get from the absid the supermodule, tower and eta/phi numbers | |
1141 | fGeom->GetCellIndex(celllist[i],iSupMod,iTower,iIphi,iIeta); | |
1142 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
1143 | ||
1144 | eCell = fCells->GetCellAmplitude(celllist[i]); | |
1145 | ||
1146 | if (totalFClusterEnergy > 0 && eCell > 0) { | |
1147 | w = TMath::Max( 0., fW0 + TMath::Log( eCell / totalFClusterEnergy )); | |
1148 | etai=(Double_t)ieta; | |
1149 | phii=(Double_t)iphi; | |
1150 | if (w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); | |
1151 | } | |
1152 | else | |
1153 | AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, fCluster->E())); | |
1154 | }// cell loop | |
1155 | ||
1156 | //Normalize to the weigth and set shower shape parameters | |
1157 | if (wtot > 0 && nstat > 1) { | |
1158 | d /= wtot; | |
1159 | dxx /= wtot; | |
1160 | dzz /= wtot; | |
1161 | dxz /= wtot; | |
1162 | dxx -= xmean * xmean; | |
1163 | dzz -= zmean * zmean; | |
1164 | dxz -= xmean * zmean; | |
1165 | m02=(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )); | |
1166 | m20=(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )); | |
1167 | } | |
1168 | else{ | |
1169 | d=0.; | |
1170 | m02=0; | |
1171 | m20=0; | |
1172 | } | |
1173 | ||
1174 | if (d>=0) | |
1175 | dispersion=TMath::Sqrt(d); | |
1176 | else | |
1177 | dispersion=0; | |
1178 | } |