AnalysisDUT.cpp 39.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
30
31
32
    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");
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
    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);
85
    hClusterChargeMapAssoc = new TProfile2D("clusterChargeMapAssoc",
86
                                            "Charge map for associated clusters; cluster charge [e]; #entries",
87
                                            m_detector->nPixels().X(),
88
89
                                            -0.5,
                                            m_detector->nPixels().X() - 0.5,
90
                                            m_detector->nPixels().Y(),
91
92
                                            -0.5,
                                            m_detector->nPixels().Y() - 0.5,
93
94
                                            0,
                                            500);
95
96
    hClusterChargeVsColAssoc =
        new TProfile("clusterChargeVsColAssoc",
97
                     "cluster charge vs. column for assoc clusters;cluster column;mean cluster charge [e]",
98
99
100
101
102
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     0,
                     100);
103
    hClusterChargeVsRowAssoc = new TProfile("clusterChargeVsRowAssoc",
Jens Kroeger's avatar
Jens Kroeger committed
104
105
                                            "cluster charge vs. row for assoc clusters;cluster row;mean cluster charge [e]",
                                            m_detector->nPixels().Y(),
106
                                            -0.5,
Jens Kroeger's avatar
Jens Kroeger committed
107
                                            m_detector->nPixels().Y() - 0.5,
108
109
                                            0,
                                            100);
110

111
112
    hSeedChargeVsColAssoc =
        new TProfile("seedChargeVsColAssoc",
113
                     "seed pixel charge vs. column for assoc clusters;cluster column;mean seed pixel charge [e]",
114
115
116
117
118
119
120
121
122
123
124
125
126
                     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);
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    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);
143

Lennart Huth's avatar
Lennart Huth committed
144
    hTrackZPosDUT = new TH1F("globalTrackZPosOnDUT",
Simon Spannagel's avatar
Simon Spannagel committed
145
                             "Global z-position of track on the DUT; global z of track intersection [mm]; #entries ",
146
147
148
                             400,
                             m_detector->displacement().z() - 10,
                             m_detector->displacement().z() + 10);
149
    // Per-pixel histograms
150
    hHitMapAssoc = new TH2F("hitMapAssoc",
151
                            "Hit map of pixels from associated clusters; hit column; hit row",
152
                            m_detector->nPixels().X(),
153
154
                            -0.5,
                            m_detector->nPixels().X() - 0.5,
155
                            m_detector->nPixels().Y(),
156
157
                            -0.5,
                            m_detector->nPixels().Y() - 0.5);
158
159
160
161
162
    hPixelRawValueAssoc = new TH1F("pixelRawValueAssoc",
                                   "Charge distribution of pixels from associated clusters;pixel raw value;#entries",
                                   1024,
                                   -0.5,
                                   1023.5);
163
    hPixelRawValueMapAssoc = new TProfile2D("pixelRawValueMapAssoc",
164
                                            "Charge map of pixels from associated clusters;pixel raw values;# entries",
165
                                            m_detector->nPixels().X(),
166
167
                                            -0.5,
                                            m_detector->nPixels().X() - 0.5,
168
                                            m_detector->nPixels().Y(),
169
170
                                            -0.5,
                                            m_detector->nPixels().Y() - 0.5,
171
172
                                            0,
                                            255);
173

174
    associatedTracksVersusTime =
175
        new TH1F("associatedTracksVersusTime", "Associated tracks over time;time [s];# associated tracks", 300000, 0, 300);
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    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,
        0.,
        0.2);

    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);
219
220
221
222
223
224

    clusterChargeAssoc = new TH1F("clusterChargeAssociated",
                                  "Charge distribution of associated clusters;cluster charge [e];# entries",
                                  10000,
                                  0,
                                  10000);
225
226
227
228
229
    seedChargeAssoc = new TH1F("seedChargeAssociated",
                               "Charge distribution of seed pixels for associated clusters;seed pixel charge [e];# entries",
                               10000,
                               0,
                               10000);
230
231
    clusterSizeAssoc = new TH1F(
        "clusterSizeAssociated", "Size distribution of associated clusters;cluster size; # entries", 30, -0.5, 29.5);
232
233
234
235
236
237
    clusterSizeAssocNorm =
        new TH1F("clusterSizeAssociatedNormalized",
                 "Normalized size distribution of associated clusters;cluster size;# entries (normalized)",
                 30,
                 0,
                 30);
238
239
240
241
242
243
244
245
246
247
    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);
248

249
    // In-pixel studies:
Jin Zhang's avatar
Jin Zhang committed
250
251
    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"));
252
    std::string mod_axes = "in-pixel x_{track} [#mum];in-pixel y_{track} [#mum];";
253

254
    // cut flow histogram
Jin Zhang's avatar
Jin Zhang committed
255
    std::string title = m_detector->getName() + ": number of tracks discarded by different cuts;cut type;tracks";
256
    hCutHisto = new TH1F("hCutHisto", title.c_str(), 4, 1, 5);
257
258
259
260
    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");
261

262
    title = "Resolution in X;" + mod_axes + "MAD(#Deltax) [#mum]";
263
264
265
266
267
268
269
270
    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.);
271

272
    title = "Resolution in Y;" + mod_axes + "MAD(#Deltay) [#mum]";
273
274
275
276
277
278
279
280
    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.);
281

282
    title = "Resolution;" + mod_axes + "MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
283
284
285
286
287
288
289
290
    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.);
291

292
    title = "Mean cluster charge map;" + mod_axes + "<cluster charge> [ke]";
293
294
295
296
297
298
299
300
301
302
    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);
303

304
    title = "Most probable cluster charge map, Moyal approx.;" + mod_axes + "cluster charge MPV [ke]";
305
306
307
308
309
310
311
312
313
314
    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);
315

316
    title = "Seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
317
318
319
320
321
322
323
324
325
326
    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);
327

328
    title = "Mean cluster size map;" + mod_axes + "<pixels/cluster>";
329
330
331
332
333
334
335
336
337
338
    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);
339

340
    title = "1-pixel cluster map;" + mod_axes + "clusters";
341
342
343
344
345
346
347
348
    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.);
349

350
    title = "2-pixel cluster map;" + mod_axes + "clusters";
351
352
353
354
355
356
357
358
    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.);
359

360
    title = "3-pixel cluster map;" + mod_axes + "clusters";
361
362
363
364
365
366
367
368
    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.);
369

370
    title = "4-pixel cluster map;" + mod_axes + "clusters";
371
372
373
374
375
376
377
378
    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.);
379

380
381
    residualsTime = new TH1F("residualsTime",
                             "Time residual;time_{track}-time_{hit} [ns];#entries",
382
                             n_timebins_,
383
384
                             -(n_timebins_ + 1) / 2. * time_binning_,
                             (n_timebins_ - 1) / 2. * time_binning_);
385
386

    hTrackCorrelationX =
387
        new TH1F("hTrackCorrelationX", "Track residual X, all clusters;x_{track}-x_{hit} [mm];# entries", 4000, -10., 10.);
388
    hTrackCorrelationY =
389
390
391
392
393
394
395
396
397
398
399
        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);
400
401
    hTrackCorrelationPosVsCorrelationTime =
        new TH2F("hTrackCorrelationPosVsCorrelationTime",
402
                 "Track time vs. distance residual;time_{track}-time_{hit} [ns];|pos_{track}-pos_{hit}| [mm];# entries",
403
404
405
406
407
408
                 20000,
                 -5000,
                 5000,
                 2100,
                 -1.,
                 10.);
409
410

    residualsTimeVsTime = new TH2F("residualsTimeVsTime",
411
                                   "Time residual vs. time;time [s];time_{track}-time_{hit} [ns];# entries",
412
413
414
                                   20000,
                                   0,
                                   200,
415
                                   n_timebins_,
416
417
                                   -(n_timebins_ + 1) / 2. * time_binning_,
                                   (n_timebins_ - 1) / 2. * time_binning_);
418
419
420
    residualsTimeVsTot =
        new TH2F("residualsTimeVsTot",
                 "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel ToT [lsb];# entries",
421
                 n_timebins_,
422
423
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
424
425
426
                 1000,
                 -0.5,
                 999.5);
427

428
429
430
    residualsTimeVsCol =
        new TH2F("residualsTimeVsCol",
                 "Time residual vs. pixel charge;time_{track} - time_{hit} [ns];seed pixel column;# entries",
431
                 n_timebins_,
432
433
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
434
435
436
437
438
439
                 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",
440
                                  n_timebins_,
441
442
                                  -(n_timebins_ + 1) / 2. * time_binning_,
                                  (n_timebins_ - 1) / 2. * time_binning_,
443
444
445
446
                                  m_detector->nPixels().X(),
                                  -0.5,
                                  m_detector->nPixels().X() - 0.5);

447
448
449
    residualsTimeVsSignal =
        new TH2F("residualsTimeVsSignal",
                 "Time residual vs. cluster charge;cluster charge [e];time_{track}-time_{hit} [mm];# entries",
450
                 n_timebins_,
451
452
                 -(n_timebins_ + 1) / 2. * time_binning_,
                 (n_timebins_ - 1) / 2. * time_binning_,
453
454
455
                 20000,
                 0,
                 100000);
456
457

    hAssociatedTracksGlobalPosition =
458
        new TH2F("hAssociatedTracksGlobalPosition",
459
                 "Map of associated track positions (global);global intercept x [mm];global intercept y [mm]",
460
461
462
                 400,
                 -20,
                 20,
463
464
465
                 200,
                 -10,
                 10);
466
467
468
469
470
471
472
473
474
    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);
475
476
    hUnassociatedTracksGlobalPosition =
        new TH2F("hUnassociatedTracksGlobalPosition",
477
                 "Map of not associated track positions (global); global intercept x [mm]; global intercept y [mm]",
478
479
480
481
482
483
                 200,
                 -10,
                 10,
                 200,
                 -10,
                 10);
Jens Kroeger's avatar
Jens Kroeger committed
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525

    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);
526
527
}

528
StatusCode AnalysisDUT::run(const std::shared_ptr<Clipboard>& clipboard) {
529
530

    // Get the telescope tracks from the clipboard
531
    auto tracks = clipboard->getData<Track>();
532
533

    // Loop over all tracks
534
    for(auto& track : tracks) {
535
536
        // Flags to select clusters and tracks
        bool has_associated_cluster = false;
537
538
539
        LOG(DEBUG) << "Looking at next track";

        // Cut on the chi2/ndof
540
        if(track->getChi2ndof() > chi2_ndof_cut_) {
541
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
542
            hCutHisto->Fill(1);
543
            num_tracks_++;
544
545
546
547
            continue;
        }

        // Check if it intercepts the DUT
548
        auto globalIntercept = m_detector->getIntercept(track.get());
549
        auto localIntercept = m_detector->globalToLocal(globalIntercept);
550

551
        if(!m_detector->hasIntercept(track.get(), 0.5)) {
552
            LOG(DEBUG) << " - track outside DUT area";
553
            hCutHisto->Fill(2);
554
            num_tracks_++;
555
556
557
            continue;
        }

558
        // Check that track is within region of interest using winding number algorithm
559
        if(!m_detector->isWithinROI(track.get())) {
560
            continue;
561
562
563
        }

        // Check that it doesn't go through/near a masked pixel
564
        if(m_detector->hitMasked(track.get(), 1.)) {
565
            LOG(DEBUG) << " - track close to masked pixel";
566
            hCutHisto->Fill(3);
567
            num_tracks_++;
568
569
570
            continue;
        }

Simon Spannagel's avatar
Simon Spannagel committed
571
        // Get the event:
572
        auto event = clipboard->getEvent();
Simon Spannagel's avatar
Simon Spannagel committed
573

574
        // Discard tracks which are very close to the frame edges
575
        if(fabs(track->timestamp() - event->end()) < time_cut_frameedge_) {
576
            // Late edge - eventEnd points to the end of the frame`
577
            LOG(DEBUG) << " - track close to end of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
578
579
                       << Units::display(fabs(track->timestamp() - event->end()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
580
            hCutHisto->Fill(4);
581
            num_tracks_++;
582
            continue;
583
        } else if(fabs(track->timestamp() - event->start()) < time_cut_frameedge_) {
584
            // Early edge - eventStart points to the beginning of the frame
585
            LOG(DEBUG) << " - track close to start of readout frame: "
Simon Spannagel's avatar
Simon Spannagel committed
586
587
                       << Units::display(fabs(track->timestamp() - event->start()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
588
            hCutHisto->Fill(4);
589
            num_tracks_++;
590
591
592
            continue;
        }

593
        // Calculate in-pixel position of track in microns
594
595
596
        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"));
597

598
        // Loop over all associated DUT clusters:
599
        for(auto assoc_cluster : track->getAssociatedClusters(m_detector->getName())) {
600
            LOG(DEBUG) << " - Looking at next associated DUT cluster";
601

602
            // if closest cluster should be used continue if current associated cluster is not the closest one
603
            if(use_closest_cluster_ && track->getClosestCluster(m_detector->getName()) != assoc_cluster) {
604
                continue;
605
            }
606
            has_associated_cluster = true;
607

608
            hTrackZPosDUT->Fill(track->getState(m_detector->getName()).z());
609
            // Check distance between track and cluster
610
            ROOT::Math::XYZPoint intercept = track->getIntercept(assoc_cluster->global().z());
611
612
613
614
615
            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();
616
617
618
619
620
            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
621
622
            hTrackCorrelationX->Fill(static_cast<double>(Units::convert(xdistance, "mm")));
            hTrackCorrelationY->Fill(static_cast<double>(Units::convert(ydistance, "mm")));
623
624
625
            hTrackCorrelationTime->Fill(tdistance);
            hTrackCorrelationPos->Fill(posDiff);
            hTrackCorrelationPosVsCorrelationTime->Fill(track->timestamp() - assoc_cluster->timestamp(), posDiff);
626

627
628
629
            hClusterMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row());
            hClusterSizeMapAssoc->Fill(
                assoc_cluster->column(), assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
630
631
            hClusterSizeVsColAssoc->Fill(assoc_cluster->column(), static_cast<double>(assoc_cluster->size()));
            hClusterSizeVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
632
633
            hClusterWidthColVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->columnWidth()));
            hClusterWidthRowVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->rowWidth()));
634
635
636
637

            // 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...
638
            auto cluster_charge = assoc_cluster->charge();
639
            auto normalized_charge = cluster_charge * norm;
640
641

            // clusterChargeAssoc->Fill(normalized_charge);
642
            clusterChargeAssoc->Fill(cluster_charge);
643
            seedChargeAssoc->Fill(assoc_cluster->getSeedPixel()->charge());
644
            hClusterChargeMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row(), cluster_charge);
645
646
647
648
            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());
649
650
            hClusterChargeVsRowAssoc_2D->Fill(assoc_cluster->row(), cluster_charge);
            hSeedChargeVsRowAssoc_2D->Fill(assoc_cluster->row(), assoc_cluster->getSeedPixel()->charge());
651
652

            // Fill per-pixel histograms
653
            for(auto& pixel : assoc_cluster->pixels()) {
654
655
656
657
                hHitMapAssoc->Fill(pixel->column(), pixel->row());
                hPixelRawValueAssoc->Fill(pixel->raw());
                hPixelRawValueMapAssoc->Fill(pixel->column(), pixel->row(), pixel->raw());
            }
658

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

661
            // Residuals
662
663
664
665
            residualsX->Fill(static_cast<double>(Units::convert(xdistance, "um")));
            residualsY->Fill(static_cast<double>(Units::convert(ydistance, "um")));
            residualsPos->Fill(static_cast<double>(Units::convert(posDiff, "um")));
            residualsPosVsresidualsTime->Fill(tdistance, static_cast<double>(Units::convert(posDiff, "um")));
666

667
            if(assoc_cluster->columnWidth() == 1) {
668
                residualsX1pix->Fill(static_cast<double>(Units::convert(xdistance, "um")));
669
670
            }
            if(assoc_cluster->rowWidth() == 1) {
671
                residualsY1pix->Fill(static_cast<double>(Units::convert(ydistance, "um")));
672
            }
673
674

            if(assoc_cluster->columnWidth() == 2) {
675
                residualsX2pix->Fill(static_cast<double>(Units::convert(xdistance, "um")));
676
677
            }
            if(assoc_cluster->rowWidth() == 2) {
678
                residualsY2pix->Fill(static_cast<double>(Units::convert(ydistance, "um")));
679
            }
680

681
            if(assoc_cluster->columnWidth() == 3) {
682
                residualsX3pix->Fill(static_cast<double>(Units::convert(xdistance, "um")));
683
684
            }
            if(assoc_cluster->rowWidth() == 3) {
685
                residualsY3pix->Fill(static_cast<double>(Units::convert(ydistance, "um")));
686
687
688
            }

            if(assoc_cluster->columnWidth() == 4) {
689
                residualsX4pix->Fill(static_cast<double>(Units::convert(xdistance, "um")));
690
691
            }
            if(assoc_cluster->rowWidth() == 4) {
692
                residualsY4pix->Fill(static_cast<double>(Units::convert(ydistance, "um")));
693
694
695
            }

            if(assoc_cluster->columnWidth() >= 5) {
696
                residualsXatLeast5pix->Fill(static_cast<double>(Units::convert(xdistance, "um")));
697
698
            }
            if(assoc_cluster->rowWidth() >= 5) {
699
                residualsYatLeast5pix->Fill(static_cast<double>(Units::convert(ydistance, "um")));
700
701
            }

702
703
            // Time residuals
            residualsTime->Fill(tdistance);
704
            residualsTimeVsTime->Fill(static_cast<double>(Units::convert(track->timestamp(), "s")), tdistance);
705
            residualsTimeVsTot->Fill(tdistance, assoc_cluster->getSeedPixel()->raw());
706
707
            residualsTimeVsCol->Fill(tdistance, assoc_cluster->getSeedPixel()->column());
            residualsTimeVsRow->Fill(tdistance, assoc_cluster->getSeedPixel()->row());
708
            residualsTimeVsSignal->Fill(tdistance, cluster_charge);
709

710
711
            clusterSizeAssoc->Fill(static_cast<double>(assoc_cluster->size()));
            clusterSizeAssocNorm->Fill(static_cast<double>(assoc_cluster->size()));
Simon Spannagel's avatar
Simon Spannagel committed
712
713
            clusterWidthRowAssoc->Fill(static_cast<double>(assoc_cluster->rowWidth()));
            clusterWidthColAssoc->Fill(static_cast<double>(assoc_cluster->columnWidth()));
714
715

            // Fill in-pixel plots: (all as function of track position within pixel cell)
716
717
718
719
720
721
            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());

Jens Kroeger's avatar
Jens Kroeger committed
722
723
724
725
726
            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!
                    }
727
728
                    pxTimeMinusSeedTime->Fill(static_cast<double>(
                        Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")));
Jens Kroeger's avatar
Jens Kroeger committed
729
                    pxTimeMinusSeedTime_vs_pxCharge->Fill(
730
731
                        static_cast<double>(
                            Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
732
733
734
735
                        px->charge());

                    if(assoc_cluster->size() == 2) {
                        pxTimeMinusSeedTime_vs_pxCharge_2px->Fill(
736
737
                            static_cast<double>(
                                Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
738
739
740
                            px->charge());
                    } else if(assoc_cluster->size() == 3) {
                        pxTimeMinusSeedTime_vs_pxCharge_3px->Fill(
741
742
                            static_cast<double>(
                                Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
743
744
745
                            px->charge());
                    } else if(assoc_cluster->size() == 4) {
                        pxTimeMinusSeedTime_vs_pxCharge_4px->Fill(
746
747
                            static_cast<double>(
                                Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
748
749
750
751
752
                            px->charge());
                    }
                }
            }

753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
            // 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)));

769
770
            hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
            hAssociatedTracksLocalPosition->Fill(m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept));
771
        }
772
773
        if(!has_associated_cluster) {
            hUnassociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
774
        }
775
        num_tracks_++;
776
777
    }
    // Return value telling analysis to keep running
778
    return StatusCode::Success;
779
}
780

781
void AnalysisDUT::finalize(const std::shared_ptr<ReadonlyClipboard>&) {
782
    hCutHisto->Scale(1 / double(num_tracks_));
783
784
    clusterSizeAssocNorm->Scale(1 / clusterSizeAssoc->Integral());
}