In case of recalibrating all cells, check not only the energy recalibration flag...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx
index 8b3552a..db11d25 100644 (file)
@@ -37,7 +37,7 @@
 #include <TH2F.h>
 #include <TArrayI.h>
 #include <TArrayF.h>
-#include "TObjArray.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.), fNonLinearThreshold(30),
-  fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
-  fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
-  fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
-  fAODFilterMask(32),
-  fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0), 
-  fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kTRUE), fCutEtaPhiSeparate(kFALSE), fCutR(0.1), fCutEta(0.02), fCutPhi(0.04), fMass(0.139), fStep(50),
-  fRejectExoticCluster(kFALSE),
-  fTrackCutsType(kTPCOnlyCut), fCutMinTrackPt(0), fCutMinNClusterTPC(-1), fCutMinNClusterITS(-1), fCutMaxChi2PerClusterTPC(1e10), fCutMaxChi2PerClusterITS(1e10),
-  fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
-  fCutMaxDCAToVertexXY(1e10), fCutMaxDCAToVertexZ(1e10),fCutDCAToVertex2D(kFALSE),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.
@@ -107,47 +111,57 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   //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();
-  
-  InitTrackCuts();
   fPIDUtils              = new AliEMCALPIDUtils();
+
+  InitTrackCuts();
 }
 
 //______________________________________________________________________
 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) 
-: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction), 
-  fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), fNonLinearThreshold(reco.fNonLinearThreshold),
-  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),
-  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),
-  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),
-  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]  ; }
 
 }
 
@@ -160,24 +174,40 @@ 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;
@@ -200,12 +230,6 @@ 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.fResidualEta){
     // assign or copy construct
@@ -264,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;
+  }
  
   delete fMatchedTrackIndex   ; 
   delete fMatchedClusterIndex ; 
   delete fResidualEta         ; 
   delete fResidualPhi         ; 
+  delete fPIDUtils            ;
 
+  InitTrackCuts();
 }
 
 //_______________________________________________________________
@@ -389,11 +420,36 @@ Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
     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 ;
 
-  if(cluster->GetNCells()<(1+cluster->E()/3.))
-    return kTRUE;
-  else
-    return kFALSE;
 }
 
 //__________________________________________________
@@ -405,7 +461,6 @@ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
     return 0;
   }
   
-  
   Float_t energy = cluster->E();
   
   switch (fNonLinearityFunction) {
@@ -676,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(){
@@ -690,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();
        
@@ -701,8 +779,10 @@ 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!");
@@ -716,34 +796,160 @@ 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)
 {
@@ -785,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);
@@ -866,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);
@@ -1051,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){
@@ -1170,17 +1395,7 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  Al
       AliESDtrack *esdTrack = esdevent->GetTrack(itr);
       if(!esdTrack || !IsAccepted(esdTrack)) continue;
       if(esdTrack->Pt()<fCutMinTrackPt) continue;
-      const AliESDfriendTrack*  friendTrack = esdTrack->GetFriendTrack();
-      if(friendTrack && friendTrack->GetTPCOut())
-      {
-        //Use TPC Out as starting point if it is available
-        trackParam=  const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
-      }
-      else
-      {
-        //Otherwise use TPC inner
-        trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
-      }
+      trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
     }
     
     //If the input event is AOD, the starting point for extrapolation is at vertex
@@ -1320,12 +1535,7 @@ Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event
   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++)
@@ -1379,18 +1589,14 @@ 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, 0.8, -1)) return kFALSE; 
+  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]);
 
-  Float_t clsPhi = (Float_t)clsPosVec.Phi();
-  if(clsPhi<0) clsPhi+=2*TMath::Pi();
-  Float_t trkPhi = (Float_t)trkPosVec.Phi();
-  if(trkPhi<0) trkPhi+=2*TMath::Pi();
-  tmpPhi = clsPhi-trkPhi;  // track cluster matching
+  // track cluster matching
+  tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
 
   return kTRUE;
@@ -1565,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];
@@ -1614,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
@@ -1628,6 +1840,8 @@ Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
   else 
     return kTRUE;
 }
+
+
 //__________________________________________________
 void AliEMCALRecoUtils::InitTrackCuts()
 {
@@ -1639,7 +1853,7 @@ void AliEMCALRecoUtils::InitTrackCuts()
     {
     case kTPCOnlyCut:
       {
-       AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()\n"));
+       AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
        //TPC
        SetMinNClustersTPC(70);
        SetMaxChi2PerClusterTPC(4);
@@ -1657,7 +1871,7 @@ void AliEMCALRecoUtils::InitTrackCuts()
     
     case kGlobalCut:
       {
-       AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)\n"));
+       AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
        //TPC
        SetMinNClustersTPC(70);
        SetMaxChi2PerClusterTPC(4);
@@ -1672,6 +1886,15 @@ void AliEMCALRecoUtils::InitTrackCuts()
 
        break;
       }
+
+    case kLooseCut:
+      {
+       AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
+       SetMinNClustersTPC(50);
+       SetAcceptKinkDaughters(kTRUE);
+
+       break;
+      }
     }
 }
 
@@ -1721,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));
@@ -1746,6 +1969,6 @@ void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
       }
     }
   }
-   fTimeCorrectionFactorsSet = kTRUE;
+   fRunCorrectionFactorsSet = kTRUE;
 }