MBaselineModule.cc

Go to the documentation of this file.

/* Diana Reconstruction program
 *
 * Author: Laura Kogler 29/01/07
 * $Id: MBaselineModule.cc 544 2006-12-11 17:57:43Z vignatim $
 *
 * Class MBaselineModule: Calculate baseline parameters
 *
*/
 
 
#include "MBaselineModule.hh"
#include "QVector.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QRunData.hh"
#include "QBaselineData.hh"
#include "QRunDataHandle.hh"
#include "QTypes.hh"
#include <math.h>
#include <gsl/gsl_fit.h>
 
using namespace Diana;
 
REGISTER_MODULE(MBaselineModule)
 
// Init method is called before event loop
void MBaselineModule::Init(QEvent& ev) 
{
        
        fNumberOfBaselinePoints = GetInt("NumPoints",0,false);
    ev.Add<QBaselineData>("BaselineData");
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Require<QHeader>("DAQ","Header");
    fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
    fTomV = GetBool("ConvertTomV",true,false);
    ev.RequireByLabel<QPulse>(fPulseLabel);
        
}
 
// Doit method is called for each event, getting the event as argument
void MBaselineModule::Do( QEvent& ev) 
{
 
    const QVector& samples = ev.GetByLabel<QPulse>(fPulseLabel).GetSamples();
    const int run = ev.Get<QHeader>("DAQ","Header").GetRun();
    const int channel = ev.Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
    QBaselineData& bData = ev.Get<QBaselineData>("BaselineData");
 
    size_t nSamples = samples.Size();
    const double* SampleArray = samples.GetConstArray();
    size_t nPoints;
    if(fNumberOfBaselinePoints < 0.) 
        nPoints = nSamples;
    else if(fNumberOfBaselinePoints == 0) {
        nPoints = int(ev.Get<QPulseInfo>("DAQ","PulseInfo").GetMasterSample().GetSampleIndex()*3./4);
    } else {
        nPoints = fNumberOfBaselinePoints;
    }
    if(nPoints > nSamples) {
        Warn("Number of points (%d) exceeds number of samples (%d): setting it to %d",
                nPoints, nSamples, nSamples);
        nPoints = nSamples;
    }
    Double_t* TimeArray = new Double_t[nPoints];
    for(size_t i = 0; i < nPoints; i++)
        TimeArray[i] = (Double_t)i;
 
    fBaseline = samples.Sum(nPoints,0);
    fBaseline /=  nPoints;
    fRightBaseline = samples.Sum(nPoints,samples.Size()-1-nPoints)/nPoints;
    fRightLeftBaseline = fRightBaseline - fBaseline;
    
    double cov00, cov01, cov11;
    gsl_fit_linear (TimeArray, 1, SampleArray, 1, nPoints, 
            &fBaselineIntercept, &fBaselineSlope, &cov00, &cov01, &cov11, &fBaselineRMS);
    
 
    fBaselineRMS = sqrt(fBaselineRMS/nPoints); // this gives RMS of deviation from best fit to baseline
    double baselineFlatRMS = samples.GetRMS(nPoints,0); 
 
    fBaselineSlopeRMSWindow = fBaselineSlope/fBaselineRMS*(double)nSamples;    
    fRightBaselineInRMS = fRightBaseline/fBaselineRMS;
    fRightLeftBaselineInRMS = fRightLeftBaseline/fBaselineRMS;
 
    QRunDataHandle rHandle(run);
    GlobalData().Get("",&rHandle,"");
    const QRunData& runData = rHandle.Get();
    const QChannelRunData& chanRunData = runData.GetChannelRunData(channel);
    double ADC2mV = 1;
    if(fTomV) ADC2mV = chanRunData.fADC2mV;
 
    // Event assignement
    bData.SetBaseline(fBaseline*ADC2mV);
    bData.SetRightBaseline(fRightBaseline*ADC2mV);
    bData.SetRightLeftBaseline(fRightLeftBaseline*ADC2mV);
    bData.SetRightBaselineInRMS(fRightBaselineInRMS);
    bData.SetRightLeftBaselineInRMS(fRightLeftBaselineInRMS);
    bData.SetBaselineIntercept(fBaselineIntercept*ADC2mV);
    bData.SetBaselineSlope(fBaselineSlope*ADC2mV);
    bData.SetBaselineSlopeRMSWindow(fBaselineSlopeRMSWindow);
    bData.SetBaselineRMS(fBaselineRMS*ADC2mV);
    bData.SetBaselineFlatRMS(baselineFlatRMS*ADC2mV);
 
    Debug("Mean is %f, Intercept is %f, Slope is %f, RMS is %f",fBaseline, fBaselineIntercept,fBaselineSlope, fBaselineRMS);
    delete [] TimeArray;
}
 
// Done method is called after event loop
void MBaselineModule::Done() 
{
 
}