AnalysisTimingATLASpix.cpp 50.8 KB
Newer Older
1
2
3
/**
 * @file
 * @brief Implementation of [AnalysisEfficiency] module
4
5
 *
 * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors.
6
7
8
9
10
 * 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 "AnalysisTimingATLASpix.h"
12
13
14
15
16
17
18
19
20
21

#include "objects/Cluster.hpp"
#include "objects/Pixel.hpp"
#include "objects/Track.hpp"

#include "TF1.h"
#include "TFile.h"

using namespace corryvreckan;

22
AnalysisTimingATLASpix::AnalysisTimingATLASpix(Configuration config, std::shared_ptr<Detector> detector)
23
24
    : Module(std::move(config), detector) {

25
26
27
    // Backwards compatibilty: also allow timing_cut to be used for time_cut_abs
    m_config.setAlias("time_cut_abs", "timing_cut", true);

28
29
    using namespace ROOT::Math;
    m_detector = detector;
30
31
32
33
34
35
36
37
    if(config.count({"time_cut_rel", "time_cut_abs"}) > 1) {
        throw InvalidCombinationError(
            m_config, {"time_cut_rel", "time_cut_abs"}, "Absolute and relative time cuts are mutually exclusive.");
    } else if(m_config.has("time_cut_abs")) {
        m_timeCut = m_config.get<double>("time_cut_abs");
    } else {
        m_timeCut = m_config.get<double>("time_cut_rel", 3.0) * m_detector->getTimeResolution();
    }
38
    m_chi2ndofCut = m_config.get<double>("chi2ndof_cut", 3.);
39
    m_timeCutFrameEdge = m_config.get<double>("time_cut_frameedge", static_cast<double>(Units::convert(20, "ns")));
40
41
42
43
44
45
46

    if(m_config.has("cluster_charge_cut")) {
        m_clusterChargeCut = m_config.get<double>("cluster_charge_cut");
    }
    if(m_config.has("cluster_size_cut")) {
        m_clusterSizeCut = m_config.get<size_t>("cluster_size_cut");
    }
47
    m_highTotCut = m_config.get<int>("high_tot_cut", 40);
48
    m_lowTotCut = m_config.get<int>("low_tot_cut", 40);
49
    m_timingTailCut = m_config.get<double>("timing_tail_cut", static_cast<double>(Units::convert(20, "ns")));
50
51
52

    if(m_config.has("correction_file_row")) {
        m_correctionFile_row = m_config.get<std::string>("correction_file_row");
53
        m_correctionGraph_row = m_config.get<std::string>("correction_graph_row");
54
55
56
57
58
59
        m_pointwise_correction_row = true;
    } else {
        m_pointwise_correction_row = false;
    }
    if(m_config.has("correction_file_timewalk")) {
        m_correctionFile_timewalk = m_config.get<std::string>("correction_file_timewalk");
60
        m_correctionGraph_timewalk = m_config.get<std::string>("correction_graph_timewalk");
61
62
63
64
65
66
67
68
        m_pointwise_correction_timewalk = true;
    } else {
        m_pointwise_correction_timewalk = false;
    }

    m_calcCorrections = m_config.get<bool>("calc_corrections", false);
    m_totBinExample = m_config.get<int>("tot_bin_example", 3);

69
70
71
    m_inpixelBinSize = m_config.get<XYVector>(
        "inpixel_bin_size",
        {static_cast<double>(Units::convert(1.0, "um")), static_cast<double>(Units::convert(1.0, "um"))});
72

73
74
75
76
77
78
79
    total_tracks_uncut = 0;
    tracks_afterChi2Cut = 0;
    tracks_hasIntercept = 0;
    tracks_isWithinROI = 0;
    tracks_afterMasking = 0;
    total_tracks = 0;
    matched_tracks = 0;
80
    tracks_afterClusterChargeCut = 0;
81
82
83
    tracks_afterClusterSizeCut = 0;
}

84
void AnalysisTimingATLASpix::initialise() {
85
86
87
88
89
90

    auto pitch_x = static_cast<double>(Units::convert(m_detector->pitch().X(), "um"));
    auto pitch_y = static_cast<double>(Units::convert(m_detector->pitch().Y(), "um"));

    std::string name = "hTrackCorrelationTime";
    hTrackCorrelationTime =
91
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
92
    hTrackCorrelationTime->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
93
    hTrackCorrelationTime->GetYaxis()->SetTitle("# events");
94
95
96

    name = "hTrackCorrelationTimeAssoc";
    hTrackCorrelationTimeAssoc =
97
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
98
    hTrackCorrelationTimeAssoc->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
99
    hTrackCorrelationTimeAssoc->GetYaxis()->SetTitle("# events");
100

101
102
    name = "hTrackCorrelationTimeAssocVsTime";
    hTrackCorrelationTimeAssocVsTime = new TH2F(name.c_str(), name.c_str(), 3e3, 0, 3e3, 1e3, -5, 5);
103
    hTrackCorrelationTimeAssocVsTime->GetYaxis()->SetTitle("ts_{track} - ts_{cluster} [us]");
104
    hTrackCorrelationTimeAssocVsTime->GetXaxis()->SetTitle("time [s]");
105
    hTrackCorrelationTimeAssocVsTime->GetZaxis()->SetTitle("# events");
106

107
108
109
110
111
112
113
114
    if(m_pointwise_correction_row) {
        name = "hTrackCorrelationTime_rowCorr";
        std::string title = "hTrackCorrelationTime_rowCorr: row-by-row correction";
        hTrackCorrelationTime_rowCorr =
            new TH1F(name.c_str(), title.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
        hTrackCorrelationTime_rowCorr->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
        hTrackCorrelationTime_rowCorr->GetYaxis()->SetTitle("# events");
    }
115

116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
    if(m_pointwise_correction_timewalk) {
        name = "hTrackCorrelationTime_rowAndTWCorr";
        hTrackCorrelationTime_rowAndTWCorr =
            new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
        hTrackCorrelationTime_rowAndTWCorr->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
        hTrackCorrelationTime_rowAndTWCorr->GetYaxis()->SetTitle("# events");

        name = "hTrackCorrelationTime_rowAndTWCorr_l25";
        hTrackCorrelationTime_rowAndTWCorr_l25 =
            new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
        hTrackCorrelationTime_rowAndTWCorr_l25->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns] (if seed tot < 25lsb)");
        hTrackCorrelationTime_rowAndTWCorr_l25->GetYaxis()->SetTitle("# events");

        name = "hTrackCorrelationTime_rowAndTWCorr_l40";
        hTrackCorrelationTime_rowAndTWCorr_l40 =
            new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
        hTrackCorrelationTime_rowAndTWCorr_l40->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns] (if seed tot < 40lsb)");
        hTrackCorrelationTime_rowAndTWCorr_l40->GetYaxis()->SetTitle("# events");

        name = "hTrackCorrelationTime_rowAndTWCorr_g40";
        hTrackCorrelationTime_rowAndTWCorr_g40 =
            new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
        hTrackCorrelationTime_rowAndTWCorr_g40->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns] (if seed tot > 40lsb)");
        hTrackCorrelationTime_rowAndTWCorr_g40->GetYaxis()->SetTitle("# events");

        name = "hTrackCorrelationTime_totBin_" + std::to_string(m_totBinExample);
        hTrackCorrelationTime_example = new TH1D(name.c_str(), name.c_str(), 20000, -5000, 5000);
        hTrackCorrelationTime_example->GetXaxis()->SetTitle(
            "track time stamp - pixel time stamp [ns] (all pixels from cluster)");
        hTrackCorrelationTime_example->GetYaxis()->SetTitle("# events");
    }
147
148
149
150

    // 2D histograms:
    // column dependence
    name = "hTrackCorrelationTimeVsCol";
151
152
    hTrackCorrelationTimeVsCol = new TH2F(
        name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().X(), -0.5, m_detector->nPixels().X() - 0.5);
153
    hTrackCorrelationTimeVsCol->GetYaxis()->SetTitle("pixel column");
154
    hTrackCorrelationTimeVsCol->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns]");
155
156
    // row dependence
    name = "hTrackCorrelationTimeVsRow";
157
158
    hTrackCorrelationTimeVsRow = new TH2F(
        name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), -0.5, m_detector->nPixels().Y() - 0.5);
159
    hTrackCorrelationTimeVsRow->GetYaxis()->SetTitle("pixel row");
160
    hTrackCorrelationTimeVsRow->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns]");
161
    name = "hTrackCorrelationTimeVsRow_1px";
162
163
    hTrackCorrelationTimeVsRow_1px = new TH2F(
        name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), -0.5, m_detector->nPixels().Y() - 0.5);
164
    hTrackCorrelationTimeVsRow_1px->GetYaxis()->SetTitle("pixel row");
165
    hTrackCorrelationTimeVsRow_1px->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns] (single-pixel clusters)");
166
    name = "hTrackCorrelationTimeVsRow_npx";
167
168
    hTrackCorrelationTimeVsRow_npx = new TH2F(
        name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), -0.5, m_detector->nPixels().Y() - 0.5);
169
    hTrackCorrelationTimeVsRow_npx->GetYaxis()->SetTitle("pixel row");
170
    hTrackCorrelationTimeVsRow_npx->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns] (multi-pixel clusters)");
171
172

    // control plot: row dependence after row correction
173
174
175
176
177
178
179
180
181
182
183
184
185
    if(m_pointwise_correction_row) {
        name = "hTrackCorrelationTimeVsRow_rowCorr";
        hTrackCorrelationTimeVsRow_rowCorr = new TH2F(name.c_str(),
                                                      name.c_str(),
                                                      20000,
                                                      -5000,
                                                      5000,
                                                      m_detector->nPixels().Y(),
                                                      -0.5,
                                                      m_detector->nPixels().Y() - 0.5);
        hTrackCorrelationTimeVsRow_rowCorr->GetYaxis()->SetTitle("pixel row");
        hTrackCorrelationTimeVsRow_rowCorr->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns]");
    }
186
187
188

    // control plot: time walk dependence, not row corrected
    name = "hTrackCorrelationTimeVsTot";
189
    hTrackCorrelationTimeVsTot = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
190
    hTrackCorrelationTimeVsTot->GetYaxis()->SetTitle("seed pixel tot [lsb]");
191
    hTrackCorrelationTimeVsTot->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
192
193

    name = "hTrackCorrelationTimeVsTot_1px";
194
    hTrackCorrelationTimeVsTot_1px = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
195
    hTrackCorrelationTimeVsTot_1px->GetYaxis()->SetTitle("seed pixel tot [lsb] (if clustersize = 1)");
196
    hTrackCorrelationTimeVsTot_1px->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
197
198

    name = "hTrackCorrelationTimeVsTot_npx";
199
    hTrackCorrelationTimeVsTot_npx = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
200
    hTrackCorrelationTimeVsTot_npx->GetYaxis()->SetTitle("seed pixel tot [lsb] (if clustersize > 1)");
201
    hTrackCorrelationTimeVsTot_npx->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
202

203
    name = "hTrackCorrelationTimeVsTot_px";
204
    hTrackCorrelationTimeVsTot_px = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
205
206
207
    hTrackCorrelationTimeVsTot_px->GetYaxis()->SetTitle("pixel tot [lsb]");
    hTrackCorrelationTimeVsTot_px->GetXaxis()->SetTitle("ts_{track} - ts_{pixel} (all pixels from cluster) [ns]");

208
209
210
    name = "hClusterTimeMinusPixelTime";
    hClusterTimeMinusPixelTime = new TH1F(name.c_str(), name.c_str(), 2000, -1000, 1000);
    hClusterTimeMinusPixelTime->GetXaxis()->SetTitle(
211
        "ts_{cluster} - ts_{pixel} [ns] (all pixels from cluster (if clusterSize>1))");
212
213

    // timewalk after row correction
214
215
    if(m_pointwise_correction_row) {
        name = "hTrackCorrelationTimeVsTot_rowCorr";
216
        hTrackCorrelationTimeVsTot_rowCorr = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
217
        hTrackCorrelationTimeVsTot_rowCorr->GetYaxis()->SetTitle("seed pixel tot [lsb]");
218
219
220
        hTrackCorrelationTimeVsTot_rowCorr->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns]");

        name = "hTrackCorrelationTimeVsTot_rowCorr_1px";
221
        hTrackCorrelationTimeVsTot_rowCorr_1px = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
222
        hTrackCorrelationTimeVsTot_rowCorr_1px->GetYaxis()->SetTitle("seed pixel tot [lsb] (single-pixel clusters)");
223
224
225
        hTrackCorrelationTimeVsTot_rowCorr_1px->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");

        name = "hTrackCorrelationTimeVsTot_rowCorr_npx";
226
        hTrackCorrelationTimeVsTot_rowCorr_npx = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
227
        hTrackCorrelationTimeVsTot_rowCorr_npx->GetYaxis()->SetTitle("seed pixel tot [lsb] (multi-pixel clusters)");
228
229
        hTrackCorrelationTimeVsTot_rowCorr_npx->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
    }
230
231

    // final plots with both row and timewalk correction:
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    if(m_pointwise_correction_timewalk) {
        name = "hTrackCorrelationTimeVsRow_rowAndTWCorr";
        hTrackCorrelationTimeVsRow_rowAndTWCorr = new TH2F(name.c_str(),
                                                           name.c_str(),
                                                           20000,
                                                           -5000,
                                                           5000,
                                                           m_detector->nPixels().Y(),
                                                           -0.5,
                                                           m_detector->nPixels().Y() - 0.5);
        hTrackCorrelationTimeVsRow_rowAndTWCorr->GetYaxis()->SetTitle("row");
        hTrackCorrelationTimeVsRow_rowAndTWCorr->GetXaxis()->SetTitle("ts_{track} - ts_{seed pixel} [ns]");

        name = "hTrackCorrelationTimeVsTot_rowAndTWCorr";
246
        hTrackCorrelationTimeVsTot_rowAndTWCorr = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 64, 0, 64);
247
        hTrackCorrelationTimeVsTot_rowAndTWCorr->GetYaxis()->SetTitle("seed pixel tot [lsb]");
248
249
        hTrackCorrelationTimeVsTot_rowAndTWCorr->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
    }
250

251
    // in-pixel time resolution plots:
252
253
    auto nbins_x = static_cast<int>(std::ceil(m_detector->pitch().X() / m_inpixelBinSize.X()));
    auto nbins_y = static_cast<int>(std::ceil(m_detector->pitch().Y() / m_inpixelBinSize.Y()));
254
255
256
257
    if(nbins_x > 1e4 || nbins_y > 1e4) {
        throw InvalidValueError(m_config, "inpixel_bin_size", "Too many bins for in-pixel histograms.");
    }

258
259
    std::string title =
        "in-pixel time resolution map;in-pixel x_{track} [#mum];in-pixel y_{track} [#mum];ts_{track} - ts_{cluster} [ns]";
260
261
262
263
264
265
266
267
268
269
270
    hPixelTrackCorrelationTimeMap = new TProfile2D("pixelTrackCorrelationTimeMap",
                                                   title.c_str(),
                                                   nbins_x,
                                                   -pitch_x / 2.,
                                                   pitch_x / 2.,
                                                   nbins_y,
                                                   -pitch_y / 2.,
                                                   pitch_y / 2.,
                                                   -1 * m_timeCut,
                                                   m_timeCut);

271
    name = "hClusterSizeVsTot_Assoc";
272
273
    hClusterSizeVsTot_Assoc = new TH2F(name.c_str(), name.c_str(), 20, 0, 20, 64, 0, 64);
    hClusterSizeVsTot_Assoc->GetYaxis()->SetTitle("pixel ToT [lsb] (all pixels from cluster)");
274
275
276
    hClusterSizeVsTot_Assoc->GetXaxis()->SetTitle("clusterSize");

    hHitMapAssoc = new TH2F("hitMapAssoc",
277
                            "hitMapAssoc; x_{track} [px]; y_{track} [px];# entries",
278
                            m_detector->nPixels().X(),
279
280
                            -0.5,
                            m_detector->nPixels().X() - 0.5,
281
                            m_detector->nPixels().Y(),
282
283
284
285
286
287
288
289
290
291
                            0,
                            m_detector->nPixels().Y());
    hHitMapAssoc_highToT = new TH2F("hitMapAssoc_highToT",
                                    "hitMapAssoc_highToT; x_{track} [px]; y_{track} [px];# entries",
                                    m_detector->nPixels().X(),
                                    -0.5,
                                    m_detector->nPixels().X() - 0.5,
                                    m_detector->nPixels().Y(),
                                    -0.5,
                                    m_detector->nPixels().Y() - 0.5);
292
    hHitMapAssoc_inPixel = new TH2F("hitMapAssoc_inPixel",
293
                                    "hitMapAssoc_inPixel; in-pixel x_{track} [#mum]; in-pixel y_{track} [#mum]",
294
                                    static_cast<int>(pitch_x),
295
296
                                    -pitch_x / 2.,
                                    pitch_x / 2.,
297
                                    static_cast<int>(pitch_y),
298
299
                                    -pitch_y / 2.,
                                    pitch_y / 2.);
300
301
302
303
304
305
306
307
308
309
310
    if(m_config.has("high_tot_cut")) {
        hHitMapAssoc_inPixel_highToT =
            new TH2F("hitMapAssoc_inPixel_highToT",
                     "hitMapAssoc_inPixel_highToT;in-pixel x_{track} [#mum];in-pixel y_{track} [#mum]",
                     static_cast<int>(pitch_x),
                     -pitch_x / 2.,
                     pitch_x / 2.,
                     static_cast<int>(pitch_y),
                     -pitch_y / 2.,
                     pitch_y / 2.);
    }
311
    hClusterMapAssoc = new TH2F("hClusterMapAssoc",
312
                                "hClusterMapAssoc; x_{cluster} [px]; y_{cluster} [px];# entries",
313
                                m_detector->nPixels().X(),
314
315
                                -0.5,
                                m_detector->nPixels().X() - 0.5,
316
                                m_detector->nPixels().Y(),
317
318
                                -0.5,
                                m_detector->nPixels().Y() - 0.5);
319

320
321
322
323
324
325
326
327
328
    hTotVsRow = new TH2F("hTotVsRow",
                         "hTotVsRow;seed-pixel row;seed-pixel ToT [lsb]",
                         m_detector->nPixels().Y(),
                         0,
                         m_detector->nPixels().Y(),
                         64,
                         0,
                         64);

329
330
331
    hTotVsTime = new TH2F("hTotVsTime", "hTotVsTime", 64, 0, 64, 1e6, 0, 100);
    hTotVsTime->GetXaxis()->SetTitle("pixel ToT [lsb]");
    hTotVsTime->GetYaxis()->SetTitle("time [s]");
332
    if(m_config.has("high_tot_cut")) {
333
334
335
        hTotVsTime_highToT = new TH2F("hTotVsTime_highToT", "hTotVsTime_highToT", 64, 0, 64, 1e6, 0, 100);
        hTotVsTime_highToT->GetXaxis()->SetTitle("pixel ToT [lsb] if > high_tot_cut");
        hTotVsTime_highToT->GetYaxis()->SetTitle("time [s]");
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
    }

    // control plots for "left/right tail" and "main peak" of the track time correlation
    if(m_config.has("timing_tail_cut") && m_pointwise_correction_timewalk) {
        hInPixelMap_leftTail = new TH2F("hPixelMap_leftTail",
                                        "in-pixel track position (left tail of time residual);in-pixel x_{track} "
                                        "[#mum];in-pixel y_{track} [#mum];# entries",
                                        nbins_x,
                                        -pitch_x / 2.,
                                        pitch_x / 2.,
                                        nbins_y,
                                        -pitch_y / 2.,
                                        pitch_y / 2.);
        hInPixelMap_rightTail = new TH2F("hPixelMap_rightTail",
                                         "in-pixel track position (right tail of time residual);in-pixel x_{track} "
                                         "[#mum];in-pixel y_{track} [#mum];# entries",
                                         nbins_x,
                                         -pitch_x / 2.,
                                         pitch_x / 2.,
                                         nbins_y,
                                         -pitch_y / 2.,
                                         pitch_y / 2.);
        hInPixelMap_mainPeak = new TH2F("hPixelMap_mainPeak",
                                        "in-pixel track position (main peak of time residual);in-pixel x_{track} "
                                        "[#mum];in-pixel y_{track} [#mum];# entries",
                                        nbins_x,
                                        -pitch_x / 2.,
                                        pitch_x / 2.,
                                        nbins_y,
                                        -pitch_y / 2.,
                                        pitch_y / 2.);
        hClusterMap_leftTail =
            new TH2F("hClusterMap_leftTail",
                     "hClusterMap (left tail of time residual); x_{cluster} [px]; x_{cluster} [px];# entries",
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     m_detector->nPixels().Y(),
                     -0.5,
                     m_detector->nPixels().Y() - 0.5);
        hClusterMap_rightTail =
            new TH2F("hClusterMap_rightTail",
                     "hClusterMap (right tail of time residual); x_{cluster} [px]; x_{cluster} [px];# entries",
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     m_detector->nPixels().Y(),
                     -0.5,
                     m_detector->nPixels().Y() - 0.5);
        hClusterMap_mainPeak =
            new TH2F("hClusterMap_mainPeak",
                     "hClusterMap (main peak of time residual); x_{cluster} [px]; x_{cluster} [px];# entries",
                     m_detector->nPixels().X(),
                     -0.5,
                     m_detector->nPixels().X() - 0.5,
                     m_detector->nPixels().Y(),
                     -0.5,
                     m_detector->nPixels().Y() - 0.5);
        hClusterSize_leftTail =
            new TH1F("clusterSize_leftTail", "clusterSize (left tail of time residual);cluster size;# entries", 100, 0, 100);
        hClusterSize_rightTail = new TH1F(
            "clusterSize_rightTail", "clusterSize (right tail of time residual);cluster size;# entries", 100, 0, 100);
        hClusterSize_mainPeak =
            new TH1F("clusterSize_mainPeak", "clusterSize (main peak of time residual);cluster size;# entries", 100, 0, 100);
        hTot_leftTail =
            new TH1F("hTot_leftTail", "ToT (left tail of time residual);seed pixel ToT [lsb];# events", 2 * 64, -64, 64);
        hTot_rightTail =
            new TH1F("hTot_rightTail", "ToT (right tail of time residual);seed pixel ToT [lsb];# events", 2 * 64, -64, 64);
        hTot_mainPeak = new TH1F("hTot_mainPeak",
                                 "ToT (main peak of time residual, 1px clusters);seed pixel ToT [lsb];# events",
406
407
408
                                 2 * 64,
                                 -64,
                                 64);
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
        hTot_leftTail_1px = new TH1F("hTot_leftTail_1px",
                                     "ToT (left tail of time residual, 1px clusters);seed pixel ToT [lsb];# events",
                                     2 * 64,
                                     -64,
                                     64);
        hTot_rightTail_1px = new TH1F("hTot_rightTail_1px",
                                      "ToT (right tail of time residual, 1px clusters);seed pixel ToT [lsb];# events",
                                      2 * 64,
                                      -64,
                                      64);
        hTot_mainPeak_1px = new TH1F("hTot_mainPeak_1px",
                                     "ToT (main peak of time residual, 1px clusters);seed pixel ToT [lsb];# events",
                                     2 * 64,
                                     -64,
                                     64);
        hPixelTimestamp_leftTail = new TH1F("pixelTimestamp_leftTail",
                                            "pixelTimestamp (left tail of time residual);pixel timestamp [ms];# events",
                                            3e6,
                                            0,
                                            3e3);
        hPixelTimestamp_rightTail = new TH1F("pixelTimestamp_leftTail",
                                             "pixelTimestamp (left tail of time residual);pixel timestamp [ms];# events",
                                             3e6,
                                             0,
                                             3e3);
        hPixelTimestamp_mainPeak = new TH1F("pixelTimestamp_mainPeak",
                                            "pixelTimestamp (left tail of time residual);pixel timestamp [ms];# events",
                                            3e6,
                                            0,
                                            3e3);
    }
440
441
442
443
444

    // /////////////////////////////////////////// //
    // TGraphErrors for Timewalk & Row Correction: //
    // /////////////////////////////////////////// //

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
470
471
472
473
474
    if(m_calcCorrections) {
        gTimeCorrelationVsRow = new TGraphErrors();
        gTimeCorrelationVsRow->SetName("gTimeCorrelationVsRow");
        gTimeCorrelationVsRow->SetTitle("gTimeCorrelationVsRow");
        gTimeCorrelationVsRow->GetXaxis()->SetTitle("row");
        gTimeCorrelationVsRow->GetYaxis()->SetTitle("time correlation peak [ns]");

        int nBinsToT = hTrackCorrelationTimeVsTot_rowCorr->GetNbinsY();
        gTimeCorrelationVsTot_rowCorr = new TGraphErrors(nBinsToT);
        gTimeCorrelationVsTot_rowCorr->SetName("gTimeCorrelationVsTot_rowCorr");
        gTimeCorrelationVsTot_rowCorr->SetTitle("gTimeCorrelationVsTot_rowCorr");
        gTimeCorrelationVsTot_rowCorr->GetXaxis()->SetTitle("pixel ToT [ns]");
        gTimeCorrelationVsTot_rowCorr->GetYaxis()->SetTitle("time correlation peak [ns]");

        nBinsToT = hTrackCorrelationTimeVsTot_rowCorr_1px->GetNbinsY();
        gTimeCorrelationVsTot_rowCorr_1px = new TGraphErrors(nBinsToT);
        gTimeCorrelationVsTot_rowCorr_1px->SetName("gTimeCorrelationVsTot_rowCorr_1px");
        gTimeCorrelationVsTot_rowCorr_1px->SetTitle("gTimeCorrelationVsTot_rowCorr_1px");
        gTimeCorrelationVsTot_rowCorr_1px->GetXaxis()->SetTitle("pixel ToT [ns] (single-pixel clusters)");
        gTimeCorrelationVsTot_rowCorr_1px->GetYaxis()->SetTitle("time correlation peak [ns]");

        nBinsToT = hTrackCorrelationTimeVsTot_rowCorr_npx->GetNbinsY();
        gTimeCorrelationVsTot_rowCorr_npx = new TGraphErrors(nBinsToT);
        gTimeCorrelationVsTot_rowCorr_npx->SetName("gTimeCorrelationVsTot_rowCorr_npx");
        gTimeCorrelationVsTot_rowCorr_npx->SetTitle("gTimeCorrelationVsTot_rowCorr_npx");
        gTimeCorrelationVsTot_rowCorr_npx->GetXaxis()->SetTitle("pixel ToT [ns] (multi-pixel clusters");
        gTimeCorrelationVsTot_rowCorr_npx->GetYaxis()->SetTitle("time correlation peak [ns]");
    }

    LOG(INFO) << "Calculate corrections: = " << m_calcCorrections;
475
476
477

    if(m_pointwise_correction_row) {
        // Import TGraphErrors for row corection:
478
479
480
481
482
483
        TFile file(m_correctionFile_row.c_str());
        if(!file.IsOpen()) {
            throw InvalidValueError(m_config,
                                    "correction_file_row",
                                    "ROOT file doesn't exist. If no row correction shall be applied, remove this parameter "
                                    "from the configuration file.");
484
        }
485
486

        gRowCorr = static_cast<TGraphErrors*>(file.Get(m_correctionGraph_row.c_str()));
487
        // Check if graph exists in ROOT file:
488
        if(!gRowCorr) {
489
490
491
            throw InvalidValueError(
                m_config, "correction_graph_row", "Graph doesn't exist in ROOT file. Use full/path/to/graph.");
        }
492
493
494
495
496
    } else {
        LOG(STATUS) << "----> NO POINTWISE ROW CORRECTION!!!";
    }
    if(m_pointwise_correction_timewalk) {
        // Import TGraphErrors for timewalk corection:
497
498
499
500
501
502
503
        TFile file(m_correctionFile_timewalk.c_str());
        if(!file.IsOpen()) {
            throw InvalidValueError(m_config,
                                    "correction_file_timewalk",
                                    "ROOT file doesn't exist. If no row correction shall be applied, remove this parameter "
                                    "from the configuration file.");
        }
504
505

        gTimeWalkCorr = static_cast<TGraphErrors*>(file.Get(m_correctionGraph_timewalk.c_str()));
506
        // Check if graph exists in ROOT file:
507
        if(!gTimeWalkCorr) {
508
509
            throw InvalidValueError(
                m_config, "correction_graph_timewalk", "Graph doesn't exist in ROOT file. Use full/path/to/graph.");
510
511
512
513
514
515
        }
    } else {
        LOG(STATUS) << "----> NO POINTWISE TIMEWALK CORRECTION!!!";
    }
}

516
StatusCode AnalysisTimingATLASpix::run(std::shared_ptr<Clipboard> clipboard) {
517
518

    // Get the telescope tracks from the clipboard
519
    auto tracks = clipboard->getData<Track>();
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
    if(tracks == nullptr) {
        LOG(DEBUG) << "No tracks on the clipboard";
        return StatusCode::Success;
    }

    // Loop over all tracks
    for(auto& track : (*tracks)) {
        bool has_associated_cluster = false;
        bool is_within_roi = true;
        LOG(DEBUG) << "Looking at next track";
        total_tracks_uncut++;

        // Cut on the chi2/ndof
        if(track->chi2ndof() > m_chi2ndofCut) {
            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
            continue;
        }
        tracks_afterChi2Cut++;

        // Check if it intercepts the DUT
        if(!m_detector->hasIntercept(track, 1.)) {
            LOG(DEBUG) << " - track outside DUT area";
            continue;
        }
        tracks_hasIntercept++;

        // Check that track is within region of interest using winding number algorithm
        if(!m_detector->isWithinROI(track)) {
            LOG(DEBUG) << " - track outside ROI";
            is_within_roi = false;
        }
        tracks_isWithinROI++;

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

        // Get the event:
561
        auto event = clipboard->getEvent();
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581

        // Discard tracks which are very close to the frame edges
        if(fabs(track->timestamp() - event->end()) < m_timeCutFrameEdge) {
            // Late edge - eventEnd points to the end of the frame`
            LOG(DEBUG) << " - track close to end of readout frame: "
                       << Units::display(fabs(track->timestamp() - event->end()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
            continue;
        } else if(fabs(track->timestamp() - event->start()) < m_timeCutFrameEdge) {
            // Early edge - eventStart points to the start of the frame
            LOG(DEBUG) << " - track close to start of readout frame: "
                       << Units::display(fabs(track->timestamp() - event->start()), {"us", "ns"}) << " at "
                       << Units::display(track->timestamp(), {"us"});
            continue;
        }

        // Count this as reference track:
        total_tracks++;

        // Get the DUT clusters from the clipboard
582
        auto clusters = clipboard->getData<Cluster>(m_detector->name());
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
        if(clusters == nullptr) {
            LOG(DEBUG) << " - no DUT clusters";
        } else {

            // Loop over all DUT clusters to find matches:
            for(auto* cluster : (*clusters)) {
                LOG(DEBUG) << " - Looking at next DUT cluster";

                hTrackCorrelationTime->Fill(track->timestamp() - cluster->timestamp());

                auto associated_clusters = track->associatedClusters();
                if(std::find(associated_clusters.begin(), associated_clusters.end(), cluster) != associated_clusters.end()) {
                    LOG(DEBUG) << "Found associated cluster " << (*cluster);
                    has_associated_cluster = true;
                    matched_tracks++;

599
                    if(m_config.has("cluster_charge_cut") && cluster->charge() > m_clusterChargeCut) {
600
                        LOG(DEBUG) << " - track discarded due to clusterChargeCut";
601
602
                        continue;
                    }
603
                    tracks_afterClusterChargeCut++;
604

605
                    if(m_config.has("cluster_size_cut") && cluster->size() > m_clusterSizeCut) {
606
607
608
609
610
611
612
                        LOG(DEBUG) << " - track discarded due to clusterSizeCut";
                        continue;
                    }
                    tracks_afterClusterSizeCut++;

                    double timeDiff = track->timestamp() - cluster->timestamp();
                    hTrackCorrelationTimeAssoc->Fill(timeDiff);
613
614
                    hTrackCorrelationTimeAssocVsTime->Fill(static_cast<double>(Units::convert(cluster->timestamp(), "s")),
                                                           static_cast<double>(Units::convert(timeDiff, "us")));
615

616
                    hTrackCorrelationTimeVsTot->Fill(timeDiff, cluster->getSeedPixel()->raw());
617
618
                    hTrackCorrelationTimeVsCol->Fill(timeDiff, cluster->getSeedPixel()->column());
                    hTrackCorrelationTimeVsRow->Fill(timeDiff, cluster->getSeedPixel()->row());
619

620
621
                    hTotVsRow->Fill(cluster->getSeedPixel()->row(), cluster->getSeedPixel()->raw());

622
                    // single-pixel and multi-pixel clusters:
623
                    if(cluster->size() == 1) {
624
                        hTrackCorrelationTimeVsTot_1px->Fill(timeDiff, cluster->getSeedPixel()->raw());
625
                        hTrackCorrelationTimeVsRow_1px->Fill(timeDiff, cluster->getSeedPixel()->row());
626
                    } else {
627
                        hTrackCorrelationTimeVsTot_npx->Fill(timeDiff, cluster->getSeedPixel()->raw());
628
                        hTrackCorrelationTimeVsRow_npx->Fill(timeDiff, cluster->getSeedPixel()->row());
629
630
631
632
633
634
635
636
637
                    }

                    // Calculate in-pixel position of track in microns
                    auto globalIntercept = m_detector->getIntercept(track);
                    auto localIntercept = m_detector->globalToLocal(globalIntercept);
                    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"));
                    hHitMapAssoc_inPixel->Fill(xmod, ymod);
638
639
                    if(m_config.has("high_tot_cut") && cluster->charge() > m_highTotCut && cluster->size() == 1) {
                        hHitMapAssoc_inPixel_highToT->Fill(xmod, ymod);
640
                    }
641
                    hPixelTrackCorrelationTimeMap->Fill(xmod, ymod, timeDiff);
642
643

                    // 2D histograms: --> fill for all pixels from cluster
644
                    for(auto& pixel : cluster->pixels()) {
645

646
647
                        hTrackCorrelationTimeVsTot_px->Fill(track->timestamp() - pixel->timestamp(), pixel->raw());

648
649
650
651
652
                        // to check that cluster timestamp = earliest pixel timestamp
                        if(cluster->size() > 1) {
                            hClusterTimeMinusPixelTime->Fill(cluster->timestamp() - pixel->timestamp());
                        }

653
                        hClusterSizeVsTot_Assoc->Fill(static_cast<double>(cluster->size()), pixel->raw());
654
                        hHitMapAssoc->Fill(pixel->column(), pixel->row());
655
                        hTotVsTime->Fill(pixel->raw(), static_cast<double>(Units::convert(pixel->timestamp(), "s")));
656
657
                        if(m_config.has("high_tot_cut") && pixel->raw() > m_highTotCut) {
                            hHitMapAssoc_highToT->Fill(pixel->column(), pixel->row());
658
659
                            hTotVsTime_highToT->Fill(pixel->raw(),
                                                     static_cast<double>(Units::convert(pixel->timestamp(), "s")));
660
661
662
                        }
                    }
                    hClusterMapAssoc->Fill(cluster->column(), cluster->row());
663
664
665
666
667
668
669
670
                    if(m_config.has("timing_tail_cut") && m_pointwise_correction_timewalk) {
                        if(track->timestamp() - cluster->timestamp() < -m_timingTailCut) {
                            hInPixelMap_leftTail->Fill(xmod, ymod);
                        } else if(track->timestamp() - cluster->timestamp() > m_timingTailCut) {
                            hInPixelMap_rightTail->Fill(xmod, ymod);
                        } else {
                            hInPixelMap_mainPeak->Fill(xmod, ymod);
                        }
671
                    }
672
673
674
675
676
677
678
679
680
681
682
683

                    // !!! Have to do this in the end because it changes the cluster time and position!!!
                    // row-by-row correction using points from TGraphError directly instead of fit.

                    // point-wise correction:
                    if(m_pointwise_correction_row) {
                        correctClusterTimestamp(cluster, 0); // mode=0 --> row correction
                        hTrackCorrelationTime_rowCorr->Fill(track->timestamp() - cluster->timestamp());
                        // for(auto& pixel : (*cluster->pixels())) {
                        hTrackCorrelationTimeVsRow_rowCorr->Fill(track->timestamp() - cluster->getSeedPixel()->timestamp(),
                                                                 cluster->getSeedPixel()->row());
                        hTrackCorrelationTimeVsTot_rowCorr->Fill(track->timestamp() - cluster->getSeedPixel()->timestamp(),
684
                                                                 cluster->getSeedPixel()->raw());
685
686
                        if(cluster->size() == 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_1px->Fill(
687
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
688
689
690
                        }
                        if(cluster->size() > 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_npx->Fill(
691
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
692
693
694
695
696
697
                        }
                        //} for(auto& pixels : ...)
                    }
                    // point-wise timewalk correction on top:
                    if(m_pointwise_correction_timewalk) {
                        correctClusterTimestamp(cluster, 1); // mode=1 --> timewalk correction
698
                        hTrackCorrelationTime_rowAndTWCorr->Fill(track->timestamp() - cluster->timestamp());
699
                        if(cluster->getSeedPixel()->raw() < 25) {
700
                            hTrackCorrelationTime_rowAndTWCorr_l25->Fill(track->timestamp() - cluster->timestamp());
701
                        }
702
                        if(cluster->getSeedPixel()->raw() < 40) {
703
                            hTrackCorrelationTime_rowAndTWCorr_l40->Fill(track->timestamp() - cluster->timestamp());
704
                        }
705
                        if(cluster->getSeedPixel()->raw() > 40) {
706
                            hTrackCorrelationTime_rowAndTWCorr_g40->Fill(track->timestamp() - cluster->timestamp());
707
708
                        }

709
                        hTrackCorrelationTimeVsRow_rowAndTWCorr->Fill(
710
                            track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->row());
711
                        hTrackCorrelationTimeVsTot_rowAndTWCorr->Fill(
712
                            track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
713

714
                        // control plots to investigate "left/right tail" in time correlation:
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
                        if(m_config.has("timing_tail_cut")) {
                            if(track->timestamp() - cluster->timestamp() < -m_timingTailCut) {
                                hClusterMap_leftTail->Fill(cluster->column(), cluster->row());
                                hTot_leftTail->Fill(cluster->getSeedPixel()->raw());
                                hPixelTimestamp_leftTail->Fill(cluster->getSeedPixel()->timestamp());
                                hClusterSize_leftTail->Fill(static_cast<double>(cluster->size()));
                                if(cluster->size() == 1) {
                                    hTot_leftTail_1px->Fill(cluster->getSeedPixel()->raw());
                                }
                            } else if(track->timestamp() - cluster->timestamp() > m_timingTailCut) {
                                hClusterMap_rightTail->Fill(cluster->column(), cluster->row());
                                hTot_rightTail->Fill(cluster->getSeedPixel()->raw());
                                hPixelTimestamp_rightTail->Fill(cluster->getSeedPixel()->timestamp());
                                hClusterSize_rightTail->Fill(static_cast<double>(cluster->size()));
                                if(cluster->size() == 1) {
                                    hTot_rightTail_1px->Fill(cluster->getSeedPixel()->raw());
                                }
                            } else {
                                hClusterMap_mainPeak->Fill(cluster->column(), cluster->row());
                                hTot_mainPeak->Fill(cluster->getSeedPixel()->raw());
                                hPixelTimestamp_mainPeak->Fill(cluster->getSeedPixel()->timestamp());
                                hClusterSize_mainPeak->Fill(static_cast<double>(cluster->size()));
                                if(cluster->size() == 1) {
                                    hTot_mainPeak_1px->Fill(cluster->getSeedPixel()->raw());
                                }
740
                            }
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
                        }
                    }
                }

            } // for loop over all clusters
        }     // else (clusters != nullptr)

        LOG(DEBUG) << "is_within_roi = " << is_within_roi;
        LOG(DEBUG) << "has_associated_cluster = " << has_associated_cluster;

    } // for loop over all tracks

    return StatusCode::Success;
}

756
void AnalysisTimingATLASpix::finalise() {
757
758
759
760
761
762
763
764
765
766
767
768
769
770
    LOG(STATUS) << "Timing analysis finished for detector " << m_detector->name() << ": ";

    if(m_calcCorrections) {

        /// ROW CORRECTION ///
        std::string fitOption = "q"; // set to "" for terminal output
        int binMax = 0;
        double timePeak = 0.;
        double timePeakErr = 0.;
        int nRows = hTrackCorrelationTimeVsRow->GetNbinsY();

        for(int iBin = 0; iBin < nRows; iBin++) {
            TH1D* hTemp = hTrackCorrelationTimeVsRow->ProjectionX("timeCorrelationInOneTotBin", iBin, iBin + 1);

771
772
            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 250) { // too few entries to fit
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
                delete hTemp;
                timePeak = 0;
                timePeakErr = 0;
                continue;
            } else {
                binMax = hTemp->GetMaximumBin();
                timePeak = hTemp->GetXaxis()->GetBinCenter(binMax);

                // fitting a Gaus for a good estimate of the peak positon:
                // NOTE: initial values for Gaussian are hard-coded at the moment!
                TF1* fPeak = new TF1("fPeak", "gaus");
                fPeak->SetParameters(1, 100, 45);
                double timeInt = 50;
                hTemp->Fit("fPeak", fitOption.c_str(), "", timePeak - timeInt, timePeak + timeInt);
                fPeak = hTemp->GetFunction("fPeak");
                timePeak = fPeak->GetParameter(1);
                timePeakErr = fPeak->GetParError(1);
                delete fPeak;
                delete hTemp;
            }
            // TGraphErrors should only have as many bins as it has sensible entries
            // (If it has multiple x=0 entries, the Spline interpolation will fail.
            int nBins = gTimeCorrelationVsRow->GetN();
            LOG(STATUS) << "nBins = " << nBins << ", x = " << iBin << ", y = " << timePeak;
            gTimeCorrelationVsRow->SetPoint(nBins, iBin, timePeak);
            gTimeCorrelationVsRow->SetPointError(nBins, 0., timePeakErr);

        } // for(iBin)

        /// TIME WALK CORRECTION on top of ROW CORRECTION: ///
        fitOption = "q"; // set to "" if you want terminal output
        binMax = 0;
        timePeak = 0.;
        timePeakErr = 0.;
        int nBinsToT = hTrackCorrelationTimeVsTot_rowCorr->GetNbinsY();
        LOG(DEBUG) << "nBinsToT = " << nBinsToT;

        for(int iBin = 0; iBin < nBinsToT; iBin++) {
            TH1D* hTemp = hTrackCorrelationTimeVsTot_rowCorr->ProjectionX("timeCorrelationInOneTotBin", iBin, iBin + 1);

            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 1000) { // too few entries to fit
815
816
                timePeak = hTemp->GetMean();
                timePeakErr = 0; // hTemp->GetStdDev();
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
                delete hTemp;
            } else {
                binMax = hTemp->GetMaximumBin();
                timePeak = hTemp->GetXaxis()->GetBinCenter(binMax);

                // fitting a Gaus for a good estimate of the peak positon:
                // initial parameters are hardcoded at the moment!
                TF1* fPeak = new TF1("fPeak", "gaus");
                fPeak->SetParameters(1, 100, 45);
                double timeInt = 50;
                hTemp->Fit("fPeak", fitOption.c_str(), "", timePeak - timeInt, timePeak + timeInt);
                fPeak = hTemp->GetFunction("fPeak");
                timePeak = fPeak->GetParameter(1);
                timePeakErr = fPeak->GetParError(1);

                delete fPeak;
                delete hTemp;
            }
            gTimeCorrelationVsTot_rowCorr->SetPoint(iBin, iBin, timePeak);
            gTimeCorrelationVsTot_rowCorr->SetPointError(iBin, 0, timePeakErr);

        } // for(iBin)

        // SAME FOR SINGLE-PIXEL CLUSTERS:
        nBinsToT = hTrackCorrelationTimeVsTot_rowCorr_1px->GetNbinsY();
        for(int iBin = 0; iBin < nBinsToT; iBin++) {
            TH1D* hTemp = hTrackCorrelationTimeVsTot_rowCorr_1px->ProjectionX("timeCorrelationInOneTotBin", iBin, iBin + 1);

            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 1000) { // too few entries to fit
                delete hTemp;
                timePeak = 0;
                timePeakErr = 0;
                continue;
            } else {
                binMax = hTemp->GetMaximumBin();
                timePeak = hTemp->GetXaxis()->GetBinCenter(binMax);

                // fitting a Gaus for a good estimate of the peak positon:
                // initial parameters are hardcoded at the moment!
                TF1* fPeak = new TF1("fPeak", "gaus");
                fPeak->SetParameters(1, 100, 45);
                double timeInt = 50;
                hTemp->Fit("fPeak", fitOption.c_str(), "", timePeak - timeInt, timePeak + timeInt);
                fPeak = hTemp->GetFunction("fPeak");
                timePeak = fPeak->GetParameter(1);
                timePeakErr = fPeak->GetParError(1);

                delete fPeak;
                delete hTemp;
            }
            gTimeCorrelationVsTot_rowCorr_1px->SetPoint(iBin, iBin, timePeak);
            gTimeCorrelationVsTot_rowCorr_1px->SetPointError(iBin, 0, timePeakErr);
        } // for(iBin)

        // SAME FOR MULTI-PIXEL CLUSTERS:
        nBinsToT = hTrackCorrelationTimeVsTot_rowCorr_npx->GetNbinsY();
        for(int iBin = 0; iBin < nBinsToT; iBin++) {
            TH1D* hTemp = hTrackCorrelationTimeVsTot_rowCorr_npx->ProjectionX("timeCorrelationInOneTotBin", iBin, iBin + 1);

            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 1000) { // too few entries to fit
                delete hTemp;
                timePeak = 0;
                timePeakErr = 0;
                continue;
            } else {
                binMax = hTemp->GetMaximumBin();
                timePeak = hTemp->GetXaxis()->GetBinCenter(binMax);

                // fitting a Gaus for a good estimate of the peak positon:
                // initial parameters are hardcoded at the moment!
                TF1* fPeak = new TF1("fPeak", "gaus");
                fPeak->SetParameters(1, 100, 45);
                double timeInt = 50;
                hTemp->Fit("fPeak", fitOption.c_str(), "", timePeak - timeInt, timePeak + timeInt);
                fPeak = hTemp->GetFunction("fPeak");
                timePeak = fPeak->GetParameter(1);
                timePeakErr = fPeak->GetParError(1);

                delete fPeak;
                delete hTemp;
            }
            gTimeCorrelationVsTot_rowCorr_npx->SetPoint(iBin, iBin, timePeak);
            gTimeCorrelationVsTot_rowCorr_npx->SetPointError(iBin, 0, timePeakErr);
        } // for(iBin)

        /// END TIME WALK CORRECTION ///

906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
        // Example Slice to investigate quality of Gaussian fit:
        hTrackCorrelationTime_example = hTrackCorrelationTimeVsTot_rowCorr->ProjectionX(
            ("hTrackCorrelationTime_totBin_" + std::to_string(m_totBinExample)).c_str(),
            m_totBinExample,
            m_totBinExample + 1);

        binMax = hTrackCorrelationTime_example->GetMaximumBin();
        timePeak = hTrackCorrelationTime_example->GetXaxis()->GetBinCenter(binMax);

        TF1* fPeak = new TF1("fPeak", "gaus");
        fPeak->SetParameters(1, 100, 45);
        double timeInt = 50;
        fitOption = "q"; // set to "q" = quiet for suppressed terminial output
        hTrackCorrelationTime_example->Fit("fPeak", fitOption.c_str(), "", timePeak - timeInt, timePeak + timeInt);
        delete fPeak;

        gTimeCorrelationVsRow->Write();
        gTimeCorrelationVsTot_rowCorr->Write();
        gTimeCorrelationVsTot_rowCorr_1px->Write();
        gTimeCorrelationVsTot_rowCorr_npx->Write();
926
927
928
929
930
931
932
933
    } // if(m_calcCorrections)

    LOG(INFO) << "matched/total tracks: " << matched_tracks << "/" << total_tracks;
    LOG(INFO) << "total tracks (uncut):\t" << total_tracks_uncut;
    LOG(INFO) << "after chi2 cut:\t" << tracks_afterChi2Cut;
    LOG(INFO) << "with intercept:\t" << tracks_hasIntercept;
    LOG(INFO) << "withing ROI:\t\t" << tracks_isWithinROI;
    LOG(INFO) << "frameEdge cut:\t\t" << matched_tracks;
934
    LOG(INFO) << "after clusterTotCut:\t" << tracks_afterClusterChargeCut;
935
936
937
    LOG(INFO) << "after clusterSizeCut:\t" << tracks_afterClusterSizeCut;
}

938
void AnalysisTimingATLASpix::correctClusterTimestamp(Cluster* cluster, int mode) {
939

940
    /*
941
942
943
944
945
946
     * MODE:
     *  0 --> row correction
     *  1 --> timewalk correction
     */

    // Get the pixels on this cluster
947
948
    auto pixels = cluster->pixels();
    auto first_pixel = pixels.front();
949
    double correction = 0;
950
951

    if(mode == 0) {
952
        correction = gRowCorr->Eval(first_pixel->row());
953
    } else if(mode == 1) {
954
        correction = gTimeWalkCorr->Eval(first_pixel->raw());
955
956
957
958
959
960
961
    } else {
        LOG(ERROR) << "Mode " << mode << " does not exist!\n"
                   << "Choose\n\t0 --> row correction \n\t1-->timewalk correction";
        return;
    }

    // Initial guess for cluster timestamp:
962
    double timestamp = first_pixel->timestamp() + correction;
963
964

    // Loop over all pixels:
965
966
967
    for(auto& pixel : pixels) {
        // FIXME ugly hack
        auto px = const_cast<Pixel*>(pixel);
968
969
970
971

        if(mode == 0) {
            correction = gRowCorr->Eval(pixel->row());
        } else if(mode == 1) {
972
            correction = gTimeWalkCorr->Eval(pixel->raw());
973
974
975
976
977
        } else {
            return;
        }

        // Override pixel timestamps:
978
        px->setTimestamp(pixel->timestamp() + correction);
979
980
981
982
983
984
985
986
987

        // timestamp = earliest pixel:
        if(pixel->timestamp() < timestamp) {
            timestamp = pixel->timestamp();
        }
    }

    cluster->setTimestamp(timestamp);
}