Prong.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Prong.cxx
3 // \brief A Prong is a Cluster with defined start position and direction
4 // \version $Id: Prong.cxx,v 1.12 2012-09-21 02:32:54 greenc Exp $
5 // \author $Author: greenc $
6 // \date $Date: 2012-09-21 02:32:54 $
7 ////////////////////////////////////////////////////////////////////////
8 #include "RecoBase/Prong.h"
9 
10 #include "Geometry/Geometry.h"
11 #include "GeometryObjects/Geo.h"
13 #include "RecoBase/CellHit.h"
14 
16 
17 namespace rb
18 {
19  //......................................................................
20  // 3D
22  TVector3 start, TVector3 dir, int id)
23  : rb::Cluster(cells, id)
24  {
25  SetStart(start);
26  SetDir(dir);
27  }
28 
29  Prong::Prong(const rb::Cluster& clust,
30  TVector3 start, TVector3 dir, int id)
31  : rb::Cluster(clust)
32  {
33  SetID(id);
34  SetStart(start);
35  SetDir(dir);
36  }
37 
38  //......................................................................
39  // 2D prong
41  geo::View_t view, double v0, double z0, double dv, double dz, int id)
42  : rb::Cluster(view, cells, id)
43  {
44  SetStart(v0, z0);
45  SetDir(dv, dz);
46  }
47 
48  //......................................................................
49  // 2D prong
50  Prong::Prong(const std::vector< art::Ptr<rb::CellHit> >& cells,
51  geo::View_t view, double v0, double z0, double dv, double dz, int id)
52  : rb::Cluster(view, cells, id)
53  {
54  SetStart(v0, z0);
55  SetDir(dv, dz);
56  }
57 
58  Prong::Prong(const rb::Cluster& clust,
59  double v0, double z0, double dv, double dz, int id)
60  : rb::Cluster(clust)
61  {
62  SetID(id);
63  SetStart(v0, z0);
64  SetDir(dv, dz);
65  }
66 
67  //......................................................................
68 
70  {
71  }
72 
73  //......................................................................
74 
75  void Prong::SetStart(TVector3 start)
76  {
77  fPrecalcTotalGeV = -1; // Invalidate any cached value
78 
79  // I've removed this assert to allow for 2D prongs to set a 3D start
80  // point if they so choose.
81  // assert(fView == geo::kXorY);
82 
83  fStart = start;
84  }
85 
86  //......................................................................
87 
88  void Prong::SetStart(double v0, double z0)
89  {
90  fPrecalcTotalGeV = -1; // Invalidate any cached value
91 
92  assert(fView == geo::kX || fView == geo::kY);
93 
94  if(fView == geo::kX){
95  fStart = TVector3(v0, 0, z0);
96  }
97  else{
98  fStart = TVector3(0, v0, z0);
99  }
100  }
101 
102  //......................................................................
103 
104  void Prong::SetDir(TVector3 dir)
105  {
106  fPrecalcTotalGeV = -1; // Invalidate any cached value
107 
108  // I've removed this assert to allow for 2D prongs to set a 3D start
109  // direction if they so choose.
110  // assert(fView == geo::kXorY);
111 
112  fDir = dir.Unit(); // Ensure normalized
113  }
114 
115  //......................................................................
116 
117  void Prong::SetDir(double dv, double dz)
118  {
119  fPrecalcTotalGeV = -1; // Invalidate any cached value
120 
121  assert(fView == geo::kX || fView == geo::kY);
122 
123  if(fView == geo::kX){
124  fDir = TVector3(dv, 0, dz);
125  }
126  else{
127  fDir = TVector3(0, dv, dz);
128  }
129 
130  fDir = fDir.Unit(); // Ensure normalized
131  }
132 
133  //......................................................................
134 
135  double Prong::W(const rb::CellHit* chit) const
136  {
137 
138  //TODO: If the prong is 2D, return a position of 0
139  //for the opposite view. This is in line with what happens
140  //in the w method for a track and cluster. In the future
141  //trying to do this should fail in all cases. This patch
142  //keeps anything from breaking until there is a long-term fix.
143 
144  // Removing this line so that 2D prongs can get a sensible W value.
145  // if (chit->View() == fView) return 0;
146 
148 
149  double xyz[3];
150  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
151  const double z1 = xyz[2];
152 
153  /// \todo CJB - we're just moving the point in the XY plane to get
154  /// it on the line. Maybe we should take the route perpendicular
155  /// to the line instead. Similarly for Track
156 
157  const double dt = (z1-fStart.Z())/fDir.Z();
158 
159  if(chit->View() == geo::kX) return fStart.Y()+dt*fDir.Y();
160  if(geom->DetId() == novadaq::cnv::kTESTBEAM) return -(fStart.X()+dt*fDir.X()); // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
161  return fStart.X()+dt*fDir.X();
162 
163  /*
164  // Solve the problem of the closest approach of the trajectory to the line
165  // made by the cell by projecting down to a 2D trajectory and 2D point for
166  // the cell.
167  TVector3 start2D = fStart;
168  TVector3 dir2D = fDir;
169  TVector3 pt2D(xyz);
170  const int otherView = 1-chit->View();
171  start2D[otherView] = dir2D[otherView] = pt2D[otherView] = 0;
172 
173  // Actually find the closest point
174  TVector3 closest;
175  geo::ClosestApproach(pt2D, start2D, dir2D, closest);
176 
177  // Need to solve for the W position of the line, the exact part that we
178  // projected away above. So just figure out Z and evaluate the prong there.
179  const double dt = (closest.Z()-fStart.Z())/fDir.Z();
180 
181  return (fStart+dt*fDir)[otherView];
182  */
183  }
184 
185  //......................................................................
186  double Prong::TotalLength() const
187  {
189 
190  double ret = 0;
191 
192  const unsigned int cellMax = NCell();
193 
194  for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
195  const art::Ptr<rb::CellHit> chit = Cell(cellIdx);
196 
197  TVector3 xyz;
198  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
199 
200  const TVector3 dxyz = (chit->View() == geo::kX) ? TVector3(0, 1, 0)
201  : TVector3(1, 0, 0);
202 
203  double dist;
205  xyz, xyz+dxyz,
206  &dist);
207 
208  ret = std::max(ret, dist);
209  }
210  return ret;
211  }
212 
213  //......................................................................
214  double Prong::DistanceFromStart(double z) const
215  {
216  // A prong is a straight line so this is easy
217 
218  // For a completely vertical track any answer is nonsense
219  if(fDir.Z() == 0) return 0;
220 
221  return (z-fStart.Z())/fDir.Z();
222  }
223 
224  //......................................................................
225 
226  bool Prong::operator< (const Prong& other) const
227  {
228  return this->Start().Z() < other.Start().Z();
229  }
230 
231 } // end namespace rb
232 ////////////////////////////////////////////////////////////////////////
T max(const caf::Proxy< T > &a, T b)
virtual ~Prong()
Definition: Prong.cxx:69
double fPrecalcTotalGeV
-1 = uninitialized
Definition: Cluster.h:324
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Prong.cxx:135
TVector3 fDir
Direction at starting point.
Definition: Prong.h:112
Prong()
Definition: Prong.h:22
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
virtual double DistanceFromStart(double z) const
Definition: Prong.cxx:214
A collection of associated CellHits.
Definition: Cluster.h:47
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double dist
Definition: runWimpSim.h:113
unsigned short Cell() const
Definition: CellHit.h:40
Definition: Cand.cxx:23
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
double dz[NP][NC]
unsigned int NCell() const
Number of cells in either view.
Definition: Cluster.h:110
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
geo::View_t fView
view this cluster is in
Definition: Cluster.h:316
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Collect Geo headers and supply basic geometry functions.
void SetID(int id)
Definition: Cluster.h:74
bool operator<(const Prong &other) const
Definition: Prong.cxx:226
Perform a "2 point" Hough transform on a collection of hits.
z
Definition: test.py:28
TVector3 fStart
Start location (xyz, cm)
Definition: Prong.h:111
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
A Cluster with defined start position and direction.
Definition: Prong.h:19
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
TDirectory * dir
Definition: macro.C:5
const int cellMax
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
Encapsulate the geometry of one entire detector (near, far, ndos)