19 #include "cetlib/pow.h" 21 #include "TLorentzVector.h" 34 for (
unsigned int ip = 0; ip <
fNPlanes; ip++) {
63 double ln(0), mn(0), nn(0);
64 double phis(0), thetan(0);
69 unsigned int Cplane = 0, Iplane = 1;
73 double angleC, angleI, slopeC, slopeI, omegaC, omegaI;
79 if (omega0 == 0 || omega1 == 0) {
88 double alt_backwards = 0;
90 if (std::fabs(omega0) > (TMath::Pi() / 2.0) || std::fabs(omega1) > (TMath::Pi() / 2.0)) {
128 slopeC = std::tan(omegaC);
129 slopeI = std::tan(omegaI);
137 if (omegaC < TMath::Pi() && omegaC > 0)
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));
147 nn = (ln / (2 * std::cos(angleC))) *
148 ((1 / slopeC) + (1 / slopeI) + nfact * ((1 / slopeC) - (1 / slopeI)));
151 if (std::fabs(omegaC) > 0.01)
154 phis = std::asin(ln / TMath::Sqrt(ln * ln + nn * nn));
156 if (std::fabs(slopeC + slopeI) < 0.001)
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))
162 phis = (std::fabs(omegaC) > TMath::Pi() / 2) ? TMath::Pi() : 0;
164 if (nn < 0 && phis > 0 &&
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;
172 phi = phis * 180 / TMath::Pi();
180 thetan = -std::asin(mn / std::hypot(ln, mn, nn));
181 theta = thetan * 180 / TMath::Pi();
196 double splane, lplane;
199 if (dw0 == 0 && dw1 == 0)
return -1;
218 double tantheta = top / bottom;
230 auto const& plane =
fGeom.
Plane({0, 0, iplane});
232 mf::LogError(Form(
"Warning : no Pitch foreseen for view %d", plane.View()));
236 double const pi = TMath::Pi();
237 double const fTheta = pi / 2 - theta;
238 double const fPhi = -(phi + pi / 2);
240 double wirePitch = plane.WirePitch();
242 double cosgamma =
std::abs(std::sin(angleToVert) * std::cos(fTheta) +
243 std::cos(angleToVert) * std::sin(fTheta) * std::sin(fPhi));
245 return cosgamma > 0 ? wirePitch / cosgamma : -1.;
254 auto const& plane =
fGeom.
Plane({0, 0, iplane});
256 mf::LogError(Form(
"Warning : no Pitch foreseen for view %d", plane.View()));
260 double const fTheta = theta;
261 double const fPhi = phi;
263 double wirePitch = plane.WirePitch();
265 double cosgamma =
std::abs(std::sin(angleToVert) * std::cos(fTheta) +
266 std::cos(angleToVert) * std::sin(fTheta) * std::sin(fPhi));
268 return cosgamma > 0 ? wirePitch / cosgamma : -1.;
278 double timestart)
const 291 return Get2Dslope(endpoint->
w - startpoint->
w, endpoint->
t - startpoint->
t);
311 double timestart)
const 323 return Get2Dangle(endpoint->
w - startpoint->
w, endpoint->
t - startpoint->
t);
335 BC = ((double)dwire);
336 AC = ((double)dtime);
337 omega = std::asin(AC / std::hypot(AC, BC));
341 omega = TMath::Pi() -
std::abs(omega);
344 omega = -TMath::Pi() +
std::abs(omega);
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.));
377 TVector3
end = start + dir_vector;
382 start[2] * std::cos(alpha) - start[1] * std::sin(alpha)),
387 end[1] * std::sin(alpha)),
409 return std::hypot(point1->
w - point2->
w, point1->
t - point2->
t);
418 double radangle = TMath::Pi() * angle / 180;
419 if (std::cos(radangle) == 0)
return 9999;
439 double& timeout)
const 443 if (slope) { invslope = -1. / slope; }
445 double ort_intercept = time1 - invslope * (double)wire1;
447 if ((slope - invslope) != 0)
448 wireout = (ort_intercept - intercept) / (slope - invslope);
451 timeout = slope * wireout + intercept;
466 double intercept = startpoint->
t - slope * startpoint->
w;
482 if (slope) { invslope = -1. / slope; }
484 double ort_intercept = point1->
t - invslope * point1->
w;
486 if ((slope - invslope) != 0)
487 pointout.
w = (ort_intercept - intercept) / (slope - invslope);
489 pointout.
w = point1->
w;
491 pointout.
t = slope * pointout.
w + intercept;
506 double& timeout)
const 508 double intercept = timestart - slope * (double)wirestart;
510 return GetPointOnLine(slope, intercept, wire1, time1, wireout, timeout);
518 double ort_intercept,
520 double& timeout)
const 524 if (slope) { invslope = -1. / slope; }
528 wireout = (ort_intercept - intercept) / (slope - invslope);
529 timeout = slope * wireout + intercept;
543 double ort_intercept,
548 if (slope) invslope = -1. / slope;
550 pointonline.
w = (ort_intercept - intercept) / (slope - invslope);
551 pointonline.
t = slope * pointonline.
w + intercept;
559 for (
unsigned int i = 0; i <
fNPlanes; ++i) {
560 if (i == p0->
plane || i == p1->
plane)
continue;
576 pos.SetCoordinates(x, y, z);
597 std::cout <<
"\033[93mWarning\033[00m " 598 "\033[95m<<GeometryUtilities::GetYZ>>\033[00m" 600 <<
" 2D wire position " << p0->
w <<
" [cm] corresponds to negative wire number." 602 <<
" Forcing it to wire=0..." << std::endl
603 <<
"\033[93mWarning ends...\033[00m" << std::endl;
607 std::cout <<
"\033[93mWarning\033[00m " 608 "\033[95m<<GeometryUtilities::GetYZ>>\033[00m" 610 <<
" 2D wire position " << p0->
w <<
" [cm] exceeds max wire number " 612 <<
" Forcing it to the max wire number..." << std::endl
613 <<
"\033[93mWarning ends...\033[00m" << std::endl;
617 std::cout <<
"\033[93mWarning\033[00m " 618 "\033[95m<<GeometryUtilities::GetYZ>>\033[00m" 620 <<
" 2D wire position " << p1->
w <<
" [cm] corresponds to negative wire number." 622 <<
" Forcing it to wire=0..." << std::endl
623 <<
"\033[93mWarning ends...\033[00m" << std::endl;
627 std::cout <<
"\033[93mWarning\033[00m " 628 "\033[95m<<GeometryUtilities::GetYZ>>\033[00m" 630 <<
" 2D wire position " << p1->
w <<
" [cm] exceeds max wire number " 632 <<
" Forcing it to the max wire number..." << std::endl
633 <<
"\033[93mWarning ends...\033[00m" << std::endl;
695 unsigned int plane)
const 719 unsigned int plane)
const 721 double xyznew[3] = {(*xyz)[0], (*xyz)[1], (*xyz)[2]};
740 double dirs[3] = {0.};
745 double wirePitch = 0.;
746 double angleToVert = 0.;
748 auto const& planegeom =
fGeom.
Plane({0, 0, plane});
755 double cosgamma =
std::abs(std::sin(angleToVert) * dirs[1] + std::cos(angleToVert) * dirs[2]);
757 if (cosgamma < 1.
e-5)
760 std::cout <<
" returning 100" << std::endl;
764 return wirePitch / cosgamma;
770 theta *= (TMath::Pi() / 180);
771 phi *= (TMath::Pi() / 180);
772 dirs[0] = std::cos(theta) * std::sin(phi);
773 dirs[1] = std::sin(theta);
774 dirs[2] = std::cos(theta) * std::cos(phi);
778 std::vector<const PxHit*>& hitlistlocal,
782 double& lineslopetest)
const 786 hitlist, hitlistlocal, startHit, linearlimit, ortlimit, lineslopetest, testHit);
794 std::vector<const PxHit*>& hitlistlocal,
798 double& lineslopetest,
799 PxHit& averageHit)
const 802 hitlistlocal.clear();
803 std::vector<unsigned int> hitlistlocal_index;
805 hitlistlocal_index.clear();
808 hitlist, hitlistlocal_index, startHit, linearlimit, ortlimit, lineslopetest);
812 for (
size_t i = 0; i < hitlistlocal_index.size(); ++i) {
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;
820 if (hitlistlocal.size()) {
821 averageHit.
w = wiresum / (double)hitlistlocal.size();
822 averageHit.
t = timesum / ((double)hitlistlocal.size());
827 std::vector<unsigned int>& hitlistlocal_index,
831 double& lineslopetest)
const 834 hitlistlocal_index.clear();
835 double locintercept = startHit.
t - startHit.
w * lineslopetest;
837 for (
size_t i = 0; i < hitlist.size(); ++i) {
849 if (lindist < linearlimit && ortdist < ortlimit) { hitlistlocal_index.push_back(i); }
857 std::vector<PxHit const*>& hitlistlocal)
const 861 hitlistlocal.clear();
862 unsigned char plane = hitlist.front().plane;
865 std::map<double, const PxHit*> hitmap;
867 for (
auto const& h : hitlist) {
868 hitmap.try_emplace(h.charge, &h);
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();
877 qintegral += (*hiter).first;
878 ordered_hits.push_back((*hiter).second);
882 std::vector<size_t> hit_index(8, 0);
883 std::vector<double> hit_distance(8, 1e9);
887 std::vector<PxPoint> edges(4,
PxPoint(plane, 0, 0));
891 for (
size_t index = 0; index < ordered_hits.size(); ++index) {
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) {
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,
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;
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;
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;
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;
941 for (
size_t index = 0; index < ordered_hits.size(); ++index) {
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;
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;
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;
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;
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);
988 for (
auto& index : hit_index) {
989 candidate_polygon.push_back(index);
993 if (unique_index.size() > 8)
throw UtilException(
"Size of the polygon > 8!");
996 candidate_polygon =
PolyOverlap(ordered_hits, candidate_polygon);
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))));
1003 if (unique_index.size() > 8)
throw UtilException(
"Size of the polygon > 8!");
1007 std::vector<size_t> candidate_polygon)
const 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;
1017 for (
unsigned int j = i + 2; j < (candidate_polygon.size() - 1); j++) {
1019 if (candidate_polygon.at(i) == candidate_polygon.at(j + 1))
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;
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;
1033 candidate_polygon.at(candidate_polygon.size() - 1) = candidate_polygon.at(0);
1035 return PolyOverlap(ordered_hits, candidate_polygon);
1040 return candidate_polygon;
1048 double const Cy)
const 1050 return (Cy - Ay) * (Bx - Ax) > (By - Ay) * (Cx - Ax);
1054 unsigned int const wirein,
1055 double const timein)
const 1061 unsigned int const wirein,
1062 double const timein)
const 1064 double min_length_from_start = 99999;
1065 unsigned int ret_ind = 0;
1067 for (
unsigned int ii = 0; ii < hitlist.size(); ii++) {
1070 if (dist_mod < min_length_from_start) {
1071 min_length_from_start = dist_mod;
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.
Class def header for exception classes used in GeometryUtilities.
std::vector< double > vertangle
double PitchInView(unsigned int plane, double phi, double theta) const
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
The data type to uniquely identify a Plane.
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.
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.
double Efield(unsigned int planegap=0) const
kV/cm
detinfo::DetectorPropertiesData const & fDetProp
3-dimensional objects, potentially hits, clusters, prongs, etc.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
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.
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.
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
unsigned int NumberTimeSamples() 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.
Description of geometry of one entire detector.
Detector simulation of raw signals on wires.
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.
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.
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
double Get2Dangle(double deltawire, double deltatime) const
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
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.
constexpr Point origin()
Returns a origin position with a point of the specified type.
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
double WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.