UnivFitterKG.h

Go to the documentation of this file.

#ifndef _UnivFitterKG_UnivFitterKG_h_
#define _UnivFitterKG_UnivFitterKG_h_

#include <tls/UnivParam.h>
#include <tls/UnivParamTime.h>
#include <tls/Atmosphere.h>
#include <tls/UnivCalibConstants.h>

#include "Minuit2/FCNBase.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/FunctionMinimum.h"

#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/TimeStamp.h>

#include <tls/UnivTimeKG.h>

#include <FitterUtil.h>


namespace UnivFitterKG {

  enum InitMode {
    eMC,
    eSD,
    eFD
  };


  struct fitConstraint {
    double mean = 0;
    double sigma = 0;
  };


  struct StationFitData {
    int stationId = 0;
    double totalSignal = 0;
    double traceSum = 0;
    double Smaxbin = 0;
    double recoveredSignal = 0;
    double recoveredSignalUncertainty = 0;
    bool isSaturated = false;
    bool useUnsaturatedTrace = false;

    double muonSignal = 0;
    double muonSignalError = 0;
    std::vector<double> muonTimes; // time referred to MC core
    bool useForMuonTimeFit = false;

    std::vector<double> trace;
    double binwidth = FitterUtil::tracebinning;
    int firstFitBin = 0;
    int lastFitBin = 0;

    utl::TimeStamp startTime;
    int startBin = 0;
    int stopBin = 0;

    double t10 = 0;
    double t50 = 0;
    double t90 = 0;

    utl::Point position;

    bool isDense = false;
    double denseRho = 0;
    double densePsi = 0;

    bool useForLDFFit = false;
    bool useForShapeFit = false;
    bool useForStartTimeFit = false;

    // the fields below are set mutable so they can be modified from within FitFunction
    // the values are set after the last step of the minimization
    mutable double signalUncertainty = 0;

    mutable double LDFLikelihood = 0;
    mutable double ShapeLikelihood = 0;
    mutable double StartTimeLikelihood = 0;

    // atmosphere integrated parallel to the shower axis (see GAP 2011-087)
    mutable double DX_diffusive = 0;
    // atmosphere integrated along the straight line that connects the position of the shower maximum and the station
    mutable double DX_rectilinear = 0;
    // atmosphere integrated along the straight line that connects the position of the muonic shower maximum and the station
    mutable double DX_mu_rectilinear = 0;
    mutable double height = 0;
    mutable double planeFrontTime = 0;
    mutable double distToFirstInteraction = 0;
    mutable std::vector<double> firstParticleTimes; // referred to plane front
    mutable std::vector<double> compSignals; // estimated component signals
  };


  struct Parameters {
    bool valid = false;

    double xmax = 0;
    double xmaxUnc = 0;

    double xmaxMu = 0;
    double xmaxMuUnc = 0;

    double x0 = 0;
    double x0Unc = 0;

    double hfi = 0;

    double showersize = 0;
    double showersizeUnc = 0;

    double energy = 0;
    double energyUnc = 0;

    double Nmu = 0;
    double NmuUnc = 0;

    double TimeModelOffset = 0;
    double TimeModelOffsetUnc = 0;

    utl::Vector axis;
    double axisUncTheta = 0;
    double axisUncPhi = 0;

    utl::Point core;
    double coreUncX = 0;
    double coreUncY = 0;

    utl::CoordinateSystemPtr coreCS;   // perpendicular to the ground, x-axis (E), y-axis (N)
    utl::CoordinateSystemPtr showerCS; // aligned with the shower, x-axis (early)

    utl::TimeStamp coretime;
  };


  enum fitMode {
    eFree,
    eCoupled,
    eFixed
  };


  struct minScan1D {
    std::vector<double> x;
    std::vector<double> f;
    std::vector<double> fLDF;
    std::vector<double> fShape;
    std::vector<double> fStartTime;
  };


  struct minScan2D {
    std::vector<double> x;
    std::vector<double> y;
    std::vector<double> f;  // one-dimensional to make it work with swig, has to be reshaped to two dimensions
  };


  class Fitter;


  struct FitFunction : public ROOT::Minuit2::FCNBase {
    FitFunction(const Fitter& theEf);
    double Up() const;
    double operator()(const std::vector<double>& pars) const;

    const Fitter& ef;
  };


  struct WCDFitFunction {
    WCDFitFunction(const Fitter& theEf);
    double operator()(const std::vector<double>& pars) const;

    const Fitter& ef;
  };


  class AMIGAFitFunction {
    AMIGAFitFunction(const Fitter& theEf);
    double operator()(const std::vector<double>& pars) const;
  };


  class ASCIIFitFunction {
    double operator()(const std::vector<double>& pars) const;
  };


  class Fitter {
  public:
    minScan1D scanMinimum1D(const int parIndex, const double val, const double delta, const int nSteps);
    minScan2D scanMinimum2D(const int parIndexx, const int parIndexy, const double valx, const double valy,
                            const double deltaX, const double deltaY, const int nSteps);
    void writeMinScan();

    Fitter();
    ~Fitter();
    void reset();
    void initTimeModel(const unsigned int version);
    void clearTimeModel();
    void setRecType(const unsigned int type);

    void addStationFitData(const StationFitData& sdata);
    StationFitData& getStationFitData(const int stationId);
    bool isShapeFitCandidate(const int stationId);
    bool isLDFFitCandidate(const int stationId);
    bool isStartTimeFitCandidate(const int stationId);
    bool isCandidateStation(const int stationId);

    void setupMinuit(ROOT::Minuit2::MnUserParameters& upar);
    void init();
    bool run();
    bool isFitOk(const ROOT::Minuit2::FunctionMinimum& min);
    void printResults();

    static double getAngle(const utl::Vector& v1, const utl::Vector& v2) { return Angle(v1, v2); }

    std::unique_ptr<UnivParamNS::UnivParam> fUnivParam;
    std::unique_ptr<UnivParamTimeNS::UnivParamTime> fUnivParamTime;
    std::unique_ptr<AtmosphereNS::Atmosphere> fAtmosphere;
    std::unique_ptr<UnivCalibConstantsNS::UnivCalibConstants> fUnivCalibConstants;
    // mutable: not nice, but necessary to access timeModels from within the FitFunction
    mutable std::vector<UnivTimeKG::TimeModel*> timeModels;

    utl::CoordinateSystemPtr fSiteCS;
    utl::CoordinateSystemPtr fBaryCS;
    utl::CoordinateSystemPtr fHottestStationCS;
    utl::Point fBarycenter;

    std::vector<StationFitData> fitData;
    std::vector<double> fFittedMinuitParams;
    std::unique_ptr<FitFunction> fFunc;
    std::unique_ptr<WCDFitFunction> fWCDFunc;
    //std::unique_ptr<AMIGAFitFunction> fAMIGAFunc;
    //std::unique_ptr<ASCIIFitFunction> fASCIIFunc;

    InitMode initMode;
    int atmModel;
    std::string atmType;
    bool atmUseGDAS;
    double fSiteCSHeight;
    double fDensity;
    bool eventHasSaturation;

    double fUPFudgeFactor;
    double fStartTimeBias;

    int fNSigComps;
    int fTimeModelVersion;
    int fX0XmaxModelVersion;
    int fCalibOpt;
    int fRecMixture;

    int fMinBinsAboveShapeFitLimit;
    double fMinBinContentForShapeFit;
    double fMaxBinContentForShapeFit;
    double fMinBinContentForSatShapeFit;
    double fMaxBinContentForSatShapeFit;

    double fLowerRadiusCutLDF;
    double fUpperRadiusCutLDF;
    double fLowerRadiusCutShape;
    double fUpperRadiusCutShape;
    double fLowerRadiusCutTime;
    double fUpperRadiusCutTime;

    double fXmaxLowerLimit;
    double fXmaxUpperLimit;
    double fAxisFitLimit;
    double fCoreFitLimit;

    double fMigradTol;
    int fFailLikelihood;
    int fMigradStrategy;
    double fEnergyScale;
    double fGrammageScale;
    int fVerbosityLevel;

    bool fPrintResults;
    bool fRescaleEnergy;
    bool fApplyXmaxBiasCorrection;
    bool fUseMCEnergyScale;

    bool fCalibMode;
    bool fEarlyTraceFit;
    bool fQuantileShapeFit;
    bool fIterativeFit;
    bool fOnlySignalLikelihood;
    bool fUseTimeModelOffset;
    bool fUseNmuModel;

    bool isMCEvent;

    bool doLDFFit;
    bool doShapeFit;
    bool doPartialShapeFit;
    bool doStartTimeFit;
    bool doSaturatedStartTimeFit;
    bool addMuonLDFConstraint;
    bool doMuonTimeFit;
    bool smearFirstParticlePdf;
    bool addRecovSignalConstraint;

    bool coreFixed;
    bool axisFixed;
    bool xmaxFixed;
    fitMode xmaxMuMode;
    bool NmuFixed;
    bool energyFixed;
    bool TimeModelOffsetFixed;
    fitMode x0Mode;

    Parameters mcPars;
    Parameters sdPars;
    Parameters fdPars;
    Parameters initPars;
    Parameters fittedPars;

    fitConstraint fCoreConstraint;
    fitConstraint fCoreTimeConstraint;
    fitConstraint fThetaConstraint;
    fitConstraint fPhiConstraint;
    fitConstraint fEnergyConstraint;
    fitConstraint fXmaxConstraint;
    fitConstraint fXmaxMuConstraint;
    fitConstraint fNmuConstraint;
    fitConstraint fX0Constraint;

    mutable double fLDFLikelihood;
    mutable double fShapeLikelihood;
    mutable double fStartTimeLikelihood;

    int fNStationsLDF;
    int fNStationsShape;
    int fNStationsStartTime;
  };

}


#endif