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 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 }