BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
HTrackParameter.cxx
Go to the documentation of this file.
1#include "VertexFit/HTrackParameter.h"
2#include "CLHEP/Units/PhysicalConstants.h"
3#include "VertexFit/BField.h"
4#include "VertexFit/WTrackParameter.h"
5using namespace CLHEP;
6const double alpha = -0.00299792458;
7// const double bField = 1.0;
8
9HTrackParameter::HTrackParameter() {
10 m_trackID = -1;
11 m_partID = -1;
12 m_hel = HepVector( 5, 0 );
13 m_eHel = HepSymMatrix( 5, 0 );
14}
15
16HTrackParameter::HTrackParameter( const HTrackParameter& htrk ) {
17 m_trackID = htrk.m_trackID;
18 m_partID = htrk.m_partID;
19 m_hel = htrk.m_hel;
20 m_eHel = htrk.m_eHel;
21}
22
23HTrackParameter& HTrackParameter ::operator=( const HTrackParameter& htrk ) {
24 m_trackID = htrk.m_trackID;
25 m_partID = htrk.m_partID;
26 m_hel = htrk.m_hel;
27 m_eHel = htrk.m_eHel;
28 return ( *this );
29}
30
31//
32// for DstMdcKalTrack --> HTrackParameter
33//
34
35HTrackParameter::HTrackParameter( const HepVector h, const HepSymMatrix eH, const int trackid,
36 const int partid ) {
37
38 m_hel = HepVector( 5, 0 );
39 m_eHel = HepSymMatrix( 5, 0 );
40
41 m_trackID = trackid;
42 m_partID = partid;
43 m_hel = h;
44 m_eHel = eH;
45}
46
47//
48// for DstMdcTrack --> HTrackParameter
49//
50
51HTrackParameter::HTrackParameter( const HepVector h, const double eH[], const int trackid,
52 const int partid ) {
53
54 m_hel = HepVector( 5, 0 );
55 m_eHel = HepSymMatrix( 5, 0 );
56
57 m_trackID = trackid;
58 m_partID = partid;
59 m_hel = h;
60 int k = 0;
61 for ( int i = 0; i < 5; i++ )
62 {
63 for ( int j = 0; j < 5; j++ )
64 {
65 m_eHel[i][j] = eH[k];
66 m_eHel[j][i] = eH[k];
67 k++;
68 }
69 }
70}
71
72//
73// WTrackParameter --> HTrackParameter
74//
75
76HTrackParameter::HTrackParameter( const WTrackParameter wtrk ) {
77
78 HepVector p( 3, 0 );
79 HepVector x( 3, 0 );
80 for ( int i = 0; i < 3; i++ )
81 {
82 p[i] = wtrk.w()[i];
83 x[i] = wtrk.w()[i + 4];
84 }
85
86 HTrackParameter htrk( wtrk.charge(), p, x );
87 HepMatrix A( 5, 3, 0 );
88 HepMatrix B( 5, 3, 0 );
89
90 A = htrk.dHdx( p, x );
91 B = htrk.dHdp( p, x );
92 HepMatrix T( 5, 7, 0 );
93 for ( int i = 0; i < 5; i++ )
94 {
95 for ( int j = 0; j < 3; j++ )
96 {
97 T[i][j] = B[i][j];
98 T[i][j + 4] = A[i][j];
99 }
100 }
101
102 HepSymMatrix eH( 5, 0 );
103 eH = ( wtrk.Ew() ).similarity( T );
104 int m_trackID = -1;
105 double mass = ( wtrk.p() ).m();
106 int m_partID = -1;
107 for ( int i = 0; i < 5; i++ )
108 {
109 if ( fabs( mass - xmass( i ) ) < 0.01 )
110 {
111 m_partID = i;
112 break;
113 }
114 }
115 // htrk.setTrackID(trackID);
116 // htrk.setPartID(partID);
117 htrk.setEHel( eH );
118 m_hel = htrk.hel();
119 m_eHel = htrk.eHel();
120}
121
122//
123// radius and centers of Helix in xy plane
124//
125
126double HTrackParameter::radius() const {
127 double bField = VertexFitBField::instance()->getBFieldZ( x3() );
128 int charge = m_hel[2] > 0 ? +1 : -1;
129 double a = alpha * bField * charge;
130 double pxy = charge / m_hel[2];
131 return ( pxy / a );
132}
133
134HepPoint3D HTrackParameter::center() const {
135 double bField = VertexFitBField::instance()->getBFieldZ( x3() );
136 int charge = m_hel[2] > 0 ? +1 : -1;
137 double a = alpha * bField * charge;
138 double pxy = charge / m_hel[2];
139 double rad = pxy / a;
140 double x = ( m_hel[0] + rad ) * cos( m_hel[1] );
141 double y = ( m_hel[0] + rad ) * sin( m_hel[1] );
142 double z = 0.0;
143 return HepPoint3D( x, y, z );
144}
145
146//
147// Helix parameter from instantaneous position and momentum
148//
149
150HTrackParameter::HTrackParameter( const int charge, const HepVector p, const HepVector vx ) {
151 m_trackID = -1;
152 m_partID = -1;
153 m_hel = HepVector( 5, 0 );
154 m_eHel = HepSymMatrix( 5, 0 );
155
156 HepPoint3D xyz( vx[0], vx[1], vx[2] );
157 double bField = VertexFitBField::instance()->getBFieldZ( xyz );
158 // std::cout << "bField in HTrackParameter(charge,p,vx) = " << bField << std::endl;
159
160 double a = alpha * bField * charge;
161 double px = p[0];
162 double py = p[1];
163 double pz = p[2];
164
165 double x = vx[0];
166 double y = vx[1];
167 double z = vx[2];
168
169 double pxy = sqrt( px * px + py * py );
170 double T = sqrt( ( px + a * y ) * ( px + a * y ) + ( py - a * x ) * ( py - a * x ) );
171 double J = ( x * px + y * py ) / T * a / pxy;
172
173 double drho = ( pxy / a ) - ( T / a );
174 double phi0 =
175 fmod( ( 4 * CLHEP::pi ) + atan2( 0 - px - a * y, py - a * x ), ( 2 * CLHEP::pi ) );
176 double kappa = charge / pxy;
177 double dz = z - ( pz / a ) * asin( J );
178 double lambda = pz / pxy;
179
180 m_hel[0] = drho;
181 m_hel[1] = phi0;
182 m_hel[2] = kappa;
183 m_hel[3] = dz;
184 m_hel[4] = lambda;
185}
186
187//
188// derivative Matrix, used in Kalman Vertex Fit A= dH/dx
189//
190
191HepMatrix HTrackParameter::dHdx( const HepVector p, const HepVector vx ) {
192 HepPoint3D xyz( vx[0], vx[1], vx[2] );
193 double bField = VertexFitBField::instance()->getBFieldZ( xyz );
194 // std::cout << "bField in dHdx(p,vx) = " << bField << std::endl;
195
196 double a = alpha * bField * charge();
197 double px = p[0];
198 double py = p[1];
199 double pz = p[2];
200
201 double x = vx[0];
202 double y = vx[1];
203 // double z = vx[2];
204
205 // double pxy = sqrt(px*px+py*py);
206 double T = sqrt( ( px + a * y ) * ( px + a * y ) + ( py - a * x ) * ( py - a * x ) );
207 // double J= (x*px+y*py)/T*a/pxy;
208
209 HepMatrix m_A( 5, 3, 0 );
210
211 m_A[0][0] = ( py - a * x ) / T;
212 m_A[0][1] = 0 - ( px + a * y ) / T;
213 m_A[1][0] = 0 - a * ( px + a * y ) / T / T;
214 m_A[1][1] = 0 - a * ( py - a * x ) / T / T;
215 m_A[3][0] = 0 - ( pz / T ) * ( px + a * y ) / T;
216 m_A[3][1] = 0 - ( pz / T ) * ( py - a * x ) / T;
217 m_A[3][2] = 1;
218
219 return m_A;
220}
221
222//
223// derivative Matrix, used in Kalman Vertex Fit B= dH/dp
224//
225
226HepMatrix HTrackParameter::dHdp( const HepVector p, const HepVector vx ) {
227 HepPoint3D xyz( vx[0], vx[1], vx[2] );
228 double bField = VertexFitBField::instance()->getBFieldZ( xyz );
229
230 double a = alpha * bField * charge();
231 double px = p[0];
232 double py = p[1];
233 double pz = p[2];
234
235 double x = vx[0];
236 double y = vx[1];
237 // double z = vx[2];
238
239 double pxy = sqrt( px * px + py * py );
240 double T = sqrt( ( px + a * y ) * ( px + a * y ) + ( py - a * x ) * ( py - a * x ) );
241 double J = ( x * px + y * py ) / T * a / pxy;
242
243 HepMatrix m_B( 5, 3, 0 );
244
245 m_B[0][0] = ( px / pxy - ( px + a * y ) / T ) / a;
246 m_B[0][1] = ( py / pxy - ( py - a * x ) / T ) / a;
247 m_B[1][0] = 0 - ( py - a * x ) / T / T;
248 m_B[1][1] = ( px + a * y ) / T / T;
249 m_B[2][0] = 0 - charge() * px / pxy / pxy / pxy;
250 m_B[2][1] = 0 - charge() * py / pxy / pxy / pxy;
251 m_B[3][0] = ( pz / a ) * ( py / pxy / pxy - ( py - a * x ) / T / T );
252 m_B[3][1] = 0 - ( pz / a ) * ( px / pxy / pxy - ( px + a * y ) / T / T );
253 m_B[3][2] = 0 - asin( J ) / a;
254 m_B[4][0] = 0 - ( px / pxy ) * ( pz / pxy ) / pxy;
255 m_B[4][1] = 0 - ( py / pxy ) * ( pz / pxy ) / pxy;
256 m_B[4][2] = 1 / pxy;
257
258 return m_B;
259}
260
261//
262// HTrackParameter --> WTrackParameter
263//
264
265WTrackParameter HTrackParameter::wTrack( const double mass ) const {
266
267 WTrackParameter wtrk;
268 wtrk.setCharge( charge() );
269 HepVector w( 7, 0 );
270 HepMatrix dWdh( 7, 5, 0 );
271
272 double ptrk = p3().rho();
273 double E = sqrt( ptrk * ptrk + mass * mass );
274 double px = p3().x();
275 double py = p3().y();
276 double pz = p3().z();
277 double x0 = x3().x();
278 double y0 = x3().y();
279 double z0 = x3().z();
280 w[0] = px;
281 w[1] = py;
282 w[2] = pz;
283 w[3] = E;
284 w[4] = x0;
285 w[5] = y0;
286 w[6] = z0;
287 wtrk.setW( w );
288 dWdh[0][1] = -py;
289 dWdh[0][2] = 0 - px / kappa();
290 dWdh[1][1] = px;
291 dWdh[1][2] = 0 - py / kappa();
292 dWdh[2][2] = 0 - pz / kappa();
293 dWdh[2][4] = charge() / kappa();
294 dWdh[3][2] = 0 - ( 1 + lambda() * lambda() ) / E / kappa() / kappa() / kappa();
295 dWdh[3][4] = lambda() / E / kappa() / kappa();
296 dWdh[4][0] = cos( phi0() );
297 dWdh[4][1] = 0 - y0;
298 dWdh[5][0] = sin( phi0() );
299 dWdh[5][1] = x0;
300 dWdh[6][3] = 1;
301 HepSymMatrix Ew( 7, 0 );
302 Ew = m_eHel.similarity( dWdh );
303 wtrk.setEw( Ew );
304 return wtrk;
305}
306
307WTrackParameter HTrackParameter::wTrack() const {
308 WTrackParameter wtrk;
309 if ( m_partID > -1 && m_partID < 5 )
310 {
311 double mass = xmass( m_partID );
312 wtrk = wTrack( mass );
313 }
314 return wtrk;
315}
316
317//
318// Utilities functions
319//
320
321double HTrackParameter::xmass( int n ) const {
322 double mass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
323 if ( n < 0 || n >= 5 ) return 0.0;
324 return mass[n];
325}
326
327//
328// Intersection position of two helice in xy plane
329//
330HepPoint3D HTrackParameter::positionTwoHelix( const HTrackParameter htrk ) const {
331 double bField1 = VertexFitBField::instance()->getBFieldZ( x3() );
332 double bField2 = VertexFitBField::instance()->getBFieldZ( htrk.x3() );
333
334 HepPoint3D pos1 = x3();
335 HepPoint3D pos2 = htrk.x3();
336 double rad1 = radius();
337 double rad2 = htrk.radius();
338 HepPoint3D xc1 = center();
339 HepPoint3D xc2 = htrk.center();
340 //
341 // requires the solution of
342 // a * y**2 + b*y + c = 0
343 // and
344 // x = (rt - 2*yt * y) / (2 * xt)
345 // where
346 // xt = xc2.x() - xc1.x()
347 // yt = xc2.y() - xc1.y()
348 // rt = rad1*rad1-rad2*rad2-xc1.perp2()+xc2.perp2()
349 // a = 1 + yt*yt/(xt*xt);
350 // b = 2*xc1.x()*yt/xt-yt*rt/(xt*xt)-2*xc1.y();
351 // c = rt*rt/(4*xt*xt)-xc1.x()*rt/xt-rad1*rad1+xc1.perp2();
352 //
353
354 double xt = xc2.x() - xc1.x();
355 double yt = xc2.y() - xc1.y();
356 if ( xt == 0 && yt == 0 ) return HepPoint3D( 999, 999, 999 );
357 double rt = rad1 * rad1 - rad2 * rad2 - xc1.perp2() + xc2.perp2();
358 double x1, y1, x2, y2;
359 if ( xt != 0 )
360 {
361 double a = 1 + yt * yt / ( xt * xt );
362 if ( a == 0 ) return HepPoint3D( 999, 999, 999 );
363 double b = 2 * xc1.x() * yt / xt - yt * rt / ( xt * xt ) - 2 * xc1.y();
364 double c = rt * rt / ( 4 * xt * xt ) - xc1.x() * rt / xt - rad1 * rad1 + xc1.perp2();
365 double d = b * b - 4 * a * c;
366 if ( d < 0 ) return HepPoint3D( 999, 999, 999 );
367 d = sqrt( d );
368 // two solution of intersection points
369
370 y1 = ( -b + d ) / ( 2 * a );
371 x1 = ( rt - 2 * yt * y1 ) / ( 2 * xt );
372
373 y2 = ( -b - d ) / ( 2 * a );
374 x2 = ( rt - 2 * yt * y2 ) / ( 2 * xt );
375 }
376 else
377 {
378 double a = 1 + xt * xt / ( yt * yt );
379 if ( a == 0 ) return HepPoint3D( 999, 999, 999 );
380 double b = 2 * xc1.y() * xt / yt - xt * rt / ( yt * yt ) - 2 * xc1.x();
381 double c = rt * rt / ( 4 * yt * yt ) - xc1.y() * rt / yt - rad1 * rad1 + xc1.perp2();
382 double d = b * b - 4 * a * c;
383 if ( d < 0 ) return HepPoint3D( 999, 999, 999 );
384 d = sqrt( d );
385 // two solution of intersection points
386
387 x1 = ( -b + d ) / ( 2 * a );
388 y1 = ( rt - 2 * xt * x1 ) / ( 2 * yt );
389
390 x2 = ( -b - d ) / ( 2 * a );
391 y2 = ( rt - 2 * xt * x2 ) / ( 2 * yt );
392 }
393
394 double z1[2], z2[2], J1[2], J2[2];
395
396 double a1 = alpha * bField1 * charge();
397 double a2 = alpha * bField2 * htrk.charge();
398 // z for solotion one
399 J1[0] = a1 * ( ( x1 - x3().x() ) * p3().x() + ( y1 - x3().y() ) * p3().y() ) / p3().perp2();
400 J1[1] = a2 *
401 ( ( x1 - htrk.x3().x() ) * htrk.p3().x() + ( y1 - htrk.x3().y() ) * htrk.p3().y() ) /
402 htrk.p3().perp2();
403 z1[0] = x3().z() + p3().z() / a1 * asin( J1[0] );
404 z1[1] = htrk.x3().z() + htrk.p3().z() / a2 * asin( J1[1] );
405 // z for solotion two
406 J2[0] = a1 * ( ( x2 - x3().x() ) * p3().x() + ( y2 - x3().y() ) * p3().y() ) / p3().perp2();
407 J2[1] = a2 *
408 ( ( x2 - htrk.x3().x() ) * htrk.p3().x() + ( y2 - htrk.x3().y() ) * htrk.p3().y() ) /
409 htrk.p3().perp2();
410 z2[0] = x3().z() + p3().z() / a2 * asin( J2[0] );
411 z2[1] = htrk.x3().z() + htrk.p3().z() / a2 * asin( J2[1] );
412
413 // take the solution if delta z is small
414
415 if ( fabs( z1[0] - z1[1] ) < fabs( z2[0] - z2[1] ) )
416 { return HepPoint3D( x1, y1, 0.5 * ( z1[0] + z1[1] ) ); }
417 else { return HepPoint3D( x2, y2, 0.5 * ( z2[0] + z2[1] ) ); }
418}
419
420//
421// intersection position of helix and Plane
422//
423
424HepPoint3D HTrackParameter::positionPlane( const double zp ) const {
425 return HepPoint3D( 999, 999, 999 );
426}
427
428//
429// intersection position of helix and a cylinder
430//
431
432HepPoint3D HTrackParameter::positionCylinder( const double R ) const {
433 return HepPoint3D( 999, 999, 999 );
434}
435
436//
437// intersection position of helix and a cone
438//
439
440HepPoint3D HTrackParameter::positionCone() const { return HepPoint3D( 999, 999, 999 ); }
441
442//
443// Minimum distance of two helice
444//
445
446double HTrackParameter::minDistanceTwoHelix( const HTrackParameter G, HepPoint3D& mpos ) {
447 int ifail;
448 double h = radius();
449 double g = G.radius();
450 double phiH0 = fmod( ( 4 * CLHEP::pi ) + atan2( p3().y(), p3().x() ), ( 2 * CLHEP::pi ) );
451 double phiG0 =
452 fmod( ( 4 * CLHEP::pi ) + atan2( G.p3().y(), G.p3().x() ), ( 2 * CLHEP::pi ) );
453 double lamH = lambda();
454 double lamG = G.lambda();
455 double a = x3().x() - G.x3().x() + g * sin( phiG0 ) - h * sin( phiH0 );
456 double b = x3().y() - G.x3().y() - g * cos( phiG0 ) + h * cos( phiH0 );
457 double c1 = h * lamH * lamH;
458 double c2 = 0 - g * lamG * lamG;
459 double d1 = 0 - g * lamG * lamH;
460 double d2 = h * lamG * lamH;
461 double e1 = lamH * ( x3().z() - G.x3().z() - h * phiH0 * lamH + g * phiG0 * lamG );
462 double e2 = lamG * ( x3().z() - G.x3().z() - h * phiH0 * lamH + g * phiG0 * lamG );
463
464 HepVector phiE( 2, 0 );
465 phiE[0] = phiH0;
466 phiE[1] = phiG0;
467 for ( int iter = 0; iter < 100; iter++ )
468 {
469 HepVector z( 2, 0 );
470 HepSymMatrix Omega( 2, 0 );
471 double phiH = phiE[0];
472 double phiG = phiE[1];
473 z[0] = cos( phiH ) * ( a - g * sin( phiG ) ) + sin( phiH ) * ( b + g * cos( phiG ) ) +
474 c1 * phiH + d1 * phiG + e1;
475 z[1] = cos( phiG ) * ( a + h * sin( phiH ) ) + sin( phiG ) * ( b - h * cos( phiH ) ) +
476 c2 * phiG + d2 * phiH + e2;
477 Omega[0][0] =
478 0 - sin( phiH ) * ( a - g * sin( phiG ) ) + cos( phiH ) * ( b + g * cos( phiG ) ) + c1;
479 Omega[0][1] = -g * cos( phiH ) * cos( phiG ) - g * sin( phiH ) * sin( phiG ) + d1;
480 Omega[1][0] = h * cos( phiH ) * cos( phiG ) + h * sin( phiG ) * sin( phiH ) + d2;
481 Omega[1][1] =
482 -sin( phiG ) * ( a + h * sin( phiH ) ) + cos( phiG ) * ( b - h * cos( phiH ) ) + c2;
483 HepVector phi( 2, 0 );
484 phi = phiE - Omega.inverse( ifail ) * z;
485 if ( ( fabs( phi[0] - phiE[0] ) < 1E-3 ) && ( fabs( phi[1] - phiE[1] ) < 1E-3 ) )
486 {
487 phiE = phi;
488 // std::cout << "number of iteration = " << iter <<std::endl;
489 break;
490 }
491 phiE = phi;
492 }
493
494 double phiH = phiE[0];
495 double phiG = phiE[1];
496 HepPoint3D posH = HepPoint3D( x3().x() + h * ( sin( phiH ) - sin( phiH0 ) ),
497 x3().y() + h * ( cos( phiH0 ) - cos( phiH ) ),
498 x3().z() + lamH * ( phiH - phiH0 ) );
499 HepPoint3D posG = HepPoint3D( G.x3().x() + g * ( sin( phiG ) - sin( phiG0 ) ),
500 G.x3().y() + g * ( cos( phiG0 ) - cos( phiG ) ),
501 G.x3().z() + lamG * ( phiG - phiG0 ) );
502 double dis = ( posH - posG ).rho();
503 mpos = 0.5 * ( posH + posG );
504 return dis;
505}
506//
507// Minimum distance of helix and a line
508//
mass
double mass
Definition CosmicGenerator.cxx:128
n
const Int_t n
Definition DataBase/tau_mode.c:57
e1
Double_t e1
Definition DataBase/tau_mode.c:7
e2
Double_t e2
Definition DataBase/tau_mode.c:7
alpha
double alpha
Definition EvtPhokharaDef.hh:17
w
double w
Definition EvtPhokharaDef.hh:28
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
xmass
const double xmass[5]
Definition Gam4pikp.cxx:35
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
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition InstallArea/x86_64-el9-gcc13-opt/include/VertexFit/HTrackParameter.h:23
HTrackParameter::positionCylinder
HepPoint3D positionCylinder(const double) const
Definition HTrackParameter.cxx:432
HTrackParameter::positionTwoHelix
HepPoint3D positionTwoHelix(const HTrackParameter) const
Definition HTrackParameter.cxx:330
HTrackParameter::dHdp
HepMatrix dHdp(const HepVector p, const HepVector x)
Definition HTrackParameter.cxx:226
HTrackParameter::minDistanceTwoHelix
double minDistanceTwoHelix(const HTrackParameter G, HepPoint3D &pos)
Definition HTrackParameter.cxx:446
HTrackParameter::dHdx
HepMatrix dHdx(const HepVector p, const HepVector x)
Definition HTrackParameter.cxx:191
HTrackParameter::xmass
double xmass(const int i) const
Definition HTrackParameter.cxx:321
HTrackParameter::wTrack
WTrackParameter wTrack() const
Definition HTrackParameter.cxx:307
HTrackParameter::positionPlane
HepPoint3D positionPlane(const double) const
Definition HTrackParameter.cxx:424
HTrackParameter::center
HepPoint3D center() const
Definition HTrackParameter.cxx:134
HTrackParameter::radius
double radius() const
Definition HTrackParameter.cxx:126
HTrackParameter::positionCone
HepPoint3D positionCone() const
Definition HTrackParameter.cxx:440
VertexFitBField::getBFieldZ
double getBFieldZ(const HepPoint3D &vtx)
Definition Analysis/VertexFit/src/BField.cxx:51