VertexDT_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 ///
3 /// \file VertexDT_module.cxx
4 /// \brief Reconstruct rb::Vertex basd on input rb::HoughResult
5 /// \author gsdavies@indiana.edu, xbbu@fnal.gov, messier@indiana.edu
6 ///
7 //////////////////////////////////////////////////////////////////////
12 
13 #include <algorithm>
14 #include <utility>
15 
16 #include "TNtuple.h"
17 #include "TH1F.h"
18 #include "TH2F.h"
19 
25 
27 #include "Geometry/Geometry.h"
29 #include "GeometryObjects/Geo.h"
31 #include "RecoBase/CellHit.h"
32 #include "RecoBase/Cluster.h"
33 #include "RecoBase/FilterList.h"
34 #include "RecoBase/HoughResult.h"
35 #include "RecoBase/Vertex.h"
36 #include "Utilities/AssociationUtil.h"
37 #include "Utilities/func/MathUtil.h"
38 
39 
40 /// VertexDT finding algorithm
41 namespace vdt
42 {
43  ///
44  /// A stand-alone module for running the VertexDT algorithm
45  ///
46  class VertexDT : public art::EDProducer
47  {
48  public:
49  explicit VertexDT(fhicl::ParameterSet const &pset);
50  virtual ~VertexDT();
51 
52  void beginJob();
53  void reconfigure(const fhicl::ParameterSet& pset);
54  void produce(art::Event &evt);
55  void endJob();
56 
58 
59  private:
60 
61  void FindMinMax(std::vector<double> Vhits, std::vector<double> Zhits,
62  double* minV, double* maxV,
63  double* minZ, double* maxZ);
64 
65  bool PassPDistReq(double x, double y, double vzm, double vzb);
66 
67  bool fMakeHitNt; ///< Make hit-level debugging ntuple
68  // Tunable parameters
69  unsigned int fMinHit; ///< Slices must have at least this many hits
70  unsigned int fMinHitX; ///< Miniumm hits in x view / slice
71  unsigned int fMinHitY; ///< Minimum hit in y view / slice
72  unsigned int fMaxHoughResults; ///< Only use the top hough results
73  double fHoughPeakThr; ///< Threshold to accept a peak as good
74  double fMaxPdistHits; ///< hits' perpendicular dist to line must be less than this value
75  double fMaxPdistEndP; ///< end points perpendicular dist to line must be less than this value
76 
77  // Folder labels
78  std::string fCellHitLabel; ///< Module label for cell hits
79  std::string fSliceLabel; ///< Module label for Clusters
80  std::string fHoughResultLabel; ///< Module label for Hough lines
81 
82  TNtuple* fHitNt;
83  };
84 }
85 
86 namespace vdt
87 {
88  //......................................................................
90  fHitNt(0)
91  {
92  this->reconfigure(pset);
93  this->produces< std::vector<rb::Vertex> >();
94  this->produces<art::Assns<rb::Vertex, rb::Cluster> >();
95  }
96 
97  //......................................................................
99  {
100  if(fMakeHitNt){
102  fHitNt = tfs->make<TNtuple>("hnt","hnt","evt:run:srun:xy:z:view:chi2:chi2endpoint:chi2line");
103  }
104  }
105  //......................................................................
107  {
108  fCellHitLabel = pset.get< std::string >("CellHitLabel");
109  fSliceLabel = pset.get< std::string >("SliceLabel");
110  fHoughResultLabel = pset.get< std::string >("HoughResultLabel");
111 
112 #define SET(T,N) { f##N = pset.get<T>(#N); }
113 
114  SET(bool, MakeHitNt);
115  SET(unsigned int, MinHit);
116  SET(unsigned int, MinHitX);
117  SET(unsigned int, MinHitY);
118 
119  SET(double, HoughPeakThr);
120  SET(unsigned int, MaxHoughResults);
121 
122  SET(double, MaxPdistHits);
123  SET(double, MaxPdistEndP);
124 
125 #undef SET
126  }
127 
128  //......................................................................
130  {
131  }
132 
133 
134  //......................................................................
136  {
137  std::unique_ptr< std::vector<rb::Vertex> > vtx(new std::vector<rb::Vertex>);
138  std::unique_ptr< art::Assns<rb::Vertex, rb::Cluster> >
140 
142 
143  //
144  // Pull out the time slices and hough results for the event
145  //
147  evt.getByLabel(fSliceLabel,slice);
148 
149  // Object to get hough lines associated with a slice
151 
152  for (unsigned int i = 0; i < slice->size(); ++i) {
153  const rb::Cluster& s = (*slice)[i];
154 
155  if(s.IsNoise()) continue;
156 
157  // Scrub this cluster of isolated hits,
158  // when slicer is fixed this will no longer be necessary ???
159  rb::Cluster neighbors(this->Scrub(s.AllCells()));
160 
161  //
162  // Skip the noise cluster and require some minimum number of hits
163  //
164  if(neighbors.NXCell() < fMinHitX) continue;
165  if(neighbors.NYCell() < fMinHitY) continue;
166  if(neighbors.NCell() < fMinHit) continue;
167 
168  //
169  // Find any hough results that go with this time slice
170  //
171  const rb::HoughResult* hrx = 0;
172  const rb::HoughResult* hry = 0;
173 
174  std::vector<const rb::HoughResult*> hr = fmhr.at(i);
175 
176  for (unsigned j = 0; j < hr.size(); ++j) {
177  if(hr[j]->fView==geo::kX) {
178  hrx = hr[j];
179  }
180  if(hr[j]->fView==geo::kY) {
181  hry = hr[j];
182  }
183  }
184  //
185  // No hough results, no joy.
186  //
187  if(hrx == 0 || hry == 0) continue;
188  if(hrx->fPeak.size() < 1) continue;
189  if(hry->fPeak.size() < 1) continue;
190 
191  //
192  // Select hough lines for each view
193  //
194  unsigned int nx = 0;
195  for(unsigned j = 0; j < hrx->fPeak.size(); ++j){
196  if(hrx->fPeak[j].fH/hrx->fPeak[j].fThr >= fHoughPeakThr) ++nx;
197  } // loop XZ view hough lines
198 
199  unsigned int ny = 0;
200  for(unsigned j = 0; j < hry->fPeak.size(); ++j) {
201  if(hry->fPeak[j].fH/hry->fPeak[j].fThr >= fHoughPeakThr) ++ny;
202  } // loop YZ view hough lines
203 
204  unsigned int nsegs = 1;
205  if(nx > nsegs) nsegs = nx;
206  if(ny > nsegs) nsegs = ny;
207  if(nsegs > fMaxHoughResults) nsegs = fMaxHoughResults;
208 
209 
210  // Build segments in XZ, YZ views
211  std::vector<vdt::Segment> segXZ, segYZ;
212 
213  // Loop over all cells and sort into xhits/yhits
214  //Function MakeSegment
215  double seedX = 0.;
216  double seedY = 0.;
217  double seedZ = 0.;
218 
219  for(unsigned int idx = 0; idx < nsegs; ++idx){
220 
221  for(int view: {0,1}){
222 
223  std::vector<double> Vhits;
224  std::vector<double> Zhits;
225 
226  unsigned int ncell = (view == 0) ? neighbors.NXCell() : neighbors.NYCell();
227  for(unsigned int nIdx = 0; nIdx < ncell; ++nIdx){
228  const rb::CellHit& nhit = (view == 0) ? *neighbors.XCell(nIdx) : *neighbors.YCell(nIdx);
229 
230  TVector3 xyz;
231  const geo::CellGeo* gcell = geom->Plane(nhit.Plane())->Cell(nhit.Cell());
232  gcell->GetCenter(xyz);
233  double V = (view == 0) ? xyz[0] : xyz[1];
234 
235  const rb::HoughResult* hr = (view == 0) ? hrx : hry;
236  if(idx < hr->fPeak.size()){
237  double tmpm, tmpc;
238  hr->SlopeIcept(idx, &tmpm, &tmpc);
239  if(PassPDistReq(xyz[2], V, tmpm, tmpc)){
240  (view == 0) ? Vhits.push_back(xyz[0]) : Vhits.push_back(xyz[1]);
241  Zhits.push_back(xyz[2]);
242  }
243  }
244  } // Loop over cells in X/Y view
245 
246  double minV, minZ;
247  double maxV, maxZ;
248  double minZ_t = geom->DetLength();
249  this->FindMinMax(Vhits, Zhits, &minV, &maxV, &minZ, &maxZ);
250 
251  vdt::Segment seg ;
252  // Now we're able to flip the segments in the case
253  // of it being too vertical...
254  if(fabs(maxV - minV) > fabs(maxZ - minZ)) seg.Flip();
255 
256  for(unsigned int hIdx = 0; hIdx < Zhits.size(); ++hIdx){
257  seg.AddPoint(Zhits[hIdx], 3., Vhits[hIdx], 2.);
258  }
259  seg.SetEndPoint1(minZ, minV, 2., 40.); ///TODO: fcl params
260  seg.SetEndPoint2(maxZ, maxV, 2., 40.);
261 
262 
263  (view == 0) ? segXZ.push_back(seg) : segYZ.push_back(seg);
264 
265  if(fMakeHitNt){
266  float x[9];
267  // can't fill ntuple until after segment is complete
268  /// TODO: Make this a debug option
269  for(unsigned int hIdx = 0; hIdx < Zhits.size(); ++hIdx){
270  x[0] = evt.id().event();
271  x[1] = evt.run();
272  x[2] = evt.subRun();
273  x[3] = Vhits[hIdx];
274  x[4] = Zhits[hIdx];
275  x[5] = view;
276  x[6] = seg.Chi2(Zhits[hIdx], Vhits[hIdx]);
277  x[7] = seg.Chi2EndPoint(Zhits[hIdx], Vhits[hIdx]);
278  x[8] = seg.Chi2Line(Zhits[hIdx], Vhits[hIdx]);
279  fHitNt->Fill(x);
280  }
281  }
282 
283  if(idx == 0){
284  if(view == 0) seedX = minV;
285  if(view == 1) seedY = minV;
286  if(minZ < minZ_t) minZ_t = minZ;
287  seedZ = minZ_t;
288  }
289  } // end loop over view
290  } // MaxHoughResults loop
291 
292 
293  // TODO: Minimize 2D vertices and use their vertex as
294  // seed for 3D vertex minimisation
295  // For now let's just minimize 3D vertex
296  vdt::Minimizer3D mn3D(1, segXZ.size(), segYZ.size());
297 
298  for(unsigned int ixz = 0; ixz < segXZ.size(); ++ixz)
299  mn3D.SetXZSegment(ixz, segXZ[ixz]);
300  for(unsigned int iyz = 0; iyz < segYZ.size(); ++iyz)
301  mn3D.SetYZSegment(iyz, segYZ[iyz]);
302 
303  mn3D.SetVertex(0, seedX, seedY, seedZ);
304  mn3D.SetLambda(30.0); /// TODO: Set this in fcl
305  mn3D.Anneal(1.0E-5, 1.0, 2.0); /// TODO: fcl parameters
306 
307  double Vx, Vy, Vz;
308  double ex, ey, ez;
309  mn3D.GetVertex(0, &Vx, &Vy, &Vz, &ex, &ey, &ez);
310 
311  art::Ptr<rb::Cluster> sptr(slice, i);
312  rb::Vertex v;
313 
314  v.SetX(Vx);
315  v.SetY(Vy);
316  v.SetZ(Vz);
317  v.SetT(sptr->MeanTNS());
318 
319  mf::LogDebug("VertexDT") << "Event Vertex: " << evt.id().event() << "\n"
320  << "Vertex X: " << v.GetX() << "\n"
321  << "Vertex Y: " << v.GetY() << "\n"
322  << "Vertex Z: " << v.GetZ();
323 
324  vtx->push_back(v);
325  util::CreateAssn(*this, evt, *vtx, sptr, *vAssns);
326 
327  } // end loop on slices
328 
329  evt.put(std::move(vtx));
330  evt.put(std::move(vAssns));
331  } // end produce
332 
333  //......................................................................
334  // Function to remove hits with no neighbors from the slice
336  {
337 
338  std::vector<int> hit_numNeighbors;
339 
340  // In an ideal detector a hit must have 1 neighbor
341  // inside this distance.
342  double maxGap = 60.0;
343 
344  art::PtrVector<rb::CellHit> trimCells;
345 
346  int iPlane, fPlane, iCell, fCell;
347  double goodCh, totCh;
348 
349  iPlane = 0;
350  iCell = 0;
351  fPlane = 0;
352  fCell = 0;
353 
356 
357  for(unsigned int i = 0; i < c.size(); i++){
358  hit_numNeighbors.push_back(0);
359  }
360 
361  // Loop to find a cells closest neighbor
362  for(unsigned int i = 0; i < c.size(); i++){
363  double minDist = 9999.0;
364  int minCell = -1;
365 
366  if(hit_numNeighbors[i] > 0) continue;
367 
368  for(unsigned int j = 0; j < c.size(); j++){
369 
370  if(i == j) continue;
371  if(c[i]->View() != c[j]->View()) continue;
372  double xyz1[3], xyz2[3];
373  geom->Plane((c[i])->Plane())->Cell((c[i])->Cell())->GetCenter(xyz1);
374  geom->Plane((c[j])->Plane())->Cell((c[j])->Cell())->GetCenter(xyz2);
375  double dist = util::pythag(xyz1[0]-xyz2[0],
376  xyz1[1]-xyz2[1],
377  xyz1[2]-xyz2[2]);
378  if(dist < minDist){
379  minDist = dist;
380  minCell = j;
381  }
382  }
383 
384  if(minCell == -1) continue;
385 
386  // Define the rectangle of cells such that a cell and its
387  // nearest neighbor are at diagonal corners
388  if(c[i]->Plane() < c[minCell]->Plane()){
389  iPlane = c[i]->Plane();
390  fPlane = c[minCell]->Plane();
391  }
392  else{
393  iPlane = c[minCell]->Plane();
394  fPlane = c[i]->Plane();
395  }
396  if(c[i]->Cell() < c[minCell]->Cell()) {
397  iCell = c[i]->Cell();
398  fCell = c[minCell]->Cell();
399  }
400  else{
401  iCell = c[minCell]->Cell();
402  fCell = c[i]->Cell();
403  }
404 
405  // Within the rectangle find the fraction of cells that are good
406  goodCh = 0;
407  totCh = 0;
408  for(int c = iCell; c <= fCell; c++){
409  for(int p = iPlane; p <= fPlane; p+=2){
410  totCh++;
411  if(!(badc->IsBad(p,c))) ++goodCh;
412  }
413  }
414  if((minDist*goodCh/totCh) < maxGap){
415  hit_numNeighbors[i]++;
416  hit_numNeighbors[minCell]++;
417  }
418  }
419 
420  // Keep the hits that have at least 1 neighbor in the slice
421  for(unsigned int i = 0; i < c.size(); i++){
422  if(hit_numNeighbors[i] > 0){
423  trimCells.push_back(c[i]);
424  }
425  }
426 
427  return trimCells;
428  }
429 
430  //------------------------------------------------------------------
431  void VertexDT::FindMinMax(std::vector<double> Vhits, std::vector<double> Zhits,
432  double* minV, double* maxV,
433  double* minZ, double* maxZ)
434  {
435 
437 
438  // Function MinMax
439  double minZ_t = geom->DetLength();
440  double maxZ_t = 0.;
441  double minV_t = geom->DetHalfWidth();
442  double maxV_t = -geom->DetHalfHeight();
443  // Now calculate the minV, minZ values for this segment
444  // We're tying the max V(x/y) to the min/max Z-coordinate
445  // Optimise this in future
446  for(unsigned int hIdx = 0; hIdx < Zhits.size(); ++hIdx){
447  if(Zhits[hIdx] > maxZ_t){
448  maxZ_t = Zhits[hIdx];
449  maxV_t = Vhits[hIdx];
450  }
451  if(Zhits[hIdx] < minZ_t){
452  minZ_t = Zhits[hIdx];
453  minV_t = Vhits[hIdx];
454  }
455  }
456  *minV = minV_t;
457  *minZ = minZ_t;
458  *maxV = maxV_t;
459  *maxZ = maxZ_t;
460  }
461 
462  //......................................................................
463  bool VertexDT::PassPDistReq(double x, double y, double vzm, double vzb)
464  {
465  double pdist = 0.;
466 
467  pdist = sqrt(geo::DsqrToLine(x, y, 0.0,
468  vzb, 100.0, 100.0*vzm + vzb));
469 
470  if(pdist < fMaxPdistHits) return true;
471 
472  return false;
473  }
474 
475  //......................................................................
477  {
478  }
479 
480 } // end namespace vdt
481 
483 ////////////////////////////////////////////////////////////////////////
#define SET(T, N)
virtual ~VertexDT()
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
double fMaxPdistHits
hits&#39; perpendicular dist to line must be less than this value
SubRunNumber_t subRun() const
Definition: Event.h:72
void SetT(double t)
Definition: Vertex.h:32
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
unsigned int fMinHit
Slices must have at least this many hits.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void SetEndPoint1(double x, double y, double resin, double resex)
Definition: Segment.cxx:48
double GetX() const
Definition: Vertex.h:23
unsigned short Plane() const
Definition: CellHit.h:39
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Definition: HoughResult.h:61
art::PtrVector< rb::CellHit > Scrub(art::PtrVector< rb::CellHit > c)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
void produce(art::Event &evt)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
3D vertex "minimizer" function declarations
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
unsigned int fMaxHoughResults
Only use the top hough results.
nhit
Definition: demo1.py:25
double GetY() const
Definition: Vertex.h:24
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double GetZ() const
Definition: Vertex.h:25
double dist
Definition: runWimpSim.h:113
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
void Flip(bool flip=true)
Definition: Segment.cxx:33
void SlopeIcept(unsigned int i, double *m, double *b) const
Slope and intercepts of Hough lines.
Definition: HoughResult.cxx:62
bool PassPDistReq(double x, double y, double vzm, double vzb)
void FindMinMax(std::vector< double > Vhits, std::vector< double > Zhits, double *minV, double *maxV, double *minZ, double *maxZ)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
void AddPoint(double x, double sigx, double y, double sigy)
Definition: Segment.cxx:84
double Chi2(double x, double y)
Definition: Segment.cxx:115
Float_t E
Definition: plot.C:20
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
T get(std::string const &key) const
Definition: ParameterSet.h:231
double Chi2EndPoint(double x, double y)
Definition: Segment.cxx:155
Collect Geo headers and supply basic geometry functions.
const double j
Definition: BetheBloch.cxx:29
std::string fCellHitLabel
Module label for cell hits.
double DetHalfHeight() const
Vertex location in position and time.
Definition: View.py:1
2D vertex "minimizer" function declarations
VertexDT(fhicl::ParameterSet const &pset)
size_type size() const
Definition: PtrVector.h:308
std::string fHoughResultLabel
Module label for Hough lines.
unsigned int fMinHitX
Miniumm hits in x view / slice.
2D vertex "energy" function
void reconfigure(const fhicl::ParameterSet &pset)
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
EventNumber_t event() const
Definition: EventID.h:116
double fHoughPeakThr
Threshold to accept a peak as good.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
double DetHalfWidth() const
Data resulting from a Hough transform on the cell hit positions.
void SetX(double x)
Definition: Vertex.h:29
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void SetZ(double z)
Definition: Vertex.h:31
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
VertexDT finding algorithm.
unsigned int fMinHitY
Minimum hit in y view / slice.
void SetXZSegment(unsigned int k, Segment &s)
Definition: Minimizer3D.cxx:38
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
double fMaxPdistEndP
end points perpendicular dist to line must be less than this value
RunNumber_t run() const
Definition: Event.h:77
Summary from a Hough transform applied to a group of cell hits.
Definition: HoughResult.h:35
Encapsulate the cell geometry.
Definition: CellGeo.h:25
bool IsBad(int plane, int cell)
bool fMakeHitNt
Make hit-level debugging ntuple.
double Chi2Line(double x, double y)
Definition: Segment.cxx:136
void SetY(double y)
Definition: Vertex.h:30
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
TNtuple * fHitNt
void SetEndPoint2(double x, double y, double resin, double resex)
Definition: Segment.cxx:66
std::string fSliceLabel
Module label for Clusters.