HoughT_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 "TH2S.h"
9 
13 
14 #include "Geometry/Geometry.h"
15 #include "HoughTransform/Hough2P.h"
17 #include "RecoBase/Cluster.h"
18 #include "RecoBase/HitList.h"
19 #include "RecoBase/HoughResult.h"
20 #include "Utilities/AssociationUtil.h"
21 
22 
23 
24 /// Hough transform header
25 namespace hough {
26  ///
27  /// Run and test Hough transform algorithms
28  ///
29  class HoughT : public art::EDProducer {
30  public:
31  explicit HoughT(fhicl::ParameterSet const &p);
32  virtual ~HoughT();
33 
34  void beginJob();
35  virtual void produce(art::Event &e);
36  void reconfigure(const fhicl::ParameterSet& p);
37  private:
38  void WriteHoughHistos(const Hough2P& h);
39  bool is_good_slice(rb::Cluster const& slice, geo::View_t view) const;
40  private:
41  int fWriteHoughHistos; ///< Write the hough maps out?
42  int fWeightByQ; ///< Weight hits by charge?
43  double fRhoBinSz; ///< Size of rho bins (cm)
44  double fThetaBinSz; ///< Angular size of angle bins (degrees)
45  int fRhoSmoothWin; ///< Number of bins to smooth over
46  int fThetaSmoothWin; ///< Number of bins to smooth over
47  double fRsqr; ///< Distance cut on hit pairs
48  double fPco; ///< Number of sigma above average peak height to use as threshold cutoff in Hough space
49  unsigned int fHitMin; ///< Minimum number of hits required per view
50  std::string fSlicerLabel; ///< Label for the process that made slices
51  };
52 }
53 
54 
55 
56 namespace hough
57 {
58  //.....................................................................
60  {
61  reconfigure(p);
62  produces<std::vector<rb::HoughResult> >();
63  produces<art::Assns<rb::HoughResult, rb::Cluster> >();
64  }
65 
66  //......................................................................
67 
69  {
70  }
71 
72  //......................................................................
73 
75  {
76  fWriteHoughHistos = p.get<int> ("WriteHoughHistos");
77  fWeightByQ = p.get<int> ("WeightByQ");
78  fRhoBinSz = p.get<double> ("RhoBinSz");
79  fThetaBinSz = p.get<double> ("ThetaBinSz");
80  fRhoSmoothWin = p.get<int> ("RhoSmoothWin");
81  fThetaSmoothWin = p.get<int> ("ThetaSmoothWin");
82  fRsqr = p.get<double> ("Rsqr");
83  fPco = p.get<double> ("PeakCO");
84  fHitMin = p.get<unsigned int> ("HitMin");
85  fSlicerLabel = p.get<std::string> ("SlicerLabel");
86  }
87 
88  //......................................................................
89 
91 
92  //......................................................................
93 
95  {
96  unsigned int i;
98 
99  // Container for results to be stored in the event
100  std::unique_ptr<std::vector<rb::HoughResult> >
101  htr(new std::vector<rb::HoughResult>);
102  std::unique_ptr<art::Assns<rb::HoughResult, rb::Cluster> >
104 
105  // Pull out the time slices and shuffle them into a useful container
107  evt.getByLabel(fSlicerLabel,slice);
108 
109  for (i=0; i<slice->size(); ++i)
110  {
111  const rb::Cluster& s = (*slice)[i];
112  //
113  // Skip any noise clusters
114  //
115  if (s.IsNoise()) continue;
116 
117  //
118  // Build a stripped down hit list from the cells in the cluster
119  //
121 
122  //
123  // Construct the Hough transforms of the X and Y views
124  //
127  if(is_good_slice(s, geo::View_t::kX)) hx.Map(h);
128  if(is_good_slice(s, geo::View_t::kY)) hy.Map(h);
129 
130  /*
131  std::cout << "*=== Slice " << i << std::endl;
132  std::cout << "Found " << hx.fH.fPeak.size() << " peaks in X" << std::endl;
133  std::cout << "Found " << hy.fH.fPeak.size() << " peaks in Y" << std::endl;
134  */
135 
136  //
137  // Put the results into the event and associate the last htr element
138  // with the slice.
139  //
140  art::Ptr<rb::Cluster> slice_ptr(slice, i);
141  htr->push_back(hx.fH);
142  util::CreateAssn(*this, evt, *htr, slice_ptr, *assns);
143  htr->push_back(hy.fH);
144  util::CreateAssn(*this, evt, *htr, slice_ptr, *assns);
145 
146  //
147  // Optionally write the hough maps to the ROOT histogram file
148  //
149  if (fWriteHoughHistos)
150  {
151  if (s.NXCell() > 0)
152  this->WriteHoughHistos(hx);
153 
154  if (s.NYCell() > 0)
155  this->WriteHoughHistos(hy);
156  }
157  }
158 
159  evt.put(std::move(htr));
160  evt.put(std::move(assns));
161  }
162 
163  //......................................................................
164 
166  {
168  f->make<TH2F>(*h.fHspaceW);
169  f->make<TH2S>(*h.fHspaceI);
170  }
171 
172 
173 
175  {
176  if (slice.NCell(view) <= fHitMin)
177  return false;
178 
179  // It can happen that a slice contains several hits from one cell,
180  // but occurring at differnt times.
181  if (slice.MinPlane(view) == slice.MaxPlane(view) &&
182  slice.MinCell(view) == slice.MaxCell(view))
183  return false;
184 
185  return true;
186  }
187 
188 } // end namespace hough
189 ////////////////////////////////////////////////////////////////////////
190 
191 
192 
194 namespace hough
195 {
197 }
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 NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
int fWeightByQ
Weight hits by charge?
Perform a "2 point" Hough transform on a collection of hits.
Definition: Hough2P.cxx:20
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void WriteHoughHistos(const Hough2P &h)
IslandsTH2 * fHspaceI
Islands for x view.
Definition: Hough2P.h:77
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
virtual void produce(art::Event &e)
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
DEFINE_ART_MODULE(TestTMapFile)
int fRhoSmoothWin
Number of bins to smooth over.
void reconfigure(const fhicl::ParameterSet &p)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
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;})
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Hough transform based on 2-point combinations.
Definition: Hough2P.h:18
const XML_Char * s
Definition: expat.h:262
std::string fSlicerLabel
Label for the process that made slices.
Just the essential information required for track fitting.
Definition: HitList.h:40
unsigned int fHitMin
Minimum number of hits required per view.
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual ~HoughT()
double fRsqr
Distance cut on hit pairs.
bool is_good_slice(rb::Cluster const &slice, geo::View_t view) const
HoughT(fhicl::ParameterSet const &p)
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
::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.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
rb::HoughResult fH
Results of Hough transformation;.
Definition: Hough2P.h:78
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
void geom(int which=0)
Definition: geom.C:163
Definition: geo.h:1
double fThetaBinSz
Angular size of angle bins (degrees)
void Map(const rb::HitList &h)
Definition: Hough2P.cxx:119
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
int fThetaSmoothWin
Number of bins to smooth over.
double fPco
Number of sigma above average peak height to use as threshold cutoff in Hough space.
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Float_t e
Definition: plot.C:35
int fWriteHoughHistos
Write the hough maps out?
Encapsulate the geometry of one entire detector (near, far, ndos)
double fRhoBinSz
Size of rho bins (cm)