ElasticArmsHS_module.cc
Go to the documentation of this file.
1 ///
2 /// \file ElasticArmsHS.cxx
3 /// \brief Application of the elastic arms fitter to hough results
4 /// \author messier@indiana.edu
5 /// \version $Id: ElasticArmsHS.cxx,v 1.21 2012-10-09 01:15:52 edniner Exp $
6 ///
10 
11 #include <algorithm>
12 
13 #include "TNtuple.h"
14 #include "TH1F.h"
15 
19 #include "art_root_io/TFileService.h"
20 
23 #include "Geometry/Geometry.h"
25 #include "RecoBase/CellHit.h"
26 #include "RecoBase/Vertex.h"
27 #include "RecoBase/Cluster.h"
28 #include "RecoBase/Prong.h"
29 #include "RecoBase/FilterList.h"
30 #include "RecoBase/HoughResult.h"
33 #include "Utilities/AssociationUtil.h"
34 #include "Utilities/func/MathUtil.h"
35 
36 /// Elastic Arms vertex finding algorithm
37 namespace earms {
38  ///
39  /// A stand-alone module for running the elastic arm algorithm
40  ///
42  {
43  public:
44  explicit ElasticArmsHS(fhicl::ParameterSet const &p);
45  virtual ~ElasticArmsHS();
46 
47  void beginJob();
48  void reconfigure(const fhicl::ParameterSet& p);
49  void produce(art::Event &e);
51  void endJob();
52 
53  private:
54  void FindSeed(double xmin,
55  double xmax,
56  double ymin,
57  double ymax,
58  double zmin,
59  double zmax,
60  const rb::HoughResult& hrx,
61  const rb::HoughResult& hry,
62  ElasticArms& arms);
63  void MCVertex(art::Event& evt,
64  std::vector<double>& mcx,
65  std::vector<double>& mcy,
66  std::vector<double>& mcz);
67 
68  private:
69  int fCompareToMC; ///< Should we make a comparison to truth?
70  int fMakeSummaryNt; ///< Make a summary ntuple entry for each event
71  int fMakeHitNt; ///< Make debugging hit-level ntuple
72  unsigned int fMinHit; ///< Slices must have at least this many hits
73  unsigned int fMinHitX; ///< Miniumm hits in x view / slice
74  unsigned int fMinHitY; ///< Minimum hit in y view / slice
75  double fLambdaV; ///< Distance scale vertex to first hit
76  //
77  // Paramters which control the initial grid search
78  //
79  bool fBeamBias; ///< True = assume +z beam bias
80  ///< False = assume -z beam bias
81  double fZTruncFracLo; ///< z hit to base the window on (fraction)
82  double fZTruncFracHi; ///< z hit to base the window on (fraction)
83  double fZTruncDz1; ///< How far up stream to set window (cm)
84  double fZTruncDz2; ///< How far down stream to set window (cm)
85  double fVtxBoxBuffer; ///< Vertex must be this close to hits [cm]
86  std::vector<double> fGridSearchFz; ///< Fractional locations of z vtx seeds
87  double fGridSearchColinear; ///< Colinearity cut
88  unsigned int fGridSearchNdir; ///< Number of minimum bias directions to use
89  double fGridSearchT; ///< Temperature to use for grid search
90  double fGridSearchLambda; ///< Noise penalty for grid search
91  unsigned int fMaxHoughResults; ///< Only use the top hough results
92  unsigned int fNhoughMatch; ///< How many lines to check in the other view
93  double fHoughPeakAsymCut; ///< Which x/y hough results to pair?
94  double fHoughPeakThr; ///< Threshold to accept a peak as good
95 
96 
97  //
98  // Parameters which control the final minimization
99  //
100  double fMinimizeLambda; ///< Noise penalty during minimization
101  double fMinimizeT0; ///< Initial temperature for annealing
102  double fMinimizeT1; ///< Final temperature for annealing
103  double fAnnealFrac; ///< Fractional decrease in temperature
104  unsigned int fAnnealMaxItr; ///< Maximum number of annealing steps
105  //
106  // Diagnostics
107  //
108  TNtuple* fSumNt; ///< Summary ntuple
109  TNtuple* fHitNt; ///< Debugging hit-level ntuple
110  TH1F* fVtxDx; ///< Distribution of x(fit)-x(mc)
111  TH1F* fVtxDy; ///< Distribution of y(fit)-y(mc)
112  TH1F* fVtxDz; ///< Distribution of z(fit)-z(mc)
113 
114  //folder labels
116  std::string fHoughResultLabel; ///< Module label to output HoughResults
117  std::string fMCLabel; ///< Module label for MC
118 
119  bool fObeyPreselection; ///< Check rb::IsFiltered?
120  };
121 }
122 
123 namespace earms{
124 
125  //......................................................................
126 
128  EDProducer(p),
129  fHitNt(0),
130  fVtxDx(0),
131  fVtxDy(0),
132  fVtxDz(0),
133  fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SliceLabel")))
134  {
135  this->reconfigure(p);
136  this->produces< std::vector<rb::Vertex> >();
137  this->produces< std::vector<rb::Prong> >();
138  this->produces<art::Assns<rb::Vertex, rb::Cluster> >();
139  this->produces<art::Assns<rb::Prong, rb::Cluster> >();
140  }
141 
142  //......................................................................
143 
145  {
147 
148  if(fMakeSummaryNt) {
149  fSumNt = f->make<TNtuple>("nt",
150  "nt",
151  "run:subrun:evt:e:nhit:m:sumvia:xf:yf:zf:xmc:ymc:zmc");
152  }
153  if(fMakeHitNt) {
154  fHitNt = f->make<TNtuple>("hnt","hnt","run:subrun:evt:xy:z:view:a:mia:via");
155  }
156  if(fCompareToMC) {
157  fVtxDx = f->make<TH1F>("fVtxDx",
158  ";x(fit) - x(mc) [cm];Events",
159  100,-50.0,50.0);
160  fVtxDy = f->make<TH1F>("fVtxDy",
161  ";y(fit) - y(mc) [cm];Events",
162  100,-50.0,50.0);
163  fVtxDz = f->make<TH1F>("fVtxDz",
164  ";z(fit) - z(mc) [cm];Events",
165  100,-50.0,50.0);
166  }
167  }
168 
169  //......................................................................
170 
172  {
173 #define SET(N) { f##N = p.get<__typeof__(f##N)>(#N); }
174 
175  SET(CompareToMC);
176  SET(MakeSummaryNt);
177  SET(MakeHitNt);
178 
179  SET(MinHit);
180  SET(MinHitX);
181  SET(MinHitY);
182  SET(MaxHoughResults);
183  SET(NhoughMatch);
184  SET(AnnealMaxItr);
185  SET(GridSearchNdir);
186 
187  SET(LambdaV);
188  SET(ZTruncFracLo);
189  SET(ZTruncFracHi);
190  SET(ZTruncDz1);
191  SET(ZTruncDz2);
192  SET(VtxBoxBuffer);
193  SET(HoughPeakAsymCut);
194  SET(GridSearchColinear);
195  SET(GridSearchT);
196  SET(GridSearchLambda);
197  SET(MinimizeLambda);
198  SET(MinimizeT0);
199  SET(MinimizeT1);
200  SET(AnnealFrac);
201  SET(HoughPeakThr);
202 
203  SET(HoughResultLabel);
204  SET(MCLabel);
205 
206  SET(ObeyPreselection);
207  SET(BeamBias);
208 
209  SET(GridSearchFz);
210 #undef SET
211  }
212 
213  //......................................................................
214 
216 
217  //......................................................................
218 
220  double xmax,
221  double ymin,
222  double ymax,
223  double zmin,
224  double zmax,
225  const rb::HoughResult& hrx,
226  const rb::HoughResult& hry,
227  ElasticArms& arms)
228  {
229  unsigned int i, j;
231  xmax+fVtxBoxBuffer,
232  ymin-fVtxBoxBuffer,
233  ymax+fVtxBoxBuffer,
234  zmin-fVtxBoxBuffer,
235  zmax+fVtxBoxBuffer,
238  //
239  // Produce a list of seeds using all intersections of hough lines
240  // and points near hits in the event
241  //
242  std::vector<double> dxv;
243  std::vector<double> dyv;
244  std::vector<double> bxv;
245  std::vector<double> byv;
246  for(i = 0; i < fNhoughMatch; ++i) {
247  double tmpd, tmpb;
248  if(i < hrx.fPeak.size()) {
249  hrx.SlopeIcept(i,&tmpd,&tmpb);
250  dxv.push_back(tmpd);
251  bxv.push_back(tmpb);
252  }
253  if(i < hry.fPeak.size()) {
254  hry.SlopeIcept(i,&tmpd,&tmpb);
255  dyv.push_back(tmpd);
256  byv.push_back(tmpb);
257  }
258  }
259  gs.AddVtxPoints(arms, dxv, bxv, dyv, byv, fGridSearchFz);
260  gs.AddHoughIntersections(hrx, hry, fMaxHoughResults, fNhoughMatch);
261 
262  //
263  // Set the seed directions. Use a standard set plus combinations of
264  // the strongest hough lines.
265  //
267  unsigned int n = fMaxHoughResults;
268  if(hrx.fPeak.size() < n) n = hrx.fPeak.size();
269  if(hry.fPeak.size() < n) n = hry.fPeak.size();
270  double hx, hy;
271  double dx, bx, dy, by;
272  for(i = 0; i < n; ++i) {
273  hrx.SlopeIcept(i,&dx,&bx);
274  hx = hrx.fPeak[i].fH;
275  for(j = 0; j < n; ++j) {
276  hry.SlopeIcept(j,&dy,&by);
277  hy = hry.fPeak[j].fH;
278  //
279  // To beat down the combinatorics, only match lines that have
280  // nearly the same peak height in each view
281  //
282  if(fabs(hx-hy)/(hx+hy) < fHoughPeakAsymCut) {
283  gs.AddDirection(acos( 1/sqrt(dx*dx+dy*dy+1)),atan2(dy,dx));
284  gs.AddDirection(acos(-1/sqrt(dx*dx+dy*dy+1)),atan2(dy,dx));
285  }
286  }
287  }
288 
289  arms.SetT (fGridSearchT);
291  gs.FindBestVtx(arms);
292  }
293 
294  //....................................................................
295  // Find the primary vertex location using true Monte Carlo
296  // information. In principle there can be more than one MC
297  // interaction, so return them as lists of x,y,z positions.
298  //
300  std::vector<double>& mcx,
301  std::vector<double>& mcy,
302  std::vector<double>& mcz)
303 
304  {
305  unsigned int j;
307  evt.getByLabel(fMCLabel, mclist);
308  mcx.resize(mclist->size());
309  mcy.resize(mclist->size());
310  mcz.resize(mclist->size());
311  for (j=0; j<mclist->size(); ++j) {
312  art::Ptr<simb::MCTruth> mc(mclist,j);
313  const simb::MCNeutrino& ne(mc->GetNeutrino());
314  const simb::MCParticle& nu(ne.Nu());
315  mcx[j] = nu.Vx();
316  mcy[j] = nu.Vy();
317  mcz[j] = nu.Vz();
318  }
319  }
320 
321  //......................................................................
322 
324  {
325  unsigned int i, j;
326 
327  std::unique_ptr< std::vector<rb::Vertex> > vtx(new std::vector<rb::Vertex>);
328  std::unique_ptr< std::vector<rb::Prong> > prg(new std::vector<rb::Prong>);
329  std::unique_ptr< art::Assns<rb::Vertex, rb::Cluster> >
331  std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> >
333 
335 
336  //
337  // Pull out the time slices and hough results for the event
338  //
340  evt.getByToken(fSliceToken,slice);
341 
342  // Object to get hough lines associated with a slice
344 
345  for(i = 0; i < slice->size(); ++i) {
346  if(fObeyPreselection && rb::IsFiltered(evt, slice, i)) continue;
347  const rb::Cluster& s = (*slice)[i];
348  if(s.IsNoise()) continue;
349  // Scrub this cluster of isolated hits.
350  // When slicer is fixed this will no longer be necessary
351  rb::Cluster neighbors(this->Scrub(s.AllCells()));
352  //
353  // Skip the noise cluster and require some minimum number of hits
354  //
355  if(neighbors.NXCell() < fMinHitX) continue;
356  if(neighbors.NYCell() < fMinHitY) continue;
357  if(neighbors.NCell() < fMinHit) continue;
358 
359  //
360  // Find any hough results that go with this time slice
361  //
362  const rb::HoughResult* hrx = 0;
363  const rb::HoughResult* hry = 0;
364 
365  std::vector<const rb::HoughResult*> hr = fmhr.at(i);
366  for(j = 0; j < hr.size(); ++j) {
367  if(hr[j]->fView == geo::kX) {
368  hrx = hr[j];
369  }
370  if(hr[j]->fView == geo::kY) {
371  hry = hr[j];
372  }
373  }
374  //
375  // No hough results, no joy.
376  //
377  if(hrx == 0 || hry == 0) continue;
378  if(hrx->fPeak.size() < 1) continue;
379  if(hry->fPeak.size() < 1) continue;
380 
381  //
382  // Decide how many prongs we have
383  //
384  unsigned int nx = 0;
385  for(j = 0; j < hrx->fPeak.size(); ++j) {
386  if(hrx->fPeak[j].fH/hrx->fPeak[j].fThr >= fHoughPeakThr) ++nx;
387  }
388  unsigned int ny = 0;
389  for(j = 0; j < hry->fPeak.size(); ++j) {
390  if(hry->fPeak[j].fH/hry->fPeak[j].fThr >= fHoughPeakThr) ++ny;
391  }
392  unsigned int nprong = 1;
393  if(nx > nprong) nprong = nx;
394  if(ny > nprong) nprong = ny;
395  if(nprong > fMaxHoughResults) nprong = fMaxHoughResults;
396 
397  //
398  // Decide where to truncate the event along z
399  //
400  std::vector<double> zs;
401  for (j = 0; j < neighbors.NCell(); ++j) {
402  const rb::CellHit* h = neighbors.Cell(j).get();
403 
404  unsigned short p = h->Plane();
405  unsigned short c = h->Cell();
406 
407  double xyz[3];
408  double dxyz[3];
410  fGeo->CellInfo(p, c, &view, xyz, dxyz);
411 
412  // To form a vertex position that has all the implicit upstream biased
413  // reversed, just flip all the points around before examining them.
414  zs.push_back(fBeamBias ? xyz[2] : -xyz[2]);
415  }
416  sort(zs.begin(), zs.end());
417 
418  unsigned int iz = (fZTruncFracLo*zs.size())/100;
419  double zlo = zs[iz] - fZTruncDz1;
420  double zhi = zs[iz] + fZTruncDz2;
421 
422  unsigned int ncell = 0, ic = 0;
423  std::vector<double> zstemp;
424  for(j = 0; j < zs.size(); ++j) {
425  if(zs[j] >= zlo && zs[j] <= zhi){
426  ++ncell;
427  zstemp.push_back(zs[j]);
428  }
429  }
430 
431  if(ncell == 0) {
432  zstemp.clear();
433  zlo = zs[0] - fZTruncDz1;
434  zhi = zs[zs.size()-1] + fZTruncDz2;
435  for(j = 0; j < zs.size(); ++j) {
436  if(zs[j] >= zlo && zs[j] <= zhi){
437  ++ncell;
438  zstemp.push_back(zs[j]);
439  }
440  }
441  }
442 
443  //
444  // Construct the arms and fill with cell hits
445  //
446  unsigned int nhitx = 0, nhity = 0;
447  ElasticArms arms(ncell, nprong, fGridSearchLambda, fLambdaV);
448  double minZ = fBeamBias ? fGeo->DetLength() : 0;
449  double maxZ = fBeamBias ? 0 : fGeo->DetLength();
450  double minX = fGeo->DetHalfWidth();
451  double maxX = -fGeo->DetHalfWidth();
452  double minY = fGeo->DetHalfHeight();
453  double maxY = -fGeo->DetHalfHeight();
454  for(j = 0; j < neighbors.NCell(); ++j) {
455  const rb::CellHit* h = neighbors.Cell(j).get();
456 
457  unsigned short p = h->Plane();
458  unsigned short c = h->Cell();
459 
460  double xyz[3];
461  double dxyz[3];
463  fGeo->CellInfo(p, c, &view, xyz, dxyz);
464 
465  if(xyz[2] >= zlo && xyz[2] <= zhi) {
466  if(view == geo::kX) {
467  arms.SetHit(ic++, xyz[2], xyz[0], dxyz[2]/sqrt(12.), view);
468  ++nhitx;
469  // For the purpose of finding the Hough vertex seed, don't flip
470  // anything (ie we have to put it back). This prevents having to edit
471  // any of that code.
472  if(xyz[2] > maxZ) maxZ = fBeamBias ? xyz[2] : -xyz[2];
473  if(xyz[2] < minZ) minZ = fBeamBias ? xyz[2] : -xyz[2];
474  if(xyz[0] > maxX) maxX = xyz[0];
475  if(xyz[0] < minX) minX = xyz[0];
476  }
477  else if(view == geo::kY) {
478  arms.SetHit(ic++, xyz[2], xyz[1], dxyz[2]/sqrt(12.), view);
479  ++nhity;
480  if(xyz[2] > maxZ) maxZ = fBeamBias ? xyz[2] : -xyz[2];
481  if(xyz[2] < minZ) minZ = fBeamBias ? xyz[2] : -xyz[2];
482  if(xyz[1] > maxY) maxY = xyz[1];
483  if(xyz[1] < minY) minY = xyz[1];
484  }
485  else abort();
486  }
487  }
488 
489  //
490  // Need at least one hit in each view to continue
491  //
492  if(nhitx == 0 || nhity == 0) continue;
493 
494  //
495  // Seed the arms before trying to minimize
496  //
497  this->FindSeed(minX, maxX, minY, maxY, minZ, maxZ, *hrx, *hry, arms);
498 
499  if(!fBeamBias){
500  // Flip the resulting seed around to match the flipped points we'll be
501  // using.
502  double vx, vy, vz;
503  arms.GetVertex(vx, vy, vz);
504  arms.SetVertex(vx, vy, -vz);
505  }
506 
507  //
508  // With the arms seeded, now try to find the best fit
509  //
510  //arms.SetVertex(415.0,-603.0,4980.0);
512  Minimizer mn(&arms);
513  mn.SetT0(fMinimizeT0);
514  mn.SetT1(fMinimizeT1);
517  mn.Fit();
518 
519  //
520  // Fill the hit-level debugging ntuple
521  //
522  if(fMakeHitNt) {
523  float x[9];
524  for(j = 0; j < arms.fN; ++j) {
525  x[0] = evt.run();
526  x[1] = evt.subRun();
527  x[2] = evt.id().event();
528  x[3] = arms.fXorY[j];
529  x[4] = arms.fZ[j];
530  x[5] = arms.fView[j];
531  for(unsigned int a = 0; a < arms.fM; ++a) {
532  x[6] = a;
533  x[7] = arms.fMia[j][a];
534  x[8] = arms.fVia[j][a];
535  fHitNt->Fill(x);
536  }
537  }
538  }
539 
540  //
541  // Push the results out
542  //
543  art::Ptr<rb::Cluster> sptr(slice, i);
544  rb::Vertex v(arms.fX0, arms.fY0, fBeamBias ? arms.fZ0 : -arms.fZ0,
545  neighbors.MinTNS());
546  vtx->push_back(v);
547  util::CreateAssn(evt, *vtx, sptr, *vAssns);
548 
549  for(j = 0; j < arms.fdXds.size(); ++j) {
550  rb::Prong p(s,
551  TVector3(arms.fX0, arms.fY0, arms.fZ0),
552  TVector3(arms.fdXds[j], arms.fdYds[j], arms.fdZds[j]));
553  prg->push_back(p);
554  util::CreateAssn(evt, *prg, sptr, *pAssns);
555  }
556 
557  //
558  // If we've been asked, make a comparison to MC truth
559  //
560  std::vector<double> xmc; // MC x vertex locations
561  std::vector<double> ymc; // MC y vertex locations
562  std::vector<double> zmc; // MC z vertex locations
563  if (fCompareToMC) {
564  this->MCVertex(evt, xmc, ymc, zmc);
565  for(j = 0; j < xmc.size(); ++j) {
566  double dx = arms.fX0 - xmc[j];
567  double dy = arms.fY0 - ymc[j];
568  double dz = arms.fZ0 - zmc[j];
569  fVtxDx->Fill(dx);
570  fVtxDy->Fill(dy);
571  fVtxDz->Fill(dz);
572  if(fabs(dx) > 4*4 || fabs(dy) > 4*4 || fabs(dz) > 4*6) {
573  std::cout << "ElasticArmsHS: Blown fit ["
574  << evt.run() << ":"
575  << evt.subRun() << ":"
576  << evt.id().event() << "] dx,dy,dz=("
577  << dx << ","
578  << dy << ","
579  << dz << ")"
580  << std::endl;
581  } // if we have a blown fit
582  } // for each entry in mclist
583  } // if making mc comparisons
584 
585  //
586  // Record a summary entry to the ntuple
587  //
588  if(fMakeSummaryNt) {
589  double sum = 0.0;
590  for(i = 0; i < arms.fN; ++i) {
591  for(j = 0; j < arms.fM; ++j) {
592  sum += arms.fVia[i][j];
593  }
594  }
595  float x[13];
596  x[ 0] = evt.run();
597  x[ 1] = evt.subRun();
598  x[ 2] = evt.id().event();
599  x[ 3] = arms.E();
600  x[ 4] = arms.fN;
601  x[ 5] = arms.fM;
602  x[ 6] = sum;
603  x[ 7] = arms.fX0;
604  x[ 8] = arms.fY0;
605  x[ 9] = arms.fZ0;
606  if (xmc.size()>0) {
607  x[10] = xmc[0];
608  x[11] = ymc[0];
609  x[12] = zmc[0];
610  }
611  else {
612  x[10] = -9e9;
613  x[11] = -9e9;
614  x[12] = -9e9;
615  }
616  fSumNt->Fill(x);
617  } // if we are making a summary ntuple
618 
619  } // Overall loop on slices
620 
621  evt.put(std::move(vtx));
622  evt.put(std::move(prg));
623  evt.put(std::move(vAssns));
624  evt.put(std::move(pAssns));
625  }
626 
627  //......................................................................
628  //function to remove hits with no neighbors from the slice
630  {
631 
632  std::vector<int> hit_numNeighbors;
633  //in an ideal detector a hit must have 1 neighbor
634  //inside this distance.
635  double maxGap = 60.0;
636  art::PtrVector<rb::CellHit> trimCells;
637  int iPlane, fPlane, iCell, fCell;
638  double goodCh, totCh;
639 
640  iPlane = 0;
641  iCell = 0;
642  fPlane = 0;
643  fCell = 0;
644 
647 
648  for(unsigned int i = 0; i < c.size(); i++){
649  hit_numNeighbors.push_back(0);
650  }
651 
652  //loop to find a cells closest neighbor
653  for(unsigned int i = 0; i < c.size(); i++){
654  double minDist = 9999.0;
655  int minCell = -1;
656  if(hit_numNeighbors[i] > 0) continue;
657  for(unsigned int j = 0; j < c.size(); j++){
658  if(i == j) continue;
659  if(c[i]->View() != c[j]->View()) continue;
660  double xyz1[3], xyz2[3];
661  geom->Plane((c[i])->Plane())->Cell((c[i])->Cell())->GetCenter(xyz1);
662  geom->Plane((c[j])->Plane())->Cell((c[j])->Cell())->GetCenter(xyz2);
663  double dist = util::pythag(xyz1[0]-xyz2[0],
664  xyz1[1]-xyz2[1],
665  xyz1[2]-xyz2[2]);
666  if(dist < minDist){
667  minDist = dist;
668  minCell = j;
669  }
670  }
671  if(minCell == -1) continue;
672 
673  //define the rectangle of cells such that a cell and its
674  //nearest neighbor are at diagonal corners
675  if(c[i]->Plane() < c[minCell]->Plane()){
676  iPlane = c[i]->Plane();
677  fPlane = c[minCell]->Plane();
678  }
679  else {
680  iPlane = c[minCell]->Plane();
681  fPlane = c[i]->Plane();
682  }
683  if(c[i]->Cell() < c[minCell]->Cell()) {
684  iCell = c[i]->Cell();
685  fCell = c[minCell]->Cell();
686  }
687  else {
688  iCell = c[minCell]->Cell();
689  fCell = c[i]->Cell();
690  }
691 
692  //within the rectangle find the fraction of cells that are good
693  goodCh = 0;
694  totCh = 0;
695  for(int c = iCell; c <= fCell; c++){
696  for(int p = iPlane; p <= fPlane; p+=2){
697  totCh++;
698  if(!(badc->IsBad(p,c))) ++goodCh;
699  }
700  }
701  if((minDist*goodCh/totCh) < maxGap){
702  hit_numNeighbors[i]++;
703  hit_numNeighbors[minCell]++;
704  }
705  }
706 
707  //keep the hits that have at least 1 neighbor in the slice
708  for(unsigned int i = 0; i < c.size(); i++){
709  if(hit_numNeighbors[i] > 0){
710  trimCells.push_back(c[i]);
711  }
712  }
713 
714  return trimCells;
715  }
716 
717 
718  //......................................................................
719 
721  {
722  if(fCompareToMC) {
723  std::cout << "ElasticArmsHS: dx | rms=" << fVtxDx->GetRMS()
724  << " | outliers="
725  << fVtxDx->Integral(0,30)+fVtxDx->Integral(70,101)
726  << std::endl;
727  std::cout << "ElasticArmsHS: dy | rms=" << fVtxDy->GetRMS()
728  << " | outliers="
729  << fVtxDy->Integral(0,30)+fVtxDy->Integral(70,101)
730  << std::endl;
731  std::cout << "ElasticArmsHS: dz | rms=" << fVtxDz->GetRMS()
732  << " | outliers="
733  << fVtxDz->Integral(0,30)+fVtxDz->Integral(70,101)
734  << std::endl;
735  }
736  }
737 
738 }
739  ////////////////////////////////////////////////////////////////////////
740 namespace earms
741 {
743 }
744 ////////////////////////////////////////////////////////////////////////
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
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.
void MCVertex(art::Event &evt, std::vector< double > &mcx, std::vector< double > &mcy, std::vector< double > &mcz)
void AddDirection(double theta, double phi)
Definition: GridSearch.cxx:131
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:77
double minY
Definition: plot_hist.C:9
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TNtuple * fSumNt
Summary ntuple.
double fGridSearchColinear
Colinearity cut.
std::map< std::string, double > xmax
bool getByToken(ProductToken< PROD > const &, Handle< PROD > &result) const
Definition: DataViewImpl.h:462
void SetAnnealMaxItr(unsigned int i)
Definition: Minimizer.cxx:58
double fMinimizeT1
Final temperature for annealing.
double fGridSearchT
Temperature to use for grid search.
double fMinimizeLambda
Noise penalty during minimization.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
void SetStandardDirections()
Definition: GridSearch.cxx:225
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Definition: HoughResult.h:61
void SetLambda(double lam)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
double fY0
Vertex y location (cm)
Definition: ElasticArms.h:127
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double fZTruncFracLo
z hit to base the window on (fraction)
double fVtxBoxBuffer
Vertex must be this close to hits [cm].
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
double fGridSearchLambda
Noise penalty for grid search.
double maxY
Definition: plot_hist.C:10
double fAnnealFrac
Fractional decrease in temperature.
const PlaneGeo * Plane(unsigned int i) const
double fZ0
Vertex z location (cm)
Definition: ElasticArms.h:128
DEFINE_ART_MODULE(TestTMapFile)
double fZTruncDz2
How far down stream to set window (cm)
T acos(T number)
Definition: d0nt_math.hpp:54
void AddVtxPoints(const ElasticArms &arms, const std::vector< double > &mx, const std::vector< double > &bx, const std::vector< double > &my, const std::vector< double > &by, const std::vector< double > &f)
Definition: GridSearch.cxx:52
unsigned int fMinHit
Slices must have at least this many hits.
std::string fHoughResultLabel
Module label to output HoughResults.
#define SET(N)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Elastic Arms vertex finding algorithm.
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
void AddHoughIntersections(const rb::HoughResult &hrx, const rb::HoughResult &hry, unsigned int nmx, unsigned int nmh, bool flipped=false)
Definition: GridSearch.cxx:139
unsigned int fM
Number of arms.
Definition: ElasticArms.h:129
int fMakeHitNt
Make debugging hit-level ntuple.
double dist
Definition: runWimpSim.h:113
Double_t ymax
Definition: plot.C:25
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
double MinTNS() const
Definition: Cluster.cxx:482
std::vector< double > fXorY
Transverse hit positions (cm)
Definition: ElasticArms.h:109
void SetHit(unsigned int i, double z, double xory, double sigma, int view)
Definition: ElasticArms.cxx:75
std::vector< double > fdZds
Track dz/ds.
Definition: ElasticArms.h:132
void SlopeIcept(unsigned int i, double *m, double *b) const
Slope and intercepts of Hough lines.
Definition: HoughResult.cxx:62
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
double dy[NP][NC]
double dx[NP][NC]
double fHoughPeakAsymCut
Which x/y hough results to pair?
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
void FindSeed(double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, const rb::HoughResult &hrx, const rb::HoughResult &hry, ElasticArms &arms)
std::void_t< T > n
const double a
unsigned int fMaxHoughResults
Only use the top hough results.
double fMinimizeT0
Initial temperature for annealing.
unsigned int fNhoughMatch
How many lines to check in the other view.
int evt
void SetT(double t)
Definition: ElasticArms.cxx:90
void GetVertex(double &x, double &y, double &z) const
SubRunNumber_t subRun() const
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
double fLambdaV
Distance scale vertex to first hit.
double DetHalfHeight() const
std::vector< double > fGridSearchFz
Fractional locations of z vtx seeds.
unsigned int fN
Data to fit.
Definition: ElasticArms.h:108
double fZTruncFracHi
z hit to base the window on (fraction)
RunNumber_t run() const
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
unsigned int fGridSearchNdir
Number of minimum bias directions to use.
size_type size() const
Definition: PtrVector.h:302
unsigned int fMinHitX
Miniumm hits in x view / slice.
double FindBestVtx(ElasticArms &arms)
Definition: GridSearch.cxx:296
void SetAnnealFrac(double f)
Definition: Minimizer.cxx:54
OStream cout
Definition: OStream.cxx:6
std::vector< double > fZ
Z locations of hits (cm)
Definition: ElasticArms.h:111
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
TH1F * fVtxDz
Distribution of z(fit)-z(mc)
std::string fMCLabel
Module label for MC.
EventNumber_t event() const
Definition: EventID.h:116
void SetT1(double T1)
Definition: Minimizer.cxx:50
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
double Vx(const int i=0) const
Definition: MCParticle.h:220
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
double fX0
Model parameters.
Definition: ElasticArms.h:126
A Cluster with defined start position and direction.
Definition: Prong.h:19
double DetHalfWidth() const
void produce(art::Event &e)
Data resulting from a Hough transform on the cell hit positions.
std::vector< double > fdXds
Track dx/ds.
Definition: ElasticArms.h:130
matrix fVia
Strength of association hit i to track a.
Definition: ElasticArms.h:122
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
bool fObeyPreselection
Check rb::IsFiltered?
void geom(int which=0)
Definition: geom.C:163
void reconfigure(const fhicl::ParameterSet &p)
TH1F * fVtxDy
Distribution of y(fit)-y(mc)
ElasticArmsHS(fhicl::ParameterSet const &p)
ProductToken< T > consumes(InputTag const &)
Definition: ModuleBase.h:55
double fZTruncDz1
How far up stream to set window (cm)
Double_t ymin
Definition: plot.C:24
art::PtrVector< rb::CellHit > Scrub(art::PtrVector< rb::CellHit > c)
std::vector< int > fView
Which view is the hit in?
Definition: ElasticArms.h:110
TH1F * fVtxDx
Distribution of x(fit)-x(mc)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
unsigned int fMinHitY
Minimum hit in y view / slice.
int fMakeSummaryNt
Make a summary ntuple entry for each event.
TNtuple * fHitNt
Debugging hit-level ntuple.
std::vector< double > fdYds
Track dy/ds.
Definition: ElasticArms.h:131
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
T const * get() const
Definition: Ptr.h:149
Summary from a Hough transform applied to a group of cell hits.
Definition: HoughResult.h:35
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
Event generator information.
Definition: MCNeutrino.h:18
bool IsBad(int plane, int cell)
double fHoughPeakThr
Threshold to accept a peak as good.
void SetT0(double T0)
Definition: Minimizer.cxx:46
EventID id() const
void SetVertex(double x, double y, double z)
Encapsulate the geometry of one entire detector (near, far, ndos)
unsigned int fAnnealMaxItr
Maximum number of annealing steps.
matrix fMia
Distance^2 hit i to track a (sigma^2)
Definition: ElasticArms.h:117
T atan2(T number)
Definition: d0nt_math.hpp:72
int fCompareToMC
Should we make a comparison to truth?
enum BeamMode string