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 //////////////////////////////////////////////////////////////////
8 
11 #include "Geometry/Geometry.h"
14 #include "RecoBase/CellHit.h"
15 #include "RecoBase/HitMap.h"
16 #include "RecoBase/Cluster.h"
17 #include "RecoBase/Prong.h"
18 #include "RecoBase/WeightedProng.h"
19 #include "RecoBase/Vertex.h"
20 #include "RecoBase/FilterList.h"
21 #include "RecoBase/PID.h"
22 #include "Utilities/func/MathUtil.h"
23 #include "Utilities/AssociationUtil.h"
24 
30 
31 #include "CVN/func/Result.h"
32 
33 #include "TVector3.h"
34 #include "TTree.h"
35 #include "TH1D.h"
36 #include <cmath>
37 
41 
42 namespace nerd
43 {
45  {
46  template<class T> using Atom = fhicl::Atom<T>;
47  template<class T> using Table = fhicl::Table<T>;
49  using Name = fhicl::Name;
50 
51  // Parameters controling the splitting of clusters
55 
56  Atom <bool> ObeyPreselection{ Name("ObeyPreselection"),
57  Comment("Only run on presel slices") };
59  Comment("Specific presel used") };
60 
62  Comment("Instance label 2D") };
64  Comment("Instance label 2D") };
66  Comment("Instance labels to store results in")};
68  Comment("Instance labels to store associations in")};
69 
71  Comment("Instance label 3D") };
73  Comment("Instance label 3D") };
75  Comment("Instance labels to store results in")};
77  Comment("Instance labels to store associations in")};
78 
80  Comment("Instance label orphaned prong") };
82  Comment("Instance label orphaned prong") };
83 
85  Comment("Instance Label for X-View Scores") };
87  Comment("Instance Label for X-View Scores") };
89  Comment("Instance Label for Y-View Scores") };
91  Comment("Instance Label for Y-View Scores") };
93  Comment("Instance Label for 2D Scores") };
95  Comment("Instance Label for 2D Scores") };
96 
98  Comment("Minimum # of hits") };
100  Comment("Minimum # of hits X") };
102  Comment("Minimum # of hits Y") };
104  Comment("CVN Score Label") };
106  Comment("CVN cut for evaluation") };
107 
108  Atom <float> DistanceThreshold{ Name("DistanceThreshold"),
109  Comment("Maximum distance in cm to include unclustered hits") };
110 
111 
115 
116  };
117 
118 
119  class NERDProng : public art::EDProducer
120  {
121  public:
122  // Allows 'nova --print-description' to work
124 
125  explicit NERDProng(const Parameters& params);
126  virtual ~NERDProng();
127 
128  void beginJob();
129  void produce(art::Event &evt);
130 
131  private:
133 
134  // An object to run the segmentation graph
136 
137  /// An object of ViewMatchAlg to do the 3D view matching
139 
140  /// An object of WeightSharedHitsAlg to do the hit sharing
142 
144 
145  /// \brief Builds prongs from tensorflow cluster + PID results
146  /// \param prongs A pointer vector for the resulting 2d prongs
147  /// \param slice The input rb::Cluster object to use in clustering
148  /// \param vert The vertex for associations
149  void MakeProngs( std::vector<rb::Prong>* prongs,
150  std::vector<std::vector<rb::PID>>* pids,
151  art::Ptr<rb::Cluster> slice,
152  cvn::PixelMap pixelmap,
153  TVector3 vert);
154 
155  /// \brief Builds clusters from a pmResult
156  std::vector<rb::Cluster> MakeCluster(art::Ptr<rb::Cluster> slice,
157  nerd::PixelPIDMaps pmResult,
158  cvn::PixelMap pixelmap,
159  geo::View_t view);
160 
161  /// \brief Takes a 2D cluster and turns into a prong for view matching
162  void AddProngs(rb::Cluster cluster,
163  TVector3 vert,
164  std::vector<rb::Prong>* prongs);
165 
166  /// \brief Compute the distance between two hits
167  double HitToHitDistance(art::Ptr<rb::CellHit> hit1, art::Ptr<rb::CellHit> hit2);
168 
169  // \brief Min dist between a hit and a physics slice
170  double MinHitClusterDist(art::Ptr<rb::CellHit> hit, rb::Cluster clus);
171 
172  // \brief Find the nearest cluster to a hit
173  std::pair<int, double> FindNearestCluster(art::Ptr<rb::CellHit> hit, std::vector<rb::Cluster> clus);
174 
175  /// \brief Function to find hits that are not contained in any prong
176  art::PtrVector<rb::CellHit> FindOrphanedHits(art::PtrVector<rb::CellHit> c, std::vector<rb::Prong> prongs);
177 
178  /// \brief Check if a slice passes the cvncut
179  bool PassCVNCut(const cvn::Result& result, float threshold);
180 
181  };
182 }
183 ////////////////////////////////////////////////////////////////////////////////
184 namespace nerd
185 {
186 
188  : fSliceToken( consumes<std::vector<rb::Cluster>>(params().SliceLabel()) ),
189  fNERDEval(new NERDEval(params().NERDEvalPSet())),
190  fViewMatchAlg(new ViewMatchAlg(params().ViewMatch())),
191  fWeightSharedHitsAlg(new WeightSharedHitsAlg(params().WeightSharedHits())),
192  fParams(params())
193  {
194  this->produces< std::vector<rb::Prong> >(fParams.PngInst2D());
195  this->produces< std::vector<rb::Prong> >(fParams.PngInst3D());
196  this->produces< std::vector<rb::Prong> >(fParams.PngInstOrph());
197  this->produces< std::vector<rb::WeightedProng> >(fParams.PngInst3DWt());
198  this->produces< std::vector<rb::WeightedProng> >(fParams.PngInst2DWt());
199 
200  this->produces< art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssn2D());
201  this->produces< art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssn3D());
202  this->produces< art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssnOrph());
203  this->produces< art::Assns<rb::WeightedProng, rb::Vertex> >(fParams.PngAssn3DWt());
204  this->produces< art::Assns<rb::WeightedProng, rb::Vertex> >(fParams.PngAssn2DWt());
205 
206  this->produces< art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssn2D());
207  this->produces< art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssn3D());
208  this->produces< art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssnOrph());
209  this->produces< art::Assns<rb::WeightedProng, rb::Cluster> >(fParams.PngAssn3DWt());
210  this->produces< art::Assns<rb::WeightedProng, rb::Cluster> >(fParams.PngAssn2DWt());
211  this->produces< art::Assns<rb::Prong, rb::WeightedProng> >(fParams.PngAssn3DWt());
212  this->produces< art::Assns<rb::Prong, rb::WeightedProng> >(fParams.PngAssn2DWt());
213 
214  this->produces< std::vector<rb::PID> >(fParams.ScoreInst2D());
215  this->produces< std::vector<rb::PID> >(fParams.ScoreInst3DX());
216  this->produces< std::vector<rb::PID> >(fParams.ScoreInst3DY());
217 
218  this->produces< art::Assns<rb::PID,rb::Prong> >(fParams.ScoreAssn2D());
219  this->produces< art::Assns<rb::PID,rb::Prong> >(fParams.ScoreAssn3DX());
220  this->produces< art::Assns<rb::PID,rb::Prong> >(fParams.ScoreAssn3DY());
221  }
222 
223  //......................................................................
224 
226  {
227  }
228 
229  //......................................................................
230 
232  {
233  if(fNERDEval) delete fNERDEval;
234  if(fViewMatchAlg) delete fViewMatchAlg;
236 
237  }
238 
239  //......................................................................
240 
242  {
244 
245  // New products
246  // 2D Clusters
247  std::unique_ptr< std::vector<rb::Prong> > png2D(new std::vector<rb::Prong>);
248 
249  std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> >
250  png2DSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
251  std::unique_ptr< art::Assns<rb::Prong, rb::Vertex> >
252  png2DAssns(new art::Assns<rb::Prong, rb::Vertex>);
253 
254  std::unique_ptr< std::vector<rb::PID> >
255  PIDCol2(new std::vector<rb::PID>);
256 
257  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
258  PID2Assns(new art::Assns<rb::PID, rb::Prong>);
259 
260  // 3D Clusters
261  std::unique_ptr< std::vector<rb::Prong> > png3D(new std::vector<rb::Prong>);
262 
263  std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> >
264  png3DSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
265  std::unique_ptr< art::Assns<rb::Prong, rb::Vertex> >
266  png3DAssns(new art::Assns<rb::Prong, rb::Vertex>);
267 
268  std::unique_ptr< std::vector<rb::PID> >
269  PIDColX(new std::vector<rb::PID>);
270  std::unique_ptr< std::vector<rb::PID> >
271  PIDColY(new std::vector<rb::PID>);
272 
273  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
274  PIDXAssns(new art::Assns<rb::PID, rb::Prong>);
275  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
276  PIDYAssns(new art::Assns<rb::PID, rb::Prong>);
277 
278  // 2D Weighted Prong Clusters
279  std::unique_ptr< std::vector<rb::WeightedProng> > png2DWt(new std::vector<rb::WeightedProng>);
280  std::unique_ptr< art::Assns<rb::WeightedProng, rb::Vertex> > png2DWtAssns(new art::Assns<rb::WeightedProng, rb::Vertex>);
281  std::unique_ptr< art::Assns<rb::WeightedProng, rb::Cluster> > png2DWtSlcAssns(new art::Assns<rb::WeightedProng, rb::Cluster>);
282 
283  //3D Weighted Prong Clusters
284  std::unique_ptr< std::vector<rb::WeightedProng> > png3DWt(new std::vector<rb::WeightedProng>);
285  std::unique_ptr< art::Assns<rb::WeightedProng, rb::Vertex> > png3DWtAssns(new art::Assns<rb::WeightedProng, rb::Vertex>);
286  std::unique_ptr< art::Assns<rb::WeightedProng, rb::Cluster> > png3DWtSlcAssns(new art::Assns<rb::WeightedProng, rb::Cluster>);
287 
288  std::unique_ptr<art::Assns<rb::Prong, rb::WeightedProng> > png3Dpng3DWtAssns(new art::Assns<rb::Prong, rb::WeightedProng>);
289  std::unique_ptr<art::Assns<rb::Prong, rb::WeightedProng> > png2Dpng2DWtAssns(new art::Assns<rb::Prong, rb::WeightedProng>);
290 
291 
292  // Orphaned Energy
293  std::unique_ptr< std::vector<rb::Prong> > pngOrph(new std::vector<rb::Prong>);
294  std::unique_ptr<art::Assns<rb::Prong, rb::Vertex> > pngOrphAssns(new art::Assns<rb::Prong, rb::Vertex>);
295  std::unique_ptr<art::Assns<rb::Prong, rb::Cluster> > pngOrphSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
296 
297 
298 
299  // Get from evt
301  evt.getByToken(fSliceToken, sliceHandle);
302 
303  art::FindManyP<cvn::PixelMap> fmPixMap( sliceHandle,
304  evt,
306 
307  art::FindManyP<rb::Vertex> fmVtx( sliceHandle,
308  evt,
309  fParams.VertexLabel());
310 
311  art::FindManyP<cvn::Result> fmCVN( sliceHandle,
312  evt,
313  fParams.CVNLabel());
314 
315  // Local
316  std::vector<rb::Prong> prongs;
317  std::vector<std::vector<rb::PID>> pids;
318 
319  std::vector<size_t> png3DIndx;
320  std::vector<size_t> png2DIndx;
321  std::vector<size_t> png3DWtIndx;
322  std::vector<size_t> png2DWtIndx;
323 
324  // Loop over slices
325  for ( unsigned int iSlice=0; iSlice<sliceHandle->size(); ++iSlice ) {
326  //std::cout<<"Looking at slice "<<iSlice<<std::endl;
327 
328  std::vector<rb::WeightedProng> pngWt;
329 
330  // Preselection
331  if( fParams.ObeyPreselection() &&
332  rb::IsFiltered(evt,sliceHandle,iSlice,fParams.PreselLabels()) )
333  continue;
334 
335  // CVN selection
336  std::vector<art::Ptr<cvn::Result>> cvnresult = fmCVN.at(iSlice);
337  if(!cvnresult.empty()){
338  if (!PassCVNCut(*cvnresult[0], fParams.CVNCut()))
339  continue;
340  }
341 
342  art::Ptr<rb::Cluster> sliceptr(sliceHandle, iSlice);
343 
344  if ( sliceptr->IsNoise() ) continue;
345  if ( sliceptr->NXCell() < fParams.MinHitX() ) continue;
346  if ( sliceptr->NYCell() < fParams.MinHitY() ) continue;
347  if ( sliceptr->NCell() < fParams.MinHit() ) continue;
348 
349  // Find any verticies that go with this slice
350  std::vector<art::Ptr<rb::Vertex> > vtx = fmVtx.at(iSlice);
351  if( vtx.size() != 1 ) continue;
352 
353  TVector3 vert = vtx[0]->GetXYZ();
354  // Same vertex check as FuzzyK
355  if( std::isnan(vert.X()) ||
356  std::isinf(vert.X()) || std::abs(vert.X())>10000) continue;
357  if( std::isnan(vert.Y()) ||
358  std::isinf(vert.Y()) || std::abs(vert.Y())>10000) continue;
359  if( std::isnan(vert.Z()) ||
360  std::isinf(vert.Z()) || std::abs(vert.Z())>10000) continue;
361 
362  // Find the PixelMap that goes with this slice
363  const std::vector<art::Ptr<cvn::PixelMap> > pixMaps = fmPixMap.at(iSlice);
364  if( pixMaps.size() != 1 ) continue;
365 
366  art::Ptr<cvn::PixelMap> pixmap = pixMaps[0];
367  // make the 2D prong candidates for this slice
368  prongs.clear();
369  pids.clear();
370 
371  this->MakeProngs(&prongs, &pids, sliceptr, *pixmap, vert);
372 
373  // perform the view matching with this set of 2D prongs
374  //std::cout<<"Beginning Matching"<<std::endl;
375  fViewMatchAlg->LoadProngs(prongs, pids);
377 
378  //std::cout<<"2D Prongs: "<<prongs.size()<<std::endl;
379  //std::cout<<"3D Prongs: "<<fViewMatchAlg->fMatchedProngs.size()<<std::endl;
380  //std::cout<<"Unmatched Prongs: "<<fViewMatchAlg->fUnmatchedProngs.size()<<std::endl;
381 
382  //std::cout<<"Preparing 3D prongs"<<std::endl;
383 
384  // loop over 3D matched prongs to push into the event
385  for( unsigned int iPng = 0; iPng < fViewMatchAlg->fMatchedProngs.size(); ++iPng ){
387 
388  //move the prong start point from the elastic arms vertex to the closest hit
389  TVector3 start;
390  double minZ = png.MinZ();
391  double maxZ = png.MaxZ();
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  start[0] = (png.Dir().X()/png.Dir().Z())*(minZ-png.Start().Z()) + png.Start().X();
404  start[1] = (png.Dir().Y()/png.Dir().Z())*(minZ-png.Start().Z()) + png.Start().Y();
405  start[2] = minZ;
406 
407  png.SetStart(start);
408 
409  png3D->push_back(png);
410 
411  util::CreateAssn(*this, evt, *png3D, vtx[0], *png3DAssns, UINT_MAX, fParams.PngInst3D());
412  util::CreateAssn(*this, evt, *png3D, sliceptr, *png3DSlcAssns, UINT_MAX, fParams.PngInst3D());
413 
414  // handle weighted prongs
415  pngWt.push_back(png);
416 
417  // Prepare for associating the pid scores
418  std::vector<rb::PID> pidX(fViewMatchAlg->fMatchedXScores[iPng]);
419  std::vector<rb::PID> pidY(fViewMatchAlg->fMatchedYScores[iPng]);
420 
421  for(unsigned int i = 0; i < pidX.size(); ++i){
422  PIDColX->push_back(pidX[i]);
423  PIDColY->push_back(pidY[i]);
424 
425  util::CreateAssn(*this, evt, *PIDXAssns, *PIDColX, *png3D, UINT_MAX, UINT_MAX,
427  util::CreateAssn(*this, evt, *PIDYAssns, *PIDColY, *png3D, UINT_MAX, UINT_MAX,
429  }
430 
431  } // 3d prongs
432 
433  // track the indices so png:pngWt associations can be made
434  if (fViewMatchAlg->fMatchedProngs.size() > 0){
435  for(unsigned int idx=(*png3D).size()-fViewMatchAlg->fMatchedProngs.size();idx<(*png3D).size();++idx){
436  png3DIndx.push_back(idx);
437  }
438  }
439 
440  //std::cout<<"Preparing 2D prongs"<<std::endl;
441 
442  //Loop over the 2d prongs for vtx related calculations
443  for( unsigned int iPng=0; iPng < fViewMatchAlg->fUnmatchedProngs.size(); ++iPng ){
445 
446  double xyStart = 0.0;
447  double zStart = 0.0;
448  double minZ = png.MinZ();
449  double maxZ = png.MaxZ();
450 
451  if ( (png.Dir().Z() > 0.0) && (minZ < png.Start().Z()) )
452  minZ = png.Start().Z();
453  if ( (png.Dir().Z() < 0.0) && (maxZ > png.Start().Z()) )
454  maxZ = png.Start().Z();
455 
456  if ( png.Dir().Z() < 0.0 ){
457  double temp = maxZ;
458  maxZ = minZ;
459  minZ = temp;
460  }
461 
462  TVector3 start;
463 
464  if( png.View() == geo::kX ){
465 
466  xyStart = (png.Dir().X()/png.Dir().Z()) * (minZ-png.Start().Z()) + png.Start().X();
467  zStart = minZ;
468 
469  if( (png.MinPlane() == png.MaxPlane()) ){
470  if ( png.Dir().X() > 0.0 ) xyStart = png.MinV(geo::kX);
471  else if ( png.Dir().X() < 0.0 ) xyStart = png.MaxV(geo::kX);
472  }
473  start[0] = xyStart;
474  start[1] = vert.Y();
475  start[2] = zStart;
476  }//kX
477 
478  if( png.View() == geo::kY ){
479 
480  xyStart = (png.Dir().Y()/png.Dir().Z()) * (minZ-png.Start().Z()) + png.Start().Y();
481  zStart = minZ;
482 
483  if( (png.MinPlane() == png.MaxPlane()) ){
484  if ( png.Dir().Y() > 0.0 ) xyStart = png.MinV(geo::kY);
485  else if ( png.Dir().Y() < 0.0 ) xyStart = png.MaxV(geo::kY);
486  }
487  start[0] = vert.X();
488  start[1] = xyStart;
489  start[2] = zStart;
490  }//kY
491 
492  png.SetStart(start);
493 
494  png2D->push_back(png);
495 
496  // handle weighted prongs
497  pngWt.push_back(png);
498 
499  util::CreateAssn(*this, evt, *png2D, vtx[0],
500  *png2DAssns, UINT_MAX, fParams.PngInst2D());
501  util::CreateAssn(*this, evt, *png2D, sliceptr,
502  *png2DSlcAssns, UINT_MAX, fParams.PngInst2D());
503 
504  std::vector<rb::PID> pid2(fViewMatchAlg->fUnmatchedScores[iPng]);
505 
506  for(unsigned int i = 0; i < pid2.size(); ++i){
507  PIDCol2->push_back(pid2[i]);
508 
509  util::CreateAssn(*this, evt, *PID2Assns, *PIDCol2, *png2D, UINT_MAX, UINT_MAX,
511  }
512 
513  } // 2d prongs
514 
515  // prepare for png:pngWt assns
516  if (fViewMatchAlg->fUnmatchedProngs.size() > 0){
517  for(unsigned int idx=(*png2D).size()-fViewMatchAlg->fUnmatchedProngs.size();idx<(*png2D).size();++idx){
518  png2DIndx.push_back(idx);
519  }
520  }
521 
522  // do hit sharing / weighting for pngWt
524  for(unsigned int k=0; k<pngWt.size(); ++k){
525  if(pngWt[k].Is3D()) {
526  png3DWt->push_back(pngWt[k]);
527  util::CreateAssn(*this, evt, *png3DWt, vtx[0], *png3DWtAssns, UINT_MAX, fParams.PngInst3DWt());
528  util::CreateAssn(*this, evt, *png3DWt, sliceptr, *png3DWtSlcAssns, UINT_MAX, fParams.PngInst3DWt());
529 
530  }
531  else {
532  png2DWt->push_back(pngWt[k]);
533  util::CreateAssn(*this, evt, *png2DWt, vtx[0], *png2DWtAssns, UINT_MAX, fParams.PngInst2DWt());
534  util::CreateAssn(*this, evt, *png2DWt, sliceptr, *png2DWtSlcAssns, UINT_MAX, fParams.PngInst2DWt());
535  }
536  }
537 
538  // prepare for png:pngWt assns
539  if (fViewMatchAlg->fMatchedProngs.size() > 0){
540  for(unsigned int idx=(*png3DWt).size()-fViewMatchAlg->fMatchedProngs.size();idx<(*png3DWt).size();++idx){
541  png3DWtIndx.push_back(idx);
542  }
543  }
544 
545  // prepare for png:pngWt assns
546  if (fViewMatchAlg->fUnmatchedProngs.size() > 0){
547  for(unsigned int idx=(*png2DWt).size()-fViewMatchAlg->fUnmatchedProngs.size();idx<(*png2DWt).size();++idx){
548  png2DWtIndx.push_back(idx);
549  }
550  }
551 
552 
553  //std::cout<<"Preparing orphaned energy"<<std::endl;
554 
555  // Stolen from fuzzyk
556  // handle orphaned hits here
557  // artificially set the orphaned prong start to be the vertex, direction to be entirely z
558  TVector3 start(vert.X(), vert.Y(), vert.Z());
559  TVector3 dir(0.0, 0.0, 1.0);
560 
561  rb::Prong pngO(FindOrphanedHits(sliceptr->AllCells(), prongs), start, dir);
562  pngOrph->push_back(pngO);
563 
564  util::CreateAssn(*this, evt, *pngOrph, vtx[0], *pngOrphAssns, UINT_MAX, fParams.PngInstOrph());
565  util::CreateAssn(*this, evt, *pngOrph, sliceptr, *pngOrphSlcAssns, UINT_MAX, fParams.PngInstOrph());
566 
567  } //slices
568 
569  // move the prong-to-weighted-prong associations here and model it like the fake vertex associations
570  for(unsigned int i=0; i<png3DIndx.size(); ++i){
571  util::CreateAssn(*this, evt, *png3Dpng3DWtAssns, *png3D, *png3DWt, png3DWtIndx[i], png3DIndx[i], fParams.PngInst3D(), fParams.PngInst3DWt());
572  }
573 
574  for(unsigned int i=0; i<png2DIndx.size(); ++i){
575  util::CreateAssn(*this, evt, *png2Dpng2DWtAssns, *png2D, *png2DWt, png2DWtIndx[i], png2DIndx[i], fParams.PngInst2D(), fParams.PngInst2DWt());
576  }
577 
578  //std::cout<<"Move it move it"<<std::endl;
579 
580  // Put on the event
581  evt.put(std::move(png2D),fParams.PngInst2D());
582  evt.put(std::move(png3D),fParams.PngInst3D());
583  evt.put(std::move(pngOrph),fParams.PngInstOrph());
584  evt.put(std::move(png2DWt),fParams.PngInst2DWt());
585  evt.put(std::move(png3DWt),fParams.PngInst3DWt());
586 
587  evt.put(std::move(png2DAssns),fParams.PngAssn2D());
588  evt.put(std::move(png3DAssns),fParams.PngAssn3D());
589  evt.put(std::move(pngOrphAssns),fParams.PngAssnOrph());
590  evt.put(std::move(png3DWtAssns),fParams.PngAssn3DWt());
591  evt.put(std::move(png2DWtAssns),fParams.PngAssn2DWt());
592 
593  evt.put(std::move(png2DSlcAssns),fParams.PngAssn2D());
594  evt.put(std::move(png3DSlcAssns),fParams.PngAssn3D());
595  evt.put(std::move(pngOrphSlcAssns),fParams.PngAssnOrph());
596  evt.put(std::move(png3DWtSlcAssns),fParams.PngAssn3DWt());
597  evt.put(std::move(png2DWtSlcAssns),fParams.PngAssn2DWt());
598 
599  evt.put(std::move(PIDColX),fParams.ScoreInst3DX());
600  evt.put(std::move(PIDColY),fParams.ScoreInst3DY());
601 
602  evt.put(std::move(PIDXAssns),fParams.ScoreAssn3DX());
603  evt.put(std::move(PIDYAssns),fParams.ScoreAssn3DY());
604 
605  evt.put(std::move(png3Dpng3DWtAssns),fParams.PngAssn3DWt());
606  evt.put(std::move(png2Dpng2DWtAssns),fParams.PngAssn2DWt());
607 
608  evt.put(std::move(PIDCol2),fParams.ScoreInst2D());
609 
610  evt.put(std::move(PID2Assns),fParams.ScoreAssn2D());
611 
612 
613  } //producer
614 
615  //......................................................................
616  /// Creates the 2D prongs in both views
617  void NERDProng::MakeProngs(std::vector<rb::Prong>* prongs,
618  std::vector<std::vector<rb::PID>>* pids,
619  art::Ptr<rb::Cluster> slice,
620  cvn::PixelMap pixelmap,
621  TVector3 vert)
622  {
623  // Run the network prediction and get an ouput mask
624  nerd::PixelPIDMaps pmResult = fNERDEval->predict(pixelmap);
625 
626  // Create clusters from the result
627  std::vector<rb::Cluster> NERDClusterX = MakeCluster(slice, pmResult, pixelmap, geo::kX);
628  std::vector<rb::Cluster> NERDClusterY = MakeCluster(slice, pmResult, pixelmap, geo::kY);
629 
630  // Add X
631  std::vector<rb::PID> temp;
632  auto xlabels = pmResult.GetViewLabel(geo::kX);
633  for ( unsigned int iPng = 0; iPng < NERDClusterX.size(); ++iPng ){
634  if ( NERDClusterX[iPng].NCell() > 0 ){
635  this->AddProngs(NERDClusterX[iPng], vert, prongs);
636 
637  temp.clear();
638  for(unsigned int i = 0; i < xlabels[iPng].size(); ++i)
639  temp.push_back({GetPDGByNType((NType)i), xlabels[iPng][i]});
640  (*pids).push_back(temp);
641  }
642  }
643 
644  // Add Y
645  auto ylabels = pmResult.GetViewLabel(geo::kY);
646  for ( unsigned int iPng = 0; iPng < NERDClusterY.size(); ++iPng ){
647  if ( NERDClusterY[iPng].NCell() > 0 ){
648  this->AddProngs(NERDClusterY[iPng], vert, prongs);
649 
650  temp.clear();
651  for(unsigned int i = 0; i < ylabels[iPng].size(); ++i)
652  temp.push_back({GetPDGByNType((NType)i), ylabels[iPng][i]});
653  (*pids).push_back(temp);
654  }
655  }
656 
657  }
658 
659  //......................................................................
660  /// Make clusters from a pmResult
661  std::vector<rb::Cluster> NERDProng::MakeCluster(art::Ptr<rb::Cluster> slice,
662  nerd::PixelPIDMaps pmResult,
663  cvn::PixelMap pixelmap,
665  {
666  std::vector<rb::Cluster> NERDCluster;
667 
668  cvn::Boundary bound = pixelmap.Bound();
669 
670  unsigned int Nprongs2D = pmResult.NClusters(view);
671 
672  // Initialize empty clusters in the appropriate view
673  for(unsigned int i = 0; i < Nprongs2D; ++i)
674  NERDCluster.push_back(rb::Cluster(view));
675 
676  std::vector<art::Ptr<rb::CellHit>> Unclustered;
677 
678  //auto fMask = view==geo::kX ? pmResult.fMaskX : pmResult.fMaskY;
679  //auto fLabel= view==geo::kX ? pmResult.fLabelX : pmResult.fLabelY;
680  auto fMask = pmResult.GetViewMask(view);
681  auto fLabel = pmResult.GetViewLabel(view);
682 
683  bool isclustered = false;
684 
685  // Loop over cellhits, add them to corresponding prong
686  unsigned int NCells = view==geo::kX ? slice->NXCell() : slice->NYCell();
687  for ( unsigned int hitId = 0; hitId < NCells; ++hitId ){
688 
689  art::Ptr< rb::CellHit > iHit = view==geo::kX ?
690  slice->XCell(hitId) :
691  slice->YCell(hitId) ;
692 
693 
694  // Skip over hits outside the pixelmap
695  if(!bound.IsWithin(iHit->Plane(), iHit->Cell(), iHit->View())){
696  Unclustered.push_back(iHit);
697  continue;
698  }
699 
700  unsigned int pmHitIdx = pixelmap.GlobalToIndexSingle( iHit->Plane(),
701  iHit->Cell(),
702  iHit->View());
703 
704  isclustered = false;
705 
706  for ( unsigned int pngIdx = 0; pngIdx < Nprongs2D; ++pngIdx ){
707  if ( fMask[pmHitIdx][pngIdx] == true ) {
708  isclustered = true;
709  NERDCluster[pngIdx].Add(iHit);
710  }
711  }
712 
713  // if a hit doesn't make into a cluster we handle it separately
714  if(!isclustered){
715  Unclustered.push_back(iHit);
716  }
717 
718  }// hits
719 
720  float maxScore = -1;
722 
723  std::vector<NType> pids;
724 
725  // Determine the identity of each cluster
726  for(unsigned int i = 0; i < fLabel.size(); ++i){
727  maxScore = -1;
728  pid = kBackground;
729  for(unsigned int j = 0; j < fLabel[0].size(); ++j){
730  if(fLabel[i][j] > maxScore){
731  maxScore = fLabel[i][j];
732  pid = (NType)j;
733  }
734  }
735  pids.push_back(pid);
736  }
737 
738  // TODO: This should be something more sophisticated
739  // Compute the distance of each hit to the nearest cluster
740  std::vector<std::pair<art::Ptr<rb::CellHit>, double>> UnclusteredDist;
741  for(unsigned int ix = 0; ix < Unclustered.size(); ++ix){
742  std::pair<int, double> dist = this->FindNearestCluster(Unclustered[ix], NERDCluster);
743  UnclusteredDist.push_back({Unclustered[ix], dist.second});
744  }
745 
746  // Sort hits by distance to nearest cluster
747  std::sort(UnclusteredDist.begin(), UnclusteredDist.end(), [](auto a, auto b) { return a.second < b.second; });
748 
749  // If hit is within threshold add it to the nearest cluster
750  for(unsigned int i = 0; i < UnclusteredDist.size(); ++i){
751  std::pair<int, double> dist = this->FindNearestCluster(UnclusteredDist[i].first, NERDCluster);
752 
753  switch(pids[dist.first]){
754  case kElectron:
755  case kPhoton:
756  if(dist.second < fParams.DistanceThreshold())
757  NERDCluster[dist.first].Add(UnclusteredDist[i].first);
758  break;
759 
760  case kMuon:
761  case kPion:
762  case kProton:
763  if(dist.second < fParams.DistanceThreshold())
764  NERDCluster[dist.first].Add(UnclusteredDist[i].first);
765  break;
766  default:
767  break;
768  }
769  }
770 
771  return NERDCluster;
772  }
773 
774  //......................................................................
775  /// Turn a cluster of hits into a prong
777  TVector3 vert,
778  std::vector<rb::Prong>* prongs)
779  {
780  // Vertex coordinates in this view
781  double zstart = vert.Z();
782  double xystart = cluster.View() == geo::kX ? vert.X() : vert.Y();
783 
784  // Mean cluster coordinates
785  double zclus = cluster.MeanZ();
786  double xyclus = cluster.MeanV(cluster.View());
787 
788  // Prong direction from vertex
789  double zdelta = zclus-zstart;
790  double xydelta = xyclus-xystart;
791 
792  // unit vector
793  double hippo = sqrt(zdelta*zdelta + xydelta*xydelta);
794  zdelta /= hippo;
795  xydelta /= hippo;
796 
797  (*prongs).push_back(rb::Prong(cluster, xystart, zstart, xydelta, zdelta));
798 
799  }
800 
801  //......................................................................
802  /// Compute distance between two cells
804  {
806 
807  // Only calculate if both hits are in the same view
808  if (hit1->View() != hit2->View()) return 100000.;
809  else{
810  double xyz1[3], xyz2[3];
811 
812  unsigned short plane1_idx = hit1->Plane();
813  unsigned short cell1_idx = hit1->Cell();
814  const geo::PlaneGeo* plane1 = fGeo->Plane(plane1_idx);
815  const geo::CellGeo* cell1 = plane1->Cell(cell1_idx);
816  cell1->GetCenter(xyz1);
817 
818  unsigned short plane2_idx = hit2->Plane();
819  unsigned short cell2_idx = hit2->Cell();
820  const geo::PlaneGeo* plane2 = fGeo->Plane(plane2_idx);
821  const geo::CellGeo* cell2 = plane2->Cell(cell2_idx);
822  cell2->GetCenter(xyz2);
823  // Agnostic to whether hits are in X or Y views
824  return util::pythag(xyz1[0]-xyz2[0], xyz1[1]-xyz2[1], xyz1[2]-xyz2[2]);
825  }
826  }
827 
828  //......................................................................
829  // Min dist between a hit and a cluster
831  double minDist = 10000;
832  for (unsigned int i = 0; i < clus.NCell(); i++){
833  double curDist = HitToHitDistance(hit, clus.Cell(i));
834  if (curDist < minDist){
835  minDist = curDist;
836  }
837  }
838  return minDist;
839  }
840 
841  //......................................................................
842  // Find the nearest cluster to a given hit
843  std::pair<int, double> NERDProng::FindNearestCluster(art::Ptr<rb::CellHit> hit, std::vector<rb::Cluster> clus){
844  double minDist = 10000;
845  double minIdx = -1;
846  for (unsigned int i = 0; i < clus.size(); i++){
847  double curDist = MinHitClusterDist(hit, clus[i]);
848  if (curDist < minDist){
849  minDist = curDist;
850  minIdx = i;
851  }
852  }
853  // nearest cluster {index, distance}
854  return {minIdx, minDist};
855  }
856 
857  //......................................................................
858  /// Check if a cvn result passes a score threshold
860  {
861  double cosid = result.fOutput[3];
862 
863  return cosid < threshold;
864  }
865 
866 
867  // function to find orphaned hits
869  {
870  art::PtrVector<rb::CellHit> orphanHits;
871  // loop over all hits in the slice
872  for (unsigned int i=0; i<c.size(); i++){
873  bool inprong = false;
874 
875  // loop over all prongs
876  for (unsigned int j=0; j<prongs.size(); ++j){
877  // check to see if this hit is in the prong
878  for(unsigned int k = 0; k < prongs[j].NCell(); ++k) {
879  if (*c[i] == *prongs[j].Cell(k)) {
880  inprong = true;
881  break; // if hit in prong, break out of loops
882  }
883  }// end hits in prong
884 
885  if (inprong) break;
886 
887  }// end prongs
888 
889  if (!inprong) { // if this slice hit was in no prongs, add to orphan hit prong
890  orphanHits.push_back(c[i]);
891  }
892  }//end hits in slice
893 
894  return orphanHits;
895  }
896 
897 }// namespace nerd
898 ////////////////////////////////////////////////////////////////////////
899 
900 namespace nerd
901 {
903 }
904 ////////////////////////////////////////////////////////////////////////
double MinV(geo::View_t view) const
Definition: Cluster.cxx:454
Table< nerd::WeightSharedHitsParams > WeightSharedHits
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
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
std::vector< std::vector< bool > > GetViewMask(geo::View_t view)
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
Atom< std::string > PngInst3DWt
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.
Definition: NERDEvaluate.cxx:9
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
Atom< std::string > PngInst2DWt
double dist
Definition: runWimpSim.h:113
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< std::string > PngAssn2DWt
Atom< unsigned int > MinHitY
std::vector< std::vector< float > > GetViewLabel(geo::View_t view)
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
int evt
WeightSharedHitsAlg * fWeightSharedHitsAlg
An object of WeightSharedHitsAlg to do the hit sharing.
void GetWeights(std::vector< rb::WeightedProng > &prong)
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.
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: structs.h:12
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 > PngAssn3DWt
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
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)