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