SliceMergeViews_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file MakeSlice_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: MakeSlice_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 "Utilities/func/MathUtil.h"
22 
23 #include <cmath>
24 #include <string>
25 #include <algorithm>
26 
27 
28 
29 // MergeViews header
30 namespace slicemergeviews {
31 
32  struct PtInfo {
33  double tns;
34  double plane;
35  double cell;
36  double isol;
37  double dens;
38  int nn;
39  };
40 
41  struct SliceZT {
42  double minZ;
43  double maxZ;
44  double T;
45  };
46 
48  {
49  template<class T> using Atom = fhicl::Atom<T>;
51  using Name = fhicl::Name;
52 
53  Atom<double> minCell
54  {
55  Name("MinCell"),
56  Comment("Min number of hits in each view to make merged slice")
57  };
58  Atom<double> zScale
59  {
60  Name("ZScale"),
61  Comment("Scale for adding in difference in start / end Z")
62  };
63  Atom<double> tScale
64  {
65  Name("TScale"),
66  Comment("Scale for adding in difference in Z extent of 2D slices")
67  };
68  Atom<double> tolerance
69  {
70  Name("Tolerance"),
71  Comment("Tolerance in matching X and Y slices.")
72  };
73  };
74 
76  {
77  template<class T> using Atom = fhicl::Atom<T>;
78  template<class T> using Table = fhicl::Table<T>;
80  using Name = fhicl::Name;
81 
82  Atom<std::string> slice2DInput
83  {
84  Name("Slice2DInput"),
85  Comment("Label for 2D slicer input")
86  };
88  {
89  Name("nd"),
90  Comment("Parameters for ND.")
91  };
93  {
94  Name("fd"),
95  Comment("Parameters for FD.")
96  };
98  {
99  Name("tb"),
100  Comment("Parameters for TB.")
101  };
102  };
103 
105  public:
107 
108  explicit SliceMergeViews(const Parameters& params);
109  virtual ~SliceMergeViews();
110 
111  void beginRun(art::Run& run);
112  void produce(art::Event& evt);
113  void endJob();
114 
115  void TrimNoise(rb::Cluster *cNoise, std::vector<rb::Cluster> *c2D);
116  double GetZTDist(SliceZT cX, SliceZT cY);
117  std::vector<SliceZT> GetSliceZT(std::vector<rb::Cluster> s);
118 
119  std::vector<rb::Cluster> DoMerge(rb::Cluster *cNoise,
120  std::vector<rb::Cluster> *cX,
121  std::vector<rb::Cluster> *cY);
122 
123  protected:
124 
127 
128  };
129 
130 }
131 
132 
133 
134 
135 
136 namespace slicemergeviews
137 {
138  //----------------------------------------------------------------------
140  : fParams(params())
141  {
142  produces< std::vector<rb::Cluster> >();
143  }
144 
145  //----------------------------------------------------------------------
147  { }
148 
149  //----------------------------------------------------------------------
151  {
153  switch (geom->DetId()) {
156  break;
159  break;
162  break;
163  default:
164  assert(0 && "Unknown detector");
165  }
166  }
167 
168  //----------------------------------------------------------------------
170  {
171  }
172 
173 
174 
175 
176  // Add hits from small clusters into noise slice and delete
178  std::vector<rb::Cluster> *c2D)
179  {
180  double counter = 0;
181  for (int i = c2D->size()-1; i >= 0; i--){
182  if (c2D->at(i).NXCell() < fDetectorParams.minCell() &&
183  c2D->at(i).NYCell() < fDetectorParams.minCell()){
184  for (int j = 0; j < (int)c2D->at(i).NCell(); j++){
185  cNoise->Add(c2D->at(i).Cell(j));
186  }
187  c2D->erase(c2D->begin()+i);
188  }
189  counter++;
190  }
191  }
192 
193  std::vector<SliceZT> SliceMergeViews::GetSliceZT(std::vector<rb::Cluster> s)
194  {
195  std::vector<SliceZT> info(s.size());
196  for (unsigned int i=0; i<s.size(); ++i){
197  info[i].minZ = s[i].MinZ();
198  info[i].maxZ = s[i].MaxZ();
199  info[i].T = s[i].MeanTNS();
200  }
201  return info;
202  }
203 
205  {
206  double dz1= (cX.maxZ - cY.maxZ)/fDetectorParams.zScale();
207  double dz2= (cX.minZ - cY.minZ)/fDetectorParams.zScale();
208  double dt = (cX.T - cY.T)/fDetectorParams.tScale();
209  return sqrt(dz1*dz1 + dz2*dz2 + dt*dt);
210  }
211 
212  std::vector<rb::Cluster> SliceMergeViews::DoMerge(rb::Cluster *cNoise,
213  std::vector<rb::Cluster> *cX,
214  std::vector<rb::Cluster> *cY)
215  {
216  TrimNoise(cNoise, cX);
217  TrimNoise(cNoise, cY);
218 
219  std::vector<rb::Cluster> clusts;
220  std::vector<SliceZT> cXInfo = GetSliceZT(*cX);
221  std::vector<SliceZT> cYInfo = GetSliceZT(*cY);
222 
223  while (cX->size() > 0 && cY->size() > 0){
224  double dzt[cX->size()][cY->size()];
225  int minX = 0, minY = 0;
226  for (int i = 0; i < (int)cX->size(); i++){
227  for (int j = 0; j < (int)cY->size(); j++){
228  dzt[i][j] = GetZTDist(cXInfo[i], cYInfo[j]);
229  if (dzt[i][j] < dzt[minX][minY]){
230  minX = i; minY = j;
231  }
232  }
233  }
234  if (dzt[minX][minY] < fDetectorParams.tolerance()){
235  for (int j = 0; j < (int)cY->at(minY).NCell(); j++){
236  cX->at(minX).Add(cY->at(minY).Cell(j));
237  }
238  clusts.push_back(cX->at(minX));
239  cX->erase(cX->begin()+minX);
240  cY->erase(cY->begin()+minY);
241  cXInfo.erase(cXInfo.begin()+minX);
242  cYInfo.erase(cYInfo.begin()+minY);
243  continue;
244  }
245  else break;
246  // if (cX->size() == 0 || cY->size() == 0) break;
247  }
248  while (cX->size() > 0){
249  for (int i = 0; i < (int)cX->at(0).NCell(); i++){
250  //clusts.push_back(cX->at(0));
251  cNoise->Add(cX->at(0).Cell(i));
252  }
253  cX->erase(cX->begin());
254  cXInfo.erase(cXInfo.begin());
255  }
256  while (cY->size() > 0){
257  for (int i = 0; i < (int)cY->at(0).NCell(); i++){
258  //clusts.push_back(cY->at(0));
259  cNoise->Add(cY->at(0).Cell(i));
260  }
261  cY->erase(cY->begin());
262  cYInfo.erase(cYInfo.begin());
263  }
264  return clusts;
265  }
266 
267  //----------------------------------------------------------------------
269  {
271  evt.getByLabel(fParams.slice2DInput(), slice2Dcol);
272 
273  std::vector<rb::Cluster> clustsX;
274  std::vector<rb::Cluster> clustsY;
275 
276  std::vector<rb::Cluster> clusters;
277  clusters.push_back(slice2Dcol->at(0));
278 
279  for (int i = 1; i < int(slice2Dcol->size()); i++){
280  slice2Dcol->at(i).NXCell() > 0 ?
281  clustsX.push_back(slice2Dcol->at(i)) :
282  clustsY.push_back(slice2Dcol->at(i));
283  }
284 
285 
286  std::unique_ptr< std::vector<rb::Cluster> >
287  clustercol(new std::vector<rb::Cluster>);
288 
289  std::vector<rb::Cluster> curclusts =
290  DoMerge(&clusters[0], &clustsX, &clustsY);
291  for (int j = 0; j < (int)curclusts.size(); j++){
292  clusters.push_back(curclusts[j]);
293  }
294 
295  // Tag the noise cluster as such.
296  clusters[0].SetNoise(true);
297 
298  // Always put in the noise cluster even if it is empty.
299  for (int i = 0; i < (int)clusters.size(); i++){
300  // if (clusters[i].NXCell() < 3) continue;
301  // if (clusters[i].NYCell() < 3) continue;
302  clustercol->push_back(clusters[i]);
303  }
304 
305  //
306  // Before saving the clusters, put the hits inside them into a
307  // standard order
308  //
309  for (unsigned int i=0; i<clustercol->size(); ++i) {
310  (*clustercol)[i].StandardSort();
311  }
312 
313  // put the slices into the event
314  evt.put(std::move(clustercol));
315 
316  } // end Produce function
317 
318 } // end namespace mergeviews
319 /////////////////////////////////////////////////////////////////////////
320 
321 namespace slicemergeviews
322 {
324 }
Table< SliceMergeViewsDetectorParams > fdParams
const XML_Char XML_Encoding * info
Definition: expat.h:530
double minY
Definition: plot_hist.C:9
double GetZTDist(SliceZT cX, SliceZT cY)
SliceMergeViews(const Parameters &params)
Table< SliceMergeViewsDetectorParams > ndParams
T sqrt(T number)
Definition: d0nt_math.hpp:156
A collection of associated CellHits.
Definition: Cluster.h:47
Table< SliceMergeViewsDetectorParams > tbParams
DEFINE_ART_MODULE(TestTMapFile)
void TrimNoise(rb::Cluster *cNoise, std::vector< rb::Cluster > *c2D)
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
const XML_Char * s
Definition: expat.h:262
Definition: Cand.cxx:23
std::vector< SliceZT > GetSliceZT(std::vector< rb::Cluster > s)
Far Detector at Ash River, MN.
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
Near Detector in the NuMI cavern.
const double j
Definition: BetheBloch.cxx:29
SliceMergeViewsDetectorParams fDetectorParams
Definition: run.py:1
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
std::vector< rb::Cluster > DoMerge(rb::Cluster *cNoise, std::vector< rb::Cluster > *cX, std::vector< rb::Cluster > *cY)
Encapsulate the geometry of one entire detector (near, far, ndos)