]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/EMCAL/AliEMCALClusterParams.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Move GenerateFixedBinArray to AliAnalysisTaskEmcal and change data type from Float...
[u/mrichter/AliRoot.git] / PWG / EMCAL / AliEMCALClusterParams.cxx
CommitLineData
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 20using std::endl;
21using std::cout;
22using std::cerr;
58c00d93 23
24ClassImp(AliEMCALClusterParams)
25
26//____________________________________________________________________________
27AliEMCALClusterParams::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//____________________________________________________________________________
41Double_t AliEMCALClusterParams::GetPe() const
42{
43 Double_t pe=fTrack->Pt()/fCluster->E();
44 return pe;
45}
46
47//____________________________________________________________________________
48Int_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//____________________________________________________________________________
72void 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//____________________________________________________________________________
95void 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//____________________________________________________________________________
105void 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//____________________________________________________________________________
145Double_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//_____________________________________________________________________________
158Double_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//____________________________________________________________________________
191Double_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//____________________________________________________________________________
298Double_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//____________________________________________________________________________
334Double_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//____________________________________________________________________________
369Double_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//____________________________________________________________________________
404Double_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//____________________________________________________________________________
417void AliEMCALClusterParams::GetEllipseParameters(Double_t &param1, Double_t &param2) 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//____________________________________________________________________________
465Double_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//____________________________________________________________________________
505void 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//____________________________________________________________________________
548Double_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//_____________________________________________________________________________
561Double_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//____________________________________________________________________________
598Double_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//____________________________________________________________________________
705Double_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//____________________________________________________________________________
739Double_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//____________________________________________________________________________
818Double_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//____________________________________________________________________________
897Double_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//____________________________________________________________________________
910void AliEMCALClusterParams::GetWeightedEllipseParameters(Double_t &param1, Double_t &param2) 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//____________________________________________________________________________
975Double_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//____________________________________________________________________________
1067void 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}