Commit 210a81ae authored by Jens Kroeger's avatar Jens Kroeger
Browse files

resolve merge conflicts after merging master (after time resolution MR was merged)

parents 9602bf7a 02ba5a48
......@@ -286,6 +286,7 @@ All supported rotations are extrinsic active rotations, i.e. the vector itself i
\item The \parameter{number_of_pixels} parameter represents a two-dimensional vector with the number of pixels in the active matrix in the column and row direction, respectively.
\item The \parameter{pixel_pitch} is a two-dimensional vector defining the size of a single pixel.
\item The intrinsic resolution of the detector can be specified using the \parameter{resolution} parameter, a two-dimensional vector holding the position resolution for the column and row directions.
\item The intrinsic time resolution of the detector should be specified using the \parameter{time_resolution} parameter with units of time. This can be used to apply detector specific time cuts in modules.
\item The \parameter{time_offset} can be used to shift the individual detector time frames of reference to e.g.\ account for time of flight effects between different detector planes by adding a fixed offset.
\item Pixels to be masked in the offline analysis can be placed in a separate file specified by the \parameter{mask_file} parameter explained in detail in Section~\ref{sec:masking}.
\item A region of interest in the given detector can be defined using the \parameter{roi} parameter. More details on this functionality can be found in Section~\ref{sec:roi}.
......@@ -300,6 +301,7 @@ orientation = 9deg, 9deg, 0deg
orientation_mode = "xyz"
pixel_pitch = 55um, 55um
position = 0um, 0um, 10mm
time_resolution = 20ns
type = "Timepix3"
[016_CP_PS]
......@@ -309,6 +311,7 @@ orientation = -0.02deg, 0.0deg, -0.015deg
orientation_mode = "xyz"
pixel_pitch = 25um, 25um
position = -0.9mm, 0.21mm, 106.0mm
time_resolution = 1ms
role = "dut"
type = "CLICpix2"
......@@ -318,6 +321,7 @@ orientation = -9deg, 9deg, 0deg
orientation_mode = "xyz"
pixel_pitch = 55um, 55um
position = 0um, 0um, 204mm
time_resolution = 20ns
role = "reference"
type = "Timepix3"
\end{minted}
......
......@@ -84,7 +84,8 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {
m_detectorType = config.get<std::string>("type");
std::transform(m_detectorType.begin(), m_detectorType.end(), m_detectorType.begin(), ::tolower);
m_timingOffset = config.get<double>("time_offset", 0.0);
m_timeOffset = config.get<double>("time_offset", 0.0);
m_timeResolution = config.get<double>("time_resolution");
this->initialise();
......@@ -92,9 +93,10 @@ Detector::Detector(const Configuration& config) : m_role(DetectorRole::NONE) {
<< Units::display(m_pitch, {"mm", "um"});
LOG(TRACE) << " Position: " << Units::display(m_displacement, {"mm", "um"});
LOG(TRACE) << " Orientation: " << Units::display(m_orientation, {"deg"}) << " (" << m_orientation_mode << ")";
if(m_timingOffset > 0.) {
LOG(TRACE) << "Timing offset: " << m_timingOffset;
if(m_timeOffset > 0.) {
LOG(TRACE) << "Time offset: " << m_timeOffset;
}
LOG(TRACE) << " Time resolution: " << Units::display(m_timeResolution, {"ms", "us"});
if(!isAuxiliary()) {
if(config.has("mask_file")) {
......@@ -250,10 +252,12 @@ Configuration Detector::getConfiguration() const {
config.set("orientation_mode", m_orientation_mode);
config.set("orientation", m_orientation, {"deg"});
if(m_timingOffset != 0.) {
config.set("time_offset", m_timingOffset, {"ns", "us", "ms", "s"});
if(m_timeOffset != 0.) {
config.set("time_offset", m_timeOffset, {"ns", "us", "ms", "s"});
}
config.set("time_resolution", m_timeResolution, {"ns", "us", "ms", "s"});
// material budget
if(m_materialBudget > 0.0) {
config.set("material_budget", m_materialBudget);
......
......@@ -137,9 +137,15 @@ namespace corryvreckan {
/**
* @brief Get detector time offset from global clock, can be used to correct for constant shifts or time of flight
* @return Timing offset fo respective detector
* @return Time offset of respective detector
*/
double timingOffset() const { return m_timingOffset; }
double timeOffset() const { return m_timeOffset; }
/**
* @brief Get detector time resolution, used for timing cuts during clustering, track formation, etc.
* @return Time resolutiom of respective detector
*/
double getTimeResolution() const { return m_timeResolution; }
/**
* @brief Update detector position in the world
......@@ -281,7 +287,8 @@ namespace corryvreckan {
XYVector m_pitch{};
XYVector m_spatial_resolution{};
ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<int>> m_nPixels{};
double m_timingOffset;
double m_timeOffset;
double m_timeResolution;
double m_materialBudget;
std::vector<std::vector<int>> m_roi{};
......
......@@ -6,7 +6,14 @@ using namespace std;
Clustering4D::Clustering4D(Configuration config, std::shared_ptr<Detector> detector)
: Module(std::move(config), detector), m_detector(detector) {
timingCut = m_config.get<double>("timing_cut", Units::get<double>(100, "ns"));
if(config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
throw InvalidCombinationError(
m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
} else if(m_config.has("time_cut_abs")) {
timeCut = m_config.get<double>("time_cut_abs");
} else {
timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution();
}
neighbourRadiusRow = m_config.get<int>("neighbour_radius_row", 1);
neighbourRadiusCol = m_config.get<int>("neighbour_radius_col", 1);
chargeWeighting = m_config.get<bool>("charge_weighting", true);
......@@ -31,6 +38,9 @@ void Clustering4D::initialise() {
clusterTimes = new TH1F("clusterTimes", title.c_str(), 3e6, 0, 3e9);
title = m_detector->name() + " Cluster multiplicity;clusters;events";
clusterMultiplicity = new TH1F("clusterMultiplicity", title.c_str(), 50, 0, 50);
// Get resolution in time of detector and calculate time cut to be applied
LOG(DEBUG) << "Time cut to be applied for " << m_detector->name() << " is "
<< Units::display(timeCut, {"ns", "us", "ms"});
}
// Sort function for pixels from low to high times
......@@ -84,9 +94,9 @@ StatusCode Clustering4D::run(std::shared_ptr<Clipboard> clipboard) {
for(size_t iNeighbour = (iP + 1); iNeighbour < totalPixels; iNeighbour++) {
Pixel* neighbour = (*pixels)[iNeighbour];
// Check if they are compatible in time with the cluster pixels
if((neighbour->timestamp() - clusterTime) > timingCut)
if(abs(neighbour->timestamp() - clusterTime) > timeCut)
break;
// if(!closeInTime(neighbour,cluster)) break;
// Check if they have been used
if(used[neighbour])
continue;
......@@ -157,7 +167,7 @@ bool Clustering4D::closeInTime(Pixel* neighbour, Cluster* cluster) {
for(auto& px : pixels) {
double timeDifference = abs(neighbour->timestamp() - px->timestamp());
if(timeDifference < timingCut)
if(timeDifference < timeCut)
CloseInTime = true;
}
return CloseInTime;
......@@ -195,8 +205,6 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) {
column_sum_chargeweighted += (pixel->column() * pixel->charge());
row_sum_chargeweighted += (pixel->row() * pixel->charge());
LOG(DEBUG) << "- cluster col, row: " << column << "," << row;
if(pixel->timestamp() < timestamp) {
timestamp = pixel->timestamp();
}
......@@ -224,6 +232,8 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) {
// Calculate global cluster position
auto positionGlobal = m_detector->localToGlobal(positionLocal);
LOG(DEBUG) << "- cluster col, row, time : " << column << "," << row << ","
<< Units::display(timestamp, {"ns", "us", "ms"});
// Set the cluster parameters
cluster->setColumn(column);
cluster->setRow(row);
......
......@@ -40,7 +40,7 @@ namespace corryvreckan {
TH1F* clusterTimes;
TH1F* clusterMultiplicity;
double timingCut;
double timeCut;
int neighbourRadiusRow;
int neighbourRadiusCol;
bool chargeWeighting;
......
......@@ -14,7 +14,8 @@ Thus, the arithmetic mean is safer.
Split clusters can be recovered using a larger search radius for neighbouring pixels.
### Parameters
* `timing_cut`: The maximum value of the time difference between two pixels for them to be associated in a cluster. Default value is `100ns`.
* `time_cut_rel`: Number of standard deviations the `time_resolution` of the detector plane will be multiplied by. This value is then used as the maximum time difference allowed between pixels for association to a cluster. 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 pixels for association to a cluster. Absolute and relative time cuts are mutually exclusive. No default value.
* `neighbour_radius_col`: Search radius for neighbouring pixels in column direction, defaults to `1` (do not allow split clusters)
* `neighbour_radius_row`: Search radius for neighbouring pixels in row direction, defaults to `1` (do not allow split clusters)
* `charge_weighting`: If true, calculate a charge-weighted mean for the cluster centre. If false, calculate the simple arithmetic mean. Defaults to `true`.
......@@ -34,5 +35,5 @@ For each detector the following plots are produced:
### Usage
```toml
[Clustering4D]
timing_cut = 200ns
time_cut_rel = 3.0
```
......@@ -6,7 +6,17 @@ using namespace std;
DUTAssociation::DUTAssociation(Configuration config, std::shared_ptr<Detector> detector)
: Module(std::move(config), detector), m_detector(detector) {
timingCut = m_config.get<double>("timing_cut", Units::get<double>(200, "ns"));
// timing cut, relative (x * time_resolution) or absolute:
if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
throw InvalidCombinationError(
m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
} else if(m_config.has("time_cut_abs")) {
timeCut = m_config.get<double>("time_cut_abs");
} else {
timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution();
}
// 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.");
......@@ -15,10 +25,9 @@ DUTAssociation::DUTAssociation(Configuration config, std::shared_ptr<Detector> d
} else {
spatialCut = m_config.get<double>("spatial_cut_rel", 3.0) * m_detector->getSpatialResolution();
}
useClusterCentre = m_config.get<bool>("use_cluster_centre", false);
LOG(DEBUG) << "timing_cut = " << Units::display(timingCut, {"ms", "us", "ns"});
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;
}
......@@ -82,6 +91,7 @@ void DUTAssociation::initialise() {
// Nr of associated clusters per track
title = m_detector->name() + ": number of associated clusters per track;associated clusters;events";
hNoAssocCls = new TH1F("no_assoc_cls", title.c_str(), 10, 0, 10);
LOG(DEBUG) << "DUT association time cut = " << Units::display(timeCut, {"ms", "ns"});
}
StatusCode DUTAssociation::run(std::shared_ptr<Clipboard> clipboard) {
......@@ -164,7 +174,7 @@ StatusCode DUTAssociation::run(std::shared_ptr<Clipboard> clipboard) {
}
// Check if the cluster is close in time
if(std::abs(cluster->timestamp() - track->timestamp()) > timingCut) {
if(std::abs(cluster->timestamp() - track->timestamp()) > timeCut) {
LOG(DEBUG) << "Discarding DUT cluster with time difference "
<< Units::display(std::abs(cluster->timestamp() - track->timestamp()), {"ms", "s"});
hCutHisto->Fill(2);
......
......@@ -27,7 +27,7 @@ namespace corryvreckan {
private:
std::shared_ptr<Detector> m_detector;
double timingCut;
double timeCut;
ROOT::Math::XYVector spatialCut;
bool useClusterCentre;
......
......@@ -16,7 +16,8 @@ The other option is to compare the distance between the cluster centre and the t
### Parameters
* `spatial_cut`: Maximum spatial distance in local coordinates in x- and y-direction allowed between cluster and track for association with the DUT. Expects two values for the two coordinates, defaults to twice the pixel pitch.
* `timing_cut`: Maximum time difference allowed between cluster and track for association for the DUT. Default value is `200ns`.
* `time_cut_rel`: Number of standard deviations the `time_resolution` of the detector plane will be multiplied by. This value is then used as the maximum time difference allowed between a DUT cluster and track for association. 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 DUT cluster and track for association. Absolute and relative time cuts are mutually exclusive. No default value.
* `use_cluster_centre`: If set true, the cluster centre will be compared to the track position for the spatial cut. If false, the nearest pixel in the cluster will be used. Defaults to `false`.
### Plots produced
......@@ -35,7 +36,7 @@ The other option is to compare the distance between the cluster centre and the t
```toml
[DUTAssociation]
spatial_cut = 100um, 50um
timing_cut = 200ns
time_cut_rel = 3.0
use_cluster_centre = false
```
......@@ -353,7 +353,7 @@ bool EventLoaderATLASpix::read_caribou_data() { // return false when reaching eo
long long hit_ts = static_cast<long long>(readout_ts_) + ts_diff;
// Convert the timestamp to nanoseconds:
double timestamp = m_clockCycle * static_cast<double>(hit_ts) + m_detector->timingOffset();
double timestamp = m_clockCycle * static_cast<double>(hit_ts) + m_detector->timeOffset();
// If this pixel is masked, do not save it
if(m_detector->masked(col, row)) {
......
......@@ -268,7 +268,7 @@ StatusCode EventLoaderCLICpix2::run(std::shared_ptr<Clipboard> clipboard) {
toa = cp2_pixel->GetTOA();
// Convert ToA form 100MHz clk into ns and sutract from shutterStopTime. Then add configured detector time
// offset
timestamp = shutterStopTime - static_cast<double>(toa) / 0.1 + m_detector->timingOffset();
timestamp = shutterStopTime - static_cast<double>(toa) / 0.1 + m_detector->timeOffset();
if(!discardZeroToT || tot > 0) {
hPixelToA->Fill(toa);
hTimeWalk->Fill(toa, tot);
......
......@@ -259,8 +259,8 @@ Event::Position EventLoaderEUDAQ2::is_within_event(std::shared_ptr<Clipboard> cl
}
// Read time from EUDAQ2 event and convert from picoseconds to nanoseconds:
double event_start = static_cast<double>(evt->GetTimeBegin()) / 1000 + m_detector->timingOffset();
double event_end = static_cast<double>(evt->GetTimeEnd()) / 1000 + m_detector->timingOffset();
double event_start = static_cast<double>(evt->GetTimeBegin()) / 1000 + m_detector->timeOffset();
double event_end = static_cast<double>(evt->GetTimeEnd()) / 1000 + m_detector->timeOffset();
// Store the original position of the event before adjusting its length:
double event_timestamp = event_start;
LOG(DEBUG) << "event_start = " << Units::display(event_start, "us")
......@@ -351,9 +351,9 @@ std::shared_ptr<PixelVector> EventLoaderEUDAQ2::get_pixel_data(std::shared_ptr<e
double ts;
if(plane.GetTimestamp(i) == 0) {
// If the plane timestamp is not defined, we place all pixels in the center of the EUDAQ2 event:
ts = static_cast<double>(evt->GetTimeBegin() + evt->GetTimeEnd()) / 2 / 1000 + m_detector->timingOffset();
ts = static_cast<double>(evt->GetTimeBegin() + evt->GetTimeEnd()) / 2 / 1000 + m_detector->timeOffset();
} else {
ts = static_cast<double>(plane.GetTimestamp(i)) / 1000 + m_detector->timingOffset();
ts = static_cast<double>(plane.GetTimestamp(i)) / 1000 + m_detector->timeOffset();
}
if(col >= m_detector->nPixels().X() || row >= m_detector->nPixels().Y()) {
......
......@@ -61,7 +61,7 @@ The decoder promises to
* return `true` only when there is a fully decoded event available and `false` in all other cases.
* not return any event before a possible T0 signal in the data.
* return the smallest possible granularity of data in time either as even or as sub-events within one event.
* always return valid event time stamps. If the device does not have timestamps, it should return zero for the beginning of the event and have a valid trigger number set.
* always return valid event timestamps. If the device does not have timestamps, it should return zero for the beginning of the event and have a valid trigger number set.
* provide the detector type via the `GetDetectorType()` function in the decoded StandardEvent.
### Configuring EUDAQ2 Event Converters
......
......@@ -556,7 +556,7 @@ bool EventLoaderTimepix3::loadData(std::shared_ptr<Clipboard> clipboard,
}
// Convert final timestamp into ns and add the timing offset (in nano seconds) from the detectors file (if any)
const double timestamp = static_cast<double>(time) / (4096. / 25.) + m_detector->timingOffset();
const double timestamp = static_cast<double>(time) / (4096. / 25.) + m_detector->timeOffset();
auto position = event->getTimestampPosition(timestamp);
// Ignore pixel data if it is before the "eventStart" read from the clipboard storage:
......
......@@ -10,8 +10,9 @@ This module collects `pixel` and `cluster` objects from the clipboard and create
### Parameters
* `make_correlations`: Boolean to change if correlation plots should be outputted. Default value is `false`.
* `do_timing_cut`: Boolean to switch on/off the cut on cluster times for correlations. Defaults to `false`.
* `timing_cut`: maximum time difference between clusters to be taken into account. Only used if `do_timing_cut` is set to `true`, defaults to `100ns`.
* `do_time_cut`: Boolean to switch on/off the cut on cluster times for correlations. Defaults to `false`.
* `time_cut_rel`: Number of standard deviations the `time_resolution` of the detector plane will be multiplied by. This value is then used as the maximum time difference for cluster correlation if `do_time_cut = true`. A relative time cut is applied by default when `do_time_cut = true`. 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 for cluster correlation if `do_time_cut = true`. Absolute and relative time cuts are mutually exclusive. No default value.
* `correlation_vs_time`: Enable plotting of spatial and time correlation as a function of time. Default value is `false` because of the time required to fill the histogram with many bins.
### Plots produced
......
......@@ -7,8 +7,16 @@ TestAlgorithm::TestAlgorithm(Configuration config, std::shared_ptr<Detector> det
: Module(std::move(config), detector), m_detector(detector) {
makeCorrelations = m_config.get<bool>("make_correlations", false);
timingCut = m_config.get<double>("timing_cut", Units::get<double>(100, "ns"));
do_timing_cut_ = m_config.get<bool>("do_timing_cut", false);
do_time_cut_ = m_config.get<bool>("do_time_cut", false);
if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
throw InvalidCombinationError(
m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
} else if(m_config.has("time_cut_abs")) {
timeCut = m_config.get<double>("time_cut_abs");
} else {
timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution();
}
m_corr_vs_time = m_config.get<bool>("correlation_vs_time", false);
}
......@@ -56,8 +64,7 @@ void TestAlgorithm::initialise() {
// time correlation plot range should cover length of events. nanosecond binning.
title = m_detector->name() + "Reference cluster time stamp - cluster time stamp;t_{ref}-t [ns];events";
correlationTime =
new TH1F("correlationTime", title.c_str(), static_cast<int>(2. * timingCut), -1 * timingCut, timingCut);
correlationTime = new TH1F("correlationTime", title.c_str(), static_cast<int>(2. * timeCut), -1 * timeCut, timeCut);
if(m_corr_vs_time) {
title = m_detector->name() + " Correlation X versus time;t [s];x_{ref}-x [mm];events";
......@@ -75,14 +82,14 @@ void TestAlgorithm::initialise() {
3e3,
0,
3e3,
static_cast<int>(2. * timingCut),
-1 * timingCut,
timingCut);
static_cast<int>(2. * timeCut),
-1 * timeCut,
timeCut);
}
title = m_detector->name() + "Reference pixel time stamp - pixel time stamp;t_{ref}-t [ns];events";
correlationTime_px =
new TH1F("correlationTime_px", title.c_str(), static_cast<int>(2. * timingCut), -1 * timingCut, timingCut);
new TH1F("correlationTime_px", title.c_str(), static_cast<int>(2. * timeCut), -1 * timeCut, timeCut);
title = m_detector->name() + "Reference cluster time stamp - cluster time stamp;t_{ref}-t [1/40MHz];events";
correlationTimeInt = new TH1F("correlationTimeInt", title.c_str(), 8000, -40000, 40000);
......@@ -218,7 +225,7 @@ StatusCode TestAlgorithm::run(std::shared_ptr<Clipboard> clipboard) {
long long int timeDifferenceInt = static_cast<long long int>(timeDifference / 25);
// Correlation plots
if(abs(timeDifference) < timingCut || !do_timing_cut_) {
if(abs(timeDifference) < timeCut || !do_time_cut_) {
correlationX->Fill(refCluster->global().x() - cluster->global().x());
correlationX2D->Fill(cluster->global().x(), refCluster->global().x());
correlationX2Dlocal->Fill(cluster->column(), refCluster->column());
......@@ -237,7 +244,7 @@ StatusCode TestAlgorithm::run(std::shared_ptr<Clipboard> clipboard) {
<< ", Time cluster: " << Units::display(cluster->timestamp(), {"ns", "us"});
if(m_corr_vs_time) {
if(abs(timeDifference) < timingCut || !do_timing_cut_) {
if(abs(timeDifference) < timeCut || !do_time_cut_) {
correlationXVsTime->Fill(static_cast<double>(Units::convert(cluster->timestamp(), "s")),
refCluster->global().x() - cluster->global().x());
correlationYVsTime->Fill(static_cast<double>(Units::convert(cluster->timestamp(), "s")),
......
......@@ -53,8 +53,8 @@ namespace corryvreckan {
// Parameters which can be set by user
bool makeCorrelations;
double timingCut;
bool do_timing_cut_;
double timeCut;
bool do_time_cut_;
bool m_corr_vs_time;
};
} // namespace corryvreckan
......
......@@ -9,7 +9,8 @@ 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.
### Parameters
* `timing_cut`: Maximum time difference allowed between clusters for association. Default value is `200ns`.
* `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_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.
* `spatial_cut`: Maximum spatial distance in the XY plane allowed between clusters for association for the telescope planes. Default value is `0.2mm`.
* `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`.
......@@ -37,10 +38,10 @@ For each detector the following plots are produced:
### Usage
```toml
[BasicTracking]
[Tracking4D]
min_hits_on_track = 4
spatial_cut = 300um
timing_cut = 200ns
time_cut_rel = 3.0
exclude_dut = true
require_detectors = "ExampleDetector_0", "ExampleDetector_1"
```
......@@ -9,7 +9,21 @@ using namespace std;
Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
: Module(std::move(config), std::move(detectors)) {
timingCut = m_config.get<double>("timing_cut", Units::get<double>(200, "ns"));
// timing cut, relative (x * time_resolution) or absolute:
if(m_config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
throw InvalidCombinationError(
m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
} else if(m_config.has("time_cut_abs")) {
double time_cut_abs_ = m_config.get<double>("time_cut_abs");
for(auto& detector : get_detectors()) {
time_cuts_[detector] = time_cut_abs_;
}
} else {
double time_cut_rel_ = m_config.get<double>("time_cut_rel", 3.0);
for(auto& detector : get_detectors()) {
time_cuts_[detector] = detector->getTimeResolution() * time_cut_rel_;
}
}
minHitsOnTrack = m_config.get<size_t>("min_hits_on_track", 6);
excludeDUT = m_config.get<bool>("exclude_dut", true);
requireDetectors = m_config.getArray<std::string>("require_detectors", {""});
......@@ -114,6 +128,7 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
LOG(DEBUG) << "Picked up " << tempClusters->size() << " clusters from " << detectorID;
if(firstDetector && !detector->isDUT()) {
referenceClusters = tempClusters;
time_cut_reference_ = time_cuts_[detector];
seedPlane = detector->name();
LOG(DEBUG) << "Seed plane is " << seedPlane;
firstDetector = false;
......@@ -176,10 +191,20 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
LOG(DEBUG) << "Searching for neighbouring cluster on device " << detectorID;
LOG(DEBUG) << "- cluster time is " << Units::display(cluster->timestamp(), {"ns", "us", "s"});
Cluster* closestCluster = nullptr;
// 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());
auto neighbours = trees[detectorID]->getAllClustersInTimeWindow(cluster, timingCut);
// For default configuration, comparing time cuts calculated from the time resolution of the current detector and
// the first plane in Z,
// and taking the maximal value as the cut in time for track-cluster association
// If an absolute cut is to be used, then time_cut_reference_=time_cuts_[det]= time_cut_abs parameter
double timeCut = std::max(time_cut_reference_, time_cuts_[det]);
LOG(TRACE) << "Reference calcuated time cut = " << Units::display(time_cut_reference_, {"ns", "us", "s"})
<< "; detector plane " << detectorID
<< " calculated time cut = " << Units::display(time_cuts_[det], {"ns", "us", "s"});
LOG(DEBUG) << "Using timing cut of " << Units::display(timeCut, {"ns", "us", "s"});
auto neighbours = trees[detectorID]->getAllClustersInTimeWindow(cluster, timeCut);
LOG(DEBUG) << "- found " << neighbours.size() << " neighbours";
......
......@@ -42,10 +42,11 @@ namespace corryvreckan {
std::map<std::string, TH1F*> residualsYwidth3;
// Cuts for tracking
double timingCut;
double time_cut_reference_;
size_t minHitsOnTrack;
bool excludeDUT;
std::vector<std::string> requireDetectors;
std::map<std::shared_ptr<Detector>, double> time_cuts_;
std::map<std::shared_ptr<Detector>, XYVector> spatial_cuts_;
std::string timestampFrom;
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment