RdTimeCalibration.cc

Go to the documentation of this file.

#include "RdTimeCalibration.h"
#include <fwk/VModule.h>
#include <fwk/CoordinateSystemRegistry.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 <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Station.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/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/AugerException.h>
#include <utl/CoordinateSystem.h>
#include <utl/Point.h>
#include <utl/Vector.h>
//#include <utl/UTCDateTime.h>

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

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

const double speedOfLight = kSpeedOfLight / 1.00025; // for adjusting refractive index

namespace RdTimeCalibration {

  RdTimeCalibration::RdTimeCalibration() :
    fInfoLevel(eDebug),
    //fRunTimeExpectation(map<int, vector<double> >()),
    //fSourcePosition(vector<double>()),
    fReferenceAntennaID(145),
    fReferenceChannel(1),
    fAdjustTime(false),
    fFirstRun(true),
    fZone(19),
    fBand('H'),
    fEllipsoidName("WGS84")
  { }


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

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

    topBranch.GetChild("SourcePosition").GetData(fSourcePosition);
    topBranch.GetChild("ReferenceChannel").GetData(fReferenceChannel);
    topBranch.GetChild("ReferenceAntennaID").GetData(fReferenceAntennaID);
    topBranch.GetChild("AdjustTime").GetData(fAdjustTime);

    //check if the SourcePosition has the right format
    if (fSourcePosition.size() < 3) {
      WARNING("Incorrect choice of position: One or more coordinates are missing!");
      return eBreakLoop;
    }
    cout << fSourcePosition[1] << endl;

    // get reference coordinate system of detector (usually PampaAmarilla)
    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    const ReferenceEllipsoid::EllipsoidID ellipsoid = ReferenceEllipsoid::GetEllipsoidIDFromString(fEllipsoidName);
    const Point posSource(fSourcePosition[0], fSourcePosition[1], fSourcePosition[2], referenceCS);
    const UTMPoint utmPosSource(posSource, ellipsoid);
    ostringstream message;
    message << "Position of source in UTM: "
               "Easting " << utmPosSource.GetEasting() << ", "
               "Northing " << utmPosSource.GetNorthing() << ", "
               "Height" << utmPosSource.GetHeight();
    INFO(message);

    return eSuccess;
  }


  VModule::ResultFlag
  RdTimeCalibration::Run(evt::Event& event)
  {
    ofstream ofile;
    INFO("RdTimeCalibration::Run()");
    ostringstream message;

    // Check if there are events at all
    if (!event.HasREvent()) {
      WARNING("RdTimeCalibration::No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();

    //only take events if RefAnt is in
    if (!rEvent.HasStation(fReferenceAntennaID)) {
      WARNING("RdTimeCalibration::Reference Antenne is not in this event!");
      return eContinueLoop;
    }
    //get the detector description
    det::Detector& det = det::Detector::GetInstance();
    const rdet::RDetector& radioDet = det.GetRDetector();

    // get reference coordinate system of detector (usually PampaAmarilla)
    const CoordinateSystemPtr referenceCS = det.GetReferenceCoordinateSystem();

    if (fFirstRun) {
      if (fInfoLevel == eDebug) {
        INFO("Entering branch for 'first run'.");
      } 
      //get data of the chosen Reference Antenna
      const rdet::Station& referenceAntenna = radioDet.GetStation(fReferenceAntennaID);
      const utl::Point position = referenceAntenna.GetPosition();
      const Point sourcePosition(fSourcePosition[0], fSourcePosition[1], fSourcePosition[2], referenceCS);
      
      const double referenceDelay = (position - sourcePosition).GetMag() / speedOfLight;

      //calculate the distance, delay and relative delay
      for (rdet::RDetector::StationIterator rdIt = radioDet.StationsBegin();
           rdIt != radioDet.StationsEnd(); ++rdIt) {

        //position of stations
        const Point position = rdIt->GetPosition();

        vector<double> distanceDelayRelativeDelay;
        distanceDelayRelativeDelay.reserve(3);
        distanceDelayRelativeDelay.push_back((position - sourcePosition).GetMag());
        distanceDelayRelativeDelay.push_back(distanceDelayRelativeDelay[0] / speedOfLight);
        distanceDelayRelativeDelay.push_back(distanceDelayRelativeDelay[1] - referenceDelay);

        fRunTimeExpectation[rdIt->GetId()] = distanceDelayRelativeDelay;
      }
      fFirstRun = false;
    }
    // match time stamps of stations (if different)
    // this functionalty should be moved to its own module later
    MatchStationTimeStamps(rEvent);

    //cross correlation
    map<int, double> constTimeOffset;
    const Channel& channelRef = rEvent.GetStation(fReferenceAntennaID).GetChannel(fReferenceChannel);
    const ChannelTimeSeries& timeSeriesRef = channelRef.GetChannelTimeSeries();
    const ChannelTimeSeries::SizeType referenceSampleLength = channelRef.GetChannelTimeSeries().GetSize();
    const double binSize = channelRef.GetChannelTimeSeries().GetBinning();

    //Begin and End of the Sum for the crosscorrelation
    const unsigned int sumIndex[2] = { (unsigned int)((0.5 - 0.15)*referenceSampleLength), (unsigned int)((0.5 + 0.15)*referenceSampleLength) };

    if (fInfoLevel == eDebug) {
      INFO("Iterating over stations");
    } 
    for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
         sIt != rEvent.CandidateStationsEnd(); ++sIt) {
      Station& station = *sIt;

      //building the name for debug output files (file for every Station, name is Station ID)
      /*string s = boost::lexical_cast<string>(sIt->GetId());
      ofile.open(s.c_str(), ios::app);
      */

      //cout << "Station " << station.GetId() << " channel" << flush;

      for (Station::ConstChannelIterator cIt = station.ChannelsBegin();
           cIt != station.ChannelsEnd(); ++cIt) {
        const Channel& channel = *cIt;
        if (channel.GetId() != fReferenceChannel)
            continue;

        //cout << ' ' << channel.GetId() << endl;    

        const ChannelTimeSeries& timeSeries = channel.GetChannelTimeSeries();
        
        /*ofile <<  "Event-Zeit: " << (UTCDateTime(event.GetHeader().GetTime())).GetInAugerFormat()
              <<" For Station " << station.GetId()
              << "(Expected: " << fRunTimeExpectation[station.GetId()][2] << ")";
        */

        const int offsetStart = round(fRunTimeExpectation[station.GetId()][2] / binSize) - round(1500. / binSize);
        const int offsetStop = round(fRunTimeExpectation[station.GetId()][2] / binSize) + round(1500. / binSize);
        double offsetOfMax = 0;
        double max = 0;
        
        //cout << "offsetStart = " << offsetStart << " \t offsetStop = " << offsetStop << endl;
        //cout << "Size of timeSeries: " << timeSeries.GetSize() << "\t lowest sample needed = " << offsetStart + sumIndex[0] << "\t highest sample needed = " << offsetStop + sumIndex[1] << endl;
        
        if ((offsetStop + sumIndex[1]) >= timeSeries.GetSize()) {
          station.SetRejected(eBeaconCorrectionFailed);
          message.str("");
          message << "Trace of station " << station.GetId() << " has only "
                  << timeSeries.GetSize() << " samples, but a mimimum of "
                  << offsetStop + sumIndex[1] + 1 << " samples required - station rejected!";
          WARNING(message);
          continue;
        }

        if ((offsetStart < 0) && (std::abs(offsetStart) > sumIndex[0])) {
          ERROR("Trace to short!");
          return eBreakLoop;
        } 
          
        for (int offset = offsetStart; offset <= offsetStop; ++offset) {
          double cc = 0;
          //double refNorm = 0;
          //double norm = 0;
          for (unsigned int j = sumIndex[0]; j <= sumIndex[1]; ++j) {
            cc += timeSeriesRef[j] * timeSeries[std::abs(int(j) + offset)];
            //only necessary if you want calculate the normalized CC for debugging
            /*refNorm +=  pow((channelRef.GetChannelTimeSeries()[j] / (micro*volt/meter)), 2);
            norm += pow((channel.GetChannelTimeSeries()[j + offset] / (micro*volt/meter)), 2);
            */
          }

          if (cc > max &&
              offset > offsetStart + round(1100.0 / binSize) &&
              offset < offsetStop - round(1100.0 / binSize)) {
            max = cc;
            //normalized CC for debugging
            //NCC = CC/sqrt(RefNorm * Norm);
            offsetOfMax = offset * binSize;
            //cout << "Set new maximum for offset * binSize = " << offset << " * " << binSize << endl;
          }
        }
        //uncomment this for debugging output
        /*ofile << OffsetOfMax << " " << NCC << " "
              << OffsetOfMax - fRunTimeExpectation[sIt->GetId()][2] << endl << endl;
        */
        constTimeOffset[station.GetId()] = offsetOfMax - fRunTimeExpectation[sIt->GetId()][2];
        //cout << "Calculation for station " << station.GetId() << "   (offsetOfMax - fRunTimeExpectation[sIt->GetId()][2] = constTimeOffset): " 
        //     << offsetOfMax << " - " << fRunTimeExpectation[sIt->GetId()][2] << " = " << constTimeOffset[station.GetId()] << endl;
        //ofile.close();
      }
    }

    //shift the start time of station by the constant Offset
    if (fAdjustTime) {
      //move to a own function?
      for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
           sIt != rEvent.CandidateStationsEnd(); ++sIt) {
        Station& station = *sIt;
        station.GetRecData().SetParameter(
          eTraceStartTime,
          station.GetRecData().GetParameter(eTraceStartTime) - constTimeOffset[station.GetId()]
        );

        message.str("");
        message << "Shifted trace start time of Station " << station.GetId() << " by "
                << constTimeOffset[station.GetId()] << " ns";
        INFO(message);
      }
    }

    return eSuccess;
  }


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

    INFO("RdTimeCalibration::Finish()");

    return eSuccess;
  }


  void
  RdTimeCalibration::MatchStationTimeStamps(revt::REvent& rEvent)
    const
  {
    if (fInfoLevel == eDebug) {
      INFO("Matching station time stamps");
    }  
    // store time stamp of the first station as reference
    const 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;
      const double stationTimeStamp = station.GetRecData().GetParameter(eTraceStartTime);

      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);
            //cout << "RdTimeCalibration: Station " << channel.GetStationId() << " channel " << channel.GetId() << " shifted by " << stationTimeStamp-firstTimeStamp << "ns." << endl;
        }
        station.GetRecData().SetParameter(eTraceStartTime, firstTimeStamp);
      }
    }
  }

}