#include "Timepix3Clustering.h" using namespace corryvreckan; using namespace std; Timepix3Clustering::Timepix3Clustering(Configuration config, std::vector detectors) : Algorithm(std::move(config), std::move(detectors)) { timingCut = m_config.get("timingCut", 0.0000001); // 100 ns } void Timepix3Clustering::initialise() { timingCutInt = (timingCut * 4096. * 40000000.); 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(), 200, 0, 1000); name = "clusterPositionGlobal_" + detector->name(); clusterPositionGlobal[detector->name()] = new TH2F(name.c_str(), name.c_str(), 400, -10., 10., 400, -10., 10.); } } // Sort function for pixels from low to high times bool sortByTime(Pixel* pixel1, Pixel* pixel2) { return (pixel1->m_timestamp < pixel2->m_timestamp); } StatusCode Timepix3Clustering::run(Clipboard* clipboard) { // Loop over all Timepix3 and for each device perform the clustering for(auto& detector : get_detectors()) { // Check if they are a Timepix3 if(detector->type() != "Timepix3") continue; // Get the pixels Pixels* pixels = (Pixels*)clipboard->get(detector->name(), "pixels"); if(pixels == NULL) { LOG(DEBUG) << "Detector " << detector->name() << " does not have any pixels on the clipboard"; continue; } LOG(DEBUG) << "Picked up " << pixels->size() << " pixels for device " << detector->name(); // if(pixels->size() > 500.){ // LOG(DEBUG) <<"Skipping large event with "<size()<<" pixels // for device "<name(); // continue; // } // Sort the pixels from low to high timestamp std::sort(pixels->begin(), pixels->end(), sortByTime); int totalPixels = pixels->size(); // Make the cluster storage Clusters* deviceClusters = new Clusters(); // Keep track of which pixels are used map used; // Start to cluster for(int iP = 0; iP < pixels->size(); iP++) { Pixel* pixel = (*pixels)[iP]; // Check if pixel is used if(used[pixel]) continue; // Make the new cluster object Cluster* cluster = new Cluster(); LOG(DEBUG) << "==== New cluster"; // Keep adding hits to the cluster until no more are found cluster->addPixel(pixel); long long int clusterTime = pixel->m_timestamp; used[pixel] = true; LOG(DEBUG) << "Adding pixel: " << pixel->m_row << "," << pixel->m_column; int nPixels = 0; while(cluster->size() != nPixels) { nPixels = cluster->size(); // Loop over all pixels for(int iNeighbour = (iP + 1); iNeighbour < totalPixels; iNeighbour++) { Pixel* neighbour = (*pixels)[iNeighbour]; // Check if they are compatible in time with the cluster pixels if((neighbour->m_timestamp - clusterTime) > timingCutInt) break; // if(!closeInTime(neighbour,cluster)) break; // Check if they have been used if(used[neighbour]) continue; // Check if they are touching cluster pixels if(!touching(neighbour, cluster)) continue; // Add to cluster cluster->addPixel(neighbour); clusterTime = neighbour->m_timestamp; used[neighbour] = true; LOG(DEBUG) << "Adding pixel: " << neighbour->m_row << "," << neighbour->m_column; } } // Finalise the cluster and save it calculateClusterCentre(cluster); deviceClusters->push_back(cluster); } // Put the clusters on the clipboard if(deviceClusters->size() > 0) clipboard->put(detector->name(), "clusters", (TestBeamObjects*)deviceClusters); LOG(DEBUG) << "Made " << deviceClusters->size() << " clusters for device " << detector->name(); } for(auto& detector : get_detectors()) { // Get the clusters Clusters* clusters = (Clusters*)clipboard->get(detector->name(), "clusters"); if(clusters == NULL) { LOG(DEBUG) << "Detector " << detector->name() << " does not have any clusters on the clipboard"; continue; } // Loop over all clusters and fill histograms for(auto& cluster : (*clusters)) { // Fill cluster histograms clusterSize[detector->name()]->Fill(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()); } } return Success; } // Check if a pixel touches any of the pixels in a cluster bool Timepix3Clustering::touching(Pixel* neighbour, Cluster* cluster) { bool Touching = false; Pixels* pixels = cluster->pixels(); for(int iPix = 0; iPix < pixels->size(); iPix++) { if(abs((*pixels)[iPix]->m_row - neighbour->m_row) <= 1 && abs((*pixels)[iPix]->m_column - neighbour->m_column) <= 1) { Touching = true; break; } } return Touching; } // Check if a pixel is close in time to the pixels of a cluster bool Timepix3Clustering::closeInTime(Pixel* neighbour, Cluster* cluster) { bool CloseInTime = false; Pixels* pixels = cluster->pixels(); for(int iPix = 0; iPix < pixels->size(); iPix++) { long long int timeDifference = abs(neighbour->m_timestamp - (*pixels)[iPix]->m_timestamp); if(timeDifference < timingCutInt) CloseInTime = true; } return CloseInTime; } void Timepix3Clustering::calculateClusterCentre(Cluster* cluster) { // Empty variables to calculate cluster position double row(0), column(0), tot(0); long long int timestamp; // Get the pixels on this cluster Pixels* pixels = cluster->pixels(); string detectorID = (*pixels)[0]->m_detectorID; timestamp = (*pixels)[0]->m_timestamp; // Loop over all pixels for(int pix = 0; pix < pixels->size(); pix++) { double pixelToT = (*pixels)[pix]->m_adc; if(pixelToT == 0) { LOG(DEBUG) << "Pixel with ToT 0!"; pixelToT = 1; } tot += pixelToT; row += ((*pixels)[pix]->m_row * pixelToT); column += ((*pixels)[pix]->m_column * pixelToT); if((*pixels)[pix]->m_timestamp < timestamp) timestamp = (*pixels)[pix]->m_timestamp; } // Row and column positions are tot-weighted row /= tot; column /= tot; auto detector = get_detector(detectorID); // Create object with local cluster position PositionVector3D> positionLocal( detector->pitchX() * (column - detector->nPixelsX() / 2), detector->pitchY() * (row - detector->nPixelsY() / 2), 0); // Calculate global cluster position PositionVector3D> positionGlobal = *(detector->m_localToGlobal) * positionLocal; // Set the cluster parameters cluster->setRow(row); cluster->setColumn(column); cluster->setTot(tot); cluster->setError(0.004); cluster->setTimestamp(timestamp); cluster->setDetectorID(detectorID); cluster->setClusterCentre(positionGlobal.X(), positionGlobal.Y(), positionGlobal.Z()); cluster->setClusterCentreLocal(positionLocal.X(), positionLocal.Y(), positionLocal.Z()); } void Timepix3Clustering::finalise() {}