RdBeamTimeOptimizer.cc

Go to the documentation of this file.

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

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/BeamQuantities.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.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/Triple.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/CoordinateSystemPtr.h>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>

#include <Minuit2/MnMinimize.h>
#include <Minuit2/VariableMetricMinimizer.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnPrint.h>

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


//Just an abbreviation of repeated code from UserModule.cc skeleton
#define PRINT(...) {stringstream msg; msg << __VA_ARGS__; INFO(msg.str());}

namespace RdBeamTimeOptimizer {
  RdBeamTimeOptimizer::RdBeamTimeOptimizer()
  {
  }

  RdBeamTimeOptimizer::~RdBeamTimeOptimizer()
  {
  }

  VModule::ResultFlag
  RdBeamTimeOptimizer::Init()
  {
    INFO("RdBeamTimeOptimizer::Init()");

    Branch conf =
      CentralConfig::GetInstance()->GetTopBranch("RdBeamTimeOptimizer");

    conf.GetChild("startBin").GetData(startbin);
    conf.GetChild("stopBin").GetData(stopbin);

    string confstring;
    conf.GetChild("waveModel").GetData(confstring);
    waveModel = (confstring == "ePlane") ? (RdWaveModel*) new RdPlaneWaveModel() :
      (RdWaveModel*) new RdSphericalWaveModel();

//  const Branch txtOut = conf.GetChild("asciiOutput");
//     if (txtOut) {
//       txtOut.GetChild("outputFileName").GetData(outfile);
//       txtOut.GetChild("outputGridType").GetData(confstring);
//       gridtype = (confstring == "ePolar") ? polar :
//         (confstring == "eCarthesian") ? carthesian : none;
//     } else {
//       gridtype = none; outfile = "";
//     }

    return eSuccess;
  }

  void RdBeamTimeOptimizer::shiftsToDelays(REvent& rEvent, TimeStamp meanTime, vector<double> &times) {
    if ((uint)rEvent.GetNumberOfStations() != times.size()) {
      throw utl::OutOfBoundException("Wrong number of delays supplied.");}

    REvent::StationIterator stat = rEvent.StationsBegin();
    vector<double>::iterator time = times.begin();
    for (; stat != rEvent.StationsEnd(); ++stat, ++time) {
      (*time) += ( stat -> GetTraceStartTime() - meanTime).GetInterval();
    }
  }

  void RdBeamTimeOptimizer::delaysToShifts(REvent& rEvent, TimeStamp meanTime, vector<double> &times) {
    if ((uint)rEvent.GetNumberOfStations() != times.size()) {
      throw utl::OutOfBoundException("Wrong number of timeshifts supplied.");}

    REvent::StationIterator stat = rEvent.StationsBegin();
    vector<double>::iterator time = times.begin();
    for (; stat != rEvent.StationsEnd(); ++stat, ++time) {
      (*time) -= ( stat -> GetTraceStartTime() - meanTime ).GetInterval();
    }
  }

  vector<double> RdBeamTimeOptimizer::delaysFromRRec(REvent& rEvent,  evt::ShowerRRecData& rrec){
    //#warning TH: replace Core Position with other reference?
    Point skyPos = rrec.GetCorePosition() + rrec.GetAxis();
    waveModel->setOrigin(rrec.GetCorePosition());
    waveModel->setSkyPos (skyPos);
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    vector<double> delays;
    REvent::StationIterator stat = rEvent.StationsBegin();

    for (; stat != rEvent.StationsEnd(); ++stat) {
      const Point pos = rDetector.GetStation(*stat).GetPosition();
      delays.push_back( waveModel->delay(pos) );
    }
    return delays;
  }

  vector<double> RdBeamTimeOptimizer::shiftsFromRRec(REvent& rEvent,  evt::ShowerRRecData& rrec){
    //#warning TH: replace Core Position with other reference?
    Point skyPos = rrec.GetCorePosition() + rrec.GetAxis();
    waveModel->setOrigin(rrec.GetCorePosition());
    waveModel->setSkyPos (skyPos);
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    vector<double> shifts;
    REvent::StationIterator stat = rEvent.StationsBegin();

    double middleOfTrace;
    {
      int start = stat->GetConstStationTimeSeries().GetStart(),
        stop = stat->GetConstStationTimeSeries().GetStop();
      if ((startbin != -1) && (startbin > start)) start = startbin;
      if ((stopbin != -1) && (stopbin < stop)) stop = stopbin;
      middleOfTrace = (stop+start) * stat->GetConstStationTimeSeries().GetBinning() /2.0;
    }

    //#warning TH: TraceStartTime to be replaced by new logic (ParameterStorage), replace core as reference?
    for (; stat != rEvent.StationsEnd(); ++stat) {
      const Point pos = rDetector.GetStation(*stat).GetPosition();
      shifts.push_back( waveModel->delay(pos) - (stat -> GetTraceStartTime() - rrec.GetCoreTime()).GetInterval() + middleOfTrace );
    }
    return shifts;
  }

  // this helper function is not functionally closed. it requires an out trace of the right
  // dimension. a more selfcontained design would require more allocating and work.
  // this is really just a helper to make beamtrace more readable, so don't overengineer.
  inline void RdBeamTimeOptimizer::efieldtopower (const StationTimeSeries& in, TraceD& out, int start, int stop) const {
    StationTimeSeries::ConstIterator it = in.Begin() + start;
    for (int i = 0; i < (stop-start); i++, it++)
      out[i] = it->GetMag2();
  }

  TraceD RdBeamTimeOptimizer::powertrace (REvent& rEvent,  evt::ShowerRRecData& /*rrec*/,
                                        const vector<double>& shifts) const {

    static FFTDataContainerAlgorithm fft;

    //#warning TH: Code assumes that all traces have same binning and size!

    // trust that all traces are of same binning&size; if not, we are bound to fail anyway,
    // if traces have to be adjusted, this could be handled elsewhere beforehand.
    REvent::StationIterator it = rEvent.StationsBegin();

    int start = it->GetConstStationTimeSeries().GetStart(),
      stop = it->GetConstStationTimeSeries().GetStop();
    if ((startbin != -1) && (startbin > start)) start = startbin;
    if ((stopbin != -1) && (stopbin < stop)) stop = stopbin;

    TraceD sum(stop-start, it->GetConstStationTimeSeries().GetBinning()),
      addend(stop-start, it->GetConstStationTimeSeries().GetBinning());

    vector<double>::const_iterator shiftIt = shifts.begin();
    uint num = 0;
    for (; it != rEvent.StationsEnd(); ++it, ++shiftIt) {
      StationFFTDataContainer data = it->GetFFTDataContainer();

      const TimeInterval shift = (*shiftIt);

      fft.ShiftTimeSeries(data, -shift);

      const StationTimeSeries& ts = data.GetTimeSeries();
      efieldtopower (ts, addend, start, stop);
      sum += addend;
      ++ num;
    }

    return sum * (1.0 / num) ; //sqrt (sum/rEvent.GetNumberOfStations())
  }




  // this helper function is not functionally closed. it requires an out trace of the right
  // dimension. a more selfcontained design would require more allocating and work.
  // this is really just a helper to make beamtrace more readable, so don't overengineer.
  inline void RdBeamTimeOptimizer::efieldproduct (const StationTimeSeries& in1, const StationTimeSeries& in2, TraceD& out, int start, int stop) const {
    StationTimeSeries::ConstIterator it1 = in1.Begin() + start, it2 = in2.Begin() + start;
    uint n = it1->GetSize(); //should be 3 ;-)
    for (int i = 0; i < (stop-start); ++i, ++it1, ++it2) {
      double prod = 0.0;
      for (uint j = 0; j < n; ++j)
        prod += it1->At(j) * it2->At(j);
      out[i] = prod;
    }
  }
  TraceD RdBeamTimeOptimizer::crosscorr (REvent& rEvent, evt::ShowerRRecData& /*rrec*/,
                                        const vector<double>& shifts) const {
    static FFTDataContainerAlgorithm fft;

    // trust that all traces are of same binning&size; if not, we are bound to fail anyway,
    // if traces have to be adjusted, this could be handled elsewhere beforehand.
    REvent::StationIterator it = rEvent.StationsBegin();

    int start = it->GetConstStationTimeSeries().GetStart(),
      stop = it->GetConstStationTimeSeries().GetStop();
    if ((startbin != -1) && (startbin > start)) start = startbin;
    if ((stopbin != -1) && (stopbin < stop)) stop = stopbin;

    TraceD sum(stop-start, it->GetConstStationTimeSeries().GetBinning()),
      addend(stop-start, it->GetConstStationTimeSeries().GetBinning());

    uint n = rEvent.GetNumberOfStations();
    StationFFTDataContainer stationData[n];
    vector<double>::const_iterator shiftIt = shifts.begin();
    for (uint i=0; it != rEvent.StationsEnd(); ++it, ++shiftIt, ++i) {
      stationData[i] = it->GetFFTDataContainer();

      TimeInterval shift = (*shiftIt);
      fft.ShiftTimeSeries(stationData[i], -shift);
     }

    uint num = 0;
    for (uint i=0; i < n; ++i)
      for (uint j=i+1; j < n; ++j) {
        efieldproduct (stationData[i].GetTimeSeries(),stationData[j].GetTimeSeries(), addend, start, stop);
        sum += addend;
        ++num;
      }

    return sum * (1.0/num);
  }

  TraceD
  RdBeamTimeOptimizer::xtrace (TraceD ccTrace, double ccMean, TraceD pwNormal) {
    TraceD result (ccTrace);
    TraceD::Iterator it = result.Begin();
    TraceD::ConstIterator itPW = pwNormal.Begin();
    for (;it != result.End(); it++, itPW++) {
        (*it) *= abs ( (*it - ccMean) / (*itPW) );
      }
    return result;
  }


  /* findPeak should probably use a very fancy algorithm.
     For now a simple max-finder should do the trick.
     There is no model for the peak shape anyway. */
  evt::BeamPeak
  RdBeamTimeOptimizer::findPeak (const TraceD& trace) {
    TraceD::ConstIterator max = max_element( trace.Begin(), trace.End() );
    const double maxval = *max;
    const int maxi = (max - trace.Begin());
    const double maxtime = maxi * trace.GetBinning();

    static TraceAlgorithm algo;
    //#warning TH: replace call to RMS with RootMeanSquare everywhere
    const double rms = algo.RMS ( trace, trace.GetStart(), trace.GetStop() -1 );
    const double offset = algo.Mean ( trace, trace.GetStart(), trace.GetStop() -1);

    peakTime = maxtime;
    return evt::BeamPeak ( maxval, maxtime, rms, offset);
  }

  void
  RdBeamTimeOptimizer::writeRecData(REvent &rEvent, std::vector<double> &shifts){

    std::vector<double>::const_iterator shift = shifts.begin();
    for (REvent::StationIterator stat = rEvent.StationsBegin();
        stat != rEvent.StationsEnd(); ++stat, ++shift){
      if (!stat->HasRecData()) {stat->MakeRecData();}
      stat->GetRecData().SetSignalTime(stat -> GetTraceStartTime() + ((*shift) + peakTime));
      //#warning TH: change to new logic of signal times (ParameterStorage)
    }

  }


  VModule::ResultFlag
  RdBeamTimeOptimizer::Run(evt::Event& event)
  {
    DEBUGLOG("RdBeamTimeOptimizer::Run()");

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

    revt::REvent &rEvent = event.GetREvent();

    stringstream fMessage;
    fMessage << "Radio event found with "
        << rEvent.GetNumberOfStations()
        << " stations!";

    INFO(fMessage.str());

    // Read event header
    const Header& rHeader = rEvent.GetHeader();
    fMessage.str("");
    fMessage << "Header ID "
        << rHeader.GetId();
    //             << " and timestamp: "
    //             << rHeader.GetTime();
    INFO(fMessage.str());

    const uint nStat = rEvent.GetNumberOfStations();
    // Unfortunately the "simple" minuit API does not work for me ...
    //    std::vector <double> startTimes (nStat,0.0);
    //    std::vector <double> startErrors (nStat, 100.0 * nanosecond);

    using namespace ROOT::Minuit2;
    MnUserParameters upar;

    if (!event.HasRecShower() || !event.GetRecShower().HasRRecShower()) {
      WARNING("No RRecShower found!");
      WARNING("Falling back to: All time shifts 0.0 initially." );

      for (int i=0; i<nStat; i++)
        { stringstream parname; parname << "shift " <<i;
        upar.Add(parname.str(), 0.0, 10000.0 * nanosecond);
        }
      //      upar.Fix(0);
    }
    else {
      INFO("RRecShower found.\nCalculating start time shift values from wave model "
        "with Axis and CorePosition.");

      evt::ShowerRRecData &rrec = event.GetRecShower().GetRRecShower();
      std::vector <double> startTimes = shiftsFromRRec(rEvent, rrec);

      int i=0;
      for (std::vector<double>::const_iterator it=startTimes.begin();
        it!=startTimes.end();it++,i++)
        { stringstream parname; parname << "shift " <<i;
        upar.Add(parname.str(), (*it), 10000.0 * nanosecond);
        }
    }

    {
      stringstream msg;
      msg<< "Start parameters (time shifts): ";
      vector<double> params = upar.Params();
       copy (params.begin(),
             params.end(),
             ostream_iterator<double> (msg, "\t"));
      INFO (msg.str());
    }

    INFO ("Maximizing power peak ...");
    {
      powerFCN fcn (&rEvent,this);
      MnMinimize m( static_cast<const FCNBase&>(fcn), upar);
      FunctionMinimum fmin = m(10000);

      upar = fmin.UserParameters();
      cout << fmin;
    }

    {
      stringstream msg;
      msg<< "Shifts after optimizing power peak: ";
      vector<double> params = upar.Params();
      copy (params.begin(),
        params.end(),
        ostream_iterator<double> (msg, "\t"));
      INFO (msg.str());
    }

    INFO ("Maximizing cross corellation peak");
    {
      crossFCN fcn (&rEvent,this);
      MnMinimize m( static_cast<const FCNBase&>(fcn), upar);
      FunctionMinimum fmin = m(10000);

      upar = fmin.UserParameters();
    }

    {
      stringstream msg;
      msg<< "Shifts after optimizing cross correlation peak: ";
      vector<double> params = upar.Params();
      copy (params.begin(),
        params.end(),
        ostream_iterator<double> (msg, "\t"));
      INFO (msg.str());

      INFO ("Writing Params to StationRecData.");
      writeRecData(rEvent, params);
    }

    //    MnMinimize m( static_cast<const FCNBase&>(*this), startTimes, startErrors, 1);

    //    FunctionMinimum fmin = m(10000);

    //    cout << fmin;

//     DEBUGLOG ("calculating beam trace");

//     TraceD cc = crosscorr (rEvent, rrec);
//     TraceD power = powertrace (rEvent, rrec);

//     evt::BeamPeak ccPeak = findPeak (cc);
//     evt::BeamPeak pwPeak = findPeak (power);

//     TraceD x = xtrace(cc, ccPeak.GetOffset(), power-pwPeak.GetOffset());
//     evt::BeamPeak xPeak = findPeak (x);

//     CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create( rrec.GetCorePosition() );
//     if (gridtype == polar) {
//       const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
//       WriteASCIIPolar (cc, power, x, ccPeak, pwPeak, xPeak, polarCoords);
//     }
//     else if (gridtype == carthesian) {
//       const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
//       WriteASCII_XYZ (cc, power, x, ccPeak, pwPeak, xPeak, xyz);
//     }

//     if (rrec.HasCurrentBeamMap()) {
//       DEBUGLOG ("filling BeamQuantities into map");

//       evt::BeamQuantities quants (rrec.GetAxis().GetTheta(coreCS) ,
//                                   rrec.GetAxis().GetPhi(coreCS) ,
//                                   waveModel->curvature(rrec.GetAxis(),coreCS),
//                                   waveModel->coneAngle(rrec.GetAxis(),coreCS),
//                                   ccPeak, pwPeak, xPeak);

//       evt::BeamMap map = rrec.GetCurrentBeamMap();
//       map.AddQuantities (quants);
//     }
//     if (!rrec.HasBeamResultQuantities()
//         || (rrec.GetBeamResultQuantities().GetCCPeak().GetHeight() / rrec.GetBeamResultQuantities().GetCCPeak().GetRMS()
//             < ccPeak.GetHeight() / ccPeak.GetRMS() )) {
//       evt::BeamQuantities resultQuants (rrec.GetAxis().GetTheta(coreCS) ,
//                                         rrec.GetAxis().GetPhi(coreCS) ,
//                                         waveModel->curvature(rrec.GetAxis(),coreCS),
//                                         waveModel->coneAngle(rrec.GetAxis(),coreCS),
//                                         ccPeak, pwPeak, xPeak);
//       rrec.SetBeamResultQuantities( resultQuants );
//     }

    return eSuccess;
  }

  VModule::ResultFlag
  RdBeamTimeOptimizer::Finish()
  {
    INFO("RdBeamTimeOptimizer::Finish()");
    delete (waveModel);

    return eSuccess;
  }

  void
  RdBeamTimeOptimizer::WriteASCIISpectrum(const revt::Station& station,
                                const string& filename) const
  {
    // Get the frequency spectrum
    const revt::StationFrequencySpectrum& spectrum =
      station.GetConstStationFrequencySpectrum();

    // open file for output
    ofstream outfile;
    outfile.open(filename.c_str(),ios::out);

    // write the frequency axis and the amplitudes to the file
    for (StationFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
      outfile << station.GetFrequencyOfBin(i)/megahertz
              << " "
              << abs(spectrum[i][0])/((micro*volt)/(meter*megahertz))
              << " "
              << abs(spectrum[i][1])/((micro*volt)/(meter*megahertz))
              << " "
              << abs(spectrum[i][2])/((micro*volt)/(meter*megahertz))
              << endl;

    outfile.close();
  }

  void RdBeamTimeOptimizer::WriteASCIIPolar (const TraceD &power, const TraceD &cc, const TraceD &x,
                                const evt::BeamPeak &pwPeak, const evt::BeamPeak & ccPeak, const evt::BeamPeak &xPeak,
                                const utl::Triple &coords) const {

    utl::Triple rThetaPhi = boost::make_tuple( coords.get<0>() / km, coords.get<1>() / degree, coords.get<2>() / degree );

    WriteTrace3D((power-pwPeak.GetOffset()) / pwPeak.GetRMS(),
                "pw_"+outfile, rThetaPhi);

    WriteTrace3D((cc-ccPeak.GetOffset()) / ccPeak.GetRMS(),
                "cc_"+outfile, rThetaPhi);

    WriteTrace3D((x-xPeak.GetOffset()) / xPeak.GetRMS(),
                "x_"+outfile, rThetaPhi);

    // uncomment to write out RMS values:
    // writeValue3D(ccPeak.GetRMS()
    //                     / ((micro*volt/meter)*(micro*volt/meter)),
    //                     "cc_RMS.dat", rThetaPhi);
    // writeValue3D(pwPeak.GetRMS()
    //                     /((micro*volt/meter)*(micro*volt/meter)),
    //                     "pw_RMS.dat", rThetaPhi);
    // writeValue3D(xPeak.GetRMS()
    //                     /((micro*volt/meter)*(micro*volt/meter)),
    //                     "x_RMS.dat", rThetaPhi);
  }

  void RdBeamTimeOptimizer::WriteASCII_XYZ (const TraceD &power, const TraceD &cc, const TraceD &x,
                                const evt::BeamPeak &pwPeak,
                                const evt::BeamPeak &ccPeak,
                                const evt::BeamPeak &xPeak,
                                const utl::Triple &coords) const {

    utl::Triple xyz = boost::make_tuple( coords.get<0>() / km, coords.get<1>() / km, coords.get<2>() / km );

    WriteTrace3D((power /*-pwPeak.GetOffset()*/ ) / pwPeak.GetRMS(),
                "pw_"+outfile, xyz);

    WriteTrace3D((cc /*-ccPeak.GetOffset()*/) / ccPeak.GetRMS(),
                "cc_"+outfile, xyz);

    WriteTrace3D((x /*-xPeak.GetOffset()*/) / xPeak.GetRMS(),
                "x_"+outfile, xyz);
  }

  void
  RdBeamTimeOptimizer::WriteTrace3D(const TraceD& trace,
                        const string& filename,
                        const utl::Triple &coords) const
  {
    // open file for output
    ofstream outfile;
    outfile.open(filename.c_str(),ios::out|ios::app);

    // write the time axis and the amplitudes to the file
    for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i)
      outfile << coords.get<0>() << " "
        << coords.get<1>() << " "
        << coords.get<2>() << " "
        << trace.GetBinning()*static_cast<double>(i+startbin)/nanosecond
              << " "
              << trace[i]
              << endl;
    outfile << endl;

    outfile.close();
  }

  void writeValue3D (double val, const string& filename,
                const utl::Triple &coords) {
    ofstream outfile;
    outfile.open(filename.c_str(),ios::out|ios::app);
    outfile << coords.get<0> () << " "
        << coords.get<1> () << " "
                << coords.get<2> () << " "
        << val << endl;
    outfile.close();
  }


 // void
 //  RdBeamTimeOptimizer::WriteTracePolar(const TraceD& trace,
 //                                 const string& filename,
 //                                 const Vector& skyVec ) const
 //  {
 //    // open file for output
 //    ofstream outfile;
 //    outfile.open(filename.c_str(),ios::out|ios::app);

 //    // write the time axis and the amplitudes to the file
 //    for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i)
 //      outfile << skyVec.GetTheta(skyVec.GetCoordinateSystem()) / M_PI * 180 << " "
 //               << skyVec.GetPhi(skyVec.GetCoordinateSystem()) / M_PI * 180 << " "

 //              << " "
 //              << trace[i] // / ((micro*volt/meter)*(micro*volt/meter))
 //              << endl;
 //    outfile << endl;

 //    outfile.close();
 //  }

}

#undef PRINT