diff --git a/inc/TRestDetectorReadoutPlane.h b/inc/TRestDetectorReadoutPlane.h index bc5df878..c3a9e8d8 100644 --- a/inc/TRestDetectorReadoutPlane.h +++ b/inc/TRestDetectorReadoutPlane.h @@ -81,6 +81,8 @@ class TRestDetectorReadoutPlane { void SetRotation(Double_t radians); + void SetAxisX(const TVector3& axis); + // Getters /// Returns an integer with the plane id number. inline Int_t GetID() const { return fId; } @@ -154,6 +156,11 @@ class TRestDetectorReadoutPlane { Int_t GetModuleIDFromPosition(Double_t x, Double_t y, Double_t z); + TVector2 GetPositionInPlane(const TVector3& point) const; + Double_t GetDistanceToPlane(const TVector3& point) const; + + TVector3 GetPositionInWorld(const TVector2& point, Double_t height = 0) const; + void Draw(); void Print(Int_t DetailLevel = 0); diff --git a/src/TRestDetectorReadoutPlane.cxx b/src/TRestDetectorReadoutPlane.cxx index e86236f2..9f291022 100644 --- a/src/TRestDetectorReadoutPlane.cxx +++ b/src/TRestDetectorReadoutPlane.cxx @@ -449,25 +449,25 @@ void TRestDetectorReadoutPlane::GetBoundaries(double& xmin, double& xmax, double } void TRestDetectorReadoutPlane::UpdateAxes() { // idempotent - const TVector3 originalNormal = {0, 0, 1}; + const TVector3 zUnit = {0, 0, 1}; fAxisX = {1, 0, 0}; fAxisY = {0, 1, 0}; constexpr double tolerance = 1E-6; - // Check if fNormal is different from the original normal - if ((fNormal - originalNormal).Mag2() < tolerance) { + // Check if fNormal is different from (0,0,1) + if ((fNormal - zUnit).Mag2() < tolerance) { // do nothing - } else if ((fNormal + originalNormal).Mag2() < tolerance) { - // normal vector is opposite to the original normal (0,0,-1), we must also flip the axes + } else if ((fNormal + zUnit).Mag2() < tolerance) { + // normal vector is opposite to (0,0,1), we must also flip the axes fAxisX = {0, -1, 0}; fAxisY = {-1, 0, 0}; } else { // Calculate the rotation axis by taking the cross product between the original normal and fNormal - TVector3 rotationAxis = originalNormal.Cross(fNormal); + TVector3 rotationAxis = zUnit.Cross(fNormal); // Calculate the rotation angle using the dot product between the original normal and fNormal - double rotationAngle = acos(originalNormal.Dot(fNormal) / (originalNormal.Mag() * fNormal.Mag())); + double rotationAngle = acos(zUnit.Dot(fNormal) / (zUnit.Mag() * fNormal.Mag())); // Rotate the axes around the rotation axis by the rotation angle fAxisX.Rotate(rotationAngle, rotationAxis); @@ -504,6 +504,36 @@ void TRestDetectorReadoutPlane::UpdateAxes() { // idempotent } void TRestDetectorReadoutPlane::SetRotation(Double_t radians) { - fRotation = radians; + // sets fRotation modulo 2pi + fRotation = TVector2::Phi_0_2pi(radians); UpdateAxes(); } + +TVector2 TRestDetectorReadoutPlane::GetPositionInPlane(const TVector3& point) const { + // Given a point in space, returns the position of the point in the plane using the plane's axes + // The position is returned in the plane's local coordinates (in mm) + return {fAxisX.Dot(point - fPosition), + fAxisY.Dot(point - fPosition)}; // dot product between vectors is the projection +} + +TVector3 TRestDetectorReadoutPlane::GetPositionInWorld(const TVector2& point, Double_t height) const { + return fPosition + point.X() * fAxisX + point.Y() * fAxisY + height * fNormal; +} + +Double_t TRestDetectorReadoutPlane::GetDistanceToPlane(const TVector3& point) const { + return (point - fPosition).Dot(fNormal); +} + +void TRestDetectorReadoutPlane::SetAxisX(const TVector3& axis) { + const TVector3 axisInPlane = (axis - axis.Dot(fNormal) * fNormal).Unit(); + if (axisInPlane.Mag2() < 1E-6) { + RESTError << "TRestDetectorReadoutPlane::SetAxisX() : " + << "The X-axis vector must not be parallel to the normal vector." << RESTendl; + exit(1); + } + + // compute the angle between the new axis and the old axis + const Double_t angle = fAxisX.Angle(axisInPlane); + + SetRotation(fRotation - angle); +} diff --git a/test/src/TRestDetectorReadout.cxx b/test/src/TRestDetectorReadout.cxx index add1536f..9bc23846 100644 --- a/test/src/TRestDetectorReadout.cxx +++ b/test/src/TRestDetectorReadout.cxx @@ -54,3 +54,43 @@ TEST(TRestDetectorReadout, Axes) { EXPECT_TRUE(AreEqual(plane.GetAxisY(), {0, 0.99995, -0.0099995})); } } + +TEST(TRestDetectorReadout, SetAxisX) { + TRestDetectorReadoutPlane plane; + plane.SetNormal({0, 0.5, 1}); + plane.SetPosition({10.0, 5.0, 50.0}); + plane.SetRotation(TMath::DegToRad() * 32.0); + plane.SetHeight(10.0); + plane.PrintMetadata(); + // default values + + // axis will be projected into the plane + const TVector3 newXAxis = {1.0, -2.0, 5.0}; + + plane.SetAxisX(newXAxis); + + const TVector3 axisX = plane.GetAxisX(); + const TVector3 newXAxisProjected = + (newXAxis - plane.GetNormal() * newXAxis.Dot(plane.GetNormal())).Unit(); + cout << "axisX: " << axisX.X() << ", " << axisX.Y() << ", " << axisX.Z() << endl; + cout << "axisXProjected: " << newXAxisProjected.X() << ", " << newXAxisProjected.Y() << ", " + << newXAxisProjected.Z() << endl; + EXPECT_TRUE(AreEqual(plane.GetAxisX(), newXAxisProjected)); +} + +TEST(TRestDetectorReadout, Module) { + TRestDetectorReadoutPlane plane; + plane.SetNormal({0, 0, 1}); + plane.SetPosition({10.0, 5.0, 50.0}); + plane.SetRotation(TMath::DegToRad() * 90.0); + plane.SetHeight(10.0); + + plane.PrintMetadata(); + + const TVector3 worldPoint = {0, 0, 100.0}; + const TVector2 local = plane.GetPositionInPlane(worldPoint); + cout << "local: " << local.X() << ", " << local.Y() << endl; + + const auto distance = plane.GetDistanceToPlane(worldPoint); + cout << "distance: " << distance << endl; +}