RdStationInterpolatorStarShape.cc

Go to the documentation of this file.

#include "RdStationInterpolatorStarShape.h"
#include <fwk/CentralConfig.h>
#include <fwk/VModule.h>

#include <rdet/RDetector.h>
#include <rdet/Station.h>
#include <utl/CoordinateSystem.h>
#include <utl/Point.h>
#include <evt/RadioSimulation.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/UTCDateTime.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <utl/RadioGeometryUtilities.h>
#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <revt/REvent.h>
#include <revt/EventTrigger.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationTriggerData.h>
#include <revt/StationSimData.h>
#include <utl/StringCompare.h>
#include <cmath>
#include <algorithm>
#include <set>

#include <utl/FFTDataContainer.h>
#include <string>
#include <evt/SimRadioPulse.h>
#include <utl/MathConstants.h>
#include <complex>
#include <map>
#include <vector>
#include <iomanip>
#include <fstream>
#include <math.h>
#include <boost/utility.hpp>
#include <TRandom2.h>
#include <boost/algorithm/string.hpp>
#include <stdlib.h>
#include <rdet/Channel.h>

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "
using namespace std;
using namespace evt;
using namespace det;
using namespace utl;
using namespace revt;
using namespace fwk;

namespace RdStationInterpolatorStarShape {

  RdStationInterpolatorStarShape::RdStationInterpolatorStarShape() :
  fMaximumDistance(5.*meter), // not yet used since not yet tested, up to which distance the interpolation works properly
  fTracePaddingAtFront(2000.*nanosecond),
  fMinimumTraceLength(12000.*nanosecond),
  fEventRate(0.0001*hertz),
  fEventTimeOffset(0.0),
  fUseRateTiming(false),
  fIncludedStationIds(vector<int>()),
  fExcludedStationIds(vector<int>()),
  fInfoLevel(-1)
  {
  }

  RdStationInterpolatorStarShape::~RdStationInterpolatorStarShape()
  {
  }

  fwk::VModule::ResultFlag RdStationInterpolatorStarShape::Init()
  {
    // Read in the configurations of the xml file
    OUT(eFinal)<<"RdStationInterpolatorStarShape::Init()"<<endl;
    string tmpstring;
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationInterpolatorStarShape");
    topBranch.GetChild("MaximumAllowedDistance").GetData(fMaximumDistance);
    topBranch.GetChild("TracePaddingAtFront").GetData(fTracePaddingAtFront);
    topBranch.GetChild("MinimumTraceLength").GetData(fMinimumTraceLength);
    topBranch.GetChild("SimEventRate").GetData(fEventRate);
    topBranch.GetChild("FirstEventTime").GetData(fFirstEventTime);
    topBranch.GetChild("UseRateTiming").GetData(tmpstring);
    fUseRateTiming = utl::StringEquivalent(tmpstring,"yes");
    topBranch.GetChild("IncludedStationIds").GetData(fIncludedStationIds);
    topBranch.GetChild("ExcludedStationIds").GetData(fExcludedStationIds);
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    ostringstream info;
    if (fUseRateTiming) {
      info << "\n Ignoring event times provided by simulations! Using rate timing instead.\n"
           << " First event date " << UTCDateTime(fFirstEventTime) << "\n"
           << " Events have a rate of " << fEventRate/hertz << " Hz\n";
    }
    info << "\n Padding time series data at front by "
         << fTracePaddingAtFront/nanosecond
         << " ns and at end to a total length of at least "
         << fMinimumTraceLength/nanosecond
         << " ns."
         << endl;
    OUT(eFinal) <<info.str() << endl;

    return eSuccess;
  }//Init

  fwk::VModule::ResultFlag RdStationInterpolatorStarShape::Run(evt::Event& theevent)
  {
    OUT(eIntermediate) <<"RdStationInterpolatorStarShape::Run" << endl;

    evt::ShowerSimData& simshow = theevent.GetSimShower();
    evt::RadioSimulation& rsim = simshow.GetRadioSimulation();

    const utl::CoordinateSystemPtr corsys = rsim.GetLocalCoordinateSystem();
    BuildSimShower(simshow,corsys);
    Detector& Det = Detector::GetInstance();
    TimeStamp eventtime = rsim.GetEventTime();

    // overwrite event time if fUseRateTiming is activated
    if (fUseRateTiming)
      eventtime = fFirstEventTime+fEventTimeOffset;

    // check if the event time has an uninitialized (illegal) value
    if (eventtime == TimeStamp(0,0)) {
      WARNING("RdStationInterpolatorStarShape: Event time of RadioSimulation is not available. Aborting ...");
      return eFailure;
    }

    Det.Update(eventtime);
    const rdet::RDetector& RadioDet = Det.GetRDetector();

    if(!theevent.HasREvent())
        theevent.MakeREvent();

    revt::REvent& therev = theevent.GetREvent();

    if (!rsim.GoToFirstSimRadioPulse()) {
      WARNING("RdStationInterpolatorStarShape: No SimRadioPulse in RadioSimulation. Aborting ...");
      return eFailure;
    }

    therev.GetHeader().SetId(rsim.GetEventNumber());
    therev.GetHeader().SetRunNumber(rsim.GetRunNumber());
    therev.GetHeader().SetTime(eventtime);
    simshow.MakeTimeStamp(eventtime);

    // Add trigger information to the event, by default events are assumed as externally triggered
    therev.MakeTrigger();
    therev.GetTrigger().SetExternalTrigger(true);

    // determine the spacings in radius in the shower plane and in azimuth
    vector<SimPulseEntry> simPulses;
    set<double> radiusValues;
    set<double> azimuthValues;
    set<double> altitudeValues;
    bool ok = true;
    do {
      const SimRadioPulse& srpulse = rsim.GetNextSimRadioPulse(ok);
      // leave the loop if there were no further SimRadioPulses
      if (!ok)
        break;

      // register simPulses
      double radius = RadioGeometryUtilities::GetDistanceToAxis(simshow.GetDirection(), simshow.GetPosition(),
        srpulse.GetLocation());

      double azimuth = (srpulse.GetLocation()-simshow.GetPosition()).GetPhi(corsys);

      double altitude = srpulse.GetLocation().GetZ(corsys);

      simPulses.push_back(SimPulseEntry(radius, azimuth, &srpulse));
      radiusValues.insert(long(radius*1000+0.5)/1000.); // store rounded radius value
      azimuthValues.insert(long(azimuth*1000+0.5)/1000.); // store rounded azimuth value
      altitudeValues.insert(long(altitude*1000+0.5)/1000.); // store rounded altitude value

    } while (ok); //do

    // check that we actually have a star-shape and count rings and rays
    if (radiusValues.size()*azimuthValues.size() != simPulses.size()) {
      WARNING("RdStationInterpolatorStarShape: Simulation does not have star shape. Aborting ...");
      return eFailure;
    }

    // output some info
    std::cout << "\nSimulation has " << radiusValues.size() << " rings, " << azimuthValues.size()
              << " angle bins and spans heights from " << (*altitudeValues.begin()) / meter
              << " m to " << (*altitudeValues.rbegin()) / meter << " m above core.\n\n" << endl;

    //fixme TH: should check for maximum height difference and abort if too big

    // sort the list of simpulses in rings with increasing azimuths
    sort(simPulses.begin(), simPulses.end());

    // now loop over radio stations and interpolate them from the sim-pulses
    for (rdet::RDetector::StationIterator rdIt = RadioDet.StationsBegin(); rdIt != RadioDet.StationsEnd(); ++rdIt)
    {
      utl::Point Position = rdIt->GetPosition();
      int statID = rdIt->GetId();

      OUT(eObscure) << "Station " << statID;

      // if station id is listed in fExcludedStationIds, skip this simulated radio pulse (it will not be re-used for the second-closest antenna station!)
      if (not fExcludedStationIds.empty())
        if ( find(fExcludedStationIds.begin(), fExcludedStationIds.end(), statID) != fExcludedStationIds.end() ) {
          OUT(eObscure) << " -> not associated as it is in the exclude list." << endl;
          continue;
        }

      // if there is a list of fIncludedStationIds skip this simulated radio pulse if station id is not in the list
      if (not fIncludedStationIds.empty())
        if ( find(fIncludedStationIds.begin(), fIncludedStationIds.end(), statID) == fIncludedStationIds.end() ) {
          OUT(eObscure) <<" -> station not associated as it is not in the include list." << endl;
          continue;
        }

      //find the nearest neighbors of simulated pulses to the station (station ID, timeseries)
      double stationRadius = RadioGeometryUtilities::GetDistanceToAxis(simshow.GetDirection(), simshow.GetPosition(),
        Position);

      double stationAzimuth = (Position-simshow.GetPosition()).GetPhi(corsys);

      double stationAltitude = (Position.GetZ(corsys));

      OUT(eObscure) << " associated at radius " << stationRadius/meter
                    << " m with azimuth " << stationAzimuth/deg
                    << " degree at altitude " << stationAltitude / meter
                    << " m\n";

      // special treatment if station is inside innermost ring or outside outermost ring
      if (stationRadius >= (*radiusValues.rbegin())) {
//      cout << "Station is outside outermost ring, skip it.\n";
        continue;
      }
      else if (stationRadius <= (*radiusValues.begin())) {
          cout << "\nStation " << statID << " is in innermost ring, skip it.\n\n";
          continue;
          //fixme TH: devise scheme to interpolate in innermost ring
      }

      // find the four adjacent sim radio pulses
      long indexHighHigh = distance(radiusValues.begin(),radiusValues.upper_bound(stationRadius))*azimuthValues.size()+distance(azimuthValues.begin(),azimuthValues.upper_bound(stationAzimuth));
      long indexHighLow = indexHighHigh-1;
      // check if we have to apply a special treatment because of wraparound for periodic azimuth values
      if (indexHighHigh >= long(simPulses.size()) ||
          fabs(simPulses.at(indexHighHigh).fDistanceFromShowerAxis-simPulses.at(indexHighLow).fDistanceFromShowerAxis)/simPulses.at(indexHighHigh).fDistanceFromShowerAxis > 1.0e-3)
      {
        indexHighHigh -= azimuthValues.size();
      }
      long indexLowHigh = indexHighHigh-azimuthValues.size();
      long indexLowLow = indexHighLow-azimuthValues.size();

/*
      std::cout << "radius " << simPulses.at(indexLowLow).fDistanceFromShowerAxis << " azimuth " << simPulses.at(indexLowLow).fObserverAzimuthAngle << " number " << indexLowLow << "\n";
      std::cout << "radius " << simPulses.at(indexLowHigh).fDistanceFromShowerAxis << " azimuth " << simPulses.at(indexLowHigh).fObserverAzimuthAngle << " number " << indexLowHigh << "\n";
      std::cout << "radius " << simPulses.at(indexHighLow).fDistanceFromShowerAxis << " azimuth " << simPulses.at(indexHighLow).fObserverAzimuthAngle << " number " << indexHighLow << "\n";
      std::cout << "radius " << simPulses.at(indexHighHigh).fDistanceFromShowerAxis << " azimuth " << simPulses.at(indexHighHigh).fObserverAzimuthAngle << " number " << indexHighHigh << "\n";
*/

      //get the design frequencies from the channels of the station and check if they are the same for all channels
      double designLowerFreq = 0.0;
      double designLowerFreqFirstCh = 0.;
      double designUpperFreq = 0.0;
      double designUpperFreqFirstCh = 0.;
      for(rdet::Station::ChannelIterator cIt = rdIt->ChannelsBegin(); cIt!= rdIt->ChannelsEnd(); ++cIt){
        const rdet::Channel& cDetChannel = RadioDet.GetStation(rdIt->GetId()).GetChannel(cIt->GetId());
        designLowerFreq = cDetChannel.GetDesignLowerFreq();
        designUpperFreq = cDetChannel.GetDesignUpperFreq();
        if (cIt == rdIt->ChannelsBegin()) {
          designLowerFreqFirstCh = designLowerFreq;
          designUpperFreqFirstCh = designUpperFreq;
        }
        if(designLowerFreq != designLowerFreqFirstCh || designUpperFreq != designUpperFreqFirstCh) {
          WARNING("RdStationInterpolatorStarShape: design frequency different for different channels");
          return eFailure;
        }
      }

/*
      // do not limit interpolation to design frequency band
      designLowerFreqFirstCh=0.;
      designUpperFreqFirstCh=1.e12;
*/

      // now interpolate between the four positions
      const double starttime1 = simPulses.at(indexLowLow).fRadioPulse->GetStartTime();
      const double radius1 = simPulses.at(indexLowLow).fDistanceFromShowerAxis;
      const double azimuth1 = simPulses.at(indexLowLow).fObserverAzimuthAngle;
      const double binning1 = simPulses.at(indexLowLow).fRadioPulse->GetEfieldTimeSeries().GetBinning();
      // padding of zeroes to end of timeseries
      StationFFTDataContainer simFFTDataContainer1;
      simFFTDataContainer1.GetTimeSeries() = PadTimeSeries(*simPulses.at(indexLowLow).fRadioPulse);

      const double starttime2 = simPulses.at(indexHighLow).fRadioPulse->GetStartTime();;
      const double radius2 = simPulses.at(indexHighLow).fDistanceFromShowerAxis;
      const double binning2 = simPulses.at(indexHighLow).fRadioPulse->GetEfieldTimeSeries().GetBinning();
      StationFFTDataContainer simFFTDataContainer2;
      simFFTDataContainer2.GetTimeSeries() = PadTimeSeries(*simPulses.at(indexHighLow).fRadioPulse);

      const double starttime3 = simPulses.at(indexLowHigh).fRadioPulse->GetStartTime();;
      const double azimuth2 = simPulses.at(indexLowHigh).fObserverAzimuthAngle;
      const double binning3 = simPulses.at(indexLowHigh).fRadioPulse->GetEfieldTimeSeries().GetBinning();
      StationFFTDataContainer simFFTDataContainer3;
      simFFTDataContainer3.GetTimeSeries() = PadTimeSeries(*simPulses.at(indexLowHigh).fRadioPulse);

      const double starttime4 = simPulses.at(indexHighHigh).fRadioPulse->GetStartTime();;
      const double binning4 = simPulses.at(indexHighHigh).fRadioPulse->GetEfieldTimeSeries().GetBinning();
      StationFFTDataContainer simFFTDataContainer4;
      simFFTDataContainer4.GetTimeSeries() = PadTimeSeries(*simPulses.at(indexHighHigh).fRadioPulse);

      // check for same binning of timeseries. Must have the same binning
      if ( (fabs(binning2 - binning1)/binning1 > 1.e-4)
        || (fabs(binning3 - binning1)/binning1 > 1.e-4)
        || (fabs(binning4 - binning1)/binning1 > 1.e-4) )
      {
            WARNING("RdStationInterpolatorStarShape: simtimeseries don't have the same binning. Aborting ...");
            return eFailure;
      }

      if ( simFFTDataContainer1.GetNyquistZone() != 1 || simFFTDataContainer2.GetNyquistZone() != 1 ) {
        WARNING("RdStationInterpolatorStarShape: simtimeseries not in first Nyquist zone.");
        return eFailure;
      }

      if ( simFFTDataContainer1.GetConstFrequencySpectrum().GetSize() != simFFTDataContainer2.GetConstFrequencySpectrum().GetSize() ) {
        WARNING("RdStationInterpolatorStarshape: simtimeseries don't have the same length.");
        return eFailure;
      }

      // for each of the azimuth bins interpolate linearly between the two rings
      const revt::StationFFTDataContainer InterFFTDataContainerLeft = Interpolate(simFFTDataContainer1, simFFTDataContainer2, radius1, radius2, stationRadius, designLowerFreqFirstCh, designUpperFreqFirstCh);
      const double starttimeLeft = starttime1*(radius2-stationRadius)/(radius2-radius1)+starttime2*(stationRadius-radius1)/(radius2-radius1);
      const revt::StationFFTDataContainer InterFFTDataContainerRight = Interpolate(simFFTDataContainer3, simFFTDataContainer4, radius1, radius2, stationRadius, designLowerFreqFirstCh, designUpperFreqFirstCh);
      const double starttimeRight = starttime3*(radius2-stationRadius)/(radius2-radius1)+starttime4*(stationRadius-radius1)/(radius2-radius1);

      // now interpolate the resulting spectra in azimuth
      const revt::StationFFTDataContainer InterFFTDataContainer = Interpolate(InterFFTDataContainerLeft, InterFFTDataContainerRight, azimuth1, azimuth2, stationAzimuth, designLowerFreqFirstCh, designUpperFreqFirstCh);
      const double InterStartTime = starttimeLeft*(azimuth2-stationAzimuth)/(azimuth2-azimuth1)+starttimeRight*(stationAzimuth-azimuth1)/(azimuth2-azimuth1);
      const revt::StationTimeSeries InterTimeSeries = InterFFTDataContainer.GetConstTimeSeries();

      //fill SimRadioPulse into the associated station

      therev.MakeStation(statID); // Create a new Station, already checks if the station exists
      revt::Station& Stat = therev.GetStation(statID);

      // create data structures
      Stat.MakeGPSData();     // should check before if exists?
      Stat.MakeTriggerData(); // should check before if exists?
      Stat.MakeSimData();     // should check before if exists?

      revt::StationTimeSeries& StatTimeSeries = Stat.GetStationTimeSeries() ;  // reference to the Station time series

      // check if there is data already and throw a warning
      if(StatTimeSeries.GetSize() != 0)
      {
        WARNING( "RDStationInterpolator: StatTimeSeries not empty!" );
        StatTimeSeries.ClearTrace();
      }

      // fill the simtimeseries into the Station TimeSeries
      StatTimeSeries.SetBinning(binning1);

      long numpaddedsamples = static_cast<long>(fTracePaddingAtFront/binning1+0.5);

      // interpolate and set start time
      Stat.SetRawTraceStartTime(therev.GetHeader().GetTime() + TimeInterval(InterStartTime - numpaddedsamples*binning1));

      // fill with the interpolated trace
      for (revt::StationTimeSeries::ConstIterator iter=InterTimeSeries.Begin(); iter!=InterTimeSeries.End(); ++iter)
      {
        StatTimeSeries.PushBack(*iter);
      }

    }// for (rdet::RDetector::StationIterator...)

    return eSuccess;
  } // Run

  fwk::VModule::ResultFlag RdStationInterpolatorStarShape::Finish()
  {
    OUT(eIntermediate)<<"RdStationInterpolatorStarShape::Finish()"<<endl;
    return eSuccess;
  } // Finish




  revt::StationTimeSeries RdStationInterpolatorStarShape::PadTimeSeries(const SimRadioPulse& simtimeseries)
  {
    revt::StationTimeSeries PaddedTimeSeries;

    Vector3D V3DZero;
    V3DZero = 0.0, 0.0, 0.0;

    double padbinning = simtimeseries.GetBinning();
    PaddedTimeSeries.SetBinning(padbinning);

    long numpaddedsamples = static_cast<long>(fTracePaddingAtFront/padbinning+0.5);

    for (long i=0; i<numpaddedsamples; ++i)
      PaddedTimeSeries.PushBack(V3DZero);

    for (revt::StationTimeSeries::ConstIterator iter=(simtimeseries.GetEfieldTimeSeries()).Begin(); iter!=(simtimeseries.GetEfieldTimeSeries()).End(); ++iter)
      PaddedTimeSeries.PushBack(*iter);

    // now pad the station trace at the end up to the given minimum trace length
    while (fMinimumTraceLength - PaddedTimeSeries.GetSize()*padbinning > 1.e-6)
        PaddedTimeSeries.PushBack(V3DZero);

    // if number of samples is uneven after padding, pad one more sample (otherwise it will again be clipped at the next FFT)
    if (PaddedTimeSeries.GetSize()/2*2 != PaddedTimeSeries.GetSize())
        PaddedTimeSeries.PushBack(V3DZero);

    return PaddedTimeSeries;
  }



  revt::StationFFTDataContainer RdStationInterpolatorStarShape::Interpolate( const revt::StationFFTDataContainer& simData1, const revt::StationFFTDataContainer& simData2, double x1, double x2, double x, double designlowerfreq, double designupperfreq)
  {
    const double specbinning1 = (simData1.GetConstFrequencySpectrum()).GetBinning();

    const StationFrequencySpectrum::SizeType size1 = (simData1.GetConstFrequencySpectrum()).GetSize();
    const StationFrequencySpectrum& freq1 = (simData1.GetConstFrequencySpectrum());

    // unused    const StationFrequencySpectrum::SizeType size2 = (simData2.GetConstFrequencySpectrum()).GetSize();
    const StationFrequencySpectrum& freq2 = (simData2.GetConstFrequencySpectrum());

    TraceV3C V_interspectrum;
    V_interspectrum.SetBinning(specbinning1);
    Vector3C element;

    Vector3D phaseElement1;
    Vector3D phaseElement2;
    Vector3D phaseElementInter;

    Vector3D jumps1;
    jumps1 = 0.0,0.0,0.0;
    Vector3D jumps2;
    jumps2 = 0.0,0.0,0.0;

    for(StationFrequencySpectrum::SizeType i=0 ; i < size1 ; i++)
    {
      for(unsigned int j = 0; j < 3; j++) // loop over x,y and z component
      {
        if( i >= designlowerfreq/specbinning1 && i <= designupperfreq/specbinning1 ) {
          // linearly interpolate amplitudes
          double interamplitude = (abs(freq1[i][j])*(x2-x)/(x2-x1) + abs(freq2[i][j])*(x-x1)/(x2-x1));

          // unwrap phases of initial signals
          if(i > 1) {

            if( arg(freq1[i][j]) < arg(freq1[i-1][j]) && jumps1[j] >= 0 )
            {
              if( fabs(arg(freq1[i][j])-arg(freq1[i-1][j])) > fabs(arg(freq1[i][j]) + utl::kPi - arg(freq1[i-1][j])) )
              {
                ++jumps1[j];
              }
            }
            else {
              if( fabs(arg(freq1[i][j])-arg(freq1[i-1][j])) > fabs(arg(freq1[i][j]) - utl::kPi - arg(freq1[i-1][j])) )
              {
                --jumps1[j];
              }
            }

            if( arg(freq2[i][j]) < arg(freq2[i-1][j]) && jumps2[j] >= 0 )
            {
              if( fabs(arg(freq2[i][j])-arg(freq2[i-1][j])) > fabs(arg(freq2[i][j]) + utl::kPi - arg(freq2[i-1][j])) )
              {
                ++jumps2[j];
              }
            }
            else {
              if( fabs(arg(freq2[i][j])-arg(freq2[i-1][j])) > fabs(arg(freq2[i][j]) - utl::kPi - arg(freq2[i-1][j])) )
              {
                --jumps2[j];
              }
            }

          }

          phaseElement1[j] = arg(freq1[i][j]) + jumps1[j]*2*utl::kPi;
          phaseElement2[j] = arg(freq2[i][j]) + jumps2[j]*2*utl::kPi;

          // interpolate phases
          phaseElementInter[j] = (phaseElement1[j]*(x2-x)/(x2-x1) + phaseElement2[j]*(x-x1)/(x2-x1));

          // interpolated spectrum
          element[j]=polar(interamplitude,phaseElementInter[j]);
        }
        else{
          element[j] = 0.0;
        }
      }

      V_interspectrum.PushBack(element);
    }

    revt::StationFFTDataContainer interpolatedFFTDataContainer;
    interpolatedFFTDataContainer.GetFrequencySpectrum() = V_interspectrum;

    return interpolatedFFTDataContainer;
  }  //Interpolate



  void RdStationInterpolatorStarShape::BuildSimShower(evt::ShowerSimData& theshower,const utl::CoordinateSystemPtr& coord) {
    // This function will build a Sim shower with correct data
    // Inspired by the EventGenerator Module

    if (!theshower.HasPosition())
      theshower.MakeGeometry(Point(0, 0, 0, coord));

  }

} // nameSpace