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