MAveragePulses.cc

Go to the documentation of this file.

#include "MAveragePulses.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QVector.hh"
#include "QGlobalHandle.hh"
#include "QAveragePulseHandle.hh"
#include "QAverageVector.hh"
#include "QBaseType.hh"
#include "QRunDataHandle.hh"
#include <math.h>
 
using std::cout;
using std::endl;
 
using namespace Diana;
REGISTER_MODULE(MAveragePulses)
 
// Init method is called before event loop
void MAveragePulses::Init(QEvent& ev) 
{
    // get parameters from config file
 
    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    ev.RequireByLabel<QPulse>(fPulseLabel);
 
 
    fValidityKind = GetString("ValidityKind","dataset",false);
    if(fValidityKind != "dataset" && fValidityKind != "run") {
      Panic("ValidityKind must be set to dataset or run");
    }
 
    fValidityStart = GetInt("ValidityStart",Q_INT_DEFAULT,false);
    fValidityEnd = GetInt("ValidityEnd",Q_INT_DEFAULT,false);
 
    fAlignmentOn = GetBool("AlignmentOn",false);
    fOutput = GetString("Output","avgpulses.root");
    fSubtractBaseline = GetBool("SubtractBaseline",false);
    if(fSubtractBaseline) fNumBaselinePoints = GetInt("BaselineNumPoints",0);
    if(fAlignmentOn) { 
        fMaxShift = GetInt("MaxShift",10);
        fFractionalShift = GetBool("FractionalShift",true);
    }
 
    fCurrentRun = -1;
    ev.Require<QHeader>("DAQ","Header");
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.RequireByLabel<QPulse>(fPulseLabel);
}
 
// Do method is called for each event, getting the event as argument
void MAveragePulses::Do( QEvent& ev)
{
    const QHeader& header = ev.Get<QHeader>("DAQ","Header");
    const QVector& pulse = ev.GetByLabel<QPulse>(fPulseLabel).GetSamples();
    const int run = header.GetRun();
 
    int chanId = ev.Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
 
    if (run != fCurrentRun)
    {
        fCurrentRun = run;
        // store run number into SourceRuns vector
        if(GetIteration()== 1) fSRuns.push_back(fCurrentRun);
    }
 
    if( GetIteration() == 1 && fChannelData.find(chanId) == fChannelData.end()) {
        //const QVector& pulse = ev.GetByLabel<QPulse>(fPulseLabel).GetSamples();
        size_t fNumSamples = pulse.Size();
        ChannelData thisChanData;
        thisChanData.fAvgPulse.Resize(fNumSamples);
        thisChanData.fAvgPulse.Initialize(0.);
        thisChanData.fNumPulses = 0;
        if(fAlignmentOn) { 
            thisChanData.fPeakPositions.Resize(fNumSamples);
            thisChanData.fPeakPositions.Initialize(0.);
            thisChanData.fAlignPos = -1.;
        }
        fChannelData[chanId] = thisChanData;
    }
 
 
    ChannelData& thisChanData = fChannelData[chanId];
    if(fAlignmentOn) 
    {
        //Alignment on
        double maxPos = pulse.GetMaxIndex();
        if(GetIteration() == 1) 
        {
            thisChanData.fPeakPositions[(int)maxPos] += 1;
            thisChanData.fNumPulses++;
        } 
        else 
        {
            // use maxPos to align pulses: standard procedure
            float delay = maxPos;
 
            // integer shift & fractional shift
            float stepreal = delay - thisChanData.fAlignPos;
            int stepint = (int) stepreal;
            stepreal = stepreal - stepint;
 
            QVector shifted = pulse;
 
            // first integer shift
            // maximum shift allowed. If exceeded, event is discarded.    
 
            if(stepint < fMaxShift && stepint > -fMaxShift)
            {
                // differentiate
                shifted.Differentiate();
 
                // mean derivative 
                int derivMeanSize = int(0.02*pulse.Size());
                QVector beginDeriv(derivMeanSize);
                QVector endDeriv(derivMeanSize);
 
                for(int i = 0; i < derivMeanSize; i++)
                { 
                    beginDeriv[i] = shifted[i];
                    endDeriv[i] = shifted[shifted.Size() - 1 -i];
                }
                double derivMeanBegin = beginDeriv.Sum(derivMeanSize)/derivMeanSize;    
                double derivMeanEnd = endDeriv.Sum(derivMeanSize)/derivMeanSize;
 
                // shift vector and sum mean derivative to the end or subtract mean derivative to the beginning
                shifted.Shift(-stepint);
                if (stepint >0)
                {
                    for(int i = 0; i < stepint; i++)
                    { 
                        shifted[shifted.Size() - stepint +i] += derivMeanEnd;
                    }
                }
                else if (stepint <0)
                {
                    for(int i = 0; i < stepint; i++)
                    { 
                        shifted[i] -= derivMeanBegin;
                    }
                }
 
                // integrate
                shifted.Integrate(pulse[0]);
                if (fFractionalShift)
                {
                    shifted.ShiftReal(-stepreal);
                    if (stepreal >0) shifted[shifted.Size() - 1] = shifted[shifted.Size() - 2] + derivMeanEnd;
                    else if (stepreal <0) shifted[0] =  shifted[1] - derivMeanBegin;
                }
                thisChanData.fAvgPulse += shifted;
                thisChanData.fNumPulses++;
            }
        }
    }
    else
    {
        // fAlignmentOn = false
        thisChanData.fAvgPulse += pulse;
        thisChanData.fNumPulses++;
    }
 
    return;
}
 
// Done method is called after event loop
void MAveragePulses::Done()
{
    std::map<int,ChannelData>::iterator chanIter = fChannelData.begin();
    while(chanIter != fChannelData.end()) {
        ChannelData& thisChanData = chanIter->second;
        int chanId = chanIter->first;
        
        if(chanIter!=fChannelData.begin()){
        --chanIter;
        int prevChan = (chanId-1);
        int prevEntry = (chanIter->first);
        ++chanIter;
            if(prevEntry != (prevChan)) {
                Error("Missing channel %d - Could not build AP",(prevChan));}
        }
        if(chanIter!= --fChannelData.end()){
        ++chanIter;
        int nextChan = (chanId+1);
        int nextEntry = (chanIter->first);
        --chanIter;
                if(nextEntry != (nextChan)) {
                Error("Missing channel %d - Could not build AP",(nextChan));}
        }
        
                
 
        if (thisChanData.fNumPulses != 0) {
            if(fAlignmentOn && GetIteration() == 1) {
                // alignment positions is the mode of the distribution peak distribution
                thisChanData.fAlignPos = thisChanData.fPeakPositions.GetMaxIndex() ;
                Debug("Channel[%d]: Number of events: %d, PeakMode %f",
                        chanId,thisChanData.fNumPulses, thisChanData.fAlignPos);
                thisChanData.fNumPulses = 0;
            } else {
                Debug("Channel[%d]: Number of events: %d", chanId, thisChanData.fNumPulses);
                thisChanData.fAvgPulse /= thisChanData.fNumPulses;
 
                // dump 
                std::string algo;
                if(fAlignmentOn) algo = "Aligned";
                else algo = "Normal";
                // Calculate and subtract baseline from average pulse
                if (fSubtractBaseline)
                {
                    double AvgBaseline = 0;
                    QVector AvgBaselineVector;
                    int nbsamples;
                    if(fNumBaselinePoints <=0) {
                        nbsamples = int (thisChanData.fAvgPulse.GetMaxIndex()*3./4);
                    } else {
                        nbsamples = fNumBaselinePoints;
                    }
                    for (int k = 0; k < nbsamples; ++k) 
                    {
                        AvgBaseline += thisChanData.fAvgPulse[k];
                    }   
                    AvgBaseline /= nbsamples;
                    AvgBaselineVector.Resize(thisChanData.fAvgPulse.Size());
                    AvgBaselineVector.Initialize(AvgBaseline);
                    thisChanData.fAvgPulse-= AvgBaselineVector;
                    char buf[128]; snprintf(buf,128,"_%.2d",chanId); 
                }
 
                // set number of events and source runs of QAverageVector
 
                thisChanData.fAvgPulse.SetSourceRuns(fSRuns); 
                thisChanData.fAvgPulse.SetNumEvents(thisChanData.fNumPulses); 
 
                // Fill Output
 
                QAveragePulseHandle aphandle(chanId);   
                if(fValidityKind=="dataset"){
                    // get dataset
                    GlobalHandle<QInt> dHandle("Dataset");
                    GlobalData().Get("",&dHandle,"");
                    int dataset = dHandle.Get();
                    aphandle.SetDataset(dataset);
                }
                else if (fValidityKind=="run"){
                    if (fSRuns.size()==1){
                        aphandle.SetRun(fSRuns[0]);
                    }
                    else Error("ValidityKind set to run is not permitted when reading multiple runs.");
                }
                else Error("ValidityKind can be set only to run or dataset. Check cfg.");
 
                aphandle.Set(thisChanData.fAvgPulse);
                GlobalData().Set(&aphandle, fOutput);                  
            }
        }
        else {
                Error("Could not build average pulse for channel %i: , 0 events passing the selection cuts",chanId);
        }       
        chanIter++;
    }
 
    if(fAlignmentOn && GetIteration() == 1) SetRunAgain(true);
    else SetRunAgain(false);
 
}