RdSphericalFit.cc

Go to the documentation of this file.

#include "RdSphericalFit.h"
#include "SphericalFitFunction.h"

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

// utilities
#include <utl/ErrorLogger.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Minou.h>
#include <utl/AnalyticCoordinateTransformator.h>

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

// revent
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>
#include <revt/StationRRecDataQuantities.h>

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

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

#include <TMath.h>
#include <TMatrixD.h>

#include <string>
#include <vector>


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


namespace RdSphericalFit {

  vector<double>
  GetTimeResiduals(const vector<StationFitData>& vStationFitData,
                   const Vector& showerAxis, const double tau0)
  {
    vector<double> expectedTimes;
    vector<double> timeResiduals;

    for (const auto& d : vStationFitData) {
      const double expTime = (d.antennaPos - showerAxis).GetMag() / kSpeedOfLight;
      expectedTimes.push_back(expTime);
    }

    const double t0 = TMath::Mean(expectedTimes.begin(), expectedTimes.end());

    for (unsigned int i = 0, n = vStationFitData.size(); i < n; ++i) {
      const double residual = (vStationFitData[i].signalTime - tau0) - (expectedTimes[i] - t0);
      timeResiduals.push_back(residual);
    }

    return timeResiduals;
  }


  VModule::ResultFlag
  RdSphericalFit::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdSphericalFit");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("minNumberOfStations").GetData(fMinNumberOfStations);
    topBranch.GetChild("useInitialGeometry").GetData(fUseInitialGeometry);
    topBranch.GetChild("thetaVarC").GetData(fThetaVarC);
    topBranch.GetChild("thetaMin").GetData(fThetaMin);
    topBranch.GetChild("thetaMax").GetData(fThetaMax);
    topBranch.GetChild("rMax").GetData(fRadiusMax);

    const string yesnostr = topBranch.GetChild("allowMultipleCalls").Get<string>();
    fLockParameterStorage = (yesnostr == "no");

    return eSuccess;
  }


  VModule::ResultFlag
  RdSphericalFit::Run(Event& event)
  {
    // nothing to do if there is no REvent
    if (!event.HasREvent()) {
      INFOFinal("eContinueLoop: No REvent, so no reconstruction");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();
    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();
    const int fNumberOfStations = rEvent.GetNumberOfSignalStations();

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

    ostringstream info;

    if (fNumberOfStations < fMinNumberOfStations) {
      info.str("");
      info << "Event contains only " << fNumberOfStations
           << " station(s), not enough for spherical fit. Event is skipped";
      INFOFinal(info);
      return eContinueLoop;
    }

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

    // get coordinate origin depending on xml settings
    Vector initialShowerAxis;
    Point initialCore;

    if (fUseInitialGeometry == "SD") {
      if (!event.GetRecShower().HasSRecShower()) {
        WARNING("Event has no SRecShower. Event is skipped");
        return eContinueLoop;
      }
      initialShowerAxis = event.GetRecShower().GetSRecShower().GetAxis();
      initialCore = event.GetRecShower().GetSRecShower().GetCorePosition();
    } else if (fUseInitialGeometry == "MC") {
      if (!event.HasSimShower()) {
        WARNING("Event has no simulated shower. Event is skipped");
        return eContinueLoop;
      }
      initialShowerAxis = -event.GetSimShower().GetDirection();
      initialCore = event.GetSimShower().GetPosition();
    } else if (fUseInitialGeometry == "RD") {
      if (!event.GetRecShower().HasRRecShower()) {
        WARNING("Event has no HasRRecShower. Event is skipped");
        return eContinueLoop;
      }
      initialShowerAxis = showerrrec.GetAxis();
      initialCore = showerrrec.GetCorePosition();
    } else if (fUseInitialGeometry == "Reference") {
      if (!event.GetRecShower().HasRRecShower()) {
        WARNING("Event has no HasRRecShower. Event is skipped");
        return eContinueLoop;
      }
      initialShowerAxis = event.GetRecShower().GetRRecShower().GetReferenceAxis(event);
      initialCore = event.GetRecShower().GetRRecShower().GetReferenceCorePosition(event);
    }

    const CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(initialCore);
    const double initialZenith = initialShowerAxis.GetTheta(localCS);
    const double initialAzimuth = initialShowerAxis.GetPhi(localCS);

    // initiate atm + tables (which model is used is decieded in the config)
    const double avgShowerXmax = 750*g/cm2;
    const atm::Atmosphere& atmo = detector.GetAtmosphere();
    atmo.InitSlantProfileModel(initialCore, initialShowerAxis, 10*g/cm2);
    const atm::ProfileResult& distanceFromDepth = atmo.EvaluateDistanceVsSlantDepth();
    // jep the result is negative, convention thing...
    const double initialRadius = abs(distanceFromDepth.Y(avgShowerXmax));

    info.str("");
    info << "Event has " << fNumberOfStations << " antennas, "
         << "initial geometry is " << fUseInitialGeometry << "\n"
         << "Zenith = " << initialZenith / deg << ", azimuth = "
         << initialAzimuth / deg << ", radius = " << initialRadius / km;
    INFOIntermediate(info);

    // absolute event time for both the test and measured times
    // choosen as mean of all measured/test times
    double tau0 = 0;
    vector<StationFitData> vStationFitData;

    for (const auto& station : rEvent.SignalStationsRange()) {
      const Vector sPos = rDetector.GetStation(station).GetPosition() - initialCore;
      const StationRecData& sRec = station.GetRecData();
      const double t = sRec.GetParameter(eSignalTime);
      const double t_err = sRec.GetParameterError(eSignalTime);

      StationFitData stationFitData;
      stationFitData.signalTime = t;
      stationFitData.signalTimeErr = t_err;
      stationFitData.antennaPos = sPos;
      vStationFitData.push_back(stationFitData);

      tau0 += t;
    }
    tau0 /= vStationFitData.size();

    // Initalize parameters (WARNING: Take care of correct ordering)
    // name, value, step, min, max, fixed
    // min = max = 0 -> no limits
    vector<Minou::ParameterDef> parameters;
    parameters.emplace_back("radius", initialRadius, 100, 0, fRadiusMax, false);
    parameters.emplace_back("zenith", initialZenith, 0.1, fThetaMin, fThetaMax, false);
    parameters.emplace_back("azimuth", initialAzimuth, 0.01, 0, 0, false);
    parameters.emplace_back("gamma", 1, 1e-3, 0, 0, true);

    // do the fit
    SphericalFitFunction fitFunction(vStationFitData, parameters, localCS, tau0);
    Minou::Minimizer<SphericalFitFunction> min(fitFunction);
    if (fInfoLevel >= eInfoIntermediate)
      min.SetVerbose(true);

    if (min.Minimize(500)) {
      INFOFinal("Fit falied, skipping event");
      showerrrec.SetParameter(eFitSuccess, 0, fLockParameterStorage);
      return eContinueLoop;
    }

    // Redo fit if zenith above threshold allowing variation of speed of light, cf. GAP2011-92
    // (Reconstruction of the Wave Front of Radio Signals at AERA)
    // unphysical situations can arise when the signal is coming from directions close to the horizon
    // not implemented yet.
    /*
    if (min.GetParameter(1) > fThetaVarC)
      parameters.emplace_back("gamma", 1, 1e-4, 0.8, 1.2, false);
    */

    // get fit results
    const auto tmpFitResult = min.GetParameters();
    const auto bestParameterError = min.GetParametersAndErrors();
    const TMatrixD pcov = min.GetCovarianceMatrix();

    const double radius = bestParameterError.at(0).first;
    const double theta = bestParameterError.at(1).first;
    const double phi = bestParameterError.at(2).first;
    const double gamma = bestParameterError.at(3).first;

    const double radiusErr = bestParameterError.at(0).second;
    const double thetaErr = bestParameterError.at(1).second;
    const double phiErr = bestParameterError.at(2).second;
    const double gammaErr = bestParameterError.at(3).second;

    // the spherical wave fit can be performed with just 4 stations. The number of degrees
    // of freedom is 0, though. To be able to compute chiSquare/NDF the NDF is set to 1.
    const int ndf = max(fitFunction.GetNDF(), 1);
    const double chi2 = fitFunction.GetChi2(tmpFitResult);

    info.str("");
    info << "Spherical wavefront fit succesfull \n"
            "\ttheta = " << theta / deg  << " +/- " << thetaErr << "\n"
            "\tphi = " << phi / deg  << " +/- " << phiErr / deg << "\n"
            "\tR = " << radius / km  << " +/- " << radiusErr / km << "\n"
            "\tchi2 / ndf = " << chi2 << " / " << ndf << " = " << chi2/ndf;
    INFOFinal(info);

    // save results
    const Vector axis(radius, theta, phi, localCS, Vector::kSpherical);
    showerrrec.SetParameter(eShowerAxisX, axis.GetX(localCS), fLockParameterStorage);
    showerrrec.SetParameter(eShowerAxisY, axis.GetY(localCS), fLockParameterStorage);
    showerrrec.SetParameter(eShowerAxisZ, axis.GetZ(localCS), fLockParameterStorage);

    showerrrec.SetParameter(eGamma, gamma, fLockParameterStorage);
    showerrrec.SetParameterError(eGamma, gammaErr, fLockParameterStorage);

    showerrrec.SetParameter(eFitSuccess, 1, fLockParameterStorage);
    showerrrec.SetParameter(eWaveFitChi2, chi2, fLockParameterStorage);
    showerrrec.SetParameter(eWaveFitNDF, ndf, fLockParameterStorage);

    double tempCov =
      AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 0, radius, theta, phi,
        pcov(0,0), pcov(1,1), pcov(2,2), pcov(0,1), pcov(0,2), pcov(1,2));
    showerrrec.SetParameterCovariance(eShowerAxisX, eShowerAxisX, tempCov, fLockParameterStorage);
    tempCov =
      AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(1, 1, radius, theta, phi,
        pcov(0,0), pcov(1,1), pcov(2,2), pcov(0,1), pcov(0,2), pcov(1,2));
    showerrrec.SetParameterCovariance(eShowerAxisY, eShowerAxisY, tempCov, fLockParameterStorage);
    tempCov =
      AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(2, 2, radius, theta, phi,
        pcov(0,0), pcov(1,1), pcov(2,2), pcov(0,1), pcov(0,2), pcov(1,2));
    showerrrec.SetParameterCovariance(eShowerAxisZ, eShowerAxisZ, tempCov, fLockParameterStorage);
    tempCov =
      AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 1, radius, theta, phi,
        pcov(0,0), pcov(1,1), pcov(2,2), pcov(0,1), pcov(0,2), pcov(1,2));
    showerrrec.SetParameterCovariance(eShowerAxisX, eShowerAxisY, tempCov, fLockParameterStorage);
    tempCov =
      AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 2, radius, theta, phi,
        pcov(0,0), pcov(1,1), pcov(2,2), pcov(0,1), pcov(0,2), pcov(1,2));
    showerrrec.SetParameterCovariance(eShowerAxisX, eShowerAxisZ, tempCov, fLockParameterStorage);
    tempCov =
      AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(1, 2, radius, theta, phi,
        pcov(0,0), pcov(1,1), pcov(2,2), pcov(0,1), pcov(0,2), pcov(1,2));
    showerrrec.SetParameterCovariance(eShowerAxisY, eShowerAxisZ, tempCov, fLockParameterStorage);

    // set time residuals and arrival direction of signal stations
    const vector<double> timeResiduals = GetTimeResiduals(vStationFitData, axis, tau0);

    int iStation = 0;
    for (auto& station : rEvent.SignalStationsRange()) {
      StationRecData& sRec = station.GetRecData();

      sRec.SetParameter(eTimeResidual, timeResiduals[iStation], fLockParameterStorage);
      sRec.SetParameterError(eTimeResidual, sRec.GetParameterError(eSignalTime), fLockParameterStorage);

      // for spherical wave use direction to point source aka axis
      // therefore, emission point is transfered to local coordniate system of station
      const CoordinateSystemPtr stationCS = rDetector.GetStation(station).GetLocalCoordinateSystem();
      sRec.SetParameter(eSignalArrivalDirectionX, axis.GetX(stationCS), fLockParameterStorage);
      sRec.SetParameter(eSignalArrivalDirectionY, axis.GetY(stationCS), fLockParameterStorage);
      sRec.SetParameter(eSignalArrivalDirectionZ, axis.GetZ(stationCS), fLockParameterStorage);

      ++iStation;
    }

    return eSuccess;
  }


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

}