EMCAL/AliEMCALRecoUtils.cxx

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