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 cellsAmpFraction[i]*=energy*calibData->GetADCchannel(iSupMod,ieta,iphi);
94 SetCellsAmplitudeFraction(cellsAmpFraction);
97 //____________________________________________________________________________
98 void AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
99 //If recalibrated - recalculate all cluster parameters
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());
110 //____________________________________________________________________________
111 void AliEMCALAodCluster::EvalEnergy(){
113 if(!fRecalibrated) // no need to recalibrate
117 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
118 energy+=GetCellAmplitudeFraction(iDigit) ;
120 //printf("EvalEnergy: e %f\n", energy);
125 ////____________________________________________________________________________
126 //void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
127 // //apply nonlinearity correction same as in AliEMCALPID.
128 // SetE(pid->GetCalibratedEnergy(E())) ;
131 //____________________________________________________________________________
132 void AliEMCALAodCluster::EvalPID(){
134 //re-evaluate identification parameters
135 // pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
136 // pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
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);
147 //____________________________________________________________________________
148 void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
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
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") ;
184 for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
186 //Check if this maximum at 0 is true!!
187 idMax = GetCellAbsId(iDigit);
188 dist = TmaxInCm(Double_t(GetCellAmplitudeFraction(iDigit)),0);
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);
196 Double_t ei = GetCellAmplitudeFraction(iDigit) ;
197 if (E() > 0 && ei > 0) {
199 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
205 for(Int_t i = 0; i < 3; i++ ) clXYZ[i] += (w*xyzi[i]);
208 dxx += w * etai * etai ;
210 dzz += w * phii * phii ;
212 dxz += w * etai * phii ;
216 AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
219 //Normalize to the weight
221 for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
226 AliError(Form("Wrong weight %f\n", wtot));
228 //Put cluster position in the global system
230 emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);
232 SetPosition(0, gpos[0]) ;
233 SetPosition(1, gpos[1]) ;
234 SetPosition(2, gpos[2]) ;
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);
242 Double_t ei=GetCellAmplitudeFraction(iDigit) ;
243 if (E() > 0 && ei > 0) {
245 if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
248 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
251 AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
254 //Normalize to the weigth and set shower shape parameters
255 if (wtot > 0 && nstat > 1) {
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 ));
273 SetDispersion(TMath::Sqrt(d)) ;
279 //_____________________________________________________________________
280 Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
283 // key = 0(gamma, default)
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
291 tmax = TMath::Log(e) + ca;
292 if (key==0) tmax += 0.5;
294 tmax *= x0; // convert to cm