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