AliPhysics  4549ede (4549ede)
AliEMCALClusterParams.cxx
Go to the documentation of this file.
1 // $Id$
2 
3 //_________________________________________________________________________
4 // Class for easier handling of track-cluster electron PID.
5 //
6 // Author: Tomas Aronsson (Yale)
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"
20 using std::endl;
21 using std::cout;
22 using std::cerr;
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 &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 //____________________________________________________________________________
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 &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 //____________________________________________________________________________
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 }
AliEMCALClusterParams::GetWeightedDispersionMax
Double_t GetWeightedDispersionMax() const
Definition: AliEMCALClusterParams.cxx:897
Double_t
double Double_t
Definition: External.C:58
AliEMCALClusterParams::GetWeightedCentroid
void GetWeightedCentroid(Double_t &xback, Double_t &yback, Double_t &rback) const
Definition: AliEMCALClusterParams.cxx:505
AliEMCALClusterParams::GetDispersionY
Double_t GetDispersionY() const
Definition: AliEMCALClusterParams.cxx:369
AliEMCALClusterParams::GetWeightedRfactor
Double_t GetWeightedRfactor() const
Definition: AliEMCALClusterParams.cxx:561
AliEMCALClusterParams::GetWeightedEllipseParameters
void GetWeightedEllipseParameters(Double_t &param1, Double_t &param2) const
Definition: AliEMCALClusterParams.cxx:910
AliEMCALClusterParams::fCluster
AliESDCaloCluster * fCluster
Definition: AliEMCALClusterParams.h:55
AliEMCALClusterParams::GetWeightedDispersionX
Double_t GetWeightedDispersionX() const
Definition: AliEMCALClusterParams.cxx:739
AliEMCALClusterParams::GetWeightedDispersionY
Double_t GetWeightedDispersionY() const
Definition: AliEMCALClusterParams.cxx:818
AliEMCALClusterParams::GetWeightedDispersion
Double_t GetWeightedDispersion(Double_t &dispersionback) const
Definition: AliEMCALClusterParams.cxx:975
Int_t
int Int_t
Definition: External.C:63
AliEMCALClusterParams::GetCentroid
void GetCentroid(Double_t &xback, Double_t &yback, Double_t &rback) const
Definition: AliEMCALClusterParams.cxx:105
Float_t
float Float_t
Definition: External.C:68
AliEMCALClusterParams::RecalculateClusterShowerShapeParameters
void RecalculateClusterShowerShapeParameters(Double_t &m02, Double_t &m20, Double_t &dispersion) const
Definition: AliEMCALClusterParams.cxx:1067
AliEMCALClusterParams::GetEllipseParameters
void GetEllipseParameters(Double_t &param1, Double_t &param2) const
Definition: AliEMCALClusterParams.cxx:417
AliEMCALClusterParams::GetWeightedR
Double_t GetWeightedR(Double_t x, Double_t y) const
Definition: AliEMCALClusterParams.cxx:548
AliEMCALClusterParams.h
AliEMCALClusterParams::ElectronfractionWeighted
Double_t ElectronfractionWeighted(Double_t r, Double_t tce) const
Definition: AliEMCALClusterParams.cxx:598
AliEMCALClusterParams::GetR
Double_t GetR(Double_t x, Double_t y) const
Definition: AliEMCALClusterParams.cxx:145
AliEMCALClusterParams::GetDispersionMax
Double_t GetDispersionMax() const
Definition: AliEMCALClusterParams.cxx:404
AliEMCALClusterParams::GetWeightedKfactor
Double_t GetWeightedKfactor() const
Definition: AliEMCALClusterParams.cxx:705
UShort_t
unsigned short UShort_t
Definition: External.C:28
AliEMCALClusterParams::fGeom
AliEMCALGeometry * fGeom
Definition: AliEMCALClusterParams.h:56
AliEMCALClusterParams::ElectronFraction
Double_t ElectronFraction(Double_t r, Double_t tce) const
Definition: AliEMCALClusterParams.cxx:191
AliEMCALClusterParams::GetDispersionX
Double_t GetDispersionX() const
Definition: AliEMCALClusterParams.cxx:334