LDFFitFunction.h

Go to the documentation of this file.

#ifndef _RdHASLDFFitter_LDFFitFunction_h_
#define _RdHASLDFFitter_LDFFitFunction_h_

#include "RdHASLDFFitter.h"

#include <det/Detector.h>

#include <utl/Minou.h>
#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/RadioGeometryUtilities.h>

#include <fwk/LocalCoordinateSystem.h>

// Atmosphere
#include <atm/ProfileResult.h>
#include <atm/InclinedAtmosphericProfile.h>

#include <vector>

#include <TMath.h>

using namespace utl;


namespace RdHASLDFFitter {

  struct FitConfig {
    bool fUseStationsWithoutSignal = false;
    bool fFitLog = false;
    bool fUseParam = true;
    bool fUseSoftExponent = true;
  };


  struct StationFitData {
    int fId = 0;
    Point fPosition;  // = Point(0, 0, 0, det::Detector::GetInstance().GetReferenceCoordinateSystem());
    double fEnergyFluenceVxB = 0;
    double fEnergyFluenceVxBErr = 0;
    double fEnergyFluenceVxVxB = 0;
    double fEnergyFluenceVxVxBErr = 0;
    mutable double fEnergyFluenceVxBPredict = 0;
    mutable double fEnergyFluenceGeomagnetic = 0;
    mutable double fEnergyFluenceGeomagneticErr = 0;
    mutable double fCEarlyLate = 0;
    mutable double fCeFraction = 0;
    mutable bool fRejected = false;
  };


  struct ShowerFitData {
    double fSineAlpha = 0;
    double fZenith = 0;
    double fDistanceTo750 = 0;
    Vector fRefShowerAxis;
    Vector fMagneticField;
    Point fRefShowerCore;
    atm::ProfileResult fDensityFromHeight;
    atm::ProfileResult fHeightFromDistance;
    atm::ProfileResult fRefracFromHeight;
  };


  class LDFFitFunction : public Minou::Base {

  public:
    LDFFitFunction(const std::vector<StationFitData>& stationData,
                   ShowerFitData showerData,
                   const std::vector<Minou::ParameterDef>& pars,
                   const FitConfig fitconfig = FitConfig(),
                   const double egeoCorr = 1) :
      fStationData(stationData),
      fShowerData(showerData),
      fFitConfig(fitconfig),
      fTransformationRefCore(
      RadioGeometryUtilities(
        fShowerData.fRefShowerAxis,
        fwk::LocalCoordinateSystem::Create(fShowerData.fRefShowerCore),
      fShowerData.fMagneticField)),
      fEgeoCorr(egeoCorr)
    {
      GetParameterDefs() = pars;
    }

    double operator()(const std::vector<double>& pars) const;

    // return number of degree of freedom.
    int GetNDF();

    // returns final lost
    double GetChi2(const std::vector<double>& pars);

    // predicts (asymmetric) energy fluence in vxB
    double FVxBPrediction(const double xvBvvB, const double yvBvvB, const double cEarlyLate) const;

    // symmetrized signal: Predicts geomagnetic energy fluence from (measured) energy fluence in vxB
    double FGeomagneticParam(const double xvBvvB, const double yvBvvB,
                             const double efvxB, const double cEarlyLate) const;

    // Determine geomagnetic energy fluence from polarisation and pulse position in shower plane (wo any param.)
    double FGeomagneticPos(const double xvBvvB, const double yvBvvB, const double cEarlyLate,
                           const StationFitData& station) const;

      // Determine geomagnetic energy fluence error from polarisation and  pulse position in shower plane (wo any param.)
    double FGeomagneticPosErr(const double xvBvvB, const double yvBvvB, const double efGeoPos,
                              const StationFitData station) const;

    // Param of charge-excees fraction as function of the (early-late corredted) axis distrance
    // and the distance to Xmax (resp. density at Xmax)
    double GetChargeExcessFraction(const double axisDistance) const;

    // updates member variables with given set of parameter
    void GetPrediction(const std::vector<double>& pars);

    // implementation of the (param.) geomagnetic LDF for the different models
    virtual double FEgeo(const double r) const = 0;

    // implementation of the integral over the geomagnetic LDF (shape)
    virtual double GetNormalization() const = 0;

    // Updates parameterized parameter in each iteration of the fit
    virtual void UpdateParameterParam(const std::vector<double>& pars) const = 0;

    // return parametrized charge-excess fraction which is defined as a = sin(alpha) ** 2 * f_ce / f_geo
    static double ChargeExcessFractionParamICRC19(const double axisDistance,
                                                  const double distanceToXmax,
                                                  const double densityXmax);

    // return parametrized charge-excess fraction which is defined as a = sin(alpha) ** 2 * f_ce / f_geo
    static double ChargeExcessFractionParamICRC21(const double axisDistance,
                                                  const double distanceToXmax,
                                                  const double densityXmax);

    // return factor c_geo in: f_geo = f_vxB * c_geo
    static double GeomagneticEmissionFactor(const double ceFraction,
                                            const double sineAlpha,
                                            const double cosAzimuth);

    // return factor c_geo in: f_geo = f_vxB * c_geo
    static double GeomagneticEmissionFactor(const double ceFraction,
                                            const double sineAlpha,
                                            const utl::Vector& showerAxis,
                                            const utl::Vector& magneticFieldVector,
                                            const utl::Point& showerCore,
                                            const utl::Point& stationPosition);

  protected:
    // lost function (chi^2)
    double GetLoss(const double model, const double data, const double uncertainty) const;

    // lost function with log
    double GetLossLog(const double model, const double data, const double uncertainty) const;


    //void TestStations(const std::vector<double>& pars, const double chi2Mean);

    // Get transformation for updated core position
    RadioGeometryUtilities GetFittedRadioCoreTransformation() const;

    const std::vector<StationFitData>& fStationData;
    ShowerFitData fShowerData;
    const FitConfig fFitConfig;
    RadioGeometryUtilities fTransformationRefCore;

    const double fEgeoCorr = 1; // To correct an overall bias (used in GetNormalization())

    // generic fit parameter
    mutable double fEgeo = 0;
    mutable double fRnorm = 0;
    mutable double fNorm = 1;

    mutable double fCoreX = 0;
    mutable double fCoreY = 0;

    mutable double fDistanceToXmax = 0; // should be given in meter
    mutable double fDensityAtXmax = 0;  // not fitted but calculated from fDistanceToXmax
    mutable double fHeight = 0;  // not fitted but calculated from fDistanceToXmax

  };


  class LDFFitFunctionPoly3 : public LDFFitFunction {

  public:
    using LDFFitFunction::LDFFitFunction;

    // Returns integral over FABCD form 0 to fRnorm * 2pi
    double GetNormalization() const override;

    // LDF model (of the shape) of the geomagnetic emission
    static double FABCD(const double r, const double r0, const double a,
                        const double b, const double c, const double d);

    /* LDF model of the geomagnetic emission:
       FEgeo = E_geo * FABCD / GetNormalization() */
    double FEgeo(const double r) const override;

    // param of C': C = exp(C') (with C' beeing the fit parameter)
    static
    double
    GetC(const double dxmax)
    {
      const double dx = dxmax / 1000;
      return dx * (6.283649656834273e-05 * dx - 0.02574916964033989) + 20.865090319962185 / dx + 1.6373732045578846;  // conversion in km
    }

    // param of D': D = exp(D') (with D' beeing the fit parameter)
    static
    double
    GetD(const double dxmax)
    {
      const double dx = dxmax / 1000;
      return dx * (9.113744732999866e-05 * dx - 0.037475161541008134) + 33.83363364863988 / dx + 1.9430065360834818;  // conversion in km
    }

    // simple parameterisation of the cherenkov radius. Systematically lower than CalculateR0
    static
    double
    CherenkovRadius(const double densitymax, const double dxmax)
    {
      const double cherenkovAngleInDeg = 0.24905864 * std::log(densitymax) + 0.92165234;
      return std::tan(cherenkovAngleInDeg * utl::deg) * dxmax;
    }

  private:
    // specific fit parameter

    const double fB = 0;  // fixed
    mutable double fC = 0;
    mutable double fD = 0;
    const double fr0 = 800;  // fixed

    void
    UpdateParameterParam(const std::vector<double>& pars)
      const override
    {
      // Fit parameter
      fEgeo = pars[0];
      fDistanceToXmax = pars[1];
      fCoreX = pars[2];
      fCoreY = pars[3];

      // parametrized parameter
      fC = GetC(fDistanceToXmax);
      fD = GetD(fDistanceToXmax);

      // calculate density from distance
      // the -1 is needed because of convention of the table (opposit direction of the init vector)
      const double height = fShowerData.fHeightFromDistance.Y(-fDistanceToXmax);
      fDensityAtXmax = fShowerData.fDensityFromHeight.Y(height) / (kg / m3);

      fRnorm = 5 * CherenkovRadius(fDensityAtXmax, fDistanceToXmax);
      fNorm = GetNormalization();  // integral over FABCD from 0 to fRnorm * 2pi
    }

  };


  class LDFFitFunctionGaussSigmoid : public LDFFitFunction {

  public:
    using LDFFitFunction::LDFFitFunction;

    // Returns integral over FGaussSigmoid form 0 to fRnorm * 2pi
    double GetNormalization() const override;

    // Model of the geomagnetic LDF (shape) with soft exponent
    static double FGaussSigmoidSoft(const double r, const double r0, const double sig,
                                    const double p, const double arel, const double slope, const double r02);

    // Model of the geomagnetic LDF (shape) with hard exponent
    static double FGaussSigmoidHard(const double r, const double r0, const double sig,
                                    const double p, const double arel, const double slope, const double r02);

    // Model of the geomagnetic LDF (shape)
    double FGaussSigmoid(const double r, const double r0, const double sig,
                         const double p, const double arel, const double slope, const double r02) const
    { return (fFitConfig.fUseSoftExponent) ? FGaussSigmoidSoft(r, r0, sig, p, arel, slope, r02)
                                           : FGaussSigmoidHard(r, r0, sig, p, arel, slope, r02); }

    /* LDF model of the geomagnetic emission:
       FEgeo = E_geo * FGaussSigmoid / GetNormalization() */
    double FEgeo(const double r) const override;

    /* Returns the cherenkov radius calculated with simple model: Radius of a cone with its apex at Xmax
       and an opening angle equal the cherenkov angle at that height given by alpha = arccos(1 / n(h)) */
    double CalculateR0() const;

    //New, "hard p", /wo d750
    static double GetSigHard(const double dmax) {
      return 0.16460107 * std::pow(dmax - 5000, 0.69675388) + 39.46602294;
    }

    static double GetPHard(const double dmax) {
      return -4.91997573e-01 * exp(-3.88166824e-05 * dmax) + 1.85788594e+00;
    }

    static double GetR02Hard(const double dmax) {
      return 5.48029929e-01 + 7.06584317e-07 * dmax + 5.06592690e+07 / dmax / dmax;
    }

    static double GetArelHard(const double dmax) {
      return 7.57300343e-01 + 7.76904953e-07 * dmax + 1.63545932e+07 / dmax / dmax;
    }

    /*// New p soft, only dmax, with d750
    double GetSig(const double dmax) const {
      return 1.39506207e-01 * std::pow(dmax - 5000, 7.09861340e-01) +
             5.39164661e+01 - 1.84052852e-03 * (dmax - fShowerData.fDistanceTo750);
    }

    static double GetPSoft(const double dmax) {
      return 1.63989873e+02 * exp(-2.75777834e-05 * dmax) + 6.81576551e+01;
    }

    static double GetR02(const double dmax) {
      return 5.41470307e-01 + 7.18512356e-07 * dmax + 6.62411621e+07 / dmax / dmax;
    }

    double GetArel(const double dmax) const {
      return 7.52546909e-01 + 7.92709635e-07 * dmax + 2.23641574e+07 / dmax / dmax
        + 2.99174458e-06 * (dmax - fShowerData.fDistanceTo750);
    }*/

    // New p soft, only dmax, without d750
    static double GetSigSoft(const double dmax)
    { return 0.13176183 * std::pow(dmax - 5000, 0.71437054) + 56.30941015; }

    static double GetPSoft(const double dmax)
    { return 1.54942405e+02 * exp(-2.50177591e-05 * dmax) + 6.49133050e+01; }

    static double GetR02Soft(const double dmax)
    { return 5.51709206e-01 + 6.87661913e-07 * dmax + 6.61581738e+07 / dmax / dmax; }

    static double GetArelSoft(const double dmax)
    { return 7.57449757e-01 + 7.68399447e-07 * dmax + 1.98026585e+07 / dmax / dmax; }

    double GetSig(const double dmax) const
    { return (fFitConfig.fUseSoftExponent) ? GetSigSoft(dmax) : GetSigHard(dmax); }

    double GetP(const double dmax) const
    { return (fFitConfig.fUseSoftExponent) ? GetPSoft(dmax) : GetPHard(dmax); }

    double GetR02(const double dmax) const
    { return (fFitConfig.fUseSoftExponent) ? GetR02Soft(dmax) : GetR02Hard(dmax); }

    double GetArel(const double dmax) const
    { return (fFitConfig.fUseSoftExponent) ? GetArelSoft(dmax) : GetArelHard(dmax); }


  private:
    // specific fit parameter
    mutable double fR0 = 0;
    mutable double fSig = 0;
    mutable double fP = 0;
    mutable double fR02 = 0;
    mutable double fArel = 0;
    const double fSlope = 5;  // fixed

    void
    UpdateParameterParam(const std::vector<double>& pars)
      const override
    {
      // Fit parameter
      fEgeo = pars[0];
      fDistanceToXmax = pars[1];
      fCoreX = pars[2];
      fCoreY = pars[3];

      // calculate density from distance
      // the -1 is needed because of convention of the table (opposit direction of the init vector)
      fHeight = fShowerData.fHeightFromDistance.Y(-fDistanceToXmax);
      fDensityAtXmax = fShowerData.fDensityFromHeight.Y(fHeight) / (kg / m3);

      // parametrized parameter
      fSig = GetSig(fDistanceToXmax);
      fP = GetP(fDistanceToXmax);
      fR02 = GetR02(fDistanceToXmax);
      fArel = GetArel(fDistanceToXmax);

      fR0 = CalculateR0();  // needs fDensityAtXmax to be updated
      fRnorm = 5 * fR0;
      fNorm = GetNormalization();  // needs fRnorm to be updated
    }
  };

}


#endif