Commit 5cedc587 authored by Morag Jean Williams's avatar Morag Jean Williams
Browse files

Merge branch 'improve_clustering4d' into 'master'

Improve Clustering - allow arithmetic or charge-weighted mean for cluster centre

See merge request !114
parents a7a0b683 bb5a4ddc
Pipeline #1083373 passed with stages
in 17 minutes and 56 seconds
......@@ -7,8 +7,9 @@ Clustering4D::Clustering4D(Configuration config, std::shared_ptr<Detector> detec
: Module(std::move(config), detector), m_detector(detector) {
timingCut = m_config.get<double>("timing_cut", Units::get<double>(100, "ns"));
neighbour_radius_row = m_config.get<int>("neighbour_radius_row", 1);
neighbour_radius_col = m_config.get<int>("neighbour_radius_col", 1);
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);
}
void Clustering4D::initialise() {
......@@ -138,7 +139,7 @@ bool Clustering4D::touching(Pixel* neighbour, Cluster* cluster) {
int row_distance = abs(pixel->row() - neighbour->row());
int col_distance = abs(pixel->column() - neighbour->column());
if(row_distance <= neighbour_radius_row && col_distance <= neighbour_radius_col) {
if(row_distance <= neighbourRadiusRow && col_distance <= neighbourRadiusCol) {
if(row_distance > 1 || col_distance > 1) {
cluster->setSplit(true);
}
......@@ -169,6 +170,9 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) {
LOG(DEBUG) << "== Making cluster centre";
// Empty variables to calculate cluster position
double column(0), row(0), charge(0);
double column_sum(0), column_sum_chargeweighted(0);
double row_sum(0), row_sum_chargeweighted(0);
bool found_charge_zero = false;
// Get the pixels on this cluster
Pixels* pixels = cluster->pixels();
......@@ -178,21 +182,38 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) {
// Loop over all pixels
for(auto& pixel : (*pixels)) {
// If charge == 0 (use epsilon to avoid errors in floating-point arithmetics):
if(pixel->charge() < std::numeric_limits<double>::epsilon()) {
// apply arithmetic mean if a pixel has zero charge
found_charge_zero = true;
}
charge += pixel->charge();
column += (pixel->column() * pixel->charge());
row += (pixel->row() * pixel->charge());
// We need both column_sum and column_sum_chargeweighted
// as we don't know a priori if there will be a pixel with
// charge==0 such that we have to fall back to the arithmetic mean.
column_sum += (pixel->column() * pixel->charge());
row_sum += (pixel->row() * pixel->charge());
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();
}
}
// Column and row positions are charge-weighted
// If charge == 0 (use epsilon to avoid errors in floating-point arithmetics)
// calculate simple arithmetic mean
column /= (charge > std::numeric_limits<double>::epsilon() ? charge : 1);
row /= (charge > std::numeric_limits<double>::epsilon() ? charge : 1);
LOG(DEBUG) << "- cluster col, row: " << column << "," << row;
if(chargeWeighting && !found_charge_zero) {
// Charge-weighted centre-of-gravity for cluster centre:
// (here it's safe to divide by the charge as it cannot be zero due to !found_charge_zero)
column = column_sum_chargeweighted / charge;
row = row_sum_chargeweighted / charge;
} else {
// Arithmetic cluster centre:
column = column_sum / static_cast<double>(cluster->size());
row = row_sum / static_cast<double>(cluster->size());
}
if(detectorID != m_detector->name()) {
// Should never happen...
......
......@@ -41,8 +41,9 @@ namespace corryvreckan {
TH1F* clusterMultiplicity;
double timingCut;
int neighbour_radius_row;
int neighbour_radius_col;
int neighbourRadiusRow;
int neighbourRadiusCol;
bool chargeWeighting;
};
} // namespace corryvreckan
#endif // CLUSTERING4D_H
......@@ -5,7 +5,11 @@
**Status**: Functional
### Description
This module performs clustering for detectors with valid individual hit timestamps. The clustering method is a charge-weighted centre of gravity calculation, using a positional cut and a timing cut on proximity. If the pixel information is binary (i.e. no valid charge-equivalent information is available), the arithmetic mean is calculated for the position.
This module performs clustering for detectors with valid individual hit timestamps.
The clustering method is either an arithmetic mean or a a charge-weighted centre-of-gravity calculation, using a positional cut and a timing cut on proximity.
If the pixel information is binary (i.e. no valid charge-equivalent information is available), the arithmetic mean is calculated for the position.
Also, if one pixel of a cluster has charge zero, the arithmetic mean is calculated even if charge-weighting is selected because it is assumed that the zero-reading is false and does not to represent a low charge but an unknown value.
Thus, the arithmetic mean is safer.
Split clusters can be recovered using a larger search radius for neighbouring pixels.
......@@ -13,6 +17,7 @@ Split clusters can be recovered using a larger search radius for neighbouring pi
* `timing_cut`: The maximum value of the time difference between two pixels for them to be associated in a cluster. Default value is `100ns`.
* `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`.
### Plots produced
For each detector the following plots are produced:
......
......@@ -8,6 +8,7 @@ ClusteringSpatial::ClusteringSpatial(Configuration config, std::shared_ptr<Detec
: Module(std::move(config), detector), m_detector(detector) {
useTriggerTimestamp = m_config.get<bool>("use_trigger_timestamp", false);
chargeWeighting = m_config.get<bool>("charge_weighting", true);
}
void ClusteringSpatial::initialise() {
......@@ -176,6 +177,9 @@ void ClusteringSpatial::calculateClusterCentre(Cluster* cluster) {
LOG(DEBUG) << "== Making cluster centre";
// Empty variables to calculate cluster position
double column(0), row(0), charge(0);
double column_sum(0), column_sum_chargeweighted(0);
double row_sum(0), row_sum_chargeweighted(0);
bool found_charge_zero = false;
// Get the pixels on this cluster
Pixels* pixels = cluster->pixels();
......@@ -184,18 +188,34 @@ void ClusteringSpatial::calculateClusterCentre(Cluster* cluster) {
// Loop over all pixels
for(auto& pixel : (*pixels)) {
// If charge == 0 (use epsilon to avoid errors in floating-point arithmetics):
if(pixel->charge() < std::numeric_limits<double>::epsilon()) {
// apply arithmetic mean if a pixel has zero charge
found_charge_zero = true;
}
charge += pixel->charge();
column += (pixel->column() * pixel->charge());
row += (pixel->row() * pixel->charge());
// We need both column_sum and column_sum_chargeweighted
// as we don't know a priori if there will be a pixel with
// charge==0 such that we have to fall back to the arithmetic mean.
column_sum += (pixel->column() * pixel->charge());
row_sum += (pixel->row() * pixel->charge());
column_sum_chargeweighted += (pixel->column() * pixel->charge());
row_sum_chargeweighted += (pixel->row() * pixel->charge());
LOG(DEBUG) << "- pixel col, row: " << pixel->column() << "," << pixel->row();
}
// Column and row positions are charge-weighted
// If charge == 0 (use epsilon to avoid errors in floating-point arithmetics)
// calculate simple arithmetic mean
column /= (charge > std::numeric_limits<double>::epsilon() ? charge : 1);
row /= (charge > std::numeric_limits<double>::epsilon() ? charge : 1);
if(chargeWeighting && !found_charge_zero) {
// Charge-weighted centre-of-gravity for cluster centre:
// (here it's safe to divide by the charge as it cannot be zero due to !found_charge_zero)
column = column_sum_chargeweighted / charge;
row = row_sum_chargeweighted / charge;
} else {
// Arithmetic cluster centre:
column = column_sum / static_cast<double>(cluster->size());
row = row_sum / static_cast<double>(cluster->size());
}
LOG(DEBUG) << "- cluster col, row: " << column << "," << row << " at time "
<< Units::display(cluster->timestamp(), "us");
......
......@@ -38,6 +38,7 @@ namespace corryvreckan {
TH1F* clusterTimes;
bool useTriggerTimestamp;
bool chargeWeighting;
};
} // namespace corryvreckan
#endif // ClusteringSpatial_H
......@@ -5,10 +5,15 @@
**Status**: Functioning
### Description
This module clusters the input data of a detector without individual hit timestamps. The clustering method only uses positional information: charge-weighted centre of gravity calculation using touching neighbours method, and no timing information. If the pixel information is binary (i.e. no valid charge-equivalent information is available), the arithmetic mean is calculated for the position. These clusters are stored on the clipboard for each device.
This module clusters the input data of a detector without individual hit timestamps.
The clustering method only uses positional information: either charge-weighted centre-of-gravity or arithmetic mean calculation using touching neighbours method, and no timing information.
If the pixel information is binary (i.e. no valid charge-equivalent information is available), the arithmetic mean is calculated for the position.
Also, if one pixel of a cluster has charge zero, the arithmetic mean is calculated even if charge-weighting is selected because it is assumed that the zero-reading is false and does not to represent a low charge but an unknown value.
These clusters are stored on the clipboard for each device.
### Parameters
* `use_trigger_timestamp`: If true, set trigger timestamp of Corryvreckan event as cluster timestamp. If false, set pixel timestamp. Default value is `false`.
* `charge_weighting`: If true, calculate a charge-weighted mean for the cluster centre. If false, calculate the simple arithmetic mean. Defaults to `true`.
### Plots produced
For each detector the following plots are produced:
......
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