EventLoaderATLASpix.cpp 30.5 KB
Newer Older
1
#include "EventLoaderATLASpix.h"
2
3
4
5
6
#include <regex>

using namespace corryvreckan;
using namespace std;

7
8
EventLoaderATLASpix::EventLoaderATLASpix(Configuration config, std::shared_ptr<Detector> detector)
    : Module(std::move(config), detector), m_detector(detector) {
9

10
    m_inputDirectory = m_config.getPath("input_directory");
11
    m_clockCycle = m_config.get<double>("clock_cycle", Units::get<double>(6.25, "ns"));
12
13

    // Allow reading of legacy data format using the Karlsruhe readout system:
14
    m_legacyFormat = m_config.get<bool>("legacy_format", false);
15

16
17
18
    // m_clkdivendM = m_config.get<int>("clkdivend", 0.) + 1;
    m_clkdivend2M = m_config.get<int>("clkdivend2", 0.) + 1;

19
20
    m_highToTCut = m_config.get<int>("high_tot_cut", 40);

21
22
23
24
    if(m_config.has("calibration_file")) {
        m_calibrationFile = m_config.getPath("calibration_file");
    }

25
26
    // ts1Range = 0x800 * m_clkdivendM;
    ts2Range = 0x40 * m_clkdivend2M;
27
}
28

29
uint32_t EventLoaderATLASpix::gray_decode(uint32_t gray) {
30
31
32
33
34
    uint32_t bin = gray;
    while(gray >>= 1) {
        bin ^= gray;
    }
    return bin;
35
}
36

37
void EventLoaderATLASpix::initialise() {
38

39
    uint32_t datain;
40

41
42
    m_detectorBusy = false;

Tomas Vanat's avatar
Tomas Vanat committed
43
44
    // File structure is RunX/data.bin
    // Assume that the ATLASpix is the DUT (if running this algorithm)
45

46
    // Open the root directory
47
    DIR* directory = opendir(m_inputDirectory.c_str());
48
    if(directory == nullptr) {
49
        LOG(ERROR) << "Directory " << m_inputDirectory << " does not exist";
50
51
52
53
54
        return;
    }
    dirent* entry;

    // Read the entries in the folder
55
    while((entry = readdir(directory))) {
56
        // Check for the data file
57
        string filename = m_inputDirectory + "/" + entry->d_name;
58
        if(filename.find("data.bin") != string::npos) {
59
60
61
62
63
            m_filename = filename;
        }
    }

    // If no data was loaded, give a warning
64
    if(m_filename.length() == 0) {
65
        LOG(WARNING) << "No data file was found for ATLASpix in " << m_inputDirectory;
66
    } else {
67
        LOG(STATUS) << "Opened data file for ATLASpix: (dbg)" << m_filename;
68
    }
69

70
71
72
73
74
    // Open the binary data file for later
    m_file.open(m_filename.c_str(), ios::in | ios::binary);
    LOG(DEBUG) << "Opening file " << m_filename;

    // fast forward file to T0 event
75
    old_fpga_ts = 0;
76
    while(1) {
77
        m_file.read(reinterpret_cast<char*>(&datain), 4);
78
79
80
81
        if(m_file.eof()) {
            m_file.clear();
            m_file.seekg(ios::beg);
            LOG(WARNING) << "No T0 event was found in file " << m_filename
82
                         << ". Rewinding the file to the beginning and loading all events.";
83
84
85
86
87
            break;
        } else if((datain & 0xFF000000) == 0x70000000) {
            LOG(STATUS) << "Found T0 event at position " << m_file.tellg() << ". Skipping all data before this event.";
            oldpos = m_file.tellg();
            unsigned long ts3 = datain & 0x00FFFFFF;
88
89
            old_fpga_ts = (static_cast<unsigned long long>(ts3));
            int checkpos = static_cast<int>(oldpos) - 8;
90
            while(checkpos >= 0) {
91
92
                std::streampos tmppos = checkpos;
                m_file.seekg(tmppos);
93
                m_file.read(reinterpret_cast<char*>(&datain), 4);
94
95
96
                unsigned int message_type = (datain >> 24);
                // TS2
                if(message_type == 0b00100000) {
97
                    old_fpga_ts = ((static_cast<unsigned long long>(datain & 0x00FFFFFF)) << 24) | ts3;
98
                    LOG(DEBUG) << "Set old_fpga_ts to " << old_fpga_ts;
99
100
101
102
103
104
105
106
107
                    break;
                }
                // TS3
                else if(message_type == 0b01100000) {
                    if(ts3 != (datain & 0x00FFFFFF)) {
                        LOG(WARNING) << "Last FPGA timestamp " << (datain & 0x00FFFFFF) << " does not match to T0 event "
                                     << ts3 << ". Some timestamps at the begining might be corrupted.";
                    }
                }
108
                checkpos = static_cast<int>(tmppos) - 4;
109
110
111
112
113
            }
            m_file.seekg(oldpos);
            break;
        }
    }
114

115
    // Make histograms for debugging
116
    hHitMap = new TH2F("hitMap",
117
                       "hitMap; pixel column; pixel row; # events",
118
                       m_detector->nPixels().X(),
119
                       0,
120
121
                       m_detector->nPixels().X(),
                       m_detector->nPixels().Y(),
122
                       0,
123
                       m_detector->nPixels().Y());
Jens Kroeger's avatar
Jens Kroeger committed
124

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
    hHitMap_highTot = new TH2F("hitMap_highTot",
                               "hitMap_hithTot; pixel column; pixel row; # events",
                               m_detector->nPixels().X(),
                               0,
                               m_detector->nPixels().X(),
                               m_detector->nPixels().Y(),
                               0,
                               m_detector->nPixels().Y());

    hHitMap_totWeighted = new TProfile2D("hHitMap_totWeighted",
                                         "hHitMap_totWeighted; pixel column; pixel row; # events",
                                         m_detector->nPixels().X(),
                                         0,
                                         m_detector->nPixels().X(),
                                         m_detector->nPixels().Y(),
                                         0,
                                         m_detector->nPixels().Y(),
                                         0,
                                         100);

    hPixelToT = new TH1F("pixelToT", "pixelToT; pixel ToT in TS2 clock cycles; # events", 64, 0, 64);
    hPixelToT_beforeCorrection = new TH1F(
        "pixelToT_beforeCorrection", "pixelToT_beforeCorrection; pixel ToT in TS2 clock cycles; # events", 2 * 64, -64, 64);
    hPixelCharge = new TH1F("pixelCharge", "pixelCharge; pixel charge [e]; # events", 100, 0, 100);
    hPixelToA = new TH1F("pixelToA", "pixelToA; pixel ToA [ns]; # events", 100, 0, 100);
150
    hPixelMultiplicity = new TH1F("pixelMultiplicity", "Pixel Multiplicity; # pixels; # events", 200, 0, 200);
151
152
    hPixelTimes = new TH1F("hPixelTimes", "pixelTimes; hit timestamp [ms]; # events", 3e6, 0, 3e3);
    hPixelTimes_long = new TH1F("hPixelTimes_long", "pixelTimes_long; hit timestamp [s]; # events", 3e6, 0, 3e3);
153
154
155
156
157
158

    hPixelTS1 = new TH1F("pixelTS1", "pixelTS1; pixel TS1 [lsb]; # events", 2050, 0, 2050);
    hPixelTS2 = new TH1F("pixelTS2", "pixelTS2; pixel TS2 [lsb]; # events", 130, 0, 130);
    hPixelTS1bits = new TH1F("pixelTS1bits", "pixelTS1bits; pixel TS1 bit [lsb->msb]; # events", 12, 0, 12);
    hPixelTS2bits = new TH1F("pixelTS2bits", "pixelTS2bits; pixel TS2 bit [lsb->msb]; # events", 8, 0, 8);

159
160
    hTriggersPerEvent = new TH1D("hTriggersPerEvent", "hTriggersPerEvent;triggers per event;entries", 20, 0, 20);

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
    // low ToT:
    hPixelTS1_lowToT = new TH1F("pixelTS1_lowToT", "pixelTS1_lowToT; pixel TS1 [lsb]; # events", 2050, 0, 2050);
    hPixelTS2_lowToT = new TH1F("pixelTS2_lowToT", "pixelTS2_lowToT; pixel TS2 [lsb]; # events", 130, 0, 130);
    hPixelTS1bits_lowToT =
        new TH1F("pixelTS1bits_lowToT", "pixelTS1bits_lowToT; pixel TS1 bit [lsb->msb]; # events", 12, 0, 12);
    hPixelTS2bits_lowToT =
        new TH1F("pixelTS2bits_lowToT", "pixelTS2bits_lowToT; pixel TS2 bit [lsb->msb]; # events", 8, 0, 8);

    // high ToT:
    hPixelTS1_highToT = new TH1F("pixelTS1_highToT", "pixelTS1_highToT; pixel TS1 [lsb]; # events", 2050, 0, 2050);
    hPixelTS2_highToT = new TH1F("pixelTS2_highToT", "pixelTS2_highToT; pixel TS2 [lsb]; # events", 130, 0, 130);
    hPixelTS1bits_highToT =
        new TH1F("pixelTS1bits_highToT", "pixelTS1bits_highToT; pixel TS1 bit [lsb->msb]; # events", 12, 0, 12);
    hPixelTS2bits_highToT =
        new TH1F("pixelTS2bits_highToT", "pixelTS2bits_highToT; pixel TS2 bit [lsb->msb]; # events", 8, 0, 8);
176

177
178
179
180
181
    hPixelTimeEventBeginResidual = new TH1F("hPixelTimeEventBeginResidual",
                                            "hPixelTimeEventBeginResidual;pixel_ts - clipboard event begin [us]; # entries",
                                            2.1e5,
                                            -10,
                                            200);
182
183
184
185
186
187
    hPixelTimeEventBeginResidual_wide =
        new TH1F("hPixelTimeEventBeginResidual_wide",
                 "hPixelTimeEventBeginResidual_wide;pixel_ts - clipboard event begin [us]; # entries",
                 1e5,
                 -5000,
                 5000);
188
189
190
191
192
193
194
195
196
197
    hPixelTimeEventBeginResidualOverTime =
        new TH2F("hPixelTimeEventBeginResidualOverTime",
                 "hPixelTimeEventBeginResidualOverTime; pixel time [s];pixel_ts - clipboard event begin [us]",
                 3e3,
                 0,
                 3e3,
                 2.1e4,
                 -10,
                 200);

198
    std::string histTitle = "hPixelTriggerTimeResidualOverTime_0;time [s];pixel_ts - trigger_ts [us];# entries";
199
200
201
    hPixelTriggerTimeResidualOverTime =
        new TH2D("hPixelTriggerTimeResidualOverTime_0", histTitle.c_str(), 3e3, 0, 3e3, 1e4, -50, 50);

202
    // Read calibration:
203
    m_calibrationFactors.resize(static_cast<size_t>(m_detector->nPixels().X() * m_detector->nPixels().Y()), 1.0);
204
205
206
207
208
209
210
211
212
    if(!m_calibrationFile.empty()) {
        std::ifstream calibration(m_calibrationFile);
        std::string line;
        std::getline(calibration, line);

        int col, row;
        double calibfactor;
        while(getline(calibration, line)) {
            std::istringstream(line) >> col >> row >> calibfactor;
213
            m_calibrationFactors.at(static_cast<size_t>(row * 25 + col)) = calibfactor;
214
215
        }
        calibration.close();
216
217
    }

Jens Kroeger's avatar
Jens Kroeger committed
218
    LOG(INFO) << "Using clock cycle length of " << m_clockCycle << " ns." << std::endl;
219

220
221
    m_oldtoa = 0;
    m_overflowcounter = 0;
222
223
}

224
StatusCode EventLoaderATLASpix::run(std::shared_ptr<Clipboard> clipboard) {
225

226
    // Check if event frame is defined:
Simon Spannagel's avatar
Simon Spannagel committed
227
    auto event = clipboard->get_event();
228

229
230
231
    // If have reached the end of file, close it and exit program running
    if(m_file.eof()) {
        m_file.close();
232
        LOG(DEBUG) << "Returning <Failure> status, ATLASPix data file reached the end.";
233
        return StatusCode::Failure;
234
235
    }

Simon Spannagel's avatar
Simon Spannagel committed
236
237
    double start_time = event->start();
    double end_time = event->end();
238
    bool busy_at_start = m_detectorBusy;
239
240

    // Read pixel data
241
    Pixels* pixels = (m_legacyFormat ? read_legacy_data(start_time, end_time) : read_caribou_data(start_time, end_time));
242

243
244
245
246
247
    if(busy_at_start || m_detectorBusy) {
        LOG(DEBUG) << "Returning <DeadTime> status, ATLASPix is BUSY.";
        return StatusCode::DeadTime;
    }

248
249
    // Here fill some histograms for data quality monitoring:
    auto nTriggers = event->triggerList().size();
250
    LOG(DEBUG) << "nTriggers = " << nTriggers;
251
252
    hTriggersPerEvent->Fill(static_cast<double>(nTriggers));

253
254
    for(auto px : (*pixels)) {
        hHitMap->Fill(px->column(), px->row());
255
        if(px->raw() > m_highToTCut) {
256
257
            hHitMap_highTot->Fill(px->column(), px->row());
        }
258
259
        hHitMap_totWeighted->Fill(px->column(), px->row(), px->raw());
        hPixelToT->Fill(px->raw());
260
        hPixelCharge->Fill(px->charge());
261
        hPixelToA->Fill(px->timestamp());
262

263
        hPixelTimeEventBeginResidual->Fill(static_cast<double>(Units::convert(px->timestamp() - start_time, "us")));
264
        hPixelTimeEventBeginResidual_wide->Fill(static_cast<double>(Units::convert(px->timestamp() - start_time, "us")));
265
266
        hPixelTimeEventBeginResidualOverTime->Fill(static_cast<double>(Units::convert(px->timestamp(), "s")),
                                                   static_cast<double>(Units::convert(px->timestamp() - start_time, "us")));
267
268
269
270
        size_t iTrigger = 0;
        for(auto& trigger : event->triggerList()) {
            // check if histogram exists already, if not: create it
            if(hPixelTriggerTimeResidual.find(iTrigger) == hPixelTriggerTimeResidual.end()) {
271
                std::string histName = "hPixelTriggerTimeResidual_" + to_string(iTrigger);
272
                std::string histTitle = histName + ";pixel_ts - trigger_ts [us];# entries";
273
274
275
276
277
                hPixelTriggerTimeResidual[iTrigger] = new TH1D(histName.c_str(), histTitle.c_str(), 2e5, -100, 100);
            }
            // use iTrigger, not trigger ID (=trigger.first) (which is unique and continuously incrementing over the runtime)
            hPixelTriggerTimeResidual[iTrigger]->Fill(
                static_cast<double>(Units::convert(px->timestamp() - trigger.second, "us")));
278
279
280
281
282
            if(iTrigger == 0) { // fill only for 0th trigger
                hPixelTriggerTimeResidualOverTime->Fill(
                    static_cast<double>(Units::convert(px->timestamp(), "s")),
                    static_cast<double>(Units::convert(px->timestamp() - trigger.second, "us")));
            }
283
284
            iTrigger++;
        }
285
286
        hPixelTimes->Fill(static_cast<double>(Units::convert(px->timestamp(), "ms")));
        hPixelTimes_long->Fill(static_cast<double>(Units::convert(px->timestamp(), "s")));
287
288
    }

289
    hPixelMultiplicity->Fill(static_cast<double>(pixels->size()));
290

291
292
    // Put the data on the clipboard
    if(!pixels->empty()) {
293
        clipboard->put(m_detector->name(), "pixels", reinterpret_cast<Objects*>(pixels));
294
    } else {
Jens Kroeger's avatar
Jens Kroeger committed
295
        delete pixels;
296
        LOG(DEBUG) << "Returning <NoData> status, no hits found.";
297
        return StatusCode::NoData;
298
299
300
    }

    // Return value telling analysis to keep running
301
    LOG(DEBUG) << "Returning <Success> status, hits found in this event window.";
302
    return StatusCode::Success;
303
304
}

305
Pixels* EventLoaderATLASpix::read_caribou_data(double start_time, double end_time) {
306
307
    LOG(DEBUG) << "Searching for events in interval from " << Units::display(start_time, {"s", "us", "ns"}) << " to "
               << Units::display(end_time, {"s", "us", "ns"}) << ", file read position " << m_file.tellg()
308
               << ", old_fpga_ts = " << old_fpga_ts << ".";
309

310
311
312
313
    // Pixel container
    Pixels* pixels = new Pixels();

    // Read file and load data
314
    uint32_t datain;
315
    long long ts_diff; // tmp
316
317

    // Initialize all to 0 for a case that hit data come before timestamp/trigger data
318
319
320
    long long hit_ts = 0;            // 64bit value of a hit timestamp combined from readout and pixel hit timestamp
    unsigned long long atp_ts = 0;   // 16bit value of ATLASpix readout timestamp
    unsigned long long trig_cnt = 0; // 32bit value of trigger counter (in FPGA)
321
322
323
324
325
326
327
    unsigned long long readout_ts =
        old_readout_ts;                       // 64bit value of a readout timestamp combined from FPGA and ATp timestamp
    unsigned long long fpga_ts = old_fpga_ts; // 64bit value of FPGA readout timestamp
    unsigned long long fpga_ts1 = 0;          // tmp [63:48] of FPGA readout timestamp
    unsigned long long fpga_ts2 = 0;          // tmp [47:24] of FPGA readout timestamp
    unsigned long long fpga_ts3 = 0;          // tmp [23:0] of FPGA readout timestamp
    unsigned long long fpga_tsx = 0;          // tmp for FPGA readout timestamp
328
329
    bool new_ts1 = true;
    bool new_ts2 = true;
330
331
332
    bool window_end = false;
    bool keep_pointer_stored = false;
    bool keep_reading = true;
333
334

    // Repeat until input EOF:
335
    while(keep_reading) {
336
        // Read next 4-byte data from file
337
        m_file.read(reinterpret_cast<char*>(&datain), 4);
338
339
340
341
        if(m_file.eof()) {
            LOG(DEBUG) << "EOF...";
            break;
        }
342

343
344
345
346
347
348
        // LOG(DEBUG) << "Found " << (datain & 0x80000000 ? "pixel data" : "header information");

        // Check if current word is a pixel data:
        if(datain & 0x80000000) {
            // Structure: {1'b1, column_addr[5:0], row_addr[8:0], rise_timestamp[9:0], fall_timestamp[5:0]}
            // Extract pixel data
349
            long ts2 = gray_decode((datain)&0x003F);
350
            // long ts2 = gray_decode((datain>>6)&0x003F);
351
            // TS1 counter is by default half speed of TS2. By multiplying with 2 we make it equal.
352
353
            long ts1 = (gray_decode((datain >> 6) & 0x03FF)) << 1;
            // long ts1 = (gray_decode(((datain << 4) & 0x3F0) | ((datain >> 12) & 0xF)))<<1;
354
355
            int row = ((datain >> (6 + 10)) & 0x01FF);
            int col = ((datain >> (6 + 10 + 9)) & 0x003F);
356
357
            // long tot = 0;

358
            ts_diff = ts1 - static_cast<long long>(readout_ts & 0x07FF);
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378

            if(ts_diff > 0) {
                // Hit probably came before readout started and meanwhile an OVF of TS1 happened
                if(ts_diff > 0x01FF) {
                    ts_diff -= 0x0800;
                }
                // Hit probably came after readout started and is within range.
                // else {
                //    // OK...
                //}
            } else {
                // Hit probably came after readout started and after OVF of TS1.
                if(ts_diff < (0x01FF - 0x0800)) {
                    ts_diff += 0x0800;
                }
                // Hit probably came before readout started and is within range.
                // else {
                //    // OK...
                //}
            }
379

380
            hit_ts = static_cast<long long>(readout_ts) + ts_diff;
381

382
            // Convert the timestamp to nanoseconds:
383
            double timestamp = m_clockCycle * static_cast<double>(hit_ts) + m_detector->timingOffset();
384

385
            if(timestamp > end_time) {
386
387
                keep_pointer_stored = true;
                LOG(DEBUG) << "Skipping processing event, pixel is after event window ("
388
                           << Units::display(timestamp, {"s", "us", "ns"}) << " > "
389
                           << Units::display(end_time, {"s", "us", "ns"}) << ")";
390
                continue;
391
392
            }

393
394
395
396
397
398
            if(timestamp < start_time) {
                LOG(DEBUG) << "Skipping pixel hit, pixel is before event window ("
                           << Units::display(timestamp, {"s", "us", "ns"}) << " < "
                           << Units::display(start_time, {"s", "us", "ns"}) << ")";
                continue;
            }
399
400
            // this window still contains data in the event window, do not stop processing
            window_end = false;
401
402
403
404
405
406
407
            if(m_detectorBusy && (busy_readout_ts < readout_ts)) {
                LOG(WARNING) << "ATLASPix went BUSY between "
                             << Units::display((m_clockCycle * static_cast<double>(busy_readout_ts)), {"s", "us", "ns"})
                             << " and "
                             << Units::display((m_clockCycle * static_cast<double>(readout_ts)), {"s", "us", "ns"}) << ".";
                m_detectorBusy = false;
            }
408
            data_pixel_++;
409
            // If this pixel is masked, do not save it
410
            if(m_detector->masked(col, row)) {
411
412
413
                continue;
            }

414
            // calculate ToT only when pixel is good for storing (division is time consuming)
415
            int tot = static_cast<int>(ts2 - ((hit_ts % static_cast<long long>(64 * m_clkdivend2M)) / m_clkdivend2M));
416
            hPixelToT_beforeCorrection->Fill(tot);
417
            if(tot < 0) {
418
                tot += 64;
419
            }
Tomas Vanat's avatar
Tomas Vanat committed
420
            // convert ToT to nanoseconds
421
            // double tot_ns = tot * m_clockCycle;
422

423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
            hPixelTS1->Fill(static_cast<double>(ts1));
            hPixelTS2->Fill(static_cast<double>(ts2));
            if(tot < m_highToTCut) {
                hPixelTS1_lowToT->Fill(static_cast<double>(ts1));
                hPixelTS2_lowToT->Fill(static_cast<double>(ts2));
            } else {
                hPixelTS1_highToT->Fill(static_cast<double>(ts1));
                hPixelTS2_highToT->Fill(static_cast<double>(ts2));
            }

            // histogram each bit of the 10-bit TS1
            for(int i = 0; i < 12; i++) {
                hPixelTS1bits->Fill(static_cast<double>(i), static_cast<double>((ts1 >> i) & 0b1));
                if(tot < m_highToTCut) {
                    hPixelTS1bits_lowToT->Fill(static_cast<double>(i), static_cast<double>((ts1 >> i) & 0b1));
                } else {
                    hPixelTS1bits_highToT->Fill(static_cast<double>(i), static_cast<double>((ts1 >> i) & 0b1));
                }
            }
            // histogram each bit of the 6-bit TS2
            for(int i = 0; i < 8; i++) {
                hPixelTS2bits->Fill(static_cast<double>(i), static_cast<double>((ts2 >> i) & 0b1));
                if(tot < m_highToTCut) {
                    hPixelTS2bits_lowToT->Fill(static_cast<double>(i), static_cast<double>((ts1 >> i) & 0b1));
                } else {
                    hPixelTS2bits_highToT->Fill(static_cast<double>(i), static_cast<double>((ts1 >> i) & 0b1));
                }
            }

452
453
454
455
            LOG(TRACE) << "HIT: TS1: " << ts1 << "\t0x" << std::hex << ts1 << "\tTS2: " << ts2 << "\t0x" << std::hex << ts2
                       << "\tTS_FULL: " << hit_ts << "\t" << Units::display(timestamp, {"s", "us", "ns"})
                       << "\tTOT: " << tot; // << "\t" << Units::display(tot_ns, {"s", "us", "ns"});

456
457
            if(col >= m_detector->nPixels().X() || row >= m_detector->nPixels().Y()) {
                LOG(WARNING) << "Pixel address " << col << ", " << row << " is outside of pixel matrix.";
458
                continue;
459
460
            }

461
462
            // when calibration is not available, set charge = tot
            Pixel* pixel = new Pixel(m_detector->name(), col, row, tot, tot, timestamp);
463

464
            // FIXME: implement conversion from ToT to charge:
465
466
            // thres-->e: 1620e/0.15V, or 1080e/100mV
            // How to get from ToT to charge?
467
468
469
470
471
            //
            // Do something like:
            // charge = charge_calibration(tot);
            // pixel->setCharge(charge);
            // pixel->setCalibrated = true;
472

473
            LOG(DEBUG) << "PIXEL:\t" << *pixel;
474
            pixels->push_back(pixel);
475

476
        } else {
477
478
479
480
481
            // data is not hit information
            // if (keep_pointer_stored) then we will go through the data one more time and wee will count that next time.
            if(!keep_pointer_stored) {
                data_header_++;
            }
482

483
484
485
486
487
            // Decode the message content according to 8 MSBits
            unsigned int message_type = (datain >> 24);
            switch(message_type) {
            // Timestamp from ATLASpix [23:0]
            case 0b01000000:
488
489
490
491
492
493
494
495
496
                // the whole previous readout was behind the time window, we can stop processing and return
                if(window_end) {
                    LOG(TRACE) << "Rewinding to file pointer : " << oldpos;
                    m_file.seekg(oldpos);
                    // exit the while loop
                    keep_reading = false;
                    // exit case
                    break;
                }
497
498

                atp_ts = (datain >> 7) & 0x1FFFE;
499
                ts_diff = static_cast<long long>(atp_ts) - static_cast<long long>(fpga_ts & 0x1FFFF);
500
501
502
503
504
505
506
507
508
509

                if(ts_diff > 0) {
                    if(ts_diff > 0x10000) {
                        ts_diff -= 0x20000;
                    }
                } else {
                    if(ts_diff < -0x1000) {
                        ts_diff += 0x20000;
                    }
                }
510
                readout_ts = static_cast<unsigned long long>(static_cast<long long>(fpga_ts) + ts_diff);
511

512
513
514
515
516
517
518
519
520
521
522
                if(!keep_pointer_stored) {
                    // Store this position in the file in case we need to rewind:
                    LOG(TRACE) << "Storing file pointer position: " << m_file.tellg() << " and readout TS: " << readout_ts;
                    oldpos = m_file.tellg();
                    old_readout_ts = readout_ts;
                } else {
                    LOG(TRACE) << "File pointer position already stored for the first event out of the window. Current "
                                  "readout_ts = "
                               << readout_ts;
                }
                // If the readout time is after the window, mark it as a candidate for last readout in the window
523
                if((static_cast<double>(readout_ts) * m_clockCycle) > end_time) {
524
525
                    window_end = true;
                }
526
                break;
527

528
529
530
531
532
533
534
535
            // Trigger counter from FPGA [23:0] (1/4)
            case 0b00010000:
                trig_cnt = datain & 0x00FFFFFF;
                break;

            // Trigger counter from FPGA [31:24] and timestamp from FPGA [63:48] (2/4)
            case 0b00110000:
                trig_cnt |= (datain << 8) & 0xFF000000;
536
                fpga_ts1 = ((static_cast<unsigned long long>(datain) << 48) & 0xFFFF000000000000);
537
538
539
540
541
                new_ts1 = true;
                break;

            // Timestamp from FPGA [47:24] (3/4)
            case 0b00100000:
542
                fpga_tsx = ((static_cast<unsigned long long>(datain) << 24) & 0x0000FFFFFF000000);
543
544
                if((!new_ts1) && (fpga_tsx < fpga_ts2)) {
                    fpga_ts1 += 0x0001000000000000;
545
                    LOG(WARNING) << "Missing TS_FPGA_1, adding one";
546
547
548
549
550
551
552
553
554
555
556
557
                }
                new_ts1 = false;
                new_ts2 = true;
                fpga_ts2 = fpga_tsx;
                break;

            // Timestamp from FPGA [23:0] (4/4)
            case 0b01100000:
                m_identifiers["FPGA_TS"]++;
                fpga_tsx = ((datain)&0xFFFFFF);
                if((!new_ts2) && (fpga_tsx < fpga_ts3)) {
                    fpga_ts2 += 0x0000000001000000;
558
                    LOG(WARNING) << "Missing TS_FPGA_2, adding one";
559
560
561
562
563
                }
                new_ts2 = false;
                fpga_ts3 = fpga_tsx;
                fpga_ts = fpga_ts1 | fpga_ts2 | fpga_ts3;
                break;
564

565
566
567
            // BUSY was asserted due to FIFO_FULL + 24 LSBs of FPGA timestamp when it happened
            case 0b00000010:
                m_identifiers["BUSY_ASSERT"]++;
568
569
                busy_readout_ts = readout_ts;
                m_detectorBusy = true;
570
571
572
573
574
575
576
577
578
579
                break;

            // T0 received
            case 0b01110000:
                LOG(WARNING) << "Another T0 event was found in the data at position " << m_file.tellg();
                break;

            // Empty data - should not happen
            case 0b00000000:
                m_identifiers["EMPTY_DATA"]++;
Tomas Vanat's avatar
Tomas Vanat committed
580
                LOG(DEBUG) << "EMPTY_DATA";
581
582
583
584
585
586
587
588
                break;

            // Other options...
            default:
                // LOG(DEBUG) << "...Other";
                // Unknown message identifier
                if(message_type & 0b11110010) {
                    m_identifiers["UNKNOWN_MESSAGE"]++;
Tomas Vanat's avatar
Tomas Vanat committed
589
                    LOG(DEBUG) << "UNKNOWN_MESSAGE";
590
591
592
593
                } else {
                    // Buffer for chip data overflow (data that came after this word were lost)
                    if((message_type & 0b11110011) == 0b00000001) {
                        m_identifiers["BUFFER_OVERFLOW"]++;
Tomas Vanat's avatar
Tomas Vanat committed
594
                        LOG(DEBUG) << "BUFFER_OVERFLOW";
595
596
597
598
                    }
                    // SERDES lock established (after reset or after lock lost)
                    if((message_type & 0b11111110) == 0b00001000) {
                        m_identifiers["SERDES_LOCK_ESTABLISHED"]++;
Tomas Vanat's avatar
Tomas Vanat committed
599
                        LOG(DEBUG) << "SERDES_LOCK_ESTABLISHED";
600
601
602
603
                    }
                    // SERDES lock lost (data might be nonsense, including up to 2 previous messages)
                    else if((message_type & 0b11111110) == 0b00001100) {
                        m_identifiers["SERDES_LOCK_LOST"]++;
Tomas Vanat's avatar
Tomas Vanat committed
604
                        LOG(DEBUG) << "SERDES_LOCK_LOST";
605
606
607
608
                    }
                    // Unexpected data came from the chip or there was a checksum error.
                    else if((message_type & 0b11111110) == 0b00000100) {
                        m_identifiers["WEIRD_DATA"]++;
Tomas Vanat's avatar
Tomas Vanat committed
609
                        LOG(DEBUG) << "WEIRD_DATA";
610
611
612
613
614
615
616
617
618
619
620
621
622
                    }
                    // Unknown message identifier
                    else {
                        m_identifiers["UNKNOWN_MESSAGE"]++;
                        LOG(WARNING) << "UNKNOWN_MESSAGE";
                    }
                }
                break;
                // End case
            }
        }
    }
    LOG(DEBUG) << "Returning " << pixels->size() << " pixels";
623
624
625
    return pixels;
}

626
Pixels* EventLoaderATLASpix::read_legacy_data(double, double) {
627

628
    // Pixel container
629
630
631
    Pixels* pixels = new Pixels();

    // Read file and load data
632
633
    while(!m_file.eof()) {

634
        int col, row, tot, ts;
635
636
637
638
639
        unsigned long long int toa, TriggerDebugTS, dummy, bincounter;

        m_file >> col >> row >> ts >> tot >> dummy >> dummy >> bincounter >> TriggerDebugTS;

        // If this pixel is masked, do not save it
640
        if(m_detector->masked(col, row)) {
641
642
643
644
645
646
            continue;
        }

        // TOT
        if(tot <= (ts * 2 & 0x3F)) {
            tot = 64 + tot - (ts * 2 & 0x3F);
647
        } else {
648
649
650
651
            tot = tot - (ts * 2 & 0x3F);
        }

        // Apply calibration:
652
653
654
        double cal_tot = tot * m_calibrationFactors.at(static_cast<size_t>(row * 25 + col));
        LOG(TRACE) << "Hit " << row << "\t" << col << ": " << m_calibrationFactors.at(static_cast<size_t>(row * 25 + col))
                   << " * " << tot << " = " << cal_tot;
655
656
657
658

        ts &= 0xFF;
        ts *= 2; // atlaspix timestamp runs at 10MHz, multiply by to to get 20.

659
        if((bincounter & 0x1FF) < static_cast<unsigned long long>(ts)) {
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
            toa = ((bincounter & 0xFFFFFFFFFFFFFE00) - (1 << 9)) | (ts & 0x1FF);
        } else {
            toa = (bincounter & 0xFFFFFFFFFFFFFE00) | (ts & 0x1FF);
        }

        if(((toa + 10000) & 0xFFFFF000) < (m_oldtoa & 0xFFFFF000)) {
            m_overflowcounter++;
            LOG(DEBUG) << "Overflow detected " << m_overflowcounter << " " << (toa & 0xFFFFF000) << " "
                       << (m_oldtoa & 0xFFFFF000);
        } // Atlaspix only! Toa has overflow at 32 bits.

        toa += (0x100000000 * m_overflowcounter);
        m_oldtoa = toa & 0xFFFFFFFF;
        LOG(DEBUG) << "    " << row << "\t" << col << ": " << tot << " " << ts << " " << bincounter << " " << toa << " "
                   << (TriggerDebugTS - toa);

676
        // TriggerDebugTS *= 4096. / 5;              // runs with 200MHz, divide by 5 to scale counter value to 40MHz
677
678
        double toa_timestamp =
            4096. * 2 * static_cast<double>(toa); // runs with 20MHz, multiply by 2 to scale counter value to 40MHz
679

680
        // Convert TOA to nanoseconds:
681
        toa_timestamp /= (4096. * 0.04);
682

683
684
        // when calibration is not available, set charge = tot
        Pixel* pixel = new Pixel(m_detector->name(), col, row, tot, tot, toa_timestamp);
685
        pixel->setCharge(cal_tot);
686
        pixels->push_back(pixel);
687
688
    }

689
    return pixels;
690
691
}

692
void EventLoaderATLASpix::finalise() {
693

694
695
696
697
    LOG(INFO) << "Identifier distribution:";
    for(auto id : m_identifiers) {
        LOG(INFO) << "\t" << id.first << ": " << id.second;
    }
698
699

    LOG(INFO) << "Found " << data_pixel_ << " pixel data blocks and " << data_header_ << " header words";
700
}