DBSlicer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: DBSlicer
3 // Module Type: filter
4 // File: DBSlicer_module.cc
5 //
6 // Generated at Tue May 13 14:24:37 2014 by Jan Zirnstein using artmod
7 // from cetpkgsupport v1_05_03.
8 ////////////////////////////////////////////////////////////////////////
9 
17 #include "fhiclcpp/ParameterSet.h"
18 
21 
22 #include <memory>
23 #include <string>
24 #include <float.h>
25 #include <math.h>
26 
27 namespace novaddt {
28  class DBSlicer;
29 }
30 
32 
33 public:
34  explicit DBSlicer(fhicl::ParameterSet const & p);
35  virtual ~DBSlicer();
36 
37  bool filter(art::Event & e) override;
38 
39 
40 private:
41 
42  // FCL parameters
45  float fDistPen;
46  float fTimePen;
47  float fTimeRes;
48  float fEps;
49  unsigned int fMinPts;
50  int fTimeGap;
51  double fMergeKnob;
52 
53  unsigned int fSliceID;
54  std::vector<std::vector< unsigned int> > fNeighbors;
55  std::vector<int> fSliceIDs;
56 
57  // Impliment DSCAN algorithm
58  void expandSlice(unsigned int Point, std::vector<unsigned int>& Neighbors, unsigned int SliceID);
59  std::vector<unsigned int> regionQuery(unsigned int Point){ return fNeighbors[Point]; }
61  double NeighborScore(const art::Handle<novaddt::HitList>& Hits, unsigned int i, unsigned int j);
62 
63  // Functions to help split and combine slices
64  void GetSingleViewSlice(const std::unique_ptr< std::vector<HitList> >& allClusters,
65  std::unique_ptr< std::vector<HitList> >& desiredClusters, std::vector<double>& desiredTMean,
66  std::vector<double>& desiredPMean, std::vector<double>& desiredTRMS,
67  std::vector<double>& desiredPRMS, daqchannelmap::DetView_TYPE desiredView);
68 
69  void MergeSlices(const std::unique_ptr< std::vector<HitList> >& XSlices,
70  std::vector<double>& XTMean, std::vector<double>& XPMean, std::vector<double>& XTRMS,
71  std::vector<double>& XPRMS, const std::unique_ptr< std::vector<HitList> >& YSlices,
72  std::vector<double>& YTMean, std::vector<double>& YPMean, std::vector<double>& YTRMS,
73  std::vector<double>& YPRMS, std::unique_ptr< std::vector<HitList> >& CombSlices);
74 
75 };
76 
77 
79 :
80 fHLabel (p.get< std::string >("hits_label")),
81 fInstance (p.get< std::string >("instance")),
82 fDistPen (p.get< float >("Distance_Penalty" )),
83 fTimePen (p.get< float >("Timing_Penalty" )),
84 fTimeRes (p.get< float >("Timing_Resolution" )),
85 fEps (p.get< float >("Epsilon" )),
86 fMinPts (p.get< unsigned int >("MinPts" )),
87 fTimeGap (p.get< int >("TimeGap")),
88 fMergeKnob(p.get< double >("ViewMergeKnob"))
89 {
90  produces<std::vector<HitList>>("DBSCANSlicesComb");
91  produces<std::vector<HitList>>("DBSCANSlicesX");
92  produces<std::vector<HitList>>("DBSCANSlicesY");
93 }
94 
96 {
97  // Clean up dynamic memory and other resources here.
98 }
99 
100 void novaddt::DBSlicer::expandSlice(unsigned int Point, std::vector<unsigned int>& Neighbors, unsigned int SliceID)
101 {
102  // Impliment expandSlice
103  fSliceIDs[Point] = SliceID;
104  for(auto it:Neighbors){
105  if(fSliceIDs[it] == -5) {
106  fSliceIDs[it] = -1;
107  std::vector<unsigned int> MoreNeighbors = regionQuery(it);
108  if(MoreNeighbors.size() >= fMinPts) // Add the two neighborhoods
109  fNeighbors[Point].insert(fNeighbors[Point].end(), MoreNeighbors.begin(), MoreNeighbors.end());
110  // Add hit to slice if it doesn't belong to another one already
111  if(fSliceIDs[it] < 0) fSliceIDs[it] = SliceID;
112  }
113  }
114 }
115 
116 double novaddt::DBSlicer::NeighborScore(const art::Handle<novaddt::HitList>& Hits, unsigned int i, unsigned int j)
117 {
118  // Ensure no merging of views at this stage
119  if(Hits->at(i).View() != Hits->at(j).View()) return DBL_MAX;
120 
121  // Use FHICL parameter to determine how tight to cluster in space, numbers are cell dimensions in cm
122  float DistPenPlane= fDistPen/6.2;
123  float TimePenalty = fTimePen/fTimeRes;
124 
125  //double dC = abs(Hits->at(i).Cell() - Hits->at(j).Cell());
126  double dT = fabs(Hits->at(i).TDC() - Hits->at(j).TDC()) * 15.625; // Conversion to TNS, for coarse timing
127  double dP = abs(Hits->at(i).Plane() - Hits->at(j).Plane());
128 
129  double Z = dP/DistPenPlane;
130  double T = dT/TimePenalty;
131 
132  return Z*Z+T*T;
133 }
134 
136 {
137  fNeighbors.clear();
138  fSliceIDs.clear();
139  unsigned int size = Hits->size();
140 
141  // Initialize IDs to unset
142  fSliceIDs.resize(size,-5);
143 
144  // Find Neighbors
145  fNeighbors.resize(size);
146  for(unsigned int it=0; it<size; ++it) {
147  fNeighbors[it].push_back(it);
148  for(unsigned int jt=it+1; jt<size; ++jt) {
149  // If the hits are too far removed in time, don't include any more points
150  if(Hits->at(jt).TDC() - Hits->at(it).TDC() > (uint64_t) fTimeGap) break;
151  double score = NeighborScore(Hits, it, jt);
152  if(score <= fEps) {
153  fNeighbors[it].push_back(jt);
154  fNeighbors[jt].push_back(it);
155  } // end if
156  } // end second index
157  } // end first index
158 }
159 
160 void novaddt::DBSlicer::GetSingleViewSlice(const std::unique_ptr< std::vector<HitList> >& allClusters,
161  std::unique_ptr< std::vector<HitList> >& desiredClusters, std::vector<double>& desiredTMean,
162  std::vector<double>& desiredPMean, std::vector<double>& desiredTRMS,
163  std::vector<double>& desiredPRMS, daqchannelmap::DetView_TYPE desiredView)
164 {
165  for(auto it=allClusters->begin(); it!=allClusters->end(); ++it){
166  if((*it)[0].View().val != desiredView) continue;
167  desiredClusters->push_back(*it);
168  double tsum = 0.0;
169  double psum = 0.0;
170  double tsum2 = 0.0;
171  double psum2 = 0.0;
172  for(auto jt:*it){
173  tsum += jt.TDC();
174  tsum2 += jt.TDC()*jt.TDC();
175  psum += jt.Plane();
176  psum2 += jt.Plane()*jt.Plane();
177  desiredTMean.push_back(tsum/it->size());
178  desiredPMean.push_back(psum/it->size());
179  desiredTRMS.push_back(sqrt(tsum2/it->size()));
180  desiredPRMS.push_back(sqrt(psum2/it->size()));
181  }
182  }
183 }
184 
185 void novaddt::DBSlicer::MergeSlices(const std::unique_ptr< std::vector<HitList> >& XSlices,
186  std::vector<double>& XTMean, std::vector<double>& XPMean, std::vector<double>& XTRMS,
187  std::vector<double>& XPRMS, const std::unique_ptr< std::vector<HitList> >& YSlices,
188  std::vector<double>& YTMean, std::vector<double>& YPMean, std::vector<double>& YTRMS,
189  std::vector<double>& YPRMS, std::unique_ptr< std::vector<HitList> >& CombSlices)
190 {
191  // This function will have to decide which Cluster from the X View to combine with which Cluster from the Y view
192  for(unsigned int it = 0; it<XSlices->size(); ++it){
193  // Initialize a few things to keep track of, the minimum score and which Slice it corresponds to
194  double minscore = DBL_MAX;
195  unsigned int YSlice = 0;
196  // A vector to mark which YSlices we have already used to combine with XSlices
197  std::vector<unsigned int> used(YSlices->size(),0);
198  for(unsigned int jt=0; jt<YSlices->size(); ++jt){
199  double score = fabs(XTMean[it]-YTMean[jt])/sqrt(XTRMS[it]*XTRMS[it]+YTRMS[jt]*YTRMS[jt]) +
200  fabs(XPMean[it]-YPMean[jt])/sqrt(XPRMS[it]*XPRMS[it]+YPRMS[jt]*YPRMS[jt]);
201  if(score < minscore && used[jt] == 0){
202  minscore = score;
203  YSlice = jt;
204  }
205  }
206  // If there is a YSlice that matches the XSlice, put them both into the Combslices, and mark the slice as used
207  if(minscore < fMergeKnob){
208  novaddt::HitList tmphits(XSlices->at(it));
209  tmphits.insert(tmphits.end(), YSlices->at(YSlice).begin(), YSlices->at(YSlice).end());
210  CombSlices->push_back(tmphits);
211  used[YSlice] = 1;
212  }
213  }
214 }
215 
217 {
219  e.getByLabel(fHLabel, fInstance, hits);
220 
221  // Now use that product to fill the things we are actually going to put in the event,
222  // view split 2D slices and a view combined slicelist of hitlists
223  std::unique_ptr< std::vector<HitList> >productcomb(new std::vector<HitList>);
224  std::unique_ptr< std::vector<HitList> >productx(new std::vector<HitList>);
225  std::unique_ptr< std::vector<HitList> >producty(new std::vector<HitList>);
226 
227  if(hits->empty()) {
228  e.put(std::move(productcomb),"DBSCANSlicesComb");
229  e.put(std::move(productx),"DBSCANSlicesX");
230  e.put(std::move(producty),"DBSCANSlicesY");
231  return false;
232  }
233  InitializeLists(hits);
234 
235  // Meat and Potatoes of DBSCAN
236  fSliceID = 0;
237  for(unsigned int it = 0; it < hits->size(); ++it) {
238  if(fSliceIDs[it] != -5) continue;
239  // Mark hit as "visited"
240  fSliceIDs[it] = -1;
241  std::vector<unsigned int> Neighbors = regionQuery(it);
242  // Bail if the hit doesn't have enough neighbors
243  if(Neighbors.size() < fMinPts) fSliceIDs[it] = 0;
244  else {
245  // Make a new slice and figure out who goes into the slice
246  fSliceID++;
247  expandSlice(it,Neighbors,fSliceID);
248  } // end else
249  } //for all unvisited hits
250 
251  // Initialize the product
252  std::unique_ptr< std::vector<HitList> >product(new std::vector<HitList>);
253 
254  // Resize it to the number of HitLists we'll be producing
255  product->resize(fSliceID);
256  // Fill the product HitLists with the corresponding DAQHits from the original Handle
257  for(unsigned int it = 0; it<fSliceIDs.size(); ++it)
258  product->at(fSliceIDs[it]).push_back(hits->at(it));
259 
260  // Containers for storing information to be used in the merging step
261  std::vector<double> XMeanTime, XMeanPlane, XRMSTime, XRMSPlane;
262  std::vector<double> YMeanTime, YMeanPlane, YRMSTime, YRMSPlane;
263 
264  GetSingleViewSlice(product, productx, XMeanTime, XMeanPlane, XRMSTime, XRMSPlane, daqchannelmap::X_VIEW);
265  GetSingleViewSlice(product, producty, YMeanTime, YMeanPlane, YRMSTime, YRMSPlane, daqchannelmap::Y_VIEW);
266  MergeSlices(productx, XMeanTime, XMeanPlane, XRMSTime, XRMSPlane,
267  producty, YMeanTime, YMeanPlane, YRMSTime, YRMSPlane, productcomb);
268 
269  //put the hits into the event for later use
270  e.put(std::move(productcomb),"DBSCANSlicesComb");
271  e.put(std::move(productx),"DBSCANSlicesX");
272  e.put(std::move(producty),"DBSCANSlicesY");
273  return true;
274 
275 }
276 
void InitializeLists(const art::Handle< novaddt::HitList > &Hits)
DBSlicer(fhicl::ParameterSet const &p)
set< int >::iterator it
void MergeSlices(const std::unique_ptr< std::vector< HitList > > &XSlices, std::vector< double > &XTMean, std::vector< double > &XPMean, std::vector< double > &XTRMS, std::vector< double > &XPRMS, const std::unique_ptr< std::vector< HitList > > &YSlices, std::vector< double > &YTMean, std::vector< double > &YPMean, std::vector< double > &YTRMS, std::vector< double > &YPRMS, std::unique_ptr< std::vector< HitList > > &CombSlices)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int fSliceID
double NeighborScore(const art::Handle< novaddt::HitList > &Hits, unsigned int i, unsigned int j)
bool filter(art::Event &e) override
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string fInstance
void abs(TH1 *hist)
DEFINE_ART_MODULE(TestTMapFile)
std::vector< unsigned int > regionQuery(unsigned int Point)
void expandSlice(unsigned int Point, std::vector< unsigned int > &Neighbors, unsigned int SliceID)
std::vector< int > fSliceIDs
Identifier for the Y measuring view of the detector (side)
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void hits()
Definition: readHits.C:15
Identifier for the X measuring view of the detector (top)
const double j
Definition: BetheBloch.cxx:29
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T product(std::vector< T > dims)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
DetView_TYPE
Types of Detector View.
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
unsigned int fMinPts
std::vector< std::vector< unsigned int > > fNeighbors
void GetSingleViewSlice(const std::unique_ptr< std::vector< HitList > > &allClusters, std::unique_ptr< std::vector< HitList > > &desiredClusters, std::vector< double > &desiredTMean, std::vector< double > &desiredPMean, std::vector< double > &desiredTRMS, std::vector< double > &desiredPRMS, daqchannelmap::DetView_TYPE desiredView)