MNoiseAvgPowerSpectrum.cc

Go to the documentation of this file.

#include "MNoiseAvgPowerSpectrum.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QVector.hh"
#include "QVectorC.hh"
#include "QRealComplexFFT.hh"
#include "QAverageNoiseHandle.hh"
#include "QAverageVector.hh"
#include "QBaseType.hh"
#include <iostream>
#include <math.h>
 
using std::cout;
using std::endl;
 
using namespace Diana;
REGISTER_MODULE(MNoiseAvgPowerSpectrum)
 
// Init method is called before event loop
void MNoiseAvgPowerSpectrum::Init(Diana::QEvent& ev) 
{
  // get parameters from config file
  
  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);
  
  fProcessedFirstEvent = false;
  fOutput = GetString("Output","avgnoise.root");
  fWindow = GetString("WindowType","None");  
  fWindowType = QFFT::StrToWindowType(fWindow);    
  if(fWindowType == QFFT::WT_None)  {
      Warn("Window type not defined, using rectangular window");
      fWindowType = QFFT::WT_Rect;
  }
  else Info("Window type: %d (%s)", fWindowType,fWindow.c_str());
  if(fWindowType != QFFT::WT_Rect) 
      fCoherentGain = GetInt("WindowCoherentGain",2);
  else fCoherentGain = false;
 
  fPulseLabel = GetString("PulseLabel","DAQ@Pulse",false);
  fCurrentRun = -1;
  fSRuns.clear();
 
  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 MNoiseAvgPowerSpectrum::Do( Diana::QEvent& ev) 
{
  const QHeader& header = ev.Get<QHeader>("DAQ","Header");
  const int run = header.GetRun();
  
 
  // inspect first event to determine number of points in the acquired window
  if (!fProcessedFirstEvent) {
    // validity time interval for DB MV FIXME
    fBeginValidity = header.GetTime().GetTimeSec();
    fEndValidity = fBeginValidity + 1000;
    fProcessedFirstEvent = true;
  }
 
    // store run number into SourceRuns vector
    if (run != fCurrentRun)
      {
          fCurrentRun = run;
          fSRuns.push_back(fCurrentRun);
      }
 
    int channel = ev.Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
  
    const QVector& noisePulse = ev.GetByLabel<QPulse>(fPulseLabel).GetSamples();
 
    if( fNoiseAvgPowerSpectrum.find(channel) == fNoiseAvgPowerSpectrum.end()) {
        fNumFreq = noisePulse.Size();
        fNoiseAvgPowerSpectrum[channel].Resize(fNumFreq);
        fNoiseAvgPowerSpectrum[channel].Initialize(0);
        fCount[channel] = 0;
    }
 
 
    //  bool noPulse = ev->GetCountPulses().GetNumberOfPulses() == 0;
    fNumFreq = noisePulse.Size();
    double mean = noisePulse.Sum(fNumFreq)/fNumFreq;
    QVector zeroMeanPulse = noisePulse;
    QVector meanPulse(fNumFreq); meanPulse.Initialize(mean);
    zeroMeanPulse -= meanPulse;
    QRealComplexFFT transformer(fNumFreq);
 
    transformer.SetWindowType(fWindowType, fCoherentGain);
    QVectorC transformedPulse(fNumFreq);
    // normal pulse
    int err = transformer.TransformToFreq(zeroMeanPulse, transformedPulse);
    if(err != 0) Error("QFFTRealComplex::TransformToFreq: pulse");
    QVector realPart = transformedPulse.Re();
    QVector imagPart = transformedPulse.Im();
   
    // build average
    // sum for now; divide by counts in Done() method
    for (int k = 1; k < fNumFreq; k++) { 
        double magSquared = realPart[k]*realPart[k] + imagPart[k]*imagPart[k];
        fNoiseAvgPowerSpectrum[channel][k] += magSquared;
    }
   fNoiseAvgPowerSpectrum[channel][0]+=mean*mean;
        
    // increment count of number of noise samples
    ++fCount[channel];
    return ;
}
 
// Done method is called after event loop
void MNoiseAvgPowerSpectrum::Done()
{
    // average magnitude squared of noise frequency spectrum
    std::map<int,QAverageVector>::const_iterator iter = fNoiseAvgPowerSpectrum.begin();
    while (iter != fNoiseAvgPowerSpectrum.end() ) {
        int channel = iter->first;
        
        if(iter != fNoiseAvgPowerSpectrum.begin()){
        --iter; 
        int prevChan = (channel-1);
        int prevEntry = (iter->first);
        ++iter;
                if(prevEntry != (prevChan)) {
                Error("Missing channel %d - Could not build ANPS",(prevChan));}
        }
        if(iter != --fNoiseAvgPowerSpectrum.end()){
        ++iter;
        int nextChan = (channel+1);
        int nextEntry = (iter->first);
        --iter;
                if(nextEntry != (nextChan)) {
                Error("Missing channel %d - Could not build ANPS",(nextChan));}
        }
 
        if (fCount[channel] != 0) {
            fNoiseAvgPowerSpectrum[channel] /= fCount[channel];
           
            // set number of events and source runs of QAverageVector
            
            fNoiseAvgPowerSpectrum[channel].SetSourceRuns(fSRuns); 
            fNoiseAvgPowerSpectrum[channel].SetNumEvents(fCount[channel]); 
     
            // Fill Output
 
            QAverageNoiseHandle npsHandle(channel); 
            if(fValidityKind=="dataset"){
              GlobalHandle<QInt> dHandle("Dataset");
              GlobalData().Get("",&dHandle,"");
              int dataset = dHandle.Get();
              npsHandle.SetDataset(dataset);
            }
            else if (fValidityKind=="run"){
              if (fSRuns.size()==1){
                npsHandle.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.");
 
            npsHandle.Set(fNoiseAvgPowerSpectrum[channel]);
            GlobalData().Set(&npsHandle,fOutput);
            
        }
        else {
                Error("Could not build average nps for channel %i: , 0 events passing the selection cuts",channel);
        }
        iter++;
    }
}