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 
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  fHitNt(0),
129  fVtxDx(0),
130  fVtxDy(0),
131  fVtxDz(0),
132  fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SliceLabel")))
133  {
134  this->reconfigure(p);
135  this->produces< std::vector<rb::Vertex> >();
136  this->produces< std::vector<rb::Prong> >();
137  this->produces<art::Assns<rb::Vertex, rb::Cluster> >();
138  this->produces<art::Assns<rb::Prong, rb::Cluster> >();
139  }
140 
141  //......................................................................
142 
144  {
146 
147  if(fMakeSummaryNt) {
148  fSumNt = f->make<TNtuple>("nt",
149  "nt",
150  "run:subrun:evt:e:nhit:m:sumvia:xf:yf:zf:xmc:ymc:zmc");
151  }
152  if(fMakeHitNt) {
153  fHitNt = f->make<TNtuple>("hnt","hnt","run:subrun:evt:xy:z:view:a:mia:via");
154  }
155  if(fCompareToMC) {
156  fVtxDx = f->make<TH1F>("fVtxDx",
157  ";x(fit) - x(mc) [cm];Events",
158  100,-50.0,50.0);
159  fVtxDy = f->make<TH1F>("fVtxDy",
160  ";y(fit) - y(mc) [cm];Events",
161  100,-50.0,50.0);
162  fVtxDz = f->make<TH1F>("fVtxDz",
163  ";z(fit) - z(mc) [cm];Events",
164  100,-50.0,50.0);
165  }
166  }
167 
168  //......................................................................
169 
171  {
172 #define SET(N) { f##N = p.get<__typeof__(f##N)>(#N); }
173 
174  SET(CompareToMC);
175  SET(MakeSummaryNt);
176  SET(MakeHitNt);
177 
178  SET(MinHit);
179  SET(MinHitX);
180  SET(MinHitY);
181  SET(MaxHoughResults);
182  SET(NhoughMatch);
183  SET(AnnealMaxItr);
184  SET(GridSearchNdir);
185 
186  SET(LambdaV);
187  SET(ZTruncFracLo);
188  SET(ZTruncFracHi);
189  SET(ZTruncDz1);
190  SET(ZTruncDz2);
191  SET(VtxBoxBuffer);
192  SET(HoughPeakAsymCut);
193  SET(GridSearchColinear);
194  SET(GridSearchT);
195  SET(GridSearchLambda);
196  SET(MinimizeLambda);
197  SET(MinimizeT0);
198  SET(MinimizeT1);
199  SET(AnnealFrac);
200  SET(HoughPeakThr);
201 
202  SET(HoughResultLabel);
203  SET(MCLabel);
204 
205  SET(ObeyPreselection);
206  SET(BeamBias);
207 
208  SET(GridSearchFz);
209 #undef SET
210  }
211 
212  //......................................................................
213 
215 
216  //......................................................................
217 
219  double xmax,
220  double ymin,
221  double ymax,
222  double zmin,
223  double zmax,
224  const rb::HoughResult& hrx,
225  const rb::HoughResult& hry,
226  ElasticArms& arms)
227  {
228  unsigned int i, j;
230  xmax+fVtxBoxBuffer,
231  ymin-fVtxBoxBuffer,
232  ymax+fVtxBoxBuffer,
233  zmin-fVtxBoxBuffer,
234  zmax+fVtxBoxBuffer,
237  //
238  // Produce a list of seeds using all intersections of hough lines
239  // and points near hits in the event
240  //
241  std::vector<double> dxv;
242  std::vector<double> dyv;
243  std::vector<double> bxv;
244  std::vector<double> byv;
245  for(i = 0; i < fNhoughMatch; ++i) {
246  double tmpd, tmpb;
247  if(i < hrx.fPeak.size()) {
248  hrx.SlopeIcept(i,&tmpd,&tmpb);
249  dxv.push_back(tmpd);
250  bxv.push_back(tmpb);
251  }
252  if(i < hry.fPeak.size()) {
253  hry.SlopeIcept(i,&tmpd,&tmpb);
254  dyv.push_back(tmpd);
255  byv.push_back(tmpb);
256  }
257  }
258  gs.AddVtxPoints(arms, dxv, bxv, dyv, byv, fGridSearchFz);
259  gs.AddHoughIntersections(hrx, hry, fMaxHoughResults, fNhoughMatch);
260 
261  //
262  // Set the seed directions. Use a standard set plus combinations of
263  // the strongest hough lines.
264  //
266  unsigned int n = fMaxHoughResults;
267  if(hrx.fPeak.size() < n) n = hrx.fPeak.size();
268  if(hry.fPeak.size() < n) n = hry.fPeak.size();
269  double hx, hy;
270  double dx, bx, dy, by;
271  for(i = 0; i < n; ++i) {
272  hrx.SlopeIcept(i,&dx,&bx);
273  hx = hrx.fPeak[i].fH;
274  for(j = 0; j < n; ++j) {
275  hry.SlopeIcept(j,&dy,&by);
276  hy = hry.fPeak[j].fH;
277  //
278  // To beat down the combinatorics, only match lines that have
279  // nearly the same peak height in each view
280  //
281  if(fabs(hx-hy)/(hx+hy) < fHoughPeakAsymCut) {
282  gs.AddDirection(acos( 1/sqrt(dx*dx+dy*dy+1)),atan2(dy,dx));
283  gs.AddDirection(acos(-1/sqrt(dx*dx+dy*dy+1)),atan2(dy,dx));
284  }
285  }
286  }
287 
288  arms.SetT (fGridSearchT);
290  gs.FindBestVtx(arms);
291  }
292 
293  //....................................................................
294  // Find the primary vertex location using true Monte Carlo
295  // information. In principle there can be more than one MC
296  // interaction, so return them as lists of x,y,z positions.
297  //
299  std::vector<double>& mcx,
300  std::vector<double>& mcy,
301  std::vector<double>& mcz)
302 
303  {
304  unsigned int j;
306  evt.getByLabel(fMCLabel, mclist);
307  mcx.resize(mclist->size());
308  mcy.resize(mclist->size());
309  mcz.resize(mclist->size());
310  for (j=0; j<mclist->size(); ++j) {
311  art::Ptr<simb::MCTruth> mc(mclist,j);
312  const simb::MCNeutrino& ne(mc->GetNeutrino());
313  const simb::MCParticle& nu(ne.Nu());
314  mcx[j] = nu.Vx();
315  mcy[j] = nu.Vy();
316  mcz[j] = nu.Vz();
317  }
318  }
319 
320  //......................................................................
321 
323  {
324  unsigned int i, j;
325 
326  std::unique_ptr< std::vector<rb::Vertex> > vtx(new std::vector<rb::Vertex>);
327  std::unique_ptr< std::vector<rb::Prong> > prg(new std::vector<rb::Prong>);
328  std::unique_ptr< art::Assns<rb::Vertex, rb::Cluster> >
330  std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> >
332 
334 
335  //
336  // Pull out the time slices and hough results for the event
337  //
339  evt.getByToken(fSliceToken,slice);
340 
341  // Object to get hough lines associated with a slice
343 
344  for(i = 0; i < slice->size(); ++i) {
345  if(fObeyPreselection && rb::IsFiltered(evt, slice, i)) continue;
346  const rb::Cluster& s = (*slice)[i];
347  if(s.IsNoise()) continue;
348  // Scrub this cluster of isolated hits.
349  // When slicer is fixed this will no longer be necessary
350  rb::Cluster neighbors(this->Scrub(s.AllCells()));
351  //
352  // Skip the noise cluster and require some minimum number of hits
353  //
354  if(neighbors.NXCell() < fMinHitX) continue;
355  if(neighbors.NYCell() < fMinHitY) continue;
356  if(neighbors.NCell() < fMinHit) continue;
357 
358  //
359  // Find any hough results that go with this time slice
360  //
361  const rb::HoughResult* hrx = 0;
362  const rb::HoughResult* hry = 0;
363 
364  std::vector<const rb::HoughResult*> hr = fmhr.at(i);
365  for(j = 0; j < hr.size(); ++j) {
366  if(hr[j]->fView == geo::kX) {
367  hrx = hr[j];
368  }
369  if(hr[j]->fView == geo::kY) {
370  hry = hr[j];
371  }
372  }
373  //
374  // No hough results, no joy.
375  //
376  if(hrx == 0 || hry == 0) continue;
377  if(hrx->fPeak.size() < 1) continue;
378  if(hry->fPeak.size() < 1) continue;
379 
380  //
381  // Decide how many prongs we have
382  //
383  unsigned int nx = 0;
384  for(j = 0; j < hrx->fPeak.size(); ++j) {
385  if(hrx->fPeak[j].fH/hrx->fPeak[j].fThr >= fHoughPeakThr) ++nx;
386  }
387  unsigned int ny = 0;
388  for(j = 0; j < hry->fPeak.size(); ++j) {
389  if(hry->fPeak[j].fH/hry->fPeak[j].fThr >= fHoughPeakThr) ++ny;
390  }
391  unsigned int nprong = 1;
392  if(nx > nprong) nprong = nx;
393  if(ny > nprong) nprong = ny;
394  if(nprong > fMaxHoughResults) nprong = fMaxHoughResults;
395 
396  //
397  // Decide where to truncate the event along z
398  //
399  std::vector<double> zs;
400  for (j = 0; j < neighbors.NCell(); ++j) {
401  const rb::CellHit* h = neighbors.Cell(j).get();
402 
403  unsigned short p = h->Plane();
404  unsigned short c = h->Cell();
405 
406  double xyz[3];
407  double dxyz[3];
409  fGeo->CellInfo(p, c, &view, xyz, dxyz);
410 
411  // To form a vertex position that has all the implicit upstream biased
412  // reversed, just flip all the points around before examining them.
413  zs.push_back(fBeamBias ? xyz[2] : -xyz[2]);
414  }
415  sort(zs.begin(), zs.end());
416 
417  unsigned int iz = (fZTruncFracLo*zs.size())/100;
418  double zlo = zs[iz] - fZTruncDz1;
419  double zhi = zs[iz] + fZTruncDz2;
420 
421  unsigned int ncell = 0, ic = 0;
422  std::vector<double> zstemp;
423  for(j = 0; j < zs.size(); ++j) {
424  if(zs[j] >= zlo && zs[j] <= zhi){
425  ++ncell;
426  zstemp.push_back(zs[j]);
427  }
428  }
429 
430  if(ncell == 0) {
431  zstemp.clear();
432  zlo = zs[0] - fZTruncDz1;
433  zhi = zs[zs.size()-1] + fZTruncDz2;
434  for(j = 0; j < zs.size(); ++j) {
435  if(zs[j] >= zlo && zs[j] <= zhi){
436  ++ncell;
437  zstemp.push_back(zs[j]);
438  }
439  }
440  }
441 
442  //
443  // Construct the arms and fill with cell hits
444  //
445  unsigned int nhitx = 0, nhity = 0;
446  ElasticArms arms(ncell, nprong, fGridSearchLambda, fLambdaV);
447  double minZ = fBeamBias ? fGeo->DetLength() : 0;
448  double maxZ = fBeamBias ? 0 : fGeo->DetLength();
449  double minX = fGeo->DetHalfWidth();
450  double maxX = -fGeo->DetHalfWidth();
451  double minY = fGeo->DetHalfHeight();
452  double maxY = -fGeo->DetHalfHeight();
453  for(j = 0; j < neighbors.NCell(); ++j) {
454  const rb::CellHit* h = neighbors.Cell(j).get();
455 
456  unsigned short p = h->Plane();
457  unsigned short c = h->Cell();
458 
459  double xyz[3];
460  double dxyz[3];
462  fGeo->CellInfo(p, c, &view, xyz, dxyz);
463 
464  if(xyz[2] >= zlo && xyz[2] <= zhi) {
465  if(view == geo::kX) {
466  arms.SetHit(ic++, xyz[2], xyz[0], dxyz[2]/sqrt(12.), view);
467  ++nhitx;
468  // For the purpose of finding the Hough vertex seed, don't flip
469  // anything (ie we have to put it back). This prevents having to edit
470  // any of that code.
471  if(xyz[2] > maxZ) maxZ = fBeamBias ? xyz[2] : -xyz[2];
472  if(xyz[2] < minZ) minZ = fBeamBias ? xyz[2] : -xyz[2];
473  if(xyz[0] > maxX) maxX = xyz[0];
474  if(xyz[0] < minX) minX = xyz[0];
475  }
476  else if(view == geo::kY) {
477  arms.SetHit(ic++, xyz[2], xyz[1], dxyz[2]/sqrt(12.), view);
478  ++nhity;
479  if(xyz[2] > maxZ) maxZ = fBeamBias ? xyz[2] : -xyz[2];
480  if(xyz[2] < minZ) minZ = fBeamBias ? xyz[2] : -xyz[2];
481  if(xyz[1] > maxY) maxY = xyz[1];
482  if(xyz[1] < minY) minY = xyz[1];
483  }
484  else abort();
485  }
486  }
487 
488  //
489  // Need at least one hit in each view to continue
490  //
491  if(nhitx == 0 || nhity == 0) continue;
492 
493  //
494  // Seed the arms before trying to minimize
495  //
496  this->FindSeed(minX, maxX, minY, maxY, minZ, maxZ, *hrx, *hry, arms);
497 
498  if(!fBeamBias){
499  // Flip the resulting seed around to match the flipped points we'll be
500  // using.
501  double vx, vy, vz;
502  arms.GetVertex(vx, vy, vz);
503  arms.SetVertex(vx, vy, -vz);
504  }
505 
506  //
507  // With the arms seeded, now try to find the best fit
508  //
509  //arms.SetVertex(415.0,-603.0,4980.0);
511  Minimizer mn(&arms);
512  mn.SetT0(fMinimizeT0);
513  mn.SetT1(fMinimizeT1);
516  mn.Fit();
517 
518  //
519  // Fill the hit-level debugging ntuple
520  //
521  if(fMakeHitNt) {
522  float x[9];
523  for(j = 0; j < arms.fN; ++j) {
524  x[0] = evt.run();
525  x[1] = evt.subRun();
526  x[2] = evt.id().event();
527  x[3] = arms.fXorY[j];
528  x[4] = arms.fZ[j];
529  x[5] = arms.fView[j];
530  for(unsigned int a = 0; a < arms.fM; ++a) {
531  x[6] = a;
532  x[7] = arms.fMia[j][a];
533  x[8] = arms.fVia[j][a];
534  fHitNt->Fill(x);
535  }
536  }
537  }
538 
539  //
540  // Push the results out
541  //
542  art::Ptr<rb::Cluster> sptr(slice, i);
543  rb::Vertex v(arms.fX0, arms.fY0, fBeamBias ? arms.fZ0 : -arms.fZ0,
544  neighbors.MinTNS());
545  vtx->push_back(v);
546  util::CreateAssn(*this, evt, *vtx, sptr, *vAssns);
547 
548  for(j = 0; j < arms.fdXds.size(); ++j) {
549  rb::Prong p(s,
550  TVector3(arms.fX0, arms.fY0, arms.fZ0),
551  TVector3(arms.fdXds[j], arms.fdYds[j], arms.fdZds[j]));
552  prg->push_back(p);
553  util::CreateAssn(*this, evt, *prg, sptr, *pAssns);
554  }
555 
556  //
557  // If we've been asked, make a comparison to MC truth
558  //
559  std::vector<double> xmc; // MC x vertex locations
560  std::vector<double> ymc; // MC y vertex locations
561  std::vector<double> zmc; // MC z vertex locations
562  if (fCompareToMC) {
563  this->MCVertex(evt, xmc, ymc, zmc);
564  for(j = 0; j < xmc.size(); ++j) {
565  double dx = arms.fX0 - xmc[j];
566  double dy = arms.fY0 - ymc[j];
567  double dz = arms.fZ0 - zmc[j];
568  fVtxDx->Fill(dx);
569  fVtxDy->Fill(dy);
570  fVtxDz->Fill(dz);
571  if(fabs(dx) > 4*4 || fabs(dy) > 4*4 || fabs(dz) > 4*6) {
572  std::cout << "ElasticArmsHS: Blown fit ["
573  << evt.run() << ":"
574  << evt.subRun() << ":"
575  << evt.id().event() << "] dx,dy,dz=("
576  << dx << ","
577  << dy << ","
578  << dz << ")"
579  << std::endl;
580  } // if we have a blown fit
581  } // for each entry in mclist
582  } // if making mc comparisons
583 
584  //
585  // Record a summary entry to the ntuple
586  //
587  if(fMakeSummaryNt) {
588  double sum = 0.0;
589  for(i = 0; i < arms.fN; ++i) {
590  for(j = 0; j < arms.fM; ++j) {
591  sum += arms.fVia[i][j];
592  }
593  }
594  float x[13];
595  x[ 0] = evt.run();
596  x[ 1] = evt.subRun();
597  x[ 2] = evt.id().event();
598  x[ 3] = arms.E();
599  x[ 4] = arms.fN;
600  x[ 5] = arms.fM;
601  x[ 6] = sum;
602  x[ 7] = arms.fX0;
603  x[ 8] = arms.fY0;
604  x[ 9] = arms.fZ0;
605  if (xmc.size()>0) {
606  x[10] = xmc[0];
607  x[11] = ymc[0];
608  x[12] = zmc[0];
609  }
610  else {
611  x[10] = -9e9;
612  x[11] = -9e9;
613  x[12] = -9e9;
614  }
615  fSumNt->Fill(x);
616  } // if we are making a summary ntuple
617 
618  } // Overall loop on slices
619 
620  evt.put(std::move(vtx));
621  evt.put(std::move(prg));
622  evt.put(std::move(vAssns));
623  evt.put(std::move(pAssns));
624  }
625 
626  //......................................................................
627  //function to remove hits with no neighbors from the slice
629  {
630 
631  std::vector<int> hit_numNeighbors;
632  //in an ideal detector a hit must have 1 neighbor
633  //inside this distance.
634  double maxGap = 60.0;
635  art::PtrVector<rb::CellHit> trimCells;
636  int iPlane, fPlane, iCell, fCell;
637  double goodCh, totCh;
638 
639  iPlane = 0;
640  iCell = 0;
641  fPlane = 0;
642  fCell = 0;
643 
646 
647  for(unsigned int i = 0; i < c.size(); i++){
648  hit_numNeighbors.push_back(0);
649  }
650 
651  //loop to find a cells closest neighbor
652  for(unsigned int i = 0; i < c.size(); i++){
653  double minDist = 9999.0;
654  int minCell = -1;
655  if(hit_numNeighbors[i] > 0) continue;
656  for(unsigned int j = 0; j < c.size(); j++){
657  if(i == j) continue;
658  if(c[i]->View() != c[j]->View()) continue;
659  double xyz1[3], xyz2[3];
660  geom->Plane((c[i])->Plane())->Cell((c[i])->Cell())->GetCenter(xyz1);
661  geom->Plane((c[j])->Plane())->Cell((c[j])->Cell())->GetCenter(xyz2);
662  double dist = util::pythag(xyz1[0]-xyz2[0],
663  xyz1[1]-xyz2[1],
664  xyz1[2]-xyz2[2]);
665  if(dist < minDist){
666  minDist = dist;
667  minCell = j;
668  }
669  }
670  if(minCell == -1) continue;
671 
672  //define the rectangle of cells such that a cell and its
673  //nearest neighbor are at diagonal corners
674  if(c[i]->Plane() < c[minCell]->Plane()){
675  iPlane = c[i]->Plane();
676  fPlane = c[minCell]->Plane();
677  }
678  else {
679  iPlane = c[minCell]->Plane();
680  fPlane = c[i]->Plane();
681  }
682  if(c[i]->Cell() < c[minCell]->Cell()) {
683  iCell = c[i]->Cell();
684  fCell = c[minCell]->Cell();
685  }
686  else {
687  iCell = c[minCell]->Cell();
688  fCell = c[i]->Cell();
689  }
690 
691  //within the rectangle find the fraction of cells that are good
692  goodCh = 0;
693  totCh = 0;
694  for(int c = iCell; c <= fCell; c++){
695  for(int p = iPlane; p <= fPlane; p+=2){
696  totCh++;
697  if(!(badc->IsBad(p,c))) ++goodCh;
698  }
699  }
700  if((minDist*goodCh/totCh) < maxGap){
701  hit_numNeighbors[i]++;
702  hit_numNeighbors[minCell]++;
703  }
704  }
705 
706  //keep the hits that have at least 1 neighbor in the slice
707  for(unsigned int i = 0; i < c.size(); i++){
708  if(hit_numNeighbors[i] > 0){
709  trimCells.push_back(c[i]);
710  }
711  }
712 
713  return trimCells;
714  }
715 
716 
717  //......................................................................
718 
720  {
721  if(fCompareToMC) {
722  std::cout << "ElasticArmsHS: dx | rms=" << fVtxDx->GetRMS()
723  << " | outliers="
724  << fVtxDx->Integral(0,30)+fVtxDx->Integral(70,101)
725  << std::endl;
726  std::cout << "ElasticArmsHS: dy | rms=" << fVtxDy->GetRMS()
727  << " | outliers="
728  << fVtxDy->Integral(0,30)+fVtxDy->Integral(70,101)
729  << std::endl;
730  std::cout << "ElasticArmsHS: dz | rms=" << fVtxDz->GetRMS()
731  << " | outliers="
732  << fVtxDz->Integral(0,30)+fVtxDz->Integral(70,101)
733  << std::endl;
734  }
735  }
736 
737 }
738  ////////////////////////////////////////////////////////////////////////
739 namespace earms
740 {
742 }
743 ////////////////////////////////////////////////////////////////////////
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
SubRunNumber_t subRun() const
Definition: Event.h:72
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:76
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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
double fY0
Vertex y location (cm)
Definition: ElasticArms.h:116
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:117
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:118
int fMakeHitNt
Make debugging hit-level ntuple.
double dist
Definition: runWimpSim.h:113
Double_t ymax
Definition: plot.C:25
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
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:98
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:121
void SlopeIcept(unsigned int i, double *m, double *b) const
Slope and intercepts of Hough lines.
Definition: HoughResult.cxx:62
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:441
void FindSeed(double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, const rb::HoughResult &hrx, const rb::HoughResult &hry, ElasticArms &arms)
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
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:97
double fZTruncFracHi
z hit to base the window on (fraction)
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
Definition: View.py:1
unsigned int fGridSearchNdir
Number of minimum bias directions to use.
size_type size() const
Definition: PtrVector.h:308
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:100
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
double fX0
Model parameters.
Definition: ElasticArms.h:115
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.
T const * get() const
Definition: Ptr.h:321
std::vector< double > fdXds
Track dx/ds.
Definition: ElasticArms.h:119
matrix fVia
Strength of association hit i to track a.
Definition: ElasticArms.h:111
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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)
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:99
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.
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
std::vector< double > fdYds
Track dy/ds.
Definition: ElasticArms.h:120
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
Summary from a Hough transform applied to a group of cell hits.
Definition: HoughResult.h:35
Event generator information.
Definition: MCNeutrino.h:18
bool IsBad(int plane, int cell)
double fHoughPeakThr
Threshold to accept a peak as good.
ProductToken< T > consumes(InputTag const &)
void SetT0(double T0)
Definition: Minimizer.cxx:46
EventID id() const
Definition: Event.h:56
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:106
T atan2(T number)
Definition: d0nt_math.hpp:72
int fCompareToMC
Should we make a comparison to truth?