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) ;
72 Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
80 for(Int_t i=0; i < GetNCells(); i++){
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);
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);
94 SetCellsAmplitudeFraction(cellsAmpFraction);
98 else AliFatal("AliEMCALGeometry was not constructed\n") ;
101 //____________________________________________________________________________
102 void AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
103 //If recalibrated - recalculate all cluster parameters
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());
114 //____________________________________________________________________________
115 void AliEMCALAodCluster::EvalEnergy(){
117 if(!fRecalibrated) // no need to recalibrate
121 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
122 energy+=GetCellAmplitudeFraction(iDigit) ;
124 //printf("EvalEnergy: e %f\n", energy);
129 ////____________________________________________________________________________
130 //void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
131 // //apply nonlinearity correction same as in AliEMCALPID.
132 // SetE(pid->GetCalibratedEnergy(E())) ;
135 //____________________________________________________________________________
136 void AliEMCALAodCluster::EvalPID(){
138 //re-evaluate identification parameters
139 // pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
140 // pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
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);
151 //____________________________________________________________________________
152 void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
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
172 Double_t clXYZ[3] ={0.,0.,0.};
173 Double_t xyzi[3] ={0.,0.,0.};
182 AliEMCALGeometry * emcalgeo = AliEMCALGeometry::GetInstance(emcalGeoName) ;
185 Double_t dist = TmaxInCm(E(),0);
186 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
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);
193 Double_t ei = GetCellAmplitudeFraction(iDigit) ;
194 if (E() > 0 && ei > 0) {
196 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
202 for(Int_t i = 0; i < 3; i++ ) clXYZ[i] += (w*xyzi[i]);
205 dxx += w * etai * etai ;
207 dzz += w * phii * phii ;
209 dxz += w * etai * phii ;
213 AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
216 //Normalize to the weight
218 for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
223 AliError(Form("Wrong weight %f\n", wtot));
225 //Put cluster position in the global system
227 emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);
229 SetPositionAt(gpos[0],0) ;
230 SetPositionAt(gpos[1],1) ;
231 SetPositionAt(gpos[2],2) ;
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);
239 Double_t ei=GetCellAmplitudeFraction(iDigit) ;
240 if (E() > 0 && ei > 0) {
242 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
245 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
248 AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
251 //Normalize to the weigth and set shower shape parameters
252 if (wtot > 0 && nstat > 1) {
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 ));
270 SetDispersion(TMath::Sqrt(d)) ;
275 AliFatal("AliEMCALGeometry was not constructed\n") ;
281 //_____________________________________________________________________
282 Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
285 // key = 0(gamma, default)
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
293 tmax = TMath::Log(e) + ca;
294 if (key==0) tmax += 0.5;
296 tmax *= x0; // convert to cm