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