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

#include <EmcRecShowerPosLinShMax.h>

Inheritance diagram for EmcRecShowerPosLinShMax:

Public Member Functions

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

Constructor & Destructor Documentation

◆ EmcRecShowerPosLinShMax() [1/3]

EmcRecShowerPosLinShMax::EmcRecShowerPosLinShMax ( )
inline

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

17{}

◆ ~EmcRecShowerPosLinShMax() [1/3]

virtual EmcRecShowerPosLinShMax::~EmcRecShowerPosLinShMax ( )
inlinevirtual

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

18{}

◆ EmcRecShowerPosLinShMax() [2/3]

EmcRecShowerPosLinShMax::EmcRecShowerPosLinShMax ( )
inline

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

17{}

◆ ~EmcRecShowerPosLinShMax() [2/3]

virtual EmcRecShowerPosLinShMax::~EmcRecShowerPosLinShMax ( )
inlinevirtual

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

18{}

◆ EmcRecShowerPosLinShMax() [3/3]

EmcRecShowerPosLinShMax::EmcRecShowerPosLinShMax ( )
inline

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

17{}

◆ ~EmcRecShowerPosLinShMax() [3/3]

virtual EmcRecShowerPosLinShMax::~EmcRecShowerPosLinShMax ( )
inlinevirtual

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

18{}

Member Function Documentation

◆ Position() [1/3]

void EmcRecShowerPosLinShMax::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 20 of file EmcRecShowerPosLinShMax.cxx.

20 {
21 RecEmcFractionMap::const_iterator cit;
22 HepPoint3D pos( 0, 0, 0 );
23 HepPoint3D posMax( 0, 0, 0 );
24 HepPoint3D possum( 0, 0, 0 );
25 double w, wsum = 0;
26 double etot;
27
28 // cout<<"EmcRecShowerPosLinShMax::Position Begin"<<endl;
29
30 IEmcRecGeoSvc* iGeoSvc;
31 ISvcLocator* svcLocator = Gaudi::svcLocator();
32 StatusCode sc = svcLocator->service( "EmcRecGeoSvc", iGeoSvc );
33 if ( sc != StatusCode::SUCCESS )
34 {
35 // cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
36 }
37
38 EmcRecParameter& Para = EmcRecParameter::GetInstance();
39
40 for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
41 {
42 etot += ( cit->second.getEnergy() * cit->second.getFraction() );
43 // cout<<etot<<endl;
44 }
45 wsum = 0;
46
47 double etota = aShower.getEAll();
48 double ltot = ( log( etota / 0.0145 ) + 0.5 ) * 1.86;
49 double e55 = aShower.e5x5();
50 double len55 = ( log( etota / 0.0145 ) + 0.5 ) * 1.86;
51 Identifier seedID = aShower.getShowerId();
52 HepPoint3D seedPoint( iGeoSvc->GetCFrontCenter( seedID ) );
53 // Rs is the distance from orgin to showerCentroid, unit cm
54 // HepDouble Rs=seedPoint.mag();
55 for ( cit = aShower.Begin(); cit != aShower.End(); cit++ )
56 {
57 w = ( cit->second.getEnergy() * cit->second.getFraction() );
58
59 HepPoint3D hitPoint( iGeoSvc->GetCFrontCenter( cit->second.getCellId() ) );
60 // HepDouble Ang=seedPoint.angle(hitPoint);
61 // theDistance is the distance from orgin to showerMax
62 double theDistance = 0;
63 if ( cit->second.getCellId() == seedID )
64 {
65 // theDistance=(Rs+ltot);
66 theDistance = ( seedPoint.mag() + ltot );
67 }
68 else
69 {
70 theDistance = ( seedPoint.mag() + ltot ) / cos( seedPoint.angle( hitPoint ) );
71 // theDistance=(Rs+ltot)/cos(Ang);
72 }
73
74 // cout<<w<<endl;
75 wsum += w;
76 pos = iGeoSvc->GetCFrontCenter( cit->second.getCellId() );
77
78 posMax = ( theDistance / pos.mag() ) * pos;
79
80 possum += posMax * w;
81 }
82
83 if ( wsum <= 0 ) { cout << "[[[[[[[[[[[[[[[" << wsum; }
84 pos = possum / wsum;
85 // pos=((pos.mag()-len55)/pos.mag())*pos;
86 // aShower.setPosition(pos);
87 if ( Para.PosCorr() == 0 ) { aShower.setPosition( pos ); }
88
89 //----------------------------------------------------------------------
90 // position correction
91 RecEmcID id = aShower.getShowerId();
92 const unsigned int module = EmcID::barrel_ec( id );
93 ;
94 const unsigned int thetaModule = EmcID::theta_module( id );
95 const unsigned int phiModule = EmcID::phi_module( id );
96 HepPoint3D posCorr( pos );
97
98 if ( module == 1 )
99 { // barrel
100 unsigned int theModule;
101 if ( thetaModule > 21 )
102 {
103 theModule = 43 - thetaModule;
104 id = EmcID::crystal_id( module, theModule, phiModule );
105 posCorr.setZ( -posCorr.z() );
106 }
107 else { theModule = thetaModule; }
108
109 double IRShift, parTheta[5], parPhi[5];
110 if ( theModule == 21 ) { IRShift = 0; }
111 else if ( theModule == 20 ) { IRShift = 2.5; }
112 else { IRShift = 5.; }
113
114 EmcRecParameter& Para = EmcRecParameter::GetInstance();
115
116 for ( int i = 0; i < 5; i++ )
117 {
118
119 if ( e55 > 1.0 ) parTheta[i] = Para.BarrShLinThetaPara( theModule, i );
120 else if ( e55 <= 1.0 && e55 > 0.5 )
121 parTheta[i] = Para.BarrShLinThetaPara( theModule + 22, i );
122 else if ( e55 <= 0.5 ) parTheta[i] = Para.BarrShLinThetaPara( theModule + 44, i );
123
124 if ( e55 > 1.0 ) parPhi[i] = Para.BarrShLinPhiPara( 0, i );
125 else if ( e55 <= 1.0 && e55 > 0.5 ) parPhi[i] = Para.BarrShLinPhiPara( 1, i );
126 else if ( e55 <= 0.5 ) parPhi[i] = Para.BarrShLinPhiPara( 2, i );
127 }
128
129 HepPoint3D IR( 0, 0, IRShift );
130 HepPoint3D center01, center23; // center of point0,1 and point2,3 in crystal
131 center01 = ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
132 center23 = ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
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
140 double xIn, xOut, deltaTheta;
141 xIn = length * tan( theta01 - ( posCorr - IR ).theta() ) - d / 2;
142 xOut = xIn - ( parTheta[0] * atan( xIn * parTheta[1] - parTheta[4] ) + parTheta[2] * xIn +
143 parTheta[3] );
144 deltaTheta = atan( ( xOut + d / 2 ) / length ) - atan( ( xIn + d / 2 ) / length );
145 // cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
146 // cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
147 // + parTheta[2]*xIn + parTheta[3] <<endl;
148 /////////////////////////////////////////////////////////////////////////////////////////////////
149 double yIn, yOut, deltaPhi;
150
151 if ( phiModule == 0 && posCorr.phi() < 0 )
152 { yIn = posCorr.phi() * 180. / CLHEP::pi + 360. - 119 * 3 - 1.843; }
153 else if ( phiModule == 119 && posCorr.phi() > 0 )
154 { yIn = posCorr.phi() * 180. / CLHEP::pi - 1.843; }
155 else
156 {
157
158 yIn = posCorr.phi() < 0
159 ? posCorr.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
160 : posCorr.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
161 }
162 yOut = parPhi[0] * atan( yIn * parPhi[1] - parPhi[4] ) + parPhi[2] * yIn + parPhi[3];
163
164 // yOut =
165 // parPhi[0]+parPhi[1]*yIn+parPhi[2]*yIn*yIn+parPhi[3]*yIn*yIn*yIn+parPhi[4]*yIn*yIn*yIn*yIn+parPhi[5]*yIn*yIn*yIn*yIn*yIn;
166
167 deltaPhi = yOut * CLHEP::pi / 180.;
168
169 HepPoint3D rotateVector( 0, 1, 0 ); // y axis
170 rotateVector.rotateZ( center01.phi() );
171 posCorr.setZ( posCorr.z() - IRShift );
172 posCorr.rotate( -deltaTheta, rotateVector );
173 posCorr.setZ( posCorr.z() + IRShift );
174 // phi parameter is gotten by the average of e+ e- peak
175 // posCorr.rotateZ(-0.002994);
176 posCorr.rotateZ( -deltaPhi );
177
178 if ( thetaModule > 21 ) { posCorr.setZ( -posCorr.z() ); }
179 }
180 posCorr = ( ( posCorr.mag() - len55 ) / posCorr.mag() ) * posCorr;
181 // aShower.setPosition(posCorr);
182 if ( Para.PosCorr() == 1 ) { aShower.setPosition( posCorr ); }
183
184 ////////////////////////////
185 if ( aShower.energy() < 1e-5 ) return;
186
187 double r, theta, phi;
188 double dtheta, dphi, dx, dy, dz;
189
190 r = posCorr.mag();
191 theta = posCorr.theta();
192 phi = posCorr.phi();
193 // cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
194
195 if ( aShower.energy() > 0 )
196 {
197 dtheta = Para.SigTheta( 0 ) / sqrt( aShower.energy() ) + Para.SigTheta( 1 );
198 dphi = Para.SigPhi( 0 ) / sqrt( aShower.energy() ) + Para.SigPhi( 1 );
199 }
200 else
201 {
202 dtheta = 0;
203 dphi = 0;
204 }
205 // cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
206
207 dx = fabs( iGeoSvc->GetBarrelR() * sin( phi ) * dphi );
208 dy = fabs( iGeoSvc->GetBarrelR() * cos( phi ) * dphi );
209 if ( theta > 0 )
210 dz = fabs( iGeoSvc->GetBarrelR() * dtheta / ( sin( theta ) * sin( theta ) ) );
211 else dz = 0;
212 // cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
213
214 aShower.setDtheta( dtheta );
215 aShower.setDphi( dphi );
216
217 //----------------------------------------------------------------------
218 // Error matrix
219 HepSymMatrix matrix( 3 ); // error matrix of r, theta, phi
220 matrix[0][0] = 0;
221 matrix[1][1] = dtheta * dtheta;
222 matrix[2][2] = dphi * dphi;
223 matrix[0][1] = 0;
224 matrix[0][2] = 0;
225 matrix[1][2] = 0;
226 // cout<<"matrix:\n"<<matrix<<endl;
227
228 HepMatrix matrix1( 3, 3 ), matrix2( 3, 3 );
229 matrix1[0][0] = sin( theta ) * cos( phi );
230 matrix1[0][1] = r * cos( theta ) * cos( phi );
231 matrix1[0][2] = -r * sin( theta ) * sin( phi );
232 matrix1[1][0] = sin( theta ) * sin( phi );
233 matrix1[1][1] = r * cos( theta ) * sin( phi );
234 matrix1[1][2] = r * sin( theta ) * cos( phi );
235 matrix1[2][0] = cos( theta );
236 matrix1[2][1] = -r * sin( theta );
237 matrix1[2][2] = 0;
238 // cout<<"matrix1:\n"<<matrix1<<endl;
239
240 for ( int i = 0; i < 3; i++ )
241 {
242 for ( int j = 0; j < 3; j++ ) { matrix2[i][j] = matrix1[j][i]; }
243 }
244 // cout<<"matrix2:\n"<<matrix2<<endl;
245
246 HepMatrix matrix4( 3, 3 );
247 matrix4 = matrix1 * matrix * matrix2;
248 // cout<<"matrix4:\n"<<matrix4<<endl;
249
250 HepSymMatrix matrix5( 3 ); // error matrix of x, y, z
251 for ( int i = 0; i < 3; i++ )
252 {
253 for ( int j = 0; j <= i; j++ ) { matrix5[i][j] = matrix4[i][j]; }
254 }
255 // cout<<"matrix5:\n"<<matrix5<<endl;
256
257 aShower.setErrorMatrix( matrix5 );
258
259 if ( matrix5[0][0] > 0 ) dx = sqrt( matrix5[0][0] );
260 if ( matrix5[1][1] > 0 ) dy = sqrt( matrix5[1][1] );
261 if ( matrix5[2][2] > 0 ) dz = sqrt( matrix5[2][2] );
262}
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::PosCorr
double PosCorr() const
Definition EmcRecParameter.cxx:524
EmcRecParameter::BarrShLinPhiPara
double BarrShLinPhiPara(int n, int m) const
Definition EmcRecParameter.cxx:581
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition EmcRecParameter.cxx:537
EmcRecParameter::BarrShLinThetaPara
double BarrShLinThetaPara(int n, int m) const
Definition EmcRecParameter.cxx:577
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 EmcRecShowerPosLinShMax::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

◆ Position() [3/3]

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

Implements EmcRecShowerPosAbs.

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