ATLAS Offline Software
Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
MuonCalib::DCSLFitter Class Reference

#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 More...
 
bool fit (MuonCalibSegment &seg, HitSelection selection) const
 fit subset of the hits. More...
 
void printLevel (int level)
 set print level More...
 

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  }

◆ 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  }

◆ 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); }

◆ 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:
used
beamspotman.r
def r
Definition: beamspotman.py:676
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
IDTPM::R
float R(const U &p)
Definition: TrackParametersHelper.h:101
MuonCalib::DCSLFitter::getZ
double getZ(const Amg::Vector3D &p) const
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h:50
hist_file_dump.d
d
Definition: hist_file_dump.py:137
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:71
M_PI
#define M_PI
Definition: ActiveFraction.h:11
MuonCalib::DCSLFitter::fit
bool fit(MuonCalibSegment &seg) const
fit using all hits
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx:17
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
checkFileSG.rw
rw
Definition: checkFileSG.py:124
x
#define x
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition: getMessageSvc.cxx:20
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
python.iconfTool.models.loaders.level
level
Definition: loaders.py:20
MuonCalib::DCSLFitter::getY
double getY(const Amg::Vector3D &p) const
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h:49
lumiFormat.i
int i
Definition: lumiFormat.py:92
z
#define z
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:14
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
MuonCalib::DCSLFitter::m_debug
bool m_debug
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h:47
selection
std::string selection
Definition: fbtTestBasics.cxx:73
beamspotman.dir
string dir
Definition: beamspotman.py:623
python.grid.Add
def Add(name)
Definition: grid.py:41
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
MuonCalib::DCSLFitter::getVec
Amg::Vector3D getVec(double x, double y, double z) const
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h:51
y
#define y
h
DEBUG
#define DEBUG
Definition: page_access.h:11
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
MuonCalib::MuonCalibSegment::MdtHitPtr
std::shared_ptr< MdtCalibHitBase > MdtHitPtr
typedef for a collection of MdtCalibHitBase s
Definition: MuonCalibSegment.h:44
MuonCalib::IMdtSegmentFitter::HitSelection
std::vector< unsigned int > HitSelection
Definition: IMdtSegmentFitter.h:32
Generated on Tue Jul 2 2024 22:23:09 for ATLAS Offline Software by   doxygen 1.8.18