EMCAL/AliEMCALRecoUtils.cxx

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