BOSS 7.0.9
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 "TrkBase/TrkPoca.h"
17#include "TrkBase/HelixTraj.h"
18#include "TrkBase/TrkLineTraj.h"
19#include "CLHEP/Vector/ThreeVector.h"
20using std::cout;
21using std::endl;
22
23TrkPocaXY::TrkPocaXY(const Trajectory& traj, double flt,
24 const HepPoint3D& pt, double prec)
25 : TrkPocaBase(flt,prec) , _docaxy(-9999.0) {
26// construct a line trajectory through the given point
27 Hep3Vector zaxis(0.0,0.0,1.0);
28// length is set to 200cm (shouldn't matter)
29// TrkLineTraj zline(pt, zaxis, 200.0);//yzhang del
30 TrkLineTraj zline(pt, zaxis, 2000.0);//yzhang change
31// create a 2-traj poca with this
32 double zlen(pt.z());
33 TrkPoca zlinepoca(traj,flt,zline,zlen,prec);
34// transfer over the data
35 _status = zlinepoca.status();
36 if(status().success()){
37 _flt1 = zlinepoca.flt1();
38 _flt2 = zlinepoca.flt2();
39 _docaxy = zlinepoca.doca(); // doca should be perpendicular to the zline so 2d by construction
40 }
41}
42
43TrkPocaXY::TrkPocaXY(const Trajectory& traj1,double flt1,
44 const Trajectory& traj2,double flt2,
45 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 &&
59 ( delta1 > prec || delta2 > prec ) ){
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 // convert to m*x+q
75 if(dir1[0] == 0){
76 m1 = 1.0e12;
77 }else{
78 m1 = dir1[1]/dir1[0];
79 }
80 q1 = pos1.y()-pos1.x()*m1;
81 }else{
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 if(dir2[0] == 0){
93 m2 = 1.0e12;
94 }else{
95 m2 = dir2[1]/dir2[0];
96 }
97 q2 = pos2.y()-pos2.x()*m2;
98 }else{
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 // First the line-line case
112
113 if(curv1==0 && curv2==0){
114 // do the intersection in 2d
115 interTwoLines(m1, q1, m2, q2, xint, yint);
116 }
117
118 // next do the two circle case
119
120 if(curv1 !=0 && curv2 !=0){
121 //There are four cases
122 double cdist = sqrt((xc1-xc2)*(xc1-xc2)+(yc1-yc2)*(yc1-yc2));
123 // a - coincident centers
124 if (fabs(cdist) < 1.e-12 ) {
125 //the algorithm fails because the points have all the same distance
127 "TrkPocaXY:: the two circles have the same center...");
128 return;
129 }
130 // b - Actual intersection
131 if ( (fabs(r1-r2) < cdist) && (cdist < r1+r2 ) ) {
132 interTwoCircles(xc1,yc1,r1,xc2,yc2,r2,xint1,yint1,xint2,yint2);
133 double dist1 = (xint1-pos1b.x())*(xint1-pos1b.x()) +
134 (yint1-pos1b.y())*(yint1-pos1b.y());
135 double dist2 = (xint2-pos1b.x())*(xint2-pos1b.x()) +
136 (yint2-pos1b.y())*(yint2-pos1b.y());
137
138 //choose closest to begining of track1
139 if(dist1<dist2){
140 xint = xint1;
141 yint = yint1;
142 } else {
143 xint = xint2;
144 yint = yint2;
145 }
146 }
147
148 // c - nested circles and d - separated circles
149
150 if ( (fabs(r1-r2) > cdist) || // nested circles
151 ( cdist > (r1+r2) )) { // separated circles
152 // line going through the centers of the two circles
153
154 double m = (yc1-yc2)/(xc1-xc2); // y = m*x+q
155 double q = yc1 - m*xc1;
156
157 // intersection points between the line and the two circles
158
159 double xint1, yint1, xint2, yint2, zOfDCrossT1;
160
161 interLineCircle(m, q, xc1, yc1, r1, xint1, yint1, xint2, yint2);
162
163 double xint3, yint3, xint4, yint4 ;
164 interLineCircle(m, q, xc2, yc2, r2, xint3, yint3, xint4, yint4);
165 if (fabs(r1-r2) > cdist) { // nested circles
166 absdoca = fabs(r1-r2)-cdist;
167#ifdef MDCPATREC_DEBUG
168 cout << "ErrMsg(debugging) doing nested circles in TrkPocaXY " << endl;
169#endif
170
171 double dist1_3 = pow((xint1-xint3),2.) + pow((yint1-yint3),2.);
172 double dist2_4 = pow((xint2-xint4),2.) + pow((yint2-yint4),2.);
173
174 if(dist1_3<dist2_4) {
175 xint = 0.5*(xint1+xint3);
176 yint = 0.5*(yint1+yint3);
177 zOfDCrossT1 = (xint3-xint1)*dir1[1]-(yint3-yint1)*dir1[0];
178 } else {
179 xint = 0.5*(xint2+xint4);
180 yint = 0.5*(yint2+yint4);
181 zOfDCrossT1 = (xint4-xint2)*dir1[1]-(yint4-yint2)*dir1[0];
182 }
183 _docaxy = absdoca;
184 if( zOfDCrossT1 > 0) _docaxy = -absdoca;
185 }
186 if( cdist > (r1+r2) ) { //separated circles
187 absdoca = cdist - (r1+r2);
188#ifdef MDCPATREC_DEBUG
189 cout << "ErrMsg(debugging) doing separated circles in TrkPocaXY"<< endl;
190#endif
191
192 double dist2_3 = pow((xint2-xint3),2.) + pow((yint2-yint3),2.);
193 double dist1_4 = pow((xint1-xint4),2.) + pow((yint1-yint4),2.);
194 if (dist2_3<dist1_4){
195 xint = 0.5*(xint2+xint3);
196 yint = 0.5*(yint2+yint3);
197 zOfDCrossT1 = (xint3-xint1)*dir1[1]-(yint3-yint1)*dir1[0];
198 } else {
199 xint = 0.5*(xint1+xint4);
200 yint = 0.5*(yint1+yint4);
201 zOfDCrossT1 = (xint4-xint2)*dir1[1]-(yint4-yint2)*dir1[0];
202 }
203 _docaxy = absdoca;
204 if( zOfDCrossT1 > 0) _docaxy = -absdoca;
205 }
206 }
207 }
208
209 // Now do the line-circle case
210
211 if((curv1 == 0 && curv2 !=0) || (curv1 != 0 && curv2 == 0)) {
212 // distance between the line and the circle center
213 HepVector dirT;
214 double m,q,r,xc,yc, zOfDCrossT1;
215 if(curv1==0) {m=m1; q=q1; r=r2; xc=xc2; yc=yc2; dirT=dir2;}
216 else {m=m2; q=q2; r=r1; xc=xc1; yc=yc1; dirT=dir1;}
217
218 double dist = fabs((m*xc-yc+q)/sqrt(1+m*m));
219 if(dist <= r) {
220
221 // the intersection points
222
223 interLineCircle(m,q,xc,yc,r, xint1, yint1, xint2, yint2);
224 double dist1 = (xint1-pos1b.x())*(xint1-pos1b.x()) +
225 (yint1-pos1b.y())*(yint1-pos1b.y());
226 double dist2 = (xint2-pos1b.x())*(xint2-pos1b.x()) +
227 (yint2-pos1b.y())*(yint2-pos1b.y());
228 //choose closest to the beginning of track 1
229 if(dist1<dist2){
230 xint = xint1;
231 yint = yint1;
232 } else {
233 xint = xint2;
234 yint = yint2;
235 }
236 } else { // no intersection points
237#ifdef MDCPATREC_DEBUG
238 cout << "ErrMsg(debugging) doing separated line-circle in TrkPocaXY"<<endl;
239#endif
240
241 // line going through the circle center and perpendicular to traj1
242
243 double mperp = -1./m;
244 double qperp = yc - mperp*xc;
245
246 // intersection between this line and the two trajectories
247
248 interLineCircle(mperp, qperp, xc, yc, r, xint1,yint1,xint2,yint2);
249
250 double xint3,yint3;
251 interTwoLines(m, q, mperp, qperp, xint3, yint3);
252
253 double dist1_3 = pow((xint1-xint3),2.) + pow((yint1-yint3),2.);
254 double dist2_3 = pow((xint2-xint3),2.) + pow((yint2-yint3),2.);
255 if (dist1_3<dist2_3) {
256 xint = 0.5*(xint3 + xint1);
257 yint = 0.5*(yint3 + yint1);
258 absdoca = sqrt((xint3-xint1)*(xint3-xint1)+
259 (yint3-yint1)*(yint3-yint1));
260 zOfDCrossT1 = (xint3-xint1)*dirT[1]-(yint3-yint1)*dirT[0];
261 }
262 else {
263 xint = 0.5*(xint3 + xint2);
264 yint = 0.5*(yint3 + yint2);
265 absdoca = sqrt((xint3-xint2)*(xint3-xint2)+
266 (yint3-yint2)*(yint3-yint2));
267 zOfDCrossT1 = (xint3-xint2)*dirT[1]-(yint3-yint2)*dirT[0];
268 }
269 _docaxy = absdoca;
270 if( zOfDCrossT1 > 0) _docaxy = -absdoca;
271 }
272 }
273
274 // get the flight lengths for the intersection point
275
276 HepPoint3D intpt( xint, yint, 0.0);
277 TrkPocaXY poca1(traj1,_flt1,intpt,prec);
278 TrkPocaXY poca2(traj2,_flt2,intpt,prec);
279 _status = poca2.status();
280 if(poca1.status().success() && poca2.status().success()) {
281 delta1 = fabs(_flt1 - poca1.flt1());
282 delta2 = fabs(_flt2 - poca2.flt1());
283 _flt1 = poca1.flt1();
284 _flt2 = poca2.flt1();
285 }else{
286#ifdef MDCPATREC_WARNING
287 cout << "ErrMsg(warning)" << " poca fialure " << endl;
288#endif
289 if(poca1.status().failure()) _status=poca1.status();
290 break;
291 }
292 niter++;
293 }
294#ifdef MDCPATREC_DEBUG
295 cout <<"ErrMsg(debugging)" <<"TrkPocaXY iterations = " << niter-1
296 << " flight lengths " << _flt1 <<" " << _flt2 << endl
297 << " deltas " << delta1 <<" " << delta2 << endl;
298#endif
299 if(niter == maxiter){
300#ifdef MDCPATREC_ROUTINE
301 cout << "ErrMsg(routine)" << "TrkPocaXY:: did not converge" << endl;
302#endif
303 _status.setFailure(13, "TrkPocaXY:: did not converge");
304 }
305}
306//------------------------------------------------------------------------
307void
308TrkPocaXY::interLineCircle(const double& m, const double& q,
309 const double& xc, const double& yc, const double& r,
310 double& xint1, double& yint1,
311 double& xint2, double& yint2)
312 //-------------------------------------------------------------------------
313{
314
315 double alpha = 1 + m*m;
316
317 double beta = -xc +m*(q-yc);
318
319 double gamma = xc*xc + (q-yc)*(q-yc) -r*r;
320
321 double Delta = beta*beta - alpha*gamma;
322
323 if (Delta < 0) {
324
326 "TrkPocaXY:: the line and the circle heve no intersections...");
327 return;
328
329 }
330 else if (fabs(Delta) <1.e-12) {
331
332 xint1 = -beta/alpha;
333 xint2 = xint1;
334
335 }
336 else {
337
338 double xPlus = -beta/alpha + sqrt(beta*beta - alpha*gamma)/alpha;
339 double xMinus = -beta/alpha - sqrt(beta*beta - alpha*gamma)/alpha;
340
341 if (xPlus > xMinus) {
342 xint1 = xMinus;
343 xint2 = xPlus;
344 }
345 else {
346 xint1 = xPlus;
347 xint2 = xMinus;
348 }
349 }
350 yint1 = m*xint1 + q;
351 yint2 = m*xint2 + q;
352
353 return;
354}
355 //------------------------------------------------------------------------
356void
357TrkPocaXY::interTwoCircles
358(const double& xc1, const double& yc1, const double& r1,
359 const double& xc2, const double& yc2, const double& r2,
360 double& xint1, double& yint1, double& xint2, double& yint2)
361 //-------------------------------------------------------------------------
362{
363
364 double A = (xc1*xc1 + yc1*yc1 - r1*r1) - (xc2*xc2 + yc2*yc2 - r2*r2);
365
366 double B = -xc1 + xc2;
367
368 double C = -yc1 + yc2;
369
370 double alpha = 1 + (B*B)/(C*C);
371
372 double beta = -xc1 + B/C*(yc1+A/(2*C));
373
374 double gamma = xc1*xc1 + (yc1+A/(2*C))*(yc1+A/(2*C)) - r1*r1;
375
376 double Delta = beta*beta - alpha*gamma;
377
378 if (Delta < 0) {
379
380 _status.setFailure(14, "TrkPocaXY:: the two circles have no intersections..\
381.");
382 return;
383
384 }
385 else if (fabs(Delta) <1.e-12) {
386
387 xint1 = -beta/alpha;
388 xint2 = xint1;
389
390 }
391 else {
392 double xPlus = -beta/alpha + sqrt(beta*beta - alpha*gamma)/alpha;
393 double xMinus = -beta/alpha - sqrt(beta*beta - alpha*gamma)/alpha;
394
395 if (xPlus > xMinus) {
396 xint1 = xMinus;
397 xint2 = xPlus;
398 }
399 else {
400 xint1 = xPlus;
401 xint2 = xMinus;
402 }
403
404 }
405
406 yint1 = -(A+2*B*xint1)/(2*C);
407 yint2 = -(A+2*B*xint2)/(2*C);
408
409
410 return;
411}
412
413//------------------------------------------------------------------------
414void
415TrkPocaXY::interTwoLines
416(const double& m1, const double& q1, const double& m2, const double& q2,
417 double& xint, double& yint)
418 //-------------------------------------------------------------------------
419{
420
421 if (fabs(m1-m2) < 1.e-12) { // parallel lines
422
423 //the algorithm fails because the points have all the same distance
424
425 _status.setFailure(13, "TrkPocaXY:: the two lines are parallel...");
426 return;
427 }
428 else { // the lines have an intersection point
429
430 xint = (q2-q1)/(m1-m2);
431 yint = m1*xint + q1;
432 }
433
434return;
435}
const double alpha
****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
***************************************************************************************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
virtual HepPoint3D position(double) const =0
virtual Hep3Vector delDirect(double) const =0
virtual Hep3Vector direction(double) const =0
virtual double curvature(double) const =0
int success() const
Definition: TrkErrCode.h:62
int failure() const
Definition: TrkErrCode.h:61
void setSuccess(int i, const char *str=0)
Definition: TrkErrCode.h:84
void setFailure(int i, const char *str=0)
Definition: TrkErrCode.h:79
TrkErrCode _status
Definition: TrkPocaBase.h:42
const TrkErrCode & status() const
Definition: TrkPocaBase.h:62
double flt2() const
Definition: TrkPocaBase.h:68
double _flt2
Definition: TrkPocaBase.h:41
double flt1() const
Definition: TrkPocaBase.h:65
double _flt1
Definition: TrkPocaBase.h:40
TrkPocaXY(const Trajectory &traj, double flt, const HepPoint3D &pt, double precision=1.0e-4)
Definition: TrkPocaXY.cxx:23
double doca() const
Definition: TrkPoca.h:56