All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | Private Attributes | List of all members
icarus::crt::CRTTrackRecoAlg Class Reference

#include <CRTTrackRecoAlg.h>

Classes

struct  Config
 

Public Member Functions

 CRTTrackRecoAlg (const Config &config)
 
 CRTTrackRecoAlg (const fhicl::ParameterSet &pset)
 
 ~CRTTrackRecoAlg ()
 
void reconfigure (const Config &config)
 
vector< vector< art::Ptr
< sbn::crt::CRTHit > > > 
CreateCRTTzeros (vector< art::Ptr< sbn::crt::CRTHit >>)
 
sbn::crt::CRTTrack FillCrtTrack (sbn::crt::CRTHit hit1, sbn::crt::CRTHit hit2, bool complete)
 
vector< pair< sbn::crt::CRTHit,
vector< int > > > 
AverageHits (vector< art::Ptr< sbn::crt::CRTHit >> hits, map< art::Ptr< sbn::crt::CRTHit >, int > hitIds)
 
vector< sbn::crt::CRTHitAverageHits (vector< art::Ptr< sbn::crt::CRTHit >> hits)
 
sbn::crt::CRTHit DoAverage (vector< art::Ptr< sbn::crt::CRTHit >> hits)
 
vector< pair
< sbn::crt::CRTTrack, vector
< int > > > 
CreateTracks (vector< pair< sbn::crt::CRTHit, vector< int >>> hits)
 
vector< sbn::crt::CRTTrackCreateTracks (vector< sbn::crt::CRTHit > hits)
 
TVector3 CrossPoint (sbn::crt::CRTHit hit, TVector3 start, TVector3 diff)
 

Private Attributes

geo::GeometryCore const * fGeometryService
 
double fTimeLimit
 
double fAverageHitDistance
 
double fDistanceLimit
 
CRTHitRecoAlg hitAlg
 

Detailed Description

Definition at line 62 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.h.

Constructor & Destructor Documentation

CRTTrackRecoAlg::CRTTrackRecoAlg ( const Config config)

Definition at line 6 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

7  : hitAlg() {
8 
9  this->reconfigure(config);
10 
11  fGeometryService = lar::providerFrom<geo::Geometry>();
12 
13 }
icarus::crt::CRTTrackRecoAlg::CRTTrackRecoAlg ( const fhicl::ParameterSet &  pset)
inline

Definition at line 88 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.h.

88  :
89  CRTTrackRecoAlg(fhicl::Table<Config>(pset, {})()) {}
CRTTrackRecoAlg::~CRTTrackRecoAlg ( )

Definition at line 23 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

23 {}

Member Function Documentation

vector< pair< sbn::crt::CRTHit, vector< int > > > CRTTrackRecoAlg::AverageHits ( vector< art::Ptr< sbn::crt::CRTHit >>  hits,
map< art::Ptr< sbn::crt::CRTHit >, int >  hitIds 
)

Definition at line 118 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

119 {
120  vector<pair<sbn::crt::CRTHit, vector<int>>> returnHits;
121  vector<art::Ptr<sbn::crt::CRTHit>> aveHits;
122  vector<art::Ptr<sbn::crt::CRTHit>> spareHits;
123 
124  //if we have CRTHits
125  if (hits.size()>0){
126 
127  bool first = true;
128  TVector3 middle(0., 0., 0.);
129 
130  //loop over CRTHits
131  for (size_t i = 0; i < hits.size(); i++){
132  // Get the position of the hit
133  TVector3 pos(hits[i]->x_pos, hits[i]->y_pos, hits[i]->z_pos);
134  // If first then set average = hit pos
135  if(first){
136  middle = pos;
137  first = false;
138  }
139  // If distance from average < limit then add to average
140  if((pos-middle).Mag() < fAverageHitDistance){
141  aveHits.push_back(hits[i]);
142  }
143  // Else add to another vector
144  else{
145  spareHits.push_back(hits[i]);
146  }
147  }
148 
149  //std::cout << "size of aveHits: " << aveHits.size() << std::endl;
150  // Checking if we have Average CRTHits
151  if (aveHits.size() > 0){
152  //aveHit = DoAverage(aveHits);
153  sbn::crt::CRTHit aveHit = DoAverage(aveHits);
154  vector<int> ids;
155  for(size_t i = 0; i < aveHits.size(); i++){
156  ids.push_back(hitIds[aveHits[i]]);
157  }
158 
159  returnHits.push_back(std::make_pair(aveHit, ids));
160 
161  //Do this recursively
162  vector<pair<sbn::crt::CRTHit, vector<int>>> moreHits = AverageHits(spareHits, hitIds);
163  returnHits.insert(returnHits.end(), moreHits.begin(), moreHits.end());
164  }
165  return returnHits;
166 
167  }//endif hits
168  else { //no hits returns empty vector
169  return returnHits;
170  }
171 
172 } // CRTTrackRecoAlg::AverageHits()
sbn::crt::CRTHit DoAverage(vector< art::Ptr< sbn::crt::CRTHit >> hits)
vector< pair< sbn::crt::CRTHit, vector< int > > > AverageHits(vector< art::Ptr< sbn::crt::CRTHit >> hits, map< art::Ptr< sbn::crt::CRTHit >, int > hitIds)
vector< sbn::crt::CRTHit > CRTTrackRecoAlg::AverageHits ( vector< art::Ptr< sbn::crt::CRTHit >>  hits)

Definition at line 175 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

176 {
177  vector<sbn::crt::CRTHit> returnHits;
178  vector<art::Ptr<sbn::crt::CRTHit>> aveHits;
179  vector<art::Ptr<sbn::crt::CRTHit>> spareHits;
180 
181  if (hits.size()>0){
182  // loop over size of tx
183  bool first = true;
184  TVector3 middle(0., 0., 0.);
185  for (size_t i = 0; i < hits.size(); i++){
186  // Get the position of the hit
187  TVector3 pos(hits[i]->x_pos, hits[i]->y_pos, hits[i]->z_pos);
188  // If first then set average = hit pos
189  if(first){
190  middle = pos;
191  first = false;
192  }
193  // If distance from average < limit then add to average
194  if((pos-middle).Mag() < fAverageHitDistance){
195  aveHits.push_back(hits[i]);
196  }
197  // Else add to another vector
198  else{
199  spareHits.push_back(hits[i]);
200  }
201  }
202  sbn::crt::CRTHit aveHit = DoAverage(aveHits);
203  returnHits.push_back(aveHit);
204 
205  //Do this recursively
206  vector<sbn::crt::CRTHit> moreHits = AverageHits(spareHits);
207  returnHits.insert(returnHits.end(), moreHits.begin(), moreHits.end());
208 
209  return returnHits;
210 
211  }
212  else {
213  return returnHits;
214  }
215 } // CRTTrackRecoAlg::AverageHits()
sbn::crt::CRTHit DoAverage(vector< art::Ptr< sbn::crt::CRTHit >> hits)
vector< pair< sbn::crt::CRTHit, vector< int > > > AverageHits(vector< art::Ptr< sbn::crt::CRTHit >> hits, map< art::Ptr< sbn::crt::CRTHit >, int > hitIds)
std::vector< std::vector< art::Ptr< sbn::crt::CRTHit > > > CRTTrackRecoAlg::CreateCRTTzeros ( vector< art::Ptr< sbn::crt::CRTHit >>  hits)

Definition at line 35 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

36 {
37 
38  std::vector<std::vector<art::Ptr<sbn::crt::CRTHit>>> crtTzeroVect;
39  int iflag[2000] = {};
40 
41  // Sort CRTHits by time
42  std::sort(hits.begin(), hits.end(), [](auto& left, auto& right)->bool{
43  return left->ts0_ns < right->ts0_ns;});
44 
45  // Loop over crt hits
46  for(size_t i = 0; i<hits.size(); i++){
47  //if hit unused
48  if(iflag[i] == 0){
49  vector<art::Ptr<sbn::crt::CRTHit>> crtTzero;
50  double time_ns_A = hits[i]->ts0_ns;
51  iflag[i]=1;
52  crtTzero.push_back(hits[i]);
53 
54  // Sort into a Tzero collection
55  // Loop over all the other CRT hits
56  for(size_t j = i+1; j<hits.size(); j++){
57 
58  //if hit unused
59  if(iflag[j] == 0){
60  // If ts1_ns - ts1_ns < diff then put them in a vector
61  double time_ns_B = hits[j]->ts0_ns;
62  double diff = std::abs(time_ns_B - time_ns_A) * 1e-3; // [us]
63  if(diff < fTimeLimit){
64  iflag[j] = 1; //mark hit used
65  crtTzero.push_back(hits[j]);
66  }
67  }
68  }
69 
70  crtTzeroVect.push_back(crtTzero);
71  }//endif hit unused
72  }
73  return crtTzeroVect;
74 }//CRTTrackRecoAlg::CreateCRTTzeros
walls no right
Definition: selectors.fcl:105
T abs(T value)
walls no left
Definition: selectors.fcl:105
do i e
vector< pair< sbn::crt::CRTTrack, vector< int > > > CRTTrackRecoAlg::CreateTracks ( vector< pair< sbn::crt::CRTHit, vector< int >>>  hits)

Definition at line 262 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

263 {
264  vector<pair<sbn::crt::CRTTrack, vector<int>>> returnTracks;
265 
266  //Store list of hit pairs with distance between them
267  vector<pair<pair<size_t, size_t>, double>> hitPairDist;
268  vector<pair<size_t, size_t>> usedPairs;
269 
270  //Calculate the distance between all hits on different planes
271  for(size_t i = 0; i < hits.size(); i++){
272 
273  sbn::crt::CRTHit hit1 = hits[i].first;
274 
275  for(size_t j = 0; j < hits.size(); j++){
276 
277  sbn::crt::CRTHit hit2 = hits[j].first;
278  pair<size_t, size_t> hitPair = std::make_pair(i, j);
279  pair<size_t, size_t> rhitPair = std::make_pair(j, i);
280 
281  //Only compare hits on different taggers and don't reuse hits
282  if(hit1.tagger!=hit2.tagger && std::find(usedPairs.begin(), usedPairs.end(), rhitPair)==usedPairs.end()){
283  //Calculate the distance between hits and store
284  TVector3 pos1(hit1.x_pos, hit1.y_pos, hit1.z_pos);
285  TVector3 pos2(hit2.x_pos, hit2.y_pos, hit2.z_pos);
286  double dist = (pos1 - pos2).Mag();
287  usedPairs.push_back(hitPair);
288  hitPairDist.push_back(std::make_pair(hitPair, dist));
289  }
290  }
291  }
292 
293  //Sort map by distance
294  std::sort(hitPairDist.begin(), hitPairDist.end(), [](auto& left, auto& right){
295  return left.second > right.second;});
296 
297  //Store potential hit collections + distance along 1D hit
298  vector<pair<vector<size_t>, double>> tracks;
299  for(size_t i = 0; i < hitPairDist.size(); i++){
300 
301  size_t hit_i = hitPairDist[i].first.first;
302  size_t hit_j = hitPairDist[i].first.second;
303 
304  //Make sure bottom plane hit is always hit_i
305  if(hits[hit_j].first.tagger=="volTaggerBot_0") std::swap(hit_i, hit_j);
306  sbn::crt::CRTHit ihit = hits[hit_i].first;
307  sbn::crt::CRTHit jhit = hits[hit_j].first;
308 
309  //If the bottom plane hit is a 1D hit
310  if(ihit.x_err>100. || ihit.z_err>100.){
311 
312  double facMax = 1;
313  vector<size_t> nhitsMax;
314  double minDist = 99999;
315 
316  //Loop over the length of the 1D hit
317  for(int i = 0; i<21; i++){
318 
319  double fac = (i)/10.;
320  vector<size_t> nhits;
321  double totalDist = 0.;
322  TVector3 start(ihit.x_pos-(1.-fac)*ihit.x_err, ihit.y_pos, ihit.z_pos-(1.-fac)*ihit.z_err);
323  TVector3 end(jhit.x_pos, jhit.y_pos, jhit.z_pos);
324  TVector3 diff = start - end;
325 
326  //Loop over the rest of the hits
327  for(size_t k = 0; k < hits.size(); k++){
328 
329  if(k == hit_i || k == hit_j || hits[k].first.tagger == ihit.tagger || hits[k].first.tagger == jhit.tagger)
330  continue;
331 
332  //Calculate the distance between the track crossing point and the true hit
333  sbn::crt::CRTHit khit = hits[k].first;
334  TVector3 mid(khit.x_pos, khit.y_pos, khit.z_pos);
335  TVector3 cross = CrossPoint(khit, start, diff);
336  double dist = (cross-mid).Mag();
337 
338  //If the distance is less than some limit add the hit to the track and record the distance
339  if(dist < fDistanceLimit){
340  nhits.push_back(k);
341  totalDist += dist;
342  }
343  }
344 
345  //If the distance down the 1D hit means more hits are included and they are closer to the track record it
346  if(nhits.size()>=nhitsMax.size() && totalDist/nhits.size() < minDist){
347  nhitsMax = nhits;
348  facMax = fac;
349  minDist = totalDist/nhits.size();
350  }
351  nhits.clear();
352  }
353 
354  //Record the track candidate
355  vector<size_t> trackCand;
356  trackCand.push_back(hit_i);
357  trackCand.push_back(hit_j);
358  trackCand.insert(trackCand.end(), nhitsMax.begin(), nhitsMax.end());
359  tracks.push_back(std::make_pair(trackCand, facMax));
360  }
361 
362  //If there is no 1D hit
363  else{
364  TVector3 start(ihit.x_pos, ihit.y_pos, ihit.z_pos);
365  TVector3 end(jhit.x_pos, jhit.y_pos, jhit.z_pos);
366  TVector3 diff = start - end;
367  vector<size_t> trackCand;
368  trackCand.push_back(hit_i);
369  trackCand.push_back(hit_j);
370 
371  //Loop over all the other hits
372  for(size_t k = 0; k < hits.size(); k++){
373 
374  if(k == hit_i || k == hit_j || hits[k].first.tagger == ihit.tagger || hits[k].first.tagger == jhit.tagger)
375  continue;
376 
377  //Calculate distance to other hits not on the planes of the track hits
378  sbn::crt::CRTHit khit = hits[k].first;
379  TVector3 mid(khit.x_pos, khit.y_pos, khit.z_pos);
380  TVector3 cross = CrossPoint(khit, start, diff);
381  double dist = (cross-mid).Mag();
382 
383  //Record any within a certain distance
384  if(dist < fDistanceLimit){
385  trackCand.push_back(k);
386  }
387  }
388  tracks.push_back(std::make_pair(trackCand, 1));
389  }
390  }
391 
392  //Sort track candidates by number of hits
393  std::sort(tracks.begin(), tracks.end(), [](auto& left, auto& right){
394  return left.first.size() > right.first.size();});
395 
396  //Record used hits
397  vector<size_t> usedHits;
398 
399  //Loop over candidates
400  for(auto& track : tracks){
401 
402  size_t hit_i = track.first[0];
403  size_t hit_j = track.first[1];
404 
405  // Make sure the first hit is the top high tagger if there are only two hits
406  if(hits[hit_j].first.tagger=="volTaggerTopHigh_0")
407  std::swap(hit_i, hit_j);
408 
409  sbn::crt::CRTHit ihit = hits[hit_i].first;
410  sbn::crt::CRTHit jhit = hits[hit_j].first;
411 
412  //Check no hits in track have been used
413  bool used = false;
414 
415  //Loop over hits in track candidate
416  for(size_t i = 0; i < track.first.size(); i++){
417  //Check if any of the hits have been used
418  if(std::find(usedHits.begin(), usedHits.end(), track.first[i]) != usedHits.end())
419  used=true;
420  }
421  //If any of the hits have already been used skip this track
422  if(used)
423  continue;
424 
425  ihit.x_pos -= (1.-track.second)*ihit.x_err;
426  ihit.z_pos -= (1.-track.second)*ihit.z_err;
427 
428  //Create track
429  sbn::crt::CRTTrack crtTrack = FillCrtTrack(ihit, jhit, true);
430 
431  //If only the top two planes are hit create an incomplete/stopping track
432  if(track.first.size()==2 && ihit.tagger == "volTaggerTopHigh_0" && jhit.tagger == "volTaggerTopLow_0"){
433  crtTrack.complete = false;
434  }
435 
436  vector<int> ids;
437  for(size_t i = 0; i < track.first.size(); i++){
438  ids.insert(ids.end(), hits[i].second.begin(), hits[i].second.end());
439  }
440 
441  returnTracks.push_back(std::make_pair(crtTrack, ids));
442 
443  //Record which hits were used only if the track has more than two hits
444  //If there are multiple 2 hit tracks there is no way to distinguish between them
445  //TODO: Add charge matching for ambiguous cases
446  for(size_t i = 0; i < track.first.size(); i++){
447  if(track.first.size()>2) usedHits.push_back(track.first[i]);
448  }
449  }
450  return returnTracks;
451 
452 } // CRTTrackRecoAlg::CreateTracks()
float z_err
position uncertainty in z-direction (cm).
Definition: CRTHit.hh:43
float x_err
position uncertainty in x-direction (cm).
Definition: CRTHit.hh:39
TVector3 CrossPoint(sbn::crt::CRTHit hit, TVector3 start, TVector3 diff)
ClusterModuleLabel join with tracks
walls no right
Definition: selectors.fcl:105
process_name use argoneut_mc_hitfinder track
float z_pos
position in z-direction (cm).
Definition: CRTHit.hh:42
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
sbn::crt::CRTTrack FillCrtTrack(sbn::crt::CRTHit hit1, sbn::crt::CRTHit hit2, bool complete)
walls no left
Definition: selectors.fcl:105
bool complete
Whether or not the track is complete.
Definition: CRTTrack.hh:53
float y_pos
position in y-direction (cm).
Definition: CRTHit.hh:40
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
float x_pos
position in x-direction (cm).
Definition: CRTHit.hh:38
pdgs k
Definition: selectors.fcl:22
Vector cross(Vector const &a, Vector const &b)
Return cross product of two vectors.
std::string tagger
Name of the CRT wall (in the form of strings).
Definition: CRTHit.hh:45
vector< sbn::crt::CRTTrack > CRTTrackRecoAlg::CreateTracks ( vector< sbn::crt::CRTHit hits)

Definition at line 455 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

456 {
457  vector<sbn::crt::CRTTrack> returnTracks;
458  //Store list of hit pairs with distance between them
459  vector<pair<pair<size_t, size_t>, double>> hitPairDist;
460  vector<pair<size_t, size_t>> usedPairs;
461 
462  //Calculate the distance between all hits on different planes
463  for(size_t i = 0; i < hits.size(); i++){
464  sbn::crt::CRTHit hit1 = hits[i];
465  for(size_t j = 0; j < hits.size(); j++){
466 
467  sbn::crt::CRTHit hit2 = hits[j];
468  pair<size_t, size_t> hitPair = std::make_pair(i, j);
469  pair<size_t, size_t> rhitPair = std::make_pair(j, i);
470 
471  //Only compare hits on different taggers and don't reuse hits
472  if(hit1.tagger!=hit2.tagger && std::find(usedPairs.begin(), usedPairs.end(), rhitPair)==usedPairs.end()){
473  //Calculate the distance between hits and store
474  TVector3 pos1(hit1.x_pos, hit1.y_pos, hit1.z_pos);
475  TVector3 pos2(hit2.x_pos, hit2.y_pos, hit2.z_pos);
476  double dist = (pos1 - pos2).Mag();
477  usedPairs.push_back(hitPair);
478  hitPairDist.push_back(std::make_pair(hitPair, dist));
479  }
480  }
481  }
482 
483  //Sort map by distance
484  std::sort(hitPairDist.begin(), hitPairDist.end(), [](auto& left, auto& right){
485  return left.second > right.second;});
486 
487  //Store potential hit collections + distance along 1D hit
488  vector<pair<vector<size_t>, double>> tracks;
489  for(size_t i = 0; i < hitPairDist.size(); i++){
490 
491  size_t hit_i = hitPairDist[i].first.first;
492  size_t hit_j = hitPairDist[i].first.second;
493 
494  //Make sure bottom plane hit is always hit_i
495  if(hits[hit_j].tagger=="volTaggerBot_0")
496  std::swap(hit_i, hit_j);
497 
498  sbn::crt::CRTHit ihit = hits[hit_i];
499  sbn::crt::CRTHit jhit = hits[hit_j];
500 
501  //If the bottom plane hit is a 1D hit
502  if(ihit.x_err>100. || ihit.z_err>100.){
503 
504  double facMax = 1;
505  vector<size_t> nhitsMax;
506  double minDist = 99999;
507 
508  //Loop over the length of the 1D hit
509  for(int i = 0; i<21; i++){
510 
511  double fac = (i)/10.;
512  vector<size_t> nhits;
513  double totalDist = 0.;
514  TVector3 start(ihit.x_pos-(1.-fac)*ihit.x_err, ihit.y_pos, ihit.z_pos-(1.-fac)*ihit.z_err);
515  TVector3 end(jhit.x_pos, jhit.y_pos, jhit.z_pos);
516  TVector3 diff = start - end;
517 
518  //Loop over the rest of the hits
519  for(size_t k = 0; k < hits.size(); k++){
520 
521  if(k == hit_i || k == hit_j || hits[k].tagger == ihit.tagger || hits[k].tagger == jhit.tagger)
522  continue;
523 
524  //Calculate the distance between the track crossing point and the true hit
525  sbn::crt::CRTHit khit = hits[k];
526  TVector3 mid(khit.x_pos, khit.y_pos, khit.z_pos);
527  TVector3 cross = CrossPoint(khit, start, diff);
528  double dist = (cross-mid).Mag();
529 
530  //If the distance is less than some limit add the hit to the track and record the distance
531  if(dist < fDistanceLimit){
532  nhits.push_back(k);
533  totalDist += dist;
534  }
535  }
536  //If the distance down the 1D hit means more hits are included and they are closer to the track record it
537  if(nhits.size()>=nhitsMax.size() && totalDist/nhits.size() < minDist){
538  nhitsMax = nhits;
539  facMax = fac;
540  minDist = totalDist/nhits.size();
541  }
542  nhits.clear();
543  }
544  //Record the track candidate
545  vector<size_t> trackCand;
546  trackCand.push_back(hit_i);
547  trackCand.push_back(hit_j);
548  trackCand.insert(trackCand.end(), nhitsMax.begin(), nhitsMax.end());
549  tracks.push_back(std::make_pair(trackCand, facMax));
550  }
551  //If there is no 1D hit
552  else{
553  TVector3 start(ihit.x_pos, ihit.y_pos, ihit.z_pos);
554  TVector3 end(jhit.x_pos, jhit.y_pos, jhit.z_pos);
555  TVector3 diff = start - end;
556  vector<size_t> trackCand;
557  trackCand.push_back(hit_i);
558  trackCand.push_back(hit_j);
559 
560  //Loop over all the other hits
561  for(size_t k = 0; k < hits.size(); k++){
562 
563  if(k == hit_i || k == hit_j || hits[k].tagger == ihit.tagger || hits[k].tagger == jhit.tagger)
564  continue;
565 
566  //Calculate distance to other hits not on the planes of the track hits
567  sbn::crt::CRTHit khit = hits[k];
568  TVector3 mid(khit.x_pos, khit.y_pos, khit.z_pos);
569  TVector3 cross = CrossPoint(khit, start, diff);
570  double dist = (cross-mid).Mag();
571 
572  //Record any within a certain distance
573  if(dist < fDistanceLimit){
574  trackCand.push_back(k);
575  }
576  }
577  tracks.push_back(std::make_pair(trackCand, 1));
578  }
579  }
580 
581  //Sort track candidates by number of hits
582  std::sort(tracks.begin(), tracks.end(), [](auto& left, auto& right){
583  return left.first.size() > right.first.size();});
584 
585  //Record used hits
586  vector<size_t> usedHits;
587 
588  //Loop over candidates
589  for(auto& track : tracks){
590 
591  size_t hit_i = track.first[0];
592  size_t hit_j = track.first[1];
593 
594  // Make sure the first hit is the top high tagger if there are only two hits
595  if(hits[hit_j].tagger=="volTaggerTopHigh_0")
596  std::swap(hit_i, hit_j);
597 
598  sbn::crt::CRTHit ihit = hits[hit_i];
599  sbn::crt::CRTHit jhit = hits[hit_j];
600 
601  //Check no hits in track have been used
602  bool used = false;
603  //Loop over hits in track candidate
604  for(size_t i = 0; i < track.first.size(); i++){
605  //Check if any of the hits have been used
606  if(std::find(usedHits.begin(), usedHits.end(), track.first[i]) != usedHits.end())
607  used=true;
608  }
609  //If any of the hits have already been used skip this track
610  if(used)
611  continue;
612  ihit.x_pos -= (1.-track.second)*ihit.x_err;
613  ihit.z_pos -= (1.-track.second)*ihit.z_err;
614 
615  //Create track
616  sbn::crt::CRTTrack crtTrack = FillCrtTrack(ihit, jhit, true);
617 
618  //If only the top two planes are hit create an incomplete/stopping track
619  if(track.first.size()==2 && ihit.tagger == "volTaggerTopHigh_0" && jhit.tagger == "volTaggerTopLow_0"){
620  crtTrack.complete = false;
621  }
622 
623  returnTracks.push_back(crtTrack);
624 
625  //Record which hits were used only if the track has more than two hits
626  //If there are multiple 2 hit tracks there is no way to distinguish between them
627  //TODO: Add charge matching for ambiguous cases
628  for(size_t i = 0; i < track.first.size(); i++){
629  if(track.first.size()>2) usedHits.push_back(track.first[i]);
630  }
631  }
632 
633  return returnTracks;
634 
635 } // CRTTrackRecoAlg::CreateTracks()
float z_err
position uncertainty in z-direction (cm).
Definition: CRTHit.hh:43
float x_err
position uncertainty in x-direction (cm).
Definition: CRTHit.hh:39
TVector3 CrossPoint(sbn::crt::CRTHit hit, TVector3 start, TVector3 diff)
ClusterModuleLabel join with tracks
walls no right
Definition: selectors.fcl:105
process_name use argoneut_mc_hitfinder track
float z_pos
position in z-direction (cm).
Definition: CRTHit.hh:42
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
sbn::crt::CRTTrack FillCrtTrack(sbn::crt::CRTHit hit1, sbn::crt::CRTHit hit2, bool complete)
walls no left
Definition: selectors.fcl:105
bool complete
Whether or not the track is complete.
Definition: CRTTrack.hh:53
float y_pos
position in y-direction (cm).
Definition: CRTHit.hh:40
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
float x_pos
position in x-direction (cm).
Definition: CRTHit.hh:38
pdgs k
Definition: selectors.fcl:22
Vector cross(Vector const &a, Vector const &b)
Return cross product of two vectors.
std::string tagger
Name of the CRT wall (in the form of strings).
Definition: CRTHit.hh:45
TVector3 CRTTrackRecoAlg::CrossPoint ( sbn::crt::CRTHit  hit,
TVector3  start,
TVector3  diff 
)

Definition at line 638 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

639 {
640  TVector3 cross;
641  // Use the error to get the fixed coordinate of a tagger
642  // FIXME: can this be done better?
643  if(hit.x_err > 0.39 && hit.x_err < 0.41){
644  double xc = hit.x_pos;
645  TVector3 crossp(xc,
646  ((xc - start.X()) / (diff.X()) * diff.Y()) + start.Y(),
647  ((xc - start.X()) / (diff.X()) * diff.Z()) + start.Z());
648  cross = crossp;
649  }
650  else if(hit.y_err > 0.39 && hit.y_err < 0.41){
651  double yc = hit.y_pos;
652  TVector3 crossp(((yc - start.Y()) / (diff.Y()) * diff.X()) + start.X(),
653  yc,
654  ((yc - start.Y()) / (diff.Y()) * diff.Z()) + start.Z());
655  cross = crossp;
656  }
657  else if(hit.z_err > 0.39 && hit.z_err < 0.41){
658  double zc = hit.z_pos;
659  TVector3 crossp(((zc - start.Z()) / (diff.Z()) * diff.X()) + start.X(),
660  ((zc - start.Z()) / (diff.Z()) * diff.Y()) + start.Y(),
661  zc);
662  cross = crossp;
663  }
664 
665  return cross;
666 
667 } // CRTTrackRecoAlg::CrossPoint()
float z_err
position uncertainty in z-direction (cm).
Definition: CRTHit.hh:43
float x_err
position uncertainty in x-direction (cm).
Definition: CRTHit.hh:39
float y_err
position uncertainty in y-direction (cm).
Definition: CRTHit.hh:41
float z_pos
position in z-direction (cm).
Definition: CRTHit.hh:42
float y_pos
position in y-direction (cm).
Definition: CRTHit.hh:40
float x_pos
position in x-direction (cm).
Definition: CRTHit.hh:38
Vector cross(Vector const &a, Vector const &b)
Return cross product of two vectors.
sbn::crt::CRTHit CRTTrackRecoAlg::DoAverage ( vector< art::Ptr< sbn::crt::CRTHit >>  hits)

Definition at line 218 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

219 {
220  //std::cout << "hits inside CRTTrackRecoAlg::DoAverage:++++++++++++++++ "<< hits[0] << std::endl;
221  // Initialize values
222  std::string tagger = hits[0]->tagger;
223  double xpos = 0.;
224  double ypos = 0.;
225  double zpos = 0.;
226  double xmax = -99999; double xmin = 99999;
227  double ymax = -99999; double ymin = 99999;
228  double zmax = -99999; double zmin = 99999;
229  double ts0_ns = 0., ts1_ns = 0.;
230  int nhits = 0;
231 
232  // Loop over hits
233  for( auto& hit : hits ){
234  // Get the mean x,y,z and times
235  xpos += hit->x_pos;
236  ypos += hit->y_pos;
237  zpos += hit->z_pos;
238  ts0_ns += (double)(int)hit->ts0_ns;
239  ts1_ns += (double)(int)hit->ts1_ns;
240  // For the errors get the maximum limits
241  if(hit->x_pos + hit->x_err > xmax) xmax = hit->x_pos + hit->x_err;
242  if(hit->x_pos - hit->x_err < xmin) xmin = hit->x_pos - hit->x_err;
243  if(hit->y_pos + hit->y_err > ymax) ymax = hit->y_pos + hit->y_err;
244  if(hit->y_pos - hit->y_err < ymin) ymin = hit->y_pos - hit->y_err;
245  if(hit->z_pos + hit->z_err > zmax) zmax = hit->z_pos + hit->z_err;
246  if(hit->z_pos - hit->z_err < zmin) zmin = hit->z_pos - hit->z_err;
247  // Add all the unique IDs in the vector
248  nhits++;
249  }
250 
251  // Create a hit
252  sbn::crt::CRTHit crtHit = hitAlg.FillCRTHit(hits[0]->feb_id, hits[0]->pesmap, hits[0]->peshit,
253  (ts0_ns/nhits)*1e-3, (ts1_ns/nhits)*1e-3, hits[0]->plane, xpos/nhits, (xmax-xmin)/2,
254  ypos/nhits, (ymax-ymin)/2., zpos/nhits, (zmax-zmin)/2., tagger);
255 
256  // std::cout << "hits inside CRTTrackRecoAlg::DoAverage:++++++++++++++++ returning......... line 251" << std::endl;
257  return crtHit;
258 
259 } // CRTTrackRecoAlg::DoAverage()
CRTHit FillCRTHit(vector< uint8_t > tfeb_id, map< uint8_t, vector< pair< int, float >>> tpesmap, float peshit, uint64_t time0, Long64_t time1, int plane, double x, double ex, double y, double ey, double z, double ez, string tagger)
process_name hit
Definition: cheaterreco.fcl:51
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
do i e
sbn::crt::CRTTrack CRTTrackRecoAlg::FillCrtTrack ( sbn::crt::CRTHit  hit1,
sbn::crt::CRTHit  hit2,
bool  complete 
)

Definition at line 77 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

78 {
79  sbn::crt::CRTTrack newtr;
80  newtr.ts0_s = (hit1.ts0_s + hit2.ts0_s)/2.;
81  newtr.ts0_s_err = (uint32_t)((hit1.ts0_s - hit2.ts0_s)/2.);
82  newtr.ts0_ns_h1 = hit1.ts0_ns;
83  newtr.ts0_ns_err_h1 = hit1.ts0_ns_corr;
84  newtr.ts0_ns_h2 = hit2.ts0_ns;
85  newtr.ts0_ns_err_h2 = hit2.ts0_ns_corr;
86  newtr.ts0_ns = (uint32_t)((hit1.ts0_ns + hit2.ts0_ns)/2.);
87  newtr.ts0_ns_err = (uint16_t)(sqrt(hit1.ts0_ns_corr*hit1.ts0_ns_corr + hit2.ts0_ns_corr*hit2.ts0_ns_corr)/2.);
88  newtr.ts1_ns = (int32_t)(((double)(int)hit1.ts1_ns + (double)(int)hit2.ts1_ns)/2.);
89  newtr.ts1_ns_err = (uint16_t)(sqrt(hit1.ts0_ns_corr*hit1.ts0_ns_corr + hit2.ts0_ns_corr*hit2.ts0_ns_corr)/2.);
90  newtr.peshit = hit1.peshit+hit2.peshit;
91  newtr.x1_pos = hit1.x_pos;
92  newtr.x1_err = hit1.x_err;
93  newtr.y1_pos = hit1.y_pos;
94  newtr.y1_err = hit1.y_err;
95  newtr.z1_pos = hit1.z_pos;
96  newtr.z1_err = hit1.z_err;
97  newtr.x2_pos = hit2.x_pos;
98  newtr.x2_err = hit2.x_err;
99  newtr.y2_pos = hit2.y_pos;
100  newtr.y2_err = hit2.y_err;
101  newtr.z2_pos = hit2.z_pos;
102  newtr.z2_err = hit2.z_err;
103  float deltax = hit1.x_pos - hit2.x_pos;
104  float deltay = hit1.y_pos - hit2.y_pos;
105  float deltaz = hit1.z_pos - hit2.z_pos;
106  newtr.length = sqrt(deltax*deltax + deltay*deltay+deltaz*deltaz);
107  newtr.thetaxy = atan2(deltax,deltay);
108  newtr.phizy = atan2(deltaz,deltay);
109  newtr.plane1 = hit1.plane;
110  newtr.plane2 = hit2.plane;
111  newtr.complete = complete;
112 
113  return(newtr);
114 
115 } // CRTTrackRecoAlg::FillCrtTrack()
float z_err
position uncertainty in z-direction (cm).
Definition: CRTHit.hh:43
float x_err
position uncertainty in x-direction (cm).
Definition: CRTHit.hh:39
double ts0_ns_err_h2
T0 time error of second CRTHit.
Definition: CRTTrack.hh:51
double ts1_ns_err
Error on average T1 (nanosecond) of the two hits making the track.
Definition: CRTTrack.hh:29
double ts0_s
Average time (second) of the two hits making the track.
Definition: CRTTrack.hh:24
double ts1_ns
Timestamp T1 ([signal time w.r.t. Trigger time]), in UTC absolute time scale in nanoseconds from the ...
Definition: CRTHit.hh:34
double ts0_ns_err
Error on average T0 (nanosecond) of the two hits making the track.
Definition: CRTTrack.hh:27
int plane
Name of the CRT wall (in the form of numbers).
Definition: CRTHit.hh:36
float peshit
Total photo-electron (PE) in a crt hit.
Definition: CRTHit.hh:27
float y_err
position uncertainty in y-direction (cm).
Definition: CRTHit.hh:41
double ts0_s_err
Average time (second) spread of the two hits making the track.
Definition: CRTTrack.hh:25
float x1_pos
X position of first CRTHit.
Definition: CRTTrack.hh:33
float y1_err
Y position error of first CRTHit.
Definition: CRTTrack.hh:36
double ts0_ns_corr
[Honestly, not sure at this point, it was there since long time (BB)]
Definition: CRTHit.hh:33
double ts0_ns_h1
T0 time of first CRTHit.
Definition: CRTTrack.hh:48
uint64_t ts0_s
Second-only part of timestamp T0.
Definition: CRTHit.hh:29
float z_pos
position in z-direction (cm).
Definition: CRTHit.hh:42
double ts0_ns
Timestamp T0 (from White Rabbit), in UTC absolute time scale in nanoseconds from the Epoch...
Definition: CRTHit.hh:32
float z1_err
Z position error of first CRTHit.
Definition: CRTTrack.hh:38
float y2_err
Y position error of second CRTHit.
Definition: CRTTrack.hh:42
double ts1_ns
Average T1 (nanosecond) of the two hits making the track.
Definition: CRTTrack.hh:28
float length
Track length.
Definition: CRTTrack.hh:45
double ts0_ns
Average T0 (nanosecond) of the two hits making the track.
Definition: CRTTrack.hh:26
float peshit
Total photoelectrons for this track (sum of PEs from the two CRTHits)
Definition: CRTTrack.hh:23
double ts0_ns_err_h1
T0 time error of first CRTHit.
Definition: CRTTrack.hh:49
bool complete
Whether or not the track is complete.
Definition: CRTTrack.hh:53
float z1_pos
Z position of first CRTHit.
Definition: CRTTrack.hh:37
double ts0_ns_h2
T0 time of second CRTHit.
Definition: CRTTrack.hh:50
float y_pos
position in y-direction (cm).
Definition: CRTHit.hh:40
float phizy
Track angle on the Z-Y plane.
Definition: CRTTrack.hh:47
int plane1
Plane ID of first CRTHit.
Definition: CRTTrack.hh:30
float z2_pos
Z position of second CRTHit.
Definition: CRTTrack.hh:43
float x_pos
position in x-direction (cm).
Definition: CRTHit.hh:38
float y2_pos
Y position of second CRTHit.
Definition: CRTTrack.hh:41
float y1_pos
Y position of first CRTHit.
Definition: CRTTrack.hh:35
float x2_pos
X position of second CRTHit.
Definition: CRTTrack.hh:39
float x1_err
X position error of first CRTHit.
Definition: CRTTrack.hh:34
int plane2
Plane ID of second CRTHit.
Definition: CRTTrack.hh:31
float z2_err
Z position error of second CRTHit.
Definition: CRTTrack.hh:44
float x2_err
X position error of second CRTHit.
Definition: CRTTrack.hh:40
float thetaxy
Track angle on the X-Y plane.
Definition: CRTTrack.hh:46
void CRTTrackRecoAlg::reconfigure ( const Config config)

Definition at line 25 of file icaruscode/icaruscode/CRT/CRTUtils/CRTTrackRecoAlg.cc.

25  {
26 
27  fTimeLimit = config.TimeLimit();
28  fAverageHitDistance = config.AverageHitDistance();
29  fDistanceLimit = config.DistanceLimit();
30 
31  return;
32 
33 }

Member Data Documentation

double icarus::crt::CRTTrackRecoAlg::fAverageHitDistance
private
double icarus::crt::CRTTrackRecoAlg::fDistanceLimit
private
geo::GeometryCore const* icarus::crt::CRTTrackRecoAlg::fGeometryService
private
double icarus::crt::CRTTrackRecoAlg::fTimeLimit
private
CRTHitRecoAlg icarus::crt::CRTTrackRecoAlg::hitAlg
private

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