Segment.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 ///
3 /// \file Segment.cxx
4 /// \brief Segments are a basic unit for finding vertices
5 /// \author messier@indiana.edu
6 ///
7 ////////////////////////////////////////////////////////////////////////
9
10 #include <cmath>
11 #include <iostream>
12
13 namespace vdt
14 {
15
16  Segment::Segment(bool flip) :
17  fFlipped(flip),
18  fX1(0),
19  fY1(0),
20  fX2(0),
21  fY2(0),
22  fChi2Line0(-1),
23  fSw(0),
24  fSwx(0),
25  fSwx2(0),
26  fSwy(0),
27  fSwy2(0),
28  fSwxy(0)
29  { }
30
31  //......................................................................
32
33  void Segment::Flip(bool flip)
34  {
35  if (flip!=fFlipped) {
36  fFlipped = flip;
37  if (fSw!=0) {
38  std::cerr << __FILE__ << ":" << __LINE__
39  << " Flip called after segment initialization. "
40  << "Will lead to incorrect results!" << std::endl;
41  abort();
42  }
43  }
44  }
45
46  //......................................................................
47
48  void Segment::SetEndPoint1(double x, double y, double resin, double resex)
49  {
50  if (!fFlipped) { fX1 = x; fY1 = y; }
51  else { fX1 = y; fY1 = x; }
52
53  fdX = fX2-fX1;
54  fdY = fY2-fY1;
55
56  fResIN1sqr = resin*resin;
57  fResEX1 = resex;
58
59  double mag = sqrt(fdX*fdX+fdY*fdY);
60  fdX /= mag;
61  fdY /= mag;
62  }
63
64  //......................................................................
65
66  void Segment::SetEndPoint2(double x, double y, double resin, double resex)
67  {
68  if (!fFlipped) { fX2 = x; fY2 = y; }
69  else { fX2 = y; fY2 = x; }
70
71  fdX = fX2-fX1;
72  fdY = fY2-fY1;
73
74  fResIN2sqr = resin*resin;
75  fResEX2 = resex;
76
77  double mag = sqrt(fdX*fdX+fdY*fdY);
78  fdX /= mag;
79  fdY /= mag;
80  }
81
82  //......................................................................
83
84  void Segment::AddPoint(double xx, double sigx, double yy, double sigy)
85  {
86  double x, y, w;
87  if (!fFlipped) { x = xx; y = yy; w = 1.0/(sigy*sigy); }
88  else { x = yy; y = xx; w = 1.0/(sigx*sigx); }
89
90  double wx = w*x;
91  double wy = w*y;
92
93  fSw += w;
94  fSwx += wx;
95  fSwx2 += wx*x;
96  fSwy += wy;
97  fSwy2 += wy*y;
98  fSwxy += wx*y;
99
100  //
101  // Changing the sums invalidates this chi^2 min
102  //
103  fChi2Line0 = -1;
104  }
105
106  //......................................................................
107
108  double Segment::BestM(double x, double y)
109  {
110  return ((fSw*y-fSwy)*x-fSwx*y+fSwxy)/((fSw*x-2*fSwx)*x+fSwx2);
111  }
112
113  //......................................................................
114
115  double Segment::Chi2(double x, double y)
116  {
117  return this->Chi2Line(x,y)+this->Chi2EndPoint(x,y);
118  }
119
120  //......................................................................
121
123  {
124  double d = fSw*fSwx2-fSwx*fSwx;
125  if (d==0.0) d = 1.0E-9;
126
127  double m = (fSw*fSwxy-fSwx*fSwy)/d;
128  double b = (fSwx2*fSwy-fSwx*fSwxy)/d;
129
130  fChi2Line0 =
131  m*m*fSwx2 + 2*b*m*fSwx - 2*m*fSwxy + b*b*fSw - 2*b*fSwy + fSwy2;
132  }
133
134  //......................................................................
135
136  double Segment::Chi2Line(double xx, double yy)
137  {
138  double x, y;
139  if (!fFlipped) { x = xx; y = yy; }
140  else { x = yy; y = xx; }
141
142  if (fChi2Line0<0) this->CalcChi2Line0();
143
144  double m = this->BestM(x,y);
145  double c =
146  ( ( (x*fSw - 2*fSwx)*x + fSwx2 )*m +
147  ( (-2*y*fSw + 2*fSwy)*x + 2*y*fSwx - 2*fSwxy ) )*m +
148  (y*fSw - 2*fSwy)*y + fSwy2;
149
150  return c-fChi2Line0;
151  }
152
153  //......................................................................
154
155  double Segment::Chi2EndPoint(double xx, double yy)
156  {
157  double x, y;
158  if (!fFlipped) { x = xx; y = yy; }
159  else { x = yy; y = xx; }
160
161  double chi1;
162  double d1 = this->GetD1(x,y);
163  if (d1>0.0) chi1 = 2.0*d1/fResEX1;
164  else chi1 = d1*d1/fResIN1sqr;
165
166  double chi2;
167  double d2 = this->GetD2(x,y);
168  if (d2>0.0) chi2 = 2.0*d2/fResEX2;
169  else chi2 = d2*d2/fResIN2sqr;
170
171  if (chi1<chi2) return chi1;
172  return chi2;
173  }
174
175  //......................................................................
176
177  double Segment::GetD1(double x, double y)
178  {
179  double dx10 = fX1-x;
180  double dy10 = fY1-y;
181  return (dx10*fdX)+(dy10*fdY);
182  }
183
184  //......................................................................
185
186  double Segment::GetD2(double x, double y)
187  {
188  double dx20 = fX2-x;
189  double dy20 = fY2-y;
190  return -(dx20*fdX)-(dy20*fdY);
191  }
192
193 } // end namespace vdt
194
195 ////////////////////////////////////////////////////////////////////////
double fSwy
Definition: Segment.h:69
double BestM(double x, double y)
Definition: Segment.cxx:108
double fY2
Definition: Segment.h:56
Double_t xx
Definition: macro.C:12
double fSwxy
Definition: Segment.h:71
void SetEndPoint1(double x, double y, double resin, double resex)
Definition: Segment.cxx:48
double fX1
Definition: Segment.h:55
T sqrt(T number)
Definition: d0nt_math.hpp:156
double GetD2(double x, double y)
Definition: Segment.cxx:186
OStream cerr
Definition: OStream.cxx:7
bool fFlipped
Definition: Segment.h:54
double fSw
Definition: Segment.h:66
double fSwy2
Definition: Segment.h:70
double chi2()
double fY1
Definition: Segment.h:55
void Flip(bool flip=true)
Definition: Segment.cxx:33
double fChi2Line0
Definition: Segment.h:62
double fResIN1sqr
Definition: Segment.h:58
double fResEX2
Definition: Segment.h:61
void AddPoint(double x, double sigx, double y, double sigy)
Definition: Segment.cxx:84
double Chi2(double x, double y)
Definition: Segment.cxx:115
double fResIN2sqr
Definition: Segment.h:60
void CalcChi2Line0()
Definition: Segment.cxx:122
double Chi2EndPoint(double x, double y)
Definition: Segment.cxx:155
Float_t d
Definition: plot.C:236
double fdY
Definition: Segment.h:57
double fResEX1
Definition: Segment.h:59
double fX2
Definition: Segment.h:56
double GetD1(double x, double y)
Definition: Segment.cxx:177
double fSwx
Definition: Segment.h:67
double fSwx2
Definition: Segment.h:68
double fdX
Definition: Segment.h:57
VertexDT finding algorithm.
const hit & b
Definition: hits.cxx:21
Segment(bool flipped=false)
Definition: Segment.cxx:16
Float_t w
Definition: plot.C:20
double Chi2Line(double x, double y)
Definition: Segment.cxx:136
void SetEndPoint2(double x, double y, double resin, double resex)
Definition: Segment.cxx:66