FuzzyKVertex_module.cc
Go to the documentation of this file.
1 ///
2 /// \file FuzzyKVertex.cxx
3 /// \brief Fuzzy k-means to optimize elastic arms vertex
4 /// \author eniner@indiana.edu
5 /// \version $Id: FuzzyKVertex.cxx,v 1.15 2012-11-07 18:20:04 edniner Exp $
6 ///
10 
13 #include "Geometry/Geometry.h"
16 #include "RecoBase/CellHit.h"
17 #include "RecoBase/HitMap.h"
18 #include "RecoBase/Cluster.h"
19 #include "RecoBase/Prong.h"
20 #include "RecoBase/WeightedProng.h"
21 #include "RecoBase/Vertex.h"
22 #include "RecoBase/FilterList.h"
23 #include "Utilities/func/MathUtil.h"
24 #include "Utilities/AssociationUtil.h"
25 #include "Simulation/Simulation.h"
29 #include "Simulation/Particle.h"
30 #include "MCCheater/BackTracker.h"
31 
38 
39 #include "TVector3.h"
40 #include "TTree.h"
41 #include "TH1D.h"
42 #include <cmath>
43 
44 /// Fuzzy k-Means prong-finding algorithm
45 namespace fuzz {
47  {
48  template<class T> using Atom = fhicl::Atom<T>;
49  template<class T> using Table = fhicl::Table<T>;
51  using Name = fhicl::Name;
52 
53  Atom<unsigned int> MinHit {Name("MinHit" ), Comment("Cuts of reconstructable events")};
54  Atom<unsigned int> MinHitX{Name("MinHitX"), Comment("Cuts of reconstructable events")};
55  Atom<unsigned int> MinHitY{Name("MinHitY"), Comment("Cuts of reconstructable events")};
56 
57  Atom<int> MinClust{Name("MinClust"), Comment("Minimum number of clusters to attempt in FuzzyK 2D clustering stage")};
58  Atom<int> MaxClust{Name("MaxClust"), Comment("Maximum number of clusters to attempt in FuzzyK 2D clustering stage")};
59 
60  Atom<double> Memb{Name("Memb"), Comment("Minimum membership in a cluster")};
61 
62  Atom<double> ScrubDist{Name("ScrubDist"), Comment("All cells must have a neighbor within this distance to be used")};
63 
64  // Parameters controling the splitting of clusters
65  Atom<int> VertPlaneGap{Name("VertPlaneGap"), Comment("Maximum plane away from vertex a hit can be to use aggresive splitting")};
66  Atom<double> VertGapSize{Name("VertGapSize"), Comment("Minimum gap size, in cm, to split for hits near vertex")};
67  Atom<int> VertIntPGap{Name("VertIntPGap"), Comment("Minimum number of planes between two hits in a view at the split point in vertex region,"
68  "there can be no hit activity on that track in that view between those hits")};
69  Atom<double> NonVertGapSize{Name("NonVertGapSize"), Comment("Minimum gap size, in cm, to split for hits away from vertex region")};
70  Atom<int> NonVertIntPGap{Name("NonVertIntPGap"), Comment("Minimum number of planes between two hits in a view at the split point away from the vertex region,"
71  "there can be no hit activity on that track in that view between those hits")};
72 
73  Atom<bool> ArtificialVert{Name("ArtificialVert"), Comment("Use an artificial vertex based on MC truth")};
74  Atom<double> SmearX{Name("SmearX"), Comment("Gaussian smear around fake vertex")};
75  Atom<double> SmearY{Name("SmearY"), Comment("Gaussian smear around fake vertex")};
76  Atom<double> SmearZ{Name("SmearZ"), Comment("Gaussian smear around fake vertex")};
77 
80 
81  Atom<std::string> PngInst3D{Name("PngInst3D"), Comment("Instance labels to store results in")};
82  Atom<std::string> PngInst2D{Name("PngInst2D"), Comment("Instance labels to store results in")};
83  Atom<std::string> PngInstOrph{Name("PngInstOrph"), Comment("Instance labels to store results in")};
84  Atom<std::string> PngInst3DWt{Name("PngInst3DWt"), Comment("Instance labels to store results in")};
85  Atom<std::string> PngInst2DWt{Name("PngInst2DWt"), Comment("Instance labels to store results in")};
86 
87  Atom<std::string> PngAssn3D{Name("PngAssn3D"), Comment("Instance labels to store associations in")};
88  Atom<std::string> PngAssn2D{Name("PngAssn2D"), Comment("Instance labels to store associations in")};
89  Atom<std::string> PngAssnOrph{Name("PngAssnOrph"), Comment("Instance labels to store associations in")};
90  Atom<std::string> PngAssn3DWt{Name("PngAssn3DWt"), Comment("Instance labels to store associations in")};
91  Atom<std::string> PngAssn2DWt{Name("PngAssn2DWt"), Comment("Instance labels to store associations in")};
92 
93  Atom<bool> ObeyPreselection{Name("ObeyPreselection"), Comment("If true, a slice will be skipped if preselection did not pass")};
94  Atom<std::string> PreselectionLabels{Name("PreselectionLabels"), Comment("Specific labels for the preselection used")};
95 
96 
100  };
101 
102 
104  {
105  public:
106  // Allows 'nova --print-description' to work
108 
109  explicit FuzzyKVertex(const Parameters& params);
110  virtual ~FuzzyKVertex();
111 
112  void beginJob();
113  void produce(art::Event &e);
114 
115  private:
116 
117  /// \brief Function that runs the FuzzyKMean algorithm and produces 2d clusters for each view
118  /// \param prongs A pointer vector collecting the resulting 2d prongs after clustering
119  /// \param slice The input rb::Cluster object to use in clustering
120  /// \param vert The candidate vertex to cluster around
121  void MakeProngs(std::vector<rb::Prong>* prongs,
122  const rb::Cluster slice,
123  TVector3 vert);
124 
125  /// \brief Function for perform splitting after the 2d clustering stage and return prong candidates
126  /// \param s rb::Cluster object containing hits that are added to the prong candidates
127  /// \param fuzzy The FuzzyKMeanAlg object with the clustering results
128  /// \param cindx The index of the cluster center to turn into a prong
129  /// \param view The view, either kX or kY in which the prong takes place
130  /// \param v The vertex seed used in the clustering
131  /// \param prongs Pointer vector to the collections of prongs the new addition will be added to
132  void AddProng(const rb::Cluster s,
133  FuzzyKMeanAlg* fuzzy,
134  unsigned int cindx,
135  geo::View_t view,
136  TVector3 v,
137  std::vector<rb::Prong>* prongs);
138 
139  /// \brief Function to clean up the slice and remove hits that have no nearby neighbors to reduce noise
140  /// \param c Input is a PtrVector to the collection of cellhits in the slice.
141  /// \return A PtrVector of cellhits after the cleaning process
143 
144  /// \brief Function to find hits that are not contained in any prong
145  /// \param c A PtrVector to the collection of cellhits in the slice
146  /// \param prongs A vector of 2D prongs which have not been view matched
147  /// \return A PtrVector of cellhits from the slice which are not contained in any prong
148  art::PtrVector<rb::CellHit> FindOrphanedHits(art::PtrVector<rb::CellHit> c, std::vector<rb::Prong> prongs);
149 
150  std::set<int> FindVisibleProngs3D(const art::PtrVector<rb::CellHit>& allhits,
151  std::set<int> ids,
152  std::vector<int> *cellsX,
153  std::vector<int> *cellsY);
154  private:
156 
157  /// An object of FuzzyKMeanAlg to do the main clustering
159 
160  /// An object for recalculating cluster centers during the splitting stage
162 
163  /// An object of ViewMatchAlg to do the 3D view matching
165 
166  /// An object of WeightSharedHitsAlg to do the hit sharing
168 
170 
171  /// Tree to output debugging info for view merging process
172  //TTree* fOutTree;
173 
174  /// ntuple variables to debug view matching
175  int run;
176  int subrun;
177  int event;
178  // int xzid;
179  // int yzid;
180  // int slc;
181  // double score;
182  // double shiftd;
183  std::vector<double> xzs;
184  std::vector<double> yzs;
185  std::vector<double> xzpe;
186  std::vector<double> yzpe;
187 
188  //TH1D* png2DE;
189 
190  };
191 }
192 ////////////////////////////////////////////////////////////////////////////////////////////////
193 
194 namespace fuzz{
195 
197  : fSliceToken(consumes<std::vector<rb::Cluster>>(params().SliceLabel())),
198  fFuzzyClustAlg(new FuzzyKMeanAlg(params().FuzzyKMean())),
199  fFuzzyTestAlg(new FuzzyKMeanAlg(params().FuzzyKMean())),
200  fViewMatchAlg(new ViewMatchAlg(params().ViewMatch())),
201  fWeightSharedHitsAlg(new WeightSharedHitsAlg(params().WeightSharedHits())),
202  fParams(params())
203  {
204  this->produces< std::vector<rb::Prong> >(fParams.PngInst3D());
205  this->produces< std::vector<rb::Prong> >(fParams.PngInst2D());
206  this->produces< std::vector<rb::Prong> >(fParams.PngInstOrph());
207  this->produces< std::vector<rb::WeightedProng> >(fParams.PngInst3DWt());
208  this->produces< std::vector<rb::WeightedProng> >(fParams.PngInst2DWt());
209 
210  if(fParams.ArtificialVert()){
211  this->produces<art::Assns<rb::Vertex, rb::Prong> >(fParams.PngAssn3D());
212  this->produces<art::Assns<rb::Vertex, rb::Prong> >(fParams.PngAssn2D());
213  this->produces<art::Assns<rb::Vertex, rb::Prong> >(fParams.PngAssnOrph());
214  this->produces<art::Assns<rb::Vertex, rb::WeightedProng> >(fParams.PngAssn3DWt());
215  this->produces<art::Assns<rb::Vertex, rb::WeightedProng> >(fParams.PngAssn2DWt());
216  }
217  else{
218  this->produces<art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssn3D());
219  this->produces<art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssn2D());
220  this->produces<art::Assns<rb::Prong, rb::Vertex> >(fParams.PngAssnOrph());
221  this->produces<art::Assns<rb::WeightedProng, rb::Vertex> >(fParams.PngAssn3DWt());
222  this->produces<art::Assns<rb::WeightedProng, rb::Vertex> >(fParams.PngAssn2DWt());
223  }
224 
225  this->produces<art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssn3D());
226  this->produces<art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssn2D());
227  this->produces<art::Assns<rb::Prong, rb::Cluster> >(fParams.PngAssnOrph());
228  this->produces<art::Assns<rb::WeightedProng, rb::Cluster> >(fParams.PngAssn3DWt());
229  this->produces<art::Assns<rb::WeightedProng, rb::Cluster> >(fParams.PngAssn2DWt());
230  this->produces<art::Assns<rb::Prong, rb::WeightedProng> >(fParams.PngAssn3DWt());
231  this->produces<art::Assns<rb::Prong, rb::WeightedProng> >(fParams.PngAssn2DWt());
232 
233  //if we want to make use of a fake vertex instead of the real one
234  if(fParams.ArtificialVert()){
235  this->produces<std::vector<rb::Vertex> >("FakeVert");
236  this->produces<art::Assns<rb::Vertex,rb::Cluster> >("FakeVert");
237 
238  // get the random number seed, use a random default if not specified
239  // in the configuration file.
240  unsigned int seed = sim::GetRandomNumberSeed();
241  createEngine( seed );
242  }
243  }
244 
245  //......................................................................
246 
248  {
249  //art::ServiceHandle<art::TFileService> tfs;
250 
251  // define ta tree containing the results of the view matching for diagnostic purposes
252  //fOutTree = tfs->make<TTree>("ViewMatch","View Matching");
253 
254  //png2DE = tfs->make<TH1D>("fpng2DE", "", 20, 0, 1);
255 
256  //store view matching summary information
257  //fOutTree->Branch("xzid", &xzid, "xzid/I");
258  //fOutTree->Branch("yzid", &yzid, "yzid/I");
259  //fOutTree->Branch("slc", &slc, "slc/I");
260  //fOutTree->Branch("score", &score);
261  //fOutTree->Branch("shiftd", &shiftd);
262  //fOutTree->Branch("xzs", &xzs);
263  //fOutTree->Branch("yzs", &yzs);
264  //fOutTree->Branch("xzpe", &xzpe);
265  //fOutTree->Branch("yzpe", &yzpe);
266 
267  }
268 
269  //......................................................................
270 
272  {
273  if(fFuzzyClustAlg) delete fFuzzyClustAlg;
274  if(fFuzzyTestAlg) delete fFuzzyTestAlg;
275  if(fViewMatchAlg) delete fViewMatchAlg;
277  }
278 
279  //......................................................................
280 
282  {
283  unsigned int i;
284 
285  std::unique_ptr< std::vector<rb::Prong> > png3D(new std::vector<rb::Prong>);
286  std::unique_ptr<art::Assns<rb::Prong, rb::Vertex> > png3DAssns(new art::Assns<rb::Prong, rb::Vertex>);
287  std::unique_ptr<art::Assns<rb::Prong, rb::Cluster> > png3DSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
288 
289  std::unique_ptr< std::vector<rb::Prong> > png2D(new std::vector<rb::Prong>);
290  std::unique_ptr<art::Assns<rb::Prong, rb::Vertex> > png2DAssns(new art::Assns<rb::Prong, rb::Vertex>);
291  std::unique_ptr<art::Assns<rb::Prong, rb::Cluster> > png2DSlcAssns(new art::Assns<rb::Prong, rb::Cluster>);
292 
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  std::unique_ptr< std::vector<rb::WeightedProng> > png3DWt(new std::vector<rb::WeightedProng>);
298  std::unique_ptr<art::Assns<rb::WeightedProng, rb::Vertex> > png3DWtAssns(new art::Assns<rb::WeightedProng, rb::Vertex>);
299  std::unique_ptr<art::Assns<rb::WeightedProng, rb::Cluster> > png3DWtSlcAssns(new art::Assns<rb::WeightedProng, rb::Cluster>);
300 
301  std::unique_ptr< std::vector<rb::WeightedProng> > png2DWt(new std::vector<rb::WeightedProng>);
302  std::unique_ptr<art::Assns<rb::WeightedProng, rb::Vertex> > png2DWtAssns(new art::Assns<rb::WeightedProng, rb::Vertex>);
303  std::unique_ptr<art::Assns<rb::WeightedProng, rb::Cluster> > png2DWtSlcAssns(new art::Assns<rb::WeightedProng, rb::Cluster>);
304 
305  std::unique_ptr<art::Assns<rb::Prong, rb::WeightedProng> > png3Dpng3DWtAssns(new art::Assns<rb::Prong, rb::WeightedProng>);
306  std::unique_ptr<art::Assns<rb::Prong, rb::WeightedProng> > png2Dpng2DWtAssns(new art::Assns<rb::Prong, rb::WeightedProng>);
307 
308  //in case we want to use an artificial vertex instead of the real one, useful for testing algorithm and decoupling from ElasticArms
309  std::unique_ptr< std::vector<rb::Vertex> > aVert(new std::vector<rb::Vertex>);
310  std::unique_ptr<art::Assns<rb::Vertex, rb::Cluster> > avAssns(new art::Assns<rb::Vertex, rb::Cluster>);
311  std::unique_ptr<art::Assns<rb::Vertex, rb::Prong> > fpng3DAssns(new art::Assns<rb::Vertex, rb::Prong>);
312  std::unique_ptr<art::Assns<rb::Vertex, rb::Prong> > fpng2DAssns(new art::Assns<rb::Vertex, rb::Prong>);
313  std::unique_ptr<art::Assns<rb::Vertex, rb::Prong> > fpngOrphAssns(new art::Assns<rb::Vertex, rb::Prong>);
314  std::unique_ptr<art::Assns<rb::Vertex, rb::WeightedProng> > fpng3DWtAssns(new art::Assns<rb::Vertex, rb::WeightedProng>);
315  std::unique_ptr<art::Assns<rb::Vertex, rb::WeightedProng> > fpng2DWtAssns(new art::Assns<rb::Vertex, rb::WeightedProng>);
316 
318 
319  // Pull out the time slices and hough results for the event
321  evt.getByToken(fSliceToken, sHandle);
323  for(unsigned int i=0; i<sHandle->size(); ++i){
324  slices.push_back(art::Ptr<rb::Cluster>(sHandle, i));
325  }
326 
327  std::vector<int> bestMCId;
328  std::vector<cheat::NeutrinoWithIndex> truths;
329  std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur;
330  //if we want to use an artificial vertex based off of the true interaction vertex
331  //build a matrix to match the proper mctruth object to a slice
332  if(fParams.ArtificialVert()){
334  //get vector of neutrinos in event
335  truths = bt->allMCTruth();
336  sEffPur = bt->SlicesToMCTruthsTable(slices);
337  //get best neutrino match for each slice
338  bestMCId = bt->SliceToOrderedNuIdsByEff(truths, sEffPur);
339  }
340 
341  //object to get vertex associated with a slice
342  art::FindManyP<rb::Vertex> fmv(sHandle, evt, fParams.VertexLabel());
343 
344  //container to hold prongs
345  std::vector<rb::Prong> prongs;
346 
347  // containers needed for making associations between fake verticies and tracks.
348  // Need to keep track of the slice and the start and stop indicies of the
349  //tracks in that slice out of the entire collection
350  std::vector<size_t> vIndx3D;
351  std::vector<size_t> vIndx2D;
352  std::vector<size_t> vIndxOrph;
353  std::vector<size_t> vIndx3DWt;
354  std::vector<size_t> vIndx2DWt;
355 
356  std::vector<size_t> sIndx3D;
357  std::vector<size_t> eIndx3D;
358  std::vector<size_t> sIndx2D;
359  std::vector<size_t> eIndx2D;
360  std::vector<size_t> sIndxOrph;
361  std::vector<size_t> eIndxOrph;
362  std::vector<size_t> sIndx3DWt;
363  std::vector<size_t> eIndx3DWt;
364  std::vector<size_t> sIndx2DWt;
365  std::vector<size_t> eIndx2DWt;
366 
367  std::vector<size_t> png3DIndx;
368  std::vector<size_t> png2DIndx;
369  std::vector<size_t> png3DWtIndx;
370  std::vector<size_t> png2DWtIndx;
371 
372  //loop over slices to do the work
373  for (i=0; i<sHandle->size(); ++i) {
374 
375  std::vector<rb::WeightedProng> pngWt;
376  //obey any preselection if desired
377  if(fParams.ObeyPreselection() && rb::IsFiltered(evt,sHandle,i,fParams.PreselectionLabels())) continue;
378 
379  art::Ptr<rb::Cluster> sptr(sHandle, i);
380 
381  // Skip the noise cluster and require some minimum number of hits
382  if (sptr->IsNoise()) continue;
383 
384  // perform scrubbing
385  rb::Cluster neighbors(this->Scrub(sptr->AllCells()));
386 
387  //make sure slice is reconstructable before proceeding
388  if (neighbors.NXCell()<fParams.MinHitX()) continue;
389  if (neighbors.NYCell()<fParams.MinHitY()) continue;
390  if (neighbors.NCell() <fParams.MinHit()) continue;
391 
392  // Find any verticies that go with this time slice
393  TVector3 vert;
394  std::vector<art::Ptr<rb::Vertex> > v;
395 
396  //if a fake vertex is needed, make it now
397  if (fParams.ArtificialVert()){
398  if (bestMCId[i] == -1) continue;
399 
400  // get the RandomNumberGenerator service
402  CLHEP::HepRandomEngine &engine = rng->getEngine();
403  CLHEP::RandGaussQ gauss(engine);
404 
405  art::Ptr<simb::MCTruth> nuTru = sEffPur[i][bestMCId[i]].neutrinoInt;
406  // Go ahead and shoot gaussians
407  double Vx = gauss.fire(0.0, fParams.SmearX());
408  double Vy = gauss.fire(0.0, fParams.SmearY());
409  double Vz = gauss.fire(0.0, fParams.SmearZ());
410  double T = 0.0;
411  // Now shift to the original vertex
412  // If neutrino is set, get to the neutrino vertex
413  if(nuTru->NeutrinoSet()){
414  Vx += nuTru->GetNeutrino().Nu().Vx();
415  Vy += nuTru->GetNeutrino().Nu().Vy();
416  Vz += nuTru->GetNeutrino().Nu().Vz();
417  T += nuTru->GetNeutrino().Nu().Position().T();
418  }
419  else{
420  // Otherwise, get to the first particle vertex
421  Vx += nuTru->GetParticle(0).Vx();
422  Vy += nuTru->GetParticle(0).Vy();
423  Vz += nuTru->GetParticle(0).Vz();
424  T += nuTru->GetParticle(0).Position().T();
425  }
426 
427  //make fake vertex and put it in the event
428  rb::Vertex vtx(Vx, Vy, Vz, T);
429  aVert->push_back(vtx);
430  util::CreateAssn(*this, evt, *aVert, slices[i], *avAssns, UINT_MAX, "FakeVert");
431  vert = vtx.GetXYZ();
432  }
433  //else get the vertex associated with the slice
434  else{
435  v = fmv.at(i);
436  //NOTE: Currently only one vertex is produced per slice, if this ever changes
437  //then this logic must also change.
438  if (v.size() != 1) continue;
439  vert = v[0]->GetXYZ();
440  }
441 
442  //make sure all vertex values are real and not nan or inf. Make sure vertex is not too far outside detector(this cut is far too loose right now
443  //but for rock events a vertex outside the detector is possible, don't want to risk being too restrictive without study.
444  if(std::isnan(vert.X()) || std::isinf(vert.X()) || std::abs(vert.X())>10000) continue;
445  if(std::isnan(vert.Y()) || std::isinf(vert.Y()) || std::abs(vert.Y())>10000) continue;
446  if(std::isnan(vert.Z()) || std::isinf(vert.Z()) || std::abs(vert.Z())>10000) continue;
447 
448  //store ntuple information for view match debugging
449  run = evt.run();
450  subrun = evt.subRun();
451  event = evt.id().event();
452 
453  prongs.clear();
454  //make the 2D prong candidates for this slice
455  this->MakeProngs(&prongs, neighbors, vert);
456 
457  //perform the view matching with this set of 2D prongs
458  fViewMatchAlg->LoadProngs(prongs);
460 
461  //record view matching information into the ntuple
462  //for (unsigned int k=0; k<fViewMatchAlg->xzid.size(); ++k){
463  // xzid = -1;
464  // yzid = -1;
465  // score = -1;
466  // shiftd = -11.0;
467  // xzpe.clear();
468  // yzpe.clear();
469  // xzs.clear();
470  // yzs.clear();
471 
472  // xzid = fViewMatchAlg->xzid[k];
473  // yzid = fViewMatchAlg->yzid[k];
474  // score = fViewMatchAlg->ksscore[k];
475  // shiftd = fViewMatchAlg->ksshift[k];
476  // xzs = fViewMatchAlg->xzs[k];
477  // yzs = fViewMatchAlg->yzs[k];
478  // xzpe = fViewMatchAlg->xzpe[k];
479  // yzpe = fViewMatchAlg->yzpe[k];
480 
481  // fOutTree->Fill();
482  //}
483 
484 
485  //loop over 3D matched prongs to push into the event
486  for(unsigned int k=0; k<fViewMatchAlg->fMatchedProngs.size(); ++k){
487 
489 
490  //move the prong start point from the elastic arms vertex to the closest hit
491  TVector3 start;
492  double minZ = png.MinZ();
493  double maxZ = png.MaxZ();
494  if ((png.Dir().Z() > 0.0) &&(minZ < png.Start().Z())) minZ = png.Start().Z();
495  if ((png.Dir().Z() < 0.0) &&(maxZ > png.Start().Z())) maxZ = png.Start().Z();
496 
497  if (png.Dir().Z() < 0.0){
498  double temp = maxZ;
499  maxZ = minZ;
500  minZ = temp;
501  }
502 
503  start[0] = (fViewMatchAlg->fMatchedProngs[k].Dir().X()/fViewMatchAlg->fMatchedProngs[k].Dir().Z())*(minZ-fViewMatchAlg->fMatchedProngs[k].Start().Z())
504  + fViewMatchAlg->fMatchedProngs[k].Start().X();
505  start[1] = (fViewMatchAlg->fMatchedProngs[k].Dir().Y()/fViewMatchAlg->fMatchedProngs[k].Dir().Z())*(minZ-fViewMatchAlg->fMatchedProngs[k].Start().Z())
506  + fViewMatchAlg->fMatchedProngs[k].Start().Y();
507  start[2] = minZ;
508 
509  png.SetStart(start);
510  png3D->push_back(png);
511 
512  if (!fParams.ArtificialVert()) util::CreateAssn(*this, evt, *png3D, v[0], *png3DAssns, UINT_MAX, fParams.PngInst3D());
513  util::CreateAssn(*this, evt, *png3D, sptr, *png3DSlcAssns, UINT_MAX, fParams.PngInst3D());
514 
515  // handle weighted prongs
516  pngWt.push_back(png);
517 
518  } // end loop over matched prongs
519 
520  // track the indices so png:pngWt associations can be made
521  if (fViewMatchAlg->fMatchedProngs.size() > 0){
522  for(unsigned int idx=(*png3D).size()-fViewMatchAlg->fMatchedProngs.size();idx<(*png3D).size();++idx){
523  png3DIndx.push_back(idx);
524  }
525  }
526 
527  //if we are making a fake vertex, calculate the start and ending index for this slice's collection of prongs/tracks
528  //so that associations can be made
529  if (fParams.ArtificialVert() && fViewMatchAlg->fMatchedProngs.size() > 0){
530  vIndx3D.push_back(aVert->size()-1);
531  sIndx3D.push_back((*png3D).size()-fViewMatchAlg->fMatchedProngs.size());
532  eIndx3D.push_back((*png3D).size());
533  }
534 
535 
536  //Loop over the unmatched 2d prongs to also put them into the event
537  for(unsigned int k=0; k<fViewMatchAlg->fUnmatchedProngs.size(); ++k){
539  double xyStart = 0.0;
540  double zStart = 0.0;
541  double minZ = png.MinZ();
542  double maxZ = png.MaxZ();
543  if ((png.Dir().Z() > 0.0) &&(minZ < png.Start().Z())) minZ = png.Start().Z();
544  if ((png.Dir().Z() < 0.0) &&(maxZ > png.Start().Z())) maxZ = png.Start().Z();
545 
546  if (png.Dir().Z() < 0.0){
547  double temp = maxZ;
548  maxZ = minZ;
549  minZ = temp;
550  }
551 
552  TVector3 start;
553 
554  if (fViewMatchAlg->fUnmatchedProngs[k].View() == geo::kX){
555  xyStart = (fViewMatchAlg->fUnmatchedProngs[k].Dir().X()/fViewMatchAlg->fUnmatchedProngs[k].Dir().Z())*(minZ-fViewMatchAlg->fUnmatchedProngs[k].Start().Z())
556  + fViewMatchAlg->fUnmatchedProngs[k].Start().X();
557  zStart = minZ;
558  if ((png.MinPlane() == png.MaxPlane())){
559  if (png.Dir().X() > 0.0) {xyStart = png.MinV(geo::kX);}
560  else if (png.Dir().X() < 0.0) {xyStart = png.MaxV(geo::kX);}
561  }
562  start[0] = xyStart;
563  start[1] = vert.Y();
564  start[2] = zStart;
565  }
566 
567  if (fViewMatchAlg->fUnmatchedProngs[k].View() == geo::kY){
568  xyStart = (fViewMatchAlg->fUnmatchedProngs[k].Dir().Y()/fViewMatchAlg->fUnmatchedProngs[k].Dir().Z())*(minZ-fViewMatchAlg->fUnmatchedProngs[k].Start().Z())
569  + fViewMatchAlg->fUnmatchedProngs[k].Start().Y();
570  zStart = minZ;
571  if ((png.MinPlane() == png.MaxPlane())) {
572  if (png.Dir().Y() > 0.0) {xyStart = png.MinV(geo::kY); }
573  else if (png.Dir().Y() < 0.0) {xyStart = png.MaxV(geo::kY);}
574  }
575  start[0] = vert.X();
576  start[1] = xyStart;
577  start[2] = zStart;
578  }
579 
580  png.SetStart(start);
581  //png2DE->Fill(png.TotalGeV());
582 
583  png2D->push_back(png);
584  if (!fParams.ArtificialVert()) util::CreateAssn(*this, evt, *png2D, v[0], *png2DAssns, UINT_MAX, fParams.PngInst2D());
585  util::CreateAssn(*this, evt, *png2D, sptr, *png2DSlcAssns, UINT_MAX, fParams.PngInst2D());
586 
587  // handle weighted prongs
588  pngWt.push_back(png);
589 
590  } // end loop over unmatched 2d tracks
591 
592  // prepare for png:pngWt assns
593  if (fViewMatchAlg->fUnmatchedProngs.size() > 0){
594  for(unsigned int idx=(*png2D).size()-fViewMatchAlg->fUnmatchedProngs.size();idx<(*png2D).size();++idx){
595  png2DIndx.push_back(idx);
596  }
597  }
598 
599  //if we are making a fake vertex, prepare accordingly for the associations
601  vIndx2D.push_back(aVert->size()-1);
602  sIndx2D.push_back((*png2D).size()-fViewMatchAlg->fUnmatchedProngs.size());
603  eIndx2D.push_back((*png2D).size());
604  }
605 
606  // do hit sharing / weighting for pngWt
608  for(unsigned int k=0; k<pngWt.size(); ++k){
609  if(pngWt[k].Is3D()) {
610  png3DWt->push_back(pngWt[k]);
611  if (!fParams.ArtificialVert()) util::CreateAssn(*this, evt, *png3DWt, v[0], *png3DWtAssns, UINT_MAX, fParams.PngInst3DWt());
612  util::CreateAssn(*this, evt, *png3DWt, sptr, *png3DWtSlcAssns, UINT_MAX, fParams.PngInst3DWt());
613 
614  }
615  else {
616  png2DWt->push_back(pngWt[k]);
617  if (!fParams.ArtificialVert()) util::CreateAssn(*this, evt, *png2DWt, v[0], *png2DWtAssns, UINT_MAX, fParams.PngInst2DWt());
618  util::CreateAssn(*this, evt, *png2DWt, sptr, *png2DWtSlcAssns, UINT_MAX, fParams.PngInst2DWt());
619  }
620  }
621 
622  // prepare for png:pngWt assns
623  if (fViewMatchAlg->fMatchedProngs.size() > 0){
624  for(unsigned int idx=(*png3DWt).size()-fViewMatchAlg->fMatchedProngs.size();idx<(*png3DWt).size();++idx){
625  png3DWtIndx.push_back(idx);
626  }
627  }
628 
629  // prepare for png:pngWt assns
630  if (fViewMatchAlg->fUnmatchedProngs.size() > 0){
631  for(unsigned int idx=(*png2DWt).size()-fViewMatchAlg->fUnmatchedProngs.size();idx<(*png2DWt).size();++idx){
632  png2DWtIndx.push_back(idx);
633  }
634  }
635 
636  //if we are making a fake vertex, prepare accordingly for the associations with the prongs with weighted hits
637  if (fParams.ArtificialVert() && fViewMatchAlg->fMatchedProngs.size() > 0){
638  vIndx3DWt.push_back(aVert->size()-1);
639  sIndx3DWt.push_back((*png3DWt).size()-fViewMatchAlg->fMatchedProngs.size());
640  eIndx3DWt.push_back((*png3DWt).size());
641  }
642 
644  vIndx2DWt.push_back(aVert->size()-1);
645  sIndx2DWt.push_back((*png2DWt).size()-fViewMatchAlg->fUnmatchedProngs.size());
646  eIndx2DWt.push_back((*png2DWt).size());
647  }
648 
649  // handle orphaned hits here
650  // artificially set the orphaned prong start to be the vertex, direction to be entirely z
651  TVector3 start(vert.X(), vert.Y(), vert.Z());
652  TVector3 dir(0.0, 0.0, 1.0);
653 
654  rb::Prong pngO(FindOrphanedHits(sptr->AllCells(), prongs), start, dir);
655  pngOrph->push_back(pngO);
656 
657  if (!fParams.ArtificialVert()) util::CreateAssn(*this, evt, *pngOrph, v[0], *pngOrphAssns, UINT_MAX, fParams.PngInstOrph());
658  util::CreateAssn(*this, evt, *pngOrph, sptr, *pngOrphSlcAssns, UINT_MAX, fParams.PngInstOrph());
659 
660  //if we are making a fake vertex, prepare accordingly for the associations
661  if (fParams.ArtificialVert()){
662  vIndxOrph.push_back(aVert->size()-1);
663  sIndxOrph.push_back((*pngOrph).size()-1);
664  eIndxOrph.push_back((*pngOrph).size());
665  }
666 
667  } // Overall loop on slices
668 
669 
670  // move the prong-to-weighted-prong associations here and model it like the fake vertex associations
671  for(i=0; i<png3DIndx.size(); ++i){
672  util::CreateAssn(*this, evt, *png3Dpng3DWtAssns, *png3D, *png3DWt, png3DWtIndx[i], png3DIndx[i], fParams.PngInst3D(), fParams.PngInst3DWt());
673  }
674 
675  for(i=0; i<png2DIndx.size(); ++i){
676  util::CreateAssn(*this, evt, *png2Dpng2DWtAssns, *png2D, *png2DWt, png2DWtIndx[i], png2DIndx[i], fParams.PngInst2D(), fParams.PngInst2DWt());
677  }
678 
679  //if a fake vertex is produced, now make the associations between the vertex and tracks
680  //can include weighted prongs in each loop, sizes should be the same
681  if (fParams.ArtificialVert()){
682  for(i=0; i<vIndx3D.size(); ++i){
683  util::CreateAssn(*this, evt, *aVert, *png3D, *fpng3DAssns, sIndx3D[i],eIndx3D[i],vIndx3D[i],"FakeVert",fParams.PngInst3D());
684  util::CreateAssn(*this, evt, *aVert, *png3DWt, *fpng3DWtAssns, sIndx3DWt[i],eIndx3DWt[i],vIndx3DWt[i],"FakeVert",fParams.PngInst3DWt());
685  }
686  for(i=0; i<vIndx2D.size(); ++i){
687  util::CreateAssn(*this, evt, *aVert, *png2D, *fpng2DAssns, sIndx2D[i],eIndx2D[i],vIndx2D[i],"FakeVert",fParams.PngInst2D());
688  util::CreateAssn(*this, evt, *aVert, *png2DWt, *fpng2DWtAssns, sIndx2DWt[i],eIndx2DWt[i],vIndx2DWt[i],"FakeVert",fParams.PngInst2DWt());
689  }
690  for(i=0; i<vIndxOrph.size(); ++i){
691  util::CreateAssn(*this, evt, *aVert, *pngOrph, *fpngOrphAssns, sIndxOrph[i],eIndxOrph[i],vIndxOrph[i],"FakeVert",fParams.PngInstOrph());
692  }
693  }
694  //put products into the event
695  evt.put(std::move(png3D),fParams.PngInst3D());
696  evt.put(std::move(png3DWt),fParams.PngInst3DWt());
697 
698  if (fParams.ArtificialVert()){
699  evt.put(std::move(fpng3DAssns),fParams.PngAssn3D());
700  evt.put(std::move(fpng2DAssns),fParams.PngAssn2D());
701  evt.put(std::move(fpngOrphAssns),fParams.PngAssnOrph());
702  evt.put(std::move(fpng3DWtAssns),fParams.PngAssn3DWt());
703  evt.put(std::move(fpng2DWtAssns),fParams.PngAssn2DWt());
704  }
705  else{
706  evt.put(std::move(png3DAssns),fParams.PngAssn3D());
707  evt.put(std::move(png2DAssns),fParams.PngAssn2D());
708  evt.put(std::move(pngOrphAssns),fParams.PngAssnOrph());
709  evt.put(std::move(png3DWtAssns),fParams.PngAssn3DWt());
710  evt.put(std::move(png2DWtAssns),fParams.PngAssn2DWt());
711  }
712 
713  evt.put(std::move(png3DSlcAssns),fParams.PngAssn3D());
714  evt.put(std::move(png2DSlcAssns),fParams.PngAssn2D());
715  evt.put(std::move(pngOrphSlcAssns),fParams.PngAssnOrph());
716  evt.put(std::move(png3DWtSlcAssns),fParams.PngAssn3DWt());
717  evt.put(std::move(png2DWtSlcAssns),fParams.PngAssn2DWt());
718 
719  evt.put(std::move(png2D),fParams.PngInst2D());
720  evt.put(std::move(png2DWt),fParams.PngInst2DWt());
721  evt.put(std::move(pngOrph),fParams.PngInstOrph());
722 
723  evt.put(std::move(png3Dpng3DWtAssns),fParams.PngAssn3DWt());
724  evt.put(std::move(png2Dpng2DWtAssns),fParams.PngAssn2DWt());
725 
726  if (fParams.ArtificialVert()){
727  evt.put(std::move(aVert),"FakeVert");
728  evt.put(std::move(avAssns),"FakeVert");
729  }
730 
731  }
732 
733  //......................................................................
734 
735  /// Function to produce the 2D prongs in both views before the matching
736  void FuzzyKVertex::MakeProngs(std::vector<rb::Prong>* prongs,
737  const rb::Cluster slice,
738  TVector3 vert)
739  {
740 
741  //first cluster in the xz view
742  fFuzzyClustAlg->LoadCluster(slice);
744 
745  //loop from back to front. FuzzyKMean makes clusters with the weakest cluster seed
746  //having the last index. When producing prongs start with the weak cluster so hits can be pushed
747  //into the more dominant prongs before being dropped.
748  for (int i=fFuzzyClustAlg->fA.size()-1; i>-1; i--){
749  this->AddProng(slice, fFuzzyClustAlg, i, geo::kX, vert, prongs);
750  }
751 
752  //now cluster in the yz view
754 
755  //loop from back to front. FuzzyKMean makes clusters with the weakest cluster seed
756  //having the last index. When producing prongs start with the weak cluster so hits can be pushed
757  //into the more dominant prongs before being dropped.
758  for ( int i=fFuzzyClustAlg->fA.size()-1; i>-1; i--){
759  this->AddProng(slice, fFuzzyClustAlg, i, geo::kY, vert, prongs);
760  }
761 
762  }
763 
764  //......................................................................
765 
767  FuzzyKMeanAlg* fuzzy,
768  unsigned int cindx,
769  geo::View_t view,
770  TVector3 v,
771  std::vector<rb::Prong>* prongs)
772  {
774 
775  //initialize variables
776  double theta = fuzzy->fA[cindx];
777  double dCosXY = sin(theta);
778  double dCosZ = cos(theta);
779  double zStart = v.Z();
780  double xyStart = v[view];
781  unsigned int vertplane = 999;
782 
783  //snap vertex into detector boundaries so that it can be determined which view the vertex plane lies in
784  double vx = v.X();
785  double vy = v.Y();
786  double vz = v.Z();
787  if (v.X() < -geom->DetHalfWidth()) vx = -geom->DetHalfWidth()+15.0;
788  if (v.X() > geom->DetHalfWidth()) vx = geom->DetHalfWidth()-15.0;
789  if (v.Y() < -geom->DetHalfHeight()) vy = -geom->DetHalfHeight()+15.0;
790  if (v.Y() > geom->DetHalfHeight()) vy = geom->DetHalfHeight()-15.0;
791  if (v.Z() < 0.0) vz = 2.0;
792  if (v.Z() > geom->DetLength()) vz = geom->DetLength()-2.0;
793 
794  if ((std::abs(v.X()) < geom->DetHalfWidth()) && (std::abs(v.Y()) < geom->DetHalfHeight()) && (v.Z() > 0.0) && (v.Z() < geom->DetLength())){
795  geo::CellUniqueId cid = geom->CellId(vx, vy, vz, 1.5, 1.5, 6., 1.0);
796  if (cid != 0) vertplane = geom->getPlaneID(cid);
797  }
798 
799  //calculate distance of each hit to the vertex and sort
800  unsigned int cellcount = s.NCell(view);
801  std::vector<std::pair<float, unsigned int> > distvec;
802  for (unsigned int i=0; i<cellcount; ++i){
803  if (fuzzy->fU[cindx][i] >= fParams.Memb()){
804  distvec.push_back(std::make_pair(fuzzy->fDist[i], i));
805  }
806  }
807  sort(distvec.begin(), distvec.end());
808 
809  //loop over collection of cells in the cluster, break into smaller groups if needed
810  std::vector<std::vector<std::pair<float, unsigned int> > > cluslist;
811  std::vector<std::pair<float, unsigned int> > cluscand;
812  cluscand.push_back(distvec[0]); //first cell always put in candidate cluster
813  for (unsigned int i=1; i<distvec.size(); ++i){
814 
815  //check if hit is near the vertex. in this special case gap size for splitting is different
816  if (std::abs(1.*s.Cell(view,distvec[i-1].second)->Plane() - 1.*vertplane) <= fParams.VertPlaneGap()) {
817  if ((distvec[i].first - distvec[i-1].first > fParams.VertGapSize()) && (abs(s.Cell(view,distvec[i].second)->Plane() - s.Cell(view,distvec[i-1].second)->Plane()) >= fParams.VertIntPGap())){
818  cluslist.push_back(cluscand);
819  cluscand.clear();
820  cluscand.push_back(distvec[i]);
821  continue;
822  }
823  else {
824  cluscand.push_back(distvec[i]);
825  continue;
826  }
827  }
828  //in non vertex region use different splitting criteria
829  else if ((distvec[i].first - distvec[i-1].first >= fParams.NonVertGapSize()) &&
830  (abs(s.Cell(view,distvec[i].second)->Plane() - s.Cell(view,distvec[i-1].second)->Plane()) >= fParams.NonVertIntPGap())){
831  cluslist.push_back(cluscand);
832  cluscand.clear();
833  cluscand.push_back(distvec[i]);
834  continue;
835  }
836  //otherwise add cluster unchanged
837  else {
838  cluscand.push_back(distvec[i]);
839  continue;
840  }
841  }
842  //add final cluster to list
843  cluslist.push_back(cluscand);
844 
845  //if the list is of size 1, there is no splitting and all cells can be added to one prong
846  if (cluslist.size() == 1){
847  rb::Cluster clus(view);
848  for (unsigned int i=0; i<cluslist[0].size(); ++i){
849  unsigned int hit = cluslist[0][i].second;
850  clus.Add(s.Cell(view,hit));
851  }
852  rb::Prong p(clus, xyStart, zStart, dCosXY, dCosZ);
853  (*prongs).push_back(p);
854  }
855  //if splits occurred we must re cluster accordingly
856  else{
857  for (unsigned int i=0; i<cluslist.size(); ++i){
858  //TODO: There should be some minimum hit/energy requirement
859  //for keeping a cluster so that small things get thrown away, but want
860  //to preserve neutrons. For now no cuts are made on energy or hits.
861 
862  //Check fraction of cells with minimum membership to a different cluster.
863  //This bit is tricky since cells in the same candidate grouping may not
864  //belong to the same different cluster. In current model, if a hit has the
865  //minimum membership in another cluster let it go there.
866  for (unsigned int j=0; j<cluslist[i].size(); ++j){
867  unsigned int hit = cluslist[i][j].second;
868  for(unsigned int k=0; k<fuzzy->fA.size(); ++k){
869  if((k != cindx) && (fuzzy->fU[k][hit] >= fParams.Memb())) {fuzzy->fU[cindx][hit] = 0.0; break;}
870  }
871  }
872 
873  //now put this into a prong. We lost hits somewhere, so make a special FuzzyKMean object to calculate new angle for
874  //for this cluster.
875  rb::Cluster clus(view);
876  for (unsigned int j=0; j<cluslist[i].size(); ++j){
877  unsigned int hit = cluslist[i][j].second;
878  if (fuzzy->fU[cindx][hit] >= fParams.Memb()) clus.Add(s.Cell(view,hit));
879  }
880  if (clus.NCell() > 0){
881  fFuzzyTestAlg->LoadCluster(clus);
882  //we know all the hits belong together, just need to recalculate angle,
883  //so set min and max cluster size to 1
884  fFuzzyTestAlg->DoClustering(1,1, xyStart,zStart,view);
885  theta = fFuzzyTestAlg->fA[0];
886  dCosXY = sin(theta);
887  dCosZ = cos(theta);
888  //make sure the prong is physical before adding it to cluster
889  if (!std::isnan(xyStart) && !std::isnan(zStart) && !std::isnan(dCosXY) && !std::isnan(dCosZ) &&
890  !std::isinf(xyStart) && !std::isinf(zStart) && !std::isinf(dCosXY) && !std::isinf(dCosZ)){
891  rb::Prong p(clus, xyStart, zStart, dCosXY, dCosZ);
892  (*prongs).push_back(p);
893  }
894  }//end if
895 
896  }// end for loop
897  }//end else statement
898 
899  }
900 
901  //......................................................................
902 
903  /// Function to remove hits with no neighbors from the slice
905  {
906 
907  std::vector<int> hit_numNeighbors;
908  //in an ideal detector a hit must have 1 neighbor
909  //inside this distance.
910  double maxGap = fParams.ScrubDist();
911  art::PtrVector<rb::CellHit> trimCells;
912  int iPlane, fPlane, iCell, fCell;
913  double goodCh, totCh;
914 
915  iPlane = 0;
916  iCell = 0;
917  fPlane = 0;
918  fCell = 0;
919 
922 
923  for (unsigned int i=0; i<c.size(); i++){
924  hit_numNeighbors.push_back(0);
925  }
926 
927  //loop to find a cells closest neighbor
928  for (unsigned int i=0; i<c.size(); i++){
929  double minDist = 9999.0;
930  int minCell = -1;
931  if (hit_numNeighbors[i] > 0) continue;
932  for(unsigned int j=0; j<c.size(); ++j){
933  if (i == j) continue;
934  if (c[i]->View() != c[j]->View()) continue;
935  double xyz1[3], xyz2[3];
936  geom->Plane((c[i])->Plane())->Cell((c[i])->Cell())->GetCenter(xyz1);
937  geom->Plane((c[j])->Plane())->Cell((c[j])->Cell())->GetCenter(xyz2);
938  double dist = util::pythag(xyz1[0]-xyz2[0],
939  xyz1[1]-xyz2[1],
940  xyz1[2]-xyz2[2]);
941  if (dist < minDist){
942  minDist = dist;
943  minCell = j;
944  }
945  }
946  if (minCell == -1) continue;
947 
948  //define the rectangle of cells such that a cell and its
949  //nearest neighbor are at diagonal corners
950  if (c[i]->Plane() < c[minCell]->Plane()){
951  iPlane = c[i]->Plane();
952  fPlane = c[minCell]->Plane();
953  }
954  else {
955  iPlane = c[minCell]->Plane();
956  fPlane = c[i]->Plane();
957  }
958  if (c[i]->Cell() < c[minCell]->Cell()) {
959  iCell = c[i]->Cell();
960  fCell = c[minCell]->Cell();
961  }
962  else {
963  iCell = c[minCell]->Cell();
964  fCell = c[i]->Cell();
965  }
966 
967  //within the rectangle find the fraction of cells that are good
968  goodCh = 0;
969  totCh = 0;
970  for (int c=iCell; c<=fCell; c++){
971  for (int p=iPlane; p<=fPlane; p+=2){
972  totCh++;
973  if(!(badc->IsBad(p,c))) ++goodCh;
974  }
975  }
976  if ((minDist*goodCh/totCh) < maxGap){
977  hit_numNeighbors[i]++;
978  hit_numNeighbors[minCell]++;
979  }
980  }
981 
982  //keep the hits that have at least 1 neighbor in the slice
983  for (unsigned int i=0; i<c.size(); i++){
984  if (hit_numNeighbors[i] > 0){
985  trimCells.push_back(c[i]);
986  }
987  }
988 
989  return trimCells;
990  }
991 
992 //......................................................................
993 
994  // function to find orphaned hits
996  {
997 
998  art::PtrVector<rb::CellHit> orphanHits;
999 
1000  // loop over all hits in the slice
1001  for (unsigned int i=0; i<c.size(); i++){
1002  bool inprong = false;
1003 
1004  // loop over all prongs
1005  for (unsigned int j=0; j<prongs.size(); ++j){
1006 
1007  // check to see if this hit is in the prong
1008  for(unsigned int k = 0; k < prongs[j].NCell(); ++k) {
1009  if (*c[i] == *prongs[j].Cell(k)) {
1010  inprong = true;
1011  break; // if hit in prong, break out of loops
1012  }
1013 
1014  }// end hits in prong
1015 
1016  if (inprong) break;
1017  }// end prongs
1018  if (!inprong) { // if this slice hit was in no prongs, add to orphan hit prong
1019  orphanHits.push_back(c[i]);
1020  }
1021  }//end hits in slice
1022 
1023  return orphanHits;
1024  }
1025 
1026 
1027 //......................................................................
1028  struct Counts
1029  {
1031  {
1032  nHits[0] = nHits[1] = 0;
1033  nOnlyHits[0] = nOnlyHits[1] = 0;
1034  }
1035  int nHits[2]; // Indexed by view
1036  int nOnlyHits[2]; // Hit by only one particle
1037  };
1038 
1039 //......................................................................
1041  std::set<int> ids,
1042  std::vector<int> *cellsX,
1043  std::vector<int> *cellsY) {
1044 
1045  std::map<int, Counts> counts;
1047 
1048  // loop over all hits and check if this track id contributed to them
1049  for(unsigned int i = 0; i < allhits.size(); ++i){
1050 
1051  //get vector of track id's that contribute to this hit
1052  const art::Ptr<rb::CellHit> hit = allhits[i];
1053  const std::vector<cheat::TrackIDE> trackIDs = bt->HitToTrackIDE(hit);
1054 
1055  // loop over all the returned track ids
1056  int nTracks = 0; // number of tracks that contributed energy to the hit
1057  int loneID = -1; // Which ID it was if there was only one
1058 
1059  for(unsigned int j = 0; j < trackIDs.size(); ++j){
1060 
1061  const int trackID = trackIDs[j].trackID;
1062  const float eFrac = trackIDs[j].energyFrac;
1063  //detector energy resolution is around 10%
1064  //to be a visible track we are looking for at least 1 hit in each view
1065  //where a true particle deposited 20% (2 sigma) of the energy in that cell
1066  //this is a very minimum definition but leaves the most room for optimization
1067  if(eFrac >= 0.2){
1068  ++nTracks;
1069  loneID = trackID;
1070 
1071  // check if the current id matches the input track id
1072  if(ids.find(trackID) != ids.end()){
1073  ++counts[trackID].nHits[hit->View()];
1074  }
1075  }
1076  } // end for j
1077 
1078  if(nTracks == 1){
1079  ++counts[loneID].nOnlyHits[hit->View()];
1080  }
1081 
1082  } // end for i
1083 
1084  const int kMinHits = 1; // minimum number of hits to be considered visible
1085  const int kMinOnlyHit = 0; // minimum number of hits in which it is only particle contributing energy to be visible
1086  std::set<int> ret;
1087 
1088  for(std::set<int>::iterator it = ids.begin(); it != ids.end(); ++it){
1089 
1090  const Counts c = counts[*it];
1091  if((c.nOnlyHits[geo::kX] >= kMinOnlyHit) && (c.nOnlyHits[geo::kY] >= kMinOnlyHit)){
1092  if((c.nHits[geo::kX] >= kMinHits) && (c.nHits[geo::kY] >= kMinHits)){
1093  //std::cout<<"For track: "<<(*it)<<" xcells: "<<c.nHits[geo::kX]<<" ycells "<<c.nHits[geo::kY]<<std::endl;
1094  cellsX->push_back(c.nHits[geo::kX]);
1095  cellsY->push_back(c.nHits[geo::kY]);
1096  ret.insert(*it);
1097  }
1098  }
1099 
1100  }
1101  return ret;
1102  } // end FindVisibleProngs
1103 
1104 }
1105 ////////////////////////////////////////////////////////////////////////
1106 
1107 namespace fuzz
1108 {
1110 }
1111 ////////////////////////////////////////////////////////////////////////
double MinV(geo::View_t view) const
Definition: Cluster.cxx:454
void LoadProngs(std::vector< rb::Prong > prongs)
Function to load the vector of 2d prongs to be matched.
Atom< std::string > PreselectionLabels
Atom< std::string > PngAssn2DWt
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
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.
Atom< unsigned int > MinHitY
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
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.
set< int >::iterator it
double MinZ() const
Definition: Cluster.h:206
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
void AddProng(const rb::Cluster s, FuzzyKMeanAlg *fuzzy, unsigned int cindx, geo::View_t view, TVector3 v, std::vector< rb::Prong > *prongs)
Function for perform splitting after the 2d clustering stage and return prong candidates.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
int run
Tree to output debugging info for view merging process.
std::vector< double > fDist
Distance from vertex.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
void Matching()
Function to perform the view matching.
std::vector< int > SliceToOrderedNuIdsByEff(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Atom< unsigned int > MinHit
std::vector< NeutrinoWithIndex > allMCTruth() const
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
FuzzyKMeanAlg * fFuzzyTestAlg
An object for recalculating cluster centers during the splitting stage.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
art::PtrVector< rb::CellHit > Scrub(art::PtrVector< rb::CellHit > c)
Function to clean up the slice and remove hits that have no nearby neighbors to reduce noise...
double MaxV(geo::View_t view) const
Definition: Cluster.cxx:500
base_engine_t & createEngine(seed_t seed)
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
void abs(TH1 *hist)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< double > xzs
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
Particle class.
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
std::vector< rb::Prong > fMatchedProngs
Container for all the 3D prongs produced after matching.
Definition: ViewMatchAlg.h:77
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
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
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double dist
Definition: runWimpSim.h:113
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void MakeProngs(std::vector< rb::Prong > *prongs, const rb::Cluster slice, TVector3 vert)
Function that runs the FuzzyKMean algorithm and produces 2d clusters for each view.
const XML_Char * s
Definition: expat.h:262
ViewMatchAlg * fViewMatchAlg
An object of ViewMatchAlg to do the 3D view matching.
std::vector< double > yzpe
base_engine_t & getEngine() const
void beginJob()
std::vector< double > fA
Cluster centers.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
unsigned int seed
Definition: runWimpSim.h:102
std::vector< std::vector< double > > fU
Assign cluster i to sub. j.
void produce(art::Event &e)
Table< fuzz::WeightSharedHitsParams > WeightSharedHits
TVector3 GetXYZ() const
Definition: Vertex.cxx:45
void DoClustering(unsigned int minClust, unsigned int maxClust, double xy, double z, geo::View_t view)
A function to control all steps in clustering in 2D on a group of hits.
Atom< std::string > PngInstOrph
Table< fuzz::ViewMatchParams > ViewMatch
Utility class for making prongs from vertex.
int getPlaneID(const CellUniqueId &id) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
Atom< std::string > PngAssn2D
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
Fuzzy k-Means prong-finding algorithm.
Atom< std::string > SliceLabel
Atom< unsigned int > MinHitX
double DetHalfHeight() const
static constexpr Double_t gauss
Definition: Munits.h:360
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
Definition: View.py:1
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
Atom< std::string > PngInst2D
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
unsigned int GetRandomNumberSeed()
Definition: Simulation.cxx:13
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
double Vx(const int i=0) const
Definition: MCParticle.h:220
std::vector< double > yzs
FuzzyKVertex(const Parameters &params)
double MaxZ() const
Definition: Cluster.h:219
A Cluster with defined start position and direction.
Definition: Prong.h:19
T sin(T number)
Definition: d0nt_math.hpp:132
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Atom< std::string > PngAssn3D
Atom< std::string > PngAssn3DWt
Definition: event.h:1
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
double Vz(const int i=0) const
Definition: MCParticle.h:222
T cos(T number)
Definition: d0nt_math.hpp:78
void GetWeights(std::vector< rb::WeightedProng > &prong)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
FuzzyKMeanAlg * fFuzzyClustAlg
An object of FuzzyKMeanAlg to do the main clustering.
Atom< std::string > PngAssnOrph
std::set< int > FindVisibleProngs3D(const art::PtrVector< rb::CellHit > &allhits, std::set< int > ids, std::vector< int > *cellsX, std::vector< int > *cellsY)
double T
Definition: Xdiff_gwt.C:5
std::vector< rb::Prong > fUnmatchedProngs
Container for all the unmatched 2D prongs left over.
Definition: ViewMatchAlg.h:80
std::vector< double > xzpe
Atom< std::string > PngInst2DWt
bool NeutrinoSet() const
Definition: MCTruth.h:77
Table< fuzz::FuzzyKMeanParams > FuzzyKMean
Atom< std::string > VertexLabel
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
Int_t trackID
Definition: plot.C:84
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
bool IsBad(int plane, int cell)
Atom< std::string > PngInst3DWt
EventID id() const
Definition: Event.h:56
FuzzyKVertexParams fParams
Atom< std::string > PngInst3D
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
void LoadCluster(rb::Cluster slice)
Function to take in a cluster and build the hitlist needed to do the clustering.
WeightSharedHitsAlg * fWeightSharedHitsAlg
An object of WeightSharedHitsAlg to do the hit sharing.