AnalysisTimingATLASpix.cpp 46.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
/**
 * @file
 * @brief Implementation of [AnalysisEfficiency] module
 * @copyright Copyright (c) 2017 CERN and the Allpix Squared authors.
 * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
 * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
 * Intergovernmental Organization or submit itself to any jurisdiction.
 */

10
#include "AnalysisTimingATLASpix.h"
11
12
13
14
15
16
17
18
19
20

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

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

using namespace corryvreckan;

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

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

27
28
    using namespace ROOT::Math;
    m_detector = detector;
29
30
31
32
33
34
35
36
    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();
    }
37
    m_chi2ndofCut = m_config.get<double>("chi2ndof_cut", 3.);
38
    m_timeCutFrameEdge = m_config.get<double>("time_cut_frameedge", static_cast<double>(Units::convert(20, "ns")));
39
40
41
42
43
44
45

    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");
    }
46
    m_highTotCut = m_config.get<int>("high_tot_cut", 40);
47
    m_highChargeCut = m_config.get<double>("high_charge_cut", static_cast<double>(m_highTotCut));
48
    m_leftTailCut = m_config.get<double>("left_tail_cut", static_cast<double>(Units::convert(-10, "ns")));
49
    m_rightTailCut = m_config.get<double>("right_tail_cut", static_cast<double>(Units::convert(10, "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
    m_inpixelBinSize = m_config.get<double>("inpixel_bin_size", Units::get<double>(1.0, "um"));

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

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

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

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

99
100
101
102
    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]");
103
    hTrackCorrelationTimeAssocVsTime->GetZaxis()->SetTitle("# events");
104

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

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

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

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

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

    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)");
143
    hTrackCorrelationTime_example->GetYaxis()->SetTitle("# events");
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180

    // 2D histograms:
    // column dependence
    name = "hTrackCorrelationTimeVsCol";
    hTrackCorrelationTimeVsCol =
        new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().X(), 0, m_detector->nPixels().X());
    hTrackCorrelationTimeVsCol->GetYaxis()->SetTitle("pixel column");
    hTrackCorrelationTimeVsCol->GetXaxis()->SetTitle("track time stamp - seed pixel time stamp [ns]");
    // row dependence
    name = "hTrackCorrelationTimeVsRow";
    hTrackCorrelationTimeVsRow =
        new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
    hTrackCorrelationTimeVsRow->GetYaxis()->SetTitle("pixel row");
    hTrackCorrelationTimeVsRow->GetXaxis()->SetTitle("track time stamp - seed pixel time stamp [ns]");
    name = "hTrackCorrelationTimeVsRow_1px";
    hTrackCorrelationTimeVsRow_1px =
        new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
    hTrackCorrelationTimeVsRow_1px->GetYaxis()->SetTitle("pixel row");
    hTrackCorrelationTimeVsRow_1px->GetXaxis()->SetTitle(
        "track time stamp - seed pixel time stamp [ns] (single-pixel clusters)");
    name = "hTrackCorrelationTimeVsRow_npx";
    hTrackCorrelationTimeVsRow_npx =
        new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
    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";
    hTrackCorrelationTimeVsRow_rowCorr =
        new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
    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);
181
182
    hTrackCorrelationTimeVsTot->GetYaxis()->SetTitle("seed pixel ToT [ns]");
    hTrackCorrelationTimeVsTot->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
183
184
185
186

    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)");
187
    hTrackCorrelationTimeVsTot_1px->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
188
189
190
191

    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)");
192
    hTrackCorrelationTimeVsTot_npx->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
193
194
195
196

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

    // timewalk after row correction
    name = "hTrackCorrelationTimeVsTot_rowCorr";
    hTrackCorrelationTimeVsTot_rowCorr = new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, 512, 0, 512);
202
    hTrackCorrelationTimeVsTot_rowCorr->GetYaxis()->SetTitle("seed pixel ToT [ns]");
203
204
205
206
    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);
207
208
    hTrackCorrelationTimeVsTot_rowCorr_1px->GetYaxis()->SetTitle("seed pixel ToT [ns] (single-pixel clusters)");
    hTrackCorrelationTimeVsTot_rowCorr_1px->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
209
210
211

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

    // final plots with both row and timewalk correction:
    name = "hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr";
    hTrackCorrelationTimeVsRow_rowAndTimeWalkCorr =
        new TH2F(name.c_str(), name.c_str(), 20000, -5000, 5000, m_detector->nPixels().Y(), 0, m_detector->nPixels().Y());
    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);
224
225
    hTrackCorrelationTimeVsTot_rowAndTimeWalkCorr->GetYaxis()->SetTitle("seed pixel ToT [ns]");
    hTrackCorrelationTimeVsTot_rowAndTimeWalkCorr->GetXaxis()->SetTitle("ts_{track} - ts_{cluster} [ns]");
226

227
228
229
230
231
232
233
    // 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.");
    }

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

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

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

300
    // control plots for "left tail" and "main peak" of time correlation
301
    hInPixelMap_leftTail = new TH2F("hPixelMap_leftTail",
302
                                    "in-pixel track position (left tail of time residual);in-pixel x_{track} "
303
304
305
306
307
308
309
                                    "[#mum];in-pixel y_{track} [#mum];# entries",
                                    nbins_x,
                                    -pitch_x / 2.,
                                    pitch_x / 2.,
                                    nbins_y,
                                    -pitch_y / 2.,
                                    pitch_y / 2.);
310
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
    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,
                 m_detector->nPixels().X(),
                 m_detector->nPixels().Y(),
                 0,
                 m_detector->nPixels().Y());
    hClusterMap_rightTail =
        new TH2F("hClusterMap_rightTail",
                 "hClusterMap (right tail of time residual); x_{cluster} [px]; x_{cluster} [px]; # entries",
                 m_detector->nPixels().X(),
                 0,
                 m_detector->nPixels().X(),
                 m_detector->nPixels().Y(),
                 0,
                 m_detector->nPixels().Y());
    hClusterMap_mainPeak =
        new TH2F("hClusterMap_mainPeak",
                 "hClusterMap (main peak of time residual); x_{cluster} [px]; x_{cluster} [px]; # entries",
                 m_detector->nPixels().X(),
                 0,
                 m_detector->nPixels().X(),
                 m_detector->nPixels().Y(),
                 0,
                 m_detector->nPixels().Y());
    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 =
362
        new TH1F("hTot_leftTail", "ToT (left tail of time residual);seed pixel ToT [lsb]; # events", 2 * 64, -64, 64);
363
    hTot_rightTail =
364
        new TH1F("hTot_rightTail", "ToT (left tail of time residual);seed pixel ToT [lsb]; # events", 2 * 64, -64, 64);
365
    hTot_mainPeak =
366
        new TH1F("hTot_mainPeak", "ToT (main peak of time residual);seed pixel ToT [lsb]; # events", 2 * 64, -64, 64);
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
    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);
382
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

    // /////////////////////////////////////////// //
    // 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:
419
420
421
422
423
424
        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.");
425
        }
426
427

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

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

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

    // Get the telescope tracks from the clipboard
460
    auto tracks = clipboard->getData<Track>();
461
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
    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:
502
        auto event = clipboard->getEvent();
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522

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

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

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

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

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

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

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

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

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

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

                    // !!! 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(),
619
                                                                 cluster->getSeedPixel()->raw());
620
621
                        if(cluster->size() == 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_1px->Fill(
622
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
623
624
625
                        }
                        if(cluster->size() > 1) {
                            hTrackCorrelationTimeVsTot_rowCorr_npx->Fill(
626
                                track->timestamp() - cluster->getSeedPixel()->timestamp(), cluster->getSeedPixel()->raw());
627
628
629
630
631
632
633
                        }
                        //} 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());
634
                        if(cluster->getSeedPixel()->raw() < 25) {
635
636
                            hTrackCorrelationTime_rowAndTimeWalkCorr_l25->Fill(track->timestamp() - cluster->timestamp());
                        }
637
                        if(cluster->getSeedPixel()->raw() < 40) {
638
639
                            hTrackCorrelationTime_rowAndTimeWalkCorr_l40->Fill(track->timestamp() - cluster->timestamp());
                        }
640
                        if(cluster->getSeedPixel()->raw() > 40) {
641
642
643
644
645
646
                            hTrackCorrelationTime_rowAndTimeWalkCorr_g40->Fill(track->timestamp() - cluster->timestamp());
                        }

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

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

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

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

695
696
            // if(hTemp->GetEntries() < 500) { // too few entries to fit
            if(hTemp->GetEntries() < 250) { // too few entries to fit
697
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
                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
739
740
                timePeak = hTemp->GetMean();
                timePeakErr = 0; // hTemp->GetStdDev();
741
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
                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;
858
    LOG(INFO) << "after clusterTotCut:\t" << tracks_afterClusterChargeCut;
859
860
861
    LOG(INFO) << "after clusterSizeCut:\t" << tracks_afterClusterSizeCut;
}

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

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

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

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

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

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

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

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

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

    cluster->setTimestamp(timestamp);
}