]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSEmcRecPoint.cxx
parameters have been redistributed; Hits2SDigits etc ... introduce
[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx
index f015b40250dd9518b15b358e810199a1a6ee7d4c..bc06a2efd6b637f5d712517dcbfdb1b4bc85cbd6 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 //_________________________________________________________________________
-// Rec Point in the PHOS EM calorimeter 
-//*-- Author : Dmitri Peressounko RRC KI
-//////////////////////////////////////////////////////////////////////////////
+//  RecPoint implementation for PHOS-EMC 
+//  An EmcRecPoint is a cluster of digits   
+//           
+//*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
+
 
 // --- ROOT system ---
 #include "TPad.h"
 
 // --- Standard library ---
 
-#include <iostream>
+#include <iostream.h> 
 
 // --- AliRoot header files ---
 
 #include "AliPHOSGeometry.h"
 #include "AliPHOSEmcRecPoint.h"
 #include "AliRun.h"
+#include "AliPHOSIndexToObject.h"
 
 ClassImp(AliPHOSEmcRecPoint)
 
@@ -45,38 +50,41 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
   fMulDigit   = 0 ;  
   fAmp   = 0. ;   
   fEnergyList = new Float_t[fMaxDigit]; 
-  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
-  fDelta     =  PHOSGeom->GetCrystalSize(0) ; 
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
+  fDelta     =  phosgeom->GetCrystalSize(0) ; 
   fW0        = W0 ;          
   fLocMaxCut = LocMaxCut ; 
   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
+  
 }
 
 //____________________________________________________________________________
-AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() 
+AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint()
 {
-  // dtor 
+  // dtor
+  if ( fEnergyList )
+    delete[] fEnergyList ; 
 }
 
 //____________________________________________________________________________
-void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
+void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy)
 {
-  // adds a digit to the digits list
-  // and accumulates the total amplitude and the multiplicity 
+  // Adds a digit to the RecPoint
+  //  and accumulates the total amplitude and the multiplicity 
   
   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
     fMaxDigit*=2 ; 
-    int * tempo = new ( int[fMaxDigit] ) ; 
+    Int_t * tempo = new ( Int_t[fMaxDigit] ) ; 
     Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;
 
     Int_t index ;     
     for ( index = 0 ; index < fMulDigit ; index++ ){
-      tempo[index] = fDigitsList[index] ;
+      tempo[index]  = fDigitsList[index] ;
       tempoE[index] = fEnergyList[index] ; 
     }
     
     delete [] fDigitsList ; 
-    fDigitsList =  new ( int[fMaxDigit] ) ;
+    fDigitsList =  new ( Int_t[fMaxDigit] ) ;
  
     delete [] fEnergyList ;
     fEnergyList =  new ( Float_t[fMaxDigit] ) ;
@@ -90,28 +98,30 @@ void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
     delete [] tempoE ; 
   } // if
   
-  fDigitsList[fMulDigit]   =  (int) &digit  ; 
-  fEnergyList[fMulDigit++] = Energy ;
+  fDigitsList[fMulDigit]   = digit.GetIndexInList()  ; 
+  fEnergyList[fMulDigit]   = Energy ;
+  fMulDigit++ ; 
   fAmp += Energy ; 
 }
 
 //____________________________________________________________________________
 Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) 
 {
+  // Tells if (true) or not (false) two digits are neighbors)
   
   Bool_t aren = kFALSE ;
   
-  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
   Int_t relid1[4] ; 
-  PHOSGeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
+  phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
 
   Int_t relid2[4] ; 
-  PHOSGeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
+  phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
   
-  Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
-  Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
+  Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
+  Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
 
-  if (( ColDiff<=1 )  && ( RowDiff <= 1 ) && (ColDiff+RowDiff > 0)) 
+  if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
     aren = kTRUE ;
   
   return aren ;
@@ -120,33 +130,35 @@ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * d
 //____________________________________________________________________________
 Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
 {
+  // Compares two RecPoints according to their position in the PHOS modules
+
   Int_t rv ; 
 
   AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; 
 
  
-  Int_t PHOSMod1 = this->GetPHOSMod() ;
-  Int_t PHOSMod2 = clu->GetPHOSMod() ;
+  Int_t phosmod1 = this->GetPHOSMod() ;
+  Int_t phosmod2 = clu->GetPHOSMod() ;
 
-  TVector3 LocPos1; 
-  this->GetLocalPosition(LocPos1) ;
-  TVector3 LocPos2;  
-  clu->GetLocalPosition(LocPos2) ;  
+  TVector3 locpos1; 
+  this->GetLocalPosition(locpos1) ;
+  TVector3 locpos2;  
+  clu->GetLocalPosition(locpos2) ;  
 
-  if(PHOSMod1 == PHOSMod2 ) {
-    Int_t rowdif = (Int_t)TMath::Ceil(LocPos1.X()/fDelta)-(Int_t)TMath::Ceil(LocPos2.X()/fDelta) ;
+  if(phosmod1 == phosmod2 ) {
+    Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/fDelta)-(Int_t)TMath::Ceil(locpos2.X()/fDelta) ;
     if (rowdif> 0) 
       rv = -1 ;
      else if(rowdif < 0) 
        rv = 1 ;
-    else if(LocPos1.Z()>LocPos2.Z()) 
+    else if(locpos1.Z()>locpos2.Z()) 
       rv = -1 ;
     else 
       rv = 1 ; 
      }
 
   else {
-    if(PHOSMod1 < PHOSMod2 ) 
+    if(phosmod1 < phosmod2 ) 
       rv = -1 ;
     else 
       rv = 1 ;
@@ -158,101 +170,107 @@ Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
 //______________________________________________________________________________
 void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
 {
-//*-*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*
-//*-*                  =========================================
-//  This member function is called when a AliPHOSRecPoint is clicked with the locator
-//
-//  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
-//  and switched off when the mouse button is released.
-//
+  // Execute action corresponding to one event
+  //  This member function is called when a AliPHOSRecPoint is clicked with the locator
+  //
+  //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
+  //  and switched off when the mouse button is released.
+  //
 
   //   static Int_t pxold, pyold;
 
-   static TGraph *  DigitGraph = 0 ;
-
-   if (!gPad->IsEditable()) return;
-
-   TH2F * Histo = 0 ;
-   TCanvas * HistoCanvas ; 
-   
-   switch (event) {
-   
-   case kButton1Down: {
-     AliPHOSDigit * digit ;
-     AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
-     Int_t iDigit;
-     Int_t relid[4] ;
-     Float_t xi[fMulDigit] ;
-     Float_t zi[fMulDigit] ;
-
-     // create the histogram for the single cluster 
-     // 1. gets histogram boundaries
-     Float_t ximax = -999. ; 
-     Float_t zimax = -999. ; 
-     Float_t ximin = 999. ; 
-     Float_t zimin = 999. ;
-     for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-       digit = (AliPHOSDigit *) fDigitsList[iDigit];
-       PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-       PHOSGeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
-       if ( xi[iDigit] > ximax )
-        ximax = xi[iDigit] ; 
-       if ( xi[iDigit] < ximin )
-        ximin = xi[iDigit] ; 
-       if ( zi[iDigit] > zimax )
-        zimax = zi[iDigit] ; 
-       if ( zi[iDigit] < zimin )
-        zimin = zi[iDigit] ;     
-     }
-     ximax += PHOSGeom->GetCrystalSize(0) / 2. ;
-     zimax += PHOSGeom->GetCrystalSize(2) / 2. ;
-     ximin -= PHOSGeom->GetCrystalSize(0) / 2. ;
-     zimin -= PHOSGeom->GetCrystalSize(2) / 2. ;
-     Int_t xdim = (int)( (ximax - ximin ) / PHOSGeom->GetCrystalSize(0) + 0.5  ) ; 
-     Int_t zdim = (int)( (zimax - zimin ) / PHOSGeom->GetCrystalSize(2) + 0.5 ) ;
-     // 2. gets the histogram title
-
-     Text_t title[100] ; 
-     sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
   
-     if (!Histo) {
-       delete Histo ; 
-       Histo = 0 ; 
-     }
-     Histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
-
-     Float_t x, z ; 
-     for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-       digit = (AliPHOSDigit *) fDigitsList[iDigit];
-       PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-       PHOSGeom->RelPosInModule(relid, x, z);
-       Histo->Fill(x, z, fEnergyList[iDigit] ) ;
-     }
-
-     if (!DigitGraph) {
-       DigitGraph = new TGraph(fMulDigit,xi,zi);
-       DigitGraph-> SetMarkerStyle(5) ; 
-       DigitGraph-> SetMarkerSize(1.) ;
-       DigitGraph-> SetMarkerColor(1) ;
-       DigitGraph-> Paint("P") ;
-     }
-
-     Print() ;
-     HistoCanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
-     HistoCanvas->Draw() ; 
-     Histo->Draw("lego1") ; 
-
-     break;
-   }
-
-   case kButton1Up: 
-     if (DigitGraph) {
-       delete DigitGraph  ;
-       DigitGraph = 0 ;
-     }
-     break;
+  static TGraph *  digitgraph = 0 ;
+  
+  if (!gPad->IsEditable()) return;
+  
+  TH2F * histo = 0 ;
+  TCanvas * histocanvas ; 
+  
+  switch (event) {
+    
+  case kButton1Down: {
+    AliPHOSDigit * digit ;
+    AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
+    Int_t iDigit;
+    Int_t relid[4] ;
+    
+    const Int_t kMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ; 
+    Float_t * xi = new Float_t[kMulDigit] ; 
+    Float_t * zi = new Float_t[kMulDigit] ; 
+    
+    // create the histogram for the single cluster 
+    // 1. gets histogram boundaries
+    Float_t ximax = -999. ; 
+    Float_t zimax = -999. ; 
+    Float_t ximin = 999. ; 
+    Float_t zimin = 999. ;
+    
+    for(iDigit=0; iDigit<kMulDigit; iDigit++) {
+      digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
+      phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+      phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
+      if ( xi[iDigit] > ximax )
+       ximax = xi[iDigit] ; 
+      if ( xi[iDigit] < ximin )
+       ximin = xi[iDigit] ; 
+      if ( zi[iDigit] > zimax )
+       zimax = zi[iDigit] ; 
+      if ( zi[iDigit] < zimin )
+       zimin = zi[iDigit] ;     
+    }
+    ximax += phosgeom->GetCrystalSize(0) / 2. ;
+    zimax += phosgeom->GetCrystalSize(2) / 2. ;
+    ximin -= phosgeom->GetCrystalSize(0) / 2. ;
+    zimin -= phosgeom->GetCrystalSize(2) / 2. ;
+    Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5  ) ; 
+    Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
+    
+    // 2. gets the histogram title
+    
+    Text_t title[100] ; 
+    sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
+    
+    if (!histo) {
+      delete histo ; 
+      histo = 0 ; 
+    }
+    histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
+    
+    Float_t x, z ; 
+    for(iDigit=0; iDigit<kMulDigit; iDigit++) {
+      digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
+      phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+      phosgeom->RelPosInModule(relid, x, z);
+      histo->Fill(x, z, fEnergyList[iDigit] ) ;
+    }
+    
+    if (!digitgraph) {
+      digitgraph = new TGraph(kMulDigit,xi,zi);
+      digitgraph-> SetMarkerStyle(5) ; 
+      digitgraph-> SetMarkerSize(1.) ;
+      digitgraph-> SetMarkerColor(1) ;
+      digitgraph-> Paint("P") ;
+    }
+    
+    Print() ;
+    histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
+    histocanvas->Draw() ; 
+    histo->Draw("lego1") ; 
+    
+    delete[] xi ; 
+    delete[] zi ; 
+    
+    break;
+  }
+  
+  case kButton1Up: 
+    if (digitgraph) {
+      delete digitgraph  ;
+      digitgraph = 0 ;
+    }
+    break;
   
    }
 }
@@ -260,82 +278,195 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
 //____________________________________________________________________________
 Float_t  AliPHOSEmcRecPoint::GetDispersion() 
 {
-  Float_t D    = 0 ;
+  // Calculates the dispersion of the shower at the origine of the RecPoint
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
+  Float_t d    = 0 ;
   Float_t wtot = 0 ;
 
-  TVector3 LocPos;
-  GetLocalPosition(LocPos);
-  Float_t x = LocPos.X() ;
-  Float_t z = LocPos.Z() ;
-  //  Int_t i = GetPHOSMod() ;
+  TVector3 locpos;
+  GetLocalPosition(locpos);
+  Float_t x = locpos.X() ;
+  Float_t z = locpos.Z() ;
 
   AliPHOSDigit * digit ;
-  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
   
   Int_t iDigit;
-  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-    digit = (AliPHOSDigit *) fDigitsList[iDigit];
+  for(iDigit=0; iDigit < fMulDigit; iDigit++) {
+    digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
     Int_t relid[4] ;
     Float_t xi ;
     Float_t zi ;
-    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-    PHOSGeom->RelPosInModule(relid, xi, zi);
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xi, zi);
     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
-    D += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
+    d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
     wtot+=w ;
   }
 
-  D /= wtot ;
+  d /= wtot ;
 
-  return TMath::Sqrt(D) ;
+  return TMath::Sqrt(d) ;
+}
+//______________________________________________________________________________
+Float_t AliPHOSEmcRecPoint::CoreEnergy()
+{
+  //This function calculates energy in the core, 
+  //i.e. within radius rad = 3cm. Beyond this radius
+  //in accoradnce with shower profile energy deposition 
+  // should be less than 2%
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
+  Float_t eCore = 0 ;
+  Float_t coreRadius = 3 ;
+
+  TVector3 locpos;
+  GetLocalPosition(locpos);
+  Float_t x = locpos.X() ;
+  Float_t z = locpos.Z() ;
+
+  AliPHOSDigit * digit ;
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
+  
+  Int_t iDigit;
+  for(iDigit=0; iDigit < fMulDigit; iDigit++) {
+    digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
+    Int_t relid[4] ;
+    Float_t xi ;
+    Float_t zi ;
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xi, zi);    
+    Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(zi-z)*(zi-z)) ;
+    if(distance < coreRadius)
+      eCore += fEnergyList[iDigit] ;
+  }
+
+return eCore ;
 }
 
 //____________________________________________________________________________
 void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
 {
-  Float_t wtot = 0. ;
-  Float_t x    = 0.;
-  Float_t z    = 0.;
-  Float_t Dxx  = 0.;
-  Float_t Dzz  = 0.;
-  Float_t Dxz  = 0.;
+  // Calculates the axis of the shower ellipsoid
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
+  Double_t wtot = 0. ;
+  Double_t x    = 0.;
+  Double_t z    = 0.;
+  Double_t dxx  = 0.;
+  Double_t dzz  = 0.;
+  Double_t dxz  = 0.;
 
   AliPHOSDigit * digit ;
-  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
   Int_t iDigit;
 
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-    digit = (AliPHOSDigit *) fDigitsList[iDigit];
+    digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
     Int_t relid[4] ;
     Float_t xi ;
     Float_t zi ;
-    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-    PHOSGeom->RelPosInModule(relid, xi, zi);
-    Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
-    Dxx  += w * xi * xi ;
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xi, zi);
+    Double_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+    dxx  += w * xi * xi ;
     x    += w * xi ;
-    Dzz  += w * zi * zi ;
+    dzz  += w * zi * zi ;
     z    += w * zi ; 
-    Dxz  += w * xi * zi ; 
+    dxz  += w * xi * zi ; 
     wtot += w ;
   }
-
-  Dxx /= wtot ;
+  dxx /= wtot ;
   x   /= wtot ;
-  Dxx -= x * x ;
-  Dzz /= wtot ;
+  dxx -= x * x ;
+  dzz /= wtot ;
   z   /= wtot ;
-  Dzz -= z * z ;
-  Dxz /= wtot ;
-  Dxz -= x * z ;
+  dzz -= z * z ;
+  dxz /= wtot ;
+  dxz -= x * z ;
+
+//   //Apply correction due to non-perpendicular incidence
+//   Double_t CosX ;
+//   Double_t CosZ ;
+//   Double_t DistanceToIP= (Double_t ) ((AliPHOSGeometry *) fGeom)->GetIPtoCrystalSurface() ;
+  
+//   CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ;
+//   CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ;
+
+//   dxx = dxx/(CosX*CosX) ;
+//   dzz = dzz/(CosZ*CosZ) ;
+//   dxz = dxz/(CosX*CosZ) ;
+
 
-  lambda[0] = TMath::Sqrt( 0.5 * (Dxx + Dzz) + TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
-  lambda[1] = TMath::Sqrt( 0.5 * (Dxx + Dzz) - TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
+  lambda[0] =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
+  if(lambda[0] > 0)
+    lambda[0] = TMath::Sqrt(lambda[0]) ;
+
+  lambda[1] =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
+  if(lambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
+    lambda[1] = TMath::Sqrt(lambda[1]) ;
+  else
+    lambda[1]= 0. ;
+}
+
+//____________________________________________________________________________
+void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
+{
+  // Calculates the center of gravity in the local PHOS-module coordinates 
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
+  if( fLocPos.X() < 1000000.) { // already evaluated
+   LPos = fLocPos ;
+   return ;
+  }
+
+  Float_t wtot = 0. ;
+  Int_t relid[4] ;
+
+  Float_t x = 0. ;
+  Float_t z = 0. ;
+  
+  AliPHOSDigit * digit ;
+
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
+
+  Int_t iDigit;
+
+
+  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+    digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) );
+
+    Float_t xi ;
+    Float_t zi ;
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xi, zi);
+    Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
+    x    += xi * w ;
+    z    += zi * w ;
+    wtot += w ;
+
+  }
+
+  x /= wtot ;
+  z /= wtot ;
+  fLocPos.SetX(x)  ;
+  fLocPos.SetY(0.) ;
+  fLocPos.SetZ(z)  ;
+
+  LPos = fLocPos ;
 }
 
 //____________________________________________________________________________
 Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
 {
+  // Finds the maximum energy in the cluster
+  
   Float_t menergy = 0. ;
 
   Int_t iDigit;
@@ -349,8 +480,10 @@ Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
 }
 
 //____________________________________________________________________________
-Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) 
+Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) 
 {
+  // Calculates the multiplicity of digits with energy larger than H*energy 
+  
   Int_t multipl   = 0 ;
   Int_t iDigit ;
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
@@ -362,8 +495,13 @@ Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H)
 }
 
 //____________________________________________________________________________
-Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(int *  maxAt, Float_t * maxAtEnergy) 
+Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) 
 { 
+  // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
+  //  energy difference between two local maxima
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
   AliPHOSDigit * digit ;
   AliPHOSDigit * digitN ;
   
@@ -371,28 +509,28 @@ Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(int *  maxAt, Float_t * maxAtEner
   Int_t iDigitN ;
   Int_t iDigit ;
 
-  for(iDigit=0; iDigit<fMulDigit; iDigit++){
-    maxAt[iDigit] = fDigitsList[iDigit] ;
+  for(iDigit = 0; iDigit < fMulDigit; iDigit++){
+    maxAt[iDigit] = (Int_t) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
   }
   
-  for(iDigit=0 ; iDigit<fMulDigit; iDigit++) {   
+  for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
     if(maxAt[iDigit] != -1) {
       digit = (AliPHOSDigit *) maxAt[iDigit] ;
-         
-      for(iDigitN=0; iDigitN<fMulDigit; iDigitN++) {   
-       digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ; 
+          
+      for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {       
+       digitN = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigitN]) ) ; 
        
        if ( AreNeighbours(digit, digitN) ) {
          if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
            maxAt[iDigitN] = -1 ;
-           //but may be digit is not local max too ?
-           if(fEnergyList[iDigit]<fEnergyList[iDigitN]+fLocMaxCut) 
+           // but may be digit too is not local max ?
+           if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut) 
              maxAt[iDigit] = -1 ;
          }
          else {
            maxAt[iDigit] = -1 ;
-           //but may be digitN is not local max too ?
-           if(fEnergyList[iDigit] >fEnergyList[iDigitN]-fLocMaxCut) 
+           // but may be digitN too is not local max ?
+           if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut) 
              maxAt[iDigitN] = -1 ; 
          } 
        } // if Areneighbours
@@ -401,7 +539,7 @@ Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(int *  maxAt, Float_t * maxAtEner
   } // while digit
   
   iDigitN = 0 ;
-  for(iDigit=0; iDigit<fMulDigit; iDigit++) { 
+  for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
     if(maxAt[iDigit] != -1){
       maxAt[iDigitN] = maxAt[iDigit] ;
       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
@@ -411,55 +549,11 @@ Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(int *  maxAt, Float_t * maxAtEner
   return iDigitN ;
 }
 
-//____________________________________________________________________________
-void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
-{
-  if( fLocPos.X() < 1000000.) { //allready evaluated
-   LPos = fLocPos ;
-   return ;
-  }
-
-  Float_t wtot = 0. ;
-  Int_t relid[4] ;
-
-  Float_t x = 0. ;
-  Float_t z = 0. ;
-  
-  AliPHOSDigit * digit ;
-
-  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
-
-  Int_t iDigit;
-
-
-  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-    digit = (AliPHOSDigit *) fDigitsList[iDigit];
-
-    Float_t xi ;
-    Float_t zi ;
-    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-    PHOSGeom->RelPosInModule(relid, xi, zi);
-    Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
-    x    += xi * w ;
-    z    += zi * w ;
-    wtot += w ;
-
-  }
-
-  x /= wtot ;
-  z /= wtot ;
-  fLocPos.SetX(x)  ;
-  fLocPos.SetY(0.) ;
-  fLocPos.SetZ(z)  ;
-
-  LPos = fLocPos ;
-}
 
 // //____________________________________________________________________________
 // AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) 
 // {
-//   int * DL = Clu.GetDigitsList() ; 
+//   int * dl = Clu.GetDigitsList() ; 
  
 //  if(fDigitsList) 
 //     delete fDigitsList  ;
@@ -469,15 +563,15 @@ void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
 //   Int_t iDigit;
 
 //   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-//     digit = (AliPHOSDigit *) DL[iDigit];
+//     digit = (AliPHOSDigit *) dl[iDigit];
 //     AddDigit(*digit) ;
 //   }
 
 //   fAmp       = Clu.GetTotalEnergy() ;
 //   fGeom      = Clu.GetGeom() ;
-//   TVector3 LocPos;
-//   Clu.GetLocalPosition(LocPos) ;
-//   fLocPos    = LocPos;
+//   TVector3 locpos;
+//   Clu.GetLocalPosition(locpos) ;
+//   fLocPos    = locpos;
 //   fMulDigit       = Clu.GetMultiplicity() ;
 //   fMaxDigit  = Clu.GetMaximumMultiplicity() ;
 //   fPHOSMod   = Clu.GetPHOSMod() ;
@@ -485,7 +579,7 @@ void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
 //   fDelta     = Clu.GetDelta() ;
 //   fLocMaxCut = Clu.GetLocMaxCut() ;
   
-//   delete DL ; 
+//   delete dl ; 
  
 //   return *this ;
 // }
@@ -493,20 +587,24 @@ void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
 //____________________________________________________________________________
 void AliPHOSEmcRecPoint::Print(Option_t * option) 
 {
+  // Print the list of digits belonging to the cluster
+  
   cout << "AliPHOSEmcRecPoint: " << endl ;
 
   AliPHOSDigit * digit ; 
   Int_t iDigit;
-  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
+  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
 
   Float_t xi ;
   Float_t zi ;
   Int_t relid[4] ; 
+  
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-    digit = (AliPHOSDigit *) fDigitsList[iDigit];
-    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-    PHOSGeom->RelPosInModule(relid, xi, zi);
+    digit = please->GimeDigit( fDigitsList[iDigit] ) ; 
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xi, zi);
     cout << " Id = " << digit->GetId() ;  
     cout << "   module  = " << relid[0] ;  
     cout << "   x  = " << xi ;  
@@ -515,29 +613,29 @@ void AliPHOSEmcRecPoint::Print(Option_t * option)
   }
   cout << "       Multiplicity    = " << fMulDigit  << endl ;
   cout << "       Cluster Energy  = " << fAmp << endl ;
-  
+  cout << "       Stored at position " << GetIndexInList() << endl ; 
 }
-
 //______________________________________________________________________________
-void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
-{
-   // Stream an object of class AliPHOSEmcRecPoint.
-
-   if (R__b.IsReading()) {
-      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
-      AliPHOSRecPoint::Streamer(R__b);
-      R__b >> fDelta;
-      R__b >> fLocMaxCut;
-      R__b.ReadArray(fEnergyList);
-      R__b >> fW0;
-   } else {
-      R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
-      AliPHOSRecPoint::Streamer(R__b);
-      R__b << fDelta;
-      R__b << fLocMaxCut;
-      R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
-      R__b << fW0;
-   }
-}
-
-
+// void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
+// {
+//    // Stream an object of class AliPHOSEmcRecPoint.
+
+//    if (R__b.IsReading()) {
+//       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+//       AliPHOSRecPoint::Streamer(R__b);
+//       R__b >> fDelta;
+//       fEnergyList = new Float_t[fMulDigit] ; 
+//       R__b.ReadFastArray(fEnergyList, fMulDigit);
+//       R__b >> fLocMaxCut;
+//       R__b >> fW0;
+//    } else {
+//       R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
+//       AliPHOSRecPoint::Streamer(R__b);
+//       R__b << fDelta;
+//       R__b.WriteFastArray(fEnergyList, fMulDigit);
+//       R__b << fLocMaxCut;
+//       R__b << fW0;
+//    }
+// }