AnalysisTimingATLASpix.cpp 50.4 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
189

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

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

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

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

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

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

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

    // final plots with both row and timewalk correction:
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
    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";
        hTrackCorrelationTimeVsTot_rowAndTWCorr = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 512, 0, 512);
        hTrackCorrelationTimeVsTot_rowAndTWCorr->GetYaxis()->SetTitle("seed pixel tot [lsb]}");
        hTrackCorrelationTimeVsTot_rowAndTWCorr->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
    }
245

246
    // in-pixel time resolution plots:
247
248
    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()));
249
250
251
252
    if(nbins_x > 1e4 || nbins_y > 1e4) {
        throw InvalidValueError(m_config, "inpixel_bin_size", "Too many bins for in-pixel histograms.");
    }

253
254
    std::string title =
        "in-pixel time resolution map;in-pixel x_{track} [#mum];in-pixel y_{track} [#mum];ts_{track} - ts_{cluster} [ns]";
255
256
257
258
259
260
261
262
263
264
265
    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);

266
267
268
269
270
271
    name = "hClusterSizeVsTot_Assoc";
    hClusterSizeVsTot_Assoc = new TH2F(name.c_str(), name.c_str(), 20, 0, 20, 512, 0, 512);
    hClusterSizeVsTot_Assoc->GetYaxis()->SetTitle("pixel ToT [ns] (all pixels from cluster)");
    hClusterSizeVsTot_Assoc->GetXaxis()->SetTitle("clusterSize");

    hHitMapAssoc = new TH2F("hitMapAssoc",
272
                            "hitMapAssoc; x_{track} [px]; y_{track} [px];# entries",
273
                            m_detector->nPixels().X(),
274
275
                            -0.5,
                            m_detector->nPixels().X() - 0.5,
276
                            m_detector->nPixels().Y(),
277
278
279
280
281
282
283
284
285
286
                            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);
287
    hHitMapAssoc_inPixel = new TH2F("hitMapAssoc_inPixel",
288
                                    "hitMapAssoc_inPixel; in-pixel x_{track} [#mum]; in-pixel y_{track} [#mum]",
289
                                    static_cast<int>(pitch_x),
290
291
                                    -pitch_x / 2.,
                                    pitch_x / 2.,
292
                                    static_cast<int>(pitch_y),
293
294
                                    -pitch_y / 2.,
                                    pitch_y / 2.);
295
296
297
298
299
300
301
302
303
304
305
    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.);
    }
306
    hClusterMapAssoc = new TH2F("hClusterMapAssoc",
307
                                "hClusterMapAssoc; x_{cluster} [px]; y_{cluster} [px];# entries",
308
                                m_detector->nPixels().X(),
309
310
                                -0.5,
                                m_detector->nPixels().X() - 0.5,
311
                                m_detector->nPixels().Y(),
312
313
                                -0.5,
                                m_detector->nPixels().Y() - 0.5);
314

315
316
317
318
319
320
321
322
323
    hTotVsRow = new TH2F("hTotVsRow",
                         "hTotVsRow;seed-pixel row;seed-pixel ToT [lsb]",
                         m_detector->nPixels().Y(),
                         0,
                         m_detector->nPixels().Y(),
                         64,
                         0,
                         64);

324
325
326
    hTotVsTime = new TH2F("hTotVsTime", "hTotVsTime", 64, 0, 64, 1e6, 0, 100);
    hTotVsTime->GetXaxis()->SetTitle("pixel ToT [lsb]");
    hTotVsTime->GetYaxis()->SetTitle("time [s]");
327
328
329
330
331
332
333
334
335
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
    if(m_config.has("high_tot_cut")) {
        hTotVsTime_high = new TH2F("hTotVsTime_high", "hTotVsTime_high", 64, 0, 64, 1e6, 0, 100);
        hTotVsTime_high->GetXaxis()->SetTitle("pixel ToT [lsb] if > high_tot_cut");
        hTotVsTime_high->GetYaxis()->SetTitle("time [s]");
    }

    // 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",
401
402
403
                                 2 * 64,
                                 -64,
                                 64);
404
405
406
407
408
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
        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);
    }
435
436
437
438
439

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

440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
    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;
470
471
472

    if(m_pointwise_correction_row) {
        // Import TGraphErrors for row corection:
473
474
475
476
477
478
        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.");
479
        }
480
481

        gRowCorr = static_cast<TGraphErrors*>(file.Get(m_correctionGraph_row.c_str()));
482
        // Check if graph exists in ROOT file:
483
        if(!gRowCorr) {
484
485
486
            throw InvalidValueError(
                m_config, "correction_graph_row", "Graph doesn't exist in ROOT file. Use full/path/to/graph.");
        }
487
488
489
490
491
    } else {
        LOG(STATUS) << "----> NO POINTWISE ROW CORRECTION!!!";
    }
    if(m_pointwise_correction_timewalk) {
        // Import TGraphErrors for timewalk corection:
492
493
494
495
496
497
498
        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.");
        }
499
500

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

511
StatusCode AnalysisTimingATLASpix::run(std::shared_ptr<Clipboard> clipboard) {
512
513

    // Get the telescope tracks from the clipboard
514
    auto tracks = clipboard->getData<Track>();
515
516
517
518
519
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
    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:
556
        auto event = clipboard->getEvent();
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576

        // 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
577
        auto clusters = clipboard->getData<Cluster>(m_detector->name());
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
        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++;

594
                    if(m_config.has("cluster_charge_cut") && cluster->charge() > m_clusterChargeCut) {
595
                        LOG(DEBUG) << " - track discarded due to clusterChargeCut";
596
597
                        continue;
                    }
598
                    tracks_afterClusterChargeCut++;
599

600
                    if(m_config.has("cluster_size_cut") && cluster->size() > m_clusterSizeCut) {
601
602
603
604
605
606
607
                        LOG(DEBUG) << " - track discarded due to clusterSizeCut";
                        continue;
                    }
                    tracks_afterClusterSizeCut++;

                    double timeDiff = track->timestamp() - cluster->timestamp();
                    hTrackCorrelationTimeAssoc->Fill(timeDiff);
608
609
                    hTrackCorrelationTimeAssocVsTime->Fill(static_cast<double>(Units::convert(cluster->timestamp(), "s")),
                                                           static_cast<double>(Units::convert(timeDiff, "us")));
610

611
                    hTrackCorrelationTimeVsTot->Fill(timeDiff, cluster->getSeedPixel()->raw());
612
613
                    hTrackCorrelationTimeVsCol->Fill(timeDiff, cluster->getSeedPixel()->column());
                    hTrackCorrelationTimeVsRow->Fill(timeDiff, cluster->getSeedPixel()->row());
614

615
616
                    hTotVsRow->Fill(cluster->getSeedPixel()->row(), cluster->getSeedPixel()->raw());

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

                    // 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);
633
634
                    if(m_config.has("high_tot_cut") && cluster->charge() > m_highTotCut && cluster->size() == 1) {
                        hHitMapAssoc_inPixel_highToT->Fill(xmod, ymod);
635
                    }
636
                    hPixelTrackCorrelationTimeMap->Fill(xmod, ymod, timeDiff);
637
638

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

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

646
                        hClusterSizeVsTot_Assoc->Fill(static_cast<double>(cluster->size()), pixel->raw());
647
                        hHitMapAssoc->Fill(pixel->column(), pixel->row());
648
                        hTotVsTime->Fill(pixel->raw(), static_cast<double>(Units::convert(pixel->timestamp(), "s")));
649
650
                        if(m_config.has("high_tot_cut") && pixel->raw() > m_highTotCut) {
                            hHitMapAssoc_highToT->Fill(pixel->column(), pixel->row());
651
                            hTotVsTime_high->Fill(pixel->raw(),
652
653
654
655
                                                  static_cast<double>(Units::convert(pixel->timestamp(), "s")));
                        }
                    }
                    hClusterMapAssoc->Fill(cluster->column(), cluster->row());
656
657
658
659
660
661
662
663
                    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);
                        }
664
                    }
665
666
667
668
669
670
671
672
673
674
675
676

                    // !!! 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(),
677
                                                                 cluster->getSeedPixel()->raw());
678
679
                        if(cluster->size() == 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_1px->Fill(
680
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
681
682
683
                        }
                        if(cluster->size() > 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_npx->Fill(
684
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
685
686
687
688
689
690
                        }
                        //} for(auto& pixels : ...)
                    }
                    // point-wise timewalk correction on top:
                    if(m_pointwise_correction_timewalk) {
                        correctClusterTimestamp(cluster, 1); // mode=1 --> timewalk correction
691
                        hTrackCorrelationTime_rowAndTWCorr->Fill(track->timestamp() - cluster->timestamp());
692
                        if(cluster->getSeedPixel()->raw() < 25) {
693
                            hTrackCorrelationTime_rowAndTWCorr_l25->Fill(track->timestamp() - cluster->timestamp());
694
                        }
695
                        if(cluster->getSeedPixel()->raw() < 40) {
696
                            hTrackCorrelationTime_rowAndTWCorr_l40->Fill(track->timestamp() - cluster->timestamp());
697
                        }
698
                        if(cluster->getSeedPixel()->raw() > 40) {
699
                            hTrackCorrelationTime_rowAndTWCorr_g40->Fill(track->timestamp() - cluster->timestamp());
700
701
                        }

702
                        hTrackCorrelationTimeVsRow_rowAndTWCorr->Fill(
703
                            track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->row());
704
                        hTrackCorrelationTimeVsTot_rowAndTWCorr->Fill(
705
                            track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
706

707
                        // control plots to investigate "left/right tail" in time correlation:
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
                        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());
                                }
733
                            }
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
                        }
                    }
                }

            } // 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;
}

749
void AnalysisTimingATLASpix::finalise() {
750
751
752
753
754
755
756
757
758
759
760
761
762
763
    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);

764
765
            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 250) { // too few entries to fit
766
767
768
769
770
771
772
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
                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
808
809
                timePeak = hTemp->GetMean();
                timePeakErr = 0; // hTemp->GetStdDev();
810
811
812
813
814
815
816
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
                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 ///

899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
        // 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();
919
920
921
922
923
924
925
926
    } // 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;
927
    LOG(INFO) << "after clusterTotCut:\t" << tracks_afterClusterChargeCut;
928
929
930
    LOG(INFO) << "after clusterSizeCut:\t" << tracks_afterClusterSizeCut;
}

931
void AnalysisTimingATLASpix::correctClusterTimestamp(Cluster* cluster, int mode) {
932

933
    /*
934
935
936
937
938
939
     * MODE:
     *  0 --> row correction
     *  1 --> timewalk correction
     */

    // Get the pixels on this cluster
940
941
    auto pixels = cluster->pixels();
    auto first_pixel = pixels.front();
942
    double correction = 0;
943
944

    if(mode == 0) {
945
        correction = gRowCorr->Eval(first_pixel->row());
946
    } else if(mode == 1) {
947
        correction = gTimeWalkCorr->Eval(first_pixel->raw());
948
949
950
951
952
953
954
    } else {
        LOG(ERROR) << "Mode " << mode << " does not exist!\n"
                   << "Choose\n\t0 --> row correction \n\t1-->timewalk correction";
        return;
    }

    // Initial guess for cluster timestamp:
955
    double timestamp = first_pixel->timestamp() + correction;
956
957

    // Loop over all pixels:
958
959
960
    for(auto& pixel : pixels) {
        // FIXME ugly hack
        auto px = const_cast<Pixel*>(pixel);
961
962
963
964

        if(mode == 0) {
            correction = gRowCorr->Eval(pixel->row());
        } else if(mode == 1) {
965
            correction = gTimeWalkCorr->Eval(pixel->raw());
966
967
968
969
970
        } else {
            return;
        }

        // Override pixel timestamps:
971
        px->setTimestamp(pixel->timestamp() + correction);
972
973
974
975
976
977
978
979
980

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

    cluster->setTimestamp(timestamp);
}