19 {
20 RecEmcFractionMap::const_iterator cit;
24
25
26
27 IEmcRecGeoSvc* iGeoSvc;
28 ISvcLocator* svcLocator = Gaudi::svcLocator();
29 StatusCode sc = svcLocator->service( "EmcRecGeoSvc", iGeoSvc );
30 if ( sc != StatusCode::SUCCESS )
31 {
32
33 }
34
36
38
39
40 double e5x5 = -99999;
42 {
44 for ( cit = aShower.
Begin(); cit != aShower.
End(); cit++ )
45 {
46 w = offset + log( ( cit->second.getEnergy() /
etot ) * cit->second.getFraction() );
48 {
52 }
53 }
54 }
56 {
57 double e3x3 = aShower.
e3x3();
59 for ( cit = fracMap3x3.begin(); cit != fracMap3x3.end(); cit++ )
60 {
61 w = offset + log( ( cit->second.getEnergy() / e3x3 ) * cit->second.getFraction() );
63 {
67 }
68 }
69 }
70 else
71 {
72 e5x5 = aShower.
e5x5();
74 for ( cit = fracMap5x5.begin(); cit != fracMap5x5.end(); cit++ )
75 {
76 w = offset + log( ( cit->second.getEnergy() / e5x5 ) * cit->second.getFraction() );
78 {
82 }
83 }
84 }
85 if ( wsum <= 0 )
86 {
87
88 }
89 pos = possum / wsum;
90
91
92
93
99
100 if ( module == 1 )
101 {
102
103 double IRShift = 0, parTheta[5], parPhi[5];
104
105 unsigned int theModule;
106 if ( thetaModule > 21 )
107 {
108 theModule = 43 - thetaModule;
110 posCorr.setZ( -posCorr.z() );
111 }
112 else { theModule = thetaModule; }
113
114 if ( theModule == 21 ) { IRShift = 0; }
115 else if ( theModule == 20 ) { IRShift = 2.5; }
116 else { IRShift = 5.; }
118
119
120
122 {
123
124 for ( int i = 0; i < 5; i++ )
125 {
126 if ( theModule == 21 )
127 {
128 double par[5] = { 51.97, -0.1209, 6.175, -0.1431, -0.005497 };
129 parTheta[i] = par[i];
130 }
131 else if ( theModule == 20 )
132 {
133 double par[5] = { 56.31, -0.1135, 6.312, -1.216, -0.02978 };
134 parTheta[i] = par[i];
135 }
136 else
137 {
138 double par[5] = { 10.65, -0.2257, 2.216, -0.1562, -0.05552 };
139 parTheta[i] = par[i];
140 }
141 }
142 }
144 {
145
146 for ( int i = 0; i < 5; i++ )
147 {
148
150 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
152 else if ( e5x5 <= 0.5 ) parTheta[i] = Para.
BarrLogThetaPara( theModule + 44, i );
153
154 if ( e5x5 > 1.0 ) parPhi[i] = Para.
BarrLogPhiPara( theModule, i );
155 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
157 else if ( e5x5 <= 0.5 ) parPhi[i] = Para.
BarrLogPhiPara( theModule + 44, i );
158 }
159 }
160
164
165 double theta01, theta23, length, d;
166 theta01 = ( center01 - IR ).theta();
167 theta23 = ( center23 - IR ).theta();
168 length = ( center01 - IR ).mag();
169 d = length *
tan( theta01 - theta23 );
170
171 double xIn, xOut, deltaTheta;
172 xIn = length *
tan( theta01 - ( posCorr - IR ).theta() ) - d / 2;
173 xOut = xIn - ( parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
174 parTheta[3] );
175 deltaTheta = atan( ( xOut + d / 2 ) / length ) - atan( ( xIn + d / 2 ) / length );
176
177
178
179 double yIn, deltaY, deltaPhi;
180
181
182 if ( phiModule == 0 ) { yIn = posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843; }
183 else if ( phiModule == 119 )
184 { yIn = posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843; }
185 else
186 {
187 yIn = posCorr.phi() < 0
188 ? posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
189 : posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
190 }
191
192
193 deltaY = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
194 deltaPhi = deltaY * CLHEP::pi / 180.;
195
197 rotateVector.rotateZ( center01.phi() );
198 posCorr.setZ( posCorr.z() - IRShift );
199 posCorr.rotate( -deltaTheta, rotateVector );
200 posCorr.setZ( posCorr.z() + IRShift );
201
202 if ( Para.
MethodMode() == 0 ) { posCorr.rotateZ( -0.002994 ); }
203 else if ( Para.
MethodMode() == 1 ) { posCorr.rotateZ( -deltaPhi ); }
204
205 if ( thetaModule > 21 ) { posCorr.setZ( -posCorr.z() ); }
206
207 double lengthemc;
208 lengthemc =
abs( posCorr.z() /
cos( posCorr.theta() ) );
209
210 double posDataCorPar;
212 {
214 posCorr.setTheta( posCorr.theta() - posDataCorPar / lengthemc );
215 }
216
217
218 double posMCCorPar;
220 {
221 posMCCorPar = Para.
BarrPosMCCor( thetaModule, phiModule );
222 posCorr.setTheta( posCorr.theta() - posMCCorPar / lengthemc );
223 }
224 }
225 else
226 {
228 {
229
230 double IRShift = 0, parTheta[5], parPhi[5];
231
232 if ( module == 0 )
233 {
234
235 IRShift = 10;
236
237 for ( int i = 0; i < 5; i++ )
238 {
239
241 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
243 else if ( e5x5 <= 0.5 ) parTheta[i] = Para.
EastLogThetaPara( thetaModule + 6, i );
244
246 else if ( e5x5 <= 1.0 && e5x5 > 0.5 ) parPhi[i] = Para.
EastLogPhiPara( 1, i );
248 }
249 }
250 else if ( module == 2 )
251 {
252 IRShift = -10;
253 for ( int i = 0; i < 5; i++ )
254 {
255
257 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
259 else if ( e5x5 <= 0.5 ) parTheta[i] = Para.
WestLogThetaPara( thetaModule + 6, i );
260
262 else if ( e5x5 <= 1.0 && e5x5 > 0.5 ) parPhi[i] = Para.
WestLogPhiPara( 1, i );
264 }
265 }
266
268
269 double theta12 = -9999, theta03 = -9999, theta23 = -9999, theta14 = -9999, delta = -9999;
270 double phi01 = -9999, phi23 = -9999, phi12 = -9999, phi34 = -9999, deltaphi = -9999;
271 double xIn, deltaX, deltaTheta, yIn, deltaY, deltaPhi;
272
273 double posphi = posCorr.phi();
274 if ( phiModule == 0 ) { posphi = posphi; }
275 else if ( ( ( thetaModule == 0 || thetaModule == 1 ) && phiModule == 63 ) ||
276 ( ( thetaModule == 2 || thetaModule == 3 ) && phiModule == 79 ) ||
277 ( ( thetaModule == 4 || thetaModule == 5 ) && phiModule == 95 ) )
278 { posphi = posphi + 2 * CLHEP::pi; }
279 else { posphi = posphi < 0 ? posphi + 2 * CLHEP::pi : posphi; }
280
281 HepPoint3D center, center01, center23, center12, center03, center34, center14;
282
283 if ( ( thetaModule == 1 &&
284 ( ( phiModule + 3 ) % 4 == 0 || ( phiModule + 2 ) % 4 == 0 ) ) ||
285 ( thetaModule == 3 && ( ( phiModule + 4 ) % 5 == 0 || ( phiModule + 3 ) % 5 == 0 ||
286 ( phiModule + 2 ) % 5 == 0 ) ) )
287 {
288
289 center23 =
291 center14 =
293 center = center23;
294 theta23 = ( center23 - IR ).theta();
295 theta14 = ( center14 - IR ).theta();
296 delta = theta14 - theta23;
297 xIn = ( ( ( posCorr - IR ).theta() - theta23 ) - delta * 0.5 ) / delta;
298
299 center12 =
301 center34 =
303 phi12 = center12.phi();
304 phi34 = center34.phi();
305
306 phi12 = phi12 < 0 ? phi12 + 2 * CLHEP::pi : phi12;
307 phi34 = phi34 < 0 ? phi34 + 2 * CLHEP::pi : phi34;
308 deltaphi = phi34 - phi12;
309 yIn = ( ( posphi - phi12 ) - deltaphi * 0.5 ) / deltaphi;
310 }
311 else
312 {
313
314 center12 =
316 center03 =
318 center = center12;
319
320 theta12 = ( center12 - IR ).theta();
321 theta03 = ( center03 - IR ).theta();
322 delta = theta03 - theta12;
323 xIn = ( ( ( posCorr - IR ).theta() - theta12 ) - delta * 0.5 ) / delta;
324
325 center01 =
327 center23 =
329 phi01 = center01.phi();
330 phi23 = center23.phi();
331 phi01 = phi01 < 0 ? phi01 + 2 * CLHEP::pi : phi01;
332 phi23 = phi23 < 0 ? phi23 + 2 * CLHEP::pi : phi23;
333 deltaphi = phi23 - phi01;
334 yIn = ( ( posphi - phi01 ) - deltaphi * 0.5 ) / deltaphi;
335 }
336
337 deltaX = parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
338 parTheta[3];
339 deltaTheta = deltaX * delta;
340
341 deltaY = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
342 deltaPhi = deltaY * deltaphi;
343
345 rotateVector.rotateZ( center.phi() );
346 posCorr.setZ( posCorr.z() - IRShift );
347 posCorr.rotate( -deltaTheta, rotateVector );
348 posCorr.setZ( posCorr.z() + IRShift );
349 posCorr.rotateZ( -deltaPhi );
350
351 double lengthemc;
352 lengthemc =
abs( posCorr.z() /
cos( posCorr.theta() ) );
353 double posDataCorPar = 0.0;
355 {
356 if ( module == 0 ) posDataCorPar = Para.
EastPosDataCor( thetaModule, phiModule );
357 if ( module == 2 ) posDataCorPar = Para.
WestPosDataCor( thetaModule, phiModule );
358 posCorr.setTheta( posCorr.theta() - posDataCorPar / lengthemc );
359 }
360 }
361 }
362
363
364
366
368
369 if ( aShower.
energy() < 1e-5 )
return;
370
371 double r, theta, phi;
372 double dtheta, dphi, dx, dy, dz;
373
374 r = posCorr.mag();
375 theta = posCorr.theta();
376 phi = posCorr.phi();
377
378
379 if ( aShower.
energy() > 0 )
380 {
383 }
384 else
385 {
386 dtheta = 0;
387 dphi = 0;
388 }
389
390
393 if ( theta > 0 )
394 dz = fabs( iGeoSvc->
GetBarrelR() * dtheta / (
sin( theta ) *
sin( theta ) ) );
395 else dz = 0;
396
397
400
401
402
403 HepSymMatrix matrix( 3 );
404 matrix[0][0] = 0;
405 matrix[1][1] = dtheta * dtheta;
406 matrix[2][2] = dphi * dphi;
407 matrix[0][1] = 0;
408 matrix[0][2] = 0;
409 matrix[1][2] = 0;
410
411
412 HepMatrix matrix1( 3, 3 ), matrix2( 3, 3 );
413 matrix1[0][0] =
sin( theta ) *
cos( phi );
414 matrix1[0][1] = r *
cos( theta ) *
cos( phi );
415 matrix1[0][2] = -r *
sin( theta ) *
sin( phi );
416 matrix1[1][0] =
sin( theta ) *
sin( phi );
417 matrix1[1][1] = r *
cos( theta ) *
sin( phi );
418 matrix1[1][2] = r *
sin( theta ) *
cos( phi );
419 matrix1[2][0] =
cos( theta );
420 matrix1[2][1] = -r *
sin( theta );
421 matrix1[2][2] = 0;
422
423
424 for ( int i = 0; i < 3; i++ )
425 {
426 for ( int j = 0; j < 3; j++ ) { matrix2[i][j] = matrix1[j][i]; }
427 }
428
429
430 HepMatrix matrix4( 3, 3 );
431 matrix4 = matrix1 * matrix * matrix2;
432
433
434 HepSymMatrix matrix5( 3 );
435 for ( int i = 0; i < 3; i++ )
436 {
437 for ( int j = 0; j <= i; j++ ) { matrix5[i][j] = matrix4[i][j]; }
438 }
439
440
442
443 if ( matrix5[0][0] > 0 ) dx = sqrt( matrix5[0][0] );
444 if ( matrix5[1][1] > 0 ) dy = sqrt( matrix5[1][1] );
445 if ( matrix5[2][2] > 0 ) dz = sqrt( matrix5[2][2] );
446}
HepGeom::Point3D< double > HepPoint3D
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
double tan(const BesAngle a)
double sin(const BesAngle a)
double cos(const BesAngle a)
void setDtheta(double dt)
void setPosition(const HepPoint3D &pos)
void setErrorMatrix(const HepSymMatrix &error)
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
static unsigned int theta_module(const Identifier &id)
static unsigned int phi_module(const Identifier &id)
double LogPosOffset() const
double BarrPosMCCor(int ntheta, int nphi) const
static EmcRecParameter & GetInstance()
double EastPosDataCor(int ntheta, int nphi) const
double WestPosDataCor(int ntheta, int nphi) const
double BarrLogPhiPara(int n, int m) const
double MethodMode() const
std::string PositionMode2() const
double SigTheta(int n) const
double WestLogThetaPara(int n, int m) const
double BarrPosDataCor(int ntheta, int nphi) const
double BarrLogThetaPara(int n, int m) const
double WestLogPhiPara(int n, int m) const
double EastLogPhiPara(int n, int m) const
double EastLogThetaPara(int n, int m) const
double SigPhi(int n) const
virtual double GetBarrelR() const =0
virtual HepPoint3D GetCFrontCenter(const Identifier &id) const =0
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
RecEmcFractionMap::const_iterator End() const
RecEmcFractionMap getFractionMap5x5() const
RecEmcID getShowerId() const
RecEmcFractionMap::const_iterator Begin() const
RecEmcEnergy getEAll() const
RecEmcFractionMap getFractionMap3x3() const