LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GeometryUtilities.cxx
Go to the documentation of this file.
1 // \file GeometryUtilities.cxx
3 //
4 // \brief Functions to calculate distances and angles in 3D and 2D
5 //
6 // andrzej.szelc@yale.edu
7 //
9 
18 
19 #include "cetlib/pow.h"
20 
21 #include "TLorentzVector.h"
22 
23 #include <cmath>
24 
25 namespace util {
26 
28  detinfo::DetectorClocksData const& clockData,
29  detinfo::DetectorPropertiesData const& propData)
30  : fGeom{geom}, fClocks{clockData}, fDetProp{propData}
31  {
33  vertangle.resize(fNPlanes);
34  for (unsigned int ip = 0; ip < fNPlanes; ip++) {
35  geo::WireID const wid{0, 0, ip, 0};
36  vertangle[ip] = fGeom.Wire(wid).ThetaZ(false) - TMath::Pi() / 2.; // wire angle
37  }
38 
40  fTimeTick = sampling_rate(fClocks) / 1000.;
42 
46  }
47 
48  //-----------------------------------------------------------------------------
49  // omega0 and omega1 (calculated by CPAN in degrees):
50  // angle based on distances in wires and time - rescaled to cm.
51  // tan(angle)*fMean_wire_pitch/(fTimeTick*fDriftVelocity);
52  // as those calculated with Get2Dangle
53  // writes phi and theta in degrees.
56  int iplane1,
57  double omega0,
58  double omega1,
59  double& phi,
60  double& theta) const
61  {
62  // y, z, x coordinates
63  double ln(0), mn(0), nn(0);
64  double phis(0), thetan(0);
65 
66  // Pretend collection and induction planes.
67  // "Collection" is the plane with the vertical angle equal to zero.
68  // If both are non-zero, collection is the one with the negative angle.
69  unsigned int Cplane = 0, Iplane = 1;
70 
71  // angleC and angleI are the respective angles to vertical in C/I
72  // planes and slopeC, slopeI are the tangents of the track.
73  double angleC, angleI, slopeC, slopeI, omegaC, omegaI;
74  omegaC = kINVALID_DOUBLE;
75  omegaI = kINVALID_DOUBLE;
76 
77  // Don't know how to reconstruct these yet, so exit with error.
78  // In
79  if (omega0 == 0 || omega1 == 0) {
80  phi = 0;
81  theta = -999;
82  return -1;
83  }
84 
86 
87  // check if backwards going track
88  double alt_backwards = 0;
89 
90  if (std::fabs(omega0) > (TMath::Pi() / 2.0) || std::fabs(omega1) > (TMath::Pi() / 2.0)) {
91  alt_backwards = 1;
92  }
93 
94  if (vertangle[iplane0] == 0) {
95  // first plane is at 0 degrees
96  Cplane = iplane0;
97  Iplane = iplane1;
98  omegaC = omega0;
99  omegaI = omega1;
100  }
101  else if (vertangle[iplane1] == 0) {
102  // second plane is at 0 degrees
103  Cplane = iplane1;
104  Iplane = iplane0;
105  omegaC = omega1;
106  omegaI = omega0;
107  }
108  else if (vertangle[iplane0] != 0 && vertangle[iplane1] != 0) {
109  // both planes are at non zero degree - find the one with deg<0
110  if (vertangle[iplane1] < vertangle[iplane0]) {
111  Cplane = iplane1;
112  Iplane = iplane0;
113  omegaC = omega1;
114  omegaI = omega0;
115  }
116  else if (vertangle[iplane1] > vertangle[iplane0]) {
117  Cplane = iplane0;
118  Iplane = iplane1;
119  omegaC = omega0;
120  omegaI = omega1;
121  }
122  else {
123  // throw error - same plane.
124  return -1;
125  }
126  }
127 
128  slopeC = std::tan(omegaC);
129  slopeI = std::tan(omegaI);
130  angleC = vertangle[Cplane];
131  angleI = vertangle[Iplane];
132 
133  // 0 -1 factor depending on if one of the planes is vertical.
134  bool nfact = !(vertangle[Cplane]);
135 
136  // ln represents y, omega is 2d angle -- in first 2 quadrants y is positive.
137  if (omegaC < TMath::Pi() && omegaC > 0)
138  ln = 1;
139  else
140  ln = -1;
141 
142  // calculate x and z using y ( ln )
143  mn = (ln / (2 * std::sin(angleI))) *
144  ((std::cos(angleI) / (slopeC * std::cos(angleC))) - (1 / slopeI) +
145  nfact * (std::cos(angleI) / (std::cos(angleC) * slopeC) - 1 / slopeI));
146 
147  nn = (ln / (2 * std::cos(angleC))) *
148  ((1 / slopeC) + (1 / slopeI) + nfact * ((1 / slopeC) - (1 / slopeI)));
149 
150  // Direction angles
151  if (std::fabs(omegaC) > 0.01) // catch numeric error values
152  {
153  // phi=std::atan(ln/nn);
154  phis = std::asin(ln / TMath::Sqrt(ln * ln + nn * nn));
155 
156  if (std::fabs(slopeC + slopeI) < 0.001)
157  phis = 0;
158  else if (std::fabs(omegaI) > 0.01 &&
159  (omegaI / std::fabs(omegaI) == -omegaC / std::fabs(omegaC)) &&
160  (std::fabs(omegaC) < 1 * TMath::Pi() / 180 ||
161  std::fabs(omegaC) > 179 * TMath::Pi() / 180)) // angles have
162  phis = (std::fabs(omegaC) > TMath::Pi() / 2) ? TMath::Pi() : 0; // angles are
163 
164  if (nn < 0 && phis > 0 &&
165  !(!alt_backwards &&
166  std::fabs(phis) <
167  TMath::Pi() / 4)) // do not go back if track looks forward and phi is forward
168  phis = (TMath::Pi()) - phis;
169  else if (nn < 0 && phis < 0 && !(!alt_backwards && std::fabs(phis) < TMath::Pi() / 4))
170  phis = (-TMath::Pi()) - phis;
171 
172  phi = phis * 180 / TMath::Pi();
173  }
174  // If plane2 (collection), phi = 2d angle (omegaC in this case)
175  else {
176  phis = omegaC;
177  phi = omegaC;
178  }
179 
180  thetan = -std::asin(mn / std::hypot(ln, mn, nn));
181  theta = thetan * 180 / TMath::Pi();
182 
183  return 0;
184  }
185 
187  // Calculate theta in case phi~0
188  // returns theta in angles
191  int iplane1,
192  double dw0,
193  double dw1) const
194  {
195 
196  double splane, lplane; // plane in which the track is shorter and longer.
197  double sdw, ldw;
198 
199  if (dw0 == 0 && dw1 == 0) return -1;
200 
201  if (dw0 >= dw1) {
202  lplane = iplane0;
203  splane = iplane1;
204  ldw = dw0;
205  sdw = dw1;
206  }
207  else {
208  lplane = iplane1;
209  splane = iplane0;
210  ldw = dw1;
211  sdw = dw0;
212  }
213 
214  double top = (std::cos(vertangle[splane]) - std::cos(vertangle[lplane]) * sdw / ldw);
215  double bottom = std::tan(vertangle[lplane] * std::cos(vertangle[splane]));
216  bottom -= std::tan(vertangle[splane] * std::cos(vertangle[lplane])) * sdw / ldw;
217 
218  double tantheta = top / bottom;
219 
220  return std::atan(tantheta) * vertangle[lplane] / std::abs(vertangle[lplane]) * 180. /
221  TMath::Pi();
222  }
223 
225  // Calculate 3D pitch in beam coordinates
226  //
228  double GeometryUtilities::CalculatePitch(unsigned int iplane, double phi, double theta) const
229  {
230  auto const& plane = fGeom.Plane({0, 0, iplane});
231  if (plane.View() == geo::kUnknown || plane.View() == geo::k3D) {
232  mf::LogError(Form("Warning : no Pitch foreseen for view %d", plane.View()));
233  return -1.;
234  }
235 
236  double const pi = TMath::Pi();
237  double const fTheta = pi / 2 - theta;
238  double const fPhi = -(phi + pi / 2);
239 
240  double wirePitch = plane.WirePitch();
241  double angleToVert = 0.5 * TMath::Pi() - fGeom.WireAngleToVertical(plane.View(), plane.ID());
242  double cosgamma = std::abs(std::sin(angleToVert) * std::cos(fTheta) +
243  std::cos(angleToVert) * std::sin(fTheta) * std::sin(fPhi));
244 
245  return cosgamma > 0 ? wirePitch / cosgamma : -1.;
246  }
247 
249  // Calculate 3D pitch in polar coordinates
250  //
252  double GeometryUtilities::CalculatePitchPolar(unsigned int iplane, double phi, double theta) const
253  {
254  auto const& plane = fGeom.Plane({0, 0, iplane});
255  if (plane.View() == geo::kUnknown || plane.View() == geo::k3D) {
256  mf::LogError(Form("Warning : no Pitch foreseen for view %d", plane.View()));
257  return -1.;
258  }
259 
260  double const fTheta = theta; // KJK: Are these two variables necessary?
261  double const fPhi = phi;
262 
263  double wirePitch = plane.WirePitch();
264  double angleToVert = 0.5 * TMath::Pi() - fGeom.WireAngleToVertical(plane.View(), plane.ID());
265  double cosgamma = std::abs(std::sin(angleToVert) * std::cos(fTheta) +
266  std::cos(angleToVert) * std::sin(fTheta) * std::sin(fPhi));
267 
268  return cosgamma > 0 ? wirePitch / cosgamma : -1.;
269  }
270 
272  // Calculate 2D slope
273  // in "cm" "cm" coordinates
275  double GeometryUtilities::Get2Dslope(double wireend,
276  double wirestart,
277  double timeend,
278  double timestart) const
279  {
280 
281  return GeometryUtilities::Get2Dslope((wireend - wirestart) * fWiretoCm,
282  (timeend - timestart) * fTimetoCm);
283  }
284 
286  // Calculate 2D slope
287  // in "cm" "cm" coordinates
289  double GeometryUtilities::Get2Dslope(const PxPoint* endpoint, const PxPoint* startpoint) const
290  {
291  return Get2Dslope(endpoint->w - startpoint->w, endpoint->t - startpoint->t);
292  }
293 
295  // Calculate 2D slope
296  // in wire time coordinates coordinates
297  //
299  double GeometryUtilities::Get2Dslope(double dwire, double dtime) const
300  {
301  return std::tan(Get2Dangle(dwire, dtime)) / fWireTimetoCmCm;
302  }
303 
305  // Calculate 2D angle
306  // in "cm" "cm" coordinates
308  double GeometryUtilities::Get2Dangle(double wireend,
309  double wirestart,
310  double timeend,
311  double timestart) const
312  {
313  return Get2Dangle((wireend - wirestart) * fWiretoCm, (timeend - timestart) * fTimetoCm);
314  }
315 
317  // Calculate 2D angle
318  // in "cm" "cm" coordinates, endpoint and startpoint are assumed to be in
319  // cm/cm space
321  double GeometryUtilities::Get2Dangle(const PxPoint* endpoint, const PxPoint* startpoint) const
322  {
323  return Get2Dangle(endpoint->w - startpoint->w, endpoint->t - startpoint->t);
324  }
325 
327  // Calculate 2D angle
328  // in "cm" "cm" coordinates
330  double GeometryUtilities::Get2Dangle(double dwire, double dtime) const
331  {
332  double BC, AC;
333  double omega;
334 
335  BC = ((double)dwire); // in cm
336  AC = ((double)dtime); // in cm
337  omega = std::asin(AC / std::hypot(AC, BC));
338  if (BC < 0) // for the time being. Will check if it works for AC<0
339  {
340  if (AC > 0) {
341  omega = TMath::Pi() - std::abs(omega); //
342  }
343  else if (AC < 0) {
344  omega = -TMath::Pi() + std::abs(omega);
345  }
346  else {
347  omega = TMath::Pi();
348  }
349  }
350  return omega;
351  }
352 
353  // accepting phi and theta in degrees
354  // returning in radians.
355 
356  double GeometryUtilities::Get2DangleFrom3D(unsigned int plane, double phi, double theta) const
357  {
358  TVector3 dummyvector(std::cos(theta * TMath::Pi() / 180.) * std::sin(phi * TMath::Pi() / 180.),
359  std::sin(theta * TMath::Pi() / 180.),
360  std::cos(theta * TMath::Pi() / 180.) * std::cos(phi * TMath::Pi() / 180.));
361  return Get2DangleFrom3D(plane, dummyvector);
362  }
363 
364  // accepting TVector3
365  // returning in radians as is customary,
366 
367  double GeometryUtilities::Get2DangleFrom3D(unsigned int plane, TVector3 dir_vector) const
368  {
369  geo::PlaneID const planeid{0, 0, plane};
370  double alpha =
371  0.5 * TMath::Pi() - fGeom.WireAngleToVertical(fGeom.Plane(planeid).View(), planeid);
372  // create dummy xyz point in middle of detector and another one in unit
373  // length. calculate correspoding points in wire-time space and use the
374  // differnces between those to return 2D a angle
375 
376  TVector3 start(fGeom.DetHalfWidth(), 0., fGeom.DetLength() / 2.);
377  TVector3 end = start + dir_vector;
378 
379  // the wire coordinate is already in cm. The time needs to be converted.
380  PxPoint startp(plane,
381  (fGeom.DetHalfHeight() * std::sin(std::fabs(alpha)) +
382  start[2] * std::cos(alpha) - start[1] * std::sin(alpha)),
383  start[0]);
384 
385  PxPoint endp(plane,
386  (fGeom.DetHalfHeight() * std::sin(std::fabs(alpha)) + end[2] * std::cos(alpha) -
387  end[1] * std::sin(alpha)),
388  end[0]);
389 
390  double angle = Get2Dangle(&endp, &startp);
391 
392  return angle;
393  }
394 
396  // Calculate 2D distance
397  // in "cm" "cm" coordinates
399  double GeometryUtilities::Get2DDistance(double wire1,
400  double time1,
401  double wire2,
402  double time2) const
403  {
404  return std::hypot((wire1 - wire2) * fWiretoCm, (time1 - time2) * fTimetoCm);
405  }
406 
407  double GeometryUtilities::Get2DDistance(const PxPoint* point1, const PxPoint* point2) const
408  {
409  return std::hypot(point1->w - point2->w, point1->t - point2->t);
410  }
411 
413  // Calculate 2D distance, using 2D angle
414  // in "cm" "cm" coordinates
416  double GeometryUtilities::Get2DPitchDistance(double angle, double inwire, double wire) const
417  {
418  double radangle = TMath::Pi() * angle / 180;
419  if (std::cos(radangle) == 0) return 9999;
420  return std::abs((wire - inwire) / std::cos(radangle)) * fWiretoCm;
421  }
422 
423  //----------------------------------------------------------------------------
424  double GeometryUtilities::Get2DPitchDistanceWSlope(double slope, double inwire, double wire) const
425  {
426 
427  return std::abs(wire - inwire) * TMath::Sqrt(1 + slope * slope) * fWiretoCm;
428  }
429 
431  // Calculate wire,time coordinates of the Hit projection onto a line
432  //
435  double intercept,
436  double wire1,
437  double time1,
438  double& wireout,
439  double& timeout) const
440  {
441  double invslope = 0;
442 
443  if (slope) { invslope = -1. / slope; }
444 
445  double ort_intercept = time1 - invslope * (double)wire1;
446 
447  if ((slope - invslope) != 0)
448  wireout = (ort_intercept - intercept) / (slope - invslope);
449  else
450  wireout = wire1;
451  timeout = slope * wireout + intercept;
452 
453  return 0;
454  }
455 
457  // Calculate wire,time coordinates of the Hit projection onto a line
458  // all points are assumed to be in cm/cm space.
461  const PxPoint* startpoint,
462  const PxPoint* point1,
463  PxPoint& pointout) const
464  {
465 
466  double intercept = startpoint->t - slope * startpoint->w;
467 
468  return GetPointOnLine(slope, intercept, point1, pointout);
469  }
470 
472  // Calculate wire,time coordinates of the Hit projection onto a line
473  // all points assumed to be in cm/cm space.
476  double intercept,
477  const PxPoint* point1,
478  PxPoint& pointout) const
479  {
480  double invslope = 0;
481 
482  if (slope) { invslope = -1. / slope; }
483 
484  double ort_intercept = point1->t - invslope * point1->w;
485 
486  if ((slope - invslope) != 0)
487  pointout.w = (ort_intercept - intercept) / (slope - invslope);
488  else
489  pointout.w = point1->w;
490 
491  pointout.t = slope * pointout.w + intercept;
492 
493  return 0;
494  }
495 
497  // Calculate wire,time coordinates of the Hit projection onto a line
498  //
501  double wirestart,
502  double timestart,
503  double wire1,
504  double time1,
505  double& wireout,
506  double& timeout) const
507  {
508  double intercept = timestart - slope * (double)wirestart;
509 
510  return GetPointOnLine(slope, intercept, wire1, time1, wireout, timeout);
511  }
512 
514  // Calculate wire,time coordinates of the Hit projection onto a line
515  //
517  double intercept,
518  double ort_intercept,
519  double& wireout,
520  double& timeout) const
521  {
522  double invslope = 0;
523 
524  if (slope) { invslope = -1. / slope; }
525 
526  invslope *= fWireTimetoCmCm * fWireTimetoCmCm;
527 
528  wireout = (ort_intercept - intercept) / (slope - invslope);
529  timeout = slope * wireout + intercept;
530 
531  wireout /= fWiretoCm;
532  timeout /= fTimetoCm;
533 
534  return 0;
535  }
536 
538  // Calculate wire,time coordinates of the Hit projection onto a line
539  // slope should be in cm/cm space. PxPoint should be in cm/cm space.
542  double intercept,
543  double ort_intercept,
544  PxPoint& pointonline) const
545  {
546  double invslope = 0;
547 
548  if (slope) invslope = -1. / slope;
549 
550  pointonline.w = (ort_intercept - intercept) / (slope - invslope);
551  pointonline.t = slope * pointonline.w + intercept;
552  return 0;
553  }
554 
556  int GeometryUtilities::GetProjectedPoint(const PxPoint* p0, const PxPoint* p1, PxPoint& pN) const
557  {
558  // determine third plane number
559  for (unsigned int i = 0; i < fNPlanes; ++i) {
560  if (i == p0->plane || i == p1->plane) continue;
561  pN.plane = i;
562  }
563 
564  // Assuming there is no problem ( and we found the best pair that comes
565  // close in time ) we try to get the Y and Z coordinates for the start of
566  // the shower.
567  unsigned int chan1 = fGeom.PlaneWireToChannel(geo::WireID(0, 0, p0->plane, p0->w));
568  unsigned int chan2 = fGeom.PlaneWireToChannel(geo::WireID(0, 0, p1->plane, p1->w));
570  auto pos = fGeom.Plane(geo::PlaneID(0, 0, p0->plane)).toWorldCoords(origin);
571  double x = (p0->t - trigger_offset(fClocks)) * fTimetoCm + pos.X();
572 
573  double y, z;
574  if (!fGeom.ChannelsIntersect(chan1, chan2, y, z)) return -1;
575 
576  pos.SetCoordinates(x, y, z);
577 
578  pN = Get2DPointProjection(pos, pN.plane);
579 
580  return 0;
581  }
582 
584  int GeometryUtilities::GetYZ(const PxPoint* p0, const PxPoint* p1, double* yz) const
585  {
586  assert(p0 && p1);
587  geo::TPCID const tpcid{0, 0};
588  geo::PlaneID const plane_0{tpcid, p0->plane};
589  geo::PlaneID const plane_1{tpcid, p1->plane};
590 
591  double y, z;
592 
593  // Force to the closest wires if not in the range
594  int z0 = p0->w / fWiretoCm;
595  int z1 = p1->w / fWiretoCm;
596  if (z0 < 0) {
597  std::cout << "\033[93mWarning\033[00m "
598  "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
599  << std::endl
600  << " 2D wire position " << p0->w << " [cm] corresponds to negative wire number."
601  << std::endl
602  << " Forcing it to wire=0..." << std::endl
603  << "\033[93mWarning ends...\033[00m" << std::endl;
604  z0 = 0;
605  }
606  else if (z0 >= (int)(fGeom.Nwires(plane_0))) {
607  std::cout << "\033[93mWarning\033[00m "
608  "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
609  << std::endl
610  << " 2D wire position " << p0->w << " [cm] exceeds max wire number "
611  << (fGeom.Nwires(plane_0) - 1) << std::endl
612  << " Forcing it to the max wire number..." << std::endl
613  << "\033[93mWarning ends...\033[00m" << std::endl;
614  z0 = fGeom.Nwires(plane_0) - 1;
615  }
616  if (z1 < 0) {
617  std::cout << "\033[93mWarning\033[00m "
618  "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
619  << std::endl
620  << " 2D wire position " << p1->w << " [cm] corresponds to negative wire number."
621  << std::endl
622  << " Forcing it to wire=0..." << std::endl
623  << "\033[93mWarning ends...\033[00m" << std::endl;
624  z1 = 0;
625  }
626  if (z1 >= (int)(fGeom.Nwires(plane_1))) {
627  std::cout << "\033[93mWarning\033[00m "
628  "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
629  << std::endl
630  << " 2D wire position " << p1->w << " [cm] exceeds max wire number "
631  << (fGeom.Nwires(plane_1) - 1) << std::endl
632  << " Forcing it to the max wire number..." << std::endl
633  << "\033[93mWarning ends...\033[00m" << std::endl;
634  z1 = fGeom.Nwires(plane_1) - 1;
635  }
636 
637  unsigned int chan1 = fGeom.PlaneWireToChannel(geo::WireID(plane_0, z0));
638  unsigned int chan2 = fGeom.PlaneWireToChannel(geo::WireID(plane_1, z1));
639 
640  if (!fGeom.ChannelsIntersect(chan1, chan2, y, z)) return -1;
641 
642  yz[0] = y;
643  yz[1] = z;
644 
645  return 0;
646  }
647 
649  int GeometryUtilities::GetXYZ(const PxPoint* p0, const PxPoint* p1, double* xyz) const
650  {
652  auto const pos = fGeom.Plane(geo::PlaneID{0, 0, p0->plane}).toWorldCoords(origin);
653  double x = (p0->t) - trigger_offset(fClocks) * fTimetoCm + pos.X();
654  double yz[2];
655 
656  GetYZ(p0, p1, yz);
657 
658  xyz[0] = x;
659  xyz[1] = yz[0];
660  xyz[2] = yz[1];
661 
662  return 0;
663  }
664 
666 
667  PxPoint GeometryUtilities::Get2DPointProjection(geo::Point_t const& xyz, unsigned int plane) const
668  {
669  geo::PlaneID const planeID{0, 0, plane};
670  PxPoint pN(0, 0, 0);
672  auto pos = fGeom.Plane(planeID).toWorldCoords(origin);
673  double drifttick = (xyz.X() / fDriftVelocity) * (1. / fTimeTick);
674 
675  pos.SetY(xyz.Y());
676  pos.SetZ(xyz.Z());
677 
680 
681  pN.w = fGeom.NearestWireID(pos, planeID).Wire;
682  pN.t = drifttick - (pos.X() / fDriftVelocity) * (1. / fTimeTick) + trigger_offset(fClocks);
683  pN.plane = plane;
684 
685  return pN;
686  }
687 
689  // for now this returns the vlause in CM/CM space.
690  // this will become the default, but don't want to break the code that depends
691  // on the previous version. A.S. 03/26/14
693 
694  PxPoint GeometryUtilities::Get2DPointProjectionCM(std::vector<double> const& xyz,
695  unsigned int plane) const
696  {
697 
698  PxPoint pN(0, 0, 0);
699 
700  geo::Point_t const pos{0., xyz[1], xyz[2]};
701 
704 
705  return {plane, fGeom.NearestWireID(pos, geo::PlaneID{0, 0, plane}).Wire * fWiretoCm, xyz[0]};
706  }
707 
708  PxPoint GeometryUtilities::Get2DPointProjectionCM(double const* xyz, unsigned int plane) const
709  {
710  geo::Point_t const pos{0., xyz[1], xyz[2]};
711 
714 
715  return {plane, fGeom.NearestWireID(pos, geo::PlaneID{0, 0, plane}).Wire * fWiretoCm, xyz[0]};
716  }
717 
719  unsigned int plane) const
720  {
721  double xyznew[3] = {(*xyz)[0], (*xyz)[1], (*xyz)[2]};
722 
723  return Get2DPointProjectionCM(xyznew, plane);
724  }
725 
726  double GeometryUtilities::GetTimeTicks(double x, unsigned int plane) const
727  {
729  auto const pos = fGeom.Plane(geo::PlaneID{0, 0, plane}).toWorldCoords(origin);
730  double drifttick = (x / fDriftVelocity) * (1. / fTimeTick);
731 
732  return drifttick - (pos.X() / fDriftVelocity) * (1. / fTimeTick) + trigger_offset(fClocks);
733  }
734 
735  //----------------------------------------------------------------------
736  // provide projected wire pitch for the view // copied from track.cxx and
737  // modified
738  double GeometryUtilities::PitchInView(unsigned int plane, double phi, double theta) const
739  {
740  double dirs[3] = {0.};
741  GetDirectionCosines(phi, theta, dirs);
742 
745  double wirePitch = 0.;
746  double angleToVert = 0.;
747 
748  auto const& planegeom = fGeom.Plane({0, 0, plane});
749  wirePitch = planegeom.WirePitch();
750  angleToVert = fGeom.WireAngleToVertical(planegeom.View(), planegeom.ID()) - 0.5 * TMath::Pi();
751 
752  //(sin(angleToVert),std::cos(angleToVert)) is the direction perpendicular to
753  // wire fDir.front() is the direction of the track at the beginning of its
754  // trajectory
755  double cosgamma = std::abs(std::sin(angleToVert) * dirs[1] + std::cos(angleToVert) * dirs[2]);
756 
757  if (cosgamma < 1.e-5)
758  // throw UtilException("cosgamma is basically 0, that can't be right");
759  {
760  std::cout << " returning 100" << std::endl;
761  return 100;
762  }
763 
764  return wirePitch / cosgamma;
765  }
766 
768  void GeometryUtilities::GetDirectionCosines(double phi, double theta, double* dirs) const
769  {
770  theta *= (TMath::Pi() / 180);
771  phi *= (TMath::Pi() / 180); // working on copies, it's ok.
772  dirs[0] = std::cos(theta) * std::sin(phi);
773  dirs[1] = std::sin(theta);
774  dirs[2] = std::cos(theta) * std::cos(phi);
775  }
776 
777  void GeometryUtilities::SelectLocalHitlist(const std::vector<PxHit>& hitlist,
778  std::vector<const PxHit*>& hitlistlocal,
779  PxPoint& startHit,
780  double& linearlimit,
781  double& ortlimit,
782  double& lineslopetest) const
783  {
784  PxHit testHit;
786  hitlist, hitlistlocal, startHit, linearlimit, ortlimit, lineslopetest, testHit);
787  }
788 
792 
793  void GeometryUtilities::SelectLocalHitlist(const std::vector<PxHit>& hitlist,
794  std::vector<const PxHit*>& hitlistlocal,
795  PxPoint& startHit,
796  double& linearlimit,
797  double& ortlimit,
798  double& lineslopetest,
799  PxHit& averageHit) const
800  {
801 
802  hitlistlocal.clear();
803  std::vector<unsigned int> hitlistlocal_index;
804 
805  hitlistlocal_index.clear();
806 
808  hitlist, hitlistlocal_index, startHit, linearlimit, ortlimit, lineslopetest);
809 
810  double timesum = 0;
811  double wiresum = 0;
812  for (size_t i = 0; i < hitlistlocal_index.size(); ++i) {
813 
814  hitlistlocal.push_back((const PxHit*)(&(hitlist.at(hitlistlocal_index.at(i)))));
815  timesum += hitlist.at(hitlistlocal_index.at(i)).t;
816  wiresum += hitlist.at(hitlistlocal_index.at(i)).w;
817  }
818 
819  averageHit.plane = startHit.plane;
820  if (hitlistlocal.size()) {
821  averageHit.w = wiresum / (double)hitlistlocal.size();
822  averageHit.t = timesum / ((double)hitlistlocal.size());
823  }
824  }
825 
826  void GeometryUtilities::SelectLocalHitlistIndex(const std::vector<PxHit>& hitlist,
827  std::vector<unsigned int>& hitlistlocal_index,
828  PxPoint& startHit,
829  double& linearlimit,
830  double& ortlimit,
831  double& lineslopetest) const
832  {
833 
834  hitlistlocal_index.clear();
835  double locintercept = startHit.t - startHit.w * lineslopetest;
836 
837  for (size_t i = 0; i < hitlist.size(); ++i) {
838 
839  PxPoint hitonline;
840 
841  GetPointOnLine(lineslopetest, locintercept, (const PxHit*)(&hitlist.at(i)), hitonline);
842 
843  // calculate linear distance from start point and orthogonal distance from
844  // axis
845  double lindist = Get2DDistance((const PxPoint*)(&hitonline), (const PxPoint*)(&startHit));
846  double ortdist =
847  Get2DDistance((const PxPoint*)(&hitlist.at(i)), (const PxPoint*)(&hitonline));
848 
849  if (lindist < linearlimit && ortdist < ortlimit) { hitlistlocal_index.push_back(i); }
850  }
851  }
852 
855 
856  void GeometryUtilities::SelectPolygonHitList(std::vector<PxHit> const& hitlist,
857  std::vector<PxHit const*>& hitlistlocal) const
858  {
859  if (empty(hitlist)) { throw UtilException("Provided empty hit list!"); }
860 
861  hitlistlocal.clear();
862  unsigned char plane = hitlist.front().plane;
863 
864  // Define subset of hits to define polygon
865  std::map<double, const PxHit*> hitmap;
866  double qtotal = 0;
867  for (auto const& h : hitlist) {
868  hitmap.try_emplace(h.charge, &h);
869  qtotal += h.charge;
870  }
871  double qintegral = 0;
872  std::vector<const PxHit*> ordered_hits;
873  ordered_hits.reserve(hitlist.size());
874  for (auto hiter = hitmap.rbegin(); qintegral < qtotal * 0.95 && hiter != hitmap.rend();
875  ++hiter) {
876 
877  qintegral += (*hiter).first;
878  ordered_hits.push_back((*hiter).second);
879  }
880 
881  // Define container to hold found polygon corner PxHit index & distance
882  std::vector<size_t> hit_index(8, 0);
883  std::vector<double> hit_distance(8, 1e9);
884 
885  // Loop over hits and find corner points in the plane view
886  // Also fill corner edge points
887  std::vector<PxPoint> edges(4, PxPoint(plane, 0, 0));
888  double wire_max = fGeom.Nwires({0, 0, plane}) * fWiretoCm;
889  double time_max = fDetProp.NumberTimeSamples() * fTimetoCm;
890 
891  for (size_t index = 0; index < ordered_hits.size(); ++index) {
892 
893  if (ordered_hits.at(index)->t < 0 || ordered_hits.at(index)->w < 0 ||
894  ordered_hits.at(index)->t > time_max || ordered_hits.at(index)->w > wire_max) {
895 
896  throw UtilException(Form("Invalid wire/time (%g,%g) ... range is (0=>%g,0=>%g)",
897  ordered_hits.at(index)->w,
898  ordered_hits.at(index)->t,
899  wire_max,
900  time_max));
901  }
902 
903  double dist = 0;
904 
905  // Comparison w/ (Wire,0)
906  dist = ordered_hits.at(index)->t;
907  if (dist < hit_distance.at(1)) {
908  hit_distance.at(1) = dist;
909  hit_index.at(1) = index;
910  edges.at(0).t = ordered_hits.at(index)->t;
911  edges.at(1).t = ordered_hits.at(index)->t;
912  }
913 
914  // Comparison w/ (WireMax,Time)
915  dist = wire_max - ordered_hits.at(index)->w;
916  if (dist < hit_distance.at(3)) {
917  hit_distance.at(3) = dist;
918  hit_index.at(3) = index;
919  edges.at(1).w = ordered_hits.at(index)->w;
920  edges.at(2).w = ordered_hits.at(index)->w;
921  }
922 
923  // Comparison w/ (Wire,TimeMax)
924  dist = time_max - ordered_hits.at(index)->t;
925  if (dist < hit_distance.at(5)) {
926  hit_distance.at(5) = dist;
927  hit_index.at(5) = index;
928  edges.at(2).t = ordered_hits.at(index)->t;
929  edges.at(3).t = ordered_hits.at(index)->t;
930  }
931 
932  // Comparison w/ (0,Time)
933  dist = ordered_hits.at(index)->w;
934  if (dist < hit_distance.at(7)) {
935  hit_distance.at(7) = dist;
936  hit_index.at(7) = index;
937  edges.at(0).w = ordered_hits.at(index)->w;
938  edges.at(3).w = ordered_hits.at(index)->w;
939  }
940  }
941  for (size_t index = 0; index < ordered_hits.size(); ++index) {
942 
943  double dist = 0;
944  // Comparison w/ (0,0)
945  dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(0).t,
946  ordered_hits.at(index)->w - edges.at(0).w);
947  if (dist < hit_distance.at(0)) {
948  hit_distance.at(0) = dist;
949  hit_index.at(0) = index;
950  }
951 
952  // Comparison w/ (WireMax,0)
953  dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(1).t,
954  ordered_hits.at(index)->w - edges.at(1).w);
955  if (dist < hit_distance.at(2)) {
956  hit_distance.at(2) = dist;
957  hit_index.at(2) = index;
958  }
959 
960  // Comparison w/ (WireMax,TimeMax)
961  dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(2).t,
962  ordered_hits.at(index)->w - edges.at(2).w);
963  if (dist < hit_distance.at(4)) {
964  hit_distance.at(4) = dist;
965  hit_index.at(4) = index;
966  }
967 
968  // Comparison w/ (0,TimeMax)
969  dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(3).t,
970  ordered_hits.at(index)->w - edges.at(3).w);
971  if (dist < hit_distance.at(6)) {
972  hit_distance.at(6) = dist;
973  hit_index.at(6) = index;
974  }
975  }
976 
977  // Loop over the resulting hit indexes and append unique hits to define the
978  // polygon to the return hit list
979  std::set<size_t> unique_index;
980  std::vector<size_t> candidate_polygon;
981  candidate_polygon.reserve(9);
982  for (auto& index : hit_index) {
983  if (unique_index.find(index) == unique_index.end()) {
984  unique_index.insert(index);
985  candidate_polygon.push_back(index);
986  }
987  }
988  for (auto& index : hit_index) {
989  candidate_polygon.push_back(index);
990  break;
991  }
992 
993  if (unique_index.size() > 8) throw UtilException("Size of the polygon > 8!");
994 
995  // Untangle Polygon
996  candidate_polygon = PolyOverlap(ordered_hits, candidate_polygon);
997 
998  hitlistlocal.clear();
999  for (unsigned int i = 0; i < (candidate_polygon.size() - 1); i++) {
1000  hitlistlocal.push_back((const PxHit*)(ordered_hits.at(candidate_polygon.at(i))));
1001  }
1002  // check that polygon does not have more than 8 sides
1003  if (unique_index.size() > 8) throw UtilException("Size of the polygon > 8!");
1004  }
1005 
1006  std::vector<size_t> GeometryUtilities::PolyOverlap(std::vector<const PxHit*> ordered_hits,
1007  std::vector<size_t> candidate_polygon) const
1008  {
1009  // loop over edges
1010  for (unsigned int i = 0; i < (candidate_polygon.size() - 1); i++) {
1011  double Ax = ordered_hits.at(candidate_polygon.at(i))->w;
1012  double Ay = ordered_hits.at(candidate_polygon.at(i))->t;
1013  double Bx = ordered_hits.at(candidate_polygon.at(i + 1))->w;
1014  double By = ordered_hits.at(candidate_polygon.at(i + 1))->t;
1015  // loop over edges that have not been checked yet...
1016  // only ones furhter down in polygon
1017  for (unsigned int j = i + 2; j < (candidate_polygon.size() - 1); j++) {
1018  // avoid consecutive segments:
1019  if (candidate_polygon.at(i) == candidate_polygon.at(j + 1))
1020  continue;
1021  else {
1022  double Cx = ordered_hits.at(candidate_polygon.at(j))->w;
1023  double Cy = ordered_hits.at(candidate_polygon.at(j))->t;
1024  double Dx = ordered_hits.at(candidate_polygon.at(j + 1))->w;
1025  double Dy = ordered_hits.at(candidate_polygon.at(j + 1))->t;
1026 
1027  if ((Clockwise(Ax, Ay, Cx, Cy, Dx, Dy) != Clockwise(Bx, By, Cx, Cy, Dx, Dy)) and
1028  (Clockwise(Ax, Ay, Bx, By, Cx, Cy) != Clockwise(Ax, Ay, Bx, By, Dx, Dy))) {
1029  size_t tmp = candidate_polygon.at(i + 1);
1030  candidate_polygon.at(i + 1) = candidate_polygon.at(j);
1031  candidate_polygon.at(j) = tmp;
1032  // check that last element is still first (to close circle...)
1033  candidate_polygon.at(candidate_polygon.size() - 1) = candidate_polygon.at(0);
1034  // swapped polygon...now do recursion to make sure
1035  return PolyOverlap(ordered_hits, candidate_polygon);
1036  } // if crossing
1037  }
1038  } // second loop
1039  } // first loop
1040  return candidate_polygon;
1041  }
1042 
1043  bool GeometryUtilities::Clockwise(double const Ax,
1044  double const Ay,
1045  double const Bx,
1046  double const By,
1047  double const Cx,
1048  double const Cy) const
1049  {
1050  return (Cy - Ay) * (Bx - Ax) > (By - Ay) * (Cx - Ax);
1051  }
1052 
1053  PxHit GeometryUtilities::FindClosestHit(std::vector<PxHit> const& hitlist,
1054  unsigned int const wirein,
1055  double const timein) const
1056  {
1057  return hitlist[FindClosestHitIndex(hitlist, wirein, timein)];
1058  }
1059 
1060  unsigned int GeometryUtilities::FindClosestHitIndex(std::vector<PxHit> const& hitlist,
1061  unsigned int const wirein,
1062  double const timein) const
1063  {
1064  double min_length_from_start = 99999;
1065  unsigned int ret_ind = 0;
1066 
1067  for (unsigned int ii = 0; ii < hitlist.size(); ii++) {
1068  PxHit const& hit = hitlist[ii];
1069  double const dist_mod = Get2DDistance(wirein, timein, hit.w, hit.t);
1070  if (dist_mod < min_length_from_start) {
1071  min_length_from_start = dist_mod;
1072  ret_ind = ii;
1073  }
1074  }
1075 
1076  return ret_ind;
1077  }
1078 
1079 } // namespace
Float_t x
Definition: compare.C:6
double Get2Dslope(double deltawire, double deltatime) const
double Get2DDistance(double wire1, double time1, double wire2, double time2) const
detinfo::DetectorClocksData const & fClocks
Namespace for general, non-LArSoft-specific utilities.
Definition: PIDAAlg.h:26
Class def header for exception classes used in GeometryUtilities.
std::vector< double > vertangle
double PitchInView(unsigned int plane, double phi, double theta) const
Unknown view.
Definition: geo_types.h:142
Float_t y
Definition: compare.C:6
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
Double_t z
Definition: plot.C:276
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
double Temperature() const
In kelvin.
constexpr auto abs(T v)
Returns the absolute value of the argument.
void GetDirectionCosines(double phi, double theta, double *dirs) const
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
Float_t tmp
Definition: plot.C:35
int GetPointOnLine(double slope, double intercept, double wire1, double time1, double &wireout, double &timeout) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
WireGeo const & Wire(WireID const &wireid) const
Returns the specified wire.
PxPoint Get2DPointProjection(double const *xyz, unsigned int plane) const
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:248
double Efield(unsigned int planegap=0) const
kV/cm
detinfo::DetectorPropertiesData const & fDetProp
3-dimensional objects, potentially hits, clusters, prongs, etc.
Definition: geo_types.h:141
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
Access the description of detector geometry.
int GetYZ(const PxPoint *p0, const PxPoint *p1, double *yz) const
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:166
void SelectLocalHitlistIndex(const std::vector< PxHit > &hitlist, std::vector< unsigned int > &hitlistlocal_index, PxPoint &startHit, double &linearlimit, double &ortlimit, double &lineslopetest) const
int GetProjectedPoint(const PxPoint *p0, const PxPoint *p1, PxPoint &pN) const
geo::Point_t toWorldCoords(LocalPoint_t const &local) const
Transform point from local plane frame to world frame.
Definition: PlaneGeo.h:1225
double t
Definition: PxUtils.h:10
double Get2DPitchDistance(double angle, double inwire, double wire) const
int Get3DaxisN(int iplane0, int iplane1, double omega0, double omega1, double &phi, double &theta) const
geo::GeometryCore const & fGeom
std::vector< size_t > PolyOverlap(std::vector< const PxHit * > ordered_hits, std::vector< size_t > candidate_polygon) const
double GetTimeTicks(double x, unsigned int plane) const
void SelectLocalHitlist(const std::vector< PxHit > &hitlist, std::vector< const PxHit * > &hitlistlocal, PxPoint &startHit, double &linearlimit, double &ortlimit, double &lineslopetest) const
int GetPointOnLineWSlopes(double slope, double intercept, double ort_intercept, double &wireout, double &timeout) const
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
Detector simulation of raw signals on wires.
double w
Definition: PxUtils.h:9
Encapsulate the geometry of a wire.
bool Clockwise(double Ax, double Ay, double Bx, double By, double Cx, double Cy) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
void SelectPolygonHitList(const std::vector< PxHit > &hitlist, std::vector< const PxHit * > &hitlistlocal) const
int GetXYZ(const PxPoint *p0, const PxPoint *p1, double *xyz) const
double CalculatePitch(unsigned int iplane0, double phi, double theta) const
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
Encapsulate the construction of a single detector plane.
Contains all timing reference information for the detector.
PxHit FindClosestHit(std::vector< PxHit > const &hitlist, unsigned int wirein, double timein) const
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double Get2DangleFrom3D(unsigned int plane, double phi, double theta) const
double CalculatePitchPolar(unsigned int iplane0, double phi, double theta) const
GeometryUtilities(geo::GeometryCore const &geom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &propData)
double Get3DSpecialCaseTheta(int iplane0, int iplane1, double dw0, double dw1) const
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
int trigger_offset(DetectorClocksData const &data)
unsigned int FindClosestHitIndex(std::vector< PxHit > const &hitlist, unsigned int wirein, double timein) const
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
constexpr double kINVALID_DOUBLE
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
double Get2DPitchDistanceWSlope(double slope, double inwire, double wire) const
Float_t e
Definition: plot.C:35
double Get2Dangle(double deltawire, double deltatime) const
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
PxPoint Get2DPointProjectionCM(std::vector< double > const &xyz, unsigned int plane) const
geo::Point3DBase_t< PlaneGeoCoordinatesTag > LocalPoint_t
Type of points in the local GDML wire plane frame.
Definition: PlaneGeo.h:106
unsigned int plane
Definition: PxUtils.h:11
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:229
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:378