AnalysisCLICpix.cpp 29.5 KB
Newer Older
1
#include "AnalysisCLICpix.h"
2
#include "objects/Cluster.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
3
#include "objects/Pixel.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
4
#include "objects/Track.hpp"
5

6
using namespace corryvreckan;
7
using namespace std;
8

9
10
AnalysisCLICpix::AnalysisCLICpix(Configuration config, std::shared_ptr<Detector> detector)
    : Module(std::move(config), detector), m_detector(detector) {
11

12
13
    m_associationCut = m_config.get<double>("association_cut", static_cast<double>(Units::convert(100, "um")));
    m_proximityCut = m_config.get<double>("proximity_cut", static_cast<double>(Units::convert(125, "um")));
14
    timepix3Telescope = m_config.get<bool>("timepix3Telescope", false);
15
16
}

Simon Spannagel's avatar
Simon Spannagel committed
17
18
19
20
template <typename T> std::string convertToString(T number) {
    std::ostringstream ss;
    ss << number;
    return ss.str();
21
}
22

23
void AnalysisCLICpix::initialise() {
Simon Spannagel's avatar
Simon Spannagel committed
24
25
26
27
28
    // Initialise member variables
    m_eventNumber = 0;
    m_triggerNumber = 0;
    m_lostHits = 0.;

29
    int maxcounter = 256;
30

Simon Spannagel's avatar
Simon Spannagel committed
31
    // Cluster/pixel histograms
32
33
    hHitPixels = new TH2F("hHitPixels",
                          "hHitPixels",
34
                          m_detector->nPixels().X(),
35
                          0,
36
37
                          m_detector->nPixels().X(),
                          m_detector->nPixels().Y(),
38
                          0,
39
40
41
                          m_detector->nPixels().Y());
    hColumnHits = new TH1F("hColumnHits", "hColumnHits", m_detector->nPixels().X(), 0, m_detector->nPixels().X());
    hRowHits = new TH1F("hRowHits", "hRowHits", m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
42
43
    hPixelToTMap = new TProfile2D("pixelToTMap",
                                  "pixelToTMap",
44
                                  m_detector->nPixels().X(),
45
                                  0,
46
47
                                  m_detector->nPixels().X(),
                                  m_detector->nPixels().Y(),
48
                                  0,
49
                                  m_detector->nPixels().Y(),
50
51
                                  0,
                                  maxcounter - 1);
Simon Spannagel's avatar
Simon Spannagel committed
52
53
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
    hClusterSizeAll = new TH1F("hClusterSizeAll", "hClusterSizeAll", 20, 0, 20);
    hClusterTOTAll = new TH1F("hClusterTOTAll", "hClusterTOTAll", 50, 0, 50);
    hClustersPerEvent = new TH1F("hClustersPerEvent", "hClustersPerEvent", 500, 0, 500);
    hClustersVersusEventNo = new TH1F("hClustersVersusEventNo", "hClustersVersusEventNo", 60000, 0, 60000);
    hGlobalClusterPositions = new TH2F("hGlobalClusterPositions", "hGlobalClusterPositions", 200, -2.0, 2.0, 300, -1., 2);

    hClusterWidthRow = new TH1F("hClusterWidthRow", "hClusterWidthRow", 20, 0, 20);
    hClusterWidthCol = new TH1F("hClusterWidthCol", "hClusterWidthCol", 20, 0, 20);

    // Track histograms
    hGlobalTrackDifferenceX = new TH1F("hGlobalTrackDifferenceX", "hGlobalTrackDifferenceX", 1000, -10., 10.);
    hGlobalTrackDifferenceY = new TH1F("hGlobalTrackDifferenceY", "hGlobalTrackDifferenceY", 1000, -10., 10.);

    hGlobalResidualsX = new TH1F("hGlobalResidualsX", "hGlobalResidualsX", 600, -0.3, 0.3);
    hGlobalResidualsY = new TH1F("hGlobalResidualsY", "hGlobalResidualsY", 600, -0.3, 0.3);
    hAbsoluteResiduals = new TH1F("hAbsoluteResiduals", "hAbsoluteResiduals", 600, 0., 0.600);
    hGlobalResidualsXversusX =
        new TH2F("hGlobalResidualsXversusX", "hGlobalResidualsXversusX", 600, -3., 3., 600, -0.3, 0.3);
    hGlobalResidualsXversusY =
        new TH2F("hGlobalResidualsXversusY", "hGlobalResidualsXversusY", 600, -3., 3., 600, -0.3, 0.3);
    hGlobalResidualsYversusY =
        new TH2F("hGlobalResidualsYversusY", "hGlobalResidualsYversusY", 600, -3., 3., 600, -0.3, 0.3);
    hGlobalResidualsYversusX =
        new TH2F("hGlobalResidualsYversusX", "hGlobalResidualsYversusX", 600, -3., 3., 600, -0.3, 0.3);

    hGlobalResidualsXversusColWidth =
        new TH2F("hGlobalResidualsXversusColWidth", "hGlobalResidualsXversusColWidth", 30, 0, 30, 600, -0.3, 0.3);
    hGlobalResidualsXversusRowWidth =
        new TH2F("hGlobalResidualsXversusRowWidth", "hGlobalResidualsXversusRowWidth", 30, 0, 30, 600, -0.3, 0.3);
    hGlobalResidualsYversusColWidth =
        new TH2F("hGlobalResidualsYversusColWidth", "hGlobalResidualsYversusColWidth", 30, 0, 30, 600, -0.3, 0.3);
    hGlobalResidualsYversusRowWidth =
        new TH2F("hGlobalResidualsYversusRowWidth", "hGlobalResidualsYversusRowWidth", 30, 0, 30, 600, -0.3, 0.3);

    hTrackInterceptRow = new TH1F("hTrackInterceptRow", "hTrackInterceptRow", 660, -1., 65.);
    hTrackInterceptCol = new TH1F("hTrackInterceptCol", "hTrackInterceptCol", 660, -1., 65.);

    hAbsoluteResidualMap = new TH2F("hAbsoluteResidualMap", "hAbsoluteResidualMap", 50, 0, 50, 25, 0, 25);
    hXresidualVersusYresidual =
        new TH2F("hXresidualVersusYresidual", "hXresidualVersusYresidual", 600, -0.3, 0.3, 600, -0.3, 0.3);

    hAssociatedClustersPerEvent = new TH1F("hAssociatedClustersPerEvent", "hAssociatedClustersPerEvent", 50, 0, 50);
    ;
    hAssociatedClustersVersusEventNo =
        new TH1F("hAssociatedClustersVersusEventNo", "hAssociatedClustersVersusEventNo", 60000, 0, 60000);
    hAssociatedClustersVersusTriggerNo =
        new TH1F("hAssociatedClustersVersusTriggerNo", "hAssociatedClustersVersusTriggerNo", 50, 0, 50);
99
    hAssociatedClusterRow =
100
101
102
        new TH1F("hAssociatedClusterRow", "hAssociatedClusterRow", m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
    hAssociatedClusterColumn = new TH1F(
        "hAssociatedClusterColumn", "hAssociatedClusterColumn", m_detector->nPixels().X(), 0, m_detector->nPixels().X());
Simon Spannagel's avatar
Simon Spannagel committed
103
104
105
106
107
108
109
110
111
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
138
139
140
141
142
143
144
145
146
147
148
    hFrameEfficiency = new TH1F("hFrameEfficiency", "hFrameEfficiency", 6000, 0, 6000);
    hFrameTracks = new TH1F("hFrameTracks", "hFrameTracks", 6000, 0, 6000);
    hFrameTracksAssociated = new TH1F("hFrameTracksAssociated", "hFrameTracksAssociated", 6000, 0, 6000);

    hClusterSizeAssociated = new TH1F("hClusterSizeAssociated", "hClusterSizeAssociated", 20, 0, 20);
    hClusterWidthRowAssociated = new TH1F("hClusterWidthRowAssociated", "hClusterWidthRowAssociated", 20, 0, 20);
    hClusterWidthColAssociated = new TH1F("hClusterWidthColAssociated", "hClusterWidthColAssociated", 20, 0, 20);

    hClusterTOTAssociated = new TH1F("hClusterTOTAssociated", "hClusterTOTAssociated", 50, 0, 50);
    hClusterTOTAssociated1pix = new TH1F("hClusterTOTAssociated1pix", "hClusterTOTAssociated1pix", 50, 0, 50);
    hClusterTOTAssociated2pix = new TH1F("hClusterTOTAssociated2pix", "hClusterTOTAssociated2pix", 50, 0, 50);
    hClusterTOTAssociated3pix = new TH1F("hClusterTOTAssociated3pix", "hClusterTOTAssociated3pix", 50, 0, 50);
    hClusterTOTAssociated4pix = new TH1F("hClusterTOTAssociated4pix", "hClusterTOTAssociated4pix", 50, 0, 50);

    hPixelResponseX = new TH1F("hPixelResponseX", "hPixelResponseX", 600, -0.3, 0.3);
    hPixelResponseY = new TH1F("hPixelResponseY", "hPixelResponseY", 600, -0.3, 0.3);

    hPixelResponseGlobalX = new TH1F("hPixelResponseGlobalX", "hPixelResponseGlobalX", 600, -0.3, 0.3);
    hPixelResponseGlobalY = new TH1F("hPixelResponseGlobalY", "hPixelResponseGlobalY", 600, -0.3, 0.3);

    hPixelResponseXOddCol = new TH1F("hPixelResponseXOddCol", "hPixelResponseXOddCol", 600, -0.3, 0.3);
    hPixelResponseXEvenCol = new TH1F("hPixelResponseXEvenCol", "hPixelResponseXEvenCol", 600, -0.3, 0.3);
    hPixelResponseYOddCol = new TH1F("hPixelResponseYOddCol", "hPixelResponseYOddCol", 600, -0.3, 0.3);
    hPixelResponseYEvenCol = new TH1F("hPixelResponseYEvenCol", "hPixelResponseYEvenCol", 600, -0.3, 0.3);

    hEtaDistributionX = new TH2F("hEtaDistributionX", "hEtaDistributionX", 25, 0, 25, 25, 0, 25);
    hEtaDistributionY = new TH2F("hEtaDistributionY", "hEtaDistributionY", 25, 0, 25, 25, 0, 25);

    hResidualsLocalRow2pix = new TH1F("hResidualsLocalRow2pix", "hResidualsLocalRow2pix", 600, -0.3, 0.3);
    hResidualsLocalCol2pix = new TH1F("hResidualsLocalCol2pix", "hResidualsLocalCol2pix", 600, -0.3, 0.3);
    hClusterTOTRow2pix = new TH1F("hClusterTOTRow2pix", "hClusterTOTRow2pix", 50, 0, 50);
    hClusterTOTCol2pix = new TH1F("hClusterTOTCol2pix", "hClusterTOTCol2pix", 50, 0, 50);
    hClusterTOTRatioRow2pix = new TH1F("hClusterTOTRatioRow2pix", "hClusterTOTRatioRow2pix", 100, 0, 1);
    hClusterTOTRatioCol2pix = new TH1F("hClusterTOTRatioCol2pix", "hClusterTOTRatioCol2pix", 100, 0, 1);
    hPixelTOTRow2pix = new TH1F("hPixelTOTRow2pix", "hPixelTOTRow2pix", 50, 0, 50);
    hPixelTOTCol2pix = new TH1F("hPixelTOTCol2pix", "hPixelTOTCol2pix", 50, 0, 50);

    // Maps
    hTrackIntercepts = new TH2F("hTrackIntercepts", "hTrackIntercepts", 200, -2.0, 2.0, 300, -1., 2);
    hTrackInterceptsAssociated =
        new TH2F("hTrackInterceptsAssociated", "hTrackInterceptsAssociated", 200, -2.0, 2.0, 300, -1., 2);
    hGlobalAssociatedClusterPositions =
        new TH2F("hGlobalAssociatedClusterPositions", "hGlobalAssociatedClusterPositions", 200, -2.0, 2.0, 300, -1., 2);
    hTrackInterceptsPixel = new TH2F("hTrackInterceptsPixel", "hTrackInterceptsPixel", 50, 0, 50, 25, 0, 25);
    hTrackInterceptsPixelAssociated =
        new TH2F("hTrackInterceptsPixelAssociated", "hTrackInterceptsPixelAssociated", 50, 0, 50, 25, 0, 25);
149
150
    hTrackInterceptsChip = new TH2F("hTrackInterceptsChip",
                                    "hTrackInterceptsChip",
151
                                    m_detector->nPixels().X() + 1,
152
                                    -0.5,
153
154
                                    m_detector->nPixels().X() + 0.5,
                                    m_detector->nPixels().Y() + 1,
155
                                    -0.5,
156
                                    m_detector->nPixels().Y() + 0.5);
157
158
    hTrackInterceptsChipAssociated = new TH2F("hTrackInterceptsChipAssociated",
                                              "hTrackInterceptsChipAssociated",
159
                                              m_detector->nPixels().X() + 1,
160
                                              -0.5,
161
162
                                              m_detector->nPixels().X() + 0.5,
                                              m_detector->nPixels().Y() + 1,
163
                                              -0.5,
164
                                              m_detector->nPixels().Y() + 0.5);
165
166
    hTrackInterceptsChipUnassociated = new TH2F("hTrackInterceptsChipUnassociated",
                                                "hTrackInterceptsChipUnassociated",
167
                                                m_detector->nPixels().X() + 1,
168
                                                -0.5,
169
170
                                                m_detector->nPixels().X() + 0.5,
                                                m_detector->nPixels().Y() + 1,
171
                                                -0.5,
172
                                                m_detector->nPixels().Y() + 0.5);
173
174
    hTrackInterceptsChipLost = new TH2F("hTrackInterceptsChipLost",
                                        "hTrackInterceptsChipLost",
175
                                        m_detector->nPixels().X() + 1,
176
                                        -0.5,
177
178
                                        m_detector->nPixels().X() + 0.5,
                                        m_detector->nPixels().Y() + 1,
179
                                        -0.5,
180
                                        m_detector->nPixels().Y() + 0.5);
Simon Spannagel's avatar
Simon Spannagel committed
181
182

    hPixelEfficiencyMap = new TH2F("hPixelEfficiencyMap", "hPixelEfficiencyMap", 50, 0, 50, 25, 0, 25);
183
184
    hChipEfficiencyMap = new TH2F("hChipEfficiencyMap",
                                  "hChipEfficiencyMap",
185
                                  m_detector->nPixels().X() + 1,
186
                                  -0.5,
187
188
                                  m_detector->nPixels().X() + 0.5,
                                  m_detector->nPixels().Y() + 1,
189
                                  -0.5,
190
                                  m_detector->nPixels().Y() + 0.5);
Simon Spannagel's avatar
Simon Spannagel committed
191
192
193
194
195
196
197
198
199
200
201
202
203
    hGlobalEfficiencyMap = new TH2F("hGlobalEfficiencyMap", "hGlobalEfficiencyMap", 200, -2.0, 2.0, 300, -1., 2);

    hInterceptClusterSize1 = new TH2F("hInterceptClusterSize1", "hInterceptClusterSize1", 25, 0, 25, 25, 0, 25);
    hInterceptClusterSize2 = new TH2F("hInterceptClusterSize2", "hInterceptClusterSize2", 25, 0, 25, 25, 0, 25);
    hInterceptClusterSize3 = new TH2F("hInterceptClusterSize3", "hInterceptClusterSize3", 25, 0, 25, 25, 0, 25);
    hInterceptClusterSize4 = new TH2F("hInterceptClusterSize4", "hInterceptClusterSize4", 25, 0, 25, 25, 0, 25);

    m_nBinsX = 32;
    m_nBinsY = 32;
    hMapClusterSizeAssociated = new TH2F("hMapClusterSizeAssociated",
                                         "hMapClusterSizeAssociated",
                                         m_nBinsX,
                                         0,
204
                                         m_detector->nPixels().X(),
Simon Spannagel's avatar
Simon Spannagel committed
205
206
                                         m_nBinsY,
                                         0,
207
                                         m_detector->nPixels().Y());
Simon Spannagel's avatar
Simon Spannagel committed
208
209
210
211
212
213
214
215

    for(int x = 0; x < m_nBinsX; x++) {
        for(int y = 0; y < m_nBinsY; y++) {
            int id = x + y * m_nBinsX;
            std::string name = "hMapClusterTOTAssociated1pix" + convertToString(id);
            TH1F* hMapEntryClusterTOTAssociated1pix = new TH1F(name.c_str(), name.c_str(), 50, 0, 50);
            hMapClusterTOTAssociated1pix[id] = hMapEntryClusterTOTAssociated1pix;
        }
216
    }
217
218
}

219
StatusCode AnalysisCLICpix::run(std::shared_ptr<Clipboard> clipboard) {
220

Simon Spannagel's avatar
Simon Spannagel committed
221
    // Get the clicpix clusters in this event
222
    Clusters* clusters = reinterpret_cast<Clusters*>(clipboard->get(m_detector->name(), "clusters"));
223
    if(clusters == nullptr) {
224
        LOG(DEBUG) << "No clusters for " << m_detector->name() << " on the clipboard";
225
        return StatusCode::Success;
Simon Spannagel's avatar
Simon Spannagel committed
226
    }
227

Simon Spannagel's avatar
Simon Spannagel committed
228
229
230
231
    // If this is the first or last trigger don't use the data
    if((m_triggerNumber == 0 || m_triggerNumber == 19) && !timepix3Telescope) {
        m_eventNumber++;
        m_triggerNumber++;
232
        return StatusCode::Success;
Simon Spannagel's avatar
Simon Spannagel committed
233
    }
234

Simon Spannagel's avatar
Simon Spannagel committed
235
236
    // Fill the histograms that only need clusters/pixels
    fillClusterHistos(clusters);
237

Simon Spannagel's avatar
Simon Spannagel committed
238
    // Get the tracks in this event
239
240
    Tracks* tracks = reinterpret_cast<Tracks*>(clipboard->get("tracks"));
    if(tracks == nullptr) {
Simon Spannagel's avatar
Simon Spannagel committed
241
        LOG(DEBUG) << "No tracks on the clipboard";
242
        return StatusCode::Success;
243
    }
244

Simon Spannagel's avatar
Simon Spannagel committed
245
246
247
248
    // Set counters
    double nClustersAssociated = 0, nValidTracks = 0;

    // Loop over tracks, and compare with Clicpix clusters
249
    for(auto& track : (*tracks)) {
250
        if(!track) {
Simon Spannagel's avatar
Simon Spannagel committed
251
            continue;
252
        }
Simon Spannagel's avatar
Simon Spannagel committed
253
254

        // Cut on the track chi2/ndof
255
        if(track->chi2ndof() < 3.0) {
Simon Spannagel's avatar
Simon Spannagel committed
256
            continue;
257
        }
Simon Spannagel's avatar
Simon Spannagel committed
258
259
260

        // Get the track intercept with the clicpix plane (global and local
        // co-ordinates)
261
262
        PositionVector3D<Cartesian3D<double>> trackIntercept = m_detector->getIntercept(track);
        PositionVector3D<Cartesian3D<double>> trackInterceptLocal = m_detector->globalToLocal(trackIntercept);
Simon Spannagel's avatar
Simon Spannagel committed
263
264
265
266

        // Plot the difference between track intercepts and all clicpix clusters
        // Also record which cluster is the closest
        bool matched = false;
267
        Clusters trackclusters = track->clusters();
268
        Cluster* bestCluster = nullptr;
Simon Spannagel's avatar
Simon Spannagel committed
269
270
271
272
        double xresidualBest = 10000.;
        double yresidualBest = 10000.;
        double absoluteresidualBest = sqrt(xresidualBest * xresidualBest + yresidualBest * yresidualBest);

273
        for(auto& cluster : (trackclusters)) {
Simon Spannagel's avatar
Simon Spannagel committed
274
275

            // Get the distance between this cluster and the track intercept (global)
276
277
            double xcorr = cluster->global().x() - trackIntercept.X();
            double ycorr = cluster->global().y() - trackIntercept.Y();
Simon Spannagel's avatar
Simon Spannagel committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297

            // Fill histograms on correlations
            hGlobalTrackDifferenceX->Fill(xcorr);
            hGlobalTrackDifferenceY->Fill(ycorr);

            // Look if this cluster can be considered associated.
            // Cut on the distance from the track in x and y
            if(fabs(ycorr) > m_associationCut || fabs(xcorr) > m_associationCut)
                continue;

            // Found a matching cluster
            matched = true;

            // Check if this is a better match than the previously found cluster, and
            // store the information
            double absoluteresidual = sqrt(xcorr * xcorr + ycorr * ycorr);
            if(absoluteresidual <= absoluteresidualBest) {
                absoluteresidualBest = absoluteresidual;
                xresidualBest = xcorr;
                yresidualBest = ycorr;
298
                bestCluster = cluster;
Simon Spannagel's avatar
Simon Spannagel committed
299
300
            }
        }
301

Simon Spannagel's avatar
Simon Spannagel committed
302
        // Get the track intercept position along the chip
303
304
        double chipInterceptCol = m_detector->getColumn(trackInterceptLocal);
        double chipInterceptRow = m_detector->getRow(trackInterceptLocal);
Simon Spannagel's avatar
Simon Spannagel committed
305
306

        // Get the track intercept position along the pixel
307
308
        double pixelInterceptX = m_detector->inPixel(trackInterceptLocal).X();
        double pixelInterceptY = m_detector->inPixel(trackInterceptLocal).Y();
Simon Spannagel's avatar
Simon Spannagel committed
309
310
311

        // Cut on the track intercept - this makes sure that it actually went
        // through the chip
312
313
        if(chipInterceptCol < 0.5 || chipInterceptRow < 0.5 || chipInterceptCol > (m_detector->nPixels().X() - 0.5) ||
           chipInterceptRow > (m_detector->nPixels().Y() - 0.5))
Simon Spannagel's avatar
Simon Spannagel committed
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
            continue;

        // Check if the hit is near a masked pixel
        bool hitMasked = checkMasked(chipInterceptRow, chipInterceptCol);
        if(hitMasked)
            continue;

        // Check there are no other tracks nearby
        bool proximityCut = checkProximity(track, tracks);
        if(proximityCut)
            continue;

        // Now have tracks that went through the device
        hTrackIntercepts->Fill(trackIntercept.X(), trackIntercept.Y());
        hFrameTracks->Fill(m_eventNumber);
        nValidTracks++;
        hTrackInterceptsChip->Fill(chipInterceptCol, chipInterceptRow);
        hTrackInterceptsPixel->Fill(pixelInterceptX, pixelInterceptY);

        if(matched) {

            // Some test plots of pixel response function
336
            fillResponseHistos(trackIntercept.X(), trackIntercept.Y(), (bestCluster));
Simon Spannagel's avatar
Simon Spannagel committed
337
338
339

            // Add this cluster to the list of associated clusters held by this track.
            // This will allow alignment to take place
340
            track->addAssociatedCluster(bestCluster);
Simon Spannagel's avatar
Simon Spannagel committed
341
342
343
344
345

            // Now fill all of the histograms/efficiency counters that we want
            hTrackInterceptsAssociated->Fill(trackIntercept.X(), trackIntercept.Y());
            hGlobalResidualsX->Fill(xresidualBest);
            hGlobalResidualsY->Fill(yresidualBest);
346
            hGlobalAssociatedClusterPositions->Fill((bestCluster)->global().x(), (bestCluster)->global().y());
Simon Spannagel's avatar
Simon Spannagel committed
347
            nClustersAssociated++;
348
349
            hAssociatedClusterRow->Fill((bestCluster)->row());
            hAssociatedClusterColumn->Fill((bestCluster)->column());
Simon Spannagel's avatar
Simon Spannagel committed
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
            hTrackInterceptsChipAssociated->Fill(chipInterceptCol, chipInterceptRow);
            hTrackInterceptsPixelAssociated->Fill(pixelInterceptX, pixelInterceptY);
            hFrameTracksAssociated->Fill(m_eventNumber);
            hGlobalResidualsXversusX->Fill(trackIntercept.X(), xresidualBest);
            hGlobalResidualsXversusY->Fill(trackIntercept.Y(), xresidualBest);
            hGlobalResidualsYversusY->Fill(trackIntercept.Y(), yresidualBest);
            hGlobalResidualsYversusX->Fill(trackIntercept.X(), yresidualBest);
            //      hGlobalResidualsXversusColWidth->Fill((*bestCluster)->colWidth(),xresidualBest);
            //      hGlobalResidualsXversusRowWidth->Fill((*bestCluster)->rowWidth(),xresidualBest);
            //      hGlobalResidualsYversusColWidth->Fill((*bestCluster)->colWidth(),yresidualBest);
            //      hGlobalResidualsYversusRowWidth->Fill((*bestCluster)->rowWidth(),yresidualBest);
            hAbsoluteResidualMap->Fill(
                pixelInterceptX, pixelInterceptY, sqrt(xresidualBest * xresidualBest + yresidualBest * yresidualBest));
            hXresidualVersusYresidual->Fill(xresidualBest, yresidualBest);
            hAbsoluteResiduals->Fill(sqrt(xresidualBest * xresidualBest + yresidualBest * yresidualBest));
365
            hClusterTOTAssociated->Fill((bestCluster)->tot());
366
            hClusterSizeAssociated->Fill(static_cast<double>(bestCluster->size()));
Simon Spannagel's avatar
Simon Spannagel committed
367
368
369
            //      hClusterWidthColAssociated->Fill((*bestCluster)->colWidth());
            //      hClusterWidthRowAssociated->Fill((*bestCluster)->rowWidth());

370
            hMapClusterSizeAssociated->Fill(chipInterceptCol, chipInterceptRow, (bestCluster)->tot());
Simon Spannagel's avatar
Simon Spannagel committed
371

372
373
            if((bestCluster)->size() == 1) {
                hClusterTOTAssociated1pix->Fill((bestCluster)->tot());
Simon Spannagel's avatar
Simon Spannagel committed
374
                hInterceptClusterSize1->Fill(pixelInterceptX, pixelInterceptY);
375
376
                int id = static_cast<int>(floor(chipInterceptCol * m_nBinsX / m_detector->nPixels().X()) +
                                          floor(chipInterceptRow * m_nBinsY / m_detector->nPixels().Y()) * m_nBinsX);
377
                hMapClusterTOTAssociated1pix[id]->Fill((bestCluster)->tot());
Simon Spannagel's avatar
Simon Spannagel committed
378
            }
379
380
            if((bestCluster)->size() == 2) {
                hClusterTOTAssociated2pix->Fill((bestCluster)->tot());
Simon Spannagel's avatar
Simon Spannagel committed
381
382
                hInterceptClusterSize2->Fill(pixelInterceptX, pixelInterceptY);
            }
383
            if((bestCluster)->size() == 3) {
384
                hClusterTOTAssociated3pix->Fill((bestCluster)->tot());
Simon Spannagel's avatar
Simon Spannagel committed
385
386
                hInterceptClusterSize3->Fill(pixelInterceptX, pixelInterceptY);
            }
387
388
            if((bestCluster)->size() == 4) {
                hClusterTOTAssociated4pix->Fill((bestCluster)->tot());
Simon Spannagel's avatar
Simon Spannagel committed
389
390
                hInterceptClusterSize4->Fill(pixelInterceptX, pixelInterceptY);
            }
391
392
            Pixels* pixels = bestCluster->pixels();
            for(auto& pixel : (*pixels)) {
393
                hPixelToTMap->Fill(pixel->column(), pixel->row(), pixel->adc());
394
            }
395

Simon Spannagel's avatar
Simon Spannagel committed
396
397
398
399
400
401
402
403
        } else {

            // Search for lost hits. Basically loop through all pixels in all clusters
            // and see if any are close.
            // Large clusters (such as from deltas) can pull the cluster centre
            // sufficiently far from the track
            bool pixelMatch = false;
            /*      for (itCorrelate = clusters->begin(); itCorrelate != clusters->end(); ++itCorrelate) {
404
        if(pixelMatch) break;
405
        // Loop over pixels
406
407
        Pixels::const_iterator itPixel;
        for (itPixel = (*itCorrelate)->pixels().begin(); itPixel != (*itCorrelate)->pixels().end(); itPixel++) {
408

409
          // Get the pixel global position
410
          LOG(TRACE) <<"New pixel, row = "<<(*itPixel)->row()<<", column = "<<(*itPixel)->column();
411
          PositionVector3D<Cartesian3D<double> > pixelPositionLocal =
412
413
      m_detector->getLocalPosition((*itPixel)->row(),(*itPixel)->column()); PositionVector3D<Cartesian3D<double> >
      pixelPositionGlobal = *(m_detector->m_localToGlobal) * pixelPositionLocal;
414

415
416
417
          // Check if it is close to the track
          if( fabs( pixelPositionGlobal.X() - trackIntercept.X() ) > m_associationCut ||
              fabs( pixelPositionGlobal.Y() - trackIntercept.Y() ) > m_associationCut ) continue;
418

419
          pixelMatch=true;
420
          break;
421
422
423
        }
      }
      // Increment counter for number of hits found this way
Simon Spannagel's avatar
Simon Spannagel committed
424
425
426
427
*/ if(pixelMatch) {
                m_lostHits++;
                hTrackInterceptsChipLost->Fill(chipInterceptCol, chipInterceptRow);
            } //*/
428

Simon Spannagel's avatar
Simon Spannagel committed
429
430
431
            if(!pixelMatch)
                hTrackInterceptsChipUnassociated->Fill(chipInterceptCol, chipInterceptRow);
        }
432
    }
433

Simon Spannagel's avatar
Simon Spannagel committed
434
435
436
437
438
439
    // Histograms relating to this mimosa frame
    hAssociatedClustersPerEvent->Fill(nClustersAssociated);
    hAssociatedClustersVersusEventNo->Fill(m_eventNumber, nClustersAssociated);
    hAssociatedClustersVersusTriggerNo->Fill(m_triggerNumber, nClustersAssociated);
    if(nValidTracks != 0)
        hFrameEfficiency->Fill(m_eventNumber, nClustersAssociated / nValidTracks);
440

Simon Spannagel's avatar
Simon Spannagel committed
441
442
443
    // Increment event counter and trigger number
    m_eventNumber++;
    m_triggerNumber++;
444

Simon Spannagel's avatar
Simon Spannagel committed
445
    // Return value telling analysis to keep running
446
    return StatusCode::Success;
447
448
}

449
void AnalysisCLICpix::finalise() {
450

Simon Spannagel's avatar
Simon Spannagel committed
451
    LOG(DEBUG) << "Analysed " << m_eventNumber << " events";
452

Simon Spannagel's avatar
Simon Spannagel committed
453
454
    // Compare number of valid tracks, and number of clusters associated in order
    // to get global efficiency
455
456
    int nTracks = static_cast<int>(hTrackIntercepts->GetEntries());
    int nClusters = static_cast<int>(hGlobalAssociatedClusterPositions->GetEntries());
Simon Spannagel's avatar
Simon Spannagel committed
457
458
    double efficiency = nClusters / nTracks;
    double errorEfficiency = (1. / nTracks) * sqrt(nClusters * (1 - efficiency));
459

Simon Spannagel's avatar
Simon Spannagel committed
460
461
462
463
    if(nTracks == 0) {
        efficiency = 0.;
        errorEfficiency = 1.;
    }
Daniel Hynds's avatar
Daniel Hynds committed
464

Simon Spannagel's avatar
Simon Spannagel committed
465
    LOG(INFO) << "***** Clicpix efficiency calculation *****";
466
    LOG(INFO) << "***** ntracks: " << nTracks << ", nclusters " << nClusters;
Simon Spannagel's avatar
Simon Spannagel committed
467
    LOG(INFO) << "***** Efficiency: " << 100. * efficiency << " +/- " << 100. * errorEfficiency << " %";
468
    LOG(INFO) << "***** If including the " << m_lostHits << " lost pixel hits, this becomes "
Simon Spannagel's avatar
Simon Spannagel committed
469
              << 100. * (m_lostHits + nClusters) / nTracks << " %";
470
}
471
472

// Check if a track has gone through or near a masked pixel
473
bool AnalysisCLICpix::checkMasked(double chipInterceptRow, double chipInterceptCol) {
474

Simon Spannagel's avatar
Simon Spannagel committed
475
    // Get the pixel row and column number
476
477
    int rowNumber = static_cast<int>(ceil(chipInterceptRow));
    int colNumber = static_cast<int>(ceil(chipInterceptCol));
Simon Spannagel's avatar
Simon Spannagel committed
478
479

    // Check if that pixel is masked, or the neighbour to the left, right, etc...
480
    if(m_detector->masked(colNumber - 1, rowNumber - 1))
Simon Spannagel's avatar
Simon Spannagel committed
481
        return true;
482
    if(m_detector->masked(colNumber, rowNumber - 1))
Simon Spannagel's avatar
Simon Spannagel committed
483
        return true;
484
    if(m_detector->masked(colNumber + 1, rowNumber - 1))
Simon Spannagel's avatar
Simon Spannagel committed
485
        return true;
486
    if(m_detector->masked(colNumber - 1, rowNumber))
Simon Spannagel's avatar
Simon Spannagel committed
487
        return true;
488
    if(m_detector->masked(colNumber, rowNumber))
Simon Spannagel's avatar
Simon Spannagel committed
489
        return true;
490
    if(m_detector->masked(colNumber + 1, rowNumber))
Simon Spannagel's avatar
Simon Spannagel committed
491
        return true;
492
    if(m_detector->masked(colNumber - 1, rowNumber + 1))
Simon Spannagel's avatar
Simon Spannagel committed
493
        return true;
494
    if(m_detector->masked(colNumber, rowNumber + 1))
Simon Spannagel's avatar
Simon Spannagel committed
495
        return true;
496
    if(m_detector->masked(colNumber + 1, rowNumber + 1))
Simon Spannagel's avatar
Simon Spannagel committed
497
498
499
500
        return true;

    // If not masked
    return false;
501
502
}

503
504
// Check if there is another track close to the selected track.
// "Close" is defined as the intercept at the clicpix
505
bool AnalysisCLICpix::checkProximity(Track* track, Tracks* tracks) {
Simon Spannagel's avatar
Simon Spannagel committed
506
507
508

    // Get the intercept of the interested track at the dut
    bool close = false;
509
    PositionVector3D<Cartesian3D<double>> trackIntercept = m_detector->getIntercept(track);
Simon Spannagel's avatar
Simon Spannagel committed
510
511
512
513
514
515
516
517
518

    // Loop over all other tracks and check if they intercept close to the track
    // we are considering
    Tracks::iterator itTrack;
    for(itTrack = tracks->begin(); itTrack != tracks->end(); itTrack++) {

        // Get the track
        Track* track2 = (*itTrack);
        // Get the track intercept with the clicpix plane (global co-ordinates)
519
        PositionVector3D<Cartesian3D<double>> trackIntercept2 = m_detector->getIntercept(track2);
Simon Spannagel's avatar
Simon Spannagel committed
520
521
522
523
524
525
526
527
        // If track == track2 do nothing
        if(trackIntercept.X() == trackIntercept2.X() && trackIntercept.Y() == trackIntercept2.Y())
            continue;
        if(fabs(trackIntercept.X() - trackIntercept2.X()) <= m_proximityCut ||
           fabs(trackIntercept.Y() - trackIntercept2.Y()) <= m_proximityCut)
            close = true;
    }
    return close;
528
529
}

Daniel Hynds's avatar
Daniel Hynds committed
530
// Small sub-routine to fill histograms that only need clusters
531
void AnalysisCLICpix::fillClusterHistos(Clusters* clusters) {
Simon Spannagel's avatar
Simon Spannagel committed
532
533

    // Pick up column to generate unique pixel id
534
    int nCols = m_detector->nPixels().X();
Simon Spannagel's avatar
Simon Spannagel committed
535
536
537
538
539
540
541
542
543
544
545
546
547
    Clusters::iterator itc;

    // Check if this is a new clicpix frame (each frame may be in several events)
    // and
    // fill histograms
    bool newFrame = false;
    for(itc = clusters->begin(); itc != clusters->end(); ++itc) {

        // Loop over pixels and check if there are pixels not known
        Pixels* pixels = (*itc)->pixels();
        Pixels::iterator itp;
        for(itp = pixels->begin(); itp != pixels->end(); itp++) {
            // Check if this clicpix frame is still the current
548
549
            int pixelID = (*itp)->column() + nCols * (*itp)->row();
            if(m_hitPixels[pixelID] != (*itp)->adc()) {
Simon Spannagel's avatar
Simon Spannagel committed
550
551
552
553
554
                // New frame! Reset the stored pixels and trigger number
                if(!newFrame) {
                    m_hitPixels.clear();
                    newFrame = true;
                }
555
                m_hitPixels[pixelID] = (*itp)->adc();
Simon Spannagel's avatar
Simon Spannagel committed
556
557
                m_triggerNumber = 0;
            }
558
559
560
            hHitPixels->Fill((*itp)->column(), (*itp)->row());
            hColumnHits->Fill((*itp)->column());
            hRowHits->Fill((*itp)->row());
Simon Spannagel's avatar
Simon Spannagel committed
561
        }
Daniel Hynds's avatar
Daniel Hynds committed
562

Simon Spannagel's avatar
Simon Spannagel committed
563
        // Fill cluster histograms
564
        hClusterSizeAll->Fill(static_cast<double>((*itc)->size()));
Simon Spannagel's avatar
Simon Spannagel committed
565
        hClusterTOTAll->Fill((*itc)->tot());
566
        hGlobalClusterPositions->Fill((*itc)->global().x(), (*itc)->global().y());
Daniel Hynds's avatar
Daniel Hynds committed
567
    }
568

569
570
    hClustersPerEvent->Fill(static_cast<double>(clusters->size()));
    hClustersVersusEventNo->Fill(m_eventNumber, static_cast<double>(clusters->size()));
571

Simon Spannagel's avatar
Simon Spannagel committed
572
    return;
Daniel Hynds's avatar
Daniel Hynds committed
573
}
574

Simon Spannagel's avatar
Simon Spannagel committed
575
576
// Sub-routine to look at pixel response, ie. how far from the pixel is the
// track intercept for the pixel to still see charge
577
void AnalysisCLICpix::fillResponseHistos(double trackInterceptX, double trackInterceptY, Cluster* cluster) {
578

Simon Spannagel's avatar
Simon Spannagel committed
579
580
581
582
583
584
585
586
587
    // Loop over pixels in the cluster and show their distance from the track
    // intercept
    Pixels* pixels = cluster->pixels();
    Pixels::iterator itp;
    for(itp = pixels->begin(); itp != pixels->end(); itp++) {

        // Get the pixel
        Pixel* pixel = (*itp);
        // Get the pixel local then global position
588
589
590
        PositionVector3D<Cartesian3D<double>> pixelPositionLocal =
            m_detector->getLocalPosition(pixel->row(), pixel->column());
        PositionVector3D<Cartesian3D<double>> pixelPositionGlobal = m_detector->localToGlobal(pixelPositionLocal);
Simon Spannagel's avatar
Simon Spannagel committed
591
592
593
594
595

        // Fill the response histograms
        hPixelResponseX->Fill(pixelPositionGlobal.X() - trackInterceptX);
        hPixelResponseY->Fill(pixelPositionGlobal.Y() - trackInterceptY);
    }
596

Simon Spannagel's avatar
Simon Spannagel committed
597
    return;
598
}