EMCAL/AliEMCALRecoUtils.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 // Class AliEMCALRecoUtils
  21 // Some utilities to recalculate the cluster position or energy linearity
  22 //
  23 //
  24 // Author:  Gustavo Conesa (LPSC- Grenoble)
  25 //          Track matching part: Rongrong Ma (Yale)
  26
  27 ///////////////////////////////////////////////////////////////////////////////
  28 // --- standard c ---
  29
  30 // standard C++ includes
  31 //#include <Riostream.h>
  32
  33 // ROOT includes
  34 #include <TGeoManager.h>
  35 #include <TGeoMatrix.h>
  36 #include <TGeoBBox.h>
  37 #include <TH2F.h>
  38 #include <TArrayI.h>
  39 #include <TArrayF.h>
  40 #include <TObjArray.h>
  41
  42 // STEER includes
  43 #include "AliVCluster.h"
  44 #include "AliVCaloCells.h"
  45 #include "AliVEvent.h"
  46 #include "AliLog.h"
  47 #include "AliPID.h"
  48 #include "AliESDEvent.h"
  49 #include "AliAODEvent.h"
  50 #include "AliESDtrack.h"
  51 #include "AliAODTrack.h"
  52 #include "AliExternalTrackParam.h"
  53 #include "AliESDfriendTrack.h"
  54 #include "AliTrackerBase.h"
  55
  56 // EMCAL includes
  57 #include "AliEMCALRecoUtils.h"
  58 #include "AliEMCALGeometry.h"
  59 #include "AliTrackerBase.h"
  60 #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
  61 #include "AliEMCALPIDUtils.h"
  62
  63
  64 ClassImp(AliEMCALRecoUtils)
  65
  66 //______________________________________________
  67 AliEMCALRecoUtils::AliEMCALRecoUtils():
  68   fParticleType(kPhoton),                 fPosAlgo(kUnchanged),                   fW0(4.),
  69   fNonLinearityFunction(kNoCorrection),   fNonLinearThreshold(30),
  70   fSmearClusterEnergy(kFALSE),            fRandom(),
  71   fCellsRecalibrated(kFALSE),             fRecalibration(kFALSE),                 fEMCALRecalibrationFactors(),
  72   fTimeRecalibration(kFALSE),             fEMCALTimeRecalibrationFactors(),
  73   fUseRunCorrectionFactors(kFALSE),       fRunCorrectionFactorsSet(kFALSE),
  74   fRemoveBadChannels(kFALSE),             fRecalDistToBadChannels(kFALSE),        fEMCALBadChannelMap(),
  75   fNCellsFromEMCALBorder(0),              fNoEMCALBorderAtEta0(kTRUE),
  76   fRejectExoticCluster(kFALSE),           fRejectExoticCells(kFALSE),
  77   fExoticCellFraction(0.97),              fExoticCellDiffTime(1e6),               fExoticCellMinAmplitude(0.5),
  78   fPIDUtils(),                            fAODFilterMask(32),
  79   fMatchedTrackIndex(0x0),                fMatchedClusterIndex(0x0),
  80   fResidualEta(0x0), fResidualPhi(0x0),   fCutEtaPhiSum(kTRUE),                   fCutEtaPhiSeparate(kFALSE),
  81   fCutR(0.1),                             fCutEta(0.025),                         fCutPhi(0.05),
  82   fClusterWindow(100),                    fMass(0.139),
  83   fStepSurface(10.),                      fStepCluster(5.),
  84   fTrackCutsType(kLooseCut),              fCutMinTrackPt(0),                      fCutMinNClusterTPC(-1),
  85   fCutMinNClusterITS(-1),                 fCutMaxChi2PerClusterTPC(1e10),         fCutMaxChi2PerClusterITS(1e10),
  86   fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
  87   fCutMaxDCAToVertexXY(1e10),             fCutMaxDCAToVertexZ(1e10),              fCutDCAToVertex2D(kFALSE)
  88 {
  89 //
  90   // Constructor.
  91   // Initialize all constant values which have to be used
  92   // during Reco algorithm execution
  93   //
  94
  95   //Misalignment matrices
  96   for(Int_t i = 0; i < 15 ; i++) {
  97       fMisalTransShift[i] = 0.;
  98       fMisalRotShift[i]   = 0.;
  99   }
 100
 101   //Non linearity
 102   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = 0.;
 103
 104   //For kBeamTestCorrected case, but default is no correction
 105   fNonLinearityParams[0] =  0.99078;
 106   fNonLinearityParams[1] =  0.161499;
 107   fNonLinearityParams[2] =  0.655166;
 108   fNonLinearityParams[3] =  0.134101;
 109   fNonLinearityParams[4] =  163.282;
 110   fNonLinearityParams[5] =  23.6904;
 111   fNonLinearityParams[6] =  0.978;
 112
 113   //For kPi0GammaGamma case
 114   //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
 115   //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
 116   //fNonLinearityParams[2] = 1.046;
 117
 118   //Cluster energy smearing
 119   fSmearClusterEnergy   = kFALSE;
 120   fSmearClusterParam[0] = 0.07; // * sqrt E term
 121   fSmearClusterParam[1] = 0.00; // * E term
 122   fSmearClusterParam[2] = 0.00; // constant
 123
 124   //Track matching
 125   fMatchedTrackIndex     = new TArrayI();
 126   fMatchedClusterIndex   = new TArrayI();
 127   fResidualPhi           = new TArrayF();
 128   fResidualEta           = new TArrayF();
 129   fPIDUtils              = new AliEMCALPIDUtils();
 130
 131   InitTrackCuts();
 132 }
 133
 134 //______________________________________________________________________
 135 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
 136 : TNamed(reco),
 137   fParticleType(reco.fParticleType),                         fPosAlgo(reco.fPosAlgo),     fW0(reco.fW0),
 138   fNonLinearityFunction(reco.fNonLinearityFunction),         fNonLinearThreshold(reco.fNonLinearThreshold),
 139   fSmearClusterEnergy(reco.fSmearClusterEnergy),             fRandom(),
 140   fCellsRecalibrated(reco.fCellsRecalibrated),
 141   fRecalibration(reco.fRecalibration),                       fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
 142   fTimeRecalibration(reco.fTimeRecalibration),               fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
 143   fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),   fRunCorrectionFactorsSet(reco.fRunCorrectionFactorsSet),
 144   fRemoveBadChannels(reco.fRemoveBadChannels),               fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
 145   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
 146   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
 147   fRejectExoticCluster(reco.fRejectExoticCluster),           fRejectExoticCells(reco.fRejectExoticCells),
 148   fExoticCellFraction(reco.fExoticCellFraction),             fExoticCellDiffTime(reco.fExoticCellDiffTime),
 149   fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
 150   fPIDUtils(reco.fPIDUtils),                                 fAODFilterMask(reco.fAODFilterMask),
 151   fMatchedTrackIndex(  reco.fMatchedTrackIndex?  new TArrayI(*reco.fMatchedTrackIndex):0x0),
 152   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
 153   fResidualEta(        reco.fResidualEta?        new TArrayF(*reco.fResidualEta):0x0),
 154   fResidualPhi(        reco.fResidualPhi?        new TArrayF(*reco.fResidualPhi):0x0),
 155   fCutEtaPhiSum(reco.fCutEtaPhiSum),                         fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
 156   fCutR(reco.fCutR),        fCutEta(reco.fCutEta),           fCutPhi(reco.fCutPhi),
 157   fClusterWindow(reco.fClusterWindow),
 158   fMass(reco.fMass),        fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
 159   fTrackCutsType(reco.fTrackCutsType),                       fCutMinTrackPt(reco.fCutMinTrackPt),
 160   fCutMinNClusterTPC(reco.fCutMinNClusterTPC),               fCutMinNClusterITS(reco.fCutMinNClusterITS),
 161   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC),   fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
 162   fCutRequireTPCRefit(reco.fCutRequireTPCRefit),             fCutRequireITSRefit(reco.fCutRequireITSRefit),
 163   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),     fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
 164   fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),             fCutDCAToVertex2D(reco.fCutDCAToVertex2D)
 165 {
 166   //Copy ctor
 167
 168   for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i]      = reco.fMisalRotShift[i]      ;
 169                                    fMisalTransShift[i]    = reco.fMisalTransShift[i]    ; }
 170   for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
 171   for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }
 172
 173 }
 174
 175
 176 //______________________________________________________________________
 177 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
 178 {
 179   //Assignment operator
 180
 181   if(this == &reco)return *this;
 182   ((TNamed *)this)->operator=(reco);
 183
 184   for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i]    = reco.fMisalTransShift[i]    ;
 185                                    fMisalRotShift[i]      = reco.fMisalRotShift[i]      ; }
 186   for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
 187   for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }
 188
 189   fParticleType              = reco.fParticleType;
 190   fPosAlgo                   = reco.fPosAlgo;
 191   fW0                        = reco.fW0;
 192
 193   fNonLinearityFunction      = reco.fNonLinearityFunction;
 194   fNonLinearThreshold        = reco.fNonLinearThreshold;
 195   fSmearClusterEnergy        = reco.fSmearClusterEnergy;
 196
 197   fCellsRecalibrated         = reco.fCellsRecalibrated;
 198   fRecalibration             = reco.fRecalibration;
 199   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
 200
 201   fTimeRecalibration             = reco.fTimeRecalibration;
 202   fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
 203
 204   fUseRunCorrectionFactors   = reco.fUseRunCorrectionFactors;
 205   fRunCorrectionFactorsSet   = reco.fRunCorrectionFactorsSet;
 206
 207   fRemoveBadChannels         = reco.fRemoveBadChannels;
 208   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
 209   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
 210
 211   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
 212   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
 213
 214   fRejectExoticCluster       = reco.fRejectExoticCluster;
 215   fRejectExoticCells         = reco.fRejectExoticCells;
 216   fExoticCellFraction        = reco.fExoticCellFraction;
 217   fExoticCellDiffTime        = reco.fExoticCellDiffTime;
 218   fExoticCellMinAmplitude    = reco.fExoticCellMinAmplitude;
 219
 220   fPIDUtils                  = reco.fPIDUtils;
 221
 222   fAODFilterMask             = reco.fAODFilterMask;
 223
 224   fCutEtaPhiSum              = reco.fCutEtaPhiSum;
 225   fCutEtaPhiSeparate         = reco.fCutEtaPhiSeparate;
 226   fCutR                      = reco.fCutR;
 227   fCutEta                    = reco.fCutEta;
 228   fCutPhi                    = reco.fCutPhi;
 229   fClusterWindow             = reco.fClusterWindow;
 230   fMass                      = reco.fMass;
 231   fStepSurface               = reco.fStepSurface;
 232   fStepCluster               = reco.fStepCluster;
 233
 234   fTrackCutsType             = reco.fTrackCutsType;
 235   fCutMinTrackPt             = reco.fCutMinTrackPt;
 236   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
 237   fCutMinNClusterITS         = reco.fCutMinNClusterITS;
 238   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
 239   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
 240   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
 241   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
 242   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
 243   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
 244   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
 245   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
 246
 247   if(reco.fResidualEta){
 248     // assign or copy construct
 249     if(fResidualEta){
 250       *fResidualEta = *reco.fResidualEta;
 251     }
 252     else fResidualEta = new TArrayF(*reco.fResidualEta);
 253   }
 254   else{
 255     if(fResidualEta)delete fResidualEta;
 256     fResidualEta = 0;
 257   }
 258
 259   if(reco.fResidualPhi){
 260     // assign or copy construct
 261     if(fResidualPhi){
 262       *fResidualPhi = *reco.fResidualPhi;
 263     }
 264     else fResidualPhi = new TArrayF(*reco.fResidualPhi);
 265   }
 266   else{
 267     if(fResidualPhi)delete fResidualPhi;
 268     fResidualPhi = 0;
 269   }
 270
 271   if(reco.fMatchedTrackIndex){
 272     // assign or copy construct
 273     if(fMatchedTrackIndex){
 274       *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
 275     }
 276     else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
 277   }
 278   else{
 279     if(fMatchedTrackIndex)delete fMatchedTrackIndex;
 280     fMatchedTrackIndex = 0;
 281   }
 282
 283   if(reco.fMatchedClusterIndex){
 284     // assign or copy construct
 285     if(fMatchedClusterIndex){
 286       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
 287     }
 288     else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
 289   }
 290   else{
 291     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
 292     fMatchedClusterIndex = 0;
 293   }
 294
 295   return *this;
 296 }
 297
 298
 299 //_____________________________________
 300 AliEMCALRecoUtils::~AliEMCALRecoUtils()
 301 {
 302   //Destructor.
 303
 304   if(fEMCALRecalibrationFactors) {
 305     fEMCALRecalibrationFactors->Clear();
 306     delete  fEMCALRecalibrationFactors;
 307   }
 308
 309   if(fEMCALTimeRecalibrationFactors) {
 310                 fEMCALTimeRecalibrationFactors->Clear();
 311                 delete  fEMCALTimeRecalibrationFactors;
 312         }
 313
 314   if(fEMCALBadChannelMap) {
 315     fEMCALBadChannelMap->Clear();
 316     delete  fEMCALBadChannelMap;
 317   }
 318
 319   delete fMatchedTrackIndex   ;
 320   delete fMatchedClusterIndex ;
 321   delete fResidualEta         ;
 322   delete fResidualPhi         ;
 323   delete fPIDUtils            ;
 324
 325   InitTrackCuts();
 326 }
 327
 328 //_______________________________________________________________________________
 329 Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
 330                                               Float_t  & amp,    Double_t & time,
 331                                               AliVCaloCells* cells)
 332 {
 333   // Reject cell if criteria not passed and calibrate it
 334
 335   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 336
 337   if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE;
 338
 339   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 340   geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
 341   geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 342
 343   // Do not include bad channels found in analysis,
 344   if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) {
 345     return kFALSE;
 346   }
 347
 348   //Recalibrate energy
 349   amp  = cells->GetCellAmplitude(absID);
 350   if(IsRecalibrationOn())
 351     amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 352
 353
 354   // Recalibrate time
 355   time = cells->GetCellTime(absID);
 356
 357   RecalibrateCellTime(absID,bc,time);
 358
 359   return kTRUE;
 360 }
 361
 362 //_______________________________________________________________
 363 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
 364 {
 365         // Given the list of AbsId of the cluster, get the maximum cell and
 366         // check if there are fNCellsFromBorder from the calorimeter border
 367
 368   if(!cluster){
 369     AliInfo("Cluster pointer null!");
 370     return kFALSE;
 371   }
 372
 373   //If the distance to the border is 0 or negative just exit accept all clusters
 374         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
 375
 376   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;
 377   Bool_t shared = kFALSE;
 378   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
 379
 380   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
 381            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
 382
 383         if(absIdMax==-1) return kFALSE;
 384
 385         //Check if the cell is close to the borders:
 386         Bool_t okrow = kFALSE;
 387         Bool_t okcol = kFALSE;
 388
 389   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 390     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
 391                   iSM,ieta,iphi));
 392   }
 393
 394   //Check rows/phi
 395   if(iSM < 10){
 396     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
 397   }
 398   else{
 399     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
 400   }
 401
 402   //Check columns/eta
 403   if(!fNoEMCALBorderAtEta0){
 404     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
 405   }
 406   else{
 407     if(iSM%2==0){
 408       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;
 409     }
 410     else {
 411       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;
 412     }
 413   }//eta 0 not checked
 414
 415   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
 416            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
 417
 418         if (okcol && okrow) {
 419     //printf("Accept\n");
 420     return kTRUE;
 421   }
 422         else  {
 423     //printf("Reject\n");
 424     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
 425     return kFALSE;
 426   }
 427
 428 }
 429
 430
 431 //_________________________________________________________________________________________________________
 432 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, const Int_t nCells){
 433         // Check that in the cluster cells, there is no bad channel of those stored
 434         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 435
 436         if(!fRemoveBadChannels)  return kFALSE;
 437         if(!fEMCALBadChannelMap) return kFALSE;
 438
 439         Int_t icol = -1;
 440         Int_t irow = -1;
 441         Int_t imod = -1;
 442         for(Int_t iCell = 0; iCell<nCells; iCell++){
 443
 444                 //Get the column and row
 445     Int_t iTower = -1, iIphi = -1, iIeta = -1;
 446     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
 447     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
 448     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 449     if(GetEMCALChannelStatus(imod, icol, irow)){
 450       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
 451       return kTRUE;
 452     }
 453
 454         }// cell cluster loop
 455
 456         return kFALSE;
 457
 458 }
 459
 460 //_______________________________________________________________________
 461 Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
 462 {
 463   // Look to cell neighbourhood and reject if it seems exotic
 464   // Do before recalibrating the cells
 465   if(!fRejectExoticCells) return kFALSE;
 466
 467   AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 468
 469   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 470   geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
 471   geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 472
 473   //Get close cells index, energy and time, not in corners
 474   Int_t absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
 475   Int_t absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
 476   Int_t absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
 477   Int_t absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
 478
 479   Float_t  ecell  = 0, ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
 480   Double_t tcell  = 0, tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
 481   Bool_t   accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
 482
 483   accept  = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
 484
 485   if(!accept) return kTRUE; // reject this cell
 486
 487   if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
 488
 489   accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
 490   accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
 491   accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
 492   accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
 493
 494   /*
 495     printf("Cell absID %d \n",absID);
 496     printf("\t  accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
 497            accept1,accept2,accept3,accept4);
 498     printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
 499            absID,absID1,absID2,absID3,absID4);
 500     printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
 501            ecell,ecell1,ecell2,ecell3,ecell4);
 502     printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
 503            tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
 504            TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
 505   */
 506
 507   if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
 508   if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
 509   if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
 510   if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
 511
 512   Float_t eCross = ecell1+ecell2+ecell3+ecell4;
 513
 514   //printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell);
 515
 516   if(1-eCross/ecell > fExoticCellFraction) {
 517     AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",absID,ecell,eCross,1-eCross/ecell));
 518     return kTRUE;
 519   }
 520
 521   return kFALSE;
 522
 523 }
 524
 525 //_________________________________________________
 526 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster, AliVCaloCells *cells, const Int_t bc) {
 527   // Check if the cluster highest energy tower is exotic
 528
 529   if(!cluster){
 530     AliInfo("Cluster pointer null!");
 531     return kFALSE;
 532   }
 533
 534   if(!fRejectExoticCluster) return kFALSE;
 535
 536   // Get highest energy tower
 537   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 538   Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
 539   Bool_t shared = kFALSE;
 540   GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
 541
 542   return IsExoticCell(absId,cells,bc);
 543
 544 }
 545
 546 //__________________________________________________
 547 Float_t AliEMCALRecoUtils::SmearClusterEnergy(AliVCluster* cluster) {
 548
 549   //In case of MC analysis, smear energy to match resolution/calibration in real data
 550
 551   if(!cluster){
 552     AliInfo("Cluster pointer null!");
 553     return 0;
 554   }
 555
 556   Float_t energy    = cluster->E() ;
 557   Float_t rdmEnergy = energy ;
 558   if(fSmearClusterEnergy){
 559     rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
 560                                     fSmearClusterParam[1] * energy +
 561                                     fSmearClusterParam[2] );
 562     AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
 563   }
 564
 565   return rdmEnergy ;
 566
 567 }
 568
 569 //__________________________________________________
 570 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
 571 // Correct cluster energy from non linearity functions
 572
 573   if(!cluster){
 574     AliInfo("Cluster pointer null!");
 575     return 0;
 576   }
 577
 578   Float_t energy = cluster->E();
 579
 580   switch (fNonLinearityFunction) {
 581
 582     case kPi0MC:
 583     {
 584       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
 585       //Double_t fNonLinearityParams[0] = 1.014;
 586       //Double_t fNonLinearityParams[1] = -0.03329;
 587       //Double_t fNonLinearityParams[2] = -0.3853;
 588       //Double_t fNonLinearityParams[3] = 0.5423;
 589       //Double_t fNonLinearityParams[4] = -0.4335;
 590        energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 591                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 592                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 593       break;
 594     }
 595
 596     case kPi0GammaGamma:
 597     {
 598       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 599       //Double_t fNonLinearityParams[0] = 1.04;
 600       //Double_t fNonLinearityParams[1] = -0.1445;
 601       //Double_t fNonLinearityParams[2] = 1.046;
 602       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 603       break;
 604     }
 605
 606     case kPi0GammaConversion:
 607     {
 608       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 609       //fNonLinearityParams[0] = 0.139393/0.1349766;
 610       //fNonLinearityParams[1] = 0.0566186;
 611       //fNonLinearityParams[2] = 0.982133;
 612       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 613
 614       break;
 615     }
 616
 617     case kBeamTest:
 618     {
 619       //From beam test, Alexei's results, for different ZS thresholds
 620       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
 621       //fNonLinearityParams[0] = 1.007;      1.003;      1.002
 622       //fNonLinearityParams[1] = 0.894;      0.719;      0.797
 623       //fNonLinearityParams[2] = 0.246;      0.334;      0.358
 624       //Rescale the param[0] with 1.03
 625       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
 626
 627       break;
 628     }
 629
 630     case kBeamTestCorrected:
 631     {
 632       //From beam test, corrected for material between beam and EMCAL
 633       //fNonLinearityParams[0] =  0.99078
 634       //fNonLinearityParams[1] =  0.161499;
 635       //fNonLinearityParams[2] =  0.655166;
 636       //fNonLinearityParams[3] =  0.134101;
 637       //fNonLinearityParams[4] =  163.282;
 638       //fNonLinearityParams[5] =  23.6904;
 639       //fNonLinearityParams[6] =  0.978;
 640         energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
 641
 642       break;
 643     }
 644
 645     case kNoCorrection:
 646       AliDebug(2,"No correction on the energy\n");
 647       break;
 648
 649   }
 650
 651   return energy;
 652
 653 }
 654 //__________________________________________________
 655 void AliEMCALRecoUtils::InitNonLinearityParam()
 656 {
 657         //Initialising Non Linearity Parameters
 658
 659         if(fNonLinearityFunction == kPi0MC)
 660                 {
 661                         fNonLinearityParams[0] = 1.014;
 662       fNonLinearityParams[1] = -0.03329;
 663       fNonLinearityParams[2] = -0.3853;
 664       fNonLinearityParams[3] = 0.5423;
 665       fNonLinearityParams[4] = -0.4335;
 666                 }
 667
 668         if(fNonLinearityFunction == kPi0GammaGamma)
 669                 {
 670                         fNonLinearityParams[0] = 1.04;
 671                         fNonLinearityParams[1] = -0.1445;
 672                         fNonLinearityParams[2] = 1.046;
 673                 }
 674
 675         if(fNonLinearityFunction == kPi0GammaConversion)
 676                 {
 677       fNonLinearityParams[0] = 0.139393;
 678       fNonLinearityParams[1] = 0.0566186;
 679       fNonLinearityParams[2] = 0.982133;
 680                 }
 681
 682         if(fNonLinearityFunction == kBeamTest)
 683                 {
 684                         if(fNonLinearThreshold == 30)
 685                         {
 686                                 fNonLinearityParams[0] = 1.007;
 687                                 fNonLinearityParams[1] = 0.894;
 688                                 fNonLinearityParams[2] = 0.246;
 689                         }
 690                         if(fNonLinearThreshold == 45)
 691                         {
 692                                 fNonLinearityParams[0] = 1.003;
 693                                 fNonLinearityParams[1] = 0.719;
 694                                 fNonLinearityParams[2] = 0.334;
 695                         }
 696                         if(fNonLinearThreshold == 75)
 697                         {
 698                                 fNonLinearityParams[0] = 1.002;
 699                                 fNonLinearityParams[1] = 0.797;
 700                                 fNonLinearityParams[2] = 0.358;
 701                         }
 702                 }
 703
 704         if(fNonLinearityFunction == kBeamTestCorrected)
 705                 {
 706                         fNonLinearityParams[0] =  0.99078;
 707                         fNonLinearityParams[1] =  0.161499;
 708                         fNonLinearityParams[2] =  0.655166;
 709                         fNonLinearityParams[3] =  0.134101;
 710                         fNonLinearityParams[4] =  163.282;
 711                         fNonLinearityParams[5] =  23.6904;
 712                         fNonLinearityParams[6] =  0.978;
 713                 }
 714 }
 715
 716 //__________________________________________________
 717 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
 718 {
 719   //Calculate shower depth for a given cluster energy and particle type
 720
 721   // parameters
 722   Float_t x0    = 1.31;
 723   Float_t ecr   = 8;
 724   Float_t depth = 0;
 725
 726   switch ( iParticle )
 727   {
 728     case kPhoton:
 729       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 730       break;
 731
 732     case kElectron:
 733       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
 734       break;
 735
 736     case kHadron:
 737       // hadron
 738       // boxes anc. here
 739       if(gGeoManager){
 740         gGeoManager->cd("ALIC_1/XEN1_1");
 741         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
 742         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
 743         if(geoSM){
 744           TGeoVolume      *geoSMVol   = geoSM->GetVolume();
 745           TGeoShape       *geoSMShape = geoSMVol->GetShape();
 746           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
 747           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
 748           else AliFatal("Null GEANT box");
 749         }else AliFatal("NULL  GEANT node matrix");
 750       }
 751       else{//electron
 752         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
 753       }
 754
 755       break;
 756
 757     default://photon
 758       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 759   }
 760
 761   return depth;
 762
 763 }
 764
 765 //__________________________________________________
 766 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
 767                                          Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
 768 {
 769   //For a given CaloCluster gets the absId of the cell
 770   //with maximum energy deposit.
 771
 772   Double_t eMax        = -1.;
 773   Double_t eCell       = -1.;
 774   Float_t  fraction    = 1.;
 775   Float_t  recalFactor = 1.;
 776   Int_t    cellAbsId   = -1 ;
 777
 778   Int_t iTower  = -1;
 779   Int_t iIphi   = -1;
 780   Int_t iIeta   = -1;
 781   Int_t iSupMod0= -1;
 782
 783   if(!clu){
 784     AliInfo("Cluster pointer null!");
 785     absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
 786     return;
 787   }
 788
 789   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 790     cellAbsId = clu->GetCellAbsId(iDig);
 791     fraction  = clu->GetCellAmplitudeFraction(iDig);
 792     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
 793     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 794     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 795     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 796     if(iDig==0) iSupMod0=iSupMod;
 797     else if(iSupMod0!=iSupMod) {
 798       shared = kTRUE;
 799       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
 800     }
 801     if(IsRecalibrationOn()) {
 802       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 803     }
 804     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 805     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
 806     if(eCell > eMax)  {
 807       eMax  = eCell;
 808       absId = cellAbsId;
 809       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 810     }
 811   }// cell loop
 812
 813   //Get from the absid the supermodule, tower and eta/phi numbers
 814   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 815   //Gives SuperModule and Tower numbers
 816   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 817                                          iIphi, iIeta,iphi,ieta);
 818   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
 819   //printf("Max end---\n");
 820
 821 }
 822
 823 //________________________________________________________________
 824 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
 825         //Init EMCAL recalibration factors
 826         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
 827         //In order to avoid rewriting the same histograms
 828         Bool_t oldStatus = TH1::AddDirectoryStatus();
 829         TH1::AddDirectory(kFALSE);
 830
 831         fEMCALRecalibrationFactors = new TObjArray(10);
 832         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));
 833         //Init the histograms with 1
 834         for (Int_t sm = 0; sm < 10; sm++) {
 835                 for (Int_t i = 0; i < 48; i++) {
 836                         for (Int_t j = 0; j < 24; j++) {
 837                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
 838                         }
 839                 }
 840         }
 841         fEMCALRecalibrationFactors->SetOwner(kTRUE);
 842         fEMCALRecalibrationFactors->Compress();
 843
 844         //In order to avoid rewriting the same histograms
 845         TH1::AddDirectory(oldStatus);
 846 }
 847
 848 //________________________________________________________________
 849 void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors(){
 850         //Init EMCAL recalibration factors
 851         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
 852         //In order to avoid rewriting the same histograms
 853         Bool_t oldStatus = TH1::AddDirectoryStatus();
 854         TH1::AddDirectory(kFALSE);
 855
 856         fEMCALTimeRecalibrationFactors = new TObjArray(4);
 857         for (int i = 0; i < 4; i++)
 858     fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
 859                                                  Form("hAllTimeAvBC%d",i),
 860                                                  48*24*10,0.,48*24*10)          );
 861         //Init the histograms with 1
 862         for (Int_t bc = 0; bc < 4; bc++) {
 863                         for (Int_t i = 0; i < 48*24*10; i++)
 864         SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
 865   }
 866
 867         fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
 868         fEMCALTimeRecalibrationFactors->Compress();
 869
 870         //In order to avoid rewriting the same histograms
 871         TH1::AddDirectory(oldStatus);
 872 }
 873
 874 //________________________________________________________________
 875 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 876         //Init EMCAL bad channels map
 877         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
 878         //In order to avoid rewriting the same histograms
 879         Bool_t oldStatus = TH1::AddDirectoryStatus();
 880         TH1::AddDirectory(kFALSE);
 881
 882         fEMCALBadChannelMap = new TObjArray(10);
 883         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
 884         for (int i = 0; i < 10; i++) {
 885                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
 886         }
 887
 888         fEMCALBadChannelMap->SetOwner(kTRUE);
 889         fEMCALBadChannelMap->Compress();
 890
 891         //In order to avoid rewriting the same histograms
 892         TH1::AddDirectory(oldStatus);
 893 }
 894
 895 //________________________________________________________________
 896 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells, const Int_t bc){
 897         // Recalibrate the cluster energy and Time, considering the recalibration map
 898   // and the energy of the cells and time that compose the cluster.
 899   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
 900
 901   if(!cluster){
 902     AliInfo("Cluster pointer null!");
 903     return;
 904   }
 905
 906         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 907         UShort_t * index    = cluster->GetCellsAbsId() ;
 908         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 909         Int_t ncells = cluster->GetNCells();
 910
 911         //Initialize some used variables
 912         Float_t energy = 0;
 913         Int_t   absId  =-1;
 914   Int_t   icol   =-1, irow =-1, imod=1;
 915         Float_t factor = 1, frac = 0;
 916   Int_t   absIdMax = -1;
 917   Float_t emax     = 0;
 918
 919         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 920         for(Int_t icell = 0; icell < ncells; icell++){
 921                 absId = index[icell];
 922                 frac =  fraction[icell];
 923                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
 924
 925     if(!fCellsRecalibrated && IsRecalibrationOn()){
 926
 927       // Energy
 928       Int_t iTower = -1, iIphi = -1, iIeta = -1;
 929       geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 930       if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
 931       geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 932       factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
 933
 934       AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
 935                       imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
 936
 937     }
 938
 939                 energy += cells->GetCellAmplitude(absId)*factor*frac;
 940
 941     if(emax < cells->GetCellAmplitude(absId)*factor*frac){
 942       emax     = cells->GetCellAmplitude(absId)*factor*frac;
 943       absIdMax = absId;
 944     }
 945
 946         }
 947
 948   cluster->SetE(energy);
 949
 950   AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
 951
 952         // Recalculate time of cluster only for ESDs
 953   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))){
 954
 955     // Time
 956     Double_t weightedTime = 0;
 957     Double_t weight       = 0;
 958     Double_t weightTot    = 0;
 959     Double_t maxcellTime  = 0;
 960     for(Int_t icell = 0; icell < ncells; icell++){
 961       absId = index[icell];
 962       frac =  fraction[icell];
 963       if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
 964
 965       Double_t celltime = cells->GetCellTime(absId);
 966       RecalibrateCellTime(absId, bc, celltime);
 967       if(absId == absIdMax) maxcellTime = celltime;
 968
 969       if(!fCellsRecalibrated){
 970
 971         Int_t iTower = -1, iIphi = -1, iIeta = -1;
 972         geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 973         if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
 974         geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 975         factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
 976
 977         AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
 978                         imod,icol,irow,frac,factor,cells->GetCellTime(absId)));
 979
 980       }
 981
 982       weight        = GetCellWeight(cells->GetCellAmplitude(absId)*factor*frac , energy );
 983       weightTot    += weight;
 984       weightedTime += celltime * weight;
 985
 986     }
 987
 988     if(weightTot > 0)
 989       cluster->SetTOF(weightedTime/weightTot);
 990     else
 991       cluster->SetTOF(maxcellTime);
 992
 993   }
 994 }
 995
 996 //________________________________________________________________
 997 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc){
 998         // Recalibrate the cells time and energy, considering the recalibration map and the energy
 999   // of the cells that compose the cluster.
1000   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1001
1002   if(!IsRecalibrationOn() && !IsTimeRecalibrationOn()) return;
1003
1004   if(!cells){
1005     AliInfo("Cells pointer null!");
1006     return;
1007   }
1008
1009   fCellsRecalibrated = kTRUE;
1010
1011   Int_t    absId  =-1;
1012   Bool_t   accept = kFALSE;
1013   Float_t  ecell  = 0;
1014   Double_t tcell  = 0;
1015
1016   Int_t    nEMcell  = cells->GetNumberOfCells() ;
1017   for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
1018
1019     absId  = cells->GetCellNumber(iCell);
1020
1021     accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1022     if(!accept) {
1023       ecell = 0;
1024       tcell = 0;
1025     }
1026
1027     //Set new values
1028     cells->SetCell(iCell,absId,ecell, tcell);
1029
1030   }
1031
1032 }
1033
1034 //_________________________________________________________________________________________________
1035 void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime)
1036 {
1037         // Recalibrate time of cell with absID  considering the recalibration map
1038   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1039
1040   if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0){
1041
1042     celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;    ;
1043
1044   }
1045
1046 }
1047
1048 //__________________________________________________
1049 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
1050 {
1051   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1052
1053   if(!clu){
1054     AliInfo("Cluster pointer null!");
1055     return;
1056   }
1057
1058   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1059   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1060   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1061
1062 }
1063
1064 //__________________________________________________
1065 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
1066 {
1067   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1068   // The algorithm is a copy of what is done in AliEMCALRecPoint
1069
1070   Double_t eCell       = 0.;
1071   Float_t  fraction    = 1.;
1072   Float_t  recalFactor = 1.;
1073
1074   Int_t    absId   = -1;
1075   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
1076   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
1077   Float_t  weight = 0.,  totalWeight=0.;
1078   Float_t  newPos[3] = {0,0,0};
1079   Double_t pLocal[3], pGlobal[3];
1080   Bool_t shared = kFALSE;
1081
1082   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
1083   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
1084   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1085
1086   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1087
1088   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
1089
1090       absId = clu->GetCellAbsId(iDig);
1091       fraction  = clu->GetCellAmplitudeFraction(iDig);
1092       if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1093
1094     if(!fCellsRecalibrated){
1095
1096       geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1097       geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1098
1099       if(IsRecalibrationOn()) {
1100         recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1101       }
1102     }
1103
1104     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1105
1106     weight = GetCellWeight(eCell,clEnergy);
1107     totalWeight += weight;
1108
1109     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1110     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1111     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1112     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1113
1114     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1115
1116   }// cell loop
1117
1118   if(totalWeight>0){
1119     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
1120   }
1121
1122   //Float_t pos[]={0,0,0};
1123   //clu->GetPosition(pos);
1124   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1125   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1126
1127         if(iSupModMax > 1) {//sector 1
1128           newPos[0] +=fMisalTransShift[3];//-=3.093;
1129           newPos[1] +=fMisalTransShift[4];//+=6.82;
1130           newPos[2] +=fMisalTransShift[5];//+=1.635;
1131     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1132
1133         }
1134         else {//sector 0
1135           newPos[0] +=fMisalTransShift[0];//+=1.134;
1136           newPos[1] +=fMisalTransShift[1];//+=8.2;
1137           newPos[2] +=fMisalTransShift[2];//+=1.197;
1138     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1139
1140         }
1141   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1142
1143   clu->SetPosition(newPos);
1144
1145 }
1146
1147 //__________________________________________________
1148 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
1149 {
1150   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1151   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1152
1153   Double_t eCell       = 1.;
1154   Float_t  fraction    = 1.;
1155   Float_t  recalFactor = 1.;
1156
1157   Int_t absId   = -1;
1158   Int_t iTower  = -1;
1159   Int_t iIphi   = -1, iIeta   = -1;
1160         Int_t iSupMod = -1, iSupModMax = -1;
1161   Int_t iphi = -1, ieta =-1;
1162   Bool_t shared = kFALSE;
1163
1164   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
1165   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
1166   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1167
1168   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
1169   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1170   Int_t startingSM = -1;
1171
1172   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
1173     absId = clu->GetCellAbsId(iDig);
1174     fraction  = clu->GetCellAmplitudeFraction(iDig);
1175     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1176
1177     if     (iDig==0)  startingSM = iSupMod;
1178     else if(iSupMod != startingSM) areInSameSM = kFALSE;
1179
1180     eCell  = cells->GetCellAmplitude(absId);
1181
1182     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1183     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1184
1185     if(!fCellsRecalibrated){
1186
1187       if(IsRecalibrationOn()) {
1188
1189         recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1190
1191       }
1192     }
1193
1194     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1195
1196     weight = GetCellWeight(eCell,clEnergy);
1197     if(weight < 0) weight = 0;
1198     totalWeight += weight;
1199     weightedCol += ieta*weight;
1200     weightedRow += iphi*weight;
1201
1202     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1203
1204     }// cell loop
1205
1206   Float_t xyzNew[]={0.,0.,0.};
1207   if(areInSameSM == kTRUE) {
1208     //printf("In Same SM\n");
1209     weightedCol = weightedCol/totalWeight;
1210     weightedRow = weightedRow/totalWeight;
1211     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1212   }
1213   else {
1214     //printf("In Different SM\n");
1215     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1216   }
1217
1218   clu->SetPosition(xyzNew);
1219
1220 }
1221
1222 //____________________________________________________________________________
1223 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
1224
1225   //re-evaluate distance to bad channel with updated bad map
1226
1227   if(!fRecalDistToBadChannels) return;
1228
1229   if(!cluster){
1230     AliInfo("Cluster pointer null!");
1231     return;
1232   }
1233
1234         //Get channels map of the supermodule where the cluster is.
1235   Int_t absIdMax        = -1, iSupMod =-1, icolM = -1, irowM = -1;
1236   Bool_t shared = kFALSE;
1237   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
1238   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1239
1240   Int_t dRrow, dRcol;
1241         Float_t  minDist = 10000.;
1242         Float_t  dist    = 0.;
1243
1244   //Loop on tower status map
1245         for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1246                 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1247                         //Check if tower is bad.
1248                         if(hMap->GetBinContent(icol,irow)==0) continue;
1249       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
1250       //       iSupMod,icol, irow, icolM,irowM);
1251
1252       dRrow=TMath::Abs(irowM-irow);
1253       dRcol=TMath::Abs(icolM-icol);
1254       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1255                         if(dist < minDist){
1256         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1257         minDist = dist;
1258       }
1259
1260                 }
1261         }
1262
1263         //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
1264         if (shared) {
1265                 TH2D* hMap2 = 0;
1266                 Int_t iSupMod2 = -1;
1267
1268                 //The only possible combinations are (0,1), (2,3) ... (8,9)
1269                 if(iSupMod%2) iSupMod2 = iSupMod-1;
1270                 else          iSupMod2 = iSupMod+1;
1271                 hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
1272
1273                 //Loop on tower status map of second super module
1274                 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1275                         for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1276                                 //Check if tower is bad.
1277                                 if(hMap2->GetBinContent(icol,irow)==0) continue;
1278                                 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
1279           //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
1280
1281         dRrow=TMath::Abs(irow-irowM);
1282
1283         if(iSupMod%2) {
1284                                   dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1285                                 }
1286         else {
1287           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1288                                 }
1289
1290                                 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1291         if(dist < minDist) minDist = dist;
1292
1293                         }
1294                 }
1295
1296         }// shared cluster in 2 SuperModules
1297
1298   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1299   cluster->SetDistanceToBadChannel(minDist);
1300
1301 }
1302
1303 //____________________________________________________________________________
1304 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
1305
1306   //re-evaluate identification parameters with bayesian
1307
1308   if(!cluster){
1309     AliInfo("Cluster pointer null!");
1310     return;
1311   }
1312
1313         if ( cluster->GetM02() != 0)
1314     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1315
1316   Float_t pidlist[AliPID::kSPECIESN+1];
1317         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1318
1319   cluster->SetPID(pidlist);
1320
1321 }
1322
1323 //____________________________________________________________________________
1324 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
1325 {
1326   // Calculates new center of gravity in the local EMCAL-module coordinates
1327   // and tranfers into global ALICE coordinates
1328   // Calculates Dispersion and main axis
1329
1330   if(!cluster){
1331     AliInfo("Cluster pointer null!");
1332     return;
1333   }
1334
1335   Int_t nstat  = 0;
1336   Float_t wtot = 0. ;
1337   Double_t eCell       = 0.;
1338   Float_t  fraction    = 1.;
1339   Float_t  recalFactor = 1.;
1340
1341   Int_t iSupMod = -1;
1342   Int_t iTower  = -1;
1343   Int_t iIphi   = -1;
1344   Int_t iIeta   = -1;
1345   Int_t iphi    = -1;
1346   Int_t ieta    = -1;
1347   Double_t etai = -1.;
1348   Double_t phii = -1.;
1349
1350   Double_t w     = 0.;
1351   Double_t d     = 0.;
1352   Double_t dxx   = 0.;
1353   Double_t dzz   = 0.;
1354   Double_t dxz   = 0.;
1355   Double_t xmean = 0.;
1356   Double_t zmean = 0.;
1357
1358   //Loop on cells
1359   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1360
1361     //Get from the absid the supermodule, tower and eta/phi numbers
1362     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1363     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1364
1365     //Get the cell energy, if recalibration is on, apply factors
1366     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1367     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1368
1369     if(!fCellsRecalibrated){
1370
1371       if(IsRecalibrationOn()) {
1372         recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1373       }
1374
1375     }
1376
1377     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1378
1379     if(cluster->E() > 0 && eCell > 0){
1380
1381       w  = GetCellWeight(eCell,cluster->E());
1382
1383       etai=(Double_t)ieta;
1384       phii=(Double_t)iphi;
1385       if(w > 0.0) {
1386         wtot += w ;
1387         nstat++;
1388         //Shower shape
1389         dxx  += w * etai * etai ;
1390         xmean+= w * etai ;
1391         dzz  += w * phii * phii ;
1392         zmean+= w * phii ;
1393         dxz  += w * etai * phii ;
1394       }
1395     }
1396     else
1397       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1398   }//cell loop
1399
1400   //Normalize to the weight
1401   if (wtot > 0) {
1402     xmean /= wtot ;
1403     zmean /= wtot ;
1404   }
1405   else
1406     AliError(Form("Wrong weight %f\n", wtot));
1407
1408   //Calculate dispersion
1409   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1410
1411     //Get from the absid the supermodule, tower and eta/phi numbers
1412     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1413     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1414
1415     //Get the cell energy, if recalibration is on, apply factors
1416     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1417     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1418     if(IsRecalibrationOn()) {
1419       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1420     }
1421     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1422
1423     if(cluster->E() > 0 && eCell > 0){
1424
1425       w  = GetCellWeight(eCell,cluster->E());
1426
1427       etai=(Double_t)ieta;
1428       phii=(Double_t)iphi;
1429       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1430     }
1431     else
1432       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1433   }// cell loop
1434
1435   //Normalize to the weigth and set shower shape parameters
1436   if (wtot > 0 && nstat > 1) {
1437     d /= wtot ;
1438     dxx /= wtot ;
1439     dzz /= wtot ;
1440     dxz /= wtot ;
1441     dxx -= xmean * xmean ;
1442     dzz -= zmean * zmean ;
1443     dxz -= xmean * zmean ;
1444     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1445     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1446   }
1447   else{
1448     d=0. ;
1449     cluster->SetM20(0.) ;
1450     cluster->SetM02(0.) ;
1451   }
1452
1453   if (d>=0)
1454     cluster->SetDispersion(TMath::Sqrt(d)) ;
1455   else
1456     cluster->SetDispersion(0) ;
1457 }
1458
1459 //____________________________________________________________________________
1460 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  AliEMCALGeometry *geom)
1461 {
1462   //This function should be called before the cluster loop
1463   //Before call this function, please recalculate the cluster positions
1464   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1465   //Store matched cluster indexes and residuals
1466
1467   fMatchedTrackIndex->Reset();
1468   fMatchedClusterIndex->Reset();
1469   fResidualPhi->Reset();
1470   fResidualEta->Reset();
1471
1472   fMatchedTrackIndex->Set(500);
1473   fMatchedClusterIndex->Set(500);
1474   fResidualPhi->Set(500);
1475   fResidualEta->Set(500);
1476
1477   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1478   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1479
1480   TObjArray *clusterArray = 0x0;
1481   if(!clusterArr)
1482     {
1483       clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
1484       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1485         {
1486           AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1487           if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1488           clusterArray->AddAt(cluster,icl);
1489         }
1490     }
1491
1492   Int_t    matched=0;
1493   Double_t cv[21];
1494   for (Int_t i=0; i<21;i++) cv[i]=0;
1495   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1496   {
1497     AliExternalTrackParam *trackParam = 0;
1498
1499     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1500     AliESDtrack *esdTrack = 0;
1501     AliAODTrack *aodTrack = 0;
1502     if(esdevent)
1503     {
1504       esdTrack = esdevent->GetTrack(itr);
1505       if(!esdTrack) continue;
1506       if(!IsAccepted(esdTrack)) continue;
1507       if(esdTrack->Pt()<fCutMinTrackPt) continue;
1508       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
1509       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
1510       trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1511     }
1512
1513     //If the input event is AOD, the starting point for extrapolation is at vertex
1514     //AOD tracks are selected according to its filterbit.
1515     else if(aodevent)
1516     {
1517       aodTrack = aodevent->GetTrack(itr);
1518       if(!aodTrack) continue;
1519       if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1520       if(aodTrack->Pt()<fCutMinTrackPt) continue;
1521       Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
1522       if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
1523       Double_t pos[3],mom[3];
1524       aodTrack->GetXYZ(pos);
1525       aodTrack->GetPxPyPz(mom);
1526       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()));
1527       trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1528     }
1529
1530     //Return if the input data is not "AOD" or "ESD"
1531     else
1532     {
1533       printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1534       if(clusterArray)
1535         {
1536           clusterArray->Clear();
1537           delete clusterArray;
1538         }
1539       return;
1540     }
1541
1542     if(!trackParam) continue;
1543
1544     //Extrapolate the track to EMCal surface
1545     AliExternalTrackParam emcalParam(*trackParam);
1546     Float_t eta, phi;
1547     if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi))
1548       {
1549         if(aodevent && trackParam) delete trackParam;
1550         continue;
1551       }
1552
1553     if(esdevent)
1554       {
1555         esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
1556       }
1557
1558     if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
1559       {
1560         if(aodevent && trackParam) delete trackParam;
1561         continue;
1562       }
1563
1564
1565     //Find matched clusters
1566     Int_t index = -1;
1567     Float_t dEta = -999, dPhi = -999;
1568     if(!clusterArr)
1569       {
1570         index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
1571       }
1572     else
1573       {
1574         index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1575       }
1576
1577     if(index>-1)
1578     {
1579       fMatchedTrackIndex   ->AddAt(itr,matched);
1580       fMatchedClusterIndex ->AddAt(index,matched);
1581       fResidualEta         ->AddAt(dEta,matched);
1582       fResidualPhi         ->AddAt(dPhi,matched);
1583       matched++;
1584     }
1585     if(aodevent && trackParam) delete trackParam;
1586   }//track loop
1587
1588   if(clusterArray)
1589     {
1590       clusterArray->Clear();
1591       delete clusterArray;
1592     }
1593
1594   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1595
1596   fMatchedTrackIndex   ->Set(matched);
1597   fMatchedClusterIndex ->Set(matched);
1598   fResidualPhi         ->Set(matched);
1599   fResidualEta         ->Set(matched);
1600 }
1601
1602 //________________________________________________________________________________
1603 Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi)
1604 {
1605   //
1606   // This function returns the index of matched cluster to input track
1607   // Returns -1 if no match is found
1608   Int_t index = -1;
1609   Double_t phiV = track->Phi()*TMath::RadToDeg();
1610   if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
1611   AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1612   if(!trackParam) return index;
1613   AliExternalTrackParam emcalParam(*trackParam);
1614   Float_t eta, phi;
1615   if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) return index;
1616   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
1617
1618   TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
1619
1620   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1621   {
1622     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1623     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1624     clusterArr->AddAt(cluster,icl);
1625   }
1626
1627   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1628   clusterArr->Clear();
1629   delete clusterArr;
1630
1631   return index;
1632 }
1633
1634 //________________________________________________________________________________
1635 Int_t  AliEMCALRecoUtils::FindMatchedClusterInClusterArr(AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, TObjArray * clusterArr, Float_t &dEta, Float_t &dPhi)
1636 {
1637   dEta=-999, dPhi=-999;
1638   Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1639   Int_t index = -1;
1640   Float_t tmpEta=-999, tmpPhi=-999;
1641
1642   Double_t exPos[3] = {0.,0.,0.};
1643   if(!emcalParam->GetXYZ(exPos)) return index;
1644
1645   Float_t clsPos[3] = {0.,0.,0.};
1646   for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1647     {
1648       AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1649       if(!cluster || !cluster->IsEMCAL()) continue;
1650       cluster->GetPosition(clsPos);
1651       Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
1652       if(dR > fClusterWindow) continue;
1653
1654       AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
1655       if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
1656       if(fCutEtaPhiSum)
1657         {
1658           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1659           if(tmpR<dRMax)
1660             {
1661               dRMax=tmpR;
1662               dEtaMax=tmpEta;
1663               dPhiMax=tmpPhi;
1664               index=icl;
1665             }
1666         }
1667       else if(fCutEtaPhiSeparate)
1668         {
1669           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1670             {
1671               dEtaMax = tmpEta;
1672               dPhiMax = tmpPhi;
1673               index=icl;
1674             }
1675         }
1676       else
1677         {
1678           printf("Error: please specify your cut criteria\n");
1679           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1680           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1681           return index;
1682         }
1683     }
1684
1685   dEta=dEtaMax;
1686   dPhi=dPhiMax;
1687
1688   return index;
1689 }
1690
1691 //
1692 //------------------------------------------------------------------------------
1693 //
1694 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, Double_t emcalR, Double_t mass, Double_t step, Float_t &eta, Float_t &phi)
1695 {
1696   eta = -999, phi = -999;
1697   if(!trkParam) return kFALSE;
1698   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
1699   Double_t trkPos[3] = {0.,0.,0.};
1700   if(!trkParam->GetXYZ(trkPos)) return kFALSE;
1701   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1702   eta = trkPosVec.Eta();
1703   phi = trkPosVec.Phi();
1704   if(phi<0)
1705     phi += 2*TMath::Pi();
1706
1707   return kTRUE;
1708 }
1709
1710
1711 //
1712 //------------------------------------------------------------------------------
1713 //
1714 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, Float_t *clsPos, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
1715 {
1716   //
1717   //Return the residual by extrapolating a track param to a global position
1718   //
1719   tmpEta = -999;
1720   tmpPhi = -999;
1721   if(!trkParam) return kFALSE;
1722   Double_t trkPos[3] = {0.,0.,0.};
1723   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1724   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1725   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1726   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
1727   if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
1728
1729   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1730   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1731
1732   // track cluster matching
1733   tmpPhi = clsPosVec.DeltaPhi(trkPosVec);    // tmpPhi is between -pi and pi
1734   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
1735
1736   return kTRUE;
1737 }
1738
1739
1740 //
1741 //------------------------------------------------------------------------------
1742 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
1743 {
1744   //
1745   //Return the residual by extrapolating a track param to a cluster
1746   //
1747   tmpEta = -999;
1748   tmpPhi = -999;
1749   if(!cluster || !trkParam) return kFALSE;
1750
1751   Float_t clsPos[3] = {0.,0.,0.};
1752   cluster->GetPosition(clsPos);
1753
1754   return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
1755 }
1756
1757 //
1758 //------------------------------------------------------------------------------
1759 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1760 {
1761   //
1762   //Return the residual by extrapolating a track param to a clusterfStepCluster
1763   //
1764
1765   return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
1766 }
1767
1768
1769 //________________________________________________________________________________
1770 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1771 {
1772   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1773   //Get the residuals dEta and dPhi for this cluster to the closest track
1774   //Works with ESDs and AODs
1775
1776   if( FindMatchedPosForCluster(clsIndex) >= 999 )
1777   {
1778     AliDebug(2,"No matched tracks found!\n");
1779     dEta=999.;
1780     dPhi=999.;
1781     return;
1782   }
1783   dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1784   dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1785 }
1786 //________________________________________________________________________________
1787 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1788 {
1789   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1790   //Get the residuals dEta and dPhi for this track to the closest cluster
1791   //Works with ESDs and AODs
1792
1793   if( FindMatchedPosForTrack(trkIndex) >= 999 )
1794   {
1795     AliDebug(2,"No matched cluster found!\n");
1796     dEta=999.;
1797     dPhi=999.;
1798     return;
1799   }
1800   dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1801   dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1802 }
1803
1804 //__________________________________________________________
1805 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1806 {
1807   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1808   //Get the index of matched track to this cluster
1809   //Works with ESDs and AODs
1810
1811   if(IsClusterMatched(clsIndex))
1812     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1813   else
1814     return -1;
1815 }
1816
1817 //__________________________________________________________
1818 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1819 {
1820   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1821   //Get the index of matched cluster to this track
1822   //Works with ESDs and AODs
1823
1824   if(IsTrackMatched(trkIndex))
1825     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1826   else
1827     return -1;
1828 }
1829
1830 //__________________________________________________
1831 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
1832 {
1833   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1834   //Returns if the cluster has a match
1835   if(FindMatchedPosForCluster(clsIndex) < 999)
1836     return kTRUE;
1837   else
1838     return kFALSE;
1839 }
1840
1841 //__________________________________________________
1842 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
1843 {
1844   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1845   //Returns if the track has a match
1846   if(FindMatchedPosForTrack(trkIndex) < 999)
1847     return kTRUE;
1848   else
1849     return kFALSE;
1850 }
1851
1852 //__________________________________________________________
1853 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1854 {
1855   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1856   //Returns the position of the match in the fMatchedClusterIndex array
1857   Float_t tmpR = fCutR;
1858   UInt_t pos = 999;
1859
1860   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1861   {
1862     if(fMatchedClusterIndex->At(i)==clsIndex)
1863       {
1864         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1865         if(r<tmpR)
1866           {
1867             pos=i;
1868             tmpR=r;
1869             AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1870           }
1871       }
1872   }
1873   return pos;
1874 }
1875
1876 //__________________________________________________________
1877 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1878 {
1879   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1880   //Returns the position of the match in the fMatchedTrackIndex array
1881   Float_t tmpR = fCutR;
1882   UInt_t pos = 999;
1883
1884   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1885   {
1886     if(fMatchedTrackIndex->At(i)==trkIndex)
1887       {
1888         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1889         if(r<tmpR)
1890           {
1891             pos=i;
1892             tmpR=r;
1893             AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1894           }
1895       }
1896   }
1897   return pos;
1898 }
1899
1900 //___________________________________________________________________________________
1901 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom,
1902                                         AliVCaloCells* cells,const Int_t bc)
1903 {
1904   // check if the cluster survives some quality cut
1905   //
1906   //
1907   Bool_t isGood=kTRUE;
1908
1909   if(!cluster || !cluster->IsEMCAL())              return kFALSE;
1910
1911   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1912
1913   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1914
1915   if(IsExoticCluster(cluster, cells,bc))           return kFALSE;
1916
1917   return isGood;
1918 }
1919
1920 //__________________________________________________________
1921 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1922 {
1923   // Given a esd track, return whether the track survive all the cuts
1924
1925   // The different quality parameter are first
1926   // retrieved from the track. then it is found out what cuts the
1927   // track did not survive and finally the cuts are imposed.
1928
1929   UInt_t status = esdTrack->GetStatus();
1930
1931   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1932   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1933
1934   Float_t chi2PerClusterITS = -1;
1935   Float_t chi2PerClusterTPC = -1;
1936   if (nClustersITS!=0)
1937     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1938   if (nClustersTPC!=0)
1939     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1940
1941
1942   //DCA cuts
1943   if(fTrackCutsType==kGlobalCut)
1944     {
1945       Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1946       //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1947       SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1948     }
1949
1950
1951   Float_t b[2];
1952   Float_t bCov[3];
1953   esdTrack->GetImpactParameters(b,bCov);
1954   if (bCov[0]<=0 || bCov[2]<=0) {
1955     AliDebug(1, "Estimated b resolution lower or equal zero!");
1956     bCov[0]=0; bCov[2]=0;
1957   }
1958
1959   Float_t dcaToVertexXY = b[0];
1960   Float_t dcaToVertexZ = b[1];
1961   Float_t dcaToVertex = -1;
1962
1963   if (fCutDCAToVertex2D)
1964     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1965   else
1966     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1967
1968   // cut the track?
1969
1970   Bool_t cuts[kNCuts];
1971   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1972
1973   // track quality cuts
1974   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1975     cuts[0]=kTRUE;
1976   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1977     cuts[1]=kTRUE;
1978   if (nClustersTPC<fCutMinNClusterTPC)
1979     cuts[2]=kTRUE;
1980   if (nClustersITS<fCutMinNClusterITS)
1981     cuts[3]=kTRUE;
1982   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1983     cuts[4]=kTRUE;
1984   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1985     cuts[5]=kTRUE;
1986   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1987     cuts[6]=kTRUE;
1988   if (fCutDCAToVertex2D && dcaToVertex > 1)
1989     cuts[7] = kTRUE;
1990   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1991     cuts[8] = kTRUE;
1992   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1993     cuts[9] = kTRUE;
1994
1995   if(fTrackCutsType==kGlobalCut)
1996     {
1997       //Require at least one SPD point + anything else in ITS
1998       if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1999         cuts[10] = kTRUE;
2000     }
2001
2002   Bool_t cut=kFALSE;
2003   for (Int_t i=0; i<kNCuts; i++)
2004     if (cuts[i]) {cut = kTRUE;}
2005
2006     // cut the track
2007   if (cut)
2008     return kFALSE;
2009   else
2010     return kTRUE;
2011 }
2012
2013
2014 //__________________________________________________
2015 void AliEMCALRecoUtils::InitTrackCuts()
2016 {
2017   //Intilize the track cut criteria
2018   //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
2019   //Also you can customize the cuts using the setters
2020
2021   switch (fTrackCutsType)
2022     {
2023     case kTPCOnlyCut:
2024       {
2025         AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2026         //TPC
2027         SetMinNClustersTPC(70);
2028         SetMaxChi2PerClusterTPC(4);
2029         SetAcceptKinkDaughters(kFALSE);
2030         SetRequireTPCRefit(kFALSE);
2031
2032         //ITS
2033         SetRequireITSRefit(kFALSE);
2034         SetMaxDCAToVertexZ(3.2);
2035         SetMaxDCAToVertexXY(2.4);
2036         SetDCAToVertex2D(kTRUE);
2037
2038         break;
2039       }
2040
2041     case kGlobalCut:
2042       {
2043         AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2044         //TPC
2045         SetMinNClustersTPC(70);
2046         SetMaxChi2PerClusterTPC(4);
2047         SetAcceptKinkDaughters(kFALSE);
2048         SetRequireTPCRefit(kTRUE);
2049
2050         //ITS
2051         SetRequireITSRefit(kTRUE);
2052         SetMaxDCAToVertexZ(2);
2053         SetMaxDCAToVertexXY();
2054         SetDCAToVertex2D(kFALSE);
2055
2056         break;
2057       }
2058
2059     case kLooseCut:
2060       {
2061         AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2062         SetMinNClustersTPC(50);
2063         SetAcceptKinkDaughters(kTRUE);
2064
2065         break;
2066       }
2067     }
2068 }
2069
2070 //___________________________________________________
2071 void AliEMCALRecoUtils::Print(const Option_t *) const
2072 {
2073   // Print Parameters
2074
2075   printf("AliEMCALRecoUtils Settings: \n");
2076   printf("Misalignment shifts\n");
2077   for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
2078                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2079                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
2080   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
2081   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
2082
2083   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
2084
2085   printf("Matching criteria: ");
2086   if(fCutEtaPhiSum)
2087     {
2088       printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
2089     }
2090   else if(fCutEtaPhiSeparate)
2091     {
2092       printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
2093     }
2094   else
2095     {
2096       printf("Error\n");
2097       printf("please specify your cut criteria\n");
2098       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2099       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2100     }
2101
2102   printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
2103   printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
2104
2105   printf("Track cuts: \n");
2106   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
2107   printf("AOD track selection mask: %d\n",fAODFilterMask);
2108   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2109   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2110   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2111   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2112   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
2113
2114 }
2115
2116 //_____________________________________________________________________
2117 void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber){
2118   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
2119   //Do it only once and only if it is requested
2120
2121   if(!fUseRunCorrectionFactors) return;
2122   if(fRunCorrectionFactorsSet)  return;
2123
2124   AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber));
2125
2126   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
2127   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
2128
2129   SwitchOnRecalibration();
2130   for(Int_t ism = 0; ism < 4; ism++){
2131     for(Int_t icol = 0; icol < 48; icol++){
2132       for(Int_t irow = 0; irow < 24; irow++){
2133         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
2134         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
2135         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
2136         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
2137         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
2138         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
2139       }
2140     }
2141   }
2142    fRunCorrectionFactorsSet = kTRUE;
2143 }
2144