LDFFitFunction.cc

Go to the documentation of this file.

#include "LDFFitFunction.h"

#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AugerUnits.h>
#include <utl/Integrator.h>
#include <utl/AxialVector.h>

#include <fwk/LocalCoordinateSystem.h>

#include <iostream>
#include <cmath>
#include <vector>

#include <TMath.h>

using namespace utl;


namespace RdHASLDFFitter {

  double
  LDFFitFunction::GeomagneticEmissionFactor(const double ceFraction,
                                            const double sineAlpha,
                                            const double cosAzimuth)
  {
    return std::pow(1 + cosAzimuth * std::sqrt(ceFraction) / std::abs(sineAlpha), -2);
  }


  double
  LDFFitFunction::GeomagneticEmissionFactor(const double ceFraction,
                                            const double sineAlpha,
                                            const utl::Vector& showerAxis,
                                            const utl::Vector& magneticFieldVector,
                                            const utl::Point& showerCore,
                                            const utl::Point& stationPosition)
  {
    const utl::Vector vxB = Cross(-showerAxis, magneticFieldVector);
    const utl::Vector coreStation = stationPosition - showerCore;
    const utl::Vector stationInShowerPlane = coreStation - (showerAxis * coreStation) * showerAxis;
    const double cosAzimuth = CosAngle(vxB, stationInShowerPlane);

    return GeomagneticEmissionFactor(ceFraction, sineAlpha, cosAzimuth);
  }


  double
  LDFFitFunction::ChargeExcessFractionParamICRC19(const double axisDistance,
                                                  const double distanceToXmax,
                                                  const double densityXmax)
  {
    /* Returns charee excess fraction for each observer defined as a = sin(alpha) ** 2 * f_ce / f_geo
       See pos(ICRC2019)294 for more information */
    const double a = 0.37313183;
    const double b = 1.31124484;  // / meter
    const double p0 = 0.1889835;
    const double p1 = 6.71403002;  //  / (kg / m3);
    const double densityAvg = 0.4;  // * (kg / m3);

    return a * (axisDistance / distanceToXmax) *
      std::exp(b * axisDistance / 1000) * (std::exp(p1 * (densityXmax - densityAvg)) - p0);
  }


  double
  LDFFitFunction::ChargeExcessFractionParamICRC21(const double axisDistance,
                                                  const double distanceToXmax,
                                                  const double densityXmax)
  {
    // Returns charee excess fraction for each observer defined as a = sin(alpha) ** 2 * f_ce / f_geo
    const double a1 = -1.17523609e-06;
    const double a2 = 0.348154734;
    const double b = 1.6068519502678418;  // / meter
    const double p0 = 1.66965694e+01;
    const double p1 = 3.31954904e+00;
    const double p2 = -5.73577715e-03;

    const double a = a1 * distanceToXmax + a2;
    return a * (axisDistance / distanceToXmax) *
      std::exp(b * axisDistance / 1000) * (p0 * std::pow(densityXmax, p1) + p2);
  }


  double
  LDFFitFunction::GetChargeExcessFraction(const double axisDistance)
    const
  {
    /* Charge-excess fraction from parametrization.
       fDistanceToXmax has to be in meters
       fDensityAtXmax has to be in kg/m3 (which is taken care of in UpdateParameterParam) */

    double ceFraction = ChargeExcessFractionParamICRC21(axisDistance, fDistanceToXmax, fDensityAtXmax);

    /* For to low densityXmax resp. to large distanceToXmax / zenith angle ceFraction can become negative.
       Set ceFraction = 0 as resonable approximation.
       This happen for ChargeExcessFractionParamICRC19 with densityXmax <~0.15 kg/m3.
                 This can exclude the shower maxima of showers with a zenith angle ~> 82 deg.
       This happen for ChargeExcessFractionParamICRC21 with densityXmax <~0.09 kg/m3.
                 This does not exclude any more shower maxima of showers < 85 deg.
                 (With an imperfect fitted fDistanceToXmax this can still a < 0 can still occur) */

    if (ceFraction < 0)
      ceFraction = 0;

    return ceFraction;
  }


  double
  LDFFitFunction::FGeomagneticParam(const double xvBvvB, const double yvBvvB, const double efvxB, const double cEarlyLate)
    const
  {
    // xvBvvB, yvBvvB have to be in fitted core system

    // calculate axis distance
    const double offAxisDistance = std::sqrt(Sqr(xvBvvB) + Sqr(yvBvvB)) / cEarlyLate;

    const double polarAngle = atan2(yvBvvB, xvBvvB);

    // get charge excess fraction from parametrization. set 0 beyond 80 deg and if get negative
    const double ceFraction = GetChargeExcessFraction(offAxisDistance);

    // factor to get geomagnetic emission from f_vxB (f_geo = f_vxB * c_geo)
    const double cGeo = GeomagneticEmissionFactor(ceFraction, fShowerData.fSineAlpha, cos(polarAngle));

    // get symmetrized geomagnetic emission
    const double efGeoMeas = efvxB * cGeo * Sqr(cEarlyLate);

    return efGeoMeas;
  }


  double
  LDFFitFunction::FVxBPrediction(const double xvBvvB, const double yvBvvB, const double cEarlyLate)
    const
  {
    // xvBvvB, yvBvvB have to be in fitted core system

    // calculate axis distance
    const double offAxisDistance = std::sqrt(Sqr(xvBvvB) + Sqr(yvBvvB)) / cEarlyLate;

    const double polarAngle = std::atan2(yvBvvB, xvBvvB);

    // get charge excess fraction from parametrization. set 0 beyond 80 deg and if get negative
    const double ceFraction = GetChargeExcessFraction(offAxisDistance);

    // get symmetric signal prediction for the geomagnetic emission (at corrected position)
    //const double efGeoModel = FABCD(offAxisDistance, fr0, fA, fB, fC, fD);
    const double efGeoModel = FEgeo(offAxisDistance);

    // factor to get geomagnetic emission from f_vxB (f_geo = f_vxB * c_geo)
    const double cGeo = GeomagneticEmissionFactor(ceFraction, fShowerData.fSineAlpha, std::cos(polarAngle));

    // add charge-excess to predicted geomagnetic model + 'revert' early late correction
    const double efvxBPredict = efGeoModel / (cGeo * Sqr(cEarlyLate));

    return efvxBPredict;
  }


  double
  LDFFitFunction::FGeomagneticPos(const double xvBvvB, const double yvBvvB, const double cEarlyLate, const StationFitData& station)
    const
  {
    // xvBvvB, yvBvvB have to be in fitted core system
    const double polarAngle = std::atan2(yvBvvB, xvBvvB);

    // get symmetrized geomagnetic emission
    const double efGeoPos =
      Sqr((std::sqrt(station.fEnergyFluenceVxB) - (cos(polarAngle) / std::abs(sin(polarAngle))) * std::sqrt(station.fEnergyFluenceVxVxB)) * cEarlyLate);

    // if station is to close at or on the vxB axis reject its from minimization
    if (std::abs(std::sin(polarAngle)) < 0.1)
      station.fRejected = true;
    else
      station.fRejected = false;

    return efGeoPos;
  }


  double
  LDFFitFunction::FGeomagneticPosErr(const double xvBvvB, const double yvBvvB, const double efGeoPos, const StationFitData station)
    const
  {
    // xvBvvB, yvBvvB have to be in fitted core system
    const double polarAngle = std::atan2(yvBvvB, xvBvvB);

    // get symmetrized geomagnetic emission
    const double efGeoPosErr = sqrt(efGeoPos * (std::pow(station.fEnergyFluenceVxBErr, 2) / station.fEnergyFluenceVxB +
      std::pow((std::cos(polarAngle) * station.fEnergyFluenceVxVxBErr) / (std::abs(std::sin(polarAngle)) * std::sqrt(station.fEnergyFluenceVxVxB)), 2)));
    // std::cout << efGeoPosErr / efGeoPos << ' ' << (cos(polarAngle) / abs(sin(polarAngle))) << ' ' << station.fRejected << ' ' << abs(sin(polarAngle)) << std::endl;
    return efGeoPosErr;
  }


  RadioGeometryUtilities
  LDFFitFunction::GetFittedRadioCoreTransformation()
    const
  {
    // this "includes" the core fit
    const utl::Point radioCore = fTransformationRefCore.GetVectorFromShowerPlaneVxB(fCoreX, fCoreY);
    return RadioGeometryUtilities(fShowerData.fRefShowerAxis, fwk::LocalCoordinateSystem::Create(radioCore),
                                  fShowerData.fMagneticField);
  }


  double
  LDFFitFunction::operator()(const std::vector<double>& pars)
    const
  {
    // update parameter according to parametrization
    UpdateParameterParam(pars);

    // this "includes" the core fit
    utl::RadioGeometryUtilities transformationRadioCore = GetFittedRadioCoreTransformation();

    double loss = 0;
    for (const auto& s : fStationData) {

      double xvxB = 0;
      double yvxB = 0;
      double zvxB = 0;
      transformationRadioCore.GetVectorInShowerPlaneVxB(xvxB, yvxB, zvxB, s.fPosition);

      const double cEarlyLate =
        utl::RadioGeometryUtilities::GetEarlyLateCorrectionFactor(fDistanceToXmax, zvxB);
      s.fCEarlyLate = cEarlyLate;

      const double offAxisDistance = std::sqrt(Sqr(xvxB) + Sqr(yvxB)) / cEarlyLate;

      if (fFitConfig.fUseParam) {
        const double efvxBPredict = FVxBPrediction(xvxB, yvxB, cEarlyLate);
        const double efGeoParam = FGeomagneticParam(xvxB, yvxB, s.fEnergyFluenceVxB, cEarlyLate);
        const double ceFraction = GetChargeExcessFraction(offAxisDistance);

        s.fEnergyFluenceVxBPredict = efvxBPredict;
        s.fCeFraction = ceFraction;
        s.fEnergyFluenceGeomagnetic = efGeoParam;

        // Not Correct: Not considering uncertainty on fitted distance to Xmax, core, arrival direction
        // (Only stored for plotting, not used in fit)
        s.fEnergyFluenceGeomagneticErr = FGeomagneticParam(xvxB, yvxB, s.fEnergyFluenceVxBErr, cEarlyLate);

        // Calculate any quantity also for rejected stations, however ignore them in the minimization.
        if (s.fRejected)
          continue;

        loss += GetLoss(efvxBPredict, s.fEnergyFluenceVxB, s.fEnergyFluenceVxBErr);
      } else {
        const double efGeoPos = FGeomagneticPos(xvxB, yvxB, cEarlyLate, s);
        const double efGeoPosErr = FGeomagneticPosErr(xvxB, yvxB, efGeoPos, s);
        const double offAxisDistance = std::sqrt(Sqr(xvxB) + Sqr(yvxB)) / cEarlyLate;
        const double efGeoModel = FEgeo(offAxisDistance);

        s.fEnergyFluenceGeomagnetic = efGeoPos;
        s.fEnergyFluenceGeomagneticErr = efGeoPosErr;

        // Calculate any quantity also for rejected stations, however ignore them in the
        // minimization. FGeomagneticPos is rejected if abs(sin(polarAngle)) < 0.1
        if (s.fRejected)
          continue;

        loss += GetLoss(efGeoModel, efGeoPos, efGeoPosErr);
      }

    }

    return std::isnan(loss) ? std::numeric_limits<double>::max() : loss;
  }


  void
  LDFFitFunction::GetPrediction(const std::vector<double>& pars)
  {
    this->operator()(pars);
  }


  double
  LDFFitFunction::GetLossLog(const double model,
                             const double data,
                             const double uncertainty)
    const
  {
    return Sqr((std::log(model) - std::log(data)) / (uncertainty / data));  // check if that makes sence
  }


  double
  LDFFitFunction::GetLoss(const double model,
                          const double data,
                          const double uncertainty)
    const
  {
    return Sqr((model - data) / uncertainty);
  }


  double
  LDFFitFunction::GetChi2(const std::vector<double>& pars)
  {
    return this->operator()(pars);
  }


  int
  LDFFitFunction::GetNDF()
  {
    double ndf = 0;

    // loop over all stations
    for (const auto& s : fStationData) {
      if (s.fRejected)
        continue;
      ++ndf;
    }

    return ndf - GetNFreeParameters();
  }


  /*void
  LDFFitFunction::TestStations(const std::vector<double>& pars, const double chi2_mean)
  {
   double chi2_station = 0;
   for (const auto& s : fStationData) {
     chi2_station = GetChi2HASLDFModel1(s.off_axis_distance, pars, s.f_geo, s.f_geo_err);
     if (chi2_station > 10 * chi2_mean) {
       std::cout << "Station " << s.id << " doesnt fit: " << chi2_station << " > 10 * "
                 << chi2_mean << " (f= " << s.f_geo << " +- " << s.f_geo_err << ")" << std::endl;
     }
   }
   // can reject stations
  }*/


  // ----------------------- LDFFitFunctionPoly3 -------------------------//
  double
  LDFFitFunctionPoly3::GetNormalization()
    const
  {
    const double B = fB;
    const double C = fC;
    const double D = fD;
    const double r0 = fr0;

    auto LDFnorm = [&B, &C, &D, &r0](const double r)
    { return kTwoPi * r * FABCD(r, r0, 1, B, C, D); };

    return MakeIntegrator(LDFnorm).GetIntegral(0, fRnorm) * fEgeoCorr;
  }


  double
  LDFFitFunctionPoly3::FEgeo(const double r)
    const
  {
    return fEgeo * FABCD(r, fr0, 1, fB, fC, fD) / fNorm;
  }


  double
  LDFFitFunctionPoly3::FABCD(const double r, const double r0, const double a,
                             const double b, const double c, const double d)
  {
    const double x = r / r0;
    return a * std::exp(x * (-b + x * (std::exp(c) - x * std::exp(d))));
  }


  // ------------------- LDFFitFunctionGaussSigmoid ----------------------//
  double
  LDFFitFunctionGaussSigmoid::CalculateR0()
    const
  {
    const double refracIndexAtXmax = fShowerData.fRefracFromHeight.Y(fHeight);
    const double cherenkovAngle = std::acos(1 / refracIndexAtXmax);
    return std::tan(cherenkovAngle) * fDistanceToXmax;
  }


  double
  LDFFitFunctionGaussSigmoid::GetNormalization()
    const
  {
    const double r0 = fR0;
    const double sig = fSig;
    const double p = fP;
    const double r02 = fR02;
    const double arel = fArel;
    const double slope = fSlope;

    if (fFitConfig.fUseSoftExponent) {
      auto LDFnorm = [&r0, &sig, &p, &arel, &slope, &r02](const double r)
      { return kTwoPi * r * FGaussSigmoidSoft(r, r0, sig, p, arel, slope, r02); };

      return MakeIntegrator(LDFnorm).GetIntegral(0, fRnorm) * fEgeoCorr;
    } else {
      auto LDFnorm = [&r0, &sig, &p, &arel, &slope, &r02](const double r)
      { return kTwoPi * r * FGaussSigmoidHard(r, r0, sig, p, arel, slope, r02); };

      return MakeIntegrator(LDFnorm).GetIntegral(0, fRnorm) * fEgeoCorr;
    }
  }


  double
  LDFFitFunctionGaussSigmoid::FEgeo(const double r)
    const
  {
    return fEgeo * FGaussSigmoid(r, fR0, fSig, fP, fArel, fSlope, fR02) / fNorm;
  }


  double
  LDFFitFunctionGaussSigmoid::FGaussSigmoidSoft(
    const double r, const double r0, const double sig, double p,
    const double arel, const double slope, const double r02)
  {
    const double pEff = (r > r0) ? 2 * std::pow(r0 / r, p / 1000) : 2;
    return std::exp(-std::pow((r - r0) / sig, pEff)) + arel / (1 + std::exp(slope * (r / r0 - r02)));
  }

  double
  LDFFitFunctionGaussSigmoid::FGaussSigmoidHard(
    const double r, const double r0, const double sig, double p,
    const double arel, const double slope, const double r02)
  {
    const double pEff = (r > r0) ? p : 2;
    return std::exp(-std::pow((r - r0) / sig, pEff)) + arel / (1 + std::exp(slope * (r / r0 - r02)));
  }

}