All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 37 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
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:493
art::ServiceHandle< geo::Geometry const > geom
const detinfo::LArProperties * larprop
std::vector< TVector3 > trkPos
do i e
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
auto const detProp
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)
auto const detProp
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  }
double z
z position of intersection
Definition: geo_types.h:805
double ftmatch
tolerance for time matching (in ticks)
std::vector< TVector3 > trajPos
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
BEGIN_PROLOG TPC
T abs(T value)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
double fsmatch
tolerance for distance matching (in cm)
auto begin(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:573
art::ServiceHandle< geo::Geometry const > geom
const detinfo::LArProperties * larprop
double y
y position of intersection
Definition: geo_types.h:804
do i e
then echo ***************************************echo Variable FHICL_FILE_PATH not found echo You porbably haven t set up larsoft echo Try setup uboonecode vXX_XX_XX q e10
Definition: find_fhicl.sh:6
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Definition: TrackState.h:17
bool empty(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:555
auto const detProp
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  }
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
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
BEGIN_PROLOG TPC
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
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
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)
auto begin(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:573
art::ServiceHandle< geo::Geometry const > geom
recob::tracking::Plane Plane
Definition: TrackState.h:17
std::vector< TVector3 > trajDir
auto const detProp

Member Data Documentation

double trkf::CosmicTrackerAlg::fsmatch
private

tolerance for distance matching (in cm)

Definition at line 70 of file CosmicTrackerAlg.h.

int trkf::CosmicTrackerAlg::fSPTAlg
private

Definition at line 55 of file CosmicTrackerAlg.h.

double trkf::CosmicTrackerAlg::ftmatch
private

tolerance for time matching (in ticks)

Definition at line 69 of file CosmicTrackerAlg.h.

TrackTrajectoryAlg trkf::CosmicTrackerAlg::fTrackTrajectoryAlg
private

Definition at line 78 of file CosmicTrackerAlg.h.

bool trkf::CosmicTrackerAlg::fTrajOnly
private

Definition at line 58 of file CosmicTrackerAlg.h.

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

Definition at line 85 of file CosmicTrackerAlg.h.

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

Definition at line 86 of file CosmicTrackerAlg.h.

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

Definition at line 47 of file CosmicTrackerAlg.h.

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

Definition at line 48 of file CosmicTrackerAlg.h.

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

Definition at line 46 of file CosmicTrackerAlg.h.

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

Definition at line 52 of file CosmicTrackerAlg.h.

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

Definition at line 51 of file CosmicTrackerAlg.h.

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

Definition at line 82 of file CosmicTrackerAlg.h.

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

Definition at line 83 of file CosmicTrackerAlg.h.

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

Definition at line 81 of file CosmicTrackerAlg.h.


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