EMCAL/AliEMCALRecoUtils.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 // Class AliEMCALRecoUtils
  21 // Some utilities to recalculate the cluster position or energy linearity
  22 //
  23 //
  24 // Author:  Gustavo Conesa (LPSC- Grenoble)
  25 ///////////////////////////////////////////////////////////////////////////////
  26
  27 // --- standard c ---
  28
  29 // standard C++ includes
  30 //#include <Riostream.h>
  31
  32 // ROOT includes
  33 #include <TGeoManager.h>
  34 #include <TGeoMatrix.h>
  35 #include <TGeoBBox.h>
  36
  37 // STEER includes
  38 #include "AliEMCALRecoUtils.h"
  39 #include "AliEMCALGeometry.h"
  40 #include "AliVCluster.h"
  41 #include "AliVCaloCells.h"
  42 #include "AliLog.h"
  43
  44 ClassImp(AliEMCALRecoUtils)
  45
  46 //______________________________________________
  47 AliEMCALRecoUtils::AliEMCALRecoUtils():
  48   fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
  49   fPosAlgo(kUnchanged),fW0(4.),
  50   fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
  51   fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(),
  52   fNCellsFromEMCALBorder(0),fNoEMCALBorderAtEta0(kFALSE)
  53 {
  54 //
  55   // Constructor.
  56   // Initialize all constant values which have to be used
  57   // during Reco algorithm execution
  58   //
  59
  60   for(Int_t i = 0; i < 15 ; i++) {
  61       fMisalTransShift[i] = 0.;
  62       fMisalRotShift[i] = 0.;
  63   }
  64   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = 0.;
  65   //For kPi0GammaGamma case, but default is no correction
  66   fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
  67   fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
  68   fNonLinearityParams[2] = 1.046;
  69
  70 }
  71
  72 //______________________________________________________________________
  73 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
  74 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
  75   fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
  76   fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
  77   fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
  78   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0)
  79
  80 {
  81   //Copy ctor
  82
  83   for(Int_t i = 0; i < 15 ; i++) {
  84       fMisalRotShift[i] = reco.fMisalRotShift[i];
  85       fMisalTransShift[i] = reco.fMisalTransShift[i];
  86   }
  87   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
  88 }
  89
  90
  91 //______________________________________________________________________
  92 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
  93 {
  94   //Assignment operator
  95
  96   if(this == &reco)return *this;
  97   ((TNamed *)this)->operator=(reco);
  98
  99   fNonLinearityFunction  = reco.fNonLinearityFunction;
 100   fParticleType          = reco.fParticleType;
 101   fPosAlgo               = reco.fPosAlgo;
 102   fW0                    = reco.fW0;
 103   fRecalibration         = reco.fRecalibration;
 104   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
 105   fRemoveBadChannels     = reco.fRemoveBadChannels;
 106   fEMCALBadChannelMap    = reco.fEMCALBadChannelMap;
 107   fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
 108   fNoEMCALBorderAtEta0   = reco.fNoEMCALBorderAtEta0;
 109
 110   for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
 111   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
 112
 113   return *this;
 114 }
 115
 116
 117 //__________________________________________________
 118 AliEMCALRecoUtils::~AliEMCALRecoUtils()
 119 {
 120   //Destructor.
 121
 122         if(fEMCALRecalibrationFactors) {
 123                 fEMCALRecalibrationFactors->Clear();
 124                 delete  fEMCALRecalibrationFactors;
 125         }
 126
 127   if(fEMCALBadChannelMap) {
 128                 fEMCALBadChannelMap->Clear();
 129                 delete  fEMCALBadChannelMap;
 130         }
 131
 132 }
 133
 134 //_______________________________________________________________
 135 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
 136 {
 137         // Given the list of AbsId of the cluster, get the maximum cell and
 138         // check if there are fNCellsFromBorder from the calorimeter border
 139
 140   //If the distance to the border is 0 or negative just exit accept all clusters
 141         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
 142
 143   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;
 144   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi);
 145
 146   AliDebug(2,Form("AliEMCALRecoUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 147            cells->GetCellAmplitude(absIdMax), cluster->E()));
 148
 149         if(absIdMax==-1) return kFALSE;
 150
 151         //Check if the cell is close to the borders:
 152         Bool_t okrow = kFALSE;
 153         Bool_t okcol = kFALSE;
 154
 155   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 156     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
 157                   iSM,ieta,iphi));
 158   }
 159
 160   //Check rows/phi
 161   if(iSM < 10){
 162     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
 163   }
 164   else{
 165     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
 166   }
 167
 168   //Check columns/eta
 169   if(!fNoEMCALBorderAtEta0){
 170     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
 171   }
 172   else{
 173     if(iSM%2==0){
 174       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;
 175     }
 176     else {
 177       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;
 178     }
 179   }//eta 0 not checked
 180
 181   AliDebug(2,Form("AliEMCALRecoUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d",
 182            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
 183
 184         if (okcol && okrow) return kTRUE;
 185         else                return kFALSE;
 186
 187 }
 188
 189
 190 //_________________________________________________________________________________________________________
 191 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
 192         // Check that in the cluster cells, there is no bad channel of those stored
 193         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 194
 195         if(!fRemoveBadChannels)  return kFALSE;
 196         if(!fEMCALBadChannelMap) return kFALSE;
 197
 198         Int_t icol = -1;
 199         Int_t irow = -1;
 200         Int_t imod = -1;
 201         for(Int_t iCell = 0; iCell<nCells; iCell++){
 202
 203                 //Get the column and row
 204     Int_t iTower = -1, iIphi = -1, iIeta = -1;
 205     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
 206     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
 207     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 208     if(GetEMCALChannelStatus(imod, icol, irow))return kTRUE;
 209
 210         }// cell cluster loop
 211
 212         return kFALSE;
 213
 214 }
 215
 216 //__________________________________________________
 217 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
 218 // Correct cluster energy from non linearity functions
 219   Float_t energy = cluster->E();
 220
 221   switch (fNonLinearityFunction) {
 222
 223     case kPi0MC:
 224       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
 225       //Double_t par0 = 1.001;
 226       //Double_t par1 = -0.01264;
 227       //Double_t par2 = -0.03632;
 228       //Double_t par3 = 0.1798;
 229       //Double_t par4 = -0.522;
 230        energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 231                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 232                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 233       break;
 234
 235     case kPi0GammaGamma:
 236
 237       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 238       //Double_t par0 = 0.1457;
 239       //Double_t par1 = -0.02024;
 240       //Double_t par2 = 1.046;
 241       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 242       break;
 243
 244     case kPi0GammaConversion:
 245
 246       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 247       //Double_t C = 0.139393/0.1349766;
 248       //Double_t a = 0.0566186;
 249       //Double_t b = 0.982133;
 250       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 251
 252       break;
 253
 254     case kNoCorrection:
 255       AliDebug(2,"No correction on the energy\n");
 256       break;
 257
 258   }
 259
 260   return energy;
 261
 262 }
 263
 264 //__________________________________________________
 265 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
 266 {
 267   //Calculate shower depth for a given cluster energy and particle type
 268
 269   // parameters
 270   Float_t x0    = 1.23;
 271   Float_t ecr   = 8;
 272   Float_t depth = 0;
 273
 274   switch ( iParticle )
 275   {
 276     case kPhoton:
 277       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 278       break;
 279
 280     case kElectron:
 281       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
 282       break;
 283
 284     case kHadron:
 285       // hadron
 286       // boxes anc. here
 287       if(gGeoManager){
 288         gGeoManager->cd("ALIC_1/XEN1_1");
 289         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
 290         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
 291         if(geoSM){
 292           TGeoVolume      *geoSMVol   = geoSM->GetVolume();
 293           TGeoShape       *geoSMShape = geoSMVol->GetShape();
 294           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
 295           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
 296           else AliFatal("Null GEANT box");
 297         }else AliFatal("NULL  GEANT node matrix");
 298       }
 299       else{//electron
 300         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
 301       }
 302
 303       break;
 304
 305     default://photon
 306       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
 307   }
 308
 309   return depth;
 310
 311 }
 312
 313 //__________________________________________________
 314 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
 315 {
 316   //For a given CaloCluster gets the absId of the cell
 317   //with maximum energy deposit.
 318
 319   Double_t eMax        = -1.;
 320   Double_t eCell       = -1.;
 321   Float_t  fraction    = 1.;
 322   Float_t  recalFactor = 1.;
 323   Int_t    cellAbsId   = -1 ;
 324
 325   Int_t iTower  = -1;
 326   Int_t iIphi   = -1;
 327   Int_t iIeta   = -1;
 328
 329   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 330     cellAbsId = clu->GetCellAbsId(iDig);
 331     fraction  = clu->GetCellAmplitudeFraction(iDig);
 332     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 333     if(IsRecalibrationOn()) {
 334       geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 335       geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 336       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 337     }
 338     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 339
 340     if(eCell > eMax)  {
 341       eMax  = eCell;
 342       absId = cellAbsId;
 343       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 344     }
 345   }// cell loop
 346
 347   //Get from the absid the supermodule, tower and eta/phi numbers
 348   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 349   //Gives SuperModule and Tower numbers
 350   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 351                                          iIphi, iIeta,iphi,ieta);
 352
 353 }
 354
 355 //________________________________________________________________
 356 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
 357         //Init EMCAL recalibration factors
 358         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
 359         //In order to avoid rewriting the same histograms
 360         Bool_t oldStatus = TH1::AddDirectoryStatus();
 361         TH1::AddDirectory(kFALSE);
 362
 363         fEMCALRecalibrationFactors = new TObjArray(12);
 364         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));
 365         //Init the histograms with 1
 366         for (Int_t sm = 0; sm < 12; sm++) {
 367                 for (Int_t i = 0; i < 48; i++) {
 368                         for (Int_t j = 0; j < 24; j++) {
 369                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
 370                         }
 371                 }
 372         }
 373         fEMCALRecalibrationFactors->SetOwner(kTRUE);
 374         fEMCALRecalibrationFactors->Compress();
 375
 376         //In order to avoid rewriting the same histograms
 377         TH1::AddDirectory(oldStatus);
 378 }
 379
 380
 381 //________________________________________________________________
 382 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 383         //Init EMCAL bad channels map
 384         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
 385         //In order to avoid rewriting the same histograms
 386         Bool_t oldStatus = TH1::AddDirectoryStatus();
 387         TH1::AddDirectory(kFALSE);
 388
 389         fEMCALBadChannelMap = new TObjArray(12);
 390         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
 391         for (int i = 0; i < 12; i++) {
 392                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
 393         }
 394
 395         //delete hTemp;
 396
 397         fEMCALBadChannelMap->SetOwner(kTRUE);
 398         fEMCALBadChannelMap->Compress();
 399
 400         //In order to avoid rewriting the same histograms
 401         TH1::AddDirectory(oldStatus);
 402 }
 403
 404 //________________________________________________________________
 405 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
 406         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 407
 408         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 409         UShort_t * index    = cluster->GetCellsAbsId() ;
 410         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 411         Int_t ncells = cluster->GetNCells();
 412
 413         //Initialize some used variables
 414         Float_t energy = 0;
 415         Int_t absId    = -1;
 416   Int_t icol = -1, irow = -1, imod=1;
 417         Float_t factor = 1, frac = 0;
 418
 419         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 420         for(Int_t icell = 0; icell < ncells; icell++){
 421                 absId = index[icell];
 422                 frac =  fraction[icell];
 423                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
 424                 Int_t iTower = -1, iIphi = -1, iIeta = -1;
 425                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 426                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
 427                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 428                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
 429     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
 430              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
 431
 432                 energy += cells->GetCellAmplitude(absId)*factor*frac;
 433         }
 434
 435
 436                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
 437
 438         cluster->SetE(energy);
 439
 440 }
 441
 442
 443 //__________________________________________________
 444 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 445 {
 446   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 447
 448   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
 449   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
 450   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
 451
 452 }
 453
 454 //__________________________________________________
 455 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 456 {
 457   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 458   // The algorithm is a copy of what is done in AliEMCALRecPoint
 459
 460   Double_t eCell       = 0.;
 461   Float_t  fraction    = 1.;
 462   Float_t  recalFactor = 1.;
 463
 464   Int_t    absId   = -1;
 465   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
 466   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
 467   Float_t  weight = 0.,  totalWeight=0.;
 468   Float_t  newPos[3] = {0,0,0};
 469   Double_t pLocal[3], pGlobal[3];
 470
 471   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
 472
 473   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
 474   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
 475
 476   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 477     absId = clu->GetCellAbsId(iDig);
 478     fraction  = clu->GetCellAmplitudeFraction(iDig);
 479     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 480     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
 481     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 482
 483     if(IsRecalibrationOn()) {
 484       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
 485     }
 486     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 487
 488     weight = GetCellWeight(eCell,clEnergy);
 489     totalWeight += weight;
 490     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
 491     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
 492
 493     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
 494
 495   }// cell loop
 496
 497   if(totalWeight>0){
 498     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
 499   }
 500
 501   //printf("iSM %d \n",iSupMod);
 502   //Float_t pos[]={0,0,0};
 503   //clu->GetPosition(pos);
 504   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
 505
 506
 507   //printf("NewPos a: %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
 508
 509         if(iSupModMax > 1) {//sector 1
 510           newPos[0] +=fMisalTransShift[3];//-=3.093;
 511           newPos[1] +=fMisalTransShift[4];//+=6.82;
 512           newPos[2] +=fMisalTransShift[5];//+=1.635;
 513         }
 514         else {//sector 0
 515           newPos[0] +=fMisalTransShift[0];//+=1.134;
 516           newPos[1] +=fMisalTransShift[1];//+=8.2;
 517           newPos[2] +=fMisalTransShift[2];//+=1.197;
 518         }
 519
 520   clu->SetPosition(newPos);
 521
 522 }
 523
 524 //__________________________________________________
 525 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 526 {
 527   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 528   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
 529
 530   Double_t eCell       = 1.;
 531   Float_t  fraction    = 1.;
 532   Float_t  recalFactor = 1.;
 533
 534   Int_t absId   = -1;
 535   Int_t iTower  = -1;
 536   Int_t iIphi   = -1, iIeta   = -1;
 537         Int_t iSupMod = -1, iSupModMax = -1;
 538   Int_t iphi = -1, ieta =-1;
 539
 540   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
 541   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
 542   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
 543
 544   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
 545   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
 546   Int_t startingSM = -1;
 547
 548   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 549     absId = clu->GetCellAbsId(iDig);
 550     fraction  = clu->GetCellAmplitudeFraction(iDig);
 551     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 552     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 553     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 554
 555     if     (iDig==0)  startingSM = iSupMod;
 556     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 557
 558     eCell  = cells->GetCellAmplitude(absId);
 559
 560     if(IsRecalibrationOn()) {
 561       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 562     }
 563     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
 564
 565     weight = GetCellWeight(eCell,clEnergy);
 566     if(weight < 0) weight = 0;
 567     totalWeight += weight;
 568     weightedCol += ieta*weight;
 569     weightedRow += iphi*weight;
 570
 571     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
 572
 573     }// cell loop
 574
 575   Float_t xyzNew[]={0.,0.,0.};
 576   if(areInSameSM == kTRUE) {
 577     //printf("In Same SM\n");
 578     weightedCol = weightedCol/totalWeight;
 579     weightedRow = weightedRow/totalWeight;
 580     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 581   }
 582   else {
 583     //printf("In Different SM\n");
 584     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
 585   }
 586
 587   clu->SetPosition(xyzNew);
 588
 589 }
 590
 591 //__________________________________________________
 592 void AliEMCALRecoUtils::Print(const Option_t *) const
 593 {
 594   // Print Parameters
 595
 596   printf("AliEMCALRecoUtils Settings: \n");
 597   printf("Misalignment shifts\n");
 598   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,
 599                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
 600                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
 601   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
 602   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
 603
 604   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
 605
 606 }