Deprecated/RdLDFFitter/LikelihoodFunction.h

Go to the documentation of this file.

#ifndef _RdLDFFitter_LikelihoodFunctions_h_
#define _RdLDFFitter_LikelihoodFunctions_h_

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

#include <Minuit2/FCNBase.h>

#include <utl/CoordinateSystemPtr.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/Math.h>
#include <utl/PhysicalFunctions.h>

namespace RdLDFFitter {

  struct FitConfig {
    int LDFModel;
    bool useStationsWithoutSignal;
    bool useSDCoreToImproveRadioCore;   // Likelihood of SD core position will be added to combined Likelihood
    bool useChargeExcessCorrectionInLDFFit; // signal strength will be corrected for charge excess contribution in LDF fit
    bool usePolarisation;   // polarization of radio signal will be used to estimate core position (polarization Likelihood will be added to combined Likelihood)
    double chargeExcessStrength;    // value of the charge excess strength relative to geo-magnetic strength for sin(alpha) = 1
    bool fitChargeExcess;     // treat charge excess strength as free fit parameters?
    bool useScintillator;     // use scintiallator signals to fit a particle LDF

    FitConfig() :
            LDFModel(1),
            useStationsWithoutSignal(false),
            useSDCoreToImproveRadioCore(false),
            useChargeExcessCorrectionInLDFFit(true),
            usePolarisation(true),
            chargeExcessStrength(0.14),
            fitChargeExcess(false),
            useScintillator(false)
    {
    }
  };

  struct StationFitData {
    utl::Point antennaPosition;
    double signal;
    double signalError;
    utl::Vector efield;
    double lorentzAngleError;

    StationFitData() :
        antennaPosition(), signal(0), signalError(0), efield(), lorentzAngleError(0)
    {
    }
  };

  struct ScintillatorFitData {
    utl::Point scintillatorPosition;
    double signal;
    bool silent;

    ScintillatorFitData() :
      scintillatorPosition(), signal(0.), silent(false)
    {
    }
  };

  struct EventFitData {
    utl::Vector fShowerAxis;
    utl::CoordinateSystemPtr fLocalCS;

    utl::Point SDCore; // the SD core only needs to be specified if the useSDCoreToImproveRadioCore flag is set true
    double SDCoreXError;
    double SDCoreYError;
    double SDCoreXYCorrelation;

    EventFitData() :
        fShowerAxis(), fLocalCS(), SDCore(), SDCoreXError(0), SDCoreYError(0), SDCoreXYCorrelation(0)
    {
    }

  };

  class LDFLikelihoodFunction: public ROOT::Minuit2::FCNBase {

  public:
    LDFLikelihoodFunction(const FitConfig ldfconfig, const EventFitData eventData,
        const std::vector<StationFitData>& stationData,
        const std::vector<ScintillatorFitData>& scintillatorData,
        const utl::Vector magneticFieldVector);

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

    double Up() const
    {
      return fTheErrorDef;   // needs to be one because the neg. log. likelihood is defined as -2*log(Likelihood)
    }

    void SetErrorDef(double def) {
      fTheErrorDef = def;
    }

  protected:

  private:
    double 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 GetLikelihoodPolarisation(const utl::Point& core, const utl::Vector& showeraxis,
        const utl::Point& stationPosition, const utl::Vector& EField,
        const double& lorentzAngleError, const double& chargeExcessStrength) const;

    double GetSDCoreLikelihood(const double coreX, const double coreY, const double SDCoreX,
        const double SDCoreY, const double SDCoreXError, const double SDCoreYError,
        const double SDCoreXYCorrelation) const;

    double 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;

    const FitConfig fFitConfig;

    const EventFitData fEventData;
    const std::vector<StationFitData>& fStationData;
    const std::vector<ScintillatorFitData>& fScintillatorData;

    utl::Vector fMagneticFieldVector;

    double fTheErrorDef;
  };

}

#endif