gambit is hosted by Hepforge, IPPP Durham
 GAMBIT  v1.5.0-2191-ga4742ac a Global And Modular Bsm Inference Tool
PointsAndLines.hpp
Go to the documentation of this file.
1 #pragma once
2 #include <cmath>
3 #include <sstream>
4 #include <iostream>
5 #include <limits>
6
7 namespace Gambit {
8  namespace ColliderBit {
9
10
13  inline double addInQuad(const double& a, const double& b) { return sqrt(a*a + b*b); }
14
15
23  class P2 {
25
26  private:
27  double _x, _y;
29
31
32  public:
34  P2() : _x(0.), _y(0.) { }
36  P2(double x, double y) : _x(x), _y(y) { }
38  P2(const P2& other) : _x(other.getx()), _y(other.gety()) { }
40  P2& operator = (const P2& other) { _x = other.getx(); _y = other.gety(); return *this; }
42
44
45  public:
46  P2& setx(double x) { _x = x; return *this; }
47  P2& sety(double y) { _y = y; return *this; }
48  P2& setxy(double x, double y) { _x = x; _y = y; return *this; }
50
52
53  public:
54  double getx() const { return _x; }
55  double gety() const { return _y; }
61
63
64  public:
65  P2& operator += (const P2& other) { _x += other.getx(); _y += other.gety(); return *this; }
66  P2& operator -= (const P2& other) { _x -= other.getx(); _y -= other.gety(); return *this; }
67  P2& operator *= (double a) { _x *= a; _y *= a; return *this; }
68  P2& operator /= (double a) { _x /= a; _y /= a; return *this; }
70  };
71
73
74  inline P2 operator + (const P2& a, const P2& b) { P2 rtn = a; return rtn += b; }
75  inline P2 operator - (const P2& a, const P2& b) { P2 rtn = a; return rtn -= b; }
76  inline P2 operator * (const P2& a, double f) { P2 rtn = a; return rtn *= f; }
77  inline P2 operator * (double f, const P2& a) { P2 rtn = a; return rtn *= f; }
78  inline P2 operator / (const P2& a, double f) { P2 rtn = a; return rtn /= f; }
79
81  inline std::string to_str(const P2& p2) {
82  std::stringstream ss;
83  ss << "(" << p2.getx() << ", " << p2.gety() << ")";
84  return ss.str();
85  }
87  inline std::ostream& operator <<(std::ostream& ostr, const P2& p2) {
88  ostr << to_str(p2);
89  return ostr;
90  }
92
93
98  class LineSegment {
100
101  private:
102  P2 _p1, _p2, _diff;
104
106
107  public:
109  void init(double x1, double y1, double x2, double y2, double extendFrac=-1.) {
110  P2 rawpt1, rawpt2, extendEnds;
111  if (x1 > x2 or (x1 == x2 and y1 > y2)) {
112  rawpt2.setxy(x1, y1);
113  rawpt1.setxy(x2, y2);
114  } else {
115  rawpt1.setxy(x1, y1);
116  rawpt2.setxy(x2, y2);
117  }
118  if(extendFrac > 0.) {
119  extendEnds = (rawpt2 - rawpt1) * extendFrac;
120  _p2 = rawpt2 + extendEnds;
121  _p1 = rawpt1 - extendEnds;
122  } else {
123  _p2 = rawpt2;
124  _p1 = rawpt1;
125  }
126  _diff = _p2 - _p1;
127  }
128
130  void init(const P2& p1, const P2& p2, double extendFrac=-1.) {
131  init(p1.getx(), p1.gety(), p2.getx(), p2.gety(), extendFrac);
132  }
133
136  init(0., 0., 0., 0.);
137  }
138
140  LineSegment(double x1, double y1, double x2, double y2, double extendFrac=-1.) {
141  init(x1, y1, x2, y2, extendFrac);
142  }
143
145  LineSegment(const P2& p1, const P2& p2, double extendFrac=0.) {
146  init(p1, p2, extendFrac);
147  }
148
150  LineSegment(const LineSegment& other) {
151  _p1 = other.getp1();
152  _p2 = other.getp2();
153  _diff = _p2 - _p1;
154  }
155
158  _p1 = other.getp1();
159  _p2 = other.getp2();
160  _diff = _p2 - _p1;
161  return *this;
162  }
164
166
167  public:
168  const P2 getp1() const { return _p1; }
169  const P2 getp2() const { return _p2; }
170
172  double slope() const {
173  if (_p1.getx() == _p2.getx())
174  return std::numeric_limits<double>::infinity();
175  else
176  return _diff.gety() / _diff.getx();
177  }
179  double m() const { return slope(); }
180
182  double intercept() const {
183  if (_p1.getx() == _p2.getx())
184  return std::numeric_limits<double>::infinity();
185  else
186  return _p1.gety() - m() * _p1.getx();
187  }
189  double b() const { return intercept(); }
190
192  double abs() const { return _diff.abs(); }
194  double r() const { return abs(); }
196
198
199  public:
201  P2 intersectsAt(const LineSegment& other) const {
202  P2 result(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity());
203  double xintersect, yintersect;
204
205  // If the slopes are equal, they will never intersect
206  if (slope() == other.slope())
207  return result;
208
209  // If self or other has an infinite slope, change the intersect calculation
210  if (slope() == std::numeric_limits<double>::infinity()) {
211  xintersect = _p1.getx();
212  yintersect = other.m() * xintersect + other.b();
213  if (xintersect >= other.getp1().getx() and xintersect <= other.getp2().getx()
214  and yintersect >= _p1.gety() and yintersect <= _p2.gety())
215  result.setxy(xintersect, yintersect);
216  } else if (other.slope() == std::numeric_limits<double>::infinity()) {
217  xintersect = other.getp1().getx();
218  yintersect = m() * xintersect + b();
219  if (xintersect >= _p1.getx() and xintersect <= _p2.getx()
220  and yintersect >= other.getp1().gety() and yintersect <= other.getp2().gety())
221  result.setxy(xintersect, yintersect);
222  } else { // Regular intercept calculation
223  xintersect = (other.b() - b()) / (m() - other.m());
224  yintersect = m() * xintersect + b();
225  if (xintersect >= _p1.getx() and xintersect <= _p2.getx()
226  and xintersect >= other.getp1().getx() and xintersect <= other.getp2().getx())
227  result.setxy(xintersect, yintersect);
228  }
229  return result;
230  }
232  };
233
235
236  inline std::string to_str(const LineSegment& lineseg) {
238  std::stringstream ss;
239  ss << to_str(lineseg.getp1()) << " -> " << to_str(lineseg.getp2());
240  return ss.str();
241  }
243  inline std::ostream& operator <<(std::ostream& ostr, const LineSegment& lineseg) {
244  ostr << to_str(lineseg);
245  return ostr;
246  }
248
249  }
250 }
double r() const
Alias for abs.
P2 & operator+=(const P2 &other)
double b() const
Alias for intercept.
P2 operator-(const P2 &a, const P2 &b)
LineSegment(const LineSegment &other)
Copy constructor.
LineSegment(const P2 &p1, const P2 &p2, double extendFrac=0.)
Point constructor.
std::ostream & operator<<(std::ostream &os, const Cutflow &cf)
Print a Cutflow to a stream.
Definition: Cutflow.hpp:207
P2 operator+(const P2 &a, const P2 &b)
P2 operator*(const P2 &a, double f)
double m() const
Alias for slope.
P2(double x, double y)
Coordinate constructor.
LineSegment()
Default constructor.
P2 & operator-=(const P2 &other)
A simple container for a point on an xy plane.
double intercept() const
Get the y-intercept of the full line which contains this segment.
START_MODEL b
Definition: demo.hpp:270
void init(double x1, double y1, double x2, double y2, double extendFrac=-1.)
Coordinate initializer / recycler.
void init(const P2 &p1, const P2 &p2, double extendFrac=-1.)
Point initializer / recycler.
LineSegment(double x1, double y1, double x2, double y2, double extendFrac=-1.)
Coordinate constructor.
std::string to_str(const P2 &p2)
Make a string representation of the vector.
double slope() const
Get the slope of the line segment.
P2()
Default constructor.
P2 & setxy(double x, double y)
START_MODEL x2
Definition: demo.hpp:42
double r() const
Alias for abs.
A simple container for a line segment on an xy plane.
P2 intersectsAt(const LineSegment &other) const
Determine if this intersects with other and where.
#define f(x)
double abs() const
Get the length of the vector from the origin to this point.
P2 & operator=(const P2 &other)
Copy assignment operator.
double abs() const
Get the length of the vector from the origin to this point.
P2 operator/(const P2 &a, double f)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
P2(const P2 &other)
Copy constructor.