diff --git a/README.md b/README.md index 730f66f6c6a969f4024bcceba75e618047ba53d1..a7def09e7008f9035ddc56aaf343ffaf0196356e 100644 --- a/README.md +++ b/README.md @@ -55,6 +55,7 @@ The following authors, in alphabetical order, have contributed to Corryvreckan: * Matthew Daniel Buckland, University of Liverpool, @mbucklan * Carsten Daniel Burgard, DESY, @cburgard +* Eric Buschmann, CERN, @ebuschma * Manuel Colocci, CERN, @mcolocci * Jens Dopke, STFC RAL, @jdopke * Jordi Duarte-Campderros, IFCA, @duarte diff --git a/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.cpp b/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.cpp index ba9d4424ebff53c2b9333c7258fed31b8e5067d8..e4282c9b5b4c80666751fe6c5e6b83d55de80bd5 100644 --- a/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.cpp +++ b/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.cpp @@ -26,6 +26,9 @@ using namespace std; EventLoaderTimepix3::EventLoaderTimepix3(Configuration& config, std::shared_ptr<Detector> detector) : Module(config, detector), m_detector(detector), m_currentEvent(0), m_prevTime(0), m_shutterOpen(false) { + config_.setDefault<size_t>("buffer_depth", 1000); + m_buffer_depth = config_.get<size_t>("buffer_depth"); + // Take input directory from global parameters m_inputDirectory = config_.getPath("input_directory"); @@ -38,6 +41,12 @@ EventLoaderTimepix3::EventLoaderTimepix3(Configuration& config, std::shared_ptr< void EventLoaderTimepix3::initialize() { + if(m_buffer_depth < 1) { + throw InvalidValueError(config_, "buffer_depth", "Buffer depth must be larger than 0."); + } else { + LOG(INFO) << "Using buffer_depth = " << m_buffer_depth; + } + // File structure is RunX/ChipID/files.dat // Open the root directory @@ -105,6 +114,9 @@ void EventLoaderTimepix3::initialize() { m_clearedHeader = false; m_syncTimeTDC = 0; m_TDCoverflowCounter = 0; + m_prevTriggerNumber = 0; + m_triggerOverflowCounter = 0; + eof_reached = false; // Sort all files by extracting the "serial number" from the file name while ignoring the timestamp: std::sort(detector_files.begin(), detector_files.end(), [](std::string a, std::string b) { @@ -325,329 +337,370 @@ void EventLoaderTimepix3::loadCalibration(std::string path, char delim, std::vec f.close(); } -// Function to load data for a given device, into the relevant container -bool EventLoaderTimepix3::loadData(const std::shared_ptr<Clipboard>& clipboard, - PixelVector& devicedata, - SpidrSignalVector& spidrData) { - +bool EventLoaderTimepix3::decodeNextWord() { std::string detectorID = m_detector->getName(); - auto event = clipboard->getEvent(); - bool extra = false; // temp + // Check if current file is at its end and move to the next: + if((*m_file_iterator)->eof()) { + m_file_iterator++; + LOG(INFO) << "Starting to read next file for " << detectorID << ": " << (*m_file_iterator).get(); + } - LOG(DEBUG) << "Loading data for device " << detectorID; - while(1) { - // Check if current file is at its end and move to the next: - if((*m_file_iterator)->eof()) { - m_file_iterator++; - LOG(INFO) << "Starting to read next file for " << detectorID << ": " << (*m_file_iterator).get(); + // Check if the last file is finished: + if(m_file_iterator == m_files.end()) { + LOG(INFO) << "EOF for all files of " << detectorID; + eof_reached = true; + return false; + } + + // Now read the data packets. + ULong64_t pixdata = 0; + + // If we can't read data anymore, jump to begin of loop: + if(!(*m_file_iterator)->read(reinterpret_cast<char*>(&pixdata), sizeof pixdata)) { + LOG(INFO) << "No more data in current file for " << detectorID << ": " << (*m_file_iterator).get(); + return true; + } + + LOG(TRACE) << "0x" << hex << pixdata << dec << " - " << pixdata; + + // Get the header (first 4 bits) and do things depending on what it is + // 0x4 is the "heartbeat" signal, 0xA and 0xB are pixel data + const UChar_t header = static_cast<UChar_t>((pixdata & 0xF000000000000000) >> 60) & 0xF; + + // Use header 0x4 to get the long timestamps (called syncTime here) + if(header == 0x4) { + LOG(TRACE) << "Found syncTime data"; + + // The 0x4 header tells us that it is part of the timestamp + // There is a second 4-bit header that says if it is the most + // or least significant part of the timestamp + const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF; + + // This is a bug fix. There appear to be errant packets with garbage data + // - source to be tracked down. + // Between the data and the header the intervening bits should all be 0, + // check if this is the case + const UChar_t intermediateBits = ((pixdata & 0x00FF000000000000) >> 48) & 0xFF; + if(intermediateBits != 0x00) { + LOG(DEBUG) << "Detector " << detectorID << ": intermediateBits error"; + return true; } - // Check if the last file is finished: - if(m_file_iterator == m_files.end()) { - LOG(INFO) << "EOF for all files of " << detectorID; - break; + // 0x4 is the least significant part of the timestamp + if(header2 == 0x4) { + // The data is shifted 16 bits to the right, then 12 to the left in + // order to match the timestamp format (net 4 right) + m_syncTime = (m_syncTime & 0xFFFFF00000000000) + ((pixdata & 0x0000FFFFFFFF0000) >> 4); } + // 0x5 is the most significant part of the timestamp + if(header2 == 0x5) { + // The data is shifted 16 bits to the right, then 44 to the left in + // order to match the timestamp format (net 28 left) + m_syncTime = (m_syncTime & 0x00000FFFFFFFFFFF) + ((pixdata & 0x00000000FFFF0000) << 28); + if(!m_clearedHeader && static_cast<double>(m_syncTime) / (4096. * 40000000.) < 6.) { + m_clearedHeader = true; + LOG(DEBUG) << detectorID << ": Cleared header"; + } + } + } - // Now read the data packets. - ULong64_t pixdata = 0; + // In data taking during 2015 there was sometimes still data left in the buffers at the start of + // a run. For that reason we keep skipping data until this "header" data has been cleared, when + // the heart beat signal starts from a low number (~few seconds max) + if(!m_clearedHeader) { + LOG(TRACE) << "Header not cleared, skipping data block."; + return true; + } - // If we can't read data anymore, jump to begin of loop: - if(!(*m_file_iterator)->read(reinterpret_cast<char*>(&pixdata), sizeof pixdata)) { - LOG(INFO) << "No more data in current file for " << detectorID << ": " << (*m_file_iterator).get(); - continue; + // Header 0x06 and 0x07 are the start and stop signals for power pulsing + if(header == 0x0) { + // These packets should only come from the DUT. Otherwise ignore and throw warning. + // (We observed these packets a few times per run in various telescope planes in the + // November 2018 test beam.) + if(!m_detector->isDUT()) { + LOG(WARNING) + /* << "Current time: " << Units::display(event->start(), {"s", "ms", "us", "ns"}) */ + << " detector " << detectorID << " " + << "header == 0x0! (indicates power pulsing.) Ignoring this."; + return true; } + // Note that the following code is probably outdated and/or not much tested + // (Estel used her private code for her power-pulsing studies.) To be fixed! - LOG(TRACE) << "0x" << hex << pixdata << dec << " - " << pixdata; + // Get the second part of the header + const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF; - // Get the header (first 4 bits) and do things depending on what it is - // 0x4 is the "heartbeat" signal, 0xA and 0xB are pixel data - const UChar_t header = static_cast<UChar_t>((pixdata & 0xF000000000000000) >> 60) & 0xF; + // New implementation of power pulsing signals from Adrian + if(header2 == 0x6) { + LOG(TRACE) << "Found power pulsing - start"; - // Use header 0x4 to get the long timestamps (called syncTime here) - if(header == 0x4) { - LOG(TRACE) << "Found syncTime data"; + // Read time stamp and convert to nanoseconds + const double timestamp = static_cast<double>((pixdata & 0x0000000FFFFFFFFF) << 12) / (4096 * 0.04); + const uint64_t controlbits = ((pixdata & 0x00F0000000000000) >> 52) & 0xF; - // The 0x4 header tells us that it is part of the timestamp - // There is a second 4-bit header that says if it is the most - // or least significant part of the timestamp - const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF; + const uint64_t powerOn = ((controlbits & 0x2) >> 1); + const uint64_t shutterClosed = ((controlbits & 0x1)); - // This is a bug fix. There appear to be errant packets with garbage data - // - source to be tracked down. - // Between the data and the header the intervening bits should all be 0, - // check if this is the case - const UChar_t intermediateBits = ((pixdata & 0x00FF000000000000) >> 48) & 0xFF; - if(intermediateBits != 0x00) { - LOG(DEBUG) << "Detector " << detectorID << ": intermediateBits error"; - continue; - } + auto powerSignal = (powerOn ? std::make_shared<SpidrSignal>("powerOn", timestamp) + : std::make_shared<SpidrSignal>("powerOff", timestamp)); + + sorted_signals_.push(powerSignal); + LOG(DEBUG) << "Power is " << (powerOn ? "on" : "off") << " power! Time: " << Units::display(timestamp, "ns"); + + LOG(TRACE) << "Shutter closed: " << hex << shutterClosed << dec; - // 0x4 is the least significant part of the timestamp - if(header2 == 0x4) { - // The data is shifted 16 bits to the right, then 12 to the left in - // order to match the timestamp format (net 4 right) - m_syncTime = (m_syncTime & 0xFFFFF00000000000) + ((pixdata & 0x0000FFFFFFFF0000) >> 4); + auto shutterSignal = (shutterClosed ? std::make_shared<SpidrSignal>("shutterClosed", timestamp) + : std::make_shared<SpidrSignal>("shutterOpen", timestamp)); + if(!shutterClosed) { + sorted_signals_.push(shutterSignal); + m_shutterOpen = true; + LOG(TRACE) << "Have opened shutter with signal " << shutterSignal->type() << " at time " + << Units::display(timestamp, "ns"); } - // 0x5 is the most significant part of the timestamp - if(header2 == 0x5) { - // The data is shifted 16 bits to the right, then 44 to the left in - // order to match the timestamp format (net 28 left) - m_syncTime = (m_syncTime & 0x00000FFFFFFFFFFF) + ((pixdata & 0x00000000FFFF0000) << 28); - if(!m_clearedHeader && static_cast<double>(m_syncTime) / (4096. * 40000000.) < 6.) { - m_clearedHeader = true; - LOG(DEBUG) << detectorID << ": Cleared header"; - } + + if(shutterClosed && m_shutterOpen) { + sorted_signals_.push(shutterSignal); + m_shutterOpen = false; + LOG(TRACE) << "Have closed shutter with signal " << shutterSignal->type() << " at time " + << Units::display(timestamp, "ns"); } + + LOG(DEBUG) << "Shutter is " << (shutterClosed ? "closed" : "open") + << ". Time: " << Units::display(timestamp, "ns"); } - // In data taking during 2015 there was sometimes still data left in the buffers at the start of - // a run. For that reason we keep skipping data until this "header" data has been cleared, when - // the heart beat signal starts from a low number (~few seconds max) - if(!m_clearedHeader) { - LOG(TRACE) << "Header not cleared, skipping data block."; - continue; + /* + // 0x6 is power on + if(header2 == 0x6){ + const double timestamp = ((pixdata & 0x0000000FFFFFFFFF) << 12 ) / (4096 * 0.04); + auto signal = std::make_shared<SpidrSignal>("powerOn",timestamp); + spidrData.push_back(signal); + LOG(DEBUG) <<"Turned on power! Time: " << Units::display(timestamp, "ns"); + } + // 0x7 is power off + if(header2 == 0x7){ + const double timestamp = ((pixdata & 0x0000000FFFFFFFFF) << 12 ) / (4096 * 0.04); + auto signal = std::make_shared<SpidrSignal>("powerOff",timestamp); + spidrData.push_back(signal); + LOG(DEBUG) <<"Turned off power! Time: " << Units::display(timestamp, "ns"); } + */ + } - // Header 0x06 and 0x07 are the start and stop signals for power pulsing - if(header == 0x0) { - // These packets should only come from the DUT. Otherwise ignore and throw warning. - // (We observed these packets a few times per run in various telescope planes in the - // November 2018 test beam.) - if(!m_detector->isDUT()) { - LOG(WARNING) << "Current time: " << Units::display(event->start(), {"s", "ms", "us", "ns"}) << " detector " - << detectorID << " " - << "header == 0x0! (indicates power pulsing.) Ignoring this."; - continue; + // Header 0x6 indicate trigger data + if(header == 0x6) { + const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF; + if(header2 == 0xF) { + unsigned int stamp = (pixdata & 0x1E0) >> 5; + long long int timestamp_raw = static_cast<long long int>(pixdata & 0xFFFFFFFFE00) >> 9; + long long int timestamp = 0; + int triggerNumber = ((pixdata & 0xFFF00000000000) >> 44); + + int intermediate = (pixdata & 0x1F); + if(intermediate != 0) { + return true; } - // Note that the following code is probably outdated and/or not much tested - // (Estel used her private code for her power-pulsing studies.) To be fixed! - - // Get the second part of the header - const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF; - - // New implementation of power pulsing signals from Adrian - if(header2 == 0x6) { - LOG(TRACE) << "Found power pulsing - start"; - - // Read time stamp and convert to nanoseconds - const double timestamp = static_cast<double>((pixdata & 0x0000000FFFFFFFFF) << 12) / (4096 * 0.04); - const uint64_t controlbits = ((pixdata & 0x00F0000000000000) >> 52) & 0xF; - - const uint64_t powerOn = ((controlbits & 0x2) >> 1); - const uint64_t shutterClosed = ((controlbits & 0x1)); - - // Stop looking at data if the signal is after the current event window - // (and rewind the file - // reader so that we start with this signal next event) - if(event->getTimestampPosition(timestamp) == Event::Position::AFTER) { - (*m_file_iterator)->seekg(-1 * static_cast<int>(sizeof(pixdata)), std::ios_base::cur); - LOG(TRACE) << "Signal has a time beyond the current event: " << Units::display(timestamp, "ns"); - break; - } - auto powerSignal = (powerOn ? std::make_shared<SpidrSignal>("powerOn", timestamp) - : std::make_shared<SpidrSignal>("powerOff", timestamp)); - spidrData.push_back(powerSignal); - LOG(DEBUG) << "Power is " << (powerOn ? "on" : "off") << " power! Time: " << Units::display(timestamp, "ns"); + if(triggerNumber < m_prevTriggerNumber) { + m_triggerOverflowCounter++; + } - LOG(TRACE) << "Shutter closed: " << hex << shutterClosed << dec; + // if jump back in time is larger than 1 sec, overflow detected... + if((m_syncTimeTDC - timestamp_raw) > 0x1312d000) { + m_TDCoverflowCounter++; + } - auto shutterSignal = (shutterClosed ? std::make_shared<SpidrSignal>("shutterClosed", timestamp) - : std::make_shared<SpidrSignal>("shutterOpen", timestamp)); - if(!shutterClosed) { - spidrData.push_back(shutterSignal); - m_shutterOpen = true; - LOG(TRACE) << "Have opened shutter with signal " << shutterSignal->type() << " at time " - << Units::display(timestamp, "ns"); - } + timestamp = timestamp_raw + (static_cast<long long int>(m_TDCoverflowCounter) << 35); - if(shutterClosed && m_shutterOpen) { - spidrData.push_back(shutterSignal); - m_shutterOpen = false; - LOG(TRACE) << "Have closed shutter with signal " << shutterSignal->type() << " at time " - << Units::display(timestamp, "ns"); - } + double triggerTime = + (static_cast<double>(timestamp) + static_cast<double>(stamp) / 12) / (8. * 0.04); // 320 MHz clock - LOG(DEBUG) << "Shutter is " << (shutterClosed ? "closed" : "open") - << ". Time: " << Units::display(timestamp, "ns"); - } + m_syncTimeTDC = timestamp_raw; - /* - // 0x6 is power on - if(header2 == 0x6){ - const double timestamp = ((pixdata & 0x0000000FFFFFFFFF) << 12 ) / (4096 * 0.04); - auto signal = std::make_shared<SpidrSignal>("powerOn",timestamp); - spidrData.push_back(signal); - LOG(DEBUG) <<"Turned on power! Time: " << Units::display(timestamp, "ns"); - } - // 0x7 is power off - if(header2 == 0x7){ - const double timestamp = ((pixdata & 0x0000000FFFFFFFFF) << 12 ) / (4096 * 0.04); - auto signal = std::make_shared<SpidrSignal>("powerOff",timestamp); - spidrData.push_back(signal); - LOG(DEBUG) <<"Turned off power! Time: " << Units::display(timestamp, "ns"); - } - */ - } + int triggerID = triggerNumber + (m_triggerOverflowCounter << 12); + m_prevTriggerNumber = triggerNumber; - // Header 0x6 indicate trigger data - if(header == 0x6) { - const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF; - if(header2 == 0xF) { - unsigned int stamp = (pixdata & 0x1E0) >> 5; - long long int timestamp_raw = static_cast<long long int>(pixdata & 0xFFFFFFFFE00) >> 9; - long long int timestamp = 0; - // int triggerNumber = ((pixdata & 0xFFF00000000000) >> 44); - int intermediate = (pixdata & 0x1F); - if(intermediate != 0) - continue; - - // if jump back in time is larger than 1 sec, overflow detected... - if((m_syncTimeTDC - timestamp_raw) > 0x1312d000) { - m_TDCoverflowCounter++; - } - m_syncTimeTDC = timestamp_raw; - timestamp = timestamp_raw + (static_cast<long long int>(m_TDCoverflowCounter) << 35); + auto triggerSignal = std::make_shared<SpidrSignal>("trigger", triggerTime, triggerID); + sorted_signals_.push(triggerSignal); + } + } - double triggerTime = - (static_cast<double>(timestamp) + static_cast<double>(stamp) / 12) / (8. * 0.04); // 320 MHz clock - auto triggerSignal = std::make_shared<SpidrSignal>("trigger", triggerTime); - spidrData.push_back(triggerSignal); - } + // Header 0xA and 0xB indicate pixel data + if(header == 0xA || header == 0xB) { + LOG(TRACE) << "Found pixel data"; + + // Decode the pixel information from the relevant bits + const UShort_t dcol = static_cast<UShort_t>((pixdata & 0x0FE0000000000000) >> 52); + const UShort_t spix = static_cast<UShort_t>((pixdata & 0x001F800000000000) >> 45); + const UShort_t pix = static_cast<UShort_t>((pixdata & 0x0000700000000000) >> 44); + const UShort_t col = static_cast<UShort_t>(dcol + pix / 4); + const UShort_t row = static_cast<UShort_t>(spix + (pix & 0x3)); + + // Check if this pixel is masked + if(m_detector->masked(col, row)) { + LOG(DEBUG) << "Detector " << detectorID << ": pixel " << col << "," << row << " masked"; + return true; } - // Header 0xA and 0xB indicate pixel data - if(header == 0xA || header == 0xB) { - LOG(TRACE) << "Found pixel data"; + // Get the rest of the data from the pixel + // const UShort_t pixno = col * 256 + row; + const UInt_t data = static_cast<UInt_t>((pixdata & 0x00000FFFFFFF0000) >> 16); + const unsigned int tot = (data & 0x00003FF0) >> 4; + const uint64_t spidrTime(pixdata & 0x000000000000FFFF); + const uint64_t ftoa(data & 0x0000000F); + const uint64_t toa((data & 0x0FFFC000) >> 14); - // Decode the pixel information from the relevant bits - const UShort_t dcol = static_cast<UShort_t>((pixdata & 0x0FE0000000000000) >> 52); - const UShort_t spix = static_cast<UShort_t>((pixdata & 0x001F800000000000) >> 45); - const UShort_t pix = static_cast<UShort_t>((pixdata & 0x0000700000000000) >> 44); - const UShort_t col = static_cast<UShort_t>(dcol + pix / 4); - const UShort_t row = static_cast<UShort_t>(spix + (pix & 0x3)); + // Calculate the timestamp. + unsigned long long int time = + (((spidrTime << 18) + (toa << 4) + (15 - ftoa)) << 8) + (m_syncTime & 0xFFFFFC0000000000); - // Check if this pixel is masked - if(m_detector->masked(col, row)) { - LOG(DEBUG) << "Detector " << detectorID << ": pixel " << col << "," << row << " masked"; - continue; - } + // Adjusting phases for double column shift + time += ((static_cast<unsigned long long int>(col) / 2 - 1) % 16) * 256; - // Get the rest of the data from the pixel - // const UShort_t pixno = col * 256 + row; - const UInt_t data = static_cast<UInt_t>((pixdata & 0x00000FFFFFFF0000) >> 16); - const unsigned int tot = (data & 0x00003FF0) >> 4; - const uint64_t spidrTime(pixdata & 0x000000000000FFFF); - const uint64_t ftoa(data & 0x0000000F); - const uint64_t toa((data & 0x0FFFC000) >> 14); - - // Calculate the timestamp. - unsigned long long int time = - (((spidrTime << 18) + (toa << 4) + (15 - ftoa)) << 8) + (m_syncTime & 0xFFFFFC0000000000); - - // Adjusting phases for double column shift - time += ((static_cast<unsigned long long int>(col) / 2 - 1) % 16) * 256; - - // The time from the pixels has a maximum value of ~26 seconds. We compare the pixel time to the "heartbeat" - // signal (which has an overflow of ~4 years) and check if the pixel time has wrapped back around to 0 - - // If the counter overflow happens before reading the new heartbeat - // while( abs(m_syncTime-time) > 0x0000020000000000 ){ - if(!extra) { - while(static_cast<long long>(m_syncTime) - static_cast<long long>(time) > 0x0000020000000000) { - time += 0x0000040000000000; - } - } else { - while(static_cast<long long>(m_prevTime) - static_cast<long long>(time) > 0x0000020000000000) { - time += 0x0000040000000000; - } - } + // The time from the pixels has a maximum value of ~26 seconds. We compare the pixel time to the "heartbeat" + // signal (which has an overflow of ~4 years) and check if the pixel time has wrapped back around to 0 - // Convert final timestamp into ns and add the timing offset (in nano seconds) from the detectors file (if any) - const double timestamp = static_cast<double>(time) / (4096. / 25.) + m_detector->timeOffset(); - auto position = event->getTimestampPosition(timestamp); + // If the counter overflow happens before reading the new heartbeat + // while( abs(m_syncTime-time) > 0x0000020000000000 ){ + while(static_cast<long long>(m_syncTime) - static_cast<long long>(time) > 0x0000020000000000) { + time += 0x0000040000000000; + } - // Ignore pixel data if it is before the "eventStart" read from the clipboard storage: - if(position == Event::Position::BEFORE) { - LOG(TRACE) << "Skipping pixel, is before event window (" << Units::display(timestamp, {"s", "us", "ns"}) - << " < " << Units::display(event->start(), {"s", "us", "ns"}) << ")"; - continue; + // Convert final timestamp into ns and add the timing offset (in nano seconds) from the detectors file (if any) + const double timestamp = static_cast<double>(time) / (4096. / 25.) + m_detector->timeOffset(); + + pixelToT_beforecalibration->Fill(static_cast<int>(tot)); + + // Apply calibration if applyCalibration is true + if(applyCalibration && m_detector->isDUT()) { + LOG(DEBUG) << "Applying calibration to DUT"; + size_t scol = static_cast<size_t>(col); + size_t srow = static_cast<size_t>(row); + float a = vtot.at(256 * srow + scol).at(2); + float b = vtot.at(256 * srow + scol).at(3); + float c = vtot.at(256 * srow + scol).at(4); + float t = vtot.at(256 * srow + scol).at(5); + + float toa_c = vtoa.at(256 * srow + scol).at(2); + float toa_t = vtoa.at(256 * srow + scol).at(3); + float toa_d = vtoa.at(256 * srow + scol).at(4); + + // Calculating calibrated tot and toa + float fvolts = (sqrt(a * a * t * t + 2 * a * b * t + 4 * a * c - 2 * a * t * static_cast<float>(tot) + b * b - + 2 * b * static_cast<float>(tot) + static_cast<float>(tot * tot)) + + a * t - b + static_cast<float>(tot)) / + (2 * a); + double fcharge = fvolts * 1e-3 * 3e-15 * 6241.509 * 1e15; // capacitance is 3 fF or 18.7 e-/mV + + /* Note 1: fvolts is the inverse to f(x) = a*x + b - c/(x-t). Note the +/- signs! */ + /* Note 2: The capacitance is actually smaller than 3 fC, more like 2.5 fC. But there is an offset when when + * using testpulses. Multiplying the voltage value with 20 [e-/mV] is a good approximation but means one is + * over estimating the input capacitance to compensate the missing information of the offset. */ + + float t_shift = toa_c / (fvolts - toa_t) + toa_d; + timeshiftPlot->Fill(static_cast<double>(Units::convert(t_shift, "ns"))); + const double ftimestamp = timestamp - t_shift; + LOG(DEBUG) << "Time shift= " << Units::display(t_shift, {"s", "ns"}); + LOG(DEBUG) << "Timestamp calibrated = " << Units::display(ftimestamp, {"s", "ns"}); + + if(col >= m_detector->nPixels().X() || row >= m_detector->nPixels().Y()) { + LOG(WARNING) << "Pixel address " << col << ", " << row << " is outside of pixel matrix."; } + // creating new pixel object with calibrated values of tot and toa + // when calibration is not available, set charge = tot + auto pixel = std::make_shared<Pixel>(detectorID, col, row, static_cast<int>(tot), tot, ftimestamp); + pixel->setCharge(fcharge); + sorted_pixels_.push(pixel); + hHitMap->Fill(col, row); + LOG(DEBUG) << "Pixel Charge = " << fcharge << "; ToT value = " << tot; + pixelToT_aftercalibration->Fill(fcharge); + } else { + LOG(DEBUG) << "Pixel hit at " << Units::display(timestamp, {"s", "ns"}); + // creating new pixel object with non-calibrated values of tot and toa + // when calibration is not available, set charge = tot + auto pixel = std::make_shared<Pixel>(detectorID, col, row, static_cast<int>(tot), tot, timestamp); + sorted_pixels_.push(pixel); + hHitMap->Fill(col, row); + } - // Stop looking at data if the pixel is after the current event window - // (and rewind the file reader so that we start with this pixel next event) - if(position == Event::Position::AFTER) { - LOG(DEBUG) << "Stopping processing event, pixel is after " - "event window (" - << Units::display(timestamp, {"s", "us", "ns"}) << " > " - << Units::display(event->end(), {"s", "us", "ns"}) << ")"; - (*m_file_iterator)->seekg(-1 * static_cast<int>(sizeof(pixdata)), std::ios_base::cur); - break; - } + m_prevTime = time; + } - // Otherwise create a new pixel object - pixelToT_beforecalibration->Fill(static_cast<int>(tot)); - - // Apply calibration if applyCalibration is true - if(applyCalibration && m_detector->isDUT()) { - LOG(DEBUG) << "Applying calibration to DUT"; - size_t scol = static_cast<size_t>(col); - size_t srow = static_cast<size_t>(row); - float a = vtot.at(256 * srow + scol).at(2); - float b = vtot.at(256 * srow + scol).at(3); - float c = vtot.at(256 * srow + scol).at(4); - float t = vtot.at(256 * srow + scol).at(5); - - float toa_c = vtoa.at(256 * srow + scol).at(2); - float toa_t = vtoa.at(256 * srow + scol).at(3); - float toa_d = vtoa.at(256 * srow + scol).at(4); - - // Calculating calibrated tot and toa - float fvolts = (sqrt(a * a * t * t + 2 * a * b * t + 4 * a * c - 2 * a * t * static_cast<float>(tot) + - b * b - 2 * b * static_cast<float>(tot) + static_cast<float>(tot * tot)) + - a * t - b + static_cast<float>(tot)) / - (2 * a); - double fcharge = fvolts * 1e-3 * 3e-15 * 6241.509 * 1e15; // capacitance is 3 fF or 18.7 e-/mV - - /* Note 1: fvolts is the inverse to f(x) = a*x + b - c/(x-t). Note the +/- signs! */ - /* Note 2: The capacitance is actually smaller than 3 fC, more like 2.5 fC. But there is an offset when when - * using testpulses. Multiplying the voltage value with 20 [e-/mV] is a good approximation but means one is - * over estimating the input capacitance to compensate the missing information of the offset. */ - - float t_shift = toa_c / (fvolts - toa_t) + toa_d; - timeshiftPlot->Fill(static_cast<double>(Units::convert(t_shift, "ns"))); - const double ftimestamp = timestamp - t_shift; - LOG(DEBUG) << "Time shift= " << Units::display(t_shift, {"s", "ns"}); - LOG(DEBUG) << "Timestamp calibrated = " << Units::display(ftimestamp, {"s", "ns"}); - - if(col >= m_detector->nPixels().X() || row >= m_detector->nPixels().Y()) { - LOG(WARNING) << "Pixel address " << col << ", " << row << " is outside of pixel matrix."; - } - // creating new pixel object with calibrated values of tot and toa - // when calibration is not available, set charge = tot - auto pixel = std::make_shared<Pixel>(detectorID, col, row, static_cast<int>(tot), tot, ftimestamp); - pixel->setCharge(fcharge); - devicedata.push_back(pixel); - hHitMap->Fill(col, row); - LOG(DEBUG) << "Pixel Charge = " << fcharge << "; ToT value = " << tot; - pixelToT_aftercalibration->Fill(fcharge); - } else { - LOG(DEBUG) << "Pixel hit at " << Units::display(timestamp, {"s", "ns"}); - // creating new pixel object with non-calibrated values of tot and toa - // when calibration is not available, set charge = tot - auto pixel = std::make_shared<Pixel>(detectorID, col, row, static_cast<int>(tot), tot, timestamp); - devicedata.push_back(pixel); - hHitMap->Fill(col, row); - } + return true; +} - m_prevTime = time; +void EventLoaderTimepix3::fillBuffer() { + // read data from file and fill timesorted buffer + while(sorted_pixels_.size() < m_buffer_depth && !eof_reached) { + // decodeNextWord returns false when EOF is reached and true otherwise + if(!decodeNextWord()) { + LOG(TRACE) << "decodeNextWord returns false: reached EOF."; + break; } } +} + +// Function to load data for a given device, into the relevant container +bool EventLoaderTimepix3::loadData(const std::shared_ptr<Clipboard>& clipboard, + PixelVector& devicedata, + SpidrSignalVector& spidrData) { + + std::string detectorID = m_detector->getName(); + auto event = clipboard->getEvent(); + + LOG(DEBUG) << "Loading data for device " << detectorID; + fillBuffer(); // Now we have data buffered into the temporary storage. We will sort this by time, and then load // the data from one event onto it. + while(!sorted_pixels_.empty()) { + auto pixel = sorted_pixels_.top(); + + auto position = event->getTimestampPosition(pixel->timestamp()); + + if(position == Event::Position::AFTER) { + LOG(DEBUG) << "Stopping processing event, pixel is after " + "event window (" + << Units::display(pixel->timestamp(), {"s", "us", "ns"}) << " > " + << Units::display(event->end(), {"s", "us", "ns"}) << ")"; + break; + } else if(position == Event::Position::BEFORE) { + LOG(TRACE) << "Skipping pixel, is before event window (" << Units::display(pixel->timestamp(), {"s", "us", "ns"}) + << " < " << Units::display(event->start(), {"s", "us", "ns"}) << ")"; + sorted_pixels_.pop(); + } else { + devicedata.push_back(pixel); + sorted_pixels_.pop(); + } + + // Refill buffer + fillBuffer(); + } + + while(!sorted_signals_.empty()) { + auto signal = sorted_signals_.top(); + + auto position = event->getTimestampPosition(signal->timestamp()); + + if(position == Event::Position::AFTER) { + LOG(DEBUG) << "Stopping processing event, signal is after " + "event window (" + << Units::display(signal->timestamp(), {"s", "us", "ns"}) << " > " + << Units::display(event->end(), {"s", "us", "ns"}) << ")"; + break; + } else if(position == Event::Position::BEFORE) { + LOG(TRACE) << "Skipping signal, is before event window (" + << Units::display(signal->timestamp(), {"s", "us", "ns"}) << " < " + << Units::display(event->start(), {"s", "us", "ns"}) << ")"; + sorted_signals_.pop(); + } else { + spidrData.push_back(signal); + sorted_signals_.pop(); + } + } + // If no data was loaded, return false if(devicedata.empty()) { return false; diff --git a/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.h b/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.h index 8d58ad15f74e9567fc4bf610f889097568a02da8..4bce827b97fa5a120b874452a32d697da7645fbd 100644 --- a/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.h +++ b/src/modules/EventLoaderTimepix3/EventLoaderTimepix3.h @@ -14,6 +14,7 @@ #include <TCanvas.h> #include <TH1F.h> #include <TH2F.h> +#include <queue> #include <stdio.h> #include "core/module/Module.hpp" #include "objects/Pixel.hpp" @@ -49,6 +50,8 @@ namespace corryvreckan { TH2F* pixelTOAParameterT; TH1F* timeshiftPlot; + bool decodeNextWord(); + void fillBuffer(); bool loadData(const std::shared_ptr<Clipboard>& clipboard, PixelVector&, SpidrSignalVector&); void loadCalibration(std::string path, char delim, std::vector<std::vector<float>>& dat); void maskPixels(std::string); @@ -67,6 +70,8 @@ namespace corryvreckan { std::vector<std::unique_ptr<std::ifstream>> m_files; std::vector<std::unique_ptr<std::ifstream>>::iterator m_file_iterator; + bool eof_reached; + size_t m_buffer_depth; unsigned long long int m_syncTime; bool m_clearedHeader; long long int m_syncTimeTDC; @@ -76,6 +81,18 @@ namespace corryvreckan { unsigned long long int m_prevTime; bool m_shutterOpen; + int m_prevTriggerNumber; + int m_triggerOverflowCounter; + + template <typename T> struct CompareTimeGreater { + bool operator()(const std::shared_ptr<T> a, const std::shared_ptr<T> b) { + return a->timestamp() > b->timestamp(); + } + }; + + std::priority_queue<std::shared_ptr<Pixel>, PixelVector, CompareTimeGreater<Pixel>> sorted_pixels_; + std::priority_queue<std::shared_ptr<SpidrSignal>, SpidrSignalVector, CompareTimeGreater<SpidrSignal>> + sorted_signals_; }; } // namespace corryvreckan #endif // TIMEPIX3EVENTLOADER_H diff --git a/src/objects/SpidrSignal.hpp b/src/objects/SpidrSignal.hpp index 829e6715b0f1f6e4bae3a915fef9a2bfb53560ba..11a2404c753ce49818400ba496c657dfa4a22c7f 100644 --- a/src/objects/SpidrSignal.hpp +++ b/src/objects/SpidrSignal.hpp @@ -22,6 +22,8 @@ namespace corryvreckan { // Constructors and destructors SpidrSignal(){}; SpidrSignal(std::string type, double timestamp) : Object(timestamp), m_type(type){}; + SpidrSignal(std::string type, double timestamp, size_t trigger) + : Object(timestamp), m_type(type), m_triggerNumber(trigger){}; /** * @brief Static member function to obtain base class for storage on the clipboard. @@ -36,13 +38,15 @@ namespace corryvreckan { // Set properties void type(std::string type) { m_type = type; } std::string type() const { return m_type; } + size_t trigger() const { return m_triggerNumber; } protected: // Member variables std::string m_type; + size_t m_triggerNumber; // ROOT I/O class definition - update version number when you change this class! - ClassDef(SpidrSignal, 3) + ClassDef(SpidrSignal, 4) }; // Vector type declaration diff --git a/testing/test_align_dut_timepix3tel_dut_atlaspix_ebeam120.conf b/testing/test_align_dut_timepix3tel_dut_atlaspix_ebeam120.conf index 80cf4cb38847a9c279a5814b0acead821cb7bba7..505f854cd2e7caec3a4d612a6820c3a2e5b39d8f 100644 --- a/testing/test_align_dut_timepix3tel_dut_atlaspix_ebeam120.conf +++ b/testing/test_align_dut_timepix3tel_dut_atlaspix_ebeam120.conf @@ -47,4 +47,4 @@ align_position_axes = "xy" # <-- ...orientation OR position alignment! max_track_chi2ndof = 3 #DATASET timepix3tel_dut_atlaspix_ebeam120 -#PASS T(341.262um,-1.543mm,105mm) R(2.41548deg,1.66765deg,-0.46845deg) +#PASS T(341.201um,-1.54307mm,105mm) R(2.44deg,1.60279deg,-0.470055deg) diff --git a/testing/test_align_telescope_timepix3tel_dut_atlaspix_ebeam120.conf b/testing/test_align_telescope_timepix3tel_dut_atlaspix_ebeam120.conf index 543b3dd3a50f958a27a8e560f257571dae4400a5..6df4524c3ec928fad5a8d26ce93a0c951aee3fbd 100644 --- a/testing/test_align_telescope_timepix3tel_dut_atlaspix_ebeam120.conf +++ b/testing/test_align_telescope_timepix3tel_dut_atlaspix_ebeam120.conf @@ -34,4 +34,4 @@ max_track_chi2ndof = 10 #DATASET timepix3tel_dut_atlaspix_ebeam120 -#PASS T(942.734um,285.824um,0) R(10.6472deg,186.43deg,-1.35631deg) +#PASS T(942.726um,285.923um,0) R(10.6448deg,186.43deg,-1.35556deg) diff --git a/testing/test_io_read_rootobj.conf b/testing/test_io_read_rootobj.conf index 1e1f723ce44b5ea1487da680bed91bff85315e54..6cae5714b65b320aacee2c700da3fc06f511bd38 100644 --- a/testing/test_io_read_rootobj.conf +++ b/testing/test_io_read_rootobj.conf @@ -19,7 +19,7 @@ time_cut_frameedge = 10ns #DEPENDS test_io_write_rootobj.conf -#PASS Total efficiency of detector W0013_G02: 100(+0 -0.00534736)%, measured with 34446/34446 matched/total tracks +#PASS Total efficiency of detector W0013_G02: 100(+0 -0.00531864)%, measured with 34632/34632 matched/total tracks # Please note: # Pixels are masked in the [EventLoaderTimepix3] when reading in the online configuration files. diff --git a/testing/test_io_write_rootobj.conf b/testing/test_io_write_rootobj.conf index ccfdbb5c8c689cbb54831dfb14e044ec362b4393..741582913c2bdfaabc800050013344e547c5c9b5 100644 --- a/testing/test_io_write_rootobj.conf +++ b/testing/test_io_write_rootobj.conf @@ -32,4 +32,4 @@ time_cut_frameedge = 10ns #DATASET timepix3tel_ebeam120 -#PASS [F:FileWriter] Wrote 1717542 objects to 15 branches in file: +#PASS [F:FileWriter] Wrote 1722304 objects to 15 branches in file: diff --git a/testing/test_tracking_timepix3tel_ebeam120.conf b/testing/test_tracking_timepix3tel_ebeam120.conf index 14588c47943867e673022434b8312c6eef880c15..3df3c379ef59e462147a72087405d475520ad6c1 100644 --- a/testing/test_tracking_timepix3tel_ebeam120.conf +++ b/testing/test_tracking_timepix3tel_ebeam120.conf @@ -23,4 +23,4 @@ spatial_cut_abs = 200um, 200um #DATASET timepix3tel_ebeam120 -#PASS Ev: 18.8k Px: 6.23M Tr: 217.5k (11.6/ev) t = 3.7598s +#PASS Ev: 18.8k Px: 6.26M Tr: 217.4k (11.6/ev) t = 3.7598s diff --git a/testing/test_tracking_timepix3tel_ebeam120_gbl.conf b/testing/test_tracking_timepix3tel_ebeam120_gbl.conf index e86e080db0ef0c4d785edbb99490647c81ae3275..5e2590b5b1dc7a34db03f6ee20ba0285ffda06e8 100644 --- a/testing/test_tracking_timepix3tel_ebeam120_gbl.conf +++ b/testing/test_tracking_timepix3tel_ebeam120_gbl.conf @@ -24,4 +24,4 @@ volume_scattering_length = 304m #DATASET timepix3tel_ebeam120 -#PASS Ev: 18.8k Px: 6.23M Tr: 206.1k (11/ev) t = 3.7598s +#PASS Ev: 18.8k Px: 6.26M Tr: 206.0k (11/ev) t = 3.7598s diff --git a/testing/test_tracking_timepix3tel_ebeam120_multiplet.conf b/testing/test_tracking_timepix3tel_ebeam120_multiplet.conf index 0b9b88e74ecd6b2c4129e288edb6204475900cd8..320be1cd97f634c9db023152ef923a15c1456dcf 100644 --- a/testing/test_tracking_timepix3tel_ebeam120_multiplet.conf +++ b/testing/test_tracking_timepix3tel_ebeam120_multiplet.conf @@ -26,5 +26,4 @@ isolation_cut = 40um #DATASET timepix3tel_ebeam120 -#PASS Ev: 18.8k Px: 6.23M Tr: 188.5k (10/ev) t = 3.7598s - \ No newline at end of file +#PASS Ev: 18.8k Px: 6.26M Tr: 188.9k (10/ev) t = 3.7598s diff --git a/testing/test_tracking_timepix3tel_ebeam120_multiplet_gbl.conf b/testing/test_tracking_timepix3tel_ebeam120_multiplet_gbl.conf index 16e469924c883ab6636b93abbb6103d556ef63bc..c6ccb54d83ae5cdfbae468feeb3691cab5ebd47e 100644 --- a/testing/test_tracking_timepix3tel_ebeam120_multiplet_gbl.conf +++ b/testing/test_tracking_timepix3tel_ebeam120_multiplet_gbl.conf @@ -28,4 +28,4 @@ momentum=120GeV #DATASET timepix3tel_ebeam120 -#PASS Ev: 18.8k Px: 6.23M Tr: 188.5k (10/ev) t = 3.7598s \ No newline at end of file +#PASS Ev: 18.8k Px: 6.26M Tr: 188.9k (10/ev) t = 3.7598s