]>
Commit | Line | Data |
---|---|---|
0c5b726e | 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 | ||
0c5b726e | 159 | Int_t iSupMod = -1; |
160 | Int_t iTower = -1; | |
161 | Int_t iIphi = -1; | |
162 | Int_t iIeta = -1; | |
163 | Int_t iphi = -1; | |
164 | Int_t ieta = -1; | |
165 | Double_t etai = -1.; | |
166 | Double_t phii = -1.; | |
167 | ||
168 | Double_t clXYZ[3] ={0.,0.,0.}; | |
169 | Double_t xyzi[3] ={0.,0.,0.}; | |
170 | ||
171 | Double_t d = 0.; | |
172 | Double_t dxx = 0.; | |
173 | Double_t dzz = 0.; | |
174 | Double_t dxz = 0.; | |
175 | Double_t xmean = 0.; | |
176 | Double_t zmean = 0.; | |
177 | ||
178 | AliEMCALGeometry * emcalgeo = AliEMCALGeometry::GetInstance(emcalGeoName) ; | |
179 | if(!emcalgeo) | |
180 | AliFatal("AliEMCALGeometry was not constructed\n") ; | |
181 | ||
7cfcebd3 | 182 | Double_t dist = TmaxInCm(Double_t(GetCellAmplitudeFraction(0)),0); |
0c5b726e | 183 | for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) { |
7cfcebd3 | 184 | |
0c5b726e | 185 | //Get from the absid the supermodule, tower and eta/phi numbers |
186 | emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); | |
7cfcebd3 | 187 | emcalgeo->RelPosCellInSModule(GetCellAbsId(iDigit), dist, xyzi[0], xyzi[1], xyzi[2]); |
0c5b726e | 188 | emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); |
189 | ||
190 | Double_t ei = GetCellAmplitudeFraction(iDigit) ; | |
191 | if (E() > 0 && ei > 0) { | |
192 | Float_t w = ei; | |
193 | if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ; | |
194 | etai=(Double_t)ieta; | |
195 | phii=(Double_t)iphi; | |
196 | if(w > 0.0) { | |
197 | wtot += w ; | |
198 | nstat++; | |
199 | for(Int_t i = 0; i < 3; i++ ) clXYZ[i] += (w*xyzi[i]); | |
200 | ||
201 | //Shower shape | |
202 | dxx += w * etai * etai ; | |
203 | xmean+= w * etai ; | |
204 | dzz += w * phii * phii ; | |
205 | zmean+= w * phii ; | |
206 | dxz += w * etai * phii ; | |
207 | } | |
208 | } | |
209 | else | |
210 | AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E())); | |
211 | } | |
212 | ||
213 | //Normalize to the weight | |
214 | if (wtot > 0) { | |
215 | for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot; | |
216 | xmean /= wtot ; | |
217 | zmean /= wtot ; | |
218 | } | |
219 | else | |
220 | AliError(Form("Wrong weight %f\n", wtot)); | |
221 | ||
222 | //Put cluster position in the global system | |
223 | TVector3 gpos ; | |
224 | emcalgeo->GetGlobal(clXYZ, gpos, iSupMod); | |
225 | ||
226 | SetPosition(0, gpos[0]) ; | |
227 | SetPosition(1, gpos[1]) ; | |
228 | SetPosition(2, gpos[2]) ; | |
229 | ||
230 | //Calculate dispersion | |
231 | for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) { | |
232 | //Get from the absid the supermodule, tower and eta/phi numbers | |
233 | emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); | |
234 | emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); | |
235 | ||
236 | Double_t ei=GetCellAmplitudeFraction(iDigit) ; | |
237 | if (E() > 0 && ei > 0) { | |
238 | Float_t w = ei; | |
239 | if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ; | |
240 | etai=(Double_t)ieta; | |
241 | phii=(Double_t)iphi; | |
242 | if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); | |
243 | } | |
244 | else | |
245 | AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E())); | |
246 | } | |
247 | ||
248 | //Normalize to the weigth and set shower shape parameters | |
249 | if (wtot > 0 && nstat > 1) { | |
250 | d /= wtot ; | |
251 | dxx /= wtot ; | |
252 | dzz /= wtot ; | |
253 | dxz /= wtot ; | |
254 | dxx -= xmean * xmean ; | |
255 | dzz -= zmean * zmean ; | |
256 | dxz -= xmean * zmean ; | |
257 | SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ; | |
258 | SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )); | |
259 | } | |
260 | else{ | |
261 | d=0. ; | |
262 | SetM20(0.) ; | |
263 | SetM02(0.) ; | |
264 | } | |
265 | ||
266 | if (d>=0) | |
267 | SetDispersion(TMath::Sqrt(d)) ; | |
268 | else | |
269 | SetDispersion(0) ; | |
270 | ||
271 | } | |
272 | ||
273 | //_____________________________________________________________________ | |
274 | Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const | |
275 | { | |
276 | // e energy in GeV) | |
277 | // key = 0(gamma, default) | |
278 | // != 0(electron) | |
279 | static Double_t ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07) | |
280 | static Double_t x0 = 1.23; // radiation lenght (cm) | |
281 | static Double_t tmax = 0.; // position of electromagnetic shower max in cm | |
282 | ||
283 | tmax = 0.0; | |
284 | if(e>0.1) { | |
285 | tmax = TMath::Log(e) + ca; | |
286 | if (key==0) tmax += 0.5; | |
287 | else tmax -= 0.5; | |
288 | tmax *= x0; // convert to cm | |
289 | } | |
290 | return tmax; | |
291 | } | |
292 |