RdHyperbolicWavefrontFit.cc

Go to the documentation of this file.

// Author: Qader Dorosti
// email: qader.dorosti@kit.edu

#include "RdHyperbolicWavefrontFit.h"
#include "DetectorGeometry.h"
#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <utl/Branch.h>
#include <utl/String.h>
#include <utl/RadioGeometryUtilities.h>

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

#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/StationConstants.h>
#include <revt/StationRecData.h>
#include <revt/StationConstants.h>

#include <utl/TraceAlgorithm.h>

#include <sstream>
#include <algorithm>
#include <TMinuit.h>
#include <TRandom.h>

#include <math.h>

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


namespace RdHyperbolicWavefrontFit {

  // global variables
  CoordinateSystemPtr RdHyperbolicWavefrontFit::fgLocalCS;
  utl::Vector RdHyperbolicWavefrontFit::gSAxis;
  utl::Point RdHyperbolicWavefrontFit::gCore;
  REvent* RdHyperbolicWavefrontFit::fgCurrentREvent = nullptr;
  Event RdHyperbolicWavefrontFit::fgCurrentEvent;
  unsigned int RdHyperbolicWavefrontFit::nPar = 0;
  double RdHyperbolicWavefrontFit::gTime0 = 0;
  unsigned int RdHyperbolicWavefrontFit::gStation0Id = 0;


  VModule::ResultFlag
  RdHyperbolicWavefrontFit::Init()
  {
    INFODebug("RdHyperbolicWavefrontFit::Init()");

    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdHyperbolicWavefrontFit");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("UsedCore").GetData(fUsedCore);
    topBranch.GetChild("DoCoreSearch").GetData(fDoCoreSearch);

    return eSuccess;
  }


  bool
  Sorter(const std::vector<double>& x, const std::vector<double>& y)
  {
    return x[7] < y[7];
  }


  bool
  SorterSt(const pair<double, int>& x, const pair<double, int>& y)
  {
    return x.first < y.first;
  }


  VModule::ResultFlag
  RdHyperbolicWavefrontFit::Run(Event& event)
  {
    INFODebug("RdHyperbolicWavefrontFit::Run()");

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

    fgCurrentEvent = event;
    fgCurrentREvent = &event.GetREvent();

    if (!event.GetRecShower().HasRRecShower()) {
      event.GetRecShower().MakeRRecShower();
    }

    if (fUsedCore == "MC"){
      ShowerSimData& MCShower = event.GetSimShower();
      // fill the minute gloabal variables
      gCore = MCShower.GetPosition();
      gSAxis = MCShower.GetDirection();
      fgLocalCS = fwk::LocalCoordinateSystem::Create(gCore);
    }

    if (fUsedCore == "SD") {
      ShowerSRecData& showerSdRec = event.GetRecShower().GetSRecShower();
      // fill the minute gloabal variables
      gCore = showerSdRec.GetCorePosition();
      gSAxis = - showerSdRec.GetAxis();
      fgLocalCS = fwk::LocalCoordinateSystem::Create(gCore);
    }

    if (fUsedCore == "RD") {
      ShowerRRecData& showerRdRec = event.GetRecShower().GetRRecShower();
      // fill the minute gloabal variables
      gCore = showerRdRec.GetCorePosition();
      gSAxis = - showerRdRec.GetAxis();
      fgLocalCS = fwk::LocalCoordinateSystem::Create(gCore);
    }

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    //vector to store the station to find the earliest hit
    vector<pair<double, int> > stationvec;

    //loop to find the earliest station registering signal
    //and use the station as the reference station
    for (REvent::SignalStationIterator sIt = fgCurrentREvent->SignalStationsBegin();
         sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
      const double time = sIt->GetRecData().GetParameter(revt::eSignalTime);
      stationvec.push_back(make_pair(time, sIt->GetId()));
    }

    //skip event with zero stations
    if (stationvec.size() == 0){
      return eContinueLoop;
    }

    sort(stationvec.begin(), stationvec.end(), SorterSt);

    gTime0 = stationvec[0].first;
    gStation0Id = stationvec[0].second;

    //vector to store the station to find the closest one to the shower axis
    vector<pair<double, int> > stationVector;

    //loop to find the nearest station to the shower axis
    //and use the station as the reference station
    for (REvent::SignalStationIterator sIt = fgCurrentREvent->SignalStationsBegin();
         sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
      const Point& stationPosition = rDetector.GetStation(*sIt).GetPosition();
      const double distance = utl::RadioGeometryUtilities::GetDistanceToAxis(gSAxis, gCore, stationPosition);
      stationVector.push_back(make_pair(distance, sIt->GetId()));
    }

    sort(stationVector.begin(), stationVector.end(), SorterSt);

    const Point posSt = rDetector.GetStation(stationVector[0].second).GetPosition();
    const utl::CoordinateSystemPtr refStCS = fwk::LocalCoordinateSystem::Create(posSt);

    // define the centre of the core search area
    const Point core(gCore.GetX(refStCS), gCore.GetY(refStCS), 0, refStCS);
    gCore = core;
    utl::CoordinateSystemPtr FixedfgLocalCS = fwk::LocalCoordinateSystem::Create(gCore);

    //generate random numbers
    TRandom rnd;

    double r = 0;
    int Ncores = 1;
    //get the core uncertainities
    if (fDoCoreSearch && fUsedCore == "SD") {
      const ShowerSRecData& showerSdRec = event.GetRecShower().GetSRecShower();
      const double CoreErr_X = showerSdRec.GetCoreError().GetX(fgLocalCS);
      const double CoreErr_Y = showerSdRec.GetCoreError().GetY(fgLocalCS);
      r = 2 * sqrt(CoreErr_X*CoreErr_X + CoreErr_Y*CoreErr_Y);
      Ncores = 3000;
    }

    //vector storing the intermediate results
    std::vector<std::vector<double>> result;

    FitParameters fit;

    //if the SD core/direction is used,
    //generate 3000 random cores
    //and fit wavefront for each core
    for (int i = 0; i < Ncores; ++i) {
      const double phi = rnd. Uniform(0, TMath::TwoPi());
      const double radius = sqrt(rnd.Uniform(0, r*r));
      const double x = radius*TMath::Cos(phi);
      const double y = radius*TMath::Sin(phi);

      const Point core(x, y, 0, FixedfgLocalCS);

      // initialise the fit parameters
      fit.gamma = 100;
      fit.t0 = 0;
      fit.NDF = 0;
      fit.u = gSAxis.GetTheta(fgLocalCS);
      fit.v = gSAxis.GetPhi(fgLocalCS);

      //double NDF = 0;
      double weight = 0;
      for (REvent::SignalStationIterator sIt = fgCurrentREvent->SignalStationsBegin();
           sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
        ++fit.NDF;
        const double sigma_time_snr_ = 20.5 * pow(sIt->GetRecData().GetParameter(revt::eSignalToNoise), -1.03/1.1);
        const double sTimeSigma2 = sigma_time_snr_ * sigma_time_snr_;
        weight += 2 * sqrt(1 + 1/sTimeSigma2 * 5/4) - 2;
      }

      weight /= fit.NDF;

      nPar = 2;
      if (fit.NDF > 4)
        nPar = 4;

      // fit the data points
      TMinuit minuit(nPar);
      int errFlag = 0;
      double argList[10];
      argList[0] = -1;
      minuit.SetPrintLevel(-1);

      minuit.SetFCN(RdHyperbolicWavefrontFit::RdHyperbolicWavefrontFnc);
      argList[0] = 1; // chi2
      minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
      minuit.mnparm(0, "gamma", fit.gamma, 0.001, 0, 0, errFlag);
      minuit.mnparm(1, "t0", fit.t0, 0.001, 0, 0, errFlag);
      if (nPar == 4) {
        minuit.mnparm(2, "u", fit.u, 0.01, 0, 0, errFlag);
        minuit.mnparm(3, "v", fit.v, 0.01, 0, 0, errFlag);
      }
      // init fnc
      argList[0] = 1;
      minuit.mnexcm("CALI", argList, 1, errFlag);

      argList[0] = 500;
      minuit.mnexcm("MINIMIZE", argList, 1, errFlag);
      if (errFlag) {
        INFOIntermediate("  minuit minimize error: " + errFlag);
        minuit.mnexcm("MINOS", argList, 0, errFlag);
      if (errFlag) {
        INFOIntermediate("  minuit minimize error: " + errFlag);
          break;
        }
      }

      minuit.GetParameter(0, fit.gamma, fit.gammaError);
      minuit.GetParameter(1, fit.t0, fit.t0Error);
      if (nPar == 4) {
        minuit.GetParameter(2, fit.u, fit.uError);
        minuit.GetParameter(3, fit.v, fit.vError);
      }

      double edm, errdef;
      Int_t nvpar, nparx;
      minuit.mnstat(fit.minChi, edm, errdef, nvpar, nparx, fit.status);

      //fill the vector with the intermediate results
      const std::vector<double> vec({
        fit.gamma, //gamma
        fit.gammaError, // gamma error
        fit.t0, // offset
        fit.minChi/weight, // chi2
        double(fit.NDF - nvpar), // NDF
        gCore.GetX(FixedfgLocalCS), // x
        gCore.GetY(FixedfgLocalCS), // y
        fit.minChi / double(fit.NDF - nvpar) / weight, // decision parameter
        fit.u, // u
        fit.v // v
      });

      result.push_back(vec);

      gCore = core;
      fgLocalCS = fwk::LocalCoordinateSystem::Create(gCore);
    }

    //find the best fit
    sort(result.begin(), result.end(), Sorter);

    Point finalCore;
    if (result.size() > 0) {
      const auto& r0 = result.front();
      fit.gamma = r0[0];
      fit.gammaError = r0[1];
      fit.minChi = r0[3];
      fit.NDF = r0[4];
      const Point core(r0[5], r0[6], 0, FixedfgLocalCS);
      finalCore = core;
      fit.u = r0[8];
      fit.v = r0[9];
      fit.t0 = r0[2];
    }

    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    // store the fit results
    if (!event.HasRecShower()) {
      event.MakeRecShower();
    }
    if (!event.GetRecShower().HasRRecShower()) {
      event.GetRecShower().MakeRRecShower();
    }
    ShowerRRecData& showerReco = event.GetRecShower().GetRRecShower();

    const CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(finalCore);
    utl::Vector axis(1, fit.u, fit.v, localCS, Vector::kSpherical);
    axis = -axis;

    //calculates residuals and fill them in station rec
    CalculateTimeResidual(axis, fit.gamma, fit.t0);

    showerReco.SetParameter(eFitSuccess, 1);
    showerReco.SetParameter(eWaveFitNDF, fit.NDF);
    showerReco.SetParameter(eWaveFitChi2, fit.minChi);
    showerReco.SetParameter(eWaveFrontModel, ShowerRRecData::kHyperbolicWaveFit);

    //#warning TODO: Set covariances!!!

    showerReco.SetParameter(eCoreX, finalCore.GetX(referenceCS));
    showerReco.SetParameter(eCoreY, finalCore.GetY(referenceCS));
    showerReco.SetParameter(eCoreZ, finalCore.GetZ(referenceCS));
    showerReco.SetParameter(eShowerAxisX, axis.GetX(referenceCS));
    showerReco.SetParameter(eShowerAxisY, axis.GetY(referenceCS));
    showerReco.SetParameter(eShowerAxisZ, axis.GetZ(referenceCS));
    showerReco.SetParameter(eGammaHyperbolicParam, fit.gamma);
    showerReco.SetParameterError(eGammaHyperbolicParam, fit.gammaError);
    showerReco.SetParameter(eHyperbolicWfFitChi2, fit.minChi);
    showerReco.SetParameter(eHyperbolicWfFitNDF, fit.NDF);
    showerReco.SetParameter(eRecStage, ShowerRRecData::kHyperbolicWaveFit);

    return eSuccess;
  }


  VModule::ResultFlag
  RdHyperbolicWavefrontFit::Finish()
  {
    return eSuccess;
  }


  void
  RdHyperbolicWavefrontFit::RdHyperbolicWavefrontFnc(int& nPar, double* const /*grad*/, double& value, double* const par,
                                                     const int /*flag*/)
  {
    //fit parameters
    const double& gamma = par[0];
    const double& t0 = par[1];
    double theta = gSAxis.GetTheta(fgLocalCS);
    double phi = gSAxis.GetPhi(fgLocalCS);
    if (nPar == 4) {
      theta = par[2];
      phi = par[3];
    }

    const utl::Vector axis(1, theta, phi, fgLocalCS, Vector::kSpherical);

    //M-estimator
    double Mscore = 0;

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    for (REvent::SignalStationIterator sIt = fgCurrentREvent->SignalStationsBegin();
         sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {  //loop over signal stations

      //const StationRecData& sRec = sIt->GetRecData();
      const double sigma_time_snr_ = 20.5 * pow(sIt->GetRecData().GetParameter(revt::eSignalToNoise), -1.03/1.1);
      const double sTimeSigma2 = sigma_time_snr_ * sigma_time_snr_;

      // calculate Rd station distance to shower core
      const Point& stationPosition = rDetector.GetStation(*sIt).GetPosition();
      const Vector VecDistanceToCore = stationPosition - rDetector.GetStation(gStation0Id).GetPosition();
      // project the vector on shower direction
      const double projVec = axis * VecDistanceToCore;
      // calculate pulse arrival time
      const double relativePulseTime = projVec / (kSpeedOfLight / nanosecond);
      //double planeTime = relativePulseTime;
      const double signalTime = sIt->GetRecData().GetParameter(revt::eSignalTime);
      const double projTime = signalTime - gTime0 - relativePulseTime;

      const Vector a = stationPosition - gCore;
      const Vector b = (axis*a)*axis;
      const Vector perp = a - b;

      const double expectedTime = sqrt(1 +  perp.GetMag() *  perp.GetMag() / (gamma * gamma) ) * 9 - 9 + t0;
      const double A = 1;//sIt->GetRecData().GetParameter(revt::eSignal)/1.e-6;

      const double dt = projTime - expectedTime;
      Mscore += 2 * sqrt(1 + A * dt*dt / (2 * sTimeSigma2)) - 2;

    }

    value = Mscore;
  }

  void
  RdHyperbolicWavefrontFit::CalculateTimeResidual(const utl::Vector& axis, const double gamma, const double t0) const
  {
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    for(REvent::SignalStationIterator sIt = fgCurrentREvent->SignalStationsBegin(); sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) { //loop over stations
        StationRecData& sRec = sIt->GetRecData();
        const double signal_time = sRec.GetParameter(eSignalTime);
        const double sigma_signal_time = sRec.GetParameterError(eSignalTime);

        const Point& stationPosition = rDetector.GetStation(*sIt).GetPosition();
        const Vector a = stationPosition - gCore;
        const Vector b = (axis*a)*axis;
        const Vector perp = a - b;
        const double expectedTime = sqrt(1 +  perp.GetMag() *  perp.GetMag() / (gamma * gamma) ) * 9 - 9 + t0;

        sRec.SetParameter(eTimeResidual, signal_time - expectedTime);
        sRec.SetParameterError(eTimeResidual, sigma_signal_time);
    }
  }

}