Commit b1adf4d0 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Merge branch 'spatial_resolution' into 'master'

spatial resolution

See merge request !193
parents 02ba5a48 c4e455cf
Pipeline #1234234 failed with stages
in 4 minutes and 32 seconds
...@@ -67,7 +67,7 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) { ...@@ -67,7 +67,7 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {
// Auxiliary devices don't have: number_of_pixels, pixel_pitch, spatial_resolution, mask_file, region-of-interest // Auxiliary devices don't have: number_of_pixels, pixel_pitch, spatial_resolution, mask_file, region-of-interest
if(!isAuxiliary()) { if(!isAuxiliary()) {
// Intrinsic spatial resolution, no default: // Intrinsic spatial resolution, no default:
m_resolution = config.get<ROOT::Math::XYVector>("resolution"); m_spatial_resolution = config.get<ROOT::Math::XYVector>("spatial_resolution");
// Number of pixels: // Number of pixels:
m_nPixels = config.get<ROOT::Math::DisplacementVector2D<Cartesian2D<int>>>("number_of_pixels"); m_nPixels = config.get<ROOT::Math::DisplacementVector2D<Cartesian2D<int>>>("number_of_pixels");
// Size of the pixels: // Size of the pixels:
...@@ -142,7 +142,7 @@ void Detector::processMaskFile() { ...@@ -142,7 +142,7 @@ void Detector::processMaskFile() {
// Open the file with masked pixels // Open the file with masked pixels
std::ifstream inputMaskFile(m_maskfile, std::ios::in); std::ifstream inputMaskFile(m_maskfile, std::ios::in);
if(!inputMaskFile.is_open()) { if(!inputMaskFile.is_open()) {
LOG(ERROR) << "Could not open mask file " << m_maskfile; LOG(WARNING) << "Could not open mask file " << m_maskfile;
} else { } else {
int row = 0, col = 0; int row = 0, col = 0;
std::string id; std::string id;
...@@ -151,8 +151,8 @@ void Detector::processMaskFile() { ...@@ -151,8 +151,8 @@ void Detector::processMaskFile() {
if(id == "c") { if(id == "c") {
inputMaskFile >> col; inputMaskFile >> col;
if(col > nPixels().X() - 1) { if(col > nPixels().X() - 1) {
LOG(ERROR) << "Column " << col << " outside of pixel matrix, chip has only " << nPixels().X() LOG(WARNING) << "Column " << col << " outside of pixel matrix, chip has only " << nPixels().X()
<< " columns!"; << " columns!";
} }
LOG(TRACE) << "Masking column " << col; LOG(TRACE) << "Masking column " << col;
for(int r = 0; r < nPixels().Y(); r++) { for(int r = 0; r < nPixels().Y(); r++) {
...@@ -161,7 +161,7 @@ void Detector::processMaskFile() { ...@@ -161,7 +161,7 @@ void Detector::processMaskFile() {
} else if(id == "r") { } else if(id == "r") {
inputMaskFile >> row; inputMaskFile >> row;
if(row > nPixels().Y() - 1) { if(row > nPixels().Y() - 1) {
LOG(ERROR) << "Row " << col << " outside of pixel matrix, chip has only " << nPixels().Y() << " rows!"; LOG(WARNING) << "Row " << col << " outside of pixel matrix, chip has only " << nPixels().Y() << " rows!";
} }
LOG(TRACE) << "Masking row " << row; LOG(TRACE) << "Masking row " << row;
for(int c = 0; c < nPixels().X(); c++) { for(int c = 0; c < nPixels().X(); c++) {
...@@ -170,13 +170,13 @@ void Detector::processMaskFile() { ...@@ -170,13 +170,13 @@ void Detector::processMaskFile() {
} else if(id == "p") { } else if(id == "p") {
inputMaskFile >> col >> row; inputMaskFile >> col >> row;
if(col > nPixels().X() - 1 || row > nPixels().Y() - 1) { if(col > nPixels().X() - 1 || row > nPixels().Y() - 1) {
LOG(ERROR) << "Pixel " << col << " " << row << " outside of pixel matrix, chip has only " LOG(WARNING) << "Pixel " << col << " " << row << " outside of pixel matrix, chip has only "
<< nPixels().X() << " x " << nPixels().Y() << " pixels!"; << nPixels().X() << " x " << nPixels().Y() << " pixels!";
} }
LOG(TRACE) << "Masking pixel " << col << " " << row; LOG(TRACE) << "Masking pixel " << col << " " << row;
maskChannel(col, row); // Flag to mask a pixel maskChannel(col, row); // Flag to mask a pixel
} else { } else {
LOG(ERROR) << "Could not parse mask entry (id \"" << id << "\")"; LOG(WARNING) << "Could not parse mask entry (id \"" << id << "\")";
} }
} }
LOG(INFO) << m_masked.size() << " masked pixels"; LOG(INFO) << m_masked.size() << " masked pixels";
...@@ -271,7 +271,7 @@ Configuration Detector::getConfiguration() const { ...@@ -271,7 +271,7 @@ Configuration Detector::getConfiguration() const {
config.set("pixel_pitch", m_pitch, {"um"}); config.set("pixel_pitch", m_pitch, {"um"});
// Intrinsic resolution: // Intrinsic resolution:
config.set("resolution", m_resolution, {"um"}); config.set("spatial_resolution", m_spatial_resolution, {"um"});
// Pixel mask file: // Pixel mask file:
if(!m_maskfile_name.empty()) { if(!m_maskfile_name.empty()) {
......
...@@ -59,7 +59,8 @@ namespace corryvreckan { ...@@ -59,7 +59,8 @@ namespace corryvreckan {
/** /**
* @brief Detector representation in the reconstruction chain * @brief Detector representation in the reconstruction chain
* *
* Contains the detector with all its properties such as type, name, position and orientation, pitch, resolution etc. * Contains the detector with all its properties such as type, name, position and orientation, pitch, spatial resolution
* etc.
*/ */
class Detector { class Detector {
public: public:
...@@ -123,10 +124,10 @@ namespace corryvreckan { ...@@ -123,10 +124,10 @@ namespace corryvreckan {
XYVector pitch() const { return m_pitch; } XYVector pitch() const { return m_pitch; }
/** /**
* @brief Get intrinsic resolution of the detector * @brief Get intrinsic spatial resolution of the detector
* @return Intrinsic resolution in X and Y * @return Intrinsic spatial resolution in X and Y
*/ */
XYVector resolution() const { return m_resolution; } XYVector getSpatialResolution() const { return m_spatial_resolution; }
/** /**
* @brief Get number of pixels in x and y * @brief Get number of pixels in x and y
...@@ -284,7 +285,7 @@ namespace corryvreckan { ...@@ -284,7 +285,7 @@ namespace corryvreckan {
std::string m_detectorType; std::string m_detectorType;
std::string m_detectorName; std::string m_detectorName;
XYVector m_pitch{}; XYVector m_pitch{};
XYVector m_resolution{}; XYVector m_spatial_resolution{};
ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<int>> m_nPixels{}; ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<int>> m_nPixels{};
double m_timeOffset; double m_timeOffset;
double m_timeResolution; double m_timeResolution;
......
...@@ -239,8 +239,8 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) { ...@@ -239,8 +239,8 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) {
cluster->setRow(row); cluster->setRow(row);
cluster->setCharge(charge); cluster->setCharge(charge);
// Set uncertainty on position from intrinstic detector resolution: // Set uncertainty on position from intrinstic detector spatial resolution:
cluster->setError(m_detector->resolution()); cluster->setError(m_detector->getSpatialResolution());
cluster->setTimestamp(timestamp); cluster->setTimestamp(timestamp);
cluster->setDetectorID(detectorID); cluster->setDetectorID(detectorID);
......
...@@ -231,8 +231,8 @@ void ClusteringSpatial::calculateClusterCentre(Cluster* cluster) { ...@@ -231,8 +231,8 @@ void ClusteringSpatial::calculateClusterCentre(Cluster* cluster) {
cluster->setColumn(column); cluster->setColumn(column);
cluster->setCharge(charge); cluster->setCharge(charge);
// Set uncertainty on position from intrinstic detector resolution: // Set uncertainty on position from intrinstic detector spatial resolution:
cluster->setError(m_detector->resolution()); cluster->setError(m_detector->getSpatialResolution());
cluster->setDetectorID(detectorID); cluster->setDetectorID(detectorID);
cluster->setClusterCentre(positionGlobal); cluster->setClusterCentre(positionGlobal);
......
...@@ -6,6 +6,7 @@ using namespace std; ...@@ -6,6 +6,7 @@ using namespace std;
DUTAssociation::DUTAssociation(Configuration config, std::shared_ptr<Detector> detector) DUTAssociation::DUTAssociation(Configuration config, std::shared_ptr<Detector> detector)
: Module(std::move(config), detector), m_detector(detector) { : Module(std::move(config), detector), m_detector(detector) {
// timing cut, relative (x * time_resolution) or absolute:
if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) { if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
throw InvalidCombinationError( throw InvalidCombinationError(
m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive."); m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
...@@ -14,8 +15,21 @@ DUTAssociation::DUTAssociation(Configuration config, std::shared_ptr<Detector> d ...@@ -14,8 +15,21 @@ DUTAssociation::DUTAssociation(Configuration config, std::shared_ptr<Detector> d
} else { } else {
timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution(); timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution();
} }
spatialCut = m_config.get<XYVector>("spatial_cut", 2 * m_detector->pitch());
// spatial cut, relative (x * spatial_resolution) or absolute:
if(m_config.count({"spatial_cut_rel", "spatial_cut_abs"}) > 1) {
throw InvalidCombinationError(
m_config, {"spatial_cut_rel", "spatial_cut_abs"}, "Absolute and relative spatial cuts are mutually exclusive.");
} else if(m_config.has("spatial_cut_abs")) {
spatialCut = m_config.get<XYVector>("spatial_cut_abs");
} else {
spatialCut = m_config.get<double>("spatial_cut_rel", 3.0) * m_detector->getSpatialResolution();
}
useClusterCentre = m_config.get<bool>("use_cluster_centre", false); useClusterCentre = m_config.get<bool>("use_cluster_centre", false);
LOG(DEBUG) << "time_cut = " << Units::display(timeCut, {"ms", "us", "ns"});
LOG(DEBUG) << "spatial_cut = " << Units::display(spatialCut, {"um", "mm"});
LOG(DEBUG) << "use_cluster_centre = " << useClusterCentre;
} }
void DUTAssociation::initialise() { void DUTAssociation::initialise() {
......
...@@ -9,9 +9,10 @@ This module performs a basic tracking method. ...@@ -9,9 +9,10 @@ This module performs a basic tracking method.
Clusters from the first plane in Z (named the seed plane) are related to clusters close in time on the other detector planes using straight line tracks. The DUT plane can be excluded from the track finding. Clusters from the first plane in Z (named the seed plane) are related to clusters close in time on the other detector planes using straight line tracks. The DUT plane can be excluded from the track finding.
### Parameters ### Parameters
* `time_cut_rel`: Number of standard deviations the `time_resolution` of one detector plane will be multiplied by, either the `time_resolution` of the first plane in Z or the current telescope plane, whichever is largest. This calculated value is then used as the maximum time difference allowed between clusters for association to a track. This allows the time cuts between different planes to be detector appropriate. By default, a relative time cut is applied. Absolute and relative time cuts are mutually exclusive. Defaults to `3.0`. * `time_cut_rel`: Factor by which the `time_resolution` of each detector plane will be multiplied, either the `time_resolution` of the first plane in Z or the current telescope plane, whichever is largest. This calculated value is then used as the maximum time difference allowed between clusters and a track for association to the track. This allows the time cuts between different planes to be detector appropriate. By default, a relative time cut is applied. Absolute and relative time cuts are mutually exclusive. Defaults to `3.0`.
* `time_cut_abs`: Specifies an absolute value for the maximum time difference allowed between clusters for association to the track. Absolute and relative time cuts are mutually exclusive. No default value. * `time_cut_abs`: Specifies an absolute value for the maximum time difference allowed between clusters and a track for association to the track. Absolute and relative time cuts are mutually exclusive. No default value.
* `spatial_cut`: Maximum spatial distance in the XY plane allowed between clusters for association for the telescope planes. Default value is `0.2mm`. * `spatial_cut_rel`: Factor by which the `spatial_resolution` in x and y of each detector plane will be multiplied. These calculated value are defining an ellipse which is then used as the maximum distance in the XY plane allowed between clusters and a track for association to the track. This allows the spatial cuts between different planes to be detector appropriate. By default, a relative spatial cut is applied. Absolute and relative spatial cuts are mutually exclusive. Defaults to `3.0`.
* `spatial_cut_abs`: Specifies a set of absolute value (x and y) which defines an ellipse for the maximum spatial distance in the XY plane between clusters and a track for association to the track. Absolute and relative spatial cuts are mutually exclusive. No default value.
* `min_hits_on_track`: Minium number of associated clusters needed to create a track, equivalent to the minimum number of planes required for each track. Default value is `6`. * `min_hits_on_track`: Minium number of associated clusters needed to create a track, equivalent to the minimum number of planes required for each track. Default value is `6`.
* `exclude_dut`: Boolean to choose if the DUT plane is included in the track finding. Default value is `true`. * `exclude_dut`: Boolean to choose if the DUT plane is included in the track finding. Default value is `true`.
* `require_detectors`: Names of detectors which are required to have a cluster on the track. If a track does not have a cluster from all detectors listed here, it is rejected. If empty, no detector is required. Default is empty. * `require_detectors`: Names of detectors which are required to have a cluster on the track. If a track does not have a cluster from all detectors listed here, it is rejected. If empty, no detector is required. Default is empty.
......
...@@ -9,7 +9,7 @@ using namespace std; ...@@ -9,7 +9,7 @@ using namespace std;
Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detector>> detectors) Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
: Module(std::move(config), std::move(detectors)) { : Module(std::move(config), std::move(detectors)) {
// Default values for cuts // timing cut, relative (x * time_resolution) or absolute:
if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) { if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
throw InvalidCombinationError( throw InvalidCombinationError(
m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive."); m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
...@@ -24,11 +24,27 @@ Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detecto ...@@ -24,11 +24,27 @@ Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detecto
time_cuts_[detector] = detector->getTimeResolution() * time_cut_rel_; time_cuts_[detector] = detector->getTimeResolution() * time_cut_rel_;
} }
} }
spatialCut = m_config.get<double>("spatial_cut", Units::get<double>(200, "um"));
minHitsOnTrack = m_config.get<size_t>("min_hits_on_track", 6); minHitsOnTrack = m_config.get<size_t>("min_hits_on_track", 6);
excludeDUT = m_config.get<bool>("exclude_dut", true); excludeDUT = m_config.get<bool>("exclude_dut", true);
requireDetectors = m_config.getArray<std::string>("require_detectors", {""}); requireDetectors = m_config.getArray<std::string>("require_detectors", {""});
timestampFrom = m_config.get<std::string>("timestamp_from", {}); timestampFrom = m_config.get<std::string>("timestamp_from", {});
// spatial cut, relative (x * spatial_resolution) or absolute:
if(m_config.count({"spatial_cut_rel", "spatial_cut_abs"}) > 1) {
throw InvalidCombinationError(
m_config, {"spatial_cut_rel", "spatial_cut_abs"}, "Absolute and relative spatial cuts are mutually exclusive.");
} else if(m_config.has("spatial_cut_abs")) {
auto spatial_cut_abs_ = m_config.get<XYVector>("spatial_cut_abs");
for(auto& detector : get_detectors()) {
spatial_cuts_[detector] = spatial_cut_abs_;
}
} else {
// default is 3.0 * spatial_resolution
auto spatial_cut_rel_ = m_config.get<double>("spatial_cut_rel", 3.0);
for(auto& detector : get_detectors()) {
spatial_cuts_[detector] = detector->getSpatialResolution() * spatial_cut_rel_;
}
}
} }
void Tracking4D::initialise() { void Tracking4D::initialise() {
...@@ -155,22 +171,30 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) { ...@@ -155,22 +171,30 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
// Get the detector // Get the detector
auto det = get_detector(detectorID); auto det = get_detector(detectorID);
// Check if the DUT should be excluded and obey: if(trees.count(detectorID) == 0) {
if(excludeDUT && det->isDUT()) { LOG(TRACE) << "Skipping detector " << det->name() << " as it has 0 clusters.";
LOG(DEBUG) << "Skipping DUT plane.";
continue; continue;
} }
if(detectorID == seedPlane) if(detectorID == seedPlane) {
LOG(TRACE) << "Skipping seed plane " << det->name();
continue; continue;
if(trees.count(detectorID) == 0) }
// Check if the DUT should be excluded and obey:
if(excludeDUT && det->isDUT()) {
LOG(DEBUG) << "Skipping DUT plane.";
continue; continue;
}
// Get all neighbours within the timing cut // Get all neighbours within the timing cut
LOG(DEBUG) << "Searching for neighbouring cluster on " << detectorID; LOG(DEBUG) << "Searching for neighbouring cluster on device " << detectorID;
LOG(DEBUG) << "- cluster time is " << Units::display(cluster->timestamp(), {"ns", "us", "s"}); LOG(DEBUG) << "- cluster time is " << Units::display(cluster->timestamp(), {"ns", "us", "s"});
Cluster* closestCluster = nullptr; Cluster* closestCluster = nullptr;
double closestClusterDistance = spatialCut;
// Use spatial cut only as initial value (check if cluster is ellipse defined by cuts is done below):
double closestClusterDistance =
sqrt(spatial_cuts_[det].x() * spatial_cuts_[det].x() + spatial_cuts_[det].y() * spatial_cuts_[det].y());
// For default configuration, comparing time cuts calculated from the time resolution of the current detector and // For default configuration, comparing time cuts calculated from the time resolution of the current detector and
// the first plane in Z, // the first plane in Z,
// and taking the maximal value as the cut in time for track-cluster association // and taking the maximal value as the cut in time for track-cluster association
...@@ -202,8 +226,25 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) { ...@@ -202,8 +226,25 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
Cluster* newCluster = neighbours[ne]; Cluster* newCluster = neighbours[ne];
// Calculate the distance to the previous plane's cluster/intercept // Calculate the distance to the previous plane's cluster/intercept
double distance = sqrt((interceptX - newCluster->global().x()) * (interceptX - newCluster->global().x()) + double distanceX = interceptX - newCluster->global().x();
(interceptY - newCluster->global().y()) * (interceptY - newCluster->global().y())); double distanceY = interceptY - newCluster->global().y();
double distance = sqrt(distanceX * distanceX + distanceY * distanceY);
// Check if newCluster lies within ellipse defined by spatial cuts around intercept,
// following this example:
// https://www.geeksforgeeks.org/check-if-a-point-is-inside-outside-or-on-the-ellipse/
//
// ellipse defined by: x^2/a^2 + y^2/b^2 = 1: on ellipse,
// > 1: outside,
// < 1: inside
// Continue if outside of ellipse:
double norm = (distanceX * distanceX) / (spatial_cuts_[det].x() * spatial_cuts_[det].x()) +
(distanceY * distanceY) / (spatial_cuts_[det].y() * spatial_cuts_[det].y());
if(norm > 1) {
continue;
}
// If this is the closest keep it // If this is the closest keep it
if(distance < closestClusterDistance) { if(distance < closestClusterDistance) {
...@@ -295,8 +336,8 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) { ...@@ -295,8 +336,8 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
<< Units::display(track_timestamp, "us") << " to track."; << Units::display(track_timestamp, "us") << " to track.";
track->setTimestamp(track_timestamp); track->setTimestamp(track_timestamp);
} else { } else {
LOG(ERROR) << "Cannot assign timestamp to track. Use average cluster timestamp for track or set detector to " LOG(ERROR) << "Cannot assign timestamp to track. Use average cluster timestamp for track or set detector to set "
"set track timestamp. Please update the configuration file."; "track timestamp. Please update the configuration file.";
return StatusCode::Failure; return StatusCode::Failure;
} }
} }
......
...@@ -43,11 +43,11 @@ namespace corryvreckan { ...@@ -43,11 +43,11 @@ namespace corryvreckan {
// Cuts for tracking // Cuts for tracking
double time_cut_reference_; double time_cut_reference_;
double spatialCut;
size_t minHitsOnTrack; size_t minHitsOnTrack;
bool excludeDUT; bool excludeDUT;
std::vector<std::string> requireDetectors; std::vector<std::string> requireDetectors;
std::map<std::shared_ptr<Detector>, double> time_cuts_; std::map<std::shared_ptr<Detector>, double> time_cuts_;
std::map<std::shared_ptr<Detector>, XYVector> spatial_cuts_;
std::string timestampFrom; std::string timestampFrom;
}; };
} // namespace corryvreckan } // namespace corryvreckan
......
...@@ -8,7 +8,8 @@ This module performs track finding using only positional information (no timing ...@@ -8,7 +8,8 @@ This module performs track finding using only positional information (no timing
### Parameters ### Parameters
* `spatial_cut`: Cut on the maximum distance between the track and cluster for them to be considered associated. Default value is `200um`. * `spatial_cut_rel`: Factor by which the `spatial_resolution` in x and y of each detector plane will be multiplied. These calculated value are defining an ellipse which is then used as the maximum distance in the XY plane allowed between clusters and a track for association to the track. This allows the spatial cuts between different planes to be detector appropriate. By default, a relative spatial cut is applied. Absolute and relative spatial cuts are mutually exclusive. Defaults to `3.0`.
* `spatial_cut_abs`: Specifies a set of absolute value (x and y) which defines an ellipse for the maximum spatial distance in the XY plane between clusters and a track for association to the track. Absolute and relative spatial cuts are mutually exclusive. No default value.
* `min_hits_on_track`: The minimum number of planes with clusters associated to a track for it to be stored. Default value is `6`. * `min_hits_on_track`: The minimum number of planes with clusters associated to a track for it to be stored. Default value is `6`.
* `exclude_dut`: Boolean to set if the DUT should be included in the track fitting. Default value is `true`. * `exclude_dut`: Boolean to set if the DUT should be included in the track fitting. Default value is `true`.
......
...@@ -5,21 +5,38 @@ ...@@ -5,21 +5,38 @@
using namespace corryvreckan; using namespace corryvreckan;
using namespace std; using namespace std;
/*
This algorithm performs the track finding using only spatial information
(no timing). It is based on a linear extrapolation along the z axis, followed
by a nearest neighbour search, and should be well adapted to testbeam
reconstruction with a mostly colinear beam.
*/
TrackingSpatial::TrackingSpatial(Configuration config, std::vector<std::shared_ptr<Detector>> detectors) TrackingSpatial::TrackingSpatial(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
: Module(std::move(config), std::move(detectors)) { : Module(std::move(config), std::move(detectors)) {
spatialCut = m_config.get<double>("spatial_cut", Units::get<double>(200, "um"));
minHitsOnTrack = m_config.get<size_t>("min_hits_on_track", 6); minHitsOnTrack = m_config.get<size_t>("min_hits_on_track", 6);
excludeDUT = m_config.get<bool>("exclude_dut", true); excludeDUT = m_config.get<bool>("exclude_dut", true);
}
/*
This algorithm performs the track finding using only spatial information // spatial cut, relative (x * spatial_resolution) or absolute:
(no timing). It is based on a linear extrapolation along the z axis, followed if(m_config.count({"spatial_cut_rel", "spatial_cut_abs"}) > 1) {
by a nearest neighbour search, and should be well adapted to testbeam throw InvalidCombinationError(
reconstruction with a mostly colinear beam. m_config, {"spatial_cut_rel", "spatial_cut_abs"}, "Absolute and relative spatial cuts are mutually exclusive.");
} else if(m_config.has("spatial_cut_abs")) {
*/ auto spatial_cut_abs_ = m_config.get<XYVector>("spatial_cut_abs");
for(auto& detector : get_detectors()) {
spatial_cuts_[detector] = spatial_cut_abs_;
}
} else {
// default is 3.0 * spatial_resolution
auto spatial_cut_rel_ = m_config.get<double>("spatial_cut_rel", 3.0);
for(auto& detector : get_detectors()) {
spatial_cuts_[detector] = detector->getSpatialResolution() * spatial_cut_rel_;
}
}
}
void TrackingSpatial::initialise() { void TrackingSpatial::initialise() {
...@@ -138,11 +155,11 @@ StatusCode TrackingSpatial::run(std::shared_ptr<Clipboard> clipboard) { ...@@ -138,11 +155,11 @@ StatusCode TrackingSpatial::run(std::shared_ptr<Clipboard> clipboard) {
for(auto& detector : detectors) { for(auto& detector : detectors) {
auto detectorID = detector->name(); auto detectorID = detector->name();
if(trees.count(detectorID) == 0) { if(trees.count(detectorID) == 0) {
LOG(TRACE) << "Skip 0th detector."; LOG(TRACE) << "Skipping detector " << detector->name() << " as it has 0 clusters.";
continue; continue;
} }
if(detectorID == seedPlane) { if(detectorID == seedPlane) {
LOG(TRACE) << "Skip seed plane."; LOG(TRACE) << "Skipping seed plane.";
continue; continue;
} }
...@@ -153,17 +170,28 @@ StatusCode TrackingSpatial::run(std::shared_ptr<Clipboard> clipboard) { ...@@ -153,17 +170,28 @@ StatusCode TrackingSpatial::run(std::shared_ptr<Clipboard> clipboard) {
} }
// Get the closest neighbour // Get the closest neighbour
LOG(DEBUG) << "- looking for nearest cluster on device " << detectorID; LOG(DEBUG) << "Searching for nearest cluster on device " << detectorID;
Cluster* closestCluster = trees[detectorID]->getClosestNeighbour(cluster); Cluster* closestCluster = trees[detectorID]->getClosestNeighbour(cluster);
// Check if it is within the spatial window double distanceX = (cluster->global().x() - closestCluster->global().x());
double distance = sqrt((cluster->global().x() - closestCluster->global().x()) * double distanceY = (cluster->global().y() - closestCluster->global().y());
(cluster->global().x() - closestCluster->global().x()) + double distance = sqrt(distanceX * distanceX + distanceY * distanceY);
(cluster->global().y() - closestCluster->global().y()) *
(cluster->global().y() - closestCluster->global().y())); // Check if closestCluster lies within ellipse defined by spatial cuts,
// following this example:
// https://www.geeksforgeeks.org/check-if-a-point-is-inside-outside-or-on-the-ellipse/
//
// ellipse defined by: x^2/a^2 + y^2/b^2 = 1: on ellipse,
// > 1: outside,
// < 1: inside
// Continue if on or outside of ellipse:
if(distance > spatialCut) double norm = (distanceX * distanceX) / (spatial_cuts_[detector].x() * spatial_cuts_[detector].x()) +
(distanceY * distanceY) / (spatial_cuts_[detector].y() * spatial_cuts_[detector].y());
if(norm > 1) {
continue; continue;
}
// Add the cluster to the track // Add the cluster to the track
track->addCluster(closestCluster); track->addCluster(closestCluster);
......
...@@ -38,7 +38,7 @@ namespace corryvreckan { ...@@ -38,7 +38,7 @@ namespace corryvreckan {
std::map<std::string, TH1F*> residualsY; std::map<std::string, TH1F*> residualsY;
// Member variables // Member variables
double spatialCut; std::map<std::shared_ptr<Detector>, XYVector> spatial_cuts_;
size_t minHitsOnTrack; size_t minHitsOnTrack;
bool excludeDUT; bool excludeDUT;
}; };
......
...@@ -4,75 +4,76 @@ orientation = 0deg,0deg,0deg ...@@ -4,75 +4,76 @@ orientation = 0deg,0deg,0deg
orientation_mode = "xyz" orientation_mode = "xyz"
pixel_pitch = 0um, 0um pixel_pitch = 0um, 0um
position = 0mm,0mm,0mm position = 0mm,0mm,0mm
resolution = 0um,0um spatial_resolution = 0um,0um
time_resolution=1s time_resolution = 1s
type = "TLU" type = "TLU"
role = "auxiliary"
[MIMOSA26_0] [MIMOSA26_0]
mask_file = "mask_MIMOSA26_0.txt" mask_file = "maskfiles/mask_MIMOSA26_0.txt"
number_of_pixels = 1152,576 number_of_pixels = 1152,576
orientation = -0.0317992deg,-0.0306532deg,-0.356494deg