FitModels.cc

Go to the documentation of this file.

#include "FitModels.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"

#include "adst/UtilityFunctions.h"

namespace RdGlobalFit {

RdGlobalFitMinimizationCriterion::RdGlobalFitMinimizationCriterion(
    FitConfig& fitconfig, const EventFitData eventData, std::vector<StationFitData>& stationData,
    const utl::Vector magneticFieldVector, const LDFConstsTable ldfConstsTable,
    calcLDFConsts& calcLDFConsts) :
    fFitConfig(fitconfig),
    fEventData(eventData),
    fStationData(stationData),
    fTheErrorDef(1.0),
    fVerbose(false),
    fLDFConstsTable(ldfConstsTable),
    fcalcLDFConsts(calcLDFConsts)
{
  fMagneticFieldVector = magneticFieldVector;
}

void RdGlobalFitMinimizationCriterion::calcShowerCoordinates(const double arrivalDirection_phi,
                                                             const double arrivalDirection_theta,
                                                             const double core_x,
                                                             const double core_y,
                                                             const double core_z,
                                                             utl::Vector magneticFieldVector) const
{
  //Create CS with core as origin
  utl::Point dynamic_core(core_x, core_y, core_z, fEventData.fLocalCS);
  utl::CoordinateSystemPtr dynamicCS = fwk::LocalCoordinateSystem::Create(dynamic_core);
  //Calculate arrivaldirection vector in LocalSystem (needs proper calculation)
  utl::Vector dynamic_Showeraxis(sin(arrivalDirection_theta) * cos(arrivalDirection_phi),
                                 sin(arrivalDirection_theta) * sin(arrivalDirection_phi),
                                 cos(arrivalDirection_theta), fEventData.fLocalCS);
  //calculate positions in shower frame
  utl::RadioGeometryUtilities transformation = utl::RadioGeometryUtilities(dynamic_Showeraxis,
                                                                           dynamicCS,
                                                                           magneticFieldVector);
  //Calculating station-positions in shower-plane
  for (std::vector<StationFitData>::iterator sIt = fStationData.begin(), end = fStationData.end();
      sIt != end; ++sIt) {
    const utl::Point sPos(sIt->stationPosition);
    double X_VxB = 0.;
    double Y_VxB = 0.;
    double Z_VxB = 0.;
    transformation.GetVectorInShowerPlaneVxB(X_VxB, Y_VxB, Z_VxB, sPos);
    //Filling FitEventData with the calculated stationepositions in the dynamic userparameter dependend coordinatesystem
    sIt->x_vxB = X_VxB;
    sIt->y_vxB = Y_VxB;
    sIt->z_vxB = Z_VxB;
  }
}

double RdGlobalFitMinimizationCriterion::operator()(const std::vector<double>& pars) const
{
  //creating dynamic showerCS depending on Fit Parameters of direction and core and saving stationPositions in this frame in fStationData.
  calcShowerCoordinates(pars[14], pars[15], pars[1], pars[2], pars[3], fMagneticFieldVector);
  calc2dLDFConstants(pars[15]);

  double Rmax = 0;
  double DXmax = 0;

  // The following is only important if fFitConfig.fitGammaAndSigmaPlusIndependently is false.
  // Its functioning is not validated and it may be removed in future.
  // Momentarily, its only commented out.

  //  if(fFitConfig.fitShowerMaxAsRMax){
  //    DXmax = RdGlobalFitMinimizationCriterion::DXmaxFromRmax(pars[11], pars[15], 1560.);
  //    Rmax = pars[11];
  //  }
  //
  //  if(fFitConfig.fitShowerMaxAsDXMax){
  //    DXmax = pars[11];
  //  }

  double negLogLikelihood = 0;
  double twoDLDFLikelyhood = 0.;
  double arrivalTimeLikelihood = 0.;

  for (std::vector<StationFitData>::iterator sIt = fStationData.begin(), end = fStationData.end();
      sIt != end; ++sIt) {

    if (fFitConfig.fitTwoDLDF == true) {
      double stationLikelihood = calc2dLDFMinCrit(pars[0], pars[4], DXmax, pars[5], pars[6],
                                                  pars[7], pars[8], pars[9], sIt->x_vxB, sIt->y_vxB,
                                                  sIt->energyFluence, sIt->energyFluenceError, sIt);
      negLogLikelihood += stationLikelihood;
      twoDLDFLikelyhood += stationLikelihood;
    }

    if (fFitConfig.fitArrivalTime == true) {
      double stationLikelihood = calcTimeMinCrit(pars[10], Rmax, pars[12], pars[13], sIt->x_vxB,
                                                 sIt->y_vxB, sIt->z_vxB, sIt->signalTime,
                                                 sIt->signalTimeError, sIt);
      negLogLikelihood += stationLikelihood;
      arrivalTimeLikelihood += stationLikelihood;
    }
  }
  fArrivalTimeLikelihood = arrivalTimeLikelihood;
  fTwoDLDFLikelyhood = twoDLDFLikelyhood;
  return negLogLikelihood;
}

double RdGlobalFitMinimizationCriterion::sigmaPlusFromDXmax(double DXmax, double a1, double a2,
                                                            double a3, double a4) const
{
  double sigma_result;
  sigma_result = (1.0L / 6.0L)
      * (-pow(2, 2.0L / 3.0L) * pow(a1, 2)
          * pow(
              (pow(a1, 3)
                  * sqrt(
                      (4 * pow(3 * a1 * a3 - pow(a2, 2), 3)
                          + pow(27 * pow(a1, 2) * (DXmax - a4) + 9 * a1 * a2 * a3 - 2 * pow(a2, 3),
                                2)) / pow(a1, 6)) + 27 * pow(a1, 2) * (-DXmax + a4)
                  - 9 * a1 * a2 * a3 + 2 * pow(a2, 3)) / pow(a1, 3),
              2.0L / 3.0L)
          - 2 * a1 * a2
              * pow(
                  (pow(a1, 3)
                      * sqrt(
                          (4 * pow(3 * a1 * a3 - pow(a2, 2), 3)
                              + pow(
                                  27 * pow(a1, 2) * (DXmax - a4) + 9 * a1 * a2 * a3
                                      - 2 * pow(a2, 3),
                                  2)) / pow(a1, 6)) + 27 * pow(a1, 2) * (-DXmax + a4)
                      - 9 * a1 * a2 * a3 + 2 * pow(a2, 3)) / pow(a1, 3),
                  1.0L / 3.0L) + 2 * pow(2, 1.0L / 3.0L) * (3 * a1 * a3 - pow(a2, 2)))
      / (pow(a1, 2)
          * pow(
              (pow(a1, 3)
                  * sqrt(
                      (4 * pow(3 * a1 * a3 - pow(a2, 2), 3)
                          + pow(27 * pow(a1, 2) * (DXmax - a4) + 9 * a1 * a2 * a3 - 2 * pow(a2, 3),
                                2)) / pow(a1, 6)) + 27 * pow(a1, 2) * (-DXmax + a4)
                  - 9 * a1 * a2 * a3 + 2 * pow(a2, 3)) / pow(a1, 3),
              1.0L / 3.0L));
  return sigma_result;
}

double RdGlobalFitMinimizationCriterion::gammaFromRmax(double DXmax, double alpha_1, double alpha_2,
                                                       double alpha_3, double c) const
{
  double gamma_result;
  gamma_result = (1.0L / 2.0L)
      * (-alpha_2 * c
          + sqrt(c * (4 * DXmax * alpha_1 - 4 * alpha_1 * alpha_3 * c + pow(alpha_2, 2) * c)))
      / (alpha_1 * c);
  return gamma_result;
}

double RdGlobalFitMinimizationCriterion::calc2dLDFMinCrit(
    const double Aplus, const double sigmaPlus, const double DXmax, const double C2Theta,
    const double C1Theta, const double CTheta, const double C1, const double C2, const double x_vxB,
    const double y_vxB, const double signal, const double signalError,
    std::vector<StationFitData>::iterator sIt) const
{
  //Calculates for given Parameters the EnergieDensity for a given Station, compares it in some defined way to the measured Data and returns a minimization criterion.
  // (The parameters are given by the class. The class itself is called with the definition of the fit (which parameters are set fixed, which range, etc..))

  //Parameters for correlation between DXmax and sigma from J. Schulz.
  //ToDo Parameter should be set with xml.
  double a1 = 2.1 * pow(10, -4);
  double a2 = -0.081;
  double a3 = 13.5;
  double a4 = -352.;

  double sigmaPl;
  double d_xmax = DXmax;
  if (!fFitConfig.fitGammaAndSigmaPlusIndependently) {
    sigmaPl = sigmaPlusFromDXmax(d_xmax, a1, a2, a3, a4);
  } else {
    sigmaPl = sigmaPlus;
  }
  double c1Theta = C1Theta;
  double c2Theta = C2Theta;
  double c1 = C1;
  double c2 = C2;
  double cTheta = CTheta;
  c1Theta = fcalcLDFConsts.C1theta;
  c2Theta = fcalcLDFConsts.C2theta;
  c1 = fcalcLDFConsts.C3;
  c2 = fcalcLDFConsts.C4;
  cTheta = fcalcLDFConsts.A_scale;

  TF2 LDF2DFit =
      TF2("LDF2DFit",
          "[0]*TMath::Exp(-((x-([1]+[5]))**2+(y-[2])**2)/[3]**2)-[6]*[0]*TMath::Exp(-((x-([1]+[4]))**2+(y-[2])**2)/TMath::Exp([7]+[8]*[3])**2)");
  LDF2DFit.SetParameter(0, Aplus);
  LDF2DFit.SetParameter(1, 0.);
  LDF2DFit.SetParameter(2, 0.);
  LDF2DFit.SetParameter(3, sigmaPl);
  LDF2DFit.SetParameter(4, c2Theta);
  LDF2DFit.SetParameter(5, c1Theta);
  LDF2DFit.SetParameter(6, cTheta);
  LDF2DFit.SetParameter(7, c1);
  LDF2DFit.SetParameter(8, c2);

  double weight = signalError;
  double calculatedEnergieDensity = LDF2DFit.Eval(x_vxB, y_vxB);
  double minimizationCriterium = pow((signal - calculatedEnergieDensity) / weight, 2);
  sIt->fittedEnergyFluence = calculatedEnergieDensity;
  if (fVerbose) {
    std::cout << x_vxB << "\t" << y_vxB << "\t" << signal << "\t" << calculatedEnergieDensity
        << "\t" << minimizationCriterium << std::endl;
  }

  return minimizationCriterium;
}

double RdGlobalFitMinimizationCriterion::calcTimeMinCrit(
    const double gamma, const double DXmax, const double b, const double t0, const double x_vxB,
    const double y_vxB, const double z_vxB, const double signalTime, const double signalTimeError,
    std::vector<StationFitData>::iterator sIt) const
{
  TF1 arrivalTimeFit = TF1("ArrivalTimeFit",
                           "TMath::Sqrt(1 + x * x / ([0] * [0])) * [1] - [1] + [2]");
  double gammaTilde;

  double c = 10000.;
  double alpha_3 = -1.49048;
  double alpha_2 = -0.0371733;
  double alpha_1 = 0.00086932;
  if (!fFitConfig.fitGammaAndSigmaPlusIndependently) {
    gammaTilde = gammaFromRmax(DXmax, alpha_1, alpha_2, alpha_3, c);
  } else {
    gammaTilde = gamma;
  }
  arrivalTimeFit.SetParameter(0, gammaTilde);
  arrivalTimeFit.SetParameter(1, b);
  arrivalTimeFit.SetParameter(2, t0);

  double weight = signalTimeError;
  double r = sqrt(x_vxB * x_vxB + y_vxB * y_vxB);
  double calculatedArrivalTime = arrivalTimeFit.Eval(r);
  double measuredArrivalTime = signalTime;
  double timeDrift = 0;  // difference between core-time and eventtime
  double relativePulseTime = z_vxB / (utl::kSpeedOfLight / utl::nanosecond);
  double planeTime = timeDrift + relativePulseTime;
  double measuredProjTime = measuredArrivalTime - planeTime;
  double minimizationCriterium = pow((measuredProjTime - calculatedArrivalTime) / weight, 2);
  sIt->fittedProjTime = calculatedArrivalTime;
  sIt->measuredProjTime = measuredProjTime;
  if (fVerbose) {
    std::cout << x_vxB << "\t" << y_vxB << "\t" << measuredProjTime << "\t" << calculatedArrivalTime
        << "\t" << minimizationCriterium << std::endl;
  }
  return minimizationCriterium;
}

void RdGlobalFitMinimizationCriterion::calc2dLDFConstants(double zenith) const
{

  fcalcLDFConsts.C0 = 0.41;

  if ((zenith) < 10. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_0_10;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_0_10;
  } else if ((zenith) < 20. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_10_20;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_10_20;
  } else if ((zenith) < 30. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_20_30;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_20_30;
  } else if ((zenith) < 40. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_30_40;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_30_40;
  } else if ((zenith) < 50. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_40_50;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_40_50;
  } else if ((zenith) < 55. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_50_55;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_50_55;
    fcalcLDFConsts.A_scale = fLDFConstsTable.A_scale_50_55;
    fcalcLDFConsts.C0 = 0.46;
  } else {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_55;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_55;
    fcalcLDFConsts.A_scale = fLDFConstsTable.A_scale_55;
    fcalcLDFConsts.C0 = 0.71;
  }
}
}