Slicer4D_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file Slicer4D_module.cc
3 // \brief This is based on a DBSCAN density based clustering algorithm
4 // found in "A Density-Based Algorithm for Discovering Clusters in
5 // Large Spatial Databases with Noise" by M.Ester et.al. (1996)
6 // \version $Id: Slicer4D_module.cc,v 1.0 2013-05-01 17:00:00 mbaird42 Exp $
7 // \authors mbaird42@fnal.gov, messier@indiana.edu
8 ///////////////////////////////////////////////////////////////////////////
9 
15 
16 #include "Calibrator/Calibrator.h"
17 #include "Geometry/Geometry.h"
18 #include "NovaDAQConventions/DAQConventions.h"
19 #include "RecoBase/CellHit.h"
20 #include "RecoBase/Cluster.h"
21 #include "Slicer/PointManager.h"
22 
23 #include <cmath>
24 #include <string>
25 
26 
27 
28 // Slicer4D header
29 namespace slicer {
30 
31  class Slicer4D : public art::EDProducer {
32  public:
33  explicit Slicer4D(fhicl::ParameterSet const &pset);
34  virtual ~Slicer4D();
35 
36  void beginRun(art::Run& run);
37  void produce(art::Event& evt);
38  void endJob();
39 
40  protected:
41  // internal methods
42  bool GenerateCluster(unsigned int point, int ClID);
43 
44  // configuration settings
45  std::string fCellHitInput; ///< Read CellHits from this input
46  bool fDoNotSliceAnything; ///< If true, one big slice is written out
47  bool fApplyNeighborhoodCleanUp; ///< scrub loner hits from the slice?
48 
49  double fEpsilon; ///< Minimum causal separation value for neighbors
50  unsigned int fMinPts; ///< minimum number of casual neighbors for a hit to be in a cluster
51  bool fUseLowPEPen; ///< Add a causal distance penalty for PE (low PE is INCREASES the distance)
52  double fPEPen; ///< PE value for which the low PE penalty will be 1.
53  double fDistPen; ///< Distance (in cm) for which the distance penalty will be 1 (roughly one strong interaction length.)
54  double fOpVPlPen; ///< Allowed plane gap for hits in opposite views for which the distance penalty will be 1.
55  unsigned int fMinHitsPerView; ///< Minimum number of hits per view to make a cluster
56 
58 
60 
61  // global variables/objects
62  PointManager fPtMan; ///< Helper class for the 4D clustering
63  };
64 
65 }
66 
67 
68 
69 
70 
71 namespace slicer
72 {
73  //----------------------------------------------------------------------
75  fCellHitInput (pset.get< std::string >("CellHitInput") ),
76  fDoNotSliceAnything(pset.get< bool >("DoNotSliceAnything")),
77  fApplyNeighborhoodCleanUp(pset.get< bool >("ApplyNeighborhoodCleanUp")),
78  fEpsilon (0),
79  fMinPts (0),
80  fUseLowPEPen (0),
81  fPEPen (0),
82  fDistPen (0),
83  fOpVPlPen (0),
84  fMinHitsPerView (0),
85  fMaxEventHits (pset.get< int >("MaxEventHits")),
86  fPSetNDOS (pset.get<fhicl::ParameterSet>("ndos")),
87  fPSetND (pset.get<fhicl::ParameterSet>("nd")),
88  fPSetFD (pset.get<fhicl::ParameterSet>("fd")),
89  fPSetTB (pset.get<fhicl::ParameterSet>("tb"))
90  {
91  produces< std::vector<rb::Cluster> >();
92  }
93 
94  //----------------------------------------------------------------------
96  { }
97 
98  //----------------------------------------------------------------------
100  {
102 
103  fhicl::ParameterSet pset;
104  switch(geom->DetId()){
105  case novadaq::cnv::kNDOS:
106  pset = fPSetNDOS;
107  break;
109  pset = fPSetND;
110  break;
112  pset = fPSetFD;
113  break;
115  pset = fPSetTB;
116  break;
117  default:
118  assert(0 && "Unknown detector");
119  }
120 
121  fEpsilon = pset.get<double >("Epsilon");
122  fMinPts = pset.get<unsigned int>("MinPts");
123  fUseLowPEPen = pset.get<bool >("UseLowPEPen");
124  fPEPen = pset.get<double >("PEPen");
125  fDistPen = pset.get<double >("DistPen");
126  fOpVPlPen = pset.get<double >("OppViewPlanePen");
127  fMinHitsPerView = pset.get<unsigned int>("MinHitsPerView");
128  }
129 
130  //----------------------------------------------------------------------
132  {
133  // Print statements used when tuning the slicer parameters.
134  /*
135  std::cout << "\n\n========== End Job ==========\n"
136  << "End Job: (eps, minpts, PEPen, DistPen, OpVPlPen, MinHitsPerView) "
137  << fEpsilon << " "
138  << fMinPts << " "
139  << fPEPen << " "
140  << fDistPen << " "
141  << fOpVPlPen << " "
142  << fMinHitsPerView << "\n";
143  */
144  }
145 
146  //----------------------------------------------------------------------
148  {
149  // NOTE: There is currently no requirement that a slice have hits in
150  // both views.
151 
153 
154  // Get the cell hits
156  evt.getByLabel(fCellHitInput, hitcol);
157 
158  // Load the cell hits into a vector for easy use
159  std::vector<art::Ptr<rb::CellHit > > hitlist;
160  hitlist.reserve(hitcol->size());
161  for(unsigned int i = 0; i < hitcol->size(); ++i){
162  art::Ptr<rb::CellHit> hit(hitcol, i);
163  hitlist.push_back(hit);
164  }
165 
166  // Safety in case of a very large number of hits in the event, which can
167  // take a long time to slice, and the results probably won't be useful for
168  // downstream modules anyway.
169  if(fMaxEventHits > 0 && int(hitlist.size()) > fMaxEventHits){
170  mf::LogWarning("Slicer4D") << "Very large number of hits ("
171  << hitlist.size()
172  << ") in event. Higher than configured limit ("
173  << fMaxEventHits
174  << "). Will put all hits in the noise slice.";
175  rb::Cluster noise;
176  for(art::Ptr<rb::CellHit> c: hitlist) noise.Add(c);
177  noise.SetNoise(true);
178  evt.put(std::make_unique<std::vector<rb::Cluster>>(1, noise));
179  return;
180  }
181 
182  // Sort the hits by time to improve the speed.
183  // WARNING: The rest of the code assumes that the hits are
184  // time sorted!
185  rb::SortByTime(hitlist);
186 
187 
188 
189  // Initalize PointManager to help with the 4D clustering,
190  // add the hit info to the point manager, and generate the
191  // initial neighbor map.
193 
194  fPtMan.Reset(hitlist.size());
195 
196  for(unsigned int i = 0; i < hitlist.size(); ++i) {
197  fPtMan.AddHit(hitlist[i]->TNS(), cal->GetTimeRes(*hitlist[i]),
198  hitlist[i]->Plane(), hitlist[i]->Cell(), hitlist[i]->PE(),
199  hitlist[i]->View());
200  }
201 
203 
204  // OPTIONAL: Clean up the neighborhood (Get off my lawn you damn kids!)
206 
207  std::unique_ptr< std::vector<rb::Cluster> >
208  clustercol(new std::vector<rb::Cluster>);
209 
210 
211 
212 
213  //
214  // Use alternate slicing methods if requested
215  //
216  int CurrentClusterID = 1;
217 
218  if(fDoNotSliceAnything) {
219 
220  // We will only put one BIG cluster into the event. I don't know
221  // if anyone uses this feature, but I grandfathered it in from
222  // the old slicer.
223 
224  rb::Cluster cluster;
225 
226  for(unsigned int i = 0; i < hitlist.size(); ++i)
227  cluster.Add(hitlist[i]);
228 
229  // Add the cluster to clustercol
230  clustercol->push_back(cluster);
231 
232  }
233  // otherwise do the standard DBSCAN slicing
234  else {
235  for(unsigned int i = 0; i < hitlist.size(); ++i) {
236  // If the point is unclassified, attempt to expand the cluster
237  if(fPtMan.fClustIDs[i] == -1) {
238  if(this->GenerateCluster(i,CurrentClusterID))
239  CurrentClusterID++;
240  } // end if point is unclassified
241  } // end loop over all hits, i
242 
243 
244 
245  //
246  // Construct the actual slices to be put into the event.
247  //
248  std::vector<rb::Cluster> clusters(CurrentClusterID);
249 
250  for(unsigned int i = 0; i < hitlist.size(); ++i)
251  clusters[fPtMan.fClustIDs[i]].Add(hitlist[i]);
252 
253  // Move hits in non-noise clusters that are too small
254  // to the noise cluster.
255  for(unsigned int i = 1; i < clusters.size(); ++i) {
256 
257  if(clusters[i].NXCell() >= fMinHitsPerView &&
258  clusters[i].NYCell() >= fMinHitsPerView) continue;
259 
260  // This cluster is too small. Put its hits into the noise
261  // cluster and clear it.
262  for(unsigned int j = 0; j < clusters[i].NCell(); ++j) {
263  art::Ptr<rb::CellHit> hit(clusters[i].Cell(j));
264  clusters[0].Add(hit);
265  } // end for j (loop over hits)
266  clusters[i].Clear();
267  } // end for i (loop over clusters)
268 
269  // Tag the noise cluster as such.
270  clusters[0].SetNoise(true);
271 
272  // Always put in the noise cluster even if it is empty.
273  clustercol->push_back(clusters[0]);
274 
275  // Add all non-empty, non-noise clusters to clustercol
276  for(unsigned int i = 1; i < clusters.size(); ++i) {
277 
278  if(clusters[i].NCell() == 0) continue;
279 
281  if(clusters[i].NXCell() == 0) V = geo::kY;
282  if(clusters[i].NYCell() == 0) V = geo::kX;
283 
284  rb::Cluster tempcluster(clusters[i].AllCells(),V);
285  clustercol->push_back(tempcluster);
286 
287  }
288 
289  } // end else { do the DBSCAN clustering }
290 
291  //
292  // Before saving the clusters, put the hits inside them into a
293  // standard order
294  //
295  for (unsigned int i=0; i<clustercol->size(); ++i) {
296  (*clustercol)[i].StandardSort();
297  }
298 
299  // put the slices into the event
300  evt.put(std::move(clustercol));
301 
302  } // end Produce function
303 
304  //----------------------------------------------------------------------
305  bool Slicer4D::GenerateCluster(unsigned int point, int ClID)
306  {
307  // In the original DBSCAN paper, this is the function they called
308  // "ExpandCluster." How this function works is most easily understood
309  // by reading that paper (don't worry, it is pretty short...)
310 
311  // Get the list of neighbors for 'point.'
312  std::vector<unsigned int> seeds = fPtMan.regionQuery(point);
313 
314  // If 'point' is not a core point, label it as noise and return.
315  if(seeds.size() < fMinPts) {
316  fPtMan.changeClustID(point, 0);
317  return false;
318  }
319  else {
320  // Assign 'point' and the neighbors of 'point' to the cluster ClID.
321  for(unsigned int seedi = 0; seedi < seeds.size(); ++seedi)
322  fPtMan.changeClustID(seeds[seedi],ClID);
323  seeds.erase(seeds.begin());
324 
325  // Loop over the neighbors of 'point' to see if they are also core
326  // points.
327  while(seeds.size() > 0) {
328  // Get the list of neighbors for the first seed point.
329  std::vector<unsigned int> result = fPtMan.regionQuery(seeds[0]);
330 
331  // If the first seed point is also a core point, add its neighbors
332  // to the list of seeds and assigne them to the cluster ClID.
333  if(result.size() >= fMinPts) {
334  for(unsigned int resulti = 0; resulti < result.size(); ++resulti) {
335  // If the point is unassigned (-1) or previously labeled as
336  // noise (0), asign it ClID.
337  if(fPtMan.fClustIDs[result[resulti]] == -1 ||
338  fPtMan.fClustIDs[result[resulti]] == 0) {
339  // Only add this point to the list of seeds if it was
340  // previously unassigned.
341  if(fPtMan.fClustIDs[result[resulti]] == -1)
342  seeds.push_back(result[resulti]);
343  fPtMan.changeClustID(result[resulti], ClID);
344  } // end if(resulti is -1 or 0)
345 
346  } // end for resulti
347 
348  } // end if result.size() >= fMinPts
349  seeds.erase(seeds.begin());
350 
351  } // end while
352  return true;
353  } // end else
354 
355  } // end GenerateCluster function
356 
357 } // end namespace slicer
358 /////////////////////////////////////////////////////////////////////////
359 
360 namespace slicer
361 {
363 }
unsigned int fMinHitsPerView
Minimum number of hits per view to make a cluster.
double GetTimeRes(rb::CellHit const &cellhit)
fhicl::ParameterSet fPSetNDOS
void Initalize(double eps, double pepen, double distpen, double opvplpen, bool usepepen)
double fDistPen
Distance (in cm) for which the distance penalty will be 1 (roughly one strong interaction length...
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
Definition: Cluster.h:161
fhicl::ParameterSet fPSetTB
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void Reset(unsigned int n, bool skip=true)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double fPEPen
PE value for which the low PE penalty will be 1.
A collection of associated CellHits.
Definition: Cluster.h:47
Slicer4D(fhicl::ParameterSet const &pset)
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
bool GenerateCluster(unsigned int point, int ClID)
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double fEpsilon
Minimum causal separation value for neighbors.
Far Detector at Ash River, MN.
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
fhicl::ParameterSet fPSetFD
int evt
Near Detector in the NuMI cavern.
const double j
Definition: BetheBloch.cxx:29
std::vector< unsigned int > regionQuery(unsigned int point)
void AddHit(double t, double tres, unsigned int p, unsigned int c, double pe, geo::View_t v)
Definition: run.py:1
PointManager fPtMan
Helper class for the 4D clustering.
void produce(art::Event &evt)
double fOpVPlPen
Allowed plane gap for hits in opposite views for which the distance penalty will be 1...
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
Helper class for Slicer4D.
Definition: PointManager.h:17
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
unsigned int fMinPts
minimum number of casual neighbors for a hit to be in a cluster
std::string fCellHitInput
Read CellHits from this input.
assert(nhit_max >=nhit_nbins)
fhicl::ParameterSet fPSetND
bool fDoNotSliceAnything
If true, one big slice is written out.
Class to help Slicer4D manage the collection of hits.
bool fUseLowPEPen
Add a causal distance penalty for PE (low PE is INCREASES the distance)
std::vector< int > fClustIDs
Assigned cluster ID for each hit.
Definition: PointManager.h:43
void changeClustID(unsigned int point, unsigned int id)
void beginRun(art::Run &run)
Encapsulate the geometry of one entire detector (near, far, ndos)
bool fApplyNeighborhoodCleanUp
scrub loner hits from the slice?
enum BeamMode string