util::GeometryUtilities Class Reference

#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:

• lardata/lardata/Utilities/GeometryUtilities.h
• lardata/lardata/Utilities/GeometryUtilities.cxx