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
  29 // --- standard c ---
  30
  31 // standard C++ includes
  32 //#include <Riostream.h>
  33
  34 // ROOT includes
  35 #include <TGeoManager.h>
  36 #include <TGeoMatrix.h>
  37 #include <TGeoBBox.h>
  38
  39 // STEER includes
  40 #include "AliVCluster.h"
  41 #include "AliVCaloCells.h"
  42 #include "AliVEvent.h"
  43 #include "AliLog.h"
  44 #include "AliPID.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliAODTrack.h"
  49 #include "AliExternalTrackParam.h"
  50 #include "AliESDfriendTrack.h"
  51 #include "AliTrackerBase.h"
  52
  53 // EMCAL includes
  54 #include "AliEMCALRecoUtils.h"
  55 #include "AliEMCALGeometry.h"
  56 #include "AliEMCALTrack.h"
  57 #include "AliEMCALCalibTimeDepCorrection.h"
  58 #include "AliEMCALPIDUtils.h"
  59
  60 ClassImp(AliEMCALRecoUtils)
  61
  62 //______________________________________________
  63 AliEMCALRecoUtils::AliEMCALRecoUtils():
  64   fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
  65   fPosAlgo(kUnchanged),fW0(4.), fNonLinearThreshold(30),
  66   fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
  67   fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
  68   fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
  69   fAODFilterMask(32),
  70   fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
  71   fResidualZ(0x0), fResidualR(0x0), fCutR(10), fCutZ(10), fMass(0.139), fStep(1),
  72   fRejectExoticCluster(kFALSE),
  73   fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
  74   fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
  75   fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
  76   fUseTimeCorrectionFactors(kFALSE),  fTimeCorrectionFactorsSet(kFALSE)
  77 {
  78 //
  79   // Constructor.
  80   // Initialize all constant values which have to be used
  81   // during Reco algorithm execution
  82   //
  83
  84   //Misalignment matrices
  85   for(Int_t i = 0; i < 15 ; i++) {
  86       fMisalTransShift[i] = 0.;
  87       fMisalRotShift[i]   = 0.;
  88   }
  89
  90   //Non linearity
  91   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = 0.;
  92
  93   //For kBeamTestCorrected case, but default is no correction
  94   fNonLinearityParams[0] =  0.99078;
  95   fNonLinearityParams[1] =  0.161499;
  96   fNonLinearityParams[2] =  0.655166;
  97   fNonLinearityParams[3] =  0.134101;
  98   fNonLinearityParams[4] =  163.282;
  99   fNonLinearityParams[5] =  23.6904;
 100   fNonLinearityParams[6] =  0.978;
 101
 102   //For kPi0GammaGamma case
 103   //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
 104   //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
 105   //fNonLinearityParams[2] = 1.046;
 106
 107   //Track matching
 108   fMatchedTrackIndex   = new TArrayI();
 109   fMatchedClusterIndex = new TArrayI();
 110   fResidualZ           = new TArrayF();
 111   fResidualR           = new TArrayF();
 112
 113   InitTrackCuts();
 114
 115   fPIDUtils            = new AliEMCALPIDUtils();
 116
 117
 118 }
 119
 120 //______________________________________________________________________
 121 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
 122 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
 123   fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), fNonLinearThreshold(reco.fNonLinearThreshold),
 124   fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
 125   fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
 126   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
 127   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
 128   fAODFilterMask(reco.fAODFilterMask),
 129   fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
 130   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
 131   fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
 132   fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
 133   fCutR(reco.fCutR),fCutZ(reco.fCutZ),fMass(reco.fMass), fStep(reco.fStep),
 134   fRejectExoticCluster(reco.fRejectExoticCluster),
 135   fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
 136   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
 137   fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
 138   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
 139   fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
 140   fPIDUtils(reco.fPIDUtils),
 141   fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors),  fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
 142 {
 143   //Copy ctor
 144
 145   for(Int_t i = 0; i < 15 ; i++) {
 146       fMisalRotShift[i] = reco.fMisalRotShift[i];
 147       fMisalTransShift[i] = reco.fMisalTransShift[i];
 148   }
 149   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
 150
 151 }
 152
 153
 154 //______________________________________________________________________
 155 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
 156 {
 157   //Assignment operator
 158
 159   if(this == &reco)return *this;
 160   ((TNamed *)this)->operator=(reco);
 161
 162   fNonLinearityFunction      = reco.fNonLinearityFunction;
 163   fParticleType              = reco.fParticleType;
 164   fPosAlgo                   = reco.fPosAlgo;
 165   fW0                        = reco.fW0;
 166   fNonLinearThreshold        = reco.fNonLinearThreshold;
 167   fRecalibration             = reco.fRecalibration;
 168   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
 169   fRemoveBadChannels         = reco.fRemoveBadChannels;
 170   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
 171   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
 172   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
 173   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
 174
 175
 176   for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
 177   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
 178
 179   fAODFilterMask              = reco.fAODFilterMask;
 180
 181   fCutR                      = reco.fCutR;
 182   fCutZ                      = reco.fCutZ;
 183   fMass                      = reco.fMass;
 184   fStep                      = reco.fStep;
 185   fRejectExoticCluster       = reco.fRejectExoticCluster;
 186
 187   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
 188   fCutMinNClusterITS         = reco.fCutMinNClusterITS;
 189   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
 190   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
 191   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
 192   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
 193   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
 194   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
 195   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
 196   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
 197
 198   fPIDUtils                  = reco.fPIDUtils;
 199
 200   fUseTimeCorrectionFactors  = reco.fUseTimeCorrectionFactors;
 201   fTimeCorrectionFactorsSet  = reco.fTimeCorrectionFactorsSet;
 202
 203
 204   if(reco.fResidualR){
 205     // assign or copy construct
 206     if(fResidualR){
 207       *fResidualR = *reco.fResidualR;
 208     }
 209     else fResidualR = new TArrayF(*reco.fResidualR);
 210   }
 211   else{
 212     if(fResidualR)delete fResidualR;
 213     fResidualR = 0;
 214   }
 215
 216   if(reco.fResidualZ){
 217     // assign or copy construct
 218     if(fResidualZ){
 219       *fResidualZ = *reco.fResidualZ;
 220     }
 221     else fResidualZ = new TArrayF(*reco.fResidualZ);
 222   }
 223   else{
 224     if(fResidualZ)delete fResidualZ;
 225     fResidualZ = 0;
 226   }
 227
 228   if(reco.fMatchedTrackIndex){
 229     // assign or copy construct
 230     if(fMatchedTrackIndex){
 231       *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
 232     }
 233     else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
 234   }
 235   else{
 236     if(fMatchedTrackIndex)delete fMatchedTrackIndex;
 237     fMatchedTrackIndex = 0;
 238   }
 239
 240   if(reco.fMatchedClusterIndex){
 241     // assign or copy construct
 242     if(fMatchedClusterIndex){
 243       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
 244     }
 245     else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
 246   }
 247   else{
 248     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
 249     fMatchedClusterIndex = 0;
 250   }
 251
 252
 253   return *this;
 254 }
 255
 256
 257 //__________________________________________________
 258 AliEMCALRecoUtils::~AliEMCALRecoUtils()
 259 {
 260   //Destructor.
 261
 262         if(fEMCALRecalibrationFactors) {
 263                 fEMCALRecalibrationFactors->Clear();
 264                 delete  fEMCALRecalibrationFactors;
 265         }
 266
 267   if(fEMCALBadChannelMap) {
 268                 fEMCALBadChannelMap->Clear();
 269                 delete  fEMCALBadChannelMap;
 270         }
 271
 272   if(fMatchedTrackIndex)   {delete fMatchedTrackIndex;   fMatchedTrackIndex=0;}
 273   if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
 274   if(fResidualR)           {delete fResidualR;           fResidualR=0;}
 275   if(fResidualZ)           {delete fResidualZ;           fResidualZ=0;}
 276
 277 }
 278
 279 //_______________________________________________________________
 280 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
 281 {
 282         // Given the list of AbsId of the cluster, get the maximum cell and
 283         // check if there are fNCellsFromBorder from the calorimeter border
 284
 285   //If the distance to the border is 0 or negative just exit accept all clusters
 286         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
 287
 288   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;
 289   Bool_t shared = kFALSE;
 290   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
 291
 292   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
 293            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
 294
 295         if(absIdMax==-1) return kFALSE;
 296
 297         //Check if the cell is close to the borders:
 298         Bool_t okrow = kFALSE;
 299         Bool_t okcol = kFALSE;
 300
 301   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 302     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
 303                   iSM,ieta,iphi));
 304   }
 305
 306   //Check rows/phi
 307   if(iSM < 10){
 308     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
 309   }
 310   else{
 311     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
 312   }
 313
 314   //Check columns/eta
 315   if(!fNoEMCALBorderAtEta0){
 316     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
 317   }
 318   else{
 319     if(iSM%2==0){
 320       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;
 321     }
 322     else {
 323       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;
 324     }
 325   }//eta 0 not checked
 326
 327   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
 328            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
 329
 330         if (okcol && okrow) {
 331     //printf("Accept\n");
 332     return kTRUE;
 333   }
 334         else  {
 335     //printf("Reject\n");
 336     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
 337     return kFALSE;
 338   }
 339
 340 }
 341
 342
 343 //_________________________________________________________________________________________________________
 344 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
 345         // Check that in the cluster cells, there is no bad channel of those stored
 346         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 347
 348         if(!fRemoveBadChannels)  return kFALSE;
 349         if(!fEMCALBadChannelMap) return kFALSE;
 350
 351         Int_t icol = -1;
 352         Int_t irow = -1;
 353         Int_t imod = -1;
 354         for(Int_t iCell = 0; iCell<nCells; iCell++){
 355
 356                 //Get the column and row
 357     Int_t iTower = -1, iIphi = -1, iIeta = -1;
 358     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
 359     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
 360     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 361     if(GetEMCALChannelStatus(imod, icol, irow)){
 362       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
 363       return kTRUE;
 364     }
 365
 366         }// cell cluster loop
 367
 368         return kFALSE;
 369
 370 }
 371
 372 //_________________________________________________
 373 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster){
 374   // Check if the cluster has high energy  but small number of cells
 375   // The criteria comes from Gustavo's study
 376   //
 377
 378   if(cluster->GetNCells()<(1+cluster->E()/3.))
 379     return kTRUE;
 380   else
 381     return kFALSE;
 382 }
 383
 384 //__________________________________________________
 385 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
 386 // Correct cluster energy from non linearity functions
 387   Float_t energy = cluster->E();
 388
 389   switch (fNonLinearityFunction) {
 390
 391     case kPi0MC:
 392     {
 393       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
 394       //Double_t fNonLinearityParams[0] = 1.014;
 395       //Double_t fNonLinearityParams[1] = -0.03329;
 396       //Double_t fNonLinearityParams[2] = -0.3853;
 397       //Double_t fNonLinearityParams[3] = 0.5423;
 398       //Double_t fNonLinearityParams[4] = -0.4335;
 399        energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 400                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 401                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 402       break;
 403     }
 404
 405     case kPi0GammaGamma:
 406     {
 407       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 408       //Double_t fNonLinearityParams[0] = 1.04;
 409       //Double_t fNonLinearityParams[1] = -0.1445;
 410       //Double_t fNonLinearityParams[2] = 1.046;
 411       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 412       break;
 413     }
 414
 415     case kPi0GammaConversion:
 416     {
 417       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 418       //fNonLinearityParams[0] = 0.139393/0.1349766;
 419       //fNonLinearityParams[1] = 0.0566186;
 420       //fNonLinearityParams[2] = 0.982133;
 421       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 422
 423       break;
 424     }
 425
 426     case kBeamTest:
 427     {
 428       //From beam test, Alexei's results, for different ZS thresholds
 429       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
 430       //fNonLinearityParams[0] = 1.007;      1.003;      1.002
 431       //fNonLinearityParams[1] = 0.894;      0.719;      0.797
 432       //fNonLinearityParams[2] = 0.246;      0.334;      0.358
 433       //Rescale the param[0] with 1.03
 434       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
 435
 436       break;
 437     }
 438
 439     case kBeamTestCorrected:
 440     {
 441       //From beam test, corrected for material between beam and EMCAL
 442       //fNonLinearityParams[0] =  0.99078
 443       //fNonLinearityParams[1] =  0.161499;
 444       //fNonLinearityParams[2] =  0.655166;
 445       //fNonLinearityParams[3] =  0.134101;
 446       //fNonLinearityParams[4] =  163.282;
 447       //fNonLinearityParams[5] =  23.6904;
 448       //fNonLinearityParams[6] =  0.978;
 449         energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
 450
 451       break;
 452     }
 453
 454     case kNoCorrection:
 455       AliDebug(2,"No correction on the energy\n");
 456       break;
 457
 458   }
 459
 460   return energy;
 461
 462 }
 463 //__________________________________________________
 464 void AliEMCALRecoUtils::InitNonLinearityParam()
 465 {
 466         //Initialising Non Linearity Parameters
 467
 468         if(fNonLinearityFunction == kPi0MC)
 469                 {
 470                         fNonLinearityParams[0] = 1.014;
 471       fNonLinearityParams[1] = -0.03329;
 472       fNonLinearityParams[2] = -0.3853;
 473       fNonLinearityParams[3] = 0.5423;
 474       fNonLinearityParams[4] = -0.4335;
 475                 }
 476
 477         if(fNonLinearityFunction == kPi0GammaGamma)
 478                 {
 479                         fNonLinearityParams[0] = 1.04;
 480                         fNonLinearityParams[1] = -0.1445;
 481                         fNonLinearityParams[2] = 1.046;
 482                 }
 483
 484         if(fNonLinearityFunction == kPi0GammaConversion)
 485                 {
 486       fNonLinearityParams[0] = 0.139393;
 487       fNonLinearityParams[1] = 0.0566186;
 488       fNonLinearityParams[2] = 0.982133;
 489                 }
 490
 491         if(fNonLinearityFunction == kBeamTest)
 492                 {
 493                         if(fNonLinearThreshold == 30)
 494                         {
 495                                 fNonLinearityParams[0] = 1.007;
 496                                 fNonLinearityParams[1] = 0.894;
 497                                 fNonLinearityParams[2] = 0.246;
 498                         }
 499                         if(fNonLinearThreshold == 45)
 500                         {
 501                                 fNonLinearityParams[0] = 1.003;
 502                                 fNonLinearityParams[1] = 0.719;
 503                                 fNonLinearityParams[2] = 0.334;
 504                         }
 505                         if(fNonLinearThreshold == 75)
 506                         {
 507                                 fNonLinearityParams[0] = 1.002;
 508                                 fNonLinearityParams[1] = 0.797;
 509                                 fNonLinearityParams[2] = 0.358;
 510                         }
 511                 }
 512
 513         if(fNonLinearityFunction == kBeamTestCorrected)
 514                 {
 515                         fNonLinearityParams[0] =  0.99078;
 516                         fNonLinearityParams[1] =  0.161499;
 517                         fNonLinearityParams[2] =  0.655166;
 518                         fNonLinearityParams[3] =  0.134101;
 519                         fNonLinearityParams[4] =  163.282;
 520                         fNonLinearityParams[5] =  23.6904;
 521                         fNonLinearityParams[6] =  0.978;
 522                 }
 523 }
 524
 525 //__________________________________________________
 526 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
 527 {
 528   //Calculate shower depth for a given cluster energy and particle type
 529
 530   // parameters
 531   Float_t x0    = 1.31;
 532   Float_t ecr   = 8;
 533   Float_t depth = 0;
 534
 535   switch ( iParticle )
 536   {
 537     case kPhoton:
 538       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 539       break;
 540
 541     case kElectron:
 542       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
 543       break;
 544
 545     case kHadron:
 546       // hadron
 547       // boxes anc. here
 548       if(gGeoManager){
 549         gGeoManager->cd("ALIC_1/XEN1_1");
 550         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
 551         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
 552         if(geoSM){
 553           TGeoVolume      *geoSMVol   = geoSM->GetVolume();
 554           TGeoShape       *geoSMShape = geoSMVol->GetShape();
 555           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
 556           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
 557           else AliFatal("Null GEANT box");
 558         }else AliFatal("NULL  GEANT node matrix");
 559       }
 560       else{//electron
 561         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
 562       }
 563
 564       break;
 565
 566     default://photon
 567       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 568   }
 569
 570   return depth;
 571
 572 }
 573
 574 //__________________________________________________
 575 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
 576                                          Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
 577 {
 578   //For a given CaloCluster gets the absId of the cell
 579   //with maximum energy deposit.
 580
 581   Double_t eMax        = -1.;
 582   Double_t eCell       = -1.;
 583   Float_t  fraction    = 1.;
 584   Float_t  recalFactor = 1.;
 585   Int_t    cellAbsId   = -1 ;
 586
 587   Int_t iTower  = -1;
 588   Int_t iIphi   = -1;
 589   Int_t iIeta   = -1;
 590   Int_t iSupMod0= -1;
 591         //printf("---Max?\n");
 592   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 593     cellAbsId = clu->GetCellAbsId(iDig);
 594     fraction  = clu->GetCellAmplitudeFraction(iDig);
 595     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
 596     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 597     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 598     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 599     if(iDig==0) iSupMod0=iSupMod;
 600     else if(iSupMod0!=iSupMod) {
 601       shared = kTRUE;
 602       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
 603     }
 604     if(IsRecalibrationOn()) {
 605       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 606     }
 607     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 608     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
 609     if(eCell > eMax)  {
 610       eMax  = eCell;
 611       absId = cellAbsId;
 612       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 613     }
 614   }// cell loop
 615
 616   //Get from the absid the supermodule, tower and eta/phi numbers
 617   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 618   //Gives SuperModule and Tower numbers
 619   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 620                                          iIphi, iIeta,iphi,ieta);
 621   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
 622   //printf("Max end---\n");
 623
 624 }
 625
 626 //________________________________________________________________
 627 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
 628         //Init EMCAL recalibration factors
 629         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
 630         //In order to avoid rewriting the same histograms
 631         Bool_t oldStatus = TH1::AddDirectoryStatus();
 632         TH1::AddDirectory(kFALSE);
 633
 634         fEMCALRecalibrationFactors = new TObjArray(10);
 635         for (int i = 0; i < 10; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
 636         //Init the histograms with 1
 637         for (Int_t sm = 0; sm < 10; sm++) {
 638                 for (Int_t i = 0; i < 48; i++) {
 639                         for (Int_t j = 0; j < 24; j++) {
 640                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
 641                         }
 642                 }
 643         }
 644         fEMCALRecalibrationFactors->SetOwner(kTRUE);
 645         fEMCALRecalibrationFactors->Compress();
 646
 647         //In order to avoid rewriting the same histograms
 648         TH1::AddDirectory(oldStatus);
 649 }
 650
 651
 652 //________________________________________________________________
 653 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 654         //Init EMCAL bad channels map
 655         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
 656         //In order to avoid rewriting the same histograms
 657         Bool_t oldStatus = TH1::AddDirectoryStatus();
 658         TH1::AddDirectory(kFALSE);
 659
 660         fEMCALBadChannelMap = new TObjArray(10);
 661         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
 662         for (int i = 0; i < 10; i++) {
 663                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
 664         }
 665
 666         //delete hTemp;
 667
 668         fEMCALBadChannelMap->SetOwner(kTRUE);
 669         fEMCALBadChannelMap->Compress();
 670
 671         //In order to avoid rewriting the same histograms
 672         TH1::AddDirectory(oldStatus);
 673 }
 674
 675 //________________________________________________________________
 676 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
 677         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 678
 679         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 680         UShort_t * index    = cluster->GetCellsAbsId() ;
 681         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 682         Int_t ncells = cluster->GetNCells();
 683
 684         //Initialize some used variables
 685         Float_t energy = 0;
 686         Int_t absId    = -1;
 687   Int_t icol = -1, irow = -1, imod=1;
 688         Float_t factor = 1, frac = 0;
 689
 690         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 691         for(Int_t icell = 0; icell < ncells; icell++){
 692                 absId = index[icell];
 693                 frac =  fraction[icell];
 694                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
 695                 Int_t iTower = -1, iIphi = -1, iIeta = -1;
 696                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 697                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
 698                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 699                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
 700     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
 701              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
 702
 703                 energy += cells->GetCellAmplitude(absId)*factor*frac;
 704         }
 705
 706
 707                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
 708
 709         cluster->SetE(energy);
 710
 711 }
 712
 713
 714 //__________________________________________________
 715 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 716 {
 717   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 718
 719   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
 720   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
 721   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
 722
 723 }
 724
 725 //__________________________________________________
 726 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 727 {
 728   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 729   // The algorithm is a copy of what is done in AliEMCALRecPoint
 730
 731   Double_t eCell       = 0.;
 732   Float_t  fraction    = 1.;
 733   Float_t  recalFactor = 1.;
 734
 735   Int_t    absId   = -1;
 736   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
 737   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
 738   Float_t  weight = 0.,  totalWeight=0.;
 739   Float_t  newPos[3] = {0,0,0};
 740   Double_t pLocal[3], pGlobal[3];
 741   Bool_t shared = kFALSE;
 742
 743   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
 744   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
 745   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
 746
 747   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
 748
 749   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 750     absId = clu->GetCellAbsId(iDig);
 751     fraction  = clu->GetCellAmplitudeFraction(iDig);
 752     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 753     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
 754     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 755
 756     if(IsRecalibrationOn()) {
 757       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
 758     }
 759     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 760
 761     weight = GetCellWeight(eCell,clEnergy);
 762     //printf("cell energy %f, weight %f\n",eCell,weight);
 763     totalWeight += weight;
 764     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
 765     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
 766     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
 767     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
 768
 769     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
 770
 771   }// cell loop
 772
 773   if(totalWeight>0){
 774     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
 775   }
 776
 777   //Float_t pos[]={0,0,0};
 778   //clu->GetPosition(pos);
 779   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
 780   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 781
 782         if(iSupModMax > 1) {//sector 1
 783           newPos[0] +=fMisalTransShift[3];//-=3.093;
 784           newPos[1] +=fMisalTransShift[4];//+=6.82;
 785           newPos[2] +=fMisalTransShift[5];//+=1.635;
 786     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
 787
 788         }
 789         else {//sector 0
 790           newPos[0] +=fMisalTransShift[0];//+=1.134;
 791           newPos[1] +=fMisalTransShift[1];//+=8.2;
 792           newPos[2] +=fMisalTransShift[2];//+=1.197;
 793     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
 794
 795         }
 796   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 797
 798   clu->SetPosition(newPos);
 799
 800 }
 801
 802 //__________________________________________________
 803 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 804 {
 805   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 806   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
 807
 808   Double_t eCell       = 1.;
 809   Float_t  fraction    = 1.;
 810   Float_t  recalFactor = 1.;
 811
 812   Int_t absId   = -1;
 813   Int_t iTower  = -1;
 814   Int_t iIphi   = -1, iIeta   = -1;
 815         Int_t iSupMod = -1, iSupModMax = -1;
 816   Int_t iphi = -1, ieta =-1;
 817   Bool_t shared = kFALSE;
 818
 819   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
 820   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
 821   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
 822
 823   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
 824   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
 825   Int_t startingSM = -1;
 826
 827   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 828     absId = clu->GetCellAbsId(iDig);
 829     fraction  = clu->GetCellAmplitudeFraction(iDig);
 830     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 831     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 832     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 833
 834     if     (iDig==0)  startingSM = iSupMod;
 835     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 836
 837     eCell  = cells->GetCellAmplitude(absId);
 838
 839     if(IsRecalibrationOn()) {
 840       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 841     }
 842     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 843
 844     weight = GetCellWeight(eCell,clEnergy);
 845     if(weight < 0) weight = 0;
 846     totalWeight += weight;
 847     weightedCol += ieta*weight;
 848     weightedRow += iphi*weight;
 849
 850     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
 851
 852     }// cell loop
 853
 854   Float_t xyzNew[]={0.,0.,0.};
 855   if(areInSameSM == kTRUE) {
 856     //printf("In Same SM\n");
 857     weightedCol = weightedCol/totalWeight;
 858     weightedRow = weightedRow/totalWeight;
 859     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 860   }
 861   else {
 862     //printf("In Different SM\n");
 863     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 864   }
 865
 866   clu->SetPosition(xyzNew);
 867
 868 }
 869
 870 //____________________________________________________________________________
 871 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
 872
 873   //re-evaluate distance to bad channel with updated bad map
 874
 875   if(!fRecalDistToBadChannels) return;
 876
 877         //Get channels map of the supermodule where the cluster is.
 878   Int_t absIdMax        = -1, iSupMod =-1, icolM = -1, irowM = -1;
 879   Bool_t shared = kFALSE;
 880   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
 881   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
 882
 883   Int_t dRrow, dRcol;
 884         Float_t  minDist = 10000.;
 885         Float_t  dist    = 0.;
 886
 887   //Loop on tower status map
 888         for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
 889                 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
 890                         //Check if tower is bad.
 891                         if(hMap->GetBinContent(icol,irow)==0) continue;
 892       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
 893       //       iSupMod,icol, irow, icolM,irowM);
 894
 895       dRrow=TMath::Abs(irowM-irow);
 896       dRcol=TMath::Abs(icolM-icol);
 897       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
 898                         if(dist < minDist){
 899         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
 900         minDist = dist;
 901       }
 902
 903                 }
 904         }
 905
 906         //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
 907         if (shared) {
 908                 TH2D* hMap2 = 0;
 909                 Int_t iSupMod2 = -1;
 910
 911                 //The only possible combinations are (0,1), (2,3) ... (8,9)
 912                 if(iSupMod%2) iSupMod2 = iSupMod-1;
 913                 else          iSupMod2 = iSupMod+1;
 914                 hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
 915
 916                 //Loop on tower status map of second super module
 917                 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
 918                         for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
 919                                 //Check if tower is bad.
 920                                 if(hMap2->GetBinContent(icol,irow)==0) continue;
 921                                 //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",
 922           //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
 923
 924         dRrow=TMath::Abs(irow-irowM);
 925
 926         if(iSupMod%2) {
 927                                   dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
 928                                 }
 929         else {
 930           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
 931                                 }
 932
 933                                 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
 934         if(dist < minDist) minDist = dist;
 935
 936                         }
 937                 }
 938
 939         }// shared cluster in 2 SuperModules
 940
 941   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
 942   cluster->SetDistanceToBadChannel(minDist);
 943
 944 }
 945
 946 //____________________________________________________________________________
 947 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
 948
 949   //re-evaluate identification parameters with bayesian
 950
 951         if ( cluster->GetM02() != 0)
 952     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
 953
 954   Float_t pidlist[AliPID::kSPECIESN+1];
 955         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
 956
 957   cluster->SetPID(pidlist);
 958
 959 }
 960
 961 //____________________________________________________________________________
 962 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
 963 {
 964   // Calculates new center of gravity in the local EMCAL-module coordinates
 965   // and tranfers into global ALICE coordinates
 966   // Calculates Dispersion and main axis
 967
 968   Int_t nstat  = 0;
 969   Float_t wtot = 0. ;
 970   Double_t eCell       = 0.;
 971   Float_t  fraction    = 1.;
 972   Float_t  recalFactor = 1.;
 973
 974   Int_t iSupMod = -1;
 975   Int_t iTower  = -1;
 976   Int_t iIphi   = -1;
 977   Int_t iIeta   = -1;
 978   Int_t iphi    = -1;
 979   Int_t ieta    = -1;
 980   Double_t etai = -1.;
 981   Double_t phii = -1.;
 982
 983   Double_t w     = 0.;
 984   Double_t d     = 0.;
 985   Double_t dxx   = 0.;
 986   Double_t dzz   = 0.;
 987   Double_t dxz   = 0.;
 988   Double_t xmean = 0.;
 989   Double_t zmean = 0.;
 990
 991   //Loop on cells
 992   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
 993
 994     //Get from the absid the supermodule, tower and eta/phi numbers
 995     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 996     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 997
 998     //Get the cell energy, if recalibration is on, apply factors
 999     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1000     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1001     if(IsRecalibrationOn()) {
1002       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1003     }
1004     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1005
1006     if(cluster->E() > 0 && eCell > 0){
1007
1008       w  = GetCellWeight(eCell,cluster->E());
1009
1010       etai=(Double_t)ieta;
1011       phii=(Double_t)iphi;
1012       if(w > 0.0) {
1013         wtot += w ;
1014         nstat++;
1015         //Shower shape
1016         dxx  += w * etai * etai ;
1017         xmean+= w * etai ;
1018         dzz  += w * phii * phii ;
1019         zmean+= w * phii ;
1020         dxz  += w * etai * phii ;
1021       }
1022     }
1023     else
1024       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1025   }//cell loop
1026
1027   //Normalize to the weight
1028   if (wtot > 0) {
1029     xmean /= wtot ;
1030     zmean /= wtot ;
1031   }
1032   else
1033     AliError(Form("Wrong weight %f\n", wtot));
1034
1035   //Calculate dispersion
1036   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1037
1038     //Get from the absid the supermodule, tower and eta/phi numbers
1039     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1040     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1041
1042     //Get the cell energy, if recalibration is on, apply factors
1043     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1044     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1045     if(IsRecalibrationOn()) {
1046       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1047     }
1048     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1049
1050     if(cluster->E() > 0 && eCell > 0){
1051
1052       w  = GetCellWeight(eCell,cluster->E());
1053
1054       etai=(Double_t)ieta;
1055       phii=(Double_t)iphi;
1056       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1057     }
1058     else
1059       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1060   }// cell loop
1061
1062   //Normalize to the weigth and set shower shape parameters
1063   if (wtot > 0 && nstat > 1) {
1064     d /= wtot ;
1065     dxx /= wtot ;
1066     dzz /= wtot ;
1067     dxz /= wtot ;
1068     dxx -= xmean * xmean ;
1069     dzz -= zmean * zmean ;
1070     dxz -= xmean * zmean ;
1071     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1072     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1073   }
1074   else{
1075     d=0. ;
1076     cluster->SetM20(0.) ;
1077     cluster->SetM02(0.) ;
1078   }
1079
1080   if (d>=0)
1081     cluster->SetDispersion(TMath::Sqrt(d)) ;
1082   else
1083     cluster->SetDispersion(0) ;
1084 }
1085
1086 //____________________________________________________________________________
1087 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  AliEMCALGeometry *geom)
1088 {
1089   //This function should be called before the cluster loop
1090   //Before call this function, please recalculate the cluster positions
1091   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1092   //Store matched cluster indexes and residuals
1093
1094   fMatchedTrackIndex  ->Reset();
1095   fMatchedClusterIndex->Reset();
1096   fResidualZ          ->Reset();
1097   fResidualR          ->Reset();
1098
1099   fMatchedTrackIndex  ->Set(500);
1100   fMatchedClusterIndex->Set(500);
1101   fResidualZ          ->Set(500);
1102   fResidualR          ->Set(500);
1103
1104   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1105   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1106
1107   Int_t    matched=0;
1108   Double_t cv[21];
1109   for (Int_t i=0; i<21;i++) cv[i]=0;
1110   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1111   {
1112     AliExternalTrackParam *trackParam=0;
1113
1114     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1115     if(esdevent)
1116       {
1117         AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1118         if(!esdTrack || !IsAccepted(esdTrack)) continue;
1119         const AliESDfriendTrack*  friendTrack = esdTrack->GetFriendTrack();
1120         if(friendTrack && friendTrack->GetTPCOut())
1121           {
1122             //Use TPC Out as starting point if it is available
1123             trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1124           }
1125         else
1126           {
1127             //Otherwise use TPC inner
1128             trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
1129           }
1130       }
1131
1132     //If the input event is AOD, the starting point for extrapolation is at vertex
1133     //AOD tracks are selected according to its bit.
1134     else if(aodevent)
1135       {
1136         AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1137         if(!aodTrack) continue;
1138         if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1139         Double_t pos[3],mom[3];
1140         aodTrack->GetXYZ(pos);
1141         aodTrack->GetPxPyPz(mom);
1142         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()));
1143         trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1144       }
1145
1146     //Return if the input data is not "AOD" or "ESD"
1147     else
1148       {
1149         printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1150         return;
1151       }
1152
1153     if(!trackParam) continue;
1154
1155     Float_t dRMax = fCutR, dZMax=fCutZ;
1156     Int_t index = -1;
1157     if(!clusterArr){// get clusters from event
1158       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1159       {
1160         AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1161         AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1162         if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1163         Float_t tmpR=-1, tmpZ=-1;
1164         if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1165         if(tmpR<dRMax)
1166         {
1167           dRMax=tmpR;
1168           dZMax=tmpZ;
1169           index=icl;
1170         }
1171         delete trkPamTmp;
1172       }//cluster loop
1173     } else { // external cluster array, not from ESD event
1174       for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1175       {
1176         AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1177         AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1178         if(!cluster->IsEMCAL()) continue;
1179         Float_t tmpR=-1, tmpZ=-1;
1180         if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1181         if(tmpR<dRMax)
1182         {
1183           dRMax=tmpR;
1184           dZMax=tmpZ;
1185           index=icl;
1186         }
1187         delete trkPamTmp;
1188       }//cluster loop
1189     }// external list of clusters
1190
1191     if(index>-1)
1192     {
1193       fMatchedTrackIndex  ->AddAt(itr,matched);
1194       fMatchedClusterIndex->AddAt(index,matched);
1195       fResidualZ          ->AddAt(dZMax,matched);
1196       fResidualR          ->AddAt(dRMax,matched);
1197       matched++;
1198     }
1199     delete trackParam;
1200   }//track loop
1201
1202   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1203
1204   fMatchedTrackIndex  ->Set(matched);
1205   fMatchedClusterIndex->Set(matched);
1206   fResidualZ          ->Set(matched);
1207   fResidualR          ->Set(matched);
1208 }
1209
1210 //________________________________________________________________________________
1211 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1212 {
1213   //
1214   // This function returns the index of matched cluster to input track
1215   // Cut on match is dR<10cm by default. Returns -1 if no match is found
1216
1217
1218   Float_t dRMax = fCutR;
1219   Int_t index = -1;
1220
1221   AliExternalTrackParam *trackParam=0;
1222   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
1223   if(friendTrack && friendTrack->GetTPCOut())
1224     trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1225   else
1226     trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1227
1228   if(!trackParam) return index;
1229   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1230     {
1231       AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1232       AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1233       if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1234       Float_t tmpR=-1, tmpZ=-1;
1235       if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1236       if(tmpR>-1 && tmpR<dRMax)
1237         {
1238           dRMax=tmpR;
1239           index=icl;
1240         }
1241         delete trkPamTmp;
1242       }//cluster loop
1243   return index;
1244 }
1245
1246 //________________________________________________________________________________
1247 Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpR, Float_t &tmpZ)
1248 {
1249   //
1250   //Return the residual by extrapolating a track to a cluster
1251   //
1252   if(!trkParam || !cluster) return kFALSE;
1253   Float_t clsPos[3];
1254   Double_t trkPos[3];
1255   cluster->GetPosition(clsPos); //Has been recalculated
1256   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1257   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1258   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1259   trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1260   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE;
1261   trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1262   tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1263   tmpZ = clsPos[2]-trkPos[2];
1264   return kTRUE;
1265 }
1266
1267 //________________________________________________________________________________
1268 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dR, Float_t &dZ)
1269 {
1270   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1271   //Get the residuals dR and dZ for this cluster to the closest track
1272   //Works with ESDs and AODs
1273
1274   if( FindMatchedPosForCluster(clsIndex) >= 999 )
1275   {
1276     AliDebug(2,"No matched tracks found!\n");
1277     dR=999.;
1278     dZ=999.;
1279     return;
1280   }
1281   dR = fResidualR->At(FindMatchedPosForCluster(clsIndex));
1282   dZ = fResidualZ->At(FindMatchedPosForCluster(clsIndex));
1283 }
1284 //________________________________________________________________________________
1285 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dR, Float_t &dZ)
1286 {
1287   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1288   //Get the residuals dR and dZ for this track to the closest cluster
1289   //Works with ESDs and AODs
1290
1291   if( FindMatchedPosForTrack(trkIndex) >= 999 )
1292   {
1293     AliDebug(2,"No matched cluster found!\n");
1294     dR=999.;
1295     dZ=999.;
1296     return;
1297   }
1298   dR = fResidualR->At(FindMatchedPosForTrack(trkIndex));
1299   dZ = fResidualZ->At(FindMatchedPosForTrack(trkIndex));
1300 }
1301
1302 //__________________________________________________________
1303 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1304 {
1305   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1306   //Get the index of matched track to this cluster
1307   //Works with ESDs and AODs
1308
1309   if(IsClusterMatched(clsIndex))
1310     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1311   else
1312     return -1;
1313 }
1314
1315 //__________________________________________________________
1316 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1317 {
1318   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1319   //Get the index of matched cluster to this track
1320   //Works with ESDs and AODs
1321
1322   if(IsTrackMatched(trkIndex))
1323     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1324   else
1325     return -1;
1326 }
1327
1328 //__________________________________________________
1329 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1330 {
1331   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1332   //Returns if the cluster has a match
1333   if(FindMatchedPosForCluster(clsIndex) < 999)
1334     return kTRUE;
1335   else
1336     return kFALSE;
1337 }
1338
1339 //__________________________________________________
1340 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
1341 {
1342   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1343   //Returns if the track has a match
1344   if(FindMatchedPosForTrack(trkIndex) < 999)
1345     return kTRUE;
1346   else
1347     return kFALSE;
1348 }
1349
1350 //__________________________________________________________
1351 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1352 {
1353   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1354   //Returns the position of the match in the fMatchedClusterIndex array
1355   Float_t tmpR = fCutR;
1356   UInt_t pos = 999;
1357
1358   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1359   {
1360     if(fMatchedClusterIndex->At(i)==clsIndex && fResidualR->At(i)<tmpR)
1361     {
1362       pos=i;
1363       tmpR=fResidualR->At(i);
1364       AliDebug(3,Form("Matched cluster index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1365     }
1366   }
1367   return pos;
1368 }
1369
1370 //__________________________________________________________
1371 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1372 {
1373   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1374   //Returns the position of the match in the fMatchedTrackIndex array
1375   Float_t tmpR = fCutR;
1376   UInt_t pos = 999;
1377
1378   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1379   {
1380     if(fMatchedTrackIndex->At(i)==trkIndex && fResidualR->At(i)<tmpR)
1381     {
1382       pos=i;
1383       tmpR=fResidualR->At(i);
1384       AliDebug(3,Form("Matched track index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1385     }
1386   }
1387   return pos;
1388 }
1389
1390 //__________________________________________________________
1391 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1392 {
1393   // check if the cluster survives some quality cut
1394   //
1395   //
1396   Bool_t isGood=kTRUE;
1397   if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1398   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) isGood=kFALSE;
1399   if(!CheckCellFiducialRegion(geom,cluster,cells)) isGood=kFALSE;
1400   if(fRejectExoticCluster && IsExoticCluster(cluster)) isGood=kFALSE;
1401
1402   return isGood;
1403 }
1404
1405 //__________________________________________________________
1406 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1407 {
1408   // Given a esd track, return whether the track survive all the cuts
1409
1410   // The different quality parameter are first
1411   // retrieved from the track. then it is found out what cuts the
1412   // track did not survive and finally the cuts are imposed.
1413
1414   UInt_t status = esdTrack->GetStatus();
1415
1416   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1417   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1418
1419   Float_t chi2PerClusterITS = -1;
1420   Float_t chi2PerClusterTPC = -1;
1421   if (nClustersITS!=0)
1422     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1423   if (nClustersTPC!=0)
1424     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1425
1426
1427   //DCA cuts
1428   Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1429   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1430   SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1431
1432
1433   Float_t b[2];
1434   Float_t bCov[3];
1435   esdTrack->GetImpactParameters(b,bCov);
1436   if (bCov[0]<=0 || bCov[2]<=0) {
1437     AliDebug(1, "Estimated b resolution lower or equal zero!");
1438     bCov[0]=0; bCov[2]=0;
1439   }
1440
1441   Float_t dcaToVertexXY = b[0];
1442   Float_t dcaToVertexZ = b[1];
1443   Float_t dcaToVertex = -1;
1444
1445   if (fCutDCAToVertex2D)
1446     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1447   else
1448     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1449
1450   // cut the track?
1451
1452   Bool_t cuts[kNCuts];
1453   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1454
1455   // track quality cuts
1456   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1457     cuts[0]=kTRUE;
1458   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1459     cuts[1]=kTRUE;
1460   if (nClustersTPC<fCutMinNClusterTPC)
1461     cuts[2]=kTRUE;
1462   if (nClustersITS<fCutMinNClusterITS)
1463     cuts[3]=kTRUE;
1464   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1465     cuts[4]=kTRUE;
1466   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1467     cuts[5]=kTRUE;
1468   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1469     cuts[6]=kTRUE;
1470   if (fCutDCAToVertex2D && dcaToVertex > 1)
1471     cuts[7] = kTRUE;
1472   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1473     cuts[8] = kTRUE;
1474   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1475     cuts[9] = kTRUE;
1476
1477   //Require at least one SPD point + anything else in ITS
1478   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1479     cuts[10] = kTRUE;
1480
1481   Bool_t cut=kFALSE;
1482   for (Int_t i=0; i<kNCuts; i++)
1483     if (cuts[i]) {cut = kTRUE;}
1484
1485     // cut the track
1486   if (cut)
1487     return kFALSE;
1488   else
1489     return kTRUE;
1490 }
1491 //__________________________________________________
1492 void AliEMCALRecoUtils::InitTrackCuts()
1493 {
1494   //Intilize the track cut criteria
1495   //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1496   //Also you can customize the cuts using the setters
1497
1498   //TPC
1499   SetMinNClustersTPC(70);
1500   SetMaxChi2PerClusterTPC(4);
1501   SetAcceptKinkDaughters(kFALSE);
1502   SetRequireTPCRefit(kTRUE);
1503
1504   //ITS
1505   SetRequireITSRefit(kTRUE);
1506   SetMaxDCAToVertexZ(2);
1507   SetDCAToVertex2D(kFALSE);
1508   SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1509   SetMinNClustersITS();
1510 }
1511
1512 //___________________________________________________
1513 void AliEMCALRecoUtils::Print(const Option_t *) const
1514 {
1515   // Print Parameters
1516
1517   printf("AliEMCALRecoUtils Settings: \n");
1518   printf("Misalignment shifts\n");
1519   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,
1520                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1521                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1522   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1523   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1524
1525   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1526
1527   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1528   printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1529
1530   printf("Track cuts: \n");
1531   printf("AOD track selection mask: %d\n",fAODFilterMask);
1532   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1533   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1534   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1535   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1536   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1537
1538 }
1539
1540 //_____________________________________________________________________
1541 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1542   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1543   //Do it only once and only if it is requested
1544
1545   if(!fUseTimeCorrectionFactors) return;
1546   if(fTimeCorrectionFactorsSet)  return;
1547
1548   printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1549
1550   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
1551   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1552
1553   SwitchOnRecalibration();
1554   for(Int_t ism = 0; ism < 4; ism++){
1555     for(Int_t icol = 0; icol < 48; icol++){
1556       for(Int_t irow = 0; irow < 24; irow++){
1557         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1558         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1559         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1560         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1561         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1562         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1563       }
1564     }
1565   }
1566    fTimeCorrectionFactorsSet = kTRUE;
1567 }
1568