EMCAL/AliEMCALRecoUtils.cxx

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