AnalysisDUT.cpp 19.3 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
    m_timeCutFrameEdge = m_config.get<double>("time_cut_frameedge", Units::get<double>(20, "ns"));
    spatialCut = m_config.get<double>("spatial_cut", Units::get<double>(50, "um"));
14
    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
    clusterSizeAssoc = new TH1F("clusterSizeAssociated", "clusterSizeAssociated", 30, 0, 30);
90
    clusterSizeAssocNorm = new TH1F("clusterSizeAssociatedNormalized", "clusterSizeAssociatedNormalized", 30, 0, 30);
91

92
    // In-pixel studies:
93
94
    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
95
    auto mod_axes = "x_{track} mod " + std::to_string(pitch_x) + "#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;";
96
97

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

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

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

    title = "DUT cluster charge map;" + mod_axes + "<cluster charge> [ke]";
110
111
    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);
112
113

    title = "DUT cluster charge map, Moyal approx;" + mod_axes + "cluster charge MPV [ke]";
114
115
    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);
116
117

    title = "DUT seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
118
119
    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);
120
121

    title = "DUT cluster size map;" + mod_axes + "<pixels/cluster>";
122
123
    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);
124
125

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

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

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

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

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

    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);
184
185
186
187
188
189
190
191
    hAssociatedTracksLocalPosition = new TH2F("hAssociatedTracksLocalPosition",
                                              "hAssociatedTracksLocalPosition",
                                              m_detector->nPixels().X(),
                                              0,
                                              m_detector->nPixels().X(),
                                              m_detector->nPixels().Y(),
                                              0,
                                              m_detector->nPixels().Y());
192
193
}

194
StatusCode AnalysisDUT::run(std::shared_ptr<Clipboard> clipboard) {
195
196

    // Get the telescope tracks from the clipboard
197
198
    Tracks* tracks = reinterpret_cast<Tracks*>(clipboard->get("tracks"));
    if(tracks == nullptr) {
199
        LOG(DEBUG) << "No tracks on the clipboard";
200
        return StatusCode::Success;
201
202
203
204
    }

    // Loop over all tracks
    for(auto& track : (*tracks)) {
205
206
207
        // Flags to select clusters and tracks
        bool has_associated_cluster = false;
        bool is_within_roi = true;
208
209
210
211
212

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

        // Cut on the chi2/ndof
        if(track->chi2ndof() > chi2ndofCut) {
213
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
214
215
216
217
            continue;
        }

        // Check if it intercepts the DUT
218
219
        auto globalIntercept = m_detector->getIntercept(track);
        auto localIntercept = m_detector->globalToLocal(globalIntercept);
220

221
        if(!m_detector->hasIntercept(track, 0.5)) {
222
223
224
225
            LOG(DEBUG) << " - track outside DUT area";
            continue;
        }

226
        // Check that track is within region of interest using winding number algorithm
227
        if(!m_detector->isWithinROI(track)) {
228
            is_within_roi = false;
229
230
231
        }

        // Check that it doesn't go through/near a masked pixel
232
        if(m_detector->hitMasked(track, 1.)) {
233
234
235
236
            LOG(DEBUG) << " - track close to masked pixel";
            continue;
        }

Simon Spannagel's avatar
Simon Spannagel committed
237
238
239
        // Get the event:
        auto event = clipboard->get_event();

240
        // Discard tracks which are very close to the frame edges
Simon Spannagel's avatar
Simon Spannagel committed
241
        if(fabs(track->timestamp() - event->end()) < m_timeCutFrameEdge) {
242
            // Late edge - eventEnd points to the end of the frame`
243
            LOG(DEBUG) << " - track close to end of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
244
245
                       << Units::display(fabs(track->timestamp() - event->end()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
246
            continue;
Simon Spannagel's avatar
Simon Spannagel committed
247
        } else if(fabs(track->timestamp() - event->start()) < m_timeCutFrameEdge) {
248
            // Early edge - eventStart points to the beginning of the frame
249
            LOG(DEBUG) << " - track close to start of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
250
251
                       << Units::display(fabs(track->timestamp() - event->start()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
252
253
254
            continue;
        }

255
        // Calculate in-pixel position of track in microns
256
257
258
        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"));
259

260
        // Get the DUT clusters from the clipboard
261
        Clusters* clusters = reinterpret_cast<Clusters*>(clipboard->get(m_detector->name(), "clusters"));
262
        if(clusters == nullptr) {
263
264
265
266
267
268
269
270
            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
271
                ROOT::Math::XYZPoint intercept = track->intercept(cluster->global().z());
272
273

                // Correlation plots
274
275
                hTrackCorrelationX->Fill(intercept.X() - cluster->global().x());
                hTrackCorrelationY->Fill(intercept.Y() - cluster->global().y());
276
277
                hTrackCorrelationTime->Fill(track->timestamp() - cluster->timestamp());

278
                auto associated_clusters = track->associatedClusters();
Simon Spannagel's avatar
Simon Spannagel committed
279
                if(std::find(associated_clusters.begin(), associated_clusters.end(), cluster) == associated_clusters.end()) {
280
                    LOG(DEBUG) << "No associated cluster found";
281
                    hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
282
283
284
                    continue;
                }

285
                LOG(DEBUG) << "Found associated cluster";
286
287
                double xdistance = intercept.X() - cluster->global().x();
                double ydistance = intercept.Y() - cluster->global().y();
288
289
                double xabsdistance = fabs(xdistance);
                double yabsdistance = fabs(ydistance);
290
                double tdistance = track->timestamp() - cluster->timestamp();
291
292

                // We now have an associated cluster
293
                has_associated_cluster = true;
294
                // FIXME need to understand local coord of clusters - why shifted? what's normal?
295
296
297
298
299
300
301
                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());
302

303
304
                clusterTotAssoc->Fill(cluster->tot());

305
306
307
308
309
                // 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);

310
311
312
                // Fill per-pixel histograms
                for(auto& pixel : (*cluster->pixels())) {
                    hHitMapAssoc->Fill(pixel->column(), pixel->row());
313
314
315
                    if(is_within_roi) {
                        hHitMapROI->Fill(pixel->column(), pixel->row());
                    }
316
317
318
319
                    hPixelToTAssoc->Fill(pixel->tot());
                    hPixelToTMapAssoc->Fill(pixel->column(), pixel->row(), pixel->tot());
                }

320
                associatedTracksVersusTime->Fill(static_cast<double>(Units::convert(track->timestamp(), "s")));
321
322
323
324
325
326
327
328
329
330
331
332
333
334

                // 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);
                }

335
336
337
338
339
                // Time residuals
                residualsTime->Fill(tdistance);
                residualsTimeVsTime->Fill(tdistance, track->timestamp());
                residualsTimeVsSignal->Fill(tdistance, cluster->tot());

340
                clusterSizeAssoc->Fill(static_cast<double>(cluster->size()));
341
                clusterSizeAssocNorm->Fill(static_cast<double>(cluster->size()));
342

343
344
345
346
347
348
349
350
351
                // 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
352
                    npxvsxmym->Fill(xmod, ymod, static_cast<double>(cluster->size()));
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
                    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)));
                }

368
369
                track->addAssociatedCluster(cluster);
                hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
370
371
                hAssociatedTracksLocalPosition->Fill(m_detector->getColumn(localIntercept),
                                                     m_detector->getRow(localIntercept));
372
373
374
375
376
377
378

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

        // Efficiency plots:
379
380
        hGlobalEfficiencyMap->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
        hChipEfficiencyMap->Fill(
381
            m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept), has_associated_cluster);
382

383
384
        // For pixels, only look at the ROI:
        if(is_within_roi) {
385
            hPixelEfficiencyMap->Fill(xmod, ymod, has_associated_cluster);
386
        }
387
388
    }
    // Return value telling analysis to keep running
389
    return StatusCode::Success;
390
}
391
392
393
394

void AnalysisDUT::finalise() {
    clusterSizeAssocNorm->Scale(1 / clusterSizeAssoc->Integral());
}