EMCAL/AliEMCALAodCluster.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 //_________________________________________________________________________
  17 //  AliAODCaloCluster extension for EMCAL to recalculate cluster
  18 //  parameters in case of recalibration.
  19 //  Copy-paste from methods in AliEMCALRecPoint.
  20 //*--
  21 //*-- Author: Dmitri Peressounko (RRC KI) for PHOS
  22 //*-- Adapted for EMCAL: Gustavo Conesa (INFN-LNF)
  23
  24 // --- ROOT system ---
  25 #include "TVector3.h"
  26 #include "TMath.h"
  27
  28 // --- Standard library ---
  29
  30 // --- AliRoot header files ---
  31 #include "AliLog.h"
  32 #include "AliEMCALGeometry.h"
  33 #include "AliEMCALPID.h"
  34 #include "AliEMCALAodCluster.h"
  35 #include "AliEMCALCalibData.h"
  36 #include "AliAODCaloCells.h"
  37
  38 ClassImp(AliEMCALAodCluster)
  39
  40 //____________________________________________________________________________
  41 AliEMCALAodCluster::AliEMCALAodCluster() :
  42   AliAODCaloCluster(),fRecalibrated(0)
  43 {
  44   // ctor
  45 }
  46 //____________________________________________________________________________
  47 AliEMCALAodCluster::AliEMCALAodCluster(const AliAODCaloCluster & clu) :
  48   AliAODCaloCluster(clu),fRecalibrated(0)
  49 {
  50   // cpy ctor
  51 }
  52
  53 //____________________________________________________________________________
  54 AliEMCALAodCluster::~AliEMCALAodCluster()
  55 {
  56   // dtor
  57 }
  58 //____________________________________________________________________________
  59 void AliEMCALAodCluster::Recalibrate(AliEMCALCalibData * calibData, AliAODCaloCells *emcCells, TString emcalGeoName){
  60   //If not done yet, apply recalibration coefficients to energies list
  61   //NOTE that after recalibration fCellsAmpFraction contains not FRACTION but FULL energy
  62
  63   if(fRecalibrated)
  64    return ;
  65
  66   if(!calibData)
  67     return ;
  68
  69   AliEMCALGeometry * emcalgeo =  AliEMCALGeometry::GetInstance(emcalGeoName) ;
  70   if(emcalgeo){
  71
  72   Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
  73   Int_t iSupMod = -1;
  74   Int_t iTower  = -1;
  75   Int_t iIphi   = -1;
  76   Int_t iIeta   = -1;
  77   Int_t iphi    = -1;
  78   Int_t ieta    = -1;
  79
  80   for(Int_t i=0; i < GetNCells(); i++){
  81
  82     //Get from the absid the supermodule, tower and eta/phi numbers
  83     emcalgeo->GetCellIndex(GetCellAbsId(i),iSupMod,iTower,iIphi,iIeta);
  84     //Gives SuperModule and Tower numbers
  85     emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
  86                                           iIphi, iIeta,iphi,ieta);
  87
  88     Double_t energy = emcCells->GetCellAmplitude(GetCellAbsId(i)) ;
  89     //AliDebug(2,Form("Recalibrate: cell %f, calib %f, fraction %f\n",energy,calibData->GetADCchannel(iSupMod,ieta,iphi),cellsAmpFraction[i]));
  90         if(cellsAmpFraction[i]< 1e-4) cellsAmpFraction[i] = 1; //Unfolding off, took all the cell energy
  91         cellsAmpFraction[i]*=energy*calibData->GetADCchannel(iSupMod,ieta,iphi);
  92   }
  93
  94   SetCellsAmplitudeFraction(cellsAmpFraction);
  95   fRecalibrated=kTRUE;
  96
  97   }
  98   else  AliFatal("AliEMCALGeometry was not constructed\n") ;
  99
 100 }
 101 //____________________________________________________________________________
 102 void  AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
 103     //If recalibrated - recalculate all cluster parameters
 104   if(!fRecalibrated)
 105     return ;
 106   //printf("EvalAll e org %f\n",E());
 107   EvalEnergy() ; //Energy should be evaluated first
 108   //printf("EvalAll e2 %f\n",E());
 109   EvalPositionAndShowerShape(logWeight, geoname) ;
 110   //printf("EvalAll e3 %f\n",E());
 111   EvalPID() ; //Should be evaluated after energy and shower shape recalculation
 112   //printf("EvalAll e4 %f\n",E());
 113 }
 114 //____________________________________________________________________________
 115 void AliEMCALAodCluster::EvalEnergy(){
 116   //Evaluate energy
 117   if(!fRecalibrated) // no need to recalibrate
 118     return ;
 119
 120   Float_t energy=0. ;
 121   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 122     energy+=GetCellAmplitudeFraction(iDigit) ;
 123   }
 124   //printf("EvalEnergy: e %f\n", energy);
 125   SetE(energy);
 126
 127
 128 }
 129 ////____________________________________________________________________________
 130 //void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
 131 //  //apply nonlinearity correction same as in AliEMCALPID.
 132 //  SetE(pid->GetCalibratedEnergy(E())) ;
 133 //}
 134
 135 //____________________________________________________________________________
 136 void AliEMCALAodCluster::EvalPID(){
 137
 138   //re-evaluate identification parameters
 139 //  pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
 140 //  pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
 141
 142   //With bayesian
 143   AliEMCALPID *pid = new AliEMCALPID(kFALSE);
 144   pid->SetLowFluxParam(); // Need to be fixed
 145   Float_t pidlist[AliPID::kSPECIESN+1];
 146         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = pid->GetPIDFinal(i);
 147   SetPIDFromESD(pidlist);
 148
 149 }
 150
 151 //____________________________________________________________________________
 152 void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
 153 {
 154   // Calculates new center of gravity in the local EMCAL-module coordinates
 155   // and tranfers into global ALICE coordinates
 156   // Calculates Dispersion and main axis
 157   if(!fRecalibrated) // no need to recalibrate
 158     return ;
 159
 160   Int_t nstat  = 0;
 161   Float_t wtot = 0. ;
 162
 163   Int_t iSupMod = -1;
 164   Int_t iTower  = -1;
 165   Int_t iIphi   = -1;
 166   Int_t iIeta   = -1;
 167   Int_t iphi    = -1;
 168   Int_t ieta    = -1;
 169   Double_t etai = -1.;
 170   Double_t phii = -1.;
 171
 172   Double_t clXYZ[3] ={0.,0.,0.};
 173   Double_t xyzi[3]  ={0.,0.,0.};
 174
 175   Double_t d     = 0.;
 176   Double_t dxx   = 0.;
 177   Double_t dzz   = 0.;
 178   Double_t dxz   = 0.;
 179   Double_t xmean = 0.;
 180   Double_t zmean = 0.;
 181
 182   AliEMCALGeometry * emcalgeo =  AliEMCALGeometry::GetInstance(emcalGeoName) ;
 183   if(emcalgeo){
 184
 185   Double_t dist  = TmaxInCm(E(),0);
 186   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 187
 188         //Get from the absid the supermodule, tower and eta/phi numbers
 189         emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 190         emcalgeo->RelPosCellInSModule(GetCellAbsId(iDigit), dist, xyzi[0], xyzi[1], xyzi[2]);
 191         emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 192
 193     Double_t ei = GetCellAmplitudeFraction(iDigit) ;
 194     if (E() > 0 && ei > 0) {
 195                 Float_t w = ei;
 196                 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
 197                 etai=(Double_t)ieta;
 198                 phii=(Double_t)iphi;
 199                 if(w > 0.0) {
 200                         wtot += w ;
 201                         nstat++;
 202                         for(Int_t i = 0; i < 3; i++ ) clXYZ[i]    += (w*xyzi[i]);
 203
 204                         //Shower shape
 205                         dxx  += w * etai * etai ;
 206                         xmean+= w * etai ;
 207                         dzz  += w * phii * phii ;
 208                         zmean+= w * phii ;
 209                         dxz  += w * etai * phii ;
 210                 }
 211     }
 212     else
 213       AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
 214   }
 215
 216   //Normalize to the weight
 217   if (wtot > 0) {
 218                 for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
 219                 xmean /= wtot ;
 220                 zmean /= wtot ;
 221   }
 222   else
 223     AliError(Form("Wrong weight %f\n", wtot));
 224
 225   //Put cluster position in the global system
 226   TVector3 gpos ;
 227   emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);
 228
 229   SetPositionAt(gpos[0],0) ;
 230   SetPositionAt(gpos[1],1) ;
 231   SetPositionAt(gpos[2],2) ;
 232
 233   //Calculate dispersion
 234   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 235                 //Get from the absid the supermodule, tower and eta/phi numbers
 236                 emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 237                 emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 238
 239                 Double_t ei=GetCellAmplitudeFraction(iDigit) ;
 240                 if (E() > 0 && ei > 0) {
 241                         Float_t w = ei;
 242                         if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
 243                         etai=(Double_t)ieta;
 244                         phii=(Double_t)iphi;
 245                         if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
 246                 }
 247                 else
 248                         AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
 249   }
 250
 251   //Normalize to the weigth and set shower shape parameters
 252   if (wtot > 0 && nstat > 1) {
 253     d /= wtot ;
 254     dxx /= wtot ;
 255     dzz /= wtot ;
 256     dxz /= wtot ;
 257     dxx -= xmean * xmean ;
 258     dzz -= zmean * zmean ;
 259     dxz -= xmean * zmean ;
 260     SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
 261     SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 262   }
 263   else{
 264     d=0. ;
 265     SetM20(0.) ;
 266     SetM02(0.) ;
 267   }
 268
 269   if (d>=0)
 270           SetDispersion(TMath::Sqrt(d)) ;
 271   else
 272           SetDispersion(0) ;
 273   }
 274   else {
 275     AliFatal("AliEMCALGeometry was not constructed\n") ;
 276   }
 277
 278
 279 }
 280
 281 //_____________________________________________________________________
 282 Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
 283 {
 284         // e energy in GeV)
 285         // key  =  0(gamma, default)
 286         //     !=  0(electron)
 287         static Double_t ca = 4.82;  // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
 288         static Double_t x0 = 1.23;  // radiation lenght (cm)
 289         static Double_t tmax = 0.;   // position of electromagnetic shower max in cm
 290
 291         tmax = 0.0;
 292         if(e>0.1) {
 293                 tmax = TMath::Log(e) + ca;
 294                 if      (key==0) tmax += 0.5;
 295                 else             tmax -= 0.5;
 296                 tmax *= x0; // convert to cm
 297         }
 298         return tmax;
 299 }
 300