AnalysisDUT.cpp 38.7 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
    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,
190
191
        -0.25,
        199.75);
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

    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:
250
251
    auto pitch_x = m_detector->getPitch().X() * 1000.; // convert mm -> um
    auto pitch_y = m_detector->getPitch().Y() * 1000.; // convert mm -> 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
        auto inpixel = m_detector->inPixel(localIntercept);
595
596
        auto xmod = inpixel.X() * 1000.; // convert mm -> um
        auto ymod = inpixel.Y() * 1000.; // convert mm -> 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
            double xdistance = intercept.X() - assoc_cluster->global().x();
            double ydistance = intercept.Y() - assoc_cluster->global().y();
613
            double xdistance_um = xdistance * 1000.;
614
            double ydistance_um = ydistance * 1000.;
615

616
617
618
            double xabsdistance = fabs(xdistance);
            double yabsdistance = fabs(ydistance);
            double tdistance = track->timestamp() - assoc_cluster->timestamp();
619
620
621
            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()));
622
            double posDiff_um = posDiff * 1000.;
623
624

            // Correlation plots
625
626
            hTrackCorrelationX->Fill(xdistance);
            hTrackCorrelationY->Fill(ydistance);
627
628
629
            hTrackCorrelationTime->Fill(tdistance);
            hTrackCorrelationPos->Fill(posDiff);
            hTrackCorrelationPosVsCorrelationTime->Fill(track->timestamp() - assoc_cluster->timestamp(), posDiff);
630

631
632
633
            hClusterMapAssoc->Fill(assoc_cluster->column(), assoc_cluster->row());
            hClusterSizeMapAssoc->Fill(
                assoc_cluster->column(), assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
634
635
            hClusterSizeVsColAssoc->Fill(assoc_cluster->column(), static_cast<double>(assoc_cluster->size()));
            hClusterSizeVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->size()));
636
637
            hClusterWidthColVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->columnWidth()));
            hClusterWidthRowVsRowAssoc->Fill(assoc_cluster->row(), static_cast<double>(assoc_cluster->rowWidth()));
638
639
640
641

            // 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...
642
            auto cluster_charge = assoc_cluster->charge();
643
            auto normalized_charge = cluster_charge * norm;
644
645

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

            // Fill per-pixel histograms
657
            for(auto& pixel : assoc_cluster->pixels()) {
658
659
660
661
                hHitMapAssoc->Fill(pixel->column(), pixel->row());
                hPixelRawValueAssoc->Fill(pixel->raw());
                hPixelRawValueMapAssoc->Fill(pixel->column(), pixel->row(), pixel->raw());
            }
662

663
            associatedTracksVersusTime->Fill(track->timestamp() / 1e9); // convert ns -> s
664

665
            // Residuals
666
667
668
669
            residualsX->Fill(xdistance_um);
            residualsY->Fill(ydistance_um);
            residualsPos->Fill(posDiff_um);
            residualsPosVsresidualsTime->Fill(tdistance, posDiff_um);
670

671
            if(assoc_cluster->columnWidth() == 1) {
672
                residualsX1pix->Fill(xdistance_um);
673
674
            }
            if(assoc_cluster->rowWidth() == 1) {
675
                residualsY1pix->Fill(ydistance_um);
676
            }
677
678

            if(assoc_cluster->columnWidth() == 2) {
679
                residualsX2pix->Fill(xdistance_um);
680
681
            }
            if(assoc_cluster->rowWidth() == 2) {
682
                residualsY2pix->Fill(ydistance_um);
683
            }
684

685
            if(assoc_cluster->columnWidth() == 3) {
686
                residualsX3pix->Fill(xdistance_um);
687
688
            }
            if(assoc_cluster->rowWidth() == 3) {
689
                residualsY3pix->Fill(ydistance_um);
690
691
692
            }

            if(assoc_cluster->columnWidth() == 4) {
693
                residualsX4pix->Fill(xdistance_um);
694
695
            }
            if(assoc_cluster->rowWidth() == 4) {
696
                residualsY4pix->Fill(ydistance_um);
697
698
699
            }

            if(assoc_cluster->columnWidth() >= 5) {
700
                residualsXatLeast5pix->Fill(xdistance_um);
701
702
            }
            if(assoc_cluster->rowWidth() >= 5) {
703
                residualsYatLeast5pix->Fill(ydistance_um);
704
705
            }

706
707
            // Time residuals
            residualsTime->Fill(tdistance);
708
            residualsTimeVsTime->Fill(track->timestamp() / 1e9, tdistance); // convert ns -> s
709
            residualsTimeVsTot->Fill(tdistance, assoc_cluster->getSeedPixel()->raw());
710
711
            residualsTimeVsCol->Fill(tdistance, assoc_cluster->getSeedPixel()->column());
            residualsTimeVsRow->Fill(tdistance, assoc_cluster->getSeedPixel()->row());
712
            residualsTimeVsSignal->Fill(tdistance, cluster_charge);
713

714
715
            clusterSizeAssoc->Fill(static_cast<double>(assoc_cluster->size()));
            clusterSizeAssocNorm->Fill(static_cast<double>(assoc_cluster->size()));
Simon Spannagel's avatar
Simon Spannagel committed
716
717
            clusterWidthRowAssoc->Fill(static_cast<double>(assoc_cluster->rowWidth()));
            clusterWidthColAssoc->Fill(static_cast<double>(assoc_cluster->columnWidth()));
718
719

            // Fill in-pixel plots: (all as function of track position within pixel cell)
720
721
722
723
724
725
            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
726
727
728
729
730
            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!
                    }
731
732
                    pxTimeMinusSeedTime->Fill(static_cast<double>(
                        Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")));
Jens Kroeger's avatar
Jens Kroeger committed
733
                    pxTimeMinusSeedTime_vs_pxCharge->Fill(
734
735
                        static_cast<double>(
                            Units::convert(px->timestamp() - assoc_cluster->getSeedPixel()->timestamp(), "ns")),
Jens Kroeger's avatar
Jens Kroeger committed
736
737
738
739
                        px->charge());

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

757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
            // 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)));

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

785
void AnalysisDUT::finalize(const std::shared_ptr<ReadonlyClipboard>&) {
786
    hCutHisto->Scale(1 / double(num_tracks_));
787
788
    clusterSizeAssocNorm->Scale(1 / clusterSizeAssoc->Integral());
}