AnalysisDUT.cpp 43 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
    config_.setDefault<bool>("correlations", false);
28

29
30
31
32
33
    time_cut_frameedge_ = config_.get<double>("time_cut_frameedge");
    chi2_ndof_cut_ = config_.get<double>("chi2ndof_cut");
    use_closest_cluster_ = config_.get<bool>("use_closest_cluster");
    n_timebins_ = config_.get<int>("n_time_bins");
    time_binning_ = config_.get<double>("time_binning");
34
    correlations_ = config_.get<bool>("correlations");
35
36
}

37
void AnalysisDUT::initialize() {
38

39
    if(correlations_) {
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
        trackCorrelationX_beforeCuts =
            new TH1F("trackCorrelationX_beforeCuts",
                     "Track residual X (tracks before cuts), all clusters;x_{track}-x_{hit} [#mum];# entries",
                     8000,
                     -1000.5,
                     999.5);
        trackCorrelationY_beforeCuts =
            new TH1F("trackCorrelationY_beforeCuts",
                     "Track residual Y (tracks before cuts), all clusters;y_{track}-y_{hit} [#mum];# entries",
                     8000,
                     -1000.5,
                     999.5);
        trackCorrelationTime_beforeCuts =
            new TH1F("trackCorrelationTime_beforeCuts",
                     "Track time residual (tracks before cuts), all clusters;time_{track}-time_{hit} [ns];# entries",
                     n_timebins_,
                     -(n_timebins_ + 1) / 2. * time_binning_,
                     (n_timebins_ - 1) / 2. * time_binning_);

        trackCorrelationX_afterCuts =
            new TH1F("trackCorrelationX_afterCuts",
                     "Track residual X (tracks after cuts), all clusters;x_{track}-x_{hit} [#mum];# entries",
                     8000,
                     -1000.5,
                     999.5);
        trackCorrelationY_afterCuts =
            new TH1F("trackCorrelationY_afterCuts",
                     "Track residual Y (tracks after cuts), all clusters;y_{track}-y_{hit} [#mum];# entries",
                     8000,
                     -1000.5,
                     999.5);
        trackCorrelationTime_afterCuts =
            new TH1F("trackCorrelationTime_afterCuts",
                     "Track time residual (tracks after cuts), all clusters;time_{track}-time_{hit} [ns];# entries",
                     n_timebins_,
                     -(n_timebins_ + 1) / 2. * time_binning_,
                     (n_timebins_ - 1) / 2. * time_binning_);
77
78
    }

79
    hClusterMapAssoc = new TH2F("clusterMapAssoc",
80
                                "Map of associated clusters; cluster col; cluster row",
81
                                m_detector->nPixels().X(),
82
83
                                -0.5,
                                m_detector->nPixels().X() - 0.5,
84
                                m_detector->nPixels().Y(),
85
86
                                -0.5,
                                m_detector->nPixels().Y() - 0.5);
87
    hClusterSizeMapAssoc = new TProfile2D("clusterSizeMapAssoc",
88
                                          "Size map for associated clusters;cluster column;cluster row;#entries",
89
                                          m_detector->nPixels().X(),
90
91
                                          -0.5,
                                          m_detector->nPixels().X() - 0.5,
92
                                          m_detector->nPixels().Y(),
93
94
                                          -0.5,
                                          m_detector->nPixels().Y() - 0.5,
95
96
                                          0,
                                          100);
97
    hClusterSizeVsColAssoc = new TProfile("clusterSizeVsColAssoc",
98
                                          "cluster size vs. column for assoc clusters;cluster column;mean cluster size",
99
100
101
102
103
                                          m_detector->nPixels().X(),
                                          -0.5,
                                          m_detector->nPixels().X() - 0.5,
                                          0,
                                          100);
104
105
    hClusterSizeVsRowAssoc = new TProfile("clusterSizeVsRowAssoc",
                                          "cluster size vs. row for assoc clusters;cluster row;mean cluster size",
Jens Kroeger's avatar
Jens Kroeger committed
106
                                          m_detector->nPixels().Y(),
107
                                          -0.5,
Jens Kroeger's avatar
Jens Kroeger committed
108
                                          m_detector->nPixels().Y() - 0.5,
109
110
                                          0,
                                          100);
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
    hClusterWidthColVsRowAssoc =
        new TProfile("clusterWidthColVsRowAssoc",
                     "cluster column width vs. row for assoc clusters;cluster row;mean cluster column width",
                     m_detector->nPixels().Y(),
                     -0.5,
                     m_detector->nPixels().Y() - 0.5,
                     0,
                     100);
    hClusterWidthRowVsRowAssoc =
        new TProfile("clusterWidthRowVsRowAssoc",
                     "cluster row width vs. row for assoc clusters;cluster row;mean cluster row width",
                     m_detector->nPixels().Y(),
                     -0.5,
                     m_detector->nPixels().Y() - 0.5,
                     0,
                     100);
127
    hClusterChargeMapAssoc = new TProfile2D("clusterChargeMapAssoc",
128
                                            "Charge map for associated clusters; cluster charge [e]; #entries",
129
                                            m_detector->nPixels().X(),
130
131
                                            -0.5,
                                            m_detector->nPixels().X() - 0.5,
132
                                            m_detector->nPixels().Y(),
133
134
                                            -0.5,
                                            m_detector->nPixels().Y() - 0.5,
135
136
                                            0,
                                            500);
137
138
    hClusterChargeVsColAssoc =
        new TProfile("clusterChargeVsColAssoc",
139
                     "cluster charge vs. column for assoc clusters;cluster column;mean cluster charge [e]",
140
141
142
143
144
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     0,
                     100);
145
    hClusterChargeVsRowAssoc = new TProfile("clusterChargeVsRowAssoc",
Jens Kroeger's avatar
Jens Kroeger committed
146
147
                                            "cluster charge vs. row for assoc clusters;cluster row;mean cluster charge [e]",
                                            m_detector->nPixels().Y(),
148
                                            -0.5,
Jens Kroeger's avatar
Jens Kroeger committed
149
                                            m_detector->nPixels().Y() - 0.5,
150
151
                                            0,
                                            100);
152

153
154
    hSeedChargeVsColAssoc =
        new TProfile("seedChargeVsColAssoc",
155
                     "seed pixel charge vs. column for assoc clusters;cluster column;mean seed pixel charge [e]",
156
157
158
159
160
161
162
163
164
165
166
167
168
                     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);
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
    hClusterChargeVsRowAssoc_2D = new TH2F("hClusterChargeVsRowAssoc_2D",
                                           "cluster charge distribution vs. cluster row; cluster row; cluster ToT [lsb]",
                                           m_detector->nPixels().Y(),
                                           -0.5,
                                           m_detector->nPixels().Y() - 0.5,
                                           64,
                                           -0.5,
                                           63.5);
    hSeedChargeVsRowAssoc_2D = new TH2F("hSeedChargeVsRowAssoc_2D",
                                        "seed charge distribution vs. cluster row; cluster row; seed pixel ToT [lsb]",
                                        m_detector->nPixels().Y(),
                                        -0.5,
                                        m_detector->nPixels().Y() - 0.5,
                                        64,
                                        -0.5,
                                        63.5);
185

Lennart Huth's avatar
Lennart Huth committed
186
    hTrackZPosDUT = new TH1F("globalTrackZPosOnDUT",
Simon Spannagel's avatar
Simon Spannagel committed
187
                             "Global z-position of track on the DUT; global z of track intersection [mm]; #entries ",
188
189
190
                             400,
                             m_detector->displacement().z() - 10,
                             m_detector->displacement().z() + 10);
191
    // Per-pixel histograms
192
    hHitMapAssoc = new TH2F("hitMapAssoc",
193
                            "Hit map of pixels from associated clusters; hit column; hit row",
194
                            m_detector->nPixels().X(),
195
196
                            -0.5,
                            m_detector->nPixels().X() - 0.5,
197
                            m_detector->nPixels().Y(),
198
199
                            -0.5,
                            m_detector->nPixels().Y() - 0.5);
200
201
202
203
204
    hPixelRawValueAssoc = new TH1F("pixelRawValueAssoc",
                                   "Charge distribution of pixels from associated clusters;pixel raw value;#entries",
                                   1024,
                                   -0.5,
                                   1023.5);
205
    hPixelRawValueMapAssoc = new TProfile2D("pixelRawValueMapAssoc",
206
                                            "Charge map of pixels from associated clusters;pixel raw values;# entries",
207
                                            m_detector->nPixels().X(),
208
209
                                            -0.5,
                                            m_detector->nPixels().X() - 0.5,
210
                                            m_detector->nPixels().Y(),
211
212
                                            -0.5,
                                            m_detector->nPixels().Y() - 0.5,
213
214
                                            0,
                                            255);
215

216
    associatedTracksVersusTime =
217
        new TH1F("associatedTracksVersusTime", "Associated tracks over time;time [s];# associated tracks", 300000, 0, 300);
218
219
220
221
222
223
224
225
226
227
228
229
230
231
    residualsX = new TH1F("residualsX", "Residual in X;x_{track}-x_{hit}  [#mum];# entries", 4000, -500.5, 499.5);
    residualsY = new TH1F("residualsY", "Residual in Y;y_{track}-y_{hit}  [#mum];# entries", 4000, -500.5, 499.5);
    residualsPos = new TH1F("residualsPos",
                            "Absolute distance between track and hit;|pos_{track}-pos_{hit}|  [#mum];# entries",
                            4000,
                            -500.5,
                            499.5);
    residualsPosVsresidualsTime = new TH2F(
        "residualsPosVsresidualsTime",
        "Time vs. absolute position residuals;time_{track}-time_{hit} [ns];|pos_{track}-pos_{hit}| [#mum];# entries",
        n_timebins_,
        -(n_timebins_ + 1) / 2. * time_binning_,
        (n_timebins_ - 1) / 2. * time_binning_,
        800,
232
233
        -0.25,
        199.75);
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

    residualsX1pix = new TH1F(
        "residualsX1pix", "Residual for 1-pixel clusters in X;x_{track}-x_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsY1pix = new TH1F(
        "residualsY1pix", "Residual for 1-pixel clusters in Y;y_{track}-y_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsX2pix = new TH1F(
        "residualsX2pix", "Residual for 2-pixel clusters in X;x_{track}-x_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsY2pix = new TH1F(
        "residualsY2pix", "Residual for 2-pixel clusters in Y;y_{track}-y_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsX3pix = new TH1F(
        "residualsX3pix", "Residual for 3-pixel clusters in X;x_{track}-x_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsY3pix = new TH1F(
        "residualsY3pix", "Residual for 3-pixel clusters in Y;y_{track}-y_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsX4pix = new TH1F(
        "residualsX4pix", "Residual for 4-pixel clusters in X;x_{track}-x_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsY4pix = new TH1F(
        "residualsY4pix", "Residual for 4-pixel clusters in Y;y_{track}-y_{hit} [#mum];# entries", 4000, -500.5, 499.5);
    residualsXatLeast5pix = new TH1F("residualsXatLeast5pix",
                                     "Residual for >= 5-pixel clusters in X;x_{track}-x_{hit} [#mum];# entries",
                                     4000,
                                     -500.5,
                                     499.5);
    residualsYatLeast5pix = new TH1F("residualsYatLeast5pix",
                                     "Residual for >= 5-pixel clusters in Y;y_{track}-y_{hit} [#mum];# entries",
                                     4000,
                                     -500.5,
                                     499.5);
261
262
263
264
265
266

    clusterChargeAssoc = new TH1F("clusterChargeAssociated",
                                  "Charge distribution of associated clusters;cluster charge [e];# entries",
                                  10000,
                                  0,
                                  10000);
267
268
269
270
271
    seedChargeAssoc = new TH1F("seedChargeAssociated",
                               "Charge distribution of seed pixels for associated clusters;seed pixel charge [e];# entries",
                               10000,
                               0,
                               10000);
272
273
    clusterSizeAssoc = new TH1F(
        "clusterSizeAssociated", "Size distribution of associated clusters;cluster size; # entries", 30, -0.5, 29.5);
274
275
276
277
278
279
    clusterSizeAssocNorm =
        new TH1F("clusterSizeAssociatedNormalized",
                 "Normalized size distribution of associated clusters;cluster size;# entries (normalized)",
                 30,
                 0,
                 30);
280
281
282
283
284
285
286
287
288
289
    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);
290

291
    // In-pixel studies:
292
293
    auto pitch_x = m_detector->getPitch().X() * 1000.; // convert mm -> um
    auto pitch_y = m_detector->getPitch().Y() * 1000.; // convert mm -> um
294
    std::string mod_axes = "in-pixel x_{track} [#mum];in-pixel y_{track} [#mum];";
295

296
    // cut flow histogram
Jin Zhang's avatar
Jin Zhang committed
297
    std::string title = m_detector->getName() + ": number of tracks discarded by different cuts;cut type;tracks";
298
    hCutHisto = new TH1F("hCutHisto", title.c_str(), 4, 1, 5);
299
300
301
302
    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");
303

304
    title = "Resolution in X;" + mod_axes + "MAD(#Deltax) [#mum]";
305
306
307
308
309
310
311
312
    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.);
313

314
    title = "Resolution in Y;" + mod_axes + "MAD(#Deltay) [#mum]";
315
316
317
318
319
320
321
322
    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.);
323

324
    title = "Resolution;" + mod_axes + "MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
325
326
327
328
329
330
331
332
    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.);
333

334
    title = "Mean cluster charge map;" + mod_axes + "<cluster charge> [ke]";
335
336
337
338
339
340
341
342
343
344
    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);
345

346
    title = "Most probable cluster charge map, Moyal approx.;" + mod_axes + "cluster charge MPV [ke]";
347
348
349
350
351
352
353
354
355
356
    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);
357

358
    title = "Seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
359
360
361
362
363
364
365
366
367
368
    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);
369

370
    title = "Mean cluster size map;" + mod_axes + "<pixels/cluster>";
371
372
373
374
375
376
377
378
379
380
    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);
381

382
    title = "1-pixel cluster map;" + mod_axes + "clusters";
383
384
385
386
387
388
389
390
    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.);
391

392
    title = "2-pixel cluster map;" + mod_axes + "clusters";
393
394
395
396
397
398
399
400
    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.);
401

402
    title = "3-pixel cluster map;" + mod_axes + "clusters";
403
404
405
406
407
408
409
410
    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.);
411

412
    title = "4-pixel cluster map;" + mod_axes + "clusters";
413
414
415
416
417
418
419
420
    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.);
421

422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
    title = "Mean cluster charge map (1-pixel);" + mod_axes + "<cluster charge> [ke]";
    qvsxmym_1px = new TProfile2D("qvsxmym_1px",
                                 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);

    title = "Mean cluster charge map (2-pixel);" + mod_axes + "<cluster charge> [ke]";
    qvsxmym_2px = new TProfile2D("qvsxmym_2px",
                                 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);

    title = "Mean cluster charge map (3-pixel);" + mod_axes + "<cluster charge> [ke]";
    qvsxmym_3px = new TProfile2D("qvsxmym_3px",
                                 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);

    title = "Mean cluster charge map (4-pixel);" + mod_axes + "<cluster charge> [ke]";
    qvsxmym_4px = new TProfile2D("qvsxmym_4px",
                                 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);

470
471
    residualsTime = new TH1F("residualsTime",
                             "Time residual;time_{track}-time_{hit} [ns];#entries",
472
                             n_timebins_,
473
474
                             -(n_timebins_ + 1) / 2. * time_binning_,
                             (n_timebins_ - 1) / 2. * time_binning_);
475
476

    residualsTimeVsTime = new TH2F("residualsTimeVsTime",
477
                                   "Time residual vs. time;time [s];time_{track}-time_{hit} [ns];# entries",
478
479
480
                                   20000,
                                   0,
                                   200,
481
                                   n_timebins_,
482
483
                                   -(n_timebins_ + 1) / 2. * time_binning_,
                                   (n_timebins_ - 1) / 2. * time_binning_);
484
485
486
    residualsTimeVsTot =
        new TH2F("residualsTimeVsTot",
                 "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel ToT [lsb];# entries",
487
                 n_timebins_,
488
489
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
490
491
492
                 1000,
                 -0.5,
                 999.5);
493

494
495
496
    residualsTimeVsCol =
        new TH2F("residualsTimeVsCol",
                 "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel column;# entries",
497
                 n_timebins_,
498
499
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
500
501
502
503
504
505
                 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",
506
                                  n_timebins_,
507
508
                                  -(n_timebins_ + 1) / 2. * time_binning_,
                                  (n_timebins_ - 1) / 2. * time_binning_,
509
510
511
512
                                  m_detector->nPixels().X(),
                                  -0.5,
                                  m_detector->nPixels().X() - 0.5);

513
514
515
    residualsTimeVsSignal =
        new TH2F("residualsTimeVsSignal",
                 "Time residual vs. cluster charge;cluster charge [e];time_{track}-time_{hit} [mm];# entries",
516
                 n_timebins_,
517
518
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
519
520
521
                 20000,
                 0,
                 100000);
522
523

    hAssociatedTracksGlobalPosition =
524
        new TH2F("hAssociatedTracksGlobalPosition",
525
                 "Map of associated track positions (global);global intercept x [mm];global intercept y [mm]",
526
527
528
                 400,
                 -20,
                 20,
529
530
531
                 200,
                 -10,
                 10);
532
533
534
535
536
537
538
539
540
    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);
541
542
    hUnassociatedTracksGlobalPosition =
        new TH2F("hUnassociatedTracksGlobalPosition",
543
                 "Map of not associated track positions (global); global intercept x [mm]; global intercept y [mm]",
544
545
546
547
548
549
                 200,
                 -10,
                 10,
                 200,
                 -10,
                 10);
Jens Kroeger's avatar
Jens Kroeger committed
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591

    pxTimeMinusSeedTime = new TH1F("pxTimeMinusSeedTime",
                                   "pixel - seed pixel timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns];events",
                                   n_timebins_,
                                   -(n_timebins_ + 1) / 2. * time_binning_,
                                   (n_timebins_ - 1) / 2. * time_binning_);
    pxTimeMinusSeedTime_vs_pxCharge =
        new TH2F("pxTimeMinusSeedTime_vs_pxCharge",
                 "pixel - seed pixel timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns]; pixel charge [e];events",
                 n_timebins_,
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
                 256,
                 -0.5,
                 255.5);
    pxTimeMinusSeedTime_vs_pxCharge_2px =
        new TH2F("pxTimeMinusSeedTime_vs_pxCharge_2px",
                 "pixel - seed pixel timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns]; pixel charge [e];events",
                 n_timebins_,
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
                 256,
                 -0.5,
                 255.5);
    pxTimeMinusSeedTime_vs_pxCharge_3px =
        new TH2F("pxTimeMinusSeedTime_vs_pxCharge_3px",
                 "pixel - seed pixel timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns]; pixel charge [e];events",
                 n_timebins_,
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
                 256,
                 -0.5,
                 255.5);
    pxTimeMinusSeedTime_vs_pxCharge_4px =
        new TH2F("pxTimeMinusSeedTime_vs_pxCharge_4px",
                 "pixel - seed pixel timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns]; pixel charge [e];events",
                 n_timebins_,
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
                 256,
                 -0.5,
                 255.5);
592
593
}

594
StatusCode AnalysisDUT::run(const std::shared_ptr<Clipboard>& clipboard) {
595
596

    // Get the telescope tracks from the clipboard
597
    auto tracks = clipboard->getData<Track>();
598
599

    // Loop over all tracks
600
    for(auto& track : tracks) {
601
602
603
        auto globalIntercept = m_detector->getIntercept(track.get());
        auto localIntercept = m_detector->globalToLocal(globalIntercept);

604
605
        // Fill correlation plots BEFORE applying any cuts:
        if(correlations_) {
606
607
            auto clusters = clipboard->getData<Cluster>(m_detector->getName());
            for(auto& cls : clusters) {
608
609
                double xdistance_um = (globalIntercept.X() - cls->global().x()) * 1000.;
                double ydistance_um = (globalIntercept.Y() - cls->global().y()) * 1000.;
610
611
612
                trackCorrelationX_beforeCuts->Fill(xdistance_um);
                trackCorrelationY_beforeCuts->Fill(ydistance_um);
                trackCorrelationTime_beforeCuts->Fill(track->timestamp() - cls->timestamp());
613
614
615
            }
        }

616
617
        // Flags to select clusters and tracks
        bool has_associated_cluster = false;
618
619
620
        LOG(DEBUG) << "Looking at next track";

        // Cut on the chi2/ndof
621
        if(track->getChi2ndof() > chi2_ndof_cut_) {
622
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
623
            hCutHisto->Fill(1);
624
            num_tracks_++;
625
626
627
628
            continue;
        }

        // Check if it intercepts the DUT
629
        if(!m_detector->hasIntercept(track.get(), 0.5)) {
630
            LOG(DEBUG) << " - track outside DUT area";
631
            hCutHisto->Fill(2);
632
            num_tracks_++;
633
634
635
            continue;
        }

636
        // Check that track is within region of interest using winding number algorithm
637
        if(!m_detector->isWithinROI(track.get())) {
638
            continue;
639
640
641
        }

        // Check that it doesn't go through/near a masked pixel
642
        if(m_detector->hitMasked(track.get(), 1.)) {
643
            LOG(DEBUG) << " - track close to masked pixel";
644
            hCutHisto->Fill(3);
645
            num_tracks_++;
646
647
648
            continue;
        }

649
        // Fill correlation plots after applying cuts:
650
        if(correlations_) {
651
652
653
654
            auto clusters = clipboard->getData<Cluster>(m_detector->getName());
            for(auto& cls : clusters) {
                double xdistance_um = (globalIntercept.X() - cls->global().x()) * 1000.;
                double ydistance_um = (globalIntercept.Y() - cls->global().y()) * 1000.;
655
656
657
                trackCorrelationX_afterCuts->Fill(xdistance_um);
                trackCorrelationY_afterCuts->Fill(ydistance_um);
                trackCorrelationTime_afterCuts->Fill(track->timestamp() - cls->timestamp());
658
659
660
            }
        }

Simon Spannagel's avatar
Simon Spannagel committed
661
        // Get the event:
662
        auto event = clipboard->getEvent();
Simon Spannagel's avatar
Simon Spannagel committed
663

664
        // Discard tracks which are very close to the frame edges
665
        if(fabs(track->timestamp() - event->end()) < time_cut_frameedge_) {
666
            // Late edge - eventEnd points to the end of the frame`
667
            LOG(DEBUG) << " - track close to end of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
668
669
                       << Units::display(fabs(track->timestamp() - event->end()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
670
            hCutHisto->Fill(4);
671
            num_tracks_++;
672
            continue;
673
        } else if(fabs(track->timestamp() - event->start()) < time_cut_frameedge_) {
674
            // Early edge - eventStart points to the beginning of the frame
675
            LOG(DEBUG) << " - track close to start of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
676
677
                       << Units::display(fabs(track->timestamp() - event->start()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
678
            hCutHisto->Fill(4);
679
            num_tracks_++;
680
681
682
            continue;
        }

683
        // Calculate in-pixel position of track in microns
684
        auto inpixel = m_detector->inPixel(localIntercept);
685
686
        auto xmod_um = inpixel.X() * 1000.; // convert mm -> um
        auto ymod_um = inpixel.Y() * 1000.; // convert mm -> um
687

688
        // Loop over all associated DUT clusters:
689
        for(auto assoc_cluster : track->getAssociatedClusters(m_detector->getName())) {
690
            LOG(DEBUG) << " - Looking at next associated DUT cluster";
691

692
            // if closest cluster should be used continue if current associated cluster is not the closest one
693
            if(use_closest_cluster_ && track->getClosestCluster(m_detector->getName()) != assoc_cluster) {
694
                continue;
695
            }
696
            has_associated_cluster = true;
697

698
            hTrackZPosDUT->Fill(track->getState(m_detector->getName()).z());
699
            // Check distance between track and cluster
700
            ROOT::Math::XYZPoint intercept = track->getIntercept(assoc_cluster->global().z());
701
702
            double xdistance = intercept.X() - assoc_cluster->global().x();
            double ydistance = intercept.Y() - assoc_cluster->global().y();
703
            double xdistance_um = xdistance * 1000.;
704
            double ydistance_um = ydistance * 1000.;
705

706
707
708
            double xabsdistance = fabs(xdistance);
            double yabsdistance = fabs(ydistance);
            double tdistance = track->timestamp() - assoc_cluster->timestamp();
709
710
711
            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()));
712
            double posDiff_um = posDiff * 1000.;
713

714
715
716
            hClusterMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row());
            hClusterSizeMapAssoc->Fill(
                assoc_cluster->column(), assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
717
718
            hClusterSizeVsColAssoc->Fill(assoc_cluster->column(), static_cast<double>(assoc_cluster->size()));
            hClusterSizeVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
719
720
            hClusterWidthColVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->columnWidth()));
            hClusterWidthRowVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->rowWidth()));
721
722
723
724

            // 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...
725
            auto cluster_charge = assoc_cluster->charge();
726
            auto normalized_charge = cluster_charge * norm;
727
728

            // clusterChargeAssoc->Fill(normalized_charge);
729
            clusterChargeAssoc->Fill(cluster_charge);
730
            seedChargeAssoc->Fill(assoc_cluster->getSeedPixel()->charge());
731
            hClusterChargeMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row(), cluster_charge);
732
733
734
735
            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());
736
737
            hClusterChargeVsRowAssoc_2D->Fill(assoc_cluster->row(), cluster_charge);
            hSeedChargeVsRowAssoc_2D->Fill(assoc_cluster->row(), assoc_cluster->getSeedPixel()->charge());
738
739

            // Fill per-pixel histograms
740
            for(auto& pixel : assoc_cluster->pixels()) {
741
742
743
744
                hHitMapAssoc->Fill(pixel->column(), pixel->row());
                hPixelRawValueAssoc->Fill(pixel->raw());
                hPixelRawValueMapAssoc->Fill(pixel->column(), pixel->row(), pixel->raw());
            }
745

746
            associatedTracksVersusTime->Fill(track->timestamp() / 1e9); // convert ns -> s
747

748
            // Residuals
749
750
751
752
            residualsX->Fill(xdistance_um);
            residualsY->Fill(ydistance_um);
            residualsPos->Fill(posDiff_um);
            residualsPosVsresidualsTime->Fill(tdistance, posDiff_um);
753

754
            if(assoc_cluster->columnWidth() == 1) {
755
                residualsX1pix->Fill(xdistance_um);
756
757
            }
            if(assoc_cluster->rowWidth() == 1) {
758
                residualsY1pix->Fill(ydistance_um);
759
            }
760
761

            if(assoc_cluster->columnWidth() == 2) {
762
                residualsX2pix->Fill(xdistance_um);
763
764
            }
            if(assoc_cluster->rowWidth() == 2) {
765
                residualsY2pix->Fill(ydistance_um);
766
            }
767

768
            if(assoc_cluster->columnWidth() == 3) {
769
                residualsX3pix->Fill(xdistance_um);
770
771
            }
            if(assoc_cluster->rowWidth() == 3) {
772
                residualsY3pix->Fill(ydistance_um);
773
774
775
            }

            if(assoc_cluster->columnWidth() == 4) {
776
                residualsX4pix->Fill(xdistance_um);
777
778
            }
            if(assoc_cluster->rowWidth() == 4) {
779
                residualsY4pix->Fill(ydistance_um);
780
781
782
            }

            if(assoc_cluster->columnWidth() >= 5) {
783
                residualsXatLeast5pix->Fill(xdistance_um);
784
785
            }
            if(assoc_cluster->rowWidth() >= 5) {
786
                residualsYatLeast5pix->Fill(ydistance_um);
787
788
            }

789
790
            // Time residuals
            residualsTime->Fill(tdistance);
791
            residualsTimeVsTime->Fill(track->timestamp() / 1e9, tdistance); // convert ns -> s
792
            residualsTimeVsTot->Fill(tdistance, assoc_cluster->getSeedPixel()->raw());
793
794
            residualsTimeVsCol->Fill(tdistance, assoc_cluster->getSeedPixel()->column());
            residualsTimeVsRow->Fill(tdistance, assoc_cluster->getSeedPixel()->row());
795
            residualsTimeVsSignal->Fill(tdistance, cluster_charge);
796

797
798
            clusterSizeAssoc->Fill(static_cast<double>(assoc_cluster->size()));
            clusterSizeAssocNorm->Fill(static_cast<double>(assoc_cluster->size()));
Simon Spannagel's avatar
Simon Spannagel committed
799
800
            clusterWidthRowAssoc->Fill(static_cast<double>(assoc_cluster->rowWidth()));
            clusterWidthColAssoc->Fill(static_cast<double>(assoc_cluster->columnWidth()));
801
802

            // Fill in-pixel plots: (all as function of track position within pixel cell)
803
804
            qvsxmym->Fill(xmod_um, ymod_um, cluster_charge);                     // cluster charge profile
            qMoyalvsxmym->Fill(xmod_um, ymod_um, exp(-normalized_charge / 3.5)); // norm. cluster charge profile
805
806

            // mean charge of cluster seed
807
            pxqvsxmym->Fill(xmod_um, ymod_um, assoc_cluster->getSeedPixel()->charge());
808

Jens Kroeger's avatar
Jens Kroeger committed
809
810
811
812
813
            if(assoc_cluster->size() > 1) {
                for(auto& px : assoc_cluster->pixels()) {
                    if(px == assoc_cluster->getSeedPixel()) {
                        continue; // don't fill this histogram for seed pixel!
                    }
814
815
                    pxTimeMinusSeedTime->Fill(static_cast<double>(
                        Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")));
Jens Kroeger's avatar
Jens Kroeger committed
816
                    pxTimeMinusSeedTime_vs_pxCharge->Fill(
817
818
                        static_cast<double>(
                            Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
819
820
821
822
                        px->charge());

                    if(assoc_cluster->size() == 2) {
                        pxTimeMinusSeedTime_vs_pxCharge_2px->Fill(
823
824
                            static_cast<double>(
                                Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
825
826
827
                            px->charge());
                    } else if(assoc_cluster->size() == 3) {
                        pxTimeMinusSeedTime_vs_pxCharge_3px->Fill(
828
829
                            static_cast<double>(
                                Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
830
831
832
                            px->charge());
                    } else if(assoc_cluster->size() == 4) {
                        pxTimeMinusSeedTime_vs_pxCharge_4px->Fill(
833
834
                            static_cast<double>(
                                Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
835
836
837
838
839
                            px->charge());
                    }
                }
            }

840
            // mean cluster size
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
            npxvsxmym->Fill(xmod_um, ymod_um, static_cast<double>(assoc_cluster->size()));
            if(assoc_cluster->size() == 1) {
                npx1vsxmym->Fill(xmod_um, ymod_um);
                qvsxmym_1px->Fill(xmod_um, ymod_um, cluster_charge);
            }
            if(assoc_cluster->size() == 2) {
                npx2vsxmym->Fill(xmod_um, ymod_um);
                qvsxmym_2px->Fill(xmod_um, ymod_um, cluster_charge);
            }
            if(assoc_cluster->size() == 3) {
                npx3vsxmym->Fill(xmod_um, ymod_um);
                qvsxmym_3px->Fill(xmod_um, ymod_um, cluster_charge);
            }
            if(assoc_cluster->size() == 4) {
                npx4vsxmym->Fill(xmod_um, ymod_um);
                qvsxmym_4px->Fill(xmod_um, ymod_um, cluster_charge);
            }
858
859

            // residual MAD x, y, combined (sqrt(x*x + y*y))
860
861
862
            rmsxvsxmym->Fill(xmod_um, ymod_um, xabsdistance);
            rmsyvsxmym->Fill(xmod_um, ymod_um, yabsdistance);
            rmsxyvsxmym->Fill(xmod_um, ymod_um, fabs(sqrt(xdistance * xdistance + ydistance * ydistance)));
863

864
865
            hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
            hAssociatedTracksLocalPosition->Fill(m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept));
866
        }
867
868
        if(!has_associated_cluster) {
            hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
869
        }
870
        num_tracks_++;
871
872
    }
    // Return value telling analysis to keep running
873
    return StatusCode::Success;
874
}
875

876
void AnalysisDUT::finalize(const std::shared_ptr<ReadonlyClipboard>&) {
877
    hCutHisto->Scale(1 / double(num_tracks_));
878
879
    clusterSizeAssocNorm->Scale(1 / clusterSizeAssoc->Integral());
}