ProjectionMatchingAlg.h

Go to the documentation of this file.

// Class:       ProjectionMatchingAlg
// Author:      D.Stefan (Dorota.Stefan@ncbj.gov.pl) and R.Sulej (Robert.Sulej@cern.ch), May 2015
//
// Projection Matching Algorithm
// see RecoAlg/PMAlg/PmaTrack3D.h for more details.
//
//      Build 3D segments and whole tracks by matching an object 2D projections to hits, simultaneously
//      in multiple wire planes. Based on the algorithm first presented in "Precise 3D track reco..."
//      AHEP (2013) 260820, with all the tricks that we developed later and with the work for the full-event
//      topology optimization that is still under construction now (2015).
//
//      The algorithm class provides functionality to build a track from selected hits. These
//      can be detailed tracks or just simple segments (if the number of nodes to add is set to 0).
//      The parameters of optimization algorithm, fixed nodes and 3D reference points can be configured here.
//      Please, check the track making module to find a way of selecting appropriate clusteres:
//        PMAlgTrackMaker_module.cc
//

#ifndef ProjectionMatchingAlg_h
#define ProjectionMatchingAlg_h

// Framework includes
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/types/Atom.h"

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/TPCGeo.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"

#include "larreco/RecoAlg/ImagePatternAlgs/DataProvider/DataProviderAlg.h"

#include "larreco/RecoAlg/PMAlg/PmaTrack3D.h"
#include "larreco/RecoAlg/PMAlg/Utilities.h"

// ROOT & C++
#include "TH1F.h"
#include <memory>

namespace pma
{
        class ProjectionMatchingAlg;
}

class pma::ProjectionMatchingAlg
{
public:

        struct Config {
                using Name = fhicl::Name;
                using Comment = fhicl::Comment;

                fhicl::Atom<double> OptimizationEps {
                        Name("OptimizationEps"),
                        Comment("relative change of the obj.fn which stops optimization after adding a node")
                };

                fhicl::Atom<double> FineTuningEps {
                        Name("FineTuningEps"),
                        Comment("relative change of the obj.fn which stops fine-tuning of optimized track")
                };

                fhicl::Atom<double> TrkValidationDist2D {
                        Name("TrkValidationDist2D"),
                        Comment("max. distance [cm] used in the track validation in the third plane")
                };

                fhicl::Atom<double> HitTestingDist2D {
                        Name("HitTestingDist2D"),
                        Comment("max. distance [cm] used in testing compatibility of hits with the track")
                };

                fhicl::Atom<double> MinTwoViewFraction {
                        Name("MinTwoViewFraction"),
                        Comment("min. fraction of track length covered with hits from many 2D views intertwinted with each other")
                };

                fhicl::Atom<double> NodeMargin3D {
                        Name("NodeMargin3D"),
                        Comment("margin in [cm] around TPC for allowed track node positions")
                };

                fhicl::Atom<double> HitWeightU {
                        Name("HitWeightU"),
                        Comment("weights used for hits in U plane")
                };

                fhicl::Atom<double> HitWeightV {
                        Name("HitWeightV"),
                        Comment("weights used for hits in V plane")
                };

                fhicl::Atom<double> HitWeightZ {
                        Name("HitWeightZ"),
                        Comment("weights used for hits in Z plane")
                };
    };

        ProjectionMatchingAlg(const Config& config);

        ProjectionMatchingAlg(const fhicl::ParameterSet& pset) :
                ProjectionMatchingAlg(fhicl::Table<Config>(pset, {})())
        {}

        double validate_on_adc(const pma::Track3D& trk,
                const img::DataProviderAlg & adcImage, float thr) const;

        double validate_on_adc_test(const pma::Track3D& trk,
                const img::DataProviderAlg & adcImage,
                const std::vector< art::Ptr<recob::Hit> >& hits,
                TH1F * histoPassing, TH1F * histoRejected) const;

        double validate(const pma::Track3D& trk,
                const std::vector< art::Ptr<recob::Hit> >& hits) const;

        double validate(const TVector3& p0, const TVector3& p1,
                const std::vector< art::Ptr<recob::Hit> >& hits,
                unsigned int testView, unsigned int tpc, unsigned int cryo) const;

        double twoViewFraction(pma::Track3D& trk) const;

        unsigned int testHits(const pma::Track3D& trk,
                const std::vector< art::Ptr<recob::Hit> >& hits,
                double eps = 1.0) const
        { return trk.TestHits(hits, eps * fHitTestingDist2D); }

        bool isContained(const pma::Track3D& trk, float margin = 0.0F) const
        {
                return (trk.FirstElement()->SameTPC(trk.front()->Point3D(), margin) &&
                        trk.LastElement()->SameTPC(trk.back()->Point3D(), margin));
        }

        pma::Track3D* buildTrack(
                const std::vector< art::Ptr<recob::Hit> >& hits_1,
                const std::vector< art::Ptr<recob::Hit> >& hits_2 = std::vector< art::Ptr<recob::Hit> >()) const;

        pma::Track3D* buildMultiTPCTrack(const std::vector< art::Ptr<recob::Hit> >& hits) const;

        pma::Track3D* buildShowerSeg(
                const std::vector< art::Ptr<recob::Hit> >& hits, 
                const pma::Vector3D & vtx) const;

        pma::Track3D* buildSegment(
                const std::vector< art::Ptr<recob::Hit> >& hits_1,
                const std::vector< art::Ptr<recob::Hit> >& hits_2 = std::vector< art::Ptr<recob::Hit> >()) const;

        pma::Track3D* buildSegment(
                const std::vector< art::Ptr<recob::Hit> >& hits_1,
                const std::vector< art::Ptr<recob::Hit> >& hits_2,
                const TVector3& point) const;

        pma::Track3D* buildSegment(
                const std::vector< art::Ptr<recob::Hit> >& hits,
                const TVector3& point) const;

        void FilterOutSmallParts(
                double r2d,
                const std::vector< art::Ptr<recob::Hit> >& hits_in,
                std::vector< art::Ptr<recob::Hit> >& hits_out,
                const TVector2& vtx2d) const;

        void RemoveNotEnabledHits(pma::Track3D& trk) const;

        pma::Track3D* extendTrack(
                const pma::Track3D& trk,
                const std::vector< art::Ptr<recob::Hit> >& hits,
                bool add_nodes) const;

        void guideEndpoints(pma::Track3D& trk,
                const std::map< unsigned int, std::vector< art::Ptr<recob::Hit> > >& hits) const;

        void guideEndpoints(pma::Track3D& trk, pma::Track3D::ETrackEnd endpoint,
                const std::map< unsigned int, std::vector< art::Ptr<recob::Hit> > >& hits) const;

        std::vector< pma::Hit3D* > trimTrackToVolume(pma::Track3D& trk, TVector3 p0, TVector3 p1) const;

        bool alignTracks(pma::Track3D& first, pma::Track3D& second) const;

        void mergeTracks(pma::Track3D& dst, pma::Track3D& src, bool reopt) const;

        void autoFlip(pma::Track3D& trk,
                pma::Track3D::EDirection dir = Track3D::kForward,
                double thr = 0.0, unsigned int n = 0) const { trk.AutoFlip(dir, thr, n); };

        double selectInitialHits(pma::Track3D& trk, unsigned int view = geo::kZ, unsigned int* nused = 0) const;

private:

        // Helpers for guideEndpoints
        bool chkEndpointHits(int wire, int wdir, double drift_x, int view,
                unsigned int tpc, unsigned int cryo,
                const pma::Track3D& trk,
                const std::vector< art::Ptr<recob::Hit> >& hits) const;
        bool addEndpointRef(pma::Track3D& trk,
                const std::map< unsigned int, std::vector< art::Ptr<recob::Hit> > >& hits,
                std::pair<int, int> const * wires, double const * xPos,
                unsigned int tpc, unsigned int cryo) const;

        // Helpers for FilterOutSmallParts 
        bool GetCloseHits(
                double r2d, 
                const std::vector< art::Ptr<recob::Hit> >& hits_in, 
                std::vector<size_t>& used,
                std::vector< art::Ptr<recob::Hit> >& hits_out) const;

        bool Has(const std::vector<size_t>& v, size_t idx) const;

        // Make segment shorter depending on mse 
        void ShortenSeg(pma::Track3D& trk, const geo::TPCGeo& tpcgeom) const;

        // Control length of the track and number of hits which are still enabled
        bool TestTrk(pma::Track3D& trk, const geo::TPCGeo& tpcgeom) const;

        // Calculate good number of segments depending on the number of hits.
        static size_t getSegCount(size_t trk_size);


        // Parameters used in the algorithm

        double fOptimizationEps;       // relative change in the obj.function that ends optimization,
                                       // then next nodes are added or track building is finished

        double fFineTuningEps;         // relative change in the obj.function that ends final tuning

        double fTrkValidationDist2D;   // max. distance [cm] used in the track validation in the "third" plane
        double fHitTestingDist2D;      // max. distance [cm] used in testing comp. of hits with the track

        double fMinTwoViewFraction;    // min. length fraction covered with multiple 2D view hits intertwinted with each other

        // Geometry and detector properties
        geo::GeometryCore const* fGeom;
        detinfo::DetectorProperties const* fDetProp;
};

#endif