CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
EmcRecShowerPosLin Class Reference

#include <EmcRecShowerPosLin.h>

+ Inheritance diagram for EmcRecShowerPosLin:

Public Member Functions

 EmcRecShowerPosLin ()
 
virtual ~EmcRecShowerPosLin ()
 
virtual void Position (RecEmcShower &aShower)
 
- Public Member Functions inherited from EmcRecShowerPosAbs
 EmcRecShowerPosAbs ()
 
virtual ~EmcRecShowerPosAbs ()
 

Detailed Description

Definition at line 15 of file EmcRecShowerPosLin.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLin()

EmcRecShowerPosLin::EmcRecShowerPosLin ( )
inline

Definition at line 18 of file EmcRecShowerPosLin.h.

18{}

◆ ~EmcRecShowerPosLin()

virtual EmcRecShowerPosLin::~EmcRecShowerPosLin ( )
inlinevirtual

Definition at line 19 of file EmcRecShowerPosLin.h.

19{}

Member Function Documentation

◆ Position()

void EmcRecShowerPosLin::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 17 of file EmcRecShowerPosLin.cxx.

18{
19 RecEmcFractionMap::const_iterator cit;
20 HepPoint3D pos(0,0,0);
21 HepPoint3D possum(0,0,0);
22 double w,wsum;
23 double etot;
24 //cout<<"EmcRecShowerPosLin::Position Begin"<<endl;
25
26 IEmcRecGeoSvc* iGeoSvc;
27 ISvcLocator* svcLocator = Gaudi::svcLocator();
28 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
29 if(sc!=StatusCode::SUCCESS) {
30 //cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
31 }
32 etot=0;
33 //cout<<etot<<endl;
34 for(cit=aShower.Begin();cit!=aShower.End();cit++){
35 etot+=(cit->second.getEnergy()*cit->second.getFraction());
36 //cout<<etot<<endl;
37 }
38 wsum=0;
39 for(cit=aShower.Begin();cit!=aShower.End();cit++){
40 w=(cit->second.getEnergy()*cit->second.getFraction());
41 //cout<<w<<endl;
42 wsum+=w;
43 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
44 possum+=pos*w;
45 }
46
47 if(wsum<=0) {
48 cout<<"[[[[[[[[[[[[[[["<<wsum;
49 }
50
51 pos=possum/wsum;
52 //aShower.setPosition(pos);
53
54 //----------------------------------------------------------------------
55 //position correction
56 RecEmcID id = aShower.getShowerId();
57 const unsigned int module = EmcID::barrel_ec(id);;
58 const unsigned int thetaModule = EmcID::theta_module(id);
59 const unsigned int phiModule = EmcID::phi_module(id);
60 HepPoint3D posCorr(pos);
61
62 if(module==1) { //barrel
63 unsigned int theModule;
64 if(thetaModule>21) {
65 theModule = 43 - thetaModule;
66 id = EmcID::crystal_id(module,theModule,phiModule);
67 posCorr.setZ(-posCorr.z());
68 } else {
69 theModule = thetaModule;
70 }
71 //cout<<"EmcID: "<<id<<" theta: "<<theModule<<endl;
72
73 double IRShift, parTheta[5];
74 if(theModule==21) {
75 IRShift = 0;
76 double par[5] = { 1.698, -1.553, 0.9384, 0.1193, -0.01916 };
77 for(int i=0;i<5;i++) {
78 parTheta[i] = par[i];
79 }
80 } else if(theModule==20) {
81 IRShift = 2.5;
82 double par[5] = { 1.593, -1.582, 0.8881, 0.3298, -0.08856 };
83 for(int i=0;i<5;i++) {
84 parTheta[i] = par[i];
85 }
86 } else {
87 IRShift = 5.;
88 double par[5] = { 1.605, -1.702, 0.8433, 0.3072, -0.1809 };
89 for(int i=0;i<5;i++) {
90 parTheta[i] = par[i];
91 }
92 }
93
94 HepPoint3D IR(0,0,IRShift);
95 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
96 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
97 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
98 //cout<<"before correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
99
100 double theta01,theta23,length,d;
101 theta01 = (center01-IR).theta();
102 theta23 = (center23-IR).theta();
103 length = (center01-IR).mag();
104 d = length*tan(theta01-theta23); //length in x direction
105 //cout<<"theta01: "<<theta01<<"\t"<<" theta23: "<<theta23<<"\t"
106 //<<" length: "<<length<<" d: "<<d<<endl;
107
108 double xIn,xOut,deltaTheta;
109 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
110 xOut = xIn-(parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+parTheta[2]*xIn+parTheta[3]);
111 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
112 //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
113
114 double parPhi1[5] = { 0.7185, -2.539, 0.6759, 0.3472, -0.3108 }; //e-
115 double parPhi2[5] = { 0.5781, -2.917, 0.5516, -0.5745, 0.3777 }; //e+
116 double parPhi[5];// = { 0.7723, -3.1, 0.6992, -0.1441, -0.01012 }; //by photon, not used
117 for(int i=0;i<5;i++) {
118 parPhi[i] = (parPhi1[i]+parPhi2[i])/2; //average of e- and e+
119 }
120
121 double yIn,yOut,deltaPhi;
122 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
123 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
124 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
125 deltaPhi = yOut*CLHEP::pi/180.;
126 //cout<<"yIn="<<yIn<<"\tyOut="<<yOut<<"\tdeltaPhi="<<deltaPhi<<endl;
127
128 HepPoint3D rotateVector(0,1,0); //y axis
129 rotateVector.rotateZ(posCorr.phi());
130 posCorr.setZ(posCorr.z()-IRShift);
131 posCorr.rotate(-deltaTheta,rotateVector);
132 posCorr.setZ(posCorr.z()+IRShift);
133 //posCorr.rotateZ(-deltaPhi);
134 posCorr.rotateZ(-0.002994);
135 //cout<<"after correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
136
137 if(thetaModule>21) {
138 posCorr.setZ(-posCorr.z());
139 }
140 }
141
142 aShower.setPosition(posCorr);
143
144 //----------------------------------------------------------------------
145 if(aShower.energy()<1e-5) return;
146
147 double r,theta,phi;
148 double dtheta,dphi,dx,dy,dz;
149
150 r = posCorr.mag();
151 theta = posCorr.theta();
152 phi = posCorr.phi();
153 //cout<<"theta: "<<theta<<" phi: "<<phi<<endl;
154
155 EmcRecParameter& Para=EmcRecParameter::GetInstance();
156 if(aShower.energy()>0) {
157 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
158 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
159 }
160 else {
161 dtheta = 0;
162 dphi = 0;
163 }
164 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
165
166 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
167 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
168 if(theta>0)
169 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
170 else
171 dz = 0;
172 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
173
174 aShower.setDtheta(dtheta);
175 aShower.setDphi(dphi);
176
177 //----------------------------------------------------------------------
178 // Error matrix
179 HepSymMatrix matrix(3); //error matrix of r, theta, phi
180 matrix[0][0]=0;
181 matrix[1][1]=dtheta*dtheta;
182 matrix[2][2]=dphi*dphi;
183 matrix[0][1]=0;
184 matrix[0][2]=0;
185 matrix[1][2]=0;
186 //cout<<"matrix: \n"<<matrix;
187
188 HepMatrix matrix1(3,3),matrix2(3,3);
189 matrix1[0][0]=sin(theta)*cos(phi);
190 matrix1[0][1]=r*cos(theta)*cos(phi);
191 matrix1[0][2]=-r*sin(theta)*sin(phi);
192 matrix1[1][0]=sin(theta)*sin(phi);
193 matrix1[1][1]=r*cos(theta)*sin(phi);
194 matrix1[1][2]=r*sin(theta)*cos(phi);
195 matrix1[2][0]=cos(theta);
196 matrix1[2][1]=-r*sin(theta);
197 matrix1[2][2]=0;
198 //cout<<"matrix1: \n"<<matrix1;
199
200 for(int i=0;i<3;i++) {
201 for(int j=0;j<3;j++) {
202 matrix2[i][j]=matrix1[j][i];
203 }
204 }
205 //cout<<"matrix2: \n"<<matrix2;
206
207 HepMatrix matrix4(3,3);
208 matrix4=matrix1*matrix*matrix2;
209 //cout<<"matrix4: \n"<<matrix4;
210
211 HepSymMatrix matrix5(3); //error matrix of x, y, z
212 for(int i=0;i<3;i++) {
213 for(int j=0;j<=i;j++) {
214 matrix5[i][j]=matrix4[i][j];
215 }
216 }
217 //cout<<"matrix5: \n"<<matrix5;
218
219 aShower.setErrorMatrix(matrix5);
220}
tan
double tan(const BesAngle a)
Definition BesAngle.h:216
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
length
double length
Definition CgemLineFit.cxx:121
etot
Double_t etot
Definition DataBase/tau_mode.c:7
DstEmcShower::setDphi
void setDphi(double dpi)
Definition DstEmcShower.h:66
DstEmcShower::setDtheta
void setDtheta(double dt)
Definition DstEmcShower.h:65
DstEmcShower::setPosition
void setPosition(const HepPoint3D &pos)
Definition DstEmcShower.h:62
DstEmcShower::energy
double energy() const
Definition DstEmcShower.h:45
DstEmcShower::setErrorMatrix
void setErrorMatrix(const HepSymMatrix &error)
Definition DstEmcShower.h:75
EmcID::crystal_id
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
Definition EmcID.cxx:71
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition EmcID.cxx:38
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:43
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:48
EmcRecParameter::GetInstance
static EmcRecParameter & GetInstance()
Definition EmcRecParameter.cxx:115
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition EmcRecParameter.cxx:187
EmcRecParameter::SigPhi
double SigPhi(int n) const
Definition EmcRecParameter.cxx:192
IEmcRecGeoSvc::GetBarrelR
virtual double GetBarrelR() const =0
IEmcRecGeoSvc::GetCFrontCenter
virtual HepPoint3D GetCFrontCenter(const Identifier &id) const =0
IEmcRecGeoSvc::GetCrystalPoint
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
RecEmcShower::End
RecEmcFractionMap::const_iterator End() const
Definition RecEmcShower.cxx:136
RecEmcShower::getShowerId
RecEmcID getShowerId() const
Definition RecEmcShower.h:55
RecEmcShower::Begin
RecEmcFractionMap::const_iterator Begin() const
Definition RecEmcShower.cxx:131
The documentation for this class was generated from the following files: