AnalysisDUT.cpp 32.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
/**
 * @file
 * @brief Implementation of module AnalysisDUT
 *
 * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors.
 * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
 * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
 * Intergovernmental Organization or submit itself to any jurisdiction.
 */

11
#include "AnalysisDUT.h"
12

13
#include "objects/Cluster.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
14
#include "objects/Pixel.hpp"
Simon Spannagel's avatar
Simon Spannagel committed
15
#include "objects/Track.hpp"
16
17
18

using namespace corryvreckan;

19
AnalysisDUT::AnalysisDUT(Configuration& config, std::shared_ptr<Detector> detector)
20
    : Module(config, detector), m_detector(detector) {
21

22
23
24
    config_.setDefault<double>("time_cut_frameedge", Units::get<double>(20, "ns"));
    config_.setDefault<double>("chi2ndof_cut", 3.);
    config_.setDefault<bool>("use_closest_cluster", true);
25
26
    config_.setDefault<int>("n_time_bins", 20000);
    config_.setDefault<double>("time_binning", Units::get<double>(0.1, "ns"));
27
28
29

    m_timeCutFrameEdge = config_.get<double>("time_cut_frameedge");
    chi2ndofCut = config_.get<double>("chi2ndof_cut");
30
    useClosestCluster = config_.get<bool>("use_closest_cluster");
31
    nTimeBins = config_.get<int>("n_time_bins");
32
    timeBinning = config_.get<double>("time_binning");
33
34
}

35
void AnalysisDUT::initialize() {
36

37
    hClusterMapAssoc = new TH2F("clusterMapAssoc",
38
                                "Map of associated clusters; cluster col; cluster row",
39
                                m_detector->nPixels().X(),
40
41
                                -0.5,
                                m_detector->nPixels().X() - 0.5,
42
                                m_detector->nPixels().Y(),
43
44
                                -0.5,
                                m_detector->nPixels().Y() - 0.5);
45
    hClusterSizeMapAssoc = new TProfile2D("clusterSizeMapAssoc",
46
                                          "Size map for associated clusters;cluster column;cluster row;#entries",
47
                                          m_detector->nPixels().X(),
48
49
                                          -0.5,
                                          m_detector->nPixels().X() - 0.5,
50
                                          m_detector->nPixels().Y(),
51
52
                                          -0.5,
                                          m_detector->nPixels().Y() - 0.5,
53
54
                                          0,
                                          100);
55
    hClusterSizeVsColAssoc = new TProfile("clusterSizeVsColAssoc",
56
                                          "cluster size vs. column for assoc clusters;cluster column;mean cluster size",
57
58
59
60
61
                                          m_detector->nPixels().X(),
                                          -0.5,
                                          m_detector->nPixels().X() - 0.5,
                                          0,
                                          100);
62
63
    hClusterSizeVsRowAssoc = new TProfile("clusterSizeVsRowAssoc",
                                          "cluster size vs. row for assoc clusters;cluster row;mean cluster size",
Jens Kroeger's avatar
Jens Kroeger committed
64
                                          m_detector->nPixels().Y(),
65
                                          -0.5,
Jens Kroeger's avatar
Jens Kroeger committed
66
                                          m_detector->nPixels().Y() - 0.5,
67
68
                                          0,
                                          100);
69
    hClusterChargeMapAssoc = new TProfile2D("clusterChargeMapAssoc",
70
                                            "Charge map for associated clusters; cluster charge [e]; #entries",
71
                                            m_detector->nPixels().X(),
72
73
                                            -0.5,
                                            m_detector->nPixels().X() - 0.5,
74
                                            m_detector->nPixels().Y(),
75
76
                                            -0.5,
                                            m_detector->nPixels().Y() - 0.5,
77
78
                                            0,
                                            500);
79
80
    hClusterChargeVsColAssoc =
        new TProfile("clusterChargeVsColAssoc",
81
                     "cluster charge vs. column for assoc clusters;cluster column;mean cluster charge [e]",
82
83
84
85
86
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     0,
                     100);
87
    hClusterChargeVsRowAssoc = new TProfile("clusterChargeVsRowAssoc",
Jens Kroeger's avatar
Jens Kroeger committed
88
89
                                            "cluster charge vs. row for assoc clusters;cluster row;mean cluster charge [e]",
                                            m_detector->nPixels().Y(),
90
                                            -0.5,
Jens Kroeger's avatar
Jens Kroeger committed
91
                                            m_detector->nPixels().Y() - 0.5,
92
93
                                            0,
                                            100);
94

95
96
    hSeedChargeVsColAssoc =
        new TProfile("seedChargeVsColAssoc",
97
                     "seed pixel charge vs. column for assoc clusters;cluster column;mean seed pixel charge [e]",
98
99
100
101
102
103
104
105
106
107
108
109
110
111
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     0,
                     100);
    hSeedChargeVsRowAssoc =
        new TProfile("seedChargeVsRowAssoc",
                     "seed pixel charge vs. row for assoc clusters;cluster row;mean seed pixel charge [e]",
                     m_detector->nPixels().Y(),
                     -0.5,
                     m_detector->nPixels().Y() - 0.5,
                     0,
                     100);

Lennart Huth's avatar
Lennart Huth committed
112
    hTrackZPosDUT = new TH1F("globalTrackZPosOnDUT",
Simon Spannagel's avatar
Simon Spannagel committed
113
                             "Global z-position of track on the DUT; global z of track intersection [mm]; #entries ",
114
115
116
                             400,
                             m_detector->displacement().z() - 10,
                             m_detector->displacement().z() + 10);
117
    // Per-pixel histograms
118
    hHitMapAssoc = new TH2F("hitMapAssoc",
119
                            "Hit map of pixels from associated clusters; hit column; hit row",
120
                            m_detector->nPixels().X(),
121
122
                            -0.5,
                            m_detector->nPixels().X() - 0.5,
123
                            m_detector->nPixels().Y(),
124
125
                            -0.5,
                            m_detector->nPixels().Y() - 0.5);
126
127
128
129
130
    hPixelRawValueAssoc = new TH1F("pixelRawValueAssoc",
                                   "Charge distribution of pixels from associated clusters;pixel raw value;#entries",
                                   1024,
                                   -0.5,
                                   1023.5);
131
    hPixelRawValueMapAssoc = new TProfile2D("pixelRawValueMapAssoc",
132
                                            "Charge map of pixels from associated clusters;pixel raw values;# entries",
133
                                            m_detector->nPixels().X(),
134
135
                                            -0.5,
                                            m_detector->nPixels().X() - 0.5,
136
                                            m_detector->nPixels().Y(),
137
138
                                            -0.5,
                                            m_detector->nPixels().Y() - 0.5,
139
140
                                            0,
                                            255);
141

142
    associatedTracksVersusTime =
143
        new TH1F("associatedTracksVersusTime", "Associated tracks over time;time [s];# associated tracks", 300000, 0, 300);
Simon Spannagel's avatar
Simon Spannagel committed
144
    residualsX = new TH1F("residualsX", "Residual in X;x_{track}-x_{hit}  [mm];# entries", 800, -0.1, 0.1);
145
146
147
    residualsY = new TH1F("residualsY", "Residual in Y;y_{track}-y_{hit}  [mm];# entries", 800, -0.1, 0.1);
    residualsPos = new TH1F(
        "residualsPos", "Absolute distance between track and hit;|pos_{track}-pos_{hit}|  [mm];# entries", 800, -0.1, 0.1);
148
149
    residualsPosVsresidualsTime =
        new TH2F("residualsPosVsresidualsTime",
150
                 "Time vs. absolute position residuals;time_{track}-time_{hit} [ns];|pos_{track}-pos_{hit}| [mm];# entries",
151
152
153
                 nTimeBins,
                 -nTimeBins / 2. * timeBinning - timeBinning / 2.,
                 nTimeBins / 2. * timeBinning - timeBinning / 2.,
154
155
156
                 800,
                 0.,
                 0.2);
157

158
159
160
161
162
163
164
165
    residualsX1pix =
        new TH1F("residualsX1pix", "Residual for 1-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsY1pix =
        new TH1F("residualsY1pix", "Residual for 1-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsX2pix =
        new TH1F("residualsX2pix", "Residual for 2-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsY2pix =
        new TH1F("residualsY2pix", "Residual for 2-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1);
166
167
168
169
170
171
172
173
174
175
176
177
    residualsX3pix =
        new TH1F("residualsX3pix", "Residual for 3-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsY3pix =
        new TH1F("residualsY3pix", "Residual for 3-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsX4pix =
        new TH1F("residualsX4pix", "Residual for 4-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsY4pix =
        new TH1F("residualsY4pix", "Residual for 4-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsXatLeast5pix = new TH1F(
        "residualsXatLeast5pix", "Residual for >= 5-pixel clusters in X;x_{track}-x_{hit} [mm];# entries", 800, -0.1, 0.1);
    residualsYatLeast5pix = new TH1F(
        "residualsYatLeast5pix", "Residual for >= 5-pixel clusters in Y;y_{track}-y_{hit} [mm];# entries", 800, -0.1, 0.1);
178
179
180
181
182
183

    clusterChargeAssoc = new TH1F("clusterChargeAssociated",
                                  "Charge distribution of associated clusters;cluster charge [e];# entries",
                                  10000,
                                  0,
                                  10000);
184
185
186
187
188
    seedChargeAssoc = new TH1F("seedChargeAssociated",
                               "Charge distribution of seed pixels for associated clusters;seed pixel charge [e];# entries",
                               10000,
                               0,
                               10000);
189
190
    clusterSizeAssoc = new TH1F(
        "clusterSizeAssociated", "Size distribution of associated clusters;cluster size; # entries", 30, -0.5, 29.5);
191
192
193
194
195
196
    clusterSizeAssocNorm =
        new TH1F("clusterSizeAssociatedNormalized",
                 "Normalized size distribution of associated clusters;cluster size;# entries (normalized)",
                 30,
                 0,
                 30);
197
198
199
200
201
202
203
204
205
206
    clusterWidthRowAssoc = new TH1F("clusterWidthRowAssociated",
                                    "Height distribution of associated clusters (rows);cluster size row; # entries",
                                    30,
                                    -0.5,
                                    29.5);
    clusterWidthColAssoc = new TH1F("clusterWidthColAssociated",
                                    "Width distribution of associated clusters (columns);cluster size col; # entries",
                                    30,
                                    -0.5,
                                    29.5);
207

208
    // In-pixel studies:
Jin Zhang's avatar
Jin Zhang committed
209
210
    auto pitch_x = static_cast<double>(Units::convert(m_detector->getPitch().X(), "um"));
    auto pitch_y = static_cast<double>(Units::convert(m_detector->getPitch().Y(), "um"));
211
    std::string mod_axes = "in-pixel x_{track} [#mum];in-pixel y_{track} [#mum];";
212

213
    // cut flow histogram
Jin Zhang's avatar
Jin Zhang committed
214
    std::string title = m_detector->getName() + ": number of tracks discarded by different cuts;cut type;tracks";
215
    hCutHisto = new TH1F("hCutHisto", title.c_str(), 4, 1, 5);
216
217
218
219
    hCutHisto->GetXaxis()->SetBinLabel(1, "High Chi2");
    hCutHisto->GetXaxis()->SetBinLabel(2, "Outside DUT area");
    hCutHisto->GetXaxis()->SetBinLabel(3, "Close to masked pixel");
    hCutHisto->GetXaxis()->SetBinLabel(4, "Close to frame begin/end");
220

221
    title = "Resolution in X;" + mod_axes + "MAD(#Deltax) [#mum]";
222
223
224
225
226
227
228
229
    rmsxvsxmym = new TProfile2D("rmsxvsxmym",
                                title.c_str(),
                                static_cast<int>(pitch_x),
                                -pitch_x / 2.,
                                pitch_x / 2.,
                                static_cast<int>(pitch_y),
                                -pitch_y / 2.,
                                pitch_y / 2.);
230

231
    title = "Resolution in Y;" + mod_axes + "MAD(#Deltay) [#mum]";
232
233
234
235
236
237
238
239
    rmsyvsxmym = new TProfile2D("rmsyvsxmym",
                                title.c_str(),
                                static_cast<int>(pitch_x),
                                -pitch_x / 2.,
                                pitch_x / 2.,
                                static_cast<int>(pitch_y),
                                -pitch_y / 2.,
                                pitch_y / 2.);
240

241
    title = "Resolution;" + mod_axes + "MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
242
243
244
245
246
247
248
249
    rmsxyvsxmym = new TProfile2D("rmsxyvsxmym",
                                 title.c_str(),
                                 static_cast<int>(pitch_x),
                                 -pitch_x / 2.,
                                 pitch_x / 2.,
                                 static_cast<int>(pitch_y),
                                 -pitch_y / 2.,
                                 pitch_y / 2.);
250

251
    title = "Mean cluster charge map;" + mod_axes + "<cluster charge> [ke]";
252
253
254
255
256
257
258
259
260
261
    qvsxmym = new TProfile2D("qvsxmym",
                             title.c_str(),
                             static_cast<int>(pitch_x),
                             -pitch_x / 2.,
                             pitch_x / 2.,
                             static_cast<int>(pitch_y),
                             -pitch_y / 2.,
                             pitch_y / 2.,
                             0,
                             250);
262

263
    title = "Most probable cluster charge map, Moyal approx.;" + mod_axes + "cluster charge MPV [ke]";
264
265
266
267
268
269
270
271
272
273
    qMoyalvsxmym = new TProfile2D("qMoyalvsxmym",
                                  title.c_str(),
                                  static_cast<int>(pitch_x),
                                  -pitch_x / 2.,
                                  pitch_x / 2.,
                                  static_cast<int>(pitch_y),
                                  -pitch_y / 2.,
                                  pitch_y / 2.,
                                  0,
                                  250);
274

275
    title = "Seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
276
277
278
279
280
281
282
283
284
285
    pxqvsxmym = new TProfile2D("pxqvsxmym",
                               title.c_str(),
                               static_cast<int>(pitch_x),
                               -pitch_x / 2.,
                               pitch_x / 2.,
                               static_cast<int>(pitch_y),
                               -pitch_y / 2.,
                               pitch_y / 2.,
                               0,
                               250);
286

287
    title = "Mean cluster size map;" + mod_axes + "<pixels/cluster>";
288
289
290
291
292
293
294
295
296
297
    npxvsxmym = new TProfile2D("npxvsxmym",
                               title.c_str(),
                               static_cast<int>(pitch_x),
                               -pitch_x / 2.,
                               pitch_x / 2.,
                               static_cast<int>(pitch_y),
                               -pitch_y / 2.,
                               pitch_y / 2.,
                               0,
                               4.5);
298

299
    title = "1-pixel cluster map;" + mod_axes + "clusters";
300
301
302
303
304
305
306
307
    npx1vsxmym = new TH2F("npx1vsxmym",
                          title.c_str(),
                          static_cast<int>(pitch_x),
                          -pitch_x / 2.,
                          pitch_x / 2.,
                          static_cast<int>(pitch_y),
                          -pitch_y / 2.,
                          pitch_y / 2.);
308

309
    title = "2-pixel cluster map;" + mod_axes + "clusters";
310
311
312
313
314
315
316
317
    npx2vsxmym = new TH2F("npx2vsxmym",
                          title.c_str(),
                          static_cast<int>(pitch_x),
                          -pitch_x / 2.,
                          pitch_x / 2.,
                          static_cast<int>(pitch_y),
                          -pitch_y / 2.,
                          pitch_y / 2.);
318

319
    title = "3-pixel cluster map;" + mod_axes + "clusters";
320
321
322
323
324
325
326
327
    npx3vsxmym = new TH2F("npx3vsxmym",
                          title.c_str(),
                          static_cast<int>(pitch_x),
                          -pitch_x / 2.,
                          pitch_x / 2.,
                          static_cast<int>(pitch_y),
                          -pitch_y / 2.,
                          pitch_y / 2.);
328

329
    title = "4-pixel cluster map;" + mod_axes + "clusters";
330
331
332
333
334
335
336
337
    npx4vsxmym = new TH2F("npx4vsxmym",
                          title.c_str(),
                          static_cast<int>(pitch_x),
                          -pitch_x / 2.,
                          pitch_x / 2.,
                          static_cast<int>(pitch_y),
                          -pitch_y / 2.,
                          pitch_y / 2.);
338

339
340
341
342
343
    residualsTime = new TH1F("residualsTime",
                             "Time residual;time_{track}-time_{hit} [ns];#entries",
                             nTimeBins,
                             -nTimeBins / 2. * timeBinning - timeBinning / 2.,
                             nTimeBins / 2. * timeBinning - timeBinning / 2.);
344
345

    hTrackCorrelationX =
346
        new TH1F("hTrackCorrelationX", "Track residual X, all clusters;x_{track}-x_{hit} [mm];# entries", 4000, -10., 10.);
347
    hTrackCorrelationY =
348
349
350
351
352
353
354
355
356
357
358
        new TH1F("hTrackCorrelationY", "Track residual Y, all clusters;y_{track}-y_{hit} [mm];# entries", 4000, -10., 10.);
    hTrackCorrelationPos = new TH1F("hTrackCorrelationPos",
                                    "Track residual (absolute), all clusters;|pos_{track}-pos_{hit}| [mm];# entries",
                                    2100,
                                    -1.,
                                    10.);
    hTrackCorrelationTime = new TH1F("hTrackCorrelationTime",
                                     "Track time residual, all clusters;time_{track}-time_{hit} [ns];# entries",
                                     20000,
                                     -5000,
                                     5000);
359
360
    hTrackCorrelationPosVsCorrelationTime =
        new TH2F("hTrackCorrelationPosVsCorrelationTime",
361
                 "Track time vs. distance residual;time_{track}-time_{hit} [ns];|pos_{track}-pos_{hit}| [mm];# entries",
362
363
364
365
366
367
                 20000,
                 -5000,
                 5000,
                 2100,
                 -1.,
                 10.);
368
369

    residualsTimeVsTime = new TH2F("residualsTimeVsTime",
370
                                   "Time residual vs. time;time [s];time_{track}-time_{hit} [ns];# entries",
371
372
373
                                   20000,
                                   0,
                                   200,
374
375
376
                                   nTimeBins,
                                   -nTimeBins / 2. * timeBinning - timeBinning / 2.,
                                   nTimeBins / 2. * timeBinning - timeBinning / 2.);
377
378
379
    residualsTimeVsTot =
        new TH2F("residualsTimeVsTot",
                 "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel ToT [lsb];# entries",
380
381
382
383
384
385
                 nTimeBins,
                 -nTimeBins / 2. * timeBinning - timeBinning / 2.,
                 nTimeBins / 2. * timeBinning - timeBinning / 2.,
                 1000,
                 -0.5,
                 999.5);
386

387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
    residualsTimeVsCol =
        new TH2F("residualsTimeVsCol",
                 "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel column;# entries",
                 nTimeBins,
                 -nTimeBins / 2. * timeBinning - timeBinning / 2.,
                 nTimeBins / 2. * timeBinning - timeBinning / 2.,
                 m_detector->nPixels().X(),
                 -0.5,
                 m_detector->nPixels().X() - 0.5);

    residualsTimeVsRow = new TH2F("residualsTimeVsRow",
                                  "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel row;# entries",
                                  nTimeBins,
                                  -nTimeBins / 2. * timeBinning - timeBinning / 2.,
                                  nTimeBins / 2. * timeBinning - timeBinning / 2.,
                                  m_detector->nPixels().X(),
                                  -0.5,
                                  m_detector->nPixels().X() - 0.5);

406
407
408
    residualsTimeVsSignal =
        new TH2F("residualsTimeVsSignal",
                 "Time residual vs. cluster charge;cluster charge [e];time_{track}-time_{hit} [mm];# entries",
409
410
                 nTimeBins,
                 -nTimeBins / 2. * timeBinning - timeBinning / 2.,
411
412
413
414
                 nTimeBins / 2. * timeBinning - timeBinning / 2.,
                 20000,
                 0,
                 100000);
415
416

    hAssociatedTracksGlobalPosition =
417
        new TH2F("hAssociatedTracksGlobalPosition",
418
                 "Map of associated track positions (global);global intercept x [mm];global intercept y [mm]",
419
420
421
                 400,
                 -20,
                 20,
422
423
424
                 200,
                 -10,
                 10);
425
426
427
428
429
430
431
432
433
    hAssociatedTracksLocalPosition =
        new TH2F("hAssociatedTracksLocalPosition",
                 "Map of associated track positions (local);local intercept x [px];local intercept y [px]",
                 10 * m_detector->nPixels().X(),
                 -0.5,
                 m_detector->nPixels().X() - 0.5,
                 10 * m_detector->nPixels().Y(),
                 -0.5,
                 m_detector->nPixels().Y() - 0.5);
434
435
    hUnassociatedTracksGlobalPosition =
        new TH2F("hUnassociatedTracksGlobalPosition",
436
                 "Map of not associated track positions (global); global intercept x [mm]; global intercept y [mm]",
437
438
439
440
441
442
                 200,
                 -10,
                 10,
                 200,
                 -10,
                 10);
443
444
}

445
StatusCode AnalysisDUT::run(const std::shared_ptr<Clipboard>& clipboard) {
446
447

    // Get the telescope tracks from the clipboard
448
    auto tracks = clipboard->getData<Track>();
449
450

    // Loop over all tracks
451
    for(auto& track : tracks) {
452
453
        // Flags to select clusters and tracks
        bool has_associated_cluster = false;
454
455
456
        LOG(DEBUG) << "Looking at next track";

        // Cut on the chi2/ndof
457
        if(track->getChi2ndof() > chi2ndofCut) {
458
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
459
            hCutHisto->Fill(1);
460
            num_tracks++;
461
462
463
464
            continue;
        }

        // Check if it intercepts the DUT
465
        auto globalIntercept = m_detector->getIntercept(track.get());
466
        auto localIntercept = m_detector->globalToLocal(globalIntercept);
467

468
        if(!m_detector->hasIntercept(track.get(), 0.5)) {
469
            LOG(DEBUG) << " - track outside DUT area";
470
            hCutHisto->Fill(2);
471
            num_tracks++;
472
473
474
            continue;
        }

475
        // Check that track is within region of interest using winding number algorithm
476
        if(!m_detector->isWithinROI(track.get())) {
477
            continue;
478
479
480
        }

        // Check that it doesn't go through/near a masked pixel
481
        if(m_detector->hitMasked(track.get(), 1.)) {
482
            LOG(DEBUG) << " - track close to masked pixel";
483
            hCutHisto->Fill(3);
484
            num_tracks++;
485
486
487
            continue;
        }

Simon Spannagel's avatar
Simon Spannagel committed
488
        // Get the event:
489
        auto event = clipboard->getEvent();
Simon Spannagel's avatar
Simon Spannagel committed
490

491
        // Discard tracks which are very close to the frame edges
Simon Spannagel's avatar
Simon Spannagel committed
492
        if(fabs(track->timestamp() - event->end()) < m_timeCutFrameEdge) {
493
            // Late edge - eventEnd points to the end of the frame`
494
            LOG(DEBUG) << " - track close to end of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
495
496
                       << Units::display(fabs(track->timestamp() - event->end()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
497
            hCutHisto->Fill(4);
498
            num_tracks++;
499
            continue;
Simon Spannagel's avatar
Simon Spannagel committed
500
        } else if(fabs(track->timestamp() - event->start()) < m_timeCutFrameEdge) {
501
            // Early edge - eventStart points to the beginning of the frame
502
            LOG(DEBUG) << " - track close to start of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
503
504
                       << Units::display(fabs(track->timestamp() - event->start()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
505
            hCutHisto->Fill(4);
506
            num_tracks++;
507
508
509
            continue;
        }

510
        // Calculate in-pixel position of track in microns
511
512
513
        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"));
514

515
        // Loop over all associated DUT clusters:
516
        for(auto assoc_cluster : track->getAssociatedClusters(m_detector->getName())) {
517
            LOG(DEBUG) << " - Looking at next associated DUT cluster";
518

519
            // if closest cluster should be used continue if current associated cluster is not the closest one
520
521
            if(useClosestCluster && track->getClosestCluster(m_detector->getName()) != assoc_cluster) {
                continue;
522
            }
523
            has_associated_cluster = true;
524

525
            hTrackZPosDUT->Fill(track->getState(m_detector->getName()).z());
526
            // Check distance between track and cluster
527
            ROOT::Math::XYZPoint intercept = track->getIntercept(assoc_cluster->global().z());
528
529
530
531
532
            double xdistance = intercept.X() - assoc_cluster->global().x();
            double ydistance = intercept.Y() - assoc_cluster->global().y();
            double xabsdistance = fabs(xdistance);
            double yabsdistance = fabs(ydistance);
            double tdistance = track->timestamp() - assoc_cluster->timestamp();
533
534
535
536
537
538
539
540
541
542
            double posDiff =
                sqrt((intercept.X() - assoc_cluster->global().x()) * (intercept.X() - assoc_cluster->global().x()) +
                     (intercept.Y() - assoc_cluster->global().y()) * (intercept.Y() - assoc_cluster->global().y()));

            // Correlation plots
            hTrackCorrelationX->Fill(xdistance);
            hTrackCorrelationY->Fill(ydistance);
            hTrackCorrelationTime->Fill(tdistance);
            hTrackCorrelationPos->Fill(posDiff);
            hTrackCorrelationPosVsCorrelationTime->Fill(track->timestamp() - assoc_cluster->timestamp(), posDiff);
543

544
545
546
            hClusterMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row());
            hClusterSizeMapAssoc->Fill(
                assoc_cluster->column(), assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
547
548
            hClusterSizeVsColAssoc->Fill(assoc_cluster->column(), static_cast<double>(assoc_cluster->size()));
            hClusterSizeVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
549
550
551
552

            // Cluster charge normalized to path length in sensor:
            double norm = 1; // FIXME fabs(cos( turn*wt )) * fabs(cos( tilt*wt ));
            // FIXME: what does this mean? To my understanding we have the correct charge here already...
553
            auto cluster_charge = assoc_cluster->charge();
554
            auto normalized_charge = cluster_charge * norm;
555
556

            // clusterChargeAssoc->Fill(normalized_charge);
557
            clusterChargeAssoc->Fill(cluster_charge);
558
            seedChargeAssoc->Fill(assoc_cluster->getSeedPixel()->charge());
559
            hClusterChargeMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row(), cluster_charge);
560
561
562
563
            hClusterChargeVsColAssoc->Fill(assoc_cluster->column(), cluster_charge);
            hClusterChargeVsRowAssoc->Fill(assoc_cluster->row(), cluster_charge);
            hSeedChargeVsColAssoc->Fill(assoc_cluster->column(), assoc_cluster->getSeedPixel()->charge());
            hSeedChargeVsRowAssoc->Fill(assoc_cluster->row(), assoc_cluster->getSeedPixel()->charge());
564
565

            // Fill per-pixel histograms
566
            for(auto& pixel : assoc_cluster->pixels()) {
567
568
569
570
                hHitMapAssoc->Fill(pixel->column(), pixel->row());
                hPixelRawValueAssoc->Fill(pixel->raw());
                hPixelRawValueMapAssoc->Fill(pixel->column(), pixel->row(), pixel->raw());
            }
571

572
            associatedTracksVersusTime->Fill(static_cast<double>(Units::convert(track->timestamp(), "s")));
573

574
575
576
            // Residuals
            residualsX->Fill(xdistance);
            residualsY->Fill(ydistance);
577
578
            residualsPos->Fill(posDiff);
            residualsPosVsresidualsTime->Fill(tdistance, posDiff);
579

580
            if(assoc_cluster->columnWidth() == 1) {
581
                residualsX1pix->Fill(xdistance);
582
583
            }
            if(assoc_cluster->rowWidth() == 1) {
584
585
                residualsY1pix->Fill(ydistance);
            }
586
587

            if(assoc_cluster->columnWidth() == 2) {
588
                residualsX2pix->Fill(xdistance);
589
590
            }
            if(assoc_cluster->rowWidth() == 2) {
591
592
                residualsY2pix->Fill(ydistance);
            }
593

594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
            if(assoc_cluster->columnWidth() == 3) {
                residualsX3pix->Fill(xdistance);
            }
            if(assoc_cluster->rowWidth() == 3) {
                residualsY3pix->Fill(ydistance);
            }

            if(assoc_cluster->columnWidth() == 4) {
                residualsX4pix->Fill(xdistance);
            }
            if(assoc_cluster->rowWidth() == 4) {
                residualsY4pix->Fill(ydistance);
            }

            if(assoc_cluster->columnWidth() >= 5) {
                residualsXatLeast5pix->Fill(xdistance);
            }
            if(assoc_cluster->rowWidth() >= 5) {
                residualsYatLeast5pix->Fill(ydistance);
            }

615
616
            // Time residuals
            residualsTime->Fill(tdistance);
617
            residualsTimeVsTime->Fill(static_cast<double>(Units::convert(track->timestamp(), "s")), tdistance);
618
            residualsTimeVsTot->Fill(tdistance, assoc_cluster->getSeedPixel()->raw());
619
620
            residualsTimeVsCol->Fill(tdistance, assoc_cluster->getSeedPixel()->column());
            residualsTimeVsRow->Fill(tdistance, assoc_cluster->getSeedPixel()->row());
621
            residualsTimeVsSignal->Fill(tdistance, cluster_charge);
622

623
624
            clusterSizeAssoc->Fill(static_cast<double>(assoc_cluster->size()));
            clusterSizeAssocNorm->Fill(static_cast<double>(assoc_cluster->size()));
Simon Spannagel's avatar
Simon Spannagel committed
625
626
            clusterWidthRowAssoc->Fill(static_cast<double>(assoc_cluster->rowWidth()));
            clusterWidthColAssoc->Fill(static_cast<double>(assoc_cluster->columnWidth()));
627
628

            // Fill in-pixel plots: (all as function of track position within pixel cell)
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
            qvsxmym->Fill(xmod, ymod, cluster_charge);                     // 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, assoc_cluster->getSeedPixel()->charge());

            // mean cluster size
            npxvsxmym->Fill(xmod, ymod, static_cast<double>(assoc_cluster->size()));
            if(assoc_cluster->size() == 1)
                npx1vsxmym->Fill(xmod, ymod);
            if(assoc_cluster->size() == 2)
                npx2vsxmym->Fill(xmod, ymod);
            if(assoc_cluster->size() == 3)
                npx3vsxmym->Fill(xmod, ymod);
            if(assoc_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)));

651
652
            hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
            hAssociatedTracksLocalPosition->Fill(m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept));
653
        }
654
655
        if(!has_associated_cluster) {
            hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
656
        }
657
        num_tracks++;
658
659
    }
    // Return value telling analysis to keep running
660
    return StatusCode::Success;
661
}
662

663
void AnalysisDUT::finalize(const std::shared_ptr<ReadonlyClipboard>&) {
664
    hCutHisto->Scale(1 / double(num_tracks));
665
666
    clusterSizeAssocNorm->Scale(1 / clusterSizeAssoc->Integral());
}