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       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
 328       //Double_t par0 = 1.001;
 329       //Double_t par1 = -0.01264;
 330       //Double_t par2 = -0.03632;
 331       //Double_t par3 = 0.1798;
 332       //Double_t par4 = -0.522;
 333        energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 334                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 335                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 336       break;
 337
 338     case kPi0GammaGamma:
 339
 340       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 341       //Double_t par0 = 0.1457;
 342       //Double_t par1 = -0.02024;
 343       //Double_t par2 = 1.046;
 344       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 345       break;
 346
 347     case kPi0GammaConversion:
 348
 349       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 350       //Double_t C = 0.139393/0.1349766;
 351       //Double_t a = 0.0566186;
 352       //Double_t b = 0.982133;
 353       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 354
 355       break;
 356
 357     case kNoCorrection:
 358       AliDebug(2,"No correction on the energy\n");
 359       break;
 360
 361   }
 362
 363   return energy;
 364
 365 }
 366
 367 //__________________________________________________
 368 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
 369 {
 370   //Calculate shower depth for a given cluster energy and particle type
 371
 372   // parameters
 373   Float_t x0    = 1.31;
 374   Float_t ecr   = 8;
 375   Float_t depth = 0;
 376
 377   switch ( iParticle )
 378   {
 379     case kPhoton:
 380       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 381       break;
 382
 383     case kElectron:
 384       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
 385       break;
 386
 387     case kHadron:
 388       // hadron
 389       // boxes anc. here
 390       if(gGeoManager){
 391         gGeoManager->cd("ALIC_1/XEN1_1");
 392         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
 393         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
 394         if(geoSM){
 395           TGeoVolume      *geoSMVol   = geoSM->GetVolume();
 396           TGeoShape       *geoSMShape = geoSMVol->GetShape();
 397           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
 398           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
 399           else AliFatal("Null GEANT box");
 400         }else AliFatal("NULL  GEANT node matrix");
 401       }
 402       else{//electron
 403         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
 404       }
 405
 406       break;
 407
 408     default://photon
 409       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 410   }
 411
 412   return depth;
 413
 414 }
 415
 416 //__________________________________________________
 417 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
 418                                          Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
 419 {
 420   //For a given CaloCluster gets the absId of the cell
 421   //with maximum energy deposit.
 422
 423   Double_t eMax        = -1.;
 424   Double_t eCell       = -1.;
 425   Float_t  fraction    = 1.;
 426   Float_t  recalFactor = 1.;
 427   Int_t    cellAbsId   = -1 ;
 428
 429   Int_t iTower  = -1;
 430   Int_t iIphi   = -1;
 431   Int_t iIeta   = -1;
 432   Int_t iSupMod0= -1;
 433         //printf("---Max?\n");
 434   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 435     cellAbsId = clu->GetCellAbsId(iDig);
 436     fraction  = clu->GetCellAmplitudeFraction(iDig);
 437     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
 438     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 439     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 440     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 441     if(iDig==0) iSupMod0=iSupMod;
 442     else if(iSupMod0!=iSupMod) {
 443       shared = kTRUE;
 444       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
 445     }
 446     if(IsRecalibrationOn()) {
 447       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 448     }
 449     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 450     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
 451     if(eCell > eMax)  {
 452       eMax  = eCell;
 453       absId = cellAbsId;
 454       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 455     }
 456   }// cell loop
 457
 458   //Get from the absid the supermodule, tower and eta/phi numbers
 459   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 460   //Gives SuperModule and Tower numbers
 461   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 462                                          iIphi, iIeta,iphi,ieta);
 463   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
 464   //printf("Max end---\n");
 465
 466 }
 467
 468 //________________________________________________________________
 469 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
 470         //Init EMCAL recalibration factors
 471         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
 472         //In order to avoid rewriting the same histograms
 473         Bool_t oldStatus = TH1::AddDirectoryStatus();
 474         TH1::AddDirectory(kFALSE);
 475
 476         fEMCALRecalibrationFactors = new TObjArray(10);
 477         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));
 478         //Init the histograms with 1
 479         for (Int_t sm = 0; sm < 12; sm++) {
 480                 for (Int_t i = 0; i < 48; i++) {
 481                         for (Int_t j = 0; j < 24; j++) {
 482                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
 483                         }
 484                 }
 485         }
 486         fEMCALRecalibrationFactors->SetOwner(kTRUE);
 487         fEMCALRecalibrationFactors->Compress();
 488
 489         //In order to avoid rewriting the same histograms
 490         TH1::AddDirectory(oldStatus);
 491 }
 492
 493
 494 //________________________________________________________________
 495 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 496         //Init EMCAL bad channels map
 497         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
 498         //In order to avoid rewriting the same histograms
 499         Bool_t oldStatus = TH1::AddDirectoryStatus();
 500         TH1::AddDirectory(kFALSE);
 501
 502         fEMCALBadChannelMap = new TObjArray(10);
 503         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
 504         for (int i = 0; i < 12; i++) {
 505                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
 506         }
 507
 508         //delete hTemp;
 509
 510         fEMCALBadChannelMap->SetOwner(kTRUE);
 511         fEMCALBadChannelMap->Compress();
 512
 513         //In order to avoid rewriting the same histograms
 514         TH1::AddDirectory(oldStatus);
 515 }
 516
 517 //________________________________________________________________
 518 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
 519         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 520
 521         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 522         UShort_t * index    = cluster->GetCellsAbsId() ;
 523         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 524         Int_t ncells = cluster->GetNCells();
 525
 526         //Initialize some used variables
 527         Float_t energy = 0;
 528         Int_t absId    = -1;
 529   Int_t icol = -1, irow = -1, imod=1;
 530         Float_t factor = 1, frac = 0;
 531
 532         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 533         for(Int_t icell = 0; icell < ncells; icell++){
 534                 absId = index[icell];
 535                 frac =  fraction[icell];
 536                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
 537                 Int_t iTower = -1, iIphi = -1, iIeta = -1;
 538                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 539                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
 540                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 541                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
 542     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
 543              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
 544
 545                 energy += cells->GetCellAmplitude(absId)*factor*frac;
 546         }
 547
 548
 549                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
 550
 551         cluster->SetE(energy);
 552
 553 }
 554
 555
 556 //__________________________________________________
 557 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 558 {
 559   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 560
 561   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
 562   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
 563   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
 564
 565 }
 566
 567 //__________________________________________________
 568 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 569 {
 570   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 571   // The algorithm is a copy of what is done in AliEMCALRecPoint
 572
 573   Double_t eCell       = 0.;
 574   Float_t  fraction    = 1.;
 575   Float_t  recalFactor = 1.;
 576
 577   Int_t    absId   = -1;
 578   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
 579   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
 580   Float_t  weight = 0.,  totalWeight=0.;
 581   Float_t  newPos[3] = {0,0,0};
 582   Double_t pLocal[3], pGlobal[3];
 583   Bool_t shared = kFALSE;
 584
 585   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
 586   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
 587   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
 588
 589   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
 590
 591   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 592     absId = clu->GetCellAbsId(iDig);
 593     fraction  = clu->GetCellAmplitudeFraction(iDig);
 594     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 595     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
 596     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 597
 598     if(IsRecalibrationOn()) {
 599       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
 600     }
 601     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 602
 603     weight = GetCellWeight(eCell,clEnergy);
 604     //printf("cell energy %f, weight %f\n",eCell,weight);
 605     totalWeight += weight;
 606     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
 607     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
 608     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
 609     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
 610
 611     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
 612
 613   }// cell loop
 614
 615   if(totalWeight>0){
 616     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
 617   }
 618
 619   //Float_t pos[]={0,0,0};
 620   //clu->GetPosition(pos);
 621   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
 622   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 623
 624         if(iSupModMax > 1) {//sector 1
 625           newPos[0] +=fMisalTransShift[3];//-=3.093;
 626           newPos[1] +=fMisalTransShift[4];//+=6.82;
 627           newPos[2] +=fMisalTransShift[5];//+=1.635;
 628     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
 629
 630         }
 631         else {//sector 0
 632           newPos[0] +=fMisalTransShift[0];//+=1.134;
 633           newPos[1] +=fMisalTransShift[1];//+=8.2;
 634           newPos[2] +=fMisalTransShift[2];//+=1.197;
 635     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
 636
 637         }
 638   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 639
 640   clu->SetPosition(newPos);
 641
 642 }
 643
 644 //__________________________________________________
 645 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 646 {
 647   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 648   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
 649
 650   Double_t eCell       = 1.;
 651   Float_t  fraction    = 1.;
 652   Float_t  recalFactor = 1.;
 653
 654   Int_t absId   = -1;
 655   Int_t iTower  = -1;
 656   Int_t iIphi   = -1, iIeta   = -1;
 657         Int_t iSupMod = -1, iSupModMax = -1;
 658   Int_t iphi = -1, ieta =-1;
 659   Bool_t shared = kFALSE;
 660
 661   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
 662   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
 663   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
 664
 665   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
 666   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
 667   Int_t startingSM = -1;
 668
 669   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 670     absId = clu->GetCellAbsId(iDig);
 671     fraction  = clu->GetCellAmplitudeFraction(iDig);
 672     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 673     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 674     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 675
 676     if     (iDig==0)  startingSM = iSupMod;
 677     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 678
 679     eCell  = cells->GetCellAmplitude(absId);
 680
 681     if(IsRecalibrationOn()) {
 682       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 683     }
 684     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 685
 686     weight = GetCellWeight(eCell,clEnergy);
 687     if(weight < 0) weight = 0;
 688     totalWeight += weight;
 689     weightedCol += ieta*weight;
 690     weightedRow += iphi*weight;
 691
 692     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
 693
 694     }// cell loop
 695
 696   Float_t xyzNew[]={0.,0.,0.};
 697   if(areInSameSM == kTRUE) {
 698     //printf("In Same SM\n");
 699     weightedCol = weightedCol/totalWeight;
 700     weightedRow = weightedRow/totalWeight;
 701     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 702   }
 703   else {
 704     //printf("In Different SM\n");
 705     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 706   }
 707
 708   clu->SetPosition(xyzNew);
 709
 710 }
 711
 712 //____________________________________________________________________________
 713 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
 714
 715   //re-evaluate distance to bad channel with updated bad map
 716
 717   if(!fRecalDistToBadChannels) return;
 718
 719         //Get channels map of the supermodule where the cluster is.
 720   Int_t absIdMax        = -1, iSupMod =-1, icolM = -1, irowM = -1;
 721   Bool_t shared = kFALSE;
 722   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
 723   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
 724
 725   Int_t dRrow, dRcol;
 726         Float_t  minDist = 10000.;
 727         Float_t  dist    = 0.;
 728
 729   //Loop on tower status map
 730         for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
 731                 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
 732                         //Check if tower is bad.
 733                         if(hMap->GetBinContent(icol,irow)==0) continue;
 734       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
 735         //     iSupMod,icol, irow, icolM,irowM);
 736
 737       dRrow=TMath::Abs(irowM-irow);
 738       dRcol=TMath::Abs(icolM-icol);
 739       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
 740                         if(dist < minDist){
 741         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
 742         minDist = dist;
 743       }
 744
 745                 }
 746         }
 747
 748         //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
 749         if (shared) {
 750                 TH2D* hMap2 = 0;
 751                 Int_t iSupMod2 = -1;
 752
 753                 //The only possible combinations are (0,1), (2,3) ... (8,9)
 754                 if(iSupMod%2) iSupMod2 = iSupMod-1;
 755                 else          iSupMod2 = iSupMod+1;
 756                 hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
 757
 758                 //Loop on tower status map of second super module
 759                 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
 760                         for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
 761                                 //Check if tower is bad.
 762                                 if(hMap2->GetBinContent(icol,irow)==0) continue;
 763                                 //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",
 764           //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
 765
 766         dRrow=TMath::Abs(irow-irowM);
 767
 768         if(iSupMod%2) {
 769                                   dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
 770                                 }
 771         else {
 772           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
 773                                 }
 774
 775                                 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
 776         if(dist < minDist) minDist = dist;
 777
 778                         }
 779                 }
 780
 781         }// shared cluster in 2 SuperModules
 782
 783   AliDebug(2,Form("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f\n",iSupMod, icolM, irowM, minDist));
 784         cluster->SetDistanceToBadChannel(minDist);
 785
 786 }
 787
 788 //____________________________________________________________________________
 789 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
 790
 791   //re-evaluate identification parameters with bayesian
 792
 793         if ( cluster->GetM02() != 0)
 794     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
 795
 796   Float_t pidlist[AliPID::kSPECIESN+1];
 797         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
 798
 799   cluster->SetPID(pidlist);
 800
 801 }
 802
 803 //____________________________________________________________________________
 804 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
 805 {
 806   // Calculates new center of gravity in the local EMCAL-module coordinates
 807   // and tranfers into global ALICE coordinates
 808   // Calculates Dispersion and main axis
 809
 810   Int_t nstat  = 0;
 811   Float_t wtot = 0. ;
 812   Double_t eCell       = 0.;
 813   Float_t  fraction    = 1.;
 814   Float_t  recalFactor = 1.;
 815
 816   Int_t iSupMod = -1;
 817   Int_t iTower  = -1;
 818   Int_t iIphi   = -1;
 819   Int_t iIeta   = -1;
 820   Int_t iphi    = -1;
 821   Int_t ieta    = -1;
 822   Double_t etai = -1.;
 823   Double_t phii = -1.;
 824
 825   Double_t w     = 0.;
 826   Double_t d     = 0.;
 827   Double_t dxx   = 0.;
 828   Double_t dzz   = 0.;
 829   Double_t dxz   = 0.;
 830   Double_t xmean = 0.;
 831   Double_t zmean = 0.;
 832
 833   //Loop on cells
 834   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
 835
 836     //Get from the absid the supermodule, tower and eta/phi numbers
 837     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 838     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 839
 840     //Get the cell energy, if recalibration is on, apply factors
 841     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
 842     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 843     if(IsRecalibrationOn()) {
 844       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 845     }
 846     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
 847
 848     if(cluster->E() > 0 && eCell > 0){
 849
 850       w  = GetCellWeight(eCell,cluster->E());
 851
 852       etai=(Double_t)ieta;
 853       phii=(Double_t)iphi;
 854       if(w > 0.0) {
 855         wtot += w ;
 856         nstat++;
 857         //Shower shape
 858         dxx  += w * etai * etai ;
 859         xmean+= w * etai ;
 860         dzz  += w * phii * phii ;
 861         zmean+= w * phii ;
 862         dxz  += w * etai * phii ;
 863       }
 864     }
 865     else
 866       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
 867   }//cell loop
 868
 869   //Normalize to the weight
 870   if (wtot > 0) {
 871     xmean /= wtot ;
 872     zmean /= wtot ;
 873   }
 874   else
 875     AliError(Form("Wrong weight %f\n", wtot));
 876
 877   //Calculate dispersion
 878   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
 879
 880     //Get from the absid the supermodule, tower and eta/phi numbers
 881     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 882     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 883
 884     //Get the cell energy, if recalibration is on, apply factors
 885     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
 886     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 887     if(IsRecalibrationOn()) {
 888       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 889     }
 890     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
 891
 892     if(cluster->E() > 0 && eCell > 0){
 893
 894       w  = GetCellWeight(eCell,cluster->E());
 895
 896       etai=(Double_t)ieta;
 897       phii=(Double_t)iphi;
 898       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
 899     }
 900     else
 901       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
 902   }// cell loop
 903
 904   //Normalize to the weigth and set shower shape parameters
 905   if (wtot > 0 && nstat > 1) {
 906     d /= wtot ;
 907     dxx /= wtot ;
 908     dzz /= wtot ;
 909     dxz /= wtot ;
 910     dxx -= xmean * xmean ;
 911     dzz -= zmean * zmean ;
 912     dxz -= xmean * zmean ;
 913     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 914     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 915   }
 916   else{
 917     d=0. ;
 918     cluster->SetM20(0.) ;
 919     cluster->SetM02(0.) ;
 920   }
 921
 922   if (d>=0)
 923     cluster->SetDispersion(TMath::Sqrt(d)) ;
 924   else
 925     cluster->SetDispersion(0) ;
 926
 927 }
 928
 929 //__________________________________________________
 930 void AliEMCALRecoUtils::FindMatches(AliVEvent *event)
 931 {
 932   //This function should be called before the cluster loop
 933   //Before call this function, please recalculate the cluster positions
 934   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
 935   //Store matched cluster indexes and residuals
 936   //It only works with ESDs, not AODs
 937
 938   fMatchedClusterIndex->Reset();
 939   fResidualZ->Reset();
 940   fResidualR->Reset();
 941
 942   fMatchedClusterIndex->Set(100);
 943   fResidualZ->Set(100);
 944   fResidualR->Set(100);
 945
 946   Int_t    matched=0;
 947   Float_t  clsPos[3];
 948   Double_t trkPos[3];
 949   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
 950   {
 951     AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
 952     if(!track || !IsAccepted(track)) continue;
 953
 954     Float_t dRMax = fCutR, dZMax = fCutZ;
 955     Int_t index = -1;
 956     AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
 957     for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
 958     {
 959       AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
 960       if(!cluster->IsEMCAL()) continue;
 961       cluster->GetPosition(clsPos); //Has been recalculated
 962       if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
 963       emctrack->GetXYZ(trkPos);
 964       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) );
 965       Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
 966
 967       if(tmpR<dRMax)
 968             {
 969               dRMax=tmpR;
 970               dZMax=tmpZ;
 971               index=icl;
 972             }
 973
 974     }//cluser loop
 975
 976     if(index>-1)
 977     {
 978       fMatchedClusterIndex->AddAt(index,matched);
 979       fResidualZ->AddAt(dZMax,matched);
 980       fResidualR->AddAt(dRMax,matched);
 981       matched++;
 982     }
 983     delete emctrack;
 984   }//track loop
 985   fMatchedClusterIndex->Set(matched);
 986   fResidualZ->Set(matched);
 987   fResidualR->Set(matched);
 988
 989   //printf("Number of matched pairs: %d\n",matched);
 990 }
 991
 992 //__________________________________________________
 993 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
 994 {
 995   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
 996   //Get the residuals dR and dZ for this cluster
 997   //It only works with ESDs, not AODs
 998
 999   if( FindMatchedPos(index) >= 999 )
1000   {
1001     AliDebug(2,"No matched tracks found!\n");
1002     dR=999.;
1003     dZ=999.;
1004     return;
1005   }
1006   dR = fResidualR->At(FindMatchedPos(index));
1007   dZ = fResidualZ->At(FindMatchedPos(index));
1008 }
1009
1010 //__________________________________________________
1011 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1012 {
1013   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1014   //Returns if cluster has a match
1015   if(FindMatchedPos(index) < 999)
1016     return kTRUE;
1017   else
1018     return kFALSE;
1019 }
1020 //__________________________________________________
1021 UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
1022 {
1023   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1024   //Returns the position of the match in the fMatchedClusterIndex array
1025   Float_t tmpR = fCutR;
1026   UInt_t pos = 999;
1027
1028   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1029     {
1030       if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1031         {
1032           pos=i;
1033           tmpR=fResidualR->At(i);
1034         }
1035       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)));
1036     }
1037   return pos;
1038 }
1039
1040 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1041 {
1042   // Given a esd track, return whether the track survive all the cuts
1043
1044   // The different quality parameter are first
1045   // retrieved from the track. then it is found out what cuts the
1046   // track did not survive and finally the cuts are imposed.
1047
1048   UInt_t status = esdTrack->GetStatus();
1049
1050   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1051   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1052
1053   Float_t chi2PerClusterITS = -1;
1054   Float_t chi2PerClusterTPC = -1;
1055   if (nClustersITS!=0)
1056     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1057   if (nClustersTPC!=0)
1058     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1059
1060
1061   //DCA cuts
1062   Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1063   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1064   SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1065
1066
1067   Float_t b[2];
1068   Float_t bCov[3];
1069   esdTrack->GetImpactParameters(b,bCov);
1070   if (bCov[0]<=0 || bCov[2]<=0) {
1071     AliDebug(1, "Estimated b resolution lower or equal zero!");
1072     bCov[0]=0; bCov[2]=0;
1073   }
1074
1075   Float_t dcaToVertexXY = b[0];
1076   Float_t dcaToVertexZ = b[1];
1077   Float_t dcaToVertex = -1;
1078
1079   if (fCutDCAToVertex2D)
1080     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1081   else
1082     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1083
1084   // cut the track?
1085
1086   Bool_t cuts[kNCuts];
1087   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1088
1089   // track quality cuts
1090   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1091     cuts[0]=kTRUE;
1092   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1093     cuts[1]=kTRUE;
1094   if (nClustersTPC<fCutMinNClusterTPC)
1095     cuts[2]=kTRUE;
1096   if (nClustersITS<fCutMinNClusterITS)
1097     cuts[3]=kTRUE;
1098   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1099     cuts[4]=kTRUE;
1100   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1101     cuts[5]=kTRUE;
1102   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1103     cuts[6]=kTRUE;
1104   if (fCutDCAToVertex2D && dcaToVertex > 1)
1105     cuts[7] = kTRUE;
1106   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1107     cuts[8] = kTRUE;
1108   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1109     cuts[9] = kTRUE;
1110
1111   //Require at least one SPD point + anything else in ITS
1112   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1113     cuts[10] = kTRUE;
1114
1115   Bool_t cut=kFALSE;
1116   for (Int_t i=0; i<kNCuts; i++)
1117     if (cuts[i]) {cut = kTRUE;}
1118
1119     // cut the track
1120   if (cut)
1121     return kFALSE;
1122   else
1123     return kTRUE;
1124 }
1125 //__________________________________________________
1126 void AliEMCALRecoUtils::InitTrackCuts()
1127 {
1128   //Intilize the track cut criteria
1129   //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1130   //Also you can customize the cuts using the setters
1131
1132   //TPC
1133   SetMinNClustersTPC(70);
1134   SetMaxChi2PerClusterTPC(4);
1135   SetAcceptKinkDaughters(kFALSE);
1136   SetRequireTPCRefit(kTRUE);
1137
1138   //ITS
1139   SetRequireITSRefit(kTRUE);
1140   SetMaxDCAToVertexZ(2);
1141   SetDCAToVertex2D(kFALSE);
1142 }
1143
1144 //__________________________________________________
1145 void AliEMCALRecoUtils::Print(const Option_t *) const
1146 {
1147   // Print Parameters
1148
1149   printf("AliEMCALRecoUtils Settings: \n");
1150   printf("Misalignment shifts\n");
1151   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,
1152                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1153                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1154   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1155   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1156
1157   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1158
1159   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1160
1161   printf("Track cuts: \n");
1162   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1163   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1164   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1165   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1166   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1167
1168
1169 }