ATLAS Offline Software
MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "MdtCalibFitters/DCSLFitter.h"
6 
7 #include "AthenaKernel/getMessageSvc.h"
8 #include "GaudiKernel/MsgStream.h"
9 #include "cmath"
10 
11 namespace MuonCalib {
12 
13  void DCSLFitter::printLevel(int level) {
14  if (level > 0) m_debug = true;
15  }
16 
17  bool DCSLFitter::fit(MuonCalibSegment& seg) const {
18  // select all hits
19  HitSelection selection(seg.mdtHitsOnTrack(), 0);
20 
21  // call fit function
22  return fit(seg, selection);
23  }
24 
25  bool DCSLFitter::fit(MuonCalibSegment& seg, HitSelection selection) const {
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  }
307 
308 } // namespace MuonCalib
used
beamspotman.r
def r
Definition: beamspotman.py:676
getMessageSvc.h
singleton-like access to IMessageSvc via open function and helper
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::MuonCalibSegment
Definition: MuonCalibSegment.h:39
MuonCalib::DCSLFitter::fit
bool fit(MuonCalibSegment &seg) const
fit using all hits
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx:17
MuonCalib::DCSLFitter::printLevel
void printLevel(int level)
set print level
Definition: MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx:13
MuonCalib::MuonCalibSegment::position
const Amg::Vector3D & position() const
retrieve local position of segment (on station level)
Definition: MuonCalibSegment.cxx:186
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
checkFileSG.rw
rw
Definition: checkFileSG.py:124
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::MuonCalibSegment::direction
const Amg::Vector3D & direction() const
retrieve local direction of segment (on station level) retrieve the transformation from local chamber...
Definition: MuonCalibSegment.cxx:187
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
CscCalcPed - algorithm that finds the Cathode Strip Chamber pedestals from an RDO.
Definition: CscCalcPed.cxx:22
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
MuonCalib::MuonCalibSegment::mdtHitsOnTrack
unsigned int mdtHitsOnTrack() const
retrieve the number of MdtCalibHitBase s assigned to this segment
Definition: MuonCalibSegment.cxx:147
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::MuonCalibSegment::set
void set(double chi2, const Amg::Vector3D &pos, const Amg::Vector3D &dir)
Definition: MuonCalibSegment.cxx:129
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
MuonCalib::MdtCalibHitBase
Definition: MdtCalibHitBase.h:38
DEBUG
#define DEBUG
Definition: page_access.h:11
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
MuonCalib::MuonCalibSegment::setErrors
void setErrors(double error_dy0, double error_dtheta)
sets Local errors on MuonCalibSegment parameters
Definition: MuonCalibSegment.cxx:119
MuonCalib::MuonCalibSegment::mdtHOT
const MdtHitVec & mdtHOT() const
retrieve the full set of MdtCalibHitBase s assigned to this segment
Definition: MuonCalibSegment.cxx:148
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
DCSLFitter.h
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 Thu Jul 4 2024 21:13:54 for ATLAS Offline Software by   doxygen 1.8.18