AnalysisDUT.cpp 18.5 KB
Newer Older
1
#include "AnalysisDUT.h"
2

3
#include "objects/Cluster.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
4
#include "objects/Pixel.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
5
#include "objects/Track.hpp"
6
7
8

using namespace corryvreckan;

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

12
13
14
    m_timeCutFrameEdge = m_config.get<double>("time_cut_frameedge", static_cast<double>(Units::convert(20, "ns")));
    spatialCut = m_config.get<double>("spatial_cut", static_cast<double>(Units::convert(50, "um")));
    chi2ndofCut = m_config.get<double>("chi2ndof_cut", 3.);
15
16
}

17
void AnalysisDUT::initialise() {
18

19
20
    hClusterMapAssoc = new TH2F("clusterMapAssoc",
                                "clusterMapAssoc",
21
                                m_detector->nPixels().X(),
22
                                0,
23
24
                                m_detector->nPixels().X(),
                                m_detector->nPixels().Y(),
25
                                0,
26
                                m_detector->nPixels().Y());
27
28
    hClusterSizeMapAssoc = new TProfile2D("clusterSizeMapAssoc",
                                          "clusterSizeMapAssoc",
29
                                          m_detector->nPixels().X(),
30
                                          0,
31
32
                                          m_detector->nPixels().X(),
                                          m_detector->nPixels().Y(),
33
                                          0,
34
                                          m_detector->nPixels().Y(),
35
36
37
38
39
                                          0,
                                          100);

    hClusterToTMapAssoc = new TProfile2D("clusterSizeToTAssoc",
                                         "clusterToTMapAssoc",
40
                                         m_detector->nPixels().X(),
41
                                         0,
42
43
                                         m_detector->nPixels().X(),
                                         m_detector->nPixels().Y(),
44
                                         0,
45
                                         m_detector->nPixels().Y(),
46
47
                                         0,
                                         1000);
48
49

    // Per-pixel histograms
50
51
    hHitMapAssoc = new TH2F("hitMapAssoc",
                            "hitMapAssoc",
52
                            m_detector->nPixels().X(),
53
                            0,
54
55
                            m_detector->nPixels().X(),
                            m_detector->nPixels().Y(),
56
                            0,
57
                            m_detector->nPixels().Y());
58
59
    hHitMapROI = new TH2F("hitMapROI",
                          "hitMapROI",
60
                          m_detector->nPixels().X(),
61
                          0,
62
63
                          m_detector->nPixels().X(),
                          m_detector->nPixels().Y(),
64
                          0,
65
                          m_detector->nPixels().Y());
66
67
68
    hPixelToTAssoc = new TH1F("pixelToTAssoc", "pixelToTAssoc", 32, 0, 31);
    hPixelToTMapAssoc = new TProfile2D("pixelToTMapAssoc",
                                       "pixelToTMapAssoc",
69
                                       m_detector->nPixels().X(),
70
                                       0,
71
72
                                       m_detector->nPixels().X(),
                                       m_detector->nPixels().Y(),
73
                                       0,
74
                                       m_detector->nPixels().Y(),
75
76
77
78
79
80
81
82
                                       0,
                                       255);

    associatedTracksVersusTime = new TH1F("associatedTracksVersusTime", "associatedTracksVersusTime", 300000, 0, 300);
    residualsX = new TH1F("residualsX", "residualsX", 800, -0.1, 0.1);
    residualsY = new TH1F("residualsY", "residualsY", 800, -0.1, 0.1);

    residualsX1pix = new TH1F("residualsX1pix", "residualsX1pix", 400, -0.2, 0.2);
83
    residualsY1pix = new TH1F("residualsY1pix", "residualsY1pix", 400, -0.2, 0.2);
84
    residualsX2pix = new TH1F("residualsX2pix", "residualsX2pix", 400, -0.2, 0.2);
85
    residualsY2pix = new TH1F("residualsY2pix", "residualsY2pix", 400, -0.2, 0.2);
86
87

    clusterTotAssoc = new TH1F("clusterTotAssociated", "clusterTotAssociated", 10000, 0, 10000);
88
    clusterTotAssocNorm = new TH1F("clusterTotAssociatedNormalized", "clusterTotAssociatedNormalized", 10000, 0, 10000);
89
90
    clusterSizeAssoc = new TH1F("clusterSizeAssociated", "clusterSizeAssociated", 30, 0, 30);

91
    // In-pixel studies:
92
93
    auto pitch_x = static_cast<double>(Units::convert(m_detector->pitch().X(), "um"));
    auto pitch_y = static_cast<double>(Units::convert(m_detector->pitch().Y(), "um"));
Simon Spannagel's avatar
Simon Spannagel committed
94
    auto mod_axes = "x_{track} mod " + std::to_string(pitch_x) + "#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;";
95
96

    std::string title = "DUT x resolution;" + mod_axes + "MAD(#Deltax) [#mum]";
97
98
    rmsxvsxmym = new TProfile2D(
        "rmsxvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
99
100

    title = "DUT y resolution;" + mod_axes + "MAD(#Deltay) [#mum]";
101
102
    rmsyvsxmym = new TProfile2D(
        "rmsyvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
103
104

    title = "DUT resolution;" + mod_axes + "MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
105
106
    rmsxyvsxmym = new TProfile2D(
        "rmsxyvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
107
108

    title = "DUT cluster charge map;" + mod_axes + "<cluster charge> [ke]";
109
110
    qvsxmym = new TProfile2D(
        "qvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 250);
111
112

    title = "DUT cluster charge map, Moyal approx;" + mod_axes + "cluster charge MPV [ke]";
113
114
    qMoyalvsxmym = new TProfile2D(
        "qMoyalvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 250);
115
116

    title = "DUT seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
117
118
    pxqvsxmym = new TProfile2D(
        "pxqvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 250);
119
120

    title = "DUT cluster size map;" + mod_axes + "<pixels/cluster>";
121
122
    npxvsxmym = new TProfile2D(
        "npxvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 4.5);
123
124

    title = "DUT 1-pixel cluster map;" + mod_axes + "clusters";
125
126
    npx1vsxmym =
        new TH2F("npx1vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
127
128

    title = "DUT 2-pixel cluster map;" + mod_axes + "clusters";
129
130
    npx2vsxmym =
        new TH2F("npx2vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
131
132

    title = "DUT 3-pixel cluster map;" + mod_axes + "clusters";
133
134
    npx3vsxmym =
        new TH2F("npx3vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
135
136

    title = "DUT 4-pixel cluster map;" + mod_axes + "clusters";
137
138
    npx4vsxmym =
        new TH2F("npx4vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
139

140
    // Efficiency maps
141
142
143
144
145
146
147
148
149
150
    hPixelEfficiencyMap = new TProfile2D("hPixelEfficiencyMap",
                                         "hPixelEfficiencyMap",
                                         static_cast<int>(pitch_x),
                                         0,
                                         pitch_x,
                                         static_cast<int>(pitch_y),
                                         0,
                                         pitch_y,
                                         0,
                                         1);
151
152
    hChipEfficiencyMap = new TProfile2D("hChipEfficiencyMap",
                                        "hChipEfficiencyMap",
153
                                        m_detector->nPixels().X(),
154
                                        0,
155
156
                                        m_detector->nPixels().X(),
                                        m_detector->nPixels().Y(),
157
                                        0,
158
                                        m_detector->nPixels().Y(),
159
160
161
162
163
                                        0,
                                        1);
    hGlobalEfficiencyMap = new TProfile2D("hGlobalEfficiencyMap",
                                          "hGlobalEfficiencyMap",
                                          300,
164
165
                                          -1.5 * m_detector->size().X(),
                                          1.5 * m_detector->size().X(),
166
                                          300,
167
168
                                          -1.5 * m_detector->size().Y(),
                                          1.5 * m_detector->size().Y(),
169
170
                                          0,
                                          1);
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186

    residualsTime = new TH1F("residualsTime", "residualsTime", 20000, -1000, +1000);

    hTrackCorrelationX = new TH1F("hTrackCorrelationX", "hTrackCorrelationX", 4000, -10., 10.);
    hTrackCorrelationY = new TH1F("hTrackCorrelationY", "hTrackCorrelationY", 4000, -10., 10.);
    hTrackCorrelationTime = new TH1F("hTrackCorrelationTime", "hTrackCorrelationTime", 2000000, -5000, 5000);

    residualsTimeVsTime = new TH2F("residualsTimeVsTime", "residualsTimeVsTime", 20000, 0, 200, 1000, -1000, +1000);
    residualsTimeVsSignal = new TH2F("residualsTimeVsSignal", "residualsTimeVsSignal", 20000, 0, 100000, 1000, -1000, +1000);

    hAssociatedTracksGlobalPosition =
        new TH2F("hAssociatedTracksGlobalPosition", "hAssociatedTracksGlobalPosition", 200, -10, 10, 200, -10, 10);
    hUnassociatedTracksGlobalPosition =
        new TH2F("hUnassociatedTracksGlobalPosition", "hUnassociatedTracksGlobalPosition", 200, -10, 10, 200, -10, 10);
}

187
StatusCode AnalysisDUT::run(std::shared_ptr<Clipboard> clipboard) {
188
189

    // Get the telescope tracks from the clipboard
190
191
    Tracks* tracks = reinterpret_cast<Tracks*>(clipboard->get("tracks"));
    if(tracks == nullptr) {
192
        LOG(DEBUG) << "No tracks on the clipboard";
193
        return StatusCode::Success;
194
195
196
197
    }

    // Loop over all tracks
    for(auto& track : (*tracks)) {
198
199
200
        // Flags to select clusters and tracks
        bool has_associated_cluster = false;
        bool is_within_roi = true;
201
202
203
204
205

        LOG(DEBUG) << "Looking at next track";

        // Cut on the chi2/ndof
        if(track->chi2ndof() > chi2ndofCut) {
206
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
207
208
209
210
            continue;
        }

        // Check if it intercepts the DUT
211
212
        auto globalIntercept = m_detector->getIntercept(track);
        auto localIntercept = m_detector->globalToLocal(globalIntercept);
213

214
        if(!m_detector->hasIntercept(track, 1.)) {
215
216
217
218
            LOG(DEBUG) << " - track outside DUT area";
            continue;
        }

219
        // Check that track is within region of interest using winding number algorithm
220
        if(!m_detector->isWithinROI(track)) {
221
            is_within_roi = false;
222
223
224
        }

        // Check that it doesn't go through/near a masked pixel
225
        if(m_detector->hitMasked(track, 1.)) {
226
227
228
229
            LOG(DEBUG) << " - track close to masked pixel";
            continue;
        }

230
        // Discard tracks which are very close to the frame edges
231
232
        if(fabs(track->timestamp() - clipboard->get_persistent("eventEnd")) < m_timeCutFrameEdge) {
            // Late edge - eventEnd points to the end of the frame`
233
            LOG(DEBUG) << " - track close to end of readout frame: "
234
                       << Units::display(fabs(track->timestamp() - clipboard->get_persistent("eventEnd")), {"us", "ns"})
235
236
                       << " at " << Units::display(track->timestamp(), {"us"});
            continue;
237
238
        } else if(fabs(track->timestamp() - clipboard->get_persistent("eventStart")) < m_timeCutFrameEdge) {
            // Early edge - eventStart points to the beginning of the frame
239
            LOG(DEBUG) << " - track close to start of readout frame: "
240
                       << Units::display(fabs(track->timestamp() - clipboard->get_persistent("eventStart")), {"us", "ns"})
241
242
243
244
                       << " at " << Units::display(track->timestamp(), {"us"});
            continue;
        }

245
        // Calculate in-pixel position of track in microns
246
247
248
        auto inpixel = m_detector->inPixel(localIntercept);
        auto xmod = static_cast<double>(Units::convert(inpixel.X(), "um"));
        auto ymod = static_cast<double>(Units::convert(inpixel.Y(), "um"));
249

250
        // Get the DUT clusters from the clipboard
251
        Clusters* clusters = reinterpret_cast<Clusters*>(clipboard->get(m_detector->name(), "clusters"));
252
        if(clusters == nullptr) {
253
254
255
256
257
258
259
260
            LOG(DEBUG) << " - no DUT clusters";
        } else {

            // Loop over all DUT clusters to find matches:
            for(auto* cluster : (*clusters)) {
                LOG(DEBUG) << " - Looking at next DUT cluster";

                // Check distance between track and cluster
261
                ROOT::Math::XYZPoint intercept = track->intercept(cluster->global().z());
262
263

                // Correlation plots
264
265
                hTrackCorrelationX->Fill(intercept.X() - cluster->global().x());
                hTrackCorrelationY->Fill(intercept.Y() - cluster->global().y());
266
267
                hTrackCorrelationTime->Fill(track->timestamp() - cluster->timestamp());

268
                auto associated_clusters = track->associatedClusters();
Simon Spannagel's avatar
Simon Spannagel committed
269
                if(std::find(associated_clusters.begin(), associated_clusters.end(), cluster) == associated_clusters.end()) {
270
                    LOG(DEBUG) << "No associated cluster found";
271
                    hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
272
273
274
                    continue;
                }

275
                LOG(DEBUG) << "Found associated cluster";
276
277
                double xdistance = intercept.X() - cluster->global().x();
                double ydistance = intercept.Y() - cluster->global().y();
278
279
                double xabsdistance = fabs(xdistance);
                double yabsdistance = fabs(ydistance);
280
                double tdistance = track->timestamp() - cluster->timestamp();
281
282

                // We now have an associated cluster
283
                has_associated_cluster = true;
284
                // FIXME need to understand local coord of clusters - why shifted? what's normal?
285
286
287
288
289
290
291
                auto clusterLocal = m_detector->globalToLocal(cluster->global());
                hClusterMapAssoc->Fill(m_detector->getColumn(clusterLocal), m_detector->getRow(clusterLocal));
                hClusterSizeMapAssoc->Fill(m_detector->getColumn(clusterLocal),
                                           m_detector->getRow(clusterLocal),
                                           static_cast<double>(cluster->size()));
                hClusterToTMapAssoc->Fill(
                    m_detector->getColumn(clusterLocal), m_detector->getRow(clusterLocal), cluster->tot());
292

293
294
                clusterTotAssoc->Fill(cluster->tot());

295
296
297
298
299
                // Cluster charge normalized to path length in sensor:
                double norm = 1; // FIXME fabs(cos( turn*wt )) * fabs(cos( tilt*wt ));
                auto normalized_charge = cluster->tot() * norm;
                clusterTotAssocNorm->Fill(normalized_charge);

300
301
302
                // Fill per-pixel histograms
                for(auto& pixel : (*cluster->pixels())) {
                    hHitMapAssoc->Fill(pixel->column(), pixel->row());
303
304
305
                    if(is_within_roi) {
                        hHitMapROI->Fill(pixel->column(), pixel->row());
                    }
306
307
308
309
                    hPixelToTAssoc->Fill(pixel->tot());
                    hPixelToTMapAssoc->Fill(pixel->column(), pixel->row(), pixel->tot());
                }

310
                associatedTracksVersusTime->Fill(static_cast<double>(Units::convert(track->timestamp(), "s")));
311
312
313
314
315
316
317
318
319
320
321
322
323
324

                // Residuals
                residualsX->Fill(xdistance);
                residualsY->Fill(ydistance);

                if(cluster->size() == 1) {
                    residualsX1pix->Fill(xdistance);
                    residualsY1pix->Fill(ydistance);
                }
                if(cluster->size() == 2) {
                    residualsX2pix->Fill(xdistance);
                    residualsY2pix->Fill(ydistance);
                }

325
326
327
328
329
                // Time residuals
                residualsTime->Fill(tdistance);
                residualsTimeVsTime->Fill(tdistance, track->timestamp());
                residualsTimeVsSignal->Fill(tdistance, cluster->tot());

330
                clusterSizeAssoc->Fill(static_cast<double>(cluster->size()));
331

332
333
334
335
336
337
338
339
340
                // Fill in-pixel plots: (all as function of track position within pixel cell)
                if(is_within_roi) {
                    qvsxmym->Fill(xmod, ymod, cluster->tot());                     // cluster charge profile
                    qMoyalvsxmym->Fill(xmod, ymod, exp(-normalized_charge / 3.5)); // norm. cluster charge profile

                    // mean charge of cluster seed
                    pxqvsxmym->Fill(xmod, ymod, cluster->getSeedPixel()->charge());

                    // mean cluster size
341
                    npxvsxmym->Fill(xmod, ymod, static_cast<double>(cluster->size()));
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
                    if(cluster->size() == 1)
                        npx1vsxmym->Fill(xmod, ymod);
                    if(cluster->size() == 2)
                        npx2vsxmym->Fill(xmod, ymod);
                    if(cluster->size() == 3)
                        npx3vsxmym->Fill(xmod, ymod);
                    if(cluster->size() == 4)
                        npx4vsxmym->Fill(xmod, ymod);

                    // residual MAD x, y, combined (sqrt(x*x + y*y))
                    rmsxvsxmym->Fill(xmod, ymod, xabsdistance);
                    rmsyvsxmym->Fill(xmod, ymod, yabsdistance);
                    rmsxyvsxmym->Fill(xmod, ymod, fabs(sqrt(xdistance * xdistance + ydistance * ydistance)));
                }

357
358
359
360
361
362
363
364
365
                track->addAssociatedCluster(cluster);
                hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());

                // Only allow one associated cluster per track
                break;
            }
        }

        // Efficiency plots:
366
367
        hGlobalEfficiencyMap->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
        hChipEfficiencyMap->Fill(
368
            m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept), has_associated_cluster);
369

370
371
        // For pixels, only look at the ROI:
        if(is_within_roi) {
372
            hPixelEfficiencyMap->Fill(xmod, ymod, has_associated_cluster);
373
        }
374
375
376
    }

    // Return value telling analysis to keep running
377
    return StatusCode::Success;
378
}