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