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 ///////////////////////////////////////////////////////////////////////////////
  26
  27 // --- standard c ---
  28
  29 // standard C++ includes
  30 //#include <Riostream.h>
  31
  32 // ROOT includes
  33 #include <TGeoManager.h>
  34 #include <TGeoMatrix.h>
  35 #include <TGeoBBox.h>
  36
  37 // STEER includes
  38 #include "AliEMCALRecoUtils.h"
  39 #include "AliEMCALGeometry.h"
  40 #include "AliVCluster.h"
  41 #include "AliVCaloCells.h"
  42 #include "AliVEvent.h"
  43 #include "AliLog.h"
  44 #include "AliEMCALPIDUtils.h"
  45 #include "AliPID.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliEMCALTrack.h"
  49 #include "AliEMCALCalibTimeDepCorrection.h"
  50
  51 ClassImp(AliEMCALRecoUtils)
  52
  53 //______________________________________________
  54 AliEMCALRecoUtils::AliEMCALRecoUtils():
  55   fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
  56   fPosAlgo(kUnchanged),fW0(4.),
  57   fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
  58   fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
  59   fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
  60   fMatchedClusterIndex(0x0), fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
  61   fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
  62   fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
  63   fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
  64   fUseTimeCorrectionFactors(kFALSE),  fTimeCorrectionFactorsSet(kFALSE)
  65 {
  66 //
  67   // Constructor.
  68   // Initialize all constant values which have to be used
  69   // during Reco algorithm execution
  70   //
  71
  72   for(Int_t i = 0; i < 15 ; i++) {
  73       fMisalTransShift[i] = 0.;
  74       fMisalRotShift[i] = 0.;
  75   }
  76   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = 0.;
  77   //For kPi0GammaGamma case, but default is no correction
  78   fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
  79   fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
  80   fNonLinearityParams[2] = 1.046;
  81
  82   fMatchedClusterIndex = new TArrayI();
  83   fResidualZ = new TArrayF();
  84   fResidualR = new TArrayF();
  85
  86   fPIDUtils = new AliEMCALPIDUtils();
  87
  88   InitTrackCuts();
  89
  90 }
  91
  92 //______________________________________________________________________
  93 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
  94 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
  95   fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
  96   fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
  97   fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
  98   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
  99   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
 100   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
 101   fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
 102   fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
 103   fCutR(reco.fCutR),fCutZ(reco.fCutZ),
 104   fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
 105   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
 106   fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
 107   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
 108   fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
 109   fPIDUtils(reco.fPIDUtils),
 110   fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors),  fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
 111 {
 112   //Copy ctor
 113
 114   for(Int_t i = 0; i < 15 ; i++) {
 115       fMisalRotShift[i] = reco.fMisalRotShift[i];
 116       fMisalTransShift[i] = reco.fMisalTransShift[i];
 117   }
 118   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
 119
 120 }
 121
 122
 123 //______________________________________________________________________
 124 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
 125 {
 126   //Assignment operator
 127
 128   if(this == &reco)return *this;
 129   ((TNamed *)this)->operator=(reco);
 130
 131   fNonLinearityFunction      = reco.fNonLinearityFunction;
 132   fParticleType              = reco.fParticleType;
 133   fPosAlgo                   = reco.fPosAlgo;
 134   fW0                        = reco.fW0;
 135   fRecalibration             = reco.fRecalibration;
 136   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
 137   fRemoveBadChannels         = reco.fRemoveBadChannels;
 138   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
 139   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
 140   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
 141   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
 142
 143
 144   for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
 145   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
 146
 147   fCutR                      = reco.fCutR;
 148   fCutZ                      = reco.fCutZ;
 149
 150   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
 151   fCutMinNClusterITS         = reco.fCutMinNClusterITS;
 152   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
 153   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
 154   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
 155   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
 156   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
 157   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
 158   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
 159   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
 160
 161   fPIDUtils                  = reco.fPIDUtils;
 162
 163   fUseTimeCorrectionFactors  = reco.fUseTimeCorrectionFactors;
 164   fTimeCorrectionFactorsSet  = reco.fTimeCorrectionFactorsSet;
 165
 166
 167   if(reco.fResidualR){
 168     // assign or copy construct
 169     if(fResidualR){
 170       *fResidualR = *reco.fResidualR;
 171     }
 172     else fResidualR = new TArrayF(*reco.fResidualR);
 173   }
 174   else{
 175     if(fResidualR)delete fResidualR;
 176     fResidualR = 0;
 177   }
 178
 179   if(reco.fResidualZ){
 180     // assign or copy construct
 181     if(fResidualZ){
 182       *fResidualZ = *reco.fResidualZ;
 183     }
 184     else fResidualZ = new TArrayF(*reco.fResidualZ);
 185   }
 186   else{
 187     if(fResidualZ)delete fResidualZ;
 188     fResidualZ = 0;
 189   }
 190
 191
 192   if(reco.fMatchedClusterIndex){
 193     // assign or copy construct
 194     if(fMatchedClusterIndex){
 195       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
 196     }
 197     else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
 198   }
 199   else{
 200     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
 201     fMatchedClusterIndex = 0;
 202   }
 203
 204
 205   return *this;
 206 }
 207
 208
 209 //__________________________________________________
 210 AliEMCALRecoUtils::~AliEMCALRecoUtils()
 211 {
 212   //Destructor.
 213
 214         if(fEMCALRecalibrationFactors) {
 215                 fEMCALRecalibrationFactors->Clear();
 216                 delete  fEMCALRecalibrationFactors;
 217         }
 218
 219   if(fEMCALBadChannelMap) {
 220                 fEMCALBadChannelMap->Clear();
 221                 delete  fEMCALBadChannelMap;
 222         }
 223
 224   if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
 225   if(fResidualR)           {delete fResidualR;           fResidualR=0;}
 226   if(fResidualZ)           {delete fResidualZ;           fResidualZ=0;}
 227
 228 }
 229
 230 //_______________________________________________________________
 231 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
 232 {
 233         // Given the list of AbsId of the cluster, get the maximum cell and
 234         // check if there are fNCellsFromBorder from the calorimeter border
 235
 236   //If the distance to the border is 0 or negative just exit accept all clusters
 237         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
 238
 239   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;
 240   Bool_t shared = kFALSE;
 241   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
 242
 243   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
 244            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
 245
 246         if(absIdMax==-1) return kFALSE;
 247
 248         //Check if the cell is close to the borders:
 249         Bool_t okrow = kFALSE;
 250         Bool_t okcol = kFALSE;
 251
 252   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 253     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
 254                   iSM,ieta,iphi));
 255   }
 256
 257   //Check rows/phi
 258   if(iSM < 10){
 259     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
 260   }
 261   else{
 262     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
 263   }
 264
 265   //Check columns/eta
 266   if(!fNoEMCALBorderAtEta0){
 267     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
 268   }
 269   else{
 270     if(iSM%2==0){
 271       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;
 272     }
 273     else {
 274       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;
 275     }
 276   }//eta 0 not checked
 277
 278   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
 279            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
 280
 281         if (okcol && okrow) {
 282     //printf("Accept\n");
 283     return kTRUE;
 284   }
 285         else  {
 286     //printf("Reject\n");
 287     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
 288     return kFALSE;
 289   }
 290
 291 }
 292
 293
 294 //_________________________________________________________________________________________________________
 295 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
 296         // Check that in the cluster cells, there is no bad channel of those stored
 297         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 298
 299         if(!fRemoveBadChannels)  return kFALSE;
 300         if(!fEMCALBadChannelMap) return kFALSE;
 301
 302         Int_t icol = -1;
 303         Int_t irow = -1;
 304         Int_t imod = -1;
 305         for(Int_t iCell = 0; iCell<nCells; iCell++){
 306
 307                 //Get the column and row
 308     Int_t iTower = -1, iIphi = -1, iIeta = -1;
 309     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
 310     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
 311     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 312     if(GetEMCALChannelStatus(imod, icol, irow)){
 313       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
 314       return kTRUE;
 315     }
 316
 317         }// cell cluster loop
 318
 319         return kFALSE;
 320
 321 }
 322
 323 //__________________________________________________
 324 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
 325 // Correct cluster energy from non linearity functions
 326   Float_t energy = cluster->E();
 327
 328   switch (fNonLinearityFunction) {
 329
 330     case kPi0MC:
 331     {
 332       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
 333       //Double_t fNonLinearityParams[0] = 1.001;
 334       //Double_t fNonLinearityParams[1] = -0.01264;
 335       //Double_t fNonLinearityParams[2] = -0.03632;
 336       //Double_t fNonLinearityParams[3] = 0.1798;
 337       //Double_t fNonLinearityParams[4] = -0.522;
 338        energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 339                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 340                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 341       break;
 342     }
 343
 344     case kPi0GammaGamma:
 345     {
 346       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 347       //Double_t fNonLinearityParams[0] = 1.04;
 348       //Double_t fNonLinearityParams[1] = -0.1445;
 349       //Double_t fNonLinearityParams[2] = 1.046;
 350       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 351       break;
 352     }
 353
 354     case kPi0GammaConversion:
 355     {
 356       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 357       //fNonLinearityParams[0] = 0.139393/0.1349766;
 358       //fNonLinearityParams[1] = 0.0566186;
 359       //fNonLinearityParams[2] = 0.982133;
 360       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 361
 362       break;
 363     }
 364
 365     case kBeamTest:
 366     {
 367       //From beam test, Alexei's results, for different ZS thresholds
 368       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
 369       //fNonLinearityParams[0] = 1.007;      1.003;      1.002
 370       //fNonLinearityParams[1] = 0.894;      0.719;      0.797
 371       //fNonLinearityParams[2] = 0.246;      0.334;      0.358
 372       //Rescale the param[0] with 1.03
 373       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
 374
 375       break;
 376     }
 377
 378     case kNoCorrection:
 379       AliDebug(2,"No correction on the energy\n");
 380       break;
 381
 382   }
 383
 384   return energy;
 385
 386 }
 387
 388 //__________________________________________________
 389 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
 390 {
 391   //Calculate shower depth for a given cluster energy and particle type
 392
 393   // parameters
 394   Float_t x0    = 1.31;
 395   Float_t ecr   = 8;
 396   Float_t depth = 0;
 397
 398   switch ( iParticle )
 399   {
 400     case kPhoton:
 401       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 402       break;
 403
 404     case kElectron:
 405       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
 406       break;
 407
 408     case kHadron:
 409       // hadron
 410       // boxes anc. here
 411       if(gGeoManager){
 412         gGeoManager->cd("ALIC_1/XEN1_1");
 413         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
 414         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
 415         if(geoSM){
 416           TGeoVolume      *geoSMVol   = geoSM->GetVolume();
 417           TGeoShape       *geoSMShape = geoSMVol->GetShape();
 418           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
 419           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
 420           else AliFatal("Null GEANT box");
 421         }else AliFatal("NULL  GEANT node matrix");
 422       }
 423       else{//electron
 424         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
 425       }
 426
 427       break;
 428
 429     default://photon
 430       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 431   }
 432
 433   return depth;
 434
 435 }
 436
 437 //__________________________________________________
 438 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
 439                                          Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
 440 {
 441   //For a given CaloCluster gets the absId of the cell
 442   //with maximum energy deposit.
 443
 444   Double_t eMax        = -1.;
 445   Double_t eCell       = -1.;
 446   Float_t  fraction    = 1.;
 447   Float_t  recalFactor = 1.;
 448   Int_t    cellAbsId   = -1 ;
 449
 450   Int_t iTower  = -1;
 451   Int_t iIphi   = -1;
 452   Int_t iIeta   = -1;
 453   Int_t iSupMod0= -1;
 454         //printf("---Max?\n");
 455   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 456     cellAbsId = clu->GetCellAbsId(iDig);
 457     fraction  = clu->GetCellAmplitudeFraction(iDig);
 458     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
 459     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 460     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 461     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 462     if(iDig==0) iSupMod0=iSupMod;
 463     else if(iSupMod0!=iSupMod) {
 464       shared = kTRUE;
 465       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
 466     }
 467     if(IsRecalibrationOn()) {
 468       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 469     }
 470     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 471     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
 472     if(eCell > eMax)  {
 473       eMax  = eCell;
 474       absId = cellAbsId;
 475       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 476     }
 477   }// cell loop
 478
 479   //Get from the absid the supermodule, tower and eta/phi numbers
 480   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 481   //Gives SuperModule and Tower numbers
 482   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 483                                          iIphi, iIeta,iphi,ieta);
 484   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
 485   //printf("Max end---\n");
 486
 487 }
 488
 489 //________________________________________________________________
 490 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
 491         //Init EMCAL recalibration factors
 492         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
 493         //In order to avoid rewriting the same histograms
 494         Bool_t oldStatus = TH1::AddDirectoryStatus();
 495         TH1::AddDirectory(kFALSE);
 496
 497         fEMCALRecalibrationFactors = new TObjArray(10);
 498         for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
 499         //Init the histograms with 1
 500         for (Int_t sm = 0; sm < 12; sm++) {
 501                 for (Int_t i = 0; i < 48; i++) {
 502                         for (Int_t j = 0; j < 24; j++) {
 503                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
 504                         }
 505                 }
 506         }
 507         fEMCALRecalibrationFactors->SetOwner(kTRUE);
 508         fEMCALRecalibrationFactors->Compress();
 509
 510         //In order to avoid rewriting the same histograms
 511         TH1::AddDirectory(oldStatus);
 512 }
 513
 514
 515 //________________________________________________________________
 516 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 517         //Init EMCAL bad channels map
 518         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
 519         //In order to avoid rewriting the same histograms
 520         Bool_t oldStatus = TH1::AddDirectoryStatus();
 521         TH1::AddDirectory(kFALSE);
 522
 523         fEMCALBadChannelMap = new TObjArray(10);
 524         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
 525         for (int i = 0; i < 10; i++) {
 526                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
 527         }
 528
 529         //delete hTemp;
 530
 531         fEMCALBadChannelMap->SetOwner(kTRUE);
 532         fEMCALBadChannelMap->Compress();
 533
 534         //In order to avoid rewriting the same histograms
 535         TH1::AddDirectory(oldStatus);
 536 }
 537
 538 //________________________________________________________________
 539 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
 540         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 541
 542         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 543         UShort_t * index    = cluster->GetCellsAbsId() ;
 544         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 545         Int_t ncells = cluster->GetNCells();
 546
 547         //Initialize some used variables
 548         Float_t energy = 0;
 549         Int_t absId    = -1;
 550   Int_t icol = -1, irow = -1, imod=1;
 551         Float_t factor = 1, frac = 0;
 552
 553         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 554         for(Int_t icell = 0; icell < ncells; icell++){
 555                 absId = index[icell];
 556                 frac =  fraction[icell];
 557                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
 558                 Int_t iTower = -1, iIphi = -1, iIeta = -1;
 559                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 560                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
 561                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 562                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
 563     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
 564              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
 565
 566                 energy += cells->GetCellAmplitude(absId)*factor*frac;
 567         }
 568
 569
 570                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
 571
 572         cluster->SetE(energy);
 573
 574 }
 575
 576
 577 //__________________________________________________
 578 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 579 {
 580   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 581
 582   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
 583   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
 584   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
 585
 586 }
 587
 588 //__________________________________________________
 589 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 590 {
 591   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 592   // The algorithm is a copy of what is done in AliEMCALRecPoint
 593
 594   Double_t eCell       = 0.;
 595   Float_t  fraction    = 1.;
 596   Float_t  recalFactor = 1.;
 597
 598   Int_t    absId   = -1;
 599   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
 600   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
 601   Float_t  weight = 0.,  totalWeight=0.;
 602   Float_t  newPos[3] = {0,0,0};
 603   Double_t pLocal[3], pGlobal[3];
 604   Bool_t shared = kFALSE;
 605
 606   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
 607   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
 608   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
 609
 610   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
 611
 612   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 613     absId = clu->GetCellAbsId(iDig);
 614     fraction  = clu->GetCellAmplitudeFraction(iDig);
 615     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 616     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
 617     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 618
 619     if(IsRecalibrationOn()) {
 620       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
 621     }
 622     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 623
 624     weight = GetCellWeight(eCell,clEnergy);
 625     //printf("cell energy %f, weight %f\n",eCell,weight);
 626     totalWeight += weight;
 627     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
 628     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
 629     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
 630     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
 631
 632     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
 633
 634   }// cell loop
 635
 636   if(totalWeight>0){
 637     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
 638   }
 639
 640   //Float_t pos[]={0,0,0};
 641   //clu->GetPosition(pos);
 642   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
 643   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 644
 645         if(iSupModMax > 1) {//sector 1
 646           newPos[0] +=fMisalTransShift[3];//-=3.093;
 647           newPos[1] +=fMisalTransShift[4];//+=6.82;
 648           newPos[2] +=fMisalTransShift[5];//+=1.635;
 649     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
 650
 651         }
 652         else {//sector 0
 653           newPos[0] +=fMisalTransShift[0];//+=1.134;
 654           newPos[1] +=fMisalTransShift[1];//+=8.2;
 655           newPos[2] +=fMisalTransShift[2];//+=1.197;
 656     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
 657
 658         }
 659   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 660
 661   clu->SetPosition(newPos);
 662
 663 }
 664
 665 //__________________________________________________
 666 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 667 {
 668   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 669   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
 670
 671   Double_t eCell       = 1.;
 672   Float_t  fraction    = 1.;
 673   Float_t  recalFactor = 1.;
 674
 675   Int_t absId   = -1;
 676   Int_t iTower  = -1;
 677   Int_t iIphi   = -1, iIeta   = -1;
 678         Int_t iSupMod = -1, iSupModMax = -1;
 679   Int_t iphi = -1, ieta =-1;
 680   Bool_t shared = kFALSE;
 681
 682   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
 683   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
 684   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
 685
 686   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
 687   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
 688   Int_t startingSM = -1;
 689
 690   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 691     absId = clu->GetCellAbsId(iDig);
 692     fraction  = clu->GetCellAmplitudeFraction(iDig);
 693     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 694     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 695     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 696
 697     if     (iDig==0)  startingSM = iSupMod;
 698     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 699
 700     eCell  = cells->GetCellAmplitude(absId);
 701
 702     if(IsRecalibrationOn()) {
 703       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 704     }
 705     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 706
 707     weight = GetCellWeight(eCell,clEnergy);
 708     if(weight < 0) weight = 0;
 709     totalWeight += weight;
 710     weightedCol += ieta*weight;
 711     weightedRow += iphi*weight;
 712
 713     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
 714
 715     }// cell loop
 716
 717   Float_t xyzNew[]={0.,0.,0.};
 718   if(areInSameSM == kTRUE) {
 719     //printf("In Same SM\n");
 720     weightedCol = weightedCol/totalWeight;
 721     weightedRow = weightedRow/totalWeight;
 722     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 723   }
 724   else {
 725     //printf("In Different SM\n");
 726     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 727   }
 728
 729   clu->SetPosition(xyzNew);
 730
 731 }
 732
 733 //____________________________________________________________________________
 734 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
 735
 736   //re-evaluate distance to bad channel with updated bad map
 737
 738   if(!fRecalDistToBadChannels) return;
 739
 740         //Get channels map of the supermodule where the cluster is.
 741   Int_t absIdMax        = -1, iSupMod =-1, icolM = -1, irowM = -1;
 742   Bool_t shared = kFALSE;
 743   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
 744   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
 745
 746   Int_t dRrow, dRcol;
 747         Float_t  minDist = 10000.;
 748         Float_t  dist    = 0.;
 749
 750   //Loop on tower status map
 751         for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
 752                 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
 753                         //Check if tower is bad.
 754                         if(hMap->GetBinContent(icol,irow)==0) continue;
 755       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
 756       //       iSupMod,icol, irow, icolM,irowM);
 757
 758       dRrow=TMath::Abs(irowM-irow);
 759       dRcol=TMath::Abs(icolM-icol);
 760       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
 761                         if(dist < minDist){
 762         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
 763         minDist = dist;
 764       }
 765
 766                 }
 767         }
 768
 769         //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
 770         if (shared) {
 771                 TH2D* hMap2 = 0;
 772                 Int_t iSupMod2 = -1;
 773
 774                 //The only possible combinations are (0,1), (2,3) ... (8,9)
 775                 if(iSupMod%2) iSupMod2 = iSupMod-1;
 776                 else          iSupMod2 = iSupMod+1;
 777                 hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
 778
 779                 //Loop on tower status map of second super module
 780                 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
 781                         for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
 782                                 //Check if tower is bad.
 783                                 if(hMap2->GetBinContent(icol,irow)==0) continue;
 784                                 //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",
 785           //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
 786
 787         dRrow=TMath::Abs(irow-irowM);
 788
 789         if(iSupMod%2) {
 790                                   dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
 791                                 }
 792         else {
 793           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
 794                                 }
 795
 796                                 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
 797         if(dist < minDist) minDist = dist;
 798
 799                         }
 800                 }
 801
 802         }// shared cluster in 2 SuperModules
 803
 804   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
 805   cluster->SetDistanceToBadChannel(minDist);
 806
 807 }
 808
 809 //____________________________________________________________________________
 810 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
 811
 812   //re-evaluate identification parameters with bayesian
 813
 814         if ( cluster->GetM02() != 0)
 815     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
 816
 817   Float_t pidlist[AliPID::kSPECIESN+1];
 818         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
 819
 820   cluster->SetPID(pidlist);
 821
 822 }
 823
 824 //____________________________________________________________________________
 825 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
 826 {
 827   // Calculates new center of gravity in the local EMCAL-module coordinates
 828   // and tranfers into global ALICE coordinates
 829   // Calculates Dispersion and main axis
 830
 831   Int_t nstat  = 0;
 832   Float_t wtot = 0. ;
 833   Double_t eCell       = 0.;
 834   Float_t  fraction    = 1.;
 835   Float_t  recalFactor = 1.;
 836
 837   Int_t iSupMod = -1;
 838   Int_t iTower  = -1;
 839   Int_t iIphi   = -1;
 840   Int_t iIeta   = -1;
 841   Int_t iphi    = -1;
 842   Int_t ieta    = -1;
 843   Double_t etai = -1.;
 844   Double_t phii = -1.;
 845
 846   Double_t w     = 0.;
 847   Double_t d     = 0.;
 848   Double_t dxx   = 0.;
 849   Double_t dzz   = 0.;
 850   Double_t dxz   = 0.;
 851   Double_t xmean = 0.;
 852   Double_t zmean = 0.;
 853
 854   //Loop on cells
 855   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
 856
 857     //Get from the absid the supermodule, tower and eta/phi numbers
 858     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 859     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 860
 861     //Get the cell energy, if recalibration is on, apply factors
 862     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
 863     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 864     if(IsRecalibrationOn()) {
 865       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 866     }
 867     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
 868
 869     if(cluster->E() > 0 && eCell > 0){
 870
 871       w  = GetCellWeight(eCell,cluster->E());
 872
 873       etai=(Double_t)ieta;
 874       phii=(Double_t)iphi;
 875       if(w > 0.0) {
 876         wtot += w ;
 877         nstat++;
 878         //Shower shape
 879         dxx  += w * etai * etai ;
 880         xmean+= w * etai ;
 881         dzz  += w * phii * phii ;
 882         zmean+= w * phii ;
 883         dxz  += w * etai * phii ;
 884       }
 885     }
 886     else
 887       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
 888   }//cell loop
 889
 890   //Normalize to the weight
 891   if (wtot > 0) {
 892     xmean /= wtot ;
 893     zmean /= wtot ;
 894   }
 895   else
 896     AliError(Form("Wrong weight %f\n", wtot));
 897
 898   //Calculate dispersion
 899   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
 900
 901     //Get from the absid the supermodule, tower and eta/phi numbers
 902     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 903     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 904
 905     //Get the cell energy, if recalibration is on, apply factors
 906     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
 907     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 908     if(IsRecalibrationOn()) {
 909       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 910     }
 911     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
 912
 913     if(cluster->E() > 0 && eCell > 0){
 914
 915       w  = GetCellWeight(eCell,cluster->E());
 916
 917       etai=(Double_t)ieta;
 918       phii=(Double_t)iphi;
 919       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
 920     }
 921     else
 922       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
 923   }// cell loop
 924
 925   //Normalize to the weigth and set shower shape parameters
 926   if (wtot > 0 && nstat > 1) {
 927     d /= wtot ;
 928     dxx /= wtot ;
 929     dzz /= wtot ;
 930     dxz /= wtot ;
 931     dxx -= xmean * xmean ;
 932     dzz -= zmean * zmean ;
 933     dxz -= xmean * zmean ;
 934     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 935     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 936   }
 937   else{
 938     d=0. ;
 939     cluster->SetM20(0.) ;
 940     cluster->SetM02(0.) ;
 941   }
 942
 943   if (d>=0)
 944     cluster->SetDispersion(TMath::Sqrt(d)) ;
 945   else
 946     cluster->SetDispersion(0) ;
 947
 948 }
 949
 950 //__________________________________________________
 951 void AliEMCALRecoUtils::FindMatches(AliVEvent *event)
 952 {
 953   //This function should be called before the cluster loop
 954   //Before call this function, please recalculate the cluster positions
 955   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
 956   //Store matched cluster indexes and residuals
 957   //It only works with ESDs, not AODs
 958
 959   fMatchedClusterIndex->Reset();
 960   fResidualZ->Reset();
 961   fResidualR->Reset();
 962
 963   fMatchedClusterIndex->Set(100);
 964   fResidualZ->Set(100);
 965   fResidualR->Set(100);
 966
 967   Int_t    matched=0;
 968   Float_t  clsPos[3];
 969   Double_t trkPos[3];
 970   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
 971   {
 972     AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
 973     if(!track || !IsAccepted(track)) continue;
 974
 975     Float_t dRMax = fCutR, dZMax = fCutZ;
 976     Int_t index = -1;
 977     AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
 978     for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
 979     {
 980       AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
 981       if(!cluster->IsEMCAL()) continue;
 982       cluster->GetPosition(clsPos); //Has been recalculated
 983       if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
 984       emctrack->GetXYZ(trkPos);
 985       Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
 986       Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
 987
 988       if(tmpR<dRMax)
 989             {
 990               dRMax=tmpR;
 991               dZMax=tmpZ;
 992               index=icl;
 993             }
 994
 995     }//cluser loop
 996
 997     if(index>-1)
 998     {
 999       fMatchedClusterIndex->AddAt(index,matched);
1000       fResidualZ->AddAt(dZMax,matched);
1001       fResidualR->AddAt(dRMax,matched);
1002       matched++;
1003     }
1004     delete emctrack;
1005   }//track loop
1006   fMatchedClusterIndex->Set(matched);
1007   fResidualZ->Set(matched);
1008   fResidualR->Set(matched);
1009
1010   //printf("Number of matched pairs: %d\n",matched);
1011 }
1012
1013 //__________________________________________________
1014 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1015 {
1016   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1017   //Get the residuals dR and dZ for this cluster
1018   //It only works with ESDs, not AODs
1019
1020   if( FindMatchedPos(index) >= 999 )
1021   {
1022     AliDebug(2,"No matched tracks found!\n");
1023     dR=999.;
1024     dZ=999.;
1025     return;
1026   }
1027   dR = fResidualR->At(FindMatchedPos(index));
1028   dZ = fResidualZ->At(FindMatchedPos(index));
1029   //printf("dR %f, dZ %f\n",dR, dZ);
1030 }
1031
1032 //__________________________________________________
1033 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1034 {
1035   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1036   //Returns if cluster has a match
1037   if(FindMatchedPos(index) < 999)
1038     return kTRUE;
1039   else
1040     return kFALSE;
1041 }
1042 //__________________________________________________
1043 UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
1044 {
1045   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1046   //Returns the position of the match in the fMatchedClusterIndex array
1047   Float_t tmpR = fCutR;
1048   UInt_t pos = 999;
1049
1050   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1051     {
1052       if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1053         {
1054           pos=i;
1055           tmpR=fResidualR->At(i);
1056         }
1057       AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1058     }
1059   return pos;
1060 }
1061
1062 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1063 {
1064   // Given a esd track, return whether the track survive all the cuts
1065
1066   // The different quality parameter are first
1067   // retrieved from the track. then it is found out what cuts the
1068   // track did not survive and finally the cuts are imposed.
1069
1070   UInt_t status = esdTrack->GetStatus();
1071
1072   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1073   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1074
1075   Float_t chi2PerClusterITS = -1;
1076   Float_t chi2PerClusterTPC = -1;
1077   if (nClustersITS!=0)
1078     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1079   if (nClustersTPC!=0)
1080     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1081
1082
1083   //DCA cuts
1084   Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1085   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1086   SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1087
1088
1089   Float_t b[2];
1090   Float_t bCov[3];
1091   esdTrack->GetImpactParameters(b,bCov);
1092   if (bCov[0]<=0 || bCov[2]<=0) {
1093     AliDebug(1, "Estimated b resolution lower or equal zero!");
1094     bCov[0]=0; bCov[2]=0;
1095   }
1096
1097   Float_t dcaToVertexXY = b[0];
1098   Float_t dcaToVertexZ = b[1];
1099   Float_t dcaToVertex = -1;
1100
1101   if (fCutDCAToVertex2D)
1102     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1103   else
1104     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1105
1106   // cut the track?
1107
1108   Bool_t cuts[kNCuts];
1109   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1110
1111   // track quality cuts
1112   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1113     cuts[0]=kTRUE;
1114   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1115     cuts[1]=kTRUE;
1116   if (nClustersTPC<fCutMinNClusterTPC)
1117     cuts[2]=kTRUE;
1118   if (nClustersITS<fCutMinNClusterITS)
1119     cuts[3]=kTRUE;
1120   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1121     cuts[4]=kTRUE;
1122   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1123     cuts[5]=kTRUE;
1124   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1125     cuts[6]=kTRUE;
1126   if (fCutDCAToVertex2D && dcaToVertex > 1)
1127     cuts[7] = kTRUE;
1128   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1129     cuts[8] = kTRUE;
1130   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1131     cuts[9] = kTRUE;
1132
1133   //Require at least one SPD point + anything else in ITS
1134   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1135     cuts[10] = kTRUE;
1136
1137   Bool_t cut=kFALSE;
1138   for (Int_t i=0; i<kNCuts; i++)
1139     if (cuts[i]) {cut = kTRUE;}
1140
1141     // cut the track
1142   if (cut)
1143     return kFALSE;
1144   else
1145     return kTRUE;
1146 }
1147 //__________________________________________________
1148 void AliEMCALRecoUtils::InitTrackCuts()
1149 {
1150   //Intilize the track cut criteria
1151   //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1152   //Also you can customize the cuts using the setters
1153
1154   //TPC
1155   SetMinNClustersTPC(70);
1156   SetMaxChi2PerClusterTPC(4);
1157   SetAcceptKinkDaughters(kFALSE);
1158   SetRequireTPCRefit(kTRUE);
1159
1160   //ITS
1161   SetRequireITSRefit(kTRUE);
1162   SetMaxDCAToVertexZ(2);
1163   SetDCAToVertex2D(kFALSE);
1164   SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1165   SetMinNClustersITS();
1166 }
1167
1168 //__________________________________________________
1169 void AliEMCALRecoUtils::Print(const Option_t *) const
1170 {
1171   // Print Parameters
1172
1173   printf("AliEMCALRecoUtils Settings: \n");
1174   printf("Misalignment shifts\n");
1175   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,
1176                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1177                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1178   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1179   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1180
1181   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1182
1183   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1184
1185   printf("Track cuts: \n");
1186   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1187   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1188   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1189   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1190   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1191
1192
1193 }
1194
1195 //__________________________________________________
1196 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1197   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1198   //Do it only once and only if it is requested
1199
1200   if(!fUseTimeCorrectionFactors) return;
1201   if(fTimeCorrectionFactorsSet)  return;
1202
1203   printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1204
1205   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
1206   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1207
1208   SwitchOnRecalibration();
1209   for(Int_t ism = 0; ism < 4; ism++){
1210     for(Int_t icol = 0; icol < 48; icol++){
1211       for(Int_t irow = 0; irow < 24; irow++){
1212         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1213         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1214         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1215         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1216         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1217         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1218       }
1219     }
1220   }
1221    fTimeCorrectionFactorsSet = kTRUE;
1222 }
1223