fit subset of the hits.
26 {
28 if (
log.level() <= MSG::DEBUG)
log << MSG::DEBUG <<
"fit() New seg:" <<
endmsg;
29
30 int N = seg.mdtHitsOnTrack();
31
32 if (N < 2) { return false; }
33
37 } else {
38 int used(0);
39 for (
int i = 0;
i <
N; ++
i) {
41 }
42 if (used < 2) {
43 log << MSG::WARNING <<
"fit() TO FEW HITS SELECTED" <<
endmsg;
44 return false;
45 }
46 }
47
49
50 double S(0), Sy(0), Sz(0),
Zc{0}, Yc{0};
51 std::vector<double>
y(N),
x(N),
z(N);
52 std::vector<double>
r(N),
sr(N),
w(N), rw(N);
53 int ii{0}, jj{0};
54 {
56 const MdtCalibHitBase&
h = *
hit;
57 y[ii] =
getY(
h.localPosition());
58 z[ii] =
getZ(
h.localPosition());
59 r[ii] = std::abs(
h.driftRadius());
60
61 if (
h.sigma2DriftRadius() > 0) {
62 w[ii] = 1. / (
h.sigma2DriftRadius());
63 sr[ii] = std::sqrt(
h.sigma2DriftRadius());
64 } else {
67 }
68 if (
log.level() <= MSG::DEBUG)
69 log << MSG::DEBUG <<
"fit() MuCCaFitter: (" <<
x[ii] <<
"," <<
y[ii] <<
") R = " <<
r[ii] <<
" W " <<
w[ii] <<
endmsg;
70 rw[ii] =
r[ii] *
w[ii];
72 ++jj;
73 continue;
74 }
78 ++ii;
79 ++jj;
80 }
81 }
83 log << MSG::ERROR <<
"Divisor S is zero."<<
endmsg;
84 return false;
85 }
88
89 std::unique_ptr<MuCCaFitterImplementation> Fitter = std::make_unique<MuCCaFitterImplementation>();
90 if (
log.level() <= MSG::DEBUG) {
91 for (
int i = 0;
i != ii;
i++) {
92 log << MSG::DEBUG <<
"fit() MuCCaFitter hits passed to computepam Z=" <<
x[
i] <<
" phi=" <<
y[
i] <<
" zzz=" <<
z[
i]
93 <<
" drift=" <<
r[
i] <<
" error=" <<
w[
i] <<
endmsg;
94 }
95 }
96 Fitter->Computeparam3(ii,
z,
y,
r, sr);
97 if (
log.level() <= MSG::DEBUG) {
98 log << MSG::DEBUG <<
"fit() MuCCaFitter computed track a=" << Fitter->get_a() <<
" da= " << Fitter->get_da()
99 << " b= " << Fitter->get_b() << " db= " << Fitter->get_db() << " corrab= " << Fitter->get_corrab()
100 <<
" chi2f= " << Fitter->get_chi2f() <<
endmsg;
101 }
102
103 double afit = Fitter->get_a();
104 double chi2f = Fitter->get_chi2f();
105
106 std::vector<double> dist(N), ddist(N);
108 double theta = std::atan(afit);
110 double sinus = std::sin(
theta);
111 double cosin = std::cos(
theta);
112 double d = -(
getZ(pos) -
Zc) * sinus + (
getY(pos) - Yc) * cosin;
113 if (
log.level() <= MSG::DEBUG) {
114 log << MSG::DEBUG <<
"fit() MuCCaFitter>>> theta= " <<
theta <<
" sinus= " << sinus <<
" cosin= " << cosin <<
endmsg;
115 log << MSG::DEBUG <<
"fit() MuCCaFitter>>> getZ( pos )= " <<
getZ(pos) <<
" getY( pos )= " <<
getY(pos) <<
" Zc= " <<
Zc
116 <<
" Yc= " << Yc <<
endmsg;
117 }
118 double Syy(0), Szy(0), Syyzz(0), Att(0);
120 for (
int i = 0;
i <
N; ++
i) {
122 Syy += (
y[
i] - Yc) * (
y[i] - Yc) *
w[
i];
123 Szy += (
y[
i] - Yc) * (
z[i] - Zc) *
w[
i];
124 Syyzz += ((
y[
i] - Yc) - (
z[
i] -
Zc)) * ((
y[i] - Yc) + (
z[i] - Zc)) *
w[
i];
125 dis = (
y[
i] - Yc) * cosin - (
z[i] - Zc) * sinus;
126 if (dis > d) {
128 } else {
130 }
131 }
132 Att = Syy + cosin * (2 * sinus * Szy - cosin * Syyzz);
134 if (
log.level() <= MSG::DEBUG){
135 log << MSG::DEBUG <<
"fit() MuCCaFitter>>> d= " <<
d <<
" R= " <<
R <<
" S= " <<
S <<
" Att= " << Att <<
endmsg;
136 }
138 log << MSG::ERROR <<
"Divisor Att is zero."<<
endmsg;
139 return false;
140 }
141 for (
int i = 0;
i <
N; ++
i) {
143 dist[
i] = cosin * (
y[
i] - Yc) - sinus * (
z[i] - Zc) -
d;
144 double dth = -(sinus * (
y[
i] - Yc) + cosin * (
z[
i] -
Zc)) * (std::sqrt(1. / Att));
145 ddist[
i] = std::sqrt(dth * dth + (1. / S));
146 }
147 if (
log.level() <= MSG::DEBUG)
log << MSG::DEBUG <<
"fit() Transforming back to real world" <<
endmsg;
150
151 if (
log.level() <= MSG::DEBUG)
152 log << MSG::DEBUG <<
"fit() New line: position " << npos <<
" direction " << ndir <<
" chi2f " << chi2f <<
endmsg;
153
154 seg.set(chi2f / (N - 2), npos, ndir);
155
158 hit_ptr->setDistanceToTrack(dist[i], ddist[i]);
160 }
161
162 if (
log.level() <= MSG::DEBUG)
log << MSG::DEBUG <<
"fit() fit done" <<
endmsg;
163 return true;
164 }
Scalar theta() const
theta method
Amg::Vector3D getVec(double x, double y, double z) const
double getZ(const Amg::Vector3D &p) const
double getY(const Amg::Vector3D &p) const
these methods are needed to change the reference frame between the local one of the hit and one used ...
std::shared_ptr< MdtCalibHitBase > MdtHitPtr
typedef for a collection of MdtCalibHitBase s
Eigen::Matrix< double, 3, 1 > Vector3D
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
IMessageSvc * getMessageSvc(bool quiet=false)
@ Zc
Z -> cc, nC = 1 (partial containment).