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