From 612dbfc299c75e7171aef6027819bbdd9d96abef Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Thu, 18 Jun 2020 11:25:16 +0200 Subject: [PATCH 01/37] applied shearz patch from Magnus Mager --- src/modules/AlignmentMillepede/AlignmentMillepede.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp index c0645914..f99573b7 100644 --- a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp +++ b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp @@ -218,6 +218,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { std::vector fscaz(nParameters, 0.); std::vector shearx(nParameters, 0.); std::vector sheary(nParameters, 0.); + std::vector shearrz(nParameters, 0.); m_constraints.clear(); for(const auto& det : get_detectors()) { @@ -239,9 +240,10 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { shearx[i] = sz; sheary[i + nPlanes] = sz; fscaz[i + 2 * nPlanes] = sz; + shearrz[i + 5 * nPlanes] = sz; } - const std::vector constraints = {true, true, true, true, false, false, true, true, true}; + const std::vector constraints = {true, true, true, true, false, false, true, true, true, true}; // Put the constraints information in the basket. if(constraints[0] && m_dofs[0]) addConstraint(ftx, 0.0); @@ -260,6 +262,8 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { addConstraint(fry, 0.0); if(constraints[8] && m_dofs[5]) addConstraint(frz, 0.0); + if(constraints[9] && m_dofs[5]) + addConstraint(shearrz, 0.0); } //============================================================================= -- GitLab From 0798994657fa0f99ed2991ffb95c2218c292aa53 Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Thu, 18 Jun 2020 11:25:42 +0200 Subject: [PATCH 02/37] added bentpixeldetector, copied from pixeldetector --- src/core/detector/BentPixelDetector.cpp | 411 ++++++++++++++++++++++++ src/core/detector/BentPixelDetector.hpp | 205 ++++++++++++ 2 files changed, 616 insertions(+) create mode 100644 src/core/detector/BentPixelDetector.cpp create mode 100644 src/core/detector/BentPixelDetector.hpp diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp new file mode 100644 index 00000000..8874b1e0 --- /dev/null +++ b/src/core/detector/BentPixelDetector.cpp @@ -0,0 +1,411 @@ +/** @file + * @brief Implementation of the detector model + * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors. + * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md". + * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an + * Intergovernmental Organization or submit itself to any jurisdiction. + */ + +#include +#include +#include + +#include "Math/RotationX.h" +#include "Math/RotationY.h" +#include "Math/RotationZ.h" +#include "Math/RotationZYX.h" + +#include "PixelDetector.hpp" +#include "core/utils/log.h" +#include "exceptions.h" + +using namespace ROOT::Math; +using namespace corryvreckan; + +PixelDetector::PixelDetector(const Configuration& config) : Detector(config) { + + // Set detector position and direction from configuration file + SetPostionAndOrientation(config); + + // initialize transform + this->initialise(); + + // Auxiliary devices don't have: number_of_pixels, pixel_pitch, spatial_resolution, mask_file, region-of-interest + if(!isAuxiliary()) { + build_axes(config); + } +} + +void PixelDetector::build_axes(const Configuration& config) { + + m_nPixels = config.get>>("number_of_pixels"); + // Size of the pixels: + m_pitch = config.get("pixel_pitch"); + + LOG(TRACE) << "Initialized \"" << m_detectorType << "\": " << m_nPixels.X() << "x" << m_nPixels.Y() << " px, pitch of " + << Units::display(m_pitch, {"mm", "um"}); + + 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("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("roi", std::vector>()); + + if(config.has("mask_file")) { + m_maskfile_name = config.get("mask_file"); + std::string mask_file = config.getPath("mask_file", true); + LOG(DEBUG) << "Adding mask to detector \"" << config.getName() << "\", reading from " << mask_file; + set_mask_file(mask_file); + process_mask_file(); + } +} + +void PixelDetector::SetPostionAndOrientation(const Configuration& config) { + // Detector position and orientation + m_displacement = config.get("position", ROOT::Math::XYZPoint()); + m_orientation = config.get("orientation", ROOT::Math::XYZVector()); + m_orientation_mode = config.get("orientation_mode", "xyz"); + if(m_orientation_mode != "xyz" && m_orientation_mode != "zyx") { + throw InvalidValueError(config, "orientation_mode", "Invalid detector orientation mode"); + } + LOG(TRACE) << " Position: " << Units::display(m_displacement, {"mm", "um"}); + LOG(TRACE) << " Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")"; +} + +void PixelDetector::process_mask_file() { + // 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 PixelDetector::maskChannel(int chX, int chY) { + int channelID = chX + m_nPixels.X() * chY; + m_masked[channelID] = true; +} + +bool PixelDetector::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 PixelDetector::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") { + LOG(DEBUG) << "Interpreting Euler angles as XYZ rotation"; + // First angle given in the configuration file is around x, second around y, last around z: + rotations = RotationZ(m_orientation.Z()) * RotationY(m_orientation.Y()) * RotationX(m_orientation.X()); + } else if(m_orientation_mode == "zyx") { + LOG(DEBUG) << "Interpreting Euler angles as ZYX rotation"; + // First angle given in the configuration file is around z, second around y, last around x: + rotations = RotationZYX(m_orientation.x(), m_orientation.y(), m_orientation.z()); + } else if(m_orientation_mode == "zxz") { + LOG(DEBUG) << "Interpreting Euler angles as ZXZ rotation"; + // First angle given in the configuration file is around z, second around x, last around z: + rotations = EulerAngles(m_orientation.x(), m_orientation.y(), m_orientation.z()); + } else { + LOG(ERROR) << "orientation_mode should be either 'zyx', xyz' or 'zxz'"; + // FIXME: To through exception, initialise needs to be changed to "initialise(const Configuration& config)" + // throw InvalidValueError( + // config, "orientation_mode", "orientation_mode should be either 'zyx', xyz' or 'zxz'"); + } + + 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 coordinate frame and then make a vector pointing between them + m_origin = PositionVector3D>(0., 0., 0.); + m_origin = m_localToGlobal * m_origin; + PositionVector3D> localZ(0., 0., 1.); + localZ = m_localToGlobal * localZ; + m_normal = PositionVector3D>( + localZ.X() - m_origin.X(), localZ.Y() - m_origin.Y(), localZ.Z() - m_origin.Z()); +} + +// Only if detector is not auxiliary +void PixelDetector::configure_detector(Configuration& config) const { + + // Number of pixels + 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); +} + +void PixelDetector::configure_pos_and_orientation(Configuration& config) const { + config.set("position", m_displacement, {"um", "mm"}); + config.set("orientation_mode", m_orientation_mode); + config.set("orientation", m_orientation, {"deg"}); +} + +// Function to get global intercept with a track +PositionVector3D> PixelDetector::getIntercept(const Track* track) const { + + // FIXME: this is else statement can only be temporary + if(track->getType() == "GblTrack") { + return track->getState(getName()); + } else { + // Get the distance from the plane to the track initial state + double distance = (m_origin.X() - track->getState(m_detectorName).X()) * m_normal.X(); + distance += (m_origin.Y() - track->getState(m_detectorName).Y()) * m_normal.Y(); + distance += (m_origin.Z() - track->getState(m_detectorName).Z()) * m_normal.Z(); + distance /= (track->getDirection(m_detectorName).X() * m_normal.X() + + track->getDirection(m_detectorName).Y() * m_normal.Y() + + track->getDirection(m_detectorName).Z() * m_normal.Z()); + + // Propagate the track + PositionVector3D> globalIntercept( + track->getState(m_detectorName).X() + distance * track->getDirection(m_detectorName).X(), + track->getState(m_detectorName).Y() + distance * track->getDirection(m_detectorName).Y(), + track->getState(m_detectorName).Z() + distance * track->getDirection(m_detectorName).Z()); + return globalIntercept; + } +} + +PositionVector3D> PixelDetector::getLocalIntercept(const Track* track) const { + return globalToLocal(getIntercept(track)); +} + +// Function to check if a track intercepts with a plane +bool PixelDetector::hasIntercept(const Track* track, double pixelTolerance) const { + + // First, get the track intercept in global coordinates with the plane + PositionVector3D> globalIntercept = this->getIntercept(track); + + // Convert to local coordinates + PositionVector3D> 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 PixelDetector::hitMasked(const Track* track, int tolerance) const { + + // First, get the track intercept in global coordinates with the plane + PositionVector3D> globalIntercept = this->getIntercept(track); + + // Convert to local coordinates + PositionVector3D> localIntercept = this->m_globalToLocal * globalIntercept; + + // Get the row and column numbers + int row = static_cast(floor(this->getRow(localIntercept) + 0.5)); + int column = static_cast(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 PixelDetector::getRow(const PositionVector3D> localPosition) const { + return localPosition.Y() / m_pitch.Y() + static_cast(m_nPixels.Y() - 1) / 2.; +} +double PixelDetector::getColumn(const PositionVector3D> localPosition) const { + return localPosition.X() / m_pitch.X() + static_cast(m_nPixels.X() - 1) / 2.; +} + +// Function to get local position from row and column +PositionVector3D> PixelDetector::getLocalPosition(double column, double row) const { + + return PositionVector3D>(m_pitch.X() * (column - static_cast(m_nPixels.X() - 1) / 2.), + m_pitch.Y() * (row - static_cast(m_nPixels.Y() - 1) / 2.), + 0.); +} + +// Function to get in-pixel position +ROOT::Math::XYVector PixelDetector::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 PixelDetector::inPixel(const PositionVector3D> localPosition) const { + double column = getColumn(localPosition); + double row = getRow(localPosition); + return inPixel(column, row); +} + +// Check if track position is within ROI: +bool PixelDetector::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 PixelDetector::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; +} + +XYVector PixelDetector::getSize() const { + return XYVector(m_pitch.X() * m_nPixels.X(), m_pitch.Y() * m_nPixels.Y()); +} + +/* 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 PixelDetector::winding_number(std::pair probe, std::vector> 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; +} +/* 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 PixelDetector::isLeft(std::pair pt0, std::pair pt1, std::pair pt2) { + return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second)); +} diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp new file mode 100644 index 00000000..eab25bf1 --- /dev/null +++ b/src/core/detector/BentPixelDetector.hpp @@ -0,0 +1,205 @@ +/** @file + * @brief Detector model class + * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors. + * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md". + * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an + * Intergovernmental Organization or submit itself to any jurisdiction. + */ + +#ifndef CORRYVRECKAN_PLANARDETECTOR_H +#define CORRYVRECKAN_PLANARDETECTOR_H + +#include +#include +#include + +#include +#include +#include +#include "Math/Transform3D.h" +#include "Math/Vector3D.h" + +#include "Detector.hpp" +#include "core/config/Configuration.hpp" +#include "core/utils/ROOT.h" +#include "core/utils/log.h" +#include "objects/Track.hpp" + +namespace corryvreckan { + + /** + * @brief PixelDetector representation derived from Detector interface in the reconstruction chain + * + * Contains the PixelDetector with all its properties such as position and orientation, pitch, spatial resolution + * etc. + */ + class PixelDetector : public Detector { + public: + /** + * Delete default constructor + */ + PixelDetector() = delete; + + /** + * Default destructor + */ + ~PixelDetector() = default; + + /** + * @brief Constructs a detector in the geometry + * @param config Configuration object describing the detector + */ + PixelDetector(const Configuration& config); + + /** + * @brief Set position and orientation from configuration file + */ + void SetPostionAndOrientation(const Configuration& config); + + /** + * @brief Update detector position in the world + * @param displacement Vector with three position coordinates + */ + void displacement(XYZPoint displacement) override { m_displacement = displacement; } + + /** + * @brief Get position in the world + * @return Global position in Cartesian coordinates + */ + XYZPoint displacement() const override { return m_displacement; } + + /** + * @brief Get orientation in the world + * @return Vector with three rotation angles + */ + XYZVector rotation() const override { return m_orientation; } + + /** + * @brief Update detector orientation in the world + * @param rotation Vector with three rotation angles + */ + void rotation(XYZVector rotation) override { m_orientation = rotation; } + + /** + * @brief Mark a detector channel as masked + * @param chX X coordinate of the pixel to be masked + * @param chY Y coordinate of the pixel to be masked + */ + void maskChannel(int chX, int chY) override; + + /** + * @brief Check if a detector channel is masked + * @param chX X coordinate of the pixel to check + * @param chY Y coordinate of the pixel to check + * @return Mask status of the pixel in question + */ + bool masked(int chX, int chY) const override; + + // Function to get global intercept with a track + PositionVector3D> getIntercept(const Track* track) const override; + // Function to get local intercept with a track + PositionVector3D> getLocalIntercept(const Track* track) const override; + + // Function to check if a track intercepts with a plane + bool hasIntercept(const Track* track, double pixelTolerance = 0.) const override; + + // Function to check if a track goes through/near a masked pixel + bool hitMasked(const Track* track, int tolerance = 0.) const override; + + // Functions to get row and column from local position + double getRow(PositionVector3D> localPosition) const override; + + double getColumn(PositionVector3D> localPosition) const override; + + // Function to get local position from column (x) and row (y) coordinates + PositionVector3D> getLocalPosition(double column, double row) const override; + + /** + * Transformation from local (sensor) coordinates to in-pixel coordinates + * @param column Column address ranging from int_column-0.5*pitch to int_column+0.5*pitch + * @param row Row address ranging from int_column-0.5*pitch to int_column+0.5*pitch + * @return Position within a single pixel cell, given in units of length + */ + XYVector inPixel(const double column, const double row) const override; + + /** + * Transformation from local (sensor) coordinates to in-pixel coordinates + * @param localPosition Local position on the sensor + * @return Position within a single pixel cell, given in units of length + */ + XYVector inPixel(PositionVector3D> localPosition) const override; + + /** + * @brief Check whether given track is within the detector's region-of-interest + * @param track The track to be checked + * @return Boolean indicating cluster affiliation with region-of-interest + */ + bool isWithinROI(const Track* track) const override; + + /** + * @brief Check whether given cluster is within the detector's region-of-interest + * @param cluster The cluster to be checked + * @return Boolean indicating cluster affiliation with region-of-interest + */ + bool isWithinROI(Cluster* cluster) const override; + + /** + * @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 + */ + XYVector getSize() const override; + + /** + * @brief Get pitch of a single pixel + * @return Pitch of a pixel + */ + XYVector getPitch() const override { return m_pitch; } + + /** + * @brief Get intrinsic spatial resolution of the detector + * @return Intrinsic spatial resolution in X and Y + */ + XYVector getSpatialResolution() const override { return m_spatial_resolution; } + + /* + * @brief Get number of pixels in x and y + * @return Number of two dimensional pixels + */ + ROOT::Math::DisplacementVector2D> nPixels() const override { return m_nPixels; } + + private: + // Initialize coordinate transformations + void initialise() override; + + // Build axis, for devices which are not auxiliary + // Different in Pixel/Strip Detector + void build_axes(const Configuration& config) override; + + // Config detector, for devices which are not auxiliary + // Different in Pixel/Strip Detector + void configure_detector(Configuration& config) const override; + + // Config position, orientation, mode of detector + // Different in Pixel/Strip Detector + void configure_pos_and_orientation(Configuration& config) const override; + + // Functions to set and check channel masking + void process_mask_file() override; + + // Seems to be used in other coordinate + inline static int isLeft(std::pair pt0, std::pair pt1, std::pair pt2); + static int winding_number(std::pair probe, std::vector> polygon); + + // For planar detector + XYVector m_pitch{}; + XYVector m_spatial_resolution{}; + ROOT::Math::DisplacementVector2D> m_nPixels{}; + std::vector> m_roi{}; + // Displacement and rotation in x,y,z + ROOT::Math::XYZPoint m_displacement; + ROOT::Math::XYZVector m_orientation; + std::string m_orientation_mode; + }; +} // namespace corryvreckan + +#endif // CORRYVRECKAN_DETECTOR_H -- GitLab From 9e41e169eabac838e2f811a85efdb0b7475cdcbe Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Thu, 18 Jun 2020 11:51:12 +0200 Subject: [PATCH 03/37] semantic issues --- src/core/CMakeLists.txt | 1 + src/core/detector/BentPixelDetector.cpp | 48 ++++++++++++------------- src/core/detector/BentPixelDetector.hpp | 8 ++--- src/core/detector/Detector.cpp | 7 ++-- src/core/detector/Detector.hpp | 1 + 5 files changed, 35 insertions(+), 30 deletions(-) diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 96815fae..61edb3a8 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -25,6 +25,7 @@ ADD_LIBRARY(CorryvreckanCore SHARED Corryvreckan.cpp detector/Detector.cpp detector/PixelDetector.cpp + detector/BentPixelDetector.cpp detector/exceptions.cpp clipboard/Clipboard.cpp config/ConfigManager.cpp diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 8874b1e0..ec77dde6 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -15,14 +15,14 @@ #include "Math/RotationZ.h" #include "Math/RotationZYX.h" -#include "PixelDetector.hpp" +#include "BentPixelDetector.hpp" #include "core/utils/log.h" #include "exceptions.h" using namespace ROOT::Math; using namespace corryvreckan; -PixelDetector::PixelDetector(const Configuration& config) : Detector(config) { +BentPixelDetector::BentPixelDetector(const Configuration& config) : Detector(config) { // Set detector position and direction from configuration file SetPostionAndOrientation(config); @@ -36,7 +36,7 @@ PixelDetector::PixelDetector(const Configuration& config) : Detector(config) { } } -void PixelDetector::build_axes(const Configuration& config) { +void BentPixelDetector::build_axes(const Configuration& config) { m_nPixels = config.get>>("number_of_pixels"); // Size of the pixels: @@ -70,7 +70,7 @@ void PixelDetector::build_axes(const Configuration& config) { } } -void PixelDetector::SetPostionAndOrientation(const Configuration& config) { +void BentPixelDetector::SetPostionAndOrientation(const Configuration& config) { // Detector position and orientation m_displacement = config.get("position", ROOT::Math::XYZPoint()); m_orientation = config.get("orientation", ROOT::Math::XYZVector()); @@ -82,7 +82,7 @@ void PixelDetector::SetPostionAndOrientation(const Configuration& config) { LOG(TRACE) << " Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")"; } -void PixelDetector::process_mask_file() { +void BentPixelDetector::process_mask_file() { // Open the file with masked pixels std::ifstream inputMaskFile(m_maskfile, std::ios::in); if(!inputMaskFile.is_open()) { @@ -127,12 +127,12 @@ void PixelDetector::process_mask_file() { } } -void PixelDetector::maskChannel(int chX, int chY) { +void BentPixelDetector::maskChannel(int chX, int chY) { int channelID = chX + m_nPixels.X() * chY; m_masked[channelID] = true; } -bool PixelDetector::masked(int chX, int chY) const { +bool BentPixelDetector::masked(int chX, int chY) const { int channelID = chX + m_nPixels.X() * chY; if(m_masked.count(channelID) > 0) return true; @@ -140,7 +140,7 @@ bool PixelDetector::masked(int chX, int chY) const { } // Function to initialise transforms -void PixelDetector::initialise() { +void BentPixelDetector::initialise() { // Make the local to global transform, built from a displacement and // rotation @@ -180,7 +180,7 @@ void PixelDetector::initialise() { } // Only if detector is not auxiliary -void PixelDetector::configure_detector(Configuration& config) const { +void BentPixelDetector::configure_detector(Configuration& config) const { // Number of pixels config.set("number_of_pixels", m_nPixels); @@ -200,14 +200,14 @@ void PixelDetector::configure_detector(Configuration& config) const { config.setMatrix("roi", m_roi); } -void PixelDetector::configure_pos_and_orientation(Configuration& config) const { +void BentPixelDetector::configure_pos_and_orientation(Configuration& config) const { config.set("position", m_displacement, {"um", "mm"}); config.set("orientation_mode", m_orientation_mode); config.set("orientation", m_orientation, {"deg"}); } // Function to get global intercept with a track -PositionVector3D> PixelDetector::getIntercept(const Track* track) const { +PositionVector3D> BentPixelDetector::getIntercept(const Track* track) const { // FIXME: this is else statement can only be temporary if(track->getType() == "GblTrack") { @@ -230,12 +230,12 @@ PositionVector3D> PixelDetector::getIntercept(const Track* t } } -PositionVector3D> PixelDetector::getLocalIntercept(const Track* track) const { +PositionVector3D> BentPixelDetector::getLocalIntercept(const Track* track) const { return globalToLocal(getIntercept(track)); } // Function to check if a track intercepts with a plane -bool PixelDetector::hasIntercept(const Track* track, double pixelTolerance) const { +bool BentPixelDetector::hasIntercept(const Track* track, double pixelTolerance) const { // First, get the track intercept in global coordinates with the plane PositionVector3D> globalIntercept = this->getIntercept(track); @@ -258,7 +258,7 @@ bool PixelDetector::hasIntercept(const Track* track, double pixelTolerance) cons } // Function to check if a track goes through/near a masked pixel -bool PixelDetector::hitMasked(const Track* track, int tolerance) const { +bool BentPixelDetector::hitMasked(const Track* track, int tolerance) const { // First, get the track intercept in global coordinates with the plane PositionVector3D> globalIntercept = this->getIntercept(track); @@ -283,15 +283,15 @@ bool PixelDetector::hitMasked(const Track* track, int tolerance) const { } // Functions to get row and column from local position -double PixelDetector::getRow(const PositionVector3D> localPosition) const { +double BentPixelDetector::getRow(const PositionVector3D> localPosition) const { return localPosition.Y() / m_pitch.Y() + static_cast(m_nPixels.Y() - 1) / 2.; } -double PixelDetector::getColumn(const PositionVector3D> localPosition) const { +double BentPixelDetector::getColumn(const PositionVector3D> localPosition) const { return localPosition.X() / m_pitch.X() + static_cast(m_nPixels.X() - 1) / 2.; } // Function to get local position from row and column -PositionVector3D> PixelDetector::getLocalPosition(double column, double row) const { +PositionVector3D> BentPixelDetector::getLocalPosition(double column, double row) const { return PositionVector3D>(m_pitch.X() * (column - static_cast(m_nPixels.X() - 1) / 2.), m_pitch.Y() * (row - static_cast(m_nPixels.Y() - 1) / 2.), @@ -299,19 +299,19 @@ PositionVector3D> PixelDetector::getLocalPosition(double col } // Function to get in-pixel position -ROOT::Math::XYVector PixelDetector::inPixel(const double column, const double row) const { +ROOT::Math::XYVector BentPixelDetector::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 PixelDetector::inPixel(const PositionVector3D> localPosition) const { +ROOT::Math::XYVector BentPixelDetector::inPixel(const PositionVector3D> localPosition) const { double column = getColumn(localPosition); double row = getRow(localPosition); return inPixel(column, row); } // Check if track position is within ROI: -bool PixelDetector::isWithinROI(const Track* track) const { +bool BentPixelDetector::isWithinROI(const Track* track) const { // Empty region of interest: if(m_roi.empty()) { @@ -330,7 +330,7 @@ bool PixelDetector::isWithinROI(const Track* track) const { } // Check if cluster is within ROI and/or touches ROI border: -bool PixelDetector::isWithinROI(Cluster* cluster) const { +bool BentPixelDetector::isWithinROI(Cluster* cluster) const { // Empty region of interest: if(m_roi.empty()) { @@ -346,7 +346,7 @@ bool PixelDetector::isWithinROI(Cluster* cluster) const { return true; } -XYVector PixelDetector::getSize() const { +XYVector BentPixelDetector::getSize() const { return XYVector(m_pitch.X() * m_nPixels.X(), m_pitch.Y() * m_nPixels.Y()); } @@ -356,7 +356,7 @@ XYVector PixelDetector::getSize() const { * 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 PixelDetector::winding_number(std::pair probe, std::vector> polygon) { +int BentPixelDetector::winding_number(std::pair probe, std::vector> polygon) { // Two points don't make an area if(polygon.size() < 3) { LOG(DEBUG) << "No ROI given."; @@ -406,6 +406,6 @@ int PixelDetector::winding_number(std::pair probe, std::vector pt0, std::pair pt1, std::pair pt2) { +int BentPixelDetector::isLeft(std::pair pt0, std::pair pt1, std::pair pt2) { return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second)); } diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index eab25bf1..aac07111 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -33,23 +33,23 @@ namespace corryvreckan { * Contains the PixelDetector with all its properties such as position and orientation, pitch, spatial resolution * etc. */ - class PixelDetector : public Detector { + class BentPixelDetector : public Detector { public: /** * Delete default constructor */ - PixelDetector() = delete; + BentPixelDetector() = delete; /** * Default destructor */ - ~PixelDetector() = default; + ~BentPixelDetector() = default; /** * @brief Constructs a detector in the geometry * @param config Configuration object describing the detector */ - PixelDetector(const Configuration& config); + BentPixelDetector(const Configuration& config); /** * @brief Set position and orientation from configuration file diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index 8d1edef6..0b79121c 100644 --- a/src/core/detector/Detector.cpp +++ b/src/core/detector/Detector.cpp @@ -84,9 +84,12 @@ std::shared_ptr corryvreckan::Detector::factory(const Configuration& c std::transform(coordinates.begin(), coordinates.end(), coordinates.begin(), ::tolower); if(coordinates == "cartesian") { return std::make_shared(config); - } else { - throw InvalidValueError(config, "coordinates", "Coordinates can only set to be cartesian now"); } + if(coordinates == "cartesian-bent") { + return std::make_shared(config); + } + // coordinates was not found, throw an error + throw InvalidValueError(config, "coordinates", "Coordinates can only set to be cartesian now"); } double Detector::getTimeResolution() const { diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index 3ad0fbf3..60735971 100644 --- a/src/core/detector/Detector.hpp +++ b/src/core/detector/Detector.hpp @@ -373,5 +373,6 @@ namespace corryvreckan { }; } // namespace corryvreckan +#include "BentPixelDetector.hpp" #include "PixelDetector.hpp" #endif // CORRYVRECKAN_DETECTOR_H -- GitLab From 02adf80b8a683c83961ac430f5d5fa16306d270c Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Thu, 18 Jun 2020 14:06:00 +0200 Subject: [PATCH 04/37] implemented globalToLocal transformation for bent chip --- src/core/detector/BentPixelDetector.cpp | 30 +++++++++++++++++++++++++ src/core/detector/BentPixelDetector.hpp | 21 +++++++++++++++-- src/core/detector/Detector.cpp | 2 ++ src/core/detector/Detector.hpp | 6 ++--- src/core/detector/PixelDetector.hpp | 1 + 5 files changed, 54 insertions(+), 6 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index ec77dde6..94fc6a63 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -41,6 +41,9 @@ void BentPixelDetector::build_axes(const Configuration& config) { m_nPixels = config.get>>("number_of_pixels"); // Size of the pixels: m_pitch = config.get("pixel_pitch"); + m_coordinates = config.get("coordinates"); + m_y_0 = config.get("y_0"); + m_radius = config.get("radius"); LOG(TRACE) << "Initialized \"" << m_detectorType << "\": " << m_nPixels.X() << "x" << m_nPixels.Y() << " px, pitch of " << Units::display(m_pitch, {"mm", "um"}); @@ -204,6 +207,33 @@ void BentPixelDetector::configure_pos_and_orientation(Configuration& config) con config.set("position", m_displacement, {"um", "mm"}); config.set("orientation_mode", m_orientation_mode); config.set("orientation", m_orientation, {"deg"}); + config.set("coordinates", m_coordinates); + config.set("y_0", m_y_0, {"mm"}); + config.set("radius", m_radius, {"mm"}); +} + +XYZPoint BentPixelDetector::globalToLocal(XYZPoint global) const { + XYZPoint local_transformed = m_globalToLocal * global; + if(local_transformed.y() > m_y_0) { + return local_transformed; + } else { + XYZPoint local(local_transformed.x(), + m_radius * asin(local_transformed.y() / m_radius) + m_y_0, + m_radius * acos(local_transformed.z() / m_radius - 1.0) + m_y_0); + return local; + } +} + +XYZPoint BentPixelDetector::localToGlobal(XYZPoint local) const { + // transformation: cylinder starting below y_0, tangent to the flat part + if(local.y() > m_y_0) { + return m_localToGlobal * local; + } else { + XYZPoint local_transformed(local.x(), + m_radius * sin((local.y() - m_y_0) / m_radius), + m_radius * (cos((local.z() - m_y_0) / m_radius) - 1.0)); + return m_localToGlobal * local_transformed; + } } // Function to get global intercept with a track diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index aac07111..d89c6603 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -6,8 +6,8 @@ * Intergovernmental Organization or submit itself to any jurisdiction. */ -#ifndef CORRYVRECKAN_PLANARDETECTOR_H -#define CORRYVRECKAN_PLANARDETECTOR_H +#ifndef CORRYVRECKAN_BENTPIXELDETECTOR_H +#define CORRYVRECKAN_BENTPIXELDETECTOR_H #include #include @@ -114,6 +114,20 @@ namespace corryvreckan { // Function to get local position from column (x) and row (y) coordinates PositionVector3D> getLocalPosition(double column, double row) const override; + /** + * @brief Transform local coordinates of this detector into global coordinates + * @param local Local coordinates in the reference frame of this detector + * @return Global coordinates + */ + XYZPoint localToGlobal(XYZPoint local) const override; + + /** + * @brief Transform global coordinates into detector-local coordinates + * @param global Global coordinates + * @return Local coordinates in the reference frame of this detector + */ + XYZPoint globalToLocal(XYZPoint global) const override; + /** * Transformation from local (sensor) coordinates to in-pixel coordinates * @param column Column address ranging from int_column-0.5*pitch to int_column+0.5*pitch @@ -192,6 +206,8 @@ namespace corryvreckan { // For planar detector XYVector m_pitch{}; + double m_radius; + double m_y_0; XYVector m_spatial_resolution{}; ROOT::Math::DisplacementVector2D> m_nPixels{}; std::vector> m_roi{}; @@ -199,6 +215,7 @@ namespace corryvreckan { ROOT::Math::XYZPoint m_displacement; ROOT::Math::XYZVector m_orientation; std::string m_orientation_mode; + std::string m_coordinates; }; } // namespace corryvreckan diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index 0b79121c..8c217479 100644 --- a/src/core/detector/Detector.cpp +++ b/src/core/detector/Detector.cpp @@ -15,7 +15,9 @@ #include "Math/RotationZ.h" #include "Math/RotationZYX.h" +#include "BentPixelDetector.hpp" #include "Detector.hpp" +#include "PixelDetector.hpp" #include "core/utils/log.h" #include "exceptions.h" diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index 60735971..eb1dea22 100644 --- a/src/core/detector/Detector.hpp +++ b/src/core/detector/Detector.hpp @@ -280,14 +280,14 @@ namespace corryvreckan { * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; + virtual XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; + virtual XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; /** * @brief Check whether given track is within the detector's region-of-interest @@ -373,6 +373,4 @@ namespace corryvreckan { }; } // namespace corryvreckan -#include "BentPixelDetector.hpp" -#include "PixelDetector.hpp" #endif // CORRYVRECKAN_DETECTOR_H diff --git a/src/core/detector/PixelDetector.hpp b/src/core/detector/PixelDetector.hpp index 211525db..74433909 100644 --- a/src/core/detector/PixelDetector.hpp +++ b/src/core/detector/PixelDetector.hpp @@ -210,6 +210,7 @@ namespace corryvreckan { ROOT::Math::XYZPoint m_displacement; ROOT::Math::XYZVector m_orientation; std::string m_orientation_mode; + std::string m_coordinates; }; } // namespace corryvreckan -- GitLab From 29797d21705bdcda67d9fde742e0be23523726ae Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Fri, 19 Jun 2020 09:01:46 +0200 Subject: [PATCH 05/37] implemented corrected transformation, still buggy --- src/core/detector/BentPixelDetector.cpp | 26 +++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 94fc6a63..5c9ea449 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -214,24 +214,34 @@ void BentPixelDetector::configure_pos_and_orientation(Configuration& config) con XYZPoint BentPixelDetector::globalToLocal(XYZPoint global) const { XYZPoint local_transformed = m_globalToLocal * global; - if(local_transformed.y() > m_y_0) { - return local_transformed; + if(local_transformed.y() > 0) { + double row = m_nPixels.Y() - (m_y_0 - local_transformed.y()) / m_pitch.Y(); + XYZPoint local(local_transformed.x(), + m_pitch.Y() * (row - static_cast(m_nPixels.Y() - 1) / 2.), // y_0 = zero + 0); + return local; } else { + double y_distance_from_top = m_radius * asin(-local_transformed.y() / m_radius) + m_y_0; + double row = m_nPixels.Y() - y_distance_from_top / m_pitch.Y(); XYZPoint local(local_transformed.x(), - m_radius * asin(local_transformed.y() / m_radius) + m_y_0, - m_radius * acos(local_transformed.z() / m_radius - 1.0) + m_y_0); + m_pitch.Y() * (row - static_cast(m_nPixels.Y() - 1) / 2.), // y_0 = zero + 0); return local; } } XYZPoint BentPixelDetector::localToGlobal(XYZPoint local) const { // transformation: cylinder starting below y_0, tangent to the flat part - if(local.y() > m_y_0) { - return m_localToGlobal * local; + double y_distance_from_top = (m_nPixels.Y() - getRow(local)) * m_pitch.Y(); + if(y_distance_from_top < m_y_0) { + XYZPoint local_transformed(local.x(), + -y_distance_from_top + m_y_0, // y_0 = zero + local.z()); + return m_localToGlobal * local_transformed; } else { XYZPoint local_transformed(local.x(), - m_radius * sin((local.y() - m_y_0) / m_radius), - m_radius * (cos((local.z() - m_y_0) / m_radius) - 1.0)); + -m_radius * sin((y_distance_from_top - m_y_0) / m_radius), + m_radius * (cos((y_distance_from_top - m_y_0) / m_radius) - 1.0)); return m_localToGlobal * local_transformed; } } -- GitLab From c0fb0f5a1d1a43757ca5423ad8fbe6b8f807a6bf Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Fri, 19 Jun 2020 16:11:42 +0200 Subject: [PATCH 06/37] implemented getIntercept with corrected transformation --- src/core/detector/BentPixelDetector.cpp | 45 ++++++++++++++++++++++--- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 5c9ea449..70997c49 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -241,7 +241,7 @@ XYZPoint BentPixelDetector::localToGlobal(XYZPoint local) const { } else { XYZPoint local_transformed(local.x(), -m_radius * sin((y_distance_from_top - m_y_0) / m_radius), - m_radius * (cos((y_distance_from_top - m_y_0) / m_radius) - 1.0)); + m_radius * (cos((y_distance_from_top - m_y_0) / m_radius) - 1.0) + local.z()); return m_localToGlobal * local_transformed; } } @@ -262,11 +262,46 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac track->getDirection(m_detectorName).Z() * m_normal.Z()); // Propagate the track - PositionVector3D> globalIntercept( + PositionVector3D> globalPlanarIntercept( track->getState(m_detectorName).X() + distance * track->getDirection(m_detectorName).X(), track->getState(m_detectorName).Y() + distance * track->getDirection(m_detectorName).Y(), track->getState(m_detectorName).Z() + distance * track->getDirection(m_detectorName).Z()); - return globalIntercept; + + if((m_nPixels.Y() - getRow(globalToLocal(globalPlanarIntercept))) * m_pitch.Y() < m_y_0) { + return globalPlanarIntercept; + } else { // intersect with the curvature, correct + XYZPoint state(track->getState(m_detectorName).X(), + track->getState(m_detectorName).Y(), + track->getState(m_detectorName).Z()); + XYZPoint direction(track->getDirection(m_detectorName).X(), + track->getDirection(m_detectorName).Y(), + track->getDirection(m_detectorName).Z()); + XYZPoint direction_tr = m_globalToLocal * direction; + XYZPoint state_tr = m_globalToLocal * state; + // solve quadratic equation for intersect: alpha*z^2+beta*z+gamma=0 + double gamma = + state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + + 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() - m_radius * m_radius; + double beta = 2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + + 2 * state_tr.y() * direction_tr.y() / direction_tr.z(); + double alpha = 1 + direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()); + + // check if the equation has a solution + double tmp = beta * beta - 4 * alpha * gamma; + if(tmp < 0) // return a value outside the chip, none is at the moment not possible + return (localToGlobal(getLocalPosition(2 * m_nPixels.X(), 2 * m_nPixels.Y()))); + // change the intercept point with the cylinder with VZ + double z_intercept = -beta + sqrt(tmp) / (2 * alpha); + // calculate x and y from z-intercept + XYZPoint intercept_tr((z_intercept - state_tr.z()) * direction_tr.x() / direction_tr.z() + state_tr.x(), + (z_intercept - state_tr.z()) * direction_tr.y() / direction_tr.z() + state_tr.y(), + z_intercept); + XYZPoint globalIntercept = m_localToGlobal * intercept_tr; + + PositionVector3D> globalIntercept_cart( + globalIntercept.x(), globalIntercept.y(), globalIntercept.z()); + return globalIntercept_cart; + } } } @@ -281,7 +316,7 @@ bool BentPixelDetector::hasIntercept(const Track* track, double pixelTolerance) PositionVector3D> globalIntercept = this->getIntercept(track); // Convert to local coordinates - PositionVector3D> localIntercept = this->m_globalToLocal * globalIntercept; + PositionVector3D> localIntercept = globalToLocal(globalIntercept); // Get the row and column numbers double row = this->getRow(localIntercept); @@ -304,7 +339,7 @@ bool BentPixelDetector::hitMasked(const Track* track, int tolerance) const { PositionVector3D> globalIntercept = this->getIntercept(track); // Convert to local coordinates - PositionVector3D> localIntercept = this->m_globalToLocal * globalIntercept; + PositionVector3D> localIntercept = globalToLocal(globalIntercept); // Get the row and column numbers int row = static_cast(floor(this->getRow(localIntercept) + 0.5)); -- GitLab From f16563c6896206e05ae9781e2580efc42308245a Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Fri, 19 Jun 2020 17:35:38 +0200 Subject: [PATCH 07/37] fixed VZ --- src/core/detector/BentPixelDetector.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 70997c49..541954f1 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -291,7 +291,7 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac if(tmp < 0) // return a value outside the chip, none is at the moment not possible return (localToGlobal(getLocalPosition(2 * m_nPixels.X(), 2 * m_nPixels.Y()))); // change the intercept point with the cylinder with VZ - double z_intercept = -beta + sqrt(tmp) / (2 * alpha); + double z_intercept = (-beta - sqrt(tmp)) / (2 * alpha); // calculate x and y from z-intercept XYZPoint intercept_tr((z_intercept - state_tr.z()) * direction_tr.x() / direction_tr.z() + state_tr.x(), (z_intercept - state_tr.z()) * direction_tr.y() / direction_tr.z() + state_tr.y(), -- GitLab From 3977978bc7320b1cf3266379caee2fb9e45fdf61 Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Thu, 2 Jul 2020 10:36:01 +0200 Subject: [PATCH 08/37] fixed intersect bugs --- src/core/detector/BentPixelDetector.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 541954f1..69ccb479 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -276,14 +276,14 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac XYZPoint direction(track->getDirection(m_detectorName).X(), track->getDirection(m_detectorName).Y(), track->getDirection(m_detectorName).Z()); - XYZPoint direction_tr = m_globalToLocal * direction; + XYZPoint direction_tr = m_globalToLocal.Rotation() * direction; XYZPoint state_tr = m_globalToLocal * state; // solve quadratic equation for intersect: alpha*z^2+beta*z+gamma=0 double gamma = - state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + - 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() - m_radius * m_radius; - double beta = 2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + - 2 * state_tr.y() * direction_tr.y() / direction_tr.z(); + state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) - + 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() + state_tr.y() * state_tr.y(); + double beta = -2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + + 2 * state_tr.y() * direction_tr.y() / direction_tr.z() - 2 * m_radius; double alpha = 1 + direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()); // check if the equation has a solution -- GitLab From 410b76874ef6e28ea99a51cd148c3982891513f2 Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Thu, 9 Jul 2020 13:36:14 +0200 Subject: [PATCH 09/37] fixed intersect --- src/core/detector/BentPixelDetector.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 69ccb479..5f4aec3b 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -283,7 +283,7 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) - 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() + state_tr.y() * state_tr.y(); double beta = -2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + - 2 * state_tr.y() * direction_tr.y() / direction_tr.z() - 2 * m_radius; + 2 * state_tr.y() * direction_tr.y() / direction_tr.z() + 2 * m_radius; double alpha = 1 + direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()); // check if the equation has a solution @@ -291,7 +291,7 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac if(tmp < 0) // return a value outside the chip, none is at the moment not possible return (localToGlobal(getLocalPosition(2 * m_nPixels.X(), 2 * m_nPixels.Y()))); // change the intercept point with the cylinder with VZ - double z_intercept = (-beta - sqrt(tmp)) / (2 * alpha); + double z_intercept = (-beta + sqrt(tmp)) / (2 * alpha); // calculate x and y from z-intercept XYZPoint intercept_tr((z_intercept - state_tr.z()) * direction_tr.x() / direction_tr.z() + state_tr.x(), (z_intercept - state_tr.z()) * direction_tr.y() / direction_tr.z() + state_tr.y(), -- GitLab From f6b44a94ae37a17fb7b2f326038c0b1e8c894bf3 Mon Sep 17 00:00:00 2001 From: Alexander Ferk Date: Fri, 18 Sep 2020 11:16:23 +0200 Subject: [PATCH 10/37] added isNeighbor function --- src/core/detector/BentPixelDetector.cpp | 19 +++++++++++++++++++ src/core/detector/BentPixelDetector.hpp | 9 +++++++++ 2 files changed, 28 insertions(+) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 5f4aec3b..72c33744 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -421,6 +421,25 @@ bool BentPixelDetector::isWithinROI(Cluster* cluster) const { return true; } +// Check if a pixel touches any of the pixels in a cluster + +bool BentPixelDetector::isNeighbor(const std::shared_ptr& neighbor, + const std::shared_ptr& cluster, + const int neighbor_radius_row, + const int neighbor_radius_col) { + for(auto pixel : cluster->pixels()) { + int row_distance = abs(pixel->row() - neighbor->row()); + int col_distance = abs(pixel->column() - neighbor->column()); + if(row_distance <= neighbor_radius_row && col_distance <= neighbor_radius_col) { + if(row_distance > 1 || col_distance > 1) { + cluster->setSplit(true); + } + return true; + } + } + return false; +} + XYVector BentPixelDetector::getSize() const { return XYVector(m_pitch.X() * m_nPixels.X(), m_pitch.Y() * m_nPixels.Y()); } diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index d89c6603..3d55bea2 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -181,6 +181,15 @@ namespace corryvreckan { */ ROOT::Math::DisplacementVector2D> nPixels() const override { return m_nPixels; } + /** + * @brief Test whether one pixel touches the cluster + * @return true if it fulfills the condition + * @note users should define their specific clustering method in the detector class, for pixel detector, the default + * is 2D clustering + */ + virtual bool + isNeighbor(const std::shared_ptr&, const std::shared_ptr&, const int, const int) override; + private: // Initialize coordinate transformations void initialise() override; -- GitLab From 8a09eefbbff1b56ce63ddfdcd9ad517dddb9578b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 2 Mar 2021 15:54:32 +0100 Subject: [PATCH 11/37] Update bent geometry for august test beam setup - Updated localToGlobal and inverse transformations --- src/core/detector/BentPixelDetector.cpp | 32 ++++--------------------- 1 file changed, 5 insertions(+), 27 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 72c33744..3c64cc00 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -214,36 +214,14 @@ void BentPixelDetector::configure_pos_and_orientation(Configuration& config) con XYZPoint BentPixelDetector::globalToLocal(XYZPoint global) const { XYZPoint local_transformed = m_globalToLocal * global; - if(local_transformed.y() > 0) { - double row = m_nPixels.Y() - (m_y_0 - local_transformed.y()) / m_pitch.Y(); - XYZPoint local(local_transformed.x(), - m_pitch.Y() * (row - static_cast(m_nPixels.Y() - 1) / 2.), // y_0 = zero - 0); - return local; - } else { - double y_distance_from_top = m_radius * asin(-local_transformed.y() / m_radius) + m_y_0; - double row = m_nPixels.Y() - y_distance_from_top / m_pitch.Y(); - XYZPoint local(local_transformed.x(), - m_pitch.Y() * (row - static_cast(m_nPixels.Y() - 1) / 2.), // y_0 = zero - 0); - return local; - } + double xtilde = asin(local_transformed.X() / m_radius) * m_radius; + XYZPoint p{xtilde, local_transformed.Y(), local_transformed.Z() + m_radius * cos(xtilde / m_radius) - m_radius}; + return p; } XYZPoint BentPixelDetector::localToGlobal(XYZPoint local) const { - // transformation: cylinder starting below y_0, tangent to the flat part - double y_distance_from_top = (m_nPixels.Y() - getRow(local)) * m_pitch.Y(); - if(y_distance_from_top < m_y_0) { - XYZPoint local_transformed(local.x(), - -y_distance_from_top + m_y_0, // y_0 = zero - local.z()); - return m_localToGlobal * local_transformed; - } else { - XYZPoint local_transformed(local.x(), - -m_radius * sin((y_distance_from_top - m_y_0) / m_radius), - m_radius * (cos((y_distance_from_top - m_y_0) / m_radius) - 1.0) + local.z()); - return m_localToGlobal * local_transformed; - } + XYZPoint p{m_radius * sin(local.X() / m_radius), local.Y(), -m_radius * (cos(local.X() / m_radius) - 1.) + local.z()}; + return m_localToGlobal * p; } // Function to get global intercept with a track -- GitLab From afccf42762bacca96f797a01f08ca3bbfae26179 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 2 Mar 2021 15:55:35 +0100 Subject: [PATCH 12/37] Add residuals X vs residuals Y correlation plots in AnalysisDUT --- src/modules/AnalysisDUT/AnalysisDUT.cpp | 50 +++++++++++++++++++++++++ src/modules/AnalysisDUT/AnalysisDUT.h | 8 +++- 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index c559afe0..d99556d7 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -260,6 +260,49 @@ void AnalysisDUT::initialize() { 4000, -500.5, 499.5); + residualsXvsXhit = new TH2F("residualsXvsXhit", + "Residual in X vs Xhit;x_{track}-x_{hit} [#mum];x_{hit} [mm];# entries", + 4000, + -500.5, + 499.5, + 40, + -20, + 20); + residualsXvsYhit = new TH2F("residualsXvsYhit", + "Residual in X vs Yhit;x_{track}-x_{hit} [#mum];y_{hit} [mm];# entries", + 4000, + -500.5, + 499.5, + 40, + -20, + 20); + residualsYvsXhit = new TH2F("residualsYvsXhit", + "Residual in Y vs Xhit;y_{track}-y_{hit} [#mum];x_{hit} [mm];# entries", + 4000, + -500.5, + 499.5, + 40, + -20, + 20); + residualsYvsYhit = new TH2F("residualsYvsYhit", + "Residual in Y vs Yhit;y_{track}-y_{hit} [#mum];y_{hit} [mm];# entries", + 4000, + -500.5, + 499.5, + 40, + -20, + 20); + residualsXvsresidualsY = + new TH2F("residualsXvsresidualsY", + "Residual in X vs Residual in Y;x_{track}-x_{hit} [#mum];y_{track}-y_{hit} [#mum];# entries", + 4000, + -500.5, + 499.5, + 4000, + -500.5, + 499.5); + residualsXprofile = new TProfile2D("residualsXprofile", "Residual in X;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); + residualsYprofile = new TProfile2D("residualsYprofile", "Residual in Y;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); clusterChargeAssoc = new TH1F("clusterChargeAssociated", "Charge distribution of associated clusters;cluster charge [e];# entries", @@ -752,6 +795,13 @@ StatusCode AnalysisDUT::run(const std::shared_ptr& clipboard) { residualsY->Fill(ydistance_um); residualsPos->Fill(posDiff_um); residualsPosVsresidualsTime->Fill(tdistance, posDiff_um); + residualsXvsXhit->Fill(xdistance_um, assoc_cluster->global().x()); + residualsXvsYhit->Fill(xdistance_um, assoc_cluster->global().y()); + residualsYvsXhit->Fill(ydistance_um, assoc_cluster->global().x()); + residualsYvsYhit->Fill(ydistance_um, assoc_cluster->global().y()); + residualsXvsresidualsY->Fill(xdistance_um, ydistance_um); + residualsXprofile->Fill(assoc_cluster->column(), assoc_cluster->row(), xdistance_um); + residualsYprofile->Fill(assoc_cluster->column(), assoc_cluster->row(), ydistance_um); if(assoc_cluster->columnWidth() == 1) { residualsX1pix->Fill(xdistance_um); diff --git a/src/modules/AnalysisDUT/AnalysisDUT.h b/src/modules/AnalysisDUT/AnalysisDUT.h index f2202b5d..bb3479af 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -54,7 +54,13 @@ namespace corryvreckan { TProfile2D* hPixelRawValueMapAssoc; TH1F* associatedTracksVersusTime; - TH1F *residualsX, *residualsY, *residualsPos; + TH1F *residualsX, *residualsY, *residualsZ, *residualsPos; + TH2F *residualsXvsXhit, *residualsXvsYhit; + TH2F *residualsYvsXhit, *residualsYvsYhit; + TH2F *residualsXvsresidualsY; + TProfile2D* residualsXprofile; + TProfile2D* residualsYprofile; + TH2F* residualsPosVsresidualsTime; TH1F *residualsX1pix, *residualsY1pix; -- GitLab From 808d0db816ef7cd8d78d1eddbeebc29e5294e12a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 2 Mar 2021 22:38:09 +0100 Subject: [PATCH 13/37] Format additions to AnalysisDUT --- src/modules/AnalysisDUT/AnalysisDUT.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.h b/src/modules/AnalysisDUT/AnalysisDUT.h index bb3479af..d179bda9 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -57,7 +57,7 @@ namespace corryvreckan { TH1F *residualsX, *residualsY, *residualsZ, *residualsPos; TH2F *residualsXvsXhit, *residualsXvsYhit; TH2F *residualsYvsXhit, *residualsYvsYhit; - TH2F *residualsXvsresidualsY; + TH2F* residualsXvsresidualsY; TProfile2D* residualsXprofile; TProfile2D* residualsYprofile; -- GitLab From 802902f94b9bab333ae1906bb27cc06939d81f47 Mon Sep 17 00:00:00 2001 From: Mihail Bogdan Blidaru Date: Wed, 3 Mar 2021 10:06:13 +0100 Subject: [PATCH 14/37] Modified getIntercept function to match August geometry. --- src/core/detector/BentPixelDetector.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 3c64cc00..c9748fe3 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -245,7 +245,7 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac track->getState(m_detectorName).Y() + distance * track->getDirection(m_detectorName).Y(), track->getState(m_detectorName).Z() + distance * track->getDirection(m_detectorName).Z()); - if((m_nPixels.Y() - getRow(globalToLocal(globalPlanarIntercept))) * m_pitch.Y() < m_y_0) { + if(m_radius == 0) { return globalPlanarIntercept; } else { // intersect with the curvature, correct XYZPoint state(track->getState(m_detectorName).X(), @@ -257,19 +257,21 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac XYZPoint direction_tr = m_globalToLocal.Rotation() * direction; XYZPoint state_tr = m_globalToLocal * state; // solve quadratic equation for intersect: alpha*z^2+beta*z+gamma=0 + double alpha = 1 + direction_tr.x() * direction_tr.x() / (direction_tr.z() * direction_tr.z()); + double beta = 2 * state_tr.x() * direction_tr.x() / direction_tr.z() - + 2 * state_tr.z() * direction_tr.x() * direction_tr.x() / (direction_tr.z() * direction_tr.z()) - + 2 * m_radius; double gamma = - state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) - - 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() + state_tr.y() * state_tr.y(); - double beta = -2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + - 2 * state_tr.y() * direction_tr.y() / direction_tr.z() + 2 * m_radius; - double alpha = 1 + direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()); + state_tr.x() * state_tr.x() + + state_tr.z() * state_tr.z() * direction_tr.x() * direction_tr.x() / (direction_tr.z() * direction_tr.z()) - + 2 * state_tr.x() * state_tr.z() * direction_tr.x() / direction_tr.z(); // check if the equation has a solution - double tmp = beta * beta - 4 * alpha * gamma; - if(tmp < 0) // return a value outside the chip, none is at the moment not possible + double delta = beta * beta - 4 * alpha * gamma; + if(delta < 0) // return a value outside the chip, none is at the moment not possible return (localToGlobal(getLocalPosition(2 * m_nPixels.X(), 2 * m_nPixels.Y()))); // change the intercept point with the cylinder with VZ - double z_intercept = (-beta + sqrt(tmp)) / (2 * alpha); + double z_intercept = (-beta - sqrt(delta)) / (2 * alpha); // calculate x and y from z-intercept XYZPoint intercept_tr((z_intercept - state_tr.z()) * direction_tr.x() / direction_tr.z() + state_tr.x(), (z_intercept - state_tr.z()) * direction_tr.y() / direction_tr.z() + state_tr.y(), -- GitLab From e930813304fda974348e7e95dc3574f64bec82a5 Mon Sep 17 00:00:00 2001 From: p-becht <> Date: Wed, 3 Mar 2021 10:52:03 +0100 Subject: [PATCH 15/37] copied from my personal fork, shearrz constraint in Alignment Millipede applied, general implementation of bending direction and Intercept that even should allow for an angled track --- src/core/detector/BentPixelDetector.cpp | 336 ++++++++++++++---- src/core/detector/BentPixelDetector.hpp | 29 +- src/core/detector/Detector.cpp | 9 +- src/core/detector/Detector.hpp | 5 +- src/core/detector/PixelDetector.hpp | 1 - .../AlignmentMillepede/AlignmentMillepede.cpp | 45 ++- src/modules/AnalysisDUT/AnalysisDUT.h | 8 +- .../EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp | 7 +- src/objects/Event.hpp | 7 +- src/objects/exceptions.h | 14 - 10 files changed, 321 insertions(+), 140 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index c9748fe3..4c8c5097 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -10,6 +10,8 @@ #include #include +#include "Math/DisplacementVector2D.h" +#include "Math/PositionVector3D.h" #include "Math/RotationX.h" #include "Math/RotationY.h" #include "Math/RotationZ.h" @@ -41,12 +43,27 @@ void BentPixelDetector::build_axes(const Configuration& config) { m_nPixels = config.get>>("number_of_pixels"); // Size of the pixels: m_pitch = config.get("pixel_pitch"); + // Geometry related parameters m_coordinates = config.get("coordinates"); - m_y_0 = config.get("y_0"); - m_radius = config.get("radius"); + m_flat_part = config.get("flat_part", 0); + m_radius = config.get("radius", 0); + m_bent_axis = config.get("bent_axis", "column"); + std::transform(m_bent_axis.begin(), m_bent_axis.end(), m_bent_axis.begin(), ::tolower); + if(m_bent_axis == "column" || m_bent_axis == "col") { + m_bent_axis = "column"; + } else if(m_bent_axis == "row") { + m_bent_axis = "row"; + } else { + throw InvalidValueError(config, "bent_axis", "Bent axis can only be set to \"column\" or \"row\" for now"); + } LOG(TRACE) << "Initialized \"" << m_detectorType << "\": " << m_nPixels.X() << "x" << m_nPixels.Y() << " px, pitch of " << Units::display(m_pitch, {"mm", "um"}); + LOG(TRACE) << "Initialized \"" + << "BentPixelDetector" + << "\": " << m_bent_axis << " is bent axis, " + << "bending radius of " << Units::display(m_radius, {"mm"}) << ", flat part of " + << Units::display(m_flat_part, {"mm"}); 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) { @@ -68,7 +85,7 @@ void BentPixelDetector::build_axes(const Configuration& config) { m_maskfile_name = config.get("mask_file"); std::string mask_file = config.getPath("mask_file", true); LOG(DEBUG) << "Adding mask to detector \"" << config.getName() << "\", reading from " << mask_file; - set_mask_file(mask_file); + maskFile(mask_file); process_mask_file(); } } @@ -78,11 +95,16 @@ void BentPixelDetector::SetPostionAndOrientation(const Configuration& config) { m_displacement = config.get("position", ROOT::Math::XYZPoint()); m_orientation = config.get("orientation", ROOT::Math::XYZVector()); m_orientation_mode = config.get("orientation_mode", "xyz"); - if(m_orientation_mode != "xyz" && m_orientation_mode != "zyx") { - throw InvalidValueError(config, "orientation_mode", "Invalid detector orientation mode"); + + if(m_orientation_mode != "xyz" && m_orientation_mode != "zyx" && m_orientation_mode != "zxz") { + throw InvalidValueError(config, "orientation_mode", "orientation_mode should be either 'zyx', xyz' or 'zxz'"); } + + LOG(WARNING) << "BentPixelDetector!"; LOG(TRACE) << " Position: " << Units::display(m_displacement, {"mm", "um"}); LOG(TRACE) << " Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")"; + LOG(TRACE) << " m_displacement: X=" << m_displacement.X() << " Y=" << m_displacement.Y() << " Z=" << m_displacement.Z(); + LOG(TRACE) << " m_orientation: X=" << m_orientation.X() << " Y=" << m_orientation.Y() << " Z=" << m_orientation.Z(); } void BentPixelDetector::process_mask_file() { @@ -163,10 +185,7 @@ void BentPixelDetector::initialise() { // First angle given in the configuration file is around z, second around x, last around z: rotations = EulerAngles(m_orientation.x(), m_orientation.y(), m_orientation.z()); } else { - LOG(ERROR) << "orientation_mode should be either 'zyx', xyz' or 'zxz'"; - // FIXME: To through exception, initialise needs to be changed to "initialise(const Configuration& config)" - // throw InvalidValueError( - // config, "orientation_mode", "orientation_mode should be either 'zyx', xyz' or 'zxz'"); + throw InvalidSettingError(this, "orientation_mode", "orientation_mode should be either 'zyx', xyz' or 'zxz'"); } m_localToGlobal = Transform3D(rotations, translations); @@ -208,26 +227,109 @@ void BentPixelDetector::configure_pos_and_orientation(Configuration& config) con config.set("orientation_mode", m_orientation_mode); config.set("orientation", m_orientation, {"deg"}); config.set("coordinates", m_coordinates); - config.set("y_0", m_y_0, {"mm"}); config.set("radius", m_radius, {"mm"}); + config.set("flat_part", m_flat_part, {"mm"}); + config.set("bent_axis", m_bent_axis); } -XYZPoint BentPixelDetector::globalToLocal(XYZPoint global) const { - XYZPoint local_transformed = m_globalToLocal * global; - double xtilde = asin(local_transformed.X() / m_radius) * m_radius; - XYZPoint p{xtilde, local_transformed.Y(), local_transformed.Z() + m_radius * cos(xtilde / m_radius) - m_radius}; - return p; +XYZPoint BentPixelDetector::localToGlobal(XYZPoint local) const { + ROOT::Math::XYVector detector_dimensions = getSize(); + // which axis of the sensor is bent + if(m_bent_axis == "column") { + if(m_flat_part != 0 || m_radius == 0) { + LOG(WARNING) << "Column bending does not allow for a flat part at the moment or zero radius. Return cartesian " + "coordinates"; + return m_localToGlobal * local; + } + // origin of coordinate system in the middle of the flat chip (tangential point) + double col_arc_length = local.X(); + ROOT::Math::XYZPoint local_transformed( + m_radius * + sin(col_arc_length / m_radius), // >0 for col_arc_length>0 and vice versa independent of sign(m_radius) + local.y(), + local.z() - m_radius * (1.0 - cos(col_arc_length / m_radius))); // local.z() + ... for r<0 + local_transformed = m_localToGlobal * local_transformed; + LOG(TRACE) << "Transformed local point (" << local.x() << "|" << local.y() << "|" << local.z() + << ") to global point (" << local_transformed.x() << "|" << local_transformed.y() << "|" + << local_transformed.z() << "). Bent column"; + return local_transformed; + } else if(m_bent_axis == "row") { + if(m_radius == 0) { + LOG(WARNING) << "Zero radius. Return cartesian coordinates"; + return m_localToGlobal * local; + } + // origin of coordinate system in the middle of the flat chip + double row_arc_length = detector_dimensions.Y() / 2 - local.Y(); // distance from top of the chip + if(row_arc_length < m_flat_part) { // flat sensor + ROOT::Math::XYZPoint local_transformed(local.x(), -row_arc_length + detector_dimensions.Y() / 2, local.z()); + local_transformed = m_localToGlobal * local_transformed; + LOG(TRACE) << "Transformed local point (" << local.x() << "|" << local.y() << "|" << local.z() + << ") to global point (" << local_transformed.x() << "|" << local_transformed.y() << "|" + << local_transformed.z() << "). Bent row"; + return local_transformed; + } else { // bent sensor + ROOT::Math::XYZPoint local_transformed( + local.x(), + -m_radius * sin((row_arc_length - m_flat_part) / m_radius) + detector_dimensions.Y() / 2 - + m_flat_part, // works for negative radius + m_radius * (cos((row_arc_length - m_flat_part) / m_radius) - 1.0) + local.z()); // works for negative radius + local_transformed = m_localToGlobal * local_transformed; + LOG(TRACE) << "Transformed local point (" << local.x() << "|" << local.y() << "|" << local.z() + << ") to global point (" << local_transformed.x() << "|" << local_transformed.y() << "|" + << local_transformed.z() << "). Bent row"; + return local_transformed; + } + } else { + LOG(ERROR) << "m_bent_axis = " << m_bent_axis << " not implemented. Return cartesian coordinates"; + return m_localToGlobal * local; + } } -XYZPoint BentPixelDetector::localToGlobal(XYZPoint local) const { - XYZPoint p{m_radius * sin(local.X() / m_radius), local.Y(), -m_radius * (cos(local.X() / m_radius) - 1.) + local.z()}; - return m_localToGlobal * p; +XYZPoint BentPixelDetector::globalToLocal(XYZPoint global) const { + ROOT::Math::XYVector detector_dimensions = getSize(); + // transform from global to bent local + ROOT::Math::XYZPoint local_transformed = m_globalToLocal * global; + // which axis of the sensor is bent + if(m_bent_axis == "column") { + if(m_flat_part != 0 || m_radius == 0) { + LOG(WARNING) << "Column bending does not allow for a flat part at the moment. Return cartesian coordinates"; + return m_globalToLocal * global; + } + // origin of coordinate system in the middle of the chip (tangential point) + double col_arc_length = m_radius * asin(local_transformed.x() / m_radius); + ROOT::Math::XYZPoint local(col_arc_length, local_transformed.y(), 0); + LOG(TRACE) << "Transformed global point (" << global.x() << "|" << global.y() << "|" << global.z() + << ") to local point (" << local.x() << "|" << local.y() << "|" << local.z() << "). Bent column"; + return local; + } else if(m_bent_axis == "row") { + if(m_radius == 0) { + LOG(WARNING) << "Zero radius. Return cartesian coordinates"; + return m_globalToLocal * global; + } + double row_arc_length = + m_radius * asin((detector_dimensions.Y() / 2 - m_flat_part - local_transformed.Y()) / m_radius) + m_flat_part; + // origin of coordinate in the middle of the flat chip + if(row_arc_length < m_flat_part) { // flat part + ROOT::Math::XYZPoint local(local_transformed.x(), local_transformed.y(), 0); + LOG(TRACE) << "Transformed global point (" << global.x() << "|" << global.y() << "|" << global.z() + << ") to local point (" << local.x() << "|" << local.y() << "|" << local.z() << "). Bent row"; + return local; + } else { // bent part + ROOT::Math::XYZPoint local(local_transformed.x(), // not modified + detector_dimensions.Y() / 2 - row_arc_length, + 0); + LOG(TRACE) << "Transformed global point (" << global.x() << "|" << global.y() << "|" << global.z() + << ") to local point (" << local.x() << "|" << local.y() << "|" << local.z() << "). Bent row"; + return local; + } + } else { + LOG(ERROR) << "m_bent_axis = " << m_bent_axis << " not implemented. Return cartesian coordinates"; + return m_globalToLocal * global; + } } -// Function to get global intercept with a track PositionVector3D> BentPixelDetector::getIntercept(const Track* track) const { - - // FIXME: this is else statement can only be temporary + // FIXME: this else statement can only be temporary if(track->getType() == "GblTrack") { return track->getState(getName()); } else { @@ -245,44 +347,130 @@ PositionVector3D> BentPixelDetector::getIntercept(const Trac track->getState(m_detectorName).Y() + distance * track->getDirection(m_detectorName).Y(), track->getState(m_detectorName).Z() + distance * track->getDirection(m_detectorName).Z()); - if(m_radius == 0) { + // Get and transform track state and direction + PositionVector3D> state_track = track->getState(m_detectorName); + DisplacementVector3D> direction_track = track->getDirection(m_detectorName); + + // Bring track to local (transformed) coordinate system + state_track = m_globalToLocal * state_track; + direction_track = m_globalToLocal.Rotation() * direction_track; + direction_track = direction_track.Unit(); + + // From globalPlanarIntercept get intercept with bent surface of pixel detector + double intercept_parameters[2] = {-999999, -999999}; + if(m_bent_axis == "column") { + if(m_flat_part != 0 || m_radius == 0) { + LOG(WARNING) << "Column bending does not allow for a flat part at the moment or zero radius. Return " + "globalPlanarIntercept"; + return globalPlanarIntercept; + } + + // Define/initialise detector cylinder in local_transformed coordinates + ROOT::Math::PositionVector3D> state_cylinder(0, 0, -m_radius); + ROOT::Math::DisplacementVector3D> direction_cylinder( + 0, 1, 0); // cylinder axis along y (row) + + getInterceptParameters(state_track, direction_track, state_cylinder, direction_cylinder, intercept_parameters); + + } else if(m_bent_axis == "row") { + if(m_radius == 0) { + LOG(WARNING) << "Zero radius. Return globalPlanarIntercept"; + return globalPlanarIntercept; + } + + double row_arc_length = + m_radius * asin((this->getSize().Y() / 2 - m_flat_part - globalPlanarIntercept.Y()) / m_radius) + + m_flat_part; + if(row_arc_length < m_flat_part) { // flat part + return globalPlanarIntercept; + } else { // bent part + ROOT::Math::PositionVector3D> state_cylinder( + 0, this->getSize().Y() / 2 - m_flat_part, -m_radius); + ROOT::Math::DisplacementVector3D> direction_cylinder( + 1, 0, 0); // cylinder axis along x (column) + + getInterceptParameters( + state_track, direction_track, state_cylinder, direction_cylinder, intercept_parameters); + } + } else { + LOG(ERROR) << "m_bent_axis = " << m_bent_axis << " not implemented. Return globalPlanarIntercept"; return globalPlanarIntercept; - } else { // intersect with the curvature, correct - XYZPoint state(track->getState(m_detectorName).X(), - track->getState(m_detectorName).Y(), - track->getState(m_detectorName).Z()); - XYZPoint direction(track->getDirection(m_detectorName).X(), - track->getDirection(m_detectorName).Y(), - track->getDirection(m_detectorName).Z()); - XYZPoint direction_tr = m_globalToLocal.Rotation() * direction; - XYZPoint state_tr = m_globalToLocal * state; - // solve quadratic equation for intersect: alpha*z^2+beta*z+gamma=0 - double alpha = 1 + direction_tr.x() * direction_tr.x() / (direction_tr.z() * direction_tr.z()); - double beta = 2 * state_tr.x() * direction_tr.x() / direction_tr.z() - - 2 * state_tr.z() * direction_tr.x() * direction_tr.x() / (direction_tr.z() * direction_tr.z()) - - 2 * m_radius; - double gamma = - state_tr.x() * state_tr.x() + - state_tr.z() * state_tr.z() * direction_tr.x() * direction_tr.x() / (direction_tr.z() * direction_tr.z()) - - 2 * state_tr.x() * state_tr.z() * direction_tr.x() / direction_tr.z(); - - // check if the equation has a solution - double delta = beta * beta - 4 * alpha * gamma; - if(delta < 0) // return a value outside the chip, none is at the moment not possible + } + + // Select solution according to bending direction + PositionVector3D> localBentIntercept; + if(m_radius <= 0) { // for negative radius select smaller solution + if(intercept_parameters[0] < -999990) { + LOG(TRACE) << "No intercept of track and cylindrcal sensor found. Return intercept outside of sensor area"; return (localToGlobal(getLocalPosition(2 * m_nPixels.X(), 2 * m_nPixels.Y()))); - // change the intercept point with the cylinder with VZ - double z_intercept = (-beta - sqrt(delta)) / (2 * alpha); - // calculate x and y from z-intercept - XYZPoint intercept_tr((z_intercept - state_tr.z()) * direction_tr.x() / direction_tr.z() + state_tr.x(), - (z_intercept - state_tr.z()) * direction_tr.y() / direction_tr.z() + state_tr.y(), - z_intercept); - XYZPoint globalIntercept = m_localToGlobal * intercept_tr; - - PositionVector3D> globalIntercept_cart( - globalIntercept.x(), globalIntercept.y(), globalIntercept.z()); - return globalIntercept_cart; + } else { + localBentIntercept = ROOT::Math::XYZPoint(state_track.x() + intercept_parameters[0] * direction_track.x(), + state_track.y() + intercept_parameters[0] * direction_track.y(), + state_track.z() + intercept_parameters[0] * direction_track.z()); + } + } else { + if(intercept_parameters[1] < -999990) { + LOG(TRACE) << "No intercept of track and cylindrcal sensor found. Return intercept outside of sensor area"; + return (localToGlobal(getLocalPosition(2 * m_nPixels.X(), 2 * m_nPixels.Y()))); + } else { + localBentIntercept = ROOT::Math::XYZPoint(state_track.x() + intercept_parameters[1] * direction_track.x(), + state_track.y() + intercept_parameters[1] * direction_track.y(), + state_track.z() + intercept_parameters[1] * direction_track.z()); + } } + + return m_localToGlobal * localBentIntercept; + } // not GBL track +} + +void BentPixelDetector::getInterceptParameters(const PositionVector3D>& state_track, + const DisplacementVector3D>& direction_track, + const PositionVector3D>& state_cylinder, + const DisplacementVector3D>& direction_cylinder, + double (&result)[2]) const { + // Get factors for quadratic equation: alpha z^2 + beta * z + gamma + double alpha; + double beta; + double gamma; + if(direction_cylinder.X() != 0 && direction_cylinder.Y() == 0 && + direction_cylinder.Z() == 0) { // cylinder along x (column) + alpha = (direction_track.Y() * direction_track.Y()) + (direction_track.Z() * direction_track.Z()); + beta = 2 * (state_track.Y() * direction_track.Y() + state_track.Z() * direction_track.Z() - + direction_track.Y() * state_cylinder.Y() - direction_track.Z() * state_cylinder.Z()); + gamma = + ((state_track.Y() * state_track.Y()) + (state_track.Z() * state_track.Z()) + + (state_cylinder.Y() * state_cylinder.Y()) + (state_cylinder.Z() * state_cylinder.Z()) - + 2 * state_track.Y() * state_cylinder.Y() - 2 * state_track.Z() * state_cylinder.Z() - (m_radius * m_radius)); + } else if(direction_cylinder.X() == 0 && direction_cylinder.Y() != 0 && + direction_cylinder.Z() == 0) { // cylinder along y (row) + alpha = (direction_track.X() * direction_track.X()) + (direction_track.Z() * direction_track.Z()); + beta = 2 * (state_track.X() * direction_track.X() + state_track.Z() * direction_track.Z() - + direction_track.X() * state_cylinder.X() - direction_track.Z() * state_cylinder.Z()); + gamma = + ((state_track.X() * state_track.X()) + (state_track.Z() * state_track.Z()) + + (state_cylinder.X() * state_cylinder.X()) + (state_cylinder.Z() * state_cylinder.Z()) - + 2 * state_track.X() * state_cylinder.X() - 2 * state_track.Z() * state_cylinder.Z() - (m_radius * m_radius)); + } else { + LOG(ERROR) << "Unknown cylinder direction"; + result[0] = -999999; + result[1] = -999999; + return; } + + // Check if the quadratic equation has a solution + double radical = beta * beta - 4 * alpha * gamma; + if(radical < 0 || alpha == 0) { // return a value outside the chip for the track to have no intercept + result[0] = -999999; + result[1] = -999999; + return; + } + + // Solve equation + result[0] = (-beta + sqrt(radical)) / (2 * alpha); + result[1] = (-beta - sqrt(radical)) / (2 * alpha); + // sort in ascending order + std::sort(std::begin(result), std::end(result)); + return; } PositionVector3D> BentPixelDetector::getLocalIntercept(const Track* track) const { @@ -401,25 +589,6 @@ bool BentPixelDetector::isWithinROI(Cluster* cluster) const { return true; } -// Check if a pixel touches any of the pixels in a cluster - -bool BentPixelDetector::isNeighbor(const std::shared_ptr& neighbor, - const std::shared_ptr& cluster, - const int neighbor_radius_row, - const int neighbor_radius_col) { - for(auto pixel : cluster->pixels()) { - int row_distance = abs(pixel->row() - neighbor->row()); - int col_distance = abs(pixel->column() - neighbor->column()); - if(row_distance <= neighbor_radius_row && col_distance <= neighbor_radius_col) { - if(row_distance > 1 || col_distance > 1) { - cluster->setSplit(true); - } - return true; - } - } - return false; -} - XYVector BentPixelDetector::getSize() const { return XYVector(m_pitch.X() * m_nPixels.X(), m_pitch.Y() * m_nPixels.Y()); } @@ -483,3 +652,22 @@ int BentPixelDetector::winding_number(std::pair probe, std::vector pt0, std::pair pt1, std::pair pt2) { return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second)); } + +// Check if a pixel touches any of the pixels in a cluster +bool BentPixelDetector::isNeighbor(const std::shared_ptr& neighbor, + const std::shared_ptr& cluster, + const int neighbor_radius_row, + const int neighbor_radius_col) { + for(auto pixel : cluster->pixels()) { + int row_distance = abs(pixel->row() - neighbor->row()); + int col_distance = abs(pixel->column() - neighbor->column()); + + if(row_distance <= neighbor_radius_row && col_distance <= neighbor_radius_col) { + if(row_distance > 1 || col_distance > 1) { + cluster->setSplit(true); + } + return true; + } + } + return false; +} diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index 3d55bea2..f436793a 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -23,6 +23,8 @@ #include "core/config/Configuration.hpp" #include "core/utils/ROOT.h" #include "core/utils/log.h" +#include "objects/Cluster.hpp" +#include "objects/Pixel.hpp" #include "objects/Track.hpp" namespace corryvreckan { @@ -114,19 +116,21 @@ namespace corryvreckan { // Function to get local position from column (x) and row (y) coordinates PositionVector3D> getLocalPosition(double column, double row) const override; - /** + /** * @brief Transform local coordinates of this detector into global coordinates * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - XYZPoint localToGlobal(XYZPoint local) const override; + //XYZPoint localToGlobal(XYZPoint local) const override; + XYZPoint localToGlobal(XYZPoint local) const; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - XYZPoint globalToLocal(XYZPoint global) const override; + //XYZPoint globalToLocal(XYZPoint global) const override; + XYZPoint globalToLocal(XYZPoint global) const; /** * Transformation from local (sensor) coordinates to in-pixel coordinates @@ -213,10 +217,15 @@ namespace corryvreckan { inline static int isLeft(std::pair pt0, std::pair pt1, std::pair pt2); static int winding_number(std::pair probe, std::vector> polygon); - // For planar detector + // Function to facilitate the incept calculation of sraight tracks with cylindrical sensor surface + void GetInterceptParameters(const ROOT::Math::PositionVector3D> &state_track, + const ROOT::Math::DisplacementVector3D> &direction_track, + const ROOT::Math::PositionVector3D> &state_cylinder, + const ROOT::Math::DisplacementVector3D> &direction_cylinder, + double (&result)[2]) const; + + // For bent pixel detector XYVector m_pitch{}; - double m_radius; - double m_y_0; XYVector m_spatial_resolution{}; ROOT::Math::DisplacementVector2D> m_nPixels{}; std::vector> m_roi{}; @@ -224,8 +233,12 @@ namespace corryvreckan { ROOT::Math::XYZPoint m_displacement; ROOT::Math::XYZVector m_orientation; std::string m_orientation_mode; - std::string m_coordinates; + // bent geometry configuration parameters + double m_radius; + double m_flat_part; + std::string m_bent_axis; + std::string m_coordinates; }; } // namespace corryvreckan -#endif // CORRYVRECKAN_DETECTOR_H +#endif // CORRYVRECKAN_BENTPIXELDETECTOR_H diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index 8c217479..e2b8ffad 100644 --- a/src/core/detector/Detector.cpp +++ b/src/core/detector/Detector.cpp @@ -15,8 +15,8 @@ #include "Math/RotationZ.h" #include "Math/RotationZYX.h" -#include "BentPixelDetector.hpp" #include "Detector.hpp" +#include "BentPixelDetector.hpp" #include "PixelDetector.hpp" #include "core/utils/log.h" #include "exceptions.h" @@ -86,12 +86,11 @@ std::shared_ptr corryvreckan::Detector::factory(const Configuration& c std::transform(coordinates.begin(), coordinates.end(), coordinates.begin(), ::tolower); if(coordinates == "cartesian") { return std::make_shared(config); - } - if(coordinates == "cartesian-bent") { + } else if (coordinates == "cartesian-bent") { return std::make_shared(config); + } else { + throw InvalidValueError(config, "coordinates", "Coordinates can only be set to \"cartesian\" or \"cartesian-bent\" for now"); } - // coordinates was not found, throw an error - throw InvalidValueError(config, "coordinates", "Coordinates can only set to be cartesian now"); } double Detector::getTimeResolution() const { diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index eb1dea22..3ad0fbf3 100644 --- a/src/core/detector/Detector.hpp +++ b/src/core/detector/Detector.hpp @@ -280,14 +280,14 @@ namespace corryvreckan { * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - virtual XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; + XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - virtual XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; + XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; /** * @brief Check whether given track is within the detector's region-of-interest @@ -373,4 +373,5 @@ namespace corryvreckan { }; } // namespace corryvreckan +#include "PixelDetector.hpp" #endif // CORRYVRECKAN_DETECTOR_H diff --git a/src/core/detector/PixelDetector.hpp b/src/core/detector/PixelDetector.hpp index 74433909..211525db 100644 --- a/src/core/detector/PixelDetector.hpp +++ b/src/core/detector/PixelDetector.hpp @@ -210,7 +210,6 @@ namespace corryvreckan { ROOT::Math::XYZPoint m_displacement; ROOT::Math::XYZVector m_orientation; std::string m_orientation_mode; - std::string m_coordinates; }; } // namespace corryvreckan diff --git a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp index f99573b7..5d417eae 100644 --- a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp +++ b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp @@ -218,7 +218,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { std::vector fscaz(nParameters, 0.); std::vector shearx(nParameters, 0.); std::vector sheary(nParameters, 0.); - std::vector shearrz(nParameters, 0.); + std::vector shearrz(nParameters, 0.); m_constraints.clear(); for(const auto& det : get_detectors()) { @@ -240,30 +240,39 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { shearx[i] = sz; sheary[i + nPlanes] = sz; fscaz[i + 2 * nPlanes] = sz; - shearrz[i + 5 * nPlanes] = sz; + shearrz[i + 5 * nPlanes] = sz; } const std::vector constraints = {true, true, true, true, false, false, true, true, true, true}; // Put the constraints information in the basket. - if(constraints[0] && m_dofs[0]) - addConstraint(ftx, 0.0); - if(constraints[1] && m_dofs[0]) - addConstraint(shearx, 0.); - if(constraints[2] && m_dofs[1]) - addConstraint(fty, 0.0); - if(constraints[3] && m_dofs[1]) - addConstraint(sheary, 0.); - if(constraints[4] && m_dofs[2]) - addConstraint(ftz, 0.0); + if(constraints[0] && m_dofs[0]) { + addConstraint(ftx, 0.0); + } + if(constraints[1] && m_dofs[0]) { + addConstraint(shearx, 0.); + } + if(constraints[2] && m_dofs[1]) { + addConstraint(fty, 0.0); + } + if(constraints[3] && m_dofs[1]) { + addConstraint(sheary, 0.); + } + if(constraints[4] && m_dofs[2]) { + addConstraint(ftz, 0.0); + } // if (constraints[5] && m_dofs[2]) addConstraint(fscaz, 0.0); - if(constraints[6] && m_dofs[3]) - addConstraint(frx, 0.0); - if(constraints[7] && m_dofs[4]) - addConstraint(fry, 0.0); - if(constraints[8] && m_dofs[5]) + if(constraints[6] && m_dofs[3]) { + addConstraint(frx, 0.0); + } + if(constraints[7] && m_dofs[4]) { + addConstraint(fry, 0.0); + } + if(constraints[8] && m_dofs[5]) { addConstraint(frz, 0.0); - if(constraints[9] && m_dofs[5]) + } + if(constraints[9] && m_dofs[6]) { addConstraint(shearrz, 0.0); + } } //============================================================================= diff --git a/src/modules/AnalysisDUT/AnalysisDUT.h b/src/modules/AnalysisDUT/AnalysisDUT.h index d179bda9..f2202b5d 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -54,13 +54,7 @@ namespace corryvreckan { TProfile2D* hPixelRawValueMapAssoc; TH1F* associatedTracksVersusTime; - TH1F *residualsX, *residualsY, *residualsZ, *residualsPos; - TH2F *residualsXvsXhit, *residualsXvsYhit; - TH2F *residualsYvsXhit, *residualsYvsYhit; - TH2F* residualsXvsresidualsY; - TProfile2D* residualsXprofile; - TProfile2D* residualsYprofile; - + TH1F *residualsX, *residualsY, *residualsPos; TH2F* residualsPosVsresidualsTime; TH1F *residualsX1pix, *residualsY1pix; diff --git a/src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp b/src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp index 566106a8..b493e50a 100644 --- a/src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp +++ b/src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp @@ -385,7 +385,7 @@ Event::Position EventLoaderEUDAQ2::is_within_event(const std::shared_ptr " << Units::display(clipboard->getEvent()->end(), {"us", "ns"}); - } else if(position == Event::Position::DURING) { + } else { // check if event has valid trigger ID (flag = 0x10): if(evt->IsFlagTrigger()) { if(veto_triggers_ && !clipboard->getEvent()->triggerList().empty()) { @@ -566,10 +566,7 @@ StatusCode EventLoaderEUDAQ2::run(const std::shared_ptr& clipboard) { } // If this event was after the current event or if we have not enough information, stop reading: - if(current_position == Event::Position::AFTER) { - break; - } else if(current_position == Event::Position::UNKNOWN) { - event_.reset(); + if(current_position == Event::Position::AFTER || current_position == Event::Position::UNKNOWN) { break; } diff --git a/src/objects/Event.hpp b/src/objects/Event.hpp index 3e15f971..2ca45371 100644 --- a/src/objects/Event.hpp +++ b/src/objects/Event.hpp @@ -12,7 +12,6 @@ #define CORRYVRECKAN_EVENT_H 1 #include "Object.hpp" -#include "exceptions.h" namespace corryvreckan { @@ -38,11 +37,7 @@ namespace corryvreckan { * @param trigger_list Optional list of triggers assigned to this event, containing their ID and timestamps */ Event(double start, double end, std::map trigger_list = std::map()) - : Object(start), end_(end), trigger_list_(std::move(trigger_list)) { - if(end < start) { - throw InvalidEventError(typeid(*this), "Negative Event duration"); - } - }; + : Object(start), end_(end), trigger_list_(std::move(trigger_list)){}; /** * @brief Static member function to obtain base class for storage on the clipboard. diff --git a/src/objects/exceptions.h b/src/objects/exceptions.h index 52317ee4..e779f6ca 100644 --- a/src/objects/exceptions.h +++ b/src/objects/exceptions.h @@ -110,20 +110,6 @@ namespace corryvreckan { } }; - class InvalidEventError : public ObjectError { - public: - /** - * @brief InvalidEventError - * @param source - */ - - explicit InvalidEventError(const std::type_info& source, const std::string msg = "") { - error_message_ += " Event Object "; - error_message_ += corryvreckan::demangle((source.name())); - error_message_ += " invalid: "; - error_message_ += msg; - } - }; } // namespace corryvreckan #endif /* CORRYVRECKAN_OBJECT_EXCEPTIONS_H */ -- GitLab From c266dded4c554ae9b46c0d19e3dd432401f650b5 Mon Sep 17 00:00:00 2001 From: Mihail Bogdan Blidaru Date: Wed, 3 Mar 2021 12:19:19 +0100 Subject: [PATCH 16/37] Clean-up y_0. --- src/core/detector/BentPixelDetector.hpp | 36 ++++++++++++------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index f436793a..99c6c62a 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -13,10 +13,9 @@ #include #include -#include -#include -#include +#include "Math/DisplacementVector2D.h" #include "Math/Transform3D.h" +#include "Math/Vector2D.h" #include "Math/Vector3D.h" #include "Detector.hpp" @@ -116,21 +115,19 @@ namespace corryvreckan { // Function to get local position from column (x) and row (y) coordinates PositionVector3D> getLocalPosition(double column, double row) const override; - /** + /** * @brief Transform local coordinates of this detector into global coordinates * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - //XYZPoint localToGlobal(XYZPoint local) const override; - XYZPoint localToGlobal(XYZPoint local) const; + XYZPoint localToGlobal(XYZPoint local) const override; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - //XYZPoint globalToLocal(XYZPoint global) const override; - XYZPoint globalToLocal(XYZPoint global) const; + XYZPoint globalToLocal(XYZPoint global) const override; /** * Transformation from local (sensor) coordinates to in-pixel coordinates @@ -217,12 +214,13 @@ namespace corryvreckan { inline static int isLeft(std::pair pt0, std::pair pt1, std::pair pt2); static int winding_number(std::pair probe, std::vector> polygon); - // Function to facilitate the incept calculation of sraight tracks with cylindrical sensor surface - void GetInterceptParameters(const ROOT::Math::PositionVector3D> &state_track, - const ROOT::Math::DisplacementVector3D> &direction_track, - const ROOT::Math::PositionVector3D> &state_cylinder, - const ROOT::Math::DisplacementVector3D> &direction_cylinder, - double (&result)[2]) const; + // Function to facilitate the incept calculation of sraight tracks with cylindrical sensor surface + void + getInterceptParameters(const ROOT::Math::PositionVector3D>& state_track, + const ROOT::Math::DisplacementVector3D>& direction_track, + const ROOT::Math::PositionVector3D>& state_cylinder, + const ROOT::Math::DisplacementVector3D>& direction_cylinder, + double (&result)[2]) const; // For bent pixel detector XYVector m_pitch{}; @@ -233,11 +231,11 @@ namespace corryvreckan { ROOT::Math::XYZPoint m_displacement; ROOT::Math::XYZVector m_orientation; std::string m_orientation_mode; - // bent geometry configuration parameters - double m_radius; - double m_flat_part; - std::string m_bent_axis; - std::string m_coordinates; + // bent geometry configuration parameters + double m_radius; + double m_flat_part; + std::string m_bent_axis; + std::string m_coordinates; }; } // namespace corryvreckan -- GitLab From dfa7dd1242ba60190d5e2862fb52fa57893a355e Mon Sep 17 00:00:00 2001 From: p-becht <> Date: Thu, 4 Mar 2021 13:10:23 +0100 Subject: [PATCH 17/37] ran clang-format to fix formatting errors in pipeline --- .../include/BorderedBandMatrix.h | 2 +- 3rdparty/GeneralBrokenLines/include/GblData.h | 2 +- .../GeneralBrokenLines/include/GblPoint.h | 2 +- .../include/GblTrajectory.h | 4 +- .../GeneralBrokenLines/include/MilleBinary.h | 2 +- 3rdparty/GeneralBrokenLines/include/VMatrix.h | 2 +- .../src/BorderedBandMatrix.cpp | 2 +- 3rdparty/GeneralBrokenLines/src/GblData.cpp | 2 +- 3rdparty/GeneralBrokenLines/src/GblPoint.cpp | 2 +- .../GeneralBrokenLines/src/GblTrajectory.cpp | 2 +- .../GeneralBrokenLines/src/MilleBinary.cpp | 2 +- 3rdparty/GeneralBrokenLines/src/VMatrix.cpp | 2 +- src/core/detector/Detector.cpp | 9 +++-- .../AlignmentMillepede/AlignmentMillepede.cpp | 38 +++++++++---------- 14 files changed, 37 insertions(+), 36 deletions(-) diff --git a/3rdparty/GeneralBrokenLines/include/BorderedBandMatrix.h b/3rdparty/GeneralBrokenLines/include/BorderedBandMatrix.h index dacb1120..61a79fff 100644 --- a/3rdparty/GeneralBrokenLines/include/BorderedBandMatrix.h +++ b/3rdparty/GeneralBrokenLines/include/BorderedBandMatrix.h @@ -100,5 +100,5 @@ namespace gbl { VMatrix invertBand(); VMatrix bandOfAVAT(const VMatrix& anArray, const VSymMatrix& aSymArray) const; }; -} +} // namespace gbl #endif /* BORDEREDBANDMATRIX_H_ */ diff --git a/3rdparty/GeneralBrokenLines/include/GblData.h b/3rdparty/GeneralBrokenLines/include/GblData.h index 8f784288..4e12fee7 100644 --- a/3rdparty/GeneralBrokenLines/include/GblData.h +++ b/3rdparty/GeneralBrokenLines/include/GblData.h @@ -262,6 +262,6 @@ namespace gbl { } } } -} +} // namespace gbl #endif /* GBLDATA_H_ */ diff --git a/3rdparty/GeneralBrokenLines/include/GblPoint.h b/3rdparty/GeneralBrokenLines/include/GblPoint.h index 5bbbc17d..ad25c662 100644 --- a/3rdparty/GeneralBrokenLines/include/GblPoint.h +++ b/3rdparty/GeneralBrokenLines/include/GblPoint.h @@ -430,5 +430,5 @@ namespace gbl { } } } -} +} // namespace gbl #endif /* GBLPOINT_H_ */ diff --git a/3rdparty/GeneralBrokenLines/include/GblTrajectory.h b/3rdparty/GeneralBrokenLines/include/GblTrajectory.h index 3de9aba2..5dbb07ff 100644 --- a/3rdparty/GeneralBrokenLines/include/GblTrajectory.h +++ b/3rdparty/GeneralBrokenLines/include/GblTrajectory.h @@ -188,7 +188,7 @@ namespace gbl { unsigned int numOffsets; ///< Number of (points with) offsets on trajectory unsigned int numInnerTransformations; ///< Number of inner transformations to external parameters unsigned int numInnerTransOffsets; ///< Number of (points with) offsets affected by inner transformations to external - ///parameters + /// parameters unsigned int numCurvature; ///< Number of curvature parameters (0 or 1) or external parameters unsigned int numParameters; ///< Number of fit parameters unsigned int numLocals; ///< Total number of (additional) local parameters @@ -338,5 +338,5 @@ namespace gbl { : innerTransformations[0].cols() + numInnerTransformations; construct(); // construct (composed) trajectory } -} +} // namespace gbl #endif /* GBLTRAJECTORY_H_ */ diff --git a/3rdparty/GeneralBrokenLines/include/MilleBinary.h b/3rdparty/GeneralBrokenLines/include/MilleBinary.h index aa350acf..8833fb71 100644 --- a/3rdparty/GeneralBrokenLines/include/MilleBinary.h +++ b/3rdparty/GeneralBrokenLines/include/MilleBinary.h @@ -98,5 +98,5 @@ namespace gbl { std::vector doubleBuffer; ///< Double buffer bool doublePrecision; ///< Flag for storage in as *double* values }; -} +} // namespace gbl #endif /* MILLEBINARY_H_ */ diff --git a/3rdparty/GeneralBrokenLines/include/VMatrix.h b/3rdparty/GeneralBrokenLines/include/VMatrix.h index d7fe7295..5dcd273d 100644 --- a/3rdparty/GeneralBrokenLines/include/VMatrix.h +++ b/3rdparty/GeneralBrokenLines/include/VMatrix.h @@ -125,5 +125,5 @@ namespace gbl { inline double VSymMatrix::operator()(unsigned int iRow, unsigned int iCol) const { return theVec[(iRow * iRow + iRow) / 2 + iCol]; // assuming iCol <= iRow } -} +} // namespace gbl #endif /* VMATRIX_H_ */ diff --git a/3rdparty/GeneralBrokenLines/src/BorderedBandMatrix.cpp b/3rdparty/GeneralBrokenLines/src/BorderedBandMatrix.cpp index 724cbdec..bcf52a24 100644 --- a/3rdparty/GeneralBrokenLines/src/BorderedBandMatrix.cpp +++ b/3rdparty/GeneralBrokenLines/src/BorderedBandMatrix.cpp @@ -381,4 +381,4 @@ namespace gbl { } return aBand; } -} +} // namespace gbl diff --git a/3rdparty/GeneralBrokenLines/src/GblData.cpp b/3rdparty/GeneralBrokenLines/src/GblData.cpp index 6fee6135..8e19d9ac 100644 --- a/3rdparty/GeneralBrokenLines/src/GblData.cpp +++ b/3rdparty/GeneralBrokenLines/src/GblData.cpp @@ -261,4 +261,4 @@ namespace gbl { derLocal = &moreDerivatives[0]; } } -} +} // namespace gbl diff --git a/3rdparty/GeneralBrokenLines/src/GblPoint.cpp b/3rdparty/GeneralBrokenLines/src/GblPoint.cpp index fd137fac..3b839cb2 100644 --- a/3rdparty/GeneralBrokenLines/src/GblPoint.cpp +++ b/3rdparty/GeneralBrokenLines/src/GblPoint.cpp @@ -519,4 +519,4 @@ namespace gbl { } } } -} +} // namespace gbl diff --git a/3rdparty/GeneralBrokenLines/src/GblTrajectory.cpp b/3rdparty/GeneralBrokenLines/src/GblTrajectory.cpp index a904c7bb..afe91656 100644 --- a/3rdparty/GeneralBrokenLines/src/GblTrajectory.cpp +++ b/3rdparty/GeneralBrokenLines/src/GblTrajectory.cpp @@ -1518,4 +1518,4 @@ namespace gbl { itData->printData(); } } -} +} // namespace gbl diff --git a/3rdparty/GeneralBrokenLines/src/MilleBinary.cpp b/3rdparty/GeneralBrokenLines/src/MilleBinary.cpp index 39d6a94d..295e64d0 100644 --- a/3rdparty/GeneralBrokenLines/src/MilleBinary.cpp +++ b/3rdparty/GeneralBrokenLines/src/MilleBinary.cpp @@ -126,4 +126,4 @@ namespace gbl { else floatBuffer.resize(1); } -} +} // namespace gbl diff --git a/3rdparty/GeneralBrokenLines/src/VMatrix.cpp b/3rdparty/GeneralBrokenLines/src/VMatrix.cpp index 7b2ec569..6b21b015 100644 --- a/3rdparty/GeneralBrokenLines/src/VMatrix.cpp +++ b/3rdparty/GeneralBrokenLines/src/VMatrix.cpp @@ -424,4 +424,4 @@ namespace gbl { } return nrank; } -} +} // namespace gbl diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index e2b8ffad..535707d7 100644 --- a/src/core/detector/Detector.cpp +++ b/src/core/detector/Detector.cpp @@ -15,8 +15,8 @@ #include "Math/RotationZ.h" #include "Math/RotationZYX.h" -#include "Detector.hpp" #include "BentPixelDetector.hpp" +#include "Detector.hpp" #include "PixelDetector.hpp" #include "core/utils/log.h" #include "exceptions.h" @@ -86,10 +86,11 @@ std::shared_ptr corryvreckan::Detector::factory(const Configuration& c std::transform(coordinates.begin(), coordinates.end(), coordinates.begin(), ::tolower); if(coordinates == "cartesian") { return std::make_shared(config); - } else if (coordinates == "cartesian-bent") { + } else if(coordinates == "cartesian-bent") { return std::make_shared(config); - } else { - throw InvalidValueError(config, "coordinates", "Coordinates can only be set to \"cartesian\" or \"cartesian-bent\" for now"); + } else { + throw InvalidValueError( + config, "coordinates", "Coordinates can only be set to \"cartesian\" or \"cartesian-bent\" for now"); } } diff --git a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp index 5d417eae..6de0924d 100644 --- a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp +++ b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp @@ -218,7 +218,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { std::vector fscaz(nParameters, 0.); std::vector shearx(nParameters, 0.); std::vector sheary(nParameters, 0.); - std::vector shearrz(nParameters, 0.); + std::vector shearrz(nParameters, 0.); m_constraints.clear(); for(const auto& det : get_detectors()) { @@ -240,39 +240,39 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { shearx[i] = sz; sheary[i + nPlanes] = sz; fscaz[i + 2 * nPlanes] = sz; - shearrz[i + 5 * nPlanes] = sz; + shearrz[i + 5 * nPlanes] = sz; } const std::vector constraints = {true, true, true, true, false, false, true, true, true, true}; // Put the constraints information in the basket. if(constraints[0] && m_dofs[0]) { - addConstraint(ftx, 0.0); - } + addConstraint(ftx, 0.0); + } if(constraints[1] && m_dofs[0]) { - addConstraint(shearx, 0.); - } + addConstraint(shearx, 0.); + } if(constraints[2] && m_dofs[1]) { - addConstraint(fty, 0.0); - } + addConstraint(fty, 0.0); + } if(constraints[3] && m_dofs[1]) { - addConstraint(sheary, 0.); - } + addConstraint(sheary, 0.); + } if(constraints[4] && m_dofs[2]) { - addConstraint(ftz, 0.0); - } + addConstraint(ftz, 0.0); + } // if (constraints[5] && m_dofs[2]) addConstraint(fscaz, 0.0); if(constraints[6] && m_dofs[3]) { - addConstraint(frx, 0.0); - } + addConstraint(frx, 0.0); + } if(constraints[7] && m_dofs[4]) { - addConstraint(fry, 0.0); - } + addConstraint(fry, 0.0); + } if(constraints[8] && m_dofs[5]) { addConstraint(frz, 0.0); - } - if(constraints[9] && m_dofs[6]) { + } + if(constraints[9] && m_dofs[6]) { addConstraint(shearrz, 0.0); - } + } } //============================================================================= -- GitLab From f70c9268cdcd2674dd6fcdae8284d154511e7bab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 2 Mar 2021 15:55:35 +0100 Subject: [PATCH 18/37] Add residuals X vs residuals Y correlation plots in AnalysisDUT --- src/modules/AnalysisDUT/AnalysisDUT.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.h b/src/modules/AnalysisDUT/AnalysisDUT.h index f2202b5d..d179bda9 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -54,7 +54,13 @@ namespace corryvreckan { TProfile2D* hPixelRawValueMapAssoc; TH1F* associatedTracksVersusTime; - TH1F *residualsX, *residualsY, *residualsPos; + TH1F *residualsX, *residualsY, *residualsZ, *residualsPos; + TH2F *residualsXvsXhit, *residualsXvsYhit; + TH2F *residualsYvsXhit, *residualsYvsYhit; + TH2F* residualsXvsresidualsY; + TProfile2D* residualsXprofile; + TProfile2D* residualsYprofile; + TH2F* residualsPosVsresidualsTime; TH1F *residualsX1pix, *residualsY1pix; -- GitLab From e8571edd5f3777025dc192c9a5b2762b296c48d7 Mon Sep 17 00:00:00 2001 From: p-becht <> Date: Fri, 12 Mar 2021 14:16:11 +0100 Subject: [PATCH 19/37] cleaned up BentDetector, added virtual to globalToLocal and so on --- src/core/detector/Detector.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index 3ad0fbf3..68cc69d9 100644 --- a/src/core/detector/Detector.hpp +++ b/src/core/detector/Detector.hpp @@ -280,14 +280,14 @@ namespace corryvreckan { * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; + virtual XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; + virtual XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; /** * @brief Check whether given track is within the detector's region-of-interest -- GitLab From 36d30af7907bb7fc394c7c2fe3d1c181bccfcf66 Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Mon, 29 Mar 2021 14:37:05 +0200 Subject: [PATCH 20/37] AnalysisTelescope: rename histograms and clarify axis titles: px - seed --- src/modules/AnalysisTelescope/AnalysisTelescope.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp index dbcec343..e229104d 100644 --- a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp +++ b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp @@ -42,13 +42,13 @@ void AnalysisTelescope::initialize() { title = detector->getName() + " Telescope resolution Y;y_{track}-y_{MC} [mm];events"; telescopeResolutionY[detector->getName()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -0.2, 0.2); } else { - std::string title = detector->getName() + " Biased residual X (local);x_{track}-x_{cluster} [mm];events"; + std::string title = detector->getName() + " Biased residual X (local);x_{track}-x_{seed} [mm];events"; telescopeResidualsLocalX[detector->getName()] = new TH1F("residualX_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased residual Y (local);y_{track}-y_{cluster} [mm];events"; + title = detector->getName() + " Biased residual Y (local);y_{track}-y_{seed} [mm];events"; telescopeResidualsLocalY[detector->getName()] = new TH1F("residualY_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased residual X (global);x_{track}-x_{cluster} [mm];events"; + title = detector->getName() + " Biased residual X (global);x_{track}-x_{seed} [mm];events"; telescopeResidualsX[detector->getName()] = new TH1F("residualX_global", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased residual Y (global);y_{track}-y_{cluster} [mm];events"; + title = detector->getName() + " Biased residual Y (global);y_{track}-y_{seed} [mm];events"; telescopeResidualsY[detector->getName()] = new TH1F("residualY_global", title.c_str(), 400, -0.2, 0.2); title = detector->getName() + " Biased MC residual X (local);x_{track}-x_{MC} [mm];events"; -- GitLab From 321108f187df5889c29d4113a99aa18865b039e1 Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Tue, 30 Mar 2021 13:27:50 +0200 Subject: [PATCH 21/37] AnalysisTelescope: fix hist labels that were accidentally change from cluster to seed --- src/modules/AnalysisTelescope/AnalysisTelescope.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp index e229104d..dbcec343 100644 --- a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp +++ b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp @@ -42,13 +42,13 @@ void AnalysisTelescope::initialize() { title = detector->getName() + " Telescope resolution Y;y_{track}-y_{MC} [mm];events"; telescopeResolutionY[detector->getName()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -0.2, 0.2); } else { - std::string title = detector->getName() + " Biased residual X (local);x_{track}-x_{seed} [mm];events"; + std::string title = detector->getName() + " Biased residual X (local);x_{track}-x_{cluster} [mm];events"; telescopeResidualsLocalX[detector->getName()] = new TH1F("residualX_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased residual Y (local);y_{track}-y_{seed} [mm];events"; + title = detector->getName() + " Biased residual Y (local);y_{track}-y_{cluster} [mm];events"; telescopeResidualsLocalY[detector->getName()] = new TH1F("residualY_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased residual X (global);x_{track}-x_{seed} [mm];events"; + title = detector->getName() + " Biased residual X (global);x_{track}-x_{cluster} [mm];events"; telescopeResidualsX[detector->getName()] = new TH1F("residualX_global", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased residual Y (global);y_{track}-y_{seed} [mm];events"; + title = detector->getName() + " Biased residual Y (global);y_{track}-y_{cluster} [mm];events"; telescopeResidualsY[detector->getName()] = new TH1F("residualY_global", title.c_str(), 400, -0.2, 0.2); title = detector->getName() + " Biased MC residual X (local);x_{track}-x_{MC} [mm];events"; -- GitLab From a0486ae3ca596464b64268ff2b0ea898001a5aee Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Tue, 30 Mar 2021 13:41:11 +0200 Subject: [PATCH 22/37] AnalysisDUT: rename histogram titles for more clarity --- src/modules/AnalysisDUT/AnalysisDUT.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.h b/src/modules/AnalysisDUT/AnalysisDUT.h index d179bda9..528bc692 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -54,7 +54,7 @@ namespace corryvreckan { TProfile2D* hPixelRawValueMapAssoc; TH1F* associatedTracksVersusTime; - TH1F *residualsX, *residualsY, *residualsZ, *residualsPos; + TH1F *residualsX, *residualsY, *residualsPos; TH2F *residualsXvsXhit, *residualsXvsYhit; TH2F *residualsYvsXhit, *residualsYvsYhit; TH2F* residualsXvsresidualsY; -- GitLab From c1d82b2aa22b92678e23e42b77a3b47741fa568d Mon Sep 17 00:00:00 2001 From: Simon Spannagel Date: Fri, 9 Apr 2021 19:13:38 +0200 Subject: [PATCH 23/37] Clipboard: add methods to remove data from the clipboard again --- src/core/clipboard/Clipboard.hpp | 10 ++++++++++ src/core/clipboard/Clipboard.tpp | 23 +++++++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/src/core/clipboard/Clipboard.hpp b/src/core/clipboard/Clipboard.hpp index a07bbcb3..94c5f31e 100644 --- a/src/core/clipboard/Clipboard.hpp +++ b/src/core/clipboard/Clipboard.hpp @@ -178,6 +178,16 @@ namespace corryvreckan { bool append = false); /** +<<<<<<< HEAD +======= + * Helper to remove a single object from the clipboard + * @param data Existing data the object should be removed from + * @param object Object to be removed + */ + template void remove_data(std::vector>& data, const std::shared_ptr& object); + + /** +>>>>>>> Clipboard: add methods to remove data from the clipboard again * Helper to remove a set of objects from the clipboard * @param storage_element Storage element concerned * @param objects Objects to be removed diff --git a/src/core/clipboard/Clipboard.tpp b/src/core/clipboard/Clipboard.tpp index 1c4d4124..769712d6 100644 --- a/src/core/clipboard/Clipboard.tpp +++ b/src/core/clipboard/Clipboard.tpp @@ -17,7 +17,12 @@ namespace corryvreckan { } template void Clipboard::removeData(std::shared_ptr object, const std::string& key) { +<<<<<<< HEAD remove_data(data_, std::vector>{std::move(object)}, key); +======= + auto data = get_data(data_, key); + remove_data(std::move(data), std::move(object)); +>>>>>>> Clipboard: add methods to remove data from the clipboard again } template void Clipboard::removeData(std::vector>& objects, const std::string& key) { @@ -108,6 +113,7 @@ namespace corryvreckan { } template +<<<<<<< HEAD void Clipboard::remove_data(ClipboardData& storage_element, const std::vector>& objects, const std::string& key) { @@ -122,6 +128,23 @@ namespace corryvreckan { if(object_iterator != data->end()) { data->erase(object_iterator); } +======= + void Clipboard::remove_data(std::vector>& data, const std::shared_ptr& object) { + auto object_iterator = data.find(object); + if(object_iterator != data.end()) { + data.erase(object_iterator); + } + } + + template + void Clipboard::remove_data(ClipboardData& storage_element, + const std::vector>& objects, + const std::string& key) { + // Try to get existing objects: + auto data = get_data(storage_element, key); + for(const auto& object : objects) { + remove_data(data, object, key); +>>>>>>> Clipboard: add methods to remove data from the clipboard again } } -- GitLab From 040794d98240d28a8046a3561cc3bced3267bbf5 Mon Sep 17 00:00:00 2001 From: Lennart Huth Date: Wed, 21 Apr 2021 13:34:32 +0200 Subject: [PATCH 24/37] don't double include eigen3 --- 3rdparty/GeneralBrokenLines/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/3rdparty/GeneralBrokenLines/CMakeLists.txt b/3rdparty/GeneralBrokenLines/CMakeLists.txt index 1e55a345..6f9d8d2a 100644 --- a/3rdparty/GeneralBrokenLines/CMakeLists.txt +++ b/3rdparty/GeneralBrokenLines/CMakeLists.txt @@ -12,7 +12,6 @@ PROJECT(GBL VERSION "2.2.0") # Enable ROOT support FIND_PACKAGE(ROOT REQUIRED) ADD_DEFINITIONS(-DGBL_EIGEN_SUPPORT_ROOT) -FIND_PACKAGE(Eigen3 REQUIRED) # include directories INCLUDE_DIRECTORIES(SYSTEM BEFORE ./include ${ROOT_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIR} ) -- GitLab From b63d307dadbff37ef35f473cdae2efd8114189dc Mon Sep 17 00:00:00 2001 From: Lennart Huth Date: Wed, 21 Apr 2021 14:33:18 +0200 Subject: [PATCH 25/37] converting skip_time only once to ps --- src/core/clipboard/Clipboard.tpp | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/core/clipboard/Clipboard.tpp b/src/core/clipboard/Clipboard.tpp index 769712d6..1c4d4124 100644 --- a/src/core/clipboard/Clipboard.tpp +++ b/src/core/clipboard/Clipboard.tpp @@ -17,12 +17,7 @@ namespace corryvreckan { } template void Clipboard::removeData(std::shared_ptr object, const std::string& key) { -<<<<<<< HEAD remove_data(data_, std::vector>{std::move(object)}, key); -======= - auto data = get_data(data_, key); - remove_data(std::move(data), std::move(object)); ->>>>>>> Clipboard: add methods to remove data from the clipboard again } template void Clipboard::removeData(std::vector>& objects, const std::string& key) { @@ -113,7 +108,6 @@ namespace corryvreckan { } template -<<<<<<< HEAD void Clipboard::remove_data(ClipboardData& storage_element, const std::vector>& objects, const std::string& key) { @@ -128,23 +122,6 @@ namespace corryvreckan { if(object_iterator != data->end()) { data->erase(object_iterator); } -======= - void Clipboard::remove_data(std::vector>& data, const std::shared_ptr& object) { - auto object_iterator = data.find(object); - if(object_iterator != data.end()) { - data.erase(object_iterator); - } - } - - template - void Clipboard::remove_data(ClipboardData& storage_element, - const std::vector>& objects, - const std::string& key) { - // Try to get existing objects: - auto data = get_data(storage_element, key); - for(const auto& object : objects) { - remove_data(data, object, key); ->>>>>>> Clipboard: add methods to remove data from the clipboard again } } -- GitLab From 462dbb0df1ca8ee11e778858cd89e23aabf840fa Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Fri, 23 Apr 2021 11:46:30 +0200 Subject: [PATCH 26/37] AnalysisDUT: unify naming convention of member variables --- src/modules/AnalysisDUT/AnalysisDUT.cpp | 41 +++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index d99556d7..bc4477b8 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -303,6 +303,47 @@ void AnalysisDUT::initialize() { 499.5); residualsXprofile = new TProfile2D("residualsXprofile", "Residual in X;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); residualsYprofile = new TProfile2D("residualsYprofile", "Residual in Y;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); +<<<<<<< HEAD +======= + residualsY = new TH1F("residualsY", "Residual in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsYvsXhit = new TH2F( + "residualsYvsXhit", "Residual in Y;y_{track}-y_{hit} [mm];x_{hit} [mm];# entries", 800, -0.1, 0.1, 40, -20, 20); + residualsYvsYhit = new TH2F( + "residualsYvsYhit", "Residual in Y;y_{track}-y_{hit} [mm];y_{hit} [mm];# entries", 800, -0.1, 0.1, 40, -20, 20); + residualsZ = new TH1F("residualsZ", "Residual in Z;z_{track}-z_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsPos = new TH1F( + "residualsPos", "Absolute distance between track and hit;|pos_{track}-pos_{hit}| [mm];# entries", 800, -0.1, 0.1); + residualsPosVsresidualsTime = + new TH2F("residualsPosVsresidualsTime", + "Time vs. absolute position residuals;time_{track}-time_{hit} [ns];|pos_{track}-pos_{hit}| [mm];# entries", + n_timebins_, + -n_timebins_ / 2. * time_binning_ - time_binning_ / 2., + n_timebins_ / 2. * time_binning_ - time_binning_ / 2., + 800, + 0., + 0.2); + + residualsX1pix = + new TH1F("residualsX1pix", "Residual for 1-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsY1pix = + new TH1F("residualsY1pix", "Residual for 1-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsX2pix = + new TH1F("residualsX2pix", "Residual for 2-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsY2pix = + new TH1F("residualsY2pix", "Residual for 2-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsX3pix = + new TH1F("residualsX3pix", "Residual for 3-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsY3pix = + new TH1F("residualsY3pix", "Residual for 3-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsX4pix = + new TH1F("residualsX4pix", "Residual for 4-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsY4pix = + new TH1F("residualsY4pix", "Residual for 4-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsXatLeast5pix = new TH1F( + "residualsXatLeast5pix", "Residual for >= 5-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); + residualsYatLeast5pix = new TH1F( + "residualsYatLeast5pix", "Residual for >= 5-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); +>>>>>>> AnalysisDUT: unify naming convention of member variables clusterChargeAssoc = new TH1F("clusterChargeAssociated", "Charge distribution of associated clusters;cluster charge [e];# entries", -- GitLab From 04ac0198664ff9064c59f462c0b1df3d63004443 Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Fri, 23 Apr 2021 13:35:55 +0200 Subject: [PATCH 27/37] AnalysisDUT: add histograms --- src/modules/AnalysisDUT/AnalysisDUT.cpp | 27 +++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index bc4477b8..6336e6ac 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -930,6 +930,33 @@ StatusCode AnalysisDUT::run(const std::shared_ptr& clipboard) { } } + if(assoc_cluster->size() > 1) { + for(auto& px : assoc_cluster->pixels()) { + if(px == assoc_cluster->getSeedPixel()) { + continue; // don't fill this histogram for seed pixel! + } + pxTimeMinusSeedTime->Fill( + static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns"))); + pxTimeMinusSeedTime_vs_pxCharge->Fill( + static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + px->charge()); + + if(assoc_cluster->size() == 2) { + pxTimeMinusSeedTime_vs_pxCharge_2px->Fill( + static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + px->charge()); + } else if(assoc_cluster->size() == 3) { + pxTimeMinusSeedTime_vs_pxCharge_3px->Fill( + static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + px->charge()); + } else if(assoc_cluster->size() == 4) { + pxTimeMinusSeedTime_vs_pxCharge_4px->Fill( + static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + px->charge()); + } + } + } + // mean cluster size npxvsxmym->Fill(xmod_um, ymod_um, static_cast(assoc_cluster->size())); if(assoc_cluster->size() == 1) { -- GitLab From de7fcc6fde730ab6cb752f9a2fe5c0019704f6a4 Mon Sep 17 00:00:00 2001 From: Lennart Huth Date: Fri, 23 Apr 2021 16:00:47 +0200 Subject: [PATCH 28/37] reverted eigen3 change --- 3rdparty/GeneralBrokenLines/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/3rdparty/GeneralBrokenLines/CMakeLists.txt b/3rdparty/GeneralBrokenLines/CMakeLists.txt index 6f9d8d2a..3c9e535d 100644 --- a/3rdparty/GeneralBrokenLines/CMakeLists.txt +++ b/3rdparty/GeneralBrokenLines/CMakeLists.txt @@ -11,6 +11,7 @@ PROJECT(GBL VERSION "2.2.0") # Enable ROOT support FIND_PACKAGE(ROOT REQUIRED) +FIND_PACKAGE(EIGEN3 REQUIRED) ADD_DEFINITIONS(-DGBL_EIGEN_SUPPORT_ROOT) # include directories -- GitLab From 1dd7569d4bc1e26eb51108708fd5f8de131e2263 Mon Sep 17 00:00:00 2001 From: Lennart Huth Date: Fri, 23 Apr 2021 16:08:07 +0200 Subject: [PATCH 29/37] eigen3 typoe --- 3rdparty/GeneralBrokenLines/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/3rdparty/GeneralBrokenLines/CMakeLists.txt b/3rdparty/GeneralBrokenLines/CMakeLists.txt index 3c9e535d..1e55a345 100644 --- a/3rdparty/GeneralBrokenLines/CMakeLists.txt +++ b/3rdparty/GeneralBrokenLines/CMakeLists.txt @@ -11,8 +11,8 @@ PROJECT(GBL VERSION "2.2.0") # Enable ROOT support FIND_PACKAGE(ROOT REQUIRED) -FIND_PACKAGE(EIGEN3 REQUIRED) ADD_DEFINITIONS(-DGBL_EIGEN_SUPPORT_ROOT) +FIND_PACKAGE(Eigen3 REQUIRED) # include directories INCLUDE_DIRECTORIES(SYSTEM BEFORE ./include ${ROOT_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIR} ) -- GitLab From 450bf211d9efa8c60c4a59baaf79b9de23af1840 Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Fri, 23 Apr 2021 17:11:37 +0200 Subject: [PATCH 30/37] AnalysisDUT: correctly fill histograms (seed px time instead of cluster time) --- src/modules/AnalysisDUT/AnalysisDUT.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index 6336e6ac..6144f6f8 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -935,23 +935,27 @@ StatusCode AnalysisDUT::run(const std::shared_ptr& clipboard) { if(px == assoc_cluster->getSeedPixel()) { continue; // don't fill this histogram for seed pixel! } - pxTimeMinusSeedTime->Fill( - static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns"))); + pxTimeMinusSeedTime->Fill(static_cast( + Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns"))); pxTimeMinusSeedTime_vs_pxCharge->Fill( - static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + static_cast( + Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), px->charge()); if(assoc_cluster->size() == 2) { pxTimeMinusSeedTime_vs_pxCharge_2px->Fill( - static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + static_cast( + Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), px->charge()); } else if(assoc_cluster->size() == 3) { pxTimeMinusSeedTime_vs_pxCharge_3px->Fill( - static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + static_cast( + Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), px->charge()); } else if(assoc_cluster->size() == 4) { pxTimeMinusSeedTime_vs_pxCharge_4px->Fill( - static_cast(Units::convert(px->timestamp() - assoc_cluster->timestamp(), "ns")), + static_cast( + Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), px->charge()); } } -- GitLab From 8038eb54003b8f8d83d0b94d7d0a4637677a47af Mon Sep 17 00:00:00 2001 From: Simon Spannagel Date: Thu, 29 Apr 2021 16:37:57 +0200 Subject: [PATCH 31/37] Clipboard: make the data removal work properly --- src/core/clipboard/Clipboard.hpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/core/clipboard/Clipboard.hpp b/src/core/clipboard/Clipboard.hpp index 94c5f31e..a07bbcb3 100644 --- a/src/core/clipboard/Clipboard.hpp +++ b/src/core/clipboard/Clipboard.hpp @@ -178,16 +178,6 @@ namespace corryvreckan { bool append = false); /** -<<<<<<< HEAD -======= - * Helper to remove a single object from the clipboard - * @param data Existing data the object should be removed from - * @param object Object to be removed - */ - template void remove_data(std::vector>& data, const std::shared_ptr& object); - - /** ->>>>>>> Clipboard: add methods to remove data from the clipboard again * Helper to remove a set of objects from the clipboard * @param storage_element Storage element concerned * @param objects Objects to be removed -- GitLab From d6471c8d05dfe813fd6a9e46cc8518766f40ef1d Mon Sep 17 00:00:00 2001 From: Jens Kroeger Date: Mon, 17 May 2021 14:08:25 +0200 Subject: [PATCH 32/37] AnalysisDUT: only create track correlation plots when configured (+updated README) --- src/modules/AnalysisDUT/AnalysisDUT.cpp | 72 ------------------------- 1 file changed, 72 deletions(-) diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index 6144f6f8..d99556d7 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -303,47 +303,6 @@ void AnalysisDUT::initialize() { 499.5); residualsXprofile = new TProfile2D("residualsXprofile", "Residual in X;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); residualsYprofile = new TProfile2D("residualsYprofile", "Residual in Y;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); -<<<<<<< HEAD -======= - residualsY = new TH1F("residualsY", "Residual in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsYvsXhit = new TH2F( - "residualsYvsXhit", "Residual in Y;y_{track}-y_{hit} [mm];x_{hit} [mm];# entries", 800, -0.1, 0.1, 40, -20, 20); - residualsYvsYhit = new TH2F( - "residualsYvsYhit", "Residual in Y;y_{track}-y_{hit} [mm];y_{hit} [mm];# entries", 800, -0.1, 0.1, 40, -20, 20); - residualsZ = new TH1F("residualsZ", "Residual in Z;z_{track}-z_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsPos = new TH1F( - "residualsPos", "Absolute distance between track and hit;|pos_{track}-pos_{hit}| [mm];# entries", 800, -0.1, 0.1); - residualsPosVsresidualsTime = - new TH2F("residualsPosVsresidualsTime", - "Time vs. absolute position residuals;time_{track}-time_{hit} [ns];|pos_{track}-pos_{hit}| [mm];# entries", - n_timebins_, - -n_timebins_ / 2. * time_binning_ - time_binning_ / 2., - n_timebins_ / 2. * time_binning_ - time_binning_ / 2., - 800, - 0., - 0.2); - - residualsX1pix = - new TH1F("residualsX1pix", "Residual for 1-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsY1pix = - new TH1F("residualsY1pix", "Residual for 1-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsX2pix = - new TH1F("residualsX2pix", "Residual for 2-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsY2pix = - new TH1F("residualsY2pix", "Residual for 2-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsX3pix = - new TH1F("residualsX3pix", "Residual for 3-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsY3pix = - new TH1F("residualsY3pix", "Residual for 3-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsX4pix = - new TH1F("residualsX4pix", "Residual for 4-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsY4pix = - new TH1F("residualsY4pix", "Residual for 4-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsXatLeast5pix = new TH1F( - "residualsXatLeast5pix", "Residual for >= 5-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1); - residualsYatLeast5pix = new TH1F( - "residualsYatLeast5pix", "Residual for >= 5-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1); ->>>>>>> AnalysisDUT: unify naming convention of member variables clusterChargeAssoc = new TH1F("clusterChargeAssociated", "Charge distribution of associated clusters;cluster charge [e];# entries", @@ -930,37 +889,6 @@ StatusCode AnalysisDUT::run(const std::shared_ptr& clipboard) { } } - if(assoc_cluster->size() > 1) { - for(auto& px : assoc_cluster->pixels()) { - if(px == assoc_cluster->getSeedPixel()) { - continue; // don't fill this histogram for seed pixel! - } - pxTimeMinusSeedTime->Fill(static_cast( - Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns"))); - pxTimeMinusSeedTime_vs_pxCharge->Fill( - static_cast( - Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), - px->charge()); - - if(assoc_cluster->size() == 2) { - pxTimeMinusSeedTime_vs_pxCharge_2px->Fill( - static_cast( - Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), - px->charge()); - } else if(assoc_cluster->size() == 3) { - pxTimeMinusSeedTime_vs_pxCharge_3px->Fill( - static_cast( - Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), - px->charge()); - } else if(assoc_cluster->size() == 4) { - pxTimeMinusSeedTime_vs_pxCharge_4px->Fill( - static_cast( - Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")), - px->charge()); - } - } - } - // mean cluster size npxvsxmym->Fill(xmod_um, ymod_um, static_cast(assoc_cluster->size())); if(assoc_cluster->size() == 1) { -- GitLab From d5f1f0f8ef77e53ed3a930d586f4aa969d3caa68 Mon Sep 17 00:00:00 2001 From: p-becht <> Date: Wed, 3 Mar 2021 10:52:03 +0100 Subject: [PATCH 33/37] copied from my personal fork, shearrz constraint in Alignment Millipede applied, general implementation of bending direction and Intercept that even should allow for an angled track --- src/core/detector/BentPixelDetector.cpp | 61 +++++++++++++++++++ src/core/detector/BentPixelDetector.hpp | 8 ++- src/core/detector/Detector.cpp | 2 +- src/core/detector/Detector.hpp | 4 +- .../AlignmentMillepede/AlignmentMillepede.cpp | 4 +- src/modules/AnalysisDUT/AnalysisDUT.h | 1 + 6 files changed, 72 insertions(+), 8 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 4c8c5097..eee3885b 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -473,6 +473,65 @@ void BentPixelDetector::getInterceptParameters(const PositionVector3D>& state_track, + const DisplacementVector3D>& direction_track, + const PositionVector3D>& state_cylinder, + const DisplacementVector3D>& direction_cylinder, + double (&result)[2]) const { + // Build quadratic equation vectors + double gamma = direction_cylinder.x() * direction_cylinder.x() + direction_cylinder.y() * direction_cylinder.y() + + direction_cylinder.z() * direction_cylinder.z(); // in case direction_cylinder is not normalised + + double alpha = direction_track.x() * direction_cylinder.x() + direction_track.y() * direction_cylinder.y() + + direction_track.z() * direction_cylinder.z(); // direction_track dot direction_cylinder + PositionVector3D> alpha_vec(direction_track.x() - (alpha / gamma) * direction_cylinder.x(), + direction_track.y() - (alpha / gamma) * direction_cylinder.y(), + direction_track.z() - (alpha / gamma) * direction_cylinder.z()); + + alpha = state_track.x() * direction_cylinder.x() + state_track.y() * direction_cylinder.y() + + state_track.z() * direction_cylinder.z(); + double beta = state_cylinder.x() * direction_cylinder.x() + state_cylinder.y() * direction_cylinder.y() + + state_cylinder.z() * direction_cylinder.z(); + PositionVector3D> gamma_vec( + state_track.x() - state_cylinder.x() + (alpha / gamma) * direction_cylinder.x() - + (beta / gamma) * direction_cylinder.x(), + state_track.y() - state_cylinder.y() + (alpha / gamma) * direction_cylinder.y() - + (beta / gamma) * direction_cylinder.y(), + state_track.z() - state_cylinder.z() + (alpha / gamma) * direction_cylinder.z() - + (beta / gamma) * direction_cylinder.z()); + + // Get factors for quadratic equation: alpha z^2 + beta * z + gamma + alpha = alpha_vec.x() * alpha_vec.x() + alpha_vec.y() * alpha_vec.y() + + alpha_vec.z() * alpha_vec.z(); // alpha_vec dot alpha_vec + beta = 2 * (alpha_vec.x() * gamma_vec.x() + alpha_vec.y() * gamma_vec.y() + + alpha_vec.z() * gamma_vec.z()); // 2 * alpha_vec dot gamma_vec + gamma = (gamma_vec.x() * gamma_vec.x() + gamma_vec.y() * gamma_vec.y() + gamma_vec.z() * gamma_vec.z()) - + (m_radius * m_radius); // gamma_vec dot gamma_vec - radius^2 + // previous implementation (not general, depends on cylinder orientation) + // double gamma = state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * + // direction_tr.z()) - + // 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() + state_tr.y() * state_tr.y(); + // double beta = -2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + + // 2 * state_tr.y() * direction_tr.y() / direction_tr.z() + 2 * m_radius; + // double alpha = 1 + direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()); + LOG(TRACE) << "Built quadratic equation to calculate intercept parameter"; + + // Check if the quadratic equation has a solution + double radical = beta * beta - 4 * alpha * gamma; + if(radical < 0 && alpha != 0) { // return a value outside the chip for the track to have no intercept + result[0] = -999999; + result[1] = -999999; + return; + } + + // Solve equation + result[0] = (-beta + sqrt(radical)) / (2 * alpha); + result[1] = (-beta - sqrt(radical)) / (2 * alpha); + // sort in ascending order + std::sort(std::begin(result), std::end(result)); + return; +} + PositionVector3D> BentPixelDetector::getLocalIntercept(const Track* track) const { return globalToLocal(getIntercept(track)); } @@ -526,6 +585,8 @@ bool BentPixelDetector::hitMasked(const Track* track, int tolerance) const { } // Functions to get row and column from local position +// FOLLOW UP: +// might need to be modified taking into account the bent length (archlenght). Seems fine! double BentPixelDetector::getRow(const PositionVector3D> localPosition) const { return localPosition.Y() / m_pitch.Y() + static_cast(m_nPixels.Y() - 1) / 2.; } diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index 99c6c62a..2dd1188a 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -115,19 +115,21 @@ namespace corryvreckan { // Function to get local position from column (x) and row (y) coordinates PositionVector3D> getLocalPosition(double column, double row) const override; - /** + /** * @brief Transform local coordinates of this detector into global coordinates * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - XYZPoint localToGlobal(XYZPoint local) const override; + //XYZPoint localToGlobal(XYZPoint local) const override; + XYZPoint localToGlobal(XYZPoint local) const; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - XYZPoint globalToLocal(XYZPoint global) const override; + //XYZPoint globalToLocal(XYZPoint global) const override; + XYZPoint globalToLocal(XYZPoint global) const; /** * Transformation from local (sensor) coordinates to in-pixel coordinates diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index 535707d7..59bf157d 100644 --- a/src/core/detector/Detector.cpp +++ b/src/core/detector/Detector.cpp @@ -15,8 +15,8 @@ #include "Math/RotationZ.h" #include "Math/RotationZYX.h" -#include "BentPixelDetector.hpp" #include "Detector.hpp" +#include "BentPixelDetector.hpp" #include "PixelDetector.hpp" #include "core/utils/log.h" #include "exceptions.h" diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index 68cc69d9..3ad0fbf3 100644 --- a/src/core/detector/Detector.hpp +++ b/src/core/detector/Detector.hpp @@ -280,14 +280,14 @@ namespace corryvreckan { * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - virtual XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; + XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - virtual XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; + XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; /** * @brief Check whether given track is within the detector's region-of-interest diff --git a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp index 6de0924d..fec24e0d 100644 --- a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp +++ b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp @@ -218,7 +218,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { std::vector fscaz(nParameters, 0.); std::vector shearx(nParameters, 0.); std::vector sheary(nParameters, 0.); - std::vector shearrz(nParameters, 0.); + std::vector shearrz(nParameters, 0.); m_constraints.clear(); for(const auto& det : get_detectors()) { @@ -240,7 +240,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { shearx[i] = sz; sheary[i + nPlanes] = sz; fscaz[i + 2 * nPlanes] = sz; - shearrz[i + 5 * nPlanes] = sz; + shearrz[i + 5 * nPlanes] = sz; } const std::vector constraints = {true, true, true, true, false, false, true, true, true, true}; diff --git a/src/modules/AnalysisDUT/AnalysisDUT.h b/src/modules/AnalysisDUT/AnalysisDUT.h index 528bc692..57382991 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -58,6 +58,7 @@ namespace corryvreckan { TH2F *residualsXvsXhit, *residualsXvsYhit; TH2F *residualsYvsXhit, *residualsYvsYhit; TH2F* residualsXvsresidualsY; + TProfile2D* residualsXprofile; TProfile2D* residualsYprofile; -- GitLab From 820d6b820862551a9c81e8fea3089203dfdbe691 Mon Sep 17 00:00:00 2001 From: Mihail Bogdan Blidaru Date: Wed, 3 Mar 2021 12:19:19 +0100 Subject: [PATCH 34/37] Clean-up y_0. --- src/core/detector/BentPixelDetector.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index 2dd1188a..8877c73d 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -226,6 +226,7 @@ namespace corryvreckan { // For bent pixel detector XYVector m_pitch{}; + double m_radius; XYVector m_spatial_resolution{}; ROOT::Math::DisplacementVector2D> m_nPixels{}; std::vector> m_roi{}; -- GitLab From bb30b6b7f8573d9df9512afb97618cf99b8e1445 Mon Sep 17 00:00:00 2001 From: p-becht <> Date: Thu, 4 Mar 2021 13:10:23 +0100 Subject: [PATCH 35/37] ran clang-format to fix formatting errors in pipeline --- src/core/detector/BentPixelDetector.hpp | 10 +++++----- src/core/detector/Detector.cpp | 2 +- src/modules/AlignmentMillepede/AlignmentMillepede.cpp | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index 8877c73d..6a72b0c6 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -115,21 +115,21 @@ namespace corryvreckan { // Function to get local position from column (x) and row (y) coordinates PositionVector3D> getLocalPosition(double column, double row) const override; - /** + /** * @brief Transform local coordinates of this detector into global coordinates * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - //XYZPoint localToGlobal(XYZPoint local) const override; - XYZPoint localToGlobal(XYZPoint local) const; + // XYZPoint localToGlobal(XYZPoint local) const override; + XYZPoint localToGlobal(XYZPoint local) const; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - //XYZPoint globalToLocal(XYZPoint global) const override; - XYZPoint globalToLocal(XYZPoint global) const; + // XYZPoint globalToLocal(XYZPoint global) const override; + XYZPoint globalToLocal(XYZPoint global) const; /** * Transformation from local (sensor) coordinates to in-pixel coordinates diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index 59bf157d..535707d7 100644 --- a/src/core/detector/Detector.cpp +++ b/src/core/detector/Detector.cpp @@ -15,8 +15,8 @@ #include "Math/RotationZ.h" #include "Math/RotationZYX.h" -#include "Detector.hpp" #include "BentPixelDetector.hpp" +#include "Detector.hpp" #include "PixelDetector.hpp" #include "core/utils/log.h" #include "exceptions.h" diff --git a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp index fec24e0d..6de0924d 100644 --- a/src/modules/AlignmentMillepede/AlignmentMillepede.cpp +++ b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp @@ -218,7 +218,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { std::vector fscaz(nParameters, 0.); std::vector shearx(nParameters, 0.); std::vector sheary(nParameters, 0.); - std::vector shearrz(nParameters, 0.); + std::vector shearrz(nParameters, 0.); m_constraints.clear(); for(const auto& det : get_detectors()) { @@ -240,7 +240,7 @@ void AlignmentMillepede::setConstraints(const size_t nPlanes) { shearx[i] = sz; sheary[i + nPlanes] = sz; fscaz[i + 2 * nPlanes] = sz; - shearrz[i + 5 * nPlanes] = sz; + shearrz[i + 5 * nPlanes] = sz; } const std::vector constraints = {true, true, true, true, false, false, true, true, true, true}; -- GitLab From f2f4d2c83f29ea702fdde573d4456b367fd67719 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Fri, 28 May 2021 11:10:57 +0200 Subject: [PATCH 36/37] Align fixes --- src/core/detector/BentPixelDetector.cpp | 61 ------------------------- src/core/detector/BentPixelDetector.hpp | 7 +-- src/core/detector/Detector.hpp | 4 +- src/modules/AnalysisDUT/AnalysisDUT.cpp | 1 + 4 files changed, 5 insertions(+), 68 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index eee3885b..4c8c5097 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -473,65 +473,6 @@ void BentPixelDetector::getInterceptParameters(const PositionVector3D>& state_track, - const DisplacementVector3D>& direction_track, - const PositionVector3D>& state_cylinder, - const DisplacementVector3D>& direction_cylinder, - double (&result)[2]) const { - // Build quadratic equation vectors - double gamma = direction_cylinder.x() * direction_cylinder.x() + direction_cylinder.y() * direction_cylinder.y() + - direction_cylinder.z() * direction_cylinder.z(); // in case direction_cylinder is not normalised - - double alpha = direction_track.x() * direction_cylinder.x() + direction_track.y() * direction_cylinder.y() + - direction_track.z() * direction_cylinder.z(); // direction_track dot direction_cylinder - PositionVector3D> alpha_vec(direction_track.x() - (alpha / gamma) * direction_cylinder.x(), - direction_track.y() - (alpha / gamma) * direction_cylinder.y(), - direction_track.z() - (alpha / gamma) * direction_cylinder.z()); - - alpha = state_track.x() * direction_cylinder.x() + state_track.y() * direction_cylinder.y() + - state_track.z() * direction_cylinder.z(); - double beta = state_cylinder.x() * direction_cylinder.x() + state_cylinder.y() * direction_cylinder.y() + - state_cylinder.z() * direction_cylinder.z(); - PositionVector3D> gamma_vec( - state_track.x() - state_cylinder.x() + (alpha / gamma) * direction_cylinder.x() - - (beta / gamma) * direction_cylinder.x(), - state_track.y() - state_cylinder.y() + (alpha / gamma) * direction_cylinder.y() - - (beta / gamma) * direction_cylinder.y(), - state_track.z() - state_cylinder.z() + (alpha / gamma) * direction_cylinder.z() - - (beta / gamma) * direction_cylinder.z()); - - // Get factors for quadratic equation: alpha z^2 + beta * z + gamma - alpha = alpha_vec.x() * alpha_vec.x() + alpha_vec.y() * alpha_vec.y() + - alpha_vec.z() * alpha_vec.z(); // alpha_vec dot alpha_vec - beta = 2 * (alpha_vec.x() * gamma_vec.x() + alpha_vec.y() * gamma_vec.y() + - alpha_vec.z() * gamma_vec.z()); // 2 * alpha_vec dot gamma_vec - gamma = (gamma_vec.x() * gamma_vec.x() + gamma_vec.y() * gamma_vec.y() + gamma_vec.z() * gamma_vec.z()) - - (m_radius * m_radius); // gamma_vec dot gamma_vec - radius^2 - // previous implementation (not general, depends on cylinder orientation) - // double gamma = state_tr.z() * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * - // direction_tr.z()) - - // 2 * state_tr.z() * state_tr.y() * direction_tr.y() / direction_tr.z() + state_tr.y() * state_tr.y(); - // double beta = -2 * state_tr.z() * direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()) + - // 2 * state_tr.y() * direction_tr.y() / direction_tr.z() + 2 * m_radius; - // double alpha = 1 + direction_tr.y() * direction_tr.y() / (direction_tr.z() * direction_tr.z()); - LOG(TRACE) << "Built quadratic equation to calculate intercept parameter"; - - // Check if the quadratic equation has a solution - double radical = beta * beta - 4 * alpha * gamma; - if(radical < 0 && alpha != 0) { // return a value outside the chip for the track to have no intercept - result[0] = -999999; - result[1] = -999999; - return; - } - - // Solve equation - result[0] = (-beta + sqrt(radical)) / (2 * alpha); - result[1] = (-beta - sqrt(radical)) / (2 * alpha); - // sort in ascending order - std::sort(std::begin(result), std::end(result)); - return; -} - PositionVector3D> BentPixelDetector::getLocalIntercept(const Track* track) const { return globalToLocal(getIntercept(track)); } @@ -585,8 +526,6 @@ bool BentPixelDetector::hitMasked(const Track* track, int tolerance) const { } // Functions to get row and column from local position -// FOLLOW UP: -// might need to be modified taking into account the bent length (archlenght). Seems fine! double BentPixelDetector::getRow(const PositionVector3D> localPosition) const { return localPosition.Y() / m_pitch.Y() + static_cast(m_nPixels.Y() - 1) / 2.; } diff --git a/src/core/detector/BentPixelDetector.hpp b/src/core/detector/BentPixelDetector.hpp index 6a72b0c6..99c6c62a 100644 --- a/src/core/detector/BentPixelDetector.hpp +++ b/src/core/detector/BentPixelDetector.hpp @@ -120,16 +120,14 @@ namespace corryvreckan { * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - // XYZPoint localToGlobal(XYZPoint local) const override; - XYZPoint localToGlobal(XYZPoint local) const; + XYZPoint localToGlobal(XYZPoint local) const override; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - // XYZPoint globalToLocal(XYZPoint global) const override; - XYZPoint globalToLocal(XYZPoint global) const; + XYZPoint globalToLocal(XYZPoint global) const override; /** * Transformation from local (sensor) coordinates to in-pixel coordinates @@ -226,7 +224,6 @@ namespace corryvreckan { // For bent pixel detector XYVector m_pitch{}; - double m_radius; XYVector m_spatial_resolution{}; ROOT::Math::DisplacementVector2D> m_nPixels{}; std::vector> m_roi{}; diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index 3ad0fbf3..68cc69d9 100644 --- a/src/core/detector/Detector.hpp +++ b/src/core/detector/Detector.hpp @@ -280,14 +280,14 @@ namespace corryvreckan { * @param local Local coordinates in the reference frame of this detector * @return Global coordinates */ - XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; + virtual XYZPoint localToGlobal(XYZPoint local) const { return m_localToGlobal * local; }; /** * @brief Transform global coordinates into detector-local coordinates * @param global Global coordinates * @return Local coordinates in the reference frame of this detector */ - XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; + virtual XYZPoint globalToLocal(XYZPoint global) const { return m_globalToLocal * global; }; /** * @brief Check whether given track is within the detector's region-of-interest diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index d99556d7..718a3e68 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -302,6 +302,7 @@ void AnalysisDUT::initialize() { -500.5, 499.5); residualsXprofile = new TProfile2D("residualsXprofile", "Residual in X;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); + residualsYprofile = new TProfile2D("residualsYprofile", "Residual in Y;Row;Column;", 128, 0, 1024, 64, 0, 512, "s"); clusterChargeAssoc = new TH1F("clusterChargeAssociated", -- GitLab From 79b17d4b1a5ac26c3352de0eda53f0ca917d4bd0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Fri, 8 Oct 2021 13:44:15 +0200 Subject: [PATCH 37/37] Align to PixelDetector --- src/core/detector/BentPixelDetector.cpp | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/core/detector/BentPixelDetector.cpp b/src/core/detector/BentPixelDetector.cpp index 4c8c5097..fcd19633 100644 --- a/src/core/detector/BentPixelDetector.cpp +++ b/src/core/detector/BentPixelDetector.cpp @@ -81,6 +81,7 @@ void BentPixelDetector::build_axes(const Configuration& config) { // region of interest: m_roi = config.getMatrix("roi", std::vector>()); + m_maskfile_name = ""; if(config.has("mask_file")) { m_maskfile_name = config.get("mask_file"); std::string mask_file = config.getPath("mask_file", true); @@ -100,7 +101,7 @@ void BentPixelDetector::SetPostionAndOrientation(const Configuration& config) { throw InvalidValueError(config, "orientation_mode", "orientation_mode should be either 'zyx', xyz' or 'zxz'"); } - LOG(WARNING) << "BentPixelDetector!"; + LOG(TRACE) << "BentPixelDetector!"; LOG(TRACE) << " Position: " << Units::display(m_displacement, {"mm", "um"}); LOG(TRACE) << " Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")"; LOG(TRACE) << " m_displacement: X=" << m_displacement.X() << " Y=" << m_displacement.Y() << " Z=" << m_displacement.Z(); @@ -208,14 +209,18 @@ void BentPixelDetector::configure_detector(Configuration& config) const { config.set("number_of_pixels", m_nPixels); // Size of the pixels - config.set("pixel_pitch", m_pitch, {"um"}); + config.set("pixel_pitch", m_pitch, {{"um"}}); // Intrinsic resolution: - config.set("spatial_resolution", m_spatial_resolution, {"um"}); + config.set("spatial_resolution", m_spatial_resolution, {{"um"}}); // Pixel mask file: - if(!m_maskfile_name.empty()) { - config.set("mask_file", m_maskfile_name); + if(!m_maskfile.empty()) { + if(m_maskfile_name.empty()) { + config.set("mask_file", m_maskfile); + } else { + config.set("mask_file", m_maskfile_name); + } } // Region-of-interest: @@ -223,12 +228,12 @@ void BentPixelDetector::configure_detector(Configuration& config) const { } void BentPixelDetector::configure_pos_and_orientation(Configuration& config) const { - config.set("position", m_displacement, {"um", "mm"}); + config.set("position", m_displacement, {{"um", "mm"}}); config.set("orientation_mode", m_orientation_mode); - config.set("orientation", m_orientation, {"deg"}); + config.set("orientation", m_orientation, {{"deg"}}); config.set("coordinates", m_coordinates); - config.set("radius", m_radius, {"mm"}); - config.set("flat_part", m_flat_part, {"mm"}); + config.set("radius", m_radius, {{"mm"}}); + config.set("flat_part", m_flat_part, {{"mm"}}); config.set("bent_axis", m_bent_axis); } -- GitLab