DBScanAlg_DUNE35t.cxx
Go to the documentation of this file.
1 /**
2  * @file Cluster3D_module.cc
3  *
4  * @brief Producer module to create 3D clusters from input hits
5  *
6  */
7 
8 // Framework Includes
9 #include "cetlib/search_path.h"
10 #include "cetlib/cpu_timer.h"
12 
14 
15 // LArSoft includes
18 #include "lardata/RecoObjects/Cluster3D.h"
24 
25 // std includes
26 #include <string>
27 #include <functional>
28 #include <iostream>
29 #include <memory>
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 // implementation follows
33 
34 namespace lar_cluster3d {
35 
37 {
38  this->reconfigure(pset);
39 }
40 
41 //------------------------------------------------------------------------------------------------------------------------------------------
42 
44 {
45 }
46 
47 //------------------------------------------------------------------------------------------------------------------------------------------
48 
50 {
51  m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
52  m_minPairPts = pset.get<size_t>("MinPairPts", 2 );
53  m_timeAdvanceGap = pset.get<double>("TimeAdvanceGap", 50.);
54  m_numSigmaPeakTime = pset.get<double>("NumSigmaPeakTime", 5.);
55  m_EpsMaxDist = pset.get<double>("EpsilonDistanceDBScan", 5.);
56 
58 
59  m_geometry = &*geometry;
60 
61  m_timeVector.resize(NUMTIMEVALUES, 0.);
62 
63 
64 }
65 
68  reco::HitPairListPtr& curCluster,
69  size_t minPts) const
70 {
71  // This is the main inside loop for the DBScan based clustering algorithm
72  //
73  // Add the current hit to the current cluster
74  epsVecItr->second.first.setInCluster();
75  curCluster.push_back(epsVecItr->first);
76 
77  // Get the list of points in this hit's epsilon neighborhood
78  // Note this is a copy so we can modify locally
79  EpsPairNeighborhoodList epsNeighborhoodList = epsVecItr->second.second;
80 
81  while(!epsNeighborhoodList.empty())
82  {
83  // Dereference the point so we can see in the debugger...
84  const reco::ClusterHit3D* neighborPtr = epsNeighborhoodList.front();
85 
86  // Use that to look up the iterator in the general neighborhood map
87  EpsPairNeighborhoodMapVec::iterator curPtEpsVecItr = epsNeighborhoodMapVec.begin();
88 
89  std::advance(curPtEpsVecItr, neighborPtr->getID());
90 
91  // If we've not been here before then take action...
92  if (!curPtEpsVecItr->second.first.visited())
93  {
94  curPtEpsVecItr->second.first.setVisited();
95 
96  // If this epsilon neighborhood of this point is large enough then add its points to our list
97  if (curPtEpsVecItr->second.first.getCount() >= minPts)
98  {
99  // Plan is to loop through the hits in this point's neighborhood and add them to our list provided
100  // they have not already been added, or are part of a cluster, etc.
101  // So, get the list of points in the neighborhood
102  EpsPairNeighborhoodList& locList = curPtEpsVecItr->second.second;
103 
104  // And loop through them...
105  for(EpsPairNeighborhoodList::iterator hitItr = locList.begin(); hitItr != locList.end(); hitItr++)
106  {
107  epsNeighborhoodList.push_back(*hitItr);
108  }
109  }
110  }
111 
112  // If the point is not yet in a cluster then we now add
113  if (!curPtEpsVecItr->second.first.inCluster())
114  {
115  curPtEpsVecItr->second.first.setInCluster();
116  curCluster.push_back(curPtEpsVecItr->first);
117  }
118 
119  epsNeighborhoodList.pop_front();
120  }
121 
122  return;
123 }
124 
125 bool SetPositionOrder(const std::unique_ptr<reco::ClusterHit3D>& left, const std::unique_ptr<reco::ClusterHit3D>& right)
126 {
127  // The positions in the Y-Z plane are quantized so take advantage of that for ordering
128  // First check that we are in the same "bin" in the z direction
129  if (left->getHits().back()->getHit().WireID().Wire == right->getHits().back()->getHit().WireID().Wire) // These hits are "on the same w wire"
130  {
131  // We can use the U wires as a proxy for ordering hits in increasing Y
132  // where we remember that as Y increases the U wire number decreases
133  if (left->getHits().front()->getHit().WireID().Wire == right->getHits().front()->getHit().WireID().Wire)
134  {
135  // In the time direction we are not quantized at a large enough level to check
136  return left->getX() < right->getX();
137  }
138 
139  return left->getHits().front()->getHit().WireID().Wire > right->getHits().front()->getHit().WireID().Wire;
140  }
141 
142  return left->getHits().back()->getHit().WireID().Wire < right->getHits().back()->getHit().WireID().Wire;
143 }
144 
145 
146 //This function sorts by Z first. It is the default, and is also used if we have particles passing through the EW
147 //scintillation counters on the exterior of the detector
148 bool SetPositionOrderZ( const std::unique_ptr<reco::ClusterHit3D>& left, const std::unique_ptr<reco::ClusterHit3D>& right )
149 {
150  //Sort by Z first, then by Y, then by X.
151  if( left->getZ() == right->getZ() ){
152  if( left->getY() == right->getY() ){
153  return left->getX() < right->getX();
154  }
155  return left->getY() < right->getY();
156  }
157  return left->getZ() < right->getZ();
158 }
159 
160 
161 //This function sorts by X first. It is used if we have particles passing through the NS scintillation counters on the
162 //exterior of the detector.
163 bool SetPositionOrderX( const std::unique_ptr<reco::ClusterHit3D>& left, const std::unique_ptr<reco::ClusterHit3D>& right )
164 {
165  //Sort by X first, then by Z, then by Y.
166  if( left->getX() == right->getX() ){
167  if( left->getZ() == right->getZ() ){
168  return left->getY() < right->getY();
169  }
170  return left->getZ() < right->getZ();
171  }
172  return left->getX() < right->getX();
173 }
174 
175 //This function sorts by Y first. It is used if we have particles passing through the top scintillation counters above
176 //the detector.
177 bool SetPositionOrderY( const std::unique_ptr<reco::ClusterHit3D>& left, const std::unique_ptr<reco::ClusterHit3D>& right )
178 {
179  //Sort by Y first, then by Z, then by X.
180  if( left->getY() == right->getY() ){
181  if( left->getZ() == right->getZ() ){
182  return left->getX() < right->getX();
183  }
184  return left->getZ() < right->getZ();
185  }
186  return left->getY() < right->getY();
187 }
188 
189 
190 
192 {
193  return (*left)->getAvePeakTime() < (*right)->getAvePeakTime();
194 }
195 
197 {
199  {
200  // Watch out for the case where two clusters can have the same number of hits!
201  if (left->second.size() == right->second.size())
202  return left->first < right->first;
203 
204  return left->second.size() > right->second.size();
205  }
206 };
207 
209  ViewToWireToHitSetMap& viewToWiretoHitSetMap,
210  HitPairList& hitPairList,
211  reco::HitPairClusterMap& hitPairClusterMap)
212 {
213  /**
214  * @brief Driver for processing input 2D hits, transforming to 3D hits and building lists
215  * of associated 3D hits (candidate 3D clusters)
216  */
217  cet::cpu_timer theClockMakeHits;
218  cet::cpu_timer theClockBuildNeighborhood;
219  cet::cpu_timer theClockDBScan;
220 
221  m_timeVector.resize(NUMTIMEVALUES, 0.);
222 
223  if (m_enableMonitoring) theClockMakeHits.start();
224 
225 
226  // The first task is to take the lists of input 2D hits (a map of view to sorted lists of 2D hits)
227  // and then to build a list of 3D hits to be used in downstream processing
228  size_t numHitPairs = BuildHitPairMap(viewToHitVectorMap, viewToWiretoHitSetMap, hitPairList);
229 
230  //Monitoring activity
231  if (m_enableMonitoring)
232  {
233  theClockMakeHits.stop();
234  theClockBuildNeighborhood.start();
235  }
236 
237  // The container of pairs and those in each pair's epsilon neighborhood
238  EpsPairNeighborhoodMapVec epsPairNeighborhoodMapVec;
239  epsPairNeighborhoodMapVec.resize(numHitPairs, EpsPairNeighborhoodMapPair(0,EpsPairNeighborhoodPair())); //<-- initialize too!
240 
241  // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
242  // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
243  // The following call does this work
244  BuildNeighborhoodMap(hitPairList, epsPairNeighborhoodMapVec);
245 
246  //Monitoring activity
247  if (m_enableMonitoring)
248  {
249  theClockBuildNeighborhood.stop();
250  theClockDBScan.start();
251  }
252 
253  // With the neighborhood built we can now form "clusters" with DBScan
254  int pairClusterIdx(0);
255 
256  // Clear the cluster list just for something to do here...
257  hitPairClusterMap.clear();
258 
259  // Ok, here we go!
260  // We can simply iterate over the map we have just built to loop through the hits "simply"
261  for(EpsPairNeighborhoodMapVec::iterator epsPairVecItr = epsPairNeighborhoodMapVec.begin();
262  epsPairVecItr != epsPairNeighborhoodMapVec.end();
263  epsPairVecItr++)
264  {
265  // Skip the null entries (they were filtered out)
266  if (!epsPairVecItr->first) continue;
267 
268  // If this hit has been "visited" already then skip
269  if (epsPairVecItr->second.first.visited()) continue;
270 
271  // We are now visiting it so mark it as so
272  epsPairVecItr->second.first.setVisited();
273 
274  // Check that density is sufficient
275  if (epsPairVecItr->second.first.getCount() < m_minPairPts)
276  {
277  epsPairVecItr->second.first.setNoise();
278  }
279  else
280  {
281  // "Create" a new cluster and get a reference to it
282  reco::HitPairListPtr& curCluster = hitPairClusterMap[pairClusterIdx++];
283 
284  // expand the cluster
285  expandCluster(epsPairNeighborhoodMapVec, epsPairVecItr, curCluster, m_minPairPts);
286  }
287  }
288 
289  //Monitoring activity
290  if (m_enableMonitoring)
291  {
292  theClockDBScan.stop();
293 
294  m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time();
295  m_timeVector[BUILDHITTOHITMAP] = theClockBuildNeighborhood.accumulated_real_time();
296  m_timeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
297  }
298 
299  mf::LogDebug("Cluster3D") << ">>>>> DBScan done, found " << hitPairClusterMap.size() << " clusters" << std::endl;
300 
301  return;
302 }
303 
304 
305 //This is second-least restrictive: U+W+V required, now with a corrected version of nearestwireid
306 size_t DBScanAlg_DUNE35t::BuildHitPairMap(ViewToHitVectorMap& viewToHitVectorMap, ViewToWireToHitSetMap& viewToWireToHitSetMap, HitPairList& hitPairList) const
307 {
308  /**
309  * @brief Given input 2D hits, build out the lists of possible 3D hits
310  *
311  * The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
312  * However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
313  * be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
314  * and then look for the match in the V plane. In the event we don't find the match in the V plane then we
315  * will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.
316  */
317 
318 
319  // Should we set a minimum total charge for a hit?
320  size_t totalNumHits(0);
321  size_t hitPairCntr(0);
322  double minCharge[] = {0., 0., 0.};
323 
324  //-********************************************************************************
325  // First job is to find all pairs of hits in different views which are "consistent"
326  // Start loops over rows in HitVectorMap
327  ViewToHitVectorMap::iterator mapItrU = viewToHitVectorMap.find(geo::kU);
328 
329  if (mapItrU != viewToHitVectorMap.end())
330  {
331  HitVector& hitVectorU = mapItrU->second;
332 
333  totalNumHits += hitVectorU.size();
334 
335  ViewToHitVectorMap::iterator mapItrW = viewToHitVectorMap.find(geo::kW);
336 
337  if (mapItrW != viewToHitVectorMap.end())
338  {
339  HitVector& hitVectorW = mapItrW->second;
340 
341  totalNumHits += hitVectorW.size();
342 
343  if (viewToHitVectorMap.find(geo::kV) != viewToHitVectorMap.end())
344  totalNumHits += viewToHitVectorMap[geo::kV].size();
345 
346  // Take advantage that hits are sorted in "start time order"
347  // Set the inner loop iterator before starting loop over outer hits
348  HitVector::iterator hitVectorWStartItr = hitVectorW.begin();
349 
350  // Now we loop over the hits in these two layers
351  for (HitVector::iterator hitItrU = hitVectorU.begin(); hitItrU != hitVectorU.end(); hitItrU++)
352  {
353 
354  const reco::ClusterHit2D* hitPtrU = *hitItrU;
355 
356  if (hitPtrU->getHit().Integral() < minCharge[hitPtrU->getHit().View()]) continue;
357 
358 
359  // This will be used in each loop so dereference the peak time here
360  double hitUPeakTime = hitPtrU->getTimeTicks();
361 
362  // Inner loop is over hits within the time range of the outer hit
363  for (HitVector::iterator hitItrW = hitVectorWStartItr; hitItrW != hitVectorW.end(); hitItrW++)
364  {
365  const reco::ClusterHit2D* hitPtrW = *hitItrW;
366 
367  //Check to make sure that this hit is in the same TPC as the previous hit
368  if( hitPtrW->getHit().WireID().TPC != hitPtrU->getHit().WireID().TPC ) continue;
369 
370  if (hitPtrW->getHit().Integral() < minCharge[hitPtrW->getHit().View()]) continue;
371 
372  // Hits are sorted in "peak time order" which we can take advantage of to try to speed
373  // the loops. Basically, we can compare the peak time for the outer loop hit against
374  // the current hit's peak time and if outside the range of interest we can take action.
375  // The range of interest is eyeballed but is meant to account for large pulse widths
376  // Start by dereferencing the inner hit peak time
377  double hitWPeakTime = hitPtrW->getTimeTicks();
378 
379  // If the outer loop's peak time is well past the current hit's then
380  // we should advance the inner loop's start iterator and keep going
381  if (hitUPeakTime > hitWPeakTime + m_timeAdvanceGap)
382  {
383 
384 
385  hitVectorWStartItr++;
386  continue;
387  }
388 
389  // If the inner loop hit start time is past the end of the outer's end time, then break out loop
390  if (hitUPeakTime + m_timeAdvanceGap < hitWPeakTime) break;
391 
392  // We have a candidate hit pair combination, try to make a hit
393  reco::ClusterHit3D pair = makeHitPair(hitPtrU, hitPtrW);
394 
395 
396  // The sign of success here is that the average peak time of the combined hits is > 0
397  // (note that when hits are combined the first window offset is accounted for)
398  if (pair.getAvePeakTime() > 0.)
399  {
400  //bool hitInVNotFound(true);
401 
402 
403  //I need to know the plane ID to find the nearest wire in this TPC. In order to do that, I need to
404  //loop through the viewToHitVector map and find the plane ID corresponding to a hit in this plane and this TPC
405  size_t thePlane = 9999;
406  for( HitVector::iterator hitVect_iter = viewToHitVectorMap[geo::kV].begin(); hitVect_iter != viewToHitVectorMap[geo::kV].end(); ++hitVect_iter ){
407  if( (*(hitVect_iter))->getHit().WireID().TPC == hitPtrW->getHit().WireID().TPC ){
408  thePlane = (*(hitVect_iter))->getHit().WireID().Plane;
409  break;
410  }
411  }
412  geo::PlaneID thePlaneID(0,hitPtrW->getHit().WireID().TPC,thePlane);
413 
414  // Recover the WireID nearest in the V plane to the position of the pair
415  //REL Use plane ids here to avoid TPC ambiguity introduced by just using view type
416  const geo::WireID wireIDV = NearestWireID_mod(pair.getPosition(), thePlaneID);
417 
418  //Some debug printing
419  if( wireIDV.TPC != hitPtrW->getHit().WireID().TPC )
420  std::cout << "TPC of third (nearest) wire is not the same as the TPC of the first two." << std::endl;
421  if( hitPtrU->getHit().WireID().TPC != hitPtrW->getHit().WireID().TPC )
422  std::cout << "TPC of wire 1 is not TPC of wire 2." << std::endl;
423 
424  // We believe the code that returns the ID is offset by one
425  //Get the hits associated with the nearest wire
426  WireToHitSetMap::iterator wireToHitSetMapVItr = viewToWireToHitSetMap[geo::kV].find(wireIDV.Wire);
427 
428  if (wireToHitSetMapVItr != viewToWireToHitSetMap[geo::kV].end())
429  {
430  const reco::ClusterHit2D* hit2DV = FindBestMatchingHit(wireToHitSetMapVItr->second, pair, m_numSigmaPeakTime*pair.getSigmaPeakTime());
431 
432 
433  // If a V hit found then it should be straightforward to make the triplet
434 
435  if (hit2DV){
436  if( hitPtrW->getHit().WireID().TPC == hit2DV->getHit().WireID().TPC ){
437 
438  reco::ClusterHit3D triplet = makeHitTriplet(pair, hit2DV);
439 
440  if (triplet.getAvePeakTime() > 0.)
441  {
442 
443  triplet.setID(hitPairCntr++);
444  hitPairList.emplace_back(std::unique_ptr<reco::ClusterHit3D>(new reco::ClusterHit3D(triplet)));
445  // goodBool = true;
446  }
447 
448  //hitInVNotFound = false;
449  }
450  }
451  }
452 
453  }
454  }
455  }
456  }
457  }
458  hitPairList.sort(SetPositionOrderZ);
459  return hitPairList.size();
460 
461 }
462 
463 
464 
465 
466 
467 //REL modified
469  EpsPairNeighborhoodMapVec& epsPairNeighborhoodMapVec )const
470 
471 {
472  size_t consistentPairsCnt = 0;
473  size_t pairsChecked = 0;
474 
475  for (HitPairList::const_iterator pairItrO = hitPairList.begin(); pairItrO != hitPairList.end(); pairItrO++){
476 
477  const reco::ClusterHit3D* hitPairO = (*pairItrO).get();
478  const size_t hitPairOID = hitPairO->getID();
479 
480  // Need to initialize the "first" part of the vector pseudo map
481  epsPairNeighborhoodMapVec[hitPairOID].first = hitPairO;
482 
483  // Get reference to the list for this hit so we don't look it up inside the loop
484  EpsPairNeighborhoodPair& hitPairOPair(epsPairNeighborhoodMapVec[hitPairOID].second);
485 
486  HitPairList::const_iterator pairItrI = pairItrO;
487 
488  std::map<int, std::pair<double, const reco::ClusterHit3D*> > bestTripletMap;
489 
490  // Get the X,Y, and Z for this triplet
491  double pairO_X = hitPairO->getPosition()[0];
492  double pairO_Y = hitPairO->getPosition()[1];
493  double pairO_Z = hitPairO->getPosition()[2];
494 
495  // Set maximums
496  double maxDist = m_EpsMaxDist; //REL 4.0; //Was 2, then 4 REL
497 
498  //Loop over second iterator
499  while (++pairItrI != hitPairList.end())
500  {
501  const reco::ClusterHit3D* hitPairI = (*pairItrI).get();
502 
503  //OLD WAY
504  // Note that the hits have been sorted by z, then y and then in x
505  // Translation: hits are sorted by W wire, then in Y (increasing u, decreasing v)
506  // and then in time.
507 
508  //NEW WAY
509  // Hits have been sorted by Z position, then by Y, then by X (each in ascending order).
510  // There is no reference to wire number here, since we need to consider clustering
511  // across TPCs.
512 
513  // Get the X,Y, and Z for this triplet
514  double pairI_X = hitPairI->getPosition()[0];
515  double pairI_Y = hitPairI->getPosition()[1];
516  double pairI_Z = hitPairI->getPosition()[2];
517 
518  //Sorting allows us to break this loop after the maxDist is surpassed in z
519  if( pairI_Z - pairO_Z > maxDist ) break;
520 
521  //Form the distance so we can check ranges
522  double distance = pow(
523  pow(pairO_X-pairI_X,2) +
524  pow(pairO_Y-pairI_Y,2) +
525  pow(pairO_Z-pairI_Z,2),0.5);
526 
527  // If we have passed the 3d distance then we are done with the loop
528  if( distance > maxDist )continue;
529 
530 
531  // This is the tight constraint on the hits
532  if (consistentPairs(hitPairO, hitPairI))
533  {
534  double bestBin = distance;
535  double bestDist = 10000.;
536 
537  if (bestTripletMap.find(bestBin) != bestTripletMap.end()) //Look at a given case of error (a given bin) and pull out the best (shortest) distance
538  {
539  bestDist = bestTripletMap[bestBin].first;
540  }
541 
542  double newDist = fabs(hitPairI->getX() - hitPairO->getX());
543 
544  // This is an attempt to "prefer" triplets over pairs
545  if (hitPairI->getHits().size() < 3) newDist += 25.;
546 
547  if (newDist < bestDist) bestTripletMap[bestBin] = std::pair<double, const reco::ClusterHit3D*>(newDist, hitPairI); //Select the pair of hits that is closest in time
548 
549  }
550 
551  bestTripletMap[distance] = std::pair<double,const reco::ClusterHit3D*>(distance,hitPairI);
552 
553  }
554 
555  for(const auto& bestMapItr : bestTripletMap)
556  {
557  const reco::ClusterHit3D* hitPairI(bestMapItr.second.second);
558 
559  hitPairOPair.first.incrementCount();
560  hitPairOPair.second.emplace_back(hitPairI);
561 
562  epsPairNeighborhoodMapVec[hitPairI->getID()].second.first.incrementCount();
563  epsPairNeighborhoodMapVec[hitPairI->getID()].second.second.emplace_back(hitPairO);
564 
565  consistentPairsCnt++;
566  }
567  }
568 
569  mf::LogDebug("Cluster3D") << "Consistent pairs: " << consistentPairsCnt << " of " << pairsChecked << " checked." << std::endl;
570 
571  return consistentPairsCnt;
572 
573 }
574 
575 
576 
577 
579 {
580  // Most likely this will fail
581  bool consistent(false);
582 
583  // Strategy: We assume that our mission is to check only in the time direction
584  // The idea is to check each 2D hit individually but to account for
585  // skipped wires if necessary. The criterion applied here is that 2
586  // of 3 times in the views must match the "tight" tolerance with the
587  // largest deltaT must still in "in range"
588  // Note that the input values can be either pairs of hits in U-W plane
589  // or a triplet of U-V-W 2D hits
590  // Start with U plane
591  std::vector<bool> matchedHits = {false, false, false};
592  double maxDeltaT(0.);
593  double maxAllowedDeltaT(0.);
594 
595  const recob::Hit& pair1UHit(pair1->getHits().front()->getHit());
596  const recob::Hit& pair2UHit(pair2->getHits().front()->getHit());
597 
598  int maxDeltaWires = pair1->getHits().size() < pair2->getHits().size()
599  ? int(pair1->getHits().size())
600  : int(pair2->getHits().size());
601 
602  int pair1UWire = pair1UHit.WireID().Wire;
603  double pair1UPeakTime = pair1UHit.PeakTime();
604  double pair1ULimit = pair1UHit.RMS() * 2.;
605 
606  int pair2UWire = pair2UHit.WireID().Wire;
607  double pair2UPeakTime = pair2UHit.PeakTime();
608  double pair2ULimit = pair2UHit.RMS() * 2.;
609 
610  int deltaUWire = fabs(pair1UWire - pair2UWire);
611  double limitTotalU = pair1ULimit + pair2ULimit;
612 
613  // Do we need to constrain this?
614  limitTotalU = std::min(limitTotalU, 100.);
615 
616  // Special conditions
617  if (deltaUWire == 0) limitTotalU *= 1.5; //3.;
618  else if (deltaUWire > 1) limitTotalU *= 3.; //2.;
619 
620  // If we are happy with U wires then go to the next step
621  if (deltaUWire < maxDeltaWires && fabs(pair1UPeakTime - pair2UPeakTime) < limitTotalU) matchedHits[0] = true;
622 
623  maxDeltaT = std::max(maxDeltaT, fabs(pair1UPeakTime - pair2UPeakTime));
624  maxAllowedDeltaT = std::max(maxAllowedDeltaT, limitTotalU);
625 
626  if (deltaUWire < maxDeltaWires)
627  {
628  // Now do the W plane which will always be the "back" element
629  const recob::Hit& pair1WHit(pair1->getHits().back()->getHit());
630  const recob::Hit& pair2WHit(pair2->getHits().back()->getHit());
631 
632  int pair1WWire = pair1WHit.WireID().Wire;
633  double pair1WPeakTime = pair1WHit.PeakTime();
634  double pair1WLimit = pair1WHit.RMS() * 2.;
635 
636  int pair2WWire = pair2WHit.WireID().Wire;
637  double pair2WPeakTime = pair2WHit.PeakTime();
638  double pair2WLimit = pair2WHit.RMS() * 2.;
639 
640  int deltaWWire = fabs(pair1WWire - pair2WWire);
641  double limitTotalW = pair1WLimit + pair2WLimit;
642 
643  // Do we need to constrain this?
644  limitTotalW = std::min(limitTotalW, 100.);
645 
646  // Special conditions
647  if (deltaWWire == 0) limitTotalW *= 1.5; //3.;
648  else if (deltaWWire > 1) limitTotalW *= 3.0; //2.;
649 
650  // If we are happy with U wires then go to the next step
651  if (deltaWWire < maxDeltaWires && fabs(pair1WPeakTime - pair2WPeakTime) < limitTotalW) matchedHits[1] = true;
652 
653  maxDeltaT = std::max(maxDeltaT, fabs(pair1WPeakTime - pair2WPeakTime));
654  maxAllowedDeltaT = std::max(maxAllowedDeltaT, limitTotalW);
655 
656  if (deltaWWire < maxDeltaWires)
657  {
658  // If we are triplet then keep going
659  if (pair1->getHits().size() > 2 && pair2->getHits().size() > 2)
660  {
661  // Now do the W plane which will always be the "back" element
662  const recob::Hit& pair1VHit(pair1->getHits()[1]->getHit());
663  const recob::Hit& pair2VHit(pair2->getHits()[1]->getHit());
664 
665  int pair1VWire = pair1VHit.WireID().Wire;
666  double pair1VPeakTime = pair1VHit.PeakTime();
667  double pair1VLimit = pair1VHit.RMS() * 2.;
668 
669  int pair2VWire = pair2VHit.WireID().Wire;
670  double pair2VPeakTime = pair2VHit.PeakTime();
671  double pair2VLimit = pair2VHit.RMS() * 2.;
672 
673  int deltaVWire = fabs(pair1VWire - pair2VWire);
674  double limitTotalV = pair1VLimit + pair2VLimit;
675 
676  // Do we need to constrain this?
677  limitTotalV = std::min(limitTotalV, 100.);
678 
679  // Special conditions
680  if (deltaVWire == 0) limitTotalV *= 2.; //4.;
681  else if (deltaVWire > 1) limitTotalV *= 4.; //2.;
682 
683  // If we are happy with U wires then go to the next step
684  if (deltaVWire < maxDeltaWires)
685  {
686  if (fabs(pair1VPeakTime - pair2VPeakTime) < limitTotalV) matchedHits[2] = true;
687 
688  maxDeltaT = std::max(maxDeltaT, fabs(pair1VPeakTime - pair2VPeakTime));
689  maxAllowedDeltaT = std::max(maxAllowedDeltaT, limitTotalV);
690  }
691  else
692  {
693  // If wires are not in range then we zap the total vector
694  matchedHits[0] = false;
695  matchedHits[1] = false;
696  }
697  }
698  }
699  }
700 
701  // Now count the number of matches
702  int numMatches = std::count(matchedHits.begin(), matchedHits.end(), true);
703 
704  if (numMatches > 2) consistent = true;
705  else if (numMatches > 1)
706  {
707  maxAllowedDeltaT = std::min(5.*maxAllowedDeltaT, 50.);
708 
709  if (maxDeltaT < maxAllowedDeltaT) consistent = true;
710  }
711 
712  return consistent;
713 }
714 
715 
717  const reco::ClusterHit2D* hit2,
718  double hitWidthSclFctr,
719  size_t hitPairCntr) const
720 {
721  // No matter what happens, we will return a HitPair object
722  unsigned statusBits(0);
723  double position[] = {0.,0.,0.};
724  double totalCharge(0.);
725  double avePeakTime(0.);
726  double deltaPeakTime(0.);
727  double sigmaPeakTime(0.);
728  double overlapFraction(0.);
729  double hitDocaToAxis(9999.);
730  double hitArclenToPoca(0.);
731 
732  std::vector<const reco::ClusterHit2D*> hitVector;
733 
734  hitVector.push_back(hit1);
735  hitVector.push_back(hit2);
736 
737  // We assume in this routine that we are looking at hits in different views
738  // The first mission is to check that the wires intersect
739  const geo::WireID& hit1WireID = hit1->getHit().WireID();
740  const geo::WireID& hit2WireID = hit2->getHit().WireID();
741 
742  geo::WireIDIntersection widIntersect;
743 
744  if (m_geometry->WireIDsIntersect(hit1WireID, hit2WireID, widIntersect))
745  {
746  // Wires intersect so now we can check the timing
747  // We'll need the plane offsets to proceed, these can be inferred from the "time ticks" stored in the 2D hit
748  double hit1Offset = hit1->getHit().PeakTime() - hit1->getTimeTicks();
749  double hit2Offset = hit2->getHit().PeakTime() - hit2->getTimeTicks();
750 
751  // Basically, require that the hit times "overlap"
752  // Check here is that they are inconsistent
753  double hit1Peak = hit1->getHit().PeakTime() - hit1Offset;
754  double hit1Sigma = hit1->getHit().RMS() * 2. / 2.355;
755 
756  double hit2Peak = hit2->getHit().PeakTime() - hit2Offset;
757  double hit2Sigma = hit2->getHit().RMS() * 2. / 2.355;
758 
759 // double hit1Width = 2.*0.5*hit1->getHit().RMS() * 2.;
760 // double hit2Width = 2.*0.5*hit2->getHit().RMS() * 2.;
761  double hit1Width = hitWidthSclFctr * hit1->getHit().RMS() * 2.;
762  double hit2Width = hitWidthSclFctr * hit2->getHit().RMS() * 2.;
763 
764  // Check hit times are consistent
765  if (fabs(hit1Peak - hit2Peak) <= (hit1Width + hit2Width))
766  {
767  double maxUpper = std::min(hit1Peak+hit1Width,hit2Peak+hit2Width);
768  double minLower = std::max(hit1Peak-hit1Width,hit2Peak-hit2Width);
769  double overlap = maxUpper - minLower;
770 
771  overlapFraction = 0.5 * overlap / std::min(hit1Width,hit2Width);
772 
773  if (overlapFraction > 0.)
774  {
775  avePeakTime = 0.5 * (hit1Peak + hit2Peak);
776  deltaPeakTime = fabs(hit1Peak - hit2Peak);
777  sigmaPeakTime = sqrt(hit1Sigma*hit1Sigma + hit2Sigma*hit2Sigma);
778  totalCharge = hit1->getHit().Integral() + hit2->getHit().Integral();
779 
780  sigmaPeakTime = std::max(sigmaPeakTime, 0.1);
781 
782  double xPositionHit1(hit1->getXPosition());
783  double xPositionHit2(hit2->getXPosition());
784  double xPosition = 0.5 * (xPositionHit1 + xPositionHit2);
785 
786  position[0] = xPosition;
787  position[1] = widIntersect.y;
788  position[2] = widIntersect.z;
789 
790  // If to here then we need to sort out the hit pair code telling what views are used
791  statusBits = (unsigned(hit1->getHit().View()) + 1) * unsigned(hit2->getHit().View());
792  }
793  }
794  }
795 
796  // Create the 3D cluster hit
797  reco::ClusterHit3D hitPair(hitPairCntr,
798  statusBits,
799  position,
800  totalCharge,
801  avePeakTime,
802  deltaPeakTime,
803  sigmaPeakTime,
804  overlapFraction,
805  hitDocaToAxis,
806  hitArclenToPoca,
807  hitVector);
808  // Send it back
809  return hitPair;
810 }
811 
812 
814  const reco::ClusterHit2D* hitV) const
815 {
816  // No matter what happens, we will return a HitPair object
817  unsigned statusBits(0x8); // This one is easy... 0x8 means triplet
818  double position[] = {0.,0.,0.};
819  double totalCharge(0.);
820  double avePeakTime(0.);
821  double deltaPeakTime(1000.);
822  double sigmaPeakTime(0.);
823  double overlapFraction(0.);
824  double hitDocaToAxis(9999.);
825  double hitArclenToPoca(0.);
826 
827  std::vector<const reco::ClusterHit2D*> hitVector;
828 
829  // Recover all the hits involved
830  const reco::ClusterHit2D* hitU(pairUW.getHits()[0]);
831  const reco::ClusterHit2D* hitW(pairUW.getHits()[1]);
832 
833  // Let's do a quick consistency check on the input hits to make sure we are in range...
834  //double uvDeltaT = std::max(pairUV.getDeltaPeakTime(), 2.); // Set a lower limit
835  //double uvSigma = pairUV.getSigmaPeakTime();
836  //double wSigma = (hitW->getHit().RMS() * 2) / 2.355;
837  double uwDeltaT = 1.;
838  double uwSigma = 2.355 * pairUW.getSigmaPeakTime();
839  double vSigma = 2. * hitV->getHit().RMS();
840 
841  // Require the W hit to be "in range" with the UV Pair
842  if (fabs(hitV->getTimeTicks() - pairUW.getAvePeakTime()) < uwDeltaT * (uwSigma + vSigma))
843  {
844  // Use the existing code to see the U and W hits are willing to pair with the V hit
845  reco::ClusterHit3D pairUV = makeHitPair(hitU, hitV, 1.); //3.);
846  reco::ClusterHit3D pairVW = makeHitPair(hitW, hitV, 1.); //3.);
847 
848  // If pairs made here then we can make a triplet
849  // This condition means one of the two must be good
850  if (pairUV.getAvePeakTime() > 0. || pairVW.getAvePeakTime() > 0.)
851  {
852  double numGoodPairs(1.);
853 
854  // Position is simply the average
855  position[0] = pairUW.getPosition()[0];
856  position[1] = pairUW.getPosition()[1];
857  position[2] = pairUW.getPosition()[2];
858 
859  // Ave, delta and sigmas
860  avePeakTime = pairUW.getAvePeakTime();
861  deltaPeakTime = pairUW.getDeltaPeakTime();
862  sigmaPeakTime = pairUW.getSigmaPeakTime() * pairUW.getSigmaPeakTime();
863 
864  // Overlap fraction... hmmm....
865  overlapFraction = pairUW.getOverlapFraction();
866 
867  // Finally, total charge
868  totalCharge = pairUW.getTotalCharge();
869 
870  // Is the first W pair good?
871  if (pairUV.getAvePeakTime() > 0.)
872  {
873  position[0] += pairUV.getPosition()[0];
874  position[1] += pairUV.getPosition()[1];
875  position[2] += pairUV.getPosition()[2];
876  avePeakTime += pairUV.getAvePeakTime();
877  deltaPeakTime = std::max(deltaPeakTime, pairUV.getDeltaPeakTime());
878  sigmaPeakTime += pairUV.getSigmaPeakTime() * pairUV.getSigmaPeakTime();
879  overlapFraction = std::max(overlapFraction, pairUV.getOverlapFraction());
880  totalCharge += pairUV.getTotalCharge();
881  numGoodPairs++;
882  }
883 
884  // Is the second W pair good?
885  if (pairVW.getAvePeakTime() > 0.)
886  {
887  position[0] += pairVW.getPosition()[0];
888  position[1] += pairVW.getPosition()[1];
889  position[2] += pairVW.getPosition()[2];
890  avePeakTime += pairVW.getAvePeakTime();
891  deltaPeakTime = std::max(deltaPeakTime, pairVW.getDeltaPeakTime());
892  sigmaPeakTime += pairUV.getSigmaPeakTime() * pairVW.getSigmaPeakTime();
893  overlapFraction = std::max(overlapFraction, pairVW.getOverlapFraction());
894  totalCharge += pairVW.getTotalCharge();
895  numGoodPairs++;
896  }
897 
898  // Get averages
899  position[0] /= numGoodPairs;
900  position[1] /= numGoodPairs;
901  position[2] /= numGoodPairs;
902  avePeakTime /= numGoodPairs;
903  sigmaPeakTime = sqrt(sigmaPeakTime);
904 
905  // Remember our hits!
906  hitVector.push_back(hitU);
907  hitVector.push_back(hitV);
908  hitVector.push_back(hitW);
909  }
910  }
911 
912  // Create the 3D cluster hit
913  reco::ClusterHit3D hitTriplet(0,
914  statusBits,
915  position,
916  totalCharge,
917  avePeakTime,
918  deltaPeakTime,
919  sigmaPeakTime,
920  overlapFraction,
921  hitDocaToAxis,
922  hitArclenToPoca,
923  hitVector);
924  // Send it back
925  return hitTriplet;
926 }
927 
928 
929 const reco::ClusterHit2D* DBScanAlg_DUNE35t::FindBestMatchingHit(const Hit2DSet& hit2DSet, const reco::ClusterHit3D& pair, double pairDeltaTimeLimits) const
930 {
931  static const double minCharge(0.);
932 
933  const reco::ClusterHit2D* bestVHit(0);
934 
935  double pairAvePeakTime(pair.getAvePeakTime());
936 
937  // Idea is to loop through the input set of hits and look for the best combination
938  for (const auto& hit2D : hit2DSet)
939  {
940  if (hit2D->getHit().Integral() < minCharge) continue;
941 
942  double hitVPeakTime(hit2D->getTimeTicks());
943  double deltaPeakTime(pairAvePeakTime-hitVPeakTime);
944 
945  if (deltaPeakTime > pairDeltaTimeLimits) continue;
946 
947  // if (deltaPeakTime < -pairDeltaTimeLimits) break;
948  //REL
949  if (deltaPeakTime < -pairDeltaTimeLimits) continue;
950 
951  pairDeltaTimeLimits = fabs(deltaPeakTime);
952  bestVHit = hit2D;
953  }
954 
955  return bestVHit;
956 }
957 
958 int DBScanAlg_DUNE35t::FindNumberInRange(const Hit2DSet& hit2DSet, const reco::ClusterHit3D& pair, double range) const
959 {
960  static const double minCharge(0.);
961 
962  int numberInRange(0);
963  double pairAvePeakTime(pair.getAvePeakTime());
964 
965  // Idea is to loop through the input set of hits and look for the best combination
966  for (const auto& hit2D : hit2DSet)
967  {
968  if (hit2D->getHit().Integral() < minCharge) continue;
969 
970  double hitVPeakTime(hit2D->getTimeTicks());
971  double deltaPeakTime(pairAvePeakTime-hitVPeakTime);
972 
973  if (deltaPeakTime > range) continue;
974 
975  if (deltaPeakTime < -range) break;
976 
977  numberInRange++;
978  }
979 
980  return numberInRange;
981 }
982 
983 
985 {
986  geo::WireID wireID(0,view,0,0);
987 
988  // Embed the call to the geometry's services nearest wire id method in a try-catch block
989  try {
990  wireID = m_geometry->NearestWireID(position, view);
991  }
992  catch(std::exception& exc)
993  {
994  // This can happen, almost always because the coordinates are **just** out of range
995  mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - " << exc.what() << std::endl;
996 
997  // Assume extremum for wire number depending on z coordinate
998  if (position[2] < 0.5 * m_geometry->DetLength()) wireID.Wire = 0;
999  else wireID.Wire = m_geometry->Nwires(view) - 1;
1000  }
1001 
1002  return wireID;
1003 }
1004 
1005 
1007 {
1008  geo::WireID wireID(0,0,0,0);
1009 
1010  // Embed the call to the geometry's services nearest wire id method in a try-catch block
1011  try {
1012  wireID = m_geometry->NearestWireID(position, thePlaneID.Plane, thePlaneID.TPC, thePlaneID.Cryostat);
1013  }
1014  catch(std::exception& exc)
1015  {
1016  // This can happen, almost always because the coordinates are **just** out of range
1017  mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - " << exc.what() << std::endl;
1018  std::cout << "Exception caught finding nearest wire." << std::endl;
1019 
1020 
1021  // Assume extremum for wire number depending on z coordinate
1022  if (position[2] < 0.5 * m_geometry->DetLength()) wireID.Wire = 0;
1023  else wireID.Wire = m_geometry->Nwires(thePlaneID.Plane) - 1; //Could be a problem here.Testing for now...
1024 
1025  std::cout << "WireID: " << wireID.Wire << std::endl;
1026 
1027  }
1028 
1029  return wireID;
1030 }
1031 
1032 }
1033 
1034 
1035 
1036 
1037 
1038  //size_t DBScanAlg_DUNE35t::BuildHitPairMap(ViewToHitVectorMap& viewToHitVectorMap, ViewToWireToHitSetMap& viewToWireToHitSetMap, HitPairList& hitPairList) const
1039  //{
1040  /**
1041  * @brief Given input 2D hits, build out the lists of possible 3D hits
1042  *
1043  * The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
1044  * However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
1045  * be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
1046  * and then look for the match in the V plane. In the event we don't find the match in the V plane then we
1047  * will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.
1048  */
1049  /*
1050  // Should we set a minimum total charge for a hit?
1051  size_t totalNumHits(0);
1052  size_t hitPairCntr(0);
1053  double minCharge[] = {0., 0., 0.};
1054 
1055  //o*********************************************************************************
1056  // First job is to find all pairs of hits in different views which are "consistent"
1057  // Start loops over rows in HitVectorMap
1058  ViewToHitVectorMap::iterator mapItrU = viewToHitVectorMap.find(geo::kU);
1059 
1060  if (mapItrU != viewToHitVectorMap.end())
1061  {
1062  HitVector& hitVectorU = mapItrU->second;
1063 
1064  totalNumHits += hitVectorU.size();
1065 
1066  ViewToHitVectorMap::iterator mapItrW = viewToHitVectorMap.find(geo::kW);
1067 
1068  if (mapItrW != viewToHitVectorMap.end())
1069  {
1070  HitVector& hitVectorW = mapItrW->second;
1071 
1072  totalNumHits += hitVectorW.size();
1073 
1074  if (viewToHitVectorMap.find(geo::kV) != viewToHitVectorMap.end())
1075  totalNumHits += viewToHitVectorMap[geo::kV].size();
1076 
1077  // Take advantage that hits are sorted in "start time order"
1078  // Set the inner loop iterator before starting loop over outer hits
1079  HitVector::iterator hitVectorWStartItr = hitVectorW.begin();
1080 
1081  // Now we loop over the hits in these two layers
1082  for (HitVector::iterator hitItrU = hitVectorU.begin(); hitItrU != hitVectorU.end(); hitItrU++)
1083  {
1084  const reco::ClusterHit2D* hitPtrU = *hitItrU;
1085 
1086  if (hitPtrU->getHit().Integral() < minCharge[hitPtrU->getHit().View()]) continue;
1087 
1088  // This will be used in each loop so dereference the peak time here
1089  double hitUPeakTime = hitPtrU->getTimeTicks();
1090 
1091  // Inner loop is over hits within the time range of the outer hit
1092  for (HitVector::iterator hitItrW = hitVectorWStartItr; hitItrW != hitVectorW.end(); hitItrW++)
1093  {
1094  const reco::ClusterHit2D* hitPtrW = *hitItrW;
1095 
1096  if (hitPtrW->getHit().Integral() < minCharge[hitPtrW->getHit().View()]) continue;
1097 
1098  // Hits are sorted in "peak time order" which we can take advantage of to try to speed
1099  // the loops. Basically, we can compare the peak time for the outer loop hit against
1100  // the current hit's peak time and if outside the range of interest we can take action.
1101  // The range of interest is eyeballed but is meant to account for large pulse widths
1102  // Start by dereferencing the inner hit peak time
1103  double hitWPeakTime = hitPtrW->getTimeTicks();
1104 
1105  // If the outer loop's peak time is well past the current hit's then
1106  // we should advance the inner loop's start iterator and keep going
1107  if (hitUPeakTime > hitWPeakTime + m_timeAdvanceGap)
1108  {
1109  hitVectorWStartItr++;
1110  continue;
1111  }
1112 
1113  // If the inner loop hit start time is past the end of the outer's end time, then break out loop
1114  if (hitUPeakTime + m_timeAdvanceGap < hitWPeakTime) break;
1115 
1116  // We have a candidate hit pair combination, try to make a hit
1117  reco::ClusterHit3D pair = makeHitPair(hitPtrU, hitPtrW);
1118 
1119  // The sign of success here is that the average peak time of the combined hits is > 0
1120  // (note that when hits are combined the first window offset is accounted for)
1121  if (pair.getAvePeakTime() > 0.)
1122  {
1123  bool hitInVNotFound(true);
1124 
1125  // Recover the WireID nearest in the V plane to the position of the pair
1126  const geo::WireID wireIDV = NearestWireID(pair.getPosition(), geo::kV);
1127 
1128  // We believe the code that returns the ID is offset by one
1129  WireToHitSetMap::iterator wireToHitSetMapVItr = viewToWireToHitSetMap[geo::kV].find(wireIDV.Wire);
1130 
1131  if (wireToHitSetMapVItr != viewToWireToHitSetMap[geo::kV].end())
1132  {
1133  const reco::ClusterHit2D* hit2DV = FindBestMatchingHit(wireToHitSetMapVItr->second, pair, m_numSigmaPeakTime*pair.getSigmaPeakTime());
1134 
1135  // If a V hit found then it should be straightforward to make the triplet
1136  if (hit2DV)
1137  {
1138  reco::ClusterHit3D triplet = makeHitTriplet(pair, hit2DV);
1139 
1140  if (triplet.getAvePeakTime() > 0.)
1141  {
1142  triplet.setID(hitPairCntr++);
1143  hitPairList.emplace_back(std::unique_ptr<reco::ClusterHit3D>(new reco::ClusterHit3D(triplet)));
1144  }
1145 
1146  hitInVNotFound = false;
1147  }
1148  }
1149 
1150  // NO V hit found, let's consider allowing for the pair to be kept
1151  // WARNING: ugly code coming since I've not thought about how to do this more elegantly
1152  if (hitInVNotFound)
1153  {
1154  // Idea will be to search either side of the wire in question looking to see how big a "gap"
1155  // to any hits on adjacent wires that are "in time"
1156  int lowSideGap(0);
1157  int nLowSideHits(0);
1158  int hiSideGap(0);
1159  int nHiSideHits(0);
1160 
1161  // Idea is to check up to 3 wires either side of the target wire
1162  // Do this in two pieces since my brain is not seeing the easy way at the moment
1163  for(int loopIdx = 0; loopIdx < 3; loopIdx++)
1164  {
1165  // Set and check wire index
1166  int wireIdx = wireIDV.Wire - loopIdx - 1;
1167 
1168  // Make sure we are in range
1169  if (wireIdx < 0) break;
1170 
1171  // Find the set of hits along the current wire
1172  WireToHitSetMap::iterator wireToHitSetMapVItr = viewToWireToHitSetMap[geo::kV].find(wireIdx);
1173 
1174  // Check there is something there and if there is find the closest hit
1175  if (wireToHitSetMapVItr != viewToWireToHitSetMap[geo::kV].end())
1176  {
1177  // Set a range that increases as we move further away in wires
1178  double range = 3. * double(loopIdx + 1) * pair.getSigmaPeakTime();
1179  int numInRange = FindNumberInRange(wireToHitSetMapVItr->second, pair, range);
1180 
1181  // Increment count of wires with hits near target wire
1182 // nLowSideHits++;
1183 
1184  // If hits were found in range then we are done with this loop
1185  if (numInRange > 0)
1186  {
1187  nLowSideHits++;
1188  break;
1189  }
1190  // Otherwise it is a gap and we do the accounting
1191  else lowSideGap++;
1192  }
1193  }
1194 
1195  // Second piece tries to check the high side
1196  for(int loopIdx = 0; loopIdx < 3; loopIdx++)
1197  {
1198  int wireIdx = wireIDV.Wire + loopIdx + 1;
1199 
1200  if (wireIdx >= int(m_geometry->Nwires(geo::kV))) break;
1201 
1202  // Find the set of hits along the current wire
1203  WireToHitSetMap::iterator wireToHitSetMapVItr = viewToWireToHitSetMap[geo::kV].find(wireIdx);
1204 
1205  // Check there is something there and if there is find the closest hit
1206  if (wireToHitSetMapVItr != viewToWireToHitSetMap[geo::kV].end())
1207  {
1208  double range = 3. * double(loopIdx + 1) * pair.getSigmaPeakTime();
1209  int numInRange = FindNumberInRange(wireToHitSetMapVItr->second, pair, range);
1210 
1211  // Again, increment counter
1212 // nHiSideHits++;
1213 
1214  // If a V hit found then check how close to wire
1215  if (numInRange)
1216  {
1217  nHiSideHits++;
1218  break;
1219  }
1220  else hiSideGap++;
1221  }
1222  }
1223 
1224  // Condition to keep the hit pair instead of requiring a triplet is that there is a partner hit
1225  // on one side and a gap of no more than 2 wires on the other
1226  if ((nLowSideHits > 0 && nHiSideHits > 0) && (lowSideGap == 0 || hiSideGap == 0))
1227  {
1228  pair.setID(hitPairCntr++);
1229  hitPairList.emplace_back(std::unique_ptr<reco::ClusterHit3D>(new reco::ClusterHit3D(pair)));
1230  }
1231  }
1232  }
1233  }
1234  }
1235  }
1236  }
1237 
1238  // Return the hit pair list but sorted by z and y positions (faster traversal in next steps)
1239  hitPairList.sort(SetPositionOrder);
1240 
1241  // Where are we?
1242  mf::LogDebug("Cluster3D") << "Total number hits: " << totalNumHits << std::endl;
1243  mf::LogDebug("Cluster3D") << "Created a total of " << hitPairList.size() << " hit pairs, counted: " << hitPairCntr << std::endl;
1244 
1245  return hitPairList.size();
1246 
1247 
1248  */
1249 
1250 
1251 
1252 //size_t DBScanAlg_DUNE35t::BuildNeighborhoodMap(HitPairList& hitPairList, EpsPairNeighborhoodMapVec& epsPairNeighborhoodMapVec) const
1253 //{
1254  /**
1255  * @brief build out the epsilon neighborhood map to be used by DBScan
1256  */
1257  /*
1258  size_t consistentPairsCnt(0);
1259  size_t pairsChecked(0);
1260 
1261  //o**********************************************************************************
1262  // Given the list of pairs of hits which are consistent with each other, build out the
1263  // epsilon neighbor maps
1264  // Now we go through the pair list and basically repeat the same exercise as above
1265  // The following assumes that the HitPairList is ordered
1266  // a) in increasing Z for hits which are not on the "same W wire",
1267  // b) in increasing U (Y) for hits on the same W wire
1268 
1269  for (HitPairList::const_iterator pairItrO = hitPairList.begin(); pairItrO != hitPairList.end(); pairItrO++)
1270  {
1271  const reco::ClusterHit3D* hitPairO = (*pairItrO).get();
1272  const size_t hitPairOID = hitPairO->getID();
1273 
1274  // Need to initialize the "first" part of the vector pseudo map
1275  epsPairNeighborhoodMapVec[hitPairOID].first = hitPairO;
1276 
1277  // Get reference to the list for this hit so we don't look it up inside the loop
1278  EpsPairNeighborhoodPair& hitPairOPair(epsPairNeighborhoodMapVec[hitPairOID].second);
1279 
1280  HitPairList::const_iterator pairItrI = pairItrO;
1281 
1282  std::map<int, std::pair<double, const reco::ClusterHit3D*> > bestTripletMap;
1283 
1284  // Get the wire numbers for this triplet
1285  int pairOWireU(hitPairO->getHits().front()->getHit().WireID().Wire);
1286  int pairOWireV(0);
1287  int pairOWireW(hitPairO->getHits().back()->getHit().WireID().Wire);
1288 
1289  if (hitPairO->getHits().size() > 2)
1290  pairOWireV = hitPairO->getHits()[1]->getHit().WireID().Wire;
1291  else
1292  {
1293  const geo::WireID wireIDV = NearestWireID(hitPairO->getPosition(), geo::kV);
1294 
1295  pairOWireV = wireIDV.Wire;
1296  }
1297 
1298  // Set maximums
1299  int maxDeltaW(2);
1300  int maxDeltaU(2);
1301  int maxDeltaV(2);
1302  int maxSumAbsUV(4);
1303 
1304  while (++pairItrI != hitPairList.end())
1305  {
1306  const reco::ClusterHit3D* hitPairI = (*pairItrI).get();
1307 
1308  // Note that the hits have been sorted by z, then y and then in x
1309  // Translation: hits are sorted by W wire, then in Y (increasing u, decreasing v)
1310  // and then in time.
1311  int pairIWireU(hitPairI->getHits().front()->getHit().WireID().Wire);
1312  int pairIWireV(0);
1313  int pairIWireW(hitPairI->getHits().back()->getHit().WireID().Wire);
1314 
1315  if (hitPairI->getHits().size() > 2)
1316  pairIWireV = hitPairI->getHits()[1]->getHit().WireID().Wire;
1317  else
1318  {
1319  const geo::WireID wireIDV = NearestWireID(hitPairI->getPosition(), geo::kV);
1320 
1321  pairIWireV = wireIDV.Wire;
1322  }
1323 
1324  // Form the wire number differences so we can check ranges
1325  int deltaU(pairIWireU - pairOWireU);
1326  int deltaV(pairIWireV - pairOWireV);
1327  int deltaW(pairIWireW - pairOWireW);
1328 
1329  // If we have passed the max delta W then we are done with the loop
1330  if (deltaW > maxDeltaW) break;
1331 
1332  // Check limits on U,V differences and if past then continue to next
1333  if (fabs(deltaU) > maxDeltaU || fabs(deltaV) > maxDeltaV) continue;
1334 
1335  // Special case
1336  if (fabs(deltaU) + fabs(deltaV) > maxSumAbsUV) continue;
1337 
1338  // Keep count...
1339  pairsChecked++;
1340 
1341  // This is the tight constraint on the hits
1342  if (consistentPairs(hitPairO, hitPairI))
1343  {
1344  int bestBin = 100 * deltaW + 10 * deltaU + deltaV;
1345  double bestDist = 10000.;
1346 
1347  if (bestTripletMap.find(bestBin) != bestTripletMap.end())
1348  {
1349  bestDist = bestTripletMap[bestBin].first;
1350  }
1351 
1352  double newDist = fabs(hitPairI->getX() - hitPairO->getX());
1353 
1354  // This is an attempt to "prefer" triplets over pairs
1355  if (hitPairI->getHits().size() < 3) newDist += 25.;
1356 
1357  if (newDist < bestDist) bestTripletMap[bestBin] = std::pair<double, const reco::ClusterHit3D*>(newDist, hitPairI);
1358 
1359  // Check limits
1360  if (deltaW == 0 && deltaU == -1 && deltaV == 1) maxSumAbsUV = 2;
1361  else if (deltaW == 1 && ((deltaU == 0 && deltaV == 1) || (deltaU == 1 && deltaV == 0))) maxDeltaW = 1;
1362  }
1363  }
1364 
1365  for(const auto& bestMapItr : bestTripletMap)
1366  {
1367  const reco::ClusterHit3D* hitPairI(bestMapItr.second.second);
1368 
1369  hitPairOPair.first.incrementCount();
1370  hitPairOPair.second.emplace_back(hitPairI);
1371 
1372  epsPairNeighborhoodMapVec[hitPairI->getID()].second.first.incrementCount();
1373  epsPairNeighborhoodMapVec[hitPairI->getID()].second.second.emplace_back(hitPairO);
1374 
1375  consistentPairsCnt++;
1376  }
1377  }
1378 
1379  mf::LogDebug("Cluster3D") << "Consistent pairs: " << consistentPairsCnt << " of " << pairsChecked << " checked." << std::endl;
1380 
1381  return consistentPairsCnt;
1382 }
1383 */
1384 
1385 
1386 
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
intermediate_table::iterator iterator
std::map< int, HitPairListPtr > HitPairClusterMap
Definition: Cluster3D.h:338
double z
z position of intersection
Definition: geo_types.h:805
std::vector< EpsPairNeighborhoodMapPair > EpsPairNeighborhoodMapVec
float getTimeTicks() const
Definition: Cluster3D.h:76
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:158
Planes which measure V.
Definition: geo_types.h:130
geo::WireID WireID() const
Definition: Hit.h:233
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:220
void ClusterHitsDBScan(ViewToHitVectorMap &viewToHitVectorMap, ViewToWireToHitSetMap &viewToWiretoHitSetMap, HitPairList &hitPairList, reco::HitPairClusterMap &hitPairClusterMap)
Given a set of recob hits, run DBscan to form 3D clusters.
constexpr T pow(T x)
Definition: pow.h:72
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
std::map< geo::View_t, WireToHitSetMap > ViewToWireToHitSetMap
size_t BuildNeighborhoodMap(HitPairList &hitPairList, EpsPairNeighborhoodMapVec &epsPairNeighborhoodMapVec) const
Given an input HitPairList, build out the map of nearest neighbors.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
float getTotalCharge() const
Definition: Cluster3D.h:162
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:232
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
intermediate_table::const_iterator const_iterator
float getSigmaPeakTime() const
Definition: Cluster3D.h:165
std::pair< DBScanParams, EpsPairNeighborhoodList > EpsPairNeighborhoodPair
reco::ClusterHit3D makeHitPair(const reco::ClusterHit2D *hit1, const reco::ClusterHit2D *hit2, double hitWidthSclFctr=1., size_t hitPairCntr=0) const
Make a HitPair object by checking two hits.
float getOverlapFraction() const
Definition: Cluster3D.h:167
bool SetPeakHitPairIteratorOrder(const HitPairList::iterator &left, const HitPairList::iterator &right)
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
art framework interface to geometry description
const recob::Hit * getHit() const
Definition: Cluster3D.h:78
float getY() const
Definition: Cluster3D.h:160
bool SetPositionOrderY(const std::unique_ptr< reco::ClusterHit3D > &left, const std::unique_ptr< reco::ClusterHit3D > &right)
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
float getAvePeakTime() const
Definition: Cluster3D.h:163
const reco::ClusterHit2D * FindBestMatchingHit(const Hit2DSet &hit2DSet, const reco::ClusterHit3D &pair, double pairDeltaTimeLimits) const
A utility routine for finding a 2D hit closest in time to the given pair.
Planes which measure U.
Definition: geo_types.h:129
geo::WireID NearestWireID(const double *position, const geo::View_t &view) const
Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range...
reco::ClusterHit3D makeHitTriplet(const reco::ClusterHit3D &pair, const reco::ClusterHit2D *hit2) const
Make a HitPair object by checking two hits.
bool SetPositionOrderZ(const std::unique_ptr< reco::ClusterHit3D > &left, const std::unique_ptr< reco::ClusterHit3D > &right)
std::map< geo::View_t, HitVector > ViewToHitVectorMap
std::list< std::unique_ptr< reco::ClusterHit3D > > HitPairList
T get(std::string const &key) const
Definition: ParameterSet.h:271
size_t BuildHitPairMap(ViewToHitVectorMap &viewToHitVectorMap, ViewToWireToHitSetMap &viewToWiretoHitSetMap, HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
size_t getID() const
Definition: Cluster3D.h:156
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
float getXPosition() const
Definition: Cluster3D.h:75
This algorithm will create and then cluster 3D hits using DBScan.
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
static int max(int a, int b)
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
size_t m_minPairPts
Data members to follow.
Encapsulate the geometry of a wire.
std::set< const reco::ClusterHit2D *, Hit2DSetCompare > Hit2DSet
std::vector< reco::ClusterHit2D * > HitVector
What follows are several highly useful typedefs which we want to expose to the outside world...
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
bool operator()(const reco::HitPairClusterMap::iterator &left, const reco::HitPairClusterMap::iterator &right)
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
void setID(const size_t &id) const
Definition: Cluster3D.h:179
Encapsulate the construction of a single detector plane.
void reconfigure(fhicl::ParameterSet const &pset)
double y
y position of intersection
Definition: geo_types.h:804
bool SetPositionOrderX(const std::unique_ptr< reco::ClusterHit3D > &left, const std::unique_ptr< reco::ClusterHit3D > &right)
DBScanAlg_DUNE35t(fhicl::ParameterSet const &pset)
Constructor.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
bool SetPositionOrder(const std::unique_ptr< reco::ClusterHit3D > &left, const std::unique_ptr< reco::ClusterHit3D > &right)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::pair< const reco::ClusterHit3D *, EpsPairNeighborhoodPair > EpsPairNeighborhoodMapPair
void expandCluster(EpsPairNeighborhoodMapVec &nMap, EpsPairNeighborhoodMapVec::iterator vecItr, reco::HitPairListPtr &cluster, size_t minPts) const
the main routine for DBScan
float getDeltaPeakTime() const
Definition: Cluster3D.h:164
geo::WireID NearestWireID(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
Access the description of detector geometry.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:131
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
int FindNumberInRange(const Hit2DSet &hit2DSet, const reco::ClusterHit3D &pair, double range) const
A utility routine for returning the number of 2D hits from the list in a given range.
double accumulated_real_time() const
Definition: cpu_timer.cc:66
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
virtual ~DBScanAlg_DUNE35t()
Destructor.
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:171
geo::WireID NearestWireID_mod(const double *position, const geo::PlaneID &thePlaneID) const
void start()
Definition: cpu_timer.cc:83
float getZ() const
Definition: Cluster3D.h:161
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
Encapsulate the construction of a single detector plane.
std::list< const reco::ClusterHit3D * > EpsPairNeighborhoodList
float getX() const
Definition: Cluster3D.h:159
bool consistentPairs(const reco::ClusterHit3D *pair1, const reco::ClusterHit3D *pair2) const
The bigger question: are two pairs of hits consistent?