47
48
49
50
51
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
63 if ( niter == 0 ) pos1b = pos1;
71 double m1,
m2, q1, q2;
72 double r1, r2, xc1, xc2, yc1, yc2, sr1, sr2;
73 if ( curv1 == 0 )
74 {
75
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
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;
110
111
112
113 if ( curv1 == 0 && curv2 == 0 )
114 {
115
116 interTwoLines(
m1, q1,
m2, q2, xint, yint );
117 }
118
119
120
121 if ( curv1 != 0 && curv2 != 0 )
122 {
123
124 double cdist = sqrt( ( xc1 - xc2 ) * ( xc1 - xc2 ) + ( yc1 - yc2 ) * ( yc1 - yc2 ) );
125
126 if ( fabs( cdist ) < 1.e-12 )
127 {
128
129 _status.setFailure( 12,
"TrkPocaXY:: the two circles have the same center..." );
130 return;
131 }
132
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
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
155
156 if ( ( fabs( r1 - r2 ) > cdist ) ||
157 ( cdist > ( r1 + r2 ) ) )
158 {
159
160
161 double m = ( yc1 - yc2 ) / ( xc1 - xc2 );
162 double q = yc1 - m * xc1;
163
164
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 {
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 {
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
225
226 if ( ( curv1 == 0 && curv2 != 0 ) || ( curv1 != 0 && curv2 == 0 ) )
227 {
228
229 HepVector dirT;
230 double m,
q, r, xc, yc, zOfDCrossT1;
231 if ( curv1 == 0 )
232 {
235 r = r2;
236 xc = xc2;
237 yc = yc2;
238 dirT = dir2;
239 }
240 else
241 {
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
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
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 {
275#ifdef MDCPATREC_DEBUG
276 cout << "ErrMsg(debugging) doing separated line-circle in TrkPocaXY" << endl;
277#endif
278
279
280
281 double mperp = -1. / m;
282 double qperp = yc - mperp * xc;
283
284
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
315
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 <<
" "
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
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)