WCHitFinderAlg.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 // This algorithm takes in Beamline raw digits information for the //
3 // NOvA Testbeam Wire Chambers and returns a vector of "good_hits" that //
4 // can be used in finding good tracks later, using a track building //
5 // algorithm of your choice. //
6 // In effect this is just the hit finding parts of the code in //
7 // WCTrackReco_module.cc and WCTrackAlg.cxx //
8 // Author: Fan Gao, fag16@pitt.edu //
9 // Based on LArIAT code lariatsoft/LArIATRecoAlg/WCHitFinderAlg.cxx //
10 // by Greg Pulliam, gkpullia@syr.edu //
11 //////////////////////////////////////////////////////////////////////////
12 
13 // framework
14 
15 // nova
17 
18 // root
19 
20 // stl
21 #include <iostream>
22 #include <cmath>
23 #include <cstdlib>
24 #include <string>
25 
26 
27 
28 namespace beamlinereco {
29  //-------------------------------------------------------------------
30  //Constructor
32  {
33  this->reconfigure(pset);
34  }
35  //-------------------------------------------------------------------
36  //Destructor
38  {
39 
40  }
41  //-------------------------------------------------------------------
43  {
44  fTime_bin_scaling = pset.get<double>("TimeBinScaling", 1.0/1208.0);
45  fWire_scaling = pset.get<double>("WireScaling", 1.0/64.0);
46  fGoodHitAveragingEps = pset.get<double>("GoodHitAveragingEps"); // 2.0
47  fDBSCANEpsilon = pset.get<float >("DBSCANEpsilon", 1.0/16.0);
48  fDBSCANMinHits = pset.get<int >("DBSCANMinHits"); // 1
49 
50  fTrack_Type = 9999;
51 
52  return;
53  }
54  //----------------------------------------------------------------------------------
55  int WCHitFinderAlg::getTrackType(std::vector<std::vector<WCHitList> > & good_hits) //HIT but need to pull parts of shouldSkipTrigger
56  {
57  int fTrack_Type = -1;
58  //Check to see if there is only a single hit on one of the WCAxes
59  bool lonely_hit_bool = false; //Single hit on at most 7 WCAxes, at least 1
60  bool unique_hit_bool = true; //Single hit on all WCAxes
61  bool missing_hit_bool = false; //Missing hit on 1 or more WCAxes
62  for (size_t iWC = 0; iWC < 4 ; ++iWC) {
63  for (size_t iAx = 0; iAx < 2 ; ++iAx) {
64  //std::cout<<"IAx= "<<iAx<<" iWC= "<<iWC<<" size= "<<good_hits.at(iWC).at(iAx).hits.size()<<std::endl;
65  if (good_hits.at(iWC).at(iAx).hits.size() == 0) missing_hit_bool = true;
66  if (good_hits.at(iWC).at(iAx).hits.size() == 1) lonely_hit_bool = true;
67  if (good_hits.at(iWC).at(iAx).hits.size() > 1) unique_hit_bool = false;
68  }
69  }
70  if (missing_hit_bool) fTrack_Type = 0;
71  else if (lonely_hit_bool && unique_hit_bool) fTrack_Type = 1;
72  else if (lonely_hit_bool && !unique_hit_bool) fTrack_Type = 2;
73  else if (!lonely_hit_bool && !unique_hit_bool) fTrack_Type = 3;
74  else { std::cout << "Unknown track condition. Check me." << std::endl; }
75 
76  return fTrack_Type;
77  }
78 
79  //---------------------------------------------------------------------------------
80  //Main function called for each trigger
82  std::vector<std::vector<WCHitList> >& good_hits,
83  bool verbose) {
84 
85  bVerbose = verbose;
86 
87  unsigned int numWCs = fGeo->NumWCs();
88  unsigned int numWCPlaneWires = fGeo->NumWCPlaneWires();
89  float centerWire = ((float)numWCPlaneWires/2.)-1;
90 
91  // create a set of buffers with 1st-Dim length of 8 (for 8 wire chamber axes: 1X, 1Y, 2X, 2Y, etc.)
92  std::vector<float> temp_buffer;
93  std::vector<std::vector<float> > hit_time_bin_vects;
94  std::vector<std::vector<float> > hit_channel_vects;
95  for (unsigned int iWCAx = 0; iWCAx < numWCs*2; ++iWCAx) {
96  hit_time_bin_vects.push_back(temp_buffer);
97  hit_channel_vects.push_back(temp_buffer);
98  }
99 
100  // collect WC digits from each trigger and fill their time and channel info into vectors
101  for (size_t iWC = 0; iWC < wcs.size(); ++iWC) {
102 
103  art::Ptr<rawdata::RawBeamlineWC> digit = wcs.at(iWC);
104 
105  std::vector<rawdata::RawBeamlineWC::WCPulse> xPulses = digit->XPulses();
106  std::vector<rawdata::RawBeamlineWC::WCPulse> yPulses = digit->YPulses();
107 
108  // put info of xPulses into vectors
109  for (size_t x_pulse = 0; x_pulse < xPulses.size(); ++x_pulse) {
110  hit_time_bin_vects.at(2*iWC).push_back(xPulses[x_pulse].Time);
111  hit_channel_vects.at(2*iWC).push_back(centerWire-xPulses[x_pulse].Channel); //Make wire number to be zero in the middle (-64 to 63)
112  }
113 
114  // put info of yPulses into vectors
115  for (size_t y_pulse = 0; y_pulse < yPulses.size(); ++y_pulse) {
116  hit_time_bin_vects.at(2*iWC+1).push_back(yPulses[y_pulse].Time);
117  hit_channel_vects.at(2*iWC+1).push_back(centerWire-yPulses[y_pulse].Channel); //Make wire number to be zero in the middle (-64 to 63)
118  }
119 
120  }
121 
122  // Initializing the good hit arrays to a default state - these clear for every event
123  // Have 2-dimensional array of hitlists: 1st Dim: WC, 2nd Dim: Axis
124  WCHitList hitList; // A vector of x-pulses or y-pulses on one WC coming from one event
125  std::vector<WCHitList> hitListAxis;
126  for (unsigned int iAx = 0; iAx < 2; ++iAx)
127  hitListAxis.push_back(hitList);
128  for (unsigned int iWC = 0; iWC < fGeo->NumWCs(); ++iWC)
129  good_hits.push_back(hitListAxis);
130 
131  //Initialize cluster time/channel buffers that are only defined for one trigger value
132  std::vector<std::vector<float> > cluster_time_buffer;
133  std::vector<std::vector<float> > cluster_wire_buffer;
134 
135  for (unsigned int iWCAx = 0; iWCAx < fGeo->NumWCs()*2 ; ++iWCAx) {
136  cluster_time_buffer.push_back(temp_buffer);
137  cluster_wire_buffer.push_back(temp_buffer);
138  }
139 
140  //Create a vector of clusters with DBSCAN, where each entry into cluster_time_buffer/cluster_wire_buffer is a different cluster
141  createClusters(hit_time_bin_vects, hit_channel_vects, cluster_time_buffer, cluster_wire_buffer);
142 
143  //Finding the hits to be used in momentum reconstruction
144  //Note here that only one hit may be accepted from a cluster. The idea is that a particle passing through the MWPC causes noise
145  //that is spatially and temporally clustered around the initial hit of the particle.
146  //In addition, if sufficiently close together, two or more good hits can be averaged
147  findGoodHits(cluster_time_buffer, cluster_wire_buffer, good_hits);
148  getTrackType(good_hits);
149  }
150 
151  //=====================================================================
152  //Use DBSCAN to create clusters of hits that are spatially and temporally close
153  void WCHitFinderAlg::createClusters(std::vector<std::vector<float> > hit_time_buffer,
154  std::vector<std::vector<float> > hit_wire_buffer,
155  std::vector<std::vector<float> > & cluster_time_buffer,
156  std::vector<std::vector<float> > & cluster_wire_buffer)
157  {
158  //Sanity check
159  if (bVerbose) { std::cout << "Length of the hit time buffer: " << hit_time_buffer.size() << std::endl;
160  std::cout << "Length of the hit wire buffer: " << hit_wire_buffer.size() << std::endl;
161  }
162 
163  //Parameter for clusters being too large
164  size_t max_hits = 10;
165 
166  //Loop through WC Axes
167  for (unsigned int iWCAx = 0; iWCAx < fGeo->NumWCs()*2; ++iWCAx) {
168 
169  //Create hit and scaled hit vectors for use in DBSCAN's clustering
170  //WCHitList hits;
171  WCHitList scaled_hits;
172  createHitAndScaledHitVectors(iWCAx, hit_time_buffer, hit_wire_buffer, scaled_hits); //(fill the above vectors with scaled hits)
173 
174  std::vector<WCHitList> cluster_list;
175  if (scaled_hits.hits.size() != 0)
176  run_DBSCAN(iWCAx, scaled_hits, cluster_list);
177 
178  //Loop through clusters and see if they are too big (pancake-like cross-talk)
179  for (size_t iClust = 0; iClust < cluster_list.size(); ++iClust) {
180  if (cluster_list.at(iClust).hits.size() > max_hits) { continue; }
181  float wire = 9999;
182  float time = 9998;
183  findLowestTimeHitInCluster (cluster_list.at(iClust), wire, time);
184  cluster_time_buffer.at(iWCAx).push_back(time/fTime_bin_scaling);
185  cluster_wire_buffer.at(iWCAx).push_back(wire/fWire_scaling);
186  }
187  }
188  }
189 
190  //=====================================================================
191  //Create hit vectors for convenient use in DBSCAN's clustering
193  const std::vector<std::vector<float> > hit_time_buffer,
194  const std::vector<std::vector<float> > hit_wire_buffer,
195  WCHitList & scaled_hit_list)
196  {
197  //Sanity Check
198  if (bVerbose) {
199  if (hit_time_buffer.size() != hit_wire_buffer.size()) { std::cout << "Error: vector size mismatch." << std::endl; }
200  if (hit_time_buffer.at(WCAx_number).size() != hit_wire_buffer.at(WCAx_number).size()) { std::cout << "Error: sub-vector size mismatch." << std::endl; }
201  }
202 
203  //For each element in the hit time buffer (for each hit)
204  for (size_t iHit = 0; iHit < hit_time_buffer.at(WCAx_number).size(); ++iHit) {
205  //Create a scaled hit and fill it with time, wire number, and whether it has been visited
206  //Also with a hit index and cluster index (-1 means noise cluster)
207  WCHit scaled_hit;
208  scaled_hit.time = hit_time_buffer.at(WCAx_number).at(iHit)*fTime_bin_scaling;
209  scaled_hit.wire = hit_wire_buffer.at(WCAx_number).at(iHit)*fWire_scaling;
210  scaled_hit.isVisited = false;
211  scaled_hit.hit_index = iHit;
212  scaled_hit.cluster_index = -1;
213  scaled_hit_list.hits.push_back(scaled_hit);
214  }
215  }
216 
217  //=====================================================================
218  //function that runs dbscan algorithm on the hit wire/time set
219  //Output is vector of same length as # of hits with labels
220  void WCHitFinderAlg::run_DBSCAN(int WCAx_number,
221  WCHitList scaled_hits,
222  std::vector<WCHitList> & cluster_list)
223  {
224  //Parameters for algorithim
225  float epsilon = fDBSCANEpsilon;
226  size_t min_hits = fDBSCANMinHits;
227 
228  //Create matrix of neighborhoods of hits, given epsilon
229  std::vector<WCHitList> neighborhood_matrix = createNeighborhoodMatrix(scaled_hits, epsilon);
230 
231  //Cluster counter
232  int cluster_counter = 0;
233  for (size_t iSH = 0; iSH < scaled_hits.hits.size(); ++iSH) {
234  //If hit has been visited, ignore it
235  if (scaled_hits.hits.at(iSH).isVisited == true) { continue; }
236  scaled_hits.hits.at(iSH).isVisited = true;
237  WCHitList neighbor_hits = regionQuery(neighborhood_matrix, scaled_hits.hits.at(iSH), scaled_hits, epsilon);
238 
239  //If there aren't enough nearest neighbors, count this hit as noise for now
240  if (neighbor_hits.hits.size() < min_hits) { scaled_hits.hits.at(iSH).cluster_index = -1; }
241 
242  //Otherwise, create a new cluster from it
243  else {
244  cluster_counter++;
245  scaled_hits.hits.at(iSH).cluster_index = cluster_counter-1;
246  expandCluster(neighborhood_matrix, scaled_hits.hits.at(iSH), scaled_hits, neighbor_hits, epsilon, min_hits, cluster_counter-1);
247  }
248  }
249 
250  //Extract cluster info from scaled_hits and put into clusters
251  for (int iClust = 0; iClust < cluster_counter; ++iClust) {
252  WCHitList cluster;
253  cluster_list.push_back(cluster);
254  }
255  for (size_t iSH = 0; iSH < scaled_hits.hits.size(); ++iSH) {
256  if (scaled_hits.hits.at(iSH).cluster_index == -1) {
257  continue;
258  }
259  cluster_list.at(scaled_hits.hits.at(iSH).cluster_index).hits.push_back(scaled_hits.hits.at(iSH));
260  }
261  }
262 
263  //=====================================================================
264  //Creating a matrix of neighbors for use in efficient DBSCAN operation
265  std::vector<WCHitList> WCHitFinderAlg::createNeighborhoodMatrix(WCHitList scaled_hits,
266  float epsilon)
267  {
268  //Final vector
269  std::vector<WCHitList> neighborhoods_vector;
270 
271  //Loop through all hits
272  for (size_t iHit = 0; iHit < scaled_hits.hits.size(); ++iHit) {
273  //Create a hit list representing all hits (including itself) that are within epsilon of that hit
274  WCHitList neighborhood_hits;
275  for (size_t iSubHit = 0; iSubHit < scaled_hits.hits.size(); ++iSubHit) {
276  float distance = pow(pow(scaled_hits.hits.at(iHit).wire-scaled_hits.hits.at(iSubHit).wire,2) + pow(scaled_hits.hits.at(iHit).time-scaled_hits.hits.at(iSubHit).time,2),0.5);
277  if (distance < epsilon) { neighborhood_hits.hits.push_back(scaled_hits.hits.at(iSubHit)); }
278  }
279  neighborhoods_vector.push_back(neighborhood_hits);
280  }
281  return neighborhoods_vector;
282  }
283 
284  //=====================================================================
285  //DBSCAN function - find hits within epsilon of the central hit
286  WCHitList WCHitFinderAlg::regionQuery(std::vector<WCHitList> neighborhood_matrix,
287  WCHit the_hit,
288  WCHitList & scaled_hits,
289  float epsilon)
290  {
291  //Create a final hit list to return
292  WCHitList neighbor_hit_list;
293 
294  for (size_t iNH = 0; iNH < neighborhood_matrix.at(the_hit.hit_index).hits.size(); ++iNH) {
295  //if (scaled_hits.hits.at(neighborhood_matrix.at(the_hit.hit_index).hits.at(iNH).hit_index).isVisited == true) { continue; }
296  neighbor_hit_list.hits.push_back(scaled_hits.hits.at(neighborhood_matrix.at(the_hit.hit_index).hits.at(iNH).hit_index));
297  }
298  return neighbor_hit_list;
299  }
300 
301  //=====================================================================
302  //DBSCAN function - once a cluster is seeded, reach out from it and find
303  //other hits that are also in this cluster
304  void WCHitFinderAlg::expandCluster(std::vector<WCHitList> neighborhood_matrix,
305  WCHit the_hit,
306  WCHitList & scaled_hits,
307  WCHitList neighbor_hits,
308  float epsilon,
309  size_t min_hits,
310  int cluster_index)
311  {
312  //Loop through neighbors (size of neighbors can increase as new hits append the neighbor hits list)
313  for (size_t iNB = 0; iNB < neighbor_hits.hits.size(); ++iNB) {
314  if (neighbor_hits.hits.at(iNB).isVisited == false) {
315  //Need to do each setting for the neighbor hit and the scaled hits cluster (scaled hits retains all info)
316  neighbor_hits.hits.at(iNB).isVisited = true;
317  scaled_hits.hits.at(neighbor_hits.hits.at(iNB).hit_index).isVisited = true;
318  WCHitList next_neighbors = regionQuery(neighborhood_matrix, neighbor_hits.hits.at(iNB), scaled_hits, epsilon);
319  //If there are enough next-neighbors for this neighbor, append the neighbor hits list
320  if (next_neighbors.hits.size() >= min_hits) {
321  for (size_t iNN = 0; iNN < next_neighbors.hits.size(); ++iNN) {
322  //If the hit already exists in the neighbors vector, then continue
323  bool isAlreadyFound = false;
324  for (size_t iHit2 = 0; iHit2 < neighbor_hits.hits.size(); ++iHit2) {
325  if (next_neighbors.hits.at(iNN).hit_index == neighbor_hits.hits.at(iHit2).hit_index) {
326  isAlreadyFound = true;
327  }
328  }
329  if (isAlreadyFound == false) {
330  //if (scaled_hits.hits.at(next_neighbors.hits.at(iNN).hit_index).isVisited == false) {
331  neighbor_hits.hits.push_back(next_neighbors.hits.at(iNN));
332  }
333  }
334  }
335  }
336 
337  //std::cout << "Scaled hits size: " << scaled_hits.hits.size() << ", searched hit index: " << neighbor_hits.hits.at(iNB).hit_index << std::endl;
338 
339  //If this neighbor hit is not yet part of a cluster, add it
340  if (neighbor_hits.hits.at(iNB).cluster_index == -1) {
341  neighbor_hits.hits.at(iNB).cluster_index = cluster_index;
342  scaled_hits.hits.at(neighbor_hits.hits.at(iNB).hit_index).cluster_index = cluster_index;
343  }
344 
345  }
346  }
347 
348  //=====================================================================
349  //Finding the one hit with the first time to represent the true
350  //passing point of the particle
352  float & wire,
353  float & time)
354  {
355  float lowest_time = 9997;
356  float lowest_time_index = 9996;
357  for (size_t iHit = 0; iHit < cluster.hits.size(); ++iHit) {
358  if (cluster.hits.at(iHit).time < lowest_time) {
359  lowest_time = cluster.hits.at(iHit).time;
360  lowest_time_index = iHit;
361  }
362  }
363  wire = cluster.hits.at(lowest_time_index).wire;
364  time = cluster.hits.at(lowest_time_index).time;
365  }
366 
367  //=====================================================================
368  //NOTE: BE CAREFUL ABOUT REFERENCING STRUCTS - MIGHT NOT REFERENCE STUFF INSIDE? CHECK
369  void WCHitFinderAlg::findGoodHits(std::vector<std::vector<float> > cluster_time_buffer,
370  std::vector<std::vector<float> > cluster_wire_buffer,
371  std::vector<std::vector<WCHitList> > & good_hits)
372  {
373  //Loop through wire chamber axes (remember, 0-7)
374  for (unsigned int iWCAx = 0; iWCAx < fGeo->NumWCs()*2 ; ++iWCAx) {
375 
376  //Sanity check
377  if (bVerbose) {
378  if (cluster_time_buffer.at(iWCAx).size() != cluster_wire_buffer.at(iWCAx).size()) {
379  std::cout << "Cluster wire/time buffer size mismatch! Error!" << std::endl;
380  return;
381  }
382  }
383  size_t number_clusters = cluster_time_buffer.at(iWCAx).size();
384 
385  //std::cout << "number clusters WCAx: " << iWCAx << " is " << number_clusters << std::endl;
386 
387  //Loop through clusters
388  for (size_t iClust = 0; iClust < number_clusters; ++iClust) {
389  float time = float(cluster_time_buffer.at(iWCAx).at(iClust));
390  float wire = float(cluster_wire_buffer.at(iWCAx).at(iClust));
391  //Convert back to wire chamber and axis for good hits.
392  int wire_chamber = int(iWCAx/2);
393  int axis = iWCAx % 2;
394  finalizeGoodHits(wire,time,good_hits.at(wire_chamber).at(axis));
395  }
396  }
397 
398  if (bVerbose) {
399  std::cout << "Number of good hits in each wire plane: ";
400  for (int iWC = 0; iWC < 4 ; ++iWC) {
401  for (int iAx = 0; iAx < 2 ; ++iAx) {
402  std::cout << good_hits.at(iWC).at(iAx).hits.size() << ", ";
403  }
404  }
405  std::cout << std::endl;
406  }
407  }
408 
409  //=====================================================================
410  //Finalize the hits and place them into the final good hit list
412  float time,
413  WCHitList & finalGoodHitList)
414  {
415  //Loop through the existing good hit list
416  for (size_t iHit = 0; iHit < finalGoodHitList.hits.size(); ++iHit) {
417  WCHit hit = finalGoodHitList.hits.at(iHit);
418  //If there are good hits close enough to each other, average them and get rid of the old hit
419  if (fabs(hit.wire-wire) < fGoodHitAveragingEps && fabs(hit.time-time) < fGoodHitAveragingEps) {
420  float average_wire = (hit.wire+wire)/2;
421  float average_time = (hit.time+time)/2;
422  finalGoodHitList.hits.erase(finalGoodHitList.hits.begin()+iHit);
423  if (bVerbose) { std::cout<<"FINALIZING!!! Wire= "<<wire<<" changing to average wire= "<<average_wire<<" and time= "<<time<<" changing to average time "<<average_time<<std::endl; }
424  wire = average_wire;
425  time = average_time;
426  }
427  }
428  //Now that we have averaged, push back the final good hit list with the (possibly) averaged hit
429  WCHit finalGoodHit;
430  finalGoodHit.wire = wire;
431  finalGoodHit.time = time;
432  finalGoodHitList.hits.push_back(finalGoodHit);
433  }
434 
435 }
std::vector< WCPulse > YPulses() const
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void createHits(const std::vector< art::Ptr< rawdata::RawBeamlineWC > > &wcs, std::vector< std::vector< WCHitList > > &good_hits, bool verbose)
unsigned int NumWCs() const
Number of WCs in beamline.
constexpr T pow(T x)
Definition: pow.h:75
void run_DBSCAN(int WCAx_number, WCHitList scaled_hits, std::vector< WCHitList > &cluster_list)
std::vector< WCPulse > XPulses() const
void finalizeGoodHits(float wire, float time, WCHitList &finalGoodHitList)
void createHitAndScaledHitVectors(int WCAx_number, const std::vector< std::vector< float > > hit_time_buffer, const std::vector< std::vector< float > > hit_wire_buffer, WCHitList &scaled_hit_list)
void findLowestTimeHitInCluster(WCHitList cluster, float &wire, float &time)
unsigned distance(const T &t1, const T &t2)
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
std::vector< WCHit > hits
void findGoodHits(std::vector< std::vector< float > > cluster_time_buffer, std::vector< std::vector< float > > cluster_wire_buffer, std::vector< std::vector< WCHitList > > &good_hits)
void hits()
Definition: readHits.C:15
T get(std::string const &key) const
Definition: ParameterSet.h:231
void expandCluster(std::vector< WCHitList > neighborhood_matrix, WCHit the_hit, WCHitList &scaled_hits, WCHitList neighbor_hits, float epsilon, size_t min_hits, int cluster_index)
void createClusters(std::vector< std::vector< float > > hit_time_buffer, std::vector< std::vector< float > > hit_wire_buffer, std::vector< std::vector< float > > &cluster_time_buffer, std::vector< std::vector< float > > &cluster_wire_buffer)
WCHitList regionQuery(std::vector< WCHitList > neighborhood_matrix, WCHit the_hit, WCHitList &scaled_hits, float epsilon)
OStream cout
Definition: OStream.cxx:6
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
std::vector< WCHitList > createNeighborhoodMatrix(WCHitList scaled_hits, float epsilon)
WCHitFinderAlg(const fhicl::ParameterSet &pset)
Definition: event.h:1
double epsilon
unsigned int NumWCPlaneWires() const
Number of wires on each WC plane.
void reconfigure(const fhicl::ParameterSet &pset)
Definition: fwd.h:28
int getTrackType(std::vector< std::vector< WCHitList > > &good_hits)