Timepix3Clustering.cpp 8.32 KB
Newer Older
1
2
#include "Timepix3Clustering.h"

3
using namespace corryvreckan;
4
using namespace std;
5

6
7
Timepix3Clustering::Timepix3Clustering(Configuration config, std::vector<Detector*> detectors)
    : Algorithm(std::move(config), std::move(detectors)) {
8
    timingCut = m_config.get<double>("timingCut", 0.0000001); // 100 ns
9
10
}

11
void Timepix3Clustering::initialise() {
12

Morag Williams's avatar
Morag Williams committed
13
14
15
16
17
    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);
Daniel Hynds's avatar
Daniel Hynds committed
18
19
20
21
      	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);
Morag Williams's avatar
Morag Williams committed
22
23
24
25
26
        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.);
    }
27
28
}

29
// Sort function for pixels from low to high times
Simon Spannagel's avatar
Simon Spannagel committed
30
31
bool sortByTime(Pixel* pixel1, Pixel* pixel2) {
    return (pixel1->m_timestamp < pixel2->m_timestamp);
32
33
}

Simon Spannagel's avatar
Simon Spannagel committed
34
StatusCode Timepix3Clustering::run(Clipboard* clipboard) {
35

Simon Spannagel's avatar
Simon Spannagel committed
36
    // Loop over all Timepix3 and for each device perform the clustering
37
    for(auto& detector : get_detectors()) {
38

Simon Spannagel's avatar
Simon Spannagel committed
39
        // Check if they are a Timepix3
40
        if(detector->type() != "Timepix3")
Simon Spannagel's avatar
Simon Spannagel committed
41
            continue;
42

Simon Spannagel's avatar
Simon Spannagel committed
43
        // Get the pixels
44
        Pixels* pixels = (Pixels*)clipboard->get(detector->name(), "pixels");
Simon Spannagel's avatar
Simon Spannagel committed
45
        if(pixels == NULL) {
46
            LOG(DEBUG) << "Detector " << detector->name() << " does not have any pixels on the clipboard";
Simon Spannagel's avatar
Simon Spannagel committed
47
48
            continue;
        }
49
        LOG(DEBUG) << "Picked up " << pixels->size() << " pixels for device " << detector->name();
Simon Spannagel's avatar
Simon Spannagel committed
50
51
52

        //    if(pixels->size() > 500.){
        //      LOG(DEBUG) <<"Skipping large event with "<<pixels->size()<<" pixels
53
        //      for device "<<detector->name();
Simon Spannagel's avatar
Simon Spannagel committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
        //      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<Pixel*, bool> 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);
112
        }
113

Simon Spannagel's avatar
Simon Spannagel committed
114
115
        // Put the clusters on the clipboard
        if(deviceClusters->size() > 0)
116
117
            clipboard->put(detector->name(), "clusters", (TestBeamObjects*)deviceClusters);
        LOG(DEBUG) << "Made " << deviceClusters->size() << " clusters for device " << detector->name();
118
    }
119

120
    for(auto& detector : get_detectors()) {
Morag Williams's avatar
Morag Williams committed
121
122
123
124
125
126
127
128
129
130
        // 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());
Daniel Hynds's avatar
Daniel Hynds committed
131
132
          	clusterWidthRow[detector->name()]->Fill(cluster->rowWidth());
          	clusterWidthColumn[detector->name()]->Fill(cluster->columnWidth());
Morag Williams's avatar
Morag Williams committed
133
134
135
            clusterTot[detector->name()]->Fill(cluster->tot());
            clusterPositionGlobal[detector->name()]->Fill(cluster->globalX(), cluster->globalY());
        }
136
    }
Morag Williams's avatar
Morag Williams committed
137

Simon Spannagel's avatar
Simon Spannagel committed
138
    return Success;
139
140
141
}

// Check if a pixel touches any of the pixels in a cluster
Simon Spannagel's avatar
Simon Spannagel committed
142
bool Timepix3Clustering::touching(Pixel* neighbour, Cluster* cluster) {
143

Simon Spannagel's avatar
Simon Spannagel committed
144
145
146
    bool Touching = false;
    Pixels* pixels = cluster->pixels();
    for(int iPix = 0; iPix < pixels->size(); iPix++) {
147

Simon Spannagel's avatar
Simon Spannagel committed
148
149
150
151
152
        if(abs((*pixels)[iPix]->m_row - neighbour->m_row) <= 1 &&
           abs((*pixels)[iPix]->m_column - neighbour->m_column) <= 1) {
            Touching = true;
            break;
        }
153
    }
Simon Spannagel's avatar
Simon Spannagel committed
154
    return Touching;
155
156
157
}

// Check if a pixel is close in time to the pixels of a cluster
Simon Spannagel's avatar
Simon Spannagel committed
158
bool Timepix3Clustering::closeInTime(Pixel* neighbour, Cluster* cluster) {
159

Simon Spannagel's avatar
Simon Spannagel committed
160
    bool CloseInTime = false;
161

Simon Spannagel's avatar
Simon Spannagel committed
162
163
    Pixels* pixels = cluster->pixels();
    for(int iPix = 0; iPix < pixels->size(); iPix++) {
164

Simon Spannagel's avatar
Simon Spannagel committed
165
166
167
168
169
        long long int timeDifference = abs(neighbour->m_timestamp - (*pixels)[iPix]->m_timestamp);
        if(timeDifference < timingCutInt)
            CloseInTime = true;
    }
    return CloseInTime;
170
171
}

Simon Spannagel's avatar
Simon Spannagel committed
172
void Timepix3Clustering::calculateClusterCentre(Cluster* cluster) {
173

Simon Spannagel's avatar
Simon Spannagel committed
174
175
176
    // Empty variables to calculate cluster position
    double row(0), column(0), tot(0);
    long long int timestamp;
177

Simon Spannagel's avatar
Simon Spannagel committed
178
179
180
181
    // Get the pixels on this cluster
    Pixels* pixels = cluster->pixels();
    string detectorID = (*pixels)[0]->m_detectorID;
    timestamp = (*pixels)[0]->m_timestamp;
182

Simon Spannagel's avatar
Simon Spannagel committed
183
184
    // Loop over all pixels
    for(int pix = 0; pix < pixels->size(); pix++) {
Simon Spannagel's avatar
Simon Spannagel committed
185
186
187
188
189
        double pixelToT = (*pixels)[pix]->m_adc;
        if(pixelToT == 0) {
            LOG(DEBUG) << "Pixel with ToT 0!";
            pixelToT = 1;
        }
190
191
192
        tot += pixelToT;
        row += ((*pixels)[pix]->m_row * pixelToT);
        column += ((*pixels)[pix]->m_column * pixelToT);
Simon Spannagel's avatar
Simon Spannagel committed
193
194
195
196
197
198
199
        if((*pixels)[pix]->m_timestamp < timestamp)
            timestamp = (*pixels)[pix]->m_timestamp;
    }
    // Row and column positions are tot-weighted
    row /= tot;
    column /= tot;

200
201
    auto detector = get_detector(detectorID);

Simon Spannagel's avatar
Simon Spannagel committed
202
203
    // Create object with local cluster position
    PositionVector3D<Cartesian3D<double>> positionLocal(
204
        detector->pitchX() * (column - detector->nPixelsX() / 2), detector->pitchY() * (row - detector->nPixelsY() / 2), 0);
Simon Spannagel's avatar
Simon Spannagel committed
205
    // Calculate global cluster position
206
    PositionVector3D<Cartesian3D<double>> positionGlobal = *(detector->m_localToGlobal) * positionLocal;
Simon Spannagel's avatar
Simon Spannagel committed
207
208
209
210
211
212
213
214
215
216

    // 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());
217
}
Simon Spannagel's avatar
Simon Spannagel committed
218

Morag Williams's avatar
Morag Williams committed
219
void Timepix3Clustering::finalise() {}