#include <GeometryUtilities.h>

Public Member Functions
	GeometryUtilities (geo::GeometryCore const &geom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &propData)

Int_t	Get3DaxisN (Int_t iplane0, Int_t iplane1, Double_t omega0, Double_t omega1, Double_t &phi, Double_t &theta) const

Double_t	CalculatePitch (UInt_t iplane0, Double_t phi, Double_t theta) const

Double_t	CalculatePitchPolar (UInt_t iplane0, Double_t phi, Double_t theta) const

Double_t	Get3DSpecialCaseTheta (Int_t iplane0, Int_t iplane1, Double_t dw0, Double_t dw1) const

Double_t	Get2Dangle (Double_t deltawire, Double_t deltatime) const

Double_t	Get2Dangle (Double_t wireend, Double_t wirestart, Double_t timeend, Double_t timestart) const

double	Get2Dangle (const util::PxPoint endpoint, const util::PxPoint startpoint) const

double	Get2DangleFrom3D (unsigned int plane, double phi, double theta) const

double	Get2DangleFrom3D (unsigned int plane, TVector3 dir_vector) const

Double_t	Get2Dslope (Double_t deltawire, Double_t deltatime) const

Double_t	Get2Dslope (Double_t wireend, Double_t wirestart, Double_t timeend, Double_t timestart) const

double	Get2Dslope (const util::PxPoint endpoint, const util::PxPoint startpoint) const

Double_t	Get2DDistance (Double_t wire1, Double_t time1, Double_t wire2, Double_t time2) const

double	Get2DDistance (const util::PxPoint point1, const util::PxPoint point2) const

Double_t	Get2DPitchDistance (Double_t angle, Double_t inwire, Double_t wire) const

Double_t	Get2DPitchDistanceWSlope (Double_t slope, Double_t inwire, Double_t wire) const

Int_t	GetPointOnLine (Double_t slope, Double_t intercept, Double_t wire1, Double_t time1, Double_t &wireout, Double_t &timeout) const

Int_t	GetPointOnLine (Double_t slope, Double_t wirestart, Double_t timestart, Double_t wire1, Double_t time1, Double_t &wireout, Double_t &timeout) const

int	GetPointOnLine (Double_t slope, const util::PxPoint startpoint, const util::PxPoint point1, util::PxPoint &pointout) const

int	GetPointOnLine (double slope, double intercept, const util::PxPoint *point1, util::PxPoint &pointout) const

Int_t	GetPointOnLineWSlopes (Double_t slope, Double_t intercept, Double_t ort_intercept, Double_t &wireout, Double_t &timeout) const

Int_t	GetPointOnLineWSlopes (double slope, double intercept, double ort_intercept, util::PxPoint &pointonline) const

PxPoint	Get2DPointProjection (Double_t *xyz, Int_t plane) const

PxPoint	Get2DPointProjectionCM (std::vector< double > xyz, int plane) const

PxPoint	Get2DPointProjectionCM (double *xyz, int plane) const

PxPoint	Get2DPointProjectionCM (TLorentzVector *xyz, int plane) const

Double_t	GetTimeTicks (Double_t x, Int_t plane) const

Int_t	GetProjectedPoint (const PxPoint p0, const PxPoint p1, PxPoint &pN) const

util::PxHit	FindClosestHit (std::vector< util::PxHit > const &hitlist, unsigned int wirein, double timein) const

unsigned int	FindClosestHitIndex (std::vector< util::PxHit > const &hitlist, unsigned int wirein, double timein) const

Int_t	GetYZ (const PxPoint p0, const PxPoint p1, Double_t *yz) const

Int_t	GetXYZ (const PxPoint p0, const PxPoint p1, Double_t *xyz) const

Double_t	PitchInView (UInt_t plane, Double_t phi, Double_t theta) const

void	GetDirectionCosines (Double_t phi, Double_t theta, Double_t *dirs) const

void	SelectLocalHitlist (const std::vector< util::PxHit > &hitlist, std::vector< const util::PxHit * > &hitlistlocal, util::PxPoint &startHit, Double_t &linearlimit, Double_t &ortlimit, Double_t &lineslopetest) const

void	SelectLocalHitlist (const std::vector< util::PxHit > &hitlist, std::vector< const util::PxHit * > &hitlistlocal, util::PxPoint &startHit, Double_t &linearlimit, Double_t &ortlimit, Double_t &lineslopetest, util::PxHit &averageHit) const

void	SelectLocalHitlistIndex (const std::vector< util::PxHit > &hitlist, std::vector< unsigned int > &hitlistlocal_index, util::PxPoint &startHit, Double_t &linearlimit, Double_t &ortlimit, Double_t &lineslopetest) const

void	SelectPolygonHitList (const std::vector< util::PxHit > &hitlist, std::vector< const util::PxHit * > &hitlistlocal) const

std::vector< size_t >	PolyOverlap (std::vector< const util::PxHit * > ordered_hits, std::vector< size_t > candidate_polygon) const

bool	Clockwise (double Ax, double Ay, double Bx, double By, double Cx, double Cy) const

Double_t	TimeToCm () const

Double_t	WireToCm () const

Double_t	WireTimeToCmCm () const

UInt_t	Nplanes () const

Private Attributes
geo::GeometryCore const &	fGeom

detinfo::DetectorClocksData const &	fClocks

detinfo::DetectorPropertiesData const &	fDetProp

std::vector< Double_t >	vertangle

Double_t	fWirePitch

Double_t	fTimeTick

Double_t	fDriftVelocity

UInt_t	fNPlanes

Double_t	fWiretoCm

Double_t	fTimetoCm

Double_t	fWireTimetoCmCm

Detailed Description

Definition at line 35 of file GeometryUtilities.h.

Constructor & Destructor Documentation

util::GeometryUtilities::GeometryUtilities	(	geo::GeometryCore const &	geom,
		detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	propData
	)

Definition at line 26 of file GeometryUtilities.cxx.

     : fGeom{geom}, fClocks{clockData}, fDetProp{propData}
   {
     fNPlanes = fGeom.Nplanes();
     vertangle.resize(fNPlanes);
     for (UInt_t ip = 0; ip < fNPlanes; ip++)
       vertangle[ip] = fGeom.Plane(ip).Wire(0).ThetaZ(false) - TMath::Pi() / 2.; // wire angle
 
     fWirePitch = fGeom.WirePitch();
     fTimeTick = sampling_rate(fClocks) / 1000.;
     fDriftVelocity = fDetProp.DriftVelocity(fDetProp.Efield(), fDetProp.Temperature());
 
     fWiretoCm = fWirePitch;
     fTimetoCm = fTimeTick * fDriftVelocity;
     fWireTimetoCmCm = fTimetoCm / fWirePitch;
   }

Member Function Documentation

Double_t util::GeometryUtilities::CalculatePitch	(	UInt_t	iplane0,
		Double_t	phi,
		Double_t	theta
	)		const

Definition at line 225 of file GeometryUtilities.cxx.

   {
 
     Double_t pitch = -1.;
 
     if (fGeom.Plane(iplane).View() == geo::kUnknown || fGeom.Plane(iplane).View() == geo::k3D) {
       mf::LogError(Form("Warning :  no Pitch foreseen for view %d", fGeom.Plane(iplane).View()));
       return pitch;
     }
     else {
 
       Double_t pi = TMath::Pi();
       Double_t fTheta = pi / 2 - theta;
       Double_t fPhi = -(phi + pi / 2);
 
       for (UInt_t i = 0; i < fGeom.Nplanes(); ++i) {
         if (i == iplane) {
           Double_t wirePitch = fGeom.WirePitch(i);
           Double_t angleToVert =
             0.5 * TMath::Pi() - fGeom.WireAngleToVertical(fGeom.Plane(i).View());
           Double_t cosgamma =
             TMath::Abs(TMath::Sin(angleToVert) * TMath::Cos(fTheta) +
                        TMath::Cos(angleToVert) * TMath::Sin(fTheta) * TMath::Sin(fPhi));
 
           if (cosgamma > 0) pitch = wirePitch / cosgamma;
         } // end if the correct view
       }   // end loop over planes
     }     // end if a reasonable view
 
     return pitch;
   }

Double_t util::GeometryUtilities::CalculatePitchPolar	(	UInt_t	iplane0,
		Double_t	phi,
		Double_t	theta
	)		const

Definition at line 262 of file GeometryUtilities.cxx.

   {
 
     Double_t pitch = -1.;
 
     if (fGeom.Plane(iplane).View() == geo::kUnknown || fGeom.Plane(iplane).View() == geo::k3D) {
       mf::LogError(Form("Warning :  no Pitch foreseen for view %d", fGeom.Plane(iplane).View()));
       return pitch;
     }
     else {
 
       Double_t fTheta = theta;
       Double_t fPhi = phi;
 
       for (UInt_t i = 0; i < fGeom.Nplanes(); ++i) {
         if (i == iplane) {
           Double_t wirePitch = fGeom.WirePitch(i);
           Double_t angleToVert =
             0.5 * TMath::Pi() - fGeom.WireAngleToVertical(fGeom.Plane(i).View());
           Double_t cosgamma =
             TMath::Abs(TMath::Sin(angleToVert) * TMath::Cos(fTheta) +
                        TMath::Cos(angleToVert) * TMath::Sin(fTheta) * TMath::Sin(fPhi));
 
           if (cosgamma > 0) pitch = wirePitch / cosgamma;
         } // end if the correct view
       }   // end loop over planes
     }     // end if a reasonable view
 
     return pitch;
   }

bool util::GeometryUtilities::Clockwise	(	double	Ax,
		double	Ay,
		double	Bx,
		double	By,
		double	Cx,
		double	Cy
	)		const

Definition at line 1116 of file GeometryUtilities.cxx.

   {
     return (Cy - Ay) * (Bx - Ax) > (By - Ay) * (Cx - Ax);
   }

util::PxHit util::GeometryUtilities::FindClosestHit	(	std::vector< util::PxHit > const &	hitlist,
		unsigned int	wirein,
		double	timein
	)		const

Definition at line 1127 of file GeometryUtilities.cxx.

   {
     return hitlist[FindClosestHitIndex(hitlist, wirein, timein)];
   }

unsigned int util::GeometryUtilities::FindClosestHitIndex	(	std::vector< util::PxHit > const &	hitlist,
		unsigned int	wirein,
		double	timein
	)		const

Definition at line 1135 of file GeometryUtilities.cxx.

   {
     double min_length_from_start = 99999;
     unsigned int ret_ind = 0;
 
     for (unsigned int ii = 0; ii < hitlist.size(); ii++) {
       util::PxHit const& hit = hitlist[ii];
       double const dist_mod = Get2DDistance(wirein, timein, hit.w, hit.t);
       if (dist_mod < min_length_from_start) {
         min_length_from_start = dist_mod;
         ret_ind = ii;
       }
     }
 
     return ret_ind;
   }

Double_t util::GeometryUtilities::Get2Dangle	(	Double_t	deltawire,
		Double_t	deltatime
	)		const

Definition at line 361 of file GeometryUtilities.cxx.

   {
     Double_t BC, AC;
     Double_t omega;
 
     BC = ((Double_t)dwire); // in cm
     AC = ((Double_t)dtime); // in cm
     omega = std::asin(AC / std::hypot(AC, BC));
     if (BC < 0) // for the time being. Will check if it works for AC<0
     {
       if (AC > 0) {
         omega = TMath::Pi() - std::abs(omega); //
       }
       else if (AC < 0) {
         omega = -TMath::Pi() + std::abs(omega);
       }
       else {
         omega = TMath::Pi();
       }
     }
     return omega;
   }

Double_t util::GeometryUtilities::Get2Dangle	(	Double_t	wireend,
		Double_t	wirestart,
		Double_t	timeend,
		Double_t	timestart
	)		const

Definition at line 335 of file GeometryUtilities.cxx.

   {
 
     return Get2Dangle((wireend - wirestart) * fWiretoCm, (timeend - timestart) * fTimetoCm);
   }

Double_t util::GeometryUtilities::Get2Dangle	(	const util::PxPoint *	endpoint,
		const util::PxPoint *	startpoint
	)		const

Definition at line 350 of file GeometryUtilities.cxx.

   {
     return Get2Dangle(endpoint->w - startpoint->w, endpoint->t - startpoint->t);
   }

double util::GeometryUtilities::Get2DangleFrom3D	(	unsigned int	plane,
		double	phi,
		double	theta
	)		const

Definition at line 388 of file GeometryUtilities.cxx.

   {
     TVector3 dummyvector(cos(theta * TMath::Pi() / 180.) * sin(phi * TMath::Pi() / 180.),
                          sin(theta * TMath::Pi() / 180.),
                          cos(theta * TMath::Pi() / 180.) * cos(phi * TMath::Pi() / 180.));
 
     return Get2DangleFrom3D(plane, dummyvector);
   }

double util::GeometryUtilities::Get2DangleFrom3D	(	unsigned int	plane,
		TVector3	dir_vector
	)		const

Definition at line 401 of file GeometryUtilities.cxx.

   {
     double alpha = 0.5 * TMath::Pi() - fGeom.WireAngleToVertical(fGeom.Plane(plane).View());
     // create dummy  xyz point in middle of detector and another one in unit
     // length. calculate correspoding points in wire-time space and use the
     // differnces between those to return 2D a angle
 
     TVector3 start(fGeom.DetHalfWidth(), 0., fGeom.DetLength() / 2.);
     TVector3 end = start + dir_vector;
 
     // the wire coordinate is already in cm. The time needs to be converted.
     util::PxPoint startp(
       plane,
       (fGeom.DetHalfHeight() * sin(fabs(alpha)) + start[2] * cos(alpha) - start[1] * sin(alpha)),
       start[0]);
 
     util::PxPoint endp(
       plane,
       (fGeom.DetHalfHeight() * sin(fabs(alpha)) + end[2] * cos(alpha) - end[1] * sin(alpha)),
       end[0]);
 
     double angle = Get2Dangle(&endp, &startp);
 
     return angle;
   }

Double_t util::GeometryUtilities::Get2DDistance	(	Double_t	wire1,
		Double_t	time1,
		Double_t	wire2,
		Double_t	time2
	)		const

Definition at line 432 of file GeometryUtilities.cxx.

   {
 
     return std::hypot((wire1 - wire2) * fWiretoCm, (time1 - time2) * fTimetoCm);
   }

double util::GeometryUtilities::Get2DDistance	(	const util::PxPoint *	point1,
		const util::PxPoint *	point2
	)		const

Definition at line 442 of file GeometryUtilities.cxx.

   {
     return std::hypot(point1->w - point2->w, point1->t - point2->t);
   }

Double_t util::GeometryUtilities::Get2DPitchDistance	(	Double_t	angle,
		Double_t	inwire,
		Double_t	wire
	)		const

Definition at line 452 of file GeometryUtilities.cxx.

   {
     Double_t radangle = TMath::Pi() * angle / 180;
     if (std::cos(radangle) == 0) return 9999;
     return std::abs((wire - inwire) / std::cos(radangle)) * fWiretoCm;
   }

Double_t util::GeometryUtilities::Get2DPitchDistanceWSlope	(	Double_t	slope,
		Double_t	inwire,
		Double_t	wire
	)		const

Definition at line 461 of file GeometryUtilities.cxx.

   {
 
     return std::abs(wire - inwire) * TMath::Sqrt(1 + slope * slope) * fWiretoCm;
   }

PxPoint util::GeometryUtilities::Get2DPointProjection	(	Double_t *	xyz,
		Int_t	plane
	)		const

Todo:: : this should use the cryostat and tpc as well in the NearestWire method

Definition at line 715 of file GeometryUtilities.cxx.

   {
 
     PxPoint pN(0, 0, 0);
     const double origin[3] = {0.};
     Double_t pos[3];
     fGeom.Plane(plane).LocalToWorld(origin, pos);
     Double_t drifttick = (xyz[0] / fDriftVelocity) * (1. / fTimeTick);
 
     pos[1] = xyz[1];
     pos[2] = xyz[2];
 
     ///\todo: this should use the cryostat and tpc as well in the NearestWire
     /// method
 
     pN.w = fGeom.NearestWire(pos, plane);
     pN.t = drifttick - (pos[0] / fDriftVelocity) * (1. / fTimeTick) + trigger_offset(fClocks);
     pN.plane = plane;
 
     return pN;
   }

PxPoint util::GeometryUtilities::Get2DPointProjectionCM	(	std::vector< double >	xyz,
		int	plane
	)		const

Todo:: : this should use the cryostat and tpc as well in the NearestWire method

Definition at line 744 of file GeometryUtilities.cxx.

   {
 
     PxPoint pN(0, 0, 0);
 
     Double_t pos[3]{0., xyz[1], xyz[2]};
 
     ///\todo: this should use the cryostat and tpc as well in the NearestWire
     /// method
 
     pN.w = fGeom.NearestWire(pos, plane) * fWiretoCm;
     pN.t = xyz[0];
     pN.plane = plane;
 
     return pN;
   }

PxPoint util::GeometryUtilities::Get2DPointProjectionCM	(	double *	xyz,
		int	plane
	)		const

Todo:: : this should use the cryostat and tpc as well in the NearestWire method

Definition at line 762 of file GeometryUtilities.cxx.

   {
 
     PxPoint pN(0, 0, 0);
 
     Double_t pos[3]{0., xyz[1], xyz[2]};
 
     ///\todo: this should use the cryostat and tpc as well in the NearestWire
     /// method
 
     pN.w = fGeom.NearestWire(pos, plane) * fWiretoCm;
     pN.t = xyz[0];
     pN.plane = plane;
 
     return pN;
   }

PxPoint util::GeometryUtilities::Get2DPointProjectionCM	(	TLorentzVector *	xyz,
		int	plane
	)		const

Definition at line 780 of file GeometryUtilities.cxx.

   {
     double xyznew[3] = {(*xyz)[0], (*xyz)[1], (*xyz)[2]};
 
     return Get2DPointProjectionCM(xyznew, plane);
   }

Double_t util::GeometryUtilities::Get2Dslope	(	Double_t	deltawire,
		Double_t	deltatime
	)		const

Definition at line 325 of file GeometryUtilities.cxx.

   {
     return tan(Get2Dangle(dwire, dtime)) / fWireTimetoCmCm;
   }

Double_t util::GeometryUtilities::Get2Dslope	(	Double_t	wireend,
		Double_t	wirestart,
		Double_t	timeend,
		Double_t	timestart
	)		const

Definition at line 298 of file GeometryUtilities.cxx.

   {
 
     return GeometryUtilities::Get2Dslope((wireend - wirestart) * fWiretoCm,
                                          (timeend - timestart) * fTimetoCm);
   }

double util::GeometryUtilities::Get2Dslope	(	const util::PxPoint *	endpoint,
		const util::PxPoint *	startpoint
	)		const

Definition at line 313 of file GeometryUtilities.cxx.

   {
     return Get2Dslope(endpoint->w - startpoint->w, endpoint->t - startpoint->t);
   }

Int_t util::GeometryUtilities::Get3DaxisN	(	Int_t	iplane0,
		Int_t	iplane1,
		Double_t	omega0,
		Double_t	omega1,
		Double_t &	phi,
		Double_t &	theta
	)		const

Definition at line 53 of file GeometryUtilities.cxx.

   {
     // y, z, x coordinates
     Double_t ln(0), mn(0), nn(0);
     Double_t phis(0), thetan(0);
 
     // Pretend collection and induction planes.
     // "Collection" is the plane with the vertical angle equal to zero.
     // If both are non-zero, collection is the one with the negative angle.
     UInt_t Cplane = 0, Iplane = 1;
 
     // angleC and angleI are the respective angles to vertical in C/I
     // planes and slopeC, slopeI are the tangents of the track.
     Double_t angleC, angleI, slopeC, slopeI, omegaC, omegaI;
     omegaC = kINVALID_DOUBLE;
     omegaI = kINVALID_DOUBLE;
 
     // Don't know how to reconstruct these yet, so exit with error.
     // In
     if (omega0 == 0 || omega1 == 0) {
       phi = 0;
       theta = -999;
       return -1;
     }
 
     //////insert check for existence of planes.
 
     // check if backwards going track
     Double_t alt_backwards = 0;
 
     if (fabs(omega0) > (TMath::Pi() / 2.0) || fabs(omega1) > (TMath::Pi() / 2.0)) {
       alt_backwards = 1;
     }
 
     if (vertangle[iplane0] == 0) {
       // first plane is at 0 degrees
       Cplane = iplane0;
       Iplane = iplane1;
       omegaC = omega0;
       omegaI = omega1;
     }
     else if (vertangle[iplane1] == 0) {
       // second plane is at 0 degrees
       Cplane = iplane1;
       Iplane = iplane0;
       omegaC = omega1;
       omegaI = omega0;
     }
     else if (vertangle[iplane0] != 0 && vertangle[iplane1] != 0) {
       // both planes are at non zero degree - find the one with deg<0
       if (vertangle[iplane1] < vertangle[iplane0]) {
         Cplane = iplane1;
         Iplane = iplane0;
         omegaC = omega1;
         omegaI = omega0;
       }
       else if (vertangle[iplane1] > vertangle[iplane0]) {
         Cplane = iplane0;
         Iplane = iplane1;
         omegaC = omega0;
         omegaI = omega1;
       }
       else {
         // throw error - same plane.
         return -1;
       }
     }
 
     slopeC = tan(omegaC);
     slopeI = tan(omegaI);
     angleC = vertangle[Cplane];
     angleI = vertangle[Iplane];
 
     // 0 -1 factor depending on if one of the planes is vertical.
     bool nfact = !(vertangle[Cplane]);
 
     // ln represents y, omega is 2d angle -- in first 2 quadrants y is positive.
     if (omegaC < TMath::Pi() && omegaC > 0)
       ln = 1;
     else
       ln = -1;
 
     // calculate x and z using y ( ln )
     mn = (ln / (2 * sin(angleI))) * ((cos(angleI) / (slopeC * cos(angleC))) - (1 / slopeI) +
                                      nfact * (cos(angleI) / (cos(angleC) * slopeC) - 1 / slopeI));
 
     nn = (ln / (2 * cos(angleC))) *
          ((1 / slopeC) + (1 / slopeI) + nfact * ((1 / slopeC) - (1 / slopeI)));
 
     // Direction angles
     if (fabs(omegaC) > 0.01) // catch numeric error values
     {
       // phi=atan(ln/nn);
       phis = asin(ln / TMath::Sqrt(ln * ln + nn * nn));
 
       if (fabs(slopeC + slopeI) < 0.001)
         phis = 0;
       else if (fabs(omegaI) > 0.01 && (omegaI / fabs(omegaI) == -omegaC / fabs(omegaC)) &&
                (fabs(omegaC) < 1 * TMath::Pi() / 180 ||
                 fabs(omegaC) > 179 * TMath::Pi() / 180))           // angles have
         phis = (fabs(omegaC) > TMath::Pi() / 2) ? TMath::Pi() : 0; // angles are
 
       if (nn < 0 && phis > 0 &&
           !(!alt_backwards &&
             fabs(phis) <
               TMath::Pi() / 4)) // do not go back if track looks forward and phi is forward
         phis = (TMath::Pi()) - phis;
       else if (nn < 0 && phis < 0 && !(!alt_backwards && fabs(phis) < TMath::Pi() / 4))
         phis = (-TMath::Pi()) - phis;
 
       phi = phis * 180 / TMath::Pi();
     }
     // If plane2 (collection), phi = 2d angle (omegaC in this case)
     else {
       phis = omegaC;
       phi = omegaC;
     }
 
     thetan = -asin(mn / std::hypot(ln, mn, nn));
     theta = thetan * 180 / TMath::Pi();
 
     return 0;
   }

Double_t util::GeometryUtilities::Get3DSpecialCaseTheta	(	Int_t	iplane0,
		Int_t	iplane1,
		Double_t	dw0,
		Double_t	dw1
	)		const

Definition at line 187 of file GeometryUtilities.cxx.

   {
 
     Double_t splane, lplane; // plane in which the track is shorter and longer.
     Double_t sdw, ldw;
 
     if (dw0 == 0 && dw1 == 0) return -1;
 
     if (dw0 >= dw1) {
       lplane = iplane0;
       splane = iplane1;
       ldw = dw0;
       sdw = dw1;
     }
     else {
       lplane = iplane1;
       splane = iplane0;
       ldw = dw1;
       sdw = dw0;
     }
 
     Double_t top = (std::cos(vertangle[splane]) - std::cos(vertangle[lplane]) * sdw / ldw);
     Double_t bottom = tan(vertangle[lplane] * std::cos(vertangle[splane]));
     bottom -= tan(vertangle[splane] * std::cos(vertangle[lplane])) * sdw / ldw;
 
     Double_t tantheta = top / bottom;
 
     return atan(tantheta) * vertangle[lplane] / std::abs(vertangle[lplane]) * 180. / TMath::Pi();
   }

void util::GeometryUtilities::GetDirectionCosines	(	Double_t	phi,
		Double_t	theta,
		Double_t *	dirs
	)		const

Definition at line 834 of file GeometryUtilities.cxx.

   {
     theta *= (TMath::Pi() / 180);
     phi *= (TMath::Pi() / 180); // working on copies, it's ok.
     dirs[0] = TMath::Cos(theta) * TMath::Sin(phi);
     dirs[1] = TMath::Sin(theta);
     dirs[2] = TMath::Cos(theta) * TMath::Cos(phi);
   }

Int_t util::GeometryUtilities::GetPointOnLine	(	Double_t	slope,
		Double_t	intercept,
		Double_t	wire1,
		Double_t	time1,
		Double_t &	wireout,
		Double_t &	timeout
	)		const

Definition at line 472 of file GeometryUtilities.cxx.

   {
     Double_t invslope = 0;
 
     if (slope) { invslope = -1. / slope; }
 
     Double_t ort_intercept = time1 - invslope * (Double_t)wire1;
 
     if ((slope - invslope) != 0)
       wireout = (ort_intercept - intercept) / (slope - invslope);
     else
       wireout = wire1;
     timeout = slope * wireout + intercept;
 
     return 0;
   }

Int_t util::GeometryUtilities::GetPointOnLine	(	Double_t	slope,
		Double_t	wirestart,
		Double_t	timestart,
		Double_t	wire1,
		Double_t	time1,
		Double_t &	wireout,
		Double_t &	timeout
	)		const

Definition at line 541 of file GeometryUtilities.cxx.

   {
     Double_t intercept = timestart - slope * (Double_t)wirestart;
 
     return GetPointOnLine(slope, intercept, wire1, time1, wireout, timeout);
   }

int util::GeometryUtilities::GetPointOnLine	(	Double_t	slope,
		const util::PxPoint *	startpoint,
		const util::PxPoint *	point1,
		util::PxPoint &	pointout
	)		const

Definition at line 499 of file GeometryUtilities.cxx.

   {
 
     double intercept = startpoint->t - slope * startpoint->w;
 
     return GetPointOnLine(slope, intercept, point1, pointout);
   }

int util::GeometryUtilities::GetPointOnLine	(	double	slope,
		double	intercept,
		const util::PxPoint *	point1,
		util::PxPoint &	pointout
	)		const

Definition at line 515 of file GeometryUtilities.cxx.

   {
     double invslope = 0;
 
     if (slope) { invslope = -1. / slope; }
 
     double ort_intercept = point1->t - invslope * point1->w;
 
     if ((slope - invslope) != 0)
       pointout.w = (ort_intercept - intercept) / (slope - invslope);
     else
       pointout.w = point1->w;
 
     pointout.t = slope * pointout.w + intercept;
 
     return 0;
   }

Int_t util::GeometryUtilities::GetPointOnLineWSlopes	(	Double_t	slope,
		Double_t	intercept,
		Double_t	ort_intercept,
		Double_t &	wireout,
		Double_t &	timeout
	)		const

Definition at line 559 of file GeometryUtilities.cxx.

   {
     Double_t invslope = 0;
 
     if (slope) { invslope = -1. / slope; }
 
     invslope *= fWireTimetoCmCm * fWireTimetoCmCm;
 
     wireout = (ort_intercept - intercept) / (slope - invslope);
     timeout = slope * wireout + intercept;
 
     wireout /= fWiretoCm;
     timeout /= fTimetoCm;
 
     return 0;
   }

Int_t util::GeometryUtilities::GetPointOnLineWSlopes	(	double	slope,
		double	intercept,
		double	ort_intercept,
		util::PxPoint &	pointonline
	)		const

Definition at line 585 of file GeometryUtilities.cxx.

   {
     Double_t invslope = 0;
 
     if (slope) invslope = -1. / slope;
 
     pointonline.w = (ort_intercept - intercept) / (slope - invslope);
     pointonline.t = slope * pointonline.w + intercept;
     return 0;
   }

Int_t util::GeometryUtilities::GetProjectedPoint	(	const PxPoint *	p0,
		const PxPoint *	p1,
		PxPoint &	pN
	)		const

Definition at line 601 of file GeometryUtilities.cxx.

   {
 
     // determine third plane number
     for (UInt_t i = 0; i < fNPlanes; ++i) {
       if (i == p0->plane || i == p1->plane) continue;
       pN.plane = i;
     }
 
     // Assuming there is no problem ( and we found the best pair that comes
     // close in time ) we try to get the Y and Z coordinates for the start of
     // the shower.
     UInt_t chan1 = fGeom.PlaneWireToChannel(p0->plane, p0->w);
     UInt_t chan2 = fGeom.PlaneWireToChannel(p1->plane, p1->w);
     const double origin[3] = {0.};
     Double_t pos[3] = {0.};
     fGeom.Plane(p0->plane).LocalToWorld(origin, pos);
     Double_t x = (p0->t - trigger_offset(fClocks)) * fTimetoCm + pos[0];
 
     Double_t y, z;
     if (!fGeom.ChannelsIntersect(chan1, chan2, y, z)) return -1;
 
     pos[0] = x;
     pos[1] = y;
     pos[2] = z;
 
     pN = Get2DPointProjection(pos, pN.plane);
 
     return 0;
   }

Double_t util::GeometryUtilities::GetTimeTicks	(	Double_t	x,
		Int_t	plane
	)		const

Definition at line 788 of file GeometryUtilities.cxx.

   {
 
     const double origin[3] = {0.};
     Double_t pos[3];
     fGeom.Plane(plane).LocalToWorld(origin, pos);
     Double_t drifttick = (x / fDriftVelocity) * (1. / fTimeTick);
 
     return drifttick - (pos[0] / fDriftVelocity) * (1. / fTimeTick) + trigger_offset(fClocks);
   }

Int_t util::GeometryUtilities::GetXYZ	(	const PxPoint *	p0,
		const PxPoint *	p1,
		Double_t *	xyz
	)		const

Definition at line 695 of file GeometryUtilities.cxx.

   {
     const double origin[3] = {0.};
     Double_t pos[3] = {0.};
     fGeom.Plane(p0->plane).LocalToWorld(origin, pos);
     Double_t x = (p0->t) - trigger_offset(fClocks) * fTimetoCm + pos[0];
     double yz[2];
 
     GetYZ(p0, p1, yz);
 
     xyz[0] = x;
     xyz[1] = yz[0];
     xyz[2] = yz[1];
 
     return 0;
   }

Int_t util::GeometryUtilities::GetYZ	(	const PxPoint *	p0,
		const PxPoint *	p1,
		Double_t *	yz
	)		const

Definition at line 634 of file GeometryUtilities.cxx.

   {
     Double_t y, z;
 
     // Force to the closest wires if not in the range
     int z0 = p0->w / fWiretoCm;
     int z1 = p1->w / fWiretoCm;
     if (z0 < 0) {
       std::cout << "\033[93mWarning\033[00m "
                    "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
                 << std::endl
                 << " 2D wire position " << p0->w << " [cm] corresponds to negative wire number."
                 << std::endl
                 << " Forcing it to wire=0..." << std::endl
                 << "\033[93mWarning ends...\033[00m" << std::endl;
       z0 = 0;
     }
     else if (z0 >= (int)(fGeom.Nwires(p0->plane))) {
       std::cout << "\033[93mWarning\033[00m "
                    "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
                 << std::endl
                 << " 2D wire position " << p0->w << " [cm] exceeds max wire number "
                 << (fGeom.Nwires(p0->plane) - 1) << std::endl
                 << " Forcing it to the max wire number..." << std::endl
                 << "\033[93mWarning ends...\033[00m" << std::endl;
       z0 = fGeom.Nwires(p0->plane) - 1;
     }
     if (z1 < 0) {
       std::cout << "\033[93mWarning\033[00m "
                    "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
                 << std::endl
                 << " 2D wire position " << p1->w << " [cm] corresponds to negative wire number."
                 << std::endl
                 << " Forcing it to wire=0..." << std::endl
                 << "\033[93mWarning ends...\033[00m" << std::endl;
       z1 = 0;
     }
     if (z1 >= (int)(fGeom.Nwires(p1->plane))) {
       std::cout << "\033[93mWarning\033[00m "
                    "\033[95m<<GeometryUtilities::GetYZ>>\033[00m"
                 << std::endl
                 << " 2D wire position " << p1->w << " [cm] exceeds max wire number "
                 << (fGeom.Nwires(p0->plane) - 1) << std::endl
                 << " Forcing it to the max wire number..." << std::endl
                 << "\033[93mWarning ends...\033[00m" << std::endl;
       z1 = fGeom.Nwires(p1->plane) - 1;
     }
 
     UInt_t chan1 = fGeom.PlaneWireToChannel(p0->plane, z0);
     UInt_t chan2 = fGeom.PlaneWireToChannel(p1->plane, z1);
 
     if (!fGeom.ChannelsIntersect(chan1, chan2, y, z)) return -1;
 
     yz[0] = y;
     yz[1] = z;
 
     return 0;
   }

UInt_t util::GeometryUtilities::Nplanes ( ) const

inline

Definition at line 195 of file GeometryUtilities.h.

     {
       return fNPlanes;
     }

Double_t util::GeometryUtilities::PitchInView	(	UInt_t	plane,
		Double_t	phi,
		Double_t	theta
	)		const

Todo:: switch to using new Geometry::WireAngleToVertical(geo::View_t)

Todo:: and Geometry::WirePitch(geo::View_t) methods

Definition at line 803 of file GeometryUtilities.cxx.

   {
     Double_t dirs[3] = {0.};
     GetDirectionCosines(phi, theta, dirs);
 
     /// \todo switch to using new Geometry::WireAngleToVertical(geo::View_t)
     /// \todo and Geometry::WirePitch(geo::View_t) methods
     Double_t wirePitch = 0.;
     Double_t angleToVert = 0.;
 
     wirePitch = fGeom.WirePitch(plane);
     angleToVert = fGeom.WireAngleToVertical(fGeom.Plane(plane).View()) - 0.5 * TMath::Pi();
 
     //(sin(angleToVert),std::cos(angleToVert)) is the direction perpendicular to
     // wire fDir.front() is the direction of the track at the beginning of its
     // trajectory
     Double_t cosgamma =
       TMath::Abs(TMath::Sin(angleToVert) * dirs[1] + TMath::Cos(angleToVert) * dirs[2]);
 
     if (cosgamma < 1.e-5)
     // throw UtilException("cosgamma is basically 0, that can't be right");
     {
       std::cout << " returning 100" << std::endl;
       return 100;
     }
 
     return wirePitch / cosgamma;
   }

std::vector< size_t > util::GeometryUtilities::PolyOverlap	(	std::vector< const util::PxHit * >	ordered_hits,
		std::vector< size_t >	candidate_polygon
	)		const

Definition at line 1078 of file GeometryUtilities.cxx.

   {
     // loop over edges
     for (unsigned int i = 0; i < (candidate_polygon.size() - 1); i++) {
       double Ax = ordered_hits.at(candidate_polygon.at(i))->w;
       double Ay = ordered_hits.at(candidate_polygon.at(i))->t;
       double Bx = ordered_hits.at(candidate_polygon.at(i + 1))->w;
       double By = ordered_hits.at(candidate_polygon.at(i + 1))->t;
       // loop over edges that have not been checked yet...
       // only ones furhter down in polygon
       for (unsigned int j = i + 2; j < (candidate_polygon.size() - 1); j++) {
         // avoid consecutive segments:
         if (candidate_polygon.at(i) == candidate_polygon.at(j + 1))
           continue;
         else {
           double Cx = ordered_hits.at(candidate_polygon.at(j))->w;
           double Cy = ordered_hits.at(candidate_polygon.at(j))->t;
           double Dx = ordered_hits.at(candidate_polygon.at(j + 1))->w;
           double Dy = ordered_hits.at(candidate_polygon.at(j + 1))->t;
 
           if ((Clockwise(Ax, Ay, Cx, Cy, Dx, Dy) != Clockwise(Bx, By, Cx, Cy, Dx, Dy)) and
               (Clockwise(Ax, Ay, Bx, By, Cx, Cy) != Clockwise(Ax, Ay, Bx, By, Dx, Dy))) {
             size_t tmp = candidate_polygon.at(i + 1);
             candidate_polygon.at(i + 1) = candidate_polygon.at(j);
             candidate_polygon.at(j) = tmp;
             // check that last element is still first (to close circle...)
             candidate_polygon.at(candidate_polygon.size() - 1) = candidate_polygon.at(0);
             // swapped polygon...now do recursion to make sure
             return PolyOverlap(ordered_hits, candidate_polygon);
           } // if crossing
         }
       } // second loop
     }   // first loop
     return candidate_polygon;
   }

void util::GeometryUtilities::SelectLocalHitlist	(	const std::vector< util::PxHit > &	hitlist,
		std::vector< const util::PxHit * > &	hitlistlocal,
		util::PxPoint &	startHit,
		Double_t &	linearlimit,
		Double_t &	ortlimit,
		Double_t &	lineslopetest
	)		const

Definition at line 844 of file GeometryUtilities.cxx.

   {
     util::PxHit testHit;
     SelectLocalHitlist(
       hitlist, hitlistlocal, startHit, linearlimit, ortlimit, lineslopetest, testHit);
   }

void util::GeometryUtilities::SelectLocalHitlist	(	const std::vector< util::PxHit > &	hitlist,
		std::vector< const util::PxHit * > &	hitlistlocal,
		util::PxPoint &	startHit,
		Double_t &	linearlimit,
		Double_t &	ortlimit,
		Double_t &	lineslopetest,
		util::PxHit &	averageHit
	)		const

Definition at line 861 of file GeometryUtilities.cxx.

   {
 
     hitlistlocal.clear();
     std::vector<unsigned int> hitlistlocal_index;
 
     hitlistlocal_index.clear();
 
     SelectLocalHitlistIndex(
       hitlist, hitlistlocal_index, startHit, linearlimit, ortlimit, lineslopetest);
 
     double timesum = 0;
     double wiresum = 0;
     for (size_t i = 0; i < hitlistlocal_index.size(); ++i) {
 
       hitlistlocal.push_back((const util::PxHit*)(&(hitlist.at(hitlistlocal_index.at(i)))));
       timesum += hitlist.at(hitlistlocal_index.at(i)).t;
       wiresum += hitlist.at(hitlistlocal_index.at(i)).w;
     }
 
     averageHit.plane = startHit.plane;
     if (hitlistlocal.size()) {
       averageHit.w = wiresum / (double)hitlistlocal.size();
       averageHit.t = timesum / ((double)hitlistlocal.size());
     }
   }

void util::GeometryUtilities::SelectLocalHitlistIndex	(	const std::vector< util::PxHit > &	hitlist,
		std::vector< unsigned int > &	hitlistlocal_index,
		util::PxPoint &	startHit,
		Double_t &	linearlimit,
		Double_t &	ortlimit,
		Double_t &	lineslopetest
	)		const

Definition at line 895 of file GeometryUtilities.cxx.

   {
 
     hitlistlocal_index.clear();
     double locintercept = startHit.t - startHit.w * lineslopetest;
 
     for (size_t i = 0; i < hitlist.size(); ++i) {
 
       util::PxPoint hitonline;
 
       GetPointOnLine(lineslopetest, locintercept, (const util::PxHit*)(&hitlist.at(i)), hitonline);
 
       // calculate linear distance from start point and orthogonal distance from
       // axis
       Double_t lindist =
         Get2DDistance((const util::PxPoint*)(&hitonline), (const util::PxPoint*)(&startHit));
       Double_t ortdist =
         Get2DDistance((const util::PxPoint*)(&hitlist.at(i)), (const util::PxPoint*)(&hitonline));
 
       if (lindist < linearlimit && ortdist < ortlimit) { hitlistlocal_index.push_back(i); }
     }
   }

void util::GeometryUtilities::SelectPolygonHitList	(	const std::vector< util::PxHit > &	hitlist,
		std::vector< const util::PxHit * > &	hitlistlocal
	)		const

Definition at line 927 of file GeometryUtilities.cxx.

   {
     if (empty(hitlist)) { throw UtilException("Provided empty hit list!"); }
 
     hitlistlocal.clear();
     unsigned char plane = hitlist.front().plane;
 
     // Define subset of hits to define polygon
     std::map<double, const util::PxHit*> hitmap;
     double qtotal = 0;
     for (auto const& h : hitlist) {
       hitmap.try_emplace(h.charge, &h);
       qtotal += h.charge;
     }
     double qintegral = 0;
     std::vector<const util::PxHit*> ordered_hits;
     ordered_hits.reserve(hitlist.size());
     for (auto hiter = hitmap.rbegin(); qintegral < qtotal * 0.95 && hiter != hitmap.rend();
          ++hiter) {
 
       qintegral += (*hiter).first;
       ordered_hits.push_back((*hiter).second);
     }
 
     // Define container to hold found polygon corner PxHit index & distance
     std::vector<size_t> hit_index(8, 0);
     std::vector<double> hit_distance(8, 1e9);
 
     // Loop over hits and find corner points in the plane view
     // Also fill corner edge points
     std::vector<util::PxPoint> edges(4, PxPoint(plane, 0, 0));
     double wire_max = fGeom.Nwires(plane) * fWiretoCm;
     double time_max = fDetProp.NumberTimeSamples() * fTimetoCm;
 
     for (size_t index = 0; index < ordered_hits.size(); ++index) {
 
       if (ordered_hits.at(index)->t < 0 || ordered_hits.at(index)->w < 0 ||
           ordered_hits.at(index)->t > time_max || ordered_hits.at(index)->w > wire_max) {
 
         throw UtilException(Form("Invalid wire/time (%g,%g) ... range is (0=>%g,0=>%g)",
                                  ordered_hits.at(index)->w,
                                  ordered_hits.at(index)->t,
                                  wire_max,
                                  time_max));
       }
 
       double dist = 0;
 
       // Comparison w/ (Wire,0)
       dist = ordered_hits.at(index)->t;
       if (dist < hit_distance.at(1)) {
         hit_distance.at(1) = dist;
         hit_index.at(1) = index;
         edges.at(0).t = ordered_hits.at(index)->t;
         edges.at(1).t = ordered_hits.at(index)->t;
       }
 
       // Comparison w/ (WireMax,Time)
       dist = wire_max - ordered_hits.at(index)->w;
       if (dist < hit_distance.at(3)) {
         hit_distance.at(3) = dist;
         hit_index.at(3) = index;
         edges.at(1).w = ordered_hits.at(index)->w;
         edges.at(2).w = ordered_hits.at(index)->w;
       }
 
       // Comparison w/ (Wire,TimeMax)
       dist = time_max - ordered_hits.at(index)->t;
       if (dist < hit_distance.at(5)) {
         hit_distance.at(5) = dist;
         hit_index.at(5) = index;
         edges.at(2).t = ordered_hits.at(index)->t;
         edges.at(3).t = ordered_hits.at(index)->t;
       }
 
       // Comparison w/ (0,Time)
       dist = ordered_hits.at(index)->w;
       if (dist < hit_distance.at(7)) {
         hit_distance.at(7) = dist;
         hit_index.at(7) = index;
         edges.at(0).w = ordered_hits.at(index)->w;
         edges.at(3).w = ordered_hits.at(index)->w;
       }
     }
     for (size_t index = 0; index < ordered_hits.size(); ++index) {
 
       double dist = 0;
       // Comparison w/ (0,0)
       dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(0).t,
                                  ordered_hits.at(index)->w - edges.at(0).w);
       if (dist < hit_distance.at(0)) {
         hit_distance.at(0) = dist;
         hit_index.at(0) = index;
       }
 
       // Comparison w/ (WireMax,0)
       dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(1).t,
                                  ordered_hits.at(index)->w - edges.at(1).w);
       if (dist < hit_distance.at(2)) {
         hit_distance.at(2) = dist;
         hit_index.at(2) = index;
       }
 
       // Comparison w/ (WireMax,TimeMax)
       dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(2).t,
                                  ordered_hits.at(index)->w - edges.at(2).w);
       if (dist < hit_distance.at(4)) {
         hit_distance.at(4) = dist;
         hit_index.at(4) = index;
       }
 
       // Comparison w/ (0,TimeMax)
       dist = cet::sum_of_squares(ordered_hits.at(index)->t - edges.at(3).t,
                                  ordered_hits.at(index)->w - edges.at(3).w);
       if (dist < hit_distance.at(6)) {
         hit_distance.at(6) = dist;
         hit_index.at(6) = index;
       }
     }
 
     // Loop over the resulting hit indexes and append unique hits to define the
     // polygon to the return hit list
     std::set<size_t> unique_index;
     std::vector<size_t> candidate_polygon;
     candidate_polygon.reserve(9);
     for (auto& index : hit_index) {
       if (unique_index.find(index) == unique_index.end()) {
         unique_index.insert(index);
         candidate_polygon.push_back(index);
       }
     }
     for (auto& index : hit_index) {
       candidate_polygon.push_back(index);
       break;
     }
 
     if (unique_index.size() > 8) throw UtilException("Size of the polygon > 8!");
 
     // Untangle Polygon
     candidate_polygon = PolyOverlap(ordered_hits, candidate_polygon);
 
     hitlistlocal.clear();
     for (unsigned int i = 0; i < (candidate_polygon.size() - 1); i++) {
       hitlistlocal.push_back((const util::PxHit*)(ordered_hits.at(candidate_polygon.at(i))));
     }
     // check that polygon does not have more than 8 sides
     if (unique_index.size() > 8) throw UtilException("Size of the polygon > 8!");
   }

Double_t util::GeometryUtilities::TimeToCm ( ) const

inline

Definition at line 180 of file GeometryUtilities.h.

     {
       return fTimetoCm;
     }

Double_t util::GeometryUtilities::WireTimeToCmCm ( ) const

inline

Definition at line 190 of file GeometryUtilities.h.

     {
       return fWireTimetoCmCm;
     }

Double_t util::GeometryUtilities::WireToCm ( ) const

inline

Definition at line 185 of file GeometryUtilities.h.

     {
       return fWiretoCm;
     }

Member Data Documentation

detinfo::DetectorClocksData const& util::GeometryUtilities::fClocks

private

Definition at line 202 of file GeometryUtilities.h.

detinfo::DetectorPropertiesData const& util::GeometryUtilities::fDetProp

private

Definition at line 203 of file GeometryUtilities.h.

Double_t util::GeometryUtilities::fDriftVelocity

private

Definition at line 208 of file GeometryUtilities.h.

geo::GeometryCore const& util::GeometryUtilities::fGeom

private

Definition at line 201 of file GeometryUtilities.h.

UInt_t util::GeometryUtilities::fNPlanes

private

Definition at line 209 of file GeometryUtilities.h.

Double_t util::GeometryUtilities::fTimeTick

private

Definition at line 207 of file GeometryUtilities.h.

Double_t util::GeometryUtilities::fTimetoCm

private

Definition at line 211 of file GeometryUtilities.h.

Double_t util::GeometryUtilities::fWirePitch

private

Definition at line 206 of file GeometryUtilities.h.

Double_t util::GeometryUtilities::fWireTimetoCmCm

private

Definition at line 212 of file GeometryUtilities.h.

Double_t util::GeometryUtilities::fWiretoCm

private

Definition at line 210 of file GeometryUtilities.h.

std::vector<Double_t> util::GeometryUtilities::vertangle

private

Definition at line 205 of file GeometryUtilities.h.

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

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation