Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DHoughFind.h
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DHoughFind.h
4 // Created: Wed Oct 31 05:59:26 EDT 2007
5 // Creator: davidl (on Darwin Amelia.local 8.10.1 i386)
6 //
7 
8 #ifndef _DHoughFind_
9 #define _DHoughFind_
10 
11 #include <vector>
12 using std::vector;
13 
14 #include <JANA/jerror.h>
15 #include <JANA/JObject.h>
16 using namespace jana;
17 
18 #include <TH2.h>
19 
20 #include <DVector2.h>
21 
22 class DHoughFind:public JObject{
23  public:
24  DHoughFind();
25  DHoughFind(double xmin, double xmax, double ymin, double ymax, unsigned int Nbinsx, unsigned int Nbinsy);
26  virtual ~DHoughFind();
27  virtual const char* className(void) const {return static_className();}
28  static const char* static_className(void){return "DHoughFind";}
29 
30  inline double signof(double a){return a<0.0 ? -1.0:1.0;}
31 
32  void SetLimits(double xmin, double xmax, double ymin, double ymax, unsigned int Nbinsx, unsigned int Nbinsy);
33  void ResetHist(void);
34  DVector2 GetMaxBinLocation(void);
35  double GetMaxBinContent(void);
36  double GetSigmaX(void);
37  double GetSigmaY(void);
38  DVector2 Find(void);
39  DVector2 Find(const vector<DVector2> &points);
40 
41  void Fill(double x, double sigmax, double y, double sigmay);
42  static DVector2 GetMaxBinLocation(vector<const DHoughFind*> &houghs); // does not look at "this" object!
43  void Add(const DHoughFind* hough);
44 
45  void AddPoint(const DVector2 &point);
46  void AddPoint(const double &x, const double &y);
47  void AddPoints(const vector<DVector2> &points);
48  unsigned int GetNPoints(void){return points.size();}
49  vector<DVector2>GetPoints(void){return points;};
50  void ClearPoints(void);
51  void PrintHist(void);
52 
53  inline void FindIndexes(const DVector2 &pos, int &ix, int &iy);
54  inline double FindBeta(double xlo, double ylo, double widthx, double widthy, DVector2 &pos, DVector2 &step);
55  inline double FindBeta(const DVector2 &a, const DVector2 &b, const DVector2 &c, const DVector2 &d);
56 
57  TH2D* MakeIntoRootHist(string hname);
58 
59  protected:
60  vector<DVector2> points;
61 
62  double xmin, xmax, ymin, ymax;
63  unsigned int Nbinsx, Nbinsy;
64  double bin_widthx, bin_widthy, bin_size;
66  unsigned int imax_binx, imax_biny;
68 
69  double *hist;
70 
76 
77  private:
78 
79 };
80 
81 
82 // The following functions are inlined for speed
83 
84 //---------------------------------
85 // FindIndexes
86 //---------------------------------
87 inline void DHoughFind::FindIndexes(const DVector2 &pos, int &ix, int &iy)
88 {
89  ix = (int)floor((pos.X()-xmin)/bin_widthx);
90  iy = (int)floor((pos.Y()-ymin)/bin_widthy);
91 }
92 
93 //---------------------------------
94 // FindBeta
95 //---------------------------------
96 inline double DHoughFind::FindBeta(double xlo, double ylo, double widthx, double widthy, DVector2 &pos, DVector2 &step)
97 {
98  //DVector2 a0(xlo, ylo);
99  //DVector2 xdir(1.0, 0.0);
100  //DVector2 ydir(0.0, 1.0);
101  //DVector2 stepdir = step/step.Mod();
102  a0.Set(xlo, ylo);
103  xdir.Set(1.0, 0.0);
104  ydir.Set(0.0, 1.0);
105  stepdir = step/step.Mod();
106 
107  //vector<double> beta(4);
108  double beta[4];
109  //start = a0;
110  start.Set(xlo, ylo);
111  beta[0] = FindBeta(start, xdir, pos, stepdir);
112  //start = a0+widthx*xdir;
113  start.Set(xlo+widthx, ylo);
114  beta[1] = FindBeta(start, ydir, pos, stepdir);
115  //start = a0+widthx*xdir+widthy*ydir;
116  start.Set(xlo+widthx, ylo+widthy);
117  beta[2] = FindBeta(start, -1.0*xdir, pos, stepdir);
118  //start = a0+widthy*ydir;
119  start.Set(xlo, ylo+widthy);
120  beta[3] = FindBeta(start, -1.0*ydir, pos, stepdir);
121 
122  // Here, we want to choose the closest boundary of the bin which is not
123  // the boundary the point "pos" lies on. The values of beta give the
124  // distances to the lines which define the boundaries, but not clipped
125  // to the bin itself.
126 
127  // Loop over beta values
128  DVector2 bin_center(xlo+widthx/2.0, ylo+widthy/2.0);
129  double min_dist_to_center = 1.0E6;
130  double min_beta = 1.0E6;
131  for(unsigned int i=0; i<4; i++){
132  if(fabs(beta[i])<1.0E-6*bin_size)continue; // This is most likely due to our starting pos being on a boundary already. Ignore it.
133 
134  DVector2 delta_center = pos + beta[i]*stepdir - bin_center;
135  if(delta_center.Mod() < min_dist_to_center){
136  min_dist_to_center = delta_center.Mod();
137  min_beta = beta[i];
138  }
139  }
140 
141  return min_beta;
142 }
143 
144 //---------------------------------
145 // FindBeta
146 //---------------------------------
147 inline double DHoughFind::FindBeta(const DVector2 &a, const DVector2 &b, const DVector2 &c, const DVector2 &d)
148 {
149  /// Given the 2-D vectors a,b,c, and d find the scaler value "beta" such that
150  ///
151  /// a + alpha*b = c + beta*d
152  ///
153  /// It vectors b and d should be unit vectors.
154  /// The value of alpha is not returned, but can be calculated via:
155  ///
156  /// alpha = beta*(d.b) - b.(a-c)
157  ///
158  /// If the vectors b and d are parallel, then the return value will be inf
159  /// or whatever division by zero returns on the system.
160 
161  double ax = a.X();
162  double ay = a.Y();
163  double bx = b.X();
164  double by = b.Y();
165  double cx = c.X();
166  double cy = c.Y();
167  double dx = d.X();
168  double dy = d.Y();
169 
170  double k = bx*dx + by*dy;
171  double alphax = ax - cx;
172  double alphay = ay - cy;
173  double betax = dx - k*bx;
174  double betay = dy - k*by;
175 
176  return (betax*alphax + betay*alphay)/(1.0-k*k);
177 
178  //return (d*(a-c) - (d*b)*(b*(a-c)))/(1.0-pow(b*d, 2.0));
179 }
180 
181 #endif // _DHoughFind_
182 
double FindBeta(double xlo, double ylo, double widthx, double widthy, DVector2 &pos, DVector2 &step)
Definition: DHoughFind.h:96
unsigned int GetNPoints(void)
Definition: DHoughFind.h:48
unsigned int imax_biny
Definition: DHoughFind.h:66
DVector2 start
Definition: DHoughFind.h:75
double bin_widthy
Definition: DHoughFind.h:64
TVector2 DVector2
Definition: DVector2.h:9
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
DVector2 stepdir
Definition: DHoughFind.h:74
axes Fill(100, 100)
#define c
double * hist
Definition: DHoughFind.h:69
#define y
virtual const char * className(void) const
Definition: DHoughFind.h:27
double ymin
Definition: DHoughFind.h:62
double signof(double a)
Definition: DHoughFind.h:30
DVector2 xdir
Definition: DHoughFind.h:72
vector< DVector2 > GetPoints(void)
Definition: DHoughFind.h:49
DVector2 ydir
Definition: DHoughFind.h:73
vector< DVector2 > points
Definition: DHoughFind.h:60
bool max_bin_valid
Definition: DHoughFind.h:65
void FindIndexes(const DVector2 &pos, int &ix, int &iy)
Definition: DHoughFind.h:87
Double_t ymin
Definition: bcal_hist_eff.C:89
locHist_TaggerEnergy Add(locHist_TaggerEnergyAcc,-1.)
DVector2 a0
Definition: DHoughFind.h:71
double max_bin_content
Definition: DHoughFind.h:67
Double_t ymax
Definition: bcal_hist_eff.C:91
static const char * static_className(void)
Definition: DHoughFind.h:28
unsigned int Nbinsy
Definition: DHoughFind.h:63