]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSEmcRecPoint.cxx
Include AliDecayer.h instead of GenTypeDefs.h
[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx
index df8047056ea31f32670eb6eab24a250172bcd42f..ac68121bdc5338720df6ed16622cd7835b5ef826 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 "AliGenerator.h"
 #include "AliPHOSGeometry.h"
 #include "AliPHOSEmcRecPoint.h"
 #include "AliRun.h"
+#include "AliPHOSIndexToObject.h"
 
 ClassImp(AliPHOSEmcRecPoint)
 
@@ -45,38 +51,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) ; 
+  fGeom = AliPHOSGeometry::GetInstance() ;
+  fDelta     =  ((AliPHOSGeometry *) fGeom)->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,14 +99,16 @@ 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 ) 
+Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const
 {
+  // Tells if (true) or not (false) two digits are neighbors)
   
   Bool_t aren = kFALSE ;
   
@@ -118,14 +129,16 @@ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * d
 }
 
 //____________________________________________________________________________
-Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
+Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const
 {
+  // 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 phosmod1 = GetPHOSMod() ;
   Int_t phosmod2 = clu->GetPHOSMod() ;
 
   TVector3 locpos1; 
@@ -158,7 +171,7 @@ Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
 //______________________________________________________________________________
 void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
 {
-  //Execute action corresponding to one event
+  // 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    
@@ -167,98 +180,109 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
 
   //   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;
   
    }
 }
 
 //____________________________________________________________________________
-Float_t  AliPHOSEmcRecPoint::GetDispersion() 
+Float_t  AliPHOSEmcRecPoint::GetDispersion() const
 {
+  // Calculates the dispersion of the shower at the origine of the RecPoint
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
   Float_t d    = 0 ;
   Float_t wtot = 0 ;
 
@@ -266,14 +290,13 @@ Float_t  AliPHOSEmcRecPoint::GetDispersion()
   GetLocalPosition(locpos);
   Float_t x = locpos.X() ;
   Float_t z = locpos.Z() ;
-  //  Int_t i = GetPHOSMod() ;
 
   AliPHOSDigit * digit ;
   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 ;
@@ -288,29 +311,69 @@ Float_t  AliPHOSEmcRecPoint::GetDispersion()
 
   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 ;
   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 ) ) ;
+    Double_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
     dxx  += w * xi * xi ;
     x    += w * xi ;
     dzz  += w * zi * zi ;
@@ -318,7 +381,6 @@ void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
     dxz  += w * xi * zi ; 
     wtot += w ;
   }
-
   dxx /= wtot ;
   x   /= wtot ;
   dxx -= x * x ;
@@ -328,13 +390,109 @@ void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
   dxz /= wtot ;
   dxz -= x * z ;
 
-  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 ) ) ;
+//   //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] =  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::EvalAll(void )
+{
+  AliPHOSRecPoint::EvalAll() ;
+  EvalLocalPosition() ;
+}
+//____________________________________________________________________________
+void AliPHOSEmcRecPoint::EvalLocalPosition(void )
+{
+  // Calculates the center of gravity in the local PHOS-module coordinates 
+
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
+  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 ;
+
+  // Correction for the depth of the shower starting point (TDR p 127)  
+  Float_t para = 0.925 ; 
+  Float_t parb = 6.52 ; 
+
+  Float_t xo,yo,zo ; //Coordinates of the origin
+  gAlice->Generator()->GetOrigin(xo,yo,zo) ;
+
+  Float_t phi = phosgeom->GetPHOSAngle(relid[0]) ;
+
+  //Transform to the local ref.frame
+  Float_t xoL,yoL ;
+  xoL = xo*TMath::Cos(phi)-yo*TMath::Sin(phi) ;
+  yoL = xo*TMath::Sin(phi)+yo*TMath::Cos(phi) ;
+  
+  Float_t radius = TMath::Sqrt((xoL-x)*(xoL-x)+
+                               (phosgeom->GetIPtoCrystalSurface()-yoL)*(phosgeom->GetIPtoCrystalSurface()-yoL)+
+                               (zo-z)*(zo-z));
+  Float_t incidencephi = TMath::ATan((x-xoL ) / radius) ; 
+  Float_t incidencetheta = TMath::ATan((z-zo) / radius) ;
+  Float_t depthx =  ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencephi) ; 
+  Float_t depthz =  ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencetheta) ; 
+  
+
+  fLocPos.SetX(x - depthx)  ;
+  fLocPos.SetY(0.) ;
+  fLocPos.SetZ(z - depthz)  ;
+
 }
 
 //____________________________________________________________________________
-Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
+Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) const
 {
+  // Finds the maximum energy in the cluster
+  
   Float_t menergy = 0. ;
 
   Int_t iDigit;
@@ -348,8 +506,10 @@ Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
 }
 
 //____________________________________________________________________________
-Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) 
+Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) const
 {
+  // Calculates the multiplicity of digits with energy larger than H*energy 
+  
   Int_t multipl   = 0 ;
   Int_t iDigit ;
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
@@ -361,8 +521,13 @@ Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H)
 }
 
 //____________________________________________________________________________
-Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) 
+Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) const
 { 
+  // 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,15 +536,15 @@ Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEn
   Int_t iDigit ;
 
   for(iDigit = 0; iDigit < fMulDigit; iDigit++){
-    maxAt[iDigit] = fDigitsList[iDigit] ;
+    maxAt[iDigit] = (Int_t) ( please->GimeDigit(fDigitsList[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] ; 
+       digitN = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigitN]) ) ; 
        
        if ( AreNeighbours(digit, digitN) ) {
          if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
@@ -410,50 +575,6 @@ Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEn
   return iDigitN ;
 }
 
-//____________________________________________________________________________
-void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
-{
-  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 *) 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) 
@@ -472,7 +593,7 @@ void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
 //     AddDigit(*digit) ;
 //   }
 
-//   fAmp       = Clu.GetTotalEnergy() ;
+//   fAmp       = Clu.GetEnergy() ;
 //   fGeom      = Clu.GetGeom() ;
 //   TVector3 locpos;
 //   Clu.GetLocalPosition(locpos) ;
@@ -492,6 +613,8 @@ 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 ; 
@@ -501,9 +624,11 @@ void AliPHOSEmcRecPoint::Print(Option_t * option)
   Float_t xi ;
   Float_t zi ;
   Int_t relid[4] ; 
+  
+  AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
+
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
-    digit = (AliPHOSDigit *) fDigitsList[iDigit];
+    digit = please->GimeDigit( fDigitsList[iDigit] ) ; 
     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
     phosgeom->RelPosInModule(relid, xi, zi);
     cout << " Id = " << digit->GetId() ;  
@@ -514,29 +639,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;
+//    }
+// }