RdStationBeamFormer.cc

Go to the documentation of this file.

#include "RdStationBeamFormer.h"

#include <fwk/CentralConfig.h>
#include <fwk/VModule.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
// #include <utl/config.h>
// #include <utl/TraceAlgorithm.h>
// #include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>

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

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>

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

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

namespace RdStationBeamFormer {

  RdStationBeamFormer::RdStationBeamFormer():
  fInfoLevel(0),
  fWaveFront(0),
  fRefractionIndex(1),
  fSpeedOfLight(kSpeedOfLight)
  {
  }


  RdStationBeamFormer::~RdStationBeamFormer()
  {
  }


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


    // Read in the configurations of the xml file
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdStationBeamFormer");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("WaveFront").GetData(fWaveFront);
    topBranch.GetChild("RefractionIndex").GetData(fRefractionIndex);
    fSpeedOfLight = kSpeedOfLight/fRefractionIndex;

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationBeamFormer::Run(evt::Event& event)
  {
    if (fInfoLevel >= eObscure)
      INFO("RdStationBeamFormer::Run()");

    // Check if there is an radio event at all
    if(!event.HasREvent()) {
      if (fInfoLevel >= eFinal)
        WARNING("No radio event found!");
      return eContinueLoop;
    }

    // get the event and the detector (for antenna positions)
    REvent& rEvent = event.GetREvent();

    // Check if there are at least 2 stations (otherwise beam forming does not make sense)
    if(rEvent.GetNumberOfStations() < 2) {
      if (fInfoLevel >= eFinal)
        WARNING("At least 2 stations are needed for beam forming!");
      return eContinueLoop;
    }
    // get detector description (for antenna positions)
    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();

    // Look if reconstruction information is available
    if (!event.HasRecShower()) {
      ERROR("No shower reconstruction found, but needed as input.");
      return eContinueLoop;
    }
    const ShowerRecData& showerrec = event.GetRecShower();
  
    // for now, work only with the RD reconstruction
    // later also input from FD or SD could be used to get core and direction for the beam forming
    if (!showerrec.HasRRecShower()) {
      ERROR("No radio reconstruction (e.g. from RDPlaneFit) found, but needed as input.");
      return eContinueLoop;
    }
    const ShowerRRecData& showerrrec = showerrec.GetRRecShower();

    // get core position and axis
    const Point& core = showerrrec.GetCorePosition();
    const Vector& axis = showerrrec.GetAxis();

    stringstream fMessage;
     fMessage << "Output of core and axis: (work in progress) ";
    if (fInfoLevel >= eDebug)
      INFO(fMessage.str());

    // store time stamp of the first station as reference:
    // the time stamps of different stations might be different and this "additional delay"
    // has to be taken into account for the beam forming
    TimeStamp firstTimeStamp  = rEvent.StationsBegin()->GetTraceStartTime();

    // loop through all stations and calculate the delays for beam forming
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
           sIt != rEvent.StationsEnd(); ++sIt) {
      // get the time stamp and store eventual differences to correct for them later
      Station& station = *sIt;
      TimeStamp stationTimeStamp = station.GetTraceStartTime();
      TimeInterval timeStampDelay = stationTimeStamp-firstTimeStamp;

      // get the position of each station and calculate the distance from the shower plane
      // distance to shower plane = (antenna position - core) * (axis)
      const rdet::Station& dStation = rDetector.GetStation(*sIt);
      const Vector vStationPosition = (dStation.GetPosition() - core);
      const double distance = vStationPosition * axis;

      // convert the distance in a geometrical delay to shift the data
      const TimeInterval geomDelay = distance / fSpeedOfLight;
      fMessage.str("");
      fMessage << "Distance of Station: "
               << distance/meter
               << " m; \t geom. Delay: "
               << geomDelay/nanosecond
               << " ns; \t TimeStampDelay: "
               << timeStampDelay/nanosecond
               << " ns";
      if (fInfoLevel >= eDebug)
        INFO(fMessage.str());

      const TimeInterval totalDelay =  timeStampDelay + geomDelay;    
    }

    return eSuccess;
  }


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

    if (fInfoLevel >= eFinal)
      INFO("RdStationBeamFormer::Finish()");

    return eSuccess;
  }

  void 
  RdStationBeamFormer::timeShift(revt::Station& station,
                          const utl::TimeInterval& shift) const
  {
    // Get the frequency spectrum
    revt::StationFrequencySpectrum& spectrum =
      station.GetStationFrequencySpectrum();

    // multiply phase gradient to each value of the spectrum
    for (StationFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
      spectrum[i] *=  
        exp(- complex<double>(0,1) * kTwoPi * station.GetFrequencyOfBin(i) * shift.GetInterval());
  }


}