LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
pma::PMAlgStitching Class Reference

#include "PMAlgStitching.h"

Classes

struct  Config
 

Public Member Functions

 PMAlgStitching (const pma::PMAlgStitching::Config &config)
 
void StitchTracksCPA (const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &tracks)
 
void StitchTracksAPA (const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &tracks)
 

Private Member Functions

void StitchTracks (const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &tracks, bool isCPA)
 
double GetOptimalStitchShift (TVector3 &pos1, TVector3 &pos2, TVector3 &dir1, TVector3 &dir2, double &shift) const
 
double GetTrackPairDelta (TVector3 &pos1, TVector3 &pos2, TVector3 &dir1, TVector3 &dir2) const
 
void GetTPCXOffsets ()
 
double GetTPCOffset (unsigned int tpc, unsigned int cryo, bool isCPA)
 

Private Attributes

std::map< geo::TPCID, double > fTPCXOffsetsAPA
 
std::map< geo::TPCID, double > fTPCXOffsetsCPA
 
double fStitchingThreshold
 
unsigned int fNodesFromEnd
 

Detailed Description

Definition at line 29 of file PMAlgStitching.h.

Constructor & Destructor Documentation

pma::PMAlgStitching::PMAlgStitching ( const pma::PMAlgStitching::Config config)

Definition at line 21 of file PMAlgStitching.cxx.

References fNodesFromEnd, fStitchingThreshold, GetTPCXOffsets(), pma::PMAlgStitching::Config::NodesFromEnd, and pma::PMAlgStitching::Config::StitchingThreshold.

22 {
23 
24  // Set parameters from the config.
26  fNodesFromEnd = config.NodesFromEnd();
27 
28  // Get CPA and APA positions.
30 }
unsigned int fNodesFromEnd
fhicl::Atom< int > StitchingThreshold
fhicl::Atom< unsigned int > NodesFromEnd

Member Function Documentation

double pma::PMAlgStitching::GetOptimalStitchShift ( TVector3 &  pos1,
TVector3 &  pos2,
TVector3 &  dir1,
TVector3 &  dir2,
double &  shift 
) const
private

Definition at line 304 of file PMAlgStitching.cxx.

References GetTrackPairDelta().

Referenced by StitchTracks().

309 {
310 
311  double stepSize = 0.1;
312  double minShift = shift - (50. * stepSize);
313  double maxShift = shift + (50. * stepSize);
314  double bestShift = 99999;
315  double bestScore = 99999;
316 
317  for (shift = minShift; shift <= maxShift; shift += stepSize) {
318  TVector3 newPos1 = pos1;
319  TVector3 newPos2 = pos2;
320  newPos1.SetX(pos1.X() - shift);
321  newPos2.SetX(pos2.X() + shift);
322  double thisScore = GetTrackPairDelta(newPos1, newPos2, dir1, dir2);
323  if (thisScore < bestScore) {
324  bestShift = shift;
325  bestScore = thisScore;
326  }
327  }
328 
329  shift = bestShift;
330  return bestScore;
331 }
double GetTrackPairDelta(TVector3 &pos1, TVector3 &pos2, TVector3 &dir1, TVector3 &dir2) const
double pma::PMAlgStitching::GetTPCOffset ( unsigned int  tpc,
unsigned int  cryo,
bool  isCPA 
)
private

Definition at line 398 of file PMAlgStitching.cxx.

References fTPCXOffsetsAPA, and fTPCXOffsetsCPA.

Referenced by StitchTracks().

399 {
400 
401  geo::TPCID thisTPCID(cryo, tpc);
402  double offset = 0.0;
403  if (isCPA) { offset = fTPCXOffsetsCPA[thisTPCID]; }
404  else {
405  offset = fTPCXOffsetsAPA[thisTPCID];
406  }
407  return offset;
408 }
std::map< geo::TPCID, double > fTPCXOffsetsAPA
std::map< geo::TPCID, double > fTPCXOffsetsCPA
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
void pma::PMAlgStitching::GetTPCXOffsets ( )
private

Definition at line 357 of file PMAlgStitching.cxx.

References fTPCXOffsetsAPA, fTPCXOffsetsCPA, and Get.

Referenced by PMAlgStitching().

358 {
359  auto const* geom = lar::providerFrom<geo::Geometry>();
360  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
361 
362  // Loop over each TPC and store the require information
363  for (geo::TPCGeo const& aTPC : geom->Iterate<geo::TPCGeo>()) {
364  auto const& tID = aTPC.ID();
365 
366  // Loop over the 3 possible readout planes to find the x position
367  // KJK: Then why does the loop go from 0 through 3 (i.e. 4 planes)?
368  unsigned int plane = 0;
369  bool hasPlane = false;
370  for (; plane < 4; ++plane) {
371  hasPlane = wireReadoutGeom.HasPlane({tID, plane});
372  if (hasPlane) { break; }
373  }
374 
375  if (!hasPlane) { continue; }
376 
377  // Get the x position of the readout plane
378  double xAnode = wireReadoutGeom.FirstPlane(tID).GetCenter().X();
379  fTPCXOffsetsAPA.insert(std::make_pair(tID, xAnode));
380 
381  // For the cathode, we have to try a little harder. Firstly, find the min and max x
382  // values for the TPC.
383  auto const center = aTPC.GetCenter();
384  double tpcDim[3] = {aTPC.HalfWidth(), aTPC.HalfHeight(), 0.5 * aTPC.Length()};
385  double xmin = center.X() - tpcDim[0];
386  double xmax = center.X() + tpcDim[0];
387  double xCathode = 0.;
388  // Now check which is further from the APA and use that as the cathode value.
389  if (fabs(xmin - xAnode) > fabs(xmax - xAnode)) { xCathode = xmin; }
390  else {
391  xCathode = xmax;
392  }
393  fTPCXOffsetsCPA.insert(std::make_pair(tID, xCathode));
394  }
395 }
Geometry information for a single TPC.
Definition: TPCGeo.h:33
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::map< geo::TPCID, double > fTPCXOffsetsAPA
std::map< geo::TPCID, double > fTPCXOffsetsCPA
double pma::PMAlgStitching::GetTrackPairDelta ( TVector3 &  pos1,
TVector3 &  pos2,
TVector3 &  dir1,
TVector3 &  dir2 
) const
private

Definition at line 334 of file PMAlgStitching.cxx.

Referenced by GetOptimalStitchShift().

338 {
339 
340  double delta = -999.;
341 
342  // Calculate number of steps to reach the merge point in x.
343  double steps1 = (pos2.X() - pos1.X()) / dir1.X();
344  double steps2 = (pos1.X() - pos2.X()) / dir2.X();
345 
346  // Extrapolate each vector to the other's extrapolation point
347  TVector3 trk1Merge = pos1 + steps1 * dir1;
348  TVector3 trk2Merge = pos2 + steps2 * dir2;
349 
350  // Find the difference between each vector and the extrapolation
351  delta = (trk1Merge - pos2).Mag() + (trk2Merge - pos1).Mag();
352 
353  return delta;
354 }
void pma::PMAlgStitching::StitchTracks ( const detinfo::DetectorClocksData clockData,
const detinfo::DetectorPropertiesData detProp,
pma::TrkCandidateColl tracks,
bool  isCPA 
)
private

Definition at line 53 of file PMAlgStitching.cxx.

References pma::Track3D::ApplyDriftShiftInTree(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), pma::Track3D::Flip(), fNodesFromEnd, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), fStitchingThreshold, pma::TrkCandidateColl::getCandidateIndex(), pma::TrkCandidateColl::getCandidateTreeId(), GetOptimalStitchShift(), pma::Track3D::GetRoot(), GetTPCOffset(), pma::TrkCandidateColl::merge(), pma::Track3D::Nodes(), pma::TrkCandidateColl::setTreeOriginAtFront(), pma::TrkCandidateColl::size(), t1, t2, and pma::TrkCandidateColl::tracks().

Referenced by StitchTracksAPA(), and StitchTracksCPA().

57 {
58 
59  unsigned int minTrkLength = 2 * fNodesFromEnd + 3;
60  // Special case for fNodesFromEnd = 0
61  if (minTrkLength < 6) minTrkLength = 6;
62 
63  // Loop over the track collection
64  unsigned int t = 0;
65  while (t < tracks.size()) {
66 
67  pma::Track3D* t1 = tracks[t].Track();
68  if (t1->Nodes().size() < minTrkLength) {
69  ++t;
70  continue;
71  }
72 
73  // Look through the following tracks for one to stitch
74  pma::Track3D* bestTrkMatch = 0x0;
75 
76  // Don't use the very end points of the tracks in case of scatter or distortion.
77  TVector3 trk1Front = t1->Nodes()[(fNodesFromEnd)]->Point3D();
78  TVector3 trk1Back = t1->Nodes()[t1->Nodes().size() - 1 - fNodesFromEnd]->Point3D();
79  TVector3 trk1FrontDir = (trk1Front - t1->Nodes()[(fNodesFromEnd + 1)]->Point3D()).Unit();
80  TVector3 trk1BackDir =
81  (trk1Back - t1->Nodes()[t1->Nodes().size() - 1 - (fNodesFromEnd + 1)]->Point3D()).Unit();
82 
83  // For stitching, we need to consider both ends of the track.
84  double offsetFront1 = GetTPCOffset(t1->FrontTPC(), t1->FrontCryo(), isCPA);
85  double offsetBack1 = GetTPCOffset(t1->BackTPC(), t1->BackCryo(), isCPA);
86 
87  bool isBestFront1 = false;
88  bool isBestFront2 = false;
89  double xBestShift = 0;
90  double frontShift1 = t1->Nodes()[0]->Point3D().X() - offsetFront1;
91  double backShift1 = t1->Nodes()[t1->Nodes().size() - 1]->Point3D().X() - offsetBack1;
92 
93  double bestMatchScore = 99999;
94 
95  for (unsigned int u = t + 1; u < tracks.size(); ++u) {
96 
97  pma::Track3D* t2 = tracks[u].Track();
98  if (t2->Nodes().size() < minTrkLength) continue;
99 
100  // Don't use the very end points of the tracks in case of scatter or distortion.
101  TVector3 trk2Front = t2->Nodes()[(fNodesFromEnd)]->Point3D();
102  TVector3 trk2Back = t2->Nodes()[t2->Nodes().size() - 1 - fNodesFromEnd]->Point3D();
103  TVector3 trk2FrontDir = (trk2Front - t2->Nodes()[(fNodesFromEnd + 1)]->Point3D()).Unit();
104  TVector3 trk2BackDir =
105  (trk2Back - t2->Nodes()[t2->Nodes().size() - 1 - (fNodesFromEnd + 1)]->Point3D()).Unit();
106 
107  // For stitching, we need to consider both ends of the track.
108  double offsetFront2 = GetTPCOffset(t2->FrontTPC(), t2->FrontCryo(), isCPA);
109  double offsetBack2 = GetTPCOffset(t2->BackTPC(), t2->BackCryo(), isCPA);
110 
111  // If the points to match are in the same TPC, then don't bother.
112  // Remember we have 4 points to consider here.
113  bool giveUp = false;
114  geo::TPCID tpc1(0, 0);
115  geo::TPCID tpc2(0, 0);
116  // Front-to-front
117  tpc1 = geo::TPCID(t1->FrontCryo(), t1->FrontTPC());
118  tpc2 = geo::TPCID(t2->FrontCryo(), t2->FrontTPC());
119  if (tpc1 == tpc2) giveUp = true;
120  // Front-to-back
121  tpc2 = geo::TPCID(t2->BackCryo(), t2->BackTPC());
122  if (tpc1 == tpc2) giveUp = true;
123  // Back-to-front
124  tpc1 = geo::TPCID(t1->BackCryo(), t1->BackTPC());
125  tpc2 = geo::TPCID(t2->FrontCryo(), t2->FrontTPC());
126  if (tpc1 == tpc2) giveUp = true;
127  // Back-to-back
128  tpc2 = geo::TPCID(t2->BackCryo(), t2->BackTPC());
129  if (tpc1 == tpc2) giveUp = true;
130 
131  // If the tracks have one end in the same TPC, give up.
132  if (giveUp) continue;
133 
134  // Also check that these tpcs do meet at the stitching surface (not a problem for protoDUNE).
135  double surfaceGap = 10.0;
136  bool carryOn[4] = {true, true, true, true};
137  if (fabs(offsetFront1 - offsetFront2) > surfaceGap) carryOn[0] = false;
138  if (fabs(offsetFront1 - offsetBack2) > surfaceGap) carryOn[1] = false;
139  if (fabs(offsetBack1 - offsetFront2) > surfaceGap) carryOn[2] = false;
140  if (fabs(offsetBack1 - offsetBack2) > surfaceGap) carryOn[3] = false;
141 
142  // Loop over the four options
143  for (int i = 0; i < 4; ++i) {
144 
145  if (!carryOn[i]) continue;
146 
147  TVector3 t1Pos;
148  TVector3 t2Pos;
149  TVector3 t1Dir;
150  TVector3 t2Dir;
151  double xShift1;
152  if (i < 2) {
153  t1Pos = trk1Front;
154  t1Dir = trk1FrontDir;
155  xShift1 = frontShift1;
156  }
157  else {
158  t1Pos = trk1Back;
159  t1Dir = trk1BackDir;
160  xShift1 = backShift1;
161  }
162  if (i % 2 == 0) {
163  t2Pos = trk2Front;
164  t2Dir = trk2FrontDir;
165  }
166  else {
167  t2Pos = trk2Back;
168  t2Dir = trk2BackDir;
169  }
170 
171  // Make sure the x directions point towards eachother (could be an issue for matching a short track)
172  if (t1Dir.X() * t2Dir.X() > 0) { continue; }
173 
174  double score = GetOptimalStitchShift(t1Pos, t2Pos, t1Dir, t2Dir, xShift1);
175 
176  if (score < fStitchingThreshold && score < bestMatchScore) {
177 
178  bestTrkMatch = t2;
179  xBestShift = xShift1;
180  bestMatchScore = score;
181  if (i < 2) { isBestFront1 = true; }
182  else {
183  isBestFront1 = false;
184  }
185  if (i % 2 == 0) { isBestFront2 = true; }
186  else {
187  isBestFront2 = false;
188  }
189  mf::LogInfo("pma::PMAlgStitcher")
190  << "Tracks " << t << " and " << u << " matching score = " << score << std::endl
191  << " - " << t1Pos.X() << ", " << t1Pos.Y() << ", " << t1Pos.Z() << " :: " << t1Dir.X()
192  << ", " << t1Dir.Y() << ", " << t1Dir.Z() << std::endl
193  << " - " << t2Pos.X() << ", " << t2Pos.Y() << ", " << t2Pos.Z() << " :: " << t2Dir.X()
194  << ", " << t2Dir.Y() << ", " << t2Dir.Z() << std::endl
195  << " - " << t1->FrontCryo() << ", " << t1->FrontTPC() << " :: " << t1->BackCryo()
196  << ", " << t1->BackTPC() << std::endl
197  << " - " << t2->FrontCryo() << ", " << t2->FrontTPC() << " :: " << t2->BackCryo()
198  << ", " << t2->BackTPC() << std::endl
199  << " - " << isBestFront1 << " :: " << isBestFront2 << std::endl;
200  } // End successful match if
201  } // Loop over matching options
202  } // Loop over track 2
203 
204  // If we found a match, do something about it.
205  if (bestTrkMatch != 0x0) {
206 
207  bool flip1 = false;
208  bool flip2 = false;
209  bool reverse = false;
210 
211  // Front-to-front match
212  if (isBestFront1 && isBestFront2) { flip1 = true; }
213  // Front-to-back match
214  else if (isBestFront1 && !isBestFront2) {
215  reverse = true;
216  }
217  // Back-to-back match
218  else if (!isBestFront1 && !isBestFront2) {
219  flip2 = true;
220  }
221  // Back-to-front match (do nothing)
222 
223  int tid1 = tracks.getCandidateTreeId(t1);
224  int tid2 = tracks.getCandidateTreeId(bestTrkMatch);
225 
226  bool canMerge = true;
227  if ((tid1 < 0) || (tid2 < 0)) {
228  throw cet::exception("pma::PMAlgStitching")
229  << "Track not found in the collection." << std::endl;
230  }
231  if (flip1) {
232 
233  std::vector<pma::Track3D*> newTracks;
234  if (t1->Flip(detProp, newTracks)) {
235  mf::LogInfo("pma::PMAlgStitching") << "Track 1 flipped.";
236  }
237  else {
238  mf::LogInfo("pma::PMAlgStitching") << "Unable to flip Track 1.";
239  canMerge = false;
240  }
241  for (const auto ts : newTracks) { // there may be a new track even if
242  // entire flip was not possible
243  tracks.tracks().emplace_back(ts, -1, tid1);
244  }
245  }
246  if (flip2) {
247 
248  std::vector<pma::Track3D*> newTracks;
249  if (bestTrkMatch->Flip(detProp, newTracks)) {
250  mf::LogInfo("pma::PMAlgStitching") << "Track 2 flipped.";
251  }
252  else {
253  mf::LogInfo("pma::PMAlgStitching") << "Unable to flip Track 1.";
254  canMerge = false;
255  }
256  for (const auto ts : newTracks) { // there may be a new track even if entire flip
257  // was not possible
258  tracks.tracks().emplace_back(ts, -1, tid2);
259  }
260  }
261 
262  t1->GetRoot()->ApplyDriftShiftInTree(clockData, detProp, -xBestShift);
263  bestTrkMatch->GetRoot()->ApplyDriftShiftInTree(clockData, detProp, +xBestShift);
264 
265  if (canMerge) {
266  mf::LogInfo("pma::PMAlgStitching") << "Merging tracks...";
267  int idx1 = tracks.getCandidateIndex(t1);
268  int idx2 = tracks.getCandidateIndex(bestTrkMatch);
269  if (reverse) // merge current track to another track, do not increase the outer
270  // loop index t (next after current track jumps in at t)
271  {
272  if (tracks.setTreeOriginAtFront(detProp, t1)) {
273  tracks.merge((size_t)idx2, (size_t)idx1);
274  }
275  else {
276  mf::LogWarning("pma::PMAlgStitching") << " could not merge.";
277  ++t;
278  }
279  }
280  else // merge to the current track, do not increase the outer loop index t (maybe
281  // something else will match to the extended track)
282  {
283  if (tracks.setTreeOriginAtFront(detProp, bestTrkMatch)) {
284  tracks.merge((size_t)idx1, (size_t)idx2);
285  }
286  else {
287  mf::LogWarning("pma::PMAlgStitching") << " could not merge.";
288  ++t;
289  }
290  }
291  mf::LogInfo("pma::PMAlgStitching") << "...done";
292  }
293  else {
294  ++t;
295  } // track matched, but not merged, go to the next track in the outer loop
296  }
297  else {
298  ++t;
299  } // no track matched, go to the next track in the outer loop
300  }
301 }
size_t size() const
double GetTPCOffset(unsigned int tpc, unsigned int cryo, bool isCPA)
TTree * t1
Definition: plottest35.C:26
int getCandidateTreeId(pma::Track3D const *candidate) const
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void ApplyDriftShiftInTree(const detinfo::DetectorClocksData &clockData, detinfo::DetectorPropertiesData const &detProp, double dx, bool skipFirst=false)
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
unsigned int BackTPC() const
Definition: PmaTrack3D.h:121
unsigned int fNodesFromEnd
unsigned int BackCryo() const
Definition: PmaTrack3D.h:122
double GetOptimalStitchShift(TVector3 &pos1, TVector3 &pos2, TVector3 &dir1, TVector3 &dir2, double &shift) const
bool setTreeOriginAtFront(detinfo::DetectorPropertiesData const &detProp, pma::Track3D *trk)
pma::Track3D * GetRoot()
int getCandidateIndex(pma::Track3D const *candidate) const
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:118
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:119
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
void merge(size_t idx1, size_t idx2)
TTree * t2
Definition: plottest35.C:36
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:527
std::vector< TrkCandidate > const & tracks() const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void pma::PMAlgStitching::StitchTracksAPA ( const detinfo::DetectorClocksData clockData,
const detinfo::DetectorPropertiesData detProp,
pma::TrkCandidateColl tracks 
)

Definition at line 42 of file PMAlgStitching.cxx.

References pma::TrkCandidateColl::size(), and StitchTracks().

Referenced by pma::PMAlgTracker::build().

45 {
46  mf::LogInfo("pma::PMAlgStitching") << "Passed " << tracks.size() << " tracks for APA stitching.";
47  StitchTracks(clockData, detProp, tracks, false);
48 }
size_t size() const
void StitchTracks(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &tracks, bool isCPA)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void pma::PMAlgStitching::StitchTracksCPA ( const detinfo::DetectorClocksData clockData,
const detinfo::DetectorPropertiesData detProp,
pma::TrkCandidateColl tracks 
)

Definition at line 33 of file PMAlgStitching.cxx.

References pma::TrkCandidateColl::size(), and StitchTracks().

Referenced by pma::PMAlgTracker::build().

36 {
37  mf::LogInfo("pma::PMAlgStitching") << "Passed " << tracks.size() << " tracks for CPA stitching.";
38  StitchTracks(clockData, detProp, tracks, true);
39 }
size_t size() const
void StitchTracks(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, pma::TrkCandidateColl &tracks, bool isCPA)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo

Member Data Documentation

unsigned int pma::PMAlgStitching::fNodesFromEnd
private

Definition at line 81 of file PMAlgStitching.h.

Referenced by PMAlgStitching(), and StitchTracks().

double pma::PMAlgStitching::fStitchingThreshold
private

Definition at line 79 of file PMAlgStitching.h.

Referenced by PMAlgStitching(), and StitchTracks().

std::map<geo::TPCID, double> pma::PMAlgStitching::fTPCXOffsetsAPA
private

Definition at line 75 of file PMAlgStitching.h.

Referenced by GetTPCOffset(), and GetTPCXOffsets().

std::map<geo::TPCID, double> pma::PMAlgStitching::fTPCXOffsetsCPA
private

Definition at line 76 of file PMAlgStitching.h.

Referenced by GetTPCOffset(), and GetTPCXOffsets().


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