20 RecEmcFractionMap::const_iterator cit;
28 ISvcLocator* svcLocator = Gaudi::svcLocator();
29 StatusCode sc = svcLocator->service(
"EmcRecGeoSvc", iGeoSvc );
30 if ( sc != StatusCode::SUCCESS )
45 for ( cit = aShower.
Begin(); cit != aShower.
End(); cit++ )
47 w = offset + log( ( cit->second.getEnergy() /
etot ) * cit->second.getFraction() );
58 double e3x3 = aShower.
e3x3();
60 for ( cit = fracMap3x3.begin(); cit != fracMap3x3.end(); cit++ )
62 w = offset + log( ( cit->second.getEnergy() / e3x3 ) * cit->second.getFraction() );
73 e5x5 = aShower.
e5x5();
75 for ( cit = fracMap5x5.begin(); cit != fracMap5x5.end(); cit++ )
77 w = offset + log( ( cit->second.getEnergy() / e5x5 ) * cit->second.getFraction() );
85 Rc = iGeoSvc->
GetCCenter( cit->second.getCellId() );
87 pos = R + ( Rc - R ) / ( Rc - R ).mag() * NX0 * 1.86;
100 HepPoint3D posFront = pos - pos / pos.mag() * NX0 * 1.86;
116 double IRShift, parTheta[5], parPhi[5];
118 for (
int i = 0; i < 5; i++ )
122 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
127 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
139 double theta01, theta23, dtheta;
141 theta01 = ( center01 ).theta();
142 theta23 = ( center23 ).theta();
144 dtheta = theta01 - theta23;
145 double xIn, xOut, deltaTheta, deltaX;
146 xIn = ( ( ( posCorr ).theta() - theta23 ) - dtheta / 2 ) / dtheta;
148 deltaX = parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
150 deltaTheta = deltaX * dtheta;
155 double phi12 = -9999, phi03 = -9999, deltaphi = -9999;
157 phi03 = center03.phi();
158 phi12 = center12.phi();
161 posPhi = posCorr.phi();
163 if ( phiModule == 0 )
169 else if ( phiModule == 119 )
171 posPhi = posPhi + 2 * CLHEP::pi;
172 phi03 = phi03 + 2 * CLHEP::pi;
173 phi12 = phi12 + 2 * CLHEP::pi;
177 posPhi = posPhi < 0 ? posPhi + 2 * CLHEP::pi : posPhi;
178 phi03 = phi03 < 0 ? phi03 + 2 * CLHEP::pi : phi03;
179 phi12 = phi12 < 0 ? phi12 + 2 * CLHEP::pi : phi12;
182 deltaphi = phi12 - phi03;
184 double yIn, deltaY, deltaPhi;
185 yIn = ( ( posPhi - phi03 ) - deltaphi * 0.5 ) / deltaphi;
187 deltaY = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
188 deltaPhi = deltaY * deltaphi;
208 rotateVector.rotateZ( center23.phi() );
210 posCorr.rotate( -deltaTheta, rotateVector );
211 posCorr.rotateZ( -deltaPhi );
235 double IRShift = 0, parTheta[5], parPhi[5];
242 for (
int i = 0; i < 5; i++ )
246 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
255 else if ( module == 2 )
258 for (
int i = 0; i < 5; i++ )
262 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
274 double theta12 = -9999, theta03 = -9999, theta23 = -9999, theta14 = -9999, delta = -9999;
275 double phi01 = -9999, phi23 = -9999, phi12 = -9999, phi34 = -9999, deltaphi = -9999;
276 double xIn, deltaX, deltaTheta, yIn, deltaY, deltaPhi;
278 double posphi = posCorr.phi();
279 if ( phiModule == 0 ) { posphi = posphi; }
280 else if ( ( ( thetaModule == 0 || thetaModule == 1 ) && phiModule == 63 ) ||
281 ( ( thetaModule == 2 || thetaModule == 3 ) && phiModule == 79 ) ||
282 ( ( thetaModule == 4 || thetaModule == 5 ) && phiModule == 95 ) )
283 { posphi = posphi + 2 * CLHEP::pi; }
284 else { posphi = posphi < 0 ? posphi + 2 * CLHEP::pi : posphi; }
286 HepPoint3D center, center01, center23, center12, center03, center34, center14;
288 if ( ( thetaModule == 1 &&
289 ( ( phiModule + 3 ) % 4 == 0 || ( phiModule + 2 ) % 4 == 0 ) ) ||
290 ( thetaModule == 3 && ( ( phiModule + 4 ) % 5 == 0 || ( phiModule + 3 ) % 5 == 0 ||
291 ( phiModule + 2 ) % 5 == 0 ) ) )
299 theta23 = ( center23 - IR ).theta();
300 theta14 = ( center14 - IR ).theta();
301 delta = theta14 - theta23;
302 xIn = ( ( ( posCorr - IR ).theta() - theta23 ) - delta * 0.5 ) / delta;
308 phi12 = center12.phi();
309 phi34 = center34.phi();
311 if ( phiModule == 0 )
316 else if ( ( ( thetaModule == 0 || thetaModule == 1 ) && phiModule == 63 ) ||
317 ( ( thetaModule == 2 || thetaModule == 3 ) && phiModule == 79 ) ||
318 ( ( thetaModule == 4 || thetaModule == 5 ) && phiModule == 95 ) )
320 phi12 = phi12 + 2 * CLHEP::pi;
321 phi34 = phi34 + 2 * CLHEP::pi;
325 phi12 = phi12 < 0 ? phi12 + 2 * CLHEP::pi : phi12;
326 phi34 = phi34 < 0 ? phi34 + 2 * CLHEP::pi : phi34;
329 deltaphi = phi34 - phi12;
330 yIn = ( ( posphi - phi12 ) - deltaphi * 0.5 ) / deltaphi;
341 theta12 = ( center12 - IR ).theta();
342 theta03 = ( center03 - IR ).theta();
343 delta = theta03 - theta12;
344 xIn = ( ( ( posCorr - IR ).theta() - theta12 ) - delta * 0.5 ) / delta;
350 phi01 = center01.phi();
351 phi23 = center23.phi();
353 if ( phiModule == 0 )
358 else if ( ( ( thetaModule == 0 || thetaModule == 1 ) && phiModule == 63 ) ||
359 ( ( thetaModule == 2 || thetaModule == 3 ) && phiModule == 79 ) ||
360 ( ( thetaModule == 4 || thetaModule == 5 ) && phiModule == 95 ) )
362 phi01 = phi01 + 2 * CLHEP::pi;
363 phi23 = phi23 + 2 * CLHEP::pi;
367 phi01 = phi01 < 0 ? phi01 + 2 * CLHEP::pi : phi01;
368 phi23 = phi23 < 0 ? phi23 + 2 * CLHEP::pi : phi23;
371 deltaphi = phi23 - phi01;
372 yIn = ( ( posphi - phi01 ) - deltaphi * 0.5 ) / deltaphi;
375 deltaX = parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
377 deltaTheta = deltaX * delta;
379 deltaY = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
380 deltaPhi = deltaY * deltaphi;
383 rotateVector.rotateZ( center.phi() );
384 posCorr.setZ( posCorr.z() - IRShift );
385 posCorr.rotate( -deltaTheta, rotateVector );
386 posCorr.setZ( posCorr.z() + IRShift );
387 posCorr.rotateZ( -deltaPhi );
423 if ( aShower.
energy() < 1e-5 )
return;
425 double r, theta, phi;
426 double dtheta, dphi, dx, dy, dz;
429 theta = posCorr.theta();
433 if ( aShower.
energy() > 0 )
448 dz = fabs( iGeoSvc->
GetBarrelR() * dtheta / (
sin( theta ) *
sin( theta ) ) );
457 HepSymMatrix matrix( 3 );
459 matrix[1][1] = dtheta * dtheta;
460 matrix[2][2] = dphi * dphi;
466 HepMatrix matrix1( 3, 3 ), matrix2( 3, 3 );
467 matrix1[0][0] =
sin( theta ) *
cos( phi );
468 matrix1[0][1] = r *
cos( theta ) *
cos( phi );
469 matrix1[0][2] = -r *
sin( theta ) *
sin( phi );
470 matrix1[1][0] =
sin( theta ) *
sin( phi );
471 matrix1[1][1] = r *
cos( theta ) *
sin( phi );
472 matrix1[1][2] = r *
sin( theta ) *
cos( phi );
473 matrix1[2][0] =
cos( theta );
474 matrix1[2][1] = -r *
sin( theta );
478 for (
int i = 0; i < 3; i++ )
480 for (
int j = 0; j < 3; j++ ) { matrix2[i][j] = matrix1[j][i]; }
484 HepMatrix matrix4( 3, 3 );
485 matrix4 = matrix1 * matrix * matrix2;
488 HepSymMatrix matrix5( 3 );
489 for (
int i = 0; i < 3; i++ )
491 for (
int j = 0; j <= i; j++ ) { matrix5[i][j] = matrix4[i][j]; }
497 if ( matrix5[0][0] > 0 ) dx = sqrt( matrix5[0][0] );
498 if ( matrix5[1][1] > 0 ) dy = sqrt( matrix5[1][1] );
499 if ( matrix5[2][2] > 0 ) dz = sqrt( matrix5[2][2] );