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

Constructor & Destructor Documentation

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

Definition at line 20 of file PMAlgStitching.cxx.

21 {
22 
23  // Set parameters from the config.
25  fNodesFromEnd = config.NodesFromEnd();
26 
27  // Get CPA and APA positions.
29 }
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 311 of file PMAlgStitching.cxx.

316 {
317 
318  double stepSize = 0.1;
319  double minShift = shift - (50. * stepSize);
320  double maxShift = shift + (50. * stepSize);
321  double bestShift = 99999;
322  double bestScore = 99999;
323 
324  for (shift = minShift; shift <= maxShift; shift += stepSize) {
325  TVector3 newPos1 = pos1;
326  TVector3 newPos2 = pos2;
327  newPos1.SetX(pos1.X() - shift);
328  newPos2.SetX(pos2.X() + shift);
329  double thisScore = GetTrackPairDelta(newPos1, newPos2, dir1, dir2);
330  if (thisScore < bestScore) {
331  bestShift = shift;
332  bestScore = thisScore;
333  }
334  }
335 
336  shift = bestShift;
337  return bestScore;
338 }
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 411 of file PMAlgStitching.cxx.

412 {
413 
414  geo::TPCID thisTPCID(cryo, tpc);
415  double offset = 0.0;
416  if (isCPA) { offset = fTPCXOffsetsCPA[thisTPCID]; }
417  else {
418  offset = fTPCXOffsetsAPA[thisTPCID];
419  }
420  return offset;
421 }
std::map< geo::TPCID, double > fTPCXOffsetsAPA
std::map< geo::TPCID, double > fTPCXOffsetsCPA
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
void pma::PMAlgStitching::GetTPCXOffsets ( )
private

Definition at line 366 of file PMAlgStitching.cxx.

367 {
368 
369  // Grab hold of the geometry
370  auto const* geom = lar::providerFrom<geo::Geometry>();
371 
372  // Loop over each TPC and store the require information
373  for (geo::TPCID const& tID : geom->IterateTPCIDs()) {
374 
375  geo::TPCGeo const& aTPC = geom->TPC(tID);
376 
377  // Loop over the 3 possible readout planes to find the x position
378  unsigned int plane = 0;
379  bool hasPlane = false;
380  for (; plane < 4; ++plane) {
381  hasPlane = aTPC.HasPlane(plane);
382  if (hasPlane) { break; }
383  }
384 
385  if (!hasPlane) { continue; }
386 
387  // Get the x position of the readout plane
388  double xAnode = aTPC.PlaneLocation(plane)[0];
389  fTPCXOffsetsAPA.insert(std::make_pair(tID, xAnode));
390 
391  // For the cathode, we have to try a little harder. Firstly, find the
392  // min and max x values for the TPC.
393  double origin[3] = {0.};
394  double center[3] = {0.};
395  aTPC.LocalToWorld(origin, center);
396  double tpcDim[3] = {aTPC.HalfWidth(), aTPC.HalfHeight(), 0.5 * aTPC.Length()};
397  double xmin = center[0] - tpcDim[0];
398  double xmax = center[0] + tpcDim[0];
399  double xCathode = 0.;
400  // Now check which is further from the APA and use that as the cathode value.
401  if (fabs(xmin - xAnode) > fabs(xmax - xAnode)) { xCathode = xmin; }
402  else {
403  xCathode = xmax;
404  }
405  fTPCXOffsetsCPA.insert(std::make_pair(tID, xCathode));
406  }
407 }
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
Definition: TPCGeo.h:175
Geometry information for a single TPC.
Definition: TPCGeo.h:38
std::map< geo::TPCID, double > fTPCXOffsetsAPA
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
std::map< geo::TPCID, double > fTPCXOffsetsCPA
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
def center(depos, point)
Definition: depos.py:117
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
const double * PlaneLocation(unsigned int p) const
Definition: TPCGeo.cxx:382
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:107
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
double pma::PMAlgStitching::GetTrackPairDelta ( TVector3 &  pos1,
TVector3 &  pos2,
TVector3 &  dir1,
TVector3 &  dir2 
) const
private

Definition at line 342 of file PMAlgStitching.cxx.

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

Definition at line 55 of file PMAlgStitching.cxx.

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

Definition at line 43 of file PMAlgStitching.cxx.

46 {
47  mf::LogInfo("pma::PMAlgStitching") << "Passed " << tracks.size() << " tracks for APA stitching.";
48  StitchTracks(clockData, detProp, tracks, false);
49 }
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.

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 79 of file PMAlgStitching.h.

double pma::PMAlgStitching::fStitchingThreshold
private

Definition at line 77 of file PMAlgStitching.h.

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

Definition at line 73 of file PMAlgStitching.h.

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

Definition at line 74 of file PMAlgStitching.h.


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