EMCAL/AliEMCALRecoUtils.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 // Class AliEMCALRecoUtils
  21 // Some utilities to recalculate the cluster position or energy linearity
  22 //
  23 //
  24 // Author:  Gustavo Conesa (LPSC- Grenoble)
  25 //          Track matching part: Rongrong Ma (Yale)
  26
  27 ///////////////////////////////////////////////////////////////////////////////
  28 // --- standard c ---
  29
  30 // standard C++ includes
  31 //#include <Riostream.h>
  32
  33 // ROOT includes
  34 #include <TGeoManager.h>
  35 #include <TGeoMatrix.h>
  36 #include <TGeoBBox.h>
  37 #include <TH2F.h>
  38 #include <TArrayI.h>
  39 #include <TArrayF.h>
  40 #include <TObjArray.h>
  41
  42 // STEER includes
  43 #include "AliVCluster.h"
  44 #include "AliVCaloCells.h"
  45 #include "AliLog.h"
  46 #include "AliPID.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAODEvent.h"
  49 #include "AliESDtrack.h"
  50 #include "AliAODTrack.h"
  51 #include "AliExternalTrackParam.h"
  52 #include "AliESDfriendTrack.h"
  53 #include "AliTrackerBase.h"
  54
  55 // EMCAL includes
  56 #include "AliEMCALRecoUtils.h"
  57 #include "AliEMCALGeometry.h"
  58 #include "AliTrackerBase.h"
  59 #include "AliEMCALPIDUtils.h"
  60
  61 ClassImp(AliEMCALRecoUtils)
  62
  63 //_____________________________________
  64 AliEMCALRecoUtils::AliEMCALRecoUtils():
  65   fParticleType(0),                       fPosAlgo(0),                            fW0(0),
  66   fNonLinearityFunction(0),               fNonLinearThreshold(0),
  67   fSmearClusterEnergy(kFALSE),            fRandom(),
  68   fCellsRecalibrated(kFALSE),             fRecalibration(kFALSE),                 fEMCALRecalibrationFactors(),
  69   fTimeRecalibration(kFALSE),             fEMCALTimeRecalibrationFactors(),       fUseRunCorrectionFactors(kFALSE),
  70   fRemoveBadChannels(kFALSE),             fRecalDistToBadChannels(kFALSE),        fEMCALBadChannelMap(),
  71   fNCellsFromEMCALBorder(0),              fNoEMCALBorderAtEta0(kTRUE),
  72   fRejectExoticCluster(kFALSE),           fRejectExoticCells(kFALSE),
  73   fExoticCellFraction(0),                 fExoticCellDiffTime(0),                 fExoticCellMinAmplitude(0),
  74   fPIDUtils(),                            fAODFilterMask(0),                      fAODHybridTracks(0),
  75   fMatchedTrackIndex(0x0),                fMatchedClusterIndex(0x0),
  76   fResidualEta(0x0), fResidualPhi(0x0),   fCutEtaPhiSum(kFALSE),                  fCutEtaPhiSeparate(kFALSE),
  77   fCutR(0),                               fCutEta(0),                             fCutPhi(0),
  78   fClusterWindow(0),                      fMass(0),
  79   fStepSurface(0),                        fStepCluster(0),
  80   fITSTrackSA(kFALSE),
  81   fTrackCutsType(0),                      fCutMinTrackPt(0),                      fCutMinNClusterTPC(0),
  82   fCutMinNClusterITS(0),                  fCutMaxChi2PerClusterTPC(0),            fCutMaxChi2PerClusterITS(0),
  83   fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
  84   fCutMaxDCAToVertexXY(0),                fCutMaxDCAToVertexZ(0),                 fCutDCAToVertex2D(kFALSE),
  85   fCutRequireITSStandAlone(kFALSE),       fCutRequireITSpureSA(kFALSE)
  86 {
  87 //
  88   // Constructor.
  89   // Initialize all constant values which have to be used
  90   // during Reco algorithm execution
  91   //
  92
  93   // Init parameters
  94   InitParameters();
  95
  96   //Track matching
  97   fMatchedTrackIndex     = new TArrayI();
  98   fMatchedClusterIndex   = new TArrayI();
  99   fResidualPhi           = new TArrayF();
 100   fResidualEta           = new TArrayF();
 101   fPIDUtils              = new AliEMCALPIDUtils();
 102
 103 }
 104
 105 //______________________________________________________________________
 106 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
 107 : TNamed(reco),
 108   fParticleType(reco.fParticleType),                         fPosAlgo(reco.fPosAlgo),     fW0(reco.fW0),
 109   fNonLinearityFunction(reco.fNonLinearityFunction),         fNonLinearThreshold(reco.fNonLinearThreshold),
 110   fSmearClusterEnergy(reco.fSmearClusterEnergy),             fRandom(),
 111   fCellsRecalibrated(reco.fCellsRecalibrated),
 112   fRecalibration(reco.fRecalibration),                       fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
 113   fTimeRecalibration(reco.fTimeRecalibration),               fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
 114   fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),
 115   fRemoveBadChannels(reco.fRemoveBadChannels),               fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
 116   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
 117   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
 118   fRejectExoticCluster(reco.fRejectExoticCluster),           fRejectExoticCells(reco.fRejectExoticCells),
 119   fExoticCellFraction(reco.fExoticCellFraction),             fExoticCellDiffTime(reco.fExoticCellDiffTime),
 120   fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
 121   fPIDUtils(reco.fPIDUtils),                                 fAODFilterMask(reco.fAODFilterMask),
 122   fAODHybridTracks(reco.fAODHybridTracks),
 123   fMatchedTrackIndex(  reco.fMatchedTrackIndex?  new TArrayI(*reco.fMatchedTrackIndex):0x0),
 124   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
 125   fResidualEta(        reco.fResidualEta?        new TArrayF(*reco.fResidualEta):0x0),
 126   fResidualPhi(        reco.fResidualPhi?        new TArrayF(*reco.fResidualPhi):0x0),
 127   fCutEtaPhiSum(reco.fCutEtaPhiSum),                         fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
 128   fCutR(reco.fCutR),        fCutEta(reco.fCutEta),           fCutPhi(reco.fCutPhi),
 129   fClusterWindow(reco.fClusterWindow),
 130   fMass(reco.fMass),        fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
 131   fITSTrackSA(reco.fITSTrackSA),
 132   fTrackCutsType(reco.fTrackCutsType),                       fCutMinTrackPt(reco.fCutMinTrackPt),
 133   fCutMinNClusterTPC(reco.fCutMinNClusterTPC),               fCutMinNClusterITS(reco.fCutMinNClusterITS),
 134   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC),   fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
 135   fCutRequireTPCRefit(reco.fCutRequireTPCRefit),             fCutRequireITSRefit(reco.fCutRequireITSRefit),
 136   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),     fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
 137   fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),             fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
 138   fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone),   fCutRequireITSpureSA(reco.fCutRequireITSpureSA)
 139 {
 140   //Copy ctor
 141
 142   for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i]      = reco.fMisalRotShift[i]      ;
 143                                    fMisalTransShift[i]    = reco.fMisalTransShift[i]    ; }
 144   for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
 145   for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }
 146
 147 }
 148
 149
 150 //______________________________________________________________________
 151 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
 152 {
 153   //Assignment operator
 154
 155   if(this == &reco)return *this;
 156   ((TNamed *)this)->operator=(reco);
 157
 158   for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i]    = reco.fMisalTransShift[i]    ;
 159                                    fMisalRotShift[i]      = reco.fMisalRotShift[i]      ; }
 160   for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
 161   for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }
 162
 163   fParticleType              = reco.fParticleType;
 164   fPosAlgo                   = reco.fPosAlgo;
 165   fW0                        = reco.fW0;
 166
 167   fNonLinearityFunction      = reco.fNonLinearityFunction;
 168   fNonLinearThreshold        = reco.fNonLinearThreshold;
 169   fSmearClusterEnergy        = reco.fSmearClusterEnergy;
 170
 171   fCellsRecalibrated         = reco.fCellsRecalibrated;
 172   fRecalibration             = reco.fRecalibration;
 173   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
 174
 175   fTimeRecalibration             = reco.fTimeRecalibration;
 176   fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
 177
 178   fUseRunCorrectionFactors   = reco.fUseRunCorrectionFactors;
 179
 180   fRemoveBadChannels         = reco.fRemoveBadChannels;
 181   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
 182   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
 183
 184   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
 185   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
 186
 187   fRejectExoticCluster       = reco.fRejectExoticCluster;
 188   fRejectExoticCells         = reco.fRejectExoticCells;
 189   fExoticCellFraction        = reco.fExoticCellFraction;
 190   fExoticCellDiffTime        = reco.fExoticCellDiffTime;
 191   fExoticCellMinAmplitude    = reco.fExoticCellMinAmplitude;
 192
 193   fPIDUtils                  = reco.fPIDUtils;
 194
 195   fAODFilterMask             = reco.fAODFilterMask;
 196   fAODHybridTracks           = reco.fAODHybridTracks;
 197
 198   fCutEtaPhiSum              = reco.fCutEtaPhiSum;
 199   fCutEtaPhiSeparate         = reco.fCutEtaPhiSeparate;
 200   fCutR                      = reco.fCutR;
 201   fCutEta                    = reco.fCutEta;
 202   fCutPhi                    = reco.fCutPhi;
 203   fClusterWindow             = reco.fClusterWindow;
 204   fMass                      = reco.fMass;
 205   fStepSurface               = reco.fStepSurface;
 206   fStepCluster               = reco.fStepCluster;
 207   fITSTrackSA                = reco.fITSTrackSA;
 208
 209   fTrackCutsType             = reco.fTrackCutsType;
 210   fCutMinTrackPt             = reco.fCutMinTrackPt;
 211   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
 212   fCutMinNClusterITS         = reco.fCutMinNClusterITS;
 213   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
 214   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
 215   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
 216   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
 217   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
 218   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
 219   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
 220   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
 221   fCutRequireITSStandAlone   = reco.fCutRequireITSStandAlone;
 222   fCutRequireITSpureSA       = reco.fCutRequireITSpureSA;
 223   if(reco.fResidualEta)
 224   {
 225     // assign or copy construct
 226     if(fResidualEta)
 227     {
 228       *fResidualEta = *reco.fResidualEta;
 229     }
 230     else
 231     {
 232       fResidualEta = new TArrayF(*reco.fResidualEta);
 233     }
 234   }
 235   else
 236   {
 237     if(fResidualEta)delete fResidualEta;
 238     fResidualEta = 0;
 239   }
 240
 241   if(reco.fResidualPhi)
 242   {
 243     // assign or copy construct
 244     if(fResidualPhi)
 245     {
 246       *fResidualPhi = *reco.fResidualPhi;
 247     }
 248     else
 249     {
 250       fResidualPhi = new TArrayF(*reco.fResidualPhi);
 251     }
 252   }
 253   else
 254   {
 255     if(fResidualPhi)delete fResidualPhi;
 256     fResidualPhi = 0;
 257   }
 258
 259   if(reco.fMatchedTrackIndex)
 260   {
 261     // assign or copy construct
 262     if(fMatchedTrackIndex)
 263     {
 264       *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
 265     }
 266     else
 267     {
 268       fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
 269     }
 270   }
 271   else
 272   {
 273     if(fMatchedTrackIndex)delete fMatchedTrackIndex;
 274     fMatchedTrackIndex = 0;
 275   }
 276
 277   if(reco.fMatchedClusterIndex)
 278   {
 279     // assign or copy construct
 280     if(fMatchedClusterIndex)
 281     {
 282       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
 283     }
 284     else
 285     {
 286       fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
 287     }
 288   }
 289   else
 290   {
 291     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
 292     fMatchedClusterIndex = 0;
 293   }
 294
 295   return *this;
 296 }
 297
 298
 299 //_____________________________________
 300 AliEMCALRecoUtils::~AliEMCALRecoUtils()
 301 {
 302   //Destructor.
 303
 304   if(fEMCALRecalibrationFactors)
 305   {
 306     fEMCALRecalibrationFactors->Clear();
 307     delete fEMCALRecalibrationFactors;
 308   }
 309
 310   if(fEMCALTimeRecalibrationFactors)
 311   {
 312     fEMCALTimeRecalibrationFactors->Clear();
 313     delete fEMCALTimeRecalibrationFactors;
 314   }
 315
 316   if(fEMCALBadChannelMap)
 317   {
 318     fEMCALBadChannelMap->Clear();
 319     delete fEMCALBadChannelMap;
 320   }
 321
 322   delete fMatchedTrackIndex   ;
 323   delete fMatchedClusterIndex ;
 324   delete fResidualEta         ;
 325   delete fResidualPhi         ;
 326   delete fPIDUtils            ;
 327
 328   InitTrackCuts();
 329 }
 330
 331 //_______________________________________________________________________________
 332 Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
 333                                               Float_t  & amp,    Double_t & time,
 334                                               AliVCaloCells* cells)
 335 {
 336   // Reject cell if criteria not passed and calibrate it
 337
 338   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 339
 340   if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE;
 341
 342   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 343
 344   if(!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
 345   {
 346     // cell absID does not exist
 347     amp=0; time = 1.e9;
 348     return kFALSE;
 349   }
 350
 351   geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 352
 353   // Do not include bad channels found in analysis,
 354   if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi))
 355   {
 356     return kFALSE;
 357   }
 358
 359   //Recalibrate energy
 360   amp  = cells->GetCellAmplitude(absID);
 361   if(!fCellsRecalibrated && IsRecalibrationOn())
 362     amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 363
 364
 365   // Recalibrate time
 366   time = cells->GetCellTime(absID);
 367
 368   RecalibrateCellTime(absID,bc,time);
 369
 370   return kTRUE;
 371 }
 372
 373 //_____________________________________________________________________________
 374 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom,
 375                                                   const AliVCluster* cluster,
 376                                                   AliVCaloCells* cells)
 377 {
 378   // Given the list of AbsId of the cluster, get the maximum cell and
 379   // check if there are fNCellsFromBorder from the calorimeter border
 380
 381   if(!cluster)
 382   {
 383     AliInfo("Cluster pointer null!");
 384     return kFALSE;
 385   }
 386
 387   //If the distance to the border is 0 or negative just exit accept all clusters
 388   if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
 389
 390   Int_t absIdMax  = -1, iSM =-1, ieta = -1, iphi = -1;
 391   Bool_t shared = kFALSE;
 392   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
 393
 394   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
 395                   absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
 396
 397   if(absIdMax==-1) return kFALSE;
 398
 399   //Check if the cell is close to the borders:
 400   Bool_t okrow = kFALSE;
 401   Bool_t okcol = kFALSE;
 402
 403   if(iSM < 0 || iphi < 0 || ieta < 0 )
 404   {
 405     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
 406                   iSM,ieta,iphi));
 407   }
 408
 409   //Check rows/phi
 410   if(iSM < 10)
 411   {
 412     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
 413   }
 414   else if (iSM >=10 && ( ( geom->GetEMCGeometry()->GetGeoName()).Contains("12SMV1")))
 415   {
 416     if(iphi >= fNCellsFromEMCALBorder && iphi < 8-fNCellsFromEMCALBorder) okrow =kTRUE; //1/3 sm case
 417   }
 418   else
 419   {
 420     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; // half SM case
 421   }
 422
 423   //Check columns/eta
 424   if(!fNoEMCALBorderAtEta0)
 425   {
 426     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
 427   }
 428   else
 429   {
 430     if(iSM%2==0)
 431     {
 432       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;
 433     }
 434     else
 435     {
 436       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;
 437     }
 438   }//eta 0 not checked
 439
 440   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
 441                   fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
 442
 443   if (okcol && okrow)
 444   {
 445     //printf("Accept\n");
 446     return kTRUE;
 447   }
 448   else
 449   {
 450     //printf("Reject\n");
 451     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
 452     return kFALSE;
 453   }
 454
 455 }
 456
 457
 458 //_______________________________________________________________________________
 459 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom,
 460                                                     const UShort_t* cellList,
 461                                                     const Int_t nCells)
 462 {
 463   // Check that in the cluster cells, there is no bad channel of those stored
 464   // in fEMCALBadChannelMap or fPHOSBadChannelMap
 465
 466   if(!fRemoveBadChannels)  return kFALSE;
 467   if(!fEMCALBadChannelMap) return kFALSE;
 468
 469   Int_t icol = -1;
 470   Int_t irow = -1;
 471   Int_t imod = -1;
 472   for(Int_t iCell = 0; iCell<nCells; iCell++)
 473   {
 474     //Get the column and row
 475     Int_t iTower = -1, iIphi = -1, iIeta = -1;
 476     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
 477     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
 478     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 479     if(GetEMCALChannelStatus(imod, icol, irow))
 480     {
 481       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
 482       return kTRUE;
 483     }
 484
 485   }// cell cluster loop
 486
 487   return kFALSE;
 488 }
 489
 490
 491 //___________________________________________________________________________
 492 Float_t AliEMCALRecoUtils::GetECross(const Int_t absID, const Double_t tcell,
 493                                      AliVCaloCells* cells, const Int_t bc)
 494 {
 495   //Calculate the energy in the cross around the energy given cell
 496
 497   AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 498
 499   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 500   geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
 501   geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 502
 503   //Get close cells index, energy and time, not in corners
 504
 505   Int_t absID1 = -1;
 506   Int_t absID2 = -1;
 507
 508   if( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
 509   if( iphi > 0 )                                absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
 510
 511   // In case of cell in eta = 0 border, depending on SM shift the cross cell index
 512
 513   Int_t absID3 = -1;
 514   Int_t absID4 = -1;
 515
 516   if     ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
 517   {
 518     absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
 519     absID4 = geom-> GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1);
 520   }
 521   else if( ieta == 0 && imod%2 )
 522   {
 523     absID3 = geom-> GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
 524     absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
 525   }
 526   else
 527   {
 528     if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
 529       absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
 530     if( ieta > 0 )
 531       absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
 532   }
 533
 534   //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
 535
 536   Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
 537   Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
 538   Bool_t   accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
 539
 540   accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
 541   accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
 542   accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
 543   accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
 544
 545   /*
 546    printf("Cell absID %d \n",absID);
 547    printf("\t  accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
 548    accept1,accept2,accept3,accept4);
 549    printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
 550    absID,absID1,absID2,absID3,absID4);
 551    printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
 552    ecell,ecell1,ecell2,ecell3,ecell4);
 553    printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
 554    tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
 555    TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
 556    */
 557
 558   if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
 559   if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
 560   if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
 561   if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
 562
 563   return ecell1+ecell2+ecell3+ecell4;
 564
 565 }
 566
 567 //_____________________________________________________________________________________________
 568 Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
 569 {
 570   // Look to cell neighbourhood and reject if it seems exotic
 571   // Do before recalibrating the cells
 572
 573   if(!fRejectExoticCells) return kFALSE;
 574
 575   Float_t  ecell  = 0;
 576   Double_t tcell  = 0;
 577   Bool_t   accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
 578
 579   if(!accept) return kTRUE; // reject this cell
 580
 581   if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
 582
 583   Float_t eCross = GetECross(absID,tcell,cells,bc);
 584
 585   if(1-eCross/ecell > fExoticCellFraction)
 586   {
 587     AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
 588                     absID,ecell,eCross,1-eCross/ecell));
 589     return kTRUE;
 590   }
 591
 592   return kFALSE;
 593 }
 594
 595 //___________________________________________________________________
 596 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
 597                                           AliVCaloCells *cells,
 598                                           const Int_t bc)
 599 {
 600   // Check if the cluster highest energy tower is exotic
 601
 602   if(!cluster)
 603   {
 604     AliInfo("Cluster pointer null!");
 605     return kFALSE;
 606   }
 607
 608   if(!fRejectExoticCluster) return kFALSE;
 609
 610   // Get highest energy tower
 611   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 612   Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
 613   Bool_t shared = kFALSE;
 614   GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
 615
 616   return IsExoticCell(absId,cells,bc);
 617
 618 }
 619
 620 //_______________________________________________________________________
 621 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
 622 {
 623   //In case of MC analysis, smear energy to match resolution/calibration in real data
 624
 625   if(!cluster)
 626   {
 627     AliInfo("Cluster pointer null!");
 628     return 0;
 629   }
 630
 631   Float_t energy    = cluster->E() ;
 632   Float_t rdmEnergy = energy ;
 633   if(fSmearClusterEnergy)
 634   {
 635     rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
 636                                     fSmearClusterParam[1] * energy +
 637                                     fSmearClusterParam[2] );
 638     AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
 639   }
 640
 641   return rdmEnergy;
 642 }
 643
 644 //____________________________________________________________________________
 645 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
 646 {
 647   // Correct cluster energy from non linearity functions
 648
 649   if(!cluster)
 650   {
 651     AliInfo("Cluster pointer null!");
 652     return 0;
 653   }
 654
 655   Float_t energy = cluster->E();
 656
 657   if(energy < 0.05)
 658   {
 659     // Clusters with less than 50 MeV or negative are not possible
 660     AliInfo(Form("Too Low Cluster energy!, E = %f < 0.05 GeV",energy));
 661     return 0;
 662   }
 663
 664   switch (fNonLinearityFunction)
 665   {
 666
 667     case kPi0MC:
 668     {
 669       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
 670       //fNonLinearityParams[0] = 1.014;
 671       //fNonLinearityParams[1] =-0.03329;
 672       //fNonLinearityParams[2] =-0.3853;
 673       //fNonLinearityParams[3] = 0.5423;
 674       //fNonLinearityParams[4] =-0.4335;
 675        energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 676                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 677                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 678       break;
 679     }
 680
 681     case kPi0MCv2:
 682     {
 683       //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
 684       //fNonLinearityParams[0] = 3.11111e-02;
 685       //fNonLinearityParams[1] =-5.71666e-02;
 686       //fNonLinearityParams[2] = 5.67995e-01;
 687
 688       energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
 689       break;
 690     }
 691
 692     case kPi0MCv3:
 693     {
 694       //Same as beam test corrected, change parameters
 695       //fNonLinearityParams[0] =  9.81039e-01
 696       //fNonLinearityParams[1] =  1.13508e-01;
 697       //fNonLinearityParams[2] =  1.00173e+00;
 698       //fNonLinearityParams[3] =  9.67998e-02;
 699       //fNonLinearityParams[4] =  2.19381e+02;
 700       //fNonLinearityParams[5] =  6.31604e+01;
 701       //fNonLinearityParams[6] =  1;
 702       energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
 703
 704       break;
 705     }
 706
 707
 708     case kPi0GammaGamma:
 709     {
 710       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 711       //fNonLinearityParams[0] = 1.04;
 712       //fNonLinearityParams[1] = -0.1445;
 713       //fNonLinearityParams[2] = 1.046;
 714       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 715       break;
 716     }
 717
 718     case kPi0GammaConversion:
 719     {
 720       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 721       //fNonLinearityParams[0] = 0.139393/0.1349766;
 722       //fNonLinearityParams[1] = 0.0566186;
 723       //fNonLinearityParams[2] = 0.982133;
 724       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 725
 726       break;
 727     }
 728
 729     case kBeamTest:
 730     {
 731       //From beam test, Alexei's results, for different ZS thresholds
 732       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
 733       //fNonLinearityParams[0] = 1.007;      1.003;      1.002
 734       //fNonLinearityParams[1] = 0.894;      0.719;      0.797
 735       //fNonLinearityParams[2] = 0.246;      0.334;      0.358
 736       //Rescale the param[0] with 1.03
 737       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
 738
 739       break;
 740     }
 741
 742     case kBeamTestCorrected:
 743     {
 744       //From beam test, corrected for material between beam and EMCAL
 745       //fNonLinearityParams[0] =  0.99078
 746       //fNonLinearityParams[1] =  0.161499;
 747       //fNonLinearityParams[2] =  0.655166;
 748       //fNonLinearityParams[3] =  0.134101;
 749       //fNonLinearityParams[4] =  163.282;
 750       //fNonLinearityParams[5] =  23.6904;
 751       //fNonLinearityParams[6] =  0.978;
 752         energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
 753
 754       break;
 755     }
 756
 757     case kNoCorrection:
 758       AliDebug(2,"No correction on the energy\n");
 759       break;
 760
 761   }
 762
 763   return energy;
 764 }
 765
 766 //__________________________________________________
 767 void AliEMCALRecoUtils::InitNonLinearityParam()
 768 {
 769   //Initialising Non Linearity Parameters
 770
 771   if(fNonLinearityFunction == kPi0MC)
 772   {
 773     fNonLinearityParams[0] = 1.014;
 774     fNonLinearityParams[1] = -0.03329;
 775     fNonLinearityParams[2] = -0.3853;
 776     fNonLinearityParams[3] = 0.5423;
 777     fNonLinearityParams[4] = -0.4335;
 778   }
 779
 780   if(fNonLinearityFunction == kPi0MCv2)
 781   {
 782     fNonLinearityParams[0] = 3.11111e-02;
 783     fNonLinearityParams[1] =-5.71666e-02;
 784     fNonLinearityParams[2] = 5.67995e-01;
 785   }
 786
 787   if(fNonLinearityFunction == kPi0MCv3)
 788   {
 789     fNonLinearityParams[0] =  9.81039e-01;
 790     fNonLinearityParams[1] =  1.13508e-01;
 791     fNonLinearityParams[2] =  1.00173e+00;
 792     fNonLinearityParams[3] =  9.67998e-02;
 793     fNonLinearityParams[4] =  2.19381e+02;
 794     fNonLinearityParams[5] =  6.31604e+01;
 795     fNonLinearityParams[6] =  1;
 796   }
 797
 798   if(fNonLinearityFunction == kPi0GammaGamma)
 799   {
 800     fNonLinearityParams[0] = 1.04;
 801     fNonLinearityParams[1] = -0.1445;
 802     fNonLinearityParams[2] = 1.046;
 803   }
 804
 805   if(fNonLinearityFunction == kPi0GammaConversion)
 806   {
 807     fNonLinearityParams[0] = 0.139393;
 808     fNonLinearityParams[1] = 0.0566186;
 809     fNonLinearityParams[2] = 0.982133;
 810   }
 811
 812   if(fNonLinearityFunction == kBeamTest)
 813   {
 814     if(fNonLinearThreshold == 30)
 815     {
 816       fNonLinearityParams[0] = 1.007;
 817       fNonLinearityParams[1] = 0.894;
 818       fNonLinearityParams[2] = 0.246;
 819     }
 820     if(fNonLinearThreshold == 45)
 821     {
 822       fNonLinearityParams[0] = 1.003;
 823       fNonLinearityParams[1] = 0.719;
 824       fNonLinearityParams[2] = 0.334;
 825     }
 826     if(fNonLinearThreshold == 75)
 827     {
 828       fNonLinearityParams[0] = 1.002;
 829       fNonLinearityParams[1] = 0.797;
 830       fNonLinearityParams[2] = 0.358;
 831     }
 832   }
 833
 834   if(fNonLinearityFunction == kBeamTestCorrected)
 835   {
 836     fNonLinearityParams[0] =  0.99078;
 837     fNonLinearityParams[1] =  0.161499;
 838     fNonLinearityParams[2] =  0.655166;
 839     fNonLinearityParams[3] =  0.134101;
 840     fNonLinearityParams[4] =  163.282;
 841     fNonLinearityParams[5] =  23.6904;
 842     fNonLinearityParams[6] =  0.978;
 843   }
 844 }
 845
 846 //_________________________________________________________
 847 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy,
 848                                      const Int_t iParticle,
 849                                      const Int_t iSM) const
 850 {
 851   //Calculate shower depth for a given cluster energy and particle type
 852
 853   // parameters
 854   Float_t x0    = 1.31;
 855   Float_t ecr   = 8;
 856   Float_t depth = 0;
 857   Float_t arg   = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV
 858
 859   switch ( iParticle )
 860   {
 861     case kPhoton:
 862       if (arg < 1)
 863         depth = 0;
 864       else
 865         depth = x0 * (TMath::Log(arg) + 0.5);
 866       break;
 867
 868     case kElectron:
 869       if (arg < 1)
 870         depth = 0;
 871       else
 872         depth = x0 * (TMath::Log(arg) - 0.5);
 873       break;
 874
 875     case kHadron:
 876       // hadron
 877       // boxes anc. here
 878       if(gGeoManager)
 879       {
 880         gGeoManager->cd("ALIC_1/XEN1_1");
 881         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
 882         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
 883         if(geoSM)
 884         {
 885           TGeoVolume      *geoSMVol   = geoSM->GetVolume();
 886           TGeoShape       *geoSMShape = geoSMVol->GetShape();
 887           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
 888           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
 889           else AliFatal("Null GEANT box");
 890         }
 891         else AliFatal("NULL  GEANT node matrix");
 892       }
 893       else
 894       {//electron
 895         if (arg < 1)
 896           depth = 0;
 897         else
 898           depth = x0 * (TMath::Log(arg) - 0.5);
 899       }
 900
 901       break;
 902
 903     default://photon
 904       if (arg < 1)
 905         depth = 0;
 906       else
 907         depth = x0 * (TMath::Log(arg) + 0.5);
 908   }
 909
 910   return depth;
 911 }
 912
 913 //____________________________________________________________________
 914 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
 915                                          AliVCaloCells* cells,
 916                                          const AliVCluster* clu,
 917                                          Int_t  & absId,
 918                                          Int_t  & iSupMod,
 919                                          Int_t  & ieta,
 920                                          Int_t  & iphi,
 921                                          Bool_t & shared)
 922 {
 923   //For a given CaloCluster gets the absId of the cell
 924   //with maximum energy deposit.
 925
 926   Double_t eMax        = -1.;
 927   Double_t eCell       = -1.;
 928   Float_t  fraction    = 1.;
 929   Float_t  recalFactor = 1.;
 930   Int_t    cellAbsId   = -1 ;
 931
 932   Int_t iTower  = -1;
 933   Int_t iIphi   = -1;
 934   Int_t iIeta   = -1;
 935   Int_t iSupMod0= -1;
 936
 937   if(!clu)
 938   {
 939     AliInfo("Cluster pointer null!");
 940     absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
 941     return;
 942   }
 943
 944   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
 945   {
 946     cellAbsId = clu->GetCellAbsId(iDig);
 947     fraction  = clu->GetCellAmplitudeFraction(iDig);
 948     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
 949     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 950     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 951     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 952     if     (iDig==0)
 953     {
 954       iSupMod0=iSupMod;
 955     }
 956     else if(iSupMod0!=iSupMod)
 957     {
 958       shared = kTRUE;
 959       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
 960     }
 961     if(!fCellsRecalibrated && IsRecalibrationOn())
 962     {
 963       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 964     }
 965     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 966     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
 967     if(eCell > eMax)
 968     {
 969       eMax  = eCell;
 970       absId = cellAbsId;
 971       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 972     }
 973   }// cell loop
 974
 975   //Get from the absid the supermodule, tower and eta/phi numbers
 976   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 977   //Gives SuperModule and Tower numbers
 978   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 979                                          iIphi, iIeta,iphi,ieta);
 980   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
 981   //printf("Max end---\n");
 982 }
 983
 984 //______________________________________
 985 void AliEMCALRecoUtils::InitParameters()
 986 {
 987   // Initialize data members with default values
 988
 989   fParticleType = kPhoton;
 990   fPosAlgo      = kUnchanged;
 991   fW0           = 4.5;
 992
 993   fNonLinearityFunction = kNoCorrection;
 994   fNonLinearThreshold   = 30;
 995
 996   fExoticCellFraction     = 0.97;
 997   fExoticCellDiffTime     = 1e6;
 998   fExoticCellMinAmplitude = 0.5;
 999
1000   fAODFilterMask   = 32;
1001   fAODHybridTracks = kFALSE;
1002
1003   fCutEtaPhiSum      = kTRUE;
1004   fCutEtaPhiSeparate = kFALSE;
1005
1006   fCutR   = 0.05;
1007   fCutEta = 0.025;
1008   fCutPhi = 0.05;
1009
1010   fClusterWindow = 100;
1011   fMass          = 0.139;
1012
1013   fStepSurface   = 20.;
1014   fStepCluster   = 5.;
1015   fTrackCutsType = kLooseCut;
1016
1017   fCutMinTrackPt     = 0;
1018   fCutMinNClusterTPC = -1;
1019   fCutMinNClusterITS = -1;
1020
1021   fCutMaxChi2PerClusterTPC  = 1e10;
1022   fCutMaxChi2PerClusterITS  = 1e10;
1023
1024   fCutRequireTPCRefit     = kFALSE;
1025   fCutRequireITSRefit     = kFALSE;
1026   fCutAcceptKinkDaughters = kFALSE;
1027
1028   fCutMaxDCAToVertexXY = 1e10;
1029   fCutMaxDCAToVertexZ  = 1e10;
1030   fCutDCAToVertex2D    = kFALSE;
1031
1032   fCutRequireITSStandAlone = kFALSE; //MARCEL
1033   fCutRequireITSpureSA     = kFALSE; //Marcel
1034
1035   //Misalignment matrices
1036   for(Int_t i = 0; i < 15 ; i++)
1037   {
1038     fMisalTransShift[i] = 0.;
1039     fMisalRotShift[i]   = 0.;
1040   }
1041
1042   //Non linearity
1043   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = 0.;
1044
1045   //For kBeamTestCorrected case, but default is no correction
1046   fNonLinearityParams[0] =  0.99078;
1047   fNonLinearityParams[1] =  0.161499;
1048   fNonLinearityParams[2] =  0.655166;
1049   fNonLinearityParams[3] =  0.134101;
1050   fNonLinearityParams[4] =  163.282;
1051   fNonLinearityParams[5] =  23.6904;
1052   fNonLinearityParams[6] =  0.978;
1053
1054   //For kPi0GammaGamma case
1055   //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
1056   //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
1057   //fNonLinearityParams[2] = 1.046;
1058
1059   //Cluster energy smearing
1060   fSmearClusterEnergy   = kFALSE;
1061   fSmearClusterParam[0] = 0.07; // * sqrt E term
1062   fSmearClusterParam[1] = 0.00; // * E term
1063   fSmearClusterParam[2] = 0.00; // constant
1064 }
1065
1066 //_____________________________________________________
1067 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
1068 {
1069   //Init EMCAL recalibration factors
1070   AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1071   //In order to avoid rewriting the same histograms
1072   Bool_t oldStatus = TH1::AddDirectoryStatus();
1073   TH1::AddDirectory(kFALSE);
1074
1075   fEMCALRecalibrationFactors = new TObjArray(12);
1076   for (int i = 0; i < 12; i++)
1077     fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1078                                              Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
1079   //Init the histograms with 1
1080   for (Int_t sm = 0; sm < 12; sm++)
1081   {
1082     for (Int_t i = 0; i < 48; i++)
1083     {
1084       for (Int_t j = 0; j < 24; j++)
1085       {
1086         SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1087       }
1088     }
1089   }
1090
1091   fEMCALRecalibrationFactors->SetOwner(kTRUE);
1092   fEMCALRecalibrationFactors->Compress();
1093
1094   //In order to avoid rewriting the same histograms
1095   TH1::AddDirectory(oldStatus);
1096 }
1097
1098 //_________________________________________________________
1099 void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1100 {
1101   //Init EMCAL recalibration factors
1102   AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1103   //In order to avoid rewriting the same histograms
1104   Bool_t oldStatus = TH1::AddDirectoryStatus();
1105   TH1::AddDirectory(kFALSE);
1106
1107   fEMCALTimeRecalibrationFactors = new TObjArray(4);
1108   for (int i = 0; i < 4; i++)
1109     fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1110                                                  Form("hAllTimeAvBC%d",i),
1111                                                  48*24*12,0.,48*24*12)          );
1112   //Init the histograms with 1
1113   for (Int_t bc = 0; bc < 4; bc++)
1114   {
1115     for (Int_t i = 0; i < 48*24*12; i++)
1116       SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
1117   }
1118
1119   fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1120   fEMCALTimeRecalibrationFactors->Compress();
1121
1122   //In order to avoid rewriting the same histograms
1123   TH1::AddDirectory(oldStatus);
1124 }
1125
1126 //____________________________________________________
1127 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1128 {
1129   //Init EMCAL bad channels map
1130   AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1131   //In order to avoid rewriting the same histograms
1132   Bool_t oldStatus = TH1::AddDirectoryStatus();
1133   TH1::AddDirectory(kFALSE);
1134
1135   fEMCALBadChannelMap = new TObjArray(12);
1136   //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1137   for (int i = 0; i < 12; i++)
1138   {
1139     fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1140   }
1141
1142   fEMCALBadChannelMap->SetOwner(kTRUE);
1143   fEMCALBadChannelMap->Compress();
1144
1145   //In order to avoid rewriting the same histograms
1146   TH1::AddDirectory(oldStatus);
1147 }
1148
1149 //____________________________________________________________________________
1150 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1151                                                  AliVCluster * cluster,
1152                                                  AliVCaloCells * cells,
1153                                                  const Int_t bc)
1154 {
1155   // Recalibrate the cluster energy and Time, considering the recalibration map
1156   // and the energy of the cells and time that compose the cluster.
1157   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1158
1159   if(!cluster)
1160   {
1161     AliInfo("Cluster pointer null!");
1162     return;
1163   }
1164
1165   //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1166   UShort_t * index    = cluster->GetCellsAbsId() ;
1167   Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1168   Int_t ncells = cluster->GetNCells();
1169
1170   //Initialize some used variables
1171   Float_t energy = 0;
1172   Int_t   absId  =-1;
1173   Int_t   icol   =-1, irow =-1, imod=1;
1174   Float_t factor = 1, frac = 0;
1175   Int_t   absIdMax = -1;
1176   Float_t emax     = 0;
1177
1178   //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1179   for(Int_t icell = 0; icell < ncells; icell++)
1180   {
1181     absId = index[icell];
1182     frac =  fraction[icell];
1183     if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1184
1185     if(!fCellsRecalibrated && IsRecalibrationOn())
1186     {
1187       // Energy
1188       Int_t iTower = -1, iIphi = -1, iIeta = -1;
1189       geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1190       if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
1191       geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1192       factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1193
1194       AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1195                       imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1196
1197     }
1198
1199     energy += cells->GetCellAmplitude(absId)*factor*frac;
1200
1201     if(emax < cells->GetCellAmplitude(absId)*factor*frac)
1202     {
1203       emax     = cells->GetCellAmplitude(absId)*factor*frac;
1204       absIdMax = absId;
1205     }
1206   }
1207
1208   AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1209
1210   cluster->SetE(energy);
1211
1212   // Recalculate time of cluster
1213   Double_t timeorg = cluster->GetTOF();
1214
1215   Double_t time = cells->GetCellTime(absIdMax);
1216   if(!fCellsRecalibrated && IsTimeRecalibrationOn())
1217     RecalibrateCellTime(absIdMax,bc,time);
1218
1219   cluster->SetTOF(time);
1220
1221   AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1222 }
1223
1224 //_____________________________________________________________
1225 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells,
1226                                          const Int_t bc)
1227 {
1228   // Recalibrate the cells time and energy, considering the recalibration map and the energy
1229   // of the cells that compose the cluster.
1230   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1231
1232   if(!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn()) return;
1233
1234   if(!cells)
1235   {
1236     AliInfo("Cells pointer null!");
1237     return;
1238   }
1239
1240   Short_t  absId  =-1;
1241   Bool_t   accept = kFALSE;
1242   Float_t  ecell  = 0;
1243   Double_t tcell  = 0;
1244   Double_t ecellin = 0;
1245   Double_t tcellin = 0;
1246   Int_t  mclabel = -1;
1247   Double_t efrac = 0;
1248
1249   Int_t nEMcell  = cells->GetNumberOfCells() ;
1250   for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1251   {
1252     cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1253
1254     accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
1255     if(!accept)
1256     {
1257       ecell = 0;
1258       tcell = -1;
1259     }
1260
1261     //Set new values
1262     cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1263   }
1264
1265   fCellsRecalibrated = kTRUE;
1266 }
1267
1268 //_______________________________________________________________________________________________________
1269 void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime) const
1270 {
1271   // Recalibrate time of cell with absID  considering the recalibration map
1272   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1273
1274   if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0)
1275   {
1276     celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;    ;
1277   }
1278 }
1279
1280 //______________________________________________________________________________
1281 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1282                                                    AliVCaloCells* cells,
1283                                                    AliVCluster* clu)
1284 {
1285   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1286
1287   if(!clu)
1288   {
1289     AliInfo("Cluster pointer null!");
1290     return;
1291   }
1292
1293   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1294   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1295   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1296 }
1297
1298 //_____________________________________________________________________________________________
1299 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom,
1300                                                                   AliVCaloCells* cells,
1301                                                                   AliVCluster* clu)
1302 {
1303   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1304   // The algorithm is a copy of what is done in AliEMCALRecPoint
1305
1306   Double_t eCell       = 0.;
1307   Float_t  fraction    = 1.;
1308   Float_t  recalFactor = 1.;
1309
1310   Int_t    absId   = -1;
1311   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
1312   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
1313   Float_t  weight = 0.,  totalWeight=0.;
1314   Float_t  newPos[3] = {0,0,0};
1315   Double_t pLocal[3], pGlobal[3];
1316   Bool_t shared = kFALSE;
1317
1318   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
1319   if (clEnergy <= 0)
1320     return;
1321   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
1322   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1323
1324   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1325
1326   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1327   {
1328     absId = clu->GetCellAbsId(iDig);
1329     fraction  = clu->GetCellAmplitudeFraction(iDig);
1330     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1331
1332     if (!fCellsRecalibrated)
1333     {
1334       geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1335       geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1336
1337       if(IsRecalibrationOn())
1338       {
1339         recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1340       }
1341     }
1342
1343     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1344
1345     weight = GetCellWeight(eCell,clEnergy);
1346     totalWeight += weight;
1347
1348     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1349     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1350     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1351     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1352
1353     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1354   }// cell loop
1355
1356   if(totalWeight>0)
1357   {
1358     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
1359   }
1360
1361   //Float_t pos[]={0,0,0};
1362   //clu->GetPosition(pos);
1363   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1364   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1365
1366   if(iSupModMax > 1) //sector 1
1367   {
1368     newPos[0] +=fMisalTransShift[3];//-=3.093;
1369     newPos[1] +=fMisalTransShift[4];//+=6.82;
1370     newPos[2] +=fMisalTransShift[5];//+=1.635;
1371     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1372   } else //sector 0
1373   {
1374     newPos[0] +=fMisalTransShift[0];//+=1.134;
1375     newPos[1] +=fMisalTransShift[1];//+=8.2;
1376     newPos[2] +=fMisalTransShift[2];//+=1.197;
1377     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1378   }
1379   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1380
1381   clu->SetPosition(newPos);
1382 }
1383
1384 //____________________________________________________________________________________________
1385 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom,
1386                                                                  AliVCaloCells* cells,
1387                                                                  AliVCluster* clu)
1388 {
1389   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1390   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1391
1392   Double_t eCell       = 1.;
1393   Float_t  fraction    = 1.;
1394   Float_t  recalFactor = 1.;
1395
1396   Int_t absId   = -1;
1397   Int_t iTower  = -1;
1398   Int_t iIphi   = -1, iIeta   = -1;
1399   Int_t iSupMod = -1, iSupModMax = -1;
1400   Int_t iphi = -1, ieta =-1;
1401   Bool_t shared = kFALSE;
1402
1403   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
1404
1405   if (clEnergy <= 0)
1406     return;
1407   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
1408   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1409
1410   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
1411   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1412   Int_t startingSM = -1;
1413
1414   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1415   {
1416     absId = clu->GetCellAbsId(iDig);
1417     fraction  = clu->GetCellAmplitudeFraction(iDig);
1418     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1419
1420     if     (iDig==0)  startingSM = iSupMod;
1421     else if(iSupMod != startingSM) areInSameSM = kFALSE;
1422
1423     eCell  = cells->GetCellAmplitude(absId);
1424
1425     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1426     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1427
1428     if (!fCellsRecalibrated)
1429     {
1430       if(IsRecalibrationOn())
1431       {
1432         recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1433       }
1434     }
1435
1436     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1437
1438     weight = GetCellWeight(eCell,clEnergy);
1439     if(weight < 0) weight = 0;
1440     totalWeight += weight;
1441     weightedCol += ieta*weight;
1442     weightedRow += iphi*weight;
1443
1444     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1445   }// cell loop
1446
1447   Float_t xyzNew[]={0.,0.,0.};
1448   if(areInSameSM == kTRUE)
1449   {
1450     //printf("In Same SM\n");
1451     weightedCol = weightedCol/totalWeight;
1452     weightedRow = weightedRow/totalWeight;
1453     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1454   }
1455   else
1456   {
1457     //printf("In Different SM\n");
1458     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1459   }
1460
1461   clu->SetPosition(xyzNew);
1462 }
1463
1464 //___________________________________________________________________________________________
1465 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
1466                                                                AliVCaloCells* cells,
1467                                                                AliVCluster * cluster)
1468 {
1469   //re-evaluate distance to bad channel with updated bad map
1470
1471   if(!fRecalDistToBadChannels) return;
1472
1473   if(!cluster)
1474   {
1475     AliInfo("Cluster pointer null!");
1476     return;
1477   }
1478
1479   //Get channels map of the supermodule where the cluster is.
1480   Int_t absIdMax  = -1, iSupMod =-1, icolM = -1, irowM = -1;
1481   Bool_t shared = kFALSE;
1482   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
1483   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1484
1485   Int_t dRrow, dRcol;
1486   Float_t  minDist = 10000.;
1487   Float_t  dist    = 0.;
1488
1489   //Loop on tower status map
1490   for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1491   {
1492     for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1493     {
1494       //Check if tower is bad.
1495       if(hMap->GetBinContent(icol,irow)==0) continue;
1496       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
1497       //       iSupMod,icol, irow, icolM,irowM);
1498
1499       dRrow=TMath::Abs(irowM-irow);
1500       dRcol=TMath::Abs(icolM-icol);
1501       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1502       if(dist < minDist)
1503       {
1504         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1505         minDist = dist;
1506       }
1507     }
1508   }
1509
1510   //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
1511   if (shared)
1512   {
1513     TH2D* hMap2 = 0;
1514     Int_t iSupMod2 = -1;
1515
1516     //The only possible combinations are (0,1), (2,3) ... (8,9)
1517     if(iSupMod%2) iSupMod2 = iSupMod-1;
1518     else          iSupMod2 = iSupMod+1;
1519     hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
1520
1521     //Loop on tower status map of second super module
1522     for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1523     {
1524       for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1525       {
1526         //Check if tower is bad.
1527         if(hMap2->GetBinContent(icol,irow)==0) continue;
1528         //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",
1529         //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
1530         dRrow=TMath::Abs(irow-irowM);
1531
1532         if(iSupMod%2)
1533         {
1534           dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1535         } else
1536         {
1537           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1538         }
1539
1540         dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1541         if(dist < minDist) minDist = dist;
1542       }
1543     }
1544   }// shared cluster in 2 SuperModules
1545
1546   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1547   cluster->SetDistanceToBadChannel(minDist);
1548 }
1549
1550 //__________________________________________________________________
1551 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1552 {
1553   //re-evaluate identification parameters with bayesian
1554
1555   if(!cluster)
1556   {
1557     AliInfo("Cluster pointer null!");
1558     return;
1559   }
1560
1561   if ( cluster->GetM02() != 0)
1562     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1563
1564   Float_t pidlist[AliPID::kSPECIESCN+1];
1565   for(Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1566
1567   cluster->SetPID(pidlist);
1568 }
1569
1570 //___________________________________________________________________________________________________________________
1571 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1572                                                                 AliVCaloCells* cells,
1573                                                                 AliVCluster * cluster,
1574                                                                 Float_t & l0,   Float_t & l1,
1575                                                                 Float_t & disp, Float_t & dEta, Float_t & dPhi,
1576                                                                 Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
1577 {
1578   // Calculates new center of gravity in the local EMCAL-module coordinates
1579   // and tranfers into global ALICE coordinates
1580   // Calculates Dispersion and main axis
1581
1582   if(!cluster)
1583   {
1584     AliInfo("Cluster pointer null!");
1585     return;
1586   }
1587
1588   Double_t eCell       = 0.;
1589   Float_t  fraction    = 1.;
1590   Float_t  recalFactor = 1.;
1591
1592   Int_t    iSupMod = -1;
1593   Int_t    iTower  = -1;
1594   Int_t    iIphi   = -1;
1595   Int_t    iIeta   = -1;
1596   Int_t    iphi    = -1;
1597   Int_t    ieta    = -1;
1598   Double_t etai    = -1.;
1599   Double_t phii    = -1.;
1600
1601   Int_t    nstat   = 0 ;
1602   Float_t  wtot    = 0.;
1603   Double_t w       = 0.;
1604   Double_t etaMean = 0.;
1605   Double_t phiMean = 0.;
1606
1607   //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added
1608   // and to check if the cluster is between 2 SM in eta
1609   Int_t   iSM0   = -1;
1610   Bool_t  shared = kFALSE;
1611   Float_t energy = 0;
1612
1613   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1614   {
1615     //Get from the absid the supermodule, tower and eta/phi numbers
1616     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1617     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1618
1619     //Check if there are cells of different SM
1620     if     (iDigit == 0   ) iSM0 = iSupMod;
1621     else if(iSupMod!= iSM0) shared = kTRUE;
1622
1623     //Get the cell energy, if recalibration is on, apply factors
1624     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1625     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1626
1627     if(IsRecalibrationOn())
1628     {
1629       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1630     }
1631
1632     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1633
1634     energy += eCell;
1635
1636   }//cell loop
1637
1638   //Loop on cells
1639   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1640   {
1641     //Get from the absid the supermodule, tower and eta/phi numbers
1642     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1643     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1644
1645     //Get the cell energy, if recalibration is on, apply factors
1646     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1647     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1648
1649     if (!fCellsRecalibrated)
1650     {
1651       if(IsRecalibrationOn())
1652       {
1653         recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1654       }
1655     }
1656
1657     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1658
1659     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1660     // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
1661     if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
1662
1663     if(cluster->E() > 0 && eCell > 0)
1664     {
1665       w  = GetCellWeight(eCell,cluster->E());
1666
1667       etai=(Double_t)ieta;
1668       phii=(Double_t)iphi;
1669
1670       if(w > 0.0)
1671       {
1672         wtot += w ;
1673         nstat++;
1674         //Shower shape
1675         sEta     += w * etai * etai ;
1676         etaMean  += w * etai ;
1677         sPhi     += w * phii * phii ;
1678         phiMean  += w * phii ;
1679         sEtaPhi  += w * etai * phii ;
1680       }
1681     }
1682     else
1683       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1684   }//cell loop
1685
1686   //Normalize to the weight
1687   if (wtot > 0)
1688   {
1689     etaMean /= wtot ;
1690     phiMean /= wtot ;
1691   }
1692   else
1693     AliError(Form("Wrong weight %f\n", wtot));
1694
1695   //Calculate dispersion
1696   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1697   {
1698     //Get from the absid the supermodule, tower and eta/phi numbers
1699     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1700     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1701
1702     //Get the cell energy, if recalibration is on, apply factors
1703     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1704     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1705     if (IsRecalibrationOn())
1706     {
1707       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1708     }
1709     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1710
1711     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1712     // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
1713     if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
1714
1715     if(cluster->E() > 0 && eCell > 0)
1716     {
1717       w  = GetCellWeight(eCell,cluster->E());
1718
1719       etai=(Double_t)ieta;
1720       phii=(Double_t)iphi;
1721       if(w > 0.0)
1722       {
1723         disp +=  w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
1724         dEta +=  w * (etai-etaMean)*(etai-etaMean) ;
1725         dPhi +=  w * (phii-phiMean)*(phii-phiMean) ;
1726       }
1727     }
1728     else
1729       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1730   }// cell loop
1731
1732   //Normalize to the weigth and set shower shape parameters
1733   if (wtot > 0 && nstat > 1)
1734   {
1735     disp    /= wtot ;
1736     dEta    /= wtot ;
1737     dPhi    /= wtot ;
1738     sEta    /= wtot ;
1739     sPhi    /= wtot ;
1740     sEtaPhi /= wtot ;
1741
1742     sEta    -= etaMean * etaMean ;
1743     sPhi    -= phiMean * phiMean ;
1744     sEtaPhi -= etaMean * phiMean ;
1745
1746     l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1747     l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1748   }
1749   else
1750   {
1751     l0   = 0. ;
1752     l1   = 0. ;
1753     dEta = 0. ; dPhi = 0. ; disp    = 0. ;
1754     sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
1755   }
1756
1757 }
1758
1759 //____________________________________________________________________________________________
1760 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1761                                                                 AliVCaloCells* cells,
1762                                                                 AliVCluster * cluster)
1763 {
1764   // Calculates new center of gravity in the local EMCAL-module coordinates
1765   // and tranfers into global ALICE coordinates
1766   // Calculates Dispersion and main axis and puts them into the cluster
1767
1768   Float_t l0   = 0., l1   = 0.;
1769   Float_t disp = 0., dEta = 0., dPhi    = 0.;
1770   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1771
1772   AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
1773                                                              dEta, dPhi, sEta, sPhi, sEtaPhi);
1774
1775   cluster->SetM02(l0);
1776   cluster->SetM20(l1);
1777   if(disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
1778
1779 }
1780
1781 //____________________________________________________________________________
1782 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
1783                                     TObjArray * clusterArr,
1784                                     const AliEMCALGeometry *geom)
1785 {
1786   //This function should be called before the cluster loop
1787   //Before call this function, please recalculate the cluster positions
1788   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1789   //Store matched cluster indexes and residuals
1790
1791   fMatchedTrackIndex  ->Reset();
1792   fMatchedClusterIndex->Reset();
1793   fResidualPhi->Reset();
1794   fResidualEta->Reset();
1795
1796   fMatchedTrackIndex  ->Set(1000);
1797   fMatchedClusterIndex->Set(1000);
1798   fResidualPhi->Set(1000);
1799   fResidualEta->Set(1000);
1800
1801   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1802   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1803
1804   // init the magnetic field if not already on
1805   if(!TGeoGlobalMagField::Instance()->GetField())
1806   {
1807     AliInfo("Init the magnetic field\n");
1808     if     (esdevent)
1809     {
1810       esdevent->InitMagneticField();
1811     }
1812     else if(aodevent)
1813     {
1814       Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
1815       Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
1816       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
1817       TGeoGlobalMagField::Instance()->SetField(field);
1818     }
1819     else
1820     {
1821       AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1822     }
1823
1824   } // Init mag field
1825
1826   if (esdevent) {
1827     UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
1828     UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
1829     Bool_t desc1 = (mask1 >> 3) & 0x1;
1830     Bool_t desc2 = (mask2 >> 3) & 0x1;
1831     if (desc1==0 || desc2==0) {
1832 //       AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
1833 //       mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
1834 //       mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
1835       fITSTrackSA=kTRUE;
1836     }
1837   }
1838
1839   TObjArray *clusterArray = 0x0;
1840   if(!clusterArr)
1841     {
1842       clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
1843       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1844   {
1845     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1846     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1847     clusterArray->AddAt(cluster,icl);
1848   }
1849     }
1850
1851   Int_t    matched=0;
1852   Double_t cv[21];
1853   for (Int_t i=0; i<21;i++) cv[i]=0;
1854   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1855   {
1856     AliExternalTrackParam *trackParam = 0;
1857
1858     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1859     AliESDtrack *esdTrack = 0;
1860     AliAODTrack *aodTrack = 0;
1861     if(esdevent)
1862     {
1863       esdTrack = esdevent->GetTrack(itr);
1864       if(!esdTrack) continue;
1865       if(!IsAccepted(esdTrack)) continue;
1866       if(esdTrack->Pt()<fCutMinTrackPt) continue;
1867       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
1868       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
1869       if(!fITSTrackSA)
1870         trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());  // if TPC Available
1871       else
1872         trackParam =  new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
1873     }
1874
1875     //If the input event is AOD, the starting point for extrapolation is at vertex
1876     //AOD tracks are selected according to its filterbit.
1877     else if(aodevent)
1878     {
1879       aodTrack = aodevent->GetTrack(itr);
1880       if(!aodTrack) continue;
1881
1882       //Check mask if not hybrid
1883       if(!fAODHybridTracks && !aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1884
1885       //Check hybrid
1886       if( fAODHybridTracks && !aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
1887
1888       if(aodTrack->Pt()<fCutMinTrackPt) continue;
1889
1890       Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
1891       if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
1892       Double_t pos[3],mom[3];
1893       aodTrack->GetXYZ(pos);
1894       aodTrack->GetPxPyPz(mom);
1895       AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
1896
1897       trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1898     }
1899
1900     //Return if the input data is not "AOD" or "ESD"
1901     else
1902     {
1903       printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1904       if(clusterArray)
1905   {
1906     clusterArray->Clear();
1907     delete clusterArray;
1908   }
1909       return;
1910     }
1911
1912     if(!trackParam) continue;
1913
1914     //Extrapolate the track to EMCal surface
1915     AliExternalTrackParam emcalParam(*trackParam);
1916     Float_t eta, phi;
1917     if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi))
1918       {
1919   if(aodevent && trackParam) delete trackParam;
1920   if(fITSTrackSA && trackParam) delete trackParam;
1921   continue;
1922       }
1923
1924 //    if(esdevent)
1925 //      {
1926 //  esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
1927 //      }
1928
1929     if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
1930       {
1931   if(aodevent && trackParam) delete trackParam;
1932   if(fITSTrackSA && trackParam) delete trackParam;
1933   continue;
1934       }
1935
1936
1937     //Find matched clusters
1938     Int_t index = -1;
1939     Float_t dEta = -999, dPhi = -999;
1940     if(!clusterArr)
1941       {
1942   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
1943       }
1944     else
1945       {
1946   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1947       }
1948
1949     if(index>-1)
1950     {
1951       fMatchedTrackIndex   ->AddAt(itr,matched);
1952       fMatchedClusterIndex ->AddAt(index,matched);
1953       fResidualEta         ->AddAt(dEta,matched);
1954       fResidualPhi         ->AddAt(dPhi,matched);
1955       matched++;
1956     }
1957     if(aodevent && trackParam) delete trackParam;
1958     if(fITSTrackSA && trackParam) delete trackParam;
1959   }//track loop
1960
1961   if(clusterArray)
1962     {
1963       clusterArray->Clear();
1964       delete clusterArray;
1965     }
1966
1967   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1968
1969   fMatchedTrackIndex   ->Set(matched);
1970   fMatchedClusterIndex ->Set(matched);
1971   fResidualPhi         ->Set(matched);
1972   fResidualEta         ->Set(matched);
1973 }
1974
1975 //________________________________________________________________________________
1976 Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
1977                                                    const AliVEvent *event,
1978                                                    const AliEMCALGeometry *geom,
1979                                                    Float_t &dEta, Float_t &dPhi)
1980 {
1981   //
1982   // This function returns the index of matched cluster to input track
1983   // Returns -1 if no match is found
1984   Int_t index = -1;
1985   Double_t phiV = track->Phi()*TMath::RadToDeg();
1986   if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
1987   AliExternalTrackParam *trackParam = 0;
1988   if(!fITSTrackSA)
1989     trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());  // If TPC
1990   else
1991     trackParam = new AliExternalTrackParam(*track);
1992
1993   if(!trackParam) return index;
1994   AliExternalTrackParam emcalParam(*trackParam);
1995   Float_t eta, phi;
1996
1997   if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) {
1998         if(fITSTrackSA) delete trackParam;
1999         return index;
2000   }
2001   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){
2002         if(fITSTrackSA) delete trackParam;
2003         return index;
2004   }
2005
2006   TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2007
2008   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
2009   {
2010     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
2011     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2012     clusterArr->AddAt(cluster,icl);
2013   }
2014
2015   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2016   clusterArr->Clear();
2017   delete clusterArr;
2018   if(fITSTrackSA) delete trackParam;
2019
2020   return index;
2021 }
2022
2023 //_______________________________________________________________________________________________
2024 Int_t  AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
2025                                                          AliExternalTrackParam *trkParam,
2026                                                          const TObjArray * clusterArr,
2027                                                          Float_t &dEta, Float_t &dPhi)
2028 {
2029   // Find matched cluster in array
2030
2031   dEta=-999, dPhi=-999;
2032   Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2033   Int_t index = -1;
2034   Float_t tmpEta=-999, tmpPhi=-999;
2035
2036   Double_t exPos[3] = {0.,0.,0.};
2037   if(!emcalParam->GetXYZ(exPos)) return index;
2038
2039   Float_t clsPos[3] = {0.,0.,0.};
2040   for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
2041     {
2042       AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2043       if(!cluster || !cluster->IsEMCAL()) continue;
2044       cluster->GetPosition(clsPos);
2045       Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
2046       if(dR > fClusterWindow) continue;
2047
2048       AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
2049       if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
2050       if(fCutEtaPhiSum)
2051         {
2052           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2053           if(tmpR<dRMax)
2054       {
2055         dRMax=tmpR;
2056         dEtaMax=tmpEta;
2057         dPhiMax=tmpPhi;
2058         index=icl;
2059       }
2060         }
2061       else if(fCutEtaPhiSeparate)
2062         {
2063           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
2064       {
2065         dEtaMax = tmpEta;
2066         dPhiMax = tmpPhi;
2067         index=icl;
2068       }
2069         }
2070       else
2071         {
2072           printf("Error: please specify your cut criteria\n");
2073           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2074           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2075           return index;
2076         }
2077     }
2078
2079   dEta=dEtaMax;
2080   dPhi=dPhiMax;
2081
2082   return index;
2083 }
2084
2085 //------------------------------------------------------------------------------------
2086 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
2087                                                          const Double_t emcalR,
2088                                                          const Double_t mass,
2089                                                          const Double_t step,
2090                                                          Float_t &eta,
2091                                                          Float_t &phi)
2092 {
2093   //Extrapolate track to EMCAL surface
2094
2095   eta = -999, phi = -999;
2096   if(!trkParam) return kFALSE;
2097   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
2098   Double_t trkPos[3] = {0.,0.,0.};
2099   if(!trkParam->GetXYZ(trkPos)) return kFALSE;
2100   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2101   eta = trkPosVec.Eta();
2102   phi = trkPosVec.Phi();
2103   if(phi<0)
2104     phi += 2*TMath::Pi();
2105
2106   return kTRUE;
2107 }
2108
2109 //-----------------------------------------------------------------------------------
2110 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2111                                                      const Float_t *clsPos,
2112                                                      Double_t mass,
2113                                                      Double_t step,
2114                                                      Float_t &tmpEta,
2115                                                      Float_t &tmpPhi)
2116 {
2117   //
2118   //Return the residual by extrapolating a track param to a global position
2119   //
2120   tmpEta = -999;
2121   tmpPhi = -999;
2122   if(!trkParam) return kFALSE;
2123   Double_t trkPos[3] = {0.,0.,0.};
2124   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2125   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2126   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
2127   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
2128   if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
2129
2130   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2131   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2132
2133   // track cluster matching
2134   tmpPhi = clsPosVec.DeltaPhi(trkPosVec);    // tmpPhi is between -pi and pi
2135   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2136
2137   return kTRUE;
2138 }
2139
2140 //----------------------------------------------------------------------------------
2141 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
2142                                                     const AliVCluster *cluster,
2143                                                     const Double_t mass,
2144                                                     const Double_t step,
2145                                                     Float_t &tmpEta,
2146                                                     Float_t &tmpPhi)
2147 {
2148   //
2149   //Return the residual by extrapolating a track param to a cluster
2150   //
2151   tmpEta = -999;
2152   tmpPhi = -999;
2153   if(!cluster || !trkParam) return kFALSE;
2154
2155   Float_t clsPos[3] = {0.,0.,0.};
2156   cluster->GetPosition(clsPos);
2157
2158   return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2159 }
2160
2161 //---------------------------------------------------------------------------------
2162 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
2163                                                     const AliVCluster *cluster,
2164                                                     Float_t &tmpEta,
2165                                                     Float_t &tmpPhi)
2166 {
2167   //
2168   //Return the residual by extrapolating a track param to a clusterfStepCluster
2169   //
2170
2171   return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
2172 }
2173
2174 //_______________________________________________________________________
2175 void AliEMCALRecoUtils::GetMatchedResiduals(const Int_t clsIndex,
2176                                             Float_t &dEta, Float_t &dPhi)
2177 {
2178   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2179   //Get the residuals dEta and dPhi for this cluster to the closest track
2180   //Works with ESDs and AODs
2181
2182   if( FindMatchedPosForCluster(clsIndex) >= 999 )
2183   {
2184     AliDebug(2,"No matched tracks found!\n");
2185     dEta=999.;
2186     dPhi=999.;
2187     return;
2188   }
2189   dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2190   dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
2191 }
2192
2193 //______________________________________________________________________________________________
2194 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
2195 {
2196   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2197   //Get the residuals dEta and dPhi for this track to the closest cluster
2198   //Works with ESDs and AODs
2199
2200   if( FindMatchedPosForTrack(trkIndex) >= 999 )
2201   {
2202     AliDebug(2,"No matched cluster found!\n");
2203     dEta=999.;
2204     dPhi=999.;
2205     return;
2206   }
2207   dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2208   dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
2209 }
2210
2211 //__________________________________________________________
2212 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2213 {
2214   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2215   //Get the index of matched track to this cluster
2216   //Works with ESDs and AODs
2217
2218   if(IsClusterMatched(clsIndex))
2219     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2220   else
2221     return -1;
2222 }
2223
2224 //__________________________________________________________
2225 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
2226 {
2227   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2228   //Get the index of matched cluster to this track
2229   //Works with ESDs and AODs
2230
2231   if(IsTrackMatched(trkIndex))
2232     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
2233   else
2234     return -1;
2235 }
2236
2237 //______________________________________________________________
2238 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
2239 {
2240   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2241   //Returns if the cluster has a match
2242   if(FindMatchedPosForCluster(clsIndex) < 999)
2243     return kTRUE;
2244   else
2245     return kFALSE;
2246 }
2247
2248 //____________________________________________________________
2249 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
2250 {
2251   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2252   //Returns if the track has a match
2253   if(FindMatchedPosForTrack(trkIndex) < 999)
2254     return kTRUE;
2255   else
2256     return kFALSE;
2257 }
2258
2259 //______________________________________________________________________
2260 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
2261 {
2262   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2263   //Returns the position of the match in the fMatchedClusterIndex array
2264   Float_t tmpR = fCutR;
2265   UInt_t pos = 999;
2266
2267   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
2268   {
2269     if(fMatchedClusterIndex->At(i)==clsIndex)
2270     {
2271       Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
2272       if(r<tmpR)
2273       {
2274         pos=i;
2275         tmpR=r;
2276         AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2277                         fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
2278       }
2279     }
2280   }
2281   return pos;
2282 }
2283
2284 //____________________________________________________________________
2285 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2286 {
2287   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2288   //Returns the position of the match in the fMatchedTrackIndex array
2289   Float_t tmpR = fCutR;
2290   UInt_t pos = 999;
2291
2292   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
2293   {
2294     if(fMatchedTrackIndex->At(i)==trkIndex)
2295     {
2296       Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
2297       if(r<tmpR)
2298       {
2299         pos=i;
2300         tmpR=r;
2301         AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2302                         fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
2303       }
2304     }
2305   }
2306   return pos;
2307 }
2308
2309 //__________________________________________________________________________
2310 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
2311                                         const AliEMCALGeometry *geom,
2312                                         AliVCaloCells* cells,const Int_t bc)
2313 {
2314   // check if the cluster survives some quality cut
2315   //
2316   //
2317   Bool_t isGood=kTRUE;
2318
2319   if(!cluster || !cluster->IsEMCAL())              return kFALSE;
2320
2321   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
2322
2323   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
2324
2325   if(IsExoticCluster(cluster, cells,bc))           return kFALSE;
2326
2327   return isGood;
2328 }
2329
2330 //__________________________________________________________
2331 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2332 {
2333   // Given a esd track, return whether the track survive all the cuts
2334
2335   // The different quality parameter are first
2336   // retrieved from the track. then it is found out what cuts the
2337   // track did not survive and finally the cuts are imposed.
2338
2339   UInt_t status = esdTrack->GetStatus();
2340
2341   Int_t nClustersITS = esdTrack->GetITSclusters(0);
2342   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2343
2344   Float_t chi2PerClusterITS = -1;
2345   Float_t chi2PerClusterTPC = -1;
2346   if (nClustersITS!=0)
2347     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2348   if (nClustersTPC!=0)
2349     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
2350
2351
2352   //DCA cuts
2353   if(fTrackCutsType==kGlobalCut)
2354     {
2355       Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2356       //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2357       SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2358     }
2359
2360
2361   Float_t b[2];
2362   Float_t bCov[3];
2363   esdTrack->GetImpactParameters(b,bCov);
2364   if (bCov[0]<=0 || bCov[2]<=0)
2365   {
2366     AliDebug(1, "Estimated b resolution lower or equal zero!");
2367     bCov[0]=0; bCov[2]=0;
2368   }
2369
2370   Float_t dcaToVertexXY = b[0];
2371   Float_t dcaToVertexZ = b[1];
2372   Float_t dcaToVertex = -1;
2373
2374   if (fCutDCAToVertex2D)
2375     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2376   else
2377     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2378
2379   // cut the track?
2380
2381   Bool_t cuts[kNCuts];
2382   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2383
2384   // track quality cuts
2385   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2386     cuts[0]=kTRUE;
2387   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2388     cuts[1]=kTRUE;
2389   if (nClustersTPC<fCutMinNClusterTPC)
2390     cuts[2]=kTRUE;
2391   if (nClustersITS<fCutMinNClusterITS)
2392     cuts[3]=kTRUE;
2393   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
2394     cuts[4]=kTRUE;
2395   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
2396     cuts[5]=kTRUE;
2397   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2398     cuts[6]=kTRUE;
2399   if (fCutDCAToVertex2D && dcaToVertex > 1)
2400     cuts[7] = kTRUE;
2401   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2402     cuts[8] = kTRUE;
2403   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2404     cuts[9] = kTRUE;
2405
2406   if(fTrackCutsType==kGlobalCut)
2407     {
2408       //Require at least one SPD point + anything else in ITS
2409       if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
2410   cuts[10] = kTRUE;
2411     }
2412
2413       // ITS
2414   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
2415     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
2416       // TPC tracks
2417       cuts[11] = kTRUE;
2418     }else{
2419       // ITS standalone tracks
2420       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
2421         if(status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
2422       }else if(fCutRequireITSpureSA){
2423         if(!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
2424       }
2425     }
2426   }
2427
2428   Bool_t cut=kFALSE;
2429   for (Int_t i=0; i<kNCuts; i++)
2430     if (cuts[i]) { cut = kTRUE ; }
2431
2432     // cut the track
2433   if (cut)
2434     return kFALSE;
2435   else
2436     return kTRUE;
2437 }
2438
2439 //_____________________________________
2440 void AliEMCALRecoUtils::InitTrackCuts()
2441 {
2442   //Intilize the track cut criteria
2443   //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
2444   //Also you can customize the cuts using the setters
2445
2446   switch (fTrackCutsType)
2447   {
2448     case kTPCOnlyCut:
2449     {
2450       AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2451       //TPC
2452       SetMinNClustersTPC(70);
2453       SetMaxChi2PerClusterTPC(4);
2454       SetAcceptKinkDaughters(kFALSE);
2455       SetRequireTPCRefit(kFALSE);
2456
2457       //ITS
2458       SetRequireITSRefit(kFALSE);
2459       SetMaxDCAToVertexZ(3.2);
2460       SetMaxDCAToVertexXY(2.4);
2461       SetDCAToVertex2D(kTRUE);
2462
2463       break;
2464     }
2465
2466     case kGlobalCut:
2467     {
2468       AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2469       //TPC
2470       SetMinNClustersTPC(70);
2471       SetMaxChi2PerClusterTPC(4);
2472       SetAcceptKinkDaughters(kFALSE);
2473       SetRequireTPCRefit(kTRUE);
2474
2475       //ITS
2476       SetRequireITSRefit(kTRUE);
2477       SetMaxDCAToVertexZ(2);
2478       SetMaxDCAToVertexXY();
2479       SetDCAToVertex2D(kFALSE);
2480
2481       break;
2482     }
2483
2484     case kLooseCut:
2485     {
2486       AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2487       SetMinNClustersTPC(50);
2488       SetAcceptKinkDaughters(kTRUE);
2489
2490       break;
2491     }
2492
2493     case kITSStandAlone:
2494     {
2495       AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
2496       SetRequireITSRefit(kTRUE);
2497       SetRequireITSStandAlone(kTRUE);
2498       SetITSTrackSA(kTRUE);
2499       break;
2500     }
2501
2502   }
2503 }
2504
2505
2506 //________________________________________________________________________
2507 void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
2508 {
2509   // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
2510
2511   Int_t nTracks = event->GetNumberOfTracks();
2512   for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
2513   {
2514     AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
2515     if (!track)
2516     {
2517       AliWarning(Form("Could not receive track %d", iTrack));
2518       continue;
2519     }
2520
2521     Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
2522     track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
2523     /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2524     AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2525     if (esdtrack) {
2526       if(matchClusIndex != -1)
2527         esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2528       else
2529         esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2530     } else {
2531       AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
2532       if(matchClusIndex != -1)
2533         aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2534       else
2535         aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2536     }
2537
2538   }
2539   AliDebug(2,"Track matched to closest cluster");
2540 }
2541
2542 //_________________________________________________________________________
2543 void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
2544 {
2545   // Checks if EMC clusters are matched to ESD track.
2546   // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2547
2548   for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
2549   {
2550     AliVCluster *cluster = event->GetCaloCluster(iClus);
2551     if (!cluster->IsEMCAL())
2552       continue;
2553
2554     Int_t nTracks = event->GetNumberOfTracks();
2555     TArrayI arrayTrackMatched(nTracks);
2556
2557     // Get the closest track matched to the cluster
2558     Int_t nMatched = 0;
2559     Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
2560     if (matchTrackIndex != -1)
2561     {
2562       arrayTrackMatched[nMatched] = matchTrackIndex;
2563       nMatched++;
2564     }
2565
2566     // Get all other tracks matched to the cluster
2567     for(Int_t iTrk=0; iTrk<nTracks; ++iTrk)
2568     {
2569       AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
2570       if(iTrk == matchTrackIndex) continue;
2571       if(track->GetEMCALcluster() == iClus)
2572       {
2573         arrayTrackMatched[nMatched] = iTrk;
2574         ++nMatched;
2575       }
2576     }
2577
2578     //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2579
2580     arrayTrackMatched.Set(nMatched);
2581     AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2582     if (esdcluster)
2583       esdcluster->AddTracksMatched(arrayTrackMatched);
2584     else if (nMatched>0) {
2585       AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2586       if (aodcluster)
2587         aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2588     }
2589
2590     Float_t eta= -999, phi = -999;
2591     if (matchTrackIndex != -1)
2592       GetMatchedResiduals(iClus, eta, phi);
2593     cluster->SetTrackDistance(phi, eta);
2594   }
2595
2596   AliDebug(2,"Cluster matched to tracks");
2597 }
2598
2599 //___________________________________________________
2600 void AliEMCALRecoUtils::Print(const Option_t *) const
2601 {
2602   // Print Parameters
2603
2604   printf("AliEMCALRecoUtils Settings: \n");
2605   printf("Misalignment shifts\n");
2606   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,
2607                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2608                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
2609   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
2610   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
2611
2612   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
2613
2614   printf("Matching criteria: ");
2615   if(fCutEtaPhiSum)
2616     {
2617       printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
2618     }
2619   else if(fCutEtaPhiSeparate)
2620     {
2621       printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
2622     }
2623   else
2624     {
2625       printf("Error\n");
2626       printf("please specify your cut criteria\n");
2627       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2628       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2629     }
2630
2631   printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
2632   printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
2633
2634   printf("Track cuts: \n");
2635   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
2636   printf("AOD track selection mask: %d\n",fAODFilterMask);
2637   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2638   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2639   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2640   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2641   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
2642 }
2643