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