]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALRecPoint.cxx
Added data members to store energy deposit per primary, parent track in EMCALDigit...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecPoint.cxx
index 686d7f433eb0ca001b6dafbe3b885210ea140e86..7ef32f40f9ebd8de5d49f8e04144c6051a1c55a3 100644 (file)
 //*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04
 
 // --- ROOT system ---
-#include <Riostream.h>
+class Riostream;
 #include <TPad.h>
-#include <TGraph.h>
-#include <TPaveText.h>
+class TGraph;
+class TPaveText;
 #include <TClonesArray.h>
 #include <TMath.h>
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
-#include "AliGenerator.h"
+//#include "AliGenerator.h"
+class AliGenerator;
 #include "AliRunLoader.h"
 #include "AliRun.h"
-#include "AliEMCAL.h"
+class AliEMCAL;
 #include "AliEMCALLoader.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCALHit.h"
@@ -45,52 +46,113 @@ ClassImp(AliEMCALRecPoint)
 
 //____________________________________________________________________________
 AliEMCALRecPoint::AliEMCALRecPoint()
-  : AliRecPoint()
+  : AliRecPoint(),
+    fGeomPtr(0),
+    fClusterType(-1),
+    fCoreEnergy(0),
+    fDispersion(0),
+    fEnergyList(0),
+    fTimeList(0),
+    fAbsIdList(0),
+    fTime(0.),
+    fCoreRadius(10),  //HG check this
+    fDETracksList(0),
+    fMulParent(0),
+    fMaxParent(0),
+    fParentsList(0),
+    fDEParentsList(0),
+    fSuperModuleNumber(0)
 {
   // ctor
-  fClusterType = -1;
   fMaxTrack = 0 ;
-  fMulDigit   = 0 ;  
-  fMaxParent = 0;
-  fMulParent = 0;
+  fMulDigit = 0 ;
   fAmp   = 0. ;   
-  fCoreEnergy = 0 ; 
-  fEnergyList = 0 ;
-  fTimeList = 0 ;
-  fAbsIdList  = 0;
-  fParentsList = 0;
-  fTime = 0. ;
-  //  fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
-  fCoreRadius = 10;        //HG Check this
 
   AliRunLoader *rl = AliRunLoader::GetRunLoader();
-  fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
-  //fGeom = AliEMCALGeometry::GetInstance();
-  fGeom->GetTransformationForSM(); // Global <-> Local
+  if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
+    fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  else
+    fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
+  fGeomPtr->GetTransformationForSM(); // Global <-> Local
 }
 
 //____________________________________________________________________________
-AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) : AliRecPoint(opt)
+AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) 
+  : AliRecPoint(opt),
+    fGeomPtr(0),
+    fClusterType(-1),
+    fCoreEnergy(0),
+    fDispersion(0),
+    fEnergyList(0),
+    fTimeList(0),
+    fAbsIdList(0),
+    fTime(-1.),
+    fCoreRadius(10),  //HG check this
+    fDETracksList(0),
+    fMulParent(0),
+    fMaxParent(1000),
+    fParentsList(0),
+    fDEParentsList(0),
+    fSuperModuleNumber(0)
 {
   // ctor
-  fClusterType = -1;
   fMaxTrack = 1000 ;
-  fMaxParent = 1000;
   fMulDigit   = 0 ; 
-  fMulParent = 0; 
   fAmp   = 0. ;   
-  fCoreEnergy = 0 ; 
-  fEnergyList = 0 ;
-  fTimeList = 0 ;
-  fAbsIdList  = 0;
+  fDETracksList = new Float_t[fMaxTrack];
   fParentsList = new Int_t[fMaxParent];
-  fTime = -1. ;
-  //fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
-  fCoreRadius = 10;        //HG Check this
-  //fGeom = AliEMCALGeometry::GetInstance();
+  fDEParentsList = new Float_t[fMaxParent];
+  for (Int_t i = 0; i < fMaxTrack; i++)
+    fDETracksList[i] = 0;
+  for (Int_t i = 0; i < fMaxParent; i++)
+    fDEParentsList[i] = 0;
+
   AliRunLoader *rl = AliRunLoader::GetRunLoader();
-  fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
-  fGeom->GetTransformationForSM(); // Global <-> Local
+  if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
+    fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  else
+    fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
+  fGeomPtr->GetTransformationForSM(); // Global <-> Local
+}
+
+//____________________________________________________________________________
+AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp) 
+  : AliRecPoint(rp),
+    fGeomPtr(rp.fGeomPtr),
+    fClusterType(rp.fClusterType),
+    fCoreEnergy(rp.fCoreEnergy),
+    fDispersion(rp.fDispersion),
+    fEnergyList(0),
+    fTimeList(0),
+    fAbsIdList(0),
+    fTime(rp.fTime),
+    fCoreRadius(rp.fCoreRadius),
+    fDETracksList(0),
+    fMulParent(rp.fMulParent),
+    fMaxParent(rp.fMaxParent),
+    fParentsList(0),
+    fDEParentsList(0),
+    fSuperModuleNumber(rp.fSuperModuleNumber)
+{
+  //copy ctor
+  fLambda[0] = rp.fLambda[0];
+  fLambda[1] = rp.fLambda[1];
+
+  fEnergyList = new Float_t[rp.fMaxDigit];
+  fTimeList = new Float_t[rp.fMaxDigit];
+  fAbsIdList = new Int_t[rp.fMaxDigit];
+  for(Int_t i = 0; i < rp.fMulDigit; i++) {
+    fEnergyList[i] = rp.fEnergyList[i];
+    fTimeList[i] = rp.fTimeList[i];
+    fAbsIdList[i] = rp.fAbsIdList[i];
+  }
+  fDETracksList = new Float_t[rp.fMaxTrack];
+  for(Int_t i = 0; i < rp.fMulTrack; i++) fDETracksList[i] = rp.fDETracksList[i];
+  fParentsList = new Int_t[rp.fMaxParent];
+  for(Int_t i = 0; i < rp.fMulParent; i++) fParentsList[i] = rp.fParentsList[i];
+  fDEParentsList = new Float_t[rp.fMaxParent];
+  for(Int_t i = 0; i < rp.fMulParent; i++) fDEParentsList[i] = rp.fDEParentsList[i];
+
 }
 //____________________________________________________________________________
 AliEMCALRecPoint::~AliEMCALRecPoint()
@@ -102,8 +164,12 @@ AliEMCALRecPoint::~AliEMCALRecPoint()
     delete[] fTimeList ; 
   if ( fAbsIdList )
     delete[] fAbsIdList ; 
+   if ( fDETracksList)
+    delete[] fDETracksList;
    if ( fParentsList)
     delete[] fParentsList;
+   if ( fDEParentsList)
+    delete[] fDEParentsList;
 }
 
 //____________________________________________________________________________
@@ -118,7 +184,7 @@ void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
     fTimeList =  new Float_t[fMaxDigit]; 
   if(fAbsIdList == 0) {
     fAbsIdList =  new Int_t[fMaxDigit];
-    fSuperModuleNumber = fGeom->GetSuperModuleNumber(digit.GetId());
+    fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId());
   }
 
   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
@@ -176,18 +242,18 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d
   // A neighbour is defined as being two digits which share a corner
   
   static Bool_t areNeighbours = kFALSE ;
-  static Int_t nSupMod=0, nTower=0, nIphi=0, nIeta=0;
-  static int nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0;
+  static Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
+  static int nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
   static Int_t relid1[2] , relid2[2] ; // ieta, iphi
   static Int_t rowdiff=0, coldiff=0;
 
   areNeighbours = kFALSE ;
 
-  fGeom->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta);
-  fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, relid1[0],relid1[1]);
+  fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nModule,nIphi,nIeta);
+  fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
 
-  fGeom->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
-  fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, relid2[0],relid2[1]);
+  fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
+  fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
   
   rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
   coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
@@ -362,28 +428,53 @@ void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
 void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
 {
   // Calculates the dispersion of the shower at the origin of the RecPoint
+  // in cell units - Nov 16,2006
 
-  Double_t d = 0., wtot = 0., w = 0., xyzi[3], diff=0.;
-  Int_t iDigit=0, nstat=0, i=0;
+  Double_t d = 0., wtot = 0., w = 0.;
+  Int_t iDigit=0, nstat=0;
   AliEMCALDigit * digit ;
-
-  // Calculates the centre of gravity in the local EMCAL-module coordinates   
-  if (!fLocPos.Mag()) 
-    EvalLocalPosition(logWeight, digits) ;
   
-  // Calculates the dispersion in coordinates 
+  // Calculates the dispersion in cell units 
+  Double_t etai, phii, etaMean=0.0, phiMean=0.0; 
+  int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
+  int iphi=0, ieta=0;
+  // Calculate mean values
   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
 
-    fGeom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
-    w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+    if (fAmp>0 && fEnergyList[iDigit]>0) {
+      fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
+      fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
+      etai=(Double_t)ieta;
+      phii=(Double_t)iphi;
+      w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+
+      if(w>0.0) {
+        phiMean += phii*w;
+        etaMean += etai*w;
+        wtot    += w;
+      }
+    }
+  }
+  if (wtot>0) {
+    phiMean /= wtot ;
+    etaMean /= wtot ;
+  } else AliError(Form("Wrong weight %f\n", wtot));
 
-    if(w>0.0) {
-      wtot += w;
-      nstat++;
-      for(i=0; i<3; i++ ) {
-        diff = xyzi[i] - double(fLocPos[i]);
-        d   += w * diff*diff;
+  // Calculate dispersion
+  for(iDigit=0; iDigit < fMulDigit; iDigit++) {
+    digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
+
+    if (fAmp>0 && fEnergyList[iDigit]>0) {
+      fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
+      fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
+      etai=(Double_t)ieta;
+      phii=(Double_t)iphi;
+      w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+
+      if(w>0.0) {
+        nstat++;
+        d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
       }
     }
   }
@@ -407,11 +498,11 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit
   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
 
-    fGeom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
+    fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
     // printf(" Id %i : Local x,y,z %f %f %f \n", digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
 
-    if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
-    else                w = fEnergyList[iDigit]; // just energy
+    if(logWeight > 0.0)  w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
+    else  w = fEnergyList[iDigit]; // just energy
 
     if(w>0.0) {
       wtot += w ;
@@ -461,9 +552,10 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
   // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius
   // in accordance with shower profile the energy deposition 
   // should be less than 2%
+  // Unfinished - Nov 15,2006
+  // Distance is calculate in (phi,eta) units
 
   AliEMCALDigit * digit ;
-  const Float_t kDeg2Rad = TMath::DegToRad() ;
 
   Int_t iDigit;
 
@@ -471,16 +563,16 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
     EvalLocalPosition(logWeight, digits);
   }
   
+  Double_t phiPoint = fLocPos.Phi(), etaPoint = fLocPos.Eta();
+  Double_t eta, phi, distance;
   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
     
-    Float_t ieta = 0.0;
-    Float_t iphi = 0.0;
-    //fGeom->PosInAlice(digit->GetId(), ieta, iphi);
-    fGeom->EtaPhiFromIndex(digit->GetId(),ieta, iphi) ;
-    iphi = iphi * kDeg2Rad;
+    eta = phi = 0.0;
+    fGeomPtr->EtaPhiFromIndex(digit->GetId(),eta, phi) ;
+    phi = phi * TMath::DegToRad();
   
-    Float_t distance = TMath::Sqrt((ieta-fLocPos.X())*(ieta-fLocPos.X())+(iphi-fLocPos.Y())*(iphi-fLocPos.Y())) ;
+    distance = TMath::Sqrt((eta-etaPoint)*(eta-etaPoint)+(phi-phiPoint)*(phi-phiPoint));
     if(distance < fCoreRadius)
       fCoreEnergy += fEnergyList[iDigit] ;
   }
@@ -490,40 +582,44 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
 void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
 {
   // Calculates the axis of the shower ellipsoid in eta and phi
+  // in cell units
+
+  static TString gn(fGeomPtr->GetName());
 
-  Double_t wtot = 0. ;
+  Double_t wtot = 0.;
   Double_t x    = 0.;
   Double_t z    = 0.;
   Double_t dxx  = 0.;
   Double_t dzz  = 0.;
   Double_t dxz  = 0.;
 
-  const Float_t kDeg2Rad = TMath::DegToRad();
-  AliEMCALDigit * digit ;
+  AliEMCALDigit * digit = 0;
 
-  TString gn(fGeom->GetName());
-
-  Int_t iDigit;
-  
-  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+  Double_t etai , phii, w; 
+  int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
+  int iphi=0, ieta=0;
+  for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
-    Float_t etai = 0. ;
-    Float_t phii = 0. ; 
-    if(gn.Contains("SHISH")) { // have to be change - Feb 28, 2006
-    //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion
-    // for this geometry does not exist
-      int nSupMod=0, nTower=0, nIphi=0, nIeta=0;
-      int iphi=0, ieta=0;
-      fGeom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
-      fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
-      etai=(Float_t)ieta;
-      phii=(Float_t)iphi;
-    } else {
-      fGeom->EtaPhiFromIndex(digit->GetId(), etai, phii);
-      phii = phii * kDeg2Rad;
+    etai = phii = 0.; 
+    if(gn.Contains("SHISH")) { 
+    // Nov 15,2006 - use cell numbers as coordinates
+    // Copied for shish-kebab geometry, ieta,iphi is cast as double as eta,phi
+    // We can use the eta,phi(or coordinates) of cell
+      nSupMod = nModule = nIphi = nIeta = iphi = ieta = 0;
+
+      fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
+      fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
+      etai=(Double_t)ieta;
+      phii=(Double_t)iphi;
+    } else { // 
+      fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii);
+      phii = phii * TMath::DegToRad();
     }
 
-    Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+    w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+    // fAmp summed amplitude of digits, i.e. energy of recpoint 
+    // Gives smaller value of lambda than log weight  
+    // w = fEnergyList[iDigit] / fAmp; // Nov 16, 2006 - try just energy
 
     dxx  += w * etai * etai ;
     x    += w * etai ;
@@ -569,21 +665,19 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
 //______________________________________________________________________________
 void  AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
 {
-  // Constructs the list of primary particles (tracks) which have contributed to this RecPoint
+  // Constructs the list of primary particles (tracks) which 
+  // have contributed to this RecPoint and calculate deposited energy 
+  // for each track
   
   AliEMCALDigit * digit ;
-  Int_t * tempo    = new Int_t[fMaxTrack] ;
+  Int_t * primArray = new Int_t[fMaxTrack] ;
+  Float_t * dEPrimArray = new Float_t[fMaxTrack] ;
 
   Int_t index ;  
   for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
     digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ; 
     Int_t nprimaries = digit->GetNprimary() ;
     if ( nprimaries == 0 ) continue ;
-    Int_t * newprimaryarray = new Int_t[nprimaries] ;
-    Int_t ii ; 
-    for ( ii = 0 ; ii < nprimaries ; ii++)
-      newprimaryarray[ii] = digit->GetPrimary(ii+1) ; 
-
     Int_t jndex ;
     for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
       if ( fMulTrack > fMaxTrack ) {
@@ -591,29 +685,34 @@ void  AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
        Error("GetNprimaries", "increase fMaxTrack ")  ;
        break ;
       }
-      Int_t newprimary = newprimaryarray[jndex] ;
+      Int_t newPrimary = digit->GetPrimary(jndex+1);
+      Float_t dEPrimary = digit->GetDEPrimary(jndex+1);
       Int_t kndex ;
       Bool_t already = kFALSE ;
       for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored
-       if ( newprimary == tempo[kndex] ){
+       if ( newPrimary == primArray[kndex] ){
          already = kTRUE ;
+         dEPrimArray[kndex] += dEPrimary; 
          break ;
        }
       } // end of check
       if ( !already && (fMulTrack < fMaxTrack)) { // store it
-       tempo[fMulTrack] = newprimary ; 
+       primArray[fMulTrack] = newPrimary ; 
+       dEPrimArray[fMulTrack] = dEPrimary ; 
        fMulTrack++ ;
       } // store it
     } // all primaries in digit
-    delete [] newprimaryarray ; 
   } // all digits
 
-  
-  fTracksList = new Int_t[fMulTrack] ;
-  for(index = 0; index < fMulTrack; index++)
-   fTracksList[index] = tempo[index] ;
-  delete [] tempo ;
+  Int_t *sortIdx = new Int_t[fMulTrack];
+  TMath::Sort(fMulTrack,dEPrimArray,sortIdx); 
+  for(index = 0; index < fMulTrack; index++) {
+    fTracksList[index] = primArray[sortIdx[index]] ;    
+    fDETracksList[index] = dEPrimArray[sortIdx[index]] ;
+  }
+  delete [] sortIdx;
+  delete [] primArray ;
+  delete [] dEPrimArray ;
 
 }
 
@@ -623,17 +722,14 @@ void  AliEMCALRecPoint::EvalParents(TClonesArray * digits)
   // Constructs the list of parent particles (tracks) which have contributed to this RecPoint
  
   AliEMCALDigit * digit ;
-  Int_t * tempo    = new Int_t[fMaxParent] ;
+  Int_t * parentArray = new Int_t[fMaxTrack] ;
+  Float_t * dEParentArray = new Float_t[fMaxTrack] ;
 
   Int_t index ;  
   for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
     digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ; 
     Int_t nparents = digit->GetNiparent() ;
     if ( nparents == 0 ) continue ;
-    Int_t * newparentarray = new Int_t[nparents] ;
-    Int_t ii ; 
-    for ( ii = 0 ; ii < nparents ; ii++)
-      newparentarray[ii] = digit->GetIparent(ii+1) ; 
 
     Int_t jndex ;
     for ( jndex = 0 ; jndex < nparents ; jndex++ ) { // all primaries in digit
@@ -642,31 +738,37 @@ void  AliEMCALRecPoint::EvalParents(TClonesArray * digits)
        Error("GetNiparent", "increase fMaxParent")  ;
        break ;
       }
-      Int_t newparent = newparentarray[jndex] ;
+      Int_t newParent = digit->GetIparent(jndex+1) ;
+      Float_t newdEParent = digit->GetDEParent(jndex+1) ;
       Int_t kndex ;
       Bool_t already = kFALSE ;
       for ( kndex = 0 ; kndex < fMulParent ; kndex++ ) { //check if not already stored
-       if ( newparent == tempo[kndex] ){
+       if ( newParent == parentArray[kndex] ){
+         dEParentArray[kndex] += newdEParent;
          already = kTRUE ;
          break ;
        }
       } // end of check
       if ( !already && (fMulTrack < fMaxTrack)) { // store it
-       tempo[fMulParent] = newparent ; 
+       parentArray[fMulParent] = newParent ; 
+       dEParentArray[fMulParent] = newdEParent ; 
        fMulParent++ ;
       } // store it
     } // all parents in digit
-    delete [] newparentarray ; 
   } // all digits
 
   if (fMulParent>0) {
-    fParentsList = new Int_t[fMulParent] ;
-    for(index = 0; index < fMulParent; index++)
-      fParentsList[index] = tempo[index] ;
+    Int_t *sortIdx = new Int_t[fMulParent];
+    TMath::Sort(fMulParent,dEParentArray,sortIdx); 
+    for(index = 0; index < fMulParent; index++) {
+      fParentsList[index] = parentArray[sortIdx[index]] ;      
+      fDEParentsList[index] = dEParentArray[sortIdx[index]] ;
+    }
+    delete [] sortIdx;
   }
  
-  delete [] tempo ;
-
+  delete [] parentArray;
+  delete [] dEParentArray;
 }
 
 //____________________________________________________________________________
@@ -685,7 +787,7 @@ void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
   // returns the position of the cluster in the global reference system of ALICE
   // These are now the Cartesian X, Y and Z
   //  cout<<" geom "<<geom<<endl;
-  fGeom->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
+  fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
 }
 
 //____________________________________________________________________________
@@ -776,59 +878,11 @@ Int_t  AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit **  maxAt, Float_t *
 Int_t AliEMCALRecPoint::GetPrimaryIndex() const  
 {
   // Get the primary track index in TreeK which deposits the most energy 
-  // in Digits which forms RecPoint. Kinematics, Hits and Digits must be 
-  // loaded before the call of the method.
-
-  AliRunLoader *rl = AliRunLoader::GetRunLoader(); 
-  if (!rl) 
-    AliError(Form(" No Runloader ")) ; 
-  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
-    (rl->GetDetectorLoader("EMCAL"));
-
-  // Get the list of digits forming this RecPoint
-  Int_t  nDigits   = fMulDigit   ;
-  Int_t *digitList = fDigitsList ;
-  
-  // Find the digit with maximum amplitude
-  AliEMCALDigit *digit = 0;
-  TClonesArray *digits = emcalLoader->Digits();
-  Int_t maxAmp = 0;
-  Int_t bestDigitIndex = -1;
-  for (Int_t iDigit=0; iDigit<nDigits; iDigit++) {
-    digit = static_cast<AliEMCALDigit *>(digits->At(digitList[iDigit]));
-    if (digit->GetAmp() > maxAmp) {
-      maxAmp = digit->GetAmp();
-      bestDigitIndex = iDigit;
-    }
-  }
+  // in Digits which forms RecPoint. 
 
-  digit = static_cast<AliEMCALDigit *>(digits->At(digitList[bestDigitIndex]));  
-
-  // Get the list of hits producing this digit,
-  // find which hit has deposited more energy 
-  // and find the primary track.
-
-  AliEMCALHit *hit = 0;
-  TClonesArray *hits = emcalLoader->Hits();
-
-  Double_t maxedep  =  0;
-  Int_t    maxtrack = -1;
-  Int_t    nHits    = hits ->GetEntries();
-  Int_t    id       = digit->GetId();
-  for (Int_t iHit=0; iHit<nHits; iHit++) {
-    hit = static_cast<AliEMCALHit*> (hits->At(iHit)) ;
-    if(hit->GetId() == id){
-      Double_t edep  = hit->GetEnergy();
-      Int_t    track = hit->GetIparent();//Primary();
-      if(edep > maxedep){
-       maxedep  = edep;
-       maxtrack = track;
-      }
-    }
-  }
-  if (maxtrack != -1) return maxtrack; 
-  return -12345;                       // no track found :(
+  if (fMulTrack)
+    return fTracksList[0];
+  return -12345;
 }
 
 //____________________________________________________________________________
@@ -923,3 +977,11 @@ void AliEMCALRecPoint::Print(Option_t *) const
   message += "       Stored at position %d" ; 
   Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ;  
 }
+
+Double_t  AliEMCALRecPoint::GetPointEnergy() const
+{
+  static double e;
+  e=0.0;
+  for(int ic=0; ic<GetMultiplicity(); ic++) e += double(fEnergyList[ic]);
+  return e;
+}