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