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