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