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

Rename SpatialClustering -> ClusteringSpatial, make DETECTOR_MODULE

parent a457b199
Pipeline #572780 passed with stages
in 4 minutes and 32 seconds
# Define module and return the generated name as MODULE_NAME
CORRYVRECKAN_GLOBAL_MODULE(MODULE_NAME)
CORRYVRECKAN_DETECTOR_MODULE(MODULE_NAME)
# Add source files to library
CORRYVRECKAN_MODULE_SOURCES(${MODULE_NAME}
SpatialClustering.cpp
# ADD SOURCE FILES HERE...
ClusteringSpatial.cpp
)
# Provide standard install target
......
#include "ClusteringSpatial.h"
#include "objects/Pixel.h"
using namespace corryvreckan;
using namespace std;
ClusteringSpatial::ClusteringSpatial(Configuration config, std::shared_ptr<Detector> detector)
: Module(std::move(config), detector), m_detector(detector) {}
/*
This algorithm clusters the input data, at the moment only if the detector type
is defined as Timepix1. The clustering is based on touching neighbours, and
uses
no timing information whatsoever. It is based on filling a map and checking
neighbours in the neighbouring row and column positions.
*/
void ClusteringSpatial::initialise() {
// Cluster plots
std::string title = m_detector->name() + " Cluster size;cluster size;events";
clusterSize = new TH1F("clusterSize", title.c_str(), 100, 0, 100);
title = m_detector->name() + " Cluster Width - Rows;cluster width [rows];events";
clusterWidthRow = new TH1F("clusterWidthRow", title.c_str(), 25, 0, 25);
title = m_detector->name() + " Cluster Width - Columns;cluster width [columns];events";
clusterWidthColumn = new TH1F("clusterWidthColumn", title.c_str(), 100, 0, 100);
title = m_detector->name() + " Cluster Charge;cluster charge [e];events";
clusterTot = new TH1F("clusterTot", title.c_str(), 10000, 0, 100000);
title = m_detector->name() + " Cluster Position (Global);x [mm];y [mm];events";
clusterPositionGlobal = new TH2F("clusterPositionGlobal",
title.c_str(),
400,
-m_detector->size().X() / 1.5,
m_detector->size().X() / 1.5,
400,
-m_detector->size().Y() / 1.5,
m_detector->size().Y() / 1.5);
}
StatusCode ClusteringSpatial::run(Clipboard* clipboard) {
// Get the pixels
Pixels* pixels = reinterpret_cast<Pixels*>(clipboard->get(m_detector->name(), "pixels"));
if(pixels == nullptr) {
LOG(DEBUG) << "Detector " << m_detector->name() << " does not have any pixels on the clipboard";
return Success;
}
// Make the cluster container and the maps for clustering
Objects* deviceClusters = new Objects();
map<Pixel*, bool> used;
map<int, map<int, Pixel*>> hitmap;
bool addedPixel;
// Get the device dimensions
int nRows = m_detector->nPixelsY();
int nCols = m_detector->nPixelsX();
// Fill the hitmap with pixels
for(auto& pixel : (*pixels)) {
hitmap[pixel->row()][pixel->column()] = pixel;
if(used[pixel]) {
continue;
}
// New pixel => new cluster
Cluster* cluster = new Cluster();
cluster->addPixel(pixel);
cluster->setTimestamp(pixel->timestamp());
used[pixel] = true;
addedPixel = true;
// Somewhere to store found neighbours
Pixels neighbours;
// Now we check the neighbours and keep adding more hits while there are
// connected pixels
while(addedPixel) {
addedPixel = false;
for(int row = (pixel->row() - 1); row <= (pixel->row() + 1); row++) {
// If out of bounds for row
if(row < 0 || row >= nRows)
continue;
for(int col = (pixel->column() - 1); col <= (pixel->column() + 1); col++) {
// If out of bounds for column
if(col < 0 || col >= nCols)
continue;
// If no pixel in this position, or is already in a cluster, do
// nothing
if(!hitmap[row][col])
continue;
if(used[hitmap[row][col]])
continue;
// Otherwise add the pixel to the cluster and store it as a found
// neighbour
cluster->addPixel(hitmap[row][col]);
used[hitmap[row][col]] = true;
neighbours.push_back(hitmap[row][col]);
}
}
// If we have neighbours that have not yet been checked, continue
// looking for more pixels
if(neighbours.size() > 0) {
addedPixel = true;
pixel = neighbours.back();
neighbours.pop_back();
}
}
// Finalise the cluster and save it
calculateClusterCentre(cluster);
// Fill cluster histograms
clusterSize->Fill(static_cast<double>(cluster->size()));
clusterWidthRow->Fill(cluster->rowWidth());
clusterWidthColumn->Fill(cluster->columnWidth());
clusterTot->Fill(cluster->tot());
clusterPositionGlobal->Fill(cluster->globalX(), cluster->globalY());
deviceClusters->push_back(cluster);
}
clipboard->put(m_detector->name(), "clusters", deviceClusters);
LOG(DEBUG) << "Put " << deviceClusters->size() << " clusters on the clipboard for detector " << m_detector->name()
<< ". From " << pixels->size() << " pixels";
// Return value telling analysis to keep running
return Success;
}
/*
Function to calculate the centre of gravity of a cluster.
Sets the local and global cluster positions as well.
*/
void ClusteringSpatial::calculateClusterCentre(Cluster* cluster) {
LOG(DEBUG) << "== Making cluster centre";
// Empty variables to calculate cluster position
double row(0), column(0), tot(0);
// Get the pixels on this cluster
Pixels* pixels = cluster->pixels();
string detectorID = (*pixels)[0]->detectorID();
LOG(DEBUG) << "- cluster has " << (*pixels).size() << " pixels";
// Loop over all pixels
for(auto& pixel : (*pixels)) {
tot += pixel->adc();
row += (pixel->row() * pixel->adc());
column += (pixel->column() * pixel->adc());
LOG(DEBUG) << "- pixel row, col: " << pixel->row() << "," << pixel->column();
}
// Row and column positions are tot-weighted
row /= (tot > 0 ? tot : 1);
column /= (tot > 0 ? tot : 1);
LOG(DEBUG) << "- cluster row, col: " << row << "," << column;
// Create object with local cluster position
PositionVector3D<Cartesian3D<double>> positionLocal(m_detector->pitch().X() * (column - m_detector->nPixelsX() / 2.),
m_detector->pitch().Y() * (row - m_detector->nPixelsY() / 2.),
0);
// Calculate global cluster position
PositionVector3D<Cartesian3D<double>> positionGlobal = m_detector->localToGlobal(positionLocal);
// Set the cluster parameters
cluster->setRow(row);
cluster->setColumn(column);
cluster->setTot(tot);
// Set uncertainty on position from intrinstic detector resolution:
cluster->setError(m_detector->resolution());
cluster->setDetectorID(detectorID);
cluster->setClusterCentre(positionGlobal.X(), positionGlobal.Y(), positionGlobal.Z());
cluster->setClusterCentreLocal(positionLocal.X(), positionLocal.Y(), positionLocal.Z());
}
#ifndef SpatialClustering_H
#define SpatialClustering_H 1
#ifndef ClusteringSpatial_H
#define ClusteringSpatial_H 1
#include <TCanvas.h>
#include <TH1F.h>
......@@ -11,30 +11,28 @@
namespace corryvreckan {
/** @ingroup Modules
*/
class SpatialClustering : public Module {
class ClusteringSpatial : public Module {
public:
// Constructors and destructors
SpatialClustering(Configuration config, std::vector<std::shared_ptr<Detector>> detectors);
~SpatialClustering() {}
ClusteringSpatial(Configuration config, std::shared_ptr<Detector> detector);
~ClusteringSpatial() {}
// Functions
void initialise();
StatusCode run(Clipboard* clipboard);
void finalise();
void calculateClusterCentre(std::shared_ptr<Detector> detector, Cluster*);
private:
// Cluster histograms
std::map<std::string, TH1F*> clusterSize;
std::map<std::string, TH1F*> clusterWidthRow;
std::map<std::string, TH1F*> clusterWidthColumn;
std::map<std::string, TH1F*> clusterTot;
std::map<std::string, TH2F*> clusterPositionGlobal;
std::shared_ptr<Detector> m_detector;
// Member variables
int m_eventNumber;
void calculateClusterCentre(Cluster*);
// Cluster histograms
TH1F* clusterSize;
TH1F* clusterWidthRow;
TH1F* clusterWidthColumn;
TH1F* clusterTot;
TH2F* clusterPositionGlobal;
};
} // namespace corryvreckan
#endif // SpatialClustering_H
#endif // ClusteringSpatial_H
#include "SpatialClustering.h"
#include "objects/Pixel.h"
using namespace corryvreckan;
using namespace std;
SpatialClustering::SpatialClustering(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
: Module(std::move(config), std::move(detectors)) {}
/*
This algorithm clusters the input data, at the moment only if the detector type
is defined as Timepix1. The clustering is based on touching neighbours, and
uses
no timing information whatsoever. It is based on filling a map and checking
neighbours in the neighbouring row and column positions.
*/
void SpatialClustering::initialise() {
for(auto& detector : get_detectors()) {
// Cluster plots
string name = "clusterSize_" + detector->name();
clusterSize[detector->name()] = new TH1F(name.c_str(), name.c_str(), 25, 0, 25);
name = "clusterWidthRow_" + detector->name();
clusterWidthRow[detector->name()] = new TH1F(name.c_str(), name.c_str(), 25, 0, 25);
name = "clusterWidthColumn_" + detector->name();
clusterWidthColumn[detector->name()] = new TH1F(name.c_str(), name.c_str(), 25, 0, 25);
name = "clusterTot_" + detector->name();
clusterTot[detector->name()] = new TH1F(name.c_str(), name.c_str(), 10000, 0, 100000);
name = "clusterPositionGlobal_" + detector->name();
clusterPositionGlobal[detector->name()] = new TH2F(name.c_str(),
name.c_str(),
400,
-detector->size().X() / 1.5,
detector->size().X() / 1.5,
400,
-detector->size().Y() / 1.5,
detector->size().Y() / 1.5);
}
// Initialise member variables
m_eventNumber = 0;
}
StatusCode SpatialClustering::run(Clipboard* clipboard) {
// Loop over all detectors of this algorithm:
for(auto& detector : get_detectors()) {
LOG(TRACE) << "Executing loop for detector " << detector->name();
// Get the pixels
Pixels* pixels = reinterpret_cast<Pixels*>(clipboard->get(detector->name(), "pixels"));
if(pixels == nullptr) {
LOG(DEBUG) << "Detector " << detector->name() << " does not have any pixels on the clipboard";
continue;
}
// Make the cluster container and the maps for clustering
Objects* deviceClusters = new Objects();
map<Pixel*, bool> used;
map<int, map<int, Pixel*>> hitmap;
bool addedPixel;
// Get the device dimensions
int nRows = detector->nPixelsY();
int nCols = detector->nPixelsX();
// Fill the hitmap with pixels
for(auto& pixel : (*pixels)) {
hitmap[pixel->row()][pixel->column()] = pixel;
if(used[pixel]) {
continue;
}
// New pixel => new cluster
Cluster* cluster = new Cluster();
cluster->addPixel(pixel);
cluster->setTimestamp(pixel->timestamp());
used[pixel] = true;
addedPixel = true;
// Somewhere to store found neighbours
Pixels neighbours;
// Now we check the neighbours and keep adding more hits while there are
// connected pixels
while(addedPixel) {
addedPixel = false;
for(int row = (pixel->row() - 1); row <= (pixel->row() + 1); row++) {
// If out of bounds for row
if(row < 0 || row >= nRows)
continue;
for(int col = (pixel->column() - 1); col <= (pixel->column() + 1); col++) {
// If out of bounds for column
if(col < 0 || col >= nCols)
continue;
// If no pixel in this position, or is already in a cluster, do
// nothing
if(!hitmap[row][col])
continue;
if(used[hitmap[row][col]])
continue;
// Otherwise add the pixel to the cluster and store it as a found
// neighbour
cluster->addPixel(hitmap[row][col]);
used[hitmap[row][col]] = true;
neighbours.push_back(hitmap[row][col]);
}
}
// If we have neighbours that have not yet been checked, continue
// looking for more pixels
if(neighbours.size() > 0) {
addedPixel = true;
pixel = neighbours.back();
neighbours.pop_back();
}
}
// Finalise the cluster and save it
calculateClusterCentre(detector, cluster);
// Fill cluster histograms
clusterSize[detector->name()]->Fill(static_cast<double>(cluster->size()));
clusterWidthRow[detector->name()]->Fill(cluster->rowWidth());
clusterWidthColumn[detector->name()]->Fill(cluster->columnWidth());
clusterTot[detector->name()]->Fill(cluster->tot());
clusterPositionGlobal[detector->name()]->Fill(cluster->globalX(), cluster->globalY());
deviceClusters->push_back(cluster);
}
clipboard->put(detector->name(), "clusters", deviceClusters);
LOG(DEBUG) << "Put " << deviceClusters->size() << " clusters on the clipboard for detector " << detector->name()
<< ". From " << pixels->size() << " pixels";
}
// Increment event counter
m_eventNumber++;
// Return value telling analysis to keep running
return Success;
}
void SpatialClustering::finalise() {
LOG(DEBUG) << "Analysed " << m_eventNumber << " events";
}
/*
Function to calculate the centre of gravity of a cluster.
Sets the local and global cluster positions as well.
*/
void SpatialClustering::calculateClusterCentre(std::shared_ptr<Detector> detector, Cluster* cluster) {
LOG(DEBUG) << "== Making cluster centre";
// Empty variables to calculate cluster position
double row(0), column(0), tot(0);
// Get the pixels on this cluster
Pixels* pixels = cluster->pixels();
string detectorID = (*pixels)[0]->detectorID();
LOG(DEBUG) << "- cluster has " << (*pixels).size() << " pixels";
// Loop over all pixels
for(auto& pixel : (*pixels)) {
tot += pixel->adc();
row += (pixel->row() * pixel->adc());
column += (pixel->column() * pixel->adc());
LOG(DEBUG) << "- pixel row, col: " << pixel->row() << "," << pixel->column();
}
// Row and column positions are tot-weighted
row /= (tot > 0 ? tot : 1);
column /= (tot > 0 ? tot : 1);
LOG(DEBUG) << "- cluster row, col: " << row << "," << column;
// Create object with local cluster position
PositionVector3D<Cartesian3D<double>> positionLocal(detector->pitch().X() * (column - detector->nPixelsX() / 2.),
detector->pitch().Y() * (row - detector->nPixelsY() / 2.),
0);
// Calculate global cluster position
PositionVector3D<Cartesian3D<double>> positionGlobal = detector->localToGlobal(positionLocal);
// Set the cluster parameters
cluster->setRow(row);
cluster->setColumn(column);
cluster->setTot(tot);
// Set uncertainty on position from intrinstic detector resolution:
cluster->setError(detector->resolution());
cluster->setDetectorID(detectorID);
cluster->setClusterCentre(positionGlobal.X(), positionGlobal.Y(), positionGlobal.Z());
cluster->setClusterCentreLocal(positionLocal.X(), positionLocal.Y(), positionLocal.Z());
}
Supports Markdown
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