PWG/EMCAL/AliEMCALClusterParams.cxx

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