Track3D.cxx
Go to the documentation of this file.
1 #include "Monopole/Track3D.h"
2 
3 #include "Monopole/Constants.h"
4 #include "Geometry/Geometry.h"
5 
7 // #include <canvas/Utilities/Exception.h>
8 #include <art/Utilities/Exception.h>
9 
10 #include <TVector3.h>
11 
12 #include <algorithm>
13 #include <iostream>
14 #include <map>
15 
16 
17 
19 {
20 }
21 
22 
23 
24 mono::Track3D::Track3D(Track2D const& x_track, Track2D const& y_track) :
25  default_value_ (-9999),
38  x_track_ (x_track),
39  y_track_ (y_track)
40 {
41  x0_ = x_track.intercept;
42  y0_ = y_track.intercept;
43  z0_ = 0.0;
44 
45  dxdz_ = x_track.slope;
46  dydz_ = y_track.slope;
47 
48  bool xt_fit_valid =
50  bool yt_fit_valid =
52 
53  if (xt_fit_valid && yt_fit_valid)
54  {
57  beta_ = velocity / Constants::SPEED_OF_LIGHT;
58  }
59 
60 
61  /*
62  Hit Matching
63 
64  Note that this can be done here or at the beginning of the constructor,
65  but it does not affect the x_track and y_track themselves, which are used
66  for the fit.
67  */
68  for (auto const& x_hit : x_track.cluster.AllCells())
69  {
70  bool has_match = false;
71  for (auto const& y_hit : y_track.cluster.AllCells())
72  if (abs(int(x_hit->Plane()) - int(y_hit->Plane())) <= 1)
73  has_match = true;
74 
75  if (has_match)
76  cluster_.Add(x_hit);
77  }
78 
79  for (auto const& y_hit : y_track.cluster.AllCells())
80  {
81  bool has_match = false;
82  for (auto const& x_hit : x_track.cluster.AllCells())
83  if (abs(int(x_hit->Plane()) - int(y_hit->Plane())) <= 1)
84  has_match = true;
85 
86  if (has_match)
87  cluster_.Add(y_hit);
88  }
89 
92 }
93 
94 
95 
96 TVector3 mono::Track3D::start() const
97 {
98  return TVector3(x0_, y0_, z0_);
99 }
100 
101 
102 
103 TVector3 mono::Track3D::velocity() const
104 {
105  return TVector3(dxdt_, dydt_, dzdt_);
106 }
107 
108 
109 
110 double mono::Track3D::dxdz() const
111 {
112  return dxdz_;
113 }
114 
115 
116 
117 double mono::Track3D::dydz() const
118 {
119  return dydz_;
120 }
121 
122 
123 
125 {
126  if (view == geo::View_t::kX)
127  return dzdt_x_track_;
128 
129  if (view == geo::View_t::kY)
130  return dzdt_y_track_;
131 
132  return default_value_;
133 }
134 
135 
136 
137 double mono::Track3D::r2_xt() const
138 {
139  return r2_xt_;
140 }
141 
142 
143 
144 double mono::Track3D::r2_yt() const
145 {
146  return r2_yt_;
147 }
148 
149 
150 
152 {
153  if (view == geo::View_t::kX)
154  return r2_zt_x_track_;
155 
156  if (view == geo::View_t::kY)
157  return r2_zt_y_track_;
158 
159  return default_value_;
160 }
161 
162 
163 
165 {
166  return max_time_gap_xz_;
167 }
168 
169 
170 
172 {
173  return max_time_gap_yz_;
174 }
175 
176 
177 
178 double mono::Track3D::beta() const
179 {
180  return beta_;
181 }
182 
183 
184 
186 {
187  return cluster_;
188 }
189 
190 
191 
193 (Track2D const& track,
194  double & dvdt, double & dzdt,
195  double & r2_vt, double & r2_zt) const
196 {
197  dvdt = default_value_;
198  dzdt = default_value_;
199  double N = track.cluster.NCell();
200  // we need at least two hits to get a non-zero variance
201  if (N <= 2)
202  return false;
203 
205 
206  std::map<std::string, double> mean = { {"v", 0.0}, {"z", 0.0}, {"t", 0.0} };
207  for (auto const& hit : track.cluster.AllCells())
208  {
209  double xyz[3];
210  geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
211  TVector3 hit_vector(xyz);
212 
213  mean.at("v") += geometry->CellTpos(hit->Plane(), hit->Cell()) / N;
214  mean.at("z") += hit_vector.z() / N;
215  mean.at("t") += hit->TNS() / N;
216  }
217 
218  std::map<std::string, double> variance =
219  { {"t", 0.0}, {"v", 0.0}, {"z", 0.0}, {"vt", 0.0}, {"zt", 0.0} };
220  for (auto const& hit : track.cluster.AllCells())
221  {
222  double xyz[3];
223  geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
224  TVector3 hit_vector(xyz);
225  double v = geometry->CellTpos(hit->Plane(), hit->Cell());
226  double z = hit_vector.z();
227  double t = hit->TNS();
228 
229  variance.at("t") += get_variance(t, mean.at("t")) / (N - 1);
230  variance.at("v") += get_variance(v, mean.at("v")) / (N - 1);
231  variance.at("z") += get_variance(z, mean.at("z")) / (N - 1);
232  variance.at("vt") +=
233  get_covariance(v, mean.at("v"), t, mean.at("t")) / (N - 1);
234  variance.at("zt") +=
235  get_covariance(z, mean.at("z"), t, mean.at("t")) / (N - 1);
236  }
237 
238  if (variance.at("t") == 0 || variance.at("v") == 0 || variance.at("z") == 0)
239  return false;
240  // throw art::Exception(art::errors::InvalidNumber)
241  // << "Monopole/Track3D.cxx\n"
242  // << "Variance is zero. Aborting.\n" << std::endl;
243 
244 
245  //
246  // Set all of the results:
247  //
248  dvdt = variance.at("vt") / variance.at("t");
249  dzdt = variance.at("zt") / variance.at("t");
250 
251  // The intercept is currently not used.
252  // double v0 = mean.at("v") - dvdt * mean.at("t");
253  // double z0 = mean.at("z") - dzdt * mean.at("t");
254 
255  // Calculate the squared correlation coefficienct (r^2):
256  r2_vt =
257  std::pow(variance.at("vt"), 2) / (variance.at("v") * variance.at("t"));
258  r2_zt =
259  std::pow(variance.at("zt"), 2) / (variance.at("z") * variance.at("t"));
260 
261  return true;
262 }
263 
264 
265 
267 {
268  double result = default_value_;
269 
270  std::set<double> times;
271  for (auto const& hit : hits)
272  times.insert(hit->TNS());
273 
274  std::adjacent_find
275  (times.begin(), times.end(), [&](double object, double next_object)
276  {
277  double gap = next_object - object;
278  result = std::max(gap, result);
279 
280  // This return is meaningless, but a boolean return is required by
281  // adjacent find. We choose false, so the loop goes until the end.
282  return false;
283  });
284 
285  return result;
286 }
287 
288 
289 
291 (double const& a, double const& mean_a) const
292 {
293  return (a - mean_a) * (a - mean_a);
294 }
295 
296 
297 
299 (double const& a, double const& mean_a,
300  double const& b, double const& mean_b) const
301 {
302  return (a - mean_a) * (b - mean_b);
303 }
304 
305 
306 
307 std::ostream& mono::operator<<(std::ostream& os, mono::Track3D const& t)
308 {
310 
311  os << "Track3D: dxdz, dydz = " << t.dxdz() << ", " << t.dydz() << "\n";
312 
313  for (auto const& hit : t.cluster().AllCells())
314  {
316  if (hit->View() == geo::View_t::kX)
317  view = "XZ";
318  else if (hit->View() == geo::View_t::kY)
319  view = "YZ";
320  else
321  view = "BAD";
322 
323  double xyz[3];
324  geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
325  TVector3 hit_vector(xyz);
326  double v = geometry->CellTpos(hit->Plane(), hit->Cell());
327  double z = hit_vector.z();
328  double t = hit->TNS();
329 
330  os << "Track3D: View " << view
331  << " Hit: t, v, z = " << t << ", " << v << ", " << z << "\n";
332  }
333 
334  return os;
335 }
double time_gap(art::PtrVector< rb::CellHit > const &hits) const
Definition: Track3D.cxx:266
T max(const caf::Proxy< T > &a, T b)
double max_time_gap_yz() const
Definition: Track3D.cxx:171
double CellTpos(unsigned int ip, unsigned int ic, double w=0.0) const
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
double max_time_gap_yz_
Definition: Track3D.h:67
const art::PtrVector< rb::CellHit > & XCells() const
Get all cells from the x-view.
Definition: Cluster.h:124
double get_variance(double const &a, double const &mean_a) const
Definition: Track3D.cxx:291
double z0_
Definition: Track3D.h:61
double r2_xt_
Definition: Track3D.h:65
double beta() const
Definition: Track3D.cxx:178
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
constexpr T pow(T x)
Definition: pow.h:75
double dydz() const
Definition: Track3D.cxx:117
double x0_
Definition: Track3D.h:61
Definition: event.h:19
double y0_
Definition: Track3D.h:61
A collection of associated CellHits.
Definition: Cluster.h:47
void abs(TH1 *hist)
const double SPEED_OF_LIGHT
Definition: Constants.h:8
double dzdt_y_track_
Definition: Track3D.h:64
bool time_fit(Track2D const &track, double &dvdt, double &dzdt, double &r2_vt, double &r2_zt) const
Definition: Track3D.cxx:193
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double intercept
Definition: Track2D.h:27
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double dxdz() const
Definition: Track3D.cxx:110
double dydz_
Definition: Track3D.h:62
TVector3 start() const
Definition: Track3D.cxx:96
double dzdt_x_track_
Definition: Track3D.h:64
Definition: novas.h:112
std::ostream & operator<<(std::ostream &os, Track3D const &t)
Definition: Track3D.cxx:307
void hits()
Definition: readHits.C:15
const double a
const art::PtrVector< rb::CellHit > & YCells() const
Get all cells from the x-view.
Definition: Cluster.h:126
rb::Cluster cluster() const
Definition: Track3D.cxx:185
double dxdz_
Definition: Track3D.h:62
double max_time_gap_xz() const
Definition: Track3D.cxx:164
TVector3 velocity() const
Definition: Track3D.cxx:103
double r2_zt_y_track_
Definition: Track3D.h:66
z
Definition: test.py:28
double dzdt(geo::View_t view) const
Definition: Track3D.cxx:124
double max_time_gap_xz_
Definition: Track3D.h:67
double r2_zt_x_track_
Definition: Track3D.h:66
Track2D y_track_
Definition: Track3D.h:70
double dxdt_
Definition: Track3D.h:63
rb::Cluster cluster_
Definition: Track3D.h:71
double dzdt_
Definition: Track3D.h:63
double r2_yt() const
Definition: Track3D.cxx:144
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double slope
Definition: Track2D.h:27
double r2_yt_
Definition: Track3D.h:65
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
double r2_zt(geo::View_t view) const
Definition: Track3D.cxx:151
rb::Cluster cluster
Definition: Track2D.h:26
Definition: geo.h:1
const hit & b
Definition: hits.cxx:21
Track2D x_track_
Definition: Track3D.h:70
double default_value_
Definition: Track3D.h:59
double r2_xt() const
Definition: Track3D.cxx:137
double get_covariance(double const &a, double const &mean_a, double const &b, double const &mean_b) const
Definition: Track3D.cxx:299
double variance(Eigen::VectorXd x)
double dydt_
Definition: Track3D.h:63
double beta_
Definition: Track3D.h:68
Encapsulate the geometry of one entire detector (near, far, ndos)