CalcBeaconRefPhase.cc

Go to the documentation of this file.

#include "CalcBeaconRefPhase.h"
#include <fwk/CentralConfig.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/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/AugerException.h>

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

#include <complex>
#include <cmath>
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>

#include <mysql/mysql.h>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace revt;


namespace CalcBeaconRefPhase {

CalcBeaconRefPhase::CalcBeaconRefPhase():
    fBeaconFrequencies(vector<double>()),
    fOutputRefPhaseFile(""),
    fDBServer(""),
    fReferenceStation(145),
    fReferenceChannel(1),
    fStationWithZeroPhase(101),
    fConvergenceLimit(0.001),
    fEventCounter(0),
    fMaxNumberEvents(true),
    fPhaseDifferences(std::map<double ,std::map<int ,vector<double> > >() )
{
}


CalcBeaconRefPhase::~CalcBeaconRefPhase()
{
}


VModule::ResultFlag
CalcBeaconRefPhase::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("CalcBeaconRefPhase::Init()");

    // Read in the configurations of the xml file
    Branch topBranch =
        CentralConfig::GetInstance()->GetTopBranch("CalcBeaconRefPhase");
    topBranch.GetChild("infoLevel").GetData(fInfoLevel);
    if (fInfoLevel < eNone || fInfoLevel > eDebug) {
        ERROR("<infoLevel> out of bounds");
        return eFailure;
    }

    topBranch.GetChild("BeaconFrequencies").GetData(fBeaconFrequencies);
    topBranch.GetChild("OutputRefPhaseFile").GetData(fOutputRefPhaseFile);
    topBranch.GetChild("DBServer").GetData(fDBServer);
    topBranch.GetChild("ReferenceStation").GetData(fReferenceStation);
    topBranch.GetChild("ReferenceChannel").GetData(fReferenceChannel);
    topBranch.GetChild("StationWithZeroPhase").GetData(fStationWithZeroPhase);
    topBranch.GetChild("MaxNumberEvents").GetData(fMaxNumberEvents);

    return eSuccess;
}


VModule::ResultFlag
CalcBeaconRefPhase::Run(evt::Event& event)
{
  stringstream message;
  message.str("CalcBeaconRefPhase::Run()");
  INFO(message.str());

  // check if beacon frequencies have to be read from database, or are provided in XML file
  // in the first case, fBeaconFrequencies is empty, i.e. size()==0  
  if(fBeaconFrequencies.size() < 1) {
    fBeaconFrequencies = det::Detector::GetInstance().GetRDetector().GetBeaconFrequencies();
  }

  
  // obtain phases of event
  std::map< double, std::map<int, double> > phases = map< double, map<int, double> >();


  // Check if there are events at all
  if(!event.HasREvent()) {
    WARNING("PhaseFinder::No radio event found!");
    return eContinueLoop;
  }
  
  REvent& rEvent = event.GetREvent();    

  // match time stamps of stations (if different)
  // this functionalty should be moved to its own module later
  matchStationTimeStamps(rEvent);
  
  
  // get phases of reference station
  map<double, double> refStationPhases = map<double, double> ();
  Channel& refChannel = rEvent.GetStation(fReferenceStation).GetChannel(fReferenceChannel);
  const revt::ChannelFrequencySpectrum refSpectrum = refChannel.GetChannelFrequencySpectrum();

  // check if spectrum has the correct lenght, such that the frequency bins are at the expected values
  // for which the beacon has been optimized
  if ( ((refSpectrum.GetSize()-1) % 1024) != 0 ) {
    message.str("");
    message << "Trace of reference station " << fReferenceStation << " must consist of 2048 samples or an integer multiple of it!";
    ERROR(message.str());
    return eFailure;
  }
  
  for (unsigned int beaconFreqN=0; beaconFreqN < fBeaconFrequencies.size(); ++beaconFreqN) {
    // check if frequency is inside the range (take care that the order of the refSpectrum is unknown)
    if( (fBeaconFrequencies[beaconFreqN] < min(refChannel.GetFrequencyOfBin(0),refChannel.GetFrequencyOfBin(refSpectrum.GetSize()-1)))
         || (fBeaconFrequencies[beaconFreqN] >  max(refChannel.GetFrequencyOfBin(0),refChannel.GetFrequencyOfBin(refSpectrum.GetSize()-1))) ) {
      WARNING("Beacon frequency must be inside of the spectrum!");
      return eContinueLoop;
    }

    // find nearest frequency bin:
    // bin = (f_beacon-f_0) * #lastBin / (f_1 - f_0)
    double bin = (fBeaconFrequencies[beaconFreqN]-refChannel.GetFrequencyOfBin(0)) * (refSpectrum.GetSize()-1)
                / (refChannel.GetFrequencyOfBin(refSpectrum.GetSize()-1) - refChannel.GetFrequencyOfBin(0));

    // get the phase
    refStationPhases[ fBeaconFrequencies[beaconFreqN] ] = arg(refSpectrum[round(bin)]);
  }  
 
 
  // loop through stations and for every station through every channel
  for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
        sIt != rEvent.CandidateStationsEnd(); ++sIt) {
    Station& station = *sIt;
    
    // check, if this station is already present in the event counter and increase counter
    if (fEventsPerStation.find(station.GetId())==fEventsPerStation.end()) {
      fEventsPerStation[station.GetId()]=1;
    } else {
      ++fEventsPerStation[station.GetId()];
    }

    Channel& channel = station.GetChannel(fReferenceChannel);

    if (fInfoLevel>=eObscure) {
      message.str("");
      message << "Calculating the phases of the beacon frequencies for channel "
              << channel.GetId()
              << " of station "
              << channel.GetStationId()
              << ".";
      INFO(message.str());
    }

    // Read the frequency spectrum of this channel
    const revt::ChannelFrequencySpectrum spectrum = channel.GetChannelFrequencySpectrum();

    if ( ((refSpectrum.GetSize()-1) % 1024) != 0 ) {
      message.str("");
      message << "Trace of station " << station.GetId() << " does not consist of 2048 samples or an integer multiple of it!"
              << "It has " << (spectrum.GetSize()-1)*2 << " samples.";
       WARNING(message.str());
      //return eFailure;
    }


    // get the amplitude spectrum for SNR calculation
    //Trace<double> amplitudeSpectrum;
    //for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
    //  amplitudeSpectrum.PushBack(abs(spectrum[i]));

    // loop through beacon frequencies and find corresponding frequency in the spectrum
    for (unsigned int beaconFreqN=0; beaconFreqN < fBeaconFrequencies.size(); ++beaconFreqN) {
      if (fInfoLevel>=eDebug) {
        message.str("");
        message << "Beacon frequency: "
                << fBeaconFrequencies[beaconFreqN]/megahertz
                << " MHz";
        INFO(message.str());

        message.str("");
        message << "Spectrum reaches from "
                << channel.GetFrequencyOfBin(0)/megahertz
                << " MHz to "
                << channel.GetFrequencyOfBin(spectrum.GetSize()-1)/megahertz
                << " MHz, binning = "
                << spectrum.GetBinning()/megahertz
                << " MHz.";
        INFO(message.str());
      }

      // check if frequency is inside the range (take care that the order of the spectrum is unknown)
      if( (fBeaconFrequencies[beaconFreqN] < min(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1)))
              || (fBeaconFrequencies[beaconFreqN] >  max(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1))) ) {
        WARNING("Beacon frequency must be inside of the spectrum!");
        return eContinueLoop;
      }

      // find nearest frequency bin:
      // bin = (f_beacon-f_0) * #lastBin / (f_1 - f_0)
      double bin = (fBeaconFrequencies[beaconFreqN]-channel.GetFrequencyOfBin(0)) * (spectrum.GetSize()-1)
                  / (channel.GetFrequencyOfBin(spectrum.GetSize()-1) - channel.GetFrequencyOfBin(0));

      // get the phase
      double phase = arg(spectrum[round(bin)]);
      
      if (fInfoLevel>=eDebug) {
        message.str("");
        message << "Looking at frequency: "
                << channel.GetFrequencyOfBin(round(bin)) / megahertz
                << " MHz "
                << "Phase: " << phase;
        INFO(message.str());
      }
      
      // store phase
      fPhaseDifferences[fBeaconFrequencies[beaconFreqN]][station.GetId()].push_back(phase - refStationPhases[fBeaconFrequencies[beaconFreqN]]);
    }
  }  
    
  // check if maximum number of events was reached
  ++fEventCounter;
  if (fEventCounter >= fMaxNumberEvents) {
    return eBreakLoop;
  }

  return eSuccess;
}


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

  INFO("CalcBeaconRefPhase::Finish()");
  INFO("");
  INFO("Calculating reference phases...");
  INFO("");
    
  message.str("");
  message << "\nPlease check manually, if the beacon frequency IDs are correct (i.e. the same as in the database):\n";
  for (unsigned int beaconFreqN = 0; beaconFreqN < fBeaconFrequencies.size(); ++beaconFreqN)
    message << beaconFreqN+1 << ": " << fBeaconFrequencies[beaconFreqN] << "\n";
  message << endl;
  WARNING(message.str());

  INFO("");
  
  // loop through beacon frequencies and calculate the reference phases
  for (unsigned int beaconFreqN=0; beaconFreqN < fBeaconFrequencies.size(); ++beaconFreqN) {
    // Calculate the reference phases iteratively by shifting the phases, until the mean becomes 0
    // the total shift then corresponds to the reference phase
    // In the map there are the reference phases and the corresponding RMS for each station
    map<int, pair<double, double> > referencePhases = calculateReferencePhases(fPhaseDifferences[fBeaconFrequencies[beaconFreqN]]);
   
    // shift all reference phases such that the desired station has a reference phase of zero
    if (referencePhases.find(fStationWithZeroPhase) == referencePhases.end()) {
      message.str("");
      message << "fStationWithZeroPhase (=" << fStationWithZeroPhase << ") not found in data!";
      ERROR(message.str());
      return eFailure;
    }

    double zeroPhaseStationPhase = referencePhases[fStationWithZeroPhase].first;

    for (map<int, pair<double, double> >::iterator it = referencePhases.begin(); it != referencePhases.end() ; ++it) {
      // substract phase of station which finally should have a reference phase of 0
      it->second.first -= zeroPhaseStationPhase;
      // confine phase between +/- pi
      while ( it->second.first > kPi ) {
        it->second.first -= kTwoPi;
      }
      while ( it->second.first < -kPi ) {
        it->second.first += kTwoPi;
      }

      // add all station IDs to the list of existing IDs at this occaision
      fExistingStationIDs.insert( it->first);
      //cout << "Result for station " << it->first << " \tRefPhase = " << it->second.first << " \tRMS = " << it->second.second << endl;
    }

    // store reference phases in global map
    fReferencePhases.push_back(referencePhases);
  }

  // Write values to file, for debug output
  if (fOutputRefPhaseFile != "") {
    writeReferencePhasesFile();
    writeFileForDatabase();
  }

  
  message.str("");
  message << "Station ID \t #events";
  for (vector<double>::iterator it=fBeaconFrequencies.begin(); it!=fBeaconFrequencies.end(); ++it) {
    message << " \t" << *it/megahertz << " MHz";
  }
  INFO(message.str());

  // print number of events per frequency and station
  for (set<unsigned int>::iterator it = fExistingStationIDs.begin(); it != fExistingStationIDs.end(); ++it) {
    message.str("");
    message << "   " << *it << "\t " << fEventsPerStation[*it];
    
    for (unsigned int beaconFreqN=0; beaconFreqN < fBeaconFrequencies.size(); ++beaconFreqN) {
      message << "\t " << fReferencePhases[beaconFreqN][*it].first << " +/- " << fReferencePhases[beaconFreqN][*it].second;
    }  
    INFO(message.str());
  }

  INFO("");

  
  
  return eSuccess;
}


void
CalcBeaconRefPhase::matchStationTimeStamps(revt::REvent& rEvent) const
{
  // store time stamp of the first station as reference
  double firstTimeStamp  = rEvent.CandidateStationsBegin()->GetRecData().GetParameter(eTraceStartTime);
  // loop through all stations:
  // shift traces if the time stamp does not match the one of the first station
  for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
          sIt != rEvent.CandidateStationsEnd(); ++sIt) {
    Station& station = *sIt;
    double stationTimeStamp = station.GetRecData().GetParameter(eTraceStartTime);
    cout << station.GetId() << ": \t" << stationTimeStamp << endl;
    if (stationTimeStamp != firstTimeStamp) {
      // correct for the difference and set new time stamp
      for (Station::ChannelIterator cIt = station.ChannelsBegin();
              cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;

        FFTDataContainerAlgorithm::ShiftTimeSeries(channel.GetFFTDataContainer(), stationTimeStamp-firstTimeStamp);
      }
      station.GetRecData().SetParameter(eTraceStartTime,firstTimeStamp);
    }
  }
}


map<int, pair<double, double> >
CalcBeaconRefPhase::calculateReferencePhases(std::map<int,std::vector<double> > phaseDifferences)
{
  map<int, pair<double, double> > referencePhases; // map to store results: reference phase and its RMS

  // loop through all stations
  for(map<int, vector<double> >::iterator it=phaseDifferences.begin(); it!=phaseDifferences.end(); ++it) {
    double mean = 0.;
    double totalShift = 0.;
    double rms = 0.;

    // shift all phase differences of each station, until the remaining mean is (almost) 0
    // the total shift corresponds to the reference phase
    // Taking directly the mean does not work, because the phase angles could extend over +/- pi.
    do {
      double phasesum = 0.;
      double rmssum = 0.;

      // shift phases by mean of last iterationstep
      // loop through all phase differences of the current station
      for(vector<double>::iterator vpos=it->second.begin(); vpos!=it->second.end(); ++vpos) {
        *vpos = *vpos - mean;
        if ( *vpos > kPi ) {
            *vpos -= kTwoPi;
        } else if ( *vpos < -kPi ) {
            *vpos += kTwoPi;
        }
        phasesum += *vpos;
        rmssum += (*vpos * *vpos);
      }
      mean = phasesum/double(it->second.size());
      rms = sqrt(rmssum/double(it->second.size()));

      totalShift += mean; //store the sum of shifts in the very first element of the vector
      //cout << "Station " << it->first << " \tReference phase: " << totalShift << " \tRMS: " << rms << " \tMean: " << mean << endl;
    } while ( abs(mean) >=  fConvergenceLimit);

    // check if final result is inbetween +/- pi
    while ( totalShift > kPi ) {
      totalShift -= kTwoPi;
    }
    while ( totalShift < -kPi ) {
      totalShift += kTwoPi;
    }

    // store calculated reference phase and its RMS
    referencePhases[it->first] = pair<double,double>(totalShift, rms);
  }

  return referencePhases;
}


void
CalcBeaconRefPhase::writeFileForDatabase(void)
{
  string startDate = "2012-08-15";
  
  // Create file name containing frequency
  stringstream outFilenameStringStream;
  outFilenameStringStream << fOutputRefPhaseFile << "_ScriptForDatabase";
  string outFilename(outFilenameStringStream.str());

  stringstream message;
  message.str("");
  message << "\nWriting script for database to file: " << outFilename;
  INFO(message.str());

  ofstream outFile;
  outFile.open(outFilename.c_str());

  // loop through all station IDs and beacon frequencies
  for (set<unsigned int>::iterator it = fExistingStationIDs.begin(); it != fExistingStationIDs.end(); ++it) {
    for (unsigned int nBeaconFreq = 0; nBeaconFreq < fBeaconFrequencies.size(); ++nBeaconFreq) {
      outFile << "update BeaconRefPhase set Stop='" << startDate << "' where BeaconFreq_id=" << nBeaconFreq+1 
              << " and RStationId=" << *it << " and Stop>'" << startDate << "';" << endl;
      outFile << "INSERT INTO BeaconRefPhase (BeaconFreq_id, RStationId,  RefPhase, Start) VALUES ('"
              << nBeaconFreq+1 << "', '" << *it << "', '" <<  fReferencePhases[nBeaconFreq][*it].first
              << "', '" << startDate << "' );" << endl;
    }
    outFile << endl;
  }

  outFile.close();
}


void
CalcBeaconRefPhase::writeReferencePhasesFile(void)
{
  // Create file name containing frequency
  stringstream outFilenameStringStream;
  outFilenameStringStream << fOutputRefPhaseFile;
  string outFilename(outFilenameStringStream.str());

  stringstream message;
  message.str("");
  message << "\nWriting reference phases to file: " << outFilename;
  INFO(message.str());

  ofstream outFile;
  outFile.open(outFilename.c_str());

  // write file header
  outFile << "Station";
  for (unsigned int beaconFreqN = 0; beaconFreqN < fBeaconFrequencies.size(); ++beaconFreqN)
      outFile << " RefPhase[rad]_for_" << fBeaconFrequencies[beaconFreqN]/megahertz << "MHz";
  outFile << endl;

  // loop through all station IDs and beacon frequencies
  for (set<unsigned int>::iterator it = fExistingStationIDs.begin(); it != fExistingStationIDs.end(); ++it) {
    outFile << *it;
    for (unsigned int nBeaconFreq = 0; nBeaconFreq < fBeaconFrequencies.size(); ++nBeaconFreq) {
      outFile << " " << fReferencePhases[nBeaconFreq][*it].first;
    }
    outFile << endl;
  }

  outFile.close();
}

}