diff --git a/3rdparty/GeneralBrokenLines/include/BorderedBandMatrix.h b/3rdparty/GeneralBrokenLines/include/BorderedBandMatrix.h index dacb11202ebebe6a937aad1bc0142c12ce29754c..61a79fffa6508a864e443b7b30c765f86707a677 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 8f7842880c5f05c7c59500fdf7e8a7695e2e1b1f..4e12fee7144ee5a8ed0a25bb9a8ccfb79654432e 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 5bbbc17d92b95580fc4b1ac07ae0ce9fbdd79594..ad25c662ea2599a0e546c09f893d926202a78508 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 3de9aba2934c264205cc8f558cb7a29d8a5b4849..5dbb07ffb7a74f7affe2b4713d43f4ca5437c5dd 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 aa350acfd32df3b76c6b9861830b0e7b50d37593..8833fb71f6af5ba6ff7dba2e9833c0d10d2ebeaa 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 d7fe7295494bdd46ce219d756cd634ef6c96bd69..5dcd273d1e8426cd20e89a4bf881f347e54e30e8 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 724cbdec13a6bb69d1480ad36c4a493470786fd1..bcf52a24d49be3dd7fb9c5fc863c9729ad0bf731 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 6fee613507b8455dc717bf382c99761eb7f27908..8e19d9ac79e2fb6f24d81a5f97c254c0fb2239e4 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 fd137fac20ed420334113c5e55fe1c4b585ffa66..3b839cb2f6fa60e3f08d2b02a822820bb7b02d22 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 a904c7bbc3a04102460942a57669be3e89814734..afe91656568cc9df2d87c2d1f759352734fe9435 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 39d6a94dc1db356b5743be3c4aea97c26257e726..295e64d0dcebde202112fbd7b21a8dd00db9dced 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 7b2ec569d7c477fc2344202f453211f2f7d1d8c6..6b21b01593874032c611f1133a0ee3597a9c062e 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/CMakeLists.txt b/src/core/CMakeLists.txt index 96815fae7941d3b27ccfe9f21a082763de0fef1c..61edb3a8511513b08e195326f5cd10c83e4feeb6 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 new file mode 100644 index 0000000000000000000000000000000000000000..fcd196338e1e0eaca0a9cb34749f46088b64e091 --- /dev/null +++ b/src/core/detector/BentPixelDetector.cpp @@ -0,0 +1,678 @@ +/** @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/DisplacementVector2D.h" +#include "Math/PositionVector3D.h" +#include "Math/RotationX.h" +#include "Math/RotationY.h" +#include "Math/RotationZ.h" +#include "Math/RotationZYX.h" + +#include "BentPixelDetector.hpp" +#include "core/utils/log.h" +#include "exceptions.h" + +using namespace ROOT::Math; +using namespace corryvreckan; + +BentPixelDetector::BentPixelDetector(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 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_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) { + 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>()); + + m_maskfile_name = ""; + 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; + maskFile(mask_file); + process_mask_file(); + } +} + +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()); + m_orientation_mode = config.get("orientation_mode", "xyz"); + + 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(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(); + LOG(TRACE) << " m_orientation: X=" << m_orientation.X() << " Y=" << m_orientation.Y() << " Z=" << m_orientation.Z(); +} + +void BentPixelDetector::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 BentPixelDetector::maskChannel(int chX, int chY) { + int channelID = chX + m_nPixels.X() * chY; + m_masked[channelID] = true; +} + +bool BentPixelDetector::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 BentPixelDetector::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 { + throw InvalidSettingError(this, "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 BentPixelDetector::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.empty()) { + if(m_maskfile_name.empty()) { + config.set("mask_file", m_maskfile); + } else { + config.set("mask_file", m_maskfile_name); + } + } + + // Region-of-interest: + config.setMatrix("roi", m_roi); +} + +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"}}); + config.set("coordinates", m_coordinates); + config.set("radius", m_radius, {{"mm"}}); + config.set("flat_part", m_flat_part, {{"mm"}}); + config.set("bent_axis", m_bent_axis); +} + +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::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; + } +} + +PositionVector3D> BentPixelDetector::getIntercept(const Track* track) const { + // FIXME: this 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> 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()); + + // 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; + } + + // 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()))); + } 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 { + return globalToLocal(getIntercept(track)); +} + +// Function to check if a track intercepts with a plane +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); + + // Convert to local coordinates + PositionVector3D> localIntercept = 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 BentPixelDetector::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 = 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 BentPixelDetector::getRow(const PositionVector3D> localPosition) const { + return localPosition.Y() / m_pitch.Y() + static_cast(m_nPixels.Y() - 1) / 2.; +} +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> 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.), + 0.); +} + +// Function to get in-pixel position +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 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 BentPixelDetector::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 BentPixelDetector::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 BentPixelDetector::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 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."; + 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 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)); +} + +// 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 new file mode 100644 index 0000000000000000000000000000000000000000..99c6c62a1d8228c8aa39232f9bd0fc72fed32a0a --- /dev/null +++ b/src/core/detector/BentPixelDetector.hpp @@ -0,0 +1,242 @@ +/** @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_BENTPIXELDETECTOR_H +#define CORRYVRECKAN_BENTPIXELDETECTOR_H + +#include +#include +#include + +#include "Math/DisplacementVector2D.h" +#include "Math/Transform3D.h" +#include "Math/Vector2D.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/Cluster.hpp" +#include "objects/Pixel.hpp" +#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 BentPixelDetector : public Detector { + public: + /** + * Delete default constructor + */ + BentPixelDetector() = delete; + + /** + * Default destructor + */ + ~BentPixelDetector() = default; + + /** + * @brief Constructs a detector in the geometry + * @param config Configuration object describing the detector + */ + BentPixelDetector(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; + + /** + * @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 + * @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; } + + /** + * @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; + + // 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); + + // 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{}; + 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; + // bent geometry configuration parameters + double m_radius; + double m_flat_part; + std::string m_bent_axis; + std::string m_coordinates; + }; +} // namespace corryvreckan + +#endif // CORRYVRECKAN_BENTPIXELDETECTOR_H diff --git a/src/core/detector/Detector.cpp b/src/core/detector/Detector.cpp index 8d1edef67d3d22c30a2a17677f2b3ded980ed9a6..535707d72fbfdc37b57e341751aee51f12d4ebbc 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" @@ -84,8 +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") { + return std::make_shared(config); } else { - throw InvalidValueError(config, "coordinates", "Coordinates can only set to be cartesian now"); + throw InvalidValueError( + config, "coordinates", "Coordinates can only be set to \"cartesian\" or \"cartesian-bent\" for now"); } } diff --git a/src/core/detector/Detector.hpp b/src/core/detector/Detector.hpp index 3ad0fbf394662878e4c480be02ee346594af28d0..68cc69d9cc7d2e589eb6dfb4f0319dcf562a94e1 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/AlignmentMillepede/AlignmentMillepede.cpp b/src/modules/AlignmentMillepede/AlignmentMillepede.cpp index c0645914b91232347e42fb4e2dc20df7b06ac0a4..6de0924dc23b5dd5d938620804f2be4451ae952f 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,27 +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; } - 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]) + if(constraints[0] && m_dofs[0]) { addConstraint(ftx, 0.0); - if(constraints[1] && m_dofs[0]) + } + if(constraints[1] && m_dofs[0]) { addConstraint(shearx, 0.); - if(constraints[2] && m_dofs[1]) + } + if(constraints[2] && m_dofs[1]) { addConstraint(fty, 0.0); - if(constraints[3] && m_dofs[1]) + } + if(constraints[3] && m_dofs[1]) { addConstraint(sheary, 0.); - if(constraints[4] && m_dofs[2]) + } + 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]) + if(constraints[6] && m_dofs[3]) { addConstraint(frx, 0.0); - if(constraints[7] && m_dofs[4]) + } + if(constraints[7] && m_dofs[4]) { addConstraint(fry, 0.0); - if(constraints[8] && m_dofs[5]) + } + if(constraints[8] && m_dofs[5]) { addConstraint(frz, 0.0); + } + if(constraints[9] && m_dofs[6]) { + addConstraint(shearrz, 0.0); + } } //============================================================================= diff --git a/src/modules/AnalysisDUT/AnalysisDUT.cpp b/src/modules/AnalysisDUT/AnalysisDUT.cpp index c559afe0f49e30a03a24d4848b9f4b0ca1a03bf2..718a3e68eb30e94d9a61bb863b9e5646a1c99737 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.cpp +++ b/src/modules/AnalysisDUT/AnalysisDUT.cpp @@ -260,6 +260,50 @@ 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 +796,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 f2202b5d18675f07aaa2339c04ace5169bbe63ab..57382991b6ebdf0911286af23ef437b206cc5446 100644 --- a/src/modules/AnalysisDUT/AnalysisDUT.h +++ b/src/modules/AnalysisDUT/AnalysisDUT.h @@ -55,6 +55,13 @@ namespace corryvreckan { TH1F* associatedTracksVersusTime; TH1F *residualsX, *residualsY, *residualsPos; + TH2F *residualsXvsXhit, *residualsXvsYhit; + TH2F *residualsYvsXhit, *residualsYvsYhit; + TH2F* residualsXvsresidualsY; + + TProfile2D* residualsXprofile; + TProfile2D* residualsYprofile; + TH2F* residualsPosVsresidualsTime; TH1F *residualsX1pix, *residualsY1pix; diff --git a/src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp b/src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp index 566106a87f7795cb4939038f84a4b80367dacb66..b493e50a9b77e693393ce42b1e20dbe503ebd96c 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 3e15f971e286a985152cfc9e5c9d5c2905627e72..2ca4537124f2fd1986e10dd49d512910765ac496 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 52317ee421404e3ee5e0bd759a2a27686e0b6389..e779f6cacb41ada939a55ce75a01b97d2fcb3435 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 */