Clustering4D.cpp 7.83 KB
Newer Older
1
#include "Clustering4D.h"
2

3
using namespace corryvreckan;
4
using namespace std;
5

6
Clustering4D::Clustering4D(Configuration config, std::shared_ptr<Detector> detector)
7
    : Module(std::move(config), detector), m_detector(detector) {
8

9
    timingCut = m_config.get<double>("timing_cut", static_cast<double>(Units::convert(100, "ns"))); // 100 ns
10
11
    neighbour_radius_row = m_config.get<int>("neighbour_radius_row", 1);
    neighbour_radius_col = m_config.get<int>("neighbour_radius_col", 1);
12
13
}

14
void Clustering4D::initialise() {
15

16
    // Cluster plots
17
18
19
20
21
22
23
24
25
26
    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, -10., 10., 400, -10., 10.);
27
28
}

29
// Sort function for pixels from low to high times
30
bool Clustering4D::sortByTime(Pixel* pixel1, Pixel* pixel2) {
31
    return (pixel1->timestamp() < pixel2->timestamp());
32
33
}

34
StatusCode Clustering4D::run(std::shared_ptr<Clipboard> clipboard) {
35

36
    // Get the pixels
37
    Pixels* pixels = reinterpret_cast<Pixels*>(clipboard->get(m_detector->name(), "pixels"));
38
    if(pixels == nullptr) {
39
        LOG(DEBUG) << "Detector " << m_detector->name() << " does not have any pixels on the clipboard";
40
        return StatusCode::Success;
41
    }
42
    LOG(DEBUG) << "Picked up " << pixels->size() << " pixels for device " << m_detector->name();
43
44
45

    // Sort the pixels from low to high timestamp
    std::sort(pixels->begin(), pixels->end(), sortByTime);
46
    size_t totalPixels = pixels->size();
47
48
49
50
51
52

    // Make the cluster storage
    Clusters* deviceClusters = new Clusters();

    // Keep track of which pixels are used
    map<Pixel*, bool> used;
53

54
    // Start to cluster
55
    for(size_t iP = 0; iP < pixels->size(); iP++) {
56
57
58
        Pixel* pixel = (*pixels)[iP];

        // Check if pixel is used
59
        if(used[pixel]) {
Simon Spannagel's avatar
Simon Spannagel committed
60
            continue;
61
        }
62

63
        if(pixel->adc() == 0.) {
Simon Spannagel's avatar
Simon Spannagel committed
64
            continue;
65
        }
66
67
68
69
70
71
72
73
74
75

        // 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);
        double clusterTime = pixel->timestamp();
        used[pixel] = true;
        LOG(DEBUG) << "Adding pixel: " << pixel->row() << "," << pixel->column();
76
        size_t nPixels = 0;
77
78
79
80
        while(cluster->size() != nPixels) {

            nPixels = cluster->size();
            // Loop over all pixels
81
            for(size_t iNeighbour = (iP + 1); iNeighbour < totalPixels; iNeighbour++) {
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
                Pixel* neighbour = (*pixels)[iNeighbour];
                // Check if they are compatible in time with the cluster pixels
                if((neighbour->timestamp() - clusterTime) > timingCut)
                    break;
                //          if(!closeInTime(neighbour,cluster)) break;
                // Check if they have been used
                if(used[neighbour])
                    continue;

                if(neighbour->adc() == 0.)
                    continue;

                // Check if they are touching cluster pixels
                if(!touching(neighbour, cluster))
                    continue;

                // Add to cluster
                cluster->addPixel(neighbour);
                clusterTime = neighbour->timestamp();
                used[neighbour] = true;
                LOG(DEBUG) << "Adding pixel: " << neighbour->row() << "," << neighbour->column() << " time "
                           << Units::display(neighbour->timestamp(), {"ns", "us", "s"});
Simon Spannagel's avatar
Simon Spannagel committed
104
            }
105
        }
Simon Spannagel's avatar
Simon Spannagel committed
106

107
108
        // Finalise the cluster and save it
        calculateClusterCentre(cluster);
109

110
        // Fill cluster histograms
111
        clusterSize->Fill(static_cast<double>(cluster->size()));
112
113
114
        clusterWidthRow->Fill(cluster->rowWidth());
        clusterWidthColumn->Fill(cluster->columnWidth());
        clusterTot->Fill(cluster->tot());
115
        clusterPositionGlobal->Fill(cluster->global().x(), cluster->global().y());
116

117
118
        deviceClusters->push_back(cluster);
    }
119

120
121
    // Put the clusters on the clipboard
    if(deviceClusters->size() > 0) {
122
        clipboard->put(m_detector->name(), "clusters", reinterpret_cast<Objects*>(deviceClusters));
123
    }
124
    LOG(DEBUG) << "Made " << deviceClusters->size() << " clusters for device " << m_detector->name();
Morag Williams's avatar
Morag Williams committed
125

126
    return StatusCode::Success;
127
128
129
}

// Check if a pixel touches any of the pixels in a cluster
130
bool Clustering4D::touching(Pixel* neighbour, Cluster* cluster) {
131

Simon Spannagel's avatar
Simon Spannagel committed
132
    bool Touching = false;
133

134
    for(auto pixel : (*cluster->pixels())) {
135
136
        int row_distance = abs(pixel->row() - neighbour->row());
        int col_distance = abs(pixel->column() - neighbour->column());
137
138
139
140
141

        if(row_distance <= neighbour_radius_row && col_distance <= neighbour_radius_col) {
            if(row_distance > 1 || col_distance > 1) {
                cluster->setSplit(true);
            }
Simon Spannagel's avatar
Simon Spannagel committed
142
143
144
            Touching = true;
            break;
        }
145
    }
Simon Spannagel's avatar
Simon Spannagel committed
146
    return Touching;
147
148
149
}

// Check if a pixel is close in time to the pixels of a cluster
150
bool Clustering4D::closeInTime(Pixel* neighbour, Cluster* cluster) {
151

Simon Spannagel's avatar
Simon Spannagel committed
152
    bool CloseInTime = false;
153

Simon Spannagel's avatar
Simon Spannagel committed
154
    Pixels* pixels = cluster->pixels();
155
    for(size_t iPix = 0; iPix < pixels->size(); iPix++) {
156

157
        double timeDifference = abs(neighbour->timestamp() - (*pixels)[iPix]->timestamp());
158
        if(timeDifference < timingCut)
Simon Spannagel's avatar
Simon Spannagel committed
159
160
161
            CloseInTime = true;
    }
    return CloseInTime;
162
163
}

164
void Clustering4D::calculateClusterCentre(Cluster* cluster) {
165

Simon Spannagel's avatar
Simon Spannagel committed
166
167
    // Empty variables to calculate cluster position
    double row(0), column(0), tot(0);
168

Simon Spannagel's avatar
Simon Spannagel committed
169
170
    // Get the pixels on this cluster
    Pixels* pixels = cluster->pixels();
171
    string detectorID = (*pixels)[0]->detectorID();
172
    double timestamp = (*pixels)[0]->timestamp();
173

Simon Spannagel's avatar
Simon Spannagel committed
174
    // Loop over all pixels
175
176
    for(auto& pixel : (*pixels)) {
        double pixelToT = pixel->adc();
Simon Spannagel's avatar
Simon Spannagel committed
177
178
179
180
        if(pixelToT == 0) {
            LOG(DEBUG) << "Pixel with ToT 0!";
            pixelToT = 1;
        }
181
        tot += pixelToT;
182
183
184
185
        row += (pixel->row() * pixelToT);
        column += (pixel->column() * pixelToT);
        if(pixel->timestamp() < timestamp)
            timestamp = pixel->timestamp();
Simon Spannagel's avatar
Simon Spannagel committed
186
187
    }
    // Row and column positions are tot-weighted
188
189
    row /= (tot > 0 ? tot : 1);
    column /= (tot > 0 ? tot : 1);
190
191
192
193
194

    if(detectorID != m_detector->name()) {
        // Should never happen...
        return;
    }
195

Simon Spannagel's avatar
Simon Spannagel committed
196
    // Create object with local cluster position
197
198
    PositionVector3D<Cartesian3D<double>> positionLocal(m_detector->pitch().X() * (column - m_detector->nPixels().X() / 2),
                                                        m_detector->pitch().Y() * (row - m_detector->nPixels().Y() / 2),
Simon Spannagel's avatar
Simon Spannagel committed
199
                                                        0);
Simon Spannagel's avatar
Simon Spannagel committed
200
    // Calculate global cluster position
201
    PositionVector3D<Cartesian3D<double>> positionGlobal = m_detector->localToGlobal(positionLocal);
Simon Spannagel's avatar
Simon Spannagel committed
202
203
204
205
206

    // Set the cluster parameters
    cluster->setRow(row);
    cluster->setColumn(column);
    cluster->setTot(tot);
207
208

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

Simon Spannagel's avatar
Simon Spannagel committed
211
212
    cluster->setTimestamp(timestamp);
    cluster->setDetectorID(detectorID);
213
214
    cluster->setClusterCentre(positionGlobal);
    cluster->setClusterCentreLocal(positionLocal);
215
}