RdStationSignalInterpolator.cc

Go to the documentation of this file.

#include "RdStationSignalInterpolator.h"

#include <fwk/CentralConfig.h>
#include <utl/ErrorLogger.h>

#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/FFTDataContainer.h>
#include <utl/Point.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/StringCompare.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>

#include <evt/Event.h>
#include <evt/RadioSimulation.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerRRecDataQuantities.h>

#include <revt/REvent.h>
#include <revt/EventTrigger.h>
#include <revt/Station.h>
#include <revt/Header.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/EventTrigger.h>

#include <rdet/RDetector.h>
#include <rdet/Station.h>

#include <atm/AerosolDB.h>
#include <atm/AerosolZone.h>
#include <atm/ProfileResult.h>
#include <atm/MonthlyAvgDBProfileModel.h>
#include <atm/InclinedAtmosphericProfile.h>

using namespace std;
using namespace evt;
using namespace det;
using namespace utl;
using namespace fwk;
using namespace revt;


namespace RdStationSignalInterpolator {

  VModule::ResultFlag
  RdStationSignalInterpolator::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationSignalInterpolator");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    topBranch.GetChild("MinimumStationEnergy").GetData(fMinimumEnergy);

    fUseEarlyLate =
      utl::StringEquivalent(topBranch.GetChild("UseEarlyLateCorrection").Get<string>(), "yes");

    // FS:
    if (fUseEarlyLate)
      ERROR("Early Late Correction will not be used within in the interpolation of "
            "the signals (even if configured). It is not properly implemented.");

    fRequireSDStationTrigger = utl::StringEquivalent(topBranch.GetChild("RequireSDStationTrigger").Get<string>(), "yes");

    topBranch.GetChild("InterpolationMethod").GetData(fInterpolationMethod);

    ostringstream info;
    info << "\n\t Interpolate only stations with energy fluence above " << fMinimumEnergy << " eV."
            "\n\t Interpolation method is " << fInterpolationMethod << "; (0: nearest neighbour, 1: bilinear interpolation,"
            " 2: bicubic interpolation)"
            "\n\t Require SD T2 trigger: " << fRequireSDStationTrigger << "\n";
    INFO(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSignalInterpolator::Run(evt::Event& event)
  {
    INFOFinal("");

    evt::ShowerSimData& simshow = event.GetSimShower();
    Detector& det = Detector::GetInstance();
    const rdet::RDetector& radioDet = det.GetRDetector();

    // At the moment: created Station IDs at StarShape station positions (not the actual Antenna Positions!) start at 10 000
    const int maxStationID = 9999;

    REvent& rEvent = event.GetREvent();
    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    const Vector simShowerAxis = -simshow.GetDirection();
    const Vector magneticField = showerrrec.GetMagneticFieldVector();
    const utl::CoordinateSystemPtr corsys = simshow.GetLocalCoordinateSystem();
    RadioGeometryUtilities transformation = RadioGeometryUtilities(simShowerAxis, corsys, magneticField);

    // Init some stuff needed:
    vector<StarShapeEntry> starShapeStations;
    set<double> radiusValues;
    set<double> azimuthValues;
    set<double> azimuthValuesStarShape;
    set<double> altitudeValues;

    // TODO: Better selection criterium for StarShape pattern identification possible.
    // Careful: it is possible for all arms not to have the same amount of stations!

    ostringstream info;
    // Loop over all star-shaped grid stations -> fill Vector
    int sstationcount = 0;
    for (auto rdIt = rEvent.CandidateStationsBegin(); rdIt != rEvent.CandidateStationsEnd(); ++rdIt) {
      const int statID = rdIt->GetId();
      if (statID < maxStationID)
        continue;

      Station& station = *rdIt;
      const Point position = radioDet.GetStation(station).GetPosition();

      // due to the choice of coordinate system the station position will be evaluated relative to the core position
      double x_vxB = 0;
      double y_vxB = 0;
      double z_vxB = 0;
      transformation.GetVectorInShowerPlaneVxB(x_vxB, y_vxB, z_vxB, position);

      if (!station.HasRecData()) {
        WARNING("Station has no RecData! Please call RdEventInitializer first!");
        return eFailure;
      }

      if (!station.GetRecData().HasParameter(eSignalEnergyFluenceMag)) {
        info.str("");
        info << "Station doesn't have rec. signal..skip. ID: " << statID;
        INFOIntermediate(info);
        continue;
      }

      const double radius = sqrt(pow(x_vxB, 2) + pow(y_vxB, 2));
      const double azimuth = atan2(y_vxB, x_vxB);
      const double altitude = position.GetZ(corsys);
      const unsigned int sizeAzi = azimuthValues.size();

      starShapeStations.emplace_back(radius, azimuth, statID);
      radiusValues.insert(round(radius * 10) / 10); // store rounded radius value
      azimuthValues.insert(round(azimuth * 10000) / 10000);
      altitudeValues.insert(round(altitude)); // store rounded altitude value

      if (azimuthValues.size() == sizeAzi)
        azimuthValuesStarShape.insert(round(azimuth * 1000) / 1000);

      ++sstationcount;
    }

    if (starShapeStations.empty()) {
      ERROR("No antennas on a star. Skip...");
      return eContinueLoop;
    }

    info.str("");
    info << "Simulation has " << radiusValues.size() << " rings, " << azimuthValuesStarShape.size()
         << " angle bins [" << azimuthValues.size() << " angles in total],and spans heights from " << (*altitudeValues.begin()) / meter
         << " m to " << (*altitudeValues.rbegin()) / meter << " m above core.";
    INFOIntermediate(info);

    if (azimuthValuesStarShape.size() != 8) {
      if (azimuthValues.size() == 8) {
        azimuthValuesStarShape = azimuthValues;
      } else {
        WARNING("RdStationInterpolatorStarShape: Simulation does not have 8 arms. Aborting ...");
        return eContinueLoop;
      }
    }
    info.str("");
    info << "Number of star shaped stations: " << sstationcount;
    INFOIntermediate(info);

    //if (fUseEarlyLate) {
    //  //     The following is needed for the early late correction:
    //  //    Early late correction is applied when transforming to the shower-plane coordinate system and it is inversely
    //  //    applied after the interpolation process.
    //  //  Needed for early late correction:
    //  const auto& gh = simshow.GetGHParameters();
    //  const double shower_xmax = gh.GetXMax();

    //  Point simCore = simshow.GetPosition();  // MC core
    //  const atm::Atmosphere& theAtm = Det.GetAtmosphere();

    //  theAtm.InitSlantProfileModel(simCore, simShowerAxis, 5.*g/cm2);
    //  const atm::ProfileResult& distanceVsdepth = theAtm.EvaluateDistanceVsSlantDepth();
    //  // shower maximum
    //  const double distance_to_xmax = abs(distanceVsdepth.Y(shower_xmax));
    //  const Point showerMax = simCore + simShowerAxis * distance_to_xmax;
    //  info.str("");
    //  info << "Applying the early late correction (on energy fluence and distance to axis)"
    //       << "\nDepth of Xmax: " << shower_xmax / (g / cm2)
    //       << "\nDistance of Xmax: " << distance_to_xmax / meter
    //       << "\nDistance of Xmax (sim event): " << simshow.GetDistanceOfShowerMaximum() / meter;
    //  INFOIntermediate(info);
    // }

    // The following enums are interpolated. After that, SNR is calculated based on the formula in the StationSignalReconstructor
    std::vector<StationRRecDataQuantities> eNumNames = {
      eSignalEnergyFluenceMag,
      eSignalEnergyFluenceVxB,
      eSignalEnergyFluenceVxVxB,
      eNoiseEnergyFluenceMag,
      eNoiseEnergyFluenceVxB,
      eNoiseEnergyFluenceVxVxB,
      eSignal,
      eNoise,
      eSignalTime,
      eSignalToNoise
    };

    std::vector<StationRRecDataQuantities> eNumNamesError = {
      eSignalEnergyFluenceMag,
      eSignalEnergyFluenceVxB,
      eSignalEnergyFluenceVxVxB,
      eSignal,
      eNoise,
      eSignalTime,
      eSignalToNoise
    };

    // sort the list of StarShapeStations in arms with same azimuth first, and after with increasing radius
    sort(starShapeStations.begin(), starShapeStations.end());

    // now loop over radio (upgrade) detector stations and interpolate them from the virtual stations.
    unsigned int n_interpolated_stations = 0;
    for (auto rdIt = radioDet.StationsBegin(); rdIt != radioDet.StationsEnd(); ++rdIt) {

      const Point pos = rdIt->GetPosition();
      const int statID = rdIt->GetId();

      if (!rEvent.HasStation(statID))
        rEvent.MakeStation(statID); // Create a new Station, already checks if the station exists

      revt::Station& revtStation = rEvent.GetStation(statID);

      // skip station on star
      if (statID >= maxStationID) {
        revtStation.SetExcludedReason(revt::eManuallyExcluded);
        continue;
      }

      double x_vxB = 0;
      double y_vxB = 0;
      double z_vxB = 0;
      transformation.GetVectorInShowerPlaneVxB(x_vxB, y_vxB, z_vxB, pos);

      //double earlyLateCorrectionFactor = 1;
      //if (fUseEarlyLate){
      //  // checked with distance_to_xmax / (distance_to_xmax + Z_VxB);
      //  const double earlyLateCorrectionFactor =
      //    RadioGeometryUtilities::GetEarlyLateCorrectionFactor(simCore, Pos, showerMax, simShowerAxis);
      //  info.str("");
      //  info << "c early late: " << c_early_late << endl;
      //  INFOIntermediate(info);
      //}

      // Early Late correction (distance)
      //x_vxB /= earlyLateCorrectionFactor;
      //y_vxB /= earlyLateCorrectionFactor;

      const double stationRadius = sqrt(pow(x_vxB, 2) + pow(y_vxB, 2));
      const double stationAzimuth = atan2(y_vxB, x_vxB);

      // skip stations outside the star
      if (stationRadius > *radiusValues.rbegin()) {
        revtStation.SetExcludedReason(revt::eManuallyExcluded);
        continue;
      }

      // skip stations inside the innermost ring.. TODO: Fix Interpolation in innermost ring.
      if (stationRadius <= *radiusValues.begin()) {
        revtStation.SetExcludedReason(revt::eManuallyExcluded);
        continue;
      }

      info.str("");
      info << "StationID: " << statID << ", Station position: r = "
           << stationRadius << ", phi = " << stationAzimuth;
      INFODebug(info);

      // FS: Comment following line because it seems redundant
      //sort(StarShapeStations.begin(), StarShapeStations.end());

      // create data structures
      if (!revtStation.HasGPSData())
        revtStation.MakeGPSData();
      if (!revtStation.HasTriggerData())
        revtStation.MakeTriggerData();
      if (!revtStation.HasSimData())
        revtStation.MakeSimData();
      if (!revtStation.HasRecData())
        revtStation.MakeRecData();

      vector<double> interpolatedValues;
      vector<double> interpolatedValuesError;
      if (fInterpolationMethod == eNN) {
        interpolatedValues =
          NearestNeighbourInterpolation(stationRadius, stationAzimuth, azimuthValuesStarShape,
                                        starShapeStations, eNumNames, rEvent, false);
        interpolatedValuesError =
          NearestNeighbourInterpolation(stationRadius, stationAzimuth, azimuthValuesStarShape,
                                        starShapeStations, eNumNamesError, rEvent, true);
      } else if (fInterpolationMethod == eBilinear) {
        interpolatedValues =
          BilinearInterpolation(stationRadius, stationAzimuth, azimuthValuesStarShape,
                                starShapeStations, eNumNames, rEvent, false);
        interpolatedValuesError =
          BilinearInterpolation(stationRadius, stationAzimuth, azimuthValuesStarShape,
                                starShapeStations, eNumNamesError, rEvent, true);
      } else if (fInterpolationMethod == eBicubic) {
        interpolatedValues =
          BicubicInterpolation(stationRadius, stationAzimuth, azimuthValuesStarShape,
                               starShapeStations, eNumNames, rEvent, false);
        interpolatedValuesError =
          BicubicInterpolation(stationRadius, stationAzimuth, azimuthValuesStarShape,
                               starShapeStations, eNumNamesError, rEvent, true);
      }

      if (fRequireSDStationTrigger) {

        if (!event.HasSEvent()) {
          ERROR("SEvent not available but required!");
          return eFailure;
        }

        const sevt::SEvent& sEvent = event.GetSEvent();

        info.str("");
        info << "SD station " << statID;
        if (sEvent.HasStation(statID)) {
          info << " is in the SEvent, ";
          if (sEvent.GetStation(statID).HasTriggerData()) {
            const sevt::StationTriggerData& trig = sEvent.GetStation(statID).GetTriggerData();
            if (trig.IsT2() || trig.IsT1()){
              info << " and has trigger (T1 or T2). Continue interpolation.";
            } else {
              info << " but doesn't have T1 or T2 trigger data. skip.";
              revtStation.SetExcludedReason(revt::eStationWithoutTrigger);
              INFOIntermediate(info);
              continue;
            }
          } else {
            info << " but doesn't have any trigger data. skip.";
            revtStation.SetExcludedReason(revt::eStationWithoutTrigger);
            INFOIntermediate(info);
            continue;
          }
        } else {
          info << " is not in the SEvent. skip.";
          revtStation.SetExcludedReason(revt::eStationWithoutTrigger);
          INFOIntermediate(info);
          continue;
        }
        INFOIntermediate(info);

      }

      info.str("");
      if (interpolatedValues[0] > fMinimumEnergy && interpolatedValues[0] < std::numeric_limits<double>::max()) {
        info << "Station " << statID << ", station radius: " << stationRadius
             << ", station azimuth: " << stationAzimuth << "--> EnergyFluence: "
             << interpolatedValues[0] << " +/- " << interpolatedValuesError[0];
        INFOIntermediate(info);

        // check if there is data already and throw a warning... This should not happen.
        if (revtStation.GetStationTimeSeries().GetSize())
          WARNING("RDStationSignalInterpolator: StationTimeSeries not empty!");

        info.str("");
        int i = 0;
        for (vector<double>::const_iterator it = interpolatedValues.begin(); it != interpolatedValues.end(); ++it, ++i) {
          info << "Interpolated Value: " << *it << ", enum: " << eNumNames[i] << '\n';
          const double intVal = *it;
          //if (fUseEarlyLate)
          //  intVal /= pow(EnergyFluenceCorrectionEarlyLateEffect, 2);
          revtStation.GetRecData().SetParameter(eNumNames[i], intVal);
        }

        i = 0;
        for (vector<double>::const_iterator it = interpolatedValuesError.begin(); it != interpolatedValuesError.end(); ++it, ++i) {
          info << "Interpolated Value Error: " << *it << ", enum: " << eNumNamesError[i] << '\n';
          revtStation.GetRecData().SetParameterError(eNumNamesError[i], *it);
        }
        revtStation.SetSignal();
        revtStation.GetRecData().SetPulseFound(true);

        info << "Station will be used for Reconstruction. ID is " << statID;
        INFODebug(info);
      } else {
        //revtStation.SetSignal();
        revtStation.SetRejectedReason(eBadSNR);
        info << "Interpolated Energy Fluence of Station " << statID << " is too low(" << interpolatedValues[0] << ") -> skip\n"
             << "RejectionStatus, badSNR: " << revtStation.GetRejectedReason();
        INFODebug(info);
        continue;
      }

      // calculating SNR and Errors (10% on the value, copied from StationSignalReconstructor
      const double peakAmplitude = revtStation.GetRecData().GetParameter(eSignal);
      const double rmsNoise = revtStation.GetRecData().GetParameter(eNoise);
      const double signalToNoise = peakAmplitude / rmsNoise;
      info.str("");
      //info << "interpolated SNR: " << revtStation.GetRecData().GetParameter(eSignalToNoise) << ", ";
      info << "Calculated Signal to Noise ratio: " <<  Sqr(signalToNoise);
      //info << "interpolatedrrors: " << revtStation.GetRecData().GetParameterError(eSignalToNoise)
      //<< ", now calculate it from Signal and Noise: "
      //<<  0.1*Sqr(SignalToNoise) << ", " << endl;
      revtStation.GetRecData().SetParameter(eSignalToNoise, Sqr(signalToNoise));
      revtStation.GetRecData().SetParameterError(eSignalToNoise, 0.1 * Sqr(signalToNoise));
      // hard coded 10%, defined like this in the StationSignalReconstructor
      INFOIntermediate(info);

      ++n_interpolated_stations;
    }

    info.str("");
    info << "Interpolated (and accepted) " << n_interpolated_stations << " signals to a corresponding station on a regular grid.";
    INFOFinal(info);

    return eSuccess;
  }


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

    INFO("RdStationSignalInterpolator::Finish()");

    return eSuccess;
  }


  double
  RdStationSignalInterpolator::PolarDistance(const double r1, const double phi1, const double r2, const double phi2)
  {
    return sqrt(pow(r1, 2) + pow(r2, 2) - 2*r1*r2*cos(phi1 - phi2));
  }


  void
  RdStationSignalInterpolator::MoveItemToBack(std::vector<StarShapeEntry>& v, const size_t itemIndex)
  {
    const auto it = v.begin() + itemIndex;
    std::rotate(it, it + 1, v.end());
  }


  vector<StarShapeEntry>
  RdStationSignalInterpolator::SortForInterpolation(const std::vector<StarShapeEntry>& list)
  {
    auto starShapeList = list;  // copy
    /*ostringstream debug;
    debug.str("");
    debug << starShapeList[0].fObserverAzimuthAngle
          << ", " << starShapeList[4].fObserverAzimuthAngle
          << ", " << starShapeList[8].fObserverAzimuthAngle
          << ", " << starShapeList[12].fObserverAzimuthAngle;*/
    const double az1 = starShapeList[0].fObserverAzimuthAngle;
    const double az2 = starShapeList[4].fObserverAzimuthAngle;
    const double az3 = starShapeList[8].fObserverAzimuthAngle;
    const double az4 = starShapeList[12].fObserverAzimuthAngle;
    //vector<StarShapeEntry> manipulatedStarShapeList = starShapeList;
    const double d1 = az2 - az1;
    const double d2 = az3 - az2;
    const double d3 = az4 - az3;
    //const double d4 = az1 - az4;
    //debug << "\ndistance: \n" << d1 << endl << d2 << endl << d3 << endl << d4 << endl;

    // find the jump from +Pi to -Pi:
    if (fabs(d1) > M_PI*2/7) {
      //debug << "jump between 1 and 2: " << az1 << ", " << az2;
      for (unsigned int i = 0; i < 4; ++i) {
        MoveItemToBack(starShapeList, 0);
      }
      //debug << "Azi, Radius: " << endl;
      /*for (vector<StarShapeEntry>::const_iterator it = starShapeList.begin(); it != starShapeList.end(); ++it) {
        debug << it->fObserverAzimuthAngle << ", " << it ->fDistanceFromShowerAxis << endl;
      }*/
    }
    if (fabs(d2) > M_PI*2/7) {
      //debug << "jump between 2 and 3: " << az2 << ", " << az3;
      for (unsigned int i = 0; i < 8; ++i) {
        MoveItemToBack(starShapeList, 0);
      }
      //debug << "Azi, Radius: " << endl;
      /*for (vector<StarShapeEntry>::const_iterator it = starShapeList.begin(); it != starShapeList.end(); ++it) {
        debug << it->fObserverAzimuthAngle << ", " << it ->fDistanceFromShowerAxis << endl;
      }*/
    }
    if (fabs(d3) > M_PI*2/7) {
      //debug << "jump between 3 and 4: " << az3 << ", " << az4;
      for (unsigned int i = 0; i < 12; ++i) {
        MoveItemToBack(starShapeList, 0);
      }
      /*debug << "Azi, Radius: " << endl;
      for(vector<StarShapeEntry>::const_iterator it = StarShapeList.begin(); it != StarShapeList.end(); it++) {
          debug << it->fObserverAzimuthAngle << ", " << it ->fDistanceFromShowerAxis << endl;
      }*/
    }
    /*if (fabs(d4) > M_PI*2/7) {
      debug << "all good.." << endl;
    }*/
    //INFODebug(debug);
    return starShapeList;
  }


  //  returns the summand for the cubic interpolation in one dimension. X being the desired position between the xi. vi are the value at the positions xi!
  double
  RdStationSignalInterpolator::LagrangianFactorCubic(const double X,
                                                     const double x1, const double x2, const double x3, const double x4,
                                                     const double v1, const double v2, const double v3, const double v4)
  {
    const double value =
      (X - x2) / (x1 - x2) * (X - x3) / (x1 - x3) * (X - x4) / (x1 - x4) * v1 +
      (X - x1) / (x2 - x1) * (X - x3) / (x2 - x3) * (X - x4) / (x2 - x4) * v2 +
      (X - x1) / (x3 - x1) * (X - x2) / (x3 - x2) * (X - x4) / (x3 - x4) * v3 +
      (X - x1) / (x4 - x1) * (X - x2) / (x4 - x2) * (X - x3) / (x4 - x2) * v4;
    /*ostringstream debug;
    debug << "interpolated value: " << value << endl;
    INFODebug(debug);*/
    return value;
  }


  // definition of findIndexOfSurroundingNN, this gives the indices of the four closest stations (or rather: the two closest stations for the two closest azimuth angles)
  vector<int>
  RdStationSignalInterpolator::FindIndexOfSurroundingNN(const double radius, const double azimuth,
                                                        const std::set<double>& azimuthListStarShape,
                                                        const std::vector<StarShapeEntry>& starShapeList)
  {
    //ostringstream debug;
    bool ok = true;
    double highAzi = -10;
    double lowAzi = -10;
    double roundedAzi;
    vector<int> listNN;
    for (unsigned int tmppul = 0; tmppul < starShapeList.size()-1; ++tmppul) {
      roundedAzi = round(starShapeList.at(tmppul).fObserverAzimuthAngle*1000) / 1000;
      if (azimuthListStarShape.find(roundedAzi) == azimuthListStarShape.end())
        continue;
      if (azimuth > roundedAzi)
        lowAzi = roundedAzi;
      else {
        highAzi = roundedAzi;
        break;
      }
    }  // for loop: find the Azimuth angle of the arm on the right (HighAzi) and the left (LowAzi) of the starshaped station list.

    //debug << "Higher Azimuth:" << HighAzi << ", Lower Azimuth: " << LowAzi << endl;
    if (azimuth < *azimuthListStarShape.begin() || azimuth > *azimuthListStarShape.rbegin()) {
      lowAzi = *azimuthListStarShape.rbegin();
      highAzi = *azimuthListStarShape.begin();
    }
    if (azimuthListStarShape.size() > 8 || lowAzi == -10 || highAzi == -10)
      return { -1 };
    //debug << "Higher Azimuth:" << HighAzi << ", Lower Azimuth: " << LowAzi << endl;
    //INFODebug(debug);
    int indexHighHigh = -1;
    int indexHighLow = -1;
    int indexLowHigh = -1;
    int indexLowLow = -1;
    for (unsigned int tmppul = 0; tmppul < starShapeList.size()-1; ++tmppul) {  // find closest radii of lower Azimuth
      if ((fabs(lowAzi - starShapeList.at(tmppul).fObserverAzimuthAngle) >= 0.01) ||
          (fabs(lowAzi - starShapeList.at(tmppul+1).fObserverAzimuthAngle) >= 0.01))
        continue;

      if (starShapeList.at(tmppul).fDistanceFromShowerAxis < radius &&
          starShapeList.at(tmppul+1).fDistanceFromShowerAxis > radius) {
        indexLowLow = tmppul;
        indexLowHigh = tmppul+1;
        break;
      }
      if (tmppul == starShapeList.size()-1)
        ok = false; //No Pulses with lower Azi with according radii
    }
    if (!ok || indexLowHigh == -1)
      return { -1 };
    for (unsigned int tmppul = 0; tmppul < starShapeList.size()-1; ++tmppul) {  // find closest radii of higher Azimuth
      if ((fabs(highAzi - starShapeList.at(tmppul).fObserverAzimuthAngle) >= 0.01) ||
          (fabs(highAzi - starShapeList.at(tmppul+1).fObserverAzimuthAngle) >= 0.01))
        continue;

      if (starShapeList.at(tmppul).fDistanceFromShowerAxis < radius &&
          starShapeList.at(tmppul+1).fDistanceFromShowerAxis > radius) {
        indexHighLow = tmppul;
        indexHighHigh = tmppul+1;
        break;
      }
      if (tmppul == starShapeList.size()-1)
        ok = false; //No Pulses with higher Azi with according radii
    }
    if (!ok || indexHighHigh == -1)
      return { -1 };
    listNN = { indexHighHigh, indexHighLow, indexLowHigh, indexLowLow };
    return listNN;
  }


  vector<double>
  RdStationSignalInterpolator::NearestNeighbourInterpolation(const double radius, const double azimuth,
                                                             const std::set<double>& azimuthListStarShape,
                                                             const std::vector<StarShapeEntry>& starShapeList,
                                                             const std::vector<StationRRecDataQuantities>& eNumNames,
                                                             REvent& rEvent, const bool errors)
  {
    //ostringstream final, debug;
    double r0 = radius;
    double phi0 = azimuth;
    vector<int> nn = FindIndexOfSurroundingNN(radius, azimuth, azimuthListStarShape, starShapeList);
    vector <double> interpolatedValues;
    //debug << "NN-Size: " << NN.size() << endl;
    if (nn.size() != 4 || eNumNames.size() == 0)
      return { 0 };

    double tmpDist = PolarDistance(starShapeList.at(nn[0]).fDistanceFromShowerAxis,
                                   starShapeList.at(nn[0]).fObserverAzimuthAngle, r0, phi0);
    int nnStationID = starShapeList.at(nn[0]).fStationID;
    for (unsigned int i = 0; i < nn.size(); ++i) {
      if (nn[i] == -1)
        return { 0 };
      if (PolarDistance(starShapeList.at(nn[i]).fDistanceFromShowerAxis,
                        starShapeList.at(nn[i]).fObserverAzimuthAngle, r0, phi0) < tmpDist) {
        tmpDist = PolarDistance(starShapeList.at(nn[i]).fDistanceFromShowerAxis,
                                starShapeList.at(nn[i]).fObserverAzimuthAngle, r0, phi0);
        nnStationID = starShapeList.at(nn[i]).fStationID;
      }
      //debug << "nn ID: " << nnStationID << endl;
    }
    //debug << "nn Station ID: " << nnStationID << endl;
    //INFODebug(debug);
    if (errors) {
      for (vector<StationRRecDataQuantities>::const_iterator it = eNumNames.begin(); it != eNumNames.end(); ++it) {
        if (rEvent.GetStation(nnStationID).GetRecData().HasParameterError(*it))
          interpolatedValues.push_back(rEvent.GetStation(nnStationID).GetRecData().GetParameterError(*it));
        else {
          /*final << "eNum Error for " << *it << " not set, skip interpolation for Station" << endl;
          INFOFinal(final);*/
          return { 0 };
        }
      }
    } else {
      for (vector<StationRRecDataQuantities>::const_iterator it = eNumNames.begin(); it != eNumNames.end(); ++it) {
        if (rEvent.GetStation(nnStationID).GetRecData().HasParameter(*it))
          interpolatedValues.push_back(rEvent.GetStation(nnStationID).GetRecData().GetParameter(*it));
        else {
          /*final << "eNum for " << *it << " not set, skip interpolation for Station" << endl;
          INFOFinal(final);*/
          return { 0 };
        }
      }
    }
    return interpolatedValues;
  }


  vector<double>
  RdStationSignalInterpolator::BilinearInterpolation(const double radius, const double azimuth,
                                                     const std::set<double>& azimuthListStarShape,
                                                     const std::vector<StarShapeEntry>& starShapeList,
                                                     const std::vector<StationRRecDataQuantities>& eNumNames,
                                                     REvent& rEvent, const bool errors)
  {
    //ostringstream final,debug;
    double r0 = radius;
    double phi0 = azimuth;
    vector<int> nn = FindIndexOfSurroundingNN(radius, azimuth, azimuthListStarShape, starShapeList);
    vector<int> nStationID;
    vector <double> interpolatedValues;
    //debug << "nn-Size: " << nn.size() << endl;
    if (nn.size() != 4 || eNumNames.size() == 0)
      return { 0 };
    for (unsigned int i = 0; i < nn.size(); ++i) {
      if (nn[i] == -1)
        return { 0 };
      nStationID.push_back(starShapeList.at(nn[i]).fStationID);
      //debug << "nn: " <<nn[i] << endl;
    }
    //INFODebug(debug);
    double r1 = starShapeList.at(nn[0]).fDistanceFromShowerAxis;
    double r2 = starShapeList.at(nn[1]).fDistanceFromShowerAxis;
    double r3 = starShapeList.at(nn[2]).fDistanceFromShowerAxis;
    double r4 = starShapeList.at(nn[3]).fDistanceFromShowerAxis;
    double phi1 = starShapeList.at(nn[0]).fObserverAzimuthAngle;
    double phi2 = starShapeList.at(nn[2]).fObserverAzimuthAngle;
    if (phi1 < phi2) {
      phi1 += 2*M_PI;
      if (phi0 < 0)
        phi0 += 2*M_PI;
    }
    if (errors) {
      for (vector<StationRRecDataQuantities>::const_iterator it = eNumNames.begin(); it != eNumNames.end(); ++it) {
        if (rEvent.GetStation(nStationID[0]).GetRecData().HasParameterError(*it) &&
            rEvent.GetStation(nStationID[1]).GetRecData().HasParameterError(*it) &&
            rEvent.GetStation(nStationID[2]).GetRecData().HasParameterError(*it) &&
            rEvent.GetStation(nStationID[3]).GetRecData().HasParameterError(*it)) {
          double par1 = rEvent.GetStation(nStationID[0]).GetRecData().GetParameterError(*it);
          double par2 = rEvent.GetStation(nStationID[1]).GetRecData().GetParameterError(*it);
          double par3 = rEvent.GetStation(nStationID[2]).GetRecData().GetParameterError(*it);
          double par4 = rEvent.GetStation(nStationID[3]).GetRecData().GetParameterError(*it);
          interpolatedValues.push_back((phi2 - phi0) / (phi2 - phi1) * ((r2 - r0) / (r2-r1) * par1 + (r0 - r1) / (r2 - r1) * par2) +
                                       ((phi0 - phi1) / (phi2 - phi1) * ((r4 - r0) / (r4 - r3)*par3 + (r0 - r3) / (r4 - r3)*par4)));
        } else {
          /*final << "eNum Error for " << *it << " not set, skip interpolation for Station" << endl;
          INFOFinal(final);*/
          return { 0 };
        }
      }
    } else {
      for (vector<StationRRecDataQuantities>::const_iterator it = eNumNames.begin(); it != eNumNames.end(); ++it) {
        if (rEvent.GetStation(nStationID[0]).GetRecData().HasParameter(*it) &&
            rEvent.GetStation(nStationID[1]).GetRecData().HasParameter(*it) &&
            rEvent.GetStation(nStationID[2]).GetRecData().HasParameter(*it) &&
            rEvent.GetStation(nStationID[3]).GetRecData().HasParameter(*it)) {
          double par1 = rEvent.GetStation(nStationID[0]).GetRecData().GetParameter(*it);
          double par2 = rEvent.GetStation(nStationID[1]).GetRecData().GetParameter(*it);
          double par3 = rEvent.GetStation(nStationID[2]).GetRecData().GetParameter(*it);
          double par4 = rEvent.GetStation(nStationID[3]).GetRecData().GetParameter(*it);
          interpolatedValues.push_back((phi2 - phi0) / (phi2 - phi1) * (( r2 - r0) / (r2 - r1) * par1 + (r0 - r1) / (r2 - r1) * par2)
                  + ((phi0 - phi1) / (phi2 - phi1) * ((r4 - r0) / (r4 - r3) * par3 + (r0 - r3) / (r4 - r3) * par4)));
        } else {
          /*final << "eNum for " << *it << " not set, skip interpolation for Station" << endl;
          INFOFinal(final);*/
          return { 0 };
        }
      }
    }
    return interpolatedValues;
  }


  // bicubic interpolation at the station position
  vector<double>
  RdStationSignalInterpolator::BicubicInterpolation(const double radius, const double azimuth,
                                                    const std::set<double>& azimuthListStarShape,
                                                    const std::vector<StarShapeEntry>& starShapeList,
                                                    const std::vector<StationRRecDataQuantities>& eNumNames,
                                                    REvent& rEvent, const bool errors)
  {
    ostringstream intermediate;
    ostringstream debug;
    double r0 = radius;
    double phi0 = azimuth;
    vector<StarShapeEntry> interpolationStationList;
    vector<int> tmpNN = FindIndexOfSurroundingNN(radius, azimuth, azimuthListStarShape, starShapeList);
    vector<int> nStationID;
    if (tmpNN.size() != 4 || eNumNames.size() == 0)
      return { 0 };
    for (unsigned int i = 0; i < tmpNN.size(); ++i) {
      if (tmpNN[i] == -1)
        return { 0 };
      nStationID.push_back(starShapeList.at(tmpNN[i]).fStationID);
      //debug << "NN: " <<tmpNN[i] << endl;
    }
    double r1 = starShapeList.at(tmpNN[0]).fDistanceFromShowerAxis;
    double r2 = starShapeList.at(tmpNN[1]).fDistanceFromShowerAxis;
    double r3 = starShapeList.at(tmpNN[2]).fDistanceFromShowerAxis;
    double r4 = starShapeList.at(tmpNN[3]).fDistanceFromShowerAxis;
    double phi1 = starShapeList.at(tmpNN[0]).fObserverAzimuthAngle;
    double phi2 = starShapeList.at(tmpNN[1]).fObserverAzimuthAngle;
    double phi3 = starShapeList.at(tmpNN[2]).fObserverAzimuthAngle;
    double phi4 = starShapeList.at(tmpNN[3]).fObserverAzimuthAngle;
    //debug << "radius 1-4: " << r1 << ", " << r2 << ", " << r3 << ", " << r4 << endl;
    //debug << "phi 1-4: " << phi1 << ", " << phi2 << ", " << phi3 << ", " << phi4 << endl;

    // get the 16 closest stations
    vector<int> nnuu = FindIndexOfSurroundingNN(r1+10, phi1+0.1, azimuthListStarShape, starShapeList);
    vector<int> nnud = FindIndexOfSurroundingNN(r2-10, phi2+0.1, azimuthListStarShape, starShapeList);
    vector<int> nndu = FindIndexOfSurroundingNN(r3+10, phi3-0.1, azimuthListStarShape, starShapeList);
    vector<int> nndd = FindIndexOfSurroundingNN(r4-10, phi4-0.1, azimuthListStarShape, starShapeList);
    /*for (int i = 0; i < 4; ++i) {
      debug << "tmpNN: " << tmpNN[i] << ", " << "UU: " << NNuu[i] << ", " << "UD: "
            << NNud[i] << ", " << "DU: " << NNdu[i] << ", " << "DD: " << NNdd[i] << endl;
    }*/
    if (nnuu[1] == nnud[0] || nnuu[1] == nnud[1]) {
      nnuu = FindIndexOfSurroundingNN(r1+30, phi1+0.1, azimuthListStarShape, starShapeList);
      nnud = FindIndexOfSurroundingNN(r2-30, phi2+0.1, azimuthListStarShape, starShapeList);
      //debug << "No perfect star shape => nnuu and nnud changed!" << endl;
      if (nnuu[1] == nnud[0] || nnuu[1] == nnud[1]) {
        /*debug << "... but it didn't work." << endl;
        INFODebug(debug);*/
        return { -1 };
      }
    }
    if (nndd[2] == nndu[3] || nndd[3] == nndu[3]) {
      nndd = FindIndexOfSurroundingNN(r3+30, phi3-0.1, azimuthListStarShape, starShapeList);
      nndu = FindIndexOfSurroundingNN(r4-30, phi4-0.1, azimuthListStarShape, starShapeList);
      debug << "No perfect star shape => NNdd and NNdu changed!" << endl;
      if (nndd[3] == nndu[2] || nndd[3] == nndu[3]){
        debug << "... but it didn't work." << endl;
        INFODebug(debug);
        return{-1};
      }
    }

    nStationID = { };
    if (nnuu.size() != 4 || nnud.size() != 4 || nndu.size() != 4 || nndd.size() != 4)
      return { 0 };
    for (unsigned int i = 0; i < nnuu.size(); ++i) { // fill in grid of 16 stations
      if (nnuu[i] == -1 || nnud[i] == -1 || nndu[i] == -1 || nndd[i] == -1)
        return { 0 };
      nStationID.push_back(starShapeList.at(tmpNN[i]).fStationID);
      interpolationStationList.push_back(starShapeList.at(nnuu[i]));
      interpolationStationList.push_back(starShapeList.at(nnud[i]));
      interpolationStationList.push_back(starShapeList.at(nndu[i]));
      interpolationStationList.push_back(starShapeList.at(nndd[i]));
      //debug << StarShapeList.at(tmpNN[i]).fStationID << "  ";
    }
    sort(interpolationStationList.begin(), interpolationStationList.end());

    // now manipulate:
    vector<StarShapeEntry> interpolationStationListOrdered = SortForInterpolation(interpolationStationList);

    // now begin interpolation:
    // if first azimuth is positive, then the "sortForInterpolation" might have put the lowest angles at the end of the vector.
    // Therefor 2Pi must be added to them to account for correct "skipping" from 2Pi to -2Pi.

    phi1 = interpolationStationListOrdered[0].fObserverAzimuthAngle;
    phi2 = interpolationStationListOrdered[4].fObserverAzimuthAngle;
    phi3 = interpolationStationListOrdered[8].fObserverAzimuthAngle;
    phi4 = interpolationStationListOrdered[12].fObserverAzimuthAngle;
    if (interpolationStationListOrdered[0].fObserverAzimuthAngle > 0) {
      if (phi0 < 0)
        phi0 += 2*M_PI;
      if (phi2 < 0)
        phi2 += 2*M_PI;
      if (phi3 < 0)
        phi3 += 2*M_PI;
      if (phi4 < 0)
        phi4 += 2*M_PI;
    }

    //debug << "check azimuth differences: " << phi2-phi1 << ", " << phi3-phi2 << ", " << phi4-phi3 << ", " << phi1-phi4 << endl;
    //debug << "azimuth: " << phi0 << ", phis: " << phi1 << ", " << phi2 << ", " << phi3 << ", " << phi4 << endl;
    //INFODebug(debug);
    // now the bicubic interpolation starts:
    double rTmp[16]; // store the radius values of all 16 positions (correctly ordered azimuth!)
    double bicTmpValues[4]; // calculated lagrangian factors
    vector <double> interpolatedValues;
    if (errors) {
      for (vector<StationRRecDataQuantities>::const_iterator it = eNumNames.begin(); it != eNumNames.end(); ++it) {
        for (unsigned int i = 0; i < 16; ++i) {
          rTmp[i] = interpolationStationListOrdered[i].fDistanceFromShowerAxis;
          if (!rEvent.GetStation(interpolationStationListOrdered[i].fStationID).GetRecData().HasParameterError(*it)) {
            intermediate << "eNum Error for " << *it << " not set for at least one simulated position, skip interpolation for Station" << endl;
            INFOIntermediate(intermediate);
            return { 0 };
          }
        }
        for (unsigned int i = 0; i < 4; ++i) {
          bicTmpValues[i] =
            LagrangianFactorCubic(r0, rTmp[4*i], rTmp[4*i+1], rTmp[4*i+2], rTmp[4*i+3],
                                  rEvent.GetStation(interpolationStationListOrdered[4*i].fStationID).GetRecData().GetParameterError(*it),
                                  rEvent.GetStation(interpolationStationListOrdered[4*i+1].fStationID).GetRecData().GetParameterError(*it),
                                  rEvent.GetStation(interpolationStationListOrdered[4*i+2].fStationID).GetRecData().GetParameterError(*it),
                                  rEvent.GetStation(interpolationStationListOrdered[4*i+3].fStationID).GetRecData().GetParameterError(*it));
        }
        double iValue = LagrangianFactorCubic(phi0, phi1, phi2, phi3, phi4, bicTmpValues[0], bicTmpValues[1], bicTmpValues[2], bicTmpValues[3]);
        interpolatedValues.push_back(iValue);
      }
    } else {
      for (vector<StationRRecDataQuantities>::const_iterator it = eNumNames.begin(); it != eNumNames.end(); ++it) {
        for (unsigned int i = 0; i < 16; ++i) {
          rTmp[i] = interpolationStationListOrdered[i].fDistanceFromShowerAxis;
          if (!rEvent.GetStation(interpolationStationListOrdered[i].fStationID).GetRecData().HasParameter(*it)) {
            intermediate << "eNum for " << *it << " not set for at least one simulated position, skip interpolation for Station" << endl;
            INFOIntermediate(intermediate);
            return { 0 };
          }
        }
        for (unsigned int i = 0; i < 4; ++i) {
          bicTmpValues[i] =
            LagrangianFactorCubic(r0, rTmp[4*i], rTmp[4*i+1], rTmp[4*i+2], rTmp[4*i+3],
                                  rEvent.GetStation(interpolationStationListOrdered[4*i].fStationID).GetRecData().GetParameter(*it),
                                  rEvent.GetStation(interpolationStationListOrdered[4*i+1].fStationID).GetRecData().GetParameter(*it),
                                  rEvent.GetStation(interpolationStationListOrdered[4*i+2].fStationID).GetRecData().GetParameter(*it),
                                  rEvent.GetStation(interpolationStationListOrdered[4*i+3].fStationID).GetRecData().GetParameter(*it));
        }
        double iValue = LagrangianFactorCubic(phi0, phi1, phi2, phi3, phi4, bicTmpValues[0], bicTmpValues[1], bicTmpValues[2], bicTmpValues[3]);
        interpolatedValues.push_back(iValue);
        /*debug.str("");
        debug << "interpolated Value: " << iValue;*/
      }
    }
    //INFODebug(debug);
    return interpolatedValues;
  }

}