diff --git a/inc/TRestDetectorExperimentalReadout.h b/inc/TRestDetectorExperimentalReadout.h new file mode 100644 index 00000000..37d81d44 --- /dev/null +++ b/inc/TRestDetectorExperimentalReadout.h @@ -0,0 +1,69 @@ +/** + * @file TRestDetectorExperimentalReadout.h + * @brief Defines the TRestDetectorExperimentalReadout class. + */ + +#ifndef REST_TRESTDETECTOREXPERIMENTALREADOUT_H +#define REST_TRESTDETECTOREXPERIMENTALREADOUT_H + +#include + +#include + +#include "TRestDetectorExperimentalReadoutModule.h" + +/** + * @class TRestDetectorExperimentalReadout + * @brief Represents an experimental readout detector. + * @inherits TNamed + * @author Luis Obis + */ +class TRestDetectorExperimentalReadout : public TNamed { + private: + std::vector fModules; /**< The modules of the detector. */ + + public: + /** + * @brief Constructs a TRestDetectorExperimentalReadout object. + */ + TRestDetectorExperimentalReadout() = default; + + /** + * @brief Destructor for TRestDetectorExperimentalReadout. + */ + ~TRestDetectorExperimentalReadout() override = default; + + /** + * @brief Returns the number of modules in the detector. + * @return The number of modules. + */ + size_t GetNumberOfModules() { return fModules.size(); } + + /** + * @brief Adds a module to the detector. + * @param module The module to add. + */ + void AddModule(const TRestDetectorExperimentalReadoutModule& module) { fModules.push_back(module); } + + /** + * @brief Returns the modules of the detector. + * @return The modules. + */ + std::vector GetModules() const; + + /** + * @brief Returns the modules that contain a given point. + * @param point The point to check. + * @return The modules containing the point. + */ + std::vector GetModulesForPoint( + const TVector3& point) const; + + std::vector GetPixelsForPoint(const TVector3& point) const; + + std::vector GetChannelsForPoint(const TVector3& point) const; + + ClassDef(TRestDetectorExperimentalReadout, 1); +}; + +#endif // REST_TRESTDETECTOREXPERIMENTALREADOUT_H diff --git a/inc/TRestDetectorExperimentalReadoutChannel.h b/inc/TRestDetectorExperimentalReadoutChannel.h new file mode 100644 index 00000000..c0388764 --- /dev/null +++ b/inc/TRestDetectorExperimentalReadoutChannel.h @@ -0,0 +1,30 @@ +/** + * @file TRestDetectorExperimentalReadoutModule.h + * @brief Defines the TRestDetectorExperimentalReadoutModule class. + */ + +#ifndef REST_TRESTDETECTOREXPERIMENTALREADOUTCHANNEL_H +#define REST_TRESTDETECTOREXPERIMENTALREADOUTCHANNEL_H + +#include + +#include "TRestDetectorExperimentalReadoutPixel.h" + +/** + * @class TRestDetectorExperimentalReadoutChannel + * @brief Represents a channel in an experimental readout detector. + * @author Luis Obis + */ +class TRestDetectorExperimentalReadoutChannel { + private: + std::vector fPixels; /**< The pixels of the channel. */ + unsigned short fChannel; /**< The channel of the pixel. */ + + public: + TRestDetectorExperimentalReadoutChannel() = default; + ~TRestDetectorExperimentalReadoutChannel() = default; + + ClassDef(TRestDetectorExperimentalReadoutChannel, 1); +}; + +#endif // REST_TRESTDETECTOREXPERIMENTALREADOUTCHANNEL_H diff --git a/inc/TRestDetectorExperimentalReadoutModule.h b/inc/TRestDetectorExperimentalReadoutModule.h new file mode 100644 index 00000000..0b25393e --- /dev/null +++ b/inc/TRestDetectorExperimentalReadoutModule.h @@ -0,0 +1,139 @@ +/** + * @file TRestDetectorExperimentalReadoutModule.h + * @brief Defines the TRestDetectorExperimentalReadoutModule class. + */ + +#ifndef REST_TRESTDETECTOREXPERIMENTALREADOUTMODULE_H +#define REST_TRESTDETECTOREXPERIMENTALREADOUTMODULE_H + +#include + +#include "TRestDetectorExperimentalReadoutPixel.h" + +class KDTree; + +/** + * @class TRestDetectorExperimentalReadoutModule + * @brief Represents a module in an experimental readout detector. + * @author Luis Obis + */ +class TRestDetectorExperimentalReadoutModule { + private: + std::vector fPixels; /**< The pixels of the module. */ + + TVector3 fPosition = {0, 0, 0}; /**< The position of the module in 3D space. */ + TVector3 fNormal = {0, 0, 1}; /**< The normal vector to the module. */ + double fHeight = 1.0; /**< The height of the module (drift distance). */ + std::pair fCoordinateAxes = {{1, 0, 0}, + {0, 1, 0}}; /**< The coordinate axes of the module. */ + double fRotation = + 0; /**< The rotation angle of the module with respect to the normal vector in degrees. */ + // auxiliary data structures + KDTree* fKDTree = nullptr; /**< The KDTree used for spatial queries. */ + double fSearchRadius = 0; /**< The search radius for spatial queries. */ + std::vector fConvexHull; /**< The convex hull of the module. */ + std::map> + fChannelToPixels; /**< A map from channel to pixels. */ + /** + * @brief Builds the KDTree for the module. + */ + void BuildKDTree(); + + /** + * @brief Computes the convex hull of the module. + * @return The convex hull of the module. + */ + std::vector ComputeConvexHull(); + + void UpdateAxes(); + + public: + // needed for root + TRestDetectorExperimentalReadoutModule() = default; + virtual ~TRestDetectorExperimentalReadoutModule() { + // delete fKDTree; + } + + void SetPosition(const TVector3& position) { fPosition = position; } + void SetNormal(const TVector3& normal); + void SetHeight(double height); + void SetRotation(double rotationInDegrees); + + TVector3 GetPosition() const { return fPosition; } + TVector3 GetNormal() const { return fNormal; } + double GetHeight() const { return fHeight; } + std::pair GetCoordinateAxes() const { return fCoordinateAxes; } + double GetRotation() const { return fRotation; } + /** + * @brief Returns the number of pixels in the module. + * @return The number of pixels. + */ + size_t GetNumberOfPixels() const { return fPixels.size(); } + + /** + * @brief Returns the pixels of the module. + * @return The pixels. + */ + std::vector GetPixels() const { return fPixels; } + + /** + * @brief Sets the pixels of the module. + * @param pixels The pixels to set. + */ + void SetPixels(const std::vector& pixels); + + /** + * @brief Returns the convex hull of the module. + * @return The convex hull. + */ + std::vector GetConvexHull() const { return fConvexHull; } + + /** + * @brief Returns the pixels that contain a given point. + * @param point The point to check. + * @return The pixels containing the point. + */ + std::vector GetPixelsForPoint(const TVector2& point) const; + std::vector GetPixelsForPoint(const TVector3& point) const; + std::vector GetPixelsForChannel( + unsigned short channel) const; + + /** + * @brief Returns the Z coordinate of a given point inside the module. + * @param point The point to check. + * @return The Z coordinate of the point. + */ + double GetZ(const TVector3& point) const; + + /** + * @brief Checks if a given point is inside the module in the module's local XY plane. + * @param point The point to check. + * @return True if the point is inside the module, false otherwise. + */ + bool IsInside(const TVector2& point) const; + + /** + * @brief Checks if a given point is inside the module in 3D space. + * @param point The point to check. + * @return True if the point is inside the module, false otherwise. + */ + bool IsInside(const TVector3& point) const; + + /** + * @brief Transforms a point from 3D space to the module's local XY plane coordinates. + * @param point The point to transform. + * @return The transformed point in the local XY plane. + */ + TVector2 TransformToModuleCoordinates(const TVector3& point) const; + + /** + * @brief Transforms a point from the module's local XY plane coordinates to 3D space. + * @param point The point to transform. + * @return The transformed point in 3D space. + */ + TVector3 TransformToAbsoluteCoordinates(const TVector2& point) const; + + ClassDef(TRestDetectorExperimentalReadoutModule, 1); +}; + +#endif // REST_TRESTDETECTOREXPERIMENTALREADOUTMODULE_H diff --git a/inc/TRestDetectorExperimentalReadoutPixel.h b/inc/TRestDetectorExperimentalReadoutPixel.h new file mode 100644 index 00000000..97dea0c5 --- /dev/null +++ b/inc/TRestDetectorExperimentalReadoutPixel.h @@ -0,0 +1,151 @@ +/** + * @file TRestDetectorExperimentalReadoutPixel.h + * @brief Defines the TRestDetectorExperimentalReadoutPixel class. + */ + +#ifndef REST_TRESTDETECTOREXPERIMENTALREADOUTPIXEL_H +#define REST_TRESTDETECTOREXPERIMENTALREADOUTPIXEL_H + +#include +#include + +#include +#include +#include +#include + +/** + + * @class TRestDetectorExperimentalReadoutPixel + * @brief Represents a pixel in an experimental readout detector. + * @author Luis Obis + +*/ +class TRestDetectorExperimentalReadoutPixel { + private: + std::vector fVertices; /**< The vertices of the pixel. */ + TVector2 fCenter; /**< The center of the pixel in mm. A pixel is guaranteed to be inside a circle with + center fCenter and radius fRadius. */ + double fRadius; /**< The radius of the pixel in mm. A pixel is guaranteed to be inside a circle with + center fCenter and radius fRadius. */ + + unsigned short fChannel; /**< The channel of the pixel. */ + /** + + @brief Initializes the vertices of the pixel. + @param vertices The vertices of the pixel. + */ + void InitializeVertices(const std::vector& vertices); + /** + + @brief Calculates the vertices of a rectangular pixel. + @param center The center of the pixel. + @param size The size of the pixel (width, height). + @return The vertices of the rectangular pixel. + */ + static std::vector GetRectangularVertices(const TVector2& center, + const std::pair& size); + + public: + /** + + @brief Constructs a TRestDetectorExperimentalReadoutPixel object with given vertices. + @param vertices The vertices of the pixel. + */ + explicit TRestDetectorExperimentalReadoutPixel(const std::vector& vertices, + unsigned short channel = 0); + /** + + @brief Constructs a rectangular TRestDetectorExperimentalReadoutPixel object. + @param center The center of the pixel. + @param size The size of the pixel (width, height). + */ + TRestDetectorExperimentalReadoutPixel(const TVector2& center, const std::pair& size, + unsigned short channel = 0); + /** + + @brief Constructs a square TRestDetectorExperimentalReadoutPixel object. + @param center The center of the pixel. + @param size The size of the pixel (width and height). + */ + TRestDetectorExperimentalReadoutPixel(const TVector2& center, double size, unsigned short channel = 0); + // needed for root + TRestDetectorExperimentalReadoutPixel() = default; + virtual ~TRestDetectorExperimentalReadoutPixel() = default; + + /** + + @brief Returns the center of the pixel. + @return The center of the pixel. + */ + TVector2 GetCenter() const { return fCenter; } + /** + + @brief Returns the radius of the pixel. + @return The radius of the pixel. + */ + double GetRadius() const { return fRadius; } + /** + + @brief Returns the vertices of the pixel. + @return The vertices of the pixel. + */ + std::vector GetVertices() const { return fVertices; } + /** + + @brief Checks if a given point is inside the pixel. + @param point The point to check. + @return True if the point is inside the pixel, false otherwise. + */ + bool IsInside(const TVector2& point) const; + + unsigned short GetChannel() const { return fChannel; } + + ClassDef(TRestDetectorExperimentalReadoutPixel, 1); +}; + +namespace readout { + +/** + +@brief Calculates the triple product of three vectors. +@param A The first vector. +@param B The second vector. +@param C The third vector. +@return The cross product of the vectors. +*/ +double crossProduct(const TVector2& A, const TVector2& B, const TVector2& C); +/** + +@brief Calculates the squared distance between two points. +@param A The first point. +@param B The second point. +@return The squared distance between the points. +*/ +double squaredDistance(const TVector2& A, const TVector2& B); +/** + +@brief Comparator function to compare two points based on their polar angle with respect to the anchor point. +@param A The first point. +@param B The second point. +@param anchor The anchor point. +@return True if point A is smaller than point B, false otherwise. +*/ +bool comparePoints(const TVector2& A, const TVector2& B, const TVector2& anchor); +/** + +@brief Finds the convex hull of a set of points using the Graham's scan algorithm. +@param _points The input points. +@return The convex hull of the points. +*/ +std::vector findConvexHull(const std::vector& _points); +/** + +@brief Checks if a given point is inside a convex hull. +@param point The point to check. +@param convexHull The convex hull. +@return True if the point is inside the convex hull, false otherwise. +*/ +bool IsPointInsideConvexHull(const TVector2& point, const std::vector& convexHull); +} // namespace readout +#endif // REST_TRESTDETECTOREXPERIMENTALREADOUTPIXEL_H diff --git a/src/TRestDetectorExperimentalReadout.cxx b/src/TRestDetectorExperimentalReadout.cxx new file mode 100644 index 00000000..02d817bc --- /dev/null +++ b/src/TRestDetectorExperimentalReadout.cxx @@ -0,0 +1,43 @@ +#include "TRestDetectorExperimentalReadout.h" + +using namespace std; + +ClassImp(TRestDetectorExperimentalReadout); + +vector TRestDetectorExperimentalReadout::GetModulesForPoint( + const TVector3& point) const { + vector modules; + for (const auto& module : fModules) { + if (module.IsInside(point)) { + modules.push_back(&module); + } + } + return modules; +} + +vector TRestDetectorExperimentalReadout::GetModules() const { + vector modules; + modules.reserve(fModules.size()); + for (const auto& module : fModules) { + modules.push_back(&module); + } + return modules; +} + +vector TRestDetectorExperimentalReadout::GetPixelsForPoint( + const TVector3& point) const { + vector pixels; + for (const auto& module : GetModulesForPoint(point)) { + module->GetPixelsForPoint(point); + } +} + +vector TRestDetectorExperimentalReadout::GetChannelsForPoint(const TVector3& point) const { + vector channels; + + for (const auto& pixel : GetPixelsForPoint(point)) { + channels.push_back(pixel->GetChannel()); + } + + return channels; +} diff --git a/src/TRestDetectorExperimentalReadoutChannel.cxx b/src/TRestDetectorExperimentalReadoutChannel.cxx new file mode 100644 index 00000000..f113732d --- /dev/null +++ b/src/TRestDetectorExperimentalReadoutChannel.cxx @@ -0,0 +1,6 @@ +#include "TRestDetectorExperimentalReadoutChannel.h" + +using namespace std; +using namespace readout; + +ClassImp(TRestDetectorExperimentalReadoutChannel); diff --git a/src/TRestDetectorExperimentalReadoutModule.cxx b/src/TRestDetectorExperimentalReadoutModule.cxx new file mode 100644 index 00000000..a3356e7c --- /dev/null +++ b/src/TRestDetectorExperimentalReadoutModule.cxx @@ -0,0 +1,244 @@ +#include "TRestDetectorExperimentalReadoutModule.h" + +using namespace std; +using namespace readout; + +ClassImp(TRestDetectorExperimentalReadoutModule); + +void TRestDetectorExperimentalReadoutModule::SetPixels( + const std::vector& pixels) { + fPixels = pixels; + fConvexHull = ComputeConvexHull(); + + // build map of pixels by channel (a pixel channel cannot be updated so this will remain current) + fChannelToPixels.clear(); + for (auto& pixel : fPixels) { + fChannelToPixels[pixel.GetChannel()].push_back(&pixel); + } + + // search radius is the maximum radius of all pixels + double searchRadius = 0; + for (auto& pixel : fPixels) { + if (pixel.GetRadius() > searchRadius) { + searchRadius = pixel.GetRadius(); + } + } + fSearchRadius = searchRadius; + + BuildKDTree(); +} + +vector TRestDetectorExperimentalReadoutModule::ComputeConvexHull() { + std::vector vertices; + + for (auto& pixel : fPixels) { + for (auto& vertex : pixel.GetVertices()) { + vertices.push_back(vertex); + } + } + + return findConvexHull(vertices); +} + +// KD Tree + +struct KDNode { + int index; // Index of the point in the original collection + KDNode* left; + KDNode* right; + + KDNode(int _index) : index(_index), left(nullptr), right(nullptr) {} +}; + +class KDTree { + private: + KDNode* root; + std::vector& points; + + KDNode* buildTree(std::vector& indices, int start, int end, int depth) { + if (start > end) { + return nullptr; + } + + int axis = depth % 2; // Split alternately on x and y coordinates + int mid = (start + end) / 2; + + // Sort the indices based on the current splitting axis + if (axis == 0) { + std::sort(indices.begin() + start, indices.begin() + end + 1, + [this](int a, int b) { return points[a].X() < points[b].X(); }); + } else { + std::sort(indices.begin() + start, indices.begin() + end + 1, + [this](int a, int b) { return points[a].Y() < points[b].Y(); }); + } + + KDNode* node = new KDNode(indices[mid]); + + node->left = buildTree(indices, start, mid - 1, depth + 1); + node->right = buildTree(indices, mid + 1, end, depth + 1); + + return node; + } + + void query(KDNode* node, const TVector2& target, double distance, std::vector& result) { + if (node == nullptr) { + return; + } + + double dist = (points[node->index] - target).Mod(); + if (dist <= distance) { + result.push_back(node->index); + } + + int axis = points[node->index].X() < target.X() + ? 0 + : 1; // Determine the axis to search based on target position + + if (axis == 0) { + query(node->left, target, distance, result); + + if (target.X() + distance > points[node->index].X()) { + query(node->right, target, distance, result); + } + } else { + query(node->right, target, distance, result); + + if (target.X() - distance < points[node->index].X()) { + query(node->left, target, distance, result); + } + } + } + + public: + KDTree(std::vector& _points) : points(_points) { + std::vector indices(points.size()); + for (size_t i = 0; i < points.size(); ++i) { + indices[i] = i; + } + root = buildTree(indices, 0, int(points.size()) - 1, 0); + } + + std::vector queryIndices(const TVector2& target, double distance) { + std::vector result; + query(root, target, distance, result); + return result; + } +}; + +void TRestDetectorExperimentalReadoutModule::BuildKDTree() { + std::vector centers; + centers.reserve(fPixels.size()); + for (const auto& pixel : fPixels) { + centers.push_back(pixel.GetCenter()); + } + + // auto tree = KDTree(centers); + + // delete fKDTree; + fKDTree = new KDTree(centers); +} + +std::vector +TRestDetectorExperimentalReadoutModule::GetPixelsForPoint(const TVector2& point) const { + auto indices = fKDTree->queryIndices(point, fSearchRadius); + cout << "Found " << indices.size() << " pixels in point" << endl; + + std::vector pixels; + pixels.reserve(indices.size()); + for (auto& index : indices) { + pixels.push_back(&fPixels[index]); + } + + return pixels; +} + +bool TRestDetectorExperimentalReadoutModule::IsInside(const TVector2& point) const { + return IsPointInsideConvexHull(point, fConvexHull); +} + +bool TRestDetectorExperimentalReadoutModule::IsInside(const TVector3& point) const { + const double z = GetZ(point); + if (z < 0 || z > fHeight) { + return false; + } + return IsInside(TransformToModuleCoordinates(point)); +} + +TVector2 TRestDetectorExperimentalReadoutModule::TransformToModuleCoordinates(const TVector3& point) const { + return { + (point - fPosition).Dot(fCoordinateAxes.first), + (point - fPosition).Dot(fCoordinateAxes.second), + }; +} + +TVector3 TRestDetectorExperimentalReadoutModule::TransformToAbsoluteCoordinates(const TVector2& point) const { + return fPosition + point.X() * fCoordinateAxes.first + point.Y() * fCoordinateAxes.second; +} + +double TRestDetectorExperimentalReadoutModule::GetZ(const TVector3& point) const { + (point - fPosition).Dot(fNormal); +} + +void TRestDetectorExperimentalReadoutModule::SetNormal(const TVector3& normal) { + fNormal = normal.Unit(); + UpdateAxes(); +} + +void TRestDetectorExperimentalReadoutModule::SetRotation(double rotationInDegrees) { + fRotation = rotationInDegrees; + UpdateAxes(); +} + +void TRestDetectorExperimentalReadoutModule::SetHeight(double height) { + if (height <= 0) { + cerr << "TRestDetectorExperimentalReadoutModule::SetHeight: Height must be positive" << endl; + return; + } + fHeight = height; +} + +void TRestDetectorExperimentalReadoutModule::UpdateAxes() { + // idempotent + + fCoordinateAxes.first = {1, 0, 0}; + fCoordinateAxes.second = {0, 1, 0}; + + // Check if fNormal is different from the original normal + const TVector3 originalNormal = {0, 0, 1}; + if (fNormal != originalNormal) { + // Calculate the rotation axis by taking the cross product between the original normal and fNormal + TVector3 rotationAxis = originalNormal.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())); + + // Rotate the axes around the rotation axis by the rotation angle + fCoordinateAxes.first.Rotate(rotationAngle, rotationAxis); + fCoordinateAxes.second.Rotate(rotationAngle, rotationAxis); + } + + // rotate around normal by rotation angle + fCoordinateAxes.first.Rotate(fRotation * TMath::DegToRad(), fNormal); + fCoordinateAxes.second.Rotate(fRotation * TMath::DegToRad(), fNormal); +} + +std::vector +TRestDetectorExperimentalReadoutModule::GetPixelsForPoint(const TVector3& point) const { + const TVector2 pointInModuleCoordinates = TransformToModuleCoordinates(point); + return GetPixelsForPoint(pointInModuleCoordinates); +} + +std::vector +TRestDetectorExperimentalReadoutModule::GetPixelsForChannel(unsigned short channel) const { + std::vector pixels; + // if channel is not in map, return + if (fChannelToPixels.find(channel) == fChannelToPixels.end()) { + return pixels; + } + + pixels.reserve(fChannelToPixels.at(channel).size()); + for (const auto& pixel : fChannelToPixels.at(channel)) { + pixels.push_back(pixel); + } + return pixels; +} diff --git a/src/TRestDetectorExperimentalReadoutPixel.cxx b/src/TRestDetectorExperimentalReadoutPixel.cxx new file mode 100644 index 00000000..c27df273 --- /dev/null +++ b/src/TRestDetectorExperimentalReadoutPixel.cxx @@ -0,0 +1,147 @@ +#include "TRestDetectorExperimentalReadoutPixel.h" + +using namespace std; +using namespace readout; + +ClassImp(TRestDetectorExperimentalReadoutPixel); + +TRestDetectorExperimentalReadoutPixel::TRestDetectorExperimentalReadoutPixel( + const std::vector& vertices, unsigned short channel) + : fChannel(channel) { + InitializeVertices(vertices); +} + +TRestDetectorExperimentalReadoutPixel::TRestDetectorExperimentalReadoutPixel( + const TVector2& center, const std::pair& size, unsigned short channel) + : TRestDetectorExperimentalReadoutPixel(GetRectangularVertices(center, size), channel) { + // rectangular pixel +} + +TRestDetectorExperimentalReadoutPixel::TRestDetectorExperimentalReadoutPixel(const TVector2& center, + double size, + unsigned short channel) + : TRestDetectorExperimentalReadoutPixel(GetRectangularVertices(center, {size, size}), channel) { + // square pixel +} + +void TRestDetectorExperimentalReadoutPixel::InitializeVertices(const std::vector& vertices) { + // at least three non-aligned vertices + if (vertices.size() < 3) { + throw std::invalid_argument("Pixel must have at least three vertices"); + } + // compute convex hull + std::vector convexHull = findConvexHull(vertices); + fVertices = convexHull; + + // compute center and radius + TVector2 center(0, 0); + for (auto& vertex : fVertices) { + center += vertex; + } + center *= 1. / double(fVertices.size()); + fCenter = center; + + // the radius is the maximum distance from the center to any vertex + double radius = 0; + for (auto& vertex : fVertices) { + double distance = (vertex - center).Mod(); + if (distance > radius) { + radius = distance; + } + } + fRadius = radius; +} + +std::vector TRestDetectorExperimentalReadoutPixel::GetRectangularVertices( + const TVector2& center, const std::pair& size) { + // rectangular pixel + std::vector vertices = {{center.X() - size.first / 2, center.Y() - size.second / 2}, + {center.X() + size.first / 2, center.Y() - size.second / 2}, + {center.X() + size.first / 2, center.Y() + size.second / 2}, + {center.X() - size.first / 2, center.Y() + size.second / 2}}; + + return vertices; +} + +bool TRestDetectorExperimentalReadoutPixel::IsInside(const TVector2& point) const { + return readout::IsPointInsideConvexHull(point, fVertices); +} + +// tools + +namespace readout { +// Helper function to calculate the cross product of two vectors +double crossProduct(const TVector2& A, const TVector2& B, const TVector2& C) { + return (B.X() - A.X()) * (C.Y() - A.Y()) - (B.Y() - A.Y()) * (C.X() - A.X()); +} + +// Helper function to calculate the squared distance between two points +double squaredDistance(const TVector2& A, const TVector2& B) { + double dx = B.X() - A.X(); + double dy = B.Y() - A.Y(); + return dx * dx + dy * dy; +} + +// Comparator function to compare two points based on their polar angle with respect to the anchor point +bool comparePoints(const TVector2& A, const TVector2& B, const TVector2& anchor) { + double cross = crossProduct(anchor, A, B); + if (cross == 0) { + // If two points have the same polar angle, choose the one that is closer to the anchor point + return squaredDistance(anchor, A) < squaredDistance(anchor, B); + } + return cross > 0; +} + +// Function to find the convex hull of a set of points using the Graham's scan algorithm +std::vector findConvexHull(const std::vector& _points) { + std::vector points = _points; + // Find the bottommost point as the anchor point + TVector2 anchor = points[0]; + for (const TVector2& point : points) { + if (point.Y() < anchor.Y() || (point.Y() == anchor.Y() && point.X() < anchor.X())) { + anchor = point; + } + } + + // Sort the points based on their polar angle with respect to the anchor point + std::sort(points.begin(), points.end(), + [&](const TVector2& A, const TVector2& B) { return comparePoints(A, B, anchor); }); + + // Build the convex hull + std::vector convexHull; + convexHull.push_back(points[0]); + convexHull.push_back(points[1]); + int hullSize = 2; + + for (int i = 2; i < points.size(); i++) { + while (hullSize >= 2 && + crossProduct(convexHull[hullSize - 2], convexHull[hullSize - 1], points[i]) <= 0) { + convexHull.pop_back(); + hullSize--; + } + convexHull.push_back(points[i]); + hullSize++; + } + + return convexHull; +} + +bool IsPointInsideConvexHull(const TVector2& point, const std::vector& convexHull) { + // Iterate through each edge of the convex hull + for (size_t i = 0; i < convexHull.size(); ++i) { + const TVector2& vertex1 = convexHull[i]; + const TVector2& vertex2 = convexHull[(i + 1) % convexHull.size()]; + + // Calculate the cross product of the point with the current edge + double crossProductValue = readout::crossProduct(vertex1, vertex2, point); + + // If the cross product is negative for any edge, the point is outside the convex hull + if (crossProductValue < 0.0) { + return false; + } + } + + // If the point is not outside any edge, it is inside the convex hull + return true; +} +} // namespace readout diff --git a/test/src/TRestDetectorReadout.cxx b/test/src/TRestDetectorReadout.cxx new file mode 100644 index 00000000..e521b3a5 --- /dev/null +++ b/test/src/TRestDetectorReadout.cxx @@ -0,0 +1,108 @@ + +#include +#include +#include +#include +#include + +#include + +namespace fs = std::filesystem; + +using namespace std; + +const auto filesPath = fs::path(__FILE__).parent_path().parent_path() / "files"; + +std::tuple GetGraphData(std::vector vertices) { + const auto n = vertices.size() + 1; + + double x[n]; + double y[n]; + + for (int i = 0; i < vertices.size(); i++) { + x[i] = vertices[i].X(); + y[i] = vertices[i].Y(); + } + + x[n - 1] = vertices[0].X(); + y[n - 1] = vertices[0].Y(); + + return {n, x, y}; +} + +TEST(TRestDetectorReadout, Initialize) { + auto readout = TRestDetectorExperimentalReadout(); + + auto module = TRestDetectorExperimentalReadoutModule(); + readout.AddModule(module); + + double pixelSize = 1.0; + int pixelsPerSide = 25; + + std::vector pixels; + for (int i = 0; i < pixelsPerSide; i++) { + for (int j = 0; j < pixelsPerSide; j++) { + TVector2 center = {i * pixelSize, j * pixelSize}; + pixels.emplace_back(center, pixelSize); + } + } + + module.SetPixels(pixels); + + module.GetPixelsForPoint({0, 0}); + return; + + auto canvas = new TCanvas("canvas", "Shape Canvas", 800, 600); + + // first draw the convex hull for the module + const auto hullVertices = module.GetConvexHull(); + + TGraph* graph; + { + const auto& vertices = hullVertices; + const auto n = vertices.size() + 1; + double x[n]; + double y[n]; + + for (int i = 0; i < n - 1; i++) { + x[i] = vertices[i].X(); + y[i] = vertices[i].Y(); + } + x[n - 1] = vertices[0].X(); + y[n - 1] = vertices[0].Y(); + graph = new TGraph(n, x, y); + } + graph->Draw(); + + // print pixel info + for (int i = 0; i < module.GetNumberOfPixels(); i++) { + const auto pixel = module.GetPixels()[i]; + // cout << "Pixel X: " << pixel.GetCenter().X() << " Y: " << pixel.GetCenter().Y() << endl; + const auto vertices = pixel.GetVertices(); + + const auto n = vertices.size() + 1; + double x[n]; + double y[n]; + + for (int i = 0; i < n - 1; i++) { + x[i] = vertices[i].X(); + y[i] = vertices[i].Y(); + } + x[n - 1] = vertices[0].X(); + y[n - 1] = vertices[0].Y(); + + for (int i = 0; i < n; i++) { + // cout << "i: " << i << " X: " << x[i] << " Y: " << y[i] << endl; + } + + auto graph = new TGraph(vertices.size() + 1, x, y); + + graph->SetLineColor(kRed); + graph->SetLineWidth(2); + + graph->Draw("SAME"); + canvas->Update(); + } + gPad->SetFixedAspectRatio(); + canvas->Update(); +}