NERDProng_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////
2 /// \file NERDProng_module.cc
3 /// \brief Produces 2D prongs from segmentation.
4 /// NERD - Neutrino Evt Reco with Deep learning ;)
5 /// \author grohmc@fnal.gov, psihas@fnal.gov
6 //////////////////////////////////////////////////////////////////
7 
10 #include "Geometry/Geometry.h"
13 #include "RecoBase/CellHit.h"
14 #include "RecoBase/HitMap.h"
15 #include "RecoBase/Cluster.h"
16 #include "RecoBase/Prong.h"
17 #include "RecoBase/Vertex.h"
18 #include "RecoBase/FilterList.h"
19 #include "RecoBase/PID.h"
20 #include "Utilities/func/MathUtil.h"
21 #include "Utilities/AssociationUtil.h"
22 
28 
29 #include "CVN/func/Result.h"
30 
31 #include "TVector3.h"
32 #include "TTree.h"
33 #include "TH1D.h"
34 #include <cmath>
35 
39 
40 namespace nerd
41 {
43  {
44  template<class T> using Atom = fhicl::Atom<T>;
45  template<class T> using Table = fhicl::Table<T>;
47  using Name = fhicl::Name;
48 
49  // Parameters controling the splitting of clusters
53 
54  Atom <bool> ObeyPreselection{ Name("ObeyPreselection"),
55  Comment("Only run on presel slices") };
57  Comment("Specific presel used") };
58 
60  Comment("Instance label 2D") };
62  Comment("Instance label 2D") };
63 
65  Comment("Instance label 3D") };
67  Comment("Instance label 3D") };
68 
70  Comment("Instance label orphaned prong") };
72  Comment("Instance label orphaned prong") };
73 
75  Comment("Instance Label for X-View Scores") };
77  Comment("Instance Label for X-View Scores") };
79  Comment("Instance Label for Y-View Scores") };
81  Comment("Instance Label for Y-View Scores") };
83  Comment("Instance Label for 2D Scores") };
85  Comment("Instance Label for 2D Scores") };
86 
88  Comment("Minimum # of hits") };
90  Comment("Minimum # of hits X") };
92  Comment("Minimum # of hits Y") };
94  Comment("CVN Score Label") };
95  Atom <float> CVNCut{ Name("CVNCut"),
96  Comment("CVN cut for evaluation") };
97 
98  Atom <float> DistanceThreshold{ Name("DistanceThreshold"),
99  Comment("Maximum distance in cm to include unclustered hits") };
100 
101 
104  };
105 
106 
107  class NERDProng : public art::EDProducer
108  {
109  public:
110  // Allows 'nova --print-description' to work
112 
113  explicit NERDProng(const Parameters& params);
114  virtual ~NERDProng();
115 
116  void beginJob();
117  void produce(art::Event &evt);
118 
119  private:
121 
122  // An object to run the segmentation graph
124 
125  /// An object of ViewMatchAlg to do the 3D view matching
127 
129 
130  /// \brief Builds prongs from tensorflow cluster + PID results
131  /// \param prongs A pointer vector for the resulting 2d prongs
132  /// \param slice The input rb::Cluster object to use in clustering
133  /// \param vert The vertex for associations
134  void MakeProngs( std::vector<rb::Prong>* prongs,
135  std::vector<std::vector<rb::PID>>* pids,
136  art::Ptr<rb::Cluster> slice,
137  cvn::PixelMap pixelmap,
138  TVector3 vert);
139 
140  /// \brief Builds clusters from a pmResult
141  std::vector<rb::Cluster> MakeCluster(art::Ptr<rb::Cluster> slice,
142  nerd::PixelPIDMaps pmResult,
143  cvn::PixelMap pixelmap,
144  geo::View_t view);
145 
146  /// \brief Takes a 2D cluster and turns into a prong for view matching
147  void AddProngs(rb::Cluster cluster,
148  TVector3 vert,
149  std::vector<rb::Prong>* prongs);
150 
151  /// \brief Compute the distance between two hits
152  double HitToHitDistance(art::Ptr<rb::CellHit> hit1, art::Ptr<rb::CellHit> hit2);
153 
154  // \brief Min dist between a hit and a physics slice
155  double MinHitClusterDist(art::Ptr<rb::CellHit> hit, rb::Cluster clus);
156 
157  // \brief Find the nearest cluster to a hit
158  std::pair<int, double> FindNearestCluster(art::Ptr<rb::CellHit> hit, std::vector<rb::Cluster> clus);
159 
160  /// \brief Function to find hits that are not contained in any prong
161  art::PtrVector<rb::CellHit> FindOrphanedHits(art::PtrVector<rb::CellHit> c, std::vector<rb::Prong> prongs);
162 
163  /// \brief Check if a slice passes the cvncut
164  bool PassCVNCut(const cvn::Result& result, float threshold);
165 
166  };
167 }
168 ////////////////////////////////////////////////////////////////////////////////
169 namespace nerd
170 {
171 
173  : fSliceToken( consumes<std::vector<rb::Cluster>>(params().SliceLabel()) ),
174  fNERDEval(new NERDEval(params().NERDEvalPSet())),
175  fViewMatchAlg(new ViewMatchAlg(params().ViewMatch())),
176  fParams(params())
177  {
178  this->produces< std::vector<rb::Prong> >(fParams.PngInst2D());
179  this->produces< std::vector<rb::Prong> >(fParams.PngInst3D());
180  this->produces< std::vector<rb::Prong> >(fParams.PngInstOrph());
181 
182  this->produces< art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssn2D());
183  this->produces< art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssn3D());
184  this->produces< art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssnOrph());
185 
186  this->produces< art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssn2D());
187  this->produces< art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssn3D());
188  this->produces< art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssnOrph());
189 
190  this->produces< std::vector<rb::PID> >(fParams.ScoreInst2D());
191  this->produces< std::vector<rb::PID> >(fParams.ScoreInst3DX());
192  this->produces< std::vector<rb::PID> >(fParams.ScoreInst3DY());
193 
194  this->produces< art::Assns<rb::PID,rb::Prong> >(fParams.ScoreAssn2D());
195  this->produces< art::Assns<rb::PID,rb::Prong> >(fParams.ScoreAssn3DX());
196  this->produces< art::Assns<rb::PID,rb::Prong> >(fParams.ScoreAssn3DY());
197  }
198 
199  //......................................................................
200 
202  {
203  }
204 
205  //......................................................................
206 
208  {
209  if(fNERDEval) delete fNERDEval;
210  if(fViewMatchAlg) delete fViewMatchAlg;
211  }
212 
213  //......................................................................
214 
216  {
218 
219  // New products
220  // 2D Clusters
221  std::unique_ptr< std::vector<rb::Prong> > png2D(new std::vector<rb::Prong>);
222 
223  std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> >
224  png2DSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
225  std::unique_ptr< art::Assns<rb::Prong, rb::Vertex> >
226  png2DAssns(new art::Assns<rb::Prong, rb::Vertex>);
227 
228  std::unique_ptr< std::vector<rb::PID> >
229  PIDCol2(new std::vector<rb::PID>);
230 
231  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
232  PID2Assns(new art::Assns<rb::PID, rb::Prong>);
233 
234  // 3D Clusters
235  std::unique_ptr< std::vector<rb::Prong> > png3D(new std::vector<rb::Prong>);
236 
237  std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> >
238  png3DSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
239  std::unique_ptr< art::Assns<rb::Prong, rb::Vertex> >
240  png3DAssns(new art::Assns<rb::Prong, rb::Vertex>);
241 
242  std::unique_ptr< std::vector<rb::PID> >
243  PIDColX(new std::vector<rb::PID>);
244  std::unique_ptr< std::vector<rb::PID> >
245  PIDColY(new std::vector<rb::PID>);
246 
247  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
248  PIDXAssns(new art::Assns<rb::PID, rb::Prong>);
249  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
250  PIDYAssns(new art::Assns<rb::PID, rb::Prong>);
251 
252  // Orphaned Energy
253  std::unique_ptr< std::vector<rb::Prong> > pngOrph(new std::vector<rb::Prong>);
254  std::unique_ptr<art::Assns<rb::Prong, rb::Vertex> > pngOrphAssns(new art::Assns<rb::Prong, rb::Vertex>);
255  std::unique_ptr<art::Assns<rb::Prong, rb::Cluster> > pngOrphSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
256 
257  // Get from evt
259  evt.getByToken(fSliceToken, sliceHandle);
260 
261  art::FindManyP<cvn::PixelMap> fmPixMap( sliceHandle,
262  evt,
264 
265  art::FindManyP<rb::Vertex> fmVtx( sliceHandle,
266  evt,
267  fParams.VertexLabel());
268 
269  art::FindManyP<cvn::Result> fmCVN( sliceHandle,
270  evt,
271  fParams.CVNLabel());
272 
273  // Local
274  std::vector<rb::Prong> prongs;
275  std::vector<std::vector<rb::PID>> pids;
276 
277  // Loop over slices
278  for ( unsigned int iSlice=0; iSlice<sliceHandle->size(); ++iSlice ) {
279  std::cout<<"Looking at slice "<<iSlice<<std::endl;
280 
281  // Preselection
282  if( fParams.ObeyPreselection() &&
283  rb::IsFiltered(evt,sliceHandle,iSlice,fParams.PreselLabels()) )
284  continue;
285 
286  // CVN selection
287  std::vector<art::Ptr<cvn::Result>> cvnresult = fmCVN.at(iSlice);
288  if(!cvnresult.empty()){
289  if (!PassCVNCut(*cvnresult[0], fParams.CVNCut()))
290  continue;
291  }
292 
293  art::Ptr<rb::Cluster> sliceptr(sliceHandle, iSlice);
294 
295  if ( sliceptr->IsNoise() ) continue;
296  if ( sliceptr->NXCell() < fParams.MinHitX() ) continue;
297  if ( sliceptr->NYCell() < fParams.MinHitY() ) continue;
298  if ( sliceptr->NCell() < fParams.MinHit() ) continue;
299 
300  // Find any verticies that go with this slice
301  std::vector<art::Ptr<rb::Vertex> > vtx = fmVtx.at(iSlice);
302  if( vtx.size() != 1 ) continue;
303 
304  TVector3 vert = vtx[0]->GetXYZ();
305  // Same vertex check as FuzzyK
306  if( std::isnan(vert.X()) ||
307  std::isinf(vert.X()) || std::abs(vert.X())>10000) continue;
308  if( std::isnan(vert.Y()) ||
309  std::isinf(vert.Y()) || std::abs(vert.Y())>10000) continue;
310  if( std::isnan(vert.Z()) ||
311  std::isinf(vert.Z()) || std::abs(vert.Z())>10000) continue;
312 
313  // Find the PixelMap that goes with this slice
314  const std::vector<art::Ptr<cvn::PixelMap> > pixMaps = fmPixMap.at(iSlice);
315  if( pixMaps.size() != 1 ) continue;
316 
317  art::Ptr<cvn::PixelMap> pixmap = pixMaps[0];
318  // make the 2D prong candidates for this slice
319  prongs.clear();
320  pids.clear();
321 
322  this->MakeProngs(&prongs, &pids, sliceptr, *pixmap, vert);
323 
324  // perform the view matching with this set of 2D prongs
325  std::cout<<"Beginning Matching"<<std::endl;
326  fViewMatchAlg->LoadProngs(prongs, pids);
328 
329  std::cout<<"2D Prongs: "<<prongs.size()<<std::endl;
330  std::cout<<"3D Prongs: "<<fViewMatchAlg->fMatchedProngs.size()<<std::endl;
331  std::cout<<"Unmatched Prongs: "<<fViewMatchAlg->fUnmatchedProngs.size()<<std::endl;
332 
333  std::cout<<"Preparing 3D prongs"<<std::endl;
334 
335  // loop over 3D matched prongs to push into the event
336  for( unsigned int iPng = 0; iPng < fViewMatchAlg->fMatchedProngs.size(); ++iPng ){
338 
339  //move the prong start point from the elastic arms vertex to the closest hit
340  TVector3 start;
341  double minZ = png.MinZ();
342  double maxZ = png.MaxZ();
343  if ( (png.Dir().Z() > 0.0) && (minZ < png.Start().Z()) )
344  minZ = png.Start().Z();
345  if ( (png.Dir().Z() < 0.0) && (maxZ > png.Start().Z()) )
346  maxZ = png.Start().Z();
347 
348  if ( png.Dir().Z() < 0.0 ){
349  double temp = maxZ;
350  maxZ = minZ;
351  minZ = temp;
352  }
353 
354  start[0] = (png.Dir().X()/png.Dir().Z())*(minZ-png.Start().Z()) + png.Start().X();
355  start[1] = (png.Dir().Y()/png.Dir().Z())*(minZ-png.Start().Z()) + png.Start().Y();
356  start[2] = minZ;
357 
358  png.SetStart(start);
359 
360  png3D->push_back(png);
361 
362  util::CreateAssn(*this, evt, *png3D, vtx[0], *png3DAssns, UINT_MAX, fParams.PngInst3D());
363  util::CreateAssn(*this, evt, *png3D, sliceptr, *png3DSlcAssns, UINT_MAX, fParams.PngInst3D());
364 
365  // Prepare for associating the pid scores
366  std::vector<rb::PID> pidX(fViewMatchAlg->fMatchedXScores[iPng]);
367  std::vector<rb::PID> pidY(fViewMatchAlg->fMatchedYScores[iPng]);
368 
369  for(unsigned int i = 0; i < pidX.size(); ++i){
370  PIDColX->push_back(pidX[i]);
371  PIDColY->push_back(pidY[i]);
372 
373  util::CreateAssn(*this, evt, *PIDXAssns, *PIDColX, *png3D, UINT_MAX, UINT_MAX,
375  util::CreateAssn(*this, evt, *PIDYAssns, *PIDColY, *png3D, UINT_MAX, UINT_MAX,
377  }
378 
379  } // 3d prongs
380 
381  std::cout<<"Preparing 2D prongs"<<std::endl;
382 
383  //Loop over the 2d prongs for vtx related calculations
384  for( unsigned int iPng=0; iPng < fViewMatchAlg->fUnmatchedProngs.size(); ++iPng ){
386 
387  double xyStart = 0.0;
388  double zStart = 0.0;
389  double minZ = png.MinZ();
390  double maxZ = png.MaxZ();
391 
392  if ( (png.Dir().Z() > 0.0) && (minZ < png.Start().Z()) )
393  minZ = png.Start().Z();
394  if ( (png.Dir().Z() < 0.0) && (maxZ > png.Start().Z()) )
395  maxZ = png.Start().Z();
396 
397  if ( png.Dir().Z() < 0.0 ){
398  double temp = maxZ;
399  maxZ = minZ;
400  minZ = temp;
401  }
402 
403  TVector3 start;
404 
405  if( png.View() == geo::kX ){
406 
407  xyStart = (png.Dir().X()/png.Dir().Z()) * (minZ-png.Start().Z()) + png.Start().X();
408  zStart = minZ;
409 
410  if( (png.MinPlane() == png.MaxPlane()) ){
411  if ( png.Dir().X() > 0.0 ) xyStart = png.MinV(geo::kX);
412  else if ( png.Dir().X() < 0.0 ) xyStart = png.MaxV(geo::kX);
413  }
414  start[0] = xyStart;
415  start[1] = vert.Y();
416  start[2] = zStart;
417  }//kX
418 
419  if( png.View() == geo::kY ){
420 
421  xyStart = (png.Dir().Y()/png.Dir().Z()) * (minZ-png.Start().Z()) + png.Start().Y();
422  zStart = minZ;
423 
424  if( (png.MinPlane() == png.MaxPlane()) ){
425  if ( png.Dir().Y() > 0.0 ) xyStart = png.MinV(geo::kY);
426  else if ( png.Dir().Y() < 0.0 ) xyStart = png.MaxV(geo::kY);
427  }
428  start[0] = vert.X();
429  start[1] = xyStart;
430  start[2] = zStart;
431  }//kY
432 
433  png.SetStart(start);
434 
435  png2D->push_back(png);
436 
437  util::CreateAssn(*this, evt, *png2D, vtx[0],
438  *png2DAssns, UINT_MAX, fParams.PngInst2D());
439  util::CreateAssn(*this, evt, *png2D, sliceptr,
440  *png2DSlcAssns, UINT_MAX, fParams.PngInst2D());
441 
442  std::vector<rb::PID> pid2(fViewMatchAlg->fUnmatchedScores[iPng]);
443 
444  for(unsigned int i = 0; i < pid2.size(); ++i){
445  PIDCol2->push_back(pid2[i]);
446 
447  util::CreateAssn(*this, evt, *PID2Assns, *PIDCol2, *png2D, UINT_MAX, UINT_MAX,
449  }
450 
451  } // 2d prongs
452 
453  std::cout<<"Preparing orphaned energy"<<std::endl;
454 
455  // Stolen from fuzzyk
456  // handle orphaned hits here
457  // artificially set the orphaned prong start to be the vertex, direction to be entirely z
458  TVector3 start(vert.X(), vert.Y(), vert.Z());
459  TVector3 dir(0.0, 0.0, 1.0);
460 
461  rb::Prong pngO(FindOrphanedHits(sliceptr->AllCells(), prongs), start, dir);
462  pngOrph->push_back(pngO);
463 
464  util::CreateAssn(*this, evt, *pngOrph, vtx[0], *pngOrphAssns, UINT_MAX, fParams.PngInstOrph());
465  util::CreateAssn(*this, evt, *pngOrph, sliceptr, *pngOrphSlcAssns, UINT_MAX, fParams.PngInstOrph());
466 
467  } //slices
468 
469  std::cout<<"Move it move it"<<std::endl;
470 
471  // Put on the event
472  evt.put(std::move(png2D),fParams.PngInst2D());
473  evt.put(std::move(png3D),fParams.PngInst3D());
474  evt.put(std::move(pngOrph),fParams.PngInstOrph());
475 
476  evt.put(std::move(png2DAssns),fParams.PngAssn2D());
477  evt.put(std::move(png3DAssns),fParams.PngAssn3D());
478  evt.put(std::move(pngOrphAssns),fParams.PngAssnOrph());
479 
480  evt.put(std::move(png2DSlcAssns),fParams.PngAssn2D());
481  evt.put(std::move(png3DSlcAssns),fParams.PngAssn3D());
482  evt.put(std::move(pngOrphSlcAssns),fParams.PngAssnOrph());
483 
484  evt.put(std::move(PIDColX),fParams.ScoreInst3DX());
485  evt.put(std::move(PIDColY),fParams.ScoreInst3DY());
486 
487  evt.put(std::move(PIDXAssns),fParams.ScoreAssn3DX());
488  evt.put(std::move(PIDYAssns),fParams.ScoreAssn3DY());
489 
490  evt.put(std::move(PIDCol2),fParams.ScoreInst2D());
491 
492  evt.put(std::move(PID2Assns),fParams.ScoreAssn2D());
493 
494 
495  } //producer
496 
497  //......................................................................
498  /// Creates the 2D prongs in both views
499  void NERDProng::MakeProngs(std::vector<rb::Prong>* prongs,
500  std::vector<std::vector<rb::PID>>* pids,
501  art::Ptr<rb::Cluster> slice,
502  cvn::PixelMap pixelmap,
503  TVector3 vert)
504  {
505 
506  // Run the network prediction and get an ouput mask
507  nerd::PixelPIDMaps pmResult = fNERDEval->predict(pixelmap);
508 
509  // Create clusters from the result
510  std::vector<rb::Cluster> NERDClusterX = MakeCluster(slice, pmResult, pixelmap, geo::kX);
511  std::vector<rb::Cluster> NERDClusterY = MakeCluster(slice, pmResult, pixelmap, geo::kY);
512 
513  // Add X
514  std::vector<rb::PID> temp;
515  for ( unsigned int iPng = 0; iPng < NERDClusterX.size(); ++iPng ){
516  if ( NERDClusterX[iPng].NCell() > 0 ){
517  this->AddProngs(NERDClusterX[iPng], vert, prongs);
518 
519  temp.clear();
520  for(unsigned int i = 0; i < pmResult.fLabelX[iPng].size(); ++i)
521  temp.push_back({GetPDGByNType((NType)i), pmResult.fLabelX[iPng][i]});
522  (*pids).push_back(temp);
523  }
524  }
525 
526  // Add Y
527  for ( unsigned int iPng = 0; iPng < NERDClusterY.size(); ++iPng ){
528  if ( NERDClusterY[iPng].NCell() > 0 ){
529  this->AddProngs(NERDClusterY[iPng], vert, prongs);
530 
531  temp.clear();
532  for(unsigned int i = 0; i < pmResult.fLabelY[iPng].size(); ++i)
533  temp.push_back({GetPDGByNType((NType)i), pmResult.fLabelY[iPng][i]});
534  (*pids).push_back(temp);
535  }
536  }
537 
538  }
539 
540  //......................................................................
541  /// Make clusters from a pmResult
542  std::vector<rb::Cluster> NERDProng::MakeCluster(art::Ptr<rb::Cluster> slice,
543  nerd::PixelPIDMaps pmResult,
544  cvn::PixelMap pixelmap,
546  {
547  std::vector<rb::Cluster> NERDCluster;
548 
549  cvn::Boundary bound = pixelmap.Bound();
550 
551  unsigned int Nprongs2D = pmResult.NClusters(view);
552 
553  // Initialize empty clusters in the appropriate view
554  for(unsigned int i = 0; i < Nprongs2D; ++i)
555  NERDCluster.push_back(rb::Cluster(view));
556 
557  std::vector<art::Ptr<rb::CellHit>> Unclustered;
558 
559  auto fMask = view==geo::kX ? pmResult.fMaskX : pmResult.fMaskY;
560  auto fLabel= view==geo::kX ? pmResult.fLabelX : pmResult.fLabelY;
561 
562  bool isclustered = false;
563 
564  // Loop over cellhits, add them to corresponding prong
565  unsigned int NCells = view==geo::kX ? slice->NXCell() : slice->NYCell();
566  for ( unsigned int hitId = 0; hitId < NCells; ++hitId ){
567 
568  art::Ptr< rb::CellHit > iHit = view==geo::kX ?
569  slice->XCell(hitId) :
570  slice->YCell(hitId) ;
571 
572 
573  // Skip over hits outside the pixelmap
574  if(!bound.IsWithin(iHit->Plane(), iHit->Cell(), iHit->View())){
575  Unclustered.push_back(iHit);
576  continue;
577  }
578 
579  unsigned int pmHitIdx = pixelmap.GlobalToIndexSingle( iHit->Plane(),
580  iHit->Cell(),
581  iHit->View());
582 
583  isclustered = false;
584 
585  for ( unsigned int pngIdx = 0; pngIdx < Nprongs2D; ++pngIdx ){
586  if ( fMask[pmHitIdx][pngIdx] == true ) {
587  isclustered = true;
588  NERDCluster[pngIdx].Add(iHit);
589  }
590  }
591 
592  // if a hit doesn't make into a cluster we handle it separately
593  if(!isclustered){
594  Unclustered.push_back(iHit);
595  }
596 
597  }// hits
598 
599  float maxScore = -1;
601 
602  std::vector<NType> pids;
603 
604  // Determine the identity of each cluster
605  for(unsigned int i = 0; i < fLabel.size(); ++i){
606  maxScore = -1;
607  pid = kBackground;
608  for(unsigned int j = 0; j < fLabel[0].size(); ++j){
609  if(fLabel[i][j] > maxScore){
610  maxScore = fLabel[i][j];
611  pid = (NType)j;
612  }
613  }
614  pids.push_back(pid);
615  }
616 
617  // TODO: This should be something more sophisticated
618  // Compute the distance of each hit to the nearest cluster
619  std::vector<std::pair<art::Ptr<rb::CellHit>, double>> UnclusteredDist;
620  for(unsigned int ix = 0; ix < Unclustered.size(); ++ix){
621  std::pair<int, double> dist = this->FindNearestCluster(Unclustered[ix], NERDCluster);
622  UnclusteredDist.push_back({Unclustered[ix], dist.second});
623  }
624 
625  // Sort hits by distance to nearest cluster
626  std::sort(UnclusteredDist.begin(), UnclusteredDist.end(), [](auto a, auto b) { return a.second < b.second; });
627 
628  // If hit is within threshold add it to the nearest cluster
629  for(unsigned int i = 0; i < UnclusteredDist.size(); ++i){
630  std::pair<int, double> dist = this->FindNearestCluster(UnclusteredDist[i].first, NERDCluster);
631 
632  switch(pids[dist.first]){
633  case kElectron:
634  case kPhoton:
635  if(dist.second < fParams.DistanceThreshold())
636  NERDCluster[dist.first].Add(UnclusteredDist[i].first);
637  break;
638 
639  case kMuon:
640  case kPion:
641  case kProton:
642  if(dist.second < fParams.DistanceThreshold())
643  NERDCluster[dist.first].Add(UnclusteredDist[i].first);
644  break;
645  default:
646  break;
647  }
648  }
649 
650  return NERDCluster;
651  }
652 
653  //......................................................................
654  /// Turn a cluster of hits into a prong
656  TVector3 vert,
657  std::vector<rb::Prong>* prongs)
658  {
659  // Vertex coordinates in this view
660  double zstart = vert.Z();
661  double xystart = cluster.View() == geo::kX ? vert.X() : vert.Y();
662 
663  // Mean cluster coordinates
664  double zclus = cluster.MeanZ();
665  double xyclus = cluster.MeanV(cluster.View());
666 
667  // Prong direction from vertex
668  double zdelta = zclus-zstart;
669  double xydelta = xyclus-xystart;
670 
671  // unit vector
672  double hippo = sqrt(zdelta*zdelta + xydelta*xydelta);
673  zdelta /= hippo;
674  xydelta /= hippo;
675 
676  (*prongs).push_back(rb::Prong(cluster, xystart, zstart, xydelta, zdelta));
677 
678  }
679 
680  //......................................................................
681  /// Compute distance between two cells
683  {
685 
686  // Only calculate if both hits are in the same view
687  if (hit1->View() != hit2->View()) return 100000.;
688  else{
689  double xyz1[3], xyz2[3];
690 
691  unsigned short plane1_idx = hit1->Plane();
692  unsigned short cell1_idx = hit1->Cell();
693  const geo::PlaneGeo* plane1 = fGeo->Plane(plane1_idx);
694  const geo::CellGeo* cell1 = plane1->Cell(cell1_idx);
695  cell1->GetCenter(xyz1);
696 
697  unsigned short plane2_idx = hit2->Plane();
698  unsigned short cell2_idx = hit2->Cell();
699  const geo::PlaneGeo* plane2 = fGeo->Plane(plane2_idx);
700  const geo::CellGeo* cell2 = plane2->Cell(cell2_idx);
701  cell2->GetCenter(xyz2);
702  // Agnostic to whether hits are in X or Y views
703  return util::pythag(xyz1[0]-xyz2[0], xyz1[1]-xyz2[1], xyz1[2]-xyz2[2]);
704  }
705  }
706 
707  //......................................................................
708  // Min dist between a hit and a cluster
710  double minDist = 10000;
711  for (unsigned int i = 0; i < clus.NCell(); i++){
712  double curDist = HitToHitDistance(hit, clus.Cell(i));
713  if (curDist < minDist){
714  minDist = curDist;
715  }
716  }
717  return minDist;
718  }
719 
720  //......................................................................
721  // Find the nearest cluster to a given hit
722  std::pair<int, double> NERDProng::FindNearestCluster(art::Ptr<rb::CellHit> hit, std::vector<rb::Cluster> clus){
723  double minDist = 10000;
724  double minIdx = -1;
725  for (unsigned int i = 0; i < clus.size(); i++){
726  double curDist = MinHitClusterDist(hit, clus[i]);
727  if (curDist < minDist){
728  minDist = curDist;
729  minIdx = i;
730  }
731  }
732  // nearest cluster {index, distance}
733  return {minIdx, minDist};
734  }
735 
736  //......................................................................
737  /// Check if a cvn result passes a score threshold
739  {
740  double cosid = result.fOutput[3];
741 
742  return cosid < threshold;
743  }
744 
745 
746  // function to find orphaned hits
748  {
749  art::PtrVector<rb::CellHit> orphanHits;
750  // loop over all hits in the slice
751  for (unsigned int i=0; i<c.size(); i++){
752  bool inprong = false;
753 
754  // loop over all prongs
755  for (unsigned int j=0; j<prongs.size(); ++j){
756  // check to see if this hit is in the prong
757  for(unsigned int k = 0; k < prongs[j].NCell(); ++k) {
758  if (*c[i] == *prongs[j].Cell(k)) {
759  inprong = true;
760  break; // if hit in prong, break out of loops
761  }
762  }// end hits in prong
763 
764  if (inprong) break;
765 
766  }// end prongs
767 
768  if (!inprong) { // if this slice hit was in no prongs, add to orphan hit prong
769  orphanHits.push_back(c[i]);
770  }
771  }//end hits in slice
772 
773  return orphanHits;
774  }
775 
776 }// namespace nerd
777 ////////////////////////////////////////////////////////////////////////
778 
779 namespace nerd
780 {
782 }
783 ////////////////////////////////////////////////////////////////////////
double MinV(geo::View_t view) const
Definition: Cluster.cxx:454
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
Atom< std::string > ScoreAssn3DY
std::vector< std::vector< bool > > fMaskX
Definition: PixelPIDMaps.h:28
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
std::vector< rb::Prong > fMatchedProngs
Containers for all the 3D prongs produced after matching.
Definition: ViewMatchAlg.h:80
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
Atom< std::string > ScoreInst3DX
double MinZ() const
Definition: Cluster.h:206
std::pair< int, double > FindNearestCluster(art::Ptr< rb::CellHit > hit, std::vector< rb::Cluster > clus)
pdg code and pid value
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Atom< float > DistanceThreshold
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double MaxV(geo::View_t view) const
Definition: Cluster.cxx:500
void produce(art::Event &evt)
std::vector< rb::Cluster > MakeCluster(art::Ptr< rb::Cluster > slice, nerd::PixelPIDMaps pmResult, cvn::PixelMap pixelmap, geo::View_t view)
Builds clusters from a pmResult.
void AddProngs(rb::Cluster cluster, TVector3 vert, std::vector< rb::Prong > *prongs)
Takes a 2D cluster and turns into a prong for view matching.
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
Atom< std::string > PngInst3D
double MeanV(geo::View_t view, rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:546
void Matching()
Function to perform the view matching.
const PlaneGeo * Plane(unsigned int i) const
NERDProng(const Parameters &params)
double MinHitClusterDist(art::Ptr< rb::CellHit > hit, rb::Cluster clus)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
Defines an enumeration for nerd classification.
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
Table< ViewMatchParams > ViewMatch
std::vector< float > fOutput
Vector of outputs from neural net.
Definition: Result.h:30
float abs(float number)
Definition: d0nt_math.hpp:39
Atom< std::string > ScoreAssn3DX
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
Atom< unsigned int > MinHitX
Give every cell in the geometry a unique ID number based on the TGeo path to the node.
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
double dist
Definition: runWimpSim.h:113
std::vector< std::vector< bool > > fMaskY
Definition: PixelPIDMaps.h:29
Table< NERDEvalParams > NERDEvalPSet
Atom< std::string > PngInst2D
std::vector< std::vector< rb::PID > > fMatchedYScores
Definition: ViewMatchAlg.h:82
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Atom< bool > ObeyPreselection
unsigned short Cell() const
Definition: CellHit.h:40
Atom< std::string > PngAssnOrph
Result for CVN.
Molds the pixel maps into the format expected then runs the tensorflow model.
Atom< std::string > CVNLabel
void beginJob()
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
NERDProngParams fParams
Atom< std::string > PngAssn2D
unsigned int NClusters(geo::View_t view)
const double a
Atom< unsigned int > MinHitY
Atom< std::string > VertexLabel
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
std::vector< std::vector< float > > fLabelX
Definition: PixelPIDMaps.h:31
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
art::PtrVector< rb::CellHit > FindOrphanedHits(art::PtrVector< rb::CellHit > c, std::vector< rb::Prong > prongs)
Function to find hits that are not contained in any prong.
Boundary Bound() const
Map boundary.
Definition: PixelMap.h:51
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
double HitToHitDistance(art::Ptr< rb::CellHit > hit1, art::Ptr< rb::CellHit > hit2)
Compute the distance between two hits.
Atom< std::string > ScoreAssn2D
std::vector< rb::Prong > fUnmatchedProngs
Container for all the unmatched 2D prongs left over.
Definition: ViewMatchAlg.h:85
size_type size() const
Definition: PtrVector.h:308
bool PassCVNCut(const cvn::Result &result, float threshold)
Check if a slice passes the cvncut.
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
std::vector< std::vector< rb::PID > > fMatchedXScores
Definition: ViewMatchAlg.h:81
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
void MakeProngs(std::vector< rb::Prong > *prongs, std::vector< std::vector< rb::PID >> *pids, art::Ptr< rb::Cluster > slice, cvn::PixelMap pixelmap, TVector3 vert)
Builds prongs from tensorflow cluster + PID results.
double MaxZ() const
Definition: Cluster.h:219
Result, basic output of CVN neural net.
Definition: Result.h:15
A Cluster with defined start position and direction.
Definition: Prong.h:19
TDirectory * dir
Definition: macro.C:5
Atom< std::string > PixelMapInput
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
Atom< std::string > PreselLabels
Atom< std::string > ScoreInst2D
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:23
Atom< unsigned int > MinHit
Atom< std::string > SliceLabel
const hit & b
Definition: hits.cxx:21
unsigned int GlobalToIndexSingle(const unsigned int &plane, const unsigned int &cell, const unsigned int &view)
Take global plane, cell (detector) and return index in fPE vector.
Definition: PixelMap.cxx:147
std::vector< std::vector< float > > fLabelY
Definition: PixelPIDMaps.h:32
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
std::vector< std::vector< rb::PID > > fUnmatchedScores
Container for all the unmatched 2D scores left over.
Definition: ViewMatchAlg.h:88
Atom< std::string > PngAssn3D
Atom< std::string > PngInstOrph
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Atom< std::string > ScoreInst3DY
Encapsulate the cell geometry.
Definition: CellGeo.h:25
NType
Definition: NERDType.h:12
void LoadProngs(std::vector< rb::Prong > prongs, std::vector< std::vector< rb::PID >> scores)
Function to load the vector of 2d prongs to be matched.
bool IsWithin(const unsigned int &plane, const unsigned int &cell, const unsigned int &view=0)
Definition: Boundary.cxx:40
Encapsulate the geometry of one entire detector (near, far, ndos)
ViewMatchAlg * fViewMatchAlg
An object of ViewMatchAlg to do the 3D view matching.
PixelPIDMaps predict(cvn::PixelMap pm)