TRestEventTimeSelectionProcess.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
 
#include "TRestEventTimeSelectionProcess.h"
 
using namespace std;
 
ClassImp(TRestEventTimeSelectionProcess);
 
TRestEventTimeSelectionProcess::TRestEventTimeSelectionProcess() { Initialize(); }
 
void TRestEventTimeSelectionProcess::Initialize() {
    SetSectionName(this->ClassName());
    fEvent = nullptr;
    fFileWithTimes = "";
    fIsActiveTime = true;
    fDelimiter = ',';
    fStartEndTimes.clear();
    fTimeOffsetInSeconds = 0;
    fTimeStartMarginInSeconds = 0;
    fTimeEndMarginInSeconds = 0;
    fUseRunStartAndEndTimes = true;
    fNEventsRejected = 0;
    fNEventsSelected = 0;
    fTotalTimeInSeconds = 0;
}
 
void TRestEventTimeSelectionProcess::InitProcess() {
    // Read the file with the time ranges
    if (!fFileWithTimes.empty()) {
        fStartEndTimes = ReadFileWithTimes(fFileWithTimes, fDelimiter);
    }
    if (fUseRunStartAndEndTimes) {
        TRestRun* run = GetRunInfo();
        if (run) {
            auto startTimeStamp = run->GetStartTimestamp();
            auto endTimeStamp = run->GetEndTimestamp();
            ApplyStartRunTime(startTimeStamp);
            ApplyEndRunTime(endTimeStamp);
        } else {
            RESTWarning << "No run information available to get TRestRun start and end times." << RESTendl;
        }
    }
    fTotalTimeInSeconds = CalculateTotalTimeInSeconds();
    fNEventsRejected = 0;
    fNEventsSelected = 0;
}
 
std::vector<Interval> TRestEventTimeSelectionProcess::ReadFileWithTimes(std::string fileWithTimes,
                                                                        Char_t delimiter) {
    std::vector<Interval> startEndTimes;
    string line;
    ifstream file(fileWithTimes);
    if (file.is_open()) {
        while (getline(file, line)) {
            if (line[0] == '#') {  // understand as comment
                continue;
            }
            std::istringstream lineStream(line);
            std::string startDate, endDate;
            if (std::getline(lineStream, startDate, delimiter) &&
                std::getline(lineStream, endDate, delimiter)) {
                // check if the time format is correct. TODO: use better way to check
                // (StringToTimeStamp usually returns a negative big number if not)
                if (StringToTimeStamp(startDate) < 0 || StringToTimeStamp(endDate) < 0) {
                    continue;
                }
                TTimeStamp sts = StringToTimeStamp(startDate);
                TTimeStamp ets = StringToTimeStamp(endDate);
                startEndTimes.emplace_back(sts, ets);
            }
        }
        file.close();
    }
 
    // sort by start time and then by end time (lexicographically)
    std::sort(startEndTimes.begin(), startEndTimes.end());
 
    return startEndTimes;
}
 
Double_t TRestEventTimeSelectionProcess::CalculateTotalTimeInSeconds() {
    Double_t totalTime = 0;
    for (auto se : fStartEndTimes) {
        TTimeStamp startTime = se.first;
        TTimeStamp endTime = se.second;
        // Reduce the time by the margin in both sides
        startTime.Add(TTimeStamp(fTimeStartMarginInSeconds));
        endTime.Add(TTimeStamp(-fTimeEndMarginInSeconds));
        auto timeDiff = endTime.AsDouble() - startTime.AsDouble();
        if (timeDiff < 0) {
            RESTDebug << "End time is before start time in time range: " << se.first << " to " << se.second
                      << RESTendl;
            continue;
        }
        totalTime += endTime.AsDouble() - startTime.AsDouble();
    }
    return totalTime;
}
 
TRestEvent* TRestEventTimeSelectionProcess::ProcessEvent(TRestEvent* inputEvent) {
    fEvent = inputEvent;
 
    TTimeStamp eventTime = fEvent->GetTimeStamp();
    eventTime.Add(TTimeStamp(fTimeOffsetInSeconds));
 
    Bool_t isInsideAnyTimeRange = false;
    for (auto se : fStartEndTimes) {
        TTimeStamp startTime = se.first;
        TTimeStamp endTime = se.second;
        if (fIsActiveTime) {
            // Reduce the active time window by the margin in both sides
            startTime.Add(TTimeStamp(fTimeStartMarginInSeconds));
            endTime.Add(TTimeStamp(-fTimeEndMarginInSeconds));
        } else {
            // Increase the dead time window by the margin in both sides
            startTime.Add(TTimeStamp(-fTimeStartMarginInSeconds));
            endTime.Add(TTimeStamp(fTimeEndMarginInSeconds));
        }
 
        if (eventTime >= startTime && eventTime <= endTime) {
            isInsideAnyTimeRange = true;
            break;
        }
    }
 
    // Decide if the event is selected or rejected based on the time ranges
    // and their meaning (active or dead periods of time).
    if (fIsActiveTime) {             // time ranges represent active periods of time
        if (isInsideAnyTimeRange) {  // time is inside an active period of time
            fNEventsSelected++;
            return fEvent;
        }
    } else {                          // time ranges represent dead periods of time
        if (!isInsideAnyTimeRange) {  // time is outside all dead period of time
            fNEventsSelected++;
            return fEvent;
        }
    }
 
    // rejected events are not returned
    fNEventsRejected++;
    return nullptr;
}
 
void TRestEventTimeSelectionProcess::EndProcess() {
    // Write here the jobs to do when all the events are processed
}
 
void TRestEventTimeSelectionProcess::ApplyStartRunTime(const TTimeStamp& runStart) {
    size_t startIndex = 0;
    bool isInsideTimeRange = false;
    for (auto se : fStartEndTimes) {
        TTimeStamp s = se.first;
        TTimeStamp e = se.second;
 
        if (runStart < s) {
            isInsideTimeRange = false;
            break;
        }
 
        if (runStart >= s && runStart <= e) {
            isInsideTimeRange = true;
            break;
        }
 
        startIndex++;
    }
 
    if (isInsideTimeRange) {
        // modify the start time of the found interval
        fStartEndTimes[startIndex].first = runStart;
    }
    // remove all intervals before the run start time
    fStartEndTimes.erase(fStartEndTimes.begin(), fStartEndTimes.begin() + startIndex);
}
 
void TRestEventTimeSelectionProcess::ApplyEndRunTime(const TTimeStamp& runEnd) {
    size_t endIndex = 0;
    bool isInsideTimeRange = false;
    for (auto se : fStartEndTimes) {
        TTimeStamp s = se.first;
        TTimeStamp e = se.second;
 
        if (runEnd < s) {
            isInsideTimeRange = false;
            break;
        }
 
        if (runEnd >= s && runEnd <= e) {
            isInsideTimeRange = true;
            break;
        }
 
        endIndex++;
    }
 
    if (isInsideTimeRange) {
        // modify the end time of the found interval
        fStartEndTimes[endIndex].second = runEnd;
        endIndex++;  // to erase from the next interval
    }
    // remove all intervals after the run end time
    fStartEndTimes.erase(fStartEndTimes.begin() + endIndex, fStartEndTimes.end());
}
 
std::string TRestEventTimeSelectionProcess::GetTimeStampCut(std::string timeStampObsName, Bool_t useOffset,
                                                            Bool_t useMargins, Int_t nTimes) {
    std::string timeCut = "";
    std::string timeStampObsNameWithOffset = timeStampObsName;
    if (useOffset && fTimeOffsetInSeconds != 0) {
        timeStampObsNameWithOffset += "+" + to_string(fTimeOffsetInSeconds);
    }
    if (nTimes < 0) nTimes = fStartEndTimes.size();
    Int_t c = 0;
    for (auto se : fStartEndTimes) {
        if (c++ >= nTimes) break;
        auto startTime = se.first;
        auto endTime = se.second;
        // Reduce the time by the margin in both sides
        if (useMargins) {
            if (fIsActiveTime) {
                // Reduce the active time window by the margin in both sides
                startTime.Add(fTimeStartMarginInSeconds);
                endTime.Add(-fTimeEndMarginInSeconds);
            } else {
                // Increase the dead time window by the margin in both sides
                startTime.Add(-fTimeStartMarginInSeconds);
                endTime.Add(fTimeEndMarginInSeconds);
            }
        }
 
        if (startTime >= endTime) {
            continue;
        }
 
        // Build the cut string
        if (!timeCut.empty()) {
            if (fIsActiveTime)
                timeCut += " || ";  // inside ANY time range
            else
                timeCut += " && ";  // outside ALL time ranges
        }
        if (!fIsActiveTime) timeCut += "!";  // NOT inside the time range
        // inside the time range
        timeCut += "(";
        timeCut += "(" + timeStampObsNameWithOffset + ">=" + to_string(startTime) + ")";
        timeCut += " && ";
        timeCut += "(" + timeStampObsNameWithOffset + "<=" + to_string(endTime) + ")";
        timeCut += ")";
    }
    return timeCut;
}
void TRestEventTimeSelectionProcess::PrintMetadata() {
    BeginPrintProcess();
    std::string typeOfTime = fIsActiveTime ? "Active" : "Dead";
 
    RESTMetadata << "File with times: " << fFileWithTimes << RESTendl;
    // print periods
    RESTMetadata << "Use run start and end: " << (fUseRunStartAndEndTimes ? "true" : "false") << RESTendl;
    RESTMetadata << "Offset time: " << fTimeOffsetInSeconds << " seconds" << RESTendl;
    RESTMetadata << "Start margin time: " << fTimeStartMarginInSeconds << " seconds" << RESTendl;
    RESTMetadata << "End margin time: " << fTimeEndMarginInSeconds << " seconds" << RESTendl;
    RESTMetadata << typeOfTime << " time periods: " << RESTendl;
    for (auto se : fStartEndTimes) {
        std::string startStr = ToDateTimeString(se.first);
        std::string endStr = ToDateTimeString(se.second);
        RESTMetadata << startStr << " to " << endStr << RESTendl;
    }
 
    // Get total time in seconds
    TTimeStamp totalTime = TTimeStamp(fTotalTimeInSeconds, 0);
    if (!fStartEndTimes.empty()) {
        TTimeStamp firstTime = fStartEndTimes.front().first;
        TTimeStamp lastTime = fStartEndTimes.back().second;
        totalTime = lastTime - firstTime;
    }
 
    double fractionOfTime = fTotalTimeInSeconds / totalTime.AsDouble() * 100;
    std::string fractionOfTimeStr = StringWithPrecision(fractionOfTime, 4) + " %";
    if ((Int_t)(fTotalTimeInSeconds / 24 / 3600) != 0)  // order of days
        RESTMetadata << "Total " << typeOfTime << " time: " << fTotalTimeInSeconds / 24 / 3600 << " days"
                     << " (" << fractionOfTimeStr << ")" << RESTendl;
    else if ((Int_t)(fTotalTimeInSeconds / 3600) != 0)  // order of hours
        RESTMetadata << "Total " << typeOfTime << " time: " << fTotalTimeInSeconds / 3600 << " hours"
                     << " (" << fractionOfTimeStr << ")" << RESTendl;
    else if ((Int_t)(fTotalTimeInSeconds / 60) != 0)  // order of minutes
        RESTMetadata << "Total " << typeOfTime << " time: " << fTotalTimeInSeconds / 60 << " minutes"
                     << " (" << fractionOfTimeStr << ")" << RESTendl;
    else
        RESTMetadata << "Total " << typeOfTime << " time: " << fTotalTimeInSeconds << " seconds"
                     << " (" << fractionOfTimeStr << ")" << RESTendl;
 
    RESTMetadata << "Number of events rejected: " << fNEventsRejected << " ("
                 << fNEventsRejected * 1. / (fNEventsRejected + fNEventsSelected) * 100 << " %)" << RESTendl;
    RESTMetadata << "Number of events selected: " << fNEventsSelected << " ("
                 << fNEventsSelected * 1. / (fNEventsRejected + fNEventsSelected) * 100 << " %)" << RESTendl;
 
    EndPrintProcess();
}