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
  34 // STEER includes
  35 #include "AliEMCALRecoUtils.h"
  36 #include "AliEMCALGeoUtils.h"
  37 #include "AliVCluster.h"
  38 #include "AliVCaloCells.h"
  39 #include "AliLog.h"
  40
  41 ClassImp(AliEMCALRecoUtils)
  42
  43 //______________________________________________
  44 AliEMCALRecoUtils::AliEMCALRecoUtils():
  45   fNonLinearityFunction (kPi0GammaGamma)
  46 {
  47 //
  48   // Constructor.
  49   // Initialize all constant values which have to be used
  50   // during Reco algorithm execution
  51   //
  52
  53   for(Int_t i = 0; i < 15 ; i++) fMisalShift[i] = 0.;
  54   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = 0.;
  55   //By default kPi0GammaGamma case
  56   fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
  57   fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
  58   fNonLinearityParams[2] = 1.046;
  59
  60 }
  61
  62 //______________________________________________________________________
  63 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
  64 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction)
  65 {
  66   //Copy ctor
  67
  68   for(Int_t i = 0; i < 15 ; i++) fMisalShift[i]         = reco.fMisalShift[i];
  69   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
  70 }
  71
  72
  73 //______________________________________________________________________
  74 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
  75 {
  76   //Assignment operator
  77
  78   if(this == &reco)return *this;
  79   ((TNamed *)this)->operator=(reco);
  80
  81   fNonLinearityFunction = reco.fNonLinearityFunction;
  82   for(Int_t i = 0; i < 15 ; i++) fMisalShift[i]         = reco.fMisalShift[i];
  83   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
  84
  85   return *this;
  86 }
  87
  88
  89 //__________________________________________________
  90 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
  91 // Correct cluster energy from non linearity functions
  92   Float_t energy = cluster->E();
  93
  94   switch (fNonLinearityFunction) {
  95
  96     case kPi0MC:
  97       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
  98       //Double_t par0 = 1.001;
  99       //Double_t par1 = -0.01264;
 100       //Double_t par2 = -0.03632;
 101       //Double_t par3 = 0.1798;
 102       //Double_t par4 = -0.522;
 103        energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
 104                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
 105                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
 106       break;
 107
 108     case kPi0GammaGamma:
 109
 110       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
 111       //Double_t par0 = 0.1457;
 112       //Double_t par1 = -0.02024;
 113       //Double_t par2 = 1.046;
 114       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
 115       break;
 116
 117     case kPi0GammaConversion:
 118
 119       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
 120       //Double_t C = 0.139393/0.1349766;
 121       //Double_t a = 0.0566186;
 122       //Double_t b = 0.982133;
 123       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
 124
 125       break;
 126
 127     case kNoCorrection:
 128       AliDebug(2,"No correction on the energy\n");
 129       break;
 130
 131   }
 132
 133   return energy;
 134
 135 }
 136
 137 //__________________________________________________
 138 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeoUtils *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
 139 {
 140   //For a given CaloCluster gets the absId of the cell
 141   //with maximum energy deposit.
 142
 143   Double_t eMax        = -1.;
 144   Double_t eCell       = -1.;
 145   Int_t    cellAbsId   = -1 ;
 146   Int_t iTower  = -1;
 147   Int_t iIphi   = -1;
 148   Int_t iIeta   = -1;
 149
 150   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 151      cellAbsId = clu->GetCellAbsId(iDig);
 152      eCell     = cells->GetCellAmplitude(cellAbsId);
 153      if(eCell > eMax)  {
 154       eMax  = eCell;
 155       absId = cellAbsId;
 156       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 157     }
 158   }// cell loop
 159
 160   //Get from the absid the supermodule, tower and eta/phi numbers
 161   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
 162   //Gives SuperModule and Tower numbers
 163   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 164                                          iIphi, iIeta,iphi,ieta);
 165
 166 }
 167
 168 //__________________________________________________
 169 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeoUtils *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t iParticle)
 170 {
 171   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
 172
 173   Double_t eMax       = -1.;
 174   Double_t eCell      = -1.;
 175   Int_t    cellAbsId  = -1;
 176
 177   Int_t maxId   = -1;
 178   Int_t iTower  = -1;
 179   Int_t iIphi   = -1;
 180   Int_t iIeta   = -1;
 181         Int_t iSupMod = -1;
 182   Int_t iphi = -1, ieta =-1;
 183
 184   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
 185   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
 186   Bool_t  areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
 187   Int_t   startingSM = -1;
 188
 189   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 190     cellAbsId = clu->GetCellAbsId(iDig);
 191     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 192     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 193     if     (iDig==0)  startingSM = iSupMod;
 194     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 195
 196     eCell  = cells->GetCellAmplitude(cellAbsId);
 197
 198     weight = TMath::Log(eCell/clEnergy) + 4;
 199     if(weight < 0) weight = 0;
 200     totalWeight += weight;
 201     weightedCol += ieta*weight;
 202     weightedRow += iphi*weight;
 203
 204     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
 205
 206     if(eCell > eMax)  {
 207       eMax  = eCell;
 208       maxId = cellAbsId;
 209       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
 210     }
 211   }// cell loop
 212
 213   //Get from the absid the supermodule, tower and eta/phi numbers
 214   geom->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta);
 215   //Gives SuperModule and Tower numbers
 216   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
 217                                          iIphi, iIeta,iphi,ieta);
 218
 219   Float_t xyzNew[]={0.,0.,0.};
 220   if(areInSameSM == kTRUE) {
 221     //printf("In Same SM\n");
 222     weightedCol = weightedCol/totalWeight;
 223     weightedRow = weightedRow/totalWeight;
 224     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupMod, clEnergy, iParticle, fMisalShift, xyzNew);
 225   }
 226   else {
 227     //printf("In Different SM\n");
 228     geom->RecalculateTowerPosition(iphi,        ieta,        iSupMod, clEnergy, iParticle, fMisalShift, xyzNew);
 229   }
 230   clu->SetPosition(xyzNew);
 231
 232
 233   //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
 234
 235 }
 236
 237 //__________________________________________________
 238 void AliEMCALRecoUtils::Print(const Option_t *) const
 239 {
 240   // Print Parameters
 241
 242   printf("AliEMCALRecoUtils Settings: \n");
 243   printf("Misalignment shifts\n");
 244   for(Int_t i=0; i<5; i++) printf("\t sector %d, (dx,dy,dz)=(%f,%f,%f)\n",i, fMisalShift[i*3],fMisalShift[i*3+1],fMisalShift[i*3+2]);
 245   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
 246   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
 247
 248 }