FitModels.h

Go to the documentation of this file.

#ifndef _FitModels_h_
#define _FitModels_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>
#include <utl/Vector.h>
#include <fwk/LocalCoordinateSystem.h>

namespace RdGlobalFit {

struct FitConfig {
  bool fitTwoDLDF;
  bool fitArrivalTime;
  bool fixCore;
  bool fitGammaAndSigmaPlusIndependently;
  bool fitShowerMaxAsRMax;
  bool fitShowerMaxAsXMAX;
  bool fitShowerMaxAsDXMax;
  bool SaveCoreContourData;
  bool printMinuitScans;
  std::string useCS;
  FitConfig() :
      fitTwoDLDF(false),
      fitArrivalTime(false),
      fixCore(false),
      fitGammaAndSigmaPlusIndependently(false),
      fitShowerMaxAsRMax(false),
      fitShowerMaxAsXMAX(false),
      fitShowerMaxAsDXMax(false),
      SaveCoreContourData(false),
      printMinuitScans(false),
      useCS("MAN")
  {
  }
};

struct LDFConstsTable {
  double A_scale;
  double C3;
  double C4;
  double C1theta_0_10;
  double C2theta_0_10;
  double C1theta_10_20;
  double C2theta_10_20;
  double C1theta_20_30;
  double C2theta_20_30;
  double C1theta_30_40;
  double C2theta_30_40;
  double C1theta_40_50;
  double C2theta_40_50;
  double C1theta_50_55;
  double C2theta_50_55;
  double A_scale_50_55;
  double C1theta_55;
  double C2theta_55;
  double A_scale_55;
  LDFConstsTable() :
      A_scale(0.),
      C3(0.),
      C4(0.),
      C1theta_0_10(0.),
      C2theta_0_10(0.),
      C1theta_10_20(0.),
      C2theta_10_20(0.),
      C1theta_20_30(0.),
      C2theta_20_30(0.),
      C1theta_30_40(0.),
      C2theta_30_40(0.),
      C1theta_40_50(0.),
      C2theta_40_50(0.),
      C1theta_50_55(0.),
      C2theta_50_55(0.),
      A_scale_50_55(0.),
      C1theta_55(0.),
      C2theta_55(0.),
      A_scale_55(0.)
  {
  }
};

struct calcLDFConsts {
  double A_scale;
  double C0;
  double C3;
  double C4;
  double C1theta;
  double C2theta;
  calcLDFConsts() :
      A_scale(0.),
      C0(0.),
      C3(0.),
      C4(0.),
      C1theta(0.),
      C2theta(0.)
  {
  }
};

struct StationFitData {
  int stationId;
  utl::Point stationPosition;
  double x_vxB;
  double y_vxB;
  double z_vxB;
  double energyFluence;
  double energyFluenceError;
  double signalTime;
  double signalTimeError;
  double fittedProjTime;
  double fittedEnergyFluence;
  double measuredProjTime;

  StationFitData() :
      stationId(0),
      stationPosition(),
      x_vxB(0),
      y_vxB(0),
      z_vxB(0),
      energyFluence(0),
      energyFluenceError(0),
      signalTime(0),
      signalTimeError(0),
      fittedProjTime(0),
      fittedEnergyFluence(0),
      measuredProjTime(0)
  {
  }
};

struct FitData {
  bool isFix;
  bool hasLimits;
  std::string startValueFrom;
  double startValue;
  double startError;
  double lowLimit;
  double highLimit;

  FitData() :
      isFix(false),
      hasLimits(false),
      startValueFrom("MAN"),
      startValue(0),
      startError(0),
      lowLimit(0),
      highLimit(0)
  {
  }

  FitData(bool isfix, bool haslimits, double startvalue, double starterror, double lowlimit,
          double highlimit) :
      isFix(isfix),
      hasLimits(haslimits),
      startValue(startvalue),
      startError(starterror),
      lowLimit(lowlimit),
      highLimit(highlimit)
  {
  }
};

struct RdGlobalFitData {
  FitData Aplus;
  FitData core_x;
  FitData core_y;
  FitData core_z;
  FitData sigmaPlus;
  FitData C2Theta;
  FitData C1Theta;
  FitData CTheta;
  FitData C1;
  FitData C2;
  FitData gamma;
  FitData b;
  FitData t0;
  FitData arrivalDirection_phi;
  FitData arrivalDirection_theta;
  FitData showerMax;
  bool fitSucess;

  RdGlobalFitData() :
      Aplus(),
      core_x(),
      core_y(),
      core_z(),
      sigmaPlus(),
      C2Theta(),
      C1Theta(),
      CTheta(),
      C1(),
      C2(),
      gamma(),
      b(),
      t0(),
      arrivalDirection_phi(),
      arrivalDirection_theta(),
      showerMax(),
      fitSucess(false)
  {
  }
};

struct EventFitData {
  utl::CoordinateSystemPtr fLocalCS;

  EventFitData() :
      fLocalCS()
  {
  }

};

class RdGlobalFitMinimizationCriterion : public ROOT::Minuit2::FCNBase {

 public:
  RdGlobalFitMinimizationCriterion(FitConfig& fitconfig, const EventFitData eventData,
                                   std::vector<StationFitData>& stationData,
                                   const utl::Vector magneticFieldVector,
                                   const LDFConstsTable ldfConstsTable,
                                   calcLDFConsts& calcLDFConsts);

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

  double Up() const
  {
    return fTheErrorDef;
  }

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

  void setVerbose(bool verbose)
  {
    fVerbose = verbose;
  }

  double getArrivalTimeLikelihood()
  {
    return fArrivalTimeLikelihood;
  }

  double getTwoDLDFLikelyhood()
  {
    return fTwoDLDFLikelyhood;
  }

  std::vector<StationFitData>& getStationData()
  {
    return fStationData;
  }

 protected:

 private:
  void 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;

  double calcTimeMinCrit(const double gamma, const double DXmax, const double b, const double t0,
                         const double x_vxB, const double y_vxvxB, const double z_vxvxB,
                         const double signalTime, const double signalTimeError,
                         std::vector<StationFitData>::iterator sIt) const;

  double calc2dLDFMinCrit(const double Aplus, const double sigma_max, 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_vxvxB, const double signal, const double signalError,
                          std::vector<StationFitData>::iterator sIt) const;

  void calc2dLDFConstants(double zenith) const;

  // These functions are needed, if one wants to fit the shower maximum "directly" in the global fit.
  // Its functioning is not yet validated. May be removed at some point.
  double sigmaPlusFromDXmax(double DXmax, double a1, double a2, double a3, double a4) const;
  double gammaFromRmax(double DXmax, double alpha_1, double alpha_2, double alpha_3,
                       double c) const;

  FitConfig& fFitConfig;

  const EventFitData fEventData;
  std::vector<StationFitData>& fStationData;
  utl::Vector fMagneticFieldVector;
  double fTheErrorDef;
  bool fVerbose;
  const LDFConstsTable fLDFConstsTable;
  calcLDFConsts& fcalcLDFConsts;
  mutable double fTwoDLDFLikelyhood;
  mutable double fArrivalTimeLikelihood;
};

}

#endif