Skip to content

GitLab

Commit 222a7f2e authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Merge branch 'split-planar-detector' into 'master' 

Split planar detector

See merge request !263
parents 283e94b6 1ae38250
Pipeline #1527473 passed with stages
in 22 minutes and 49 seconds
Hide whitespace changes
Inline Side-by-side
doc/usermanual/chapters/configuration.tex
View file @ 222a7f2e
......@@ -281,6 +281,8 @@ which describes a rotation of \SI{45}{\degree} around the $Z$ axis, followed by
All supported rotations are extrinsic active rotations, i.e. the vector itself is rotated, not the coordinate system. All angles in configuration files should be specified in the order they will be applied.
\end{warning}
\item The \parameter{coordinates} parameter represents the local coordinate of detectors, the coordinates can be \textbf{cartesian}, \textbf{polar}, etc. The default \parameter{coordinates} is \textbf{cartesian}.
\item The \parameter{number_of_pixels} parameter represents a two-dimensional vector with the number of pixels in the active matrix in the column and row directions respectively.
\item The \parameter{pixel_pitch} is a two-dimensional vector defining the size of a single pixel in the column and row directions respectively.
\item The intrinsic resolution of the detector has to be specified using the \parameter{spatial_resolution} parameter, a two-dimensional vector holding the position resolution for the column and row directions. This value is used to assign the uncertainty of cluster positions. This parameter defaults to the pitch$/\sqrt{12}$ of the respective detector if not specified.
......
doc/usermanual/chapters/framework.tex
View file @ 222a7f2e
......@@ -238,6 +238,6 @@ The global coordinate system is chosen as a right-handed Cartesian system, and t
Here, the beam direction defines the positive z-axis at the origin of the x-y-plane.
The origin along the z-axis is fixed by the placement of the detectors in the geometry of the setup.
Local coordinate systems for the detectors are also right-handed Cartesian systems, with the x- and y-axes defining the sensor plane.
The origin of this coordinate system is the center of the lower left pixel in the grid, i.e.\ the pixel with indices (0,0).
This simplifies calculations in the local coordinate system as all positions can either be stated in absolute numbers or in fractions of the pixel pitch.
There are several kinds of detectors with different local coordinate system, e.g., Cartesian coordinate for pixel detector and polar coordinate for some specific strip detectors. For this reason, the \textsl{Detector} class is designed as an interface, the detector with a specific coordinate system will be generated according to the parameter in the detector configuration file. In the current version, the Cartesian coordinate is the only available one.
For the local Cartesian coordinate system, it is right-handed with the x- and y-axes defining the sensor plane. The origin of this coordinate system is the center of the lower left pixel in the grid, i.e.\ the pixel with indices (0,0). This simplifies calculations in the local coordinate system as all positions can either be stated in absolute numbers or in fractions of the pixel pitch.
src/core/CMakeLists.txt
View file @ 222a7f2e
......@@ -5,6 +5,7 @@ INCLUDE_DIRECTORIES(SYSTEM ${CORRYVRECKAN_DEPS_INCLUDE_DIRS})
ADD_LIBRARY(CorryvreckanCore SHARED
Corryvreckan.cpp
detector/Detector.cpp
detector/PixelDetector.cpp
detector/exceptions.cpp
utils/log.cpp
utils/unit.cpp
......
src/core/detector/Detector.cpp
View file @ 222a7f2e
......@@ -46,14 +46,6 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {
throw InvalidValueError(config, "role", "Auxiliary devices cannot hold any other detector role");
}
// Detector position and orientation
m_displacement = config.get<ROOT::Math::XYZPoint>("position", ROOT::Math::XYZPoint());
m_orientation = config.get<ROOT::Math::XYZVector>("orientation", ROOT::Math::XYZVector());
m_orientation_mode = config.get<std::string>("orientation_mode", "xyz");
if(m_orientation_mode != "xyz" && m_orientation_mode != "zyx") {
throw InvalidValueError(config, "orientation_mode", "Invalid detector orientation mode");
}
m_detectorName = config.getName();
// Material budget of detector, including support material
......@@ -66,49 +58,16 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {
m_materialBudget = config.get<double>("material_budget");
}
// Auxiliary devices don't have: number_of_pixels, pixel_pitch, spatial_resolution, mask_file, region-of-interest
if(!isAuxiliary()) {
// Number of pixels:
m_nPixels = config.get<ROOT::Math::DisplacementVector2D<Cartesian2D<int>>>("number_of_pixels");
// Size of the pixels:
m_pitch = config.get<ROOT::Math::XYVector>("pixel_pitch");
if(Units::convert(m_pitch.X(), "mm") >= 1 or Units::convert(m_pitch.Y(), "mm") >= 1 or
Units::convert(m_pitch.X(), "um") <= 1 or Units::convert(m_pitch.Y(), "um") <= 1) {
LOG(WARNING) << "Pixel pitch unphysical for detector " << m_detectorName << ": " << std::endl
<< Units::display(m_pitch, {"nm", "um", "mm"});
}
// Intrinsic spatial resolution, defaults to pitch/sqrt(12):
m_spatial_resolution = config.get<ROOT::Math::XYVector>("spatial_resolution", m_pitch / std::sqrt(12));
if(!config.has("spatial_resolution")) {
LOG(WARNING) << "Spatial resolution for detector '" << m_detectorName << "' not set." << std::endl
<< "Using pitch/sqrt(12) as default";
}
// region of interest:
m_roi = config.getMatrix<int>("roi", std::vector<std::vector<int>>());
}
m_detectorType = config.get<std::string>("type");
std::transform(m_detectorType.begin(), m_detectorType.end(), m_detectorType.begin(), ::tolower);
m_timeOffset = config.get<double>("time_offset", 0.0);
// Time resolution - default ot negative number, i.e. unknown. This will trigger an exception
// when calling getTimeResolution
m_timeResolution = config.get<double>("time_resolution", -1.0);
// Initialize the detector, calculate transformations etc
this->initialise();
LOG(TRACE) << "Initialized \"" << m_detectorType << "\": " << m_nPixels.X() << "x" << m_nPixels.Y() << " px, pitch of "
<< Units::display(m_pitch, {"mm", "um"});
LOG(TRACE) << " Position: " << Units::display(m_displacement, {"mm", "um"});
LOG(TRACE) << " Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")";
if(m_timeOffset > 0.) {
LOG(TRACE) << "Time offset: " << m_timeOffset;
}
// Time resolution - default to negative number, i.e. unknown. This will trigger an exception
// when calling getTimeResolution
m_timeResolution = config.get<double>("time_resolution", -1.0);
if(m_timeResolution > 0) {
LOG(TRACE) << " Time resolution: " << Units::display(m_timeResolution, {"ms", "us"});
}
......@@ -116,15 +75,16 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {
if(config.has("calibration_file")) {
m_calibrationfile = config.getPath("calibration_file");
}
}
if(!isAuxiliary()) {
if(config.has("mask_file")) {
m_maskfile_name = config.get<std::string>("mask_file");
std::string mask_file = config.getPath("mask_file");
LOG(DEBUG) << "Adding mask to detector \"" << config.getName() << "\", reading from " << mask_file;
setMaskFile(mask_file);
processMaskFile();
}
std::shared_ptr<Detector> corryvreckan::Detector::Factory(const Configuration& config) {
// default coordinate is cartesian coordinate
std::string coordinates = config.get<std::string>("coordinates", "cartesian");
std::transform(coordinates.begin(), coordinates.end(), coordinates.begin(), ::tolower);
if(coordinates == "cartesian") {
return std::make_shared<PixelDetector>(config);
} else {
throw InvalidValueError(config, "coordinates", "Coordiantes can only set to be cartesian now");
}
}
......@@ -144,10 +104,6 @@ std::string Detector::type() const {
return m_detectorType;
}
XYVector Detector::size() const {
return XYVector(m_pitch.X() * m_nPixels.X(), m_pitch.Y() * m_nPixels.Y());
}
bool Detector::isReference() const {
return static_cast<bool>(m_role & DetectorRole::REFERENCE);
}
......@@ -165,96 +121,12 @@ void Detector::setMaskFile(std::string file) {
m_maskfile = file;
}
void Detector::processMaskFile() {
// Open the file with masked pixels
std::ifstream inputMaskFile(m_maskfile, std::ios::in);
if(!inputMaskFile.is_open()) {
LOG(WARNING) << "Could not open mask file " << m_maskfile;
} else {
int row = 0, col = 0;
std::string id;
// loop over all lines and apply masks
while(inputMaskFile >> id) {
if(id == "c") {
inputMaskFile >> col;
if(col > nPixels().X() - 1) {
LOG(WARNING) << "Column " << col << " outside of pixel matrix, chip has only " << nPixels().X()
<< " columns!";
}
LOG(TRACE) << "Masking column " << col;
for(int r = 0; r < nPixels().Y(); r++) {
maskChannel(col, r);
}
} else if(id == "r") {
inputMaskFile >> row;
if(row > nPixels().Y() - 1) {
LOG(WARNING) << "Row " << col << " outside of pixel matrix, chip has only " << nPixels().Y() << " rows!";
}
LOG(TRACE) << "Masking row " << row;
for(int c = 0; c < nPixels().X(); c++) {
maskChannel(c, row);
}
} else if(id == "p") {
inputMaskFile >> col >> row;
if(col > nPixels().X() - 1 || row > nPixels().Y() - 1) {
LOG(WARNING) << "Pixel " << col << " " << row << " outside of pixel matrix, chip has only "
<< nPixels().X() << " x " << nPixels().Y() << " pixels!";
}
LOG(TRACE) << "Masking pixel " << col << " " << row;
maskChannel(col, row); // Flag to mask a pixel
} else {
LOG(WARNING) << "Could not parse mask entry (id \"" << id << "\")";
}
}
LOG(INFO) << m_masked.size() << " masked pixels";
}
}
void Detector::maskChannel(int chX, int chY) {
int channelID = chX + m_nPixels.X() * chY;
m_masked[channelID] = true;
}
bool Detector::masked(int chX, int chY) const {
int channelID = chX + m_nPixels.X() * chY;
if(m_masked.count(channelID) > 0)
return true;
return false;
}
// Function to initialise transforms
void Detector::initialise() {
// Make the local to global transform, built from a displacement and
// rotation
Translation3D translations = Translation3D(m_displacement.X(), m_displacement.Y(), m_displacement.Z());
Rotation3D rotations;
if(m_orientation_mode == "xyz") {
rotations = RotationZ(m_orientation.Z()) * RotationY(m_orientation.Y()) * RotationX(m_orientation.X());
} else if(m_orientation_mode == "zyx") {
rotations = RotationZYX(m_orientation.x(), m_orientation.y(), m_orientation.x());
}
m_localToGlobal = Transform3D(rotations, translations);
m_globalToLocal = m_localToGlobal.Inverse();
// Find the normal to the detector surface. Build two points, the origin and a unit step in z,
// transform these points to the global co-ordinate frame and then make a vector pointing between them
m_origin = PositionVector3D<Cartesian3D<double>>(0., 0., 0.);
m_origin = m_localToGlobal * m_origin;
PositionVector3D<Cartesian3D<double>> localZ(0., 0., 1.);
localZ = m_localToGlobal * localZ;
m_normal = PositionVector3D<Cartesian3D<double>>(
localZ.X() - m_origin.X(), localZ.Y() - m_origin.Y(), localZ.Z() - m_origin.Z());
}
// Function to update transforms (such as during alignment)
void Detector::update() {
this->initialise();
}
Configuration Detector::getConfiguration() const {
Configuration Detector::GetConfiguration() const {
Configuration config(name());
config.set("type", m_detectorType);
......@@ -275,243 +147,24 @@ Configuration Detector::getConfiguration() const {
config.setArray("role", roles);
}
config.set("position", m_displacement, {"um", "mm"});
config.set("orientation_mode", m_orientation_mode);
config.set("orientation", m_orientation, {"deg"});
if(m_timeOffset != 0.) {
config.set("time_offset", m_timeOffset, {"ns", "us", "ms", "s"});
}
config.set("time_resolution", m_timeResolution, {"ns", "us", "ms", "s"});
// different for PixelDetector and StripDetector
this->configurePosAndOrientation(config);
// material budget
if(m_materialBudget > std::numeric_limits<double>::epsilon()) {
config.set("material_budget", m_materialBudget);
}
// only if detector is not auxiliary:
if(!this->isAuxiliary()) {
config.set("number_of_pixels", m_nPixels);
// Size of the pixels
config.set("pixel_pitch", m_pitch, {"um"});
// Intrinsic resolution:
config.set("spatial_resolution", m_spatial_resolution, {"um"});
// Pixel mask file:
if(!m_maskfile_name.empty()) {
config.set("mask_file", m_maskfile_name);
}
// Region-of-interest:
config.setMatrix("roi", m_roi);
this->configureDetector(config);
}
return config;
}
// Function to get global intercept with a track
PositionVector3D<Cartesian3D<double>> Detector::getIntercept(const Track* track) const {
// FIXME: this is else statement can only be temporary
if(track->getType() == "gbl") {
return track->state(name());
} else {
// Get the distance from the plane to the track initial state
double distance = (m_origin.X() - track->state(m_detectorName).X()) * m_normal.X();
distance += (m_origin.Y() - track->state(m_detectorName).Y()) * m_normal.Y();
distance += (m_origin.Z() - track->state(m_detectorName).Z()) * m_normal.Z();
distance /=
(track->direction(m_detectorName).X() * m_normal.X() + track->direction(m_detectorName).Y() * m_normal.Y() +
track->direction(m_detectorName).Z() * m_normal.Z());
// Propagate the track
PositionVector3D<Cartesian3D<double>> globalIntercept(
track->state(m_detectorName).X() + distance * track->direction(m_detectorName).X(),
track->state(m_detectorName).Y() + distance * track->direction(m_detectorName).Y(),
track->state(m_detectorName).Z() + distance * track->direction(m_detectorName).Z());
return globalIntercept;
}
}
PositionVector3D<Cartesian3D<double>> Detector::getLocalIntercept(const Track* track) const {
return globalToLocal(getIntercept(track));
}
// Function to check if a track intercepts with a plane
bool Detector::hasIntercept(const Track* track, double pixelTolerance) const {
// First, get the track intercept in global co-ordinates with the plane
PositionVector3D<Cartesian3D<double>> globalIntercept = this->getIntercept(track);
// Convert to local co-ordinates
PositionVector3D<Cartesian3D<double>> localIntercept = this->m_globalToLocal * globalIntercept;
// Get the row and column numbers
double row = this->getRow(localIntercept);
double column = this->getColumn(localIntercept);
// Check if the row and column are outside of the chip
// Chip reaches from -0.5 to nPixels-0.5
bool intercept = true;
if(row < pixelTolerance - 0.5 || row > (this->m_nPixels.Y() - pixelTolerance - 0.5) || column < pixelTolerance - 0.5 ||
column > (this->m_nPixels.X() - pixelTolerance - 0.5))
intercept = false;
return intercept;
}
// Function to check if a track goes through/near a masked pixel
bool Detector::hitMasked(Track* track, int tolerance) const {
// First, get the track intercept in global co-ordinates with the plane
PositionVector3D<Cartesian3D<double>> globalIntercept = this->getIntercept(track);
// Convert to local co-ordinates
PositionVector3D<Cartesian3D<double>> localIntercept = this->m_globalToLocal * globalIntercept;
// Get the row and column numbers
int row = static_cast<int>(floor(this->getRow(localIntercept) + 0.5));
int column = static_cast<int>(floor(this->getColumn(localIntercept) + 0.5));
// Check if the pixels around this pixel are masked
bool hitmasked = false;
for(int r = (row - tolerance); r <= (row + tolerance); r++) {
for(int c = (column - tolerance); c <= (column + tolerance); c++) {
if(this->masked(c, r))
hitmasked = true;
}
}
return hitmasked;
}
// Functions to get row and column from local position
double Detector::getRow(const PositionVector3D<Cartesian3D<double>> localPosition) const {
double row = localPosition.Y() / m_pitch.Y() + static_cast<double>(m_nPixels.Y() - 1) / 2.;
return row;
}
double Detector::getColumn(const PositionVector3D<Cartesian3D<double>> localPosition) const {
double column = localPosition.X() / m_pitch.X() + static_cast<double>(m_nPixels.X() - 1) / 2.;
return column;
}
// Function to get local position from row and column
PositionVector3D<Cartesian3D<double>> Detector::getLocalPosition(double column, double row) const {
return PositionVector3D<Cartesian3D<double>>(
m_pitch.X() * (column - (m_nPixels.X() - 1) / 2.), m_pitch.Y() * (row - (m_nPixels.Y() - 1) / 2.), 0.);
}
// Function to get in-pixel position
ROOT::Math::XYVector Detector::inPixel(const double column, const double row) const {
// a pixel ranges from (col-0.5) to (col+0.5)
return XYVector(m_pitch.X() * (column - floor(column) - 0.5), m_pitch.Y() * (row - floor(row) - 0.5));
}
ROOT::Math::XYVector Detector::inPixel(const PositionVector3D<Cartesian3D<double>> localPosition) const {
double column = getColumn(localPosition);
double row = getRow(localPosition);
return inPixel(column, row);
}
// Check if track position is within ROI:
bool Detector::isWithinROI(const Track* track) const {
// Empty region of interest:
if(m_roi.empty()) {
return true;
}
// Check that track is within region of interest using winding number algorithm
auto localIntercept = this->getLocalIntercept(track);
auto coordinates = std::make_pair(this->getColumn(localIntercept), this->getRow(localIntercept));
if(winding_number(coordinates, m_roi) != 0) {
return true;
}
// Outside ROI:
return false;
}
// Check if cluster is within ROI and/or touches ROI border:
bool Detector::isWithinROI(Cluster* cluster) const {
// Empty region of interest:
if(m_roi.empty()) {
return true;
}
// Loop over all pixels of the cluster
for(auto& pixel : cluster->pixels()) {
if(winding_number(pixel->coordinates(), m_roi) == 0) {
return false;
}
}
return true;
}
/* isLeft(): tests if a point is Left|On|Right of an infinite line.
* via: http://geomalgorithms.com/a03-_inclusion.html
* Input: three points P0, P1, and P2
* Return: >0 for P2 left of the line through P0 and P1
* =0 for P2 on the line
* <0 for P2 right of the line
* See: Algorithm 1 "Area of Triangles and Polygons"
*/
int Detector::isLeft(std::pair<int, int> pt0, std::pair<int, int> pt1, std::pair<int, int> pt2) {
return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second));
}
/* Winding number test for a point in a polygon
* via: http://geomalgorithms.com/a03-_inclusion.html
* Input: x, y = a point,
* polygon = vector of vertex points of a polygon V[n+1] with V[n]=V[0]
* Return: wn = the winding number (=0 only when P is outside)
*/
int Detector::winding_number(std::pair<int, int> probe, std::vector<std::vector<int>> polygon) {
// Two points don't make an area
if(polygon.size() < 3) {
LOG(DEBUG) << "No ROI given.";
return 0;
}
int wn = 0; // the winding number counter
// loop through all edges of the polygon
// edge from V[i] to V[i+1]
for(size_t i = 0; i < polygon.size(); i++) {
auto point_this = std::make_pair(polygon.at(i).at(0), polygon.at(i).at(1));
auto point_next = (i + 1 < polygon.size() ? std::make_pair(polygon.at(i + 1).at(0), polygon.at(i + 1).at(1))
: std::make_pair(polygon.at(0).at(0), polygon.at(0).at(1)));
// start y <= P.y
if(point_this.second <= probe.second) {
// an upward crossing
if(point_next.second > probe.second) {
// P left of edge
if(isLeft(point_this, point_next, probe) > 0) {
// have a valid up intersect
++wn;
}
}
} else {
// start y > P.y (no test needed)
// a downward crossing
if(point_next.second <= probe.second) {
// P right of edge
if(isLeft(point_this, point_next, probe) < 0) {
// have a valid down intersect
--wn;
}
}
}
}
return wn;
}
src/core/detector/Detector.hpp
View file @ 222a7f2e
......@@ -57,9 +57,9 @@ namespace corryvreckan {
}
/**
* @brief Detector representation in the reconstruction chain
* @brief Detector interface in the reconstruction chain
*
* Contains the detector with all its properties such as type, name, position and orientation, pitch, spatial resolution
* Contains the detector with common properties such as type, name, coordinate
* etc.
*/
class Detector {
......@@ -69,12 +69,24 @@ namespace corryvreckan {
*/
Detector() = delete;
/**
* Default destructor
*/
virtual ~Detector() = default;
/**
* @brief Constructs a detector in the geometry
* @param config Configuration object describing the detector
*/
Detector(const Configuration& config);
/**
* @brief Factory to dynamically create detectors
* @param config Configuration object describing the detector
* @return shared_ptr that contains the real detector
*/
static std::shared_ptr<Detector> Factory(const Configuration& config);
/**
* @brief Get type of the detector
* @return Type of the detector model
......@@ -109,31 +121,35 @@ namespace corryvreckan {
* @brief Retrieve configuration object from detector, containing all (potentially updated) parameters
* @return Configuration object for this detector
*/
Configuration getConfiguration() const;
Configuration GetConfiguration() const;
/**
* @brief Get the total size of the active matrix, i.e. pitch * number of pixels in both dimensions
* @return 2D vector with the dimensions of the pixle matrix in X and Y
* @to do: this is designed for PixelDetector, find a proper interface for other Detector type
*/
XYVector size() const;
virtual XYVector size() const = 0;
/**
* @brief Get pitch of a single pixel
* @return Pitch of a pixel
* @to do: this is designed for PixelDetector, find a proper interface for other Detector type
*/
XYVector pitch() const { return m_pitch; }
virtual XYVector pitch() const = 0;
/**
* @brief Get intrinsic spatial resolution of the detector
* @return Intrinsic spatial resolution in X and Y
* @to do: this is designed for PixelDetector, find a proper interface for other Detector type
*/
XYVector getSpatialResolution() const { return m_spatial_resolution; }
virtual XYVector getSpatialResolution() const = 0;
/**
* @brief Get number of pixels in x and y
* @return Number of two dimensional pixels
* @to do: this is designed for PixelDetector, find a proper interface for other Detector type
*/
ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<int>> nPixels() const { return m_nPixels; }
virtual ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<int>> nPixels() const = 0;
/**
* @brief Get detector time offset from global clock, can be used to correct for constant shifts or time of flight
......@@ -151,31 +167,37 @@ namespace corryvreckan {
* @brief Update detector position in the world
* @param displacement Vector with three position coordinates
*/
void displacement(XYZPoint displacement) { m_displacement = displacement; }
virtual void displacement(XYZPoint displacement) = 0;
/**
* @brief Get position in the world
* @return Global position in Cartesian coordinates
*/
XYZPoint displacement() const { return m_displacement; }
virtual XYZPoint displacement() const = 0;
/**
* @brief Get orientation in the world
* @return Vector with three rotation angles
*/
XYZVector rotation() const { return m_orientation; }
virtual XYZVector rotation() const = 0;
/**
* @brief Update detector orientation in the world
* @param rotation Vector with three rotation angles
*/
void rotation(XYZVector rotation) { m_orientation = rotation; }
virtual void rotation(XYZVector rotation) = 0;
/**
* @brief Get normal vector to sensor surface
* @return Normal vector to sensor surface