BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TrkPocaXY.cxx
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// File and Version Information:
3// $Id: TrkPocaXY.cxx,v 1.3 2005/12/01 06:18:42 zhangy Exp $
4//
5// Description:
6//
7//
8// Environment:
9// Software developed for the BaBar Detector at the SLAC B-Factory.
10//
11// Author(s): Steve Schaffner, largely taken from Art Snyder
12//
13//------------------------------------------------------------------------
14#include "TrkBase/TrkPocaXY.h"
15#include "CLHEP/Vector/ThreeVector.h"
16#include "TrkBase/HelixTraj.h"
17#include "TrkBase/TrkExchangePar.h"
18#include "TrkBase/TrkLineTraj.h"
19#include "TrkBase/TrkPoca.h"
20using std::cout;
21using std::endl;
22
23TrkPocaXY::TrkPocaXY( const Trajectory& traj, double flt, const HepPoint3D& pt, double prec )
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}
43
44TrkPocaXY::TrkPocaXY( const Trajectory& traj1, double flt1, const Trajectory& traj2,
45 double flt2, double prec )
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}
353//------------------------------------------------------------------------
354void TrkPocaXY::interLineCircle( const double& m, const double& q, const double& xc,
355 const double& yc, const double& r, double& xint1,
356 double& yint1, double& xint2, double& yint2 )
357//-------------------------------------------------------------------------
358{
359
360 double alpha = 1 + m * m;
361
362 double beta = -xc + m * ( q - yc );
363
364 double gamma = xc * xc + ( q - yc ) * ( q - yc ) - r * r;
365
366 double Delta = beta * beta - alpha * gamma;
367
368 if ( Delta < 0 )
369 {
370
371 _status.setFailure( 14, "TrkPocaXY:: the line and the circle heve no intersections..." );
372 return;
373 }
374 else if ( fabs( Delta ) < 1.e-12 )
375 {
376
377 xint1 = -beta / alpha;
378 xint2 = xint1;
379 }
380 else
381 {
382
383 double xPlus = -beta / alpha + sqrt( beta * beta - alpha * gamma ) / alpha;
384 double xMinus = -beta / alpha - sqrt( beta * beta - alpha * gamma ) / alpha;
385
386 if ( xPlus > xMinus )
387 {
388 xint1 = xMinus;
389 xint2 = xPlus;
390 }
391 else
392 {
393 xint1 = xPlus;
394 xint2 = xMinus;
395 }
396 }
397 yint1 = m * xint1 + q;
398 yint2 = m * xint2 + q;
399
400 return;
401}
402//------------------------------------------------------------------------
403void TrkPocaXY::interTwoCircles( const double& xc1, const double& yc1, const double& r1,
404 const double& xc2, const double& yc2, const double& r2,
405 double& xint1, double& yint1, double& xint2, double& yint2 )
406//-------------------------------------------------------------------------
407{
408
409 double A = ( xc1 * xc1 + yc1 * yc1 - r1 * r1 ) - ( xc2 * xc2 + yc2 * yc2 - r2 * r2 );
410
411 double B = -xc1 + xc2;
412
413 double C = -yc1 + yc2;
414
415 double alpha = 1 + ( B * B ) / ( C * C );
416
417 double beta = -xc1 + B / C * ( yc1 + A / ( 2 * C ) );
418
419 double gamma = xc1 * xc1 + ( yc1 + A / ( 2 * C ) ) * ( yc1 + A / ( 2 * C ) ) - r1 * r1;
420
421 double Delta = beta * beta - alpha * gamma;
422
423 if ( Delta < 0 )
424 {
425
426 _status.setFailure( 14, "TrkPocaXY:: the two circles have no intersections..\
427." );
428 return;
429 }
430 else if ( fabs( Delta ) < 1.e-12 )
431 {
432
433 xint1 = -beta / alpha;
434 xint2 = xint1;
435 }
436 else
437 {
438 double xPlus = -beta / alpha + sqrt( beta * beta - alpha * gamma ) / alpha;
439 double xMinus = -beta / alpha - sqrt( beta * beta - alpha * gamma ) / alpha;
440
441 if ( xPlus > xMinus )
442 {
443 xint1 = xMinus;
444 xint2 = xPlus;
445 }
446 else
447 {
448 xint1 = xPlus;
449 xint2 = xMinus;
450 }
451 }
452
453 yint1 = -( A + 2 * B * xint1 ) / ( 2 * C );
454 yint2 = -( A + 2 * B * xint2 ) / ( 2 * C );
455
456 return;
457}
458
459//------------------------------------------------------------------------
460void TrkPocaXY::interTwoLines( const double& m1, const double& q1, const double& m2,
461 const double& q2, double& xint, double& yint )
462//-------------------------------------------------------------------------
463{
464
465 if ( fabs( m1 - m2 ) < 1.e-12 )
466 { // parallel lines
467
468 // the algorithm fails because the points have all the same distance
469
470 _status.setFailure( 13, "TrkPocaXY:: the two lines are parallel..." );
471 return;
472 }
473 else
474 { // the lines have an intersection point
475
476 xint = ( q2 - q1 ) / ( m1 - m2 );
477 yint = m1 * xint + q1;
478 }
479
480 return;
481}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
alpha
double alpha
Definition EvtPhokharaDef.hh:17
q
****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
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
Trajectory::position
virtual HepPoint3D position(double) const =0
Trajectory::delDirect
virtual Hep3Vector delDirect(double) const =0
Trajectory::direction
virtual Hep3Vector direction(double) const =0
Trajectory::curvature
virtual double curvature(double) const =0
TrkPocaBase::TrkPocaBase
TrkPocaBase(double flt1, double flt2, double precision)
Definition TrkPocaBase.cxx:25
TrkPocaXY::TrkPocaXY
TrkPocaXY(const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
Definition TrkPocaXY.cxx:23
m2
double double * m2
Definition qcdloop1.h:83
m1
double * m1
Definition qcdloop1.h:83