Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
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 33 of file CosmicTrackerAlg.h.

Constructor & Destructor Documentation

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

Definition at line 41 of file CosmicTrackerAlg.cxx.

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 434 of file CosmicTrackerAlg.cxx.

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

Definition at line 51 of file CosmicTrackerAlg.cxx.

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

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

83  {
84  // Track hit vectors for fitting the trajectory
85  std::array<std::vector<geo::WireID>, 3> trkWID;
86  std::array<std::vector<double>, 3> trkX;
87  std::array<std::vector<double>, 3> trkXErr;
88 
89  std::array<std::vector<art::Ptr<recob::Hit>>, 3> trajHits;
90 
91  for (size_t i = 0; i < fHits.size(); ++i) {
92  trajHits[fHits[i]->WireID().Plane].push_back(fHits[i]);
93  }
94 
95  std::vector<PlnLen> spl;
96  PlnLen plnlen;
97  for (unsigned int ipl = 0; ipl < 3; ++ipl) {
98  plnlen.index = ipl;
99  plnlen.length = trajHits[ipl].size();
100  spl.push_back(plnlen);
101  }
102  std::sort(spl.begin(), spl.end(), greaterThan1);
103  // spl[0] has the most hits and spl.back() has the least
104 
105  for (size_t ipl = 0; ipl < 3; ++ipl) {
106  if (!trajHits[ipl].size()) continue;
107  double xbeg0 = detProp.ConvertTicksToX(trajHits[spl[0].index][0]->PeakTime(),
108  trajHits[spl[0].index][0]->WireID().Plane,
109  trajHits[spl[0].index][0]->WireID().TPC,
110  trajHits[spl[0].index][0]->WireID().Cryostat);
111  double xend0 = detProp.ConvertTicksToX(trajHits[spl[0].index].back()->PeakTime(),
112  trajHits[spl[0].index].back()->WireID().Plane,
113  trajHits[spl[0].index].back()->WireID().TPC,
114  trajHits[spl[0].index].back()->WireID().Cryostat);
115  double xbeg1 = detProp.ConvertTicksToX(trajHits[ipl][0]->PeakTime(),
116  trajHits[ipl][0]->WireID().Plane,
117  trajHits[ipl][0]->WireID().TPC,
118  trajHits[ipl][0]->WireID().Cryostat);
119  double xend1 = detProp.ConvertTicksToX(trajHits[ipl].back()->PeakTime(),
120  trajHits[ipl].back()->WireID().Plane,
121  trajHits[ipl].back()->WireID().TPC,
122  trajHits[ipl].back()->WireID().Cryostat);
123  double dx1 = std::abs(xbeg0 - xbeg1) + std::abs(xend0 - xend1);
124  double dx2 = std::abs(xbeg0 - xend1) + std::abs(xend0 - xbeg1);
125  if (std::abs(xbeg1 - xend1) >
126  0.2 // this is to make sure the track is not completely flat in t,
127  // different by ~2.5 ticks
128  && dx2 < dx1) {
129  std::reverse(trajHits[ipl].begin(), trajHits[ipl].end());
130  }
131  }
132  // make the track trajectory
133  for (size_t ipl = 0; ipl < 3; ++ipl) {
134  // if (ipl == spl.back().index) continue;
135  trkWID[ipl].resize(trajHits[ipl].size());
136  trkX[ipl].resize(trajHits[ipl].size());
137  trkXErr[ipl].resize(trajHits[ipl].size());
138  for (size_t iht = 0; iht < trajHits[ipl].size(); ++iht) {
139  trkWID[ipl][iht] = trajHits[ipl][iht]->WireID();
140  trkX[ipl][iht] = detProp.ConvertTicksToX(trajHits[ipl][iht]->PeakTime(),
141  ipl,
142  trajHits[ipl][iht]->WireID().TPC,
143  trajHits[ipl][iht]->WireID().Cryostat);
144  trkXErr[ipl][iht] = 0.2 * trajHits[ipl][iht]->RMS() * trajHits[ipl][iht]->Multiplicity();
145  } // iht
146  } // ip
147  fTrackTrajectoryAlg.TrackTrajectory(trkWID, trkX, trkXErr, trajPos, trajDir);
148  if (!trajPos.size() || !trajDir.size()) {
149  mf::LogWarning("CosmicTrackerAlg") << "Failed to reconstruct trajectory points.";
150  return;
151  }
152  //remove duplicated points
153  for (auto ipos = trajPos.begin(), idir = trajDir.begin(); ipos != trajPos.end() - 1;) {
154  auto ipos1 = ipos + 1;
155  if (*ipos1 == *ipos) {
156  ipos = trajPos.erase(ipos);
157  idir = trajDir.erase(idir);
158  }
159  else {
160  ++ipos;
161  ++idir;
162  }
163  }
164  vw.clear();
165  vt.clear();
166  vtraj.clear();
167  vw.resize(geom->Ncryostats());
168  vt.resize(geom->Ncryostats());
169  vtraj.resize(geom->Ncryostats());
170  for (size_t cstat = 0; cstat < geom->Ncryostats(); ++cstat) {
171  vw[cstat].resize(geom->Cryostat(cstat).NTPC());
172  vt[cstat].resize(geom->Cryostat(cstat).NTPC());
173  vtraj[cstat].resize(geom->Cryostat(cstat).NTPC());
174  for (size_t tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
175  const geo::TPCGeo& tpcgeom = geom->Cryostat(cstat).TPC(tpc);
176  vw[cstat][tpc].resize(tpcgeom.Nplanes());
177  vt[cstat][tpc].resize(tpcgeom.Nplanes());
178  vtraj[cstat][tpc].resize(tpcgeom.Nplanes());
179  for (size_t plane = 0; plane < tpcgeom.Nplanes(); ++plane) {
180  for (size_t i = 0; i < trajPos.size(); ++i) {
181  double wirecord =
182  geom->WireCoordinate(trajPos[i].Y(), trajPos[i].Z(), plane, tpc, cstat);
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  }
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
std::vector< std::vector< std::vector< std::vector< double > > > > vw
std::vector< std::vector< std::vector< std::vector< double > > > > vt
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
Geometry information for a single TPC.
Definition: TPCGeo.h:38
std::vector< TVector3 > trajPos
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
T abs(T value)
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:75
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
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
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
recob::tracking::Plane Plane
Definition: TrackState.h:17
std::vector< TVector3 > trajDir

Member Data Documentation

double trkf::CosmicTrackerAlg::fsmatch
private

tolerance for distance matching (in cm)

Definition at line 66 of file CosmicTrackerAlg.h.

int trkf::CosmicTrackerAlg::fSPTAlg
private

Definition at line 51 of file CosmicTrackerAlg.h.

double trkf::CosmicTrackerAlg::ftmatch
private

tolerance for time matching (in ticks)

Definition at line 65 of file CosmicTrackerAlg.h.

TrackTrajectoryAlg trkf::CosmicTrackerAlg::fTrackTrajectoryAlg
private

Definition at line 74 of file CosmicTrackerAlg.h.

bool trkf::CosmicTrackerAlg::fTrajOnly
private

Definition at line 54 of file CosmicTrackerAlg.h.

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

Definition at line 81 of file CosmicTrackerAlg.h.

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

Definition at line 82 of file CosmicTrackerAlg.h.

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

Definition at line 43 of file CosmicTrackerAlg.h.

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

Definition at line 44 of file CosmicTrackerAlg.h.

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

Definition at line 42 of file CosmicTrackerAlg.h.

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

Definition at line 48 of file CosmicTrackerAlg.h.

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

Definition at line 47 of file CosmicTrackerAlg.h.

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

Definition at line 78 of file CosmicTrackerAlg.h.

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

Definition at line 79 of file CosmicTrackerAlg.h.

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

Definition at line 77 of file CosmicTrackerAlg.h.


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