TrajectoryMCSFitter.h

Go to the documentation of this file.

#ifndef TRAJECTORYMCSFITTER_H
#define TRAJECTORYMCSFITTER_H

// Framework includes
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Comment.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Sequence.h"
#include "fhiclcpp/types/Table.h"

#include "lardata/RecoObjects/TrackState.h"
#include "lardataobj/RecoBase/MCSFitResult.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/TrackTrajectory.h"
#include "lardataobj/RecoBase/Trajectory.h"

#include <array>
#include <utility>
#include <vector>

namespace trkf {
  class TrajectoryMCSFitter {
    //
  public:
    //
    struct Config {
      using Name = fhicl::Name;
      using Comment = fhicl::Comment;
      fhicl::Atom<int> pIdHypothesis{
        Name("pIdHypothesis"),
        Comment("Default particle Id Hypothesis to be used in the fit when not specified."),
        13};
      fhicl::Atom<int> minNumSegments{
        Name("minNumSegments"),
        Comment("Minimum number of segments the track is split into."),
        3};
      fhicl::Atom<double> segmentLength{
        Name("segmentLength"),
        Comment("Nominal length of track segments used in the fit."),
        14.};
      fhicl::Atom<int> minHitsPerSegment{
        Name("minHitsPerSegment"),
        Comment("Exclude segments with less hits than this value."),
        2};
      fhicl::Atom<int> nElossSteps{
        Name("nElossSteps"),
        Comment("Number of steps for computing energy loss uptream to current segment."),
        10};
      fhicl::Atom<int> eLossMode{
        Name("eLossMode"),
        Comment("Default is MPV Landau. Choose 1 for MIP (constant); 2 for Bethe-Bloch."),
        0};
      fhicl::Atom<double> pMin{Name("pMin"),
                               Comment("Minimum momentum value in likelihood scan."),
                               0.01};
      fhicl::Atom<double> pMax{Name("pMax"),
                               Comment("Maximum momentum value in likelihood scan."),
                               7.50};
      fhicl::Atom<double> pStepCoarse{
        Name("pStepCoarse"),
        Comment("Step in momentum value in initial coase likelihood scan."),
        0.01};
      fhicl::Atom<double> pStep{Name("pStep"),
                                Comment("Step in momentum value in fine grained likelihood scan."),
                                0.01};
      fhicl::Atom<double> fineScanWindow{
        Name("fineScanWindow"),
        Comment("Window size for fine grained likelihood scan around result of coarse scan."),
        0.01};
      fhicl::Sequence<double, 5> angResol{
        Name("angResol"),
        Comment(
          "Angular resolution parameters used in Highland formula. Formula is angResol[0]/(p*p) + "
          "angResol[1]/p + angResol[2] + angResol[3]*p + angResol[4]*p*p. Unit is mrad."),
        {0., 0., 3.0, 0, 0}};
      fhicl::Sequence<double, 5> hlParams{
        Name("hlParams"),
        Comment(
          "Parameters for tuning of Highland formula. Default is pdg value of 13.6. For values as "
          "in https://arxiv.org/abs/1703.0618 set to [0.1049,0.,11.0038,0,0]. Formula is "
          "hlParams[0]/(p*p) + hlParams[1]/p + hlParams[2] + hlParams[3]*p + hlParams[4]*p*p."),
        {0., 0., 13.6, 0, 0}};
      fhicl::Atom<double> segLenTolerance{
        Name("segLenTolerance"),
        Comment("Tolerance in actual segment length (lower bound)."),
        1.0};
      fhicl::Atom<bool> applySCEcorr{
        Name("applySCEcorr"),
        Comment("Flag to turn the Space Charge Effect correction on/off."),
        false};
    };
    using Parameters = fhicl::Table<Config>;
    //
    TrajectoryMCSFitter(int pIdHyp,
                        int minNSegs,
                        double segLen,
                        int minHitsPerSegment,
                        int nElossSteps,
                        int eLossMode,
                        double pMin,
                        double pMax,
                        double pStepCoarse,
                        double pStep,
                        double fineScanWindow,
                        const std::array<double, 5>& angResol,
                        const std::array<double, 5>& hlParams,
                        double segLenTolerance,
                        bool applySCEcorr)
    {
      pIdHyp_ = pIdHyp;
      minNSegs_ = minNSegs;
      segLen_ = segLen;
      minHitsPerSegment_ = minHitsPerSegment;
      nElossSteps_ = nElossSteps;
      eLossMode_ = eLossMode;
      pMin_ = pMin;
      pMax_ = pMax;
      pStepCoarse_ = pStepCoarse;
      pStep_ = pStep;
      fineScanWindow_ = fineScanWindow;
      angResol_ = angResol;
      hlParams_ = hlParams;
      segLenTolerance_ = segLenTolerance;
      applySCEcorr_ = applySCEcorr;
    }
    explicit TrajectoryMCSFitter(const Parameters& p)
      : TrajectoryMCSFitter(p().pIdHypothesis(),
                            p().minNumSegments(),
                            p().segmentLength(),
                            p().minHitsPerSegment(),
                            p().nElossSteps(),
                            p().eLossMode(),
                            p().pMin(),
                            p().pMax(),
                            p().pStepCoarse(),
                            p().pStep(),
                            p().fineScanWindow(),
                            p().angResol(),
                            p().hlParams(),
                            p().segLenTolerance(),
                            p().applySCEcorr())
    {}
    //
    recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj) const
    {
      return fitMcs(traj, pIdHyp_);
    }
    recob::MCSFitResult fitMcs(const recob::Track& track) const { return fitMcs(track, pIdHyp_); }
    recob::MCSFitResult fitMcs(const recob::Trajectory& traj) const
    {
      return fitMcs(traj, pIdHyp_);
    }
    //
    recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj, int pid) const;
    recob::MCSFitResult fitMcs(const recob::Track& track, int pid) const
    {
      return fitMcs(track.Trajectory(), pid);
    }
    recob::MCSFitResult fitMcs(const recob::Trajectory& traj, int pid) const
    {
      recob::TrackTrajectory::Flags_t flags(traj.NPoints());
      const recob::TrackTrajectory tt(traj, std::move(flags));
      return fitMcs(tt, pid);
    }
    //
    void breakTrajInSegments(const recob::TrackTrajectory& traj,
                             std::vector<size_t>& breakpoints,
                             std::vector<float>& segradlengths,
                             std::vector<float>& cumseglens) const;
    void linearRegression(const recob::TrackTrajectory& traj,
                          const size_t firstPoint,
                          const size_t lastPoint,
                          recob::tracking::Vector_t& pcdir) const;
    double mcsLikelihood(double p,
                         double theta0x,
                         std::vector<float>& dthetaij,
                         std::vector<float>& seg_nradl,
                         std::vector<float>& cumLen,
                         bool fwd,
                         int pid) const;
    //
    struct ScanResult {
    public:
      ScanResult(double ap, double apUnc, double alogL) : p(ap), pUnc(apUnc), logL(alogL) {}
      double p, pUnc, logL;
    };
    //
    const ScanResult doLikelihoodScan(std::vector<float>& dtheta,
                                      std::vector<float>& seg_nradlengths,
                                      std::vector<float>& cumLen,
                                      bool fwdFit,
                                      int pid,
                                      float detAngResol) const;
    const ScanResult doLikelihoodScan(std::vector<float>& dtheta,
                                      std::vector<float>& seg_nradlengths,
                                      std::vector<float>& cumLen,
                                      bool fwdFit,
                                      int pid,
                                      float pmin,
                                      float pmax,
                                      float pstep,
                                      float detAngResol) const;
    //
    inline double HighlandFirstTerm(const double p) const
    {
      return hlParams_[0] / (p * p) + hlParams_[1] / p + hlParams_[2] + hlParams_[3] * p +
             hlParams_[4] * p * p;
    }
    inline double DetectorAngularResolution(const double uz) const
    {
      return angResol_[0] / (uz * uz) + angResol_[1] / uz + angResol_[2] + angResol_[3] * uz +
             angResol_[4] * uz * uz;
    }
    double mass(int pid) const
    {
      if (abs(pid) == 13) { return mumass; }
      if (abs(pid) == 211) { return pimass; }
      if (abs(pid) == 321) { return kmass; }
      if (abs(pid) == 2212) { return pmass; }
      return util::kBogusD;
    }
    double energyLossBetheBloch(const double mass, const double e2) const;
    double energyLossLandau(const double mass2, const double E2, const double x) const;
    //
    double GetE(const double initial_E, const double length_travelled, const double mass) const;
    //
    int minNSegs() const { return minNSegs_; }
    double segLen() const { return segLen_; }
    double segLenTolerance() const { return segLenTolerance_; }
    //
  private:
    int pIdHyp_;
    int minNSegs_;
    double segLen_;
    int minHitsPerSegment_;
    int nElossSteps_;
    int eLossMode_;
    double pMin_;
    double pMax_;
    double pStepCoarse_;
    double pStep_;
    double fineScanWindow_;
    std::array<double, 5> angResol_;
    std::array<double, 5> hlParams_;
    double segLenTolerance_;
    bool applySCEcorr_;
  };
}

#endif