LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
SpacePointAlg.cxx
Go to the documentation of this file.
1 
11 #include <cassert>
12 #include <iostream>
13 #include <sstream>
14 #include <cmath>
15 #include <map>
16 #include <algorithm>
17 #include "cetlib_except/exception.h"
18 #include "messagefacility/MessageLogger/MessageLogger.h"
19 #include "larreco/RecoAlg/SpacePointAlg.h"
20 #include "larcore/Geometry/Geometry.h"
21 #include "larcorealg/Geometry/CryostatGeo.h"
22 #include "larcorealg/Geometry/TPCGeo.h"
23 #include "larcorealg/Geometry/PlaneGeo.h"
24 #include "larcorealg/Geometry/WireGeo.h"
25 #include "canvas/Persistency/Common/PtrVector.h"
26 #include "art/Framework/Services/Optional/TFileService.h"
27 #include "art/Framework/Principal/View.h"
28 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
29 //\todo Remove include of BackTrackerService.h once this algorithm is stripped of test for MC
30 #include "larsim/MCCheater/BackTrackerService.h"
31 #include "lardataobj/RecoBase/Hit.h"
32 #include "lardataobj/RecoBase/SpacePoint.h"
33 #include "lardata/RecoObjects/KHitTrack.h"
34 #include "lardata/RecoObjects/KHitWireX.h"
35 
36 #include "TH1F.h"
37 
38 //----------------------------------------------------------------------
39 // Constructor.
40 //
41 namespace trkf{
42 
43  SpacePointAlg::SpacePointAlg(const fhicl::ParameterSet& pset) :
44  fMaxDT(0.),
45  fMaxS(0.),
46  fMinViews(1000),
47  fEnableU(false),
48  fEnableV(false),
49  fEnableW(false),
50  fFilter(false),
51  fMerge(false),
52  fPreferColl(false),
53  fTickOffsetU(0.),
54  fTickOffsetV(0.),
55  fTickOffsetW(0.)
56  {
57  reconfigure(pset);
58  }
59 
60  //----------------------------------------------------------------------
61  // Destructor.
62  //
63  SpacePointAlg::~SpacePointAlg()
64  {
65  }
66 
67  //----------------------------------------------------------------------
68  // Update configuration parameters.
69  //
70  void SpacePointAlg::reconfigure(const fhicl::ParameterSet& pset)
71  {
72  // Get configuration parameters.
73 
74  fMaxDT = pset.get<double>("MaxDT");
75  fMaxS = pset.get<double>("MaxS");
76 
77  fMinViews = pset.get<int>("MinViews");
78 
79  fEnableU = pset.get<bool>("EnableU");
80  fEnableV = pset.get<bool>("EnableV");
81  fEnableW = pset.get<bool>("EnableW");
82  fFilter = pset.get<bool>("Filter");
83  fMerge = pset.get<bool>("Merge");
84  fPreferColl = pset.get<bool>("PreferColl");
85  fTickOffsetU = pset.get<double>("TickOffsetU", 0.);
86  fTickOffsetV = pset.get<double>("TickOffsetV", 0.);
87  fTickOffsetW = pset.get<double>("TickOffsetW", 0.);
88 
89  // Only allow one of fFilter and fMerge to be true.
90 
91  if(fFilter && fMerge)
92  throw cet::exception("SpacePointAlg") << "Filter and Merge flags are both true.\n";
93 
94  // Report.
95 
96  std::cout << "SpacePointAlg configured with the following parameters:\n"
97  << " MaxDT = " << fMaxDT << "\n"
98  << " MaxS = " << fMaxS << "\n"
99  << " MinViews = " << fMinViews << "\n"
100  << " EnableU = " << fEnableU << "\n"
101  << " EnableV = " << fEnableV << "\n"
102  << " EnableW = " << fEnableW << "\n"
103  << " Filter = " << fFilter << "\n"
104  << " Merge = " << fMerge << "\n"
105  << " PreferColl = " << fPreferColl << "\n"
106  << " TickOffsetU = " << fTickOffsetU << "\n"
107  << " TickOffsetV = " << fTickOffsetV << "\n"
108  << " TickOffsetW = " << fTickOffsetW << std::endl;
109  }
110 
111  //----------------------------------------------------------------------
112  // Print geometry and properties constants.
113  //
114  void SpacePointAlg::update() const
115  {
116  // Generate info report on first call only.
117 
118  static bool first = true;
119  bool report = first;
120  first = false;
121 
122  // Get services.
123 
124  art::ServiceHandle<geo::Geometry> geom;
125  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
126 
127  // Calculate and print geometry information.
128 
129  if(report)
130  mf::LogInfo("SpacePointAlg") << "Updating geometry constants.\n";
131 
132 
133  for(unsigned int cstat = 0; cstat < geom->Ncryostats(); ++cstat){
134 
135  // Loop over TPCs.
136 
137  unsigned int const ntpc = geom->Cryostat(cstat).NTPC();
138 
139  for(unsigned int tpc = 0; tpc < ntpc; ++tpc) {
140  const geo::TPCGeo& tpcgeom = geom->Cryostat(cstat).TPC(tpc);
141 
142  // Loop over planes.
143 
144  unsigned int const nplane = tpcgeom.Nplanes();
145 
146  for(unsigned int plane = 0; plane < nplane; ++plane) {
147  geo::PlaneID planeid(cstat, tpc, plane);
148  const geo::PlaneGeo& pgeom = tpcgeom.Plane(planeid);
149 
150  // Fill view-dependent quantities.
151 
152  geo::View_t view = pgeom.View();
153  std::string viewname = "?";
154  if(view == geo::kU) {
155  viewname = "U";
156  }
157  else if(view == geo::kV) {
158  viewname = "V";
159  }
160  else if(view == geo::kZ) {
161  viewname = "Z";
162  }
163  else
164  throw cet::exception("SpacePointAlg") << "Bad view = "
165  << view << "\n";
166 
167  std::string sigtypename = "?";
168  geo::SigType_t sigtype = geom->SignalType(planeid);
169  if(sigtype == geo::kInduction)
170  sigtypename = "Induction";
171  else if(sigtype == geo::kCollection)
172  sigtypename = "Collection";
173  else
174  throw cet::exception("SpacePointAlg") << "Bad signal type = "
175  << sigtype << "\n";
176 
177  std::string orientname = "?";
178  geo::Orient_t orient = pgeom.Orientation();
179  if(orient == geo::kVertical)
180  orientname = "Vertical";
181  else if(orient == geo::kHorizontal)
182  orientname = "Horizontal";
183  else
184  throw cet::exception("SpacePointAlg") << "Bad orientation = "
185  << orient << "\n";
186 
187  if(report) {
188  const double* xyz = tpcgeom.PlaneLocation(plane);
189  mf::LogInfo("SpacePointAlg") << "\nCryostat, TPC, Plane: "
190  << cstat << ","
191  << tpc << ", "
192  << plane << "\n"
193  << " View: " << viewname << "\n"
194  << " SignalType: " << sigtypename << "\n"
195  << " Orientation: " << orientname << "\n"
196  << " Plane location: " << xyz[0] << "\n"
197  << " Plane pitch: "
198  << tpcgeom.Plane0Pitch(plane) << "\n"
199  << " Wire angle: "
200  << tpcgeom.Plane(plane).Wire(0).ThetaZ() << "\n"
201  << " Wire pitch: " << tpcgeom.WirePitch() << "\n"
202  << " Time offset: "
203  << detprop->GetXTicksOffset(plane,tpc,cstat) << "\n";
204  }
205 
206  if(orient != geo::kVertical)
207  throw cet::exception("SpacePointAlg")
208  << "Horizontal wire geometry not implemented.\n";
209  }// end loop over planes
210  }// end loop over tpcs
211  }// end loop over cryostats
212  }
213 
214 
215 
216  //----------------------------------------------------------------------
217  // Get corrected time for the specified hit.
218  double SpacePointAlg::correctedTime(const recob::Hit& hit) const
219  {
220  // Get services.
221 
222  art::ServiceHandle<geo::Geometry> geom;
223  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
224 
225  // Correct time for trigger offset and plane-dependent time offsets.
226 
227  double t = hit.PeakTime() - detprop->GetXTicksOffset(hit.WireID().Plane,hit.WireID().TPC,hit.WireID().Cryostat);
228  if(hit.View() == geo::kU)
229  t -= fTickOffsetU;
230  else if(hit.View() == geo::kV)
231  t -= fTickOffsetV;
232  else if(hit.View() == geo::kW)
233  t -= fTickOffsetW;
234 
235  return t;
236  }
237 
238  //----------------------------------------------------------------------
239  // Spatial separation of hits (zero if two or fewer).
240  double SpacePointAlg::separation(const art::PtrVector<recob::Hit>& hits) const
241  {
242  // Get geometry service.
243 
244  art::ServiceHandle<geo::Geometry> geom;
245 
246  // Trivial case - fewer than three hits.
247 
248  if(hits.size() < 3)
249  return 0.;
250 
251  // Error case - more than three hits.
252 
253  if(hits.size() > 3) {
254  mf::LogError("SpacePointAlg") << "Method separation called with more than three htis.";
255  return 0.;
256  }
257 
258  // Got exactly three hits.
259 
260  // Calculate angles and distance of each hit from origin.
261 
262  double dist[3] = {0., 0., 0.};
263  double sinth[3] = {0., 0., 0.};
264  double costh[3] = {0., 0., 0.};
265  unsigned int cstats[3];
266  unsigned int tpcs[3];
267  unsigned int planes[3];
268 
269  for(int i=0; i<3; ++i) {
270 
271  // Get tpc, plane, wire.
272 
273  const recob::Hit& hit = *(hits[i]);
274  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
275  cstats[i] = hit.WireID().Cryostat;
276  tpcs[i] = hit.WireID().TPC;
277  planes[i] = hit.WireID().Plane;
278 
279  // Check tpc and plane errors.
280 
281  for(int j=0; j<i; ++j) {
282 
283  if(cstats[j] != hit.WireID().Cryostat) {
284  mf::LogError("SpacePointAlg") << "Method separation called with hits from multiple cryostats..";
285  return 0.;
286  }
287 
288  if(tpcs[j] != hit.WireID().TPC) {
289  mf::LogError("SpacePointAlg") << "Method separation called with hits from multiple tpcs..";
290  return 0.;
291  }
292 
293  if(planes[j] == hit.WireID().Plane ) {
294  mf::LogError("SpacePointAlg") << "Method separation called with hits from the same plane..";
295  return 0.;
296  }
297  }
298 
299  // Get angles and distance of wire.
300 
301  double hl = wgeom.HalfL();
302  double xyz[3];
303  double xyz1[3];
304  wgeom.GetCenter(xyz);
305  wgeom.GetCenter(xyz1, hl);
306  double s = (xyz1[1] - xyz[1]) / hl;
307  double c = (xyz1[2] - xyz[2]) / hl;
308  sinth[hit.WireID().Plane] = s;
309  costh[hit.WireID().Plane] = c;
310  dist[hit.WireID().Plane] = xyz[2] * s - xyz[1] * c;
311  }
312 
313  double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0]
314  +(sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1]
315  +(sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
316  return S;
317  }
318 
319  //----------------------------------------------------------------------
320  // Check hits for compatibility.
321  // Check hits pairwise for different views and maximum time difference.
322  // Check three hits for spatial compatibility.
323  bool SpacePointAlg::compatible(const art::PtrVector<recob::Hit>& hits,
324  bool useMC) const
325  {
326  // Get services.
327 
328  art::ServiceHandle<geo::Geometry> geom;
329  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
330 
331  int nhits = hits.size();
332 
333  // Fewer than two or more than three hits can never be compatible.
334 
335  bool result = nhits >= 2 && nhits <= 3;
336  bool mc_ok = true;
337  unsigned int tpc = 0;
338  unsigned int cstat = 0;
339 
340  if(result) {
341 
342  // First do pairwise tests.
343  // Do double loop over hits.
344 
345  for(int ihit1 = 0; result && ihit1 < nhits-1; ++ihit1) {
346  const recob::Hit& hit1 = *(hits[ihit1]);
347  geo::WireID hit1WireID = hit1.WireID();
348  geo::View_t view1 = hit1.View();
349 
350  double t1 = hit1.PeakTime() - detprop->GetXTicksOffset(hit1WireID.Plane,hit1WireID.TPC,hit1WireID.Cryostat);
351 
352  // If using mc information, get a collection of track ids for hit 1.
353  // If not using mc information, this section of code will trigger the
354  // insertion of a single invalid HitMCInfo object into fHitMCMap.
355 
356  const HitMCInfo& mcinfo1 = fHitMCMap[(useMC ? &hit1 : 0)];
357  const std::vector<int>& tid1 = mcinfo1.trackIDs;
358  bool only_neg1 = tid1.size() > 0 && tid1.back() < 0;
359 
360  // Loop over second hit.
361 
362  for(int ihit2 = ihit1+1; result && ihit2 < nhits; ++ihit2) {
363  const recob::Hit& hit2 = *(hits[ihit2]);
364  geo::WireID hit2WireID = hit2.WireID();
365  geo::View_t view2 = hit2.View();
366 
367  // Test for same tpc and different views.
368 
369  result = result && hit1WireID.TPC == hit2WireID.TPC && view1 != view2 && hit1WireID.Cryostat == hit2WireID.Cryostat;
370  if(result) {
371 
372  // Remember which tpc and cryostat we are in.
373 
374  tpc = hit1WireID.TPC;
375  cstat = hit1WireID.Cryostat;
376 
377  double t2 = hit2.PeakTime() - detprop->GetXTicksOffset(hit2WireID.Plane,hit2WireID.TPC,hit2WireID.Cryostat);
378 
379  // Test maximum time difference.
380 
381  result = result && std::abs(t1-t2) <= fMaxDT;
382 
383  // Test mc truth.
384 
385  if(result && useMC) {
386 
387  // Test whether hits have a common parent track id.
388 
389  const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
390  std::vector<int> tid2 = mcinfo2.trackIDs;
391  bool only_neg2 = tid2.size() > 0 && tid2.back() < 0;
392  std::vector<int>::iterator it =
393  std::set_intersection(tid1.begin(), tid1.end(),
394  tid2.begin(), tid2.end(),
395  tid2.begin());
396  tid2.resize(it - tid2.begin());
397 
398  // Hits are compatible if they have parents in common.
399  // If the only parent id in common is negative (-999),
400  // then hits are compatible only if both hits have only
401  // negative parent tracks.
402 
403  bool only_neg3 = tid2.size() > 0 && tid2.back() < 0;
404  mc_ok = tid2.size() > 0 && (!only_neg3 || (only_neg1 && only_neg2));
405  result = result && mc_ok;
406 
407  // If we are still OK, check that either hit is
408  // the nearest neighbor of the other.
409 
410  if(result) {
411  result = mcinfo1.pchit[hit2WireID.Plane] == &hit2 ||
412  mcinfo2.pchit[hit1WireID.Plane] == &hit1;
413  }
414  }
415  }
416  }
417  }
418 
419  // If there are exactly three hits, and they pass pairwise tests, check
420  // for spatial compatibility.
421 
422  if(result && nhits == 3) {
423 
424  // Loop over hits.
425 
426  double dist[3] = {0., 0., 0.};
427  double sinth[3] = {0., 0., 0.};
428  double costh[3] = {0., 0., 0.};
429 
430  for(int i=0; i<3; ++i) {
431 
432  // Get tpc, plane, wire.
433 
434  const recob::Hit& hit = *(hits[i]);
435  geo::WireID hitWireID = hit.WireID();
436 
437  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
438  if ((hitWireID.TPC != tpc) || (hitWireID.Cryostat != cstat))
439  throw cet::exception("SpacePointAlg") << "compatible(): geometry mismatch\n";
440 
441  // Get angles and distance of wire.
442 
443  double hl = wgeom.HalfL();
444  double xyz[3];
445  double xyz1[3];
446  wgeom.GetCenter(xyz);
447  wgeom.GetCenter(xyz1, hl);
448  double s = (xyz1[1] - xyz[1]) / hl;
449  double c = (xyz1[2] - xyz[2]) / hl;
450  sinth[hit.WireID().Plane] = s;
451  costh[hit.WireID().Plane] = c;
452  dist[hit.WireID().Plane] = xyz[2] * s - xyz[1] * c;
453  }
454 
455  // Do space cut.
456 
457  double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0]
458  +(sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1]
459  +(sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
460 
461  result = result && std::abs(S) < fMaxS;
462  }
463  }
464 
465  // Done.
466 
467  return result;
468  }
469 
470  //----------------------------------------------------------------------
471  // Fill one space point using a colleciton of hits.
472  // Assume points have already been tested for compatibility.
473  //
474  void SpacePointAlg::fillSpacePoint(const art::PtrVector<recob::Hit>& hits,
475  std::vector<recob::SpacePoint> &sptv,
476  int sptid) const
477  {
478  // Get services.
479 
480  art::ServiceHandle<geo::Geometry> geom;
481  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
482 
483  double timePitch=detprop->GetXTicksCoefficient();
484 
485  int nhits = hits.size();
486 
487  // Remember associated hits internally.
488 
489  if (fSptHitMap.find(sptid) != fSptHitMap.end())
490  throw cet::exception("SpacePointAlg") << "fillSpacePoint(): hit already present!\n";
491  fSptHitMap[sptid] = hits;
492 
493  // Calculate position and error matrix.
494 
495  double xyz[3] = {0., 0., 0.};
496  double errxyz[6] = {0.,
497  0., 0.,
498  0., 0., 0.};
499 
500  // Calculate x using drift times.
501  // Loop over all hits and calculate the weighted average drift time.
502  // Also calculate time variance and chisquare.
503 
504  double sumt2w = 0.;
505  double sumtw = 0.;
506  double sumw = 0.;
507 
508 
509  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
510  ihit != hits.end(); ++ihit) {
511 
512  const recob::Hit& hit = **ihit;
513  geo::WireID hitWireID = hit.WireID();
514 
515  // Correct time for trigger offset and view-dependent time offsets.
516 
517  double t0 = detprop->GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
518  double t = hit.PeakTime() - t0;
519  double et = hit.SigmaPeakTime();
520  double w = 1./(et*et);
521 
522  sumt2w += w*t*t;
523  sumtw += w*t;
524  sumw += w;
525  }
526 
527  double drift_time = 0.;
528  double var_time = 0.;
529  double chisq = 0.;
530  if(sumw != 0.) {
531  drift_time = sumtw / sumw;
532  //var_time = sumt2w / sumw - drift_time * drift_time;
533  var_time = 1. / sumw;
534  if(var_time < 0.)
535  var_time = 0.;
536  chisq = sumt2w - sumtw * drift_time;
537  if(chisq < 0.)
538  chisq = 0.;
539  }
540  xyz[0] = drift_time * timePitch;
541  errxyz[0] = var_time * timePitch * timePitch;
542 
543  // Calculate y, z using wires (need at least two hits).
544 
545  if(nhits >= 2) {
546 
547  // Calculate y and z by chisquare minimization of wire coordinates.
548 
549  double sw = 0.; // sum w_i
550  double sus = 0.; // sum w_i u_i sin_th_i
551  double suc = 0.; // sum w_i u_i cos_th_i
552  double sc2 = 0.; // sum w_i cos2_th_i
553  double ss2 = 0.; // sum w_i sin2_th_i
554  double ssc = 0.; // sum w_i sin_th_i cos_th_i
555 
556  // Loop over points.
557 
558  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
559  ihit != hits.end(); ++ihit) {
560 
561  const recob::Hit& hit = **ihit;
562  geo::WireID hitWireID = hit.WireID();
563  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
564 
565  // Calculate angle and wire coordinate in this view.
566 
567  double hl = wgeom.HalfL();
568  double cen[3];
569  double cen1[3];
570  wgeom.GetCenter(cen);
571  wgeom.GetCenter(cen1, hl);
572  double s = (cen1[1] - cen[1]) / hl;
573  double c = (cen1[2] - cen[2]) / hl;
574  double u = cen[2] * s - cen[1] * c;
575  double eu = geom->WirePitch(hitWireID.Plane, hitWireID.TPC) / std::sqrt(12.);
576  double w = 1. / (eu * eu);
577 
578  // Summations
579 
580  sw += w;
581  sus += w*u*s;
582  suc += w*u*c;
583  sc2 += w*c*c;
584  ss2 += w*s*s;
585  ssc += w*s*c;
586  }
587 
588  // Calculate y,z
589 
590  double denom = sc2 * ss2 - ssc * ssc;
591  if(denom != 0.) {
592  xyz[1] = (-suc * ss2 + sus * ssc) / denom;
593  xyz[2] = (sus * sc2 - suc * ssc) / denom;
594  errxyz[2] = ss2 / denom;
595  errxyz[4] = ssc / denom;
596  errxyz[5] = sc2 / denom;
597  }
598 
599  // Make space point.
600 
601  recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
602  sptv.push_back(spt);
603  }
604  return;
605  }
606 
625  void SpacePointAlg::fillSpacePoints(std::vector<recob::SpacePoint>& spts,
626  std::multimap<double, KHitTrack> const& trackMap) const
627  {
628  // Loop over KHitTracks.
629 
630  art::PtrVector<recob::Hit> hits;
631  art::PtrVector<recob::Hit> compatible_hits;
632  for(std::multimap<double, KHitTrack>::const_iterator it = trackMap.begin();
633  it != trackMap.end(); ++it) {
634  const KHitTrack& track = (*it).second;
635 
636  // Extrack Hit from track.
637 
638  const std::shared_ptr<const KHitBase>& hit = track.getHit();
639  const KHitWireX* phit = dynamic_cast<const KHitWireX*>(&*hit);
640  if(phit != 0) {
641  const art::Ptr<recob::Hit> prhit = phit->getHit();
642 
643  // Test this hit for compatibility.
644 
645  hits.push_back(prhit);
646  bool ok = this->compatible(hits);
647  if(!ok) {
648 
649  // The new hit is not compatible. Make a space point out of
650  // the last known compatible hits, provided there are at least
651  // two.
652 
653  if(compatible_hits.size() >= 2) {
654  this->fillSpacePoint(compatible_hits, spts, this->numHitMap());
655  compatible_hits.clear();
656  }
657 
658  // Forget about any previous hits.
659 
660  hits.clear();
661  hits.push_back(prhit);
662  }
663 
664  // Update the list of known compatible hits.
665 
666  compatible_hits = hits;
667  }
668  }
669 
670  // Maybe make one final space point.
671 
672  if(compatible_hits.size() >= 2) {
673  this->fillSpacePoint(compatible_hits, spts, this->numHitMap());
674  }
675  }
676 
677  //----------------------------------------------------------------------
678  // Fill one space point using a colleciton of hits.
679  // Assume points have already been tested for compatibility.
680  // This version assumes there can be multiple hits per view,
681  // and gives unequal weight to different hits.
682  //
683  void SpacePointAlg::
684  fillComplexSpacePoint(const art::PtrVector<recob::Hit>& hits,
685  std::vector<recob::SpacePoint>& sptv,
686  int sptid) const
687  {
688  // Get services.
689 
690  art::ServiceHandle<geo::Geometry> geom;
691  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
692 
693  // Calculate time pitch.
694 
695  double timePitch = detprop->GetXTicksCoefficient(); // cm / tick
696 
697  // Figure out which tpc we are in.
698 
699  unsigned int tpc0 = 0;
700  unsigned int cstat0 = 0;
701  int nhits = hits.size();
702  if(nhits > 0) {
703  tpc0 = hits.front()->WireID().TPC;
704  cstat0=hits.front()->WireID().Cryostat;
705  }
706 
707  // Remember associated hits internally.
708 
709  if(fSptHitMap.count(sptid) != 0)
710  throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): hit already present!\n";
711  fSptHitMap[sptid] = hits;
712 
713  // Do a preliminary scan of hits.
714  // Determine weight given to hits in each view.
715 
716  unsigned int nplanes = geom->Cryostat(cstat0).TPC(tpc0).Nplanes();
717  std::vector<int> numhits(nplanes, 0);
718  std::vector<double> weight(nplanes, 0.);
719 
720  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
721  ihit != hits.end(); ++ihit) {
722 
723  const recob::Hit& hit = **ihit;
724  geo::WireID hitWireID = hit.WireID();
725  /* // kept as assertions for performance reasons
726  if (hitWireID.Cryostat != cstat0)
727  throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): incompatible cryostat\n";
728  if (hitWireID.TPC != tpc0);
729  throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): incompatible TPC\n";
730  if (hitWireID.Plane >= nplanes);
731  throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): unknown plane\n";
732  */
733  assert(hitWireID.Cryostat == cstat0);
734  assert(hitWireID.TPC == tpc0);
735  assert(hitWireID.Plane < nplanes);
736  ++numhits[hitWireID.Plane];
737  }
738 
739  for(unsigned int plane = 0; plane < nplanes; ++plane) {
740  double np = numhits[plane];
741  if(np > 0.)
742  weight[plane] = 1. / (np*np*np);
743  }
744 
745  // Calculate position and error matrix.
746 
747  double xyz[3] = {0., 0., 0.};
748  double errxyz[6] = {0.,
749  0., 0.,
750  0., 0., 0.};
751 
752  // Calculate x using drift times.
753  // Loop over all hits and calculate the weighted average drift time.
754  // Also calculate time variance and chisquare.
755 
756  double sumt2w = 0.;
757  double sumtw = 0.;
758  double sumw = 0.;
759 
760  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
761  ihit != hits.end(); ++ihit) {
762 
763  const recob::Hit& hit = **ihit;
764  geo::WireID hitWireID = hit.WireID();
765 
766  // Correct time for trigger offset and view-dependent time offsets.
767 
768  double t0 = detprop->GetXTicksOffset(hitWireID.Plane,hitWireID.TPC,hitWireID.Cryostat);
769  double t = hit.PeakTime() - t0;
770  double et = hit.SigmaPeakTime();
771  double w = weight[hitWireID.Plane]/(et*et);
772 
773  sumt2w += w*t*t;
774  sumtw += w*t;
775  sumw += w;
776  }
777 
778  double drift_time = 0.;
779  double var_time = 0.;
780  double chisq = 0.;
781  if(sumw != 0.) {
782  drift_time = sumtw / sumw;
783  var_time = sumt2w / sumw - drift_time * drift_time;
784  if(var_time < 0.)
785  var_time = 0.;
786  chisq = sumt2w - sumtw * drift_time;
787  if(chisq < 0.)
788  chisq = 0.;
789  }
790  xyz[0] = drift_time * timePitch;
791  errxyz[0] = var_time * timePitch * timePitch;
792 
793  // Calculate y, z using wires (need at least two hits).
794 
795  if(nhits >= 2) {
796 
797  // Calculate y and z by chisquare minimization of wire coordinates.
798 
799  double sw = 0.; // sum w_i
800  double sus = 0.; // sum w_i u_i sin_th_i
801  double suc = 0.; // sum w_i u_i cos_th_i
802  double sc2 = 0.; // sum w_i cos2_th_i
803  double ss2 = 0.; // sum w_i sin2_th_i
804  double ssc = 0.; // sum w_i sin_th_i cos_th_i
805 
806  // Loop over points.
807 
808  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
809  ihit != hits.end(); ++ihit) {
810 
811  const recob::Hit& hit = **ihit;
812  geo::WireID hitWireID = hit.WireID();
813  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
814 
815  // Calculate angle and wire coordinate in this view.
816 
817  double hl = wgeom.HalfL();
818  double cen[3];
819  double cen1[3];
820  wgeom.GetCenter(cen);
821  wgeom.GetCenter(cen1, hl);
822  double s = (cen1[1] - cen[1]) / hl;
823  double c = (cen1[2] - cen[2]) / hl;
824  double u = cen[2] * s - cen[1] * c;
825  double eu = geom->WirePitch(hitWireID.Plane, hitWireID.TPC) / std::sqrt(12.);
826  double w = weight[hitWireID.Plane] / (eu * eu);
827 
828  // Summations
829 
830  sw += w;
831  sus += w*u*s;
832  suc += w*u*c;
833  sc2 += w*c*c;
834  ss2 += w*s*s;
835  ssc += w*s*c;
836  }
837 
838  // Calculate y,z
839 
840  double denom = sc2 * ss2 - ssc * ssc;
841  if(denom != 0.) {
842  xyz[1] = (-suc * ss2 + sus * ssc) / denom;
843  xyz[2] = (sus * sc2 - suc * ssc) / denom;
844  errxyz[2] = ss2 / denom;
845  errxyz[4] = ssc / denom;
846  errxyz[5] = sc2 / denom;
847  }
848 
849  // Make space point.
850 
851  recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
852  sptv.push_back(spt);
853  }
854  return;
855  }
856 
857  //----------------------------------------------------------------------
858  // Fill a vector of space points for all compatible combinations of hits
859  // from an input vector of hits (non-mc-truth version).
860  //
861  void SpacePointAlg::makeSpacePoints(const art::PtrVector<recob::Hit>& hits,
862  std::vector<recob::SpacePoint>& spts) const
863  {
864  makeSpacePoints(hits, spts, false);
865  }
866 
867  //----------------------------------------------------------------------
868  // Fill a vector of space points for all compatible combinations of hits
869  // from an input vector of hits (mc-truth version).
870  //
871  void SpacePointAlg::makeMCTruthSpacePoints(const art::PtrVector<recob::Hit>& hits,
872  std::vector<recob::SpacePoint>& spts) const
873  {
874  makeSpacePoints(hits, spts, true);
875  }
876 
877  //----------------------------------------------------------------------
878  // Fill a vector of space points for all compatible combinations of hits
879  // from an input vector of hits (general version).
880  //
881  void SpacePointAlg::makeSpacePoints(const art::PtrVector<recob::Hit>& hits,
882  std::vector<recob::SpacePoint>& spts,
883  bool useMC) const
884  {
885  // Get services.
886 
887  art::ServiceHandle<geo::Geometry> geom;
888  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
889 
890  // Clear space point to hit map.
891 
892  fSptHitMap.clear();
893 
894  // Print diagnostic information.
895 
896  update();
897 
898  // First make sure result vector is empty.
899 
900  spts.erase(spts.begin(), spts.end());
901 
902  // Statistics.
903 
904  int n2 = 0; // Number of two-hit space points.
905  int n3 = 0; // Number of three-hit space points.
906  int n2filt = 0; // Number of two-hit space points after filtering/merging.
907  int n3filt = 0; // Number of three-hit space pointe after filtering/merging.
908 
909  // Sort hits into maps indexed by [cryostat][tpc][plane][wire].
910  // If using mc information, also generate maps of sim::IDEs and mc
911  // position indexed by hit.
912 
913  std::vector< std::vector<std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit> > > > > hitmap;
914  fHitMCMap.clear();
915 
916  unsigned int ncstat = geom->Ncryostats();
917  hitmap.resize(ncstat);
918  for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
919  unsigned int ntpc = geom->Cryostat(cstat).NTPC();
920  hitmap[cstat].resize(ntpc);
921  for(unsigned int tpc = 0; tpc < ntpc; ++tpc) {
922  int nplane = geom->Cryostat(cstat).TPC(tpc).Nplanes();
923  hitmap[cstat][tpc].resize(nplane);
924  }
925  }
926 
927  for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end(); ++ihit) {
928  const art::Ptr<recob::Hit>& phit = *ihit;
929  geo::View_t view = phit->View();
930  if((view == geo::kU && fEnableU) ||
931  (view == geo::kV && fEnableV) ||
932  (view == geo::kZ && fEnableW)) {
933  geo::WireID phitWireID = phit->WireID();
934  hitmap[phitWireID.Cryostat][phitWireID.TPC][phitWireID.Plane].insert(std::make_pair(phitWireID.Wire, phit));
935  }
936  }
937 
938  // Fill mc information, including IDEs and closest neighbors
939  // of each hit.
942  if(useMC) {
943  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
944 
945  // First loop over hits and fill track ids and mc position.
946  for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
947  for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
948  int nplane = geom->Cryostat(cstat).TPC(tpc).Nplanes();
949  for(int plane = 0; plane < nplane; ++plane) {
950  for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator ihit = hitmap[cstat][tpc][plane].begin();
951  ihit != hitmap[cstat][tpc][plane].end(); ++ihit) {
952  const art::Ptr<recob::Hit>& phit = ihit->second;
953  const recob::Hit& hit = *phit;
954  HitMCInfo& mcinfo = fHitMCMap[&hit]; // Default HitMCInfo.
955 
956  // Fill default nearest neighbor information (i.e. none).
957 
958  mcinfo.pchit.resize(nplane, 0);
959  mcinfo.dist2.resize(nplane, 1.e20);
960 
961  // Get sim::IDEs for this hit.
962 
963  std::vector<sim::IDE> ides = bt_serv->HitToAvgSimIDEs(phit);
964 
965  // Get sorted track ids. for this hit.
966 
967  mcinfo.trackIDs.reserve(ides.size());
968  for(std::vector<sim::IDE>::const_iterator i = ides.begin();
969  i != ides.end(); ++i)
970  mcinfo.trackIDs.push_back(i->trackID);
971  sort(mcinfo.trackIDs.begin(), mcinfo.trackIDs.end());
972 
973  // Get position of ionization for this hit.
974 
975  try {
976  mcinfo.xyz = bt_serv->SimIDEsToXYZ(ides);
977  }
978  catch(cet::exception& x) {
979  mcinfo.xyz.clear();
980  }
981  } // end loop over ihit
982  }// end loop oer planes
983  }// end loop over TPCs
984  }// end loop over cryostats
985 
986  // Loop over hits again and fill nearest neighbor information for real.
987  for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
988  for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
989  int nplane = geom->Cryostat(cstat).TPC(tpc).Nplanes();
990  for(int plane = 0; plane < nplane; ++plane) {
991  for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator ihit = hitmap[cstat][tpc][plane].begin();
992  ihit != hitmap[cstat][tpc][plane].end(); ++ihit) {
993  const art::Ptr<recob::Hit>& phit = ihit->second;
994  const recob::Hit& hit = *phit;
995  HitMCInfo& mcinfo = fHitMCMap[&hit];
996  if(mcinfo.xyz.size() != 0) {
997  assert(mcinfo.xyz.size() == 3);
998 
999  // Fill nearest neighbor information for this hit.
1000 
1001  for(int plane2 = 0; plane2 < nplane; ++plane2) {
1002  for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator jhit = hitmap[cstat][tpc][plane2].begin();
1003  jhit != hitmap[cstat][tpc][plane2].end(); ++jhit) {
1004  const art::Ptr<recob::Hit>& phit2 = jhit->second;
1005  const recob::Hit& hit2 = *phit2;
1006  const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
1007 
1008 
1009  if(mcinfo2.xyz.size() != 0) {
1010  assert(mcinfo2.xyz.size() == 3);
1011  double dx = mcinfo.xyz[0] - mcinfo2.xyz[0];
1012  double dy = mcinfo.xyz[1] - mcinfo2.xyz[1];
1013  double dz = mcinfo.xyz[2] - mcinfo2.xyz[2];
1014  double dist2 = dx*dx + dy*dy + dz*dz;
1015  if(dist2 < mcinfo.dist2[plane2]) {
1016  mcinfo.dist2[plane2] = dist2;
1017  mcinfo.pchit[plane2] = &hit2;
1018  }
1019  }// end if mcinfo2.xyz valid
1020  }// end loop over jhit
1021  }// end loop over plane2
1022  }// end if mcinfo.xyz valid.
1023  }// end loop over ihit
1024  }// end loop over plane
1025  }// end loop over tpc
1026  }// end loop over cryostats
1027  }// end if MC
1028 
1029  // use mf::LogDebug instead of LOG_DEBUG because we reuse it in many lines
1030  // insertions are protected by mf::isDebugEnabled()
1031  mf::LogDebug debug("SpacePointAlg");
1032  if (mf::isDebugEnabled()) {
1033  debug << "Total hits = " << hits.size() << "\n\n";
1034 
1035  for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
1036  for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
1037  int nplane = hitmap[cstat][tpc].size();
1038  for(int plane = 0; plane < nplane; ++plane) {
1039  debug << "TPC, Plane: " << tpc << ", " << plane
1040  << ", hits = " << hitmap[cstat][tpc][plane].size() << "\n";
1041  }
1042  }
1043  } // end loop over cryostats
1044  } // if debug
1045 
1046  // Make empty multimap from hit pointer on preferred
1047  // (most-populated or collection) plane to space points that
1048  // include that hit (used for sorting, filtering, and
1049  // merging).
1050 
1051  typedef const recob::Hit* sptkey_type;
1052  std::multimap<sptkey_type, recob::SpacePoint> sptmap;
1053  std::set<sptkey_type> sptkeys; // Keys of multimap.
1054 
1055  // Loop over TPCs.
1056  for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
1057  for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
1058 
1059  geo::TPCID tpcid(cstat, tpc);
1060 
1061  // Sort maps in increasing order of number of hits.
1062  // This is so that we can do the outer loops over hits
1063  // over the views with fewer hits.
1064  //
1065  // If config parameter PreferColl is true, treat the colleciton
1066  // plane as if it had the most hits, regardless of how many
1067  // hits it actually has. This will force space points to be
1068  // filtered and merged with respect to the collection plane
1069  // wires. It will also force space points to be sorted by
1070  // collection plane wire.
1071 
1072  int nplane = hitmap[cstat][tpc].size();
1073  std::vector<int> index(nplane);
1074 
1075  for(int i=0; i<nplane; ++i)
1076  index[i] = i;
1077 
1078  for(int i=0; i<nplane-1; ++i) {
1079 
1080  for(int j=i+1; j<nplane; ++j) {
1081  bool icoll = fPreferColl &&
1082  geom->SignalType(geo::PlaneID(tpcid, index[i])) == geo::kCollection;
1083  bool jcoll = fPreferColl &&
1084  geom->SignalType(geo::PlaneID(tpcid, index[j])) == geo::kCollection;
1085  if((hitmap[cstat][tpc][index[i]].size() > hitmap[cstat][tpc][index[j]].size() &&
1086  !jcoll) || icoll) {
1087  int temp = index[i];
1088  index[i] = index[j];
1089  index[j] = temp;
1090  }
1091  }
1092  }// end loop over i
1093 
1094  // how many views with hits?
1095  // This will allow for the special case where we might have only 2 planes of information and
1096  // still want space points even if a three plane TPC
1097  std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>& hitsByPlaneVec = hitmap[cstat][tpc];
1098  int nViewsWithHits(0);
1099 
1100  for(int i = 0; i < nplane; i++)
1101  {
1102  if (hitsByPlaneVec[index[i]].size() > 0) nViewsWithHits++;
1103  }
1104 
1105  // If two-view space points are allowed, make a double loop
1106  // over hits and produce space points for compatible hit-pairs.
1107 
1108  if((nViewsWithHits == 2 || nplane == 2) && fMinViews <= 2) {
1109 
1110  // Loop over pairs of views.
1111  for(int i=0; i<nplane-1; ++i) {
1112  unsigned int plane1 = index[i];
1113 
1114  if (hitmap[cstat][tpc][plane1].empty()) continue;
1115 
1116  for(int j=i+1; j<nplane; ++j) {
1117  unsigned int plane2 = index[j];
1118 
1119  if (hitmap[cstat][tpc][plane2].empty()) continue;
1120 
1121  // Get angle, pitch, and offset of plane2 wires.
1122  const geo::WireGeo& wgeo2 = geom->Cryostat(cstat).TPC(tpc).Plane(plane2).Wire(0);
1123  double hl2 = wgeo2.HalfL();
1124  double xyz21[3];
1125  double xyz22[3];
1126  wgeo2.GetCenter(xyz21, -hl2);
1127  wgeo2.GetCenter(xyz22, hl2);
1128  double s2 = (xyz22[1] - xyz21[1]) / (2.*hl2);
1129  double c2 = (xyz22[2] - xyz21[2]) / (2.*hl2);
1130  double dist2 = -xyz21[1] * c2 + xyz21[2] * s2;
1131  double pitch2 = geom->WirePitch(plane2, tpc, cstat);
1132 
1133  if(!fPreferColl && hitmap[cstat][tpc][plane1].size() > hitmap[cstat][tpc][plane2].size())
1134  throw cet::exception("SpacePointAlg") << "makeSpacePoints(): hitmaps with incompatible size\n";
1135 
1136 
1137  // Loop over pairs of hits.
1138 
1139  art::PtrVector<recob::Hit> hitvec;
1140  hitvec.reserve(2);
1141 
1142  for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
1143  ihit1 = hitmap[cstat][tpc][plane1].begin();
1144  ihit1 != hitmap[cstat][tpc][plane1].end(); ++ihit1) {
1145 
1146  const art::Ptr<recob::Hit>& phit1 = ihit1->second;
1147  geo::WireID phit1WireID = phit1->WireID();
1148  const geo::WireGeo& wgeo = geom->WireIDToWireGeo(phit1WireID);
1149 
1150  // Get endpoint coordinates of this wire.
1151  // (kept as assertions for performance reasons)
1152  assert(phit1WireID.Cryostat == cstat);
1153  assert(phit1WireID.TPC == tpc);
1154  assert(phit1WireID.Plane == plane1);
1155  double hl1 = wgeo.HalfL();
1156  double xyz1[3];
1157  double xyz2[3];
1158  wgeo.GetCenter(xyz1, -hl1);
1159  wgeo.GetCenter(xyz2, hl1);
1160 
1161  // Find the plane2 wire numbers corresponding to the endpoints.
1162 
1163  double wire21 = (-xyz1[1] * c2 + xyz1[2] * s2 - dist2) / pitch2;
1164  double wire22 = (-xyz2[1] * c2 + xyz2[2] * s2 - dist2) / pitch2;
1165 
1166  int wmin = std::max(0., std::min(wire21, wire22));
1167  int wmax = std::max(0., std::max(wire21, wire22) + 1.);
1168 
1169  std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
1170  ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
1171  ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
1172 
1173  for(; ihit2 != ihit2end; ++ihit2) {
1174 
1175  const art::Ptr<recob::Hit>& phit2 = ihit2->second;
1176 
1177  // Check current pair of hits for compatibility.
1178  // By construction, hits should always have compatible views
1179  // and times, but may not have compatible mc information.
1180 
1181  hitvec.clear();
1182  hitvec.push_back(phit1);
1183  hitvec.push_back(phit2);
1184  bool ok = compatible(hitvec, useMC);
1185  if(ok) {
1186 
1187  // Add a space point.
1188 
1189  ++n2;
1190 
1191  // make a dummy vector of recob::SpacePoints
1192  // as we are filtering or merging and don't want to
1193  // add the created SpacePoint to the final collection just yet
1194  // This dummy vector will hold just one recob::SpacePoint,
1195  // which will go into the multimap and then the vector
1196  // will go out of scope.
1197 
1198  std::vector<recob::SpacePoint> sptv;
1199  fillSpacePoint(hitvec, sptv, sptmap.size());
1200  sptkey_type key = &*phit2;
1201  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
1202  sptkeys.insert(key);
1203  }
1204  }
1205  }
1206  }
1207  }
1208  }// end if fMinViews <= 2
1209 
1210  // If three-view space points are allowed, make a triple loop
1211  // over hits and produce space points for compatible triplets.
1212 
1213  if(nplane >= 3 && fMinViews <= 3) {
1214 
1215  // Loop over triplets of hits.
1216 
1217  art::PtrVector<recob::Hit> hitvec;
1218  hitvec.reserve(3);
1219 
1220  unsigned int plane1 = index[0];
1221  unsigned int plane2 = index[1];
1222  unsigned int plane3 = index[2];
1223 
1224  // Get angle, pitch, and offset of plane1 wires.
1225 
1226  const geo::WireGeo& wgeo1 = geom->Cryostat(cstat).TPC(tpc).Plane(plane1).Wire(0);
1227  double hl1 = wgeo1.HalfL();
1228  double xyz11[3];
1229  double xyz12[3];
1230  wgeo1.GetCenter(xyz11, -hl1);
1231  wgeo1.GetCenter(xyz12, hl1);
1232  double s1 = (xyz12[1] - xyz11[1]) / (2.*hl1);
1233  double c1 = (xyz12[2] - xyz11[2]) / (2.*hl1);
1234  double dist1 = -xyz11[1] * c1 + xyz11[2] * s1;
1235  double pitch1 = geom->WirePitch(plane1, tpc, cstat);
1236  const double TicksOffset1 = detprop->GetXTicksOffset(plane1,tpc,cstat);
1237 
1238  // Get angle, pitch, and offset of plane2 wires.
1239 
1240  const geo::WireGeo& wgeo2 = geom->Cryostat(cstat).TPC(tpc).Plane(plane2).Wire(0);
1241  double hl2 = wgeo2.HalfL();
1242  double xyz21[3];
1243  double xyz22[3];
1244  wgeo2.GetCenter(xyz21, -hl2);
1245  wgeo2.GetCenter(xyz22, hl2);
1246  double s2 = (xyz22[1] - xyz21[1]) / (2.*hl2);
1247  double c2 = (xyz22[2] - xyz21[2]) / (2.*hl2);
1248  double dist2 = -xyz21[1] * c2 + xyz21[2] * s2;
1249  double pitch2 = geom->WirePitch(plane2, tpc, cstat);
1250  const double TicksOffset2 = detprop->GetXTicksOffset(plane2,tpc,cstat);
1251 
1252  // Get angle, pitch, and offset of plane3 wires.
1253 
1254  const geo::WireGeo& wgeo3 = geom->Cryostat(cstat).TPC(tpc).Plane(plane3).Wire(0);
1255  double hl3 = wgeo3.HalfL();
1256  double xyz31[3];
1257  double xyz32[3];
1258  wgeo3.GetCenter(xyz31, -hl3);
1259  wgeo3.GetCenter(xyz32, hl3);
1260  double s3 = (xyz32[1] - xyz31[1]) / (2.*hl3);
1261  double c3 = (xyz32[2] - xyz31[2]) / (2.*hl3);
1262  double dist3 = -xyz31[1] * c3 + xyz31[2] * s3;
1263  double pitch3 = geom->WirePitch(plane3, tpc, cstat);
1264  const double TicksOffset3 = detprop->GetXTicksOffset(plane3,tpc,cstat);
1265 
1266  // Get sine of angle differences.
1267 
1268  double s12 = s1 * c2 - s2 * c1; // sin(theta1 - theta2).
1269  double s23 = s2 * c3 - s3 * c2; // sin(theta2 - theta3).
1270  double s31 = s3 * c1 - s1 * c3; // sin(theta3 - theta1).
1271 
1272  // Loop over hits in plane1.
1273 
1274  std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
1275  ihit1 = hitmap[cstat][tpc][plane1].begin(),
1276  ihit1end = hitmap[cstat][tpc][plane1].end();
1277  for(; ihit1 != ihit1end; ++ihit1) {
1278 
1279  unsigned int wire1 = ihit1->first;
1280  const art::Ptr<recob::Hit>& phit1 = ihit1->second;
1281  geo::WireID phit1WireID = phit1->WireID();
1282  const geo::WireGeo& wgeo = geom->WireIDToWireGeo(phit1WireID);
1283 
1284  // Get endpoint coordinates of this wire from plane1.
1285  // (kept as assertions for performance reasons)
1286  assert(phit1WireID.Cryostat == cstat);
1287  assert(phit1WireID.TPC == tpc);
1288  assert(phit1WireID.Plane == plane1);
1289  assert(phit1WireID.Wire == wire1);
1290  double hl1 = wgeo.HalfL();
1291  double xyz1[3];
1292  double xyz2[3];
1293  wgeo.GetCenter(xyz1, -hl1);
1294  wgeo.GetCenter(xyz2, hl1);
1295 
1296  // Get corrected time and oblique coordinate of first hit.
1297 
1298  double t1 = phit1->PeakTime() - TicksOffset1;
1299  double u1 = wire1 * pitch1 + dist1;
1300 
1301  // Find the plane2 wire numbers corresponding to the endpoints.
1302 
1303  double wire21 = (-xyz1[1] * c2 + xyz1[2] * s2 - dist2) / pitch2;
1304  double wire22 = (-xyz2[1] * c2 + xyz2[2] * s2 - dist2) / pitch2;
1305 
1306  int wmin = std::max(0., std::min(wire21, wire22));
1307  int wmax = std::max(0., std::max(wire21, wire22) + 1.);
1308 
1309  std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
1310  ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
1311  ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
1312 
1313  for(; ihit2 != ihit2end; ++ihit2) {
1314 
1315  int wire2 = ihit2->first;
1316  const art::Ptr<recob::Hit>& phit2 = ihit2->second;
1317 
1318  // Get corrected time of second hit.
1319 
1320  double t2 = phit2->PeakTime() - TicksOffset2;
1321 
1322  // Check maximum time difference with first hit.
1323 
1324  bool dt12ok = std::abs(t1-t2) <= fMaxDT;
1325  if(dt12ok) {
1326 
1327  // Test first two hits for compatibility before looping
1328  // over third hit.
1329 
1330  hitvec.clear();
1331  hitvec.push_back(phit1);
1332  hitvec.push_back(phit2);
1333  bool h12ok = compatible(hitvec, useMC);
1334  if(h12ok) {
1335 
1336  // Get oblique coordinate of second hit.
1337 
1338  double u2 = wire2 * pitch2 + dist2;
1339 
1340  // Predict plane3 oblique coordinate and wire number.
1341 
1342  double u3pred = (-u1*s23 - u2*s31) / s12;
1343  double w3pred = (u3pred - dist3) / pitch3;
1344  double w3delta = std::abs(fMaxS / (s12 * pitch3));
1345  int w3min = std::max(0., std::ceil(w3pred - w3delta));
1346  int w3max = std::max(0., std::floor(w3pred + w3delta));
1347 
1348  std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
1349  ihit3 = hitmap[cstat][tpc][plane3].lower_bound(w3min),
1350  ihit3end = hitmap[cstat][tpc][plane3].upper_bound(w3max);
1351 
1352  for(; ihit3 != ihit3end; ++ihit3) {
1353 
1354  int wire3 = ihit3->first;
1355  const art::Ptr<recob::Hit>& phit3 = ihit3->second;
1356 
1357  // Get corrected time of third hit.
1358 
1359  double t3 = phit3->PeakTime() - TicksOffset3;
1360 
1361  // Check time difference of third hit compared to first two hits.
1362 
1363  bool dt123ok = std::abs(t1-t3) <= fMaxDT && std::abs(t2-t3) <= fMaxDT;
1364  if(dt123ok) {
1365 
1366  // Get oblique coordinate of third hit and check spatial separation.
1367 
1368  double u3 = wire3 * pitch3 + dist3;
1369  double S = s23 * u1 + s31 * u2 + s12 * u3;
1370  bool sok = std::abs(S) <= fMaxS;
1371  if(sok) {
1372 
1373  // Test triplet for compatibility.
1374 
1375  hitvec.clear();
1376  hitvec.push_back(phit1);
1377  hitvec.push_back(phit2);
1378  hitvec.push_back(phit3);
1379  bool h123ok = compatible(hitvec, useMC);
1380  if(h123ok) {
1381 
1382  // Add a space point.
1383 
1384  ++n3;
1385 
1386  // make a dummy vector of recob::SpacePoints
1387  // as we are filtering or merging and don't want to
1388  // add the created SpacePoint to the final collection just yet
1389  // This dummy vector will hold just one recob::SpacePoint,
1390  // which will go into the multimap and then the vector
1391  // will go out of scope.
1392 
1393  std::vector<recob::SpacePoint> sptv;
1394  fillSpacePoint(hitvec, sptv, sptmap.size()-1);
1395  sptkey_type key = &*phit3;
1396  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
1397  sptkeys.insert(key);
1398  }
1399  }
1400  }
1401  }
1402  }
1403  }
1404  }
1405  }
1406  }// end if fMinViews <= 3
1407 
1408  // Do Filtering.
1409 
1410  if(fFilter) {
1411 
1412  // Transfer (some) space points from sptmap to spts.
1413 
1414  spts.reserve(spts.size() + sptkeys.size());
1415 
1416  // Loop over keys of space point map.
1417  // Space points that have the same key are candidates for filtering.
1418 
1419  for(std::set<sptkey_type>::const_iterator i = sptkeys.begin();
1420  i != sptkeys.end(); ++i) {
1421  sptkey_type key = *i;
1422 
1423  // Loop over space points corresponding to the current key.
1424  // Choose the single best space point from among this group.
1425 
1426  double best_chisq = 0.;
1427  const recob::SpacePoint* best_spt = 0;
1428 
1429  for(std::multimap<sptkey_type, recob::SpacePoint>::const_iterator j = sptmap.lower_bound(key);
1430  j != sptmap.upper_bound(key); ++j) {
1431  const recob::SpacePoint& spt = j->second;
1432  if(best_spt == 0 || spt.Chisq() < best_chisq) {
1433  best_spt = &spt;
1434  best_chisq = spt.Chisq();
1435  }
1436  }
1437 
1438  // Transfer best filtered space point to result vector.
1439 
1440  if (!best_spt)
1441  throw cet::exception("SpacePointAlg") << "makeSpacePoints(): no best point\n";
1442  spts.push_back(*best_spt);
1443  if(fMinViews <= 2)
1444  ++n2filt;
1445  else
1446  ++n3filt;
1447  }
1448  }// end if filtering
1449 
1450  // Do merging.
1451 
1452  else if(fMerge) {
1453 
1454  // Transfer merged space points from sptmap to spts.
1455 
1456  spts.reserve(spts.size() + sptkeys.size());
1457 
1458  // Loop over keys of space point map.
1459  // Space points that have the same key are candidates for merging.
1460 
1461  for(std::set<sptkey_type>::const_iterator i = sptkeys.begin();
1462  i != sptkeys.end(); ++i) {
1463  sptkey_type key = *i;
1464 
1465  // Loop over space points corresponding to the current key.
1466  // Make a collection of hits that is the union of the hits
1467  // from each candidate space point.
1468 
1469  std::multimap<sptkey_type, recob::SpacePoint>::const_iterator
1470  jSPT = sptmap.lower_bound(key), jSPTend = sptmap.upper_bound(key);
1471 
1472  art::PtrVector<recob::Hit> merged_hits;
1473  for(; jSPT != jSPTend; ++jSPT) {
1474  const recob::SpacePoint& spt = jSPT->second;
1475 
1476  // Loop over hits from this space points.
1477  // Add each hit to the collection of all hits.
1478 
1479  const art::PtrVector<recob::Hit>& spt_hits = getAssociatedHits(spt);
1480  merged_hits.reserve(merged_hits.size() + spt_hits.size()); // better than nothing, but not ideal
1481  for(art::PtrVector<recob::Hit>::const_iterator k = spt_hits.begin();
1482  k != spt_hits.end(); ++k) {
1483  const art::Ptr<recob::Hit>& hit = *k;
1484  merged_hits.push_back(hit);
1485  }
1486  }
1487 
1488  // Remove duplicates.
1489 
1490  std::sort(merged_hits.begin(), merged_hits.end());
1491  art::PtrVector<recob::Hit>::iterator it =
1492  std::unique(merged_hits.begin(), merged_hits.end());
1493  merged_hits.erase(it, merged_hits.end());
1494 
1495  // Construct a complex space points using merged hits.
1496 
1497  fillComplexSpacePoint(merged_hits, spts, sptmap.size() + spts.size()-1);
1498 
1499  if(fMinViews <= 2)
1500  ++n2filt;
1501  else
1502  ++n3filt;
1503  }
1504  }// end if merging
1505 
1506  // No filter, no merge.
1507 
1508  else {
1509 
1510  // Transfer all space points from sptmap to spts.
1511 
1512  spts.reserve(spts.size() + sptkeys.size());
1513 
1514  // Loop over space points.
1515 
1516  for(std::multimap<sptkey_type, recob::SpacePoint>::const_iterator j = sptmap.begin();
1517  j != sptmap.end(); ++j) {
1518  const recob::SpacePoint& spt = j->second;
1519  spts.push_back(spt);
1520  }
1521 
1522  // Update statistics.
1523 
1524  n2filt = n2;
1525  n3filt = n3;
1526  }
1527  }// end loop over tpcs
1528  }// end loop over cryostats
1529 
1530  if (mf::isDebugEnabled()) {
1531  debug << "\n2-hit space points = " << n2 << "\n"
1532  << "3-hit space points = " << n3 << "\n"
1533  << "2-hit filtered/merged space points = " << n2filt << "\n"
1534  << "3-hit filtered/merged space points = " << n3filt;
1535  } // if debug
1536  }
1537 
1538  //----------------------------------------------------------------------
1539  // Get hits associated with a particular space point, based on most recent
1540  // call of any make*SpacePoints method.
1541  const art::PtrVector<recob::Hit>&
1542  SpacePointAlg::getAssociatedHits(const recob::SpacePoint& spt) const
1543  {
1544  // It is an error if no hits are associated with this space point (throw exception).
1545 
1546  std::map<int, art::PtrVector<recob::Hit> >::const_iterator it =
1547  fSptHitMap.find(spt.ID());
1548  if(it == fSptHitMap.end())
1549  {
1550  mf::LogWarning("SpacePointAlg") <<"Looking for ID " << spt.ID()<< " from " << fSptHitMap.size()<<std::endl;
1551  throw cet::exception("SpacePointAlg") << "No Hits associated with space point.\n";
1552  }
1553  return (*it).second;
1554 
1555  }
1556 
1557 
1558 } // end namespace trkf
x
Float_t x
Definition: compare.C:6
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
trkf::SpacePointAlg::correctedTime
double correctedTime(const recob::Hit &hit) const
Definition: SpacePointAlg.cxx:218
geo::WireGeo
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:61
s
Float_t s
Definition: plot.C:23
art::PtrVector::reserve
void reserve(size_type n)
Definition: PtrVector.h:343
KHitTrack.h
Basic Kalman filter track class, plus one measurement on same surface.
trkf::KHitTrack::getHit
const std::shared_ptr< const KHitBase > & getHit() const
Measurement.
Definition: KHitTrack.h:53
cheat::BackTrackerService::SimIDEsToXYZ
const std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides)
Definition: BackTrackerService_service.cc:238
trkf::SpacePointAlg::HitMCInfo::pchit
std::vector< const recob::Hit * > pchit
Pointer to nearest neighbor hit (indexed by plane).
Definition: SpacePointAlg.h:179
trkf::SpacePointAlg::fHitMCMap
std::map< const recob::Hit *, HitMCInfo > fHitMCMap
Definition: SpacePointAlg.h:182
trkf::SpacePointAlg::fTickOffsetU
double fTickOffsetU
Tick offset for plane U.
Definition: SpacePointAlg.h:169
t1
TTree * t1
Definition: plottest35.C:26
CryostatGeo.h
Encapsulate the construction of a single cyostat.
trkf::SpacePointAlg::fMaxDT
double fMaxDT
Maximum time difference between planes.
Definition: SpacePointAlg.h:160
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
mf::MaybeLogger_
Definition: MessageLogger.h:114
geo::PlaneGeo::Wire
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
trkf::SpacePointAlg::fFilter
bool fFilter
Filter flag.
Definition: SpacePointAlg.h:166
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:212
geo::kV
Planes which measure V.
Definition: geo_types.h:77
recob::Hit::WireID
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:234
geo::Orient_t
enum geo::_plane_orient Orient_t
Enumerate the possible plane projections.
geo::TPCGeo::Nplanes
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:145
art::PtrVector::begin
iterator begin()
Definition: PtrVector.h:223
Hit.h
Declaration of signal hit object.
KHitWireX.h
Kalman filter wire-time measurement on a SurfWireX surface.
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:250
trkf::SpacePointAlg::HitMCInfo::dist2
std::vector< double > dist2
Distance to nearest neighbor hit (indexed by plane).
Definition: SpacePointAlg.h:180
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:37
trkf::SpacePointAlg::fEnableW
bool fEnableW
Enable flag (W).
Definition: SpacePointAlg.h:165
art::PtrVector::erase
iterator erase(iterator position)
Definition: PtrVector.h:508
trkf::SpacePointAlg::HitMCInfo::trackIDs
std::vector< int > trackIDs
Parent trackIDs.
Definition: SpacePointAlg.h:177
trkf::KHitWireX::getHit
const art::Ptr< recob::Hit > & getHit() const
Get original hit.
Definition: KHitWireX.h:47
trkf::SpacePointAlg::update
void update() const
Definition: SpacePointAlg.cxx:114
recob::SpacePoint::Chisq
Double32_t Chisq() const
Definition: SpacePoint.h:67
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:79
recob::Hit::View
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:233
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:313
geo::GeometryCore::SignalType
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
Definition: GeometryCore.cxx:537
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:214
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:1964
geo::WireGeo::ThetaZ
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:192
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:783
trkf::SpacePointAlg::getAssociatedHits
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
Definition: SpacePointAlg.cxx:1542
trkf::SpacePointAlg::makeMCTruthSpacePoints
void makeMCTruthSpacePoints(const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
Definition: SpacePointAlg.cxx:871
trkf::SpacePointAlg::separation
double separation(const art::PtrVector< recob::Hit > &hits) const
Definition: SpacePointAlg.cxx:240
trkf::SpacePointAlg::fillSpacePoints
void fillSpacePoints(std::vector< recob::SpacePoint > &spts, std::multimap< double, KHitTrack > const &trackMap) const
Fill a collection of space points.
Definition: SpacePointAlg.cxx:625
geo::kU
Planes which measure U.
Definition: geo_types.h:76
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:171
max
Int_t max
Definition: plot.C:27
geo::kHorizontal
Planes that are in the horizontal plane.
Definition: geo_types.h:87
hits
void hits()
Definition: readHits.C:15
trkf::SpacePointAlg::SpacePointAlg
SpacePointAlg(const fhicl::ParameterSet &pset)
Definition: SpacePointAlg.cxx:43
c1
TCanvas * c1
Definition: plotHisto.C:7
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:25
trkf::SpacePointAlg::fEnableU
bool fEnableU
Enable flag (U).
Definition: SpacePointAlg.h:163
c2
TCanvas * c2
Definition: plot_hist.C:75
art::PtrVector::push_back
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
geo::kInduction
Signal from induction planes.
Definition: geo_types.h:92
tca::debug
DebugStuff debug
Definition: DebugStruct.cxx:4
trkf::SpacePointAlg::fSptHitMap
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
Definition: SpacePointAlg.h:183
geo::kVertical
Planes that are in the vertical plane (e.g. ArgoNeuT).
Definition: geo_types.h:88
trkf::SpacePointAlg::HitMCInfo::xyz
std::vector< double > xyz
Location of ionization (all tracks).
Definition: SpacePointAlg.h:178
detinfo::DetectorProperties::GetXTicksCoefficient
virtual double GetXTicksCoefficient(int t, int c) const =0
geo::SigType_t
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
evd::details::begin
std::vector< evd::details::RawDigitInfo_t >::const_iterator begin(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:314
trkf::SpacePointAlg::fillSpacePoint
void fillSpacePoint(const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
Definition: SpacePointAlg.cxx:474
geo::PlaneGeo
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:78
art::PtrVector::end
iterator end()
Definition: PtrVector.h:237
geo::CryostatGeo::NTPC
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:155
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:280
mf::isDebugEnabled
bool isDebugEnabled()
Definition: MessageLogger.cc:745
geo::TPCGeo::WirePitch
double WirePitch(unsigned plane=0) const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.cxx:427
trkf::SpacePointAlg::fMaxS
double fMaxS
Maximum space separation between wires.
Definition: SpacePointAlg.h:161
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:195
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
trkf::SpacePointAlg::fPreferColl
bool fPreferColl
Sort by collection wire.
Definition: SpacePointAlg.h:168
geo::TPCGeo::Plane0Pitch
double Plane0Pitch(unsigned int p) const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.cxx:355
art::PtrVector::iterator
data_t::iterator iterator
Definition: PtrVector.h:60
art::PtrVector::size
size_type size() const
Definition: PtrVector.h:308
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
t2
TTree * t2
Definition: plottest35.C:36
trkf::SpacePointAlg::reconfigure
void reconfigure(const fhicl::ParameterSet &pset)
Definition: SpacePointAlg.cxx:70
geo::PlaneGeo::Orientation
Orient_t Orientation() const
What is the orientation of the plane.
Definition: PlaneGeo.h:174
WireGeo.h
Encapsulate the geometry of a wire.
art::PtrVector::const_iterator
data_t::const_iterator const_iterator
Definition: PtrVector.h:61
geo::WireGeo::HalfL
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:106
art::PtrVector::front
reference front()
Definition: PtrVector.h:379
trkf::SpacePointAlg::makeSpacePoints
void makeSpacePoints(const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
Definition: SpacePointAlg.cxx:861
sw
Float_t sw
Definition: plot.C:23
weight
double weight
Definition: plottest35.C:25
recob::SpacePoint::ID
ID_t ID() const
Definition: SpacePoint.h:64
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:219
trkf::SpacePointAlg::fTickOffsetV
double fTickOffsetV
Tick offset for plane V.
Definition: SpacePointAlg.h:170
PlaneGeo.h
Encapsulate the construction of a single detector plane.
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:174
min
Int_t min
Definition: plot.C:26
trkf::SpacePointAlg::fillComplexSpacePoint
void fillComplexSpacePoint(const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
Definition: SpacePointAlg.cxx:684
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:213
trkf::SpacePointAlg::fTickOffsetW
double fTickOffsetW
Tick offset for plane W.
Definition: SpacePointAlg.h:171
trkf::SpacePointAlg::fMinViews
int fMinViews
Mininum number of views per space point.
Definition: SpacePointAlg.h:162
trkf::SpacePointAlg::numHitMap
int numHitMap() const
Definition: SpacePointAlg.h:148
recob::Hit::SigmaPeakTime
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:220
geo::kW
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:78
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:299
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:49
geo::WireGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:68
trkf::SpacePointAlg::compatible
bool compatible(const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
Definition: SpacePointAlg.cxx:323
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
SpacePointAlg.h
Algorithm for generating space points from hits.
trkf::SpacePointAlg::fMerge
bool fMerge
Merge flag.
Definition: SpacePointAlg.h:167
track
Float_t track
Definition: plot.C:34
art::PtrVector::clear
void clear()
Definition: PtrVector.h:537
w
Float_t w
Definition: plot.C:23
geo::TPCGeo::PlaneLocation
const double * PlaneLocation(unsigned int p) const
Returns the coordinates of the center of the specified plane [cm].
Definition: TPCGeo.cxx:413
cheat::BackTrackerService::HitToAvgSimIDEs
const std::vector< sim::IDE > HitToAvgSimIDEs(recob::Hit const &hit)
Definition: BackTrackerService_service.cc:210
Geometry.h
art framework interface to geometry description
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
TPCGeo.h
Encapsulate the construction of a single detector plane.
geo::kCollection
Signal from collection planes.
Definition: geo_types.h:93
trkf::SpacePointAlg::fEnableV
bool fEnableV
Enable flag (V).
Definition: SpacePointAlg.h:164
geo::GeometryCore::WireIDToWireGeo
WireGeo const & WireIDToWireGeo(geo::WireID const &wireid) const
Returns the specified wire.
Definition: GeometryCore.h:3372
detinfo::DetectorProperties::GetXTicksOffset
virtual double GetXTicksOffset(int p, int t, int c) const =0