MultiHoughT_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 /// \brief Run and test Hough transform algorithms
3 /// \authors messier@indiana.edu, mibaird@indiana.edu
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "TH1F.h"
8 #include "TH2S.h"
9 
13 
16 #include "Geometry/Geometry.h"
19 #include "RecoBase/Cluster.h"
20 #include "RecoBase/CellHit.h"
21 #include "RecoBase/HitList.h"
22 #include "RecoBase/HoughResult.h"
23 #include "RecoBase/FilterList.h"
24 #include "Utilities/func/MathUtil.h"
25 #include "Utilities/AssociationUtil.h"
26 
27 
28 
29 
30 
31 // MultiHoughT header
32 namespace hough {
33  ///
34  /// Run and test Hough transform algorithms
35  ///
36  class MultiHoughT : public art::EDProducer {
37  public:
38  explicit MultiHoughT(fhicl::ParameterSet const &p);
39  virtual ~MultiHoughT();
40 
41  void beginJob();
42  virtual void produce(art::Event &e);
43  void reconfigure(const fhicl::ParameterSet& p);
44  private:
45  void WriteHoughHistos(const MultiHough2P& h);
47 
48  private:
49  int fWriteHoughHistos; ///< Write the hough maps out?
50  int fWeightByQ; ///< Weight hits by charge?
51  double fRhoBinSz; ///< Size of rho bins (cm)
52  double fThetaBinSz; ///< Angular size of angle bins (degrees)
53  int fRhoSmoothWin; ///< Number of bins to smooth over
54  int fThetaSmoothWin; ///< Number of bins to smooth over
55  double fRsqr; ///< Distance cut on hit pairs
56  double fPco; ///< Number of sigma above average peak height to use as threshold cutoff in Hough space
57  int fLoops; ///< Max number of Multi-Hough loops
58  unsigned int fHitMin; ///< Minimum number of hits required per view
59  double fRemoveDist; ///< Distance cut for removing points on a Hough line
60  std::string fSlicerLabel; ///< Label for the process that made the slices
61  bool fObeyPreselection; ///< Check rb::IsFiltered?
62 
63  };
64 }
65 
66 
67 
68 namespace hough
69 {
70  //......................................................................
72  {
73  this->reconfigure(p);
74  this->produces< std::vector<rb::HoughResult> >();
75  this->produces<art::Assns<rb::HoughResult, rb::Cluster> >();
76  }
77 
78  //......................................................................
79 
81  { }
82 
83  //......................................................................
84 
86  {
87  fWriteHoughHistos = p.get<int> ("WriteHoughHistos");
88  fWeightByQ = p.get<int> ("WeightByQ");
89  fRhoBinSz = p.get<double> ("RhoBinSz");
90  fThetaBinSz = p.get<double> ("ThetaBinSz");
91  fRhoSmoothWin = p.get<int> ("RhoSmoothWin");
92  fThetaSmoothWin = p.get<int> ("ThetaSmoothWin");
93  fRsqr = p.get<double> ("Rsqr");
94  fPco = p.get<double> ("PeakCO");
95  fLoops = p.get<int> ("MHLoops");
96  fHitMin = p.get<unsigned int> ("HitMin");
97  fRemoveDist = p.get<double> ("RemoveDist");
98  fSlicerLabel = p.get<std::string> ("SlicerLabel");
99  fObeyPreselection = p.get<bool> ("ObeyPreselection");
100  }
101 
102 
103 
104  //......................................................................
106 
107 
108 
109  //.....................................................................
110 
112  {
113  unsigned int i;
114 
115  // Container for results to be stored in the event
116  std::unique_ptr<std::vector<rb::HoughResult> >
117  htr(new std::vector<rb::HoughResult>);
118  std::unique_ptr<art::Assns<rb::HoughResult, rb::Cluster> >
120 
121  // Pull out the time slices and shuffle them into a useful container
123  evt.getByLabel(fSlicerLabel,slice);
124 
125 
126  for (i=0; i<slice->size(); ++i) {
127  if(fObeyPreselection && rb::IsFiltered(evt,slice,i)) continue;
128  const rb::Cluster& s = (*slice)[i];
129 
130  //
131  // Skip any noise clusters
132  //
133  if (s.IsNoise()) continue;
134 
135 
136  //
137  // Build a stripped down hit list from the cells in the cluster but first,
138  // scrub the list of noisy hits.
139  //
140  rb::Cluster scrubbed(this->Scrub(s.AllCells()));
141  rb::HitList h(scrubbed, fWeightByQ);
142 
143 
144  //
145  // Construct the MultiHough transforms of the X and Y views
146  //
149 
150  // require there to be more than fHitMin hits in each view to do the multi-hough
151 
152  if(scrubbed.NXCell() < fHitMin || scrubbed.NYCell() < fHitMin) continue;
153  hx.MultiMap(h);
154  hy.MultiMap(h);
155 
156 
157  // require that there must be at least fHitMin hits in each view SEPERATELY
158  /*
159  if(scrubbed.NXCell() >= fHitMin) hx.MultiMap(h);
160  if(scrubbed.NYCell() >= fHitMin) hy.MultiMap(h);
161  */
162 
163  //
164  // Put the results into the event
165  //
166  art::Ptr<rb::Cluster> sptr(slice, i);
167  htr->push_back(hx.fH);
168  // Associate most recently added houghresult with the slice
169  util::CreateAssn(*this, evt, *htr, sptr, *assns);
170  htr->push_back(hy.fH);
171  // Associate most recently added houghresult with the slice
172  util::CreateAssn(*this, evt, *htr, sptr, *assns);
173 
174  //
175  // Optionally write the hough maps to the ROOT histogram file
176  //
177  if (fWriteHoughHistos) {
178  this->WriteHoughHistos(hx);
179  this->WriteHoughHistos(hy);
180  }
181 
182  }
183 
184  evt.put(std::move(htr));
185  evt.put(std::move(assns));
186 
187  }
188 
189  //......................................................................
190 
192  {
193  // write the Hough space map
195  TH2F* hist = h.MapToTH2F();
196  if(hist) f->make<TH2F>(*hist);
197  }
198 
199  //......................................................................
200  //function to remove hits with no neighbors from the slice
202  {
203 
204  std::vector<int> hit_numNeighbors;
205  //in an ideal detector a hit must have 1 neighbor
206  //inside this distance.
207  double maxGap = 60.0;
208  art::PtrVector<rb::CellHit> trimCells;
209  int iPlane, fPlane, iCell, fCell;
210  double goodCh, totCh;
211 
212  iPlane = 0;
213  iCell = 0;
214  fPlane = 0;
215  fCell = 0;
216 
219 
220  for (unsigned int i=0; i<c.size(); i++){
221  hit_numNeighbors.push_back(0);
222  }
223 
224  //loop to find a cells closest neighbor
225  for (unsigned int i=0; i<c.size(); i++){
226  double minDist = 9999.0;
227  int minCell = -1;
228  if (hit_numNeighbors[i] > 0) continue;
229  for(unsigned int j=0; j<c.size(); j++){
230  if (i == j) continue;
231  if (c[i]->View() != c[j]->View()) continue;
232  double xyz1[3], xyz2[3];
233  geom->Plane((c[i])->Plane())->Cell((c[i])->Cell())->GetCenter(xyz1);
234  geom->Plane((c[j])->Plane())->Cell((c[j])->Cell())->GetCenter(xyz2);
235  double dist = util::pythag(xyz1[0]-xyz2[0],
236  xyz1[1]-xyz2[1],
237  xyz1[2]-xyz2[2]);
238  if (dist < minDist){
239  minDist = dist;
240  minCell = j;
241  }
242  }
243  if (minCell == -1) continue;
244 
245  //define the rectangle of cells such that a cell and its
246  //nearest neighbor are at diagonal corners
247  if (c[i]->Plane() < c[minCell]->Plane()){
248  iPlane = c[i]->Plane();
249  fPlane = c[minCell]->Plane();
250  }
251  else {
252  iPlane = c[minCell]->Plane();
253  fPlane = c[i]->Plane();
254  }
255  if (c[i]->Cell() < c[minCell]->Cell()) {
256  iCell = c[i]->Cell();
257  fCell = c[minCell]->Cell();
258  }
259  else {
260  iCell = c[minCell]->Cell();
261  fCell = c[i]->Cell();
262  }
263 
264  //within the rectangle find the fraction of cells that are good
265  goodCh = 0;
266  totCh = 0;
267  for (int c=iCell; c<=fCell; c++){
268  for (int p=iPlane; p<=fPlane; p+=2){
269  totCh++;
270  if(!(badc->IsBad(p,c))) ++goodCh;
271  }
272  }
273  if ((minDist*goodCh/totCh) < maxGap){
274  hit_numNeighbors[i]++;
275  hit_numNeighbors[minCell]++;
276  }
277  }
278 
279  //keep the hits that have at least 1 neighbor in the slice
280  for (unsigned int i=0; i<c.size(); i++){
281  if (hit_numNeighbors[i] > 0){
282  trimCells.push_back(c[i]);
283  }
284  }
285 
286  return trimCells;
287  }
288 
289 
290 } // end namespace hough
291 ////////////////////////////////////////////////////////////////////////
292 
293 
294 
296 namespace hough
297 {
299 }
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.
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
Hough transform based on 2-point combinations.
Definition: MultiHough2P.h:21
art::PtrVector< rb::CellHit > Scrub(art::PtrVector< rb::CellHit > c)
Perform a "2 point" Hough transform on a collection of hits.
Definition: Hough2P.cxx:20
double fThetaBinSz
Angular size of angle bins (degrees)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double fPco
Number of sigma above average peak height to use as threshold cutoff in Hough space.
A collection of associated CellHits.
Definition: Cluster.h:47
double fRemoveDist
Distance cut for removing points on a Hough line.
const PlaneGeo * Plane(unsigned int i) const
virtual void produce(art::Event &e)
std::string fSlicerLabel
Label for the process that made the slices.
DEFINE_ART_MODULE(TestTMapFile)
TH2F * MapToTH2F() const
int fRhoSmoothWin
Number of bins to smooth over.
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 dist
Definition: runWimpSim.h:113
void MultiMap(rb::HitList h)
bool fObeyPreselection
Check rb::IsFiltered?
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
double fRsqr
Distance cut on hit pairs.
int fThetaSmoothWin
Number of bins to smooth over.
Just the essential information required for track fitting.
Definition: HitList.h:40
rb::HoughResult fH
Results of Hough transformation;.
Definition: MultiHough2P.h:91
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double fRhoBinSz
Size of rho bins (cm)
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
int fWriteHoughHistos
Write the hough maps out?
int fWeightByQ
Weight hits by charge?
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
int fLoops
Max number of Multi-Hough loops.
Definition: View.py:1
size_type size() const
Definition: PtrVector.h:308
void WriteHoughHistos(const MultiHough2P &h)
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
MultiHoughT(fhicl::ParameterSet const &p)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
Data resulting from a Hough transform on the cell hit positions.
unsigned int fHitMin
Minimum number of hits required per view.
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
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Float_t e
Definition: plot.C:35
bool IsBad(int plane, int cell)
void reconfigure(const fhicl::ParameterSet &p)
Encapsulate the geometry of one entire detector (near, far, ndos)