Geometric3DVertexFitter.h

Go to the documentation of this file.

#ifndef GEOMETRIC3DVERTEXFITTER_H
#define GEOMETRIC3DVERTEXFITTER_H

#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Comment.h"
#include "fhiclcpp/types/Name.h"
#include "fhiclcpp/types/Table.h"

#include "lardata/RecoObjects/TrackStatePropagator.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/TrackingTypes.h"
#include "lardataobj/RecoBase/VertexAssnMeta.h"
#include "larreco/RecoAlg/VertexWrapper.h"

namespace detinfo {
  class DetectorPropertiesData;
}

#include <functional> // std::reference_wrapper<>
#include <memory>
#include <utility> // std::pair<>
#include <vector>

namespace trkf {

  using SMatrixSym22 = recob::tracking::SMatrixSym22;
  using SVector2 = recob::tracking::SVector2;
  using SMatrixSym33 = recob::tracking::SMatrixSym33;
  using SVector3 = recob::tracking::SVector3;

  class Geometric3DVertexFitter {
  public:
    struct Config {
      using Name = fhicl::Name;
      using Comment = fhicl::Comment;
      fhicl::Atom<int> debugLevel{
        Name("debugLevel"),
        Comment("Debugging level: 0 for no printouts, 1 for minimal, 2 for full.")};
      fhicl::Atom<double> sipCut{
        Name("sipCut"),
        Comment(
          "Cut on maximum impact parameter significance to use the track in the vertex fit.")};
    };

    struct TracksFromVertexSorter {
      TracksFromVertexSorter(const recob::tracking::Point_t& vtxPos) : vtxPos_(vtxPos) {}
      bool operator()(std::reference_wrapper<const recob::Track> a,
                      std::reference_wrapper<const recob::Track> b) const
      {
        return ((a.get().Trajectory().Start() - vtxPos_).Mag2() <
                (b.get().Trajectory().Start() - vtxPos_).Mag2());
      }

    private:
      const recob::tracking::Point_t& vtxPos_;
    };

    struct ParsCovsOnPlane {
      ParsCovsOnPlane(const SVector2& p1,
                      const SVector2& p2,
                      const SMatrixSym22& c1,
                      const SMatrixSym22& c2,
                      const recob::tracking::Plane& p)
        : par1(p1), par2(p2), cov1(c1), cov2(c2), plane(p)
      {}
      SVector2 par1, par2;
      SMatrixSym22 cov1, cov2;
      recob::tracking::Plane plane;
    };

    // Constructor
    Geometric3DVertexFitter(const fhicl::Table<Config>& o,
                            const fhicl::Table<TrackStatePropagator::Config>& p)
      : debugLevel(o().debugLevel()), sipCut(o().sipCut())
    {
      prop = std::make_unique<TrackStatePropagator>(p);
    }

    VertexWrapper fitPFP(detinfo::DetectorPropertiesData const& detProp,
                         size_t iPF,
                         const art::ValidHandle<std::vector<recob::PFParticle>>& inputPFParticle,
                         const std::unique_ptr<art::FindManyP<recob::Track>>& assocTracks) const;
    VertexWrapper fitTracks(detinfo::DetectorPropertiesData const& detProp,
                            const std::vector<art::Ptr<recob::Track>>& arttracks) const;
    VertexWrapper fitTracks(detinfo::DetectorPropertiesData const& detProp,
                            TrackRefVec& tracks) const;
    VertexWrapper fitTracksWithVtx(detinfo::DetectorPropertiesData const& detProp,
                                   const std::vector<art::Ptr<recob::Track>>& tracks,
                                   const recob::tracking::Point_t& vtxPos) const;
    VertexWrapper fitTracksWithVtx(detinfo::DetectorPropertiesData const& detProp,
                                   TrackRefVec& tracks,
                                   const recob::tracking::Point_t& vtxPos) const;
    VertexWrapper closestPointAlongTrack(detinfo::DetectorPropertiesData const& detProp,
                                         const recob::Track& track,
                                         const recob::Track& other) const;
    VertexWrapper fitTwoTracks(detinfo::DetectorPropertiesData const& detProp,
                               const recob::Track& tk1,
                               const recob::Track& tk2) const;

    void addTrackToVertex(detinfo::DetectorPropertiesData const& detProp,
                          VertexWrapper& vtx,
                          const recob::Track& tk) const;

    std::vector<recob::VertexAssnMeta> computeMeta(detinfo::DetectorPropertiesData const& detProp,
                                                   const VertexWrapper& vtx);
    std::vector<recob::VertexAssnMeta> computeMeta(
      detinfo::DetectorPropertiesData const& detProp,
      const VertexWrapper& vtx,
      const std::vector<art::Ptr<recob::Track>>& arttracks);
    std::vector<recob::VertexAssnMeta> computeMeta(detinfo::DetectorPropertiesData const& detProp,
                                                   const VertexWrapper& vtx,
                                                   const TrackRefVec& trks);

    double chi2(detinfo::DetectorPropertiesData const& detProp,
                const VertexWrapper& vtx,
                const recob::Track& tk) const;
    double ip(detinfo::DetectorPropertiesData const& detProp,
              const VertexWrapper& vtx,
              const recob::Track& tk) const;
    double ipErr(detinfo::DetectorPropertiesData const& detProp,
                 const VertexWrapper& vtx,
                 const recob::Track& tk) const;
    double sip(detinfo::DetectorPropertiesData const& detProp,
               const VertexWrapper& vtx,
               const recob::Track& tk) const;
    double pDist(const VertexWrapper& vtx, const recob::Track& tk) const;

    VertexWrapper unbiasedVertex(detinfo::DetectorPropertiesData const& detProp,
                                 const VertexWrapper& vtx,
                                 const recob::Track& tk) const;
    double chi2Unbiased(detinfo::DetectorPropertiesData const& detProp,
                        const VertexWrapper& vtx,
                        const recob::Track& tk) const;
    double ipUnbiased(detinfo::DetectorPropertiesData const& detProp,
                      const VertexWrapper& vtx,
                      const recob::Track& tk) const;
    double ipErrUnbiased(detinfo::DetectorPropertiesData const& detProp,
                         const VertexWrapper& vtx,
                         const recob::Track& tk) const;
    double sipUnbiased(detinfo::DetectorPropertiesData const& detProp,
                       const VertexWrapper& vtx,
                       const recob::Track& tk) const;
    double pDistUnbiased(detinfo::DetectorPropertiesData const& detProp,
                         const VertexWrapper& vtx,
                         const recob::Track& tk) const;

  private:
    std::unique_ptr<TrackStatePropagator> prop;
    int debugLevel;
    double sipCut;

    double chi2(const ParsCovsOnPlane& pcp) const;
    double ip(const ParsCovsOnPlane& pcp) const;
    double ipErr(const ParsCovsOnPlane& pcp) const;
    double sip(const ParsCovsOnPlane& pcp) const;
    ParsCovsOnPlane getParsCovsOnPlane(detinfo::DetectorPropertiesData const& detProp,
                                       const trkf::VertexWrapper& vtx,
                                       const recob::Track& tk) const;
    std::pair<TrackState, double> weightedAverageState(ParsCovsOnPlane& pcop) const
    {
      return weightedAverageState(pcop.par1, pcop.par2, pcop.cov1, pcop.cov2, pcop.plane);
    };
    std::pair<TrackState, double> weightedAverageState(SVector2& par1,
                                                       SVector2& par2,
                                                       SMatrixSym22& cov1,
                                                       SMatrixSym22& cov2,
                                                       recob::tracking::Plane& target) const;
  };

}

#endif