ClusteringSpatial.cpp 9.6 KB
Newer Older
1
#include "ClusteringSpatial.h"
Simon Spannagel's avatar
Simon Spannagel committed
2
#include "objects/Pixel.hpp"
3
4
5
6
7

using namespace corryvreckan;
using namespace std;

ClusteringSpatial::ClusteringSpatial(Configuration config, std::shared_ptr<Detector> detector)
8
    : Module(std::move(config), detector), m_detector(detector) {
9

10
    useTriggerTimestamp = m_config.get<bool>("use_trigger_timestamp", false);
11
    chargeWeighting = m_config.get<bool>("charge_weighting", true);
12
}
13
14
15
16
17
18

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);
19
    title = m_detector->name() + " Cluster seed charge;cluster seed charge [e];events";
20
    clusterSeedCharge = new TH1F("clusterSeedCharge", title.c_str(), 256, 0, 256);
21
22
23
24
    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);
Katharina Dort's avatar
Katharina Dort committed
25
26
    title = m_detector->name() + " Cluster Charge;cluster charge [e];events";
    clusterCharge = new TH1F("clusterCharge", title.c_str(), 5000, 0, 50000);
27
28
29
30
31
32
33
34
35
    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);
36
37
38
39
40
41
42
43
44
    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.);
45
46
47

    title = ";cluster timestamp [ns]; # events";
    clusterTimes = new TH1F("clusterTimes", title.c_str(), 3e6, 0, 3e9);
48
49
}

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

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

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

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

74
    for(auto pixel : (*pixels)) {
75
76
77
78
79
80
81
        if(used[pixel]) {
            continue;
        }

        // New pixel => new cluster
        Cluster* cluster = new Cluster();
        cluster->addPixel(pixel);
82
83
84

        if(useTriggerTimestamp) {
            if(!clipboard->get_event()->triggerList().empty()) {
85
86
87
                double trigger_ts = clipboard->get_event()->triggerList().begin()->second;
                LOG(DEBUG) << "Using trigger timestamp " << Units::display(trigger_ts, "us") << " as cluster timestamp.";
                cluster->setTimestamp(trigger_ts);
88
            } else {
89
90
                LOG(WARNING) << "No trigger available. Use pixel timestamp " << Units::display(pixel->timestamp(), "us")
                             << " as cluster timestamp.";
91
                cluster->setTimestamp(pixel->timestamp());
92
93
94
            }
        } else {
            // assign pixel timestamp
95
96
            LOG(DEBUG) << "Pixel has timestamp " << Units::display(pixel->timestamp(), "us")
                       << ", set as cluster timestamp. ";
97
98
99
            cluster->setTimestamp(pixel->timestamp());
        }

100
101
102
103
104
        used[pixel] = true;
        addedPixel = true;
        // Somewhere to store found neighbours
        Pixels neighbours;

105
        // Now we check the neighbours and keep adding more hits while there are connected pixels
106
107
108
        while(addedPixel) {

            addedPixel = false;
109
            for(int row = pixel->row() - 1; row <= pixel->row() + 1; row++) {
110
                // If out of bounds for row
111
                if(row < 0 || row >= nRows) {
112
                    continue;
113
114
115
                }

                for(int col = pixel->column() - 1; col <= pixel->column() + 1; col++) {
116
                    // If out of bounds for column
117
                    if(col < 0 || col >= nCols) {
118
                        continue;
119
                    }
120

121
                    // If no pixel in this position, or is already in a cluster, do nothing
122
                    if(!hitmap[col][row]) {
123
                        continue;
124
                    }
125
                    if(used[hitmap[col][row]]) {
126
                        continue;
127
                    }
128
129
130

                    // Otherwise add the pixel to the cluster and store it as a found
                    // neighbour
131
132
133
                    cluster->addPixel(hitmap[col][row]);
                    used[hitmap[col][row]] = true;
                    neighbours.push_back(hitmap[col][row]);
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
                }
            }

            // 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()));
151

152
153
        clusterWidthRow->Fill(cluster->rowWidth());
        clusterWidthColumn->Fill(cluster->columnWidth());
154
155
        clusterCharge->Fill(cluster->charge());
        clusterSeedCharge->Fill(cluster->getSeedPixel()->charge());
156
        clusterPositionGlobal->Fill(cluster->global().x(), cluster->global().y());
157
        clusterPositionLocal->Fill(cluster->local().x(), cluster->local().y());
158
        clusterTimes->Fill(static_cast<double>(Units::convert(cluster->timestamp(), "ns")));
159
        LOG(DEBUG) << "cluster local: " << cluster->local();
160
161
162
163
164
165
166
167
        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
168
    return StatusCode::Success;
169
170
171
172
173
174
175
176
177
178
}

/*
 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
179
    double column(0), row(0), charge(0);
180
    bool found_charge_zero = false;
181
182
183
184
185
186
187
188

    // 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)) {
189
190
191
192
        if(pixel->charge() < std::numeric_limits<double>::epsilon()) {
            // apply arithmetic mean if a pixel has zero charge
            found_charge_zero = true;
        }
193
        charge += pixel->charge();
194
195
196
197
198
199
200
201
202
203

        if(chargeWeighting) {
            // charge-weighted cluster centre:
            column += (pixel->column() * pixel->charge());
            row += (pixel->row() * pixel->charge());
        } else {
            // arithmetic cluster centre:
            column += pixel->column();
            row += pixel->row();
        }
204
205

        LOG(DEBUG) << "- pixel col, row: " << pixel->column() << "," << pixel->row();
206
207
    }

208
209
210
211
212
213
214
215
216
217
    if(chargeWeighting && !found_charge_zero) {
        // Charge-weighted cluster centre:
        // If charge == 0 (use epsilon to avoid errors in floating-point arithmetics).
        column /= (charge > std::numeric_limits<double>::epsilon() ? charge : 1);
        row /= (charge > std::numeric_limits<double>::epsilon() ? charge : 1);
    } else {
        // Arightmetic cluster centre:
        column /= static_cast<double>(cluster->size());
        row /= static_cast<double>(cluster->size());
    }
218

219
220
    LOG(DEBUG) << "- cluster col, row: " << column << "," << row << " at time "
               << Units::display(cluster->timestamp(), "us");
221
222

    // Create object with local cluster position
223
224
    auto positionLocal = m_detector->getLocalPosition(column, row);

225
    // Calculate global cluster position
226
    auto positionGlobal = m_detector->localToGlobal(positionLocal);
227
228
229
230

    // Set the cluster parameters
    cluster->setRow(row);
    cluster->setColumn(column);
231
    cluster->setCharge(charge);
232
233
234
235
236

    // Set uncertainty on position from intrinstic detector resolution:
    cluster->setError(m_detector->resolution());

    cluster->setDetectorID(detectorID);
237
238
    cluster->setClusterCentre(positionGlobal);
    cluster->setClusterCentreLocal(positionLocal);
239
}