MPulserFlagByRegularTiming.cc

Go to the documentation of this file.

#include "MPulserFlagByRegularTiming.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QRunData.hh"
#include "QBool.hh"
#include "QBaseType.hh"
#include "QCoincidenceTiming.hh"
#include <sstream>
 
using std::flush;
using std::map;
using std::stringstream;
using namespace Diana;
 
REGISTER_MODULE(MPulserFlagByRegularTiming)
 
void MPulserFlagByRegularTiming::Init(QEvent& ev)
{
    // we manually handle filters
    GetBool("SkipEvents",false,false);
 
    if (GetIteration() == 1) {
 
        fExpectedTimeInterval    = GetDouble("ExpectedTimeInterval", -1);
        if (fExpectedTimeInterval <= 0) {
            fExpectedTimeInterval = 305.536;
        }
 
        fTimeIntervalUncertainty = GetDouble("TimeIntervalUncertainty", 0.500);
        fTimeIntervalStepSize    = GetDouble("TimeIntervalStepSize", 0.008);
 
        int numberOfTimeAdjustments = GetInt("NumberOfTimeAdjustments", 0);
        for (int i = 1; i <= numberOfTimeAdjustments; ++i) {
            stringstream parameterName;
            parameterName << "TimeAdjustment_" << i << flush;
            double timeAdjustment = GetDouble(parameterName.str(), 0.0);
            fTimeAdjustments.push_back(timeAdjustment);
        }
        
        fTimeTolerance     = GetDouble("TimeTolerance", 1.0);
        fFalsePositiveRate = GetDouble("FalsePositiveRate", 0.001);
        
        fUseCoincidences = GetBool("UseCoincidences", false);
        
        if (fUseCoincidences) {
            fNecessaryCoincidentEvents  = GetInt("NecessaryCoincidentEvents", 2);
            fSufficientCoincidentEvents = GetInt("SufficientCoincidentEvents", 5);
            fMaxRandomCoincidentEvents  = GetInt("MaxRandomCoincidentEvents", 2);
        }
   
    } else {
        ev.Add<QBool>("IsPulser");
        ev.Add<QInt>("Goodness");
        ev.Add<QInt>("GoodnessThreshold");
        if(fUseCoincidences) ev.Require<QCoincidenceTiming>("CoincidenceTiming_PulserFlag","Timing");
    }
 
    ev.Require<QHeader>("DAQ","Header");
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
}
 
void MPulserFlagByRegularTiming::Do(QEvent& ev)
{
 
    const QHeader& header = ev.Get<QHeader>("DAQ","Header"); 
    const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo"); 
    const int& channel = pulseInfo.GetChannelId();
            
    const double time = header.GetTime().GetFromStartRunNs() / 1e9;
    
   
    if (GetIteration() == 1 && ev.Get<QBool>("SkipEvent")==false) {
        fTimingAnalyzer[channel].AddValue(time);
    }
    else if (GetIteration() == 2) {
 
        unsigned int goodness = 0;
        unsigned int goodnessThreshold = 0;
        
        if (fTimingAnalyzer.count(channel) != 0) {
            goodness = fTimingAnalyzer[channel].GetGoodness(time);           
            goodnessThreshold = fTimingAnalyzer[channel].GetGoodnessThreshold();
        }
        
        bool isPulser = (goodness > goodnessThreshold);
        
        if (fUseCoincidences) {
            const QCoincidenceTiming& coincTiming = ev.Get<QCoincidenceTiming>("CoincidenceTiming_PulserFlag","Timing");
            int numberOfCoincidentEvents_SamePulserGroup = coincTiming.fNumberOfCoincidentEvents_SameHeaterGroup;
            int numberOfCoincidentEvents = coincTiming.fNumberOfCoincidentEvents;
            
            // require at least some coincident events in the same heater group
            isPulser = isPulser
                && numberOfCoincidentEvents_SamePulserGroup
                    >= fNecessaryCoincidentEvents;
            
            // also flag events where a whole heater group fired
            isPulser = isPulser 
                || 
                (
                    numberOfCoincidentEvents_SamePulserGroup
                    >= fSufficientCoincidentEvents
                    &&
                    numberOfCoincidentEvents
                    <= numberOfCoincidentEvents_SamePulserGroup
                       + fMaxRandomCoincidentEvents
                );
        }
        
        ev.Get<QBool>("IsPulser") = isPulser;
        ev.Get<QInt>("Goodness") = goodness;
        ev.Get<QInt>("GoodnessThreshold") = goodnessThreshold;
      
    }
    
}
 
void MPulserFlagByRegularTiming::Done()
{
    // called at the end of the event loop
    if (GetIteration() == 1) {
    
        for (map<int, QTimingAnalyzer>::iterator channelIter
                 = fTimingAnalyzer.begin();
             channelIter != fTimingAnalyzer.end();
             ++channelIter) {
            
            Debug("Operating on channel %d", channelIter->first);
            
            QTimingAnalyzer& timingAnalyzer = channelIter->second;
            
            timingAnalyzer.SetExpectedPeriod(fExpectedTimeInterval);
            timingAnalyzer.SetAdjustments(fTimeAdjustments);
            timingAnalyzer.SetTolerance(fTimeTolerance);
            
            timingAnalyzer.FindPeriod(fTimeIntervalUncertainty,
                                      fTimeIntervalStepSize);
            
            Info("Best time interval for channel %d:  %f",
                 channelIter->first,
                 timingAnalyzer.GetPeriod());
            
            timingAnalyzer.ComputeGoodnessThreshold(fFalsePositiveRate);
            // argument is false positive rate under the null hypothesis
            
        }
        SetRunAgain(true);
        
    }
    else {
        SetRunAgain(false);
    }
}