AnalysisTimingATLASpix.cpp 46.6 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_highChargeCut = m_config.get<double>("high_charge_cut", static_cast<double>(m_highTotCut));
49
    m_leftTailCut = m_config.get<double>("left_tail_cut", static_cast<double>(Units::convert(-10, "ns")));
50
    m_rightTailCut = m_config.get<double>("right_tail_cut", static_cast<double>(Units::convert(10, "ns")));
51
52
53

    if(m_config.has("correction_file_row")) {
        m_correctionFile_row = m_config.get<std::string>("correction_file_row");
54
        m_correctionGraph_row = m_config.get<std::string>("correction_graph_row");
55
56
57
58
59
60
        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");
61
        m_correctionGraph_timewalk = m_config.get<std::string>("correction_graph_timewalk");
62
63
64
65
66
67
68
69
        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);

70
71
    m_inpixelBinSize = m_config.get<double>("inpixel_bin_size", Units::get<double>(1.0, "um"));

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

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

    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 =
90
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
91
    hTrackCorrelationTime->GetXaxis()->SetTitle("Track time stamp - cluster time stamp [ns]");
92
    hTrackCorrelationTime->GetYaxis()->SetTitle("# events");
93
94
95

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

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

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

    name = "hTrackCorrelationTime_rowAndTimeWalkCorr";
    hTrackCorrelationTime_rowAndTimeWalkCorr =
115
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
116
    hTrackCorrelationTime_rowAndTimeWalkCorr->GetXaxis()->SetTitle("track time stamp - cluster time stamp [ns]");
117
    hTrackCorrelationTime_rowAndTimeWalkCorr->GetYaxis()->SetTitle("# events");
118
119
120

    name = "hTrackCorrelationTime_rowAndTimeWalkCorr_l25";
    hTrackCorrelationTime_rowAndTimeWalkCorr_l25 =
121
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
122
123
    hTrackCorrelationTime_rowAndTimeWalkCorr_l25->GetXaxis()->SetTitle(
        "track time stamp - cluster time stamp [ns] (if seed tot < 25lsb)");
124
    hTrackCorrelationTime_rowAndTimeWalkCorr_l25->GetYaxis()->SetTitle("# events");
125
126
127

    name = "hTrackCorrelationTime_rowAndTimeWalkCorr_l40";
    hTrackCorrelationTime_rowAndTimeWalkCorr_l40 =
128
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
129
130
    hTrackCorrelationTime_rowAndTimeWalkCorr_l40->GetXaxis()->SetTitle(
        "track time stamp - cluster time stamp [ns] (if seed tot < 40lsb)");
131
    hTrackCorrelationTime_rowAndTimeWalkCorr_l40->GetYaxis()->SetTitle("# events");
132
133
134

    name = "hTrackCorrelationTime_rowAndTimeWalkCorr_g40";
    hTrackCorrelationTime_rowAndTimeWalkCorr_g40 =
135
        new TH1F(name.c_str(), name.c_str(), static_cast<int>(2. * m_timeCut), -1 * m_timeCut, m_timeCut);
136
137
    hTrackCorrelationTime_rowAndTimeWalkCorr_g40->GetXaxis()->SetTitle(
        "track time stamp - cluster time stamp [ns] (if seed tot > 40lsb)");
138
    hTrackCorrelationTime_rowAndTimeWalkCorr_g40->GetYaxis()->SetTitle("# events");
139
140
141
142
143

    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)");
144
    hTrackCorrelationTime_example->GetYaxis()->SetTitle("# events");
145
146
147
148

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

    // control plot: row dependence after row correction
    name = "hTrackCorrelationTimeVsRow_rowCorr";
174
175
    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);
176
177
178
179
180
181
    hTrackCorrelationTimeVsRow_rowCorr->GetYaxis()->SetTitle("pixel row");
    hTrackCorrelationTimeVsRow_rowCorr->GetXaxis()->SetTitle("track time stamp - seed pixel time stamp [ns]");

    // 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);
182
183
    hTrackCorrelationTimeVsTot->GetYaxis()->SetTitle("seed pixel ToT [ns]");
    hTrackCorrelationTimeVsTot->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
184
185
186
187

    name = "hTrackCorrelationTimeVsTot_1px";
    hTrackCorrelationTimeVsTot_1px = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 512, 0, 512);
    hTrackCorrelationTimeVsTot_1px->GetYaxis()->SetTitle("seed pixel ToT [ns] (if clustersize = 1)");
188
    hTrackCorrelationTimeVsTot_1px->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
189
190
191
192

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

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

    // timewalk after row correction
    name = "hTrackCorrelationTimeVsTot_rowCorr";
    hTrackCorrelationTimeVsTot_rowCorr = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 512, 0, 512);
203
    hTrackCorrelationTimeVsTot_rowCorr->GetYaxis()->SetTitle("seed pixel ToT [ns]");
204
205
206
207
    hTrackCorrelationTimeVsTot_rowCorr->GetXaxis()->SetTitle("track time stamp - seed pixel time stamp [ns]");

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

    name = "hTrackCorrelationTimeVsTot_rowCorr_npx";
    hTrackCorrelationTimeVsTot_rowCorr_npx = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 512, 0, 512);
213
214
    hTrackCorrelationTimeVsTot_rowCorr_npx->GetYaxis()->SetTitle("seed pixel ToT [ns] (multi-pixel clusters)");
    hTrackCorrelationTimeVsTot_rowCorr_npx->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
215
216
217

    // final plots with both row and timewalk correction:
    name = "hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr";
218
219
    hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr = new TH2F(
        name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), -0.5, m_detector->nPixels().Y() - 0.5);
220
221
222
223
224
    hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr->GetYaxis()->SetTitle("row");
    hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr->GetXaxis()->SetTitle("track time stamp - seed pixel time stamp [ns]");

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

228
229
230
231
232
233
234
    // in-pixel time resolution plots:
    auto nbins_x = static_cast<int>(std::ceil(m_detector->pitch().X() / m_inpixelBinSize));
    auto nbins_y = static_cast<int>(std::ceil(m_detector->pitch().Y() / m_inpixelBinSize));
    if(nbins_x > 1e4 || nbins_y > 1e4) {
        throw InvalidValueError(m_config, "inpixel_bin_size", "Too many bins for in-pixel histograms.");
    }

235
    title = "in-pixel time resolution map;x [px];y [px];ts_{track} - ts_{cluster} [ns]";
236
237
238
239
240
241
242
243
244
245
246
    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);

247
248
249
250
251
252
    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",
253
                            "hitMapAssoc; x_{track} [px]; x_{track} [px]; # entries",
254
                            m_detector->nPixels().X(),
255
256
                            -0.5,
                            m_detector->nPixels().X() - 0.5,
257
                            m_detector->nPixels().Y(),
258
259
                            -0.5,
                            m_detector->nPixels().Y() - 0.5);
260
261
262
    hHitMapAssoc_highCharge = new TH2F("hitMapAssoc_highCharge",
                                       "hitMapAssoc_highCharge; x_{track} [px]; x_{track} [px]; # entries",
                                       m_detector->nPixels().X(),
263
264
                                       -0.5,
                                       m_detector->nPixels().X() - 0.5,
265
                                       m_detector->nPixels().Y(),
266
267
                                       -0.5,
                                       m_detector->nPixels().Y() - 0.5);
268
    hHitMapAssoc_inPixel = new TH2F("hitMapAssoc_inPixel",
269
                                    "hitMapAssoc_inPixel; in-pixel x_{track} [#mum]; in-pixel y_{track} [#mum]",
270
                                    static_cast<int>(pitch_x),
271
272
                                    -pitch_x / 2.,
                                    pitch_x / 2.,
273
                                    static_cast<int>(pitch_y),
274
275
                                    -pitch_y / 2.,
                                    pitch_y / 2.);
276
277
    hHitMapAssoc_inPixel_highCharge =
        new TH2F("hitMapAssoc_inPixel_highCharge",
278
                 "hitMapAssoc_inPixel_highCharge;in-pixel x_{track} [#mum];in-pixel y_{track} [#mum]",
279
                 static_cast<int>(pitch_x),
280
281
                 -pitch_x / 2.,
                 pitch_x / 2.,
282
                 static_cast<int>(pitch_y),
283
284
                 -pitch_y / 2.,
                 pitch_y / 2.);
285
    hClusterMapAssoc = new TH2F("hClusterMapAssoc",
286
                                "hClusterMapAssoc; x_{cluster} [px]; x_{cluster} [px]; # entries",
287
                                m_detector->nPixels().X(),
288
289
                                -0.5,
                                m_detector->nPixels().X() - 0.5,
290
                                m_detector->nPixels().Y(),
291
292
                                -0.5,
                                m_detector->nPixels().Y() - 0.5);
293

294
295
296
    hTotVsTime = new TH2F("hTotVsTime", "hTotVsTime", 64, 0, 64, 1e6, 0, 100);
    hTotVsTime->GetXaxis()->SetTitle("pixel ToT [lsb]");
    hTotVsTime->GetYaxis()->SetTitle("time [s]");
297
    hTotVsTime_high = new TH2F("hTotVsTime_high", "hTotVsTime_high", 64, 0, 64, 1e6, 0, 100);
298
299
    hTotVsTime_high->GetXaxis()->SetTitle("pixel ToT [lsb] if > high_tot_cut");
    hTotVsTime_high->GetYaxis()->SetTitle("time [s]");
300

301
    // control plots for "left tail" and "main peak" of time correlation
302
    hInPixelMap_leftTail = new TH2F("hPixelMap_leftTail",
303
                                    "in-pixel track position (left tail of time residual);in-pixel x_{track} "
304
305
306
307
308
309
310
                                    "[#mum];in-pixel y_{track} [#mum];# entries",
                                    nbins_x,
                                    -pitch_x / 2.,
                                    pitch_x / 2.,
                                    nbins_y,
                                    -pitch_y / 2.,
                                    pitch_y / 2.);
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
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
    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 =
363
        new TH1F("hTot_leftTail", "ToT (left tail of time residual);seed pixel ToT [lsb]; # events", 2 * 64, -64, 64);
364
    hTot_rightTail =
365
        new TH1F("hTot_rightTail", "ToT (left tail of time residual);seed pixel ToT [lsb]; # events", 2 * 64, -64, 64);
366
    hTot_mainPeak =
367
        new TH1F("hTot_mainPeak", "ToT (main peak of time residual);seed pixel ToT [lsb]; # events", 2 * 64, -64, 64);
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
    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);
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419

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

    gTimeCorrelationVsRow = new TGraphErrors();
    gTimeCorrelationVsRow->SetName("gTimeCorrelationVsRow");
    gTimeCorrelationVsRow->SetTitle("gTimeCorrelationVsRow");
    gTimeCorrelationVsRow->GetXaxis()->SetTitle("row");
    gTimeCorrelationVsRow->GetYaxis()->SetTitle("time correlation peak [ns]");

    // !!!!also fix these:!!!!
    int nBinsToT = hTrackCorrelationTimeVsTot_rowCorr->GetNbinsY();
    gTimeCorrelationVsTot_rowCorr = new TGraphErrors(nBinsToT);
    gTimeCorrelationVsTot_rowCorr->SetName("gTimeCorrelationVsTot");
    gTimeCorrelationVsTot_rowCorr->SetTitle("gTimeCorrelationVsTot");
    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_1px");
    gTimeCorrelationVsTot_rowCorr_1px->SetTitle("gTimeCorrelationVsTot_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_npx");
    gTimeCorrelationVsTot_rowCorr_npx->SetTitle("gTimeCorrelationVsTot_npx");
    gTimeCorrelationVsTot_rowCorr_npx->GetXaxis()->SetTitle("pixel ToT [ns] (multi-pixel clusters");
    gTimeCorrelationVsTot_rowCorr_npx->GetYaxis()->SetTitle("time correlation peak [ns]");

    LOG(INFO) << "calcCorrections = " << m_calcCorrections;

    if(m_pointwise_correction_row) {
        // Import TGraphErrors for row corection:
420
421
422
423
424
425
        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.");
426
        }
427
428

        gRowCorr = static_cast<TGraphErrors*>(file.Get(m_correctionGraph_row.c_str()));
429
        // Check if graph exists in ROOT file:
430
        if(!gRowCorr) {
431
432
433
            throw InvalidValueError(
                m_config, "correction_graph_row", "Graph doesn't exist in ROOT file. Use full/path/to/graph.");
        }
434
435
436
437
438
    } else {
        LOG(STATUS) << "----> NO POINTWISE ROW CORRECTION!!!";
    }
    if(m_pointwise_correction_timewalk) {
        // Import TGraphErrors for timewalk corection:
439
440
441
442
443
444
445
        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.");
        }
446
447

        gTimeWalkCorr = static_cast<TGraphErrors*>(file.Get(m_correctionGraph_timewalk.c_str()));
448
        // Check if graph exists in ROOT file:
449
        if(!gTimeWalkCorr) {
450
451
            throw InvalidValueError(
                m_config, "correction_graph_timewalk", "Graph doesn't exist in ROOT file. Use full/path/to/graph.");
452
453
454
455
456
457
        }
    } else {
        LOG(STATUS) << "----> NO POINTWISE TIMEWALK CORRECTION!!!";
    }
}

458
StatusCode AnalysisTimingATLASpix::run(std::shared_ptr<Clipboard> clipboard) {
459
460

    // Get the telescope tracks from the clipboard
461
    auto tracks = clipboard->getData<Track>();
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    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:
503
        auto event = clipboard->getEvent();
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523

        // 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
524
        auto clusters = clipboard->getData<Cluster>(m_detector->name());
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
        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++;

541
                    if(m_config.has("cluster_charge_cut") && cluster->charge() > m_clusterChargeCut) {
542
                        LOG(DEBUG) << " - track discarded due to clusterChargeCut";
543
544
                        continue;
                    }
545
                    tracks_afterClusterChargeCut++;
546

547
                    if(m_config.has("cluster_size_cut") && cluster->size() > m_clusterSizeCut) {
548
549
550
551
552
553
554
                        LOG(DEBUG) << " - track discarded due to clusterSizeCut";
                        continue;
                    }
                    tracks_afterClusterSizeCut++;

                    double timeDiff = track->timestamp() - cluster->timestamp();
                    hTrackCorrelationTimeAssoc->Fill(timeDiff);
555
556
                    hTrackCorrelationTimeAssocVsTime->Fill(static_cast<double>(Units::convert(cluster->timestamp(), "s")),
                                                           static_cast<double>(Units::convert(timeDiff, "us")));
557

558
                    hTrackCorrelationTimeVsTot->Fill(timeDiff, cluster->getSeedPixel()->raw());
559
560
                    hTrackCorrelationTimeVsCol->Fill(timeDiff, cluster->getSeedPixel()->column());
                    hTrackCorrelationTimeVsRow->Fill(timeDiff, cluster->getSeedPixel()->row());
561

562
                    // single-pixel and multi-pixel clusters:
563
                    if(cluster->size() == 1) {
564
                        hTrackCorrelationTimeVsTot_1px->Fill(timeDiff, cluster->getSeedPixel()->raw());
565
                        hTrackCorrelationTimeVsRow_1px->Fill(timeDiff, cluster->getSeedPixel()->row());
566
                    } else {
567
                        hTrackCorrelationTimeVsTot_npx->Fill(timeDiff, cluster->getSeedPixel()->raw());
568
                        hTrackCorrelationTimeVsRow_npx->Fill(timeDiff, cluster->getSeedPixel()->row());
569
570
571
572
573
574
575
576
577
                    }

                    // 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);
578
579
                    if(cluster->charge() > m_highChargeCut && cluster->size() == 1) {
                        hHitMapAssoc_inPixel_highCharge->Fill(xmod, ymod);
580
                    }
581
                    hPixelTrackCorrelationTimeMap->Fill(xmod, ymod, timeDiff);
582
583

                    // 2D histograms: --> fill for all pixels from cluster
584
                    for(auto& pixel : cluster->pixels()) {
585
586
587
588
589
590

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

591
                        hClusterSizeVsTot_Assoc->Fill(static_cast<double>(cluster->size()), pixel->raw());
592
                        hHitMapAssoc->Fill(pixel->column(), pixel->row());
593
                        hTotVsTime->Fill(pixel->raw(), static_cast<double>(Units::convert(pixel->timestamp(), "s")));
594
595
596
                        if(pixel->raw() > m_highTotCut) {
                            hHitMapAssoc_highCharge->Fill(pixel->column(), pixel->row());
                            hTotVsTime_high->Fill(pixel->raw(),
597
598
599
600
                                                  static_cast<double>(Units::convert(pixel->timestamp(), "s")));
                        }
                    }
                    hClusterMapAssoc->Fill(cluster->column(), cluster->row());
601
602
                    if(track->timestamp() - cluster->timestamp() < m_leftTailCut) {
                        hInPixelMap_leftTail->Fill(xmod, ymod);
603
604
605
606
                    } else if(track->timestamp() - cluster->timestamp() > m_rightTailCut) {
                        hInPixelMap_rightTail->Fill(xmod, ymod);
                    } else {
                        hInPixelMap_mainPeak->Fill(xmod, ymod);
607
                    }
608
609
610
611
612
613
614
615
616
617
618
619

                    // !!! 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(),
620
                                                                 cluster->getSeedPixel()->raw());
621
622
                        if(cluster->size() == 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_1px->Fill(
623
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
624
625
626
                        }
                        if(cluster->size() > 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_npx->Fill(
627
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
628
629
630
631
632
633
634
                        }
                        //} for(auto& pixels : ...)
                    }
                    // point-wise timewalk correction on top:
                    if(m_pointwise_correction_timewalk) {
                        correctClusterTimestamp(cluster, 1); // mode=1 --> timewalk correction
                        hTrackCorrelationTime_rowAndTimeWalkCorr->Fill(track->timestamp() - cluster->timestamp());
635
                        if(cluster->getSeedPixel()->raw() < 25) {
636
637
                            hTrackCorrelationTime_rowAndTimeWalkCorr_l25->Fill(track->timestamp() - cluster->timestamp());
                        }
638
                        if(cluster->getSeedPixel()->raw() < 40) {
639
640
                            hTrackCorrelationTime_rowAndTimeWalkCorr_l40->Fill(track->timestamp() - cluster->timestamp());
                        }
641
                        if(cluster->getSeedPixel()->raw() > 40) {
642
643
644
645
646
647
                            hTrackCorrelationTime_rowAndTimeWalkCorr_g40->Fill(track->timestamp() - cluster->timestamp());
                        }

                        hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr->Fill(
                            track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->row());
                        hTrackCorrelationTimeVsTot_rowAndTimeWalkCorr->Fill(
648
                            track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
649

650
                        // control plots to investigate "left/right tail" in time correlation:
651
                        if(track->timestamp() - cluster->timestamp() < m_leftTailCut) {
652
                            hClusterMap_leftTail->Fill(cluster->column(), cluster->row());
653
                            hTot_leftTail->Fill(cluster->getSeedPixel()->raw());
654
                            hPixelTimestamp_leftTail->Fill(cluster->getSeedPixel()->timestamp());
655
                            hClusterSize_leftTail->Fill(static_cast<double>(cluster->size()));
656
657
658
659
660
661
                        } else if(track->timestamp() - cluster->timestamp() > m_rightTailCut) {
                            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()));
                        } else {
662
663
664
665
                            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()));
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
                        }
                    }
                }

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

681
void AnalysisTimingATLASpix::finalise() {
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    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);

696
697
            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 250) { // too few entries to fit
698
699
700
701
702
703
704
705
706
707
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
733
734
735
736
737
738
739
                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
740
741
                timePeak = hTemp->GetMean();
                timePeakErr = 0; // hTemp->GetStdDev();
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
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
808
809
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
                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 ///

    } // if(m_calcCorrections)

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

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

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

    // hTrackCorrelationTime_example->Write();
    gTimeCorrelationVsRow->Write();
    gTimeCorrelationVsTot_rowCorr->Write();
    gTimeCorrelationVsTot_rowCorr_1px->Write();
    gTimeCorrelationVsTot_rowCorr_npx->Write();

    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;
859
    LOG(INFO) << "after clusterTotCut:\t" << tracks_afterClusterChargeCut;
860
861
862
    LOG(INFO) << "after clusterSizeCut:\t" << tracks_afterClusterSizeCut;
}

863
void AnalysisTimingATLASpix::correctClusterTimestamp(Cluster* cluster, int mode) {
864

865
    /*
866
867
868
869
870
871
     * MODE:
     *  0 --> row correction
     *  1 --> timewalk correction
     */

    // Get the pixels on this cluster
872
873
    auto pixels = cluster->pixels();
    auto first_pixel = pixels.front();
874
    double correction = 0;
875
876

    if(mode == 0) {
877
        correction = gRowCorr->Eval(first_pixel->row());
878
    } else if(mode == 1) {
879
        correction = gTimeWalkCorr->Eval(first_pixel->raw());
880
881
882
883
884
885
886
    } else {
        LOG(ERROR) << "Mode " << mode << " does not exist!\n"
                   << "Choose\n\t0 --> row correction \n\t1-->timewalk correction";
        return;
    }

    // Initial guess for cluster timestamp:
887
    double timestamp = first_pixel->timestamp() + correction;
888
889

    // Loop over all pixels:
890
891
892
    for(auto& pixel : pixels) {
        // FIXME ugly hack
        auto px = const_cast<Pixel*>(pixel);
893
894
895
896

        if(mode == 0) {
            correction = gRowCorr->Eval(pixel->row());
        } else if(mode == 1) {
897
            correction = gTimeWalkCorr->Eval(pixel->raw());
898
899
900
901
902
        } else {
            return;
        }

        // Override pixel timestamps:
903
        px->setTimestamp(pixel->timestamp() + correction);
904
905
906
907
908
909
910
911
912

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

    cluster->setTimestamp(timestamp);
}