RdChannelNoisePulseCounter.cc

Go to the documentation of this file.

#include "RdChannelNoisePulseCounter.h"
#include <utl/ErrorLogger.h>

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

#include <utl/Point.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/RadioGeometryUtilities.h>
#include <fwk/CentralConfig.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <math.h>

#define INFO2(x,y) if ((x) <= fInfoLevel) INFO(y)

using namespace std;
using namespace fwk;


namespace RdChannelNoisePulseCounter {

RdChannelNoisePulseCounter::RdChannelNoisePulseCounter() :
  fInfoLevel(0),
  fSNRThreshold(10),
  fFirstEvent(true),
  fNumberOfEvents(0),
  fOutputPath("output"),
  fOutputFilename("noisecounter"),
  fMinimumStationId(100),
  fMaximumStationId(124),
  fMaximumNeighborDistance(200 * utl::meter)
{ }


VModule::ResultFlag
RdChannelNoisePulseCounter::Init()
{
  // Initialize your module here. This method
  // is called once at the beginning of the run.
  // The eSuccess flag indicates the method ended
  // successfully.  For other possible return types,
  // see the VModule documentation.

  INFO("RdChannelNoisePulseCounter::Init()");
  utl::Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelNoisePulseCounter");
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  topBranch.GetChild("SNRThreshold").GetData(fSNRThreshold);
  topBranch.GetChild("OutputPath").GetData(fOutputPath);
  topBranch.GetChild("OutputFilename").GetData(fOutputFilename);
  topBranch.GetChild("MinimumStationId").GetData(fMinimumStationId);
  topBranch.GetChild("MaximumStationId").GetData(fMaximumStationId);
  topBranch.GetChild("MaximumNeighborDistance").GetData(fMaximumNeighborDistance);

  return eSuccess;
}


double
RdChannelNoisePulseCounter::FindPulse(revt::Station& station)
{
  double NoiseSearchWindowStart = station.GetRecData().GetParameter(revt::eNoiseWindowStart);
  double NoiseSearchWindowStop = station.GetRecData().GetParameter(revt::eNoiseWindowStop);
  double SignalSearchWindowStart = station.GetRecData().GetParameter(revt::eSignalSearchWindowStart);
  double SignalSearchWindowStop = station.GetRecData().GetParameter(revt::eSignalSearchWindowStop);

  stringstream info;
  // identify active channels
  vector<revt::Station::ConstChannelIterator> usableChannels;
  for (revt::Station::ConstChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd();
      ++cIt)
    if (cIt->IsActive())
      usableChannels.push_back(cIt);

  // remove station if it has no channels (no data) - should actually be done beforehand, but just to be sure
  if (usableChannels.size() == 0) {
    info.str("");
    info << "Station ";
    info << station.GetId();
    info
        << " does not have any channels, removing it! Should treat this case beforehand, e.g. in RdChannelSelector.";
    WARNING(info.str());
    //        rEvent.RemoveStation(sIt);
    return 0;
  }

  if (usableChannels.size() != 2) {
    if (usableChannels.size() == 3) {
      info.str("");
      info << "Station ";
      info << station.GetId();
      info
          << " does have three active channels! The PreWaveFitter does only support two active channel for the moment. The first two active channels will be used.";
      WARNING(info.str());
    } else {
      info.str("");
      info << "Station ";
      info << station.GetId();
      info
          << " does not have exactly two active channels! Have you forgotten to include the RdChannelSelector?";
      ERROR(info.str());
      throw utl::DoesNotComputeException(
          "RdPreWaveFitter does not support stations with number of Channels != 2 or != 3.");
    }
  }

  const revt::Channel& channel1 = *usableChannels.at(0);
  const revt::Channel& channel2 = *usableChannels.at(1);
  const revt::ChannelTimeSeries& channelTrace1 = channel1.GetChannelTimeSeries();
  const revt::ChannelTimeSeries& channelTrace2 = channel2.GetChannelTimeSeries();

  revt::ChannelTimeSeries channelTrace;
  channelTrace.SetBinning(channelTrace1.GetBinning());
  channelTrace.ClearTrace();

  for (unsigned int i = 0; i < channelTrace1.GetSize(); ++i) {
    channelTrace.PushBack(
        sqrt((channelTrace1[i] * channelTrace1[i]) + (channelTrace2[i] * channelTrace2[i])));
  }

  unsigned int start = NoiseSearchWindowStart / channelTrace1.GetBinning();
  unsigned int stop = NoiseSearchWindowStop / channelTrace1.GetBinning();
  unsigned int startSignal = SignalSearchWindowStart / channelTrace1.GetBinning();
  unsigned int stopSignal = SignalSearchWindowStop / channelTrace1.GetBinning();
  double rms = utl::TraceAlgorithm::RootMeanSquare(channelTrace, start, stop - 1);  // stop - 1 because RootMeanSquare function loops from i=start to i<=stop
  double amplitude = utl::TraceAlgorithm::Max(channelTrace, startSignal, stopSignal);
//  delete usableChannels;
  return amplitude * amplitude / (rms * rms);
}

VModule::ResultFlag RdChannelNoisePulseCounter::Run(evt::Event& event)
{


  INFO("RdChannelNoisePulseCounter::Run()");
  // Check if there are events at all
  if (!event.HasREvent()) {
    WARNING("RdChannelNoisePulseCounter::No radio event found!");
    return eContinueLoop;
  } else {
    revt::REvent& rEvent = event.GetREvent();
    // compute nearest neighbors
    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();
    const std::vector<int>& fullStationList = rDetector.GetFullStationList();
    std::map<int, std::vector<int> > mNeighbors;  // neighbors need to be recalculated each time as detector configuration might change
    for (std::vector<int>::const_iterator sIt = fullStationList.begin(); sIt != fullStationList.end(); ++sIt) {
      const int id = *sIt;
      if( (id > fMaximumNeighborDistance) or (id < fMinimumStationId)) {  // consider only specific detector setup
        continue;
      }
      if (!mNeighbors.count(id)) {
        std::vector<int> neighbors;
        mNeighbors[id] = neighbors;
      }
      utl::Point currentPosition = rDetector.GetStation(id).GetPosition();
      for (std::vector<int>::const_iterator sIt2 = fullStationList.begin(); sIt2 != fullStationList.end(); ++sIt2) {
        const int id2 = *sIt2;
        if(id2 == id) { // do not double count first station
          continue;
        }
        if( (id > fMaximumNeighborDistance) or (id < fMinimumStationId)) {  // consider only specific detector setup
          continue;
        }
        utl::Point currentPosition2 = rDetector.GetStation(id2).GetPosition();
        if ((currentPosition - currentPosition2).GetMag() < fMaximumNeighborDistance) {
          mNeighbors[id].push_back(id2);
        }
      }
    }
    stringstream info;
    for (std::map<int, std::vector<int> >::iterator it=mNeighbors.begin(); it!=mNeighbors.end(); ++it) {
      info.str("");
      info << "Station " << it->first << " has neighbors: ";
      for(std::vector<int>::iterator itNeighbor = it->second.begin(); itNeighbor!=it->second.end(); ++itNeighbor) {
        info << *itNeighbor << " ";
      }
      INFO2(eFinal, info.str());
    }

    fNumberOfEvents++;
    std::map<int, bool> pulseFoundMap;
    // loop through stations
    for (revt::REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
      revt::Station& station = *sIt;
      const int id = station.GetId();
      double SNR = FindPulse(station);
      if(SNR == 0) {
        continue;
      }
      double SignalSearchWindowStart = station.GetRecData().GetParameter(revt::eSignalSearchWindowStart);
      double SignalSearchWindowStop = station.GetRecData().GetParameter(revt::eSignalSearchWindowStop);
      int GPSseconds = rEvent.GetHeader().GetTime().GetGPSSecond();

      if (SNR >= fSNRThreshold) {
        fNoisePulseData.SNRs.push_back(SNR);
        fNoisePulseData.GPSsecods.push_back(GPSseconds);
        fNoisePulseData.numberOfPulses++;
        fStationNoisePulseData[id].SNRs.push_back(SNR);
        fStationNoisePulseData[id].GPSsecods.push_back(GPSseconds);
        fStationNoisePulseData[id].numberOfPulses++;
        if( (id > fMaximumNeighborDistance) or (id < fMinimumStationId)) {  // consider only specific detector setup
          pulseFoundMap[id] = true;
        }
      } else {
        if( (id > fMaximumNeighborDistance) or (id < fMinimumStationId)) {  // consider only specific detector setup
          pulseFoundMap[id] = false;
        }
      }
      fNoisePulseData.totalSignalSerachWindow += (SignalSearchWindowStop - SignalSearchWindowStart);
      fNoisePulseData.numberOfTraces++;
      fStationNoisePulseData[id].totalSignalSerachWindow += (SignalSearchWindowStop
          - SignalSearchWindowStart);
      fStationNoisePulseData[id].numberOfTraces++;
    }

    // search for coincident noise pulses in multiple stations
    std::map<int, int> coincidentPulses_tmp;
    for (std::map<int, bool>::iterator it=pulseFoundMap.begin(); it!=pulseFoundMap.end(); ++it) {
      if(it->second) {
        int id = it->first;
        int nCoincidences = 1;
        for(std::vector<int>::iterator itNeighbor = mNeighbors[id].begin(); itNeighbor!=mNeighbors[id].end(); ++itNeighbor) {
          if(pulseFoundMap[*itNeighbor]) {
            nCoincidences++;
          }
        }
        if (!coincidentPulses_tmp.count(nCoincidences)) {
          coincidentPulses_tmp[nCoincidences] = 1;
        } else {
          coincidentPulses_tmp[nCoincidences]++;
        }
      }
    }
    // if any n station coincidence has been found, increase member variable by one. Multiple n station coincidences in one event are not considered, beacuse of double counting problems.
    for (std::map<int, int>::iterator it = coincidentPulses_tmp.begin(); it != coincidentPulses_tmp.end(); ++it) {
      if (!fCoincidentPulses.count(it->first)) {
        fCoincidentPulses[it->first] = 1;
      } else {
        fCoincidentPulses[it->first]++;
      }
    }
  }
  return eSuccess;
}

VModule::ResultFlag RdChannelNoisePulseCounter::Finish()
{
  // Put any termination or cleanup code here.
  // This method is called once at the end of the run.
  INFO("RdChannelNoisePulseCounter::Finish()");

  std::ofstream ofs;
  std::ofstream ofs_SNR;
  std::ofstream ofs_GPS;

  std::string outputfilename = fOutputPath;
  if(outputfilename.substr(outputfilename.size()-1, 1) != "/") {
    outputfilename.append("/");
  }
  outputfilename.append(fOutputFilename);

  ofs_SNR.open ((outputfilename+"_SNRs.out").c_str(), std::ofstream::app);
  ofs_SNR.precision(5);
  for (std::vector<double>::iterator itSNR = fNoisePulseData.SNRs.begin(); itSNR!=fNoisePulseData.SNRs.end(); ++itSNR) {
      ofs_SNR << "\t" << *itSNR;
  }
  ofs_SNR.close();
  ofs_GPS.open ((outputfilename+"_GPSSeconds.out").c_str(), std::ofstream::app);
  for (std::vector<int>::iterator itGPS = fNoisePulseData.GPSsecods.begin(); itGPS!=fNoisePulseData.GPSsecods.end(); ++itGPS) {
      ofs_GPS << "\t" << *itGPS;
  }
  ofs_GPS.close();

  ofs.open ((outputfilename+".out").c_str(), std::ofstream::app);
  ofs << "# nTraces\tnPulseFounds\ttotalSignalSerachWindow\n";
  ofs << fNoisePulseData.numberOfTraces << "\t" << fNoisePulseData.numberOfPulses << "\t" << fNoisePulseData.totalSignalSerachWindow << endl;
  ofs.close();

  ofs.open ((outputfilename+"_stationwise.out").c_str(), std::ofstream::app);
  ofs_SNR.open ((outputfilename+"_SNRs_stationwise.out").c_str(), std::ofstream::app);
  ofs_GPS.open ((outputfilename+"_GPSSeconds_stationwise.out").c_str(), std::ofstream::app);
  ofs_SNR.precision(5);
  for (std::map<int, NoisePulseData>::iterator it=fStationNoisePulseData.begin(); it!=fStationNoisePulseData.end(); ++it) {
    ofs << it->first << "\t" << it->second.numberOfTraces << "\t" << it->second.numberOfPulses << "\t" << it->second.totalSignalSerachWindow << endl;
    ofs_SNR << it->first;
    for (std::vector<double>::iterator itSNR = it->second.SNRs.begin(); itSNR!=it->second.SNRs.end(); ++itSNR) {
      ofs_SNR << "\t" << *itSNR;
    }
    ofs_SNR << endl;

    ofs_GPS << it->first;
    for (std::vector<int>::iterator itSNR = it->second.GPSsecods.begin(); itSNR!=it->second.GPSsecods.end(); ++itSNR) {
      ofs_GPS << "\t" << *itSNR;
    }
    ofs_GPS << endl;
  }
  ofs_SNR.close();
  ofs_GPS.close();
  ofs.close();

  std::ofstream ofs_CoincidentPulses;
  ofs_CoincidentPulses.open((outputfilename+"_spatial_coincidences.out").c_str(), std::ofstream::app);
  for (std::map<int, int>::iterator it = fCoincidentPulses.begin(); it != fCoincidentPulses.end(); ++it) {
    ofs_CoincidentPulses << it->first << "\t" << it->second << "\t" << fNumberOfEvents << endl;
  }
  ofs_CoincidentPulses.close();

  stringstream info;
  info.str("");
  info << "number of processed traces is " << fNoisePulseData.numberOfPulses << "/" << fNoisePulseData.numberOfTraces;
  INFO2(eDebug, info.str());

  return eSuccess;
}

}