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