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, 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 Zc = Sz / S;
90 Yc = Sy / S;
91
92 if (m_debug) {
93 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
94 log << MSG::DEBUG << "Yc " << Yc << " Zc " << Zc << endmsg;
95 }
96 //
97 // shift hits
98 //
99 Sy = 0;
100 Sz = 0;
101 double Szz(0), Syy(0), Szy(0), Syyzz(0);
102 std::vector<double> rw(N);
103 std::vector<double> ryw(N);
104 std::vector<double> rzw(N);
105
106 for (int i = 0; i < N; ++i) {
107 y[i] -= Yc;
108 z[i] -= Zc;
109
110 if (selection[i]) continue;
111
112 rw[i] = r[i] * w[i];
113 ryw[i] = rw[i] * y[i];
114 rzw[i] = rw[i] * z[i];
115
116 Szz += z[i] * z[i] * w[i];
117 Syy += y[i] * y[i] * w[i];
118 Szy += y[i] * z[i] * w[i];
119 Syyzz += (y[i] - z[i]) * (y[i] + z[i]) * w[i];
120 }
121
122 if (m_debug) {
123 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
124 log << MSG::DEBUG << " Szz " << Szz << " Syy " << Syy << " Szy " << Szy << " Syyzz " << Syyzz << endmsg;
125 }
126
127 int count(0);
128 double R(0), Ry(0), Rz(0);
129 double Att(0);
130 double Add = S;
131 double Bt(0);
132 double Bd(0);
133 double Stt(0);
134 double Sd(0);
135
136 double cosin = getZ(dir) / dir.mag();
137 double sinus = getY(dir) / dir.mag();
138
139 if (m_debug) {
140 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
141 log << MSG::DEBUG << "cos " << cosin << " sin " << sinus << endmsg;
142 }
143
144 // make sure 0 <= theta < PI
145 if (sinus < 0.0) {
146 if (m_debug) {
147 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
148 log << MSG::DEBUG << "sinus=" << sinus << endmsg;
149 }
150 sinus = -sinus;
151 cosin = -cosin;
152 } else if (sinus == 0.0 && cosin < 0.0) {
153 if (m_debug) {
154 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
155 log << MSG::DEBUG << "sinus == 0.0 && cosin < 0.0" << endmsg;
156 }
157 cosin = -cosin;
158 }
159 //
160 // calculate shift
161 //
162 double d = -(getZ(pos) - Zc) * sinus + (getY(pos) - Yc) * cosin;
163 double theta = std::acos(cosin);
164 if (m_debug) {
165 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
166 log << MSG::DEBUG << "____________INITIAL VALUES________________" << count << endmsg;
167 log << MSG::DEBUG << "Theta " << theta << " d " << d << count << endmsg;
168 }
169
170 while (count < 100) {
171 if (m_debug) {
172 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
173 log << MSG::DEBUG << "____________NEW ITERATION________________" << count << endmsg;
174 }
175 R = 0;
176 Ry = 0;
177 Rz = 0;
178 for (int i = 0; i < N; ++i) {
179 if (selection[i]) continue;
180
181 double dist = y[i] * cosin - z[i] * sinus;
182 if (dist > d) {
183 R -= rw[i];
184 Ry -= ryw[i];
185 Rz -= rzw[i];
186 if (m_debug) {
187 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
188 log << MSG::DEBUG << " < - > " << dist - d << " r " << r[i] << endmsg;
189 }
190 } else {
191 R += rw[i];
192 Ry += ryw[i];
193 Rz += rzw[i];
194 if (m_debug) {
195 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
196 log << MSG::DEBUG << " < + > " << dist - d << " r " << r[i] << endmsg;
197 }
198 }
199 }
200 Att = Syy + cosin * (2 * sinus * Szy - cosin * Syyzz);
201 Bt = -Szy + cosin * (sinus * Syyzz + 2 * cosin * Szy + Rz) + sinus * Ry;
202 Bd = -S * d + R;
203 if (Att == 0) {
204 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
205 log << MSG::WARNING << "NewtonSLDCFitter ZERO Determinant" << endmsg;
206 return false;
207 }
208 theta += Bt / Att;
209 if (theta < 0.0) theta += M_PI;
210 if (theta >= M_PI) theta -= M_PI;
211 cosin = std::cos(theta);
212 sinus = std::sqrt(1 - cosin * cosin);
213 d = R / S;
214 if (m_debug) {
215 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
216 log << MSG::DEBUG << "R " << R << " Ry " << Ry << " Rz " << Rz << endmsg;
217 log << MSG::DEBUG << "Att " << Att << " Add " << Add << " Bt " << Bt << " Bd " << Bd << endmsg;
218 log << MSG::DEBUG << "dTheta " << Bt / Att << " dD " << Bd / Add << endmsg;
219 }
220 if (std::abs(Bt / Att) < 0.001 && std::abs(Bd / Add) < 0.001) {
221 Stt = std::sqrt(1 / Att);
222 Sd = std::sqrt(1 / Add);
223
224 if (m_debug) {
225 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
226 log << MSG::DEBUG << "Fit converged after " << count << " iterations " << endmsg;
227 log << MSG::DEBUG << "Theta " << theta << "d " << d << endmsg;
228 log << MSG::DEBUG << "Errors: theta " << Stt << " d " << Sd << endmsg;
229 }
230
231 seg.setErrors(Sd, Stt);
232
233 break;
234 }
235 ++count;
236 }
237 if (m_debug) {
238 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
239 log << MSG::DEBUG << "Calculating chi2" << endmsg;
240 }
241
242 std::vector<double> yl(N);
243 std::vector<double> dyl(N);
244 std::vector<double> res(N);
245 double chi2(0);
246 if (m_debug) {
247 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
248 log << MSG::DEBUG << "contributions to chi2: " << endmsg;
249 }
250
251 // calculate predicted hit positions from track parameters
252 for (int i = 0; i < N; ++i) {
253 yl[i] = cosin * y[i] - sinus * z[i] - d;
254 double dth = -(sinus * y[i] + cosin * z[i]) * Stt;
255 dyl[i] = std::hypot(dth, Sd);
256 res[i] = std::abs(yl[i]) - r[i];
257 if (m_debug) {
258 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
259 log << MSG::DEBUG << " r_track " << yl[i] << " dr " << dyl[i] << " r_rt " << r[i] << " res " << res[i] << " pull "
260 << res[i] * res[i] * w[i] << endmsg;
261 }
262 // skip hits that are not used in fit
263 if (selection[i]) continue;
264
265 chi2 += res[i] * res[i] * w[i];
266 }
267
268 if (m_debug) {
269 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
270 log << MSG::DEBUG << "Fit complete: Chi2 tof " << chi2 / (N - 2) << " " << !(chi2 / (N - 2) > 5) << endmsg;
271 }
272 if (chi2 / (N - 2) > 5) {
273 if (m_debug) {
274 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
275 log << MSG::DEBUG << "_______NOT GOOD " << endmsg;
276 }
277 }
278
279 if (m_debug) {
280 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
281 log << MSG::DEBUG << "Transforming back to real world" << endmsg;
282 }
283
284 Amg::Vector3D ndir = getVec(dir.x(), sinus, cosin);
285 Amg::Vector3D npos = getVec(pos.x(), Yc + cosin * d, Zc - sinus * d);
286
287 if (m_debug) {
288 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
289 log << MSG::DEBUG << "New line: position " << npos << " direction " << ndir << endmsg;
290 }
291
292 seg.set(chi2 / (N - 2), npos, ndir);
293
294 int i(0);
295 for (const MuonCalibSegment::MdtHitPtr& mdt_hit : seg.mdtHOT()) {
296 mdt_hit->setDistanceToTrack(yl[i], dyl[i]);
297 ++i;
298 }
299
300 if (m_debug) {
301 MsgStream log(Athena::getMessageSvc(), "DCSLFitter");
302 log << MSG::DEBUG << "fit done" << endmsg;
303 }
304 // tracking failed if position and directins are not real numbers
305 return (std::isfinite(ndir.y()) && std::isfinite(ndir.z()) && std::isfinite(npos.y()) && std::isfinite(npos.z()));
306 }
M_PI
#define M_PI
Definition ActiveFraction.h:11
theta
Scalar theta() const
theta method
Definition AmgMatrixBasePlugin.h:75
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:63
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:525
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:146
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
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
hist_file_dump.d
d
Definition hist_file_dump.py:142

◆ 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.14.0