LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
trkf::SpacePointAlg Class Reference

#include "SpacePointAlg.h"

Classes

struct  HitMCInfo
 

Public Member Functions

 SpacePointAlg (const fhicl::ParameterSet &pset)
 
bool filter () const noexcept
 
bool merge () const noexcept
 
double maxDT () const noexcept
 
double maxS () const noexcept
 
int minViews () const noexcept
 
bool enableU () const noexcept
 
bool enableV () const noexcept
 
bool enableW () const noexcept
 
void update (detinfo::DetectorPropertiesData const &detProp) const
 
double correctedTime (detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
 
double separation (const art::PtrVector< recob::Hit > &hits) const
 
bool compatible (detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
 
void fillSpacePoint (detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
 
void fillSpacePoints (detinfo::DetectorPropertiesData const &detProp, std::vector< recob::SpacePoint > &spts, std::multimap< double, KHitTrack > const &trackMap) const
 Fill a collection of space points. More...
 
void fillComplexSpacePoint (detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
 
void makeSpacePoints (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
 
void makeMCTruthSpacePoints (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
 
const art::PtrVector< recob::Hit > & getAssociatedHits (const recob::SpacePoint &spt) const
 
void clearHitMap () const
 
int numHitMap () const
 

Private Member Functions

void makeSpacePoints (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts, bool useMC) const
 

Private Attributes

double fMaxDT
 Maximum time difference between planes. More...
 
double fMaxS
 Maximum space separation between wires. More...
 
int fMinViews
 Mininum number of views per space point. More...
 
bool fEnableU
 Enable flag (U). More...
 
bool fEnableV
 Enable flag (V). More...
 
bool fEnableW
 Enable flag (W). More...
 
bool fFilter
 Filter flag. More...
 
bool fMerge
 Merge flag. More...
 
bool fPreferColl
 Sort by collection wire. More...
 
double fTickOffsetU
 Tick offset for plane U. More...
 
double fTickOffsetV
 Tick offset for plane V. More...
 
double fTickOffsetW
 Tick offset for plane W. More...
 
std::map< const recob::Hit *, HitMCInfofHitMCMap
 
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
 

Detailed Description

Definition at line 80 of file SpacePointAlg.h.

Constructor & Destructor Documentation

trkf::SpacePointAlg::SpacePointAlg ( const fhicl::ParameterSet pset)

Definition at line 38 of file SpacePointAlg.cxx.

References fEnableU, fEnableV, fEnableW, fFilter, fMaxDT, fMaxS, fMerge, fMinViews, fPreferColl, fTickOffsetU, fTickOffsetV, fTickOffsetW, and fhicl::ParameterSet::get().

39  : fMaxDT{pset.get<double>("MaxDT")}
40  , fMaxS{pset.get<double>("MaxS")}
41  , fMinViews{pset.get<int>("MinViews")}
42  , fEnableU{pset.get<bool>("EnableU")}
43  , fEnableV{pset.get<bool>("EnableV")}
44  , fEnableW{pset.get<bool>("EnableW")}
45  , fFilter{pset.get<bool>("Filter")}
46  , fMerge{pset.get<bool>("Merge")}
47  , fPreferColl{pset.get<bool>("PreferColl")}
48  , fTickOffsetU{pset.get<double>("TickOffsetU", 0.)}
49  , fTickOffsetV{pset.get<double>("TickOffsetV", 0.)}
50  , fTickOffsetW{pset.get<double>("TickOffsetW", 0.)}
51  {
52  // Only allow one of fFilter and fMerge to be true.
53 
54  if (fFilter && fMerge)
55  throw cet::exception("SpacePointAlg") << "Filter and Merge flags are both true.\n";
56 
57  // Report.
58 
59  std::cout << "SpacePointAlg configured with the following parameters:\n"
60  << " MaxDT = " << fMaxDT << "\n"
61  << " MaxS = " << fMaxS << "\n"
62  << " MinViews = " << fMinViews << "\n"
63  << " EnableU = " << fEnableU << "\n"
64  << " EnableV = " << fEnableV << "\n"
65  << " EnableW = " << fEnableW << "\n"
66  << " Filter = " << fFilter << "\n"
67  << " Merge = " << fMerge << "\n"
68  << " PreferColl = " << fPreferColl << "\n"
69  << " TickOffsetU = " << fTickOffsetU << "\n"
70  << " TickOffsetV = " << fTickOffsetV << "\n"
71  << " TickOffsetW = " << fTickOffsetW << std::endl;
72  }
double fTickOffsetU
Tick offset for plane U.
double fMaxDT
Maximum time difference between planes.
bool fFilter
Filter flag.
bool fEnableW
Enable flag (W).
bool fEnableU
Enable flag (U).
T get(std::string const &key) const
Definition: ParameterSet.h:314
double fMaxS
Maximum space separation between wires.
bool fPreferColl
Sort by collection wire.
double fTickOffsetV
Tick offset for plane V.
double fTickOffsetW
Tick offset for plane W.
int fMinViews
Mininum number of views per space point.
bool fMerge
Merge flag.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool fEnableV
Enable flag (V).

Member Function Documentation

void trkf::SpacePointAlg::clearHitMap ( ) const
inline
bool trkf::SpacePointAlg::compatible ( detinfo::DetectorPropertiesData const &  detProp,
const art::PtrVector< recob::Hit > &  hits,
bool  useMC = false 
) const

Definition at line 255 of file SpacePointAlg.cxx.

References util::abs(), geo::CryostatID::Cryostat, larg4::dist(), fHitMCMap, fMaxDT, fMaxS, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), trkf::SpacePointAlg::HitMCInfo::pchit, recob::Hit::PeakTime(), geo::PlaneID::Plane, art::PtrVector< T >::size(), t1, t2, geo::TPCID::TPC, trkf::SpacePointAlg::HitMCInfo::trackIDs, recob::Hit::View(), recob::Hit::WireID(), and geo::GeometryCore::WireIDToWireGeo().

Referenced by fillSpacePoints(), and makeSpacePoints().

258  {
260 
261  int nhits = hits.size();
262 
263  // Fewer than two or more than three hits can never be compatible.
264 
265  bool result = nhits >= 2 && nhits <= 3;
266  bool mc_ok = true;
267  unsigned int tpc = 0;
268  unsigned int cstat = 0;
269 
270  if (result) {
271 
272  // First do pairwise tests.
273  // Do double loop over hits.
274 
275  for (int ihit1 = 0; result && ihit1 < nhits - 1; ++ihit1) {
276  const recob::Hit& hit1 = *(hits[ihit1]);
277  geo::WireID hit1WireID = hit1.WireID();
278  geo::View_t view1 = hit1.View();
279 
280  double t1 = hit1.PeakTime() -
281  detProp.GetXTicksOffset(hit1WireID.Plane, hit1WireID.TPC, hit1WireID.Cryostat);
282 
283  // If using mc information, get a collection of track ids for hit 1.
284  // If not using mc information, this section of code will trigger the
285  // insertion of a single invalid HitMCInfo object into fHitMCMap.
286 
287  const HitMCInfo& mcinfo1 = fHitMCMap[(useMC ? &hit1 : 0)];
288  const std::vector<int>& tid1 = mcinfo1.trackIDs;
289  bool only_neg1 = tid1.size() > 0 && tid1.back() < 0;
290 
291  // Loop over second hit.
292 
293  for (int ihit2 = ihit1 + 1; result && ihit2 < nhits; ++ihit2) {
294  const recob::Hit& hit2 = *(hits[ihit2]);
295  geo::WireID hit2WireID = hit2.WireID();
296  geo::View_t view2 = hit2.View();
297 
298  // Test for same tpc and different views.
299 
300  result = result && hit1WireID.TPC == hit2WireID.TPC && view1 != view2 &&
301  hit1WireID.Cryostat == hit2WireID.Cryostat;
302  if (result) {
303 
304  // Remember which tpc and cryostat we are in.
305 
306  tpc = hit1WireID.TPC;
307  cstat = hit1WireID.Cryostat;
308 
309  double t2 = hit2.PeakTime() - detProp.GetXTicksOffset(
310  hit2WireID.Plane, hit2WireID.TPC, hit2WireID.Cryostat);
311 
312  // Test maximum time difference.
313 
314  result = result && std::abs(t1 - t2) <= fMaxDT;
315 
316  // Test mc truth.
317 
318  if (result && useMC) {
319 
320  // Test whether hits have a common parent track id.
321 
322  const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
323  std::vector<int> tid2 = mcinfo2.trackIDs;
324  bool only_neg2 = tid2.size() > 0 && tid2.back() < 0;
325  std::vector<int>::iterator it = std::set_intersection(
326  tid1.begin(), tid1.end(), tid2.begin(), tid2.end(), tid2.begin());
327  tid2.resize(it - tid2.begin());
328 
329  // Hits are compatible if they have parents in common.
330  // If the only parent id in common is negative (-999),
331  // then hits are compatible only if both hits have only
332  // negative parent tracks.
333 
334  bool only_neg3 = tid2.size() > 0 && tid2.back() < 0;
335  mc_ok = tid2.size() > 0 && (!only_neg3 || (only_neg1 && only_neg2));
336  result = result && mc_ok;
337 
338  // If we are still OK, check that either hit is
339  // the nearest neighbor of the other.
340 
341  if (result) {
342  result = mcinfo1.pchit[hit2WireID.Plane] == &hit2 ||
343  mcinfo2.pchit[hit1WireID.Plane] == &hit1;
344  }
345  }
346  }
347  }
348  }
349 
350  // If there are exactly three hits, and they pass pairwise tests, check
351  // for spatial compatibility.
352 
353  if (result && nhits == 3) {
354 
355  // Loop over hits.
356 
357  double dist[3] = {0., 0., 0.};
358  double sinth[3] = {0., 0., 0.};
359  double costh[3] = {0., 0., 0.};
360 
361  for (int i = 0; i < 3; ++i) {
362 
363  // Get tpc, plane, wire.
364 
365  const recob::Hit& hit = *(hits[i]);
366  geo::WireID hitWireID = hit.WireID();
367 
368  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
369  if ((hitWireID.TPC != tpc) || (hitWireID.Cryostat != cstat))
370  throw cet::exception("SpacePointAlg") << "compatible(): geometry mismatch\n";
371 
372  // Get angles and distance of wire.
373 
374  double const hl = wgeom.HalfL();
375  auto const xyz = wgeom.GetCenter();
376  auto const xyz1 = wgeom.GetEnd();
377  double s = (xyz1.Y() - xyz.Y()) / hl;
378  double c = (xyz1.Z() - xyz.Z()) / hl;
379  sinth[hit.WireID().Plane] = s;
380  costh[hit.WireID().Plane] = c;
381  dist[hit.WireID().Plane] = xyz.Z() * s - xyz.Y() * c;
382  }
383 
384  // Do space cut.
385 
386  double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0] +
387  (sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1] +
388  (sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
389 
390  result = result && std::abs(S) < fMaxS;
391  }
392  }
393 
394  // Done.
395 
396  return result;
397  }
intermediate_table::iterator iterator
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:114
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
std::map< const recob::Hit *, HitMCInfo > fHitMCMap
TTree * t1
Definition: plottest35.C:26
double fMaxDT
Maximum time difference between planes.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:276
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Point_t GetEnd() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:231
double fMaxS
Maximum space separation between wires.
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
TTree * t2
Definition: plottest35.C:36
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:172
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double trkf::SpacePointAlg::correctedTime ( detinfo::DetectorPropertiesData const &  detProp,
const recob::Hit hit 
) const

Definition at line 150 of file SpacePointAlg.cxx.

References fTickOffsetU, fTickOffsetV, fTickOffsetW, detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::kU, geo::kV, geo::kW, recob::Hit::PeakTime(), recob::Hit::View(), and recob::Hit::WireID().

Referenced by trkf::SpacePointAna::analyze().

152  {
153  // Get services.
154 
156 
157  // Correct time for trigger offset and plane-dependent time offsets.
158 
159  double t = hit.PeakTime() - detProp.GetXTicksOffset(hit.WireID());
160  if (hit.View() == geo::kU)
161  t -= fTickOffsetU;
162  else if (hit.View() == geo::kV)
163  t -= fTickOffsetV;
164  else if (hit.View() == geo::kW)
165  t -= fTickOffsetW;
166 
167  return t;
168  }
double fTickOffsetU
Tick offset for plane U.
Planes which measure V.
Definition: geo_types.h:136
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:276
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Planes which measure U.
Definition: geo_types.h:135
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
double fTickOffsetV
Tick offset for plane V.
double fTickOffsetW
Tick offset for plane W.
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:137
bool trkf::SpacePointAlg::enableU ( ) const
inlinenoexcept
bool trkf::SpacePointAlg::enableV ( ) const
inlinenoexcept
bool trkf::SpacePointAlg::enableW ( ) const
inlinenoexcept
void trkf::SpacePointAlg::fillComplexSpacePoint ( detinfo::DetectorPropertiesData const &  detProp,
const art::PtrVector< recob::Hit > &  hits,
std::vector< recob::SpacePoint > &  sptv,
int  sptid 
) const

Definition at line 603 of file SpacePointAlg.cxx.

References geo::PlaneID::asPlaneID(), art::PtrVector< T >::begin(), geo::CryostatID::Cryostat, art::PtrVector< T >::end(), art::PtrVector< T >::front(), fSptHitMap, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), detinfo::DetectorPropertiesData::GetXTicksCoefficient(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), hits(), geo::TPCGeo::Nplanes(), recob::Hit::PeakTime(), geo::PlaneID::Plane, recob::Hit::SigmaPeakTime(), art::PtrVector< T >::size(), ssc, geo::TPCID::TPC, geo::GeometryCore::TPC(), w, weight, recob::Hit::WireID(), geo::GeometryCore::WireIDToWireGeo(), and geo::GeometryCore::WirePitch().

Referenced by makeSpacePoints().

607  {
609 
610  // Calculate time pitch.
611 
612  double timePitch = detProp.GetXTicksCoefficient(); // cm / tick
613 
614  // Figure out which tpc we are in.
615 
616  unsigned int tpc0 = 0;
617  unsigned int cstat0 = 0;
618  int nhits = hits.size();
619  if (nhits > 0) {
620  tpc0 = hits.front()->WireID().TPC;
621  cstat0 = hits.front()->WireID().Cryostat;
622  }
623 
624  // Remember associated hits internally.
625 
626  if (fSptHitMap.count(sptid) != 0)
627  throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): hit already present!\n";
628  fSptHitMap[sptid] = hits;
629 
630  // Do a preliminary scan of hits.
631  // Determine weight given to hits in each view.
632 
633  unsigned int nplanes = geom->TPC(geo::TPCID(cstat0, tpc0)).Nplanes();
634  std::vector<int> numhits(nplanes, 0);
635  std::vector<double> weight(nplanes, 0.);
636 
637  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
638  ++ihit) {
639 
640  const recob::Hit& hit = **ihit;
641  geo::WireID hitWireID = hit.WireID();
642  // kept as assertions for performance reasons
643  assert(hitWireID.Cryostat == cstat0);
644  assert(hitWireID.TPC == tpc0);
645  assert(hitWireID.Plane < nplanes);
646  ++numhits[hitWireID.Plane];
647  }
648 
649  for (unsigned int plane = 0; plane < nplanes; ++plane) {
650  double np = numhits[plane];
651  if (np > 0.) weight[plane] = 1. / (np * np * np);
652  }
653 
654  // Calculate position and error matrix.
655 
656  double xyz[3] = {0., 0., 0.};
657  double errxyz[6] = {0., 0., 0., 0., 0., 0.};
658 
659  // Calculate x using drift times.
660  // Loop over all hits and calculate the weighted average drift time.
661  // Also calculate time variance and chisquare.
662 
663  double sumt2w = 0.;
664  double sumtw = 0.;
665  double sumw = 0.;
666 
667  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
668  ++ihit) {
669 
670  const recob::Hit& hit = **ihit;
671  geo::WireID hitWireID = hit.WireID();
672 
673  // Correct time for trigger offset and view-dependent time offsets.
674 
675  double t0 = detProp.GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
676  double t = hit.PeakTime() - t0;
677  double et = hit.SigmaPeakTime();
678  double w = weight[hitWireID.Plane] / (et * et);
679 
680  sumt2w += w * t * t;
681  sumtw += w * t;
682  sumw += w;
683  }
684 
685  double drift_time = 0.;
686  double var_time = 0.;
687  double chisq = 0.;
688  if (sumw != 0.) {
689  drift_time = sumtw / sumw;
690  var_time = sumt2w / sumw - drift_time * drift_time;
691  if (var_time < 0.) var_time = 0.;
692  chisq = sumt2w - sumtw * drift_time;
693  if (chisq < 0.) chisq = 0.;
694  }
695  xyz[0] = drift_time * timePitch;
696  errxyz[0] = var_time * timePitch * timePitch;
697 
698  // Calculate y, z using wires (need at least two hits).
699 
700  if (nhits >= 2) {
701 
702  // Calculate y and z by chisquare minimization of wire coordinates.
703 
704  double sus = 0.; // sum w_i u_i sin_th_i
705  double suc = 0.; // sum w_i u_i cos_th_i
706  double sc2 = 0.; // sum w_i cos2_th_i
707  double ss2 = 0.; // sum w_i sin2_th_i
708  double ssc = 0.; // sum w_i sin_th_i cos_th_i
709 
710  // Loop over points.
711 
712  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
713  ++ihit) {
714 
715  const recob::Hit& hit = **ihit;
716  geo::WireID hitWireID = hit.WireID();
717  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
718 
719  // Calculate angle and wire coordinate in this view.
720 
721  double const hl = wgeom.HalfL();
722  auto const cen = wgeom.GetCenter();
723  auto const cen1 = wgeom.GetEnd();
724  double s = (cen1.Y() - cen.Y()) / hl;
725  double c = (cen1.Z() - cen.Z()) / hl;
726  double u = cen.Z() * s - cen.Y() * c;
727  double eu = geom->WirePitch(hitWireID.asPlaneID()) / std::sqrt(12.);
728  double w = weight[hitWireID.Plane] / (eu * eu);
729 
730  // Summations
731 
732  sus += w * u * s;
733  suc += w * u * c;
734  sc2 += w * c * c;
735  ss2 += w * s * s;
736  ssc += w * s * c;
737  }
738 
739  // Calculate y,z
740 
741  double denom = sc2 * ss2 - ssc * ssc;
742  if (denom != 0.) {
743  xyz[1] = (-suc * ss2 + sus * ssc) / denom;
744  xyz[2] = (sus * sc2 - suc * ssc) / denom;
745  errxyz[2] = ss2 / denom;
746  errxyz[4] = ssc / denom;
747  errxyz[5] = sc2 / denom;
748  }
749 
750  // Make space point.
751 
752  recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
753  sptv.push_back(spt);
754  }
755  }
code to link reconstructed objects back to the MC truth information
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:114
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
Float_t ssc
Definition: plot.C:23
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
iterator begin()
Definition: PtrVector.h:217
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Point_t GetEnd() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:231
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
void hits()
Definition: readHits.C:15
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
iterator end()
Definition: PtrVector.h:231
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:172
reference front()
Definition: PtrVector.h:373
double weight
Definition: plottest35.C:25
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:224
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void trkf::SpacePointAlg::fillSpacePoint ( detinfo::DetectorPropertiesData const &  detProp,
const art::PtrVector< recob::Hit > &  hits,
std::vector< recob::SpacePoint > &  sptv,
int  sptid 
) const

Definition at line 403 of file SpacePointAlg.cxx.

References geo::PlaneID::asPlaneID(), art::PtrVector< T >::begin(), geo::CryostatID::Cryostat, art::PtrVector< T >::end(), fSptHitMap, geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), detinfo::DetectorPropertiesData::GetXTicksCoefficient(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), hits(), recob::Hit::PeakTime(), geo::PlaneID::Plane, recob::Hit::SigmaPeakTime(), art::PtrVector< T >::size(), ssc, geo::TPCID::TPC, w, recob::Hit::WireID(), geo::GeometryCore::WireIDToWireGeo(), and geo::GeometryCore::WirePitch().

Referenced by fillSpacePoints(), and makeSpacePoints().

407  {
409 
410  double timePitch = detProp.GetXTicksCoefficient();
411 
412  int nhits = hits.size();
413 
414  // Remember associated hits internally.
415 
416  if (fSptHitMap.find(sptid) != fSptHitMap.end())
417  throw cet::exception("SpacePointAlg") << "fillSpacePoint(): hit already present!\n";
418  fSptHitMap[sptid] = hits;
419 
420  // Calculate position and error matrix.
421 
422  double xyz[3] = {0., 0., 0.};
423  double errxyz[6] = {0., 0., 0., 0., 0., 0.};
424 
425  // Calculate x using drift times.
426  // Loop over all hits and calculate the weighted average drift time.
427  // Also calculate time variance and chisquare.
428 
429  double sumt2w = 0.;
430  double sumtw = 0.;
431  double sumw = 0.;
432 
433  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
434  ++ihit) {
435 
436  const recob::Hit& hit = **ihit;
437  geo::WireID hitWireID = hit.WireID();
438 
439  // Correct time for trigger offset and view-dependent time offsets.
440 
441  double t0 = detProp.GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
442  double t = hit.PeakTime() - t0;
443  double et = hit.SigmaPeakTime();
444  double w = 1. / (et * et);
445 
446  sumt2w += w * t * t;
447  sumtw += w * t;
448  sumw += w;
449  }
450 
451  double drift_time = 0.;
452  double var_time = 0.;
453  double chisq = 0.;
454  if (sumw != 0.) {
455  drift_time = sumtw / sumw;
456  //var_time = sumt2w / sumw - drift_time * drift_time;
457  var_time = 1. / sumw;
458  if (var_time < 0.) var_time = 0.;
459  chisq = sumt2w - sumtw * drift_time;
460  if (chisq < 0.) chisq = 0.;
461  }
462  xyz[0] = drift_time * timePitch;
463  errxyz[0] = var_time * timePitch * timePitch;
464 
465  // Calculate y, z using wires (need at least two hits).
466 
467  if (nhits >= 2) {
468 
469  // Calculate y and z by chisquare minimization of wire coordinates.
470 
471  double sus = 0.; // sum w_i u_i sin_th_i
472  double suc = 0.; // sum w_i u_i cos_th_i
473  double sc2 = 0.; // sum w_i cos2_th_i
474  double ss2 = 0.; // sum w_i sin2_th_i
475  double ssc = 0.; // sum w_i sin_th_i cos_th_i
476 
477  // Loop over points.
478 
479  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
480  ++ihit) {
481 
482  const recob::Hit& hit = **ihit;
483  geo::WireID hitWireID = hit.WireID();
484  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
485 
486  // Calculate angle and wire coordinate in this view.
487 
488  double const hl = wgeom.HalfL();
489  auto const cen = wgeom.GetCenter();
490  auto const cen1 = wgeom.GetEnd();
491  double s = (cen1.Y() - cen.Y()) / hl;
492  double c = (cen1.Z() - cen.Z()) / hl;
493  double u = cen.Z() * s - cen.Y() * c;
494  double eu = geom->WirePitch(hitWireID.asPlaneID()) / std::sqrt(12.);
495  double w = 1. / (eu * eu);
496 
497  // Summations
498 
499  sus += w * u * s;
500  suc += w * u * c;
501  sc2 += w * c * c;
502  ss2 += w * s * s;
503  ssc += w * s * c;
504  }
505 
506  // Calculate y,z
507 
508  double denom = sc2 * ss2 - ssc * ssc;
509  if (denom != 0.) {
510  xyz[1] = (-suc * ss2 + sus * ssc) / denom;
511  xyz[2] = (sus * sc2 - suc * ssc) / denom;
512  errxyz[2] = ss2 / denom;
513  errxyz[4] = ssc / denom;
514  errxyz[5] = sc2 / denom;
515  }
516 
517  // Make space point.
518 
519  recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
520  sptv.push_back(spt);
521  }
522  return;
523  }
code to link reconstructed objects back to the MC truth information
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:114
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
Float_t ssc
Definition: plot.C:23
iterator begin()
Definition: PtrVector.h:217
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Point_t GetEnd() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:231
void hits()
Definition: readHits.C:15
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
iterator end()
Definition: PtrVector.h:231
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:172
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:224
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void trkf::SpacePointAlg::fillSpacePoints ( detinfo::DetectorPropertiesData const &  detProp,
std::vector< recob::SpacePoint > &  spts,
std::multimap< double, KHitTrack > const &  trackMap 
) const

Fill a collection of space points.

Fill a collection of space points.

Arguments:

spts - Collection of space points to fill. sptalg - Space point algorithm object.

This method uses the hits contained in this track to construct space points.

This method does not have any knowledge of what constitutes a good space point, except that Hits are required to be consecutive when sorted by path distance, and space points are required to pass compatibility tests used by the space point algorithm object. This method will make space points from either two or three Hits (even for three-plane detectors), if the space point algorithm is configured to allow it.

Definition at line 543 of file SpacePointAlg.cxx.

References art::PtrVector< T >::clear(), compatible(), fillSpacePoint(), trkf::KHitWireX::getHit(), trkf::KHitTrack::getHit(), hits(), numHitMap(), art::PtrVector< T >::push_back(), art::PtrVector< T >::size(), and track.

Referenced by trkf::Track3DKalmanHit::createOutputs(), and trkf::TrackKalmanCheater::produce().

546  {
547  // Loop over KHitTracks.
548 
550  art::PtrVector<recob::Hit> compatible_hits;
551  for (std::multimap<double, KHitTrack>::const_iterator it = trackMap.begin();
552  it != trackMap.end();
553  ++it) {
554  const KHitTrack& track = (*it).second;
555 
556  // Extrack Hit from track.
557 
558  const std::shared_ptr<const KHitBase>& hit = track.getHit();
559  const KHitWireX* phit = dynamic_cast<const KHitWireX*>(&*hit);
560  if (phit != 0) {
561  const art::Ptr<recob::Hit> prhit = phit->getHit();
562 
563  // Test this hit for compatibility.
564 
565  hits.push_back(prhit);
566  bool ok = this->compatible(detProp, hits);
567  if (!ok) {
568 
569  // The new hit is not compatible. Make a space point out of
570  // the last known compatible hits, provided there are at least
571  // two.
572 
573  if (compatible_hits.size() >= 2) {
574  this->fillSpacePoint(detProp, compatible_hits, spts, this->numHitMap());
575  compatible_hits.clear();
576  }
577 
578  // Forget about any previous hits.
579 
580  hits.clear();
581  hits.push_back(prhit);
582  }
583 
584  // Update the list of known compatible hits.
585 
586  compatible_hits = hits;
587  }
588  }
589 
590  // Maybe make one final space point.
591 
592  if (compatible_hits.size() >= 2) {
593  this->fillSpacePoint(detProp, compatible_hits, spts, this->numHitMap());
594  }
595  }
intermediate_table::const_iterator const_iterator
void hits()
Definition: readHits.C:15
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
bool compatible(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
int numHitMap() const
void fillSpacePoint(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
Float_t track
Definition: plot.C:35
void clear()
Definition: PtrVector.h:533
bool trkf::SpacePointAlg::filter ( ) const
inlinenoexcept

Definition at line 86 of file SpacePointAlg.h.

86 { return fFilter; }
bool fFilter
Filter flag.
const art::PtrVector< recob::Hit > & trkf::SpacePointAlg::getAssociatedHits ( const recob::SpacePoint spt) const

Definition at line 1435 of file SpacePointAlg.cxx.

References fSptHitMap, and recob::SpacePoint::ID().

Referenced by trkf::SpacePointAna::analyze(), trkf::Track3DKalmanHit::createOutputs(), trkf::SeedFinderAlgorithm::FindSeedAtEnd(), trkf::SeedFinderAlgorithm::FindSeeds(), makeSpacePoints(), trkf::TCTrack::produce(), trkf::SpacePointCheater::produce(), trkf::SpacePointFinder::produce(), and trkf::TrackKalmanCheater::produce().

1437  {
1438  // It is an error if no hits are associated with this space point (throw exception).
1439 
1440  std::map<int, art::PtrVector<recob::Hit>>::const_iterator it = fSptHitMap.find(spt.ID());
1441  if (it == fSptHitMap.end()) {
1442  mf::LogWarning("SpacePointAlg")
1443  << "Looking for ID " << spt.ID() << " from " << fSptHitMap.size() << std::endl;
1444  throw cet::exception("SpacePointAlg") << "No Hits associated with space point.\n";
1445  }
1446  return (*it).second;
1447  }
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
ID_t ID() const
Definition: SpacePoint.h:74
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void trkf::SpacePointAlg::makeMCTruthSpacePoints ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
const art::PtrVector< recob::Hit > &  hits,
std::vector< recob::SpacePoint > &  spts 
) const

Definition at line 773 of file SpacePointAlg.cxx.

References makeSpacePoints().

Referenced by trkf::SpacePointAna::analyze(), and trkf::SpacePointCheater::produce().

777  {
778  makeSpacePoints(clockData, detProp, hits, spts, true);
779  }
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
void trkf::SpacePointAlg::makeSpacePoints ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
const art::PtrVector< recob::Hit > &  hits,
std::vector< recob::SpacePoint > &  spts 
) const

Definition at line 761 of file SpacePointAlg.cxx.

Referenced by trkf::SpacePointAna::analyze(), trkf::SeedFinderAlgorithm::FindSeeds(), trkf::FeatureTracker::Get3DFeaturePoints(), makeMCTruthSpacePoints(), cluster::ClusterMatchAlg::Match_SpacePoint(), trkf::TCTrack::produce(), and trkf::SpacePointFinder::produce().

765  {
766  makeSpacePoints(clockData, detProp, hits, spts, false);
767  }
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
void trkf::SpacePointAlg::makeSpacePoints ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
const art::PtrVector< recob::Hit > &  hits,
std::vector< recob::SpacePoint > &  spts,
bool  useMC 
) const
private
Todo:
Why are we still checking on whether this is MC or not?
Todo:
Such checks should not be in reconstruction code.

Definition at line 785 of file SpacePointAlg.cxx.

References util::abs(), util::begin(), art::PtrVector< T >::begin(), c1, c2, recob::SpacePoint::Chisq(), art::PtrVector< T >::clear(), compatible(), geo::CryostatID::Cryostat, geo::GeometryCore::Cryostat(), tca::debug, trkf::SpacePointAlg::HitMCInfo::dist2, util::empty(), art::PtrVector< T >::end(), art::PtrVector< T >::erase(), fEnableU, fEnableV, fEnableW, fFilter, fHitMCMap, fillComplexSpacePoint(), fillSpacePoint(), fMaxDT, fMaxS, fMerge, fMinViews, fPreferColl, fSptHitMap, getAssociatedHits(), geo::WireGeo::GetEnd(), geo::WireGeo::GetStart(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::WireGeo::HalfL(), cheat::BackTrackerService::HitToAvgSimIDEs(), mf::isDebugEnabled(), geo::GeometryCore::Iterate(), geo::kCollection, geo::kU, geo::kV, geo::kZ, geo::GeometryCore::Ncryostats(), geo::TPCGeo::Nplanes(), geo::CryostatGeo::NTPC(), trkf::SpacePointAlg::HitMCInfo::pchit, recob::Hit::PeakTime(), geo::PlaneID::Plane, geo::GeometryCore::Plane(), art::PtrVector< T >::push_back(), art::PtrVector< T >::reserve(), geo::GeometryCore::SignalType(), cheat::BackTrackerService::SimIDEsToXYZ(), art::PtrVector< T >::size(), util::size(), t1, t2, geo::CryostatGeo::TPC(), geo::TPCID::TPC, geo::GeometryCore::TPC(), trkf::SpacePointAlg::HitMCInfo::trackIDs, update(), recob::Hit::View(), geo::PlaneGeo::Wire(), geo::WireID::Wire, recob::Hit::WireID(), geo::GeometryCore::WireIDToWireGeo(), geo::GeometryCore::WirePitch(), x, and trkf::SpacePointAlg::HitMCInfo::xyz.

790  {
792 
793  fSptHitMap.clear();
794 
795  // Print diagnostic information.
796 
797  update(detProp);
798 
799  // First make sure result vector is empty.
800 
801  spts.erase(spts.begin(), spts.end());
802 
803  // Statistics.
804 
805  int n2 = 0; // Number of two-hit space points.
806  int n3 = 0; // Number of three-hit space points.
807  int n2filt = 0; // Number of two-hit space points after filtering/merging.
808  int n3filt = 0; // Number of three-hit space pointe after filtering/merging.
809 
810  // Sort hits into maps indexed by [cryostat][tpc][plane][wire].
811  // If using mc information, also generate maps of sim::IDEs and mc
812  // position indexed by hit.
813 
814  std::vector<std::vector<std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>>> hitmap;
815  fHitMCMap.clear();
816 
817  unsigned int ncstat = geom->Ncryostats();
818  hitmap.resize(ncstat);
819  for (auto const& cryoid : geom->Iterate<geo::CryostatID>()) {
820  unsigned int cstat = cryoid.Cryostat;
821  unsigned int ntpc = geom->Cryostat(cryoid).NTPC();
822  hitmap[cstat].resize(ntpc);
823  for (unsigned int tpc = 0; tpc < ntpc; ++tpc) {
824  int nplane = geom->Cryostat(cryoid).TPC(tpc).Nplanes();
825  hitmap[cstat][tpc].resize(nplane);
826  }
827  }
828 
829  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
830  ++ihit) {
831  const art::Ptr<recob::Hit>& phit = *ihit;
832  geo::View_t view = phit->View();
833  if ((view == geo::kU && fEnableU) || (view == geo::kV && fEnableV) ||
834  (view == geo::kZ && fEnableW)) {
835  geo::WireID phitWireID = phit->WireID();
836  hitmap[phitWireID.Cryostat][phitWireID.TPC][phitWireID.Plane].insert(
837  std::make_pair(phitWireID.Wire, phit));
838  }
839  }
840 
841  // Fill mc information, including IDEs and closest neighbors
842  // of each hit.
845  if (useMC) {
847 
848  // First loop over hits and fill track ids and mc position.
849  for (auto const& id : geom->Iterate<geo::PlaneID>()) {
850  auto const [cstat, tpc, plane] = std::make_tuple(id.Cryostat, id.TPC, id.Plane);
851  int nplane = geom->TPC(id).Nplanes();
852  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit =
853  hitmap[cstat][tpc][plane].begin();
854  ihit != hitmap[cstat][tpc][plane].end();
855  ++ihit) {
856  const art::Ptr<recob::Hit>& phit = ihit->second;
857  const recob::Hit& hit = *phit;
858  HitMCInfo& mcinfo = fHitMCMap[&hit]; // Default HitMCInfo.
859 
860  // Fill default nearest neighbor information (i.e. none).
861 
862  mcinfo.pchit.resize(nplane, 0);
863  mcinfo.dist2.resize(nplane, 1.e20);
864 
865  // Get sim::IDEs for this hit.
866 
867  std::vector<sim::IDE> ides = bt_serv->HitToAvgSimIDEs(clockData, phit);
868 
869  // Get sorted track ids. for this hit.
870 
871  mcinfo.trackIDs.reserve(ides.size());
872  for (std::vector<sim::IDE>::const_iterator i = ides.begin(); i != ides.end(); ++i)
873  mcinfo.trackIDs.push_back(i->trackID);
874  sort(mcinfo.trackIDs.begin(), mcinfo.trackIDs.end());
875 
876  // Get position of ionization for this hit.
877 
878  try {
879  mcinfo.xyz = bt_serv->SimIDEsToXYZ(ides);
880  }
881  catch (cet::exception& x) {
882  mcinfo.xyz.clear();
883  }
884  } // end loop over ihit
885  } // end loop oer planes
886 
887  // Loop over hits again and fill nearest neighbor information for real.
888  for (auto const& id : geom->Iterate<geo::PlaneID>()) {
889  auto const [cstat, tpc, plane] = std::make_tuple(id.Cryostat, id.TPC, id.Plane);
890  int nplane = geom->TPC(id).Nplanes();
891  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit =
892  hitmap[cstat][tpc][plane].begin();
893  ihit != hitmap[cstat][tpc][plane].end();
894  ++ihit) {
895  const art::Ptr<recob::Hit>& phit = ihit->second;
896  const recob::Hit& hit = *phit;
897  HitMCInfo& mcinfo = fHitMCMap[&hit];
898  if (mcinfo.xyz.size() != 0) {
899  assert(mcinfo.xyz.size() == 3);
900 
901  // Fill nearest neighbor information for this hit.
902 
903  for (int plane2 = 0; plane2 < nplane; ++plane2) {
904  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator jhit =
905  hitmap[cstat][tpc][plane2].begin();
906  jhit != hitmap[cstat][tpc][plane2].end();
907  ++jhit) {
908  const art::Ptr<recob::Hit>& phit2 = jhit->second;
909  const recob::Hit& hit2 = *phit2;
910  const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
911 
912  if (mcinfo2.xyz.size() != 0) {
913  assert(mcinfo2.xyz.size() == 3);
914  double dx = mcinfo.xyz[0] - mcinfo2.xyz[0];
915  double dy = mcinfo.xyz[1] - mcinfo2.xyz[1];
916  double dz = mcinfo.xyz[2] - mcinfo2.xyz[2];
917  double dist2 = dx * dx + dy * dy + dz * dz;
918  if (dist2 < mcinfo.dist2[plane2]) {
919  mcinfo.dist2[plane2] = dist2;
920  mcinfo.pchit[plane2] = &hit2;
921  }
922  } // end if mcinfo2.xyz valid
923  } // end loop over jhit
924  } // end loop over plane2
925  } // end if mcinfo.xyz valid.
926  } // end loop over ihit
927  } // end loop over planes
928  } // end if MC
929 
930  // use mf::LogDebug instead of MF_LOG_DEBUG because we reuse it in many lines
931  // insertions are protected by mf::isDebugEnabled()
932  mf::LogDebug debug("SpacePointAlg");
933  if (mf::isDebugEnabled()) {
934  debug << "Total hits = " << hits.size() << "\n\n";
935  for (auto const& id : geom->Iterate<geo::TPCID>()) {
936  auto const [cstat, tpc] = std::make_tuple(id.Cryostat, id.TPC);
937  int nplane = hitmap[cstat][tpc].size();
938  for (int plane = 0; plane < nplane; ++plane) {
939  debug << "TPC, Plane: " << tpc << ", " << plane
940  << ", hits = " << hitmap[cstat][tpc][plane].size() << "\n";
941  }
942  } // end loop over TPCs
943  } // if debug
944 
945  // Make empty multimap from hit pointer on preferred
946  // (most-populated or collection) plane to space points that
947  // include that hit (used for sorting, filtering, and
948  // merging).
949 
950  typedef const recob::Hit* sptkey_type;
951  std::multimap<sptkey_type, recob::SpacePoint> sptmap;
952  std::set<sptkey_type> sptkeys; // Keys of multimap.
953 
954  // Loop over TPCs.
955  for (auto const& tpcid : geom->Iterate<geo::TPCID>()) {
956 
957  auto const [cstat, tpc] = std::make_tuple(tpcid.Cryostat, tpcid.TPC);
958  // Sort maps in increasing order of number of hits.
959  // This is so that we can do the outer loops over hits
960  // over the views with fewer hits.
961  //
962  // If config parameter PreferColl is true, treat the colleciton
963  // plane as if it had the most hits, regardless of how many
964  // hits it actually has. This will force space points to be
965  // filtered and merged with respect to the collection plane
966  // wires. It will also force space points to be sorted by
967  // collection plane wire.
968 
969  int nplane = hitmap[cstat][tpc].size();
970  std::vector<int> index(nplane);
971 
972  for (int i = 0; i < nplane; ++i)
973  index[i] = i;
974 
975  for (int i = 0; i < nplane - 1; ++i) {
976 
977  for (int j = i + 1; j < nplane; ++j) {
978  bool icoll =
979  fPreferColl && geom->SignalType(geo::PlaneID(tpcid, index[i])) == geo::kCollection;
980  bool jcoll =
981  fPreferColl && geom->SignalType(geo::PlaneID(tpcid, index[j])) == geo::kCollection;
982  if ((hitmap[cstat][tpc][index[i]].size() > hitmap[cstat][tpc][index[j]].size() &&
983  !jcoll) ||
984  icoll) {
985  int temp = index[i];
986  index[i] = index[j];
987  index[j] = temp;
988  }
989  }
990  } // end loop over i
991 
992  // how many views with hits?
993  // This will allow for the special case where we might have only 2 planes of information and
994  // still want space points even if a three plane TPC
995  std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>& hitsByPlaneVec =
996  hitmap[cstat][tpc];
997  int nViewsWithHits(0);
998 
999  for (int i = 0; i < nplane; i++) {
1000  if (hitsByPlaneVec[index[i]].size() > 0) nViewsWithHits++;
1001  }
1002 
1003  // If two-view space points are allowed, make a double loop
1004  // over hits and produce space points for compatible hit-pairs.
1005 
1006  if ((nViewsWithHits == 2 || nplane == 2) && fMinViews <= 2) {
1007 
1008  // Loop over pairs of views.
1009  for (int i = 0; i < nplane - 1; ++i) {
1010  unsigned int plane1 = index[i];
1011 
1012  if (hitmap[cstat][tpc][plane1].empty()) continue;
1013 
1014  for (int j = i + 1; j < nplane; ++j) {
1015  unsigned int plane2 = index[j];
1016 
1017  if (hitmap[cstat][tpc][plane2].empty()) continue;
1018 
1019  // Get angle, pitch, and offset of plane2 wires.
1020  geo::PlaneID const plane2_id{tpcid, plane2};
1021  const geo::WireGeo& wgeo2 = geom->Plane(plane2_id).Wire(0);
1022  double const hl2 = wgeo2.HalfL();
1023  auto const xyz21 = wgeo2.GetStart();
1024  auto const xyz22 = wgeo2.GetEnd();
1025  double s2 = (xyz22.Y() - xyz21.Y()) / (2. * hl2);
1026  double c2 = (xyz22.Z() - xyz21.Z()) / (2. * hl2);
1027  double dist2 = -xyz21.Y() * c2 + xyz21.Z() * s2;
1028  double pitch2 = geom->WirePitch(plane2_id);
1029 
1030  if (!fPreferColl &&
1031  hitmap[cstat][tpc][plane1].size() > hitmap[cstat][tpc][plane2].size())
1032  throw cet::exception("SpacePointAlg")
1033  << "makeSpacePoints(): hitmaps with incompatible size\n";
1034 
1035  // Loop over pairs of hits.
1036 
1038  hitvec.reserve(2);
1039 
1040  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit1 =
1041  hitmap[cstat][tpc][plane1].begin();
1042  ihit1 != hitmap[cstat][tpc][plane1].end();
1043  ++ihit1) {
1044 
1045  const art::Ptr<recob::Hit>& phit1 = ihit1->second;
1046  geo::WireID phit1WireID = phit1->WireID();
1047  const geo::WireGeo& wgeo = geom->WireIDToWireGeo(phit1WireID);
1048 
1049  // Get endpoint coordinates of this wire.
1050  // (kept as assertions for performance reasons)
1051  assert(phit1WireID.Cryostat == cstat);
1052  assert(phit1WireID.TPC == tpc);
1053  assert(phit1WireID.Plane == plane1);
1054  auto const xyz1 = wgeo.GetStart();
1055  auto const xyz2 = wgeo.GetEnd();
1056 
1057  // Find the plane2 wire numbers corresponding to the endpoints.
1058 
1059  double wire21 = (-xyz1.Y() * c2 + xyz1.Z() * s2 - dist2) / pitch2;
1060  double wire22 = (-xyz2.Y() * c2 + xyz2.Z() * s2 - dist2) / pitch2;
1061 
1062  int wmin = std::max(0., std::min(wire21, wire22));
1063  int wmax = std::max(0., std::max(wire21, wire22) + 1.);
1064 
1065  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1066  ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
1067  ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
1068 
1069  for (; ihit2 != ihit2end; ++ihit2) {
1070 
1071  const art::Ptr<recob::Hit>& phit2 = ihit2->second;
1072 
1073  // Check current pair of hits for compatibility.
1074  // By construction, hits should always have compatible views
1075  // and times, but may not have compatible mc information.
1076 
1077  hitvec.clear();
1078  hitvec.push_back(phit1);
1079  hitvec.push_back(phit2);
1080  bool ok = compatible(detProp, hitvec, useMC);
1081  if (ok) {
1082 
1083  // Add a space point.
1084 
1085  ++n2;
1086 
1087  // make a dummy vector of recob::SpacePoints
1088  // as we are filtering or merging and don't want to
1089  // add the created SpacePoint to the final collection just yet
1090  // This dummy vector will hold just one recob::SpacePoint,
1091  // which will go into the multimap and then the vector
1092  // will go out of scope.
1093 
1094  std::vector<recob::SpacePoint> sptv;
1095  fillSpacePoint(detProp, hitvec, sptv, sptmap.size());
1096  sptkey_type key = &*phit2;
1097  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
1098  sptkeys.insert(key);
1099  }
1100  }
1101  }
1102  }
1103  }
1104  } // end if fMinViews <= 2
1105 
1106  // If three-view space points are allowed, make a triple loop
1107  // over hits and produce space points for compatible triplets.
1108 
1109  if (nplane >= 3 && fMinViews <= 3) {
1110 
1111  // Loop over triplets of hits.
1112 
1114  hitvec.reserve(3);
1115 
1116  unsigned int plane1 = index[0];
1117  unsigned int plane2 = index[1];
1118  unsigned int plane3 = index[2];
1119 
1120  // Get angle, pitch, and offset of plane1 wires.
1121 
1122  geo::PlaneID const plane1_id{tpcid, plane1};
1123  const geo::WireGeo& wgeo1 = geom->Plane(plane1_id).Wire(0);
1124  double const hl1 = wgeo1.HalfL();
1125  auto const xyz11 = wgeo1.GetStart();
1126  auto const xyz12 = wgeo1.GetEnd();
1127  double s1 = (xyz12.Y() - xyz11.Y()) / (2. * hl1);
1128  double c1 = (xyz12.Z() - xyz11.Z()) / (2. * hl1);
1129  double dist1 = -xyz11.Y() * c1 + xyz11.Z() * s1;
1130  double pitch1 = geom->WirePitch(plane1_id);
1131  const double TicksOffset1 = detProp.GetXTicksOffset(plane1_id);
1132 
1133  // Get angle, pitch, and offset of plane2 wires.
1134 
1135  geo::PlaneID const plane2_id{tpcid, plane2};
1136  const geo::WireGeo& wgeo2 = geom->Plane(plane2_id).Wire(0);
1137  double const hl2 = wgeo2.HalfL();
1138  auto const xyz21 = wgeo2.GetStart();
1139  auto const xyz22 = wgeo2.GetEnd();
1140  double s2 = (xyz22.Y() - xyz21.Y()) / (2. * hl2);
1141  double c2 = (xyz22.Z() - xyz21.Z()) / (2. * hl2);
1142  double dist2 = -xyz21.Y() * c2 + xyz21.Z() * s2;
1143  double pitch2 = geom->WirePitch(plane2_id);
1144  const double TicksOffset2 = detProp.GetXTicksOffset(plane2_id);
1145 
1146  // Get angle, pitch, and offset of plane3 wires.
1147 
1148  geo::PlaneID const plane3_id{tpcid, plane3};
1149  const geo::WireGeo& wgeo3 = geom->Plane(plane3_id).Wire(0);
1150  double const hl3 = wgeo3.HalfL();
1151  auto const xyz31 = wgeo3.GetStart();
1152  auto const xyz32 = wgeo3.GetEnd();
1153  double s3 = (xyz32.Y() - xyz31.Y()) / (2. * hl3);
1154  double c3 = (xyz32.Z() - xyz31.Z()) / (2. * hl3);
1155  double dist3 = -xyz31.Y() * c3 + xyz31.Z() * s3;
1156  double pitch3 = geom->WirePitch(plane3_id);
1157  const double TicksOffset3 = detProp.GetXTicksOffset(plane3_id);
1158 
1159  // Get sine of angle differences.
1160 
1161  double s12 = s1 * c2 - s2 * c1; // sin(theta1 - theta2).
1162  double s23 = s2 * c3 - s3 * c2; // sin(theta2 - theta3).
1163  double s31 = s3 * c1 - s1 * c3; // sin(theta3 - theta1).
1164 
1165  // Loop over hits in plane1.
1166 
1167  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1168  ihit1 = hitmap[cstat][tpc][plane1].begin(),
1169  ihit1end = hitmap[cstat][tpc][plane1].end();
1170  for (; ihit1 != ihit1end; ++ihit1) {
1171 
1172  unsigned int wire1 = ihit1->first;
1173  const art::Ptr<recob::Hit>& phit1 = ihit1->second;
1174  geo::WireID phit1WireID = phit1->WireID();
1175  const geo::WireGeo& wgeo = geom->WireIDToWireGeo(phit1WireID);
1176 
1177  // Get endpoint coordinates of this wire from plane1.
1178  // (kept as assertions for performance reasons)
1179  assert(phit1WireID.Cryostat == cstat);
1180  assert(phit1WireID.TPC == tpc);
1181  assert(phit1WireID.Plane == plane1);
1182  assert(phit1WireID.Wire == wire1);
1183  auto const xyz1 = wgeo.GetStart();
1184  auto const xyz2 = wgeo.GetEnd();
1185 
1186  // Get corrected time and oblique coordinate of first hit.
1187 
1188  double t1 = phit1->PeakTime() - TicksOffset1;
1189  double u1 = wire1 * pitch1 + dist1;
1190 
1191  // Find the plane2 wire numbers corresponding to the endpoints.
1192 
1193  double wire21 = (-xyz1.Y() * c2 + xyz1.Z() * s2 - dist2) / pitch2;
1194  double wire22 = (-xyz2.Y() * c2 + xyz2.Z() * s2 - dist2) / pitch2;
1195 
1196  int wmin = std::max(0., std::min(wire21, wire22));
1197  int wmax = std::max(0., std::max(wire21, wire22) + 1.);
1198 
1199  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1200  ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
1201  ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
1202 
1203  for (; ihit2 != ihit2end; ++ihit2) {
1204 
1205  int wire2 = ihit2->first;
1206  const art::Ptr<recob::Hit>& phit2 = ihit2->second;
1207 
1208  // Get corrected time of second hit.
1209 
1210  double t2 = phit2->PeakTime() - TicksOffset2;
1211 
1212  // Check maximum time difference with first hit.
1213 
1214  bool dt12ok = std::abs(t1 - t2) <= fMaxDT;
1215  if (dt12ok) {
1216 
1217  // Test first two hits for compatibility before looping
1218  // over third hit.
1219 
1220  hitvec.clear();
1221  hitvec.push_back(phit1);
1222  hitvec.push_back(phit2);
1223  bool h12ok = compatible(detProp, hitvec, useMC);
1224  if (h12ok) {
1225 
1226  // Get oblique coordinate of second hit.
1227 
1228  double u2 = wire2 * pitch2 + dist2;
1229 
1230  // Predict plane3 oblique coordinate and wire number.
1231 
1232  double u3pred = (-u1 * s23 - u2 * s31) / s12;
1233  double w3pred = (u3pred - dist3) / pitch3;
1234  double w3delta = std::abs(fMaxS / (s12 * pitch3));
1235  int w3min = std::max(0., std::ceil(w3pred - w3delta));
1236  int w3max = std::max(0., std::floor(w3pred + w3delta));
1237 
1238  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1239  ihit3 = hitmap[cstat][tpc][plane3].lower_bound(w3min),
1240  ihit3end = hitmap[cstat][tpc][plane3].upper_bound(w3max);
1241 
1242  for (; ihit3 != ihit3end; ++ihit3) {
1243 
1244  int wire3 = ihit3->first;
1245  const art::Ptr<recob::Hit>& phit3 = ihit3->second;
1246 
1247  // Get corrected time of third hit.
1248 
1249  double t3 = phit3->PeakTime() - TicksOffset3;
1250 
1251  // Check time difference of third hit compared to first two hits.
1252 
1253  bool dt123ok = std::abs(t1 - t3) <= fMaxDT && std::abs(t2 - t3) <= fMaxDT;
1254  if (dt123ok) {
1255 
1256  // Get oblique coordinate of third hit and check spatial separation.
1257 
1258  double u3 = wire3 * pitch3 + dist3;
1259  double S = s23 * u1 + s31 * u2 + s12 * u3;
1260  bool sok = std::abs(S) <= fMaxS;
1261  if (sok) {
1262 
1263  // Test triplet for compatibility.
1264 
1265  hitvec.clear();
1266  hitvec.push_back(phit1);
1267  hitvec.push_back(phit2);
1268  hitvec.push_back(phit3);
1269  bool h123ok = compatible(detProp, hitvec, useMC);
1270  if (h123ok) {
1271 
1272  // Add a space point.
1273 
1274  ++n3;
1275 
1276  // make a dummy vector of recob::SpacePoints
1277  // as we are filtering or merging and don't want to
1278  // add the created SpacePoint to the final collection just yet
1279  // This dummy vector will hold just one recob::SpacePoint,
1280  // which will go into the multimap and then the vector
1281  // will go out of scope.
1282 
1283  std::vector<recob::SpacePoint> sptv;
1284  fillSpacePoint(detProp, hitvec, sptv, sptmap.size() - 1);
1285  sptkey_type key = &*phit3;
1286  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
1287  sptkeys.insert(key);
1288  }
1289  }
1290  }
1291  }
1292  }
1293  }
1294  }
1295  }
1296  } // end if fMinViews <= 3
1297 
1298  // Do Filtering.
1299 
1300  if (fFilter) {
1301 
1302  // Transfer (some) space points from sptmap to spts.
1303 
1304  spts.reserve(spts.size() + sptkeys.size());
1305 
1306  // Loop over keys of space point map.
1307  // Space points that have the same key are candidates for filtering.
1308 
1309  for (std::set<sptkey_type>::const_iterator i = sptkeys.begin(); i != sptkeys.end(); ++i) {
1310  sptkey_type key = *i;
1311 
1312  // Loop over space points corresponding to the current key.
1313  // Choose the single best space point from among this group.
1314 
1315  double best_chisq = 0.;
1316  const recob::SpacePoint* best_spt = 0;
1317 
1319  sptmap.lower_bound(key);
1320  j != sptmap.upper_bound(key);
1321  ++j) {
1322  const recob::SpacePoint& spt = j->second;
1323  if (best_spt == 0 || spt.Chisq() < best_chisq) {
1324  best_spt = &spt;
1325  best_chisq = spt.Chisq();
1326  }
1327  }
1328 
1329  // Transfer best filtered space point to result vector.
1330 
1331  if (!best_spt)
1332  throw cet::exception("SpacePointAlg") << "makeSpacePoints(): no best point\n";
1333  spts.push_back(*best_spt);
1334  if (fMinViews <= 2)
1335  ++n2filt;
1336  else
1337  ++n3filt;
1338  }
1339  } // end if filtering
1340 
1341  // Do merging.
1342 
1343  else if (fMerge) {
1344 
1345  // Transfer merged space points from sptmap to spts.
1346 
1347  spts.reserve(spts.size() + sptkeys.size());
1348 
1349  // Loop over keys of space point map.
1350  // Space points that have the same key are candidates for merging.
1351 
1352  for (std::set<sptkey_type>::const_iterator i = sptkeys.begin(); i != sptkeys.end(); ++i) {
1353  sptkey_type key = *i;
1354 
1355  // Loop over space points corresponding to the current key.
1356  // Make a collection of hits that is the union of the hits
1357  // from each candidate space point.
1358 
1360  sptmap.lower_bound(key),
1361  jSPTend =
1362  sptmap.upper_bound(key);
1363 
1364  art::PtrVector<recob::Hit> merged_hits;
1365  for (; jSPT != jSPTend; ++jSPT) {
1366  const recob::SpacePoint& spt = jSPT->second;
1367 
1368  // Loop over hits from this space points.
1369  // Add each hit to the collection of all hits.
1370 
1371  const art::PtrVector<recob::Hit>& spt_hits = getAssociatedHits(spt);
1372  merged_hits.reserve(merged_hits.size() +
1373  spt_hits.size()); // better than nothing, but not ideal
1375  k != spt_hits.end();
1376  ++k) {
1377  const art::Ptr<recob::Hit>& hit = *k;
1378  merged_hits.push_back(hit);
1379  }
1380  }
1381 
1382  // Remove duplicates.
1383 
1384  std::sort(merged_hits.begin(), merged_hits.end());
1386  std::unique(merged_hits.begin(), merged_hits.end());
1387  merged_hits.erase(it, merged_hits.end());
1388 
1389  // Construct a complex space points using merged hits.
1390 
1391  fillComplexSpacePoint(detProp, merged_hits, spts, sptmap.size() + spts.size() - 1);
1392 
1393  if (fMinViews <= 2)
1394  ++n2filt;
1395  else
1396  ++n3filt;
1397  }
1398  } // end if merging
1399 
1400  // No filter, no merge.
1401 
1402  else {
1403 
1404  // Transfer all space points from sptmap to spts.
1405 
1406  spts.reserve(spts.size() + sptkeys.size());
1407 
1408  // Loop over space points.
1409 
1411  j != sptmap.end();
1412  ++j) {
1413  const recob::SpacePoint& spt = j->second;
1414  spts.push_back(spt);
1415  }
1416 
1417  // Update statistics.
1418 
1419  n2filt = n2;
1420  n3filt = n3;
1421  }
1422  } // end loop over tpcs
1423 
1424  if (mf::isDebugEnabled()) {
1425  debug << "\n2-hit space points = " << n2 << "\n"
1426  << "3-hit space points = " << n3 << "\n"
1427  << "2-hit filtered/merged space points = " << n2filt << "\n"
1428  << "3-hit filtered/merged space points = " << n3filt;
1429  } // if debug
1430  }
Float_t x
Definition: compare.C:6
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:114
void reserve(size_type n)
Definition: PtrVector.h:337
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
typename data_t::iterator iterator
Definition: PtrVector.h:54
std::map< const recob::Hit *, HitMCInfo > fHitMCMap
TTree * t1
Definition: plottest35.C:26
double fMaxDT
Maximum time difference between planes.
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:465
bool fFilter
Filter flag.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Planes which measure V.
Definition: geo_types.h:136
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
iterator begin()
Definition: PtrVector.h:217
void update(detinfo::DetectorPropertiesData const &detProp) const
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
bool fEnableW
Enable flag (W).
iterator erase(iterator position)
Definition: PtrVector.h:504
constexpr auto abs(T v)
Returns the absolute value of the argument.
Double32_t Chisq() const
Definition: SpacePoint.h:86
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
Planes which measure Z direction.
Definition: geo_types.h:138
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:276
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
intermediate_table::const_iterator const_iterator
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:430
CryostatGeo const & Cryostat(CryostatID const &cryoid=cryostat_zero) const
Returns the specified cryostat.
Point_t GetStart() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:226
Point_t GetEnd() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:231
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
Planes which measure U.
Definition: geo_types.h:135
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
TCanvas * c1
Definition: plotHisto.C:7
bool fEnableU
Enable flag (U).
TCanvas * c2
Definition: plot_hist.C:75
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
DebugStuff debug
Definition: DebugStruct.cxx:4
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
iterator end()
Definition: PtrVector.h:231
void fillComplexSpacePoint(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:175
bool isDebugEnabled()
double fMaxS
Maximum space separation between wires.
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
bool fPreferColl
Sort by collection wire.
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
TTree * t2
Definition: plottest35.C:36
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:172
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:84
std::vector< sim::IDE > HitToAvgSimIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides) const
bool compatible(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
int fMinViews
Mininum number of views per space point.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
bool fMerge
Merge flag.
recob::tracking::Plane Plane
Definition: TrackState.h:17
void fillSpacePoint(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
void clear()
Definition: PtrVector.h:533
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
Signal from collection planes.
Definition: geo_types.h:152
bool fEnableV
Enable flag (V).
double trkf::SpacePointAlg::maxDT ( ) const
inlinenoexcept

Definition at line 88 of file SpacePointAlg.h.

Referenced by cluster::ClusterMatchAlg::Match_SpacePoint().

88 { return fMaxDT; }
double fMaxDT
Maximum time difference between planes.
double trkf::SpacePointAlg::maxS ( ) const
inlinenoexcept

Definition at line 89 of file SpacePointAlg.h.

89 { return fMaxS; }
double fMaxS
Maximum space separation between wires.
bool trkf::SpacePointAlg::merge ( ) const
inlinenoexcept

Definition at line 87 of file SpacePointAlg.h.

Referenced by trkf::SpacePointAna::analyze(), and trkf::SpacePointAna::bookHistograms().

87 { return fMerge; }
bool fMerge
Merge flag.
int trkf::SpacePointAlg::minViews ( ) const
inlinenoexcept

Definition at line 90 of file SpacePointAlg.h.

Referenced by trkf::SpacePointCheater::produce(), and trkf::SpacePointFinder::produce().

90 { return fMinViews; }
int fMinViews
Mininum number of views per space point.
int trkf::SpacePointAlg::numHitMap ( ) const
inline

Definition at line 152 of file SpacePointAlg.h.

Referenced by fillSpacePoints().

152 { return fSptHitMap.size(); }
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
double trkf::SpacePointAlg::separation ( const art::PtrVector< recob::Hit > &  hits) const

Definition at line 172 of file SpacePointAlg.cxx.

References geo::CryostatID::Cryostat, larg4::dist(), geo::WireGeo::GetCenter(), geo::WireGeo::GetEnd(), geo::WireGeo::HalfL(), geo::PlaneID::Plane, art::PtrVector< T >::size(), geo::TPCID::TPC, recob::Hit::WireID(), and geo::GeometryCore::WireIDToWireGeo().

Referenced by trkf::SpacePointAna::analyze().

173  {
174  // Get geometry service.
175 
177 
178  // Trivial case - fewer than three hits.
179 
180  if (hits.size() < 3) return 0.;
181 
182  // Error case - more than three hits.
183 
184  if (hits.size() > 3) {
185  mf::LogError("SpacePointAlg") << "Method separation called with more than three htis.";
186  return 0.;
187  }
188 
189  // Got exactly three hits.
190 
191  // Calculate angles and distance of each hit from origin.
192 
193  double dist[3] = {0., 0., 0.};
194  double sinth[3] = {0., 0., 0.};
195  double costh[3] = {0., 0., 0.};
196  unsigned int cstats[3];
197  unsigned int tpcs[3];
198  unsigned int planes[3];
199 
200  for (int i = 0; i < 3; ++i) {
201 
202  // Get tpc, plane, wire.
203 
204  const recob::Hit& hit = *(hits[i]);
205  const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
206  cstats[i] = hit.WireID().Cryostat;
207  tpcs[i] = hit.WireID().TPC;
208  planes[i] = hit.WireID().Plane;
209 
210  // Check tpc and plane errors.
211 
212  for (int j = 0; j < i; ++j) {
213 
214  if (cstats[j] != hit.WireID().Cryostat) {
215  mf::LogError("SpacePointAlg")
216  << "Method separation called with hits from multiple cryostats..";
217  return 0.;
218  }
219 
220  if (tpcs[j] != hit.WireID().TPC) {
221  mf::LogError("SpacePointAlg")
222  << "Method separation called with hits from multiple tpcs..";
223  return 0.;
224  }
225 
226  if (planes[j] == hit.WireID().Plane) {
227  mf::LogError("SpacePointAlg")
228  << "Method separation called with hits from the same plane..";
229  return 0.;
230  }
231  }
232 
233  // Get angles and distance of wire.
234 
235  double const hl = wgeom.HalfL();
236  auto const xyz = wgeom.GetCenter();
237  auto const xyz1 = wgeom.GetEnd();
238  double s = (xyz1.Y() - xyz.Y()) / hl;
239  double c = (xyz1.Z() - xyz.Z()) / hl;
240  sinth[hit.WireID().Plane] = s;
241  costh[hit.WireID().Plane] = c;
242  dist[hit.WireID().Plane] = xyz.Z() * s - xyz.Y() * c;
243  }
244 
245  double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0] +
246  (sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1] +
247  (sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
248  return S;
249  }
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:114
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Point_t GetEnd() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:231
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:172
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
void trkf::SpacePointAlg::update ( detinfo::DetectorPropertiesData const &  detProp) const

Definition at line 77 of file SpacePointAlg.cxx.

References detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::GeometryCore::Iterate(), geo::kCollection, geo::kHorizontal, geo::kInduction, geo::kU, geo::kV, geo::kVertical, geo::kZ, geo::GeometryCore::SignalType(), and geo::GeometryCore::TPC().

Referenced by makeSpacePoints().

78  {
79  // Generate info report on first call only.
80 
81  static bool first = true;
82  bool report = first;
83  first = false;
84 
85  // Get services.
86 
88 
89  // Calculate and print geometry information.
90 
91  if (report) mf::LogInfo("SpacePointAlg") << "Updating geometry constants.\n";
92 
93  for (auto const& plane : geom->Iterate<geo::PlaneGeo>()) {
94 
95  // Fill view-dependent quantities.
96 
97  geo::View_t view = plane.View();
98  std::string viewname = "?";
99  if (view == geo::kU) { viewname = "U"; }
100  else if (view == geo::kV) {
101  viewname = "V";
102  }
103  else if (view == geo::kZ) {
104  viewname = "Z";
105  }
106  else
107  throw cet::exception("SpacePointAlg") << "Bad view = " << view << "\n";
108 
109  std::string sigtypename = "?";
110  geo::SigType_t sigtype = geom->SignalType(plane.ID());
111  if (sigtype == geo::kInduction)
112  sigtypename = "Induction";
113  else if (sigtype == geo::kCollection)
114  sigtypename = "Collection";
115  else
116  throw cet::exception("SpacePointAlg") << "Bad signal type = " << sigtype << "\n";
117 
118  std::string orientname = "?";
119  geo::Orient_t orient = plane.Orientation();
120  if (orient == geo::kVertical)
121  orientname = "Vertical";
122  else if (orient == geo::kHorizontal)
123  orientname = "Horizontal";
124  else
125  throw cet::exception("SpacePointAlg") << "Bad orientation = " << orient << "\n";
126 
127  if (report) {
128  auto const xyz = plane.GetCenter();
129  auto const& tpcgeom = geom->TPC(plane.ID());
130  mf::LogInfo("SpacePointAlg")
131  << '\n'
132  << plane.ID() << '\n'
133  << " View: " << viewname << "\n"
134  << " SignalType: " << sigtypename << "\n"
135  << " Orientation: " << orientname << "\n"
136  << " Plane location: " << xyz.X() << "\n"
137  << " Plane pitch: " << tpcgeom.Plane0Pitch(plane.ID().Plane) << "\n"
138  << " Wire angle: " << plane.Wire(0).ThetaZ() << "\n"
139  << " Wire pitch: " << tpcgeom.WirePitch() << "\n"
140  << " Time offset: " << detProp.GetXTicksOffset(plane.ID()) << "\n";
141  }
142 
143  if (orient != geo::kVertical)
144  throw cet::exception("SpacePointAlg") << "Horizontal wire geometry not implemented.\n";
145  } // end loop over planes
146  }
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Planes which measure V.
Definition: geo_types.h:136
enum geo::_plane_orient Orient_t
Enumerate the possible plane projections.
Planes which measure Z direction.
Definition: geo_types.h:138
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
Planes which measure U.
Definition: geo_types.h:135
Planes that are in the horizontal plane.
Definition: geo_types.h:146
Signal from induction planes.
Definition: geo_types.h:151
Planes that are in the vertical plane (e.g. ArgoNeuT).
Definition: geo_types.h:147
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:78
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Signal from collection planes.
Definition: geo_types.h:152

Member Data Documentation

bool trkf::SpacePointAlg::fEnableU
private

Enable flag (U).

Definition at line 168 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fEnableV
private

Enable flag (V).

Definition at line 169 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fEnableW
private

Enable flag (W).

Definition at line 170 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fFilter
private

Filter flag.

Definition at line 171 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

std::map<const recob::Hit*, HitMCInfo> trkf::SpacePointAlg::fHitMCMap
mutableprivate

Definition at line 186 of file SpacePointAlg.h.

Referenced by compatible(), and makeSpacePoints().

double trkf::SpacePointAlg::fMaxDT
private

Maximum time difference between planes.

Definition at line 165 of file SpacePointAlg.h.

Referenced by compatible(), makeSpacePoints(), and SpacePointAlg().

double trkf::SpacePointAlg::fMaxS
private

Maximum space separation between wires.

Definition at line 166 of file SpacePointAlg.h.

Referenced by compatible(), makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fMerge
private

Merge flag.

Definition at line 172 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

int trkf::SpacePointAlg::fMinViews
private

Mininum number of views per space point.

Definition at line 167 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

bool trkf::SpacePointAlg::fPreferColl
private

Sort by collection wire.

Definition at line 173 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and SpacePointAlg().

std::map<int, art::PtrVector<recob::Hit> > trkf::SpacePointAlg::fSptHitMap
mutableprivate
double trkf::SpacePointAlg::fTickOffsetU
private

Tick offset for plane U.

Definition at line 174 of file SpacePointAlg.h.

Referenced by correctedTime(), and SpacePointAlg().

double trkf::SpacePointAlg::fTickOffsetV
private

Tick offset for plane V.

Definition at line 175 of file SpacePointAlg.h.

Referenced by correctedTime(), and SpacePointAlg().

double trkf::SpacePointAlg::fTickOffsetW
private

Tick offset for plane W.

Definition at line 176 of file SpacePointAlg.h.

Referenced by correctedTime(), and SpacePointAlg().


The documentation for this class was generated from the following files: