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     }
1205     else { // external cluster array, not from ESD event
1206       for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1207       {
1208         AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1209         if(!cluster){
1210           AliInfo("Cluster not found!!!");
1211           continue;
1212         }
1213         if(!cluster->IsEMCAL()) continue;
1214         AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1215         Float_t tmpEta=-999, tmpPhi=-999;
1216         if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1217         if(fCutEtaPhiSum)
1218         {
1219           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1220           if(tmpR<dRMax)
1221           {
1222             dRMax=tmpR;
1223             dEtaMax=tmpEta;
1224             dPhiMax=tmpPhi;
1225             index=icl;
1226           }
1227         }
1228         else if(fCutEtaPhiSeparate)
1229         {
1230           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1231           {
1232             dEtaMax = tmpEta;
1233             dPhiMax = tmpPhi;
1234             index=icl;
1235           }
1236         }
1237         else
1238         {
1239           printf("Error: please specify your cut criteria\n");
1240           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1241           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1242           if(aodevent && trackParam) delete trackParam;
1243           return;
1244         }
1245       }//cluster loop
1246     }// external list of clusters
1247
1248     if(index>-1)
1249     {
1250       fMatchedTrackIndex  ->AddAt(itr,matched);
1251       fMatchedClusterIndex->AddAt(index,matched);
1252       fResidualEta          ->AddAt(dEtaMax,matched);
1253       fResidualPhi          ->AddAt(dPhiMax,matched);
1254       matched++;
1255     }
1256     if(aodevent && trackParam) delete trackParam;
1257   }//track loop
1258
1259   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1260
1261   fMatchedTrackIndex  ->Set(matched);
1262   fMatchedClusterIndex->Set(matched);
1263   fResidualPhi          ->Set(matched);
1264   fResidualEta          ->Set(matched);
1265 }
1266
1267 //________________________________________________________________________________
1268 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1269 {
1270   //
1271   // This function returns the index of matched cluster to input track
1272   // Returns -1 if no match is found
1273
1274
1275   Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
1276   Int_t index = -1;
1277
1278   AliExternalTrackParam *trackParam=0;
1279   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
1280   if(friendTrack && friendTrack->GetTPCOut())
1281     trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1282   else
1283     trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1284
1285   if(!trackParam) return index;
1286   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1287   {
1288     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1289     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1290     AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1291     Float_t tmpEta=-999, tmpPhi=-999;
1292     if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1293     if(fCutEtaPhiSum)
1294     {
1295       Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1296       if(tmpR<dRMax)
1297             {
1298               dRMax=tmpR;
1299               dEtaMax=tmpEta;
1300               dPhiMax=tmpPhi;
1301               index=icl;
1302             }
1303     }
1304     else if(fCutEtaPhiSeparate)
1305     {
1306       if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1307             {
1308               dEtaMax = tmpEta;
1309               dPhiMax = tmpPhi;
1310               index=icl;
1311             }
1312     }
1313     else
1314     {
1315       printf("Error: please specify your cut criteria\n");
1316       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1317       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1318       return -1;
1319     }
1320   }//cluster loop
1321   return index;
1322 }
1323
1324 //________________________________________________________________________________
1325 Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1326 {
1327   //
1328   //Return the residual by extrapolating a track to a cluster
1329   //
1330   if(!trkParam || !cluster) return kFALSE;
1331   Float_t clsPos[3];
1332   Double_t trkPos[3];
1333   cluster->GetPosition(clsPos); //Has been recalculated
1334   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1335   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1336   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1337   trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1338   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return kFALSE;
1339   trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1340
1341   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1342   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1343
1344   Float_t clsPhi = (Float_t)clsPosVec.Phi();
1345   if(clsPhi<0) clsPhi+=2*TMath::Pi();
1346   Float_t trkPhi = (Float_t)trkPosVec.Phi();
1347   if(trkPhi<0) trkPhi+=2*TMath::Pi();
1348   tmpPhi = clsPhi-trkPhi;  // track cluster matching
1349   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
1350
1351   return kTRUE;
1352 }
1353
1354 //________________________________________________________________________________
1355 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1356 {
1357   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1358   //Get the residuals dEta and dPhi for this cluster to the closest track
1359   //Works with ESDs and AODs
1360
1361   if( FindMatchedPosForCluster(clsIndex) >= 999 )
1362   {
1363     AliDebug(2,"No matched tracks found!\n");
1364     dEta=999.;
1365     dPhi=999.;
1366     return;
1367   }
1368   dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1369   dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1370 }
1371 //________________________________________________________________________________
1372 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1373 {
1374   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1375   //Get the residuals dEta and dPhi for this track to the closest cluster
1376   //Works with ESDs and AODs
1377
1378   if( FindMatchedPosForTrack(trkIndex) >= 999 )
1379   {
1380     AliDebug(2,"No matched cluster found!\n");
1381     dEta=999.;
1382     dPhi=999.;
1383     return;
1384   }
1385   dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1386   dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1387 }
1388
1389 //__________________________________________________________
1390 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1391 {
1392   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1393   //Get the index of matched track to this cluster
1394   //Works with ESDs and AODs
1395
1396   if(IsClusterMatched(clsIndex))
1397     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1398   else
1399     return -1;
1400 }
1401
1402 //__________________________________________________________
1403 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1404 {
1405   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1406   //Get the index of matched cluster to this track
1407   //Works with ESDs and AODs
1408
1409   if(IsTrackMatched(trkIndex))
1410     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1411   else
1412     return -1;
1413 }
1414
1415 //__________________________________________________
1416 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
1417 {
1418   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1419   //Returns if the cluster has a match
1420   if(FindMatchedPosForCluster(clsIndex) < 999)
1421     return kTRUE;
1422   else
1423     return kFALSE;
1424 }
1425
1426 //__________________________________________________
1427 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
1428 {
1429   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1430   //Returns if the track has a match
1431   if(FindMatchedPosForTrack(trkIndex) < 999)
1432     return kTRUE;
1433   else
1434     return kFALSE;
1435 }
1436
1437 //__________________________________________________________
1438 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1439 {
1440   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1441   //Returns the position of the match in the fMatchedClusterIndex array
1442   Float_t tmpR = fCutR;
1443   UInt_t pos = 999;
1444
1445   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1446   {
1447     if(fMatchedClusterIndex->At(i)==clsIndex)
1448       {
1449         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1450         if(r<tmpR)
1451           {
1452             pos=i;
1453             tmpR=r;
1454             AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1455           }
1456       }
1457   }
1458   return pos;
1459 }
1460
1461 //__________________________________________________________
1462 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1463 {
1464   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1465   //Returns the position of the match in the fMatchedTrackIndex array
1466   Float_t tmpR = fCutR;
1467   UInt_t pos = 999;
1468
1469   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1470   {
1471     if(fMatchedTrackIndex->At(i)==trkIndex)
1472       {
1473         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1474         if(r<tmpR)
1475           {
1476             pos=i;
1477             tmpR=r;
1478             AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1479           }
1480       }
1481   }
1482   return pos;
1483 }
1484
1485 //__________________________________________________________
1486 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1487 {
1488   // check if the cluster survives some quality cut
1489   //
1490   //
1491   Bool_t isGood=kTRUE;
1492   if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1493   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1494   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1495   if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
1496
1497   return isGood;
1498 }
1499
1500 //__________________________________________________________
1501 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1502 {
1503   // Given a esd track, return whether the track survive all the cuts
1504
1505   // The different quality parameter are first
1506   // retrieved from the track. then it is found out what cuts the
1507   // track did not survive and finally the cuts are imposed.
1508
1509   UInt_t status = esdTrack->GetStatus();
1510
1511   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1512   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1513
1514   Float_t chi2PerClusterITS = -1;
1515   Float_t chi2PerClusterTPC = -1;
1516   if (nClustersITS!=0)
1517     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1518   if (nClustersTPC!=0)
1519     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1520
1521
1522   //DCA cuts
1523   Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1524   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1525   SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1526
1527
1528   Float_t b[2];
1529   Float_t bCov[3];
1530   esdTrack->GetImpactParameters(b,bCov);
1531   if (bCov[0]<=0 || bCov[2]<=0) {
1532     AliDebug(1, "Estimated b resolution lower or equal zero!");
1533     bCov[0]=0; bCov[2]=0;
1534   }
1535
1536   Float_t dcaToVertexXY = b[0];
1537   Float_t dcaToVertexZ = b[1];
1538   Float_t dcaToVertex = -1;
1539
1540   if (fCutDCAToVertex2D)
1541     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1542   else
1543     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1544
1545   // cut the track?
1546
1547   Bool_t cuts[kNCuts];
1548   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1549
1550   // track quality cuts
1551   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1552     cuts[0]=kTRUE;
1553   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1554     cuts[1]=kTRUE;
1555   if (nClustersTPC<fCutMinNClusterTPC)
1556     cuts[2]=kTRUE;
1557   if (nClustersITS<fCutMinNClusterITS)
1558     cuts[3]=kTRUE;
1559   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1560     cuts[4]=kTRUE;
1561   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1562     cuts[5]=kTRUE;
1563   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1564     cuts[6]=kTRUE;
1565   if (fCutDCAToVertex2D && dcaToVertex > 1)
1566     cuts[7] = kTRUE;
1567   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1568     cuts[8] = kTRUE;
1569   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1570     cuts[9] = kTRUE;
1571
1572   //Require at least one SPD point + anything else in ITS
1573   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1574     cuts[10] = kTRUE;
1575
1576   Bool_t cut=kFALSE;
1577   for (Int_t i=0; i<kNCuts; i++)
1578     if (cuts[i]) {cut = kTRUE;}
1579
1580     // cut the track
1581   if (cut)
1582     return kFALSE;
1583   else
1584     return kTRUE;
1585 }
1586 //__________________________________________________
1587 void AliEMCALRecoUtils::InitTrackCuts()
1588 {
1589   //Intilize the track cut criteria
1590   //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1591   //Also you can customize the cuts using the setters
1592
1593   //TPC
1594   SetMinNClustersTPC(70);
1595   SetMaxChi2PerClusterTPC(4);
1596   SetAcceptKinkDaughters(kFALSE);
1597   SetRequireTPCRefit(kTRUE);
1598
1599   //ITS
1600   SetRequireITSRefit(kTRUE);
1601   SetMaxDCAToVertexZ(2);
1602   SetDCAToVertex2D(kFALSE);
1603   SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1604   SetMinNClustersITS();
1605 }
1606
1607 //___________________________________________________
1608 void AliEMCALRecoUtils::Print(const Option_t *) const
1609 {
1610   // Print Parameters
1611
1612   printf("AliEMCALRecoUtils Settings: \n");
1613   printf("Misalignment shifts\n");
1614   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,
1615                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1616                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1617   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1618   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1619
1620   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1621
1622   printf("Matching criteria: ");
1623   if(fCutEtaPhiSum)
1624     {
1625       printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1626     }
1627   else if(fCutEtaPhiSeparate)
1628     {
1629       printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1630     }
1631   else
1632     {
1633       printf("Error\n");
1634       printf("please specify your cut criteria\n");
1635       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1636       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1637     }
1638
1639   printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1640
1641   printf("Track cuts: \n");
1642   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
1643   printf("AOD track selection mask: %d\n",fAODFilterMask);
1644   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1645   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1646   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1647   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1648   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1649
1650 }
1651
1652 //_____________________________________________________________________
1653 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1654   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1655   //Do it only once and only if it is requested
1656
1657   if(!fUseTimeCorrectionFactors) return;
1658   if(fTimeCorrectionFactorsSet)  return;
1659
1660   printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1661
1662   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
1663   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1664
1665   SwitchOnRecalibration();
1666   for(Int_t ism = 0; ism < 4; ism++){
1667     for(Int_t icol = 0; icol < 48; icol++){
1668       for(Int_t irow = 0; irow < 24; irow++){
1669         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1670         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1671         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1672         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1673         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1674         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1675       }
1676     }
1677   }
1678    fTimeCorrectionFactorsSet = kTRUE;
1679 }
1680