LArSoft  v10_04_05
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
 
geo::WireReadoutGeom const * wireReadoutGeom = &art::ServiceHandle<geo::WireReadout>()->Get()
 
const detinfo::LArPropertieslarprop
 

Detailed Description

Definition at line 34 of file CosmicTrackerAlg.h.

Constructor & Destructor Documentation

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

Definition at line 41 of file CosmicTrackerAlg.cxx.

References fhicl::ParameterSet::get().

42  {
43  fSPTAlg = pset.get<int>("SPTAlg", 0);
44  fTrajOnly = pset.get<bool>("TrajOnly", false);
45  ftmatch = pset.get<double>("TMatch");
46  fsmatch = pset.get<double>("SMatch");
47  }
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 427 of file CosmicTrackerAlg.cxx.

References 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().

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

Definition at line 50 of file CosmicTrackerAlg.cxx.

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

53  {
54  trajPos.clear();
55  trajDir.clear();
56  trajHit.clear();
57 
58  trkPos.clear();
59  trkDir.clear();
60 
61  if (fSPTAlg == 0) { TrackTrajectory(detProp, fHits); }
62  else if (fSPTAlg == 1) {
63  Track3D(clockData, detProp, fHits);
64  }
65  else {
66  throw cet::exception("CosmicTrackerAlg")
67  << "Unknown SPTAlg " << fSPTAlg << ", needs to be 0 or 1";
68  }
69 
70  if (!fTrajOnly && trajPos.size())
71  MakeSPT(clockData, detProp, fHits);
72  else {
73  trkPos = trajPos;
74  trkDir = trajDir;
75  }
76  }
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 194 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, and geo::WireID::WireID().

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

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

81  {
82  // Track hit vectors for fitting the trajectory
83  std::array<std::vector<geo::WireID>, 3> trkWID;
84  std::array<std::vector<double>, 3> trkX;
85  std::array<std::vector<double>, 3> trkXErr;
86 
87  std::array<std::vector<art::Ptr<recob::Hit>>, 3> trajHits;
88 
89  for (size_t i = 0; i < fHits.size(); ++i) {
90  trajHits[fHits[i]->WireID().Plane].push_back(fHits[i]);
91  }
92 
93  std::vector<PlnLen> spl;
94  PlnLen plnlen;
95  for (unsigned int ipl = 0; ipl < 3; ++ipl) {
96  plnlen.index = ipl;
97  plnlen.length = trajHits[ipl].size();
98  spl.push_back(plnlen);
99  }
100  std::sort(spl.begin(), spl.end(), greaterThan1);
101  // spl[0] has the most hits and spl.back() has the least
102 
103  for (size_t ipl = 0; ipl < 3; ++ipl) {
104  if (!trajHits[ipl].size()) continue;
105  double xbeg0 = detProp.ConvertTicksToX(trajHits[spl[0].index][0]->PeakTime(),
106  trajHits[spl[0].index][0]->WireID().Plane,
107  trajHits[spl[0].index][0]->WireID().TPC,
108  trajHits[spl[0].index][0]->WireID().Cryostat);
109  double xend0 = detProp.ConvertTicksToX(trajHits[spl[0].index].back()->PeakTime(),
110  trajHits[spl[0].index].back()->WireID().Plane,
111  trajHits[spl[0].index].back()->WireID().TPC,
112  trajHits[spl[0].index].back()->WireID().Cryostat);
113  double xbeg1 = detProp.ConvertTicksToX(trajHits[ipl][0]->PeakTime(),
114  trajHits[ipl][0]->WireID().Plane,
115  trajHits[ipl][0]->WireID().TPC,
116  trajHits[ipl][0]->WireID().Cryostat);
117  double xend1 = detProp.ConvertTicksToX(trajHits[ipl].back()->PeakTime(),
118  trajHits[ipl].back()->WireID().Plane,
119  trajHits[ipl].back()->WireID().TPC,
120  trajHits[ipl].back()->WireID().Cryostat);
121  double dx1 = std::abs(xbeg0 - xbeg1) + std::abs(xend0 - xend1);
122  double dx2 = std::abs(xbeg0 - xend1) + std::abs(xend0 - xbeg1);
123  if (std::abs(xbeg1 - xend1) >
124  0.2 // this is to make sure the track is not completely flat in t,
125  // different by ~2.5 ticks
126  && dx2 < dx1) {
127  std::reverse(trajHits[ipl].begin(), trajHits[ipl].end());
128  }
129  }
130  // make the track trajectory
131  for (size_t ipl = 0; ipl < 3; ++ipl) {
132  // if (ipl == spl.back().index) continue;
133  trkWID[ipl].resize(trajHits[ipl].size());
134  trkX[ipl].resize(trajHits[ipl].size());
135  trkXErr[ipl].resize(trajHits[ipl].size());
136  for (size_t iht = 0; iht < trajHits[ipl].size(); ++iht) {
137  trkWID[ipl][iht] = trajHits[ipl][iht]->WireID();
138  trkX[ipl][iht] = detProp.ConvertTicksToX(trajHits[ipl][iht]->PeakTime(),
139  ipl,
140  trajHits[ipl][iht]->WireID().TPC,
141  trajHits[ipl][iht]->WireID().Cryostat);
142  trkXErr[ipl][iht] = 0.2 * trajHits[ipl][iht]->RMS() * trajHits[ipl][iht]->Multiplicity();
143  } // iht
144  } // ip
145  fTrackTrajectoryAlg.TrackTrajectory(trkWID, trkX, trkXErr, trajPos, trajDir);
146  if (!trajPos.size() || !trajDir.size()) {
147  mf::LogWarning("CosmicTrackerAlg") << "Failed to reconstruct trajectory points.";
148  return;
149  }
150  //remove duplicated points
151  for (auto ipos = trajPos.begin(), idir = trajDir.begin(); ipos != trajPos.end() - 1;) {
152  auto ipos1 = ipos + 1;
153  if (*ipos1 == *ipos) {
154  ipos = trajPos.erase(ipos);
155  idir = trajDir.erase(idir);
156  }
157  else {
158  ++ipos;
159  ++idir;
160  }
161  }
162  vw.clear();
163  vt.clear();
164  vtraj.clear();
165  vw.resize(geom->Ncryostats());
166  vt.resize(geom->Ncryostats());
167  vtraj.resize(geom->Ncryostats());
168  for (auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
169  auto const cstat = cryostat.ID().Cryostat;
170  vw[cstat].resize(cryostat.NTPC());
171  vt[cstat].resize(cryostat.NTPC());
172  vtraj[cstat].resize(cryostat.NTPC());
173  for (size_t tpc = 0; tpc < cryostat.NTPC(); ++tpc) {
174  const geo::TPCGeo& tpcgeom = cryostat.TPC(tpc);
175  unsigned int nplanes = wireReadoutGeom->Nplanes(tpcgeom.ID());
176  vw[cstat][tpc].resize(nplanes);
177  vt[cstat][tpc].resize(nplanes);
178  vtraj[cstat][tpc].resize(nplanes);
179  for (unsigned int plane = 0; plane < nplanes; ++plane) {
180  for (size_t i = 0; i < trajPos.size(); ++i) {
181  double wirecord = wireReadoutGeom->Plane({tpcgeom.ID(), plane})
182  .WireCoordinate(geo::Point_t{0, trajPos[i].Y(), trajPos[i].Z()});
183  double tick = detProp.ConvertXToTicks(trajPos[i].X(), plane, tpc, cstat);
184  vw[cstat][tpc][plane].push_back(wirecord);
185  vt[cstat][tpc][plane].push_back(tick);
186  vtraj[cstat][tpc][plane].push_back(i);
187  } // trajPos.size()
188  } // plane
189  } // tpc
190  } // cstat
191  }
std::vector< std::vector< std::vector< std::vector< double > > > > vw
std::vector< std::vector< std::vector< std::vector< double > > > > vt
Geometry information for a single TPC.
Definition: TPCGeo.h:33
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:42
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:303
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)
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
art::ServiceHandle< geo::Geometry const > geom
geo::WireReadoutGeom const * wireReadoutGeom
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
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
recob::tracking::Plane Plane
Definition: TrackState.h:17
std::vector< TVector3 > trajDir
Float_t X
Definition: plot.C:37
TPCID const & ID() const
Returns the identifier of this TPC.
Definition: TPCGeo.h:147

Member Data Documentation

double trkf::CosmicTrackerAlg::fsmatch
private

tolerance for distance matching (in cm)

Definition at line 67 of file CosmicTrackerAlg.h.

int trkf::CosmicTrackerAlg::fSPTAlg
private

Definition at line 52 of file CosmicTrackerAlg.h.

double trkf::CosmicTrackerAlg::ftmatch
private

tolerance for time matching (in ticks)

Definition at line 66 of file CosmicTrackerAlg.h.

TrackTrajectoryAlg trkf::CosmicTrackerAlg::fTrackTrajectoryAlg
private

Definition at line 75 of file CosmicTrackerAlg.h.

bool trkf::CosmicTrackerAlg::fTrajOnly
private

Definition at line 55 of file CosmicTrackerAlg.h.

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

Definition at line 82 of file CosmicTrackerAlg.h.

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

Definition at line 84 of file CosmicTrackerAlg.h.

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

Definition at line 44 of file CosmicTrackerAlg.h.

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

Definition at line 45 of file CosmicTrackerAlg.h.

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

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

Definition at line 43 of file CosmicTrackerAlg.h.

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

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

Definition at line 49 of file CosmicTrackerAlg.h.

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

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

Definition at line 48 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 79 of file CosmicTrackerAlg.h.

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

Definition at line 80 of file CosmicTrackerAlg.h.

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

Definition at line 78 of file CosmicTrackerAlg.h.

geo::WireReadoutGeom const* trkf::CosmicTrackerAlg::wireReadoutGeom = &art::ServiceHandle<geo::WireReadout>()->Get()
private

Definition at line 83 of file CosmicTrackerAlg.h.


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