1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // AliAODCaloCluster extension for EMCAL to recalculate cluster
18 // parameters in case of recalibration.
19 // Copy-paste from methods in AliEMCALRecPoint.
21 //*-- Author: Dmitri Peressounko (RRC KI) for PHOS
22 //*-- Adapted for EMCAL: Gustavo Conesa (INFN-LNF)
24 // --- ROOT system ---
28 // --- Standard library ---
30 // --- AliRoot header files ---
32 #include "AliEMCALGeometry.h"
33 #include "AliEMCALPID.h"
34 #include "AliEMCALAodCluster.h"
35 #include "AliEMCALCalibData.h"
36 #include "AliAODCaloCells.h"
38 ClassImp(AliEMCALAodCluster)
40 //____________________________________________________________________________
41 AliEMCALAodCluster::AliEMCALAodCluster() :
42 AliAODCaloCluster(),fRecalibrated(0)
46 //____________________________________________________________________________
47 AliEMCALAodCluster::AliEMCALAodCluster(const AliAODCaloCluster & clu) :
48 AliAODCaloCluster(clu),fRecalibrated(0)
53 //____________________________________________________________________________
54 AliEMCALAodCluster::~AliEMCALAodCluster()
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
69 AliEMCALGeometry * emcalgeo = AliEMCALGeometry::GetInstance(emcalGeoName) ;
71 AliFatal("AliEMCALGeometry was not constructed\n") ;
73 Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
81 for(Int_t i=0; i < GetNCells(); i++){
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);
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 if(cellsAmpFraction[i]< 1e-4) cellsAmpFraction[i] = 1; //Unfolding off, took all the cell energy
92 cellsAmpFraction[i]*=energy*calibData->GetADCchannel(iSupMod,ieta,iphi);
95 SetCellsAmplitudeFraction(cellsAmpFraction);
98 //____________________________________________________________________________
99 void AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
100 //If recalibrated - recalculate all cluster parameters
103 //printf("EvalAll e org %f\n",E());
104 EvalEnergy() ; //Energy should be evaluated first
105 //printf("EvalAll e2 %f\n",E());
106 EvalPositionAndShowerShape(logWeight, geoname) ;
107 //printf("EvalAll e3 %f\n",E());
108 EvalPID() ; //Should be evaluated after energy and shower shape recalculation
109 //printf("EvalAll e4 %f\n",E());
111 //____________________________________________________________________________
112 void AliEMCALAodCluster::EvalEnergy(){
114 if(!fRecalibrated) // no need to recalibrate
118 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
119 energy+=GetCellAmplitudeFraction(iDigit) ;
121 //printf("EvalEnergy: e %f\n", energy);
126 ////____________________________________________________________________________
127 //void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
128 // //apply nonlinearity correction same as in AliEMCALPID.
129 // SetE(pid->GetCalibratedEnergy(E())) ;
132 //____________________________________________________________________________
133 void AliEMCALAodCluster::EvalPID(){
135 //re-evaluate identification parameters
136 // pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
137 // pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
140 AliEMCALPID *pid = new AliEMCALPID(kFALSE);
141 pid->SetLowFluxParam(); // Need to be fixed
142 Float_t pidlist[AliPID::kSPECIESN+1];
143 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = pid->GetPIDFinal(i);
144 SetPIDFromESD(pidlist);
148 //____________________________________________________________________________
149 void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
151 // Calculates new center of gravity in the local EMCAL-module coordinates
152 // and tranfers into global ALICE coordinates
153 // Calculates Dispersion and main axis
154 if(!fRecalibrated) // no need to recalibrate
169 Double_t clXYZ[3] ={0.,0.,0.};
170 Double_t xyzi[3] ={0.,0.,0.};
179 AliEMCALGeometry * emcalgeo = AliEMCALGeometry::GetInstance(emcalGeoName) ;
181 AliFatal("AliEMCALGeometry was not constructed\n") ;
183 Double_t dist = TmaxInCm(E(),0);
184 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
186 //Get from the absid the supermodule, tower and eta/phi numbers
187 emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
188 emcalgeo->RelPosCellInSModule(GetCellAbsId(iDigit), dist, xyzi[0], xyzi[1], xyzi[2]);
189 emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
191 Double_t ei = GetCellAmplitudeFraction(iDigit) ;
192 if (E() > 0 && ei > 0) {
194 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
200 for(Int_t i = 0; i < 3; i++ ) clXYZ[i] += (w*xyzi[i]);
203 dxx += w * etai * etai ;
205 dzz += w * phii * phii ;
207 dxz += w * etai * phii ;
211 AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
214 //Normalize to the weight
216 for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
221 AliError(Form("Wrong weight %f\n", wtot));
223 //Put cluster position in the global system
225 emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);
227 SetPosition(0,gpos[0]) ;
228 SetPosition(1,gpos[1]) ;
229 SetPosition(2,gpos[2]) ;
231 //Calculate dispersion
232 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
233 //Get from the absid the supermodule, tower and eta/phi numbers
234 emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
235 emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
237 Double_t ei=GetCellAmplitudeFraction(iDigit) ;
238 if (E() > 0 && ei > 0) {
240 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
243 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
246 AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
249 //Normalize to the weigth and set shower shape parameters
250 if (wtot > 0 && nstat > 1) {
255 dxx -= xmean * xmean ;
256 dzz -= zmean * zmean ;
257 dxz -= xmean * zmean ;
258 SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
259 SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
268 SetDispersion(TMath::Sqrt(d)) ;
274 //_____________________________________________________________________
275 Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
278 // key = 0(gamma, default)
280 static Double_t ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
281 static Double_t x0 = 1.23; // radiation lenght (cm)
282 static Double_t tmax = 0.; // position of electromagnetic shower max in cm
286 tmax = TMath::Log(e) + ca;
287 if (key==0) tmax += 0.5;
289 tmax *= x0; // convert to cm