ClusteringSpatial.cpp 7.82 KB
Newer Older
1
#include "ClusteringSpatial.h"
Simon Spannagel's avatar
Simon Spannagel committed
2
#include "objects/Pixel.hpp"
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

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);
29
30
    title = m_detector->name() + " Cluster Charge;cluster charge [ke];events";
    clusterTot = new TH1F("clusterTot", title.c_str(), 300, 0, 300);
31
32
33
34
35
36
37
38
39
    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);
40
41
42
43
44
45
46
47
48
49
    title = m_detector->name() + " Cluster Position (Local);x [px];y [px];events";

    clusterPositionLocal = new TH2F("clusterPositionLocal",
                                    title.c_str(),
                                    m_detector->nPixels().X(),
                                    -m_detector->nPixels().X() / 2.,
                                    m_detector->nPixels().X() / 2.,
                                    m_detector->nPixels().Y(),
                                    -m_detector->nPixels().Y() / 2.,
                                    m_detector->nPixels().Y() / 2.);
50
51
}

52
StatusCode ClusteringSpatial::run(std::shared_ptr<Clipboard> clipboard) {
53
54
55
56
57

    // 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";
58
        return StatusCode::Success;
59
60
61
62
63
64
65
66
67
    }

    // 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
68
69
    int nRows = m_detector->nPixels().Y();
    int nCols = m_detector->nPixels().X();
70

71
72
    // Pre-fill the hitmap with pixels
    for(auto pixel : (*pixels)) {
73
        hitmap[pixel->row()][pixel->column()] = pixel;
74
    }
75

76
    for(auto pixel : (*pixels)) {
77
78
79
80
81
82
83
84
85
86
87
88
89
        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;

90
        // Now we check the neighbours and keep adding more hits while there are connected pixels
91
92
93
        while(addedPixel) {

            addedPixel = false;
94
            for(int row = pixel->row() - 1; row <= pixel->row() + 1; row++) {
95
                // If out of bounds for row
96
                if(row < 0 || row >= nRows) {
97
                    continue;
98
99
100
                }

                for(int col = pixel->column() - 1; col <= pixel->column() + 1; col++) {
101
                    // If out of bounds for column
102
                    if(col < 0 || col >= nCols) {
103
                        continue;
104
                    }
105

106
                    // If no pixel in this position, or is already in a cluster, do nothing
107
                    if(!hitmap[row][col]) {
108
                        continue;
109
110
                    }
                    if(used[hitmap[row][col]]) {
111
                        continue;
112
                    }
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

                    // 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());
138
        clusterTot->Fill(cluster->tot() * 1e-3);
139
        clusterPositionGlobal->Fill(cluster->global().x(), cluster->global().y());
140
141
        clusterPositionLocal->Fill(cluster->local().x(), cluster->local().y());
        LOG(DEBUG) << "cluster local: " << cluster->local();
142
143
144
145
146
147
148
149
        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
150
    return StatusCode::Success;
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
}

/*
 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
183
184
    PositionVector3D<Cartesian3D<double>> positionLocal(m_detector->pitch().X() * (column - m_detector->nPixels().X() / 2.),
                                                        m_detector->pitch().Y() * (row - m_detector->nPixels().Y() / 2.),
185
186
187
188
189
190
191
192
193
194
195
196
197
                                                        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);
198
199
    cluster->setClusterCentre(positionGlobal);
    cluster->setClusterCentreLocal(positionLocal);
200
}