In case of recalibrating all cells, check not only the energy recalibration flag...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx
index d985579..db11d25 100644 (file)
@@ -25,7 +25,6 @@
 //          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 "AliVCluster.h"
 #include "AliEMCALRecoUtils.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCALTrack.h"
-#include "AliEMCALCalibTimeDepCorrection.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), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
-  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
-  fAODFilterMask(32),
-  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0), 
-  fResidualZ(0x0), fResidualR(0x0), fCutR(10), fCutZ(10), fMass(0.139), fStep(1),
-  fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
-  fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
-  fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
-  fUseTimeCorrectionFactors(kFALSE),  fTimeCorrectionFactorsSet(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.
@@ -103,48 +111,57 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
   //fNonLinearityParams[2] = 1.046;
 
-  //Track matching
-  fMatchedTrackIndex   = new TArrayI();
-  fMatchedClusterIndex = new TArrayI();
-  fResidualZ           = new TArrayF();
-  fResidualR           = new TArrayF();
+  //Cluster energy smearing
+  fSmearClusterEnergy   = kFALSE;
+  fSmearClusterParam[0] = 0.07; // * sqrt E term
+  fSmearClusterParam[1] = 0.00; // * E term
+  fSmearClusterParam[2] = 0.00; // constant
   
-  InitTrackCuts();
-  
-  fPIDUtils            = new AliEMCALPIDUtils();
-
+  //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),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
+: 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),
+  fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
+  fRejectExoticCluster(reco.fRejectExoticCluster),           fPIDUtils(reco.fPIDUtils), 
   fAODFilterMask(reco.fAODFilterMask),
-  fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
+  fMatchedTrackIndex(  reco.fMatchedTrackIndex?  new TArrayI(*reco.fMatchedTrackIndex):0x0),
   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
-  fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
-  fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
-  fCutR(reco.fCutR),fCutZ(reco.fCutZ),fMass(reco.fMass), fStep(reco.fStep),
-  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),
-  fPIDUtils(reco.fPIDUtils), 
-  fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors),  fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
+  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 < 7  ; 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]  ; }
 
 }
 
@@ -157,29 +174,52 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   if(this == &reco)return *this;
   ((TNamed *)this)->operator=(reco);
 
-  fNonLinearityFunction      = reco.fNonLinearityFunction;
+  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;
+
+  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;
 
-  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]; 
-
-  fAODFilterMask              = reco.fAODFilterMask;
+  fAODFilterMask             = reco.fAODFilterMask;
   
+  fCutEtaPhiSum              = reco.fCutEtaPhiSum;
+  fCutEtaPhiSeparate         = reco.fCutEtaPhiSeparate;
   fCutR                      = reco.fCutR;
-  fCutZ                      = reco.fCutZ;
+  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;
@@ -190,35 +230,29 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
-
-  fPIDUtils                  = reco.fPIDUtils;
   
-  fUseTimeCorrectionFactors  = reco.fUseTimeCorrectionFactors;
-  fTimeCorrectionFactorsSet  = reco.fTimeCorrectionFactorsSet;
-
-  
-  if(reco.fResidualR){
+  if(reco.fResidualEta){
     // assign or copy construct
-    if(fResidualR){ 
-      *fResidualR = *reco.fResidualR;
+    if(fResidualEta){ 
+      *fResidualEta = *reco.fResidualEta;
     }
-    else fResidualR = new TArrayF(*reco.fResidualR);
+    else fResidualEta = new TArrayF(*reco.fResidualEta);
   }
   else{
-    if(fResidualR)delete fResidualR;
-    fResidualR = 0;
+    if(fResidualEta)delete fResidualEta;
+    fResidualEta = 0;
   }
   
-  if(reco.fResidualZ){
+  if(reco.fResidualPhi){
     // assign or copy construct
-    if(fResidualZ){ 
-      *fResidualZ = *reco.fResidualZ;
+    if(fResidualPhi){ 
+      *fResidualPhi = *reco.fResidualPhi;
     }
-    else fResidualZ = new TArrayF(*reco.fResidualZ);
+    else fResidualPhi = new TArrayF(*reco.fResidualPhi);
   }
   else{
-    if(fResidualZ)delete fResidualZ;
-    fResidualZ = 0;
+    if(fResidualPhi)delete fResidualPhi;
+    fResidualPhi = 0;
   }
   
   if(reco.fMatchedTrackIndex){
@@ -244,8 +278,7 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
     fMatchedClusterIndex = 0;
   }
-  
-  
+   
   return *this;
 }
 
@@ -255,21 +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;
+  }
  
-  if(fMatchedTrackIndex)   {delete fMatchedTrackIndex;   fMatchedTrackIndex=0;}
-  if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
-  if(fResidualR)           {delete fResidualR;           fResidualR=0;}
-  if(fResidualZ)           {delete fResidualZ;           fResidualZ=0;}
+  delete fMatchedTrackIndex   ; 
+  delete fMatchedClusterIndex ; 
+  delete fResidualEta         ; 
+  delete fResidualPhi         ; 
+  delete fPIDUtils            ;
 
+  InitTrackCuts();
 }
 
 //_______________________________________________________________
@@ -278,6 +318,11 @@ 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;
   
@@ -337,7 +382,7 @@ Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCl
 
 
 //_________________________________________________________________________________________________________
-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
        
@@ -365,9 +410,57 @@ 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) {
@@ -445,6 +538,68 @@ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
 
 }
 //__________________________________________________
+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 
 {
   //Calculate shower depth for a given cluster energy and particle type
@@ -510,7 +665,13 @@ void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells*
   Int_t iIphi   = -1;
   Int_t iIeta   = -1;
   Int_t iSupMod0= -1;
-       //printf("---Max?\n");
+
+  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);
@@ -554,9 +715,9 @@ void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
        TH1::AddDirectory(kFALSE);
   
        fEMCALRecalibrationFactors = new TObjArray(10);
-       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));
+       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.);
@@ -570,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(){
@@ -584,9 +770,7 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
        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();
        
@@ -595,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() ;
@@ -605,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.");
@@ -669,20 +991,26 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet
   //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);
-    //printf("cell energy %f, weight %f\n",eCell,weight);
     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);
@@ -750,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);
@@ -796,6 +1131,11 @@ void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry
   
   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;
@@ -869,7 +1209,12 @@ void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry
 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());
   
@@ -887,6 +1232,11 @@ void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry
   // 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.;
@@ -920,9 +1270,15 @@ void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry
     //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);
+    
+    if(!fCellsRecalibrated){
+      
+      if(IsRecalibrationOn()) {
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+      }
+      
     }
+    
     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
     
     if(cluster->E() > 0 && eCell > 0){
@@ -1006,165 +1362,222 @@ void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry
 }
 
 //____________________________________________________________________________
-void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr, TString dataType)
+void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  AliEMCALGeometry *geom)
 {
-  //Use dataType to indicate the input event is AOD or ESD
   //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();
+  
+  fMatchedTrackIndex->Reset();
   fMatchedClusterIndex->Reset();
-  fResidualZ          ->Reset();
-  fResidualR          ->Reset();
+  fResidualPhi->Reset();
+  fResidualEta->Reset();
   
-  fMatchedTrackIndex  ->Set(500);
+  fMatchedTrackIndex->Set(500);
   fMatchedClusterIndex->Set(500);
-  fResidualZ          ->Set(500);
-  fResidualR          ->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;
-
+    AliExternalTrackParam *trackParam = 0;
+    
     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner 
-    if(dataType.Contains("ESD"))
-      {
-       AliESDtrack *esdTrack = ((AliESDEvent*)event)->GetTrack(itr);
-       if(!esdTrack || !IsAccepted(esdTrack)) continue;
-       const AliESDfriendTrack*  friendTrack = esdTrack->GetFriendTrack();
-       if(friendTrack && friendTrack->GetTPCOut())
-         {
-           //Use TPC Out as starting point if it is available
-           trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
-         }
-       else
-         {
-           //Otherwise use TPC inner
-           trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
-         }
-      }
+    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(dataType.Contains("AOD"))
-      {
-       AliAODTrack *aodTrack = ((AliAODEvent*)event)->GetTrack(itr);
-       if(!aodTrack) continue;
-       if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
-       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());
-      }
+    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 %s! Should be \"AOD\" or \"ESD\"\n",dataType.Data());
-       return;
-      }
-  
+    {
+      printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
+      return;
+    }
+    
     if(!trackParam) continue;
-
-    Float_t dRMax = fCutR, dZMax=fCutZ;
+    
+    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++)
       {
-       AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
         AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
-        if(!cluster->IsEMCAL()) continue;                              
-       Float_t tmpR=-1, tmpZ=-1;
-       if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
-        if(tmpR<dRMax)
+        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)
         {
-          dRMax=tmpR;
-          dZMax=tmpZ;
-          index=icl;
+          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;
         }
-       delete trkPamTmp;
       }//cluster loop
-    } else { // external cluster array, not from ESD event
+    } 
+    else { // external cluster array, not from ESD event
       for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
       {
-       AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
-        AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
+        AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
+        if(!cluster){ 
+          AliInfo("Cluster not found!!!");
+          continue;
+        }
         if(!cluster->IsEMCAL()) continue;
-       Float_t tmpR=-1, tmpZ=-1;
-       if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
-        if(tmpR<dRMax)
+        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)
         {
-          dRMax=tmpR;
-          dZMax=tmpZ;
-          index=icl;
+          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;
         }
-       delete trkPamTmp;
       }//cluster loop
     }// external list of clusters
-
+    
     if(index>-1)
     {
       fMatchedTrackIndex  ->AddAt(itr,matched);
       fMatchedClusterIndex->AddAt(index,matched);
-      fResidualZ          ->AddAt(dZMax,matched);
-      fResidualR          ->AddAt(dRMax,matched);
+      fResidualEta          ->AddAt(dEtaMax,matched);
+      fResidualPhi          ->AddAt(dPhiMax,matched);
       matched++;
     }
-    delete trackParam;
+    if(aodevent && trackParam) delete trackParam;
   }//track loop
   
   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
   
   fMatchedTrackIndex  ->Set(matched);
   fMatchedClusterIndex->Set(matched);
-  fResidualZ          ->Set(matched);
-  fResidualR          ->Set(matched);
+  fResidualPhi          ->Set(matched);
+  fResidualEta          ->Set(matched);
 }
 
 //________________________________________________________________________________
-Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event)
+Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
 {
   //
   // This function returns the index of matched cluster to input track
-  // Cut on match is dR<10cm by default. Returns -1 if no match is found
-
-
-  Float_t dRMax = fCutR;
+  // Returns -1 if no match is found
+  
+  
+  Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
   Int_t index = -1;
-
-  AliExternalTrackParam *trackParam=0;
-  const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
-  if(friendTrack && friendTrack->GetTPCOut())
-    trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
-  else
-    trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
-
+  
+  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)
     {
-      AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
-      AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
-      if(!cluster->IsEMCAL()) continue;                                
-      Float_t tmpR=-1, tmpZ=-1;
-      if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
-      if(tmpR>-1 && tmpR<dRMax)
-        {
-          dRMax=tmpR;
-          index=icl;
-        }
-       delete trkPamTmp;
-      }//cluster loop
+      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 &tmpR, Float_t &tmpZ)
+Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
 {
   //
   //Return the residual by extrapolating a track to a cluster
@@ -1176,47 +1589,52 @@ Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkP
   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
-  trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
-  if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE; 
+  if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.8, -1)) return kFALSE;
   trkParam->GetXYZ(trkPos); //Get the extrapolated global position
-  tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
-  tmpZ = clsPos[2]-trkPos[2];
+
+  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 &dR, Float_t &dZ)
+void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
 {
   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
-  //Get the residuals dR and dZ for this cluster to the closest track
+  //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");
-    dR=999.;
-    dZ=999.;
+    dEta=999.;
+    dPhi=999.;
     return;
   }
-  dR = fResidualR->At(FindMatchedPosForCluster(clsIndex));
-  dZ = fResidualZ->At(FindMatchedPosForCluster(clsIndex));
+  dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
+  dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
 }
 //________________________________________________________________________________
-void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dR, Float_t &dZ)
+void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
 {
   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
-  //Get the residuals dR and dZ for this track to the closest cluster
+  //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");
-    dR=999.;
-    dZ=999.;
+    dEta=999.;
+    dPhi=999.;
     return;
   }
-  dR = fResidualR->At(FindMatchedPosForTrack(trkIndex));
-  dZ = fResidualZ->At(FindMatchedPosForTrack(trkIndex));
+  dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
+  dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
 }
 
 //__________________________________________________________
@@ -1246,7 +1664,7 @@ Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
 }
 
 //__________________________________________________
-Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
+Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
 {
   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
   //Returns if the cluster has a match
@@ -1257,7 +1675,7 @@ Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
 }
 
 //__________________________________________________
-Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
+Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const 
 {
   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
   //Returns if the track has a match
@@ -1277,12 +1695,16 @@ UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
   
   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
   {
-    if(fMatchedClusterIndex->At(i)==clsIndex && fResidualR->At(i)<tmpR)
-    {
-      pos=i;
-      tmpR=fResidualR->At(i);
-      AliDebug(3,Form("Matched cluster index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(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;
 }
@@ -1297,17 +1719,36 @@ UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
   
   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
   {
-    if(fMatchedTrackIndex->At(i)==trkIndex && fResidualR->At(i)<tmpR)
-    {
-      pos=i;
-      tmpR=fResidualR->At(i);
-      AliDebug(3,Form("Matched track index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualR->At(i),fResidualZ->At(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
@@ -1330,9 +1771,12 @@ Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
 
 
   //DCA cuts
-  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
+  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];
@@ -1379,12 +1823,15 @@ Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
     cuts[9] = kTRUE;
 
-  //Require at least one SPD point + anything else in ITS
-  if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
-    cuts[10] = 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++) 
+  for (Int_t i=0; i<kNCuts; i++)
     if (cuts[i]) {cut = kTRUE;}
 
     // cut the track
@@ -1393,25 +1840,62 @@ Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
   else 
     return kTRUE;
 }
+
+
 //__________________________________________________
 void AliEMCALRecoUtils::InitTrackCuts()
 {
   //Intilize the track cut criteria
-  //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
+  //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
   //Also you can customize the cuts using the setters
   
-  //TPC
-  SetMinNClustersTPC(70);
-  SetMaxChi2PerClusterTPC(4);
-  SetAcceptKinkDaughters(kFALSE);
-  SetRequireTPCRefit(kTRUE);
-  
-  //ITS
-  SetRequireITSRefit(kTRUE);
-  SetMaxDCAToVertexZ(2);
-  SetDCAToVertex2D(kFALSE);
-  SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
-  SetMinNClustersITS();
+  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;
+      }
+    }
 }
 
 //___________________________________________________
@@ -1429,10 +1913,27 @@ void AliEMCALRecoUtils::Print(const Option_t *) const
   
   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
 
-  printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
-  printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
+  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);
@@ -1443,14 +1944,14 @@ void AliEMCALRecoUtils::Print(const Option_t *) const
 }
 
 //_____________________________________________________________________
-void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
+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(!fUseTimeCorrectionFactors) return;
-  if(fTimeCorrectionFactorsSet)  return;
+  if(!fUseRunCorrectionFactors) return;
+  if(fRunCorrectionFactorsSet)  return;
   
-  printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
+  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));
@@ -1468,6 +1969,6 @@ void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
       }
     }
   }
-   fTimeCorrectionFactorsSet = kTRUE;
+   fRunCorrectionFactorsSet = kTRUE;
 }