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