RdChannelNoisePowerAnalyser.cc

Go to the documentation of this file.

#include "RdChannelNoisePowerAnalyser.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 <revt/ChannelRRecDataQuantities.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 <utl/TimeStamp.h>
#include <utl/UTCDate.h>
#include <fwk/CentralConfig.h>

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

#include <math.h>

#include "TProfile.h"
#include "TProfile2D.h"
#include "TFile.h"
#include "TH2D.h"
#include "TMath.h"
#include "TCanvas.h"

#include <utl/LeapSeconds.h>

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

using namespace std;

namespace RdChannelNoisePowerAnalyser {

// function to calculate LST
double RdChannelNoisePowerAnalyser::GPSSecToLSTHourAuger(const unsigned long gps)
{
  const double J2000 = 2451545.0;
  const double Dtohr = 0.066666666666666;
  const double auger_longitude = -69.6;

  // time_t utc = gps + 315964819 + 16;  // GPS UTC offset, Leap Seconds have to be updated
  time_t utc = 0;
  utl::LeapSeconds::GetInstance().ConvertGPSToUnix(gps, utc);

  const tm* const t = gmtime(&utc);
  int year = t->tm_year + 1900;
  int month = t->tm_mon + 1;
  const int day = t->tm_mday;
  const int hour = t->tm_hour;
  const int min = t->tm_min;
  const int sec = t->tm_sec;

#warning can you use the utl::ModifiedJulianDate here?
  //Date to Julian Date
  if (month <= 2) {
    --year;
    month += 12;
  }
  const int a = int(floor(year / 100.));
  const int b = 2 - a + a / 4;
  double jd = floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1)) + day + b - 1524.5;
  const double fd = double(60 * 60 * hour + 60 * min + sec) / 86400;
  jd = jd + fd;

  static double sidereal_coeff[4] = { 280.46061837, 360.98564736629, 0.000387933, 38710000 };
  const double dt = jd - J2000;
  const double dt0 = dt / 36525.;
  const double theta = sidereal_coeff[0] + dt * sidereal_coeff[1]
      + dt0 * dt0 * (sidereal_coeff[2] - dt0 / sidereal_coeff[3]);  // is in degrees
      // original J. Meeus: longitude is positive westwards  *lst = mod(theta-longitude,360.) * kSTC::Dtohr;
      // with the new convention for longitudes (positive eastwards) :
  const double lst = fmod(theta + (360 + auger_longitude), 360.) * Dtohr;
  return lst;
}

RdChannelNoisePowerAnalyser::RdChannelNoisePowerAnalyser() :
    fInfoLevel(0),
    fUnbinnedASCIIOutput(false),
    fFrequencyBinWidth(1 * utl::MHz),
    fMinFrequency(30 * utl::MHz),
    fMaxFrequency(80 * utl::MHz),
    fOutputPath("output"),
    fOutputFilename("noise_amplitudes"),
    fMinimumStationId(101),
    fMaximumStationId(224),
    fMinimumChannelId(1),
    fMaximumChannelId(2),
    fLocalSiderealDay(86164.090530833 * utl::second),  // exact definition from Aoki, S., B. Guinot, G. H. Kaplan, H. Kinoshita, D. D. McCarthy and P. K. Seidelmann: "The new definition of Universal Time". Astronomy and Astrophysics 105(2), 361, 1982, see equation 19.
    fNumberOfBinsPerDay(144),
    fNFrequencyBins(0),
    fReferenceTime(utl::UTCDateTime(2011, 1, 1, 0, 0, 0, 0.)),  // initialize reference time to 1. Jan. 2011
    fRootOut(0),
    fCurrentMonth(0)
{
}

fwk::VModule::ResultFlag RdChannelNoisePowerAnalyser::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("RdChannelNoisePowerAnalyser::Init()");
  utl::Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch(
      "RdChannelNoisePowerAnalyser");
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  topBranch.GetChild("OutputPath").GetData(fOutputPath);
  topBranch.GetChild("OutputFilename").GetData(fOutputFilename);
  topBranch.GetChild("minStationId").GetData(fMinimumStationId);
  topBranch.GetChild("maxStationId").GetData(fMaximumStationId);

  // initialize data structure
  std::string outputfilename = fOutputPath;
  if (outputfilename.substr(outputfilename.size() - 1, 1) != "/") {
    outputfilename.append("/");
  }
  outputfilename.append(fOutputFilename);
  std::string rootOutputFilename = outputfilename + ".root";
  fRootOut = new TFile(rootOutputFilename.c_str(), "RECREATE");

  fNFrequencyBins = (fMaxFrequency - fMinFrequency) / fFrequencyBinWidth;
  fFrequencies.reserve(fNFrequencyBins);
  for (int i = 0; i <= fNFrequencyBins; ++i) {
    fFrequencies.push_back(fMinFrequency + fFrequencyBinWidth * i);
  }
  cout << "creating data structure...";
  for (int iStation = fMinimumStationId; iStation <= fMaximumStationId; ++iStation) {
    cout << "." << std::flush;
    for (int iChannel = fMinimumChannelId; iChannel <= fMaximumChannelId; ++iChannel) {
      for (int iFrequency = 0; iFrequency < fNFrequencyBins; ++iFrequency) {
        Key key(iStation, iChannel, iFrequency);

        // create profile histograms
        stringstream title;
        title << "station " << boost::get<0>(key) << ", channel " << boost::get<1>(key)
              << ", frequency " << fFrequencies[boost::get<2>(key)] / utl::MHz << "MHz";
        stringstream sstr;
//        sstr.str("");
//        sstr << "Amp_" << boost::get<0>(key) << "_" << boost::get<1>(key) << "_"
//             << boost::get<2>(key);
//        TProfile *hAmp = new TProfile(sstr.str().c_str(), title.str().c_str(), fNumberOfBinsPerDay,
//                                      0, 24);
//        hAmp->SetDirectory(fRootOut);
//        hAmp->GetYaxis()->SetTitle("Noise Amplitude [V/MHz]");
//        hAmp->GetXaxis()->SetTitle("LST [h]");
//        fProfileAmplitudeMap.insert(std::pair<Key, TProfile*>(key, hAmp));
//        sstr.str("");
//        sstr << "Power_" << boost::get<0>(key) << "_" << boost::get<1>(key) << "_"
//             << boost::get<2>(key);
//        TProfile *hPower = new TProfile(sstr.str().c_str(), title.str().c_str(),
//                                        fNumberOfBinsPerDay, 0, 24);
//        hPower->SetDirectory(fRootOut);
//        hPower->GetYaxis()->SetTitle("power spectral density [W/MHz]");
//        hPower->GetXaxis()->SetTitle("LST [h]");
//        fProfilePowerMap.insert(std::pair<Key, TProfile*>(key, hPower));

        // initialize 2D histograms
        sstr.str("");
        sstr << "Power2DSparse_" << boost::get<0>(key) << "_" << boost::get<1>(key) << "_"
             << boost::get<2>(key);

        Int_t bins[2] = { fNumberOfBinsPerDay, 400 };
        Double_t xmin[2] = { 0, -14 };
        Double_t xmax[2] = { 24, -6 };
        THnSparseF* hs = new THnSparseF(sstr.str().c_str(), title.str().c_str(), 2, bins, xmin,
                                        xmax);
//        hs->SetDirectory(fRootOut);
        hs->GetAxis(0)->SetTitle("LST [h]");
        hs->GetAxis(1)->SetTitle("Log (power spectral density [W/MHz])");
        fHistogramPowerMap.insert(std::pair<Key, THnSparseF*>(key, hs));



        // initialize 2D histograms linear
        sstr.str("");
        sstr << "Power2DSparseLinear_" << boost::get<0>(key) << "_" << boost::get<1>(key) << "_"
             << boost::get<2>(key);

        Int_t bins2[2] = { fNumberOfBinsPerDay, 1000 };
        Double_t xmin2[2] = { 0, 1e-12 };
        Double_t xmax2[2] = { 24, 50e-12 };
        THnSparseF* hs2 = new THnSparseF(sstr.str().c_str(), title.str().c_str(), 2, bins2, xmin2,
                                        xmax2);
        hs2->GetAxis(0)->SetTitle("LST [h]");
        hs2->GetAxis(1)->SetTitle("power spectral density [W/MHz]");
        fHistogramPowerMapLinear.insert(std::pair<Key, THnSparseF*>(key, hs2));
//        TH2D *hHist2dPower = new TH2D(sstr.str().c_str(), title.str().c_str(), fNumberOfBinsPerDay,
//                                      0, fLocalSiderealDay, 500, -7, 0);
//        hHist2dPower->SetDirectory(fRootOut);
//        hHist2dPower->GetYaxis()->SetTitle("Log (Noise Power [V^2/MHz])");
//        hHist2dPower->GetXaxis()->SetTitle("(time - 2011/01/01) mod LSD [ns]");
//        fHistogramPowerMap.insert(std::pair<Key, TH2D*>(key, hHist2dPower));

      }
    }
  }
  cout << endl;

  return eSuccess;
}

fwk::VModule::ResultFlag RdChannelNoisePowerAnalyser::Run(evt::Event& event)
{

  INFO2(eDebug, "RdChannelNoisePowerAnalyser::Run()");
  // Check if there are events at all
  if (!event.HasREvent()) {
    WARNING("RdChannelNoisePowerAnalyser::No radio event found!");
    return eContinueLoop;
  }
  revt::REvent& revent = event.GetREvent();
  // compute nearest neighbors
//  const det::Detector& detector = det::Detector::GetInstance();
//  const rdet::RDetector& rDetector = detector.GetRDetector();

  utl::TimeStamp currentTimeStamp = revent.GetHeader().GetTime();

  // initialize dynamic spectra
  const utl::UTCDateTime date = utl::UTCDateTime(currentTimeStamp);
  const int year = date.GetYear();
  const int month = date.GetMonth();
  if (fCurrentMonth != month) {
    fCurrentMonth = month;
    // initialize dynamic spectrum histogram
    stringstream sstr;
    sstr << "initializing dynamic power spectrum for month " << month << " and year " << year;
    INFO(sstr.str());
    for (int iStation = fMinimumStationId; iStation <= fMaximumStationId; ++iStation) {
      for (int iChannel = fMinimumChannelId; iChannel <= fMaximumChannelId; ++iChannel) {
        stringstream sstr;
        sstr.str("");
        KeyYearMonth key(iStation, iChannel, year, month);
        sstr << "DynamicSpectrum_" << boost::get<0>(key) << "_" << boost::get<1>(key) << "_"
             << boost::get<2>(key) << "_" << boost::get<3>(key);
        stringstream title;
        title << "station " << boost::get<0>(key) << ", channel " << boost::get<1>(key) << ", "
              << month << " " << year;
        utl::UTCDate tmpDateStart = utl::UTCDate(year, month, 1);
        int nextMonth = month + 1;
        int nextYear = year;
        if (nextMonth > 12) {
          nextMonth = 1;
          nextYear += 1;
        }
        utl::UTCDate tmpDateStop = utl::UTCDate(nextYear, nextMonth, 1);
        const int numberOfBins = NumberOfDaysInMonth(year, month) * 24 * 6;  // every ten minutes
        stringstream info;
        info.str("");
        info << "number of days per month = " << NumberOfDaysInMonth(year, month);
        INFO(info.str());
        TProfile2D *hHist2dPower = new TProfile2D(
            sstr.str().c_str(), title.str().c_str(), numberOfBins,
            utl::UTCDateTime(tmpDateStart).GetTimeStamp().GetGPSSecond(),
            utl::UTCDateTime(tmpDateStop).GetTimeStamp().GetGPSSecond(), fNFrequencyBins,
            fMinFrequency, fMaxFrequency);
        cout << "start time " << utl::UTCDateTime(tmpDateStart).GetTimeStamp().GetGPSSecond() << endl;
        cout << "stop time " << utl::UTCDateTime(tmpDateStop).GetTimeStamp().GetGPSSecond() << endl;

        hHist2dPower->SetDirectory(fRootOut);
        hHist2dPower->GetYaxis()->SetTitle("frequency [GHz]");
        hHist2dPower->GetXaxis()->SetTitle("GPS seconds");
        hHist2dPower->GetZaxis()->SetTitle("log(power spectral density [W/MHz])");
        fDynamicPowerSpectrum.insert(std::pair<KeyYearMonth, TProfile2D*>(key, hHist2dPower));
        info.str("");
        info << "initializing dynamic power spectrum for month " << month << " and year " << year
             << "station " << iStation << " and channel " << iChannel;
        INFO(info.str());
      }
    }
  }
  // determine correct LSD bin
  const unsigned long currentGPSSeconds = revent.GetHeader().GetTime().GetGPSSecond();
  // const double deltaLST = fmod(currentTimeStamp - fReferenceTime.GetTimeStamp(), fLocalSiderealDay);
  const double deltaLST = GPSSecToLSTHourAuger(currentGPSSeconds);
  std::vector<double> LSDBinBorders;
  LSDBinBorders.reserve(fNumberOfBinsPerDay + 1);
  for (int i = 0; i <= fNumberOfBinsPerDay; ++i) {
    LSDBinBorders[i] = i * 24. / (fNumberOfBinsPerDay);
  }
  int deltaLSTBin = -1;
  for (int i = 0; i < fNumberOfBinsPerDay; ++i) {
    if (deltaLST >= LSDBinBorders[i] and deltaLST < LSDBinBorders[i + 1]) {
      deltaLSTBin = i;
    }
  }
  if (deltaLSTBin == -1) {
    cout << "delta LSTBin == -1, deltaLST = " << deltaLST << endl;
  }

  std::map<int, std::vector<std::vector<double> > > noiseMap;  // keys are [stationId][channelId][frequency] amplitudes
  for (revt::REvent::ConstCandidateStationIterator iS = revent.CandidateStationsBegin();
      iS != revent.CandidateStationsEnd(); ++iS) {
    // skip stations that are out of range
    if (iS->GetId() < fMinimumStationId or iS->GetId() > fMaximumStationId) {
      stringstream info;
      info << "station " << iS->GetId() << " is out of range, skipping station.";
      INFO2(eDebug, info.str());
      continue;
    }
    std::vector<std::vector<double> > tmpChannels;
    for (revt::Station::ConstChannelIterator iC = iS->ChannelsBegin(); iC != iS->ChannelsEnd();
        ++iC) {
      if (not iC->IsActive()) {
        stringstream info;
        info << "channel " << iC->GetId() << " is inactive. Skipping channel";
        INFO2(eDebug, info.str());
        continue;
      }
      // skip channels that are out of range
      if (iC->GetId() < fMinimumChannelId or iC->GetId() > fMaximumChannelId) {
        stringstream info;
        info << "channel " << iC->GetId() << " is out of range, skipping channel";
        INFO2(eDebug, info.str());
        continue;
      }

      const double timeBinSize = iC->GetFFTDataContainer().GetConstTimeSeries().GetBinning();
      const double traceLength = iC->GetFFTDataContainer().GetConstTimeSeries().GetSize()
          * timeBinSize;

      const revt::ChannelFrequencySpectrum& spectrum = iC->GetConstChannelFrequencySpectrum();
      // double frequencyBinWidth = spectrum.GetBinning(); //unused, see below
      // double nBins1MHz = 1 * utl::MHz / frequencyBinWidth;
      // double nBins1MHz = 1. / frequencyBinWidth;  // not used in noiseAmplitude and noisePower calc anymore since 2016-12-07 glaser commit
      std::vector<double> noiseAmplitudes;
      noiseAmplitudes.reserve(fNFrequencyBins);

      // debug
//      cout << "timetrace " << endl;
//      const revt::ChannelTimeSeries& timeseries = iC->GetChannelTimeSeries();
//      int nn = iC->GetFFTDataContainer().GetConstTimeSeries().GetSize();
//      cout << "timeseries = [";
//      for (int i = 0; i < nn; ++i) {
//        cout << timeseries[i] << ", ";
//      }
//      cout << "]" << endl;
//      cout << "tracelength = " << traceLength << endl;
      // end debug

      KeyYearMonth keyYearMonth(iS->GetId(), iC->GetId(), year, month);
      for (int i = 0; i < fNFrequencyBins; ++i) {
        Key key(iS->GetId(), iC->GetId(), i);
        int binStart = iC->GetFFTDataContainer().GetBinOfFrequency(fFrequencies[i]);
        int binStop = iC->GetFFTDataContainer().GetBinOfFrequency(fFrequencies[i + 1]);
        double frequencyInterval = iC->GetFFTDataContainer().GetFrequencyOfBin(binStop)
            - iC->GetFFTDataContainer().GetFrequencyOfBin(binStart);
        double noiseAmplitude(0.);
        double noisePower(0.);
        // average over 1 MHz interval
        for (int j = binStart; j < binStop; ++j) {
          noiseAmplitude += abs(spectrum[j]);
          noisePower += abs(spectrum[j]) * abs(spectrum[j]);
        }
        // normalize to number of bins (because it can be different for different frequency ranges as the samplingfrequency is not a multiple of our frequency bin width
        noiseAmplitude *= 1. / (binStop - binStart);  // average over 1 MHz interval
        noisePower *= 1. / (binStop - binStart);
        noisePower *= 1. / frequencyInterval * (fFrequencies[i + 1] - fFrequencies[i]);  // average over 1 MHz interval

        noisePower *= 2;  // because negative frequencies are omitted in real fft
        noisePower /= traceLength;  // energy -> power
        noisePower /= (50 * utl::ohm);  // V^2 -> Watts (50 ohm impedance system)
        noisePower /= (utl::watt / utl::MHz);  // convert to units W/MHz

//        stringstream sstr;
//        sstr << "f = " << fFrequencies[i] / utl::MHz << " - " << fFrequencies[i + 1] / utl::MHz
//             << " MHz" << " delta = " << frequencyInterval / utl::MHz << " noisepower = "
//             << noisePower;
//        INFO2(eDebug, sstr.str());
//        // add info to profiles
//        fProfileAmplitudeMap[key]->Fill(deltaLST, noiseAmplitude);
//        fProfilePowerMap[key]->Fill(deltaLST, noisePower);
        double x[2] = { deltaLST, TMath::Log10(noisePower) };
        fHistogramPowerMap[key]->Fill(x);

        double x2[2] = { deltaLST, noisePower};
        fHistogramPowerMapLinear[key]->Fill(x2);

        const double currentFrequency = 0.5 * (fFrequencies[i] + fFrequencies[i + 1]);
        fDynamicPowerSpectrum[keyYearMonth]->Fill(currentGPSSeconds, currentFrequency,
                                                  TMath::Log10(noisePower));
      }
      tmpChannels.push_back(noiseAmplitudes);
    }
  }

  return eSuccess;
}

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

  fRootOut->Write();

  typedef std::map<KeyYearMonth, TProfile2D*>::iterator it_type;
  for(it_type iterator = fDynamicPowerSpectrum.begin(); iterator != fDynamicPowerSpectrum.end(); iterator++) {
    iterator->second->Write();
    /*
    cout << "writing power spectrum" << endl;
    TCanvas *c1 = new TCanvas();
    c1->cd();
    iterator->second->Draw("COLZ");
    stringstream sstr;
    sstr << iterator->second->GetName() << ".png";
    c1->SaveAs(sstr.str().c_str());
    */

      // iterator->first = key
      // iterator->second = value
      // Repeat if you also want to iterate through the second map.
  }

  // delete all histograms from memory
//  for (int iStation = fMinimumStationId; iStation <= fMaximumStationId; ++iStation) {
//    for (int iChannel = fMinimumChannelId; iChannel <= fMaximumChannelId; ++iChannel) {
//      for (int iFrequency = 0; iFrequency < fNFrequencyBins; ++iFrequency) {
//        Key key(iStation, iChannel, iFrequency);
//        delete fProfileAmplitudeMap[key];
//        delete fProfilePowerMap[key];
//      }
//    }
//  }
  for (int iStation = fMinimumStationId; iStation <= fMaximumStationId; ++iStation) {
    for (int iChannel = fMinimumChannelId; iChannel <= fMaximumChannelId; ++iChannel) {
      for (int iFrequency = 0; iFrequency < fNFrequencyBins; ++iFrequency) {
        Key key(iStation, iChannel, iFrequency);
        fHistogramPowerMap[key]->Write();
        delete fHistogramPowerMap[key];
        fHistogramPowerMapLinear[key]->Write();
        delete fHistogramPowerMapLinear[key];
//        TH2D* h2 = fHistogramPowerMap[key]->Projection(1, 0);
//        stringstream sstr;
//        sstr << "Power2D_" << boost::get<0>(key) << "_" << boost::get<1>(key) << "_"
//                     << boost::get<2>(key);
//        h2->SetName(sstr.str().c_str());
//        h2->Write();
//        delete h2;
      }
    }
  }

  return eSuccess;
}

}