PmaTrack3D.h

Go to the documentation of this file.

#ifndef PmaTrack3D_h
#define PmaTrack3D_h

#include "larreco/RecoAlg/PMAlg/PmaHit3D.h"
#include "larreco/RecoAlg/PMAlg/PmaNode3D.h"
#include "larreco/RecoAlg/PMAlg/PmaSegment3D.h"

namespace pma
{
        class Track3D;
}

class pma::Track3D
{
public:
        enum ETrackEnd { kBegin = -1, kEnd = 1 };
        enum EDirection { kForward = -1, kBackward = 1 };
        enum ETag
        {
            kNotTagged = 0, kTrackLike = 0, kEmLike = 1, kStopping = 2, kCosmic = 4,

        kGeometry_YY = 0x000100,
        kGeometry_YZ = 0x000200,
        kGeometry_ZZ = 0x000300,
        kGeometry_XX = 0x000400,
        kGeometry_XY = 0x000500,
        kGeometry_XZ = 0x000600,

        kGeometry_Y  = 0x001000,
        kGeometry_Z  = 0x002000,
        kGeometry_X  = 0x003000,

        kOutsideDrift_Partial   = 0x010000,
        kOutsideDrift_Complete  = 0x020000,
        kBeamIncompatible       = 0x030000
        };
        ETag GetTag(void) const { return fTag; }
        bool HasTagFlag(ETag value) const { return (fTag & value); }
        void SetTagFlag(ETag value) { fTag = (ETag)(fTag | value); }
        void SetTag(ETag value) { fTag = value; }
        

        Track3D(void);
        Track3D(const Track3D& src);
        ~Track3D(void);

        bool Initialize(float initEndSegW = 0.05F);

        pma::Hit3D* release_at(size_t index);
        void push_back(pma::Hit3D* hit) { hit->fParent = this; fHits.push_back(hit); }
        bool push_back(const art::Ptr< recob::Hit > & hit);
        bool erase(const art::Ptr< recob::Hit > & hit);

        pma::Hit3D* & operator [] (size_t index) { return fHits[index]; }
        pma::Hit3D* const & operator [] (size_t index) const { return fHits[index]; }
        pma::Hit3D* & front() { return fHits.front(); }
        pma::Hit3D* const & front() const { return fHits.front(); }
        pma::Hit3D* & back() { return fHits.back(); }
        pma::Hit3D* const & back() const { return fHits.back(); }
        size_t size() const { return fHits.size(); }

        int index_of(const pma::Hit3D* hit) const;
        int index_of(const pma::Node3D* n) const;

        double Length(size_t step = 1) const { return Length(0, size() - 1, step); }
        double Length(size_t start, size_t stop, size_t step = 1) const;

        double Dist2(const TVector2& p2d, unsigned int view, unsigned int tpc, unsigned int cryo) const;
        double Dist2(const TVector3& p3d) const;

    pma::Vector3D GetDirection3D(size_t index) const;

        void AddHits(const std::vector< art::Ptr<recob::Hit> >& hits);

        void RemoveHits(const std::vector< art::Ptr<recob::Hit> >& hits);

        unsigned int NHits(unsigned int view) const;
        unsigned int NEnabledHits(unsigned int view = geo::kUnknown) const;
        bool HasTwoViews(size_t nmin = 1) const;

        std::vector< unsigned int > TPCs(void) const;
        std::vector< unsigned int > Cryos(void) const;

        unsigned int FrontTPC(void) const { return fNodes.front()->TPC(); }
        unsigned int FrontCryo(void) const { return fNodes.front()->Cryo(); }

        unsigned int BackTPC(void) const { return fNodes.back()->TPC(); }
        unsigned int BackCryo(void) const { return fNodes.back()->Cryo(); }

        bool HasTPC(int tpc) const
        {
                for (auto n : fNodes) if (n->TPC() == tpc) return true;
                return false;
        }

        std::pair< TVector2, TVector2 > WireDriftRange(unsigned int view, unsigned int tpc, unsigned int cryo) const;

        bool Flip(std::vector< pma::Track3D* >& allTracks);

        void Flip(void);

        bool CanFlip(void) const;

        void AutoFlip(pma::Track3D::EDirection dir, double thr = 0.0, unsigned int n = 0);
        bool AutoFlip(std::vector< pma::Track3D* >& allTracks, pma::Track3D::EDirection dir, double thr = 0.0, unsigned int n = 0);

        double TestHitsMse(const std::vector< art::Ptr<recob::Hit> >& hits,
                bool normalized = true) const; // normalize to the number of hits

        unsigned int TestHits(const std::vector< art::Ptr<recob::Hit> >& hits,
                double dist = 0.4) const; // max acceptable distance [cm]

        int NextHit(int index, unsigned int view = geo::kZ, bool inclDisabled = false) const;
        int PrevHit(int index, unsigned int view = geo::kZ, bool inclDisabled = false) const;

        double HitDxByView(size_t index, unsigned int view) const;


        double GetRawdEdxSequence(std::map< size_t, std::vector<double> >& dedx, unsigned int view = geo::kZ,
                unsigned int skip = 0, bool inclDisabled = false) const;

        std::vector<float> DriftsOfWireIntersection(unsigned int wire, unsigned int view) const;
        size_t CompleteMissingWires(unsigned int view);

        void AddRefPoint(const TVector3& p) { fAssignedPoints.push_back(new TVector3(p)); }
        void AddRefPoint(double x, double y, double z) { fAssignedPoints.push_back(new TVector3(x, y, z)); }
        bool HasRefPoint(TVector3* p) const;

        double GetMse(unsigned int view = geo::kUnknown) const;

        double GetObjFunction(float penaltyFactor = 1.0F) const;

        double Optimize(int nNodes = -1, double eps = 0.01,
                bool selAllHits = true, bool setAllNodes = true,
                size_t selSegHits = 0, size_t selVtxHits = 0);

        void SortHitsInTree(bool skipFirst = false);
        void MakeProjectionInTree(bool skipFirst = false);
        bool UpdateParamsInTree(bool skipFirst = false);
        double GetObjFnInTree(bool skipFirst = false);
        double TuneSinglePass(bool skipFirst = false);
        double TuneFullTree(double eps = 0.001, double gmax = 50.0);

        void ApplyDriftShiftInTree(double dx, bool skipFirst = false);
  void SetT0FromDx(double dx);
        double GetT0(void) const { return fT0; }
  bool HasT0(void) const {return fT0Flag; }

        void CleanupTails(void);

        bool ShiftEndsToHits(void);

        std::vector< pma::Segment3D* > const & Segments(void) const { return fSegments; }

        pma::Segment3D* NextSegment(pma::Node3D* vtx) const;
        pma::Segment3D* PrevSegment(pma::Node3D* vtx) const;

        std::vector< pma::Node3D* > const & Nodes(void) const { return fNodes; }
        pma::Node3D* FirstElement(void) const { return fNodes.front(); }
        pma::Node3D* LastElement(void) const { return fNodes.back(); }

        void AddNode(pma::Node3D* node);
        void AddNode(TVector3 const & p3d, unsigned int tpc, unsigned int cryo)
        {
            double ds = fNodes.empty() ? 0 : fNodes.back()->GetDriftShift();
            AddNode(new pma::Node3D(p3d, tpc, cryo, false, ds));
        }
        bool AddNode(void);

        void InsertNode(
                TVector3 const & p3d, size_t at_idx,
                unsigned int tpc, unsigned int cryo);
        bool RemoveNode(size_t idx);

        pma::Track3D* Split(size_t idx, bool try_start_at_idx = true);

        bool AttachTo(pma::Node3D* vStart, bool noFlip = false);
        bool AttachBackTo(pma::Node3D* vStart);
        bool IsAttachedTo(pma::Track3D const * trk) const;

    void ExtendWith(pma::Track3D* src);

        pma::Track3D* GetRoot(void);
        bool GetBranches(std::vector< pma::Track3D const * >& branches, bool skipFirst = false) const;

        void MakeProjection(void);
        void UpdateProjection(void);
        void SortHits(void);

        unsigned int DisableSingleViewEnds(void);
        bool SelectHits(float fraction = 1.0F);
        bool SelectRndHits(size_t segmax, size_t vtxmax);
        bool SelectAllHits(void);

        float GetEndSegWeight(void) const { return fEndSegWeight; }
        void SetEndSegWeight(float value) { fEndSegWeight = value; }

        float GetPenalty(void) const { return fPenaltyFactor; }
        void SetPenalty(float value) { fPenaltyFactor = value; }

        unsigned int GetMaxHitsPerSeg(void) const { return fMaxHitsPerSeg; }
        void SetMaxHitsPerSeg(unsigned int value) { fMaxHitsPerSeg = value; }

private:
        void ClearNodes(void);
        void MakeFastProjection(void);

        bool AttachToSameTPC(pma::Node3D* vStart);
        bool AttachToOtherTPC(pma::Node3D* vStart);

        bool AttachBackToSameTPC(pma::Node3D* vStart);
        bool AttachBackToOtherTPC(pma::Node3D* vStart);

        void InternalFlip(std::vector< pma::Track3D* >& toSort);

        void UpdateHitsRadius(void);
        double AverageDist2(void) const;

        bool InitFromHits(int tpc, int cryo, float initEndSegW = 0.05F);
        bool InitFromRefPoints(int tpc, int cryo);
        void InitFromMiddle(int tpc, int cryo);

        pma::Track3D* GetNearestTrkInTree(const TVector3& p3d_cm, double& dist, bool skipFirst = false);
        pma::Track3D* GetNearestTrkInTree(const TVector2& p2d_cm,
                unsigned int view, unsigned int tpc, unsigned int cryo,
                double& dist, bool skipFirst = false);
        void ReassignHitsInTree(pma::Track3D* plRoot = 0);

        double HitDxByView(size_t index, unsigned int view, Track3D::EDirection dir, bool secondDir = false) const;

        bool GetUnconstrainedProj3D(art::Ptr<recob::Hit> hit, TVector3& p3d, double& dist2) const;

    void DeleteSegments(void);
        void RebuildSegments(void);
        bool SwapVertices(size_t v0, size_t v1);
        bool UpdateParams(void);

        bool CheckEndSegment(pma::Track3D::ETrackEnd endCode);

        std::vector< pma::Hit3D* > fHits;

        std::vector< TVector3* > fAssignedPoints;

        pma::Element3D* GetNearestElement(const TVector2& p2d, unsigned int view, int tpc = -1,
                bool skipFrontVtx = false, bool skipBackVtx = false) const;
        pma::Element3D* GetNearestElement(const TVector3& p3d) const;

        std::vector< pma::Node3D* > fNodes;
        std::vector< pma::Segment3D* > fSegments;

        unsigned int fMaxHitsPerSeg;
        float fPenaltyFactor;
        float fMaxSegStopFactor;

        unsigned int fSegStopValue;
        unsigned int fMinSegStop;
        unsigned int fMaxSegStop;

        float fSegStopFactor;
        float fPenaltyValue;
        float fEndSegWeight;
        float fHitsRadius;

        double fT0;
  bool fT0Flag;

        ETag fTag;
};

#endif