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

#include <TrkPocaXY.h>

Inheritance diagram for TrkPocaXY:

Public Member Functions

 TrkPocaXY (const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
 TrkPocaXY (const Trajectory &traj1, double flt1, const Trajectory &traj2, double flt2, double precision=1.0e-4)
 ~TrkPocaXY ()
double docaXY () const
double fltl1 () const
double fltl2 () const
 TrkPocaXY (const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
 TrkPocaXY (const Trajectory &traj1, double flt1, const Trajectory &traj2, double flt2, double precision=1.0e-4)
 ~TrkPocaXY ()
double docaXY () const
double fltl1 () const
double fltl2 () const
 TrkPocaXY (const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
 TrkPocaXY (const Trajectory &traj1, double flt1, const Trajectory &traj2, double flt2, double precision=1.0e-4)
 ~TrkPocaXY ()
double docaXY () const
double fltl1 () const
double fltl2 () const
Public Member Functions inherited from TrkPocaBase
const TrkErrCodestatus () const
double flt1 () const
double flt2 () const
double precision ()
const TrkErrCodestatus () const
double flt1 () const
double flt2 () const
double precision ()
const TrkErrCodestatus () const
double flt1 () const
double flt2 () const
double precision ()

Additional Inherited Members

Protected Member Functions inherited from TrkPocaBase
 TrkPocaBase (double flt1, double flt2, double precision)
 TrkPocaBase (double flt1, double precision)
virtual ~TrkPocaBase ()
void minimize (const Trajectory &traj1, double f1, const Trajectory &traj2, double f2)
void minimize (const Trajectory &traj1, double f1, const HepPoint3D &pt)
void stepTowardPoca (const Trajectory &traj1, const Trajectory &traj2)
void stepToPointPoca (const Trajectory &traj, const HepPoint3D &pt)
 TrkPocaBase (double flt1, double flt2, double precision)
 TrkPocaBase (double flt1, double precision)
virtual ~TrkPocaBase ()
void minimize (const Trajectory &traj1, double f1, const Trajectory &traj2, double f2)
void minimize (const Trajectory &traj1, double f1, const HepPoint3D &pt)
void stepTowardPoca (const Trajectory &traj1, const Trajectory &traj2)
void stepToPointPoca (const Trajectory &traj, const HepPoint3D &pt)
 TrkPocaBase (double flt1, double flt2, double precision)
 TrkPocaBase (double flt1, double precision)
virtual ~TrkPocaBase ()
void minimize (const Trajectory &traj1, double f1, const Trajectory &traj2, double f2)
void minimize (const Trajectory &traj1, double f1, const HepPoint3D &pt)
void stepTowardPoca (const Trajectory &traj1, const Trajectory &traj2)
void stepToPointPoca (const Trajectory &traj, const HepPoint3D &pt)
Protected Attributes inherited from TrkPocaBase
double _precision
double _flt1
double _flt2
TrkErrCode _status
Static Protected Attributes inherited from TrkPocaBase
static double _maxDist = 1.e7
static int _maxTry = 500
static double _extrapToler = 2.

Detailed Description

Constructor & Destructor Documentation

◆ TrkPocaXY() [1/6]

TrkPocaXY::TrkPocaXY ( const Trajectory & traj,
double flt,
const HepPoint3D & pt,
double precision = 1.0e-4 )

Definition at line 23 of file TrkPocaXY.cxx.

24 : TrkPocaBase( flt, prec ), _docaxy( -9999.0 ) {
25 // construct a line trajectory through the given point
26 Hep3Vector zaxis( 0.0, 0.0, 1.0 );
27 // length is set to 200cm (shouldn't matter)
28 // TrkLineTraj zline(pt, zaxis, 200.0);//yzhang del
29 TrkLineTraj zline( pt, zaxis, 2000.0 ); // yzhang change
30 // create a 2-traj poca with this
31 double zlen( pt.z() );
32 TrkPoca zlinepoca( traj, flt, zline, zlen, prec );
33 // transfer over the data
34 _status = zlinepoca.status();
35 if ( status().success() )
36 {
37 _flt1 = zlinepoca.flt1();
38 _flt2 = zlinepoca.flt2();
39 _docaxy =
40 zlinepoca.doca(); // doca should be perpendicular to the zline so 2d by construction
41 }
42}
TrkPocaBase(double flt1, double flt2, double precision)

Referenced by TrkPocaXY().

◆ TrkPocaXY() [2/6]

TrkPocaXY::TrkPocaXY ( const Trajectory & traj1,
double flt1,
const Trajectory & traj2,
double flt2,
double precision = 1.0e-4 )

Definition at line 44 of file TrkPocaXY.cxx.

46 : TrkPocaBase( flt1, flt2, prec ), _docaxy( -9999.0 ) {
47
48 // this traj-traj version starts by approximating the trejectories as
49 // either a line or a circle and treats separately the line-line,
50 // circle-circle, and line-circle cases.
51
52 _flt1 = flt1;
53 _flt2 = flt2;
54 double delta1( 10 * prec );
55 double delta2( 10 * prec );
56 unsigned niter( 0 );
57 static const unsigned maxiter( 20 );
58 while ( niter < maxiter && ( delta1 > prec || delta2 > prec ) )
59 {
60 // get positions and directions
61 HepPoint3D pos1 = traj1.position( _flt1 );
62 HepPoint3D pos1b;
63 if ( niter == 0 ) pos1b = pos1;
64 HepPoint3D pos2 = traj2.position( _flt2 );
65 Hep3Vector dir1 = traj1.direction( _flt1 );
66 Hep3Vector dir2 = traj2.direction( _flt2 );
67 Hep3Vector dd1 = traj1.delDirect( _flt1 );
68 Hep3Vector dd2 = traj2.delDirect( _flt2 );
69 double curv1 = traj1.curvature( _flt1 );
70 double curv2 = traj2.curvature( _flt2 );
71 double m1, m2, q1, q2;
72 double r1, r2, xc1, xc2, yc1, yc2, sr1, sr2;
73 if ( curv1 == 0 )
74 {
75 // convert to m*x+q
76 if ( dir1[0] == 0 ) { m1 = 1.0e12; }
77 else { m1 = dir1[1] / dir1[0]; }
78 q1 = pos1.y() - pos1.x() * m1;
79 }
80 else
81 {
82 // get circle parameters
83 r1 = ( 1 - dir1[2] * dir1[2] ) / curv1;
84 sr1 = r1;
85 if ( dir1[0] * dd1[1] - dir1[1] * dd1[0] < 0 ) sr1 = -r1;
86 double cosphi1 = dir1[0] / sqrt( dir1[0] * dir1[0] + dir1[1] * dir1[1] );
87 double sinphi1 = cosphi1 * dir1[1] / dir1[0];
88 xc1 = pos1.x() - sr1 * sinphi1;
89 yc1 = pos1.y() + sr1 * cosphi1;
90 }
91 if ( curv2 == 0 )
92 {
93 if ( dir2[0] == 0 ) { m2 = 1.0e12; }
94 else { m2 = dir2[1] / dir2[0]; }
95 q2 = pos2.y() - pos2.x() * m2;
96 }
97 else
98 {
99 r2 = ( 1 - dir2[2] * dir2[2] ) / curv2;
100 sr2 = r2;
101 if ( dir2[0] * dd2[1] - dir2[1] * dd2[0] < 0 ) sr2 = -r2;
102 double cosphi2 = dir2[0] / sqrt( dir2[0] * dir2[0] + dir2[1] * dir2[1] );
103 double sinphi2 = cosphi2 * dir2[1] / dir2[0];
104 xc2 = pos2.x() - sr2 * sinphi2;
105 yc2 = pos2.y() + sr2 * cosphi2;
106 }
107 double xint, yint, xint1, xint2, yint1, yint2, absdoca;
108 _docaxy = 0;
109 _status.setSuccess( 3 );
110
111 // First the line-line case
112
113 if ( curv1 == 0 && curv2 == 0 )
114 {
115 // do the intersection in 2d
116 interTwoLines( m1, q1, m2, q2, xint, yint );
117 }
118
119 // next do the two circle case
120
121 if ( curv1 != 0 && curv2 != 0 )
122 {
123 // There are four cases
124 double cdist = sqrt( ( xc1 - xc2 ) * ( xc1 - xc2 ) + ( yc1 - yc2 ) * ( yc1 - yc2 ) );
125 // a - coincident centers
126 if ( fabs( cdist ) < 1.e-12 )
127 {
128 // the algorithm fails because the points have all the same distance
129 _status.setFailure( 12, "TrkPocaXY:: the two circles have the same center..." );
130 return;
131 }
132 // b - Actual intersection
133 if ( ( fabs( r1 - r2 ) < cdist ) && ( cdist < r1 + r2 ) )
134 {
135 interTwoCircles( xc1, yc1, r1, xc2, yc2, r2, xint1, yint1, xint2, yint2 );
136 double dist1 = ( xint1 - pos1b.x() ) * ( xint1 - pos1b.x() ) +
137 ( yint1 - pos1b.y() ) * ( yint1 - pos1b.y() );
138 double dist2 = ( xint2 - pos1b.x() ) * ( xint2 - pos1b.x() ) +
139 ( yint2 - pos1b.y() ) * ( yint2 - pos1b.y() );
140
141 // choose closest to begining of track1
142 if ( dist1 < dist2 )
143 {
144 xint = xint1;
145 yint = yint1;
146 }
147 else
148 {
149 xint = xint2;
150 yint = yint2;
151 }
152 }
153
154 // c - nested circles and d - separated circles
155
156 if ( ( fabs( r1 - r2 ) > cdist ) || // nested circles
157 ( cdist > ( r1 + r2 ) ) )
158 { // separated circles
159 // line going through the centers of the two circles
160
161 double m = ( yc1 - yc2 ) / ( xc1 - xc2 ); // y = m*x+q
162 double q = yc1 - m * xc1;
163
164 // intersection points between the line and the two circles
165
166 double xint1, yint1, xint2, yint2, zOfDCrossT1;
167
168 interLineCircle( m, q, xc1, yc1, r1, xint1, yint1, xint2, yint2 );
169
170 double xint3, yint3, xint4, yint4;
171 interLineCircle( m, q, xc2, yc2, r2, xint3, yint3, xint4, yint4 );
172 if ( fabs( r1 - r2 ) > cdist )
173 { // nested circles
174 absdoca = fabs( r1 - r2 ) - cdist;
175#ifdef MDCPATREC_DEBUG
176 cout << "ErrMsg(debugging) doing nested circles in TrkPocaXY " << endl;
177#endif
178
179 double dist1_3 = pow( ( xint1 - xint3 ), 2. ) + pow( ( yint1 - yint3 ), 2. );
180 double dist2_4 = pow( ( xint2 - xint4 ), 2. ) + pow( ( yint2 - yint4 ), 2. );
181
182 if ( dist1_3 < dist2_4 )
183 {
184 xint = 0.5 * ( xint1 + xint3 );
185 yint = 0.5 * ( yint1 + yint3 );
186 zOfDCrossT1 = ( xint3 - xint1 ) * dir1[1] - ( yint3 - yint1 ) * dir1[0];
187 }
188 else
189 {
190 xint = 0.5 * ( xint2 + xint4 );
191 yint = 0.5 * ( yint2 + yint4 );
192 zOfDCrossT1 = ( xint4 - xint2 ) * dir1[1] - ( yint4 - yint2 ) * dir1[0];
193 }
194 _docaxy = absdoca;
195 if ( zOfDCrossT1 > 0 ) _docaxy = -absdoca;
196 }
197 if ( cdist > ( r1 + r2 ) )
198 { // separated circles
199 absdoca = cdist - ( r1 + r2 );
200#ifdef MDCPATREC_DEBUG
201 cout << "ErrMsg(debugging) doing separated circles in TrkPocaXY" << endl;
202#endif
203
204 double dist2_3 = pow( ( xint2 - xint3 ), 2. ) + pow( ( yint2 - yint3 ), 2. );
205 double dist1_4 = pow( ( xint1 - xint4 ), 2. ) + pow( ( yint1 - yint4 ), 2. );
206 if ( dist2_3 < dist1_4 )
207 {
208 xint = 0.5 * ( xint2 + xint3 );
209 yint = 0.5 * ( yint2 + yint3 );
210 zOfDCrossT1 = ( xint3 - xint1 ) * dir1[1] - ( yint3 - yint1 ) * dir1[0];
211 }
212 else
213 {
214 xint = 0.5 * ( xint1 + xint4 );
215 yint = 0.5 * ( yint1 + yint4 );
216 zOfDCrossT1 = ( xint4 - xint2 ) * dir1[1] - ( yint4 - yint2 ) * dir1[0];
217 }
218 _docaxy = absdoca;
219 if ( zOfDCrossT1 > 0 ) _docaxy = -absdoca;
220 }
221 }
222 }
223
224 // Now do the line-circle case
225
226 if ( ( curv1 == 0 && curv2 != 0 ) || ( curv1 != 0 && curv2 == 0 ) )
227 {
228 // distance between the line and the circle center
229 HepVector dirT;
230 double m, q, r, xc, yc, zOfDCrossT1;
231 if ( curv1 == 0 )
232 {
233 m = m1;
234 q = q1;
235 r = r2;
236 xc = xc2;
237 yc = yc2;
238 dirT = dir2;
239 }
240 else
241 {
242 m = m2;
243 q = q2;
244 r = r1;
245 xc = xc1;
246 yc = yc1;
247 dirT = dir1;
248 }
249
250 double dist = fabs( ( m * xc - yc + q ) / sqrt( 1 + m * m ) );
251 if ( dist <= r )
252 {
253
254 // the intersection points
255
256 interLineCircle( m, q, xc, yc, r, xint1, yint1, xint2, yint2 );
257 double dist1 = ( xint1 - pos1b.x() ) * ( xint1 - pos1b.x() ) +
258 ( yint1 - pos1b.y() ) * ( yint1 - pos1b.y() );
259 double dist2 = ( xint2 - pos1b.x() ) * ( xint2 - pos1b.x() ) +
260 ( yint2 - pos1b.y() ) * ( yint2 - pos1b.y() );
261 // choose closest to the beginning of track 1
262 if ( dist1 < dist2 )
263 {
264 xint = xint1;
265 yint = yint1;
266 }
267 else
268 {
269 xint = xint2;
270 yint = yint2;
271 }
272 }
273 else
274 { // no intersection points
275#ifdef MDCPATREC_DEBUG
276 cout << "ErrMsg(debugging) doing separated line-circle in TrkPocaXY" << endl;
277#endif
278
279 // line going through the circle center and perpendicular to traj1
280
281 double mperp = -1. / m;
282 double qperp = yc - mperp * xc;
283
284 // intersection between this line and the two trajectories
285
286 interLineCircle( mperp, qperp, xc, yc, r, xint1, yint1, xint2, yint2 );
287
288 double xint3, yint3;
289 interTwoLines( m, q, mperp, qperp, xint3, yint3 );
290
291 double dist1_3 = pow( ( xint1 - xint3 ), 2. ) + pow( ( yint1 - yint3 ), 2. );
292 double dist2_3 = pow( ( xint2 - xint3 ), 2. ) + pow( ( yint2 - yint3 ), 2. );
293 if ( dist1_3 < dist2_3 )
294 {
295 xint = 0.5 * ( xint3 + xint1 );
296 yint = 0.5 * ( yint3 + yint1 );
297 absdoca = sqrt( ( xint3 - xint1 ) * ( xint3 - xint1 ) +
298 ( yint3 - yint1 ) * ( yint3 - yint1 ) );
299 zOfDCrossT1 = ( xint3 - xint1 ) * dirT[1] - ( yint3 - yint1 ) * dirT[0];
300 }
301 else
302 {
303 xint = 0.5 * ( xint3 + xint2 );
304 yint = 0.5 * ( yint3 + yint2 );
305 absdoca = sqrt( ( xint3 - xint2 ) * ( xint3 - xint2 ) +
306 ( yint3 - yint2 ) * ( yint3 - yint2 ) );
307 zOfDCrossT1 = ( xint3 - xint2 ) * dirT[1] - ( yint3 - yint2 ) * dirT[0];
308 }
309 _docaxy = absdoca;
310 if ( zOfDCrossT1 > 0 ) _docaxy = -absdoca;
311 }
312 }
313
314 // get the flight lengths for the intersection point
315
316 HepPoint3D intpt( xint, yint, 0.0 );
317 TrkPocaXY poca1( traj1, _flt1, intpt, prec );
318 TrkPocaXY poca2( traj2, _flt2, intpt, prec );
319 _status = poca2.status();
320 if ( poca1.status().success() && poca2.status().success() )
321 {
322 delta1 = fabs( _flt1 - poca1.flt1() );
323 delta2 = fabs( _flt2 - poca2.flt1() );
324 _flt1 = poca1.flt1();
325 _flt2 = poca2.flt1();
326 }
327 else
328 {
329#ifdef MDCPATREC_WARNING
330 cout << "ErrMsg(warning)"
331 << " poca fialure " << endl;
332#endif
333 if ( poca1.status().failure() ) _status = poca1.status();
334 break;
335 }
336 niter++;
337 }
338#ifdef MDCPATREC_DEBUG
339 cout << "ErrMsg(debugging)"
340 << "TrkPocaXY iterations = " << niter - 1 << " flight lengths " << _flt1 << " "
341 << _flt2 << endl
342 << " deltas " << delta1 << " " << delta2 << endl;
343#endif
344 if ( niter == maxiter )
345 {
346#ifdef MDCPATREC_ROUTINE
347 cout << "ErrMsg(routine)"
348 << "TrkPocaXY:: did not converge" << endl;
349#endif
350 _status.setFailure( 13, "TrkPocaXY:: did not converge" );
351 }
352}
HepGeom::Point3D< double > HepPoint3D
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
virtual HepPoint3D position(double) const =0
virtual Hep3Vector delDirect(double) const =0
virtual Hep3Vector direction(double) const =0
virtual double curvature(double) const =0
TrkPocaXY(const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
Definition TrkPocaXY.cxx:23
double double * m2
Definition qcdloop1.h:83
double * m1
Definition qcdloop1.h:83

◆ ~TrkPocaXY() [1/3]

TrkPocaXY::~TrkPocaXY ( )
inline

◆ TrkPocaXY() [3/6]

TrkPocaXY::TrkPocaXY ( const Trajectory & traj,
double flt,
const HepPoint3D & pt,
double precision = 1.0e-4 )

◆ TrkPocaXY() [4/6]

TrkPocaXY::TrkPocaXY ( const Trajectory & traj1,
double flt1,
const Trajectory & traj2,
double flt2,
double precision = 1.0e-4 )

◆ ~TrkPocaXY() [2/3]

TrkPocaXY::~TrkPocaXY ( )
inline

◆ TrkPocaXY() [5/6]

TrkPocaXY::TrkPocaXY ( const Trajectory & traj,
double flt,
const HepPoint3D & pt,
double precision = 1.0e-4 )

◆ TrkPocaXY() [6/6]

TrkPocaXY::TrkPocaXY ( const Trajectory & traj1,
double flt1,
const Trajectory & traj2,
double flt2,
double precision = 1.0e-4 )

◆ ~TrkPocaXY() [3/3]

TrkPocaXY::~TrkPocaXY ( )
inline

Member Function Documentation

◆ docaXY() [1/3]

double TrkPocaXY::docaXY ( ) const
inline

Definition at line 66 of file InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkPocaXY.h.

66{ return _docaxy; }

◆ docaXY() [2/3]

double TrkPocaXY::docaXY ( ) const
inline

◆ docaXY() [3/3]

double TrkPocaXY::docaXY ( ) const
inline

◆ fltl1() [1/3]

double TrkPocaXY::fltl1 ( ) const
inline

Definition at line 47 of file InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkPocaXY.h.

47{ return flt1(); }

◆ fltl1() [2/3]

double TrkPocaXY::fltl1 ( ) const
inline

Definition at line 47 of file InstallArea/x86_64-el9-gcc13-opt/include/TrkBase/TrkPocaXY.h.

47{ return flt1(); }

◆ fltl1() [3/3]

double TrkPocaXY::fltl1 ( ) const
inline

Definition at line 47 of file Reconstruction/MdcPatRec/TrkBase/include/TrkBase/TrkPocaXY.h.

47{ return flt1(); }

◆ fltl2() [1/3]

double TrkPocaXY::fltl2 ( ) const
inline

Definition at line 48 of file InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkPocaXY.h.

48{ return flt2(); }

◆ fltl2() [2/3]

double TrkPocaXY::fltl2 ( ) const
inline

Definition at line 48 of file InstallArea/x86_64-el9-gcc13-opt/include/TrkBase/TrkPocaXY.h.

48{ return flt2(); }

◆ fltl2() [3/3]

double TrkPocaXY::fltl2 ( ) const
inline

Definition at line 48 of file Reconstruction/MdcPatRec/TrkBase/include/TrkBase/TrkPocaXY.h.

48{ return flt2(); }

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