In case of recalibrating all cells, check not only the energy recalibration flag...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx
index ce69ba8..db11d25 100644 (file)
@@ -22,8 +22,9 @@
 //
 //
 // Author:  Gustavo Conesa (LPSC- Grenoble) 
-///////////////////////////////////////////////////////////////////////////////
+//          Track matching part: Rongrong Ma (Yale)
 
+///////////////////////////////////////////////////////////////////////////////
 // --- standard c ---
 
 // standard C++ includes
 #include <TGeoManager.h>
 #include <TGeoMatrix.h>
 #include <TGeoBBox.h>
+#include <TH2F.h>
+#include <TArrayI.h>
+#include <TArrayF.h>
+#include <TObjArray.h>
 
 // STEER includes
-#include "AliEMCALRecoUtils.h"
-#include "AliEMCALGeometry.h"
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
+#include "AliVEvent.h"
 #include "AliLog.h"
+#include "AliPID.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliESDtrack.h"
+#include "AliAODTrack.h"
+#include "AliExternalTrackParam.h"
+#include "AliESDfriendTrack.h"
+#include "AliTrackerBase.h"
+
+// EMCAL includes
+#include "AliEMCALRecoUtils.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALTrack.h"
+#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
+#include "AliEMCALPIDUtils.h"
 
 ClassImp(AliEMCALRecoUtils)
   
 //______________________________________________
 AliEMCALRecoUtils::AliEMCALRecoUtils():
-  fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
-  fPosAlgo(kUnchanged),fW0(4.),
-  fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
-  fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(),
-  fNCellsFromEMCALBorder(0),fNoEMCALBorderAtEta0(kFALSE)
+  fParticleType(kPhoton),                 fPosAlgo(kUnchanged),                   fW0(4.), 
+  fNonLinearityFunction(kNoCorrection),   fNonLinearThreshold(30),
+  fSmearClusterEnergy(kFALSE),            fRandom(),
+  fCellsRecalibrated(kFALSE),             fRecalibration(kFALSE),                 fEMCALRecalibrationFactors(),
+  fTimeRecalibration(kFALSE),             fEMCALTimeRecalibrationFactors(),
+  fUseRunCorrectionFactors(kFALSE),       fRunCorrectionFactorsSet(kFALSE),
+  fRemoveBadChannels(kFALSE),             fRecalDistToBadChannels(kFALSE),        fEMCALBadChannelMap(),
+  fNCellsFromEMCALBorder(0),              fNoEMCALBorderAtEta0(kTRUE),
+  fRejectExoticCluster(kFALSE),           fPIDUtils(),                            fAODFilterMask(32),
+  fMatchedTrackIndex(0x0),                fMatchedClusterIndex(0x0), 
+  fResidualEta(0x0), fResidualPhi(0x0),   fCutEtaPhiSum(kTRUE),                   fCutEtaPhiSeparate(kFALSE), 
+  fCutR(0.1),                             fCutEta(0.025),                         fCutPhi(0.05), 
+  fMass(0.139),                           fStep(10),
+  fTrackCutsType(kLooseCut),              fCutMinTrackPt(0),                      fCutMinNClusterTPC(-1), 
+  fCutMinNClusterITS(-1),                 fCutMaxChi2PerClusterTPC(1e10),         fCutMaxChi2PerClusterITS(1e10),
+  fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
+  fCutMaxDCAToVertexXY(1e10),             fCutMaxDCAToVertexZ(1e10),              fCutDCAToVertex2D(kFALSE)
 {
 //
   // Constructor.
@@ -57,34 +88,81 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   // during Reco algorithm execution
   //
   
+  //Misalignment matrices
   for(Int_t i = 0; i < 15 ; i++) {
       fMisalTransShift[i] = 0.; 
-      fMisalRotShift[i] = 0.; 
+      fMisalRotShift[i]   = 0.; 
   }
-  for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = 0.; 
-  //For kPi0GammaGamma case, but default is no correction
-  fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
-  fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
-  fNonLinearityParams[2] = 1.046;
   
+  //Non linearity
+  for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = 0.; 
+
+  //For kBeamTestCorrected case, but default is no correction
+  fNonLinearityParams[0] =  0.99078;
+  fNonLinearityParams[1] =  0.161499;
+  fNonLinearityParams[2] =  0.655166; 
+  fNonLinearityParams[3] =  0.134101;
+  fNonLinearityParams[4] =  163.282;
+  fNonLinearityParams[5] =  23.6904;
+  fNonLinearityParams[6] =  0.978;
+  
+  //For kPi0GammaGamma case
+  //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
+  //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
+  //fNonLinearityParams[2] = 1.046;
+
+  //Cluster energy smearing
+  fSmearClusterEnergy   = kFALSE;
+  fSmearClusterParam[0] = 0.07; // * sqrt E term
+  fSmearClusterParam[1] = 0.00; // * E term
+  fSmearClusterParam[2] = 0.00; // constant
+  
+  //Track matching
+  fMatchedTrackIndex     = new TArrayI();
+  fMatchedClusterIndex   = new TArrayI();
+  fResidualPhi           = new TArrayF();
+  fResidualEta           = new TArrayF();
+  fPIDUtils              = new AliEMCALPIDUtils();
+
+  InitTrackCuts();
 }
 
 //______________________________________________________________________
 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) 
-: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction), 
-  fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), 
-  fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
-  fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
-  fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0)
-
+: TNamed(reco), 
+  fParticleType(reco.fParticleType),                         fPosAlgo(reco.fPosAlgo),     fW0(reco.fW0),
+  fNonLinearityFunction(reco.fNonLinearityFunction),         fNonLinearThreshold(reco.fNonLinearThreshold),
+  fSmearClusterEnergy(reco.fSmearClusterEnergy),             fRandom(),
+  fCellsRecalibrated(reco.fCellsRecalibrated),
+  fRecalibration(reco.fRecalibration),                       fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
+  fTimeRecalibration(reco.fTimeRecalibration),               fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
+  fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),   fRunCorrectionFactorsSet(reco.fRunCorrectionFactorsSet),
+  fRemoveBadChannels(reco.fRemoveBadChannels),               fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
+  fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
+  fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
+  fRejectExoticCluster(reco.fRejectExoticCluster),           fPIDUtils(reco.fPIDUtils), 
+  fAODFilterMask(reco.fAODFilterMask),
+  fMatchedTrackIndex(  reco.fMatchedTrackIndex?  new TArrayI(*reco.fMatchedTrackIndex):0x0),
+  fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
+  fResidualEta(        reco.fResidualEta?        new TArrayF(*reco.fResidualEta):0x0),
+  fResidualPhi(        reco.fResidualPhi?        new TArrayF(*reco.fResidualPhi):0x0),
+  fCutEtaPhiSum(reco.fCutEtaPhiSum),                         fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate), 
+  fCutR(reco.fCutR),        fCutEta(reco.fCutEta),           fCutPhi(reco.fCutPhi),
+  fMass(reco.fMass),        fStep(reco.fStep),
+  fTrackCutsType(reco.fTrackCutsType),                       fCutMinTrackPt(reco.fCutMinTrackPt), 
+  fCutMinNClusterTPC(reco.fCutMinNClusterTPC),               fCutMinNClusterITS(reco.fCutMinNClusterITS), 
+  fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC),   fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
+  fCutRequireTPCRefit(reco.fCutRequireTPCRefit),             fCutRequireITSRefit(reco.fCutRequireITSRefit),
+  fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),     fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),    
+  fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),             fCutDCAToVertex2D(reco.fCutDCAToVertex2D)
 {
   //Copy ctor
   
-  for(Int_t i = 0; i < 15 ; i++) {
-      fMisalRotShift[i] = reco.fMisalRotShift[i]; 
-      fMisalTransShift[i] = reco.fMisalTransShift[i]; 
-  } 
-  for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
+  for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i]      = reco.fMisalRotShift[i]      ; 
+                                   fMisalTransShift[i]    = reco.fMisalTransShift[i]    ; } 
+  for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
+  for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }
+
 }
 
 
@@ -96,20 +174,111 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   if(this == &reco)return *this;
   ((TNamed *)this)->operator=(reco);
 
-  fNonLinearityFunction  = reco.fNonLinearityFunction;
-  fParticleType          = reco.fParticleType;
-  fPosAlgo               = reco.fPosAlgo; 
-  fW0                    = reco.fW0;
-  fRecalibration         = reco.fRecalibration;
+  for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i]    = reco.fMisalTransShift[i]    ; 
+                                   fMisalRotShift[i]      = reco.fMisalRotShift[i]      ; }
+  for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
+  for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }   
+  
+  fParticleType              = reco.fParticleType;
+  fPosAlgo                   = reco.fPosAlgo; 
+  fW0                        = reco.fW0;
+  
+  fNonLinearityFunction      = reco.fNonLinearityFunction;
+  fNonLinearThreshold        = reco.fNonLinearThreshold;
+  fSmearClusterEnergy        = reco.fSmearClusterEnergy;
+
+  fCellsRecalibrated         = reco.fCellsRecalibrated;
+  fRecalibration             = reco.fRecalibration;
   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
-  fRemoveBadChannels     = reco.fRemoveBadChannels;
-  fEMCALBadChannelMap    = reco.fEMCALBadChannelMap;
-  fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
-  fNoEMCALBorderAtEta0   = reco.fNoEMCALBorderAtEta0;
+
+  fTimeRecalibration             = reco.fTimeRecalibration;
+  fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
+
+  fUseRunCorrectionFactors   = reco.fUseRunCorrectionFactors;
+  fRunCorrectionFactorsSet   = reco.fRunCorrectionFactorsSet;
+  
+  fRemoveBadChannels         = reco.fRemoveBadChannels;
+  fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
+  fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
+  
+  fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
+  fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
+  fRejectExoticCluster       = reco.fRejectExoticCluster;           
+
+  fPIDUtils                  = reco.fPIDUtils;
+
+  fAODFilterMask             = reco.fAODFilterMask;
+  
+  fCutEtaPhiSum              = reco.fCutEtaPhiSum;
+  fCutEtaPhiSeparate         = reco.fCutEtaPhiSeparate;
+  fCutR                      = reco.fCutR;
+  fCutEta                    = reco.fCutEta;
+  fCutPhi                    = reco.fCutPhi;
+  fMass                      = reco.fMass;
+  fStep                      = reco.fStep;
+  fRejectExoticCluster       = reco.fRejectExoticCluster;
+
+  fTrackCutsType             = reco.fTrackCutsType;
+  fCutMinTrackPt             = reco.fCutMinTrackPt;
+  fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
+  fCutMinNClusterITS         = reco.fCutMinNClusterITS; 
+  fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
+  fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
+  fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
+  fCutRequireITSRefit        = reco.fCutRequireITSRefit;
+  fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
+  fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
+  fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
+  fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
+  
+  if(reco.fResidualEta){
+    // assign or copy construct
+    if(fResidualEta){ 
+      *fResidualEta = *reco.fResidualEta;
+    }
+    else fResidualEta = new TArrayF(*reco.fResidualEta);
+  }
+  else{
+    if(fResidualEta)delete fResidualEta;
+    fResidualEta = 0;
+  }
+  
+  if(reco.fResidualPhi){
+    // assign or copy construct
+    if(fResidualPhi){ 
+      *fResidualPhi = *reco.fResidualPhi;
+    }
+    else fResidualPhi = new TArrayF(*reco.fResidualPhi);
+  }
+  else{
+    if(fResidualPhi)delete fResidualPhi;
+    fResidualPhi = 0;
+  }
   
-  for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
-  for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
+  if(reco.fMatchedTrackIndex){
+    // assign or copy construct
+    if(fMatchedTrackIndex){ 
+      *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
+    }
+    else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
+  }
+  else{
+    if(fMatchedTrackIndex)delete fMatchedTrackIndex;
+    fMatchedTrackIndex = 0;
+  }  
   
+  if(reco.fMatchedClusterIndex){
+    // assign or copy construct
+    if(fMatchedClusterIndex){ 
+      *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
+    }
+    else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
+  }
+  else{
+    if(fMatchedClusterIndex)delete fMatchedClusterIndex;
+    fMatchedClusterIndex = 0;
+  }
+   
   return *this;
 }
 
@@ -119,16 +288,28 @@ AliEMCALRecoUtils::~AliEMCALRecoUtils()
 {
   //Destructor.
        
-       if(fEMCALRecalibrationFactors) { 
-               fEMCALRecalibrationFactors->Clear();
-               delete  fEMCALRecalibrationFactors;
+  if(fEMCALRecalibrationFactors) { 
+    fEMCALRecalibrationFactors->Clear();
+    delete  fEMCALRecalibrationFactors;
+  }    
+  
+  if(fEMCALTimeRecalibrationFactors) { 
+               fEMCALTimeRecalibrationFactors->Clear();
+               delete  fEMCALTimeRecalibrationFactors;
        }       
   
   if(fEMCALBadChannelMap) { 
-               fEMCALBadChannelMap->Clear();
-               delete  fEMCALBadChannelMap;
-       }
-  
+    fEMCALBadChannelMap->Clear();
+    delete  fEMCALBadChannelMap;
+  }
+  delete fMatchedTrackIndex   ; 
+  delete fMatchedClusterIndex ; 
+  delete fResidualEta         ; 
+  delete fResidualPhi         ; 
+  delete fPIDUtils            ;
+
+  InitTrackCuts();
 }
 
 //_______________________________________________________________
@@ -137,14 +318,20 @@ Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCl
        // Given the list of AbsId of the cluster, get the maximum cell and 
        // check if there are fNCellsFromBorder from the calorimeter border
        
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return kFALSE;
+  }
+  
   //If the distance to the border is 0 or negative just exit accept all clusters
        if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
   
-  Int_t absIdMax       = -1, iSM =-1, ieta = -1, iphi = -1;  
-  GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi);
+  Int_t absIdMax       = -1, iSM =-1, ieta = -1, iphi = -1;
+  Bool_t shared = kFALSE;
+  GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
 
-  AliDebug(2,Form("AliEMCALRecoUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n", 
-           absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E()));
+  AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", 
+           absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
        
        if(absIdMax==-1) return kFALSE;
        
@@ -178,17 +365,24 @@ Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCl
     }
   }//eta 0 not checked
     
-  AliDebug(2,Form("AliEMCALRecoUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d",
+  AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
        
-       if (okcol && okrow) return kTRUE; 
-       else                return kFALSE;
+       if (okcol && okrow) {
+    //printf("Accept\n");
+    return kTRUE;
+  }
+       else  {
+    //printf("Reject\n");
+    AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
+    return kFALSE;
+  }
        
 }      
 
 
 //_________________________________________________________________________________________________________
-Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
+Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, const Int_t nCells){
        // Check that in the cluster cells, there is no bad channel of those stored 
        // in fEMCALBadChannelMap or fPHOSBadChannelMap
        
@@ -205,7 +399,10 @@ Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, USho
     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                     
-    if(GetEMCALChannelStatus(imod, icol, irow))return kTRUE;
+    if(GetEMCALChannelStatus(imod, icol, irow)){
+      AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
+      return kTRUE;
+    }
                
        }// cell cluster loop
        
@@ -213,43 +410,123 @@ Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, USho
        
 }
 
+//_________________________________________________
+Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
+  // Check if the cluster has high energy  but small number of cells
+  // The criteria comes from Gustavo's study
+  //
+  
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return kFALSE;
+  }
+  
+  Int_t nc = cluster->GetNCells() ;
+  
+  if      ( nc > 8 )                   return kFALSE ; // Good cluster, needed for 3x3 clusterizer  
+  else if ( nc < 1 + cluster->E()/3. ) return kTRUE  ; // Bad cluster
+  else                                 return kFALSE ; // Good cluster
+  
+}
+
+//__________________________________________________
+Float_t AliEMCALRecoUtils::SmearClusterEnergy(AliVCluster* cluster) {
+
+  //In case of MC analysis, smear energy to match resolution/calibration in real data
+  
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return 0;
+  }
+  
+  Float_t energy    = cluster->E() ;
+  Float_t rdmEnergy = energy ;
+  if(fSmearClusterEnergy){
+    rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
+                                    fSmearClusterParam[1] * energy +
+                                    fSmearClusterParam[2] );
+    AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
+  }
+  
+  return rdmEnergy ;
+
+}
+
 //__________________________________________________
 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
 // Correct cluster energy from non linearity functions
+  
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return 0;
+  }
+  
   Float_t energy = cluster->E();
   
   switch (fNonLinearityFunction) {
       
     case kPi0MC:
+    {
       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
-      //Double_t par0 = 1.001;
-      //Double_t par1 = -0.01264;
-      //Double_t par2 = -0.03632;
-      //Double_t par3 = 0.1798;
-      //Double_t par4 = -0.522;
-       energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
+      //Double_t fNonLinearityParams[0] = 1.014;
+      //Double_t fNonLinearityParams[1] = -0.03329;
+      //Double_t fNonLinearityParams[2] = -0.3853;
+      //Double_t fNonLinearityParams[3] = 0.5423;
+      //Double_t fNonLinearityParams[4] = -0.4335;
+       energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
       break;
-      
+    }
+     
     case kPi0GammaGamma:
-
+    {
       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
-      //Double_t par0 = 0.1457;
-      //Double_t par1 = -0.02024;
-      //Double_t par2 = 1.046;
+      //Double_t fNonLinearityParams[0] = 1.04;
+      //Double_t fNonLinearityParams[1] = -0.1445;
+      //Double_t fNonLinearityParams[2] = 1.046;
       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
       break;
+    }
       
     case kPi0GammaConversion:
-      
+    {
       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
-      //Double_t C = 0.139393/0.1349766;
-      //Double_t a = 0.0566186;
-      //Double_t b = 0.982133;
+      //fNonLinearityParams[0] = 0.139393/0.1349766;
+      //fNonLinearityParams[1] = 0.0566186;
+      //fNonLinearityParams[2] = 0.982133;
       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
       
       break;
+    }
+      
+    case kBeamTest:
+    {
+      //From beam test, Alexei's results, for different ZS thresholds
+      //                        th=30 MeV; th = 45 MeV; th = 75 MeV
+      //fNonLinearityParams[0] = 1.007;      1.003;      1.002 
+      //fNonLinearityParams[1] = 0.894;      0.719;      0.797 
+      //fNonLinearityParams[2] = 0.246;      0.334;      0.358 
+      //Rescale the param[0] with 1.03
+      energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
+      
+      break;
+    }
+      
+    case kBeamTestCorrected:
+    {
+      //From beam test, corrected for material between beam and EMCAL
+      //fNonLinearityParams[0] =  0.99078
+      //fNonLinearityParams[1] =  0.161499;
+      //fNonLinearityParams[2] =  0.655166; 
+      //fNonLinearityParams[3] =  0.134101;
+      //fNonLinearityParams[4] =  163.282;
+      //fNonLinearityParams[5] =  23.6904;
+      //fNonLinearityParams[6] =  0.978;
+        energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
+
+      break;
+    }
       
     case kNoCorrection:
       AliDebug(2,"No correction on the energy\n");
@@ -260,6 +537,67 @@ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
   return energy;
 
 }
+//__________________________________________________
+void AliEMCALRecoUtils::InitNonLinearityParam()
+{
+       //Initialising Non Linearity Parameters
+       
+       if(fNonLinearityFunction == kPi0MC)
+               {
+                       fNonLinearityParams[0] = 1.014;
+      fNonLinearityParams[1] = -0.03329;
+      fNonLinearityParams[2] = -0.3853;
+      fNonLinearityParams[3] = 0.5423;
+      fNonLinearityParams[4] = -0.4335;
+               }
+
+       if(fNonLinearityFunction == kPi0GammaGamma)
+               {
+                       fNonLinearityParams[0] = 1.04;
+                       fNonLinearityParams[1] = -0.1445;
+                       fNonLinearityParams[2] = 1.046;
+               }       
+
+       if(fNonLinearityFunction == kPi0GammaConversion)
+               {
+      fNonLinearityParams[0] = 0.139393;
+      fNonLinearityParams[1] = 0.0566186;
+      fNonLinearityParams[2] = 0.982133;
+               }       
+
+       if(fNonLinearityFunction == kBeamTest)
+               {
+                       if(fNonLinearThreshold == 30)
+                       {
+                               fNonLinearityParams[0] = 1.007; 
+                               fNonLinearityParams[1] = 0.894; 
+                               fNonLinearityParams[2] = 0.246; 
+                       }
+                       if(fNonLinearThreshold == 45)
+                       {
+                               fNonLinearityParams[0] = 1.003; 
+                               fNonLinearityParams[1] = 0.719; 
+                               fNonLinearityParams[2] = 0.334; 
+                       }
+                       if(fNonLinearThreshold == 75)
+                       {
+                               fNonLinearityParams[0] = 1.002; 
+                               fNonLinearityParams[1] = 0.797; 
+                               fNonLinearityParams[2] = 0.358; 
+                       }
+               }
+
+       if(fNonLinearityFunction == kBeamTestCorrected)
+               {
+                       fNonLinearityParams[0] =  0.99078;
+                       fNonLinearityParams[1] =  0.161499;
+                       fNonLinearityParams[2] =  0.655166; 
+                       fNonLinearityParams[3] =  0.134101;
+                       fNonLinearityParams[4] =  163.282;
+                       fNonLinearityParams[5] =  23.6904;
+                       fNonLinearityParams[6] =  0.978;
+               }
+}
 
 //__________________________________________________
 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const 
@@ -267,7 +605,7 @@ Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle
   //Calculate shower depth for a given cluster energy and particle type
 
   // parameters 
-  Float_t x0    = 1.23;
+  Float_t x0    = 1.31;
   Float_t ecr   = 8;
   Float_t depth = 0;
   
@@ -311,7 +649,8 @@ Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle
 }
 
 //__________________________________________________
-void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
+void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, 
+                                         Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
 {
   //For a given CaloCluster gets the absId of the cell 
   //with maximum energy deposit.
@@ -325,18 +664,31 @@ void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells*
   Int_t iTower  = -1;
   Int_t iIphi   = -1;
   Int_t iIeta   = -1;
-       
+  Int_t iSupMod0= -1;
+
+  if(!clu){
+    AliInfo("Cluster pointer null!");
+    absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
+    return;
+  }
+  
   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
     cellAbsId = clu->GetCellAbsId(iDig);
     fraction  = clu->GetCellAmplitudeFraction(iDig);
+    //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
+    if(iDig==0) iSupMod0=iSupMod;
+    else if(iSupMod0!=iSupMod) {
+      shared = kTRUE;
+      //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
+    }
     if(IsRecalibrationOn()) {
-      geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
-      geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
     }
     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
-    
+    //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
     if(eCell > eMax)  { 
       eMax  = eCell; 
       absId = cellAbsId;
@@ -348,7 +700,9 @@ void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells*
   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
   //Gives SuperModule and Tower numbers
   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                         iIphi, iIeta,iphi,ieta);    
+                                         iIphi, iIeta,iphi,ieta); 
+  //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
+  //printf("Max end---\n");
   
 }
 
@@ -360,10 +714,10 @@ void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
        Bool_t oldStatus = TH1::AddDirectoryStatus();
        TH1::AddDirectory(kFALSE);
   
-       fEMCALRecalibrationFactors = new TObjArray(12);
-       for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
+       fEMCALRecalibrationFactors = new TObjArray(10);
+       for (int i = 0; i < 10; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
        //Init the histograms with 1
-       for (Int_t sm = 0; sm < 12; sm++) {
+       for (Int_t sm = 0; sm < 10; sm++) {
                for (Int_t i = 0; i < 48; i++) {
                        for (Int_t j = 0; j < 24; j++) {
                                SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
@@ -377,6 +731,31 @@ void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
        TH1::AddDirectory(oldStatus);           
 }
 
+//________________________________________________________________
+void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors(){
+       //Init EMCAL recalibration factors
+       AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
+       //In order to avoid rewriting the same histograms
+       Bool_t oldStatus = TH1::AddDirectoryStatus();
+       TH1::AddDirectory(kFALSE);
+  
+       fEMCALTimeRecalibrationFactors = new TObjArray(4);
+       for (int i = 0; i < 4; i++) 
+    fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
+                                                 Form("hAllTimeAvBC%d",i),  
+                                                 48*24*10,0.,48*24*10)          );
+       //Init the histograms with 1
+       for (Int_t bc = 0; bc < 4; bc++) {
+                       for (Int_t i = 0; i < 48*24*10; i++) 
+        SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
+  }
+
+       fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
+       fEMCALTimeRecalibrationFactors->Compress();
+       
+       //In order to avoid rewriting the same histograms
+       TH1::AddDirectory(oldStatus);           
+}
 
 //________________________________________________________________
 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
@@ -386,14 +765,12 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
        Bool_t oldStatus = TH1::AddDirectoryStatus();
        TH1::AddDirectory(kFALSE);
        
-       fEMCALBadChannelMap = new TObjArray(12);
+       fEMCALBadChannelMap = new TObjArray(10);
        //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
-       for (int i = 0; i < 12; i++) {
+       for (int i = 0; i < 10; i++) {
                fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
        }
-       
-       //delete hTemp;
-       
+               
        fEMCALBadChannelMap->SetOwner(kTRUE);
        fEMCALBadChannelMap->Compress();
        
@@ -402,9 +779,16 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 }
 
 //________________________________________________________________
-void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
-       // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
+void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells, const Int_t bc){
+       // Recalibrate the cluster energy and Time, considering the recalibration map 
+  // and the energy of the cells and time that compose the cluster.
+  // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
        
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return;
+  }  
+  
        //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
        UShort_t * index    = cluster->GetCellsAbsId() ;
        Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
@@ -412,39 +796,170 @@ void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVClu
        
        //Initialize some used variables
        Float_t energy = 0;
-       Int_t absId    = -1;
-  Int_t icol = -1, irow = -1, imod=1;
+       Int_t   absId  =-1;
+  Int_t   icol   =-1, irow =-1, imod=1;
        Float_t factor = 1, frac = 0;
-       
+  Int_t   absIdMax = -1;
+  Float_t emax     = 0;
+  
        //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
        for(Int_t icell = 0; icell < ncells; icell++){
                absId = index[icell];
                frac =  fraction[icell];
                if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
-               Int_t iTower = -1, iIphi = -1, iIeta = -1; 
-               geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
-               if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
-               geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                  
-               factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
-    AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
-             imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
-               
+    
+    if(!fCellsRecalibrated && IsRecalibrationOn()){
+      
+      // Energy  
+      Int_t iTower = -1, iIphi = -1, iIeta = -1; 
+      geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
+      if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
+      geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                   
+      factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
+      
+      AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
+                      imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
+      
+    } 
+    
                energy += cells->GetCellAmplitude(absId)*factor*frac;
+    
+    if(emax < cells->GetCellAmplitude(absId)*factor*frac){
+      emax     = cells->GetCellAmplitude(absId)*factor*frac;
+      absIdMax = absId;
+    }
+
        }
        
-       
-               AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
-       
-       cluster->SetE(energy);
-       
+  cluster->SetE(energy);
+
+  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
+
+       // Recalculate time of cluster only for ESDs
+  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))){
+    
+    // Time
+    Double_t weightedTime = 0;
+    Double_t weight       = 0;
+    Double_t weightTot    = 0;
+    Double_t maxcellTime  = 0;
+    for(Int_t icell = 0; icell < ncells; icell++){
+      absId = index[icell];
+      frac =  fraction[icell];
+      if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
+
+      Double_t celltime = cells->GetCellTime(absId);
+      RecalibrateCellTime(absId, bc, celltime);
+      if(absId == absIdMax) maxcellTime = celltime;
+
+      if(!fCellsRecalibrated){
+      
+        Int_t iTower = -1, iIphi = -1, iIeta = -1; 
+        geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
+        if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
+        geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                 
+        factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
+        
+        AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
+                        imod,icol,irow,frac,factor,cells->GetCellTime(absId)));
+        
+      } 
+      
+      weight        = GetCellWeight(cells->GetCellAmplitude(absId)*factor*frac , energy );
+      weightTot    += weight;
+      weightedTime += celltime * weight;
+      
+    }
+    
+    if(weightTot > 0)
+      cluster->SetTOF(weightedTime/weightTot);
+    else 
+      cluster->SetTOF(maxcellTime);
+    
+  }
 }
 
+//________________________________________________________________
+void AliEMCALRecoUtils::RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells * cells, Int_t bc){
+       // Recalibrate the cells time and energy, considering the recalibration map and the energy 
+  // of the cells that compose the cluster.
+  // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
+
+  if(!IsRecalibrationOn() && !IsTimeRecalibrationOn()) return;
+  
+  if(!cells){
+    AliInfo("Cells pointer null!");
+    return;
+  }  
+  
+  fCellsRecalibrated = kTRUE;
+  
+  Int_t absId  =-1;
+  Int_t icol   =-1, irow  =-1, imod  = 1;
+  Int_t iTower =-1, iIeta =-1, iIphi =-1;
 
+  Int_t nEMcell = cells->GetNumberOfCells() ;
+  
+  for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
+    
+    absId = cells->GetCellNumber(iCell);
+    
+    // Energy
+    Float_t factor = 1;
+    if(IsRecalibrationOn()){
+      geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
+      if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
+      geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);   
+      factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
+               }
+    
+    Float_t cellE      = cells->GetAmplitude(iCell) * factor ;
+    
+    //Time
+    Double_t celltime = cells->GetCellTime(absId);
+    RecalibrateCellTime(absId, bc, celltime);
+    
+    //Set new values
+    cells->SetCell(iCell,cells->GetCellNumber(iCell),cellE, celltime);
+    
+  }
+  
+}
+
+//________________________________________________________________
+void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime){
+       // Recalibrate time of cell with absID  considering the recalibration map 
+  // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
+    
+  if(!fCellsRecalibrated && IsTimeRecalibrationOn()){
+//    printf("cell time org %g, ",celltime);
+
+    Double_t timeBCoffset = 0.;
+    if( bc%4 ==0 || bc%4==1) timeBCoffset = 100.*1.e-9; //in ns        
+    
+    Double_t celloffset = GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; 
+    
+//    printf("absId %d, time %f bc %d-%d: bc0 %f, bc1 %f, bc2 %f, bc3 %f \n", absId, celltime*1.e9,bc, bc%4, 
+//           GetEMCALChannelTimeRecalibrationFactor(0,absId),GetEMCALChannelTimeRecalibrationFactor(1,absId),
+//           GetEMCALChannelTimeRecalibrationFactor(2,absId),GetEMCALChannelTimeRecalibrationFactor(3,absId));
+    
+    celltime -= timeBCoffset ;
+    celltime -= celloffset   ;  
+//    printf("new %g\n",celltime);
+  }
+  
+}
+  
 //__________________________________________________
 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 {
   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
   
+  if(!clu){
+    AliInfo("Cluster pointer null!");
+    return;
+  }
+    
   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
@@ -467,29 +982,40 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet
   Float_t  weight = 0.,  totalWeight=0.;
   Float_t  newPos[3] = {0,0,0};
   Double_t pLocal[3], pGlobal[3];
-  
+  Bool_t shared = kFALSE;
+
   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
-  
-  GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
+  GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
   
+  //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
+  
   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
-    absId = clu->GetCellAbsId(iDig);
-    fraction  = clu->GetCellAmplitudeFraction(iDig);
-    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-    geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
-    geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                      
     
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
+      absId = clu->GetCellAbsId(iDig);
+      fraction  = clu->GetCellAmplitudeFraction(iDig);
+      if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    
+    if(!fCellsRecalibrated){
+
+      geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
+      geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                    
+      
+      if(IsRecalibrationOn()) {
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
+      }
     }
+    
     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
     
     weight = GetCellWeight(eCell,clEnergy);
     totalWeight += weight;
+    
     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
+    //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
-    
+    //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
+
     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
     
   }// cell loop
@@ -498,25 +1024,27 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet
     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
   }
     
-  //printf("iSM %d \n",iSupMod);
   //Float_t pos[]={0,0,0};
   //clu->GetPosition(pos);
   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
-  
-  
-  //printf("NewPos a: %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
+  //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
   
        if(iSupModMax > 1) {//sector 1
          newPos[0] +=fMisalTransShift[3];//-=3.093; 
          newPos[1] +=fMisalTransShift[4];//+=6.82;
          newPos[2] +=fMisalTransShift[5];//+=1.635;
+    //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
+
        }
        else {//sector 0
          newPos[0] +=fMisalTransShift[0];//+=1.134;
          newPos[1] +=fMisalTransShift[1];//+=8.2;
          newPos[2] +=fMisalTransShift[2];//+=1.197;
+    //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
+
        }
-  
+  //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
+
   clu->SetPosition(newPos);
   
 }  
@@ -536,9 +1064,10 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr
   Int_t iIphi   = -1, iIeta   = -1;
        Int_t iSupMod = -1, iSupModMax = -1;
   Int_t iphi = -1, ieta =-1;
-  
+  Bool_t shared = kFALSE;
+
   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
-  GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
+  GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
 
   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
@@ -549,17 +1078,24 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr
     absId = clu->GetCellAbsId(iDig);
     fraction  = clu->GetCellAmplitudeFraction(iDig);
     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-    geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
-    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);                  
-    
+
     if     (iDig==0)  startingSM = iSupMod;
     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 
     eCell  = cells->GetCellAmplitude(absId);
     
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+    geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);          
+    
+    if(!fCellsRecalibrated){
+      
+      if(IsRecalibrationOn()) {
+                
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+        
+      }
     }
+    
     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
     
     weight = GetCellWeight(eCell,clEnergy);
@@ -588,7 +1124,781 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr
   
 }
 
+//____________________________________________________________________________
+void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){           
+       
+  //re-evaluate distance to bad channel with updated bad map
+  
+  if(!fRecalDistToBadChannels) return;
+  
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return;
+  }  
+  
+       //Get channels map of the supermodule where the cluster is.
+  Int_t absIdMax       = -1, iSupMod =-1, icolM = -1, irowM = -1;
+  Bool_t shared = kFALSE;
+  GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
+  TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
+
+  Int_t dRrow, dRcol;  
+       Float_t  minDist = 10000.;
+       Float_t  dist    = 0.;
+  
+  //Loop on tower status map 
+       for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
+               for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
+                       //Check if tower is bad.
+                       if(hMap->GetBinContent(icol,irow)==0) continue;
+      //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
+      //       iSupMod,icol, irow, icolM,irowM);
+      
+      dRrow=TMath::Abs(irowM-irow);
+      dRcol=TMath::Abs(icolM-icol);
+      dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
+                       if(dist < minDist){
+        //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
+        minDist = dist;
+      }
+      
+               }
+       }
+  
+       //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
+       if (shared) {
+               TH2D* hMap2 = 0;
+               Int_t iSupMod2 = -1;
+    
+               //The only possible combinations are (0,1), (2,3) ... (8,9)
+               if(iSupMod%2) iSupMod2 = iSupMod-1;
+               else          iSupMod2 = iSupMod+1;
+               hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
+    
+               //Loop on tower status map of second super module
+               for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
+                       for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
+                               //Check if tower is bad.
+                               if(hMap2->GetBinContent(icol,irow)==0) continue;
+                               //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
+          //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
+        
+        dRrow=TMath::Abs(irow-irowM);
+        
+        if(iSupMod%2) {
+                                 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
+                               }
+        else {
+          dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
+                               }                    
+        
+                               dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
+        if(dist < minDist) minDist = dist;        
+        
+                       }
+               }
+    
+       }// shared cluster in 2 SuperModules
+  
+  AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
+  cluster->SetDistanceToBadChannel(minDist);
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){           
+       
+  //re-evaluate identification parameters with bayesian
+  
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return;
+  }
+  
+       if ( cluster->GetM02() != 0)
+    fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
+  
+  Float_t pidlist[AliPID::kSPECIESN+1];
+       for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
+  
+  cluster->SetPID(pidlist);
+       
+}
+
+//____________________________________________________________________________
+void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
+{
+  // Calculates new center of gravity in the local EMCAL-module coordinates 
+  // and tranfers into global ALICE coordinates
+  // Calculates Dispersion and main axis
+  
+  if(!cluster){
+    AliInfo("Cluster pointer null!");
+    return;
+  }
+    
+  Int_t nstat  = 0;
+  Float_t wtot = 0. ;
+  Double_t eCell       = 0.;
+  Float_t  fraction    = 1.;
+  Float_t  recalFactor = 1.;
+
+  Int_t iSupMod = -1;
+  Int_t iTower  = -1;
+  Int_t iIphi   = -1;
+  Int_t iIeta   = -1;
+  Int_t iphi    = -1;
+  Int_t ieta    = -1;
+  Double_t etai = -1.;
+  Double_t phii = -1.;
+  
+  Double_t w     = 0.;
+  Double_t d     = 0.;
+  Double_t dxx   = 0.;
+  Double_t dzz   = 0.;
+  Double_t dxz   = 0.;  
+  Double_t xmean = 0.;
+  Double_t zmean = 0.;
+    
+  //Loop on cells
+  for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
+    
+    //Get from the absid the supermodule, tower and eta/phi numbers
+    geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);        
+    
+    //Get the cell energy, if recalibration is on, apply factors
+    fraction  = cluster->GetCellAmplitudeFraction(iDigit);
+    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    
+    if(!fCellsRecalibrated){
+      
+      if(IsRecalibrationOn()) {
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+      }
+      
+    }
+    
+    eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
+    
+    if(cluster->E() > 0 && eCell > 0){
+      
+      w  = GetCellWeight(eCell,cluster->E());
+      
+      etai=(Double_t)ieta;
+      phii=(Double_t)iphi;             
+      if(w > 0.0) {
+        wtot += w ;
+        nstat++;                       
+        //Shower shape
+        dxx  += w * etai * etai ;
+        xmean+= w * etai ;
+        dzz  += w * phii * phii ;
+        zmean+= w * phii ; 
+        dxz  += w * etai * phii ; 
+      }
+    }
+    else
+      AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
+  }//cell loop
+  
+  //Normalize to the weight    
+  if (wtot > 0) {
+    xmean /= wtot ;
+    zmean /= wtot ;
+  }
+  else
+    AliError(Form("Wrong weight %f\n", wtot));
+  
+  //Calculate dispersion       
+  for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
+    
+    //Get from the absid the supermodule, tower and eta/phi numbers
+    geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
+    
+    //Get the cell energy, if recalibration is on, apply factors
+    fraction  = cluster->GetCellAmplitudeFraction(iDigit);
+    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    if(IsRecalibrationOn()) {
+      recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+    }
+    eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
+    
+    if(cluster->E() > 0 && eCell > 0){
+      
+      w  = GetCellWeight(eCell,cluster->E());
+      
+      etai=(Double_t)ieta;
+      phii=(Double_t)iphi;             
+      if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); 
+    }
+    else
+      AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
+  }// cell loop
+  
+  //Normalize to the weigth and set shower shape parameters
+  if (wtot > 0 && nstat > 1) {
+    d /= wtot ;
+    dxx /= wtot ;
+    dzz /= wtot ;
+    dxz /= wtot ;
+    dxx -= xmean * xmean ;
+    dzz -= zmean * zmean ;
+    dxz -= xmean * zmean ;
+    cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
+    cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
+  }
+  else{
+    d=0. ;
+    cluster->SetM20(0.) ;
+    cluster->SetM02(0.) ;
+  }    
+  
+  if (d>=0)
+    cluster->SetDispersion(TMath::Sqrt(d)) ;
+  else    
+    cluster->SetDispersion(0) ;
+}
+
+//____________________________________________________________________________
+void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  AliEMCALGeometry *geom)
+{
+  //This function should be called before the cluster loop
+  //Before call this function, please recalculate the cluster positions
+  //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
+  //Store matched cluster indexes and residuals
+  
+  fMatchedTrackIndex->Reset();
+  fMatchedClusterIndex->Reset();
+  fResidualPhi->Reset();
+  fResidualEta->Reset();
+  
+  fMatchedTrackIndex->Set(500);
+  fMatchedClusterIndex->Set(500);
+  fResidualPhi->Set(500);
+  fResidualEta->Set(500);
+  
+  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
+  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+  
+  Int_t    matched=0;
+  Double_t cv[21];
+  for (Int_t i=0; i<21;i++) cv[i]=0;
+  for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
+  {
+    AliExternalTrackParam *trackParam = 0;
+    
+    //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner 
+    if(esdevent)
+    {
+      AliESDtrack *esdTrack = esdevent->GetTrack(itr);
+      if(!esdTrack || !IsAccepted(esdTrack)) continue;
+      if(esdTrack->Pt()<fCutMinTrackPt) continue;
+      trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
+    }
+    
+    //If the input event is AOD, the starting point for extrapolation is at vertex
+    //AOD tracks are selected according to its bit.
+    else if(aodevent)
+    {
+      AliAODTrack *aodTrack = aodevent->GetTrack(itr);
+      if(!aodTrack) continue;
+      if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
+      if(aodTrack->Pt()<fCutMinTrackPt) continue;
+      Double_t pos[3],mom[3];
+      aodTrack->GetXYZ(pos);
+      aodTrack->GetPxPyPz(mom);
+      AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
+      trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
+    }
+    
+    //Return if the input data is not "AOD" or "ESD"
+    else
+    {
+      printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
+      return;
+    }
+    
+    if(!trackParam) continue;
+    
+    Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
+    Int_t index = -1;
+    if(!clusterArr){// get clusters from event
+      for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
+      {
+        AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
+        if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
+        AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation 
+        Float_t tmpEta=-999, tmpPhi=-999;
+        if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
+        if(fCutEtaPhiSum)
+        {
+          Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
+          if(tmpR<dRMax)
+          {
+            dRMax=tmpR;
+            dEtaMax=tmpEta;
+            dPhiMax=tmpPhi;
+            index=icl;
+          }
+        }
+        else if(fCutEtaPhiSeparate)
+        {
+          if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
+          {
+            dEtaMax = tmpEta;
+            dPhiMax = tmpPhi;
+            index=icl;
+          }
+        }
+        else
+        {
+          printf("Error: please specify your cut criteria\n");
+          printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
+          printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
+          if(aodevent && trackParam) delete trackParam;
+          return;
+        }
+      }//cluster loop
+    } 
+    else { // external cluster array, not from ESD event
+      for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
+      {
+        AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
+        if(!cluster){ 
+          AliInfo("Cluster not found!!!");
+          continue;
+        }
+        if(!cluster->IsEMCAL()) continue;
+        AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
+        Float_t tmpEta=-999, tmpPhi=-999;
+        if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
+        if(fCutEtaPhiSum)
+        {
+          Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
+          if(tmpR<dRMax)
+          {
+            dRMax=tmpR;
+            dEtaMax=tmpEta;
+            dPhiMax=tmpPhi;
+            index=icl;
+          }
+        }
+        else if(fCutEtaPhiSeparate)
+        {
+          if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
+          {
+            dEtaMax = tmpEta;
+            dPhiMax = tmpPhi;
+            index=icl;
+          }
+        }
+        else
+        {
+          printf("Error: please specify your cut criteria\n");
+          printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
+          printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
+          if(aodevent && trackParam) delete trackParam;
+          return;
+        }
+      }//cluster loop
+    }// external list of clusters
+    
+    if(index>-1)
+    {
+      fMatchedTrackIndex  ->AddAt(itr,matched);
+      fMatchedClusterIndex->AddAt(index,matched);
+      fResidualEta          ->AddAt(dEtaMax,matched);
+      fResidualPhi          ->AddAt(dPhiMax,matched);
+      matched++;
+    }
+    if(aodevent && trackParam) delete trackParam;
+  }//track loop
+  
+  AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
+  
+  fMatchedTrackIndex  ->Set(matched);
+  fMatchedClusterIndex->Set(matched);
+  fResidualPhi          ->Set(matched);
+  fResidualEta          ->Set(matched);
+}
+
+//________________________________________________________________________________
+Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
+{
+  //
+  // This function returns the index of matched cluster to input track
+  // Returns -1 if no match is found
+  
+  
+  Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
+  Int_t index = -1;
+  
+  AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
+  
+  if(!trackParam) return index;          
+  for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
+  {
+    AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
+    if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;  
+    AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
+    Float_t tmpEta=-999, tmpPhi=-999;
+    if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
+    if(fCutEtaPhiSum)
+    {
+      Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
+      if(tmpR<dRMax)
+           {
+             dRMax=tmpR;
+             dEtaMax=tmpEta;
+             dPhiMax=tmpPhi;
+             index=icl;
+           }
+    }
+    else if(fCutEtaPhiSeparate)
+    {
+      if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
+           {
+             dEtaMax = tmpEta;
+             dPhiMax = tmpPhi;
+             index=icl;
+           }
+    }
+    else
+    {
+      printf("Error: please specify your cut criteria\n");
+      printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
+      printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
+      return -1;
+    }
+  }//cluster loop
+  return index;
+}
+
+//________________________________________________________________________________
+Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
+{
+  //
+  //Return the residual by extrapolating a track to a cluster
+  //
+  if(!trkParam || !cluster) return kFALSE;
+  Float_t clsPos[3];
+  Double_t trkPos[3];
+  cluster->GetPosition(clsPos); //Has been recalculated
+  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
+  Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
+  vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
+  if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.8, -1)) return kFALSE;
+  trkParam->GetXYZ(trkPos); //Get the extrapolated global position
+
+  TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+
+  // track cluster matching
+  tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
+  tmpEta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
+
+  return kTRUE;
+}
+
+//________________________________________________________________________________
+void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
+{
+  //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
+  //Get the residuals dEta and dPhi for this cluster to the closest track
+  //Works with ESDs and AODs
+
+  if( FindMatchedPosForCluster(clsIndex) >= 999 )
+  {
+    AliDebug(2,"No matched tracks found!\n");
+    dEta=999.;
+    dPhi=999.;
+    return;
+  }
+  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
+  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
+}
+//________________________________________________________________________________
+void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
+{
+  //Given a track index as in AliESDEvent::GetTrack(trkIndex)
+  //Get the residuals dEta and dPhi for this track to the closest cluster
+  //Works with ESDs and AODs
+
+  if( FindMatchedPosForTrack(trkIndex) >= 999 )
+  {
+    AliDebug(2,"No matched cluster found!\n");
+    dEta=999.;
+    dPhi=999.;
+    return;
+  }
+  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
+  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
+}
+
+//__________________________________________________________
+Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
+{
+  //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
+  //Get the index of matched track to this cluster
+  //Works with ESDs and AODs
+  
+  if(IsClusterMatched(clsIndex))
+    return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
+  else 
+    return -1; 
+}
+
+//__________________________________________________________
+Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
+{
+  //Given a track index as in AliESDEvent::GetTrack(trkIndex)
+  //Get the index of matched cluster to this track
+  //Works with ESDs and AODs
+  
+  if(IsTrackMatched(trkIndex))
+    return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
+  else 
+    return -1; 
+}
+
 //__________________________________________________
+Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
+{
+  //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
+  //Returns if the cluster has a match
+  if(FindMatchedPosForCluster(clsIndex) < 999) 
+    return kTRUE;
+  else
+    return kFALSE;
+}
+
+//__________________________________________________
+Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const 
+{
+  //Given a track index as in AliESDEvent::GetTrack(trkIndex)
+  //Returns if the track has a match
+  if(FindMatchedPosForTrack(trkIndex) < 999) 
+    return kTRUE;
+  else
+    return kFALSE;
+}
+
+//__________________________________________________________
+UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
+{
+  //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
+  //Returns the position of the match in the fMatchedClusterIndex array
+  Float_t tmpR = fCutR;
+  UInt_t pos = 999;
+  
+  for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
+  {
+    if(fMatchedClusterIndex->At(i)==clsIndex)
+      {
+       Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
+       if(r<tmpR)
+         {
+           pos=i;
+           tmpR=r;
+           AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
+         }
+      }
+  }
+  return pos;
+}
+
+//__________________________________________________________
+UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
+{
+  //Given a track index as in AliESDEvent::GetTrack(trkIndex)
+  //Returns the position of the match in the fMatchedTrackIndex array
+  Float_t tmpR = fCutR;
+  UInt_t pos = 999;
+  
+  for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
+  {
+    if(fMatchedTrackIndex->At(i)==trkIndex)
+      {
+       Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
+       if(r<tmpR)
+         {
+           pos=i;
+           tmpR=r;
+           AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
+         }
+      }
+  }
+  return pos;
+}
+
+//__________________________________________________________
+Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
+{
+  // check if the cluster survives some quality cut
+  //
+  //
+  Bool_t isGood=kTRUE;
+  if(!cluster || !cluster->IsEMCAL()) return kFALSE;
+  if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
+  if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
+  if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
+
+  return isGood;
+}
+
+//__________________________________________________________
+Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
+{
+  // Given a esd track, return whether the track survive all the cuts
+
+  // The different quality parameter are first
+  // retrieved from the track. then it is found out what cuts the
+  // track did not survive and finally the cuts are imposed.
+
+  UInt_t status = esdTrack->GetStatus();
+
+  Int_t nClustersITS = esdTrack->GetITSclusters(0);
+  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
+
+  Float_t chi2PerClusterITS = -1;
+  Float_t chi2PerClusterTPC = -1;
+  if (nClustersITS!=0)
+    chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
+  if (nClustersTPC!=0) 
+    chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
+
+
+  //DCA cuts
+  if(fTrackCutsType==kGlobalCut)
+    {
+      Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
+      //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
+      SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
+    }
+
+
+  Float_t b[2];
+  Float_t bCov[3];
+  esdTrack->GetImpactParameters(b,bCov);
+  if (bCov[0]<=0 || bCov[2]<=0) {
+    AliDebug(1, "Estimated b resolution lower or equal zero!");
+    bCov[0]=0; bCov[2]=0;
+  }
+
+  Float_t dcaToVertexXY = b[0];
+  Float_t dcaToVertexZ = b[1];
+  Float_t dcaToVertex = -1;
+
+  if (fCutDCAToVertex2D)
+    dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
+  else
+    dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
+    
+  // cut the track?
+  
+  Bool_t cuts[kNCuts];
+  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
+  
+  // track quality cuts
+  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
+    cuts[0]=kTRUE;
+  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
+    cuts[1]=kTRUE;
+  if (nClustersTPC<fCutMinNClusterTPC)
+    cuts[2]=kTRUE;
+  if (nClustersITS<fCutMinNClusterITS) 
+    cuts[3]=kTRUE;
+  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
+    cuts[4]=kTRUE; 
+  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
+    cuts[5]=kTRUE;  
+  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
+    cuts[6]=kTRUE;
+  if (fCutDCAToVertex2D && dcaToVertex > 1)
+    cuts[7] = kTRUE;
+  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
+    cuts[8] = kTRUE;
+  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
+    cuts[9] = kTRUE;
+
+  if(fTrackCutsType==kGlobalCut)
+    {
+      //Require at least one SPD point + anything else in ITS
+      if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
+       cuts[10] = kTRUE;
+    }
+
+  Bool_t cut=kFALSE;
+  for (Int_t i=0; i<kNCuts; i++)
+    if (cuts[i]) {cut = kTRUE;}
+
+    // cut the track
+  if (cut) 
+    return kFALSE;
+  else 
+    return kTRUE;
+}
+
+
+//__________________________________________________
+void AliEMCALRecoUtils::InitTrackCuts()
+{
+  //Intilize the track cut criteria
+  //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
+  //Also you can customize the cuts using the setters
+  
+  switch (fTrackCutsType)
+    {
+    case kTPCOnlyCut:
+      {
+       AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
+       //TPC
+       SetMinNClustersTPC(70);
+       SetMaxChi2PerClusterTPC(4);
+       SetAcceptKinkDaughters(kFALSE);
+       SetRequireTPCRefit(kFALSE);
+
+       //ITS
+       SetRequireITSRefit(kFALSE);
+       SetMaxDCAToVertexZ(3.2);
+       SetMaxDCAToVertexXY(2.4);
+       SetDCAToVertex2D(kTRUE);
+
+       break;
+      }
+    
+    case kGlobalCut:
+      {
+       AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
+       //TPC
+       SetMinNClustersTPC(70);
+       SetMaxChi2PerClusterTPC(4);
+       SetAcceptKinkDaughters(kFALSE);
+       SetRequireTPCRefit(kTRUE);
+
+       //ITS
+       SetRequireITSRefit(kTRUE);
+       SetMaxDCAToVertexZ(2);
+       SetMaxDCAToVertexXY();
+       SetDCAToVertex2D(kFALSE);
+
+       break;
+      }
+
+    case kLooseCut:
+      {
+       AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
+       SetMinNClustersTPC(50);
+       SetAcceptKinkDaughters(kTRUE);
+
+       break;
+      }
+    }
+}
+
+//___________________________________________________
 void AliEMCALRecoUtils::Print(const Option_t *) const 
 {
   // Print Parameters
@@ -602,5 +1912,63 @@ void AliEMCALRecoUtils::Print(const Option_t *) const
   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
   
   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
-    
+
+  printf("Matching criteria: ");
+  if(fCutEtaPhiSum)
+    {
+      printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
+    }
+  else if(fCutEtaPhiSeparate)
+    {
+      printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
+    }
+  else
+    {
+      printf("Error\n");
+      printf("please specify your cut criteria\n");
+      printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
+      printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
+    }
+
+  printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
+
+  printf("Track cuts: \n");
+  printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
+  printf("AOD track selection mask: %d\n",fAODFilterMask);
+  printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
+  printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
+  printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
+  printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
+  printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
+
 }
+
+//_____________________________________________________________________
+void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber){
+  //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
+  //Do it only once and only if it is requested
+  
+  if(!fUseRunCorrectionFactors) return;
+  if(fRunCorrectionFactorsSet)  return;
+  
+  AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber));
+  AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
+  corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
+  
+  SwitchOnRecalibration();
+  for(Int_t ism = 0; ism < 4; ism++){
+    for(Int_t icol = 0; icol < 48; icol++){
+      for(Int_t irow = 0; irow < 24; irow++){
+        Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
+        Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
+        GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
+        //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow, 
+        //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
+        //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
+      }
+    }
+  }
+   fRunCorrectionFactorsSet = kTRUE;
+}
+