PmaTrack3D.h

Go to the documentation of this file.

#ifndef PmaTrack3D_h
#define PmaTrack3D_h

#include "TVector2.h"
#include "TVector3.h"

#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larreco/RecoAlg/PMAlg/PmaHit3D.h"
#include "larreco/RecoAlg/PMAlg/PmaNode3D.h"
#include "larreco/RecoAlg/PMAlg/Utilities.h"
namespace detinfo {
  class DetectorClocksData;
  class DetectorPropertiesData;
}

namespace pma {
  class Segment3D;
  class Track3D;
}

#include <utility>
#include <vector>

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() const noexcept { return fTag; }
  bool HasTagFlag(ETag value) const noexcept { return (fTag & value); }
  void SetTagFlag(ETag value) { fTag = (ETag)(fTag | value); }

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

  bool Initialize(detinfo::DetectorPropertiesData const& detProp, 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 detinfo::DetectorPropertiesData& detProp, 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 const* front() const { return fHits.front(); }
  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(detinfo::DetectorPropertiesData const& detProp,
               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() const;
  std::vector<unsigned int> Cryos() const;

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

  unsigned int BackTPC() const { return fNodes.back()->TPC(); }
  unsigned int BackCryo() 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(detinfo::DetectorPropertiesData const& detProp,
                                               unsigned int view,
                                               unsigned int tpc,
                                               unsigned int cryo) const;

  bool Flip(const detinfo::DetectorPropertiesData& detProp, std::vector<pma::Track3D*>& allTracks);

  void Flip();

  bool CanFlip() const;

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

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

  unsigned int TestHits(detinfo::DetectorPropertiesData const& detProp,
                        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(detinfo::DetectorPropertiesData const& detProp,
                                              unsigned int wire,
                                              unsigned int view) const;
  size_t CompleteMissingWires(detinfo::DetectorPropertiesData const& detProp, 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(const detinfo::DetectorPropertiesData& detProp,
                  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, size_t& depth);
  double GetObjFnInTree(bool skipFirst = false);
  double TuneSinglePass(bool skipFirst = false);
  double TuneFullTree(double eps = 0.001, double gmax = 50.0);

  void ApplyDriftShiftInTree(const detinfo::DetectorClocksData& clockData,
                             detinfo::DetectorPropertiesData const& detProp,
                             double dx,
                             bool skipFirst = false);
  void SetT0FromDx(const detinfo::DetectorClocksData& clockData,
                   detinfo::DetectorPropertiesData const& detProp,
                   double dx);
  double GetT0() const { return fT0; }
  bool HasT0() const noexcept { return fT0Flag; }

  void CleanupTails();

  bool ShiftEndsToHits();

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

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

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

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

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

  pma::Track3D* Split(detinfo::DetectorPropertiesData const& detProp,
                      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();
  bool GetBranches(std::vector<pma::Track3D const*>& branches, bool skipFirst = false) const;

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

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

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

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

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

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

  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();
  double AverageDist2() const;

  bool InitFromHits(detinfo::DetectorPropertiesData const& detProp,
                    int tpc,
                    int cryo,
                    float initEndSegW = 0.05F);
  bool InitFromRefPoints(detinfo::DetectorPropertiesData const& detProp, int tpc, int cryo);
  void InitFromMiddle(detinfo::DetectorPropertiesData const& detProp, 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 = nullptr);

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

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

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

  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{70};
  float fPenaltyFactor{1.0F};
  float fMaxSegStopFactor{8.0F};

  unsigned int fSegStopValue{2};
  unsigned int fMinSegStop{2};
  unsigned int fMaxSegStop{2};

  float fSegStopFactor{0.2F};
  float fPenaltyValue{0.1F};
  float fEndSegWeight{0.05F};
  float fHitsRadius{1.0F};

  double fT0{};
  bool fT0Flag{false};

  ETag fTag{kNotTagged};
};

#endif