EventLoaderEUDAQ2.cpp 29.9 KB
Newer Older
1
2
/**
 * @file
3
 * @brief Implementation of module EventLoaderEUDAQ2
4
5
 *
 * @copyright Copyright (c) 2019-2020 CERN and the Corryvreckan authors.
6
7
8
9
10
11
 * 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.
 */

#include "EventLoaderEUDAQ2.h"
12
#include "eudaq/FileReader.hh"
13
14
15

using namespace corryvreckan;

16
EventLoaderEUDAQ2::EventLoaderEUDAQ2(Configuration& config, std::shared_ptr<Detector> detector)
17
    : Module(config, detector), detector_(detector) {
18

19
20
21
    // Backwards compatibility: also allow get_tag_vectors to be used for get_tag_profiles
    config_.setAlias("get_tag_profiles", "get_tag_vectors", true);

22
    config_.setDefault<bool>("get_time_residuals", false);
23
24
    config_.setDefault<bool>("get_tag_histograms", false);
    config_.setDefault<bool>("get_tag_profiles", false);
25
26
27
28
29
    config_.setDefault<bool>("ignore_bore", true);
    config_.setDefault<double>("skip_time", 0.);
    config_.setDefault<int>("buffer_depth", 0);
    config_.setDefault<int>("shift_triggers", 0);
    config_.setDefault<bool>("inclusive", true);
30
    config_.setDefault<std::string>("eudaq_loglevel", "ERROR");
31

32
33
    filename_ = config_.getPath("file_name", true);
    get_time_residuals_ = config_.get<bool>("get_time_residuals");
34

35
36
    get_tag_histograms_ = config_.get<bool>("get_tag_histograms");
    get_tag_profiles_ = config_.get<bool>("get_tag_profiles");
37
38
39
40
41
42
    ignore_bore_ = config_.get<bool>("ignore_bore");
    skip_time_ = config_.get<double>("skip_time");
    adjust_event_times_ = config_.getMatrix<std::string>("adjust_event_times", {});
    buffer_depth_ = config_.get<int>("buffer_depth");
    shift_triggers_ = config_.get<int>("shift_triggers");
    inclusive_ = config_.get<bool>("inclusive");
43
44
45

    // Set EUDAQ log level to desired value:
    EUDAQ_LOG_LEVEL(config_.get<std::string>("eudaq_loglevel"));
46
    LOG(INFO) << "Setting EUDAQ2 log level to \"" << config_.get<std::string>("eudaq_loglevel") << "\"";
47

48
49
50
    // Prepare EUDAQ2 config object
    eudaq::Configuration cfg;

51
52
    // Forward all settings to EUDAQ
    // WARNING: the EUDAQ Configuration class is not very flexible and e.g. booleans have to be passed as 1 and 0.
53
    auto configs = config_.getAll();
54
55
56
57
58
    for(const auto& key : configs) {
        LOG(DEBUG) << "Forwarding key \"" << key.first << " = " << key.second << "\" to EUDAQ converter";
        cfg.Set(key.first, key.second);
    }

59
60
61
    // Provide the calibration file specified in the detector geometry:
    // NOTE: This should go last to allow overwriting the calibration_file key in the module config with absolute path
    auto calibration_file =
62
        config_.has("calibration_file") ? config_.getPath("calibration_file", true) : detector_->calibrationFile();
63
64
65
66
67
    if(!calibration_file.empty()) {
        LOG(DEBUG) << "Forwarding detector calibration file: " << calibration_file;
        cfg.Set("calibration_file", calibration_file);
    }

Simon Spannagel's avatar
Simon Spannagel committed
68
    // Converting the newly built configuration to a shared pointer of a const configuration object
69
70
71
    // Unfortunbately EUDAQ does not provide appropriate member functions for their configuration class to avoid this dance
    const eudaq::Configuration eu_cfg = cfg;
    eudaq_config_ = std::make_shared<const eudaq::Configuration>(eu_cfg);
72
}
73

74
void EventLoaderEUDAQ2::initialize() {
75

76
    // Declare histograms
77
    std::string title = ";EUDAQ event start time[ms];# entries";
78
    hEudaqEventStart = new TH1D("eudaqEventStart", title.c_str(), 3e6, 0, 3e3);
79

80
81
82
    title = ";EUDAQ event start time[s];# entries";
    hEudaqEventStart_long = new TH1D("eudaqEventStart_long", title.c_str(), 3e6, 0, 3e3);

83
    title = "Corryvreckan event start times (placed on clipboard); Corryvreckan event start time [ms];# entries";
84
    hClipboardEventStart = new TH1D("clipboardEventStart", title.c_str(), 3e6, 0, 3e3);
85

86
    title = "Corryvreckan event start times (placed on clipboard); Corryvreckan event start time [s];# entries";
87
88
    hClipboardEventStart_long = new TH1D("clipboardEventStart_long", title.c_str(), 3e6, 0, 3e3);

89
    title = "Corryvreckan event end times (placed on clipboard); Corryvreckan event end time [ms];# entries";
90
    hClipboardEventEnd = new TH1D("clipboardEventEnd", title.c_str(), 3e6, 0, 3e3);
91

92
93
    title = "Corryvreckan event end times (on clipboard); Corryvreckan event duration [ms];# entries";
    hClipboardEventDuration = new TH1D("clipboardEventDuration", title.c_str(), 3e6, 0, 3e3);
94

Lennart Huth's avatar
Lennart Huth committed
95
96
97
    hTriggersPerEvent =
        new TH1D("hTriggersPerEvent", "Number of triggers per event;triggers per event;entries", 20, -0.5, 19.5);

98
    title = " # events per corry event; number of events from " + detector_->getName() + " per corry event;# entries";
Lennart Huth's avatar
Lennart Huth committed
99
    hEudaqeventsPerCorry = new TH1D("hEudaqeventsPerCorryEvent", title.c_str(), 50, -.5, 49.5);
Lennart Huth's avatar
Lennart Huth committed
100
101
    title = "number of hits in corry frame vs number of eudaq frames;eudaq frames;# hits";
    hHitsVersusEUDAQ2Frames = new TH2D("hHitsVersusEUDAQ2Frames", title.c_str(), 15, -.5, 14.5, 200, -0.5, 199.5);
102
    // Create the following histograms only when detector is not auxiliary:
103
    if(!detector_->isAuxiliary()) {
104
105
106
        title = "hitmap;column;row;# events";
        hitmap = new TH2F("hitmap",
                          title.c_str(),
107
                          detector_->nPixels().X(),
108
                          -0.5,
109
110
                          detector_->nPixels().X() - 0.5,
                          detector_->nPixels().Y(),
111
                          -0.5,
112
                          detector_->nPixels().Y() - 0.5);
113
114
115
116

        title = "rawValues; column; row; raw values";
        hRawValuesMap = new TProfile2D("hRawValuesMap",
                                       title.c_str(),
117
                                       detector_->nPixels().X(),
118
                                       -0.5,
119
120
                                       detector_->nPixels().X() - 0.5,
                                       detector_->nPixels().Y(),
121
                                       -0.5,
122
                                       detector_->nPixels().Y() - 0.5);
123
124
125
126
127
128
129
130

        title = ";hit time [ms];# events";
        hPixelTimes = new TH1F("hPixelTimes", title.c_str(), 3e6, 0, 3e3);

        title = ";hit timestamp [s];# events";
        hPixelTimes_long = new TH1F("hPixelTimes_long", title.c_str(), 3e6, 0, 3e3);

        title = ";pixel raw values;# events";
131
        hPixelRawValues = new TH1F("hPixelRawValues", title.c_str(), 1024, -0.5, 1023.5);
132

133
        title = "Pixel Multiplicity per EUDAQ Event;# pixels;# events";
134
        hPixelMultiplicityPerEudaqEvent = new TH1F("hPixelMultiplicityPerEudaqEvent", title.c_str(), 1000, -0.5, 999.5);
135

136
        title = "Pixel Multiplicity per Corry Event;# pixels;# events";
137
        hPixelMultiplicityPerCorryEvent = new TH1F("hPixelMultiplicityPerCorryEvent", title.c_str(), 1000, -0.5, 999.5);
138

139
        if(get_time_residuals_) {
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
            hPixelTimeEventBeginResidual =
                new TH1F("hPixelTimeEventBeginResidual",
                         "hPixelTimeEventBeginResidual;pixel_ts - clipboard event begin [us]; # entries",
                         2.1e5,
                         -10,
                         200);

            hPixelTimeEventBeginResidual_wide =
                new TH1F("hPixelTimeEventBeginResidual_wide",
                         "hPixelTimeEventBeginResidual_wide;pixel_ts - clipboard event begin [us]; # entries",
                         1e5,
                         -5000,
                         5000);
            hPixelTimeEventBeginResidualOverTime =
                new TH2F("hPixelTimeEventBeginResidualOverTime",
                         "hPixelTimeEventBeginResidualOverTime; pixel time [s];pixel_ts - clipboard event begin [us]",
                         3e3,
                         0,
                         3e3,
                         2.1e4,
                         -10,
                         200);
162
            std::string histTitle = "hPixelTriggerTimeResidualOverTime_0;time [s];pixel_ts - trigger_ts [us];# entries";
163
164
165
            hPixelTriggerTimeResidualOverTime =
                new TH2D("hPixelTriggerTimeResidualOverTime_0", histTitle.c_str(), 3e3, 0, 3e3, 1e4, -50, 50);
        }
166
    }
167

168
169
    // open the input file with the eudaq reader
    try {
170
        reader_ = eudaq::Factory<eudaq::FileReader>::MakeUnique(eudaq::str2hash("native"), filename_);
171
    } catch(...) {
172
        LOG(ERROR) << "eudaq::FileReader could not read the input file ' " << filename_
173
                   << " '. Please verify that the path and file name are correct.";
174
        throw InvalidValueError(config_, "file_path", "Parsing error!");
175
    }
176

177
    // Check if all elements of m_adjust_event_times have a valid size of 3, if not throw error.
178
    for(auto& shift_times : adjust_event_times_) {
179
180
        if(shift_times.size() != 3) {
            throw InvalidValueError(
181
                config_,
182
                "adjust_event_times",
183
                "Parameter needs 3 values per row: [\"event type\", shift event start, shift event end]");
184
        }
185
    }
186
187
}

188
189
std::shared_ptr<eudaq::StandardEvent> EventLoaderEUDAQ2::get_next_sorted_std_event() {

190
    // Refill the event buffer if necessary:
191
    while(static_cast<int>(sorted_events_.size()) < buffer_depth_) {
192
193
194
        LOG(DEBUG) << "Filling buffer with new event.";
        // fill buffer with new std event:
        auto new_event = get_next_std_event();
195
        sorted_events_.push(new_event);
196
197
    }

198
199
200
    // get first element of queue and erase it
    auto stdevt = sorted_events_.top();
    sorted_events_.pop();
201
    return stdevt;
202
203
204
}

std::shared_ptr<eudaq::StandardEvent> EventLoaderEUDAQ2::get_next_std_event() {
205
206
207
    auto stdevt = eudaq::StandardEvent::MakeShared();
    bool decoding_failed = true;
    do {
208
209
210
        // Create new StandardEvent
        stdevt = eudaq::StandardEvent::MakeShared();

211
212
213
        // Check if we need a new full event or if we still have some in the cache:
        if(events_.empty()) {
            LOG(TRACE) << "Reading new EUDAQ event from file";
214
            auto new_event = reader_->GetNextEvent();
215
216
217
218
            if(!new_event) {
                LOG(DEBUG) << "Reached EOF";
                throw EndOfFile();
            }
219
220
            // Build buffer from all sub-events:
            events_ = new_event->GetSubEvents();
221
            // The main event might also contain data, so add it to the buffer:
Simon Spannagel's avatar
Simon Spannagel committed
222
            if(events_.empty()) {
223
                events_.push_back(new_event);
Simon Spannagel's avatar
Simon Spannagel committed
224
            }
225
        }
226
        LOG(TRACE) << "Buffer contains " << events_.size() << " (sub-) events:";
227
228
        for(auto& evt : events_) {
            LOG(TRACE) << "  (sub-) event of type " << evt->GetDescription();
229
        }
230
231

        // Retrieve first and remove from buffer:
232
233
        auto event = events_.front();
        events_.erase(events_.begin());
234

235
        // If this is a Begin-of-Run event and we should ignore it, please do so:
236
        if(event->IsBORE() && ignore_bore_) {
237
238
239
240
            LOG(DEBUG) << "Found EUDAQ2 BORE event, ignoring it";
            continue;
        }

241
242
        decoding_failed = !eudaq::StdEventConverter::Convert(event, stdevt, eudaq_config_);

243
        // Read and store tag information:
244
        if(get_tag_histograms_ || get_tag_profiles_) {
245
            retrieve_event_tags(stdevt);
246
        }
247

248
        LOG(DEBUG) << event->GetDescription() << ": EventConverter returned " << (decoding_failed ? "false" : "true");
249
250
251
    } while(decoding_failed);
    return stdevt;
}
252

253
void EventLoaderEUDAQ2::retrieve_event_tags(const eudaq::EventSPC evt) {
254
255
256
257
258
259
260
261
    auto tags = evt->GetTags();

    for(auto tag_pair : tags) {
        // Trying to convert tag value to double:
        try {
            double value = std::stod(tag_pair.second);

            // Check if histogram exists already, if not: create it
262
            if(get_tag_histograms_) {
263
264
265
266
267
268
269
                if(tagHist.find(tag_pair.first) == tagHist.end()) {
                    std::string name = "tagHist_" + tag_pair.first;
                    std::string title = tag_pair.first + ";tag value;# entries";
                    tagHist[tag_pair.first] = new TH1D(name.c_str(), title.c_str(), 1024, -512.5, 511.5);
                }
                tagHist[tag_pair.first]->Fill(value);
            }
270
            if(get_tag_profiles_) {
271
272
273
274
275
276
                if(tagProfile.find(tag_pair.first) == tagProfile.end()) {
                    std::string name = "tagProfile_" + tag_pair.first;
                    std::string title = "tag_" + tag_pair.first + ";event / 1000;tag value";
                    tagProfile[tag_pair.first] = new TProfile(name.c_str(), title.c_str(), 4e5, 0, 100);
                }
                tagProfile[tag_pair.first]->Fill(evt->GetEventN() / 1000, value, 1);
277
278
279
280
281
            }
        } catch(std::exception& e) {
        }
    }
}
282
Event::Position EventLoaderEUDAQ2::is_within_event(const std::shared_ptr<Clipboard>& clipboard,
283
                                                   std::shared_ptr<eudaq::StandardEvent> evt) const {
284

285
    // Check if this event has timestamps available:
286
    if(evt->GetTimeBegin() == 0 && evt->GetTimeEnd() == 0) {
287
        LOG(DEBUG) << evt->GetDescription() << ": Event has no timestamp, comparing trigger IDs";
288

289
        // If there is no event defined yet, there is little we can do:
290
        if(!clipboard->isEventDefined()) {
291
            LOG(DEBUG) << "No Corryvreckan event defined - cannot define without timestamps.";
292
            return Event::Position::UNKNOWN;
293
        }
294

295
        // Get position of this time frame with respect to the defined event:
296
        auto trigger_position = clipboard->getEvent()->getTriggerPosition(evt->GetTriggerN());
297
298
299
300
301
302
303
        if(trigger_position == Event::Position::BEFORE) {
            LOG(DEBUG) << "Trigger ID " << evt->GetTriggerN() << " before triggers registered in Corryvreckan event";
        } else if(trigger_position == Event::Position::AFTER) {
            LOG(DEBUG) << "Trigger ID " << evt->GetTriggerN() << " after triggers registered in Corryvreckan event";
        } else if(trigger_position == Event::Position::UNKNOWN) {
            LOG(DEBUG) << "Trigger ID " << evt->GetTriggerN() << " within Corryvreckan event range but not registered";
        } else {
304
            // Redefine EUDAQ2 event begin and end to trigger timestamp (both were zero):
305
306
            evt->SetTimeBegin(static_cast<uint64_t>(clipboard->getEvent()->getTriggerTime(evt->GetTriggerN()) * 1000));
            evt->SetTimeEnd(static_cast<uint64_t>(clipboard->getEvent()->getTriggerTime(evt->GetTriggerN()) * 1000));
307
308
309
            LOG(DEBUG) << "Trigger ID " << evt->GetTriggerN() << " found in Corryvreckan event";
        }
        return trigger_position;
310
    }
311

312
    // Read time from EUDAQ2 event and convert from picoseconds to nanoseconds:
313
314
    double event_start = static_cast<double>(evt->GetTimeBegin()) / 1000 + detector_->timeOffset();
    double event_end = static_cast<double>(evt->GetTimeEnd()) / 1000 + detector_->timeOffset();
315
316
    // Store the original position of the event before adjusting its length:
    double event_timestamp = event_start;
317
318
    LOG(DEBUG) << "event_start = " << Units::display(event_start, "us")
               << ", event_end = " << Units::display(event_end, "us");
Jens Kroeger's avatar
Jens Kroeger committed
319

320
    // If adjustment of event start/end is required:
321
322
    const auto it = std::find_if(adjust_event_times_.begin(),
                                 adjust_event_times_.end(),
323
                                 [evt](const std::vector<std::string>& x) { return x.front() == evt->GetDescription(); });
324

325
    if(it != adjust_event_times_.end()) {
326
327
328
329
330
331
332
333
334
335
        event_start += corryvreckan::from_string<double>(it->at(1));
        event_end += corryvreckan::from_string<double>(it->at(2));
        LOG(DEBUG) << "Adjusting " << it->at(0) << " event_start by "
                   << Units::display(corryvreckan::from_string<double>(it->at(1)), {"us", "ns"}) << ", event_end by "
                   << Units::display(corryvreckan::from_string<double>(it->at(2)), {"us", "ns"});
        LOG(DEBUG) << "Adjusted event: " << Units::display(event_start, {"us", "ns"}) << " - "
                   << Units::display(event_end, {"us", "ns"}) << ", length "
                   << Units::display(event_end - event_start, {"us", "ns"});
    }

336
    // Skip if later start is requested:
337
    if(event_start < skip_time_) {
338
        LOG(DEBUG) << "Event start before requested skip time: " << Units::display(event_start, {"us", "ns"}) << " < "
339
                   << Units::display(skip_time_, {"us", "ns"});
340
        return Event::Position::BEFORE;
341
342
    }

Simon Spannagel's avatar
Simon Spannagel committed
343
    // Check if an event is defined or if we need to create it:
344
    if(!clipboard->isEventDefined()) {
345
        LOG(DEBUG) << "Defining Corryvreckan event: " << Units::display(event_start, {"us", "ns"}) << " - "
Simon Spannagel's avatar
Simon Spannagel committed
346
347
                   << Units::display(event_end, {"us", "ns"}) << ", length "
                   << Units::display(event_end - event_start, {"us", "ns"});
348
        clipboard->putEvent(std::make_shared<Event>(event_start, event_end));
349
350
351
352
353
        hClipboardEventStart->Fill(static_cast<double>(Units::convert(event_start, "ms")));
        hClipboardEventStart_long->Fill(static_cast<double>(Units::convert(event_start, "s")));
        hClipboardEventEnd->Fill(static_cast<double>(Units::convert(event_end, "ms")));
        hClipboardEventDuration->Fill(
            static_cast<double>(Units::convert(clipboard->getEvent()->end() - clipboard->getEvent()->start(), "ms")));
354
    } else {
355
356
357
358
        LOG(DEBUG) << "Corryvreckan event found on clipboard: "
                   << Units::display(clipboard->getEvent()->start(), {"us", "ns"}) << " - "
                   << Units::display(clipboard->getEvent()->end(), {"us", "ns"})
                   << ", length: " << Units::display(clipboard->getEvent()->duration(), {"us", "ns"});
359
    }
360

361
    // Get position of this time frame with respect to the defined event:
362
    auto position = clipboard->getEvent()->getFramePosition(event_start, event_end, inclusive_);
363
    if(position == Event::Position::BEFORE) {
364
        LOG(DEBUG) << "Event start before Corryvreckan event: " << Units::display(event_start, {"us", "ns"}) << " < "
365
                   << Units::display(clipboard->getEvent()->start(), {"us", "ns"});
366
    } else if(position == Event::Position::AFTER) {
367
        LOG(DEBUG) << "Event end after Corryvreckan event: " << Units::display(event_end, {"us", "ns"}) << " > "
368
                   << Units::display(clipboard->getEvent()->end(), {"us", "ns"});
369
    } else {
370
371
        // check if event has valid trigger ID (flag = 0x10):
        if(evt->IsFlagTrigger()) {
372
            // Potentially shift the trigger IDs if requested
373
            auto trigger_id = static_cast<uint32_t>(static_cast<int>(evt->GetTriggerN()) + shift_triggers_);
374
            // Store potential trigger numbers, assign to center of event:
375
            clipboard->getEvent()->addTrigger(trigger_id, event_timestamp);
376
            LOG(DEBUG) << "Stored trigger ID " << trigger_id << " at " << Units::display(event_timestamp, {"us", "ns"});
377
        }
378
    }
379
380

    return position;
381
}
382

383
PixelVector EventLoaderEUDAQ2::get_pixel_data(std::shared_ptr<eudaq::StandardEvent> evt, int plane_id) const {
384

385
    PixelVector pixels;
386

387
388
389
390
    // No plane found:
    if(plane_id < 0) {
        return pixels;
    }
391

392
393
394
    auto plane = evt->GetPlane(static_cast<size_t>(plane_id));
    // Concatenate plane name according to naming convention: sensor_type + "_" + int
    auto plane_name = plane.Sensor() + "_" + std::to_string(plane.ID());
395
    auto detector_name = detector_->getName();
396
397
398
399
400
401
402
    // Convert to lower case before string comparison to avoid errors by the user:
    std::transform(plane_name.begin(), plane_name.end(), plane_name.begin(), ::tolower);
    std::transform(detector_name.begin(), detector_name.end(), detector_name.begin(), ::tolower);
    LOG(TRACE) << plane_name << " (ID " << plane_id << ") with " << plane.HitPixels() << " pixel hits";

    // Loop over all hits and add to pixels vector:
    for(unsigned int i = 0; i < plane.HitPixels(); i++) {
403

404
405
406
407
        auto col = static_cast<int>(plane.GetX(i));
        auto row = static_cast<int>(plane.GetY(i));
        auto raw = static_cast<int>(plane.GetPixel(i)); // generic pixel raw value (could be ToT, ADC, ...)

408
409
        double ts;
        if(plane.GetTimestamp(i) == 0) {
410
411
412
413
            // If the plane timestamp is not defined, we place all pixels in the center of the EUDAQ2 event.
            // Note: If the EUDAQ2 event has zero timestamps (time begin/end == 0), then the event times are
            // redefined to time begin/end == trigger timestamp in get_trigger_position (see above), i.e. in
            // that case, all pixel timestamp will be set to the corresponding trigger timestamp.
414
            ts = static_cast<double>(evt->GetTimeBegin() + evt->GetTimeEnd()) / 2 / 1000 + detector_->timeOffset();
415
        } else {
416
            ts = static_cast<double>(plane.GetTimestamp(i)) / 1000 + detector_->timeOffset();
417
418
        }

419
        if(col >= detector_->nPixels().X() || row >= detector_->nPixels().Y()) {
Lennart Huth's avatar
Lennart Huth committed
420
            LOG(WARNING) << "Pixel address " << col << ", " << row << " is outside of pixel matrix with size "
421
                         << detector_->nPixels();
422
        }
423

424
        if(detector_->masked(col, row)) {
425
            LOG(TRACE) << "Masking pixel (col, row) = (" << col << ", " << row << ")";
426
            continue;
427
        } else {
428
429
            LOG(TRACE) << "Storing (col, row, timestamp) = (" << col << ", " << row << ", "
                       << Units::display(ts, {"ns", "us", "ms"}) << ") from EUDAQ2 event data";
430
        }
Jens Kroeger's avatar
Jens Kroeger committed
431

432
        // when calibration is not available, set charge = raw
433
        auto pixel = std::make_shared<Pixel>(detector_->getName(), col, row, raw, raw, ts);
434

435
436
437
438
        hitmap->Fill(col, row);
        hPixelTimes->Fill(static_cast<double>(Units::convert(ts, "ms")));
        hPixelTimes_long->Fill(static_cast<double>(Units::convert(ts, "s")));
        hPixelRawValues->Fill(raw);
439
        hRawValuesMap->Fill(col, row, raw);
440

441
        pixels.push_back(pixel);
442
    }
443
    hPixelMultiplicityPerEudaqEvent->Fill(static_cast<int>(pixels.size()));
444
    LOG(DEBUG) << detector_->getName() << ": Plane contains " << pixels.size() << " pixels";
445

446
    return pixels;
447
448
}

449
bool EventLoaderEUDAQ2::filter_detectors(std::shared_ptr<eudaq::StandardEvent> evt, int& plane_id) const {
450
451
452
453
454
455
456
457
    // Check if the detector type matches the currently processed detector type:
    auto detector_type = evt->GetDetectorType();
    std::transform(detector_type.begin(), detector_type.end(), detector_type.begin(), ::tolower);
    // Fall back to parsing the description if not set:
    if(detector_type.empty()) {
        LOG(TRACE) << "Using fallback comparison with EUDAQ2 event description";
        auto description = evt->GetDescription();
        std::transform(description.begin(), description.end(), description.begin(), ::tolower);
458
459
        if(description.find(detector_->getType()) == std::string::npos) {
            LOG(DEBUG) << "Ignoring event because description doesn't match type " << detector_->getType() << ": "
460
461
462
                       << description;
            return false;
        }
463
    } else if(detector_type != detector_->getType()) {
464
465
466
467
468
        LOG(DEBUG) << "Ignoring event because detector type doesn't match: " << detector_type;
        return false;
    }

    // To the best of our knowledge, this is the detector we are looking for.
469
    LOG(DEBUG) << "Found matching event for detector type " << detector_->getType();
470
471
472
473
474
475
476
477
478
479
    if(evt->NumPlanes() == 0) {
        return true;
    }

    // Check if we can identify the detector itself among the planes:
    for(size_t i_plane = 0; i_plane < evt->NumPlanes(); i_plane++) {
        auto plane = evt->GetPlane(i_plane);

        // Concatenate plane name according to naming convention: sensor_type + "_" + int
        auto plane_name = plane.Sensor() + "_" + std::to_string(plane.ID());
480
        auto detector_name = detector_->getName();
481
482
483
484
485
        // Convert to lower case before string comparison to avoid errors by the user:
        std::transform(plane_name.begin(), plane_name.end(), plane_name.begin(), ::tolower);
        std::transform(detector_name.begin(), detector_name.end(), detector_name.begin(), ::tolower);

        if(detector_name == plane_name) {
486
            plane_id = static_cast<int>(i_plane);
487
            LOG(DEBUG) << "Found matching plane in event for detector " << detector_->getName();
488
            return true;
Lennart Huth's avatar
Lennart Huth committed
489
490
        } else {
            LOG(DEBUG) << "plane " << plane_name << "does not match " << detector_name;
491
492
493
494
        }
    }

    // Detector not found among planes of this event
495
    LOG(DEBUG) << "Ignoring event because no matching plane could be found for detector " << detector_->getName();
496
    return false;
497
498
}

499
StatusCode EventLoaderEUDAQ2::run(const std::shared_ptr<Clipboard>& clipboard) {
Lennart Huth's avatar
Lennart Huth committed
500
    size_t num_eudaq_events_per_corry = 0;
501

502
    PixelVector pixels;
503

504
    Event::Position current_position = Event::Position::UNKNOWN;
505
    while(1) {
506
507
508
        // Retrieve next event from file/buffer:
        if(!event_) {
            try {
509
                if(buffer_depth_ == 0) {
510
511
512
513
514
515
                    // simply get next decoded EUDAQ StandardEvent from buffer
                    event_ = get_next_std_event();
                } else {
                    // get next decoded EUDAQ StandardEvent from timesorted buffer
                    event_ = get_next_sorted_std_event();
                }
516
517
            } catch(EndOfFile&) {
                return StatusCode::EndRun;
518
519
#ifdef STDEVENTCONVERTER_EXCEPTIONS_
                // Only evaluate if defined by EUDAQ2
520
521
522
            } catch(eudaq::DataInvalid& e) {
                LOG(ERROR) << "EUDAQ2 reports invalid data: " << e.what() << std::endl << "Ending run.";
                return StatusCode::EndRun;
523
#endif
524
            }
525
        }
526

527
528
529
        // Filter out "wrong" detectors and store plane ID if found:
        int plane_id = -1;
        if(!filter_detectors(event_, plane_id)) {
530
531
532
            event_.reset();
            continue;
        }
533
534
        // Check if this event is within the currently defined Corryvreckan event:
        current_position = is_within_event(clipboard, event_);
535

536
        if(current_position == Event::Position::DURING) {
Lennart Huth's avatar
Lennart Huth committed
537
            num_eudaq_events_per_corry++;
538
            LOG(DEBUG) << "Is within current Corryvreckan event, storing data";
539
            // Store data on the clipboard
540
            auto new_pixels = get_pixel_data(event_, plane_id);
541
            hits_ += new_pixels.size();
542
            pixels.insert(pixels.end(), new_pixels.begin(), new_pixels.end());
543
        }
544

545
546
        // If this event was after the current event or if we have not enough information, stop reading:
        if(current_position == Event::Position::AFTER || current_position == Event::Position::UNKNOWN) {
547
548
            break;
        }
549

550
        // Do not fill if current_position == Event::Position::AFTER to avoid double-counting!
551
        // Converting EUDAQ2 picoseconds into Corryvreckan nanoseconds:
552
        hEudaqEventStart->Fill(static_cast<double>(event_->GetTimeBegin()) / 1e9);       // here convert from ps to ms
553
        hEudaqEventStart_long->Fill(static_cast<double>(event_->GetTimeBegin()) / 1e12); // here convert from ps to seconds
554

Simon Spannagel's avatar
Simon Spannagel committed
555
        // Reset this shared event pointer to get a new event from the stack:
556
        event_.reset();
557
    }
558

559
    auto event = clipboard->getEvent();
Simon Spannagel's avatar
Simon Spannagel committed
560
    hTriggersPerEvent->Fill(static_cast<double>(event->triggerList().size()));
561
562
563
    if(event->triggerList().size() == 0) {
        LOG(WARNING) << "Triggers on clipboard event: " << event->triggerList().size();
    }
Simon Spannagel's avatar
Simon Spannagel committed
564
    LOG(DEBUG) << "Triggers on clipboard event: " << event->triggerList().size();
565
    for(auto& trigger : event->triggerList()) {
566
        LOG(DEBUG) << "\t ID: " << trigger.first << ", time: " << std::fixed << Units::display(trigger.second, "us");
567
568
    }

569
    // histogram only exists for non-auxiliary detectors:
570
    if(!detector_->isAuxiliary()) {
571
        hPixelMultiplicityPerCorryEvent->Fill(static_cast<int>(pixels.size()));
572
573
    }

574
    // Loop over pixels for plotting
575
    if(get_time_residuals_) {
576
        for(auto& pixel : pixels) {
577
578
579
580
581
582
583
584
585
586
587
588
589
            hPixelTimeEventBeginResidual->Fill(
                static_cast<double>(Units::convert(pixel->timestamp() - event->start(), "us")));
            hPixelTimeEventBeginResidual_wide->Fill(
                static_cast<double>(Units::convert(pixel->timestamp() - event->start(), "us")));
            hPixelTimeEventBeginResidualOverTime->Fill(
                static_cast<double>(Units::convert(pixel->timestamp(), "s")),
                static_cast<double>(Units::convert(pixel->timestamp() - event->start(), "us")));

            size_t iTrigger = 0;
            for(auto& trigger : event->triggerList()) {
                // check if histogram exists already, if not: create it
                if(hPixelTriggerTimeResidual.find(iTrigger) == hPixelTriggerTimeResidual.end()) {
                    std::string histName = "hPixelTriggerTimeResidual_" + to_string(iTrigger);
590
                    std::string histTitle = histName + ";pixel_ts - trigger_ts [us];# entries";
591
592
593
594
595
                    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(
596
                    static_cast<double>(Units::convert(pixel->timestamp() - trigger.second, "us")));
597
598
599
600
601
602
                if(iTrigger == 0) { // fill only for 0th trigger
                    hPixelTriggerTimeResidualOverTime->Fill(
                        static_cast<double>(Units::convert(pixel->timestamp(), "s")),
                        static_cast<double>(Units::convert(pixel->timestamp() - trigger.second, "us")));
                }
                iTrigger++;
603
            }
604
        }
605
606
    }

Lennart Huth's avatar
Lennart Huth committed
607
608
609
    // Store the full event data on the clipboard
    hEudaqeventsPerCorry->Fill(static_cast<double>(num_eudaq_events_per_corry));
    hHitsVersusEUDAQ2Frames->Fill(static_cast<double>(num_eudaq_events_per_corry), static_cast<double>(pixels.size()));
610
    clipboard->putData(pixels, detector_->getName());
611

612
    LOG(DEBUG) << "Finished Corryvreckan event";
613
    return StatusCode::Success;
614
}
615

616
void EventLoaderEUDAQ2::finalize(const std::shared_ptr<ReadonlyClipboard>&) {
617

618
    LOG(INFO) << "Found " << hits_ << " hits in the data.";
619
}