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

#include <EmcRecShowerPosLoglin.h>

Inheritance diagram for EmcRecShowerPosLoglin:

Public Member Functions

 EmcRecShowerPosLoglin ()
virtual ~EmcRecShowerPosLoglin ()
virtual void Position (RecEmcShower &aShower)
 EmcRecShowerPosLoglin ()
virtual ~EmcRecShowerPosLoglin ()
virtual void Position (RecEmcShower &aShower)
 EmcRecShowerPosLoglin ()
virtual ~EmcRecShowerPosLoglin ()
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/EmcRecShowerPosLoglin.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLoglin() [1/3]

EmcRecShowerPosLoglin::EmcRecShowerPosLoglin ( )
inline

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

17{}

◆ ~EmcRecShowerPosLoglin() [1/3]

virtual EmcRecShowerPosLoglin::~EmcRecShowerPosLoglin ( )
inlinevirtual

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

18{}

◆ EmcRecShowerPosLoglin() [2/3]

EmcRecShowerPosLoglin::EmcRecShowerPosLoglin ( )
inline

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

17{}

◆ ~EmcRecShowerPosLoglin() [2/3]

virtual EmcRecShowerPosLoglin::~EmcRecShowerPosLoglin ( )
inlinevirtual

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

18{}

◆ EmcRecShowerPosLoglin() [3/3]

EmcRecShowerPosLoglin::EmcRecShowerPosLoglin ( )
inline

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

17{}

◆ ~EmcRecShowerPosLoglin() [3/3]

virtual EmcRecShowerPosLoglin::~EmcRecShowerPosLoglin ( )
inlinevirtual

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

18{}

Member Function Documentation

◆ Position() [1/3]

void EmcRecShowerPosLoglin::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 18 of file EmcRecShowerPosLoglin.cxx.

18 {
19 RecEmcFractionMap::const_iterator cit;
20 HepPoint3D pos( 0, 0, 0 );
21 HepPoint3D possum( 0, 0, 0 );
22 double w, w1, w2, wsum = 0;
23 double num, num0;
24
25 // cout<<"EmcRecShowerPosLoglin::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 offset=Para.LogPosOffset();
38 double a0 = 4;
39 double offset = 0.0;
40 // cout<<offset<<endl;
41 double e5x5 = -99999;
42 if ( Para.PositionMode2() == "all" )
43 {
44 double etot = aShower.getEAll();
45 for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
46 {
47 w = offset + log( ( cit->second.getEnergy() / etot ) * cit->second.getFraction() );
48 if ( w > 0 )
49 {
50 wsum += w;
51 pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
52 possum += pos * w;
53 }
54 }
55 }
56 else if ( Para.PositionMode2() == "3x3" )
57 {
58 double e3x3 = aShower.e3x3();
59 RecEmcFractionMap fracMap3x3 = aShower.getFractionMap3x3();
60 for ( cit = fracMap3x3.begin(); cit != fracMap3x3.end(); cit++ )
61 {
62 w = offset + log( ( cit->second.getEnergy() / e3x3 ) * cit->second.getFraction() );
63 if ( w > 0 )
64 {
65 wsum += w;
66 pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
67 possum += pos * w;
68 }
69 }
70 }
71 else
72 {
73 e5x5 = aShower.e5x5();
74 offset = a0 - 1.594 * exp( -2.543 * e5x5 );
75 RecEmcFractionMap fracMap5x5 = aShower.getFractionMap5x5();
76
77 num0 = 0;
78 num = 0;
79 for ( cit = fracMap5x5.begin(); cit != fracMap5x5.end(); cit++ )
80 {
81 w1 = offset + log( ( cit->second.getEnergy() / e5x5 ) * cit->second.getFraction() );
82 w2 = ( cit->second.getEnergy() / e5x5 ) * cit->second.getFraction();
83
84 num0++;
85
86 if ( w1 > 0 ) { num++; }
87 }
88
89 // cout<<"num0=="<<num0<<endl;
90 // cout<<"num==="<<num<<endl;
91 for ( cit = fracMap5x5.begin(); cit != fracMap5x5.end(); cit++ )
92 {
93 w1 = offset + log( ( cit->second.getEnergy() / e5x5 ) * cit->second.getFraction() );
94 w2 = ( cit->second.getEnergy() / e5x5 ) * cit->second.getFraction();
95
96 num0++;
97
98 if ( w1 > 0 )
99 {
100 if ( num > 1 )
101 {
102 // cout<<"w1 be used"<<endl;
103 wsum += w1;
104 pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
105 possum += pos * w1;
106 }
107 else
108 {
109 // cout<<"w2 be used"<<endl;
110 wsum += w2;
111 pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
112 possum += pos * w2;
113 }
114 }
115 }
116 }
117 if ( wsum <= 0 )
118 {
119 // cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
120 }
121 pos = possum / wsum;
122
123 // PosCorr=0 without position correction
124 // PosCorr=1 with position correction
125 if ( Para.PosCorr() == 0 ) { aShower.setPosition( pos ); }
126
127 //----------------------------------------------------------------------
128 // position correction
129 RecEmcID id = aShower.getShowerId();
130 const unsigned int module = EmcID::barrel_ec( id );
131 ;
132 const unsigned int thetaModule = EmcID::theta_module( id );
133 const unsigned int phiModule = EmcID::phi_module( id );
134 HepPoint3D posCorr( pos );
135
136 if ( module == 1 )
137 { // barrel
138 unsigned int theModule;
139 if ( thetaModule > 21 )
140 {
141 theModule = 43 - thetaModule;
142 id = EmcID::crystal_id( module, theModule, phiModule );
143 posCorr.setZ( -posCorr.z() );
144 }
145 else { theModule = thetaModule; }
146
147 double IRShift = 0, parTheta[5], parPhi[5];
148 if ( theModule == 21 ) { IRShift = 0; }
149 else if ( theModule == 20 ) { IRShift = 2.5; }
150 else { IRShift = 5.; }
151
152 EmcRecParameter& Para = EmcRecParameter::GetInstance();
153 for ( int i = 0; i < 5; i++ )
154 {
155
156 parTheta[i] = Para.BarrLoglinThetaPara( theModule, i );
157
158 parPhi[i] = Para.BarrLoglinPhiPara( 0, i );
159
160 // cout<<"Para.BarrLoglinThetaPara"<<i<<"===="<<Para.BarrLoglinThetaPara(theModule,i)<<endl;
161 // cout<<"Para.BarrLoglinPhiPara"<<i<<"===="<<Para.BarrLoglinPhiPara(0,i)<<endl;
162
163 // cout<<"theModule===="<<theModule<<endl;
164 // cout<<"partheta"<<i<<"===="<<parTheta[i]<<endl;
165 // cout<<"parphi"<<i<<"===="<<parPhi[i]<<endl;
166 }
167
168 HepPoint3D IR( 0, 0, IRShift );
169 HepPoint3D center01, center23; // center of point0,1 and point2,3 in crystal
170 center01 = ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
171 center23 = ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
172
173 double theta01, theta23, length, d;
174 theta01 = ( center01 - IR ).theta();
175 theta23 = ( center23 - IR ).theta();
176 length = ( center01 - IR ).mag();
177 d = length * tan( theta01 - theta23 ); // length in x direction
178
179 double xIn, xOut, deltaTheta;
180 xIn = length * tan( theta01 - ( posCorr - IR ).theta() ) - d / 2;
181 xOut = xIn - ( parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
182 parTheta[3] );
183 deltaTheta = atan( ( xOut + d / 2 ) / length ) - atan( ( xIn + d / 2 ) / length );
184 // cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
185 // cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
186 // + parTheta[2]*xIn + parTheta[3] <<endl;
187
188 // obtained by photon, not used, just for comparison
189 // double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
190 double yIn, yOut, deltaPhi;
191 // yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
192 // : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
193 // yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
194
195 yIn = posCorr.phi() < 0 ? posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
196 : posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
197
198 if ( phiModule == 0 && posCorr.phi() < 0 )
199 { yIn = posCorr.phi() * 180. / CLHEP::pi + 360. - 119 * 3 - 1.843; }
200 else if ( phiModule == 119 && posCorr.phi() > 0 )
201 { yIn = posCorr.phi() * 180. / CLHEP::pi - 1.843; }
202 else
203 {
204 yIn = posCorr.phi() < 0
205 ? posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
206 : posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
207 }
208 yOut = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
209 deltaPhi = yOut * CLHEP::pi / 180.;
210
211 HepPoint3D rotateVector( 0, 1, 0 ); // y axis
212 rotateVector.rotateZ( center01.phi() );
213 posCorr.setZ( posCorr.z() - IRShift );
214 posCorr.rotate( -deltaTheta, rotateVector );
215 posCorr.setZ( posCorr.z() + IRShift );
216 // phi parameter is gotten by the average of e+ e- peak
217
218 posCorr.rotateZ( -deltaPhi );
219
220 if ( thetaModule > 21 ) { posCorr.setZ( -posCorr.z() ); }
221 }
222 if ( Para.PosCorr() == 1 ) { aShower.setPosition( posCorr ); }
223 ////////////////////////////
224 if ( aShower.energy() < 1e-5 ) return;
225
226 double r, theta, phi;
227 double dtheta, dphi, dx, dy, dz;
228
229 r = posCorr.mag();
230 theta = posCorr.theta();
231 phi = posCorr.phi();
232 // cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
233
234 if ( aShower.energy() > 0 )
235 {
236 dtheta = Para.SigTheta( 0 ) / sqrt( aShower.energy() ) + Para.SigTheta( 1 );
237 dphi = Para.SigPhi( 0 ) / sqrt( aShower.energy() ) + Para.SigPhi( 1 );
238 }
239 else
240 {
241 dtheta = 0;
242 dphi = 0;
243 }
244 // cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
245
246 dx = fabs( iGeoSvc->GetBarrelR() * sin( phi ) * dphi );
247 dy = fabs( iGeoSvc->GetBarrelR() * cos( phi ) * dphi );
248 if ( theta > 0 )
249 dz = fabs( iGeoSvc->GetBarrelR() * dtheta / ( sin( theta ) * sin( theta ) ) );
250 else dz = 0;
251 // cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
252
253 aShower.setDtheta( dtheta );
254 aShower.setDphi( dphi );
255
256 //----------------------------------------------------------------------
257 // Error matrix
258 HepSymMatrix matrix( 3 ); // error matrix of r, theta, phi
259 matrix[0][0] = 0;
260 matrix[1][1] = dtheta * dtheta;
261 matrix[2][2] = dphi * dphi;
262 matrix[0][1] = 0;
263 matrix[0][2] = 0;
264 matrix[1][2] = 0;
265 // cout<<"matrix:\n"<<matrix<<endl;
266
267 HepMatrix matrix1( 3, 3 ), matrix2( 3, 3 );
268 matrix1[0][0] = sin( theta ) * cos( phi );
269 matrix1[0][1] = r * cos( theta ) * cos( phi );
270 matrix1[0][2] = -r * sin( theta ) * sin( phi );
271 matrix1[1][0] = sin( theta ) * sin( phi );
272 matrix1[1][1] = r * cos( theta ) * sin( phi );
273 matrix1[1][2] = r * sin( theta ) * cos( phi );
274 matrix1[2][0] = cos( theta );
275 matrix1[2][1] = -r * sin( theta );
276 matrix1[2][2] = 0;
277 // cout<<"matrix1:\n"<<matrix1<<endl;
278
279 for ( int i = 0; i < 3; i++ )
280 {
281 for ( int j = 0; j < 3; j++ ) { matrix2[i][j] = matrix1[j][i]; }
282 }
283 // cout<<"matrix2:\n"<<matrix2<<endl;
284
285 HepMatrix matrix4( 3, 3 );
286 matrix4 = matrix1 * matrix * matrix2;
287 // cout<<"matrix4:\n"<<matrix4<<endl;
288
289 HepSymMatrix matrix5( 3 ); // error matrix of x, y, z
290 for ( int i = 0; i < 3; i++ )
291 {
292 for ( int j = 0; j <= i; j++ ) { matrix5[i][j] = matrix4[i][j]; }
293 }
294 // cout<<"matrix5:\n"<<matrix5<<endl;
295
296 aShower.setErrorMatrix( matrix5 );
297
298 if ( matrix5[0][0] > 0 ) dx = sqrt( matrix5[0][0] );
299 if ( matrix5[1][1] > 0 ) dy = sqrt( matrix5[1][1] );
300 if ( matrix5[2][2] > 0 ) dz = sqrt( matrix5[2][2] );
301}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
etot
Double_t etot
Definition DataBase/tau_mode.c:8
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/src/ekhara/common.h:4
RecEmcFractionMap
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
Definition Emc/EmcRecEventModel/include/EmcRecEventModel/RecEmcFraction.h:55
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
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::PosCorr
double PosCorr() const
Definition EmcRecParameter.cxx:524
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition EmcRecParameter.cxx:537
EmcRecParameter::BarrLoglinThetaPara
double BarrLoglinThetaPara(int n, int m) const
Definition EmcRecParameter.cxx:555
EmcRecParameter::BarrLoglinPhiPara
double BarrLoglinPhiPara(int n, int m) const
Definition EmcRecParameter.cxx:559
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::getFractionMap5x5
RecEmcFractionMap getFractionMap5x5() const
Definition RecEmcShower.cxx:206
RecEmcShower::Begin
RecEmcFractionMap::const_iterator Begin() const
Definition RecEmcShower.cxx:112
RecEmcShower::getFractionMap3x3
RecEmcFractionMap getFractionMap3x3() const
Definition RecEmcShower.cxx:183

◆ Position() [2/3]

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

Implements EmcRecShowerPosAbs.

◆ Position() [3/3]

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

Implements EmcRecShowerPosAbs.

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