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     AliFatal("AliEMCALGeometry was not constructed\n") ;
  72
  73   Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
  74   Int_t iSupMod = -1;
  75   Int_t iTower  = -1;
  76   Int_t iIphi   = -1;
  77   Int_t iIeta   = -1;
  78   Int_t iphi    = -1;
  79   Int_t ieta    = -1;
  80
  81   for(Int_t i=0; i < GetNCells(); i++){
  82
  83     //Get from the absid the supermodule, tower and eta/phi numbers
  84     emcalgeo->GetCellIndex(GetCellAbsId(i),iSupMod,iTower,iIphi,iIeta);
  85     //Gives SuperModule and Tower numbers
  86     emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
  87                                           iIphi, iIeta,iphi,ieta);
  88
  89     Double_t energy = emcCells->GetCellAmplitude(GetCellAbsId(i)) ;
  90     AliDebug(2,Form("Recalibrate: cell %f, calib %f, fraction %f\n",energy,calibData->GetADCchannel(iSupMod,ieta,iphi),cellsAmpFraction[i]));
  91         cellsAmpFraction[i]*=energy*calibData->GetADCchannel(iSupMod,ieta,iphi);
  92   }
  93
  94   SetCellsAmplitudeFraction(cellsAmpFraction);
  95   fRecalibrated=kTRUE;
  96 }
  97 //____________________________________________________________________________
  98 void  AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
  99     //If recalibrated - recalculate all cluster parameters
 100   if(!fRecalibrated)
 101     return ;
 102   //printf("EvalAll e org %f\n",E());
 103   EvalEnergy() ; //Energy should be evaluated first
 104   //printf("EvalAll e2 %f\n",E());
 105   EvalPositionAndShowerShape(logWeight, geoname) ;
 106   //printf("EvalAll e3 %f\n",E());
 107   EvalPID() ; //Should be evaluated after energy and shower shape recalculation
 108   //printf("EvalAll e4 %f\n",E());
 109 }
 110 //____________________________________________________________________________
 111 void AliEMCALAodCluster::EvalEnergy(){
 112   //Evaluate energy
 113   if(!fRecalibrated) // no need to recalibrate
 114     return ;
 115
 116   Float_t energy=0. ;
 117   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 118     energy+=GetCellAmplitudeFraction(iDigit) ;
 119   }
 120   //printf("EvalEnergy: e %f\n", energy);
 121   SetE(energy);
 122
 123
 124 }
 125 ////____________________________________________________________________________
 126 //void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
 127 //  //apply nonlinearity correction same as in AliEMCALPID.
 128 //  SetE(pid->GetCalibratedEnergy(E())) ;
 129 //}
 130
 131 //____________________________________________________________________________
 132 void AliEMCALAodCluster::EvalPID(){
 133
 134   //re-evaluate identification parameters
 135 //  pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
 136 //  pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
 137
 138   //With bayesian
 139   AliEMCALPID *pid = new AliEMCALPID(kFALSE);
 140   pid->SetLowFluxParam(); // Need to be fixed
 141   Float_t pidlist[AliPID::kSPECIESN+1];
 142         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = pid->GetPIDFinal(i);
 143   SetPIDFromESD(pidlist);
 144
 145 }
 146
 147 //____________________________________________________________________________
 148 void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
 149 {
 150   // Calculates new center of gravity in the local EMCAL-module coordinates
 151   // and tranfers into global ALICE coordinates
 152   // Calculates Dispersion and main axis
 153   if(!fRecalibrated) // no need to recalibrate
 154     return ;
 155
 156   Int_t nstat  = 0;
 157   Float_t wtot = 0. ;
 158
 159   Int_t idMax   = -1;
 160   Int_t iSupMod = -1;
 161   Int_t iTower  = -1;
 162   Int_t iIphi   = -1;
 163   Int_t iIeta   = -1;
 164   Int_t iphi    = -1;
 165   Int_t ieta    = -1;
 166   Double_t etai = -1.;
 167   Double_t phii = -1.;
 168
 169   Double_t clXYZ[3] ={0.,0.,0.};
 170   Double_t xyzi[3]  ={0.,0.,0.};
 171
 172   Double_t d     = 0.;
 173   Double_t dxx   = 0.;
 174   Double_t dzz   = 0.;
 175   Double_t dxz   = 0.;
 176   Double_t xmean = 0.;
 177   Double_t zmean = 0.;
 178
 179   AliEMCALGeometry * emcalgeo =  AliEMCALGeometry::GetInstance(emcalGeoName) ;
 180   if(!emcalgeo)
 181     AliFatal("AliEMCALGeometry was not constructed\n") ;
 182
 183   Double_t dist = 0;
 184   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 185         if(iDigit==0) {
 186                 //Check if this maximum  at 0 is true!!
 187                 idMax = GetCellAbsId(iDigit);
 188                 dist  = TmaxInCm(Double_t(GetCellAmplitudeFraction(iDigit)),0);
 189         }
 190
 191         //Get from the absid the supermodule, tower and eta/phi numbers
 192         emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 193         emcalgeo->RelPosCellInSModule(GetCellAbsId(iDigit),idMax, dist, xyzi[0], xyzi[1], xyzi[2]);
 194         emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 195
 196     Double_t ei = GetCellAmplitudeFraction(iDigit) ;
 197     if (E() > 0 && ei > 0) {
 198                 Float_t w = ei;
 199                 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
 200                 etai=(Double_t)ieta;
 201                 phii=(Double_t)iphi;
 202                 if(w > 0.0) {
 203                         wtot += w ;
 204                         nstat++;
 205                         for(Int_t i = 0; i < 3; i++ ) clXYZ[i]    += (w*xyzi[i]);
 206
 207                         //Shower shape
 208                         dxx  += w * etai * etai ;
 209                         xmean+= w * etai ;
 210                         dzz  += w * phii * phii ;
 211                         zmean+= w * phii ;
 212                         dxz  += w * etai * phii ;
 213                 }
 214     }
 215     else
 216       AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
 217   }
 218
 219   //Normalize to the weight
 220   if (wtot > 0) {
 221                 for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
 222                 xmean /= wtot ;
 223                 zmean /= wtot ;
 224   }
 225   else
 226     AliError(Form("Wrong weight %f\n", wtot));
 227
 228   //Put cluster position in the global system
 229   TVector3 gpos ;
 230   emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);
 231
 232   SetPosition(0, gpos[0]) ;
 233   SetPosition(1, gpos[1]) ;
 234   SetPosition(2, gpos[2]) ;
 235
 236   //Calculate dispersion
 237   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 238                 //Get from the absid the supermodule, tower and eta/phi numbers
 239                 emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
 240                 emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
 241
 242                 Double_t ei=GetCellAmplitudeFraction(iDigit) ;
 243                 if (E() > 0 && ei > 0) {
 244                         Float_t w = ei;
 245                         if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
 246                         etai=(Double_t)ieta;
 247                         phii=(Double_t)iphi;
 248                         if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
 249                 }
 250                 else
 251                         AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
 252   }
 253
 254   //Normalize to the weigth and set shower shape parameters
 255   if (wtot > 0 && nstat > 1) {
 256     d /= wtot ;
 257     dxx /= wtot ;
 258     dzz /= wtot ;
 259     dxz /= wtot ;
 260     dxx -= xmean * xmean ;
 261     dzz -= zmean * zmean ;
 262     dxz -= xmean * zmean ;
 263     SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
 264     SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 265   }
 266   else{
 267     d=0. ;
 268     SetM20(0.) ;
 269     SetM02(0.) ;
 270   }
 271
 272   if (d>=0)
 273           SetDispersion(TMath::Sqrt(d)) ;
 274   else
 275           SetDispersion(0) ;
 276
 277 }
 278
 279 //_____________________________________________________________________
 280 Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
 281 {
 282         // e energy in GeV)
 283         // key  =  0(gamma, default)
 284         //     !=  0(electron)
 285         static Double_t ca = 4.82;  // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
 286         static Double_t x0 = 1.23;  // radiation lenght (cm)
 287         static Double_t tmax = 0.;   // position of electromagnetic shower max in cm
 288
 289         tmax = 0.0;
 290         if(e>0.1) {
 291                 tmax = TMath::Log(e) + ca;
 292                 if      (key==0) tmax += 0.5;
 293                 else             tmax -= 0.5;
 294                 tmax *= x0; // convert to cm
 295         }
 296         return tmax;
 297 }
 298