RdChannelNoiseGenerator.cc

Go to the documentation of this file.

#include "RdChannelNoiseGenerator.h"
#include <iostream>
#include <cmath>

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Noise.h>
#include <utl/NoiseWhite.h>
#include <utl/NoiseSaS.h>
#include <utl/NoiseCCIR670.h>
#include <utl/RandomEngine.h>
#include <utl/Reader.h>
#include <utl/StringCompare.h>
#include <utl/TraceAlgorithm.h>

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

#include <boost/random.hpp>
#include <boost/random/exponential_distribution.hpp>
#include <boost/random/uniform_real.hpp>

using namespace utl;
using namespace fwk;



namespace RdChannelNoiseGenerator {

  RdChannelNoiseGenerator::RdChannelNoiseGenerator() :
    fInfoLevel(0),
    fSaSSigmaOverride(0.),
    fSaSAlphaOverride(0.),
    fNoiseRMS(-1),
    randgen(time(nullptr)),
    interval_noise(4.75*micro*volt, 7.28*micro*volt),
    interval_h(0, 1),
    randomGenerator_sigma(randgen, interval_noise),
    randomGenerator_h(randgen, interval_h)
  { }


  VModule::ResultFlag
  RdChannelNoiseGenerator::Init()
  {
    INFO("RdChannelNoiseGenerator::Init()");

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

    RandomEngine* const prndEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);
    topBranch.GetChild("NoiseModel").GetData(fNoiseType);

    topBranch.GetChild("NoiseRMS").GetData(fNoiseRMS);

    if (StringEquivalent(fNoiseType, "WhiteNoise")) {
      double noiseTemp;
      topBranch.GetChild("WhiteNoiseTemperature").GetData(noiseTemp);
      fNoise = new NoiseWhite(noiseTemp, prndEngine);
    } else if (StringEquivalent(fNoiseType,"CCIR670")) {
      std::string atmosphereNoiseString;
      topBranch.GetChild("CCIR670AtmosphereProfile").GetData(atmosphereNoiseString);
      AtmosphereNoise atmosphereNoise = eNoAtmosphereNoise;
      if (StringEquivalent(atmosphereNoiseString,"None")) {
        atmosphereNoise = eNoAtmosphereNoise;
      } else if (StringEquivalent(atmosphereNoiseString,"Day")) {
        atmosphereNoise = eDayAtmosphereNoise;
      } else if (StringEquivalent(atmosphereNoiseString,"Night")) {
        atmosphereNoise = eNightAtmosphereNoise;
      }

      std::string industryProfileString;
      topBranch.GetChild("CCIR670IndustryProfile").GetData(industryProfileString);
      IndustryNoise industryProfile = eNoIndustryNoise;
      if (StringEquivalent(industryProfileString, "None")) {
        industryProfile = eNoIndustryNoise;
      } else if (StringEquivalent(industryProfileString, "Rural")) {
        industryProfile = eRuralIndustryNoise;
      } else if (StringEquivalent(industryProfileString, "Urban")) {
        industryProfile = eUrbanIndustryNoise;
      }

      fNoise = new NoiseCCIR670(AtmosphereNoise(atmosphereNoise), IndustryNoise(industryProfile), prndEngine);
    } else if (StringEquivalent(fNoiseType, "SaSNoise")) {
      fNoise_time_domain = new NoiseSaS(prndEngine);
      topBranch.GetChild("SaSSigmaOverride").GetData(fSaSSigmaOverride);
      topBranch.GetChild("SaSAlphaOverride").GetData(fSaSAlphaOverride);
    } else {
      WARNING("Illegal noise type defined!");
      return eFailure;
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelNoiseGenerator::Run(evt::Event& event)
  {
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }
    INFO("Run()");
    std::cout << fNoiseRMS << "   " << fNoiseType << std::endl;

    revt::REvent& rEvent = event.GetREvent();
    det::Detector& zedet = det::Detector::GetInstance();

    if (fNoiseRMS >= 0 && (fNoiseType == "WhiteNoise" || fNoiseType == "CCIR670")) {
      INFO("adding noise");
      for (revt::REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
        revt::Station& station = *sIt;
        for (revt::Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {
          revt::Channel& channel = *cIt;
          revt::ChannelFFTDataContainer noiseChannelData  = channel.GetFFTDataContainer();
          revt::ChannelFrequencySpectrum& noiseSpec = noiseChannelData.GetFrequencySpectrum();
          const rdet::Channel& channelDet=zedet.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId());
          // generate noise in frequency domain
          for (revt::ChannelFrequencySpectrum::SizeType bin = 0; bin < noiseSpec.GetSize(); ++bin) {
            double binFrequency = channel.GetFrequencyOfBin(bin);
            double antennaHeight = channelDet.GetIntegratedEffectiveAntennaHeight(binFrequency);
            noiseSpec[bin] = antennaHeight*fNoise->GetSpectralFieldNoiseAtFrequency(binFrequency, noiseSpec.GetBinning());
          }
          // normalize noise trace to desired RMS in time domain
          revt::ChannelTimeSeries& noiseTrace = noiseChannelData.GetTimeSeries();
          const double RMS = utl::TraceAlgorithm::RootMeanSquare(noiseTrace, 0, noiseTrace.GetSize() - 1);
          revt::ChannelTimeSeries& trace = channel.GetChannelTimeSeries();
          for (revt::ChannelTimeSeries::SizeType bin = 0; bin < trace.GetSize(); ++bin) {
            trace[bin] += noiseTrace[bin] * fNoiseRMS / RMS;
          }
        }
      }
    } else if (fNoiseType == "WhiteNoise" || fNoiseType == "CCIR670"){
      // loop through all stations
      for (revt::REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
        revt::Station& station = *sIt;
        for (revt::Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {
          revt::Channel& channel = *cIt;
          const rdet::Channel& channelDet=zedet.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId());
          revt::ChannelFrequencySpectrum& mySpec = channel.GetChannelFrequencySpectrum();
          for (revt::ChannelFrequencySpectrum::SizeType bin = 0; bin < mySpec.GetSize(); ++bin) {
            double binFrequency = channel.GetFrequencyOfBin(bin);
            double antennaHeight = channelDet.GetIntegratedEffectiveAntennaHeight(binFrequency);
            mySpec[bin] += antennaHeight*fNoise->GetSpectralFieldNoiseAtFrequency(binFrequency, mySpec.GetBinning());
          }

        }

      }

    } else if (fNoiseType == "SaSNoise") {
      // SaS noise see GAP-2011-017
      evt::Event fNoiseEvent(event);
      for (revt::REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
        revt::Station& station = *sIt;
        for (revt::Station::ChannelIterator cIt = station.ChannelsBegin();cIt != station.ChannelsEnd(); ++cIt) {
          revt::Channel& channel = *cIt;
          revt::Channel& noise_channel = fNoiseEvent.GetREvent().GetStation(station.GetId()).GetChannel(channel.GetId());
          double sigma = fSaSSigmaOverride ? fSaSSigmaOverride : randomGenerator_sigma();
          double alpha = fSaSAlphaOverride ? fSaSAlphaOverride : (1./91)*log(randomGenerator_h())+1.999;
          revt::ChannelTimeSeries& noiseTime = noise_channel.GetChannelTimeSeries();
          double cof0 = pow(1 - pow(0.29, 2), -0.5);
          for (revt::ChannelTimeSeries::SizeType bin = 0; bin < noiseTime.GetSize(); ++bin) {
            if (bin == 0) {
              noiseTime[bin] = sigma*fNoise_time_domain->CreateSaSNoiseInTimeDomain(alpha);
            } else {
              noiseTime[bin] = (sigma*fNoise_time_domain->CreateSaSNoiseInTimeDomain(alpha)+ 0.29*cof0*noiseTime[bin-1])/cof0;
            }
          }
          const rdet::Channel& channelDet=zedet.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId());
          revt::ChannelFrequencySpectrum& mySpec = channel.GetChannelFrequencySpectrum();
          const revt::ChannelFrequencySpectrum& noiseSpec = noise_channel.GetChannelFrequencySpectrum();
          for (revt::ChannelFrequencySpectrum::SizeType bin = 0; bin < mySpec.GetSize(); ++bin) {
            double binFrequency = channel.GetFrequencyOfBin(bin);
            double antennaHeight = channelDet.GetIntegratedEffectiveAntennaHeight(binFrequency);
            mySpec[bin] += antennaHeight*noiseSpec[bin];
          }
        }
      }
    }
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelNoiseGenerator::Finish()
  {
    INFO("RdChannelNoiseGenerator::Finish()");
    return eSuccess;
  }

}