ATLAS Offline Software
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
MuonCalib::DCSLFitter Class Reference

Straight line fitter for drift circles. More...

#include <DCSLFitter.h>

Inheritance diagram for MuonCalib::DCSLFitter:
Collaboration diagram for MuonCalib::DCSLFitter:

Public Types

typedef std::vector< unsigned int > HitSelection

Public Member Functions

 DCSLFitter ()
bool fit (MuonCalibSegment &seg) const
 fit using all hits
bool fit (MuonCalibSegment &seg, HitSelection selection) const
 fit subset of the hits.
void printLevel (int level)
 set print level

Private Member Functions

double getY (const Amg::Vector3D &p) const
double getZ (const Amg::Vector3D &p) const
Amg::Vector3D getVec (double x, double y, double z) const

Private Attributes

bool m_debug

Detailed Description

Straight line fitter for drift circles.

Author
Niels.nosp@m..Van.nosp@m..Eldi.nosp@m.k@ce.nosp@m.rn.ch

Definition at line 28 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h.

Member Typedef Documentation

◆ HitSelection

typedef std::vector<unsigned int> MuonCalib::IMdtSegmentFitter::HitSelection
inherited

Definition at line 32 of file IMdtSegmentFitter.h.

Constructor & Destructor Documentation

◆ DCSLFitter()

MuonCalib::DCSLFitter::DCSLFitter ( )
inline

Definition at line 30 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h.

30: m_debug(false) {}

Member Function Documentation

◆ fit() [1/2]

bool MuonCalib::DCSLFitter::fit ( MuonCalibSegment & seg) const
virtual

fit using all hits

Implements MuonCalib::IMdtSegmentFitter.

Definition at line 17 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx.

17 {
18 // select all hits
19 HitSelection selection(seg.mdtHitsOnTrack(), 0);
20
21 // call fit function
22 return fit(seg, std::move(selection));
23 }
MuonCalib::DCSLFitter::fit
bool fit(MuonCalibSegment &seg) const
fit using all hits
Definition MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx:17
MuonCalib::IMdtSegmentFitter::HitSelection
std::vector< unsigned int > HitSelection
Definition IMdtSegmentFitter.h:32
selection
const std::string selection
Definition fbtTestBasics.cxx:75

◆ fit() [2/2]

bool MuonCalib::DCSLFitter::fit ( MuonCalibSegment & seg,
HitSelection selection ) const
virtual

fit subset of the hits.

If the HitSelection vector contains a 0 for a given hit the hit is used else the hit is not included in the fit. The size of the HitSelection vector should be equal to the number of hits on track else no fit is performed.

Implements MuonCalib::IMdtSegmentFitter.

Definition at line 25 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx.

25 {
26 if (m_debug) {
27 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
28 log << MSG::DEBUG << "New seg: " << endmsg;
29 }
30
31 int N = seg.mdtHitsOnTrack();
32
33 if (N < 2) { return false; }
34
35 if ((int)selection.size() != N) {
36 selection.clear();
37 selection.assign(N, 0);
38 } else {
39 int used(0);
40 for (int i = 0; i < N; ++i) {
41 if (selection[i] == 0) ++used;
42 }
43 if (used < 2) {
44 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
45 log << MSG::WARNING << "TO FEW HITS SELECTED" << endmsg;
46 return false;
47 }
48 }
49
50 Amg::Vector3D pos = seg.position();
51 Amg::Vector3D dir = seg.direction();
52
53 double S(0), Sz(0), Sy(0);
54 double Zc(0), Yc(0);
55
56 std::vector<double> y(N), z(N), r(N), w(N);
57 {
58 int ii(0);
59 for (const MuonCalibSegment::MdtHitPtr& hit_ptr : seg.mdtHOT()) {
60 const MdtCalibHitBase& h = *hit_ptr;
61
62 y[ii] = getY(h.localPosition());
63 z[ii] = getZ(h.localPosition());
64 r[ii] = std::abs(h.driftRadius());
65 if (h.sigma2DriftRadius() > 0)
66 w[ii] = 1. / (h.sigma2DriftRadius());
67 else
68 w[ii] = 0.;
69 if (r[ii] < 0) {
70 r[ii] = 0.;
71 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
72 log << MSG::WARNING << "<Negative r> " << r[ii] << endmsg;
73 }
74 if (m_debug) {
75 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
76 log << MSG::DEBUG << "DC: (" << y[ii] << "," << z[ii] << ") R = " << r[ii] << " W " << w[ii] << endmsg;
77 }
78
79 if (selection[ii]) {
80 ++ii;
81 continue;
82 }
83 S += w[ii];
84 Sz += w[ii] * z[ii];
85 Sy += w[ii] * y[ii];
86 ++ii;
87 }
88 }
89 if (S == 0)[[unlikely]]{
90 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
91 log << MSG::ERROR << "Divisor is zero" << endmsg;
92 return false;
93 }
94 Zc = Sz / S;
95 Yc = Sy / S;
96
97 if (m_debug) {
98 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
99 log << MSG::DEBUG << "Yc " << Yc << " Zc " << Zc << endmsg;
100 }
101 //
102 // shift hits
103 //
104 Sy = 0;
105 Sz = 0;
106 double Szz(0), Syy(0), Szy(0), Syyzz(0);
107 std::vector<double> rw(N);
108 std::vector<double> ryw(N);
109 std::vector<double> rzw(N);
110
111 for (int i = 0; i < N; ++i) {
112 y[i] -= Yc;
113 z[i] -= Zc;
114
115 if (selection[i]) continue;
116
117 rw[i] = r[i] * w[i];
118 ryw[i] = rw[i] * y[i];
119 rzw[i] = rw[i] * z[i];
120
121 Szz += z[i] * z[i] * w[i];
122 Syy += y[i] * y[i] * w[i];
123 Szy += y[i] * z[i] * w[i];
124 Syyzz += (y[i] - z[i]) * (y[i] + z[i]) * w[i];
125 }
126
127 if (m_debug) {
128 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
129 log << MSG::DEBUG << " Szz " << Szz << " Syy " << Syy << " Szy " << Szy << " Syyzz " << Syyzz << endmsg;
130 }
131
132 int count(0);
133 double R(0), Ry(0), Rz(0);
134 double Att(0);
135 double Add = S;
136 double Bt(0);
137 double Bd(0);
138 double Stt(0);
139 double Sd(0);
140
141 double cosin = getZ(dir) / dir.mag();
142 double sinus = getY(dir) / dir.mag();
143
144 if (m_debug) {
145 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
146 log << MSG::DEBUG << "cos " << cosin << " sin " << sinus << endmsg;
147 }
148
149 // make sure 0 <= theta < PI
150 if (sinus < 0.0) {
151 if (m_debug) {
152 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
153 log << MSG::DEBUG << "sinus=" << sinus << endmsg;
154 }
155 sinus = -sinus;
156 cosin = -cosin;
157 } else if (sinus == 0.0 && cosin < 0.0) {
158 if (m_debug) {
159 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
160 log << MSG::DEBUG << "sinus == 0.0 && cosin < 0.0" << endmsg;
161 }
162 cosin = -cosin;
163 }
164 //
165 // calculate shift
166 //
167 double d = -(getZ(pos) - Zc) * sinus + (getY(pos) - Yc) * cosin;
168 double theta = std::acos(cosin);
169 if (m_debug) {
170 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
171 log << MSG::DEBUG << "____________INITIAL VALUES________________" << count << endmsg;
172 log << MSG::DEBUG << "Theta " << theta << " d " << d << count << endmsg;
173 }
174
175 while (count < 100) {
176 if (m_debug) {
177 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
178 log << MSG::DEBUG << "____________NEW ITERATION________________" << count << endmsg;
179 }
180 R = 0;
181 Ry = 0;
182 Rz = 0;
183 for (int i = 0; i < N; ++i) {
184 if (selection[i]) continue;
185
186 double dist = y[i] * cosin - z[i] * sinus;
187 if (dist > d) {
188 R -= rw[i];
189 Ry -= ryw[i];
190 Rz -= rzw[i];
191 if (m_debug) {
192 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
193 log << MSG::DEBUG << " < - > " << dist - d << " r " << r[i] << endmsg;
194 }
195 } else {
196 R += rw[i];
197 Ry += ryw[i];
198 Rz += rzw[i];
199 if (m_debug) {
200 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
201 log << MSG::DEBUG << " < + > " << dist - d << " r " << r[i] << endmsg;
202 }
203 }
204 }
205 Att = Syy + cosin * (2 * sinus * Szy - cosin * Syyzz);
206 Bt = -Szy + cosin * (sinus * Syyzz + 2 * cosin * Szy + Rz) + sinus * Ry;
207 Bd = -S * d + R;
208 if (Att == 0) {
209 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
210 log << MSG::WARNING << "NewtonSLDCFitter ZERO Determinant" << endmsg;
211 return false;
212 }
213 theta += Bt / Att;
214 if (theta < 0.0) theta += M_PI;
215 if (theta >= M_PI) theta -= M_PI;
216 cosin = std::cos(theta);
217 sinus = std::sqrt(1 - cosin * cosin);
218 d = R / S;
219 if (m_debug) {
220 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
221 log << MSG::DEBUG << "R " << R << " Ry " << Ry << " Rz " << Rz << endmsg;
222 log << MSG::DEBUG << "Att " << Att << " Add " << Add << " Bt " << Bt << " Bd " << Bd << endmsg;
223 log << MSG::DEBUG << "dTheta " << Bt / Att << " dD " << Bd / Add << endmsg;
224 }
225 if (std::abs(Bt / Att) < 0.001 && std::abs(Bd / Add) < 0.001) {
226 Stt = std::sqrt(1 / Att);
227 Sd = std::sqrt(1 / Add);
228
229 if (m_debug) {
230 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
231 log << MSG::DEBUG << "Fit converged after " << count << " iterations " << endmsg;
232 log << MSG::DEBUG << "Theta " << theta << "d " << d << endmsg;
233 log << MSG::DEBUG << "Errors: theta " << Stt << " d " << Sd << endmsg;
234 }
235
236 seg.setErrors(Sd, Stt);
237
238 break;
239 }
240 ++count;
241 }
242 if (m_debug) {
243 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
244 log << MSG::DEBUG << "Calculating chi2" << endmsg;
245 }
246
247 std::vector<double> yl(N);
248 std::vector<double> dyl(N);
249 std::vector<double> res(N);
250 double chi2(0);
251 if (m_debug) {
252 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
253 log << MSG::DEBUG << "contributions to chi2: " << endmsg;
254 }
255
256 // calculate predicted hit positions from track parameters
257 for (int i = 0; i < N; ++i) {
258 yl[i] = cosin * y[i] - sinus * z[i] - d;
259 double dth = -(sinus * y[i] + cosin * z[i]) * Stt;
260 dyl[i] = std::hypot(dth, Sd);
261 res[i] = std::abs(yl[i]) - r[i];
262 if (m_debug) {
263 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
264 log << MSG::DEBUG << " r_track " << yl[i] << " dr " << dyl[i] << " r_rt " << r[i] << " res " << res[i] << " pull "
265 << res[i] * res[i] * w[i] << endmsg;
266 }
267 // skip hits that are not used in fit
268 if (selection[i]) continue;
269
270 chi2 += res[i] * res[i] * w[i];
271 }
272
273 if (m_debug) {
274 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
275 log << MSG::DEBUG << "Fit complete: Chi2 tof " << chi2 / (N - 2) << " " << !(chi2 / (N - 2) > 5) << endmsg;
276 }
277 if (chi2 / (N - 2) > 5) {
278 if (m_debug) {
279 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
280 log << MSG::DEBUG << "_______NOT GOOD " << endmsg;
281 }
282 }
283
284 if (m_debug) {
285 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
286 log << MSG::DEBUG << "Transforming back to real world" << endmsg;
287 }
288
289 Amg::Vector3D ndir = getVec(dir.x(), sinus, cosin);
290 Amg::Vector3D npos = getVec(pos.x(), Yc + cosin * d, Zc - sinus * d);
291
292 if (m_debug) {
293 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
294 log << MSG::DEBUG << "New line: position " << npos << " direction " << ndir << endmsg;
295 }
296
297 seg.set(chi2 / (N - 2), npos, ndir);
298
299 int i(0);
300 for (const MuonCalibSegment::MdtHitPtr& mdt_hit : seg.mdtHOT()) {
301 mdt_hit->setDistanceToTrack(yl[i], dyl[i]);
302 ++i;
303 }
304
305 if (m_debug) {
306 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
307 log << MSG::DEBUG << "fit done" << endmsg;
308 }
309 // tracking failed if position and directins are not real numbers
310 return (std::isfinite(ndir.y()) && std::isfinite(ndir.z()) && std::isfinite(npos.y()) && std::isfinite(npos.z()));
311 }
M_PI
#define M_PI
Definition ActiveFraction.h:14
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:54
res
std::pair< std::vector< unsigned int >, bool > res
Definition JetGroupProductTest.cxx:11
y
#define y
z
#define z
MuonCalib::MuonCalibSegment::MdtHitPtr
std::shared_ptr< MdtCalibHitBase > MdtHitPtr
typedef for a collection of MdtCalibHitBase s
Definition MuonCalibSegment.h:43
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:792
r
int r
Definition globals.cxx:22
count
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:148
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:48
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition AnalysisMisc.h:49
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition getMessageSvc.cxx:20
FtagLargeRLabel::Zc
@ Zc
Z -> cc, nC = 1 (partial containment).
Definition FtagLargeRJetLabelEnum.h:49
MuonR4::SegmentFit::N
constexpr std::size_t N
Definition SegmentFitterEventData.cxx:37
hist_file_dump.d
d
Definition hist_file_dump.py:142
unlikely
#define unlikely(x)
Definition pythonic_coracool.cxx:9

◆ getVec()

Amg::Vector3D MuonCalib::DCSLFitter::getVec ( double x,
double y,
double z ) const
inlineprivate

Definition at line 51 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h.

51{ return Amg::Vector3D(x, z, y); }
x
#define x

◆ getY()

double MuonCalib::DCSLFitter::getY ( const Amg::Vector3D & p) const
inlineprivate

Definition at line 49 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h.

49{ return p.z(); }

◆ getZ()

double MuonCalib::DCSLFitter::getZ ( const Amg::Vector3D & p) const
inlineprivate

Definition at line 50 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h.

50{ return p.y(); }

◆ printLevel()

void MuonCalib::DCSLFitter::printLevel ( int level)
virtual

set print level

Implements MuonCalib::IMdtSegmentFitter.

Definition at line 13 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx.

13 {
14 if (level > 0) m_debug = true;
15 }

Member Data Documentation

◆ m_debug

bool MuonCalib::DCSLFitter::m_debug
private

Definition at line 47 of file MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h.

The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.17.0