ToFPositionFinder.cxx
Go to the documentation of this file.
1 #include "ToFPositionFinder.h"
2 //************************CONSTRUCTORS***********************
4 
5 }
6 
7 
8 //***********************3-CORNER-METHOD*********************
9 // Returns a posn_t (x,y) representing the position in a coordinate system where
10 // point B is the origin, point A is at (0,sidelength) and C is (sidelength,0)
11 // In other words, B is adjacent to both A and C. A is on the y-axis and C
12 // is on the x-axis.
13 // tA is the time at which the counter at corner A received a hit.
14 // Similarly for tB and tC
16 
17  //calculate the deltas
18  double dtA = tA - tB;
19  double dtC = tC - tB;
20 
21  //the solution to a quadratic...
22  double m = 4 * (dtA*dtA + dtC*dtC - s*s);
23  double n = 4 * (dtA*dtA*dtA + dtC*dtC*dtC - dtA*s*s - dtC*s*s);
24  double o = dtA*dtA*dtA*dtA + dtC*dtC*dtC*dtC + 2 * (s*s*s*s - s*s*dtA*dtA - s*s*dtC*dtC);
25  //...is the distance in ns from the middle corner
26  double B = ( -n - std::sqrt(n*n - 4*m*o) ) / ( 2*m );
27 
28  //once we know the absolute distance to one corner,
29  //we effectively know the absolute distance to any corner,
30  //so calculating x and y becomes trivial.
31  double x = ( B*B + s*s - (B+dtC)*(B+dtC) ) / (2*s);
32  double y = ( B*B + s*s - (B+dtA)*(B+dtA) ) / (2*s);
33 
34  posn_t loc = {x,y};
35  return loc;
36 }
37 
38 
39 //***********************4-CORNER-METHOD**********************
40 
41 // Get positions from each corner and transform to one coordinate system
42 std::vector<beamlinereco::posn_t> beamlinereco::ToFPositionFinder::Get4Posns(double tA, double tB, double tC, double tD, double s) {
43  //reconstruct from each corner
44  posn_t locB = Get3CornerPosition(tA, tB, tC, s);
45  posn_t locC = Get3CornerPosition(tB, tC, tD, s);
46  posn_t locD = Get3CornerPosition(tC, tD, tA, s);
47  posn_t locA = Get3CornerPosition(tD, tA, tB, s);
48 
49  //transform to coordinate system with B at origin
50  locB = RotateCoords(locB, s, 0);
51  locC = RotateCoords(locC, s, 1);
52  locD = RotateCoords(locD, s, 2);
53  locA = RotateCoords(locA, s, 3);
54 
55  return {locA,locB,locC,locD};
56 }
57 
58 
59 beamlinereco::posn_t beamlinereco::ToFPositionFinder::Get4CornerPosition(double tA, double tB, double tC, double tD, double s) {
60  std::vector<posn_t> posns = Get4Posns(tA, tB, tC, tD, s);
61  posn_t avg = PosnAvg(posns);
62  return avg;
63 }
64 
65 double beamlinereco::ToFPositionFinder::Get4CornerSpread(double tA, double tB, double tC, double tD, double s) {
66  std::vector<posn_t> posns = Get4Posns(tA,tB,tC,tD,s);
67  double spread = GetSpread(posns);
68  return spread;
69 }
70 
71 //*************************TEST-FUNCTIONS***********************
73  std::cout << "*********************RUNNING TESTS*********************\n";
74  bool pass3 = Run3CornerTests();
75  bool pass4 = Run4CornerTests();
76  std::cout << "*******************DONE RUNNING TESTS******************\n";
77  return pass3 && pass4;
78 }
79 
80 
82  std::cout << "***********RUNNING 3_CORNER TESTS**********\n";
83  bool allPassed = true;
84  std::vector<testset_t> testSets;
85  testSets = { { 100, 100, 100, 0, 1, {.5,.5} },
86  { 101, 100+std::sqrt(2), 101, 0, 1, {1,1} },
87  { 1+std::sqrt(10), 1+std::sqrt(18), 1+std::sqrt(10), 0, 4, {3,3} },
88  { 5+std::sqrt(26), 5+std::sqrt(26), 5+std::sqrt(106), 0, 10, {1,5} },
89  { 7+std::sqrt(20), 7+std::sqrt(53), 7+std::sqrt(130), 0, 11, {2,7} } };
90 
91  for(auto ts : testSets) {
92  bool thisPassed = Run3CornerTest(ts);
93  allPassed = allPassed && thisPassed;
94  }
95  std::cout << (allPassed ? "All 3-corner tests passed.\n" : "1+ 3-corner tests failed.\n");
96  std::cout << "***********DONE RUNNING 3_CORNER TESTS**********\n";
97 
98  return allPassed;
99 }
100 
101 
103  std::cout << "Parameters A:" << ts.tA << " B:" << ts.tB << " C:" << ts.tC << " s:" << ts.s << "\n";
104  posn_t calculated = Get3CornerPosition(ts.tA, ts.tB, ts.tC, ts.s);
105  std::cout << "Calculated result: (" << calculated.x << "," << calculated.y;
106  std::cout << ") Expected result: (" << ts.expected.x << "," << ts.expected.y << ")\n\n";
107  return std::abs(ts.expected.x - calculated.x) < EPSILON
108  && std::abs(ts.expected.y - calculated.y) < EPSILON;
109 }
110 
111 
113  std::cout << "***********RUNNING 4_CORNER TESTS**********\n";
114  bool allPassed = true;
115  std::vector<testset_t> testSets;
116  testSets = { { 100, 100, 100, 100, 1, {.5,.5} },
117  { 101, 100+std::sqrt(2), 101, 100, 1, {1,1} },
118  { 1+std::sqrt(10), 1+std::sqrt(18), 1+std::sqrt(10), 1+std::sqrt(2), 4, {3,3} },
119  { 5+std::sqrt(26), 5+std::sqrt(26), 5+std::sqrt(106), 5+std::sqrt(106), 10, {1,5} },
120  { 7+std::sqrt(20), 7+std::sqrt(53), 7+std::sqrt(130), 7+std::sqrt(97), 11, {2,7} } };
121 
122  for(auto ts : testSets) {
123  allPassed = Run4CornerTest(ts) && allPassed;
124  }
125  std::cout << (allPassed ? "All 4-corner tests passed.\n" : "1+ 4-corner tests failed.\n");
126  std::cout << "***********RUNNING 4_CORNER TESTS**********\n";
127  return allPassed;
128 }
129 
130 
132  std::cout << "Parameters A:" << ts.tA << " B:" << ts.tB << " C:" << ts.tC << " D:" << ts.tD << " s:" << ts.s << "\n";
133  posn_t calculated = Get4CornerPosition(ts.tA, ts.tB, ts.tC, ts.tD, ts.s);
134  std::cout << "Calculated result: (" << calculated.x << "," << calculated.y;
135  std::cout << ") Expected result: (" << ts.expected.x << "," << ts.expected.y << ")\n\n";
136  return std::abs(ts.expected.x - calculated.x) < EPSILON
137  && std::abs(ts.expected.y - calculated.y) < EPSILON;
138 }
139 
140 
141 //******************************POSN UTILITIES***************************
142 
143 //take the average of some posns (center of mass)
145  double sumX = 0; double sumY = 0;
146  for(auto p : posns) { sumX += p.x; sumY += p.y; }
147  double avgX = sumX/posns.size();
148  double avgY = sumY/posns.size();
149  return { avgX, avgY};
150 }
151 
152 //rotate the coordinate system of a posn right around the square.
153 // if x,y are given in a coordinate system where
154 // A=(s,s); B=(0,s); C=(0,0); D=(s,0)
155 // this function will transform the coords such that
156 // A=(0,s); B=(0,0); C=(s,0); D=(s,s)
158  return { s - p.y, p.x };
159 }
160 
161 //rotate the coordinate system n times.
163  for(int i = 0; i < n; i++) {
164  p = Transform90posn(p,s);
165  }
166  return p;
167 }
168 
169 //get the distance between two posns.
171  return std::sqrt( std::pow(a.x-b.x, 2) + std::pow(a.y-b.y, 2) );
172 }
173 
174 //get the spread on a list of posns (average distance to center of mass)
175 double beamlinereco::ToFPositionFinder::GetSpread(std::vector<posn_t> posns) {
176  posn_t center = PosnAvg(posns);
177  double avgDist = 0;
178  for(posn_t p : posns) {
179  avgDist += PosnDist(p,center);
180  }
181  avgDist /= posns.size();
182  return avgDist;
183 }
184 
185 //Test whether a given posn is physically possible, given the side length
187  return 0 < p.x && p.x < s && 0 < p.y && p.y < s;
188 }
189 
190 //Find the distance to the nearest of the listed posns.
191 double beamlinereco::ToFPositionFinder::GetClosestDist(posn_t p, std::vector<posn_t> others) {
192  double mindist = 0;
193  for(posn_t other : others) {
194  double thisdist = PosnDist(p, other);
195  mindist = thisdist < mindist ? thisdist : mindist;
196  }
197  return mindist;
198 }
199 
200 //Find the distance to the nearest corner of the scintillator
202  return GetClosestDist( p, { {0,0}, {0,s}, {s,0}, {s,s} } );
203 }
204 
205 //********************NICE HISTOGRAM OUTPUT*********************
206 TH2D beamlinereco::ToFPositionFinder::GetSummaryHist(TString name, TString title, double tA, double tB, double tC, double tD, double s) {
207  std::vector<posn_t> posns = Get4Posns(tA,tB,tC,tD,s);
208  posn_t avg = PosnAvg(posns);
209  TH2D hist = TH2D(name, title, HIST_N_BINS, -s, 2*s, HIST_N_BINS, -s, 2*s);
210  for(posn_t p : posns) {
211  hist.Fill(p.x,p.y);
212  hist.Fill(avg.x,avg.y);
213  }
214  hist.SetOption("colz");
215  return hist;
216 }
const XML_Char * name
Definition: expat.h:151
TH2D GetSummaryHist(TString name, TString title, double tA, double tB, double tC, double tD, double s)
double GetClosestDist(posn_t p, std::vector< posn_t > others)
double GetDistanceToNearestCorner(posn_t p, double s)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
const double EPSILON
posn_t Get3CornerPosition(double tA, double tB, double tC, double s)
float abs(float number)
Definition: d0nt_math.hpp:39
posn_t RotateCoords(posn_t p, double s, int n)
const XML_Char * s
Definition: expat.h:262
double Get4CornerSpread(double tA, double tB, double tC, double tD, double s)
double GetSpread(std::vector< posn_t > posns)
const double a
bool IsPhysical(posn_t p, double s)
int HIST_N_BINS
OStream cout
Definition: OStream.cxx:6
posn_t Get4CornerPosition(double tA, double tB, double tC, double tD, double s)
posn_t Transform90posn(posn_t p, double s)
posn_t PosnAvg(std::vector< posn_t > posns)
std::vector< posn_t > Get4Posns(double tA, double tB, double tC, double tD, double s)
const hit & b
Definition: hits.cxx:21
double PosnDist(posn_t a, posn_t b)