AnalysisCLICpix.cpp 29.4 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", Units::get<double>(100, "um"));
    m_proximityCut = m_config.get<double>("proximity_cut", Units::get<double>(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
    hClusterSizeAll = new TH1F("hClusterSizeAll", "hClusterSizeAll", 20, 0, 20);
53
    hClusterChargeAll = new TH1F("hClusterChargeAll", "hClusterChargeAll", 50, 0, 50);
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
    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
    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);

111
112
113
114
115
    hClusterChargeAssociated = new TH1F("hClusterChargeAssociated", "hClusterChargeAssociated", 50, 0, 50);
    hClusterChargeAssociated1pix = new TH1F("hClusterChargeAssociated1pix", "hClusterChargeAssociated1pix", 50, 0, 50);
    hClusterChargeAssociated2pix = new TH1F("hClusterChargeAssociated2pix", "hClusterChargeAssociated2pix", 50, 0, 50);
    hClusterChargeAssociated3pix = new TH1F("hClusterChargeAssociated3pix", "hClusterChargeAssociated3pix", 50, 0, 50);
    hClusterChargeAssociated4pix = new TH1F("hClusterChargeAssociated4pix", "hClusterChargeAssociated4pix", 50, 0, 50);
Simon Spannagel's avatar
Simon Spannagel committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

    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);
133
134
135
136
    hClusterChargeRow2pix = new TH1F("hClusterChargeRow2pix", "hClusterChargeRow2pix", 50, 0, 50);
    hClusterChargeCol2pix = new TH1F("hClusterChargeCol2pix", "hClusterChargeCol2pix", 50, 0, 50);
    hClusterChargeRatioRow2pix = new TH1F("hClusterChargeRatioRow2pix", "hClusterChargeRatioRow2pix", 100, 0, 1);
    hClusterChargeRatioCol2pix = new TH1F("hClusterChargeRatioCol2pix", "hClusterChargeRatioCol2pix", 100, 0, 1);
Simon Spannagel's avatar
Simon Spannagel committed
137
138
139
140
141
142
143
144
145
146
147
148
    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

    for(int x = 0; x < m_nBinsX; x++) {
        for(int y = 0; y < m_nBinsY; y++) {
            int id = x + y * m_nBinsX;
212
213
214
            std::string name = "hMapClusterChargeAssociated1pix" + convertToString(id);
            TH1F* hMapEntryClusterChargeAssociated1pix = new TH1F(name.c_str(), name.c_str(), 50, 0, 50);
            hMapClusterChargeAssociated1pix[id] = hMapEntryClusterChargeAssociated1pix;
Simon Spannagel's avatar
Simon Spannagel committed
215
        }
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
    auto clusters = clipboard->getData<Cluster>(m_detector->name());
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
    auto tracks = clipboard->getData<Track>();
240
    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
        auto 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(!m_detector->hasIntercept(track, 0.5)) {
            LOG(DEBUG) << " - track outside DUT area";
Simon Spannagel's avatar
Simon Spannagel committed
314
            continue;
315
        }
Simon Spannagel's avatar
Simon Spannagel committed
316
317
318

        // Check if the hit is near a masked pixel
        bool hitMasked = checkMasked(chipInterceptRow, chipInterceptCol);
319
320
        if(hitMasked) {
            LOG(DEBUG) << " - hit is masked";
Simon Spannagel's avatar
Simon Spannagel committed
321
            continue;
322
        }
Simon Spannagel's avatar
Simon Spannagel committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338

        // 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
339
            fillResponseHistos(trackIntercept.X(), trackIntercept.Y(), (bestCluster));
Simon Spannagel's avatar
Simon Spannagel committed
340
341
342

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

            // Now fill all of the histograms/efficiency counters that we want
            hTrackInterceptsAssociated->Fill(trackIntercept.X(), trackIntercept.Y());
            hGlobalResidualsX->Fill(xresidualBest);
            hGlobalResidualsY->Fill(yresidualBest);
349
            hGlobalAssociatedClusterPositions->Fill((bestCluster)->global().x(), (bestCluster)->global().y());
Simon Spannagel's avatar
Simon Spannagel committed
350
            nClustersAssociated++;
351
352
            hAssociatedClusterRow->Fill((bestCluster)->row());
            hAssociatedClusterColumn->Fill((bestCluster)->column());
Simon Spannagel's avatar
Simon Spannagel committed
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
            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));
368
            hClusterChargeAssociated->Fill((bestCluster)->charge());
369
            hClusterSizeAssociated->Fill(static_cast<double>(bestCluster->size()));
Simon Spannagel's avatar
Simon Spannagel committed
370
371
372
            //      hClusterWidthColAssociated->Fill((*bestCluster)->colWidth());
            //      hClusterWidthRowAssociated->Fill((*bestCluster)->rowWidth());

373
            hMapClusterSizeAssociated->Fill(chipInterceptCol, chipInterceptRow, (bestCluster)->charge());
Simon Spannagel's avatar
Simon Spannagel committed
374

375
            if((bestCluster)->size() == 1) {
376
                hClusterChargeAssociated1pix->Fill((bestCluster)->charge());
Simon Spannagel's avatar
Simon Spannagel committed
377
                hInterceptClusterSize1->Fill(pixelInterceptX, pixelInterceptY);
378
379
                int id = static_cast<int>(floor(chipInterceptCol * m_nBinsX / m_detector->nPixels().X()) +
                                          floor(chipInterceptRow * m_nBinsY / m_detector->nPixels().Y()) * m_nBinsX);
380
                hMapClusterChargeAssociated1pix[id]->Fill((bestCluster)->charge());
Simon Spannagel's avatar
Simon Spannagel committed
381
            }
382
            if((bestCluster)->size() == 2) {
383
                hClusterChargeAssociated2pix->Fill((bestCluster)->charge());
Simon Spannagel's avatar
Simon Spannagel committed
384
385
                hInterceptClusterSize2->Fill(pixelInterceptX, pixelInterceptY);
            }
386
            if((bestCluster)->size() == 3) {
387
                hClusterChargeAssociated3pix->Fill((bestCluster)->charge());
Simon Spannagel's avatar
Simon Spannagel committed
388
389
                hInterceptClusterSize3->Fill(pixelInterceptX, pixelInterceptY);
            }
390
            if((bestCluster)->size() == 4) {
391
                hClusterChargeAssociated4pix->Fill((bestCluster)->charge());
Simon Spannagel's avatar
Simon Spannagel committed
392
393
                hInterceptClusterSize4->Fill(pixelInterceptX, pixelInterceptY);
            }
394
395
            auto pixels = bestCluster->pixels();
            for(auto& pixel : pixels) {
396
                hPixelToTMap->Fill(pixel->column(), pixel->row(), pixel->raw());
397
            }
398

Simon Spannagel's avatar
Simon Spannagel committed
399
400
401
402
403
404
405
406
        } 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) {
407
        if(pixelMatch) break;
408
        // Loop over pixels
409
        PixelVector::const_iterator itPixel;
410
        for (itPixel = (*itCorrelate)->pixels().begin(); itPixel != (*itCorrelate)->pixels().end(); itPixel++) {
411

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

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

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

Simon Spannagel's avatar
Simon Spannagel committed
432
433
434
            if(!pixelMatch)
                hTrackInterceptsChipUnassociated->Fill(chipInterceptCol, chipInterceptRow);
        }
435
    }
436

Simon Spannagel's avatar
Simon Spannagel committed
437
438
439
440
441
442
    // 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);
443

Simon Spannagel's avatar
Simon Spannagel committed
444
445
446
    // Increment event counter and trigger number
    m_eventNumber++;
    m_triggerNumber++;
447

Simon Spannagel's avatar
Simon Spannagel committed
448
    // Return value telling analysis to keep running
449
    return StatusCode::Success;
450
451
}

452
void AnalysisCLICpix::finalise() {
453

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

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

Simon Spannagel's avatar
Simon Spannagel committed
463
464
465
466
    if(nTracks == 0) {
        efficiency = 0.;
        errorEfficiency = 1.;
    }
Daniel Hynds's avatar
Daniel Hynds committed
467

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

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

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

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

    // If not masked
    return false;
504
505
}

506
507
// Check if there is another track close to the selected track.
// "Close" is defined as the intercept at the clicpix
508
bool AnalysisCLICpix::checkProximity(Track* track, std::shared_ptr<TrackVector> tracks) {
Simon Spannagel's avatar
Simon Spannagel committed
509
510
511

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

    // Loop over all other tracks and check if they intercept close to the track
    // we are considering
516
    TrackVector::iterator itTrack;
Simon Spannagel's avatar
Simon Spannagel committed
517
518
519
520
521
    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)
522
        PositionVector3D<Cartesian3D<double>> trackIntercept2 = m_detector->getIntercept(track2);
Simon Spannagel's avatar
Simon Spannagel committed
523
524
525
526
527
528
529
530
        // 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;
531
532
}

Daniel Hynds's avatar
Daniel Hynds committed
533
// Small sub-routine to fill histograms that only need clusters
534
void AnalysisCLICpix::fillClusterHistos(std::shared_ptr<ClusterVector> clusters) {
Simon Spannagel's avatar
Simon Spannagel committed
535
536

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

    // 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
547
548
        auto pixels = (*itc)->pixels();
        for(auto& itp : pixels) {
Simon Spannagel's avatar
Simon Spannagel committed
549
            // Check if this clicpix frame is still the current
550
551
            int pixelID = itp->column() + nCols * itp->row();
            if(m_hitPixels[pixelID] != itp->raw()) {
Simon Spannagel's avatar
Simon Spannagel committed
552
553
554
555
556
                // New frame! Reset the stored pixels and trigger number
                if(!newFrame) {
                    m_hitPixels.clear();
                    newFrame = true;
                }
557
                m_hitPixels[pixelID] = itp->raw();
Simon Spannagel's avatar
Simon Spannagel committed
558
559
                m_triggerNumber = 0;
            }
560
561
562
            hHitPixels->Fill(itp->column(), itp->row());
            hColumnHits->Fill(itp->column());
            hRowHits->Fill(itp->row());
Simon Spannagel's avatar
Simon Spannagel committed
563
        }
Daniel Hynds's avatar
Daniel Hynds committed
564

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

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

Simon Spannagel's avatar
Simon Spannagel committed
574
    return;
Daniel Hynds's avatar
Daniel Hynds committed
575
}
576

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

Simon Spannagel's avatar
Simon Spannagel committed
581
582
    // Loop over pixels in the cluster and show their distance from the track
    // intercept
583
584
    auto pixels = cluster->pixels();
    for(auto& pixel : pixels) {
Simon Spannagel's avatar
Simon Spannagel committed
585
        // Get the pixel local then global position
586
        PositionVector3D<Cartesian3D<double>> pixelPositionLocal =
587
            m_detector->getLocalPosition(pixel->column(), pixel->row());
588
        PositionVector3D<Cartesian3D<double>> pixelPositionGlobal = m_detector->localToGlobal(pixelPositionLocal);
Simon Spannagel's avatar
Simon Spannagel committed
589
590
591
592
593

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

Simon Spannagel's avatar
Simon Spannagel committed
595
    return;
596
}