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

#include "CosmicTrackerAlg.h"

Public Member Functions

 CosmicTrackerAlg (fhicl::ParameterSet const &pset)
 
void SPTReco (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
 

Public Attributes

std::vector< TVector3 > trajPos
 
std::vector< TVector3 > trajDir
 
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit
 
std::vector< TVector3 > trkPos
 
std::vector< TVector3 > trkDir
 

Private Member Functions

void TrackTrajectory (detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
 
void Track3D (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
 
void MakeSPT (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
 

Private Attributes

int fSPTAlg
 
bool fTrajOnly
 
double ftmatch
 tolerance for time matching (in ticks) More...
 
double fsmatch
 tolerance for distance matching (in cm) More...
 
TrackTrajectoryAlg fTrackTrajectoryAlg
 
std::vector< std::vector< std::vector< std::vector< double > > > > vw
 
std::vector< std::vector< std::vector< std::vector< double > > > > vt
 
std::vector< std::vector< std::vector< std::vector< unsigned int > > > > vtraj
 
art::ServiceHandle< geo::Geometry const > geom
 
const detinfo::LArPropertieslarprop
 

Detailed Description

Definition at line 36 of file CosmicTrackerAlg.h.

Constructor & Destructor Documentation

trkf::CosmicTrackerAlg::CosmicTrackerAlg ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 40 of file CosmicTrackerAlg.cxx.

References fhicl::ParameterSet::get().

41  {
42  fSPTAlg = pset.get<int>("SPTAlg", 0);
43  fTrajOnly = pset.get<bool>("TrajOnly", false);
44  ftmatch = pset.get<double>("TMatch");
45  fsmatch = pset.get<double>("SMatch");
46  }
double ftmatch
tolerance for time matching (in ticks)
double fsmatch
tolerance for distance matching (in cm)

Member Function Documentation

void trkf::CosmicTrackerAlg::MakeSPT ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> &  fHits 
)
private

Definition at line 428 of file CosmicTrackerAlg.cxx.

References geo::PlaneID::asPlaneID(), geo::CryostatID::Cryostat, detinfo::DetectorPropertiesData::DriftVelocity(), e, detinfo::DetectorPropertiesData::Efield(), geo::PlaneID::Plane, detinfo::sampling_rate(), detinfo::DetectorPropertiesData::Temperature(), geo::TPCID::TPC, geo::WireID::Wire, and geo::WireID::WireID().

431  {
432  larprop = lar::providerFrom<detinfo::LArPropertiesService>();
433 
434  double timetick = sampling_rate(clockData) * 1e-3; // time sample in us
435  double Efield_drift = detProp.Efield(0); // Electric Field in the drift region in kV/cm
436  double Temperature = detProp.Temperature(); // LAr Temperature in K
437  double driftvelocity =
438  detProp.DriftVelocity(Efield_drift, Temperature); //drift velocity in the drift region (cm/us)
439  double time_pitch = driftvelocity * timetick; //time sample (cm)
440 
441  for (size_t i = 0; i < fHits.size(); ++i) {
442  int ip1 = -1;
443  int ip2 = -1;
444  int ip2p = -1;
445  int ih1 = -1;
446  int ih2 = -1;
447  int ih2p = -1;
448  double mindis1 = 1e9;
449  double mindis2 = 1e9;
450  double mindis2p = 1e9;
451  geo::WireID wireid = fHits[i]->WireID();
452  unsigned int wire = wireid.Wire;
453  unsigned int plane = wireid.Plane;
454  unsigned int tpc = wireid.TPC;
455  unsigned int cstat = wireid.Cryostat;
456  double wire_pitch = geom->WirePitch(wireid.asPlaneID()); //wire pitch in cm
457  //find the two trajectory points that enclose the hit
458  //find the nearest trajectory point first
459  for (size_t j = 0; j < vw[cstat][tpc][plane].size(); ++j) {
460  double dis = std::hypot(wire_pitch * (wire - vw[cstat][tpc][plane][j]),
461  time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]));
462  if (dis < mindis1) {
463  mindis1 = dis;
464  ip1 = vtraj[cstat][tpc][plane][j];
465  ih1 = j;
466  }
467  }
468  //find the second nearest trajectory point, prefer the point on the other side
469  for (size_t j = 0; j < vw[cstat][tpc][plane].size(); ++j) {
470  if (int(j) == ih1) continue;
471  double dis = std::hypot(wire_pitch * (wire - vw[cstat][tpc][plane][j]),
472  time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]));
473  if (dis < mindis2) {
474  mindis2 = dis;
475  ip2 = vtraj[cstat][tpc][plane][j];
476  ih2 = j;
477  }
478  TVector3 v1(wire_pitch * (wire - vw[cstat][tpc][plane][ih1]),
479  time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]),
480  0);
481  TVector3 v2(wire_pitch * (wire - vw[cstat][tpc][plane][j]),
482  time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]),
483  0);
484  if (v1.Dot(v2) < 0) {
485  if (dis < mindis2p) {
486  mindis2p = dis;
487  ip2p = vtraj[cstat][tpc][plane][j];
488  ih2p = j;
489  }
490  }
491  }
492  if (ip2p != -1) {
493  ip2 = ip2p;
494  ih2 = ih2p;
495  }
496  if (ip1 == -1 || ip2 == -1) {
497  return;
498  throw cet::exception("CosmicTrackerAlg") << "Cannot find two nearest trajectory points.";
499  }
500  TVector3 v1(wire_pitch * (wire - vw[cstat][tpc][plane][ih1]),
501  time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][ih1]),
502  0);
503  TVector3 v2(wire_pitch * (vw[cstat][tpc][plane][ih2] - vw[cstat][tpc][plane][ih1]),
504  time_pitch * (vt[cstat][tpc][plane][ih2] - vt[cstat][tpc][plane][ih1]),
505  0);
506  if (!v2.Mag()) {
507  throw cet::exception("CosmicTrackerAlg") << "Two identical trajectory points.";
508  }
509  double ratio = v1.Dot(v2) / v2.Mag2();
510  TVector3 pos = trajPos[ip1] + ratio * (trajPos[ip2] - trajPos[ip1]);
511  trkPos.push_back(pos);
512  if (trajDir[ip1].Mag()) { trkDir.push_back(trajDir[ip1]); }
513  else if (trajDir[ip2].Mag()) {
514  trkDir.push_back(trajDir[ip2]);
515  }
516  else { //only got two trajectory points?
517  trkDir.push_back((trajPos[ip2] - trajPos[ip1]).Unit());
518  }
519  } //i
520  }
std::vector< std::vector< std::vector< std::vector< double > > > > vw
std::vector< std::vector< std::vector< std::vector< double > > > > vt
std::vector< TVector3 > trajPos
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
std::vector< TVector3 > trkDir
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
std::vector< std::vector< std::vector< std::vector< unsigned int > > > > vtraj
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
art::ServiceHandle< geo::Geometry const > geom
const detinfo::LArProperties * larprop
std::vector< TVector3 > trkPos
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
Float_t e
Definition: plot.C:35
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
std::vector< TVector3 > trajDir
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::CosmicTrackerAlg::SPTReco ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> &  fHits 
)

Definition at line 49 of file CosmicTrackerAlg.cxx.

Referenced by trkf::CosmicTracker::produce().

52  {
53  trajPos.clear();
54  trajDir.clear();
55  trajHit.clear();
56 
57  trkPos.clear();
58  trkDir.clear();
59 
60  if (fSPTAlg == 0) { TrackTrajectory(detProp, fHits); }
61  else if (fSPTAlg == 1) {
62  Track3D(clockData, detProp, fHits);
63  }
64  else {
65  throw cet::exception("CosmicTrackerAlg")
66  << "Unknown SPTAlg " << fSPTAlg << ", needs to be 0 or 1";
67  }
68 
69  if (!fTrajOnly && trajPos.size())
70  MakeSPT(clockData, detProp, fHits);
71  else {
72  trkPos = trajPos;
73  trkDir = trajDir;
74  }
75  }
std::vector< TVector3 > trajPos
std::vector< TVector3 > trkDir
void TrackTrajectory(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
std::vector< TVector3 > trkPos
void MakeSPT(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
std::vector< TVector3 > trajDir
void Track3D(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &fHits)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit
void trkf::CosmicTrackerAlg::Track3D ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> &  fHits 
)
private

Definition at line 192 of file CosmicTrackerAlg.cxx.

References util::abs(), util::begin(), c1, c2, detinfo::DetectorPropertiesData::ConvertTicksToX(), detinfo::DetectorPropertiesData::DriftVelocity(), e, detinfo::DetectorPropertiesData::Efield(), util::empty(), util::end(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::PlaneID::Plane, detinfo::sampling_rate(), util::size(), detinfo::DetectorPropertiesData::Temperature(), w, geo::WireID::Wire, geo::WireID::WireID(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

195  {
196  larprop = lar::providerFrom<detinfo::LArPropertiesService>();
197 
198  //save time/hit information along track trajectory
199  std::vector<std::map<int, double>> vtimemap(3);
200  std::vector<std::map<int, art::Ptr<recob::Hit>>> vhitmap(3);
201 
202  //find hit on each wire along the fitted line
203  for (size_t ihit = 0; ihit < fHits.size(); ++ihit) { //loop over hits
204  geo::WireID hitWireID = fHits[ihit]->WireID();
205  unsigned int w = hitWireID.Wire;
206  unsigned int pl = hitWireID.Plane;
207  double time = fHits[ihit]->PeakTime();
208  time -= detProp.GetXTicksOffset(
209  fHits[ihit]->WireID().Plane, fHits[ihit]->WireID().TPC, fHits[ihit]->WireID().Cryostat);
210  vtimemap[pl][w] = time;
211  vhitmap[pl][w] = fHits[ihit];
212  }
213 
214  // Find two clusters with the most numbers of hits, and time ranges
215  int iclu1 = -1;
216  int iclu2 = -1;
217  int iclu3 = -1;
218  unsigned maxnumhits0 = 0;
219  unsigned maxnumhits1 = 0;
220 
221  std::vector<double> tmin(vtimemap.size());
222  std::vector<double> tmax(vtimemap.size());
223  for (size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
224  tmin[iclu] = 1e9;
225  tmax[iclu] = -1e9;
226  }
227 
228  for (size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
229  for (auto itime = vtimemap[iclu].begin(); itime != vtimemap[iclu].end(); ++itime) {
230  if (itime->second > tmax[iclu]) { tmax[iclu] = itime->second; }
231  if (itime->second < tmin[iclu]) { tmin[iclu] = itime->second; }
232  }
233  if (vtimemap[iclu].size() > maxnumhits0) {
234  if (iclu1 != -1) {
235  iclu2 = iclu1;
236  maxnumhits1 = maxnumhits0;
237  }
238  iclu1 = iclu;
239  maxnumhits0 = vtimemap[iclu].size();
240  }
241  else if (vtimemap[iclu].size() > maxnumhits1) {
242  iclu2 = iclu;
243  maxnumhits1 = vtimemap[iclu].size();
244  }
245  }
246 
247  std::swap(iclu1, iclu2); //now iclu1 has fewer hits than iclu2
248 
249  //find iclu3
250  for (int iclu = 0; iclu < (int)vtimemap.size(); ++iclu) {
251  if (iclu != iclu1 && iclu != iclu2) iclu3 = iclu;
252  }
253 
254  if (iclu1 != -1 && iclu2 != -1) { //at least two good clusters
255  //select hits in a common time range
256  auto ihit = vhitmap[iclu1].begin();
257  auto itime = vtimemap[iclu1].begin();
258  while (itime != vtimemap[iclu1].end()) {
259  if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) - ftmatch ||
260  itime->second > std::min(tmax[iclu1], tmax[iclu2]) + ftmatch) {
261  vtimemap[iclu1].erase(itime++);
262  vhitmap[iclu1].erase(ihit++);
263  }
264  else {
265  ++itime;
266  ++ihit;
267  }
268  }
269 
270  ihit = vhitmap[iclu2].begin();
271  itime = vtimemap[iclu2].begin();
272  while (itime != vtimemap[iclu2].end()) {
273  if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) - ftmatch ||
274  itime->second > std::min(tmax[iclu1], tmax[iclu2]) + ftmatch) {
275  vtimemap[iclu2].erase(itime++);
276  vhitmap[iclu2].erase(ihit++);
277  }
278  else {
279  ++itime;
280  ++ihit;
281  }
282  }
283 
284  //if one cluster is empty, replace it with iclu3
285  if (!vtimemap[iclu1].size()) {
286  if (iclu3 != -1) { std::swap(iclu3, iclu1); }
287  }
288  if (!vtimemap[iclu2].size()) {
289  if (iclu3 != -1) {
290  std::swap(iclu3, iclu2);
291  std::swap(iclu1, iclu2);
292  }
293  }
294  if ((!vtimemap[iclu1].size()) || (!vtimemap[iclu2].size())) return;
295 
296  bool rev = false;
297  auto times1 = vtimemap[iclu1].begin();
298  auto timee1 = vtimemap[iclu1].end();
299  --timee1;
300  auto times2 = vtimemap[iclu2].begin();
301  auto timee2 = vtimemap[iclu2].end();
302  --timee2;
303 
304  double ts1 = times1->second;
305  double te1 = timee1->second;
306  double ts2 = times2->second;
307  double te2 = timee2->second;
308 
309  //find out if we need to flip ends
310  if (std::abs(ts1 - ts2) + std::abs(te1 - te2) > std::abs(ts1 - te2) + std::abs(te1 - ts2)) {
311  rev = true;
312  }
313 
314  double timetick = sampling_rate(clockData) * 1e-3; // time sample in us
315  double Efield_drift = detProp.Efield(0); // Electric Field in the drift region in kV/cm
316  double Temperature = detProp.Temperature(); // LAr Temperature in K
317 
318  double driftvelocity = detProp.DriftVelocity(
319  Efield_drift, Temperature); //drift velocity in the drift region (cm/us)
320  double timepitch = driftvelocity * timetick;
321 
322  double wire_pitch =
323  geom->WirePitch(vhitmap[0].begin()->second->WireID().asPlaneID()); //wire pitch in cm
324 
325  //find out 2d track length for all clusters associated with track candidate
326  std::vector<double> vtracklength;
327  for (size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
328 
329  double tracklength = 0;
330  if (vtimemap[iclu].size() == 1) { tracklength = wire_pitch; }
331  else if (!vtimemap[iclu].empty()) {
332  std::map<int, double>::const_iterator iw = vtimemap[iclu].cbegin(),
333  wend = vtimemap[iclu].cend();
334  double t0 = iw->first, w0 = iw->second;
335  while (++iw != wend) {
336  tracklength += std::hypot((iw->first - w0) * wire_pitch, (iw->second - t0) * timepitch);
337  w0 = iw->first;
338  t0 = iw->second;
339  } // while
340  }
341  vtracklength.push_back(tracklength);
342  }
343 
344  std::map<int, int> maxhitsMatch;
345 
346  auto ihit1 = vhitmap[iclu1].begin();
347  for (auto itime1 = vtimemap[iclu1].begin(); itime1 != vtimemap[iclu1].end();
348  ++itime1, ++ihit1) { //loop over min-hits
349  std::vector<art::Ptr<recob::Hit>> sp_hits;
350  sp_hits.push_back(ihit1->second);
351  double hitcoord[3] = {0., 0., 0.};
352  double length1 = 0;
353  if (vtimemap[iclu1].size() == 1) { length1 = wire_pitch; }
354  else {
355  for (auto iw1 = vtimemap[iclu1].begin(); iw1 != itime1; ++iw1) {
356  auto iw2 = iw1;
357  ++iw2;
358  length1 += std::hypot((iw1->first - iw2->first) * wire_pitch,
359  (iw1->second - iw2->second) * timepitch);
360  }
361  }
362  double difference = 1e10; //distance between two matched hits
363  auto matchedtime = vtimemap[iclu2].end();
364  auto matchedhit = vhitmap[iclu2].end();
365 
366  auto ihit2 = vhitmap[iclu2].begin();
367  for (auto itime2 = vtimemap[iclu2].begin(); itime2 != vtimemap[iclu2].end();
368  ++itime2, ++ihit2) { //loop over max-hits
369  if (maxhitsMatch[itime2->first]) continue;
370  double length2 = 0;
371  if (vtimemap[iclu2].size() == 1) { length2 = wire_pitch; }
372  else {
373  for (auto iw1 = vtimemap[iclu2].begin(); iw1 != itime2; ++iw1) {
374  auto iw2 = iw1;
375  ++iw2;
376  length2 += std::hypot((iw1->first - iw2->first) * wire_pitch,
377  (iw1->second - iw2->second) * timepitch);
378  }
379  }
380  if (rev) length2 = vtracklength[iclu2] - length2;
381  length2 = vtracklength[iclu1] / vtracklength[iclu2] * length2;
382  bool timematch = std::abs(itime1->second - itime2->second) < ftmatch;
383  if (timematch && std::abs(length2 - length1) < difference) {
384  difference = std::abs(length2 - length1);
385  matchedtime = itime2;
386  matchedhit = ihit2;
387  }
388  } //loop over hits2
389  if (difference < fsmatch) {
390  hitcoord[0] = detProp.ConvertTicksToX(
391  matchedtime->second + detProp.GetXTicksOffset((matchedhit->second)->WireID().Plane,
392  (matchedhit->second)->WireID().TPC,
393  (matchedhit->second)->WireID().Cryostat),
394  (matchedhit->second)->WireID().Plane,
395  (matchedhit->second)->WireID().TPC,
396  (matchedhit->second)->WireID().Cryostat);
397 
398  hitcoord[1] = -1e10;
399  hitcoord[2] = -1e10;
400 
401  //WireID is the exact segment of the wire where the hit is on (1 out of 3 for the 35t)
402 
403  geo::WireID c1 = (ihit1->second)->WireID();
404  geo::WireID c2 = (matchedhit->second)->WireID();
405 
406  geo::WireIDIntersection tmpWIDI;
407  bool sameTpcOrNot = geom->WireIDsIntersect(c1, c2, tmpWIDI);
408 
409  if (sameTpcOrNot) {
410  hitcoord[1] = tmpWIDI.y;
411  hitcoord[2] = tmpWIDI.z;
412  }
413 
414  if (hitcoord[1] > -1e9 && hitcoord[2] > -1e9) {
415  maxhitsMatch[matchedtime->first] = 1;
416  sp_hits.push_back(matchedhit->second);
417  }
418  } //if (difference<fsmatch)
419  if (sp_hits.size() > 1) {
420  trajPos.push_back(TVector3(hitcoord));
421  trajHit.push_back(sp_hits);
422  }
423  } //loop over hits1
424  } //if (iclu1!=-1&&iclu2!=-1){//at least two good clusters
425  }
code to link reconstructed objects back to the MC truth information
double z
z position of intersection
Definition: geo_types.h:792
double ftmatch
tolerance for time matching (in ticks)
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< TVector3 > trajPos
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
intermediate_table::const_iterator const_iterator
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
TCanvas * c1
Definition: plotHisto.C:7
TCanvas * c2
Definition: plot_hist.C:75
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
double fsmatch
tolerance for distance matching (in cm)
art::ServiceHandle< geo::Geometry const > geom
const detinfo::LArProperties * larprop
double y
y position of intersection
Definition: geo_types.h:791
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
Float_t e
Definition: plot.C:35
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Definition: TrackState.h:17
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit
void trkf::CosmicTrackerAlg::TrackTrajectory ( detinfo::DetectorPropertiesData const &  detProp,
std::vector< art::Ptr< recob::Hit >> &  fHits 
)
private

Definition at line 78 of file CosmicTrackerAlg.cxx.

References util::abs(), util::begin(), detinfo::DetectorPropertiesData::ConvertTicksToX(), detinfo::DetectorPropertiesData::ConvertXToTicks(), util::end(), geo::TPCGeo::ID(), geo::TPCGeo::Nplanes(), util::size(), and X.

80  {
81  // Track hit vectors for fitting the trajectory
82  std::array<std::vector<geo::WireID>, 3> trkWID;
83  std::array<std::vector<double>, 3> trkX;
84  std::array<std::vector<double>, 3> trkXErr;
85 
86  std::array<std::vector<art::Ptr<recob::Hit>>, 3> trajHits;
87 
88  for (size_t i = 0; i < fHits.size(); ++i) {
89  trajHits[fHits[i]->WireID().Plane].push_back(fHits[i]);
90  }
91 
92  std::vector<PlnLen> spl;
93  PlnLen plnlen;
94  for (unsigned int ipl = 0; ipl < 3; ++ipl) {
95  plnlen.index = ipl;
96  plnlen.length = trajHits[ipl].size();
97  spl.push_back(plnlen);
98  }
99  std::sort(spl.begin(), spl.end(), greaterThan1);
100  // spl[0] has the most hits and spl.back() has the least
101 
102  for (size_t ipl = 0; ipl < 3; ++ipl) {
103  if (!trajHits[ipl].size()) continue;
104  double xbeg0 = detProp.ConvertTicksToX(trajHits[spl[0].index][0]->PeakTime(),
105  trajHits[spl[0].index][0]->WireID().Plane,
106  trajHits[spl[0].index][0]->WireID().TPC,
107  trajHits[spl[0].index][0]->WireID().Cryostat);
108  double xend0 = detProp.ConvertTicksToX(trajHits[spl[0].index].back()->PeakTime(),
109  trajHits[spl[0].index].back()->WireID().Plane,
110  trajHits[spl[0].index].back()->WireID().TPC,
111  trajHits[spl[0].index].back()->WireID().Cryostat);
112  double xbeg1 = detProp.ConvertTicksToX(trajHits[ipl][0]->PeakTime(),
113  trajHits[ipl][0]->WireID().Plane,
114  trajHits[ipl][0]->WireID().TPC,
115  trajHits[ipl][0]->WireID().Cryostat);
116  double xend1 = detProp.ConvertTicksToX(trajHits[ipl].back()->PeakTime(),
117  trajHits[ipl].back()->WireID().Plane,
118  trajHits[ipl].back()->WireID().TPC,
119  trajHits[ipl].back()->WireID().Cryostat);
120  double dx1 = std::abs(xbeg0 - xbeg1) + std::abs(xend0 - xend1);
121  double dx2 = std::abs(xbeg0 - xend1) + std::abs(xend0 - xbeg1);
122  if (std::abs(xbeg1 - xend1) >
123  0.2 // this is to make sure the track is not completely flat in t,
124  // different by ~2.5 ticks
125  && dx2 < dx1) {
126  std::reverse(trajHits[ipl].begin(), trajHits[ipl].end());
127  }
128  }
129  // make the track trajectory
130  for (size_t ipl = 0; ipl < 3; ++ipl) {
131  // if (ipl == spl.back().index) continue;
132  trkWID[ipl].resize(trajHits[ipl].size());
133  trkX[ipl].resize(trajHits[ipl].size());
134  trkXErr[ipl].resize(trajHits[ipl].size());
135  for (size_t iht = 0; iht < trajHits[ipl].size(); ++iht) {
136  trkWID[ipl][iht] = trajHits[ipl][iht]->WireID();
137  trkX[ipl][iht] = detProp.ConvertTicksToX(trajHits[ipl][iht]->PeakTime(),
138  ipl,
139  trajHits[ipl][iht]->WireID().TPC,
140  trajHits[ipl][iht]->WireID().Cryostat);
141  trkXErr[ipl][iht] = 0.2 * trajHits[ipl][iht]->RMS() * trajHits[ipl][iht]->Multiplicity();
142  } // iht
143  } // ip
144  fTrackTrajectoryAlg.TrackTrajectory(trkWID, trkX, trkXErr, trajPos, trajDir);
145  if (!trajPos.size() || !trajDir.size()) {
146  mf::LogWarning("CosmicTrackerAlg") << "Failed to reconstruct trajectory points.";
147  return;
148  }
149  //remove duplicated points
150  for (auto ipos = trajPos.begin(), idir = trajDir.begin(); ipos != trajPos.end() - 1;) {
151  auto ipos1 = ipos + 1;
152  if (*ipos1 == *ipos) {
153  ipos = trajPos.erase(ipos);
154  idir = trajDir.erase(idir);
155  }
156  else {
157  ++ipos;
158  ++idir;
159  }
160  }
161  vw.clear();
162  vt.clear();
163  vtraj.clear();
164  vw.resize(geom->Ncryostats());
165  vt.resize(geom->Ncryostats());
166  vtraj.resize(geom->Ncryostats());
167  for (auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
168  auto const cstat = cryostat.ID().Cryostat;
169  vw[cstat].resize(cryostat.NTPC());
170  vt[cstat].resize(cryostat.NTPC());
171  vtraj[cstat].resize(cryostat.NTPC());
172  for (size_t tpc = 0; tpc < cryostat.NTPC(); ++tpc) {
173  const geo::TPCGeo& tpcgeom = cryostat.TPC(tpc);
174  vw[cstat][tpc].resize(tpcgeom.Nplanes());
175  vt[cstat][tpc].resize(tpcgeom.Nplanes());
176  vtraj[cstat][tpc].resize(tpcgeom.Nplanes());
177  for (unsigned int plane = 0; plane < tpcgeom.Nplanes(); ++plane) {
178  for (size_t i = 0; i < trajPos.size(); ++i) {
179  double wirecord = geom->WireCoordinate(geo::Point_t{0, trajPos[i].Y(), trajPos[i].Z()},
180  geo::PlaneID{tpcgeom.ID(), plane});
181  double tick = detProp.ConvertXToTicks(trajPos[i].X(), plane, tpc, cstat);
182  vw[cstat][tpc][plane].push_back(wirecord);
183  vt[cstat][tpc][plane].push_back(tick);
184  vtraj[cstat][tpc][plane].push_back(i);
185  } // trajPos.size()
186  } // plane
187  } // tpc
188  } // cstat
189  }
geo::TPCID const & ID() const
Returns the identifier of this TPC.
Definition: TPCGeo.h:270
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
std::vector< std::vector< std::vector< std::vector< double > > > > vw
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::vector< std::vector< std::vector< std::vector< double > > > > vt
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
Geometry information for a single TPC.
Definition: TPCGeo.h:36
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< TVector3 > trajPos
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:430
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::vector< std::vector< std::vector< std::vector< unsigned int > > > > vtraj
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
TrackTrajectoryAlg fTrackTrajectoryAlg
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
void TrackTrajectory(std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
art::ServiceHandle< geo::Geometry const > geom
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
recob::tracking::Plane Plane
Definition: TrackState.h:17
std::vector< TVector3 > trajDir
Float_t X
Definition: plot.C:37

Member Data Documentation

double trkf::CosmicTrackerAlg::fsmatch
private

tolerance for distance matching (in cm)

Definition at line 69 of file CosmicTrackerAlg.h.

int trkf::CosmicTrackerAlg::fSPTAlg
private

Definition at line 54 of file CosmicTrackerAlg.h.

double trkf::CosmicTrackerAlg::ftmatch
private

tolerance for time matching (in ticks)

Definition at line 68 of file CosmicTrackerAlg.h.

TrackTrajectoryAlg trkf::CosmicTrackerAlg::fTrackTrajectoryAlg
private

Definition at line 77 of file CosmicTrackerAlg.h.

bool trkf::CosmicTrackerAlg::fTrajOnly
private

Definition at line 57 of file CosmicTrackerAlg.h.

art::ServiceHandle<geo::Geometry const> trkf::CosmicTrackerAlg::geom
private

Definition at line 84 of file CosmicTrackerAlg.h.

const detinfo::LArProperties* trkf::CosmicTrackerAlg::larprop
private

Definition at line 85 of file CosmicTrackerAlg.h.

std::vector<TVector3> trkf::CosmicTrackerAlg::trajDir

Definition at line 46 of file CosmicTrackerAlg.h.

std::vector<std::vector<art::Ptr<recob::Hit> > > trkf::CosmicTrackerAlg::trajHit

Definition at line 47 of file CosmicTrackerAlg.h.

Referenced by trkf::CosmicTracker::produce().

std::vector<TVector3> trkf::CosmicTrackerAlg::trajPos

Definition at line 45 of file CosmicTrackerAlg.h.

Referenced by trkf::CosmicTracker::produce().

std::vector<TVector3> trkf::CosmicTrackerAlg::trkDir

Definition at line 51 of file CosmicTrackerAlg.h.

Referenced by trkf::CosmicTracker::produce().

std::vector<TVector3> trkf::CosmicTrackerAlg::trkPos

Definition at line 50 of file CosmicTrackerAlg.h.

Referenced by trkf::CosmicTracker::produce().

std::vector<std::vector<std::vector<std::vector<double> > > > trkf::CosmicTrackerAlg::vt
private

Definition at line 81 of file CosmicTrackerAlg.h.

std::vector<std::vector<std::vector<std::vector<unsigned int> > > > trkf::CosmicTrackerAlg::vtraj
private

Definition at line 82 of file CosmicTrackerAlg.h.

std::vector<std::vector<std::vector<std::vector<double> > > > trkf::CosmicTrackerAlg::vw
private

Definition at line 80 of file CosmicTrackerAlg.h.


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