RdChannelDebugWriter.cc

Go to the documentation of this file.

#include "RdChannelDebugWriter.h"

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>

#include <TROOT.h>
#include <TH1F.h>
#include <TFile.h>


using namespace fwk;
using namespace utl;
using namespace revt;
using std::vector;
using std::set;
using std::string;
using std::ostringstream;
using std::ofstream;
using std::endl;
using std::ios;


namespace RdChannelDebugWriter {

  VModule::ResultFlag
  RdChannelDebugWriter::Init()
  {
    // Read in the configurations of the xml file
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelDebugWriter");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    // read in the station list as vector, and then copy it to the set class variable
    const vector<int> dummyStationIDlist = topBranch.GetChild("StationIds").Get<vector<int>>();
    fStationIds = set<int>(dummyStationIDlist.begin(), dummyStationIDlist.end());

    // read in the channel list as vector, and then copy it to the set class variable
    const vector<int> dummyChannelIDlist = topBranch.GetChild("ChannelIds").Get<vector<int>>();
    fChannelIds = set<int>(dummyChannelIDlist.begin(), dummyChannelIDlist.end());

    topBranch.GetChild("EventFraction").GetData(fEventFraction);
    topBranch.GetChild("AsciiSpectrumNameBase").GetData(fAsciiSpectrumNameBase);
    topBranch.GetChild("RootFileNameBase").GetData(fRootFileNameBase);
    topBranch.GetChild("AsciiSpectrumNameBase").GetData(fAsciiSpectrumNameBase);
    topBranch.GetChild("AsciiTraceNameBase").GetData(fAsciiTraceNameBase);

    const string unitstr = topBranch.GetChild("UseUnits").Get<string>();
    fUseUnits = (unitstr == "yes");

    ostringstream info;
    info << "fRootFileNameBase = " << fRootFileNameBase;
    INFODebug(info);

    if (!fUseUnits)
      INFOFinal("The Module will not convert to physical unit\n"
                 "  Be sure of what you are doing");

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelDebugWriter::Run(evt::Event& event)
  {
    // Check if there are events at all
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    const REvent& rEvent = event.GetREvent();

    // Check if it is the first call for this event
    if (fLastEventID != rEvent.GetHeader().GetId()) {
      ++fEventCounter;
      fStepCounter = 0;
      fLastEventID = rEvent.GetHeader().GetId();
      fLastEventGPSsecond = rEvent.GetHeader().GetTime().GetGPSSecond();
    } else {
      ++fStepCounter;
    }

    ostringstream info;
    info << "Running RdChannelDebugWriter for event number "
         << fEventCounter << " (ID = " << fLastEventID << ") - step "
         << fStepCounter << " with GPS second " << fLastEventGPSsecond;
    INFOIntermediate(info);

    if (fEventFraction > 0)
      if (((fEventCounter - 1) % fEventFraction) != 0) {
        info.str("");
        info << "Skip event, because only every " << fEventFraction
             << "th event is considered for debug output (as set in XML).";
        INFOFinal(info);
        return eSuccess;
      }

    // write debug output into root file, if requested
    if (!fRootFileNameBase.empty())
      writeRootOutput(rEvent);

    // write spectra into ASCII file, if requested
    if (!fAsciiSpectrumNameBase.empty())
      writeAsciiSpectrum(rEvent);

    // write traces into ASCII file, if requested
    if (!fAsciiTraceNameBase.empty())
      writeAsciiTrace(rEvent);

    return eSuccess;
  }


  void
  RdChannelDebugWriter::writeRootOutput(const REvent& rEvent)
    const
  {
    ostringstream count;
    count << '_' << fLastEventGPSsecond << "_Step-" << fStepCounter;
    TFile OutputFile((fRootFileNameBase + count.str() + ".root").c_str(), "recreate");

    float unitmuv = 1;
    float unitmuvpermhz = 1;
    if (fUseUnits) {
      unitmuv = micro * volt;
      unitmuvpermhz = micro * volt / megahertz;
    }

    for (const auto& station : rEvent.StationsRange()) {

      // check if this station should be considered, otherwise continue
      if (!fStationIds.empty() && fStationIds.find(station.GetId()) == fStationIds.end())
        continue;

      ostringstream stat;
      stat << "Station_" << station.GetId();

      for (const auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        // check if this channel should be considered, otherwise continue
        if (!fChannelIds.empty() && fChannelIds.find(channel.GetId()) == fChannelIds.end())
          continue;

        // Work on the frequency spectrum and time series of this channel
        const auto& spectrum = channel.GetChannelFrequencySpectrum();
        const auto& trace = channel.GetChannelTimeSeries();

        ostringstream chan;
        chan << "_Ch." << channel.GetId();

        // Time Series histo
        const double tsSize = trace.GetSize();
        const double tsStart = trace.GetStart() / nanosecond;
        const double tsStop = (trace.GetSize() - 1) / nanosecond;
        TH1F HistTS(
          (stat.str() + chan.str() + "_TimeSeries" + count.str()).c_str(),
          (stat.str() + chan.str() + "_TimeSeries" + count.str()).c_str(),
          tsSize, tsStart, tsStop
        );

        /* Phase and absolute histos (Is taken into account the different conventions,
           i.e. the karlsrhue data where the Spectrum reaches from 80 MHz to 40 MHz
           instead of 0 MHz to 200 MHz) */
        const double binFreqSize  = spectrum.GetSize();
        const double binFreqBegin = std::min(
          channel.GetFrequencyOfBin(0),
          channel.GetFrequencyOfBin(spectrum.GetSize() - 1)
        ) / megahertz;
        const double binFreqEnd = std::max(
          channel.GetFrequencyOfBin(0),
          channel.GetFrequencyOfBin(spectrum.GetSize() - 1)
        ) / megahertz;

        TH1F HistPh(
          (stat.str() + chan.str() + "_Phase" + count.str()).c_str(),
          (stat.str() + chan.str() + "_Phase" + count.str()).c_str(),
          binFreqSize, binFreqBegin, binFreqEnd
        );
        TH1F HistAbs(
          (stat.str() + chan.str() + "_Abs" + count.str()).c_str(),
          (stat.str() + chan.str() + "_Abs" + count.str()).c_str(),
          binFreqSize, binFreqBegin, binFreqEnd
        );

        int w = 0;
        for (const auto& el : trace)
          HistTS.SetBinContent(w++, el / unitmuv);

        int k = 0;
        for (const auto& el : spectrum) {
          HistPh.SetBinContent(k, std::arg(el) / deg);
          HistAbs.SetBinContent(k++, abs(el) / unitmuvpermhz);
        }

        HistTS.SetXTitle("ns");
        HistTS.SetYTitle("microvolt");
        HistPh.SetXTitle("MHz");
        HistPh.SetYTitle("degrees");
        HistAbs.SetXTitle("MHz");
        HistAbs.SetYTitle("microvolt/MHz");
        HistTS.Write();
        HistPh.Write();
        HistAbs.Write();
      }
    }
    OutputFile.Close();
  }


  void
  RdChannelDebugWriter::writeAsciiSpectrum(const REvent& rEvent)
    const
  {
    float unitmuvpermhz = 1;
    if (fUseUnits)
      unitmuvpermhz = micro * volt / megahertz;

    for (const auto& station : rEvent.StationsRange()) {

      // check if this station should be considered, otherwise continue
      if (!fStationIds.empty() && fStationIds.find(station.GetId()) == fStationIds.end())
        continue;

      for (const auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        // check if this channel should be considered, otherwise continue
        if (!fChannelIds.empty() && fChannelIds.find(channel.GetId()) == fChannelIds.end())
          continue;

        const auto& spectrum = channel.GetChannelFrequencySpectrum();

        // write spectrum to file
        ostringstream filename;
        filename << fAsciiSpectrumNameBase
                 << '_' << fLastEventGPSsecond
                 << "_s" << channel.GetStationId()
                 << "_c" << channel.GetId()
                 << "_Step-" << fStepCounter
                 << ".dat";
        ofstream outfile(filename.str().c_str());

        // write the frequency axis, the amplitude, and the phase of the spectrum to the file
        for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
          outfile << channel.GetFrequencyOfBin(i) / megahertz
                  << ' ' << abs(spectrum[i]) / unitmuvpermhz
                  << ' ' << arg(spectrum[i]) / deg << '\n';
      }
    }
  }


  void
  RdChannelDebugWriter::writeAsciiTrace(const REvent& rEvent)
    const
  {
    float unitmuv = 1;
    if (fUseUnits)
      unitmuv = micro * volt;

    for (const auto& station : rEvent.StationsRange()) {

      // check if this station should be considered, otherwise continue
      if (!fStationIds.empty() && fStationIds.find(station.GetId()) == fStationIds.end())
        continue;

      for (const auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        // check if this channel should be considered, otherwise continue
        if (!fChannelIds.empty() && fChannelIds.find(channel.GetId()) == fChannelIds.end())
          continue;

        const auto& trace = channel.GetChannelTimeSeries();

        // write trace to file
        ostringstream filename;
        filename << fAsciiTraceNameBase
                 << '_' << fLastEventGPSsecond
                 << "_s" << channel.GetStationId()
                 << "_c" << channel.GetId()
                 << "_Step-" << fStepCounter
                 << ".dat";
        ofstream outfile(filename.str().c_str());

        // write the time axis and the amplitude to the file
        for (ChannelTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i)
          outfile << trace.GetBinning() * i / nanosecond
                  << ' ' << trace[i] / unitmuv << '\n';
      }
    }
  }

}