BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
KalFitDoca.cxx
Go to the documentation of this file.
1
2
3#include "KalFitAlg/KalFitHitMdc.h"
4#include "KalFitAlg/KalFitWire.h"
5#include "KalFitAlg/helix/Helix.h"
6// #include "KalFitAlg/KalFitTrack.h"
7
8using namespace KalmanFit;
9// this approach function is used to calculate poca beteween from the helix and
10// the wire , well it also can deal with the wire sag.
11// this function transpalted from belle trasan packages.
12
13double Helix::approach( KalFitHitMdc& hit, bool doSagCorrection ) const {
14 //...Cal. dPhi to rotate...
15 // const TMDCWire & w = * link.wire();
16 HepPoint3D positionOnWire, positionOnTrack;
17 HepPoint3D pv = pivot();
18 // std::cout<<"the track pivot in approach is "<<pv<<std::endl;
19 HepVector Va = a();
20 // std::cout<<"the track parameters is "<<Va<<std::endl;
21 HepSymMatrix Ma = Ea();
22 // std::cout<<"the error matrix is "<<Ma<<std::endl;
23
24 Helix _helix( pv, Va, Ma );
25 //_helix.pivot(IP);
26 const KalFitWire& w = hit.wire();
27 Hep3Vector Wire = w.fwd() - w.bck();
28 // xyPosition(), returns middle position of a wire. z componet is 0.
29 // w.xyPosition(wp);
30 double wp[3];
31 wp[0] = 0.5 * ( w.fwd() + w.bck() ).x();
32 wp[1] = 0.5 * ( w.fwd() + w.bck() ).y();
33 wp[2] = 0.;
34 double wb[3];
35 // w.backwardPosition(wb);
36 wb[0] = w.bck().x();
37 wb[1] = w.bck().y();
38 wb[2] = w.bck().z();
39 double v[3];
40 v[0] = Wire.unit().x();
41 v[1] = Wire.unit().y();
42 v[2] = Wire.unit().z();
43 // std::cout<<"Wire.unit() is "<<Wire.unit()<<std::endl;
44
45 //...Sag correction...
46 /* if (doSagCorrection) {
47 HepVector3D dir = w.direction();
48 HepPoint3D xw(wp[0], wp[1], wp[2]);
49 HepPoint3D wireBackwardPosition(wb[0], wb[1], wb[2]);
50 w.wirePosition(link.positionOnTrack().z(),
51 xw,
52 wireBackwardPosition,
53 dir);
54 v[0] = dir.x();
55 v[1] = dir.y();
56 v[2] = dir.z();
57 wp[0] = xw.x();
58 wp[1] = xw.y();
59 wp[2] = xw.z();
60 wb[0] = wireBackwardPosition.x();
61 wb[1] = wireBackwardPosition.y();
62 wb[2] = wireBackwardPosition.z();
63 }
64 */
65 //...Cal. dPhi to rotate...
66 const HepPoint3D& xc = _helix.center();
67
68 // std::cout<<" helix center: "<<xc<<std::endl;
69
70 double xt[3];
71 _helix.x( 0., xt );
72 double x0 = -xc.x();
73 double y0 = -xc.y();
74 double x1 = wp[0] + x0;
75 double y1 = wp[1] + y0;
76 x0 += xt[0];
77 y0 += xt[1];
78 double dPhi = atan2( x0 * y1 - y0 * x1, x0 * x1 + y0 * y1 );
79
80 // std::cout<<"dPhi is "<<dPhi<<std::endl;
81
82 //...Setup...
83 double kappa = _helix.kappa();
84 double phi0 = _helix.phi0();
85
86 //...Axial case...
87 /* if (!w.stereo()) {
88 positionOnTrack = _helix.x(dPhi);
89 HepPoint3D x(wp[0], wp[1], wp[2]);
90 x.setZ(positionOnTrack.z());
91 positionOnWire = x;
92 //link.dPhi(dPhi);
93 std::cout<<" in axial wire : positionOnTrack is "<<positionOnTrack
94 <<" positionOnWire is "<<positionOnWire<<std::endl;
95 return (positionOnTrack - positionOnWire).mag();
96 }
97 */
98 double firstdfdphi = 0.;
99 static bool first = true;
100 if ( first )
101 {
102 // extern BelleTupleManager * BASF_Histogram;
103 // BelleTupleManager * m = BASF_Histogram;
104 // h_nTrial = m->histogram("TTrack::approach nTrial", 100, 0., 100.);
105 }
106 // #endif
107
108 //...Stereo case...
109 double rho = Helix::ConstantAlpha / kappa;
110 double tanLambda = _helix.tanl();
111 static HepPoint3D x;
112 double t_x[3];
113 double t_dXdPhi[3];
114 const double convergence = 1.0e-5;
115 double l;
116 unsigned nTrial = 0;
117 while ( nTrial < 100 )
118 {
119
120 // x = link.positionOnTrack(_helix->x(dPhi));
121 positionOnTrack = _helix.x( dPhi );
122 x = _helix.x( dPhi );
123 t_x[0] = x[0];
124 t_x[1] = x[1];
125 t_x[2] = x[2];
126
127 l = v[0] * t_x[0] + v[1] * t_x[1] + v[2] * t_x[2] - v[0] * wb[0] - v[1] * wb[1] -
128 v[2] * wb[2];
129
130 double rcosPhi = rho * cos( phi0 + dPhi );
131 double rsinPhi = rho * sin( phi0 + dPhi );
132 t_dXdPhi[0] = rsinPhi;
133 t_dXdPhi[1] = -rcosPhi;
134 t_dXdPhi[2] = -rho * tanLambda;
135
136 //...f = d(Distance) / d phi...
137 double t_d2Xd2Phi[2];
138 t_d2Xd2Phi[0] = rcosPhi;
139 t_d2Xd2Phi[1] = rsinPhi;
140
141 //...iw new...
142 double n[3];
143 n[0] = t_x[0] - wb[0];
144 n[1] = t_x[1] - wb[1];
145 n[2] = t_x[2] - wb[2];
146
147 double a[3];
148 a[0] = n[0] - l * v[0];
149 a[1] = n[1] - l * v[1];
150 a[2] = n[2] - l * v[2];
151 double dfdphi = a[0] * t_dXdPhi[0] + a[1] * t_dXdPhi[1] + a[2] * t_dXdPhi[2];
152
153 // #ifdef TRKRECO_DEBUG
154 if ( nTrial == 0 )
155 {
156 // break;
157 firstdfdphi = dfdphi;
158 }
159
160 //...Check bad case...
161 if ( nTrial > 3 )
162 { std::cout << " bad case, calculate approach ntrial = " << nTrial << std::endl; }
163 // #endif
164
165 //...Is it converged?...
166 if ( fabs( dfdphi ) < convergence ) break;
167
168 double dv = v[0] * t_dXdPhi[0] + v[1] * t_dXdPhi[1] + v[2] * t_dXdPhi[2];
169 double t0 =
170 t_dXdPhi[0] * t_dXdPhi[0] + t_dXdPhi[1] * t_dXdPhi[1] + t_dXdPhi[2] * t_dXdPhi[2];
171 double d2fd2phi = t0 - dv * dv + a[0] * t_d2Xd2Phi[0] + a[1] * t_d2Xd2Phi[1];
172 // + a[2] * t_d2Xd2Phi[2];
173
174 dPhi -= dfdphi / d2fd2phi;
175
176 // cout << "nTrial=" << nTrial << endl;
177 // cout << "iw f,df,dphi=" << dfdphi << "," << d2fd2phi << "," << dPhi << endl;
178
179 ++nTrial;
180 }
181 //...Cal. positions...
182 positionOnWire[0] = wb[0] + l * v[0];
183 positionOnWire[1] = wb[1] + l * v[1];
184 positionOnWire[2] = wb[2] + l * v[2];
185
186 // std::cout<<" positionOnTrack is "<<positionOnTrack
187 // <<" positionOnWire is "<<positionOnWire<<std::endl;
188 return ( positionOnTrack - positionOnWire ).mag();
189
190 // link.dPhi(dPhi);
191
192 // #ifdef TRKRECO_DEBUG
193 // h_nTrial->accumulate((float) nTrial + .5);
194 // #endif
195 // return nTrial;
196}
197
198double Helix::approach( HepPoint3D pfwd, HepPoint3D pbwd, bool doSagCorrection ) const {
199 // ...Cal. dPhi to rotate...
200 // const TMDCWire & w = * link.wire();
201 HepPoint3D positionOnWire, positionOnTrack;
202 HepPoint3D pv = pivot();
203 HepVector Va = a();
204 HepSymMatrix Ma = Ea();
205
206 Helix _helix( pv, Va, Ma );
207 Hep3Vector Wire;
208 Wire[0] = ( pfwd - pbwd ).x();
209 Wire[1] = ( pfwd - pbwd ).y();
210 Wire[2] = ( pfwd - pbwd ).z();
211 // xyPosition(), returns middle position of a wire. z componet is 0.
212 // w.xyPosition(wp);
213 double wp[3];
214 wp[0] = 0.5 * ( pfwd + pbwd ).x();
215 wp[1] = 0.5 * ( pfwd + pbwd ).y();
216 wp[2] = 0.;
217 double wb[3];
218 // w.backwardPosition(wb);
219 wb[0] = pbwd.x();
220 wb[1] = pbwd.y();
221 wb[2] = pbwd.z();
222 double v[3];
223 v[0] = Wire.unit().x();
224 v[1] = Wire.unit().y();
225 v[2] = Wire.unit().z();
226 // std::cout<<"Wire.unit() is "<<Wire.unit()<<std::endl;
227
228 // ...Sag correction...
229 /* if (doSagCorrection) {
230 HepVector3D dir = w.direction();
231 HepPoint3D xw(wp[0], wp[1], wp[2]);
232 HepPoint3D wireBackwardPosition(wb[0], wb[1], wb[2]);
233 w.wirePosition(link.positionOnTrack().z(),
234 xw,
235 wireBackwardPosition,
236 dir);
237 v[0] = dir.x();
238 v[1] = dir.y();
239 v[2] = dir.z();
240 wp[0] = xw.x();
241 wp[1] = xw.y();
242 wp[2] = xw.z();
243 wb[0] = wireBackwardPosition.x();
244 wb[1] = wireBackwardPosition.y();
245 wb[2] = wireBackwardPosition.z();
246 }
247 */
248 // ...Cal. dPhi to rotate...
249 const HepPoint3D& xc = _helix.center();
250 double xt[3];
251 _helix.x( 0., xt );
252 double x0 = -xc.x();
253 double y0 = -xc.y();
254 double x1 = wp[0] + x0;
255 double y1 = wp[1] + y0;
256 x0 += xt[0];
257 y0 += xt[1];
258 // std::cout<<" x0 is: "<<x0<<" y0 is: "<<y0<<std::endl;
259 // std::cout<<" x1 is: "<<x1<<" y1 is: "<<y1<<std::endl;
260 // std::cout<<" xt[0] is: "<<xt[0]<<" xt[1] is: "<<xt[1]<<std::endl;
261
262 double dPhi = atan2( x0 * y1 - y0 * x1, x0 * x1 + y0 * y1 );
263 // std::cout<<" x0 * y1 - y0 * x1 is: "<<(x0 * y1 - y0 * x1)<<std::endl;
264 // std::cout<<" x0 * x1 + y0 * y1 is: "<<(x0 * x1 + y0 * y1)<<std::endl;
265 // std::cout<<" before loop dPhi is "<<dPhi<<std::endl;
266 //...Setup...
267 double kappa = _helix.kappa();
268 double phi0 = _helix.phi0();
269
270 //...Axial case...
271 /* if (!w.stereo()) {
272 positionOnTrack = _helix.x(dPhi);
273 HepPoint3D x(wp[0], wp[1], wp[2]);
274 x.setZ(positionOnTrack.z());
275 positionOnWire = x;
276 //link.dPhi(dPhi);
277 std::cout<<" in axial wire : positionOnTrack is "<<positionOnTrack
278 <<" positionOnWire is "<<positionOnWire<<std::endl;
279 return (positionOnTrack - positionOnWire).mag();
280 }
281 */
282 double firstdfdphi = 0.;
283 static bool first = true;
284 if ( first )
285 {
286 // extern BelleTupleManager * BASF_Histogram;
287 // BelleTupleManager * m = BASF_Histogram;
288 // h_nTrial = m->histogram("TTrack::approach nTrial", 100, 0., 100.);
289 }
290 // #endif
291
292 //...Stereo case...
293 double rho = Helix::ConstantAlpha / kappa;
294 double tanLambda = _helix.tanl();
295 static HepPoint3D x;
296 double t_x[3];
297 double t_dXdPhi[3];
298 const double convergence = 1.0e-5;
299 double l;
300 unsigned nTrial = 0;
301 while ( nTrial < 100 )
302 {
303
304 // x = link.positionOnTrack(_helix->x(dPhi));
305 positionOnTrack = _helix.x( dPhi );
306 x = _helix.x( dPhi );
307 t_x[0] = x[0];
308 t_x[1] = x[1];
309 t_x[2] = x[2];
310
311 l = v[0] * t_x[0] + v[1] * t_x[1] + v[2] * t_x[2] - v[0] * wb[0] - v[1] * wb[1] -
312 v[2] * wb[2];
313
314 double rcosPhi = rho * cos( phi0 + dPhi );
315 double rsinPhi = rho * sin( phi0 + dPhi );
316 t_dXdPhi[0] = rsinPhi;
317 t_dXdPhi[1] = -rcosPhi;
318 t_dXdPhi[2] = -rho * tanLambda;
319
320 //...f = d(Distance) / d phi...
321 double t_d2Xd2Phi[2];
322 t_d2Xd2Phi[0] = rcosPhi;
323 t_d2Xd2Phi[1] = rsinPhi;
324
325 //...iw new...
326 double n[3];
327 n[0] = t_x[0] - wb[0];
328 n[1] = t_x[1] - wb[1];
329 n[2] = t_x[2] - wb[2];
330
331 double a[3];
332 a[0] = n[0] - l * v[0];
333 a[1] = n[1] - l * v[1];
334 a[2] = n[2] - l * v[2];
335 double dfdphi = a[0] * t_dXdPhi[0] + a[1] * t_dXdPhi[1] + a[2] * t_dXdPhi[2];
336
337 if ( nTrial == 0 )
338 {
339 // break;
340 firstdfdphi = dfdphi;
341 }
342
343 //...Check bad case...
344 if ( nTrial > 3 )
345 {
346 // std::cout<<" BAD CASE!!, calculate approach ntrial = "<<nTrial<< endl;
347 }
348 //...Is it converged?...
349 if ( fabs( dfdphi ) < convergence ) break;
350
351 double dv = v[0] * t_dXdPhi[0] + v[1] * t_dXdPhi[1] + v[2] * t_dXdPhi[2];
352 double t0 =
353 t_dXdPhi[0] * t_dXdPhi[0] + t_dXdPhi[1] * t_dXdPhi[1] + t_dXdPhi[2] * t_dXdPhi[2];
354 double d2fd2phi = t0 - dv * dv + a[0] * t_d2Xd2Phi[0] + a[1] * t_d2Xd2Phi[1];
355 // + a[2] * t_d2Xd2Phi[2];
356
357 dPhi -= dfdphi / d2fd2phi;
358 ++nTrial;
359 }
360 // std::cout<<" dPhi is: "<<dPhi<<std::endl;
361 //...Cal. positions...
362 positionOnWire[0] = wb[0] + l * v[0];
363 positionOnWire[1] = wb[1] + l * v[1];
364 positionOnWire[2] = wb[2] + l * v[2];
365
366 // std::cout<<"wb[0] is: "<<wb[0]<<" l is: "<<l<<" v[0] is: "<<v[0]<<std::endl;
367 // std::cout<<"wb[1] is: "<<wb[1]<<" v[1] is: "<<v[1]<<std::endl;
368 // std::cout<<"wb[2] is: "<<wb[2]<<" v[2] is: "<<v[2]<<std::endl;
369
370 // std::cout<<" positionOnTrack is "<<positionOnTrack
371 // <<" positionOnWire is "<<positionOnWire<<std::endl;
372
373 return ( positionOnTrack - positionOnWire ).mag();
374 // link.dPhi(dPhi);
375 // return nTrial;
376}
n
const Int_t n
Definition DataBase/tau_mode.c:57
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
w
double w
Definition EvtPhokharaDef.hh:28
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
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
KalmanFit::Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition InstallArea/x86_64-el9-gcc13-dbg/include/KalFitAlg/helix/Helix.h:183
KalmanFit::Helix::approach
double approach(KalFitHitMdc &hit, bool doSagCorrection) const
Definition KalFitDoca.cxx:13
KalmanFit::Helix::ConstantAlpha
static const double ConstantAlpha
Constant alpha for uniform field.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/KalFitAlg/helix/Helix.h:153
KalmanFit::Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/KalFitAlg/helix/Helix.h:205
KalmanFit::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Reconstruction/KalFitAlg/src/helix/Helix.cxx:121
KalmanFit::Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/KalFitAlg/helix/Helix.h:203
KalmanFit::Helix::Helix
Helix(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
Constructor with pivot, helix parameter a, and its error matrix.
Definition Reconstruction/KalFitAlg/src/helix/Helix.cxx:44
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/KalFitAlg/helix/Helix.h:185