SNSlicer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SNSlicer
3 // Module Type: producer
4 // File: SNSlicer_module.cc
5 //
6 // Generated at Wed May 10 16:28:41 2017 by Justin Vasel using artmod
7 // from cetpkgsupport v1_10_02.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // C includes
11 #include <algorithm>
12 #include <vector>
13 
14 // Framework includes
22 #include "fhiclcpp/ParameterSet.h"
24 
25 // ROOT includes
26 #include "TVector3.h"
27 
28 // NOvASoft includes
29 #include "Geometry/Geometry.h"
32 #include "RecoBase/Cluster.h"
33 #include "RecoBaseHit/CellHit.h"
34 
35 namespace sn {
36  class SNSlicer;
37 }
38 
39 class sn::SNSlicer : public art::EDProducer {
40 public:
41  explicit SNSlicer(fhicl::ParameterSet const & p);
42 
43  SNSlicer(SNSlicer const &) = delete;
44  SNSlicer(SNSlicer &&) = delete;
45  SNSlicer & operator = (SNSlicer const &) = delete;
46  SNSlicer & operator = (SNSlicer &&) = delete;
47 
49  bool ClusterIsGood(rb::Cluster* cluster);
50 
51  // Required functions.
52  void produce(art::Event & e) override;
53 
54 
55 private:
57 
58  std::string fCellHitInput; /// CellHit product label
59  unsigned int fADCThreshold; /// ADC threshold
60  unsigned int fPromptTimeThreshold; /// Min time separation
61  unsigned int fDelayedTimeThreshold; /// Min time separation
62  unsigned int fCellThreshold; /// Min cell separation
63  unsigned int fPlaneThreshold; /// Min plane separation
64  unsigned int fMinNumberHits; /// Min number of hits for cluster to be acceptable
65  unsigned int fMaxNumberHits; /// Max number of hits for cluster to be acceptable
66  bool fGroupPromptClusters; /// Attempt to conncet disparate prompt clusters?
67 };
68 
69 
71 fCellHitInput(p.get<std::string>("CellHitInput")),
72 fADCThreshold(p.get<int>("ADCThreshold")),
73 fPromptTimeThreshold(p.get<int>("PromptTimeSeparationThreshold")),
74 fDelayedTimeThreshold(p.get<int>("DelayedTimeSeparationThreshold")),
75 fCellThreshold(p.get<int>("CellSeparationThreshold")),
76 fPlaneThreshold(p.get<int>("PlaneSeparationThreshold")),
77 fMinNumberHits(p.get<unsigned int>("MinNumberHits")),
78 fMaxNumberHits(p.get<unsigned int>("MaxNumberHits")),
79 fGroupPromptClusters(p.get<bool>("GroupPromptClusters"))
80 {
81  produces<std::vector<rb::Cluster>>();
82 }
83 
85 {
86  bool kNumberHitsOutOfRange = (cluster->NCell() > fMaxNumberHits);
87  bool kMaxPlaneExtentExceeded = (cluster->ExtentPlane() > fPlaneThreshold);
88  bool kMaxXCellExtentExceeded = (cluster->NXCell()>0) ? (cluster->ExtentCell(geo::kX) > fCellThreshold) : false;
89  bool kMaxYCellExtentExceeded = (cluster->NYCell()>0) ? (cluster->ExtentCell(geo::kY) > fCellThreshold) : false;
90  bool kMaxTnsExtentExceeded = (cluster->ExtentTNS() > fPromptTimeThreshold);
91  bool kMaxTotalAdcExceeded = (cluster->TotalADC() > fADCThreshold);
92  // bool kClusterIs2D = (cluster->NXCell()==0 || cluster->NYCell()==0);
93 
94  bool kClusterIsBad = kNumberHitsOutOfRange || kMaxPlaneExtentExceeded ||
95  kMaxXCellExtentExceeded || kMaxYCellExtentExceeded ||
96  kMaxTnsExtentExceeded || kMaxTotalAdcExceeded;
97 
98  return (!kClusterIsBad);
99 }
100 
102 {
103  // Determine the view
104  unsigned int nCellX = cluster.NXCell();
105  unsigned int nCellY = cluster.NYCell();
106 
107  rb::Cluster cFinal;
108 
109  if (nCellX>0 && nCellY==0) {
110  cFinal = rb::Cluster(geo::kX, cluster.AllCells(), cluster.ID());
111  } else if (nCellX==0 && nCellY>0) {
112  cFinal = rb::Cluster(geo::kY, cluster.AllCells(), cluster.ID());
113  } else {
114  cFinal = rb::Cluster(cluster.AllCells(), cluster.ID());
115  }
116 
117  return cFinal;
118 }
119 
120 
122 {
123  // Clusters will be stored here (already-used hits are tracked as well)
124  std::vector<rb::Cluster> promptClusters;
125  std::unique_ptr<std::vector<rb::Cluster>> clusters(new std::vector<rb::Cluster>);
126 
127  // Get the cell hits
129  e.getByLabel(fCellHitInput, hitsHandle);
130 
131  // Put the hits into a vector so we can alter them
132  std::vector<art::Ptr<rb::CellHit>> hits;
133  hits.reserve(hitsHandle->size());
134  for (unsigned int idx = 0; idx < hitsHandle->size(); ++idx){
135  art::Ptr<rb::CellHit> hit(hitsHandle, idx);
136  hits.push_back(hit);
137  }
138 
139  // Sort hits by time.
140  // NOTE: The rest of the code depends on this.
141  rb::SortByTime(hits);
142 
143  bool hitConsumed[hits.size() + 1] = { false };
144 
145  /* OUTER LOOP (HITS) */
146  size_t idxCluster = 0;
147  for (unsigned int i=0; i<hits.size(); ++i) {
148  if (hitConsumed[i]) continue;
149  ++idxCluster;
150 
151  // We add this hit to the cluster, but only if it has not yet been added
152  // to any other clusters.
153  rb::Cluster* aCluster = new rb::Cluster();
154  aCluster->SetID(idxCluster);
155 
156  if (!hitConsumed[i]) {
157  aCluster->Add(hits[i]);
158  hitConsumed[i] = true;
159  }
160  else {
161  continue;
162  }
163 
164  /*
165  The inner loop continues on until a hit fails a time condition, but we want
166  to come back to the hits we missed due to a failed plane or cell condition.
167  The "latch" variable is used to perform this book-keeping.
168 
169  At the first instance a plane or cell condition fails, the latch is set to
170  the index of that hit iff the latch is currently set to zero. The latch
171  therefore does not get re-assigned on subsequent plane or cell failures,
172  only the first.
173 
174  Once the time condition fails, the hit index i is set based on the latch
175  point if it is non-zero, otherwise it is set based on j.
176 
177  The latch get reset during each iteration of the outer loop.
178 
179  Here are some simple examples. Each box represents a hit that we are looping
180  over.
181  [*] = hit passed all conditions and is added to the cluster.
182  [x] = hit failed time condition and the inner loop will terminate here.
183  [ ] = hit failed plane or cell condition
184  [L] = latch point set here
185  ^ = carat indicates where the loop will start next time
186 
187  SCENARIO 1: Several hits pass all conditions, but some intermediate hits
188  fail cell or plane condition. Time condition fails at hit 6.
189  0 1 2 3 4 5 6
190  [*] [*] [L] [ ] [*] [*] [x]
191  ^ ---> latch point at hit 2. Next outer loop will start here.
192 
193  SCENARIO 2: All hits 0-5 pass all conditions. Hit 6 fails time condition.
194  No latch point is set, so the next iteration of the outer loop will start
195  at the failure point (Hit 6).
196  0 1 2 3 4 5 6
197  [*] [*] [*] [*] [*] [*] [x]
198  ^ ---> Next outer loop will start here.
199  */
200  unsigned int latch = 0;
201 
202  /* INNER LOOP (HITS) */
203  // std::cout << std::endl;
204  for (unsigned int j=i+1; j<hits.size(); ++j) {
205  // std::cout << "Checking hit (j) " << j+1 << "/" << hits.size();
206  bool kHitsFailTimeThreshold = hits[j]->TNS() - hits[i]->TNS() > fPromptTimeThreshold;
207  bool kHitsFailPlaneThreshold = std::abs(hits[j]->Plane() - hits[i]->Plane()) > fPlaneThreshold;
208  bool kHitsFailCellThreshold = std::abs(hits[j]->Cell() - hits[i]->Cell()) > fCellThreshold;
209  bool kHitsInSameView = hits[j]->View() == hits[i]->View();
210 
211  // This hit is too far away in time from the first, and all subsequent
212  // hits willalso be so. We terminate the inner loop and set the outer loop
213  // index depending on either the latch point or the current inner loop
214  // index.
215  if (kHitsFailTimeThreshold) {
216  // std::cout << " -- SEED FOR NEXT CLUSTER" << std::endl;
217  // std::cout << "--Break--" << std::endl;
218  i = (latch==0) ? j-1 : latch-1;
219  break;
220  }
221 
222  // This hit is too far away in Planes from the first. If the latch point
223  // has not yet been set, we set it here. We move to the next iteration
224  // of the inner loop.
225  if (kHitsFailPlaneThreshold) {
226  // std::cout << " -- FAIL PLANE CONDITION." << std::endl;
227  if (latch==0) latch = j;
228  continue;
229  }
230 
231  // This hit is too far away in Cells from the first. This can only happen
232  // if both hits are in the same view. If the latch point has not yet been
233  // set, we set it here. We move to the next iteration of the inner loop.
234  if (kHitsInSameView && kHitsFailCellThreshold) {
235  // std::cout << " -- FAIL CELL CONDITION." << std::endl;
236  if (latch==0) latch = j;
237  continue;
238  }
239 
240  // If we've made it this far, this hit belongs in a cluster, but only if
241  // it hasn't been added to any other clusters yet.
242  if (!hitConsumed[j]) {
243  aCluster->Add(hits[j]);
244  if (this->ClusterIsGood(aCluster)) {
245  hitConsumed[j] = true;
246  } else {
247  aCluster->RemoveHit(hits[j]);
248  }
249  // std::cout << " -- ADDED TO CLUSTER #" << idxCluster << std::endl;
250  }
251  }
252 
253  if (this->ClusterIsGood(aCluster) && aCluster->NXCell()>0 && aCluster->NYCell()>0 && aCluster->NCell() >= fMinNumberHits) {
254  rb::Cluster fCluster = this->FinalizeCluster(*aCluster);
255 
256  promptClusters.push_back(fCluster);
257  clusters->push_back(fCluster);
258  } else {
259  hitConsumed[i] = false;
260  aCluster->RemoveHit(hits[i]);
261  }
262  delete aCluster;
263  // std::cout << std::endl;
264  }
265 
266 
267  /* Connect disparate prompt clusters */
268  // NOTE: This section never worked completely and hasn't been touched in ages.
269  // do not use without testing it first.
270  if (fGroupPromptClusters) {
271  bool clusterConsumed[promptClusters.size() + 1] = { false };
272 
273  /* OUTER LOOP (CLUSTERS) */
274  idxCluster = 0;
275  for (size_t i=0; i<promptClusters.size(); ++i) {
276  ++idxCluster;
277  std::cout << "Checking cluster (i) " << i+1 << "/" << promptClusters.size();
278  rb::Cluster* aCluster = new rb::Cluster();
279 
280  if (!clusterConsumed[i]) {
281  aCluster->Add(promptClusters[i].AllCells());
282  clusterConsumed[i] = true;
283  std::cout << " -- ADDED TO CLUSTER #" << idxCluster;
284  }
285 
286  unsigned int latch = 0;
287 
288  /* INNER LOOP (CLUSTERS) */
289  std::cout << std::endl;
290  for (size_t j=i+1; j<promptClusters.size(); ++j) {
291  std::cout << "Checking cluster (j) " << j+1 << "/" << promptClusters.size();
292 
293  // Check time
294  bool kClustersFailTimeThreshold = promptClusters[j].MeanTNS() - promptClusters[i].MeanTNS() > fDelayedTimeThreshold;
295 
296  // Check planes
297  unsigned int clusterPlanes[] = {promptClusters[i].MinPlane(), promptClusters[i].MaxPlane(), promptClusters[j].MinPlane(), promptClusters[j].MaxPlane()};
298  unsigned int* minPlane = std::min_element(std::begin(clusterPlanes), std::end(clusterPlanes));
299  unsigned int* maxPlane = std::max_element(std::begin(clusterPlanes), std::end(clusterPlanes));
300  bool kClustersFailPlaneThreshold = (*maxPlane)-(*minPlane) > fPlaneThreshold;
301 
302  // Check Cells
303  unsigned int iClusterNX = promptClusters[i].NXCell();
304  unsigned int iClusterNY = promptClusters[i].NYCell();
305  unsigned int jClusterNX = promptClusters[j].NXCell();
306  unsigned int jClusterNY = promptClusters[j].NYCell();
307 
308  geo::View_t iClusterView;
309  geo::View_t jClusterView;
310 
311  if (iClusterNX>0 && iClusterNY==0) iClusterView = geo::kX;
312  else if (iClusterNY>0 && iClusterNX==0) iClusterView = geo::kY;
313  else iClusterView = geo::kXorY;
314 
315  if (jClusterNX>0 && jClusterNY==0) jClusterView = geo::kX;
316  else if (jClusterNY>0 && jClusterNX==0) jClusterView = geo::kY;
317  else jClusterView = geo::kXorY;
318 
319  bool kBothClustersSpanBothViews = ((iClusterView == geo::kXorY) && (jClusterView == geo::kXorY));
320  bool kOneClusterSpanBothViews = ((iClusterView != geo::kXorY) != (jClusterView != geo::kXorY));
321  bool kBothClustersSpanOneView = ((iClusterView != geo::kXorY) && (jClusterView != geo::kXorY));
322  bool kClustersFailCellThreshold = false;
323 
324  // Case 1: Both clusters span both views
325  if (kBothClustersSpanBothViews) {
326  // Check X-view
327  unsigned int clusterXCells[] = {promptClusters[i].MinCell(geo::kX), promptClusters[i].MaxCell(geo::kX), promptClusters[j].MinCell(geo::kX), promptClusters[j].MaxCell(geo::kX)};
328  unsigned int* minCell = std::min_element(std::begin(clusterXCells), std::end(clusterXCells));
329  unsigned int* maxCell = std::max_element(std::begin(clusterXCells), std::end(clusterXCells));
330  bool kClustersFailXCellThreshold = (*maxCell)-(*minCell) > fCellThreshold;
331  // Check Y-view
332  unsigned int clusterYCells[] = {promptClusters[i].MinCell(geo::kY), promptClusters[i].MaxCell(geo::kY), promptClusters[j].MinCell(geo::kY), promptClusters[j].MaxCell(geo::kY)};
333  minCell = std::min_element(std::begin(clusterYCells), std::end(clusterYCells));
334  maxCell = std::max_element(std::begin(clusterYCells), std::end(clusterYCells));
335  bool kClustersFailYCellThreshold = (*maxCell)-(*minCell) > fCellThreshold;
336  kClustersFailCellThreshold = kClustersFailXCellThreshold && kClustersFailYCellThreshold;
337  }
338  // Case 2: Only one cluster spans both views
339  else if (kOneClusterSpanBothViews) {
340  geo::View_t commonView = (iClusterView == geo::kXorY) ? jClusterView : iClusterView;
341  unsigned int clusterCells[] = {promptClusters[i].MinCell(commonView), promptClusters[i].MaxCell(commonView), promptClusters[j].MinCell(commonView), promptClusters[j].MaxCell(commonView)};
342  unsigned int* minCell = std::min_element(std::begin(clusterCells), std::end(clusterCells));
343  unsigned int* maxCell = std::max_element(std::begin(clusterCells), std::end(clusterCells));
344  kClustersFailCellThreshold = (*maxCell)-(*minCell) > fCellThreshold;
345  }
346  // Case 3: Both clusters only span one view
347  else if (kBothClustersSpanOneView) {
348  // Check their cell distance if they're in the same view
349  if (iClusterView == jClusterView) {
350  unsigned int clusterCells[] = {promptClusters[i].MinCell(iClusterView), promptClusters[i].MaxCell(iClusterView), promptClusters[j].MinCell(jClusterView), promptClusters[j].MaxCell(jClusterView)};
351  unsigned int* minCell = std::min_element(std::begin(clusterCells), std::end(clusterCells));
352  unsigned int* maxCell = std::max_element(std::begin(clusterCells), std::end(clusterCells));
353  kClustersFailCellThreshold = (*maxCell)-(*minCell) > fCellThreshold;
354  }
355  }
356  // Default: We should never get here.
357  else {
358  std::cout << "I am confused about the views for clusters (i,j) = (" << iClusterView << "," << jClusterView << ")." << std::endl;
359  }
360 
361  // This cluster is too far away in time from the first, and all subsequent
362  // clusters willalso be so. We terminate the inner loop and set the outer loop
363  // index depending on either the latch point or the current inner loop
364  // index.
365  if (kClustersFailTimeThreshold) {
366  std::cout << " -- SEED FOR NEXT CLUSTER" << std::endl;
367  std::cout << "--Break--" << std::endl;
368  i = (latch==0) ? j-1 : latch-1;
369  break;
370  }
371 
372  // This cluster is too far away in Planes from the first. If the latch point
373  // has not yet been set, we set it here. We move to the next iteration
374  // of the inner loop.
375  if (kClustersFailPlaneThreshold) {
376  std::cout << " -- FAIL PLANE CONDITION." << std::endl;
377  if (latch==0) latch = j;
378  continue;
379  }
380 
381  // This cluster is too far away in Cells from the first. This can only happen
382  // if both hits are in the same view. If the latch point has not yet been
383  // set, we set it here. We move to the next iteration of the inner loop.
384  if (kClustersFailCellThreshold) {
385  std::string cellFailure;
386  if (kBothClustersSpanBothViews) cellFailure = "Both Span Both Views";
387  else if (kOneClusterSpanBothViews) cellFailure = "One Spans Both Views";
388  else if (kBothClustersSpanOneView) cellFailure = "Both Span One View";
389  std::cout << " -- FAIL CELL CONDITION (" << cellFailure << ")" << std::endl;
390  if (latch==0) latch = j;
391  continue;
392  }
393 
394  // If we've made it this far, this hit belongs in a cluster, but only if
395  // it hasn't been added to any other clusters yet.
396  if (!clusterConsumed[j]) {
397  aCluster->Add(promptClusters[j].AllCells());
398  clusterConsumed[j] = true;
399  std::cout << " -- ADDED TO CLUSTER #" << idxCluster << std::endl;
400  }
401  else std::cout << std::endl;
402  }
403  if (ClusterIsGood(aCluster) && aCluster->NXCell() > 0 && aCluster->NYCell() > 0) clusters->push_back(*aCluster);
404  std::cout << std::endl;
405  }
406  }
407 
408  // Put the clusters into the event
409  e.put(std::move(clusters));
410 }
411 
unsigned int fMaxNumberHits
Min number of hits for cluster to be acceptable.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
rb::Cluster FinalizeCluster(rb::Cluster cluster)
unsigned int fDelayedTimeThreshold
Min time separation.
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
std::string fCellHitInput
void RemoveHit(const art::Ptr< rb::CellHit > hit)
Remove hit from current cluster.
Definition: Cluster.cxx:290
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
DEFINE_ART_MODULE(TestTMapFile)
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
unsigned int ExtentCell(geo::View_t view) const
Definition: Cluster.cxx:570
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double TotalADC() const
Sum of the ADC of all the contained hits.
Definition: Cluster.cxx:360
void hits()
Definition: readHits.C:15
void produce(art::Event &e) override
unsigned int fMinNumberHits
Min plane separation.
unsigned int fPlaneThreshold
Min cell separation.
art::ServiceHandle< geo::Geometry > fGeom
SNSlicer(fhicl::ParameterSet const &p)
bool ClusterIsGood(rb::Cluster *cluster)
Remove hits from hot and cold channels.
const double j
Definition: BetheBloch.cxx:29
void SetID(int id)
Definition: Cluster.h:74
SNSlicer & operator=(SNSlicer const &)=delete
unsigned int fADCThreshold
CellHit product label.
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
unsigned int fCellThreshold
Min time separation.
const int ID() const
Definition: Cluster.h:75
Float_t e
Definition: plot.C:35
bool fGroupPromptClusters
Max number of hits for cluster to be acceptable.
double ExtentTNS() const
Definition: Cluster.h:252
Encapsulate the geometry of one entire detector (near, far, ndos)
unsigned int fPromptTimeThreshold
ADC threshold.