Deprecated/RdLDFFitter/LikelihoodFunction.cc

Go to the documentation of this file.

#include "LikelihoodFunction.h"

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

#include <Minuit2/FCNBase.h>

#include <utl/CoordinateSystemPtr.h>
#include <utl/Vector.h>
#include <utl/BasicVector.h>
#include <utl/Point.h>
#include <utl/Math.h>
#include <utl/PhysicalFunctions.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/Probability.h>
#include <utl/AugerUnits.h>

#include "TF1.h"
#include "TF2.h"
#include "TMath.h"

namespace RdLDFFitter {

  LDFLikelihoodFunction::LDFLikelihoodFunction(const FitConfig fitconfig,
      const EventFitData eventData, const std::vector<StationFitData>& stationData,
      const std::vector<ScintillatorFitData>& scintillatorData,
      const utl::Vector magneticFieldVector) :
          fFitConfig(fitconfig),
          fEventData(eventData),
          fStationData(stationData),
          fScintillatorData(scintillatorData),
          fTheErrorDef(1.0)
  {
    fMagneticFieldVector = magneticFieldVector;
  }

  /*
   A different LDF model can be implemented easily by defining a new Likelihood of Chi2 function.
   This contribution to the combined neg. log. Likelihood needs to be added in the station loop to
   the variable "negLogLikelihood". Just follow the example for LDFModel==1.
   The parameters of the LDF function needs to be added to the minimizer in RdLDFFitter.cc
   Parameters are:
     0, 1, 2 = (coreX, coreY, coreZ) in local coordinate system
     3 = a (charge excess strength)
     4, 5, 6 = Scintillator LDF parameters
     7, 8,... = optional LDF parameters
   */

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

    utl::Point core(pars[0], pars[1], pars[2], fEventData.fLocalCS);
    const double chargeExcessStrength = pars[3];

    double negLogLikelihood = 0;

    // std::cout << "station loop " << std::endl;

    // loop over all stations
    for(std::vector<StationFitData>::const_iterator sIt = fStationData.begin(), end =
        fStationData.end(); sIt != end; ++sIt) {

      double signal = sIt->signal;
      double signalError = sIt->signalError;
      // correct signal amplitude for charge excess and geo-magnetic emission
      if(fFitConfig.useChargeExcessCorrectionInLDFFit and fFitConfig.LDFModel) { // calculate correction only if LDF fit should be performed
        double correctionFactor = utl::RadioGeometryUtilities::GetSignalCorrectionFactor(core,
            fEventData.fShowerAxis, sIt->antennaPosition, fMagneticFieldVector,
            chargeExcessStrength);
        signal *= correctionFactor;
        signalError *= correctionFactor;
      }

      // chi2 LDF fit for LDF model 1
      if(fFitConfig.LDFModel == 1) {
        negLogLikelihood += GetChi2LDFModel1(pars[7], pars[8], core, sIt->antennaPosition,
            fEventData.fShowerAxis, signal, signalError);
      }

      // polarization estimator
      if(fFitConfig.usePolarisation) {
        negLogLikelihood += -2
            * log(
                GetLikelihoodPolarisation(core, fEventData.fShowerAxis, sIt->antennaPosition,
                    sIt->efield, sIt->lorentzAngleError, chargeExcessStrength));
      }

    } // loop over stations

    // add SD core Likelihood to combined Likelihood
    if(fFitConfig.useSDCoreToImproveRadioCore) {
      negLogLikelihood += -2
          * log(
              GetSDCoreLikelihood(core.GetX(fEventData.fLocalCS), core.GetY(fEventData.fLocalCS),
                  fEventData.SDCore.GetX(fEventData.fLocalCS),
                  fEventData.SDCore.GetY(fEventData.fLocalCS), fEventData.SDCoreXError,
                  fEventData.SDCoreYError, fEventData.SDCoreXYCorrelation));
    }

    // loop over scnintillators
    if(fFitConfig.useScintillator) {
      for(std::vector<ScintillatorFitData>::const_iterator sIt = fScintillatorData.begin(), end =
          fScintillatorData.end(); sIt != end; ++sIt) {
        negLogLikelihood += GetScintillatorLDFLikelihood(core, fEventData.fShowerAxis, sIt->scintillatorPosition,
                      pars[4], pars[5], pars[6], sIt->signal, sIt->silent);
      }
    }

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

//    double Up() const { return 0.5; }

  double LDFLikelihoodFunction::GetChi2LDFModel1(const double A, const double R0,
      const utl::Point& core, const utl::Point& stationPosition, const utl::Vector& showeraxis,
      const double signal, const double signalError) const
  {
    double distanceToShowerAxis = utl::RadioGeometryUtilities::GetDistanceToAxis(showeraxis,
        core, stationPosition);
    TF1 fLDF = TF1("fLDF", "[0]*exp(-x/[1])");
    fLDF.SetParameters(A, R0);
    double LDFSignal = fLDF.Eval(distanceToShowerAxis);
    return pow((signal - LDFSignal) / signalError, 2);
  }

  double LDFLikelihoodFunction::GetLikelihoodPolarisation(const utl::Point& core,
      const utl::Vector& showeraxis, const utl::Point& stationPosition, const utl::Vector& EField,
      const double& lorentzAngleError, const double& chargeExcessStrength) const
  {
//    utl::Vector chargeExcessVector = utl::RadioGeometryUtilities::GetChargeExcessVector(core,
//        showeraxis, stationPosition);
//    utl::Vector lorentzVector = utl::RadioGeometryUtilities::GetLorentzVector(showeraxis,
//        fMagneticFieldVector);
//    utl::Vector vEexpected = lorentzVector + chargeExcessStrength * chargeExcessVector;
    double delta = utl::RadioGeometryUtilities::GetAngleToEFieldExpectation2D(EField, core,
                                                                          showeraxis,
                                                                          stationPosition,
                                                                          fMagneticFieldVector,
                                                                          chargeExcessStrength,
                                                                          fEventData.fLocalCS);
//    double x = TMath::RadToDeg() * showeraxis.GetTheta(fEventData.fLocalCS);
//    double mean = 12.6091 + x * -0.323591;
//    if(x < 40) {
//      mean = 0.;
//    }
//    double sigma = pow(x, 0) * 4.33775  + pow(x, 1) * -0.0465731  + pow(x, 2) * 0.00752399  + pow(x, 3) * -0.000265992  + pow(x, 4) * -3.97691e-06  + pow(x, 5) * 2.19571e-07  + pow(x, 6) * -1.86878e-09 ;
//    mean *= TMath::DegToRad();
//    sigma *= TMath::DegToRad();

//    double delta = utl::Angle(vEexpected, EField);
    //double tmp = utl::Probability::GetNormalPDF(delta, 0, 1 * utl::deg);
//    double tmp = utl::Probability::GetNormalPDF(delta, 0, lorentzAngleError);

//    double tmp = utl::Probability::GetVonMisesPDF(delta, 0, 1 / (1 * utl::deg * 1 * utl::deg));
    double tmp = utl::Probability::GetVonMisesPDF(delta, 0.,
                                                  1 / (lorentzAngleError * lorentzAngleError));

    //double tmp = utl::Probability::GetNormalPDF(delta, mean, sigma);
    // std::cout << "Likelihood Pol: p("
    //          << delta / utl::deg << ", 0, " <<  lorentzAngleError / utl::deg << ")=" << tmp << std::endl;
    // double tmp = utl::Probability::GetFisherPDF(delta, 1 / (sigmaDelta * sigmaDelta));
    return (tmp == 0) ? std::numeric_limits<double>::min() : tmp;


//    TF1 fRayleigh = TF1("fRayleigh", "x/[0]*exp(-x*x/(2*[0]*[0]))");
//    fRayleigh.SetParameter(0, sigmaDelta);
//    double tmp = fRayleigh.Eval(delta);
//    return (tmp == 0) ? std::numeric_limits<double>::min() : tmp;
  }

  double LDFLikelihoodFunction::GetSDCoreLikelihood(const double coreX, const double coreY,
      const double SDCoreX, const double SDCoreY, const double SDCoreXError,
      const double SDCoreYError, const double SDCoreXYCorrelation) const
  {
    TF2 fLikelihood = TF2("2DimLikelihood",
        "2.*( ((x-[0])/[1])**2 + ((y-[2])/[3])**2 - 2.*[4]*((x-[0])/[1])*((y-[2])/[3]) )");
    fLikelihood.SetParameters(SDCoreX, SDCoreXError, SDCoreY, SDCoreYError, SDCoreXYCorrelation);

    return fLikelihood.Eval(coreX, coreY);
  }

  double LDFLikelihoodFunction::GetScintillatorLDFLikelihood(const utl::Point& core,
                                                             const utl::Vector& showeraxis,
                                                             const utl::Point& stationPosition,
                                                             const double N_charged_particles,
                                                             const double showerAge,
                                                             const double moliereRadius,
                                                             const double signal,
                                                             const bool silent) const
  {
    
    double distanceToShowerAxis = utl::RadioGeometryUtilities::GetDistanceToAxis(
        showeraxis, core, stationPosition);
        
    double scint_likelihood_value = 0.;
    double sum_log_signal = 0.;

    TF1 NKG = TF1("NKG","[0]*TMath::Gamma(4.5-[1])/(2.*TMath::Pi()*[2]*[2]*TMath::Gamma([1])*TMath::Gamma(4.5-(2.*[1])))*TMath::Power((x/[2]),([1]-2.))*TMath::Power((1+(x/[2])),[1]-4.5)");
    
    NKG.SetParameter(0, N_charged_particles);
    NKG.SetParameter(1, showerAge);
    NKG.SetParameter(2, moliereRadius);
    
    if(silent){
        scint_likelihood_value = -NKG.Eval(distanceToShowerAxis)+log(1+NKG.Eval(distanceToShowerAxis));
    }
    else{
        for(int i=1; i<=static_cast<int>(round(signal)); i++){
            sum_log_signal+=log(i);
        }  
        scint_likelihood_value = signal*log(NKG.Eval(distanceToShowerAxis)) 
            - NKG.Eval(distanceToShowerAxis) - sum_log_signal;
    }

    return -2. * scint_likelihood_value; // factor two to be -2*ln(L)
  }

}