RdBeamTimeOptimizer/RdWaveModel.cc

Go to the documentation of this file.

#include "RdWaveModel.h"
//#include "revt+points.h"
#include <utl/ErrorLogger.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>

#include <det/Detector.h>
#include <rdet/RDetector.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/PhysicalConstants.h>

#include <fwk/CoordinateSystemRegistry.h>
#include <cmath>

using namespace std;

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

namespace RdBeamTimeOptimizer {

  RdWaveModel::RdWaveModel()  { }
  RdWaveModel::~RdWaveModel()  { }

  static inline double norm (utl::Vector x) {
    return sqrt(x*x);
  }

  void RdWaveModel::setOrigin (utl::Point origin) {
    this->origin = origin;
  }

  double RdWaveModel::delay(utl::Point /*skyPos*/,
                            utl::Point /*antennaPos*/) const {
    cerr << "Abstract method RdWaveModel::delay not meant to be called\n";
    return NAN;
  }
  void RdWaveModel::setSkyPos(utl::Point skyPos) {
    this -> skyPos = skyPos;
  }
  double RdWaveModel::delay (utl::Point antennaPos) const {
    return delay (skyPos, antennaPos);
  }
  
  double RdSphericalWaveModel::delay(utl::Point skyPos,
                                     utl::Point antennaPos) const {
    const utl::Vector diff = skyPos - antennaPos;
    return (norm(diff) - norm(skyPos-origin)) / kSpeedOfLight;
  }
  void RdSphericalWaveModel::setSkyPos(utl::Point skyPos) {
    this -> skyPos = skyPos;
    skyVecNorm = norm(skyPos - origin);
  }
  double RdSphericalWaveModel::delay (utl::Point antennaPos) const {
    const utl::Vector diff = antennaPos - skyPos;
    return (norm(diff) - skyVecNorm) / kSpeedOfLight;
  }
    
  double RdPlaneWaveModel::delay(utl::Point skyPos,
                                 utl::Point antennaPos) const {
    const utl::Vector rho = skyPos - origin;
    const utl::Vector x = antennaPos - origin;
    return -(rho * x) / norm (rho) / kSpeedOfLight;
  }
  
  double RdWaveModel::coneAngle(utl::Vector /*skyVec*/,
                                const utl::CoordinateSystemPtr /*cs*/) {
    return 0.0;
  }

  double RdWaveModel::curvature(utl::Vector /*skyVec*/,
                                const utl::CoordinateSystemPtr /*cs*/) {
    return 0.0;
  }

  double RdSphericalWaveModel::curvature(utl::Vector skyVec,
                                         const utl::CoordinateSystemPtr cs) {
    return 1.0 / skyVec.GetR(cs);
  }
  
  // From here on only module code for testing ...
#ifdef RDDELAYS_AS_MODULE
  
  VModule::ResultFlag
  RdWaveModel::Init() {
    INFO("RdWaveModel::Init()");
    return eSuccess;
  }


  VModule::ResultFlag
  RdWaveModel::Run(evt::Event& event) {
    INFO("RdWaveModel::Run()");
    
    if(!event.HasREvent())
    {
      WARNING("RdWaveModel::No radio event found!");
      return eContinueLoop;
    }
    
    REvent& rEvent = event.GetREvent();
    
    stringstream fMessage;
    fMessage << "Radio event found with "
             << rEvent.GetNumberOfStations()
             << " stations!";
    INFO(fMessage.str());

    const Header& rHeader = rEvent.GetHeader();
    fMessage.str("");
    fMessage << "Found header with ID "
             << rHeader.GetId()
             << " and timestamp: "
             << rHeader.GetTime();
    INFO(fMessage.str());
    
    setOrigin (revt::signalBarycenter (rEvent));

    return eSuccess;
  }

  VModule::ResultFlag
  RdWaveModel::Finish()  {
    INFO("RdWaveModel::Finish()");
    return eSuccess;
  }

#endif
}