MCoincidenceTiming.cc

Go to the documentation of this file.

#include "MCoincidenceTiming.hh"
#include <algorithm>
#include "QDbDetector.hh"
#include "QError.hh"
#include "QEvent.hh"
#include "QFiltersData.hh"
#include "QMainPulse.hh"
#include "QOBData.hh"
#include "QRawEvent.hh"
#include "QRunData.hh"
 
using std::vector;
 
REGISTER_MODULE(MCoincidenceTiming)
 
MCoincidenceTiming::EventInfo::EventInfo(const QEvent *ev)
{
    const QRawEvent& raw = ev->GetRawEvent();
    const QMainPulse& pulse = raw.GetPulse();
    fChannel = pulse.GetChannelId();
    fEventNumber = raw.GetEventNumber();
    fHeaterChannel = -1;
    fPassed = ev->GetFilters().GetPassed();
    fPulseTimeOffset = ev->GetOB().GetStartTimeOffset();
    fRelative = -1;
    fTimeFromStartRunNs = raw.GetTime().GetFromStartRunNs();
}
 
void MCoincidenceTiming::EventInfo::ApplyOffset(double samplingPeriod)
{
    long long offsetNs = (long long) (fPulseTimeOffset * samplingPeriod * 1e9);
    fTimeFromStartRunNs += offsetNs;
}
 
MCoincidenceTiming::MCoincidenceTiming(QSequence* s)
    : QModule("CoincidenceTiming",s)
{
}
 
MCoincidenceTiming::~MCoincidenceTiming()
{
}
 
void MCoincidenceTiming::Init()
{
    if (GetIteration() == 1) {
        fIsRunDataLoaded = false;
        
        fUsePulseStartTimeOffset = GetBool("UsePulseStartTimeOffset", false);
        fTimeWindow = GetDouble("TimeWindow", 0.050);
        fDetectorSetupFile = GetString("DetectorSetupFile","DB");
        fFilteredEventsType = GetString("FilteredEventsType", "", false);
    }
}
 
QEvent* MCoincidenceTiming::Do(QEvent* ev)
{
    const QRawEvent& raw = ev->GetRawEvent();
    int thisRun = raw.GetRun();
    
    EventInfo thisEvent(ev);
    
    // Fill fRelatedMap with heater channels and relatives
    if (fRelatedMap.count(thisEvent.fChannel) == 0) {
        try {
            fRelatedMap[thisEvent.fChannel].fHeaterChan
                = QDbDetector::GetPulserId(thisRun, thisEvent.fChannel,
                                           fDetectorSetupFile);
        }
        catch (const QError& err) {
            Info("%s", err.GetDescription().c_str());
            fRelatedMap[thisEvent.fChannel].fHeaterChan = -1;
        }
        
        vector<int> relatives;
        try {
            relatives = QDbDetector::GetRelativeChannels(thisRun,
                                                         thisEvent.fChannel,
                                                         fDetectorSetupFile);
        }
        catch (const QError& err) {
            Info("%s", err.GetDescription().c_str());
        }
        
        if (relatives.size() > 0) {
            fRelatedMap[thisEvent.fChannel].fRelative = relatives[0];
        }
        else {
            fRelatedMap[thisEvent.fChannel].fRelative = -1;
        }
    }
    
    thisEvent.fHeaterChannel = fRelatedMap[thisEvent.fChannel].fHeaterChan; 
    thisEvent.fRelative = fRelatedMap[thisEvent.fChannel].fRelative;
    
    if (fUsePulseStartTimeOffset) {
        if (!fIsRunDataLoaded) {
            if ( const QRunData *runData = dynamic_cast<const QRunData*>(
                 SeqAuxData().GetQObject("RunData")) ) {
                double samplingFrequency = runData->GetSamplingFrequency();
                fSamplingPeriod = 0;
                if (samplingFrequency != 0) {
                    fSamplingPeriod = 1.0 / samplingFrequency;
                }
            }
            fIsRunDataLoaded = true;
        }
        thisEvent.ApplyOffset(fSamplingPeriod);
    }
    
    if (GetIteration() == 1) {
        fEventList.push_back(thisEvent);
    }
    else if (GetIteration() == 2) {
        double timeSinceEvent_SameChannel = -999;
        double timeSinceEvent_DifferentChannel = -999;
        double timeUntilEvent_SameChannel = -999;
        double timeUntilEvent_DifferentChannel = -999;
        
        int numberOfCoincidentEvents = 0;
        int numberOfCoincidentEvents_SameHeaterGroup = 0;
        
        int thisIndex = fEventListIndex[thisEvent.fEventNumber];
        int index = thisIndex;
        double timeDifference = 0;
        
        bool previousEventFound_SameChannel = false;
        bool previousEventFound_DifferentChannel = false;
        bool nextEventFound_SameChannel = false;
        bool nextEventFound_DifferentChannel = false;
        bool allFound = false;
        
        // Look backward
        index = thisIndex - 1;
        while (index >= 0 && (!allFound || timeDifference < fTimeWindow)) {
            const EventInfo& earlierEvent = fEventList[index];
            if (earlierEvent.fPassed) {
                timeDifference = (thisEvent.fTimeFromStartRunNs
                    - earlierEvent.fTimeFromStartRunNs) / 1e9;
                if (!previousEventFound_SameChannel
                    && earlierEvent.fChannel == thisEvent.fChannel) {
                    timeSinceEvent_SameChannel = timeDifference;
                    previousEventFound_SameChannel = true;
                }
                if (!previousEventFound_DifferentChannel
                    && earlierEvent.fChannel != thisEvent.fChannel
                    && earlierEvent.fChannel != thisEvent.fRelative) {
                    timeSinceEvent_DifferentChannel = timeDifference;
                    previousEventFound_DifferentChannel = true;
                }
                if (timeDifference < fTimeWindow
                    && earlierEvent.fChannel != thisEvent.fChannel
                    && earlierEvent.fChannel != thisEvent.fRelative) {
                    ++numberOfCoincidentEvents;
                    if (earlierEvent.fHeaterChannel
                        == thisEvent.fHeaterChannel) {
                        ++numberOfCoincidentEvents_SameHeaterGroup;
                    }
                }
            }
            --index;
            allFound = previousEventFound_SameChannel
                    && previousEventFound_DifferentChannel;
        }
        
        // Look forward
        index = thisIndex + 1;
        timeDifference = 0;
        allFound = false;
        while (index < (int) fEventList.size()
               && (!allFound || timeDifference < fTimeWindow)) {
            const EventInfo& laterEvent = fEventList[index];
            if (laterEvent.fPassed) {
                timeDifference = (laterEvent.fTimeFromStartRunNs
                    - thisEvent.fTimeFromStartRunNs) / 1e9;
                if (!nextEventFound_SameChannel
                    && laterEvent.fChannel == thisEvent.fChannel) {
                    timeUntilEvent_SameChannel = timeDifference;
                    nextEventFound_SameChannel = true;
                }
                if (!nextEventFound_DifferentChannel
                    && laterEvent.fChannel != thisEvent.fChannel
                    && laterEvent.fChannel != thisEvent.fRelative) {
                    timeUntilEvent_DifferentChannel = timeDifference;
                    nextEventFound_DifferentChannel = true;
                }
                if (timeDifference < fTimeWindow
                    && laterEvent.fChannel != thisEvent.fChannel
                    && laterEvent.fChannel != thisEvent.fRelative) {
                    ++numberOfCoincidentEvents;
                    if (laterEvent.fHeaterChannel
                        == thisEvent.fHeaterChannel) {
                        ++numberOfCoincidentEvents_SameHeaterGroup;
                    }
                }
            }
            ++index;
            allFound = nextEventFound_SameChannel
                    && nextEventFound_DifferentChannel;
        }
        ev->AuxData().SetDouble(
            "TimeSince" + fFilteredEventsType + "Event_SameChannel",
            timeSinceEvent_SameChannel, "save"
        );
        ev->AuxData().SetDouble(
            "TimeSince" + fFilteredEventsType + "Event_DifferentChannel",
            timeSinceEvent_DifferentChannel, "save"
        );
        ev->AuxData().SetDouble(
            "TimeUntil" + fFilteredEventsType + "Event_SameChannel",
            timeUntilEvent_SameChannel, "save"
        );
        ev->AuxData().SetDouble(
            "TimeUntil" + fFilteredEventsType + "Event_DifferentChannel",
            timeUntilEvent_DifferentChannel, "save"
        );
        ev->AuxData().SetInt(
            "NumberOfCoincident" + fFilteredEventsType + "Events",
            numberOfCoincidentEvents, "save"
        );
        ev->AuxData().SetInt(
            "NumberOfCoincident" + fFilteredEventsType
            + "Events_SameHeaterGroup",
            numberOfCoincidentEvents_SameHeaterGroup, "save"
        );
    }
 
    return ev;
}
 
void MCoincidenceTiming::Done()
{
    if (GetIteration() == 1) {
        // Sort the event list by time
        sort(fEventList.begin(), fEventList.end());
        // Build the index by event number
        for (unsigned int index = 0; index != fEventList.size(); ++index) {
            int eventNumber = fEventList[index].fEventNumber;
            fEventListIndex[eventNumber] = index;
        }
        SetRunAgain(true);
    }    
}