BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcRecShowerPosLin Class Reference

#include <EmcRecShowerPosLin.h>

Inheritance diagram for EmcRecShowerPosLin:

Public Member Functions

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

Detailed Description

Definition at line 15 of file InstallArea/x86_64-el9-gcc13-dbg/include/EmcRec/EmcRecShowerPosLin.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLin() [1/3]

EmcRecShowerPosLin::EmcRecShowerPosLin ( )
inline

Definition at line 17 of file InstallArea/x86_64-el9-gcc13-dbg/include/EmcRec/EmcRecShowerPosLin.h.

17{}

◆ ~EmcRecShowerPosLin() [1/3]

virtual EmcRecShowerPosLin::~EmcRecShowerPosLin ( )
inlinevirtual

Definition at line 18 of file InstallArea/x86_64-el9-gcc13-dbg/include/EmcRec/EmcRecShowerPosLin.h.

18{}

◆ EmcRecShowerPosLin() [2/3]

EmcRecShowerPosLin::EmcRecShowerPosLin ( )
inline

Definition at line 17 of file InstallArea/x86_64-el9-gcc13-opt/include/EmcRec/EmcRecShowerPosLin.h.

17{}

◆ ~EmcRecShowerPosLin() [2/3]

virtual EmcRecShowerPosLin::~EmcRecShowerPosLin ( )
inlinevirtual

Definition at line 18 of file InstallArea/x86_64-el9-gcc13-opt/include/EmcRec/EmcRecShowerPosLin.h.

18{}

◆ EmcRecShowerPosLin() [3/3]

EmcRecShowerPosLin::EmcRecShowerPosLin ( )
inline

Definition at line 17 of file Reconstruction/EmcRec/include/EmcRec/EmcRecShowerPosLin.h.

17{}

◆ ~EmcRecShowerPosLin() [3/3]

virtual EmcRecShowerPosLin::~EmcRecShowerPosLin ( )
inlinevirtual

Definition at line 18 of file Reconstruction/EmcRec/include/EmcRec/EmcRecShowerPosLin.h.

18{}

Member Function Documentation

◆ Position() [1/3]

void EmcRecShowerPosLin::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 19 of file EmcRecShowerPosLin.cxx.

19 {
20 RecEmcFractionMap::const_iterator cit;
21 HepPoint3D pos( 0, 0, 0 );
22 HepPoint3D possum( 0, 0, 0 );
23 double w, wsum;
24 double etot;
25 // cout<<"EmcRecShowerPosLin::Position Begin"<<endl;
26
27 IEmcRecGeoSvc* iGeoSvc;
28 ISvcLocator* svcLocator = Gaudi::svcLocator();
29 StatusCode sc = svcLocator->service( "EmcRecGeoSvc", iGeoSvc );
30 if ( sc != StatusCode::SUCCESS )
31 {
32 // cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
33 }
34
35 EmcRecParameter& Para = EmcRecParameter::GetInstance();
36
37 double e5x5 = aShower.e5x5();
38
39 etot = 0;
40
41 for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
42 { etot += ( cit->second.getEnergy() * cit->second.getFraction() ); }
43 wsum = 0;
44 for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
45 {
46 w = ( cit->second.getEnergy() * cit->second.getFraction() );
47 wsum += w;
48 pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
49 possum += pos * w;
50 }
51
52 if ( wsum <= 0 ) { cout << "[[[[[[[[[[[[[[[" << wsum; }
53
54 pos = possum / wsum;
55 // aShower.setPosition(pos);
56
57 // PosCorr=0 without position correction
58 // PosCorr=1 with position correction
59 if ( Para.PosCorr() == 0 ) { aShower.setPosition( pos ); }
60
61 //----------------------------------------------------------------------
62 // position correction
63 RecEmcID id = aShower.getShowerId();
64 const unsigned int module = EmcID::barrel_ec( id );
65 ;
66 const unsigned int thetaModule = EmcID::theta_module( id );
67 const unsigned int phiModule = EmcID::phi_module( id );
68 HepPoint3D posCorr( pos );
69
70 if ( module == 1 )
71 { // barrel
72 unsigned int theModule;
73 if ( thetaModule > 21 )
74 {
75 theModule = 43 - thetaModule;
76 id = EmcID::crystal_id( module, theModule, phiModule );
77 posCorr.setZ( -posCorr.z() );
78 }
79 else { theModule = thetaModule; }
80 // cout<<"EmcID: "<<id<<" theta: "<<theModule<<endl;
81
82 EmcRecParameter& Para = EmcRecParameter::GetInstance();
83 double IRShift, parTheta[5], parPhi[5];
84
85 if ( theModule == 21 ) { IRShift = 0; }
86 else if ( theModule == 20 ) { IRShift = 2.5; }
87 else { IRShift = 5.; }
88
89 if ( Para.MethodMode() == 0 )
90 {
91
92 for ( int i = 0; i < 5; i++ )
93 {
94 if ( theModule == 21 )
95 {
96 double par[5] = { 1.698, -1.553, 0.9384, 0.1193, -0.01916 };
97 parTheta[i] = par[i];
98 }
99 else if ( theModule == 20 )
100 {
101 double par[5] = { 1.593, -1.582, 0.8881, 0.3298, -0.08856 };
102 parTheta[i] = par[i];
103 }
104 else
105 {
106 double par[5] = { 1.605, -1.702, 0.8433, 0.3072, -0.1809 };
107 parTheta[i] = par[i];
108 }
109 }
110 }
111 else if ( Para.MethodMode() == 1 )
112 {
113
114 for ( int i = 0; i < 5; i++ )
115 {
116
117 if ( e5x5 > 1.0 ) parTheta[i] = Para.BarrLinThetaPara( theModule, i );
118 else if ( e5x5 <= 1.0 && e5x5 > 0.5 )
119 parTheta[i] = Para.BarrLinThetaPara( theModule + 22, i );
120 else if ( e5x5 <= 0.5 ) parTheta[i] = Para.BarrLinThetaPara( theModule + 44, i );
121
122 if ( e5x5 > 1.0 ) parPhi[i] = Para.BarrLinPhiPara( 0, i );
123 else if ( e5x5 <= 1.0 && e5x5 > 0.5 ) parPhi[i] = Para.BarrLinPhiPara( 1, i );
124 else if ( e5x5 <= 0.5 ) parPhi[i] = Para.BarrLinPhiPara( 2, i );
125 }
126 }
127
128 HepPoint3D IR( 0, 0, IRShift );
129 HepPoint3D center01, center23; // center of point0,1 and point2,3 in crystal
130 center01 = ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
131 center23 = ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
132 // cout<<"before correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
133
134 double theta01, theta23, length, d;
135 theta01 = ( center01 - IR ).theta();
136 theta23 = ( center23 - IR ).theta();
137 length = ( center01 - IR ).mag();
138 d = length * tan( theta01 - theta23 ); // length in x direction
139 // cout<<"theta01: "<<theta01<<"\t"<<" theta23: "<<theta23<<"\t"
140 //<<" length: "<<length<<" d: "<<d<<endl;
141
142 double xIn, xOut, deltaTheta;
143 xIn = length * tan( theta01 - ( posCorr - IR ).theta() ) - d / 2;
144 xOut = xIn - ( parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
145 parTheta[3] );
146 deltaTheta = atan( ( xOut + d / 2 ) / length ) - atan( ( xIn + d / 2 ) / length );
147 // cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
148
149 double parPhi1[5] = { 0.7185, -2.539, 0.6759, 0.3472, -0.3108 }; // e-
150 double parPhi2[5] = { 0.5781, -2.917, 0.5516, -0.5745, 0.3777 }; // e+
151
152 // double parPhi[5];// = { 0.7723, -3.1, 0.6992, -0.1441, -0.01012 }; //by photon, not
153 // used
154
155 // for(int i=0;i<5;i++) {
156 // parPhi[i] = (parPhi1[i]+parPhi2[i])/2; //average of e- and e+
157 // }
158
159 double yIn, yOut, deltaPhi;
160 // yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
161 // : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
162 if ( phiModule == 0 && posCorr.phi() < 0 )
163 { yIn = posCorr.phi() * 180. / CLHEP::pi + 360. - 119 * 3 - 1.843; }
164 else if ( phiModule == 119 && posCorr.phi() > 0 )
165 { yIn = posCorr.phi() * 180. / CLHEP::pi - 1.843; }
166 else
167 {
168
169 yIn = posCorr.phi() < 0
170 ? posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
171 : posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
172 }
173
174 yOut = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
175 deltaPhi = yOut * CLHEP::pi / 180.;
176 // cout<<"yIn="<<yIn<<"\tyOut="<<yOut<<"\tdeltaPhi="<<deltaPhi<<endl;
177
178 HepPoint3D rotateVector( 0, 1, 0 ); // y axis
179 rotateVector.rotateZ( posCorr.phi() );
180 posCorr.setZ( posCorr.z() - IRShift );
181 posCorr.rotate( -deltaTheta, rotateVector );
182 posCorr.setZ( posCorr.z() + IRShift );
183 if ( Para.MethodMode() == 0 ) { posCorr.rotateZ( -0.002994 ); }
184 else if ( Para.MethodMode() == 1 ) { posCorr.rotateZ( -deltaPhi ); }
185
186 if ( thetaModule > 21 ) { posCorr.setZ( -posCorr.z() ); }
187 }
188 else
189 { // endcap
190
191 EmcRecParameter& Para = EmcRecParameter::GetInstance();
192
193 if ( Para.MethodMode() == 1 )
194 {
195
196 double IRShift = 0;
197 HepPoint3D IR( 0, 0, IRShift );
198 double theta12 = -9999, theta03 = -9999, theta23 = -9999, theta14 = -9999, delta = -9999;
199 double phi = -9999, phi01 = -9999, phi23 = -9999, phi12 = -9999, phi34 = -9999,
200 deltaphi = -9999;
201 double xIn, xOut, deltaTheta, yIn, yOut, deltaPhi;
202 double posphi = posCorr.phi();
203
204 double parTheta[5], parPhi[5];
205
206 HepPoint3D point0, point1, point2, point3, point4;
207 point0 = iGeoSvc->GetCrystalPoint( id, 0 );
208 point1 = iGeoSvc->GetCrystalPoint( id, 1 );
209 point2 = iGeoSvc->GetCrystalPoint( id, 2 );
210 point3 = iGeoSvc->GetCrystalPoint( id, 3 );
211 point4 = iGeoSvc->GetCrystalPoint( id, 4 );
212
213 HepPoint3D center, center01, center23, center12, center03, center34, center14;
214 center01 = ( point0 + point1 ) / 2;
215 center23 = ( point2 + point3 ) / 2;
216 center12 = ( point1 + point2 ) / 2;
217 center34 = ( point3 + point4 ) / 2;
218 center03 = ( point0 + point3 ) / 2;
219 center14 = ( point1 + point4 ) / 2;
220
221 phi01 = center01.phi();
222 phi23 = center23.phi();
223 phi12 = center12.phi();
224 phi34 = center34.phi();
225
226 if ( ( thetaModule == 1 &&
227 ( ( phiModule + 3 ) % 4 == 0 || ( phiModule + 2 ) % 4 == 0 ) ) ||
228 ( thetaModule == 3 && ( ( phiModule + 4 ) % 5 == 0 || ( phiModule + 3 ) % 5 == 0 ||
229 ( phiModule + 2 ) % 5 == 0 ) ) )
230 {
231 center12 = center23;
232 center03 = center14;
233
234 center01 = center12;
235 center23 = center34;
236
237 phi01 = phi12;
238 phi23 = phi34;
239 }
240
241 if ( phiModule == 0 ) { posphi = posphi; }
242 else if ( ( ( thetaModule == 0 || thetaModule == 1 ) && phiModule == 63 ) ||
243 ( ( thetaModule == 2 || thetaModule == 3 ) && phiModule == 79 ) ||
244 ( ( thetaModule == 4 || thetaModule == 5 ) && phiModule == 95 ) )
245 { posphi = posphi + 2 * CLHEP::pi; }
246 else { posphi = posphi < 0 ? posphi + 2 * CLHEP::pi : posphi; }
247
248 phi01 = phi01 < 0 ? phi01 + 2 * CLHEP::pi : phi01;
249 phi23 = phi23 < 0 ? phi23 + 2 * CLHEP::pi : phi23;
250
251 if ( module == 0 )
252 { // east endcap
253
254 IRShift = 10;
255 center = center23;
256 phi = phi23;
257
258 for ( int i = 0; i < 5; i++ )
259 {
260 parTheta[i] = Para.EastLinThetaPara( thetaModule, i );
261 parPhi[i] = Para.EastLinPhiPara( 0, i );
262 }
263 }
264 else if ( module == 2 )
265 { // west endcap
266
267 IRShift = -10;
268 center = center01;
269 phi = phi01;
270
271 for ( int i = 0; i < 5; i++ )
272 {
273 parTheta[i] = Para.WestLinThetaPara( thetaModule, i );
274 parPhi[i] = Para.WestLinPhiPara( 0, i );
275 }
276 }
277
278 theta12 = ( center12 - IR ).theta();
279 theta03 = ( center03 - IR ).theta();
280 delta = theta03 - theta12;
281 xIn = ( ( ( posCorr - IR ).theta() - theta12 ) - delta * 0.5 ) / delta;
282 xOut = parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
283 parTheta[3];
284 deltaTheta = xOut * delta;
285
286 deltaphi = fabs( phi23 - phi01 );
287 yIn = ( ( posphi - phi ) - deltaphi * 0.5 ) / deltaphi;
288 yOut = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
289 deltaPhi = yOut * deltaphi;
290
291 HepPoint3D rotateVector( 0, 1, 0 ); // y axis
292 rotateVector.rotateZ( center.phi() );
293 posCorr.setZ( posCorr.z() - IRShift );
294 posCorr.rotate( -deltaTheta, rotateVector );
295 posCorr.setZ( posCorr.z() + IRShift );
296 posCorr.rotateZ( -deltaPhi );
297 }
298 }
299
300 // aShower.setPosition(posCorr);
301 if ( Para.PosCorr() == 1 ) { aShower.setPosition( posCorr ); }
302
303 //----------------------------------------------------------------------
304 if ( aShower.energy() < 1e-5 ) return;
305
306 double r, theta, phi;
307 double dtheta, dphi, dx, dy, dz;
308
309 r = posCorr.mag();
310 theta = posCorr.theta();
311 phi = posCorr.phi();
312 // cout<<"theta: "<<theta<<" phi: "<<phi<<endl;
313
314 // EmcRecParameter& Para=EmcRecParameter::GetInstance();
315 if ( aShower.energy() > 0 )
316 {
317 dtheta = Para.SigTheta( 0 ) / sqrt( aShower.energy() ) + Para.SigTheta( 1 );
318 dphi = Para.SigPhi( 0 ) / sqrt( aShower.energy() ) + Para.SigPhi( 1 );
319 }
320 else
321 {
322 dtheta = 0;
323 dphi = 0;
324 }
325 // cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
326
327 dx = fabs( iGeoSvc->GetBarrelR() * sin( phi ) * dphi );
328 dy = fabs( iGeoSvc->GetBarrelR() * cos( phi ) * dphi );
329 if ( theta > 0 )
330 dz = fabs( iGeoSvc->GetBarrelR() * dtheta / ( sin( theta ) * sin( theta ) ) );
331 else dz = 0;
332 // cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
333
334 aShower.setDtheta( dtheta );
335 aShower.setDphi( dphi );
336
337 //----------------------------------------------------------------------
338 // Error matrix
339 HepSymMatrix matrix( 3 ); // error matrix of r, theta, phi
340 matrix[0][0] = 0;
341 matrix[1][1] = dtheta * dtheta;
342 matrix[2][2] = dphi * dphi;
343 matrix[0][1] = 0;
344 matrix[0][2] = 0;
345 matrix[1][2] = 0;
346 // cout<<"matrix: \n"<<matrix;
347
348 HepMatrix matrix1( 3, 3 ), matrix2( 3, 3 );
349 matrix1[0][0] = sin( theta ) * cos( phi );
350 matrix1[0][1] = r * cos( theta ) * cos( phi );
351 matrix1[0][2] = -r * sin( theta ) * sin( phi );
352 matrix1[1][0] = sin( theta ) * sin( phi );
353 matrix1[1][1] = r * cos( theta ) * sin( phi );
354 matrix1[1][2] = r * sin( theta ) * cos( phi );
355 matrix1[2][0] = cos( theta );
356 matrix1[2][1] = -r * sin( theta );
357 matrix1[2][2] = 0;
358 // cout<<"matrix1: \n"<<matrix1;
359
360 for ( int i = 0; i < 3; i++ )
361 {
362 for ( int j = 0; j < 3; j++ ) { matrix2[i][j] = matrix1[j][i]; }
363 }
364 // cout<<"matrix2: \n"<<matrix2;
365
366 HepMatrix matrix4( 3, 3 );
367 matrix4 = matrix1 * matrix * matrix2;
368 // cout<<"matrix4: \n"<<matrix4;
369
370 HepSymMatrix matrix5( 3 ); // error matrix of x, y, z
371 for ( int i = 0; i < 3; i++ )
372 {
373 for ( int j = 0; j <= i; j++ ) { matrix5[i][j] = matrix4[i][j]; }
374 }
375 // cout<<"matrix5: \n"<<matrix5;
376
377 aShower.setErrorMatrix( matrix5 );
378}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
etot
Double_t etot
Definition DataBase/tau_mode.c:8
w
double w
Definition EvtPhokharaDef.hh:28
tan
double tan(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:189
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
DstEmcShower::setPosition
void setPosition(const HepPoint3D &pos)
Definition Event/DstEvent/include/DstEvent/DstEmcShower.h:62
DstEmcShower::setErrorMatrix
void setErrorMatrix(const HepSymMatrix &error)
Definition Event/DstEvent/include/DstEvent/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:63
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition EmcID.cxx:36
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:41
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:46
EmcRecParameter::GetInstance
static EmcRecParameter & GetInstance()
Definition EmcRecParameter.cxx:494
EmcRecParameter::EastLinThetaPara
double EastLinThetaPara(int n, int m) const
Definition EmcRecParameter.cxx:617
EmcRecParameter::EastLinPhiPara
double EastLinPhiPara(int n, int m) const
Definition EmcRecParameter.cxx:613
EmcRecParameter::PosCorr
double PosCorr() const
Definition EmcRecParameter.cxx:524
EmcRecParameter::MethodMode
double MethodMode() const
Definition EmcRecParameter.cxx:522
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition EmcRecParameter.cxx:537
EmcRecParameter::WestLinThetaPara
double WestLinThetaPara(int n, int m) const
Definition EmcRecParameter.cxx:621
EmcRecParameter::BarrLinPhiPara
double BarrLinPhiPara(int n, int m) const
Definition EmcRecParameter.cxx:567
EmcRecParameter::BarrLinThetaPara
double BarrLinThetaPara(int n, int m) const
Definition EmcRecParameter.cxx:563
EmcRecParameter::WestLinPhiPara
double WestLinPhiPara(int n, int m) const
Definition EmcRecParameter.cxx:615
EmcRecParameter::SigPhi
double SigPhi(int n) const
Definition EmcRecParameter.cxx:539
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:114
RecEmcShower::Begin
RecEmcFractionMap::const_iterator Begin() const
Definition RecEmcShower.cxx:112

◆ Position() [2/3]

virtual void EmcRecShowerPosLin::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

◆ Position() [3/3]

virtual void EmcRecShowerPosLin::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

The documentation for this class was generated from the following files: