Track3D.cxx
Go to the documentation of this file.
1 #include "HoughTrack/Track3D.h"
2 
3 #include "Utilities/func/MathUtil.h"
4 
5 #include "HoughTrack/Constants.h"
6 #include "Geometry/Geometry.h"
7 #include "RecoBaseHit/CellHit.h"
8 
11 
12 #include <TVector3.h>
13 
14 #include <iostream>
15 #include <map>
16 
17 
19 {
20 }
21 
22 
23 
24 htk::Track3D::Track3D(Track2D const& x_track, Track2D const& y_track) :
25  default_value_(-9999),
26  x_track_(x_track), y_track_(y_track)
27 {
28  x0_ = x_track.intercept;
29  y0_ = y_track.intercept;
30  z0_ = 0.0;
31 
32  dxdz_ = x_track.slope;
33  dydz_ = y_track.slope;
34 
35  bool xt_fit_valid = time_fit(x_track, dxdt_, dzdt_x_track_);
36  bool yt_fit_valid = time_fit(y_track, dydt_, dzdt_y_track_);
37 
38  if (xt_fit_valid && yt_fit_valid)
39  {
40  // add a check to make sure it's almost the same for the two views
41  // having a large discrepancy in the two views might be a powerful
42  // discriminant
43 
46  beta_ = velocity / Constants::SPEED_OF_LIGHT;
47  } else {
50  }
51 
52  for (auto const& x_hit : x_track.cluster.AllCells())
53  {
54  cluster_.Add(x_hit);
55  }
56 
57  for (auto const& y_hit : y_track.cluster.AllCells())
58  {
59  cluster_.Add(y_hit);
60  }
61 }
62 
63 
64 
65 TVector3 htk::Track3D::start() const
66 {
67  return TVector3(x0_, y0_, z0_);
68 }
69 
70 
71 
72 TVector3 htk::Track3D::velocity() const
73 {
74  return TVector3(dxdt_, dydt_, dzdt_);
75 }
76 
77 
78 
79 double htk::Track3D::dxdz() const
80 {
81  return dxdz_;
82 }
83 
84 
85 
86 double htk::Track3D::dydz() const
87 {
88  return dydz_;
89 }
90 
91 
92 
94 {
95  if (view == geo::View_t::kX)
96  return dzdt_x_track_;
97 
98  if (view == geo::View_t::kY)
99  return dzdt_y_track_;
100 
101  return -9999;
102 }
103 
104 
105 
106 double htk::Track3D::beta() const
107 {
108  return beta_;
109 }
110 
111 
112 
114 {
115  return cluster_;
116 }
117 
118 
119 
121 (Track2D const& track, double & dvdt, double & dzdt) const
122 {
123  dvdt = default_value_;
124  dzdt = default_value_;
125  double N = track.cluster.NCell();
126  // we need at least two hits to get a non-zero variance
127  if (N <= 2)
128  return false;
129 
131 
132  std::map<std::string, double> mean = { {"v", 0.0}, {"z", 0.0}, {"t", 0.0} };
133  for (auto const& hit : track.cluster.AllCells())
134  {
135  double xyz[3];
136  geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
137  TVector3 hit_vector(xyz);
138 
139  mean.at("v") += geometry->CellTpos(hit->Plane(), hit->Cell()) / N;
140  mean.at("z") += hit_vector.z() / N;
141  mean.at("t") += hit->TNS() / N;
142  }
143 
144  std::map<std::string, double> variance =
145  { {"t", 0.0}, {"vt", 0.0}, {"zt", 0.0} };
146  for (auto const& hit : track.cluster.AllCells())
147  {
148  double xyz[3];
149  geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
150  TVector3 hit_vector(xyz);
151  double v = geometry->CellTpos(hit->Plane(), hit->Cell());
152  double z = hit_vector.z();
153  double t = hit->TNS();
154 
155  variance.at("t") += get_variance(t, mean.at("t")) / (N - 1);
156  variance.at("vt") +=
157  get_covariance(v, mean.at("v"), t, mean.at("t")) / (N - 1);
158  variance.at("zt") +=
159  get_covariance(z, mean.at("z"), t, mean.at("t")) / (N - 1);
160  }
161 
162  if (variance.at("t") == 0)
163  {
164  // throw art::Exception(art::errors::InvalidNumber)
165  // << "Monopole/Track3D.cxx\n"
166  // << "Variance of t is zero. Aborting.\n" << std::endl;
167 
168  dvdt = 999999;
169  dzdt = 999999;
170  }
171 
172  else
173  {
174  dvdt = variance.at("vt") / variance.at("t");
175  dzdt = variance.at("zt") / variance.at("t");
176  }
177 
178  //double v0 = mean.at("v") - dvdt * mean.at("t");
179  //double z0 = mean.at("z") - dzdt * mean.at("t");
180 
181  return true;
182 }
183 
184 
185 
187 (double const& a, double const& mean_a) const
188 {
189  return (a - mean_a) * (a - mean_a);
190 }
191 
192 
193 
195 (double const& a, double const& mean_a,
196  double const& b, double const& mean_b) const
197 {
198  return (a - mean_a) * (b - mean_b);
199 }
double z0_
Definition: Track3D.h:45
double intercept
Definition: Track2D.h:26
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
rb::Cluster cluster_
Definition: Track3D.h:52
double dxdz_
Definition: Track3D.h:46
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
double slope
Definition: Track2D.h:26
const double SPEED_OF_LIGHT
Definition: Constants.h:8
double dydt_
Definition: Track3D.h:47
rb::Cluster cluster() const
Definition: Track3D.cxx:113
rb::Cluster cluster
Definition: Track2D.h:25
double dxdz() const
Definition: Track3D.cxx:79
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
double dzdt_
Definition: Track3D.h:47
double get_variance(double const &a, double const &mean_a) const
Definition: Track3D.cxx:187
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
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double default_value_
Definition: Track3D.h:43
double dzdt_x_track_
Definition: Track3D.h:48
Track2D x_track_
Definition: Track3D.h:51
double x0_
Definition: Track3D.h:45
TVector3 velocity() const
Definition: Track3D.cxx:72
const double a
z
Definition: test.py:28
double beta() const
Definition: Track3D.cxx:106
double dydz_
Definition: Track3D.h:46
double y0_
Definition: Track3D.h:45
Track2D y_track_
Definition: Track3D.h:51
double dzdt(geo::View_t view) const
Definition: Track3D.cxx:93
TVector3 start() const
Definition: Track3D.cxx:65
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
Definition: geo.h:1
const hit & b
Definition: hits.cxx:21
double beta_
Definition: Track3D.h:49
double get_covariance(double const &a, double const &mean_a, double const &b, double const &mean_b) const
Definition: Track3D.cxx:195
double variance(Eigen::VectorXd x)
double dxdt_
Definition: Track3D.h:47
Encapsulate the geometry of one entire detector (near, far, ndos)
double dydz() const
Definition: Track3D.cxx:86
double dzdt_y_track_
Definition: Track3D.h:48
bool time_fit(Track2D const &track, double &dvdt, double &dzdt) const
Definition: Track3D.cxx:121