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