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