6#include "EmcRec/EmcRecShowerPosLogShMax.h"
11#include "EmcRecGeoSvc/EmcRecGeoSvc.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/ISvcLocator.h"
15#include "EmcRec/EmcRecParameter.h"
20 RecEmcFractionMap::const_iterator cit;
28 ISvcLocator* svcLocator = Gaudi::svcLocator();
29 StatusCode sc = svcLocator->service(
"EmcRecGeoSvc",iGeoSvc);
30 if(sc!=StatusCode::SUCCESS) {
43 for(cit=aShower.
Begin();cit!=aShower.
End();cit++){
44 w=offset+log((cit->second.getEnergy()/
etot)*cit->second.getFraction());
52 double e3x3=aShower.
e3x3();
54 for(cit=fracMap3x3.begin();cit!=fracMap3x3.end();cit++) {
55 w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
65 for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
66 w=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
73 Rc=iGeoSvc->
GetCCenter(cit->second.getCellId());
75 pos = R + (Rc-R)/(Rc-R).mag()*NX0*1.86;
88 HepPoint3D posFront=pos-pos/pos.mag()*NX0*1.86;
105 double IRShift, parTheta[5],parPhi[5];
108 for(
int i=0;i<5;i++){
122 HepPoint3D center01, center23, center12,center03;
128 double theta01,theta23,dtheta;
130 theta01 = (center01).theta();
131 theta23 = (center23).theta();
133 dtheta = theta01-theta23;
134 double xIn,xOut,deltaTheta,deltaX;
135 xIn = (((posCorr).theta()-theta23)-dtheta/2)/dtheta;
137 deltaX = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
138 deltaTheta = deltaX*dtheta;
143 double phi12=-9999,phi03=-9999,deltaphi=-9999;
145 phi03 = center03.phi();
146 phi12 = center12.phi();
150 posPhi=posCorr.phi();
158 }
else if( phiModule==119){
159 posPhi = posPhi+2*CLHEP::pi;
160 phi03 = phi03+2*CLHEP::pi;
161 phi12 = phi12+2*CLHEP::pi;
164 posPhi = posPhi <0 ? posPhi+2*CLHEP::pi : posPhi;
165 phi03 = phi03 <0 ? phi03+2*CLHEP::pi : phi03;
166 phi12 = phi12 <0 ? phi12+2*CLHEP::pi : phi12;
169 deltaphi = phi12 - phi03;
172 double yIn,deltaY,deltaPhi;
173 yIn = ((posPhi - phi03) - deltaphi*0.5)/deltaphi;
175 deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
176 deltaPhi = deltaY*deltaphi;
196 rotateVector.rotateZ(center23.phi());
198 posCorr.rotate(-deltaTheta,rotateVector);
199 posCorr.rotateZ(-deltaPhi);
222 double IRShift=0, parTheta[5],parPhi[5];
228 for(
int i=0;i<5;i++) {
240 }
else if (module==2){
242 for(
int i=0;i<5;i++){
258 double theta12=-9999,theta03=-9999,theta23=-9999,theta14=-9999,
delta=-9999;
259 double phi01=-9999,phi23=-9999,phi12=-9999,phi34=-9999,deltaphi=-9999;
260 double xIn,deltaX,deltaTheta,yIn,deltaY,deltaPhi;
262 double posphi=posCorr.phi();
265 }
else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
266 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
267 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
269 posphi = posphi+2*CLHEP::pi;
271 posphi = posphi < 0 ? posphi+2*CLHEP::pi : posphi;
274 HepPoint3D center,center01, center23, center12,center03,center34,center14;
277 if( (thetaModule==1&&((phiModule+3)%4 == 0||(phiModule+2)%4 == 0))
278 ||(thetaModule==3&&((phiModule+4)%5 == 0||(phiModule+3)%5 == 0||(phiModule+2)%5 == 0)))
284 theta23 = (center23 - IR).theta();
285 theta14 = (center14 - IR).theta();
286 delta = theta14 - theta23;
287 xIn = (((posCorr-IR).theta() - theta23) -
delta*0.5)/
delta;
292 phi12 = center12.phi();
293 phi34 = center34.phi();
298 }
else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
299 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
300 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
302 phi12 = phi12+2*CLHEP::pi;
303 phi34 = phi34+2*CLHEP::pi;
305 phi12 = phi12 <0 ? phi12+2*CLHEP::pi : phi12;
306 phi34 = phi34 <0 ? phi34+2*CLHEP::pi : phi34;
309 deltaphi = phi34 - phi12;
310 yIn = ((posphi - phi12) - deltaphi*0.5)/deltaphi;
318 theta12 = (center12 - IR).theta();
319 theta03 = (center03 - IR).theta();
320 delta = theta03 - theta12;
321 xIn = (((posCorr-IR).theta() - theta12) -
delta*0.5)/
delta;
326 phi01 = center01.phi();
327 phi23 = center23.phi();
333 }
else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
334 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
335 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
337 phi01 = phi01+2*CLHEP::pi;
338 phi23 = phi23+2*CLHEP::pi;
340 phi01 = phi01 <0 ? phi01+2*CLHEP::pi : phi01;
341 phi23 = phi23 <0 ? phi23+2*CLHEP::pi : phi23;
344 deltaphi = phi23 - phi01;
345 yIn = ((posphi - phi01) - deltaphi*0.5)/deltaphi;
349 deltaX = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
350 deltaTheta = deltaX*
delta;
352 deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4]) + parPhi[2]*yIn + parPhi[3] ;
353 deltaPhi = deltaY*deltaphi;
356 rotateVector.rotateZ(center.phi());
357 posCorr.setZ(posCorr.z()-IRShift);
358 posCorr.rotate(-deltaTheta,rotateVector);
359 posCorr.setZ(posCorr.z()+IRShift);
360 posCorr.rotateZ(-deltaPhi);
397 if(aShower.
energy()<1e-5)
return;
400 double dtheta,dphi,dx,dy,dz;
403 theta = posCorr.theta();
430 HepSymMatrix matrix(3);
432 matrix[1][1]=dtheta*dtheta;
433 matrix[2][2]=dphi*dphi;
439 HepMatrix matrix1(3,3),matrix2(3,3);
440 matrix1[0][0]=
sin(theta)*
cos(phi);
441 matrix1[0][1]=r*
cos(theta)*
cos(phi);
442 matrix1[0][2]=-r*
sin(theta)*
sin(phi);
443 matrix1[1][0]=
sin(theta)*
sin(phi);
444 matrix1[1][1]=r*
cos(theta)*
sin(phi);
445 matrix1[1][2]=r*
sin(theta)*
cos(phi);
446 matrix1[2][0]=
cos(theta);
447 matrix1[2][1]=-r*
sin(theta);
451 for(
int i=0;i<3;i++) {
452 for(
int j=0;j<3;j++) {
453 matrix2[i][j]=matrix1[j][i];
458 HepMatrix matrix4(3,3);
459 matrix4=matrix1*matrix*matrix2;
462 HepSymMatrix matrix5(3);
463 for(
int i=0;i<3;i++) {
464 for(
int j=0;j<=i;j++) {
465 matrix5[i][j]=matrix4[i][j];
472 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
473 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
474 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
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 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 WestLogShMaxPhiPara(int n, int m) const
double BarrLogShMaxPhiPara(int n, int m) const
static EmcRecParameter & GetInstance()
double EastLogShMaxPhiPara(int n, int m) const
double EastLogShMaxThetaPara(int n, int m) const
std::string PositionMode2() const
double SigTheta(int n) const
double WestLogShMaxThetaPara(int n, int m) const
double SigPhi(int n) const
double BarrLogShMaxThetaPara(int n, int m) const
virtual void Position(RecEmcShower &aShower)
virtual double GetBarrelR() const =0
virtual HepPoint3D GetCFrontCenter(const Identifier &id) const =0
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
virtual HepPoint3D GetCCenter(const Identifier &id) const =0
RecEmcFractionMap::const_iterator End() const
RecEmcFractionMap getFractionMap5x5() const
RecEmcID getShowerId() const
RecEmcFractionMap::const_iterator Begin() const
RecEmcEnergy getEAll() const
RecEmcFractionMap getFractionMap3x3() const