pma::Track3D Class Reference

#include <PmaTrack3D.h>

Public Types
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 }

Public Member Functions
ETag	GetTag () const noexcept
bool	HasTagFlag (ETag value) const noexcept
void	SetTagFlag (ETag 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)
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)
pma::Hit3D const *	operator[] (size_t index) const
pma::Hit3D const *	front () const
pma::Hit3D const *	back () const
size_t	size () const
int	index_of (const pma::Hit3D *hit) const
int	index_of (const pma::Node3D *n) const
double	Length (size_t step=1) const
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
	Get trajectory direction at given hit index. More...
void	AddHits (detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits)
void	RemoveHits (const std::vector< art::Ptr< recob::Hit >> &hits)
	Remove hits; removes also hit->node/seg assignments. More...
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
unsigned int	FrontCryo () const
unsigned int	BackTPC () const
unsigned int	BackCryo () const
bool	HasTPC (int tpc) const
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
	Check if the track can be flipped without breaking any other track. More...
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
	MSE of 2D hits. More...
unsigned int	TestHits (detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, double dist=0.4) const
	Count close 2D hits. More...
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)
void	AddRefPoint (double x, double y, double z)
bool	HasRefPoint (TVector3 *p) const
double	GetMse (unsigned int view=geo::kUnknown) const
	MSE of hits weighted with hit amplidudes and wire plane coefficients. More...
double	GetObjFunction (float penaltyFactor=1.0F) const
	Objective function optimized in track reconstruction. More...
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)
	Main optimization method. More...
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)
	Function to convert dx into dT0. More...
double	GetT0 () const
bool	HasT0 () const noexcept
void	CleanupTails ()
	Cut out tails with no hits assigned. More...
bool	ShiftEndsToHits ()
std::vector< pma::Segment3D * > const &	Segments () const noexcept
pma::Segment3D *	NextSegment (pma::Node3D *vtx) const
pma::Segment3D *	PrevSegment (pma::Node3D *vtx) const
std::vector< pma::Node3D * > const &	Nodes () const noexcept
pma::Node3D *	FirstElement () const
pma::Node3D *	LastElement () const
void	AddNode (pma::Node3D *node)
void	AddNode (detinfo::DetectorPropertiesData const &detProp, TVector3 const &p3d, unsigned int tpc, unsigned int cryo)
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)
	Extend the track with everything from src, delete the src;. More...
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
void	SetEndSegWeight (float value) noexcept
float	GetPenalty () const noexcept
void	SetPenalty (float value) noexcept
unsigned int	GetMaxHitsPerSeg () const noexcept
void	SetMaxHitsPerSeg (unsigned int value) noexcept

Private Member Functions
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)
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

Private Attributes
std::vector< pma::Hit3D * >	fHits
std::vector< TVector3 * >	fAssignedPoints
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}

Detailed Description

Definition at line 40 of file PmaTrack3D.h.

Member Enumeration Documentation

enum pma::Track3D::EDirection

Enumerator
kForward
kBackward

Definition at line 43 of file PmaTrack3D.h.

{ kForward = -1, kBackward = 1 };

enum pma::Track3D::ETag

Enumerator
kNotTagged
kTrackLike
kEmLike
kStopping
kCosmic
kGeometry_YY
kGeometry_YZ
kGeometry_ZZ
kGeometry_XX
kGeometry_XY
kGeometry_XZ
kGeometry_Y
kGeometry_Z
kGeometry_X
kOutsideDrift_Partial
kOutsideDrift_Complete
kBeamIncompatible

Definition at line 44 of file PmaTrack3D.h.

            {
    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
  };

enum pma::Track3D::ETrackEnd

Enumerator
kBegin
kEnd

Definition at line 42 of file PmaTrack3D.h.

{ kBegin = -1, kEnd = 1 };

Constructor & Destructor Documentation

pma::Track3D::Track3D ( )

default

pma::Track3D::Track3D

(

const Track3D &

src

)

Definition at line 31 of file PmaTrack3D.cxx.

  : fMaxHitsPerSeg(src.fMaxHitsPerSeg)
  , fPenaltyFactor(src.fPenaltyFactor)
  , fMaxSegStopFactor(src.fMaxSegStopFactor)
  , fSegStopValue(src.fSegStopValue)
  , fMinSegStop(src.fMinSegStop)
  , fMaxSegStop(src.fMaxSegStop)
  , fSegStopFactor(src.fSegStopFactor)
  , fPenaltyValue(src.fPenaltyValue)
  , fEndSegWeight(src.fEndSegWeight)
  , fHitsRadius(src.fHitsRadius)
  , fT0(src.fT0)
  , fT0Flag(src.fT0Flag)
  , fTag(src.fTag)
{
  fHits.reserve(src.fHits.size());
  for (auto const* hit : src.fHits) {
    pma::Hit3D* h3d = new pma::Hit3D(*hit);
    h3d->fParent = this;
    fHits.push_back(h3d);
  }

  fNodes.reserve(src.fNodes.size());
  for (auto const* node : src.fNodes)
    fNodes.push_back(new pma::Node3D(*node));

  for (auto const* point : src.fAssignedPoints)
    fAssignedPoints.push_back(new TVector3(*point));

  RebuildSegments();
  MakeProjection();
}

pma::Track3D::~Track3D

(

)

Definition at line 64 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < fHits.size(); i++)
    delete fHits[i];
  for (size_t i = 0; i < fAssignedPoints.size(); i++)
    delete fAssignedPoints[i];

  for (size_t i = 0; i < fSegments.size(); i++)
    delete fSegments[i];
  for (size_t i = 0; i < fNodes.size(); i++)
    if (!fNodes[i]->NextCount() && !fNodes[i]->Prev()) delete fNodes[i];
}

Member Function Documentation

void pma::Track3D::AddHits	(	detinfo::DetectorPropertiesData const &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	hits
	)

Add hits; does not update hit->node/seg assignments nor hit projection to track, so MakeProjection() and SortHits() should be called as needed.

Definition at line 404 of file PmaTrack3D.cxx.

{
  fHits.reserve(fHits.size() + hits.size());
  for (auto const& hit : hits)
    push_back(detProp, hit);
}

void pma::Track3D::AddNode

(

pma::Node3D *

node

)

Definition at line 1265 of file PmaTrack3D.cxx.

{
  fNodes.push_back(node);
  if (fNodes.size() > 1) RebuildSegments();
}

void pma::Track3D::AddNode ( detinfo::DetectorPropertiesData const &  detProp, TVector3 const &  p3d, unsigned int  tpc, unsigned int  cryo  )

inline

Definition at line 355 of file PmaTrack3D.h.

  {
    double ds = fNodes.empty() ? 0 : fNodes.back()->GetDriftShift();
    AddNode(new pma::Node3D(detProp, p3d, tpc, cryo, false, ds));
  }

bool pma::Track3D::AddNode

(

detinfo::DetectorPropertiesData const &

detProp

)

Definition at line 1272 of file PmaTrack3D.cxx.

{
  pma::Segment3D* seg;
  pma::Segment3D* maxSeg = nullptr;

  size_t si = 0;
  while (si < fSegments.size()) {
    if (!fSegments[si]->IsFrozen()) {
      maxSeg = fSegments[si];
      break;
    }
    else
      si++;
  }
  if (!maxSeg) return false;

  unsigned int nHitsByView[3];
  unsigned int nHits, maxHits = 0;
  unsigned int vIndex = 0, segHits, maxSegHits = 0;
  float segLength, maxLength = maxSeg->Length();
  for (unsigned int i = si + 1; i < fNodes.size(); i++) {
    seg = static_cast<pma::Segment3D*>(fNodes[i]->Prev());
    if (seg->IsFrozen()) continue;

    nHitsByView[0] = seg->NEnabledHits(geo::kU);
    nHitsByView[1] = seg->NEnabledHits(geo::kV);
    nHitsByView[2] = seg->NEnabledHits(geo::kZ);
    segHits = nHitsByView[0] + nHitsByView[1] + nHitsByView[2];

    if (segHits < 15) {
      if ((nHitsByView[0] == 0) && ((nHitsByView[1] < 4) || (nHitsByView[2] < 4))) continue;
      if ((nHitsByView[1] == 0) && ((nHitsByView[0] < 4) || (nHitsByView[2] < 4))) continue;
      if ((nHitsByView[2] == 0) && ((nHitsByView[0] < 4) || (nHitsByView[1] < 4))) continue;
    }

    nHits = fNodes[i]->NEnabledHits() + seg->NEnabledHits() + fNodes[i - 1]->NEnabledHits();

    if (nHits > maxHits) {
      maxHits = nHits;
      maxLength = seg->Length();
      maxSegHits = segHits;
      maxSeg = seg;
      vIndex = i;
    }
    else if (nHits == maxHits) {
      segLength = seg->Length();
      if (segLength > maxLength) {
        maxLength = segLength;
        maxSegHits = segHits;
        maxSeg = seg;
        vIndex = i;
      }
    }
  }

  if (maxSegHits > 1) {
    maxSeg->SortHits();

    nHitsByView[0] = maxSeg->NEnabledHits(geo::kU);
    nHitsByView[1] = maxSeg->NEnabledHits(geo::kV);
    nHitsByView[2] = maxSeg->NEnabledHits(geo::kZ);

    unsigned int maxViewIdx = 2, midViewIdx = 2;
    if ((nHitsByView[2] >= nHitsByView[1]) && (nHitsByView[1] >= nHitsByView[0])) {
      maxViewIdx = 2;
      midViewIdx = 1;
    }
    else if ((nHitsByView[1] >= nHitsByView[2]) && (nHitsByView[2] >= nHitsByView[0])) {
      maxViewIdx = 1;
      midViewIdx = 2;
    }
    else if ((nHitsByView[0] >= nHitsByView[2]) && (nHitsByView[2] >= nHitsByView[1])) {
      maxViewIdx = 0;
      midViewIdx = 2;
    }
    else if ((nHitsByView[2] >= nHitsByView[0]) && (nHitsByView[0] >= nHitsByView[1])) {
      maxViewIdx = 2;
      midViewIdx = 0;
    }
    else if ((nHitsByView[0] >= nHitsByView[1]) && (nHitsByView[1] >= nHitsByView[2])) {
      maxViewIdx = 0;
      midViewIdx = 1;
    }
    else if ((nHitsByView[1] >= nHitsByView[0]) && (nHitsByView[0] >= nHitsByView[2])) {
      maxViewIdx = 1;
      midViewIdx = 0;
    }
    if (nHitsByView[midViewIdx] < 2) midViewIdx = maxViewIdx;

    if (nHitsByView[midViewIdx] < 2) {
      mf::LogVerbatim("pma::Track3D") << "AddNode(): too few hits.";
      return false;
    }

    unsigned int mHits[3] = {0, 0, 0};
    unsigned int halfIndex = (nHitsByView[midViewIdx] >> 1) - 1;
    unsigned int n = 0, i = 0, i0 = 0, i1 = 0;
    while (i < maxSeg->NHits() - 1) {
      if (maxSeg->Hit(i).IsEnabled()) {
        mHits[maxSeg->Hit(i).View2D()]++;
        if (maxSeg->Hit(i).View2D() == midViewIdx) {
          if (n == halfIndex) break;
          n++;
        }
      }
      i++;
    }

    i0 = i;
    i++;
    while ((i < maxSeg->NHits()) &&
           !((maxSeg->Hit(i).View2D() == midViewIdx) && maxSeg->Hit(i).IsEnabled())) {
      i++;
    }
    i1 = i;

    if (!nHitsByView[0]) {
      if (nHitsByView[1] && (mHits[1] < 2)) {
        mf::LogVerbatim("pma::Track3D") << "AddNode(): low Ind2 hits.";
        return false;
      }
      if (nHitsByView[2] && (mHits[2] < 2)) {
        mf::LogVerbatim("pma::Track3D") << "AddNode(): low Coll hits.";
        return false;
      }
    }

    maxSeg->SetProjection(maxSeg->Hit(i0));
    maxSeg->SetProjection(maxSeg->Hit(i1));

    unsigned int tpc = maxSeg->Hit(i0).TPC();
    unsigned int cryo = maxSeg->Hit(i0).Cryo();

    pma::Node3D* p = new pma::Node3D(detProp,
                                     (maxSeg->Hit(i0).Point3D() + maxSeg->Hit(i1).Point3D()) * 0.5,
                                     tpc,
                                     cryo,
                                     false,
                                     fNodes[vIndex]->GetDriftShift());

    fNodes.insert(fNodes.begin() + vIndex, p);

    maxSeg->AddNext(fNodes[vIndex]);

    seg = new pma::Segment3D(this, fNodes[vIndex], fNodes[vIndex + 1]);
    fSegments.insert(fSegments.begin() + vIndex, seg);

    return true;
  }
  else
    return false;
}

void pma::Track3D::AddRefPoint ( const TVector3 &  p )

inline

Definition at line 265 of file PmaTrack3D.h.

  {
    fAssignedPoints.push_back(new TVector3(p));
  }

void pma::Track3D::AddRefPoint ( double  x, double  y, double  z  )

inline

Definition at line 270 of file PmaTrack3D.h.

  {
    fAssignedPoints.push_back(new TVector3(x, y, z));
  }

void pma::Track3D::ApplyDriftShiftInTree	(	const detinfo::DetectorClocksData &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		double	dx,
		bool	skipFirst = false
	)

Adjust track tree position in the drift direction (when T0 is being corrected).

Definition at line 2334 of file PmaTrack3D.cxx.

{
  pma::Node3D* node = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(node)) node = static_cast<pma::Node3D*>(segThis->Next());
  }

  while (node) {
    node->ApplyDriftShift(dx);

    auto segThis = NextSegment(node);
    for (size_t i = 0; i < node->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(node->Next(i));
      if (seg != segThis) seg->Parent()->ApplyDriftShiftInTree(clockData, detProp, dx, true);
    }

    if (segThis)
      node = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  for (auto h : fHits) {
    h->fPoint3D[0] += dx;
  }

  for (auto p : fAssignedPoints) {
    (*p)[0] += dx;
  }

  // For T0 we need to make sure we use the total shift, not just this current
  // one in case of multiple shifts
  double newdx = fNodes.front()->GetDriftShift();

  // Now convert this newdx into T0 and store in fT0
  SetT0FromDx(clockData, detProp, newdx);
}

bool pma::Track3D::AttachBackTo

(

pma::Node3D *

vStart

)

Definition at line 1656 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.back();

  if (vtx == vStart) return true; // already connected!

  pma::Track3D* rootThis = GetRoot();
  if (vStart->Prev()) {
    pma::Track3D* rootPrev = static_cast<pma::Segment3D*>(vStart->Prev())->Parent()->GetRoot();
    if (rootPrev->IsAttachedTo(rootThis)) { return false; }
  }
  else if (vStart->NextCount()) {
    pma::Track3D* rootNext = static_cast<pma::Segment3D*>(vStart->Next(0))->Parent()->GetRoot();
    if (rootNext->IsAttachedTo(rootThis)) { return false; }
  }
  else {
    return false;
  }

  for (auto n : fNodes)
    if (n == vStart) {
      mf::LogError("pma::Track3D") << "Don't create loops!";
      return false;
    }

  if (vStart->TPC() == vtx->TPC())
    return AttachBackToSameTPC(vStart);
  else
    return AttachBackToOtherTPC(vStart);
}

bool pma::Track3D::AttachBackToOtherTPC ( pma::Node3D *  vStart )

private

Definition at line 1688 of file PmaTrack3D.cxx.

{
  if (vStart->Prev()) return false;

  mf::LogVerbatim("pma::Track3D") << "Attach to track in another TPC.";

  fNodes.push_back(vStart);

  size_t idx = fNodes.size() - 1;
  fSegments.push_back(new pma::Segment3D(this, fNodes[idx - 1], fNodes[idx]));

  return true;
}

bool pma::Track3D::AttachBackToSameTPC ( pma::Node3D *  vStart )

private

Definition at line 1703 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.back();

  if (vStart->Prev()) {
    pma::Segment3D* seg = static_cast<pma::Segment3D*>(vStart->Prev());
    pma::Track3D* tp = seg->Parent();
    if (tp->NextSegment(vStart)) {
      mf::LogError("pma::Track3D") << "Cannot attach back to inner node of other track.";
      return false;
    }

    if (tp->CanFlip())
      tp->Flip(); // flip in remote vStart, no problem
    else {
      mf::LogError("pma::Track3D") << "Flip not possible, cannot attach.";
      return false;
    }
  }
  fNodes[fNodes.size() - 1] = vStart;
  fSegments[fSegments.size() - 1]->AddNext(vStart);

  while (vtx->NextCount()) // reconnect nexts to vStart
  {
    pma::Segment3D* seg = static_cast<pma::Segment3D*>(vtx->Next(0));
    pma::Track3D* trk = seg->Parent();

    vtx->RemoveNext(seg);
    trk->fNodes[0] = vStart;
    vStart->AddNext(seg);
  }

  if (vtx->NextCount() || vtx->Prev()) {
    throw cet::exception("pma::Track3D")
      << "Something is still using disconnected vertex." << std::endl;
  }
  else
    delete vtx; // ok

  return true;
}

bool pma::Track3D::AttachTo	(	pma::Node3D *	vStart,
		bool	noFlip = false
	)

Definition at line 1545 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (vtx == vStart) return true; // already connected!

  pma::Track3D* rootThis = GetRoot();
  if (vStart->Prev()) {
    pma::Track3D* rootPrev = static_cast<pma::Segment3D*>(vStart->Prev())->Parent()->GetRoot();
    if (rootPrev->IsAttachedTo(rootThis)) { return false; }
  }
  else if (vStart->NextCount()) {
    pma::Track3D* rootNext = static_cast<pma::Segment3D*>(vStart->Next(0))->Parent()->GetRoot();
    if (rootNext->IsAttachedTo(rootThis)) { return false; }
  }
  else {
    return false;
  }

  for (auto n : fNodes)
    if (n == vStart) {
      mf::LogError("pma::Track3D") << "Don't create loops!";
      return false;
    }

  if (!noFlip && CanFlip() && (vStart->TPC() == fNodes.back()->TPC()) &&
      (pma::Dist2(vtx->Point3D(), vStart->Point3D()) >
       pma::Dist2(fNodes.back()->Point3D(), vStart->Point3D())) &&
      (fNodes.back()->NextCount() == 0)) {
    mf::LogError("pma::Track3D") << "Flip, endpoint closer to vStart.";
    Flip();
  }

  if (vStart->TPC() == vtx->TPC())
    return AttachToSameTPC(vStart);
  else
    return AttachToOtherTPC(vStart);
}

bool pma::Track3D::AttachToOtherTPC ( pma::Node3D *  vStart )

private

Definition at line 1585 of file PmaTrack3D.cxx.

{
  if (fNodes.front()->Prev()) return false;

  mf::LogVerbatim("pma::Track3D") << "Attach to track in another TPC.";

  fNodes.insert(fNodes.begin(), vStart);

  fSegments.insert(fSegments.begin(), new pma::Segment3D(this, fNodes[0], fNodes[1]));

  return true;
}

bool pma::Track3D::AttachToSameTPC ( pma::Node3D *  vStart )

private

Definition at line 1599 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (vtx->Prev()) {
    pma::Segment3D* segPrev = static_cast<pma::Segment3D*>(vtx->Prev());
    pma::Track3D* tpPrev = segPrev->Parent();
    if (tpPrev->NextSegment(vtx)) { return false; }
    else if (tpPrev->CanFlip()) {
      tpPrev->Flip();
    } // flip in local vtx, no problem
    else {
      if (vStart->Prev()) {
        pma::Segment3D* segNew = static_cast<pma::Segment3D*>(vStart->Prev());
        pma::Track3D* tpNew = segNew->Parent();
        if (tpNew->NextSegment(vStart)) { return false; }
        else if (tpNew->CanFlip()) // flip in remote vStart, no problem
        {
          tpNew->Flip();

          mf::LogVerbatim("pma::Track3D") << "Reconnect prev to vStart.";
          tpPrev->fNodes[tpPrev->fNodes.size() - 1] = vStart;
          segPrev->AddNext(vStart);
        }
        else {
          mf::LogError("pma::Track3D") << "Flips in vtx and vStart not possible, cannot attach.";
          return false;
        }
      }
      else {
        mf::LogVerbatim("pma::Track3D") << "Reconnect prev to vStart.";
        tpPrev->fNodes[tpPrev->fNodes.size() - 1] = vStart;
        segPrev->AddNext(vStart);
      }
    }
  }

  while (vtx->NextCount()) // reconnect nexts to vStart
  {
    pma::Segment3D* seg = static_cast<pma::Segment3D*>(vtx->Next(0));
    pma::Track3D* trk = seg->Parent();

    vtx->RemoveNext(seg);
    trk->fNodes[0] = vStart;
    vStart->AddNext(seg);
  }

  if (vtx->NextCount() || vtx->Prev()) // better throw here
  {
    throw cet::exception("pma::Track3D") << "Something is still using disconnected vertex.";
  }
  else
    delete vtx; // ok
  return true;
}

void pma::Track3D::AutoFlip	(	pma::Track3D::EDirection	dir,
		double	thr = 0.0,
		unsigned int	n = 0
	)

Definition at line 673 of file PmaTrack3D.cxx.

{
  unsigned int nViews = 3;
  pma::dedx_map dedxInViews[3];
  for (unsigned int i = 0; i < nViews; i++) {
    GetRawdEdxSequence(dedxInViews[i], i, 1);
  }
  unsigned int bestView = 2;
  if (dedxInViews[0].size() > 2 * dedxInViews[2].size()) bestView = 0;
  if (dedxInViews[1].size() > 2 * dedxInViews[2].size()) bestView = 1;

  std::vector<std::vector<double>> dedx;
  for (size_t i = 0; i < size(); i++) {
    auto it = dedxInViews[bestView].find(i);
    if (it != dedxInViews[bestView].end()) { dedx.push_back(it->second); }
  }
  if (!dedx.empty()) dedx.pop_back();

  float dEdxStart = 0.0F, dEdxStop = 0.0F;
  float dEStart = 0.0F, dxStart = 0.0F;
  float dEStop = 0.0F, dxStop = 0.0F;
  if (dedx.size() > 4) {
    if (!n) // use default options
    {
      if (dedx.size() > 30)
        n = 12;
      else if (dedx.size() > 20)
        n = 8;
      else if (dedx.size() > 10)
        n = 4;
      else
        n = 3;
    }

    size_t k = (dedx.size() - 2) >> 1;
    if (n > k) n = k;

    for (size_t i = 1, j = 0; j < n; i++, j++) {
      dEStart += dedx[i][5];
      dxStart += dedx[i][6];
    }
    if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;

    for (size_t i = dedx.size() - 2, j = 0; j < n; i--, j++) {
      dEStop += dedx[i][5];
      dxStop += dedx[i][6];
    }
    if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
  }
  else if (dedx.size() == 4) {
    dEStart = dedx[0][5] + dedx[1][5];
    dxStart = dedx[0][6] + dedx[1][6];
    dEStop = dedx[2][5] + dedx[3][5];
    dxStop = dedx[2][6] + dedx[3][6];
    if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;
    if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
  }
  else if (dedx.size() > 1) {
    if (dedx.front()[2] > 0.0F) dEdxStart = dedx.front()[5] / dedx.front()[6];
    if (dedx.back()[2] > 0.0F) dEdxStop = dedx.back()[5] / dedx.back()[6];
  }
  else
    return;

  if ((dir == pma::Track3D::kForward) && ((1.0 + thr) * dEdxStop < dEdxStart)) {
    mf::LogVerbatim("pma::Track3D")
      << "Auto-flip fired (1), thr: " << (1.0 + thr) << ", value: " << dEdxStart / dEdxStop;
    Flip(); // particle stops at the end of the track
  }
  if ((dir == pma::Track3D::kBackward) && (dEdxStop > (1.0 + thr) * dEdxStart)) {
    mf::LogVerbatim("pma::Track3D")
      << "Auto-flip fired (2), thr: " << (1.0 + thr) << ", value: " << dEdxStop / dEdxStart;
    Flip(); // particle stops at the front of the track
  }
}

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

Definition at line 750 of file PmaTrack3D.cxx.

{
  unsigned int nViews = 3;
  pma::dedx_map dedxInViews[3];
  for (unsigned int i = 0; i < nViews; i++) {
    GetRawdEdxSequence(dedxInViews[i], i, 1);
  }
  unsigned int bestView = 2;
  if (dedxInViews[0].size() > 2 * dedxInViews[2].size()) bestView = 0;
  if (dedxInViews[1].size() > 2 * dedxInViews[2].size()) bestView = 1;

  std::vector<std::vector<double>> dedx;
  for (size_t i = 0; i < size(); i++) {
    auto it = dedxInViews[bestView].find(i);
    if (it != dedxInViews[bestView].end()) { dedx.push_back(it->second); }
  }
  if (!dedx.empty()) dedx.pop_back();

  float dEdxStart = 0.0F, dEdxStop = 0.0F;
  float dEStart = 0.0F, dxStart = 0.0F;
  float dEStop = 0.0F, dxStop = 0.0F;
  if (dedx.size() > 4) {
    if (!n) // use default options
    {
      if (dedx.size() > 30)
        n = 12;
      else if (dedx.size() > 20)
        n = 8;
      else if (dedx.size() > 10)
        n = 4;
      else
        n = 3;
    }

    size_t k = (dedx.size() - 2) >> 1;
    if (n > k) n = k;

    for (size_t i = 1, j = 0; j < n; i++, j++) {
      dEStart += dedx[i][5];
      dxStart += dedx[i][6];
    }
    if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;

    for (size_t i = dedx.size() - 2, j = 0; j < n; i--, j++) {
      dEStop += dedx[i][5];
      dxStop += dedx[i][6];
    }
    if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
  }
  else if (dedx.size() == 4) {
    dEStart = dedx[0][5] + dedx[1][5];
    dxStart = dedx[0][6] + dedx[1][6];
    dEStop = dedx[2][5] + dedx[3][5];
    dxStop = dedx[2][6] + dedx[3][6];
    if (dxStart > 0.0F) dEdxStart = dEStart / dxStart;
    if (dxStop > 0.0F) dEdxStop = dEStop / dxStop;
  }
  else if (dedx.size() > 1) {
    if (dedx.front()[2] > 0.0F) dEdxStart = dedx.front()[5] / dedx.front()[6];
    if (dedx.back()[2] > 0.0F) dEdxStop = dedx.back()[5] / dedx.back()[6];
  }
  else
    return false;

  bool done = false;
  if ((dir == pma::Track3D::kForward) && ((1.0 + thr) * dEdxStop < dEdxStart)) {
    mf::LogVerbatim("pma::Track3D")
      << "Auto-flip fired (1), thr: " << (1.0 + thr) << ", value: " << dEdxStart / dEdxStop;
    done = Flip(detProp, allTracks); // particle stops at the end of the track
  }
  if ((dir == pma::Track3D::kBackward) && (dEdxStop > (1.0 + thr) * dEdxStart)) {
    mf::LogVerbatim("pma::Track3D")
      << "Auto-flip fired (2), thr: " << (1.0 + thr) << ", value: " << dEdxStop / dEdxStart;
    done = Flip(detProp, allTracks); // particle stops at the front of the track
  }
  return done;
}

double pma::Track3D::AverageDist2 ( ) const

private

Definition at line 3160 of file PmaTrack3D.cxx.

{
  double sum = 0.0;
  unsigned int count = 0;

  for (auto n : fNodes) {
    sum += n->SumDist2();
    count += n->NEnabledHits();
  }

  for (auto s : fSegments) {
    sum += s->SumDist2();
    count += s->NEnabledHits();
  }

  if (count) { return sum / count; }
  else {
    mf::LogError("pma::Track3D") << "0 enabled hits in AverageDist2 calculation.";
    return 0;
  }
}

pma::Hit3D const* pma::Track3D::back ( ) const

inline

Definition at line 106 of file PmaTrack3D.h.

  {
    return fHits.back();
  }

unsigned int pma::Track3D::BackCryo ( ) const

inline

Definition at line 164 of file PmaTrack3D.h.

  {
    return fNodes.back()->Cryo();
  }

unsigned int pma::Track3D::BackTPC ( ) const

inline

Definition at line 159 of file PmaTrack3D.h.

  {
    return fNodes.back()->TPC();
  }

bool pma::Track3D::CanFlip

(

)

const

Check if the track can be flipped without breaking any other track.

Definition at line 654 of file PmaTrack3D.cxx.

{
  if (!fNodes.size()) { return false; }

  pma::Node3D* n = fNodes.front();
  if (n->Prev()) {
    pma::Segment3D* s = static_cast<pma::Segment3D*>(n->Prev());
    pma::Track3D* t = s->Parent();
    if (t->NextSegment(n)) { return false; } // cannot flip if starts from middle of another track
    else {
      return t->CanFlip();
    } // check root
  }
  else {
    return true;
  }
}

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

private

Definition at line 3249 of file PmaTrack3D.cxx.

{
  unsigned int v1, v2;
  switch (endCode) {
  case pma::Track3D::kBegin:
    if (fSegments.front()->IsFrozen()) return false;
    if (fNodes.front()->NextCount() > 1) return false;
    v1 = 0;
    v2 = 1;
    break;
  case pma::Track3D::kEnd:
    if (fSegments.back()->IsFrozen()) return false;
    if (fNodes.back()->NextCount() > 1) return false;
    v1 = fNodes.size() - 1;
    v2 = fNodes.size() - 2;
    break;
  default: return false;
  }
  if (fNodes[v1]->TPC() != fNodes[v2]->TPC()) return false;

  double g1, g0 = GetObjFunction();

  if (SwapVertices(v1, v2)) RebuildSegments();
  MakeProjection();
  g1 = GetObjFunction();

  if (g1 >= g0) {
    if (SwapVertices(v1, v2)) RebuildSegments();
    MakeProjection();
    return false;
  }
  else
    return true;
}

void pma::Track3D::CleanupTails

(

)

Cut out tails with no hits assigned.

Definition at line 2438 of file PmaTrack3D.cxx.

{
  unsigned int nhits = 0;
  while (!nhits && (fNodes.size() > 2) && !fNodes.front()->IsBranching()) {
    pma::Node3D* vtx = fNodes.front();
    nhits = vtx->NHits();

    pma::Segment3D* seg = NextSegment(vtx);
    nhits += seg->NHits();

    if (!nhits) {
      if (vtx->NextCount() < 2) delete vtx;
      fNodes.erase(fNodes.begin());

      delete seg;
      fSegments.erase(fSegments.begin());

      MakeProjection();
    }
  }

  nhits = 0;
  while (!nhits && (fNodes.size() > 2) && !fNodes.back()->IsBranching()) {
    pma::Node3D* vtx = fNodes.back();
    nhits = vtx->NHits();

    pma::Segment3D* seg = PrevSegment(vtx);
    nhits += seg->NHits();

    if (!nhits) {
      if (vtx->NextCount() < 1) delete fNodes.back();
      fNodes.pop_back();

      delete seg;
      fSegments.pop_back();

      MakeProjection();
    }
  }
}

void pma::Track3D::ClearNodes ( )

private

Definition at line 114 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < fNodes.size(); i++)
    delete fNodes[i];
  fNodes.clear();
}

size_t pma::Track3D::CompleteMissingWires	(	detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	view
	)

Definition at line 1193 of file PmaTrack3D.cxx.

{
  int dPrev, dw, w, wx, wPrev, i = NextHit(-1, view);

  pma::Hit3D* hitPrev = nullptr;
  pma::Hit3D* hit = nullptr;

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

  while (i < (int)size()) {
    hitPrev = hit;
    hit = fHits[i];

    if (hit->View2D() == view) {
      w = hit->Wire();
      if (hitPrev) {
        wPrev = hitPrev->Wire();
        dPrev = (int)hitPrev->PeakTime();
        if (abs(w - wPrev) > 1) {
          if (w > wPrev)
            dw = 1;
          else
            dw = -1;
          wx = wPrev + dw;
          int k = 1;
          while (wx != w) {
            int peakTime = 0;
            unsigned int iWire = wx;
            std::vector<float> drifts = DriftsOfWireIntersection(detProp, wx, view);
            if (drifts.size()) {
              peakTime = drifts[0];
              for (size_t d = 1; d < drifts.size(); d++)
                if (fabs(drifts[d] - dPrev) < fabs(peakTime - dPrev)) peakTime = drifts[d];
            }
            else {
              mf::LogVerbatim("pma::Track3D") << "Track does not intersect with the wire " << wx;
              break;
            }
            pma::Hit3D* hmiss =
              new pma::Hit3D(detProp, iWire, view, hit->TPC(), hit->Cryo(), peakTime, 1.0, 1.0);
            missHits.push_back(hmiss);
            wx += dw;
            k++;
          }
        }
      }
    }
    else
      hit = hitPrev;

    i = NextHit(i, view, true);
    while ((i < (int)size()) && (hit->TPC() != fHits[i]->TPC())) {
      hitPrev = hit;
      hit = fHits[i];
      i = NextHit(i, view, true);
    }
  }

  if (missHits.size()) {
    for (size_t hi = 0; hi < missHits.size(); hi++) {
      push_back(missHits[hi]);
    }

    MakeProjection();
    SortHits();
  }

  return missHits.size();
}

std::vector< unsigned int > pma::Track3D::Cryos

(

)

const

Definition at line 472 of file PmaTrack3D.cxx.

{
  std::vector<unsigned int> cryo_idxs;
  for (auto hit : fHits) {
    unsigned int cryo = hit->Cryo();

    bool found = false;
    for (size_t j = 0; j < cryo_idxs.size(); j++)
      if (cryo_idxs[j] == cryo) {
        found = true;
        break;
      }

    if (!found) cryo_idxs.push_back(cryo);
  }
  return cryo_idxs;
}

void pma::Track3D::DeleteSegments ( )

private

Definition at line 2419 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < fSegments.size(); i++)
    delete fSegments[i];
  fSegments.clear();
}

unsigned int pma::Track3D::DisableSingleViewEnds

(

)

Definition at line 2750 of file PmaTrack3D.cxx.

{
  SortHits();

  unsigned int nDisabled = 0;

  std::array<int, 4> hasHits{{}};

  pma::Hit3D* nextHit = nullptr;
  int hitIndex = -1;

  bool rebuild = false, stop = false;
  int nViews = 0;

  do {
    pma::Node3D* vtx = fNodes.front();
    pma::Segment3D* seg = NextSegment(vtx);
    if (!seg) break;

    if (vtx->NPoints() + seg->NPoints() > 0) hasHits[3] = 1;

    for (size_t i = 0; i < vtx->NHits(); i++) {
      hitIndex = index_of(&(vtx->Hit(i)));
      if ((hitIndex >= 0) && (hitIndex + 1 < (int)size()))
        nextHit = fHits[hitIndex + 1];
      else
        nextHit = 0;

      if (vtx->Hit(i).IsEnabled()) hasHits[vtx->Hit(i).View2D()] = 1;
      if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
      nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
      if (nViews < 2) {
        if (vtx->Hit(i).IsEnabled()) {
          vtx->Hit(i).SetEnabled(false);
          nDisabled++;
        }
      }
    }
    for (size_t i = 0; i < seg->NHits(); i++) {
      hitIndex = index_of(&(seg->Hit(i)));
      if ((hitIndex >= 0) && (hitIndex + 1 < (int)size()))
        nextHit = fHits[hitIndex + 1];
      else
        nextHit = 0;

      if (seg->Hit(i).IsEnabled()) hasHits[seg->Hit(i).View2D()] = 1;
      if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
      nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
      if (nViews < 2) {
        if (seg->Hit(i).IsEnabled()) {
          seg->Hit(i).SetEnabled(false);
          nDisabled++;
        }
      }
    }

    if (fNodes.size() < 3) break;

    nViews = hasHits[1] + hasHits[2] + hasHits[3];
    if (hasHits[0] || (nViews > 1))
      stop = true;
    else if (!fNodes.front()->IsBranching()) {
      pma::Node3D* vtx_front = fNodes.front();
      fNodes.erase(fNodes.begin());
      delete vtx_front;
      rebuild = true;
    }

  } while (!stop);

  stop = false;
  nViews = 0;
  hasHits.fill(0);

  do {
    pma::Node3D* vtx = fNodes.back();
    pma::Segment3D* seg = PrevSegment(vtx);
    if (!seg) break;

    if (vtx->NPoints() || seg->NPoints()) hasHits[3] = 1;

    for (int i = vtx->NHits() - 1; i >= 0; i--) {
      hitIndex = index_of(&(vtx->Hit(i)));
      if ((hitIndex >= 0) && (hitIndex - 1 >= 0))
        nextHit = fHits[hitIndex - 1];
      else
        nextHit = 0;

      if (vtx->Hit(i).IsEnabled()) hasHits[vtx->Hit(i).View2D()] = 1;
      if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
      nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
      if (nViews < 2) {
        if (vtx->Hit(i).IsEnabled()) {
          vtx->Hit(i).SetEnabled(false);
          nDisabled++;
        }
      }
    }
    for (int i = seg->NHits() - 1; i >= 0; i--) {
      hitIndex = index_of(&(seg->Hit(i)));
      if ((hitIndex >= 0) && (hitIndex - 1 >= 0))
        nextHit = fHits[hitIndex - 1];
      else
        nextHit = 0;

      if (seg->Hit(i).IsEnabled()) hasHits[seg->Hit(i).View2D()] = 1;
      if (nextHit && nextHit->IsEnabled()) hasHits[nextHit->View2D()] = 1;
      nViews = hasHits[0] + hasHits[1] + hasHits[2] + hasHits[3];
      if (nViews < 2) {
        if (seg->Hit(i).IsEnabled()) {
          seg->Hit(i).SetEnabled(false);
          nDisabled++;
        }
      }
    }

    if (fNodes.size() < 3) break;

    nViews = hasHits[1] + hasHits[2] + hasHits[3];
    if (hasHits[0] || (nViews > 1))
      stop = true;
    else if (!fNodes.back()->IsBranching()) {
      pma::Node3D* vtx_back = fNodes.back();
      fNodes.pop_back();
      delete vtx_back;
      rebuild = true;
    }

  } while (!stop);

  if (rebuild) {
    RebuildSegments();
    MakeProjection();
  }

  return nDisabled;
}

double pma::Track3D::Dist2	(	const TVector2 &	p2d,
		unsigned int	view,
		unsigned int	tpc,
		unsigned int	cryo
	)		const

Definition at line 2576 of file PmaTrack3D.cxx.

{
  double dist, min_dist = 1.0e12;

  int t = (int)tpc, c = (int)cryo;
  auto n0 = fNodes.front();
  if ((n0->TPC() == t) && (n0->Cryo() == c)) {
    dist = n0->GetDistance2To(p2d, view);
    if (dist < min_dist) min_dist = dist;
  }
  auto n1 = fNodes.back();
  if ((n1->TPC() == t) && (n1->Cryo() == c)) {
    dist = n1->GetDistance2To(p2d, view);
    if (dist < min_dist) min_dist = dist;
  }

  for (auto s : fSegments)
    if ((s->TPC() == t) && (s->Cryo() == c)) {
      dist = s->GetDistance2To(p2d, view);
      if (dist < min_dist) min_dist = dist;
    }
  return min_dist;
}

double pma::Track3D::Dist2

(

const TVector3 &

p3d

)

const

Definition at line 2604 of file PmaTrack3D.cxx.

{
  using namespace ranges;
  auto to_distance2 = [&p3d](auto seg) { return seg->GetDistance2To(p3d); };
  return min(fSegments | views::transform(to_distance2));
}

std::vector< float > pma::Track3D::DriftsOfWireIntersection	(	detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	wire,
		unsigned int	view
	)		const

Definition at line 1161 of file PmaTrack3D.cxx.

{
  std::vector<float> drifts;
  for (size_t i = 0; i < fNodes.size() - 1; i++) {
    int tpc = fNodes[i]->TPC(), cryo = fNodes[i]->Cryo();
    if ((tpc != fNodes[i + 1]->TPC()) || (cryo != fNodes[i + 1]->Cryo())) continue;

    TVector2 p0 = pma::CmToWireDrift(detProp,
                                     fNodes[i]->Projection2D(view).X(),
                                     fNodes[i]->Projection2D(view).Y(),
                                     view,
                                     fNodes[i]->TPC(),
                                     fNodes[i]->Cryo());
    TVector2 p1 = pma::CmToWireDrift(detProp,
                                     fNodes[i + 1]->Projection2D(view).X(),
                                     fNodes[i + 1]->Projection2D(view).Y(),
                                     view,
                                     fNodes[i + 1]->TPC(),
                                     fNodes[i + 1]->Cryo());

    if ((p0.X() - wire) * (p1.X() - wire) <= 0.0) {
      double frac_w = (wire - p0.X()) / (p1.X() - p0.X());
      double d = p0.Y() + frac_w * (p1.Y() - p0.Y());
      drifts.push_back((float)d);
    }
  }
  return drifts;
}

bool pma::Track3D::erase

(

const art::Ptr< recob::Hit > &

hit

)

Definition at line 365 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < size(); i++) {
    if (hit == fHits[i]->fHit) {
      pma::Hit3D* h3d = fHits[i];
      fHits.erase(fHits.begin() + i);
      delete h3d;
      return true;
    }
  }
  return false;
}

void pma::Track3D::ExtendWith

(

pma::Track3D *

src

)

Extend the track with everything from src, delete the src;.

Definition at line 1746 of file PmaTrack3D.cxx.

{
  if (src->fNodes.front()->Prev()) {
    throw cet::exception("pma::Track3D")
      << "Cant extend with a track being a daughter of another track." << std::endl;
  }

  src->DeleteSegments(); // disconnect from vertices and delete
  for (size_t i = 0; i < src->fNodes.size(); ++i) {
    fNodes.push_back(src->fNodes[i]);

    size_t idx = fNodes.size() - 1;
    fSegments.push_back(new pma::Segment3D(this, fNodes[idx - 1], fNodes[idx]));
  }
  src->fNodes.clear(); // just empty the node collection

  while (src->size()) {
    push_back(src->release_at(0));
  } // move hits

  for (auto p : src->fAssignedPoints) {
    fAssignedPoints.push_back(p);
  } // move 3D reference points
  src->fAssignedPoints.clear();

  MakeProjection();
  SortHits();

  delete src;
}

pma::Node3D* pma::Track3D::FirstElement ( ) const

inline

Definition at line 343 of file PmaTrack3D.h.

  {
    return fNodes.front();
  }

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

Invert the order of hits and vertices in the track, break other tracks if needed (new tracks are added to the allTracks vector). Returns true if successful or false if any of required track flips was not possible (e.g. resulting track would be composed of hits from a single 2D projection).

Definition at line 534 of file PmaTrack3D.cxx.

{
  if (fNodes.size() < 2) { return true; }

  std::vector<pma::Track3D*> toSort;

  pma::Node3D* n = fNodes.front();
  if (n->Prev()) {
    pma::Segment3D* s = static_cast<pma::Segment3D*>(n->Prev());
    pma::Track3D* t = s->Parent();

    if (t->NextSegment(n)) // starts from middle of another track: need to break that one first
    {
      int idx = t->index_of(n);
      if (idx >= 0) {
        pma::Track3D* u = t->Split(detProp, idx, false); // u is in front of t
        if (u) {
          allTracks.push_back(u);
          if (u->Flip(detProp, allTracks)) {
            InternalFlip(toSort);
            toSort.push_back(this);
          }
          else {
            mf::LogWarning("pma::Track3D")
              << "Flip(): Could not flip after split (but new track is preserved!!).";
            return false;
          }
        }
        else {
          return false;
        } // could not flip/break associated tracks, so give up on this one
      }
      else {
        throw cet::exception("pma::Track3D") << "Node not found." << std::endl;
      }
    }
    else // flip root
    {
      if (t->Flip(detProp, allTracks)) // all ok, so can flip this track
      {
        InternalFlip(toSort);
        toSort.push_back(this);
      }
      else {
        return false;
      } // could not flip/break associated tracks, so give up on this one
    }
  }
  else // simply flip
  {
    InternalFlip(toSort);
    toSort.push_back(this);
  }

  for (size_t t = 0; t < toSort.size(); t++) {
    bool sorted = false;
    for (size_t u = 0; u < t; u++)
      if (toSort[u] == toSort[t]) {
        sorted = true;
        break;
      }
    if (!sorted) {
      toSort[t]->MakeProjection();
      toSort[t]->SortHits();
    }
  }
  return true;
}

void pma::Track3D::Flip

(

)

Invert the order of hits and vertices in the track, will fail on configuration that causes breaking another track.

Definition at line 633 of file PmaTrack3D.cxx.

{
  std::vector<pma::Track3D*> toSort;
  InternalFlip(toSort);
  toSort.push_back(this);

  for (size_t t = 0; t < toSort.size(); t++) {
    bool sorted = false;
    for (size_t u = 0; u < t; u++)
      if (toSort[u] == toSort[t]) {
        sorted = true;
        break;
      }
    if (!sorted) {
      toSort[t]->MakeProjection();
      toSort[t]->SortHits();
    }
  }
}

pma::Hit3D const* pma::Track3D::front ( ) const

inline

Definition at line 101 of file PmaTrack3D.h.

  {
    return fHits.front();
  }

unsigned int pma::Track3D::FrontCryo ( ) const

inline

Definition at line 153 of file PmaTrack3D.h.

  {
    return fNodes.front()->Cryo();
  }

unsigned int pma::Track3D::FrontTPC ( ) const

inline

Definition at line 148 of file PmaTrack3D.h.

  {
    return fNodes.front()->TPC();
  }

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

Definition at line 1794 of file PmaTrack3D.cxx.

{
  for (auto trk : branches)
    if (trk == this) { throw cet::exception("pma::Track3D") << "Track tree has loop."; }

  branches.push_back(this);

  size_t i0 = 0;
  if (skipFirst) i0 = 1;

  for (size_t i = i0; i < fNodes.size(); ++i)
    for (size_t n = 0; n < fNodes[i]->NextCount(); n++) {
      pma::Segment3D* seg = static_cast<pma::Segment3D*>(fNodes[i]->Next(n));
      if (seg && (seg->Parent() != this)) seg->Parent()->GetBranches(branches, true);
    }

  return true;
}

pma::Vector3D pma::Track3D::GetDirection3D

(

size_t

index

)

const

Get trajectory direction at given hit index.

Definition at line 379 of file PmaTrack3D.cxx.

{
  pma::Hit3D* h = fHits[index];

  for (auto s : fSegments) {
    if (s->HasHit(h)) return s->GetDirection3D();
  }
  for (auto n : fNodes) {
    if (n->HasHit(h)) return n->GetDirection3D();
  }

  auto pe = GetNearestElement(h->Point2D(), h->View2D(), h->TPC());
  if (pe) {
    mf::LogWarning("pma::Track3D")
      << "GetDirection3D(): had to update hit assignment to segment/node.";
    pe->AddHit(h);
    return pe->GetDirection3D();
  }
  else {
    throw cet::exception("pma::Track3D") << "GetDirection3D(): direction of a not assigned hit "
                                         << index << " (size: " << fHits.size() << ")" << std::endl;
  }
}

float pma::Track3D::GetEndSegWeight ( ) const

inlinenoexcept

Definition at line 396 of file PmaTrack3D.h.

  {
    return fEndSegWeight;
  }

unsigned int pma::Track3D::GetMaxHitsPerSeg ( ) const

inlinenoexcept

Definition at line 418 of file PmaTrack3D.h.

  {
    return fMaxHitsPerSeg;
  }

double pma::Track3D::GetMse

(

unsigned int

view = geo::kUnknown

)

const

MSE of hits weighted with hit amplidudes and wire plane coefficients.

Definition at line 1838 of file PmaTrack3D.cxx.

{
  double sumMse = 0.0;
  unsigned int nEnabledHits = 0;
  for (auto n : fNodes) {
    sumMse += n->SumDist2(view);
    nEnabledHits += n->NEnabledHits(view);
  }
  for (auto s : fSegments) {
    sumMse += s->SumDist2(view);
    nEnabledHits += s->NEnabledHits(view);
  }

  if (nEnabledHits)
    return sumMse / nEnabledHits;
  else
    return 0.0;
}

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

private

Definition at line 2612 of file PmaTrack3D.cxx.

{
  if (fSegments.front()->TPC() < 0) skipFrontVtx = false;
  if (fSegments.back()->TPC() < 0) skipBackVtx = false;

  if (skipFrontVtx && skipBackVtx && (fSegments.size() == 1))
    return fSegments.front(); // no need for searching...

  size_t v0 = 0, v1 = fNodes.size();
  if (skipFrontVtx) v0 = 1;
  if (skipBackVtx) --v1;

  pma::Element3D* pe_min = nullptr;
  auto min_dist = std::numeric_limits<double>::max();
  for (size_t i = v0; i < v1; i++)
    if (fNodes[i]->TPC() == tpc) {
      double dist = fNodes[i]->GetDistance2To(p2d, view);
      if (dist < min_dist) {
        min_dist = dist;
        pe_min = fNodes[i];
      }
    }
  for (auto segment : fSegments)
    if (segment->TPC() == tpc) {
      if (segment->TPC() < 0) continue; // segment between TPC's

      double const dist = segment->GetDistance2To(p2d, view);
      if (dist < min_dist) {
        min_dist = dist;
        pe_min = segment;
      }
    }
  if (!pe_min) throw cet::exception("pma::Track3D") << "Nearest element not found." << std::endl;
  return pe_min;
}

pma::Element3D * pma::Track3D::GetNearestElement ( const TVector3 &  p3d ) const

private

Definition at line 2653 of file PmaTrack3D.cxx.

{
  pma::Element3D* pe_min = fNodes.front();
  double dist, min_dist = pe_min->GetDistance2To(p3d);
  for (size_t i = 1; i < fNodes.size(); i++) {
    dist = fNodes[i]->GetDistance2To(p3d);
    if (dist < min_dist) {
      min_dist = dist;
      pe_min = fNodes[i];
    }
  }
  for (size_t i = 0; i < fSegments.size(); i++) {
    dist = fSegments[i]->GetDistance2To(p3d);
    if (dist < min_dist) {
      min_dist = dist;
      pe_min = fSegments[i];
    }
  }
  return pe_min;
}

pma::Track3D * pma::Track3D::GetNearestTrkInTree ( const TVector3 &  p3d_cm, double &  dist, bool  skipFirst = false  )

private

Definition at line 2070 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());
  }

  pma::Track3D* result = this;
  dist = sqrt(Dist2(p3d_cm));

  pma::Track3D* candidate = nullptr;
  while (vtx) {
    auto segThis = NextSegment(vtx);

    double d;
    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) {
        candidate = seg->Parent()->GetNearestTrkInTree(p3d_cm, d, true);
        if (d < dist) {
          dist = d;
          result = candidate;
        }
      }
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  return result;
}

pma::Track3D* pma::Track3D::GetNearestTrkInTree ( const TVector2 &  p2d_cm, unsigned int  view, unsigned int  tpc, unsigned int  cryo, double &  dist, bool  skipFirst = false  )

private

double pma::Track3D::GetObjFnInTree

(

bool

skipFirst = false

)

Definition at line 2252 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());
  }

  double g = 0.0;
  while (vtx) {
    auto segThis = NextSegment(vtx);

    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) g += seg->Parent()->GetObjFnInTree(true);
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  return g + GetObjFunction();
}

double pma::Track3D::GetObjFunction

(

float

penaltyFactor = 1.0F

)

const

Objective function optimized in track reconstruction.

Definition at line 1858 of file PmaTrack3D.cxx.

{
  double sum = 0.0;
  float p = penaltyFactor * fPenaltyValue;
  for (size_t i = 0; i < fNodes.size(); i++) {
    sum += fNodes[i]->GetObjFunction(p, fEndSegWeight);
  }
  return sum / fNodes.size();
}

float pma::Track3D::GetPenalty ( ) const

inlinenoexcept

Definition at line 407 of file PmaTrack3D.h.

  {
    return fPenaltyFactor;
  }

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

Sequence of <hit_index, (wire, drift, X, Y, Z, dE, dx, range)> values for the track, hits tagged as outliers are skipped by default. Results are pushed into the dedx vector given in the function arguments:

  hit (segment middle if many hits) 2D projection in view:
    dedx[n][0] = wire;
    dedx[n][1] = drift;

  hit (segment middle if many hits) 3D position [cm]:
    dedx[n][2] = X;
    dedx[n][3] = Y;
    dedx[n][4] = Z;

    dedx[n][5] = dE [now ADC], energy assigned to the segment;

    dedx[n][6] = dx [cm], length of the segment.

    dedx[n][7] = range, total length to the track endpoint;

  Parameters:
    dedx  - vector to store results (empty at the begining);
    view  - view (U, V or Z) from which dedx is created;
    skip  - number of hits to skip at the begining (first hit has poorly estimated segment
            length so it can be convenient to set skip=1 and handle first hit charge manually);
    inclDisabled - if true then artificial hits added with CompleteMissingWires() are used,
                   otherwise only true hits found in ADC are used.

  Return value: sum of ADC's of hits skipped at the begining.  

Definition at line 1049 of file PmaTrack3D.cxx.

{
  dedx.clear();

  if (!size()) return 0.0;

  size_t step = 1;

  pma::Hit3D* hit = nullptr;

  double rv, dr, dR, dq, dEq, qSkipped = 0.0;

  size_t j = NextHit(-1, view, inclDisabled), s = skip;
  if (j >= size()) return 0.0F; // no charged hits at all
  while (j < size())            // look for the first hit index
  {
    hit = fHits[j];
    dq = hit->SummedADC();
    if (s) {
      qSkipped += dq;
      s--;
    }
    else
      break;

    j = NextHit(j, view, inclDisabled);
  }

  size_t jmax = PrevHit(size(), view, inclDisabled);

  std::vector<size_t> indexes;
  TVector3 p0(0., 0., 0.), p1(0., 0., 0.);
  TVector2 c0(0., 0.), c1(0., 0.);
  while (j <= jmax) {
    indexes.clear(); // prepare to collect hit indexes used for this dE/dx entry

    indexes.push_back(j);
    hit = fHits[j];

    p0 = hit->Point3D();
    p1 = hit->Point3D();

    c0.Set(hit->Wire(), hit->PeakTime());
    c1.Set(hit->Wire(), hit->PeakTime());

    dEq = hit->SummedADC(); // [now it is ADC sum]

    dr = HitDxByView(j,
                     view,
                     pma::Track3D::kForward); // protection against hits on the same position
    rv = HitDxByView(j, view);                // dx seen by j-th hit
    fHits[j]->fDx = rv;
    dR = rv;

    size_t m = 1; // number of hits with charge > 0
    while (((m < step) || (dR < 0.1) || (dr == 0.0)) && (j <= jmax)) {
      j = NextHit(j, view); // just next, even if tagged as outlier
      if (j > jmax) break;  // no more hits in this view

      hit = fHits[j];
      if (!inclDisabled && !hit->IsEnabled()) {
        if (dr == 0.0)
          continue;
        else
          break;
      }
      indexes.push_back(j);

      p1 = hit->Point3D();

      c1.Set(hit->Wire(), hit->PeakTime());

      dq = hit->SummedADC();

      dEq += dq;

      dr = HitDxByView(j, view, pma::Track3D::kForward);
      rv = HitDxByView(j, view);
      fHits[j]->fDx = rv;
      dR += rv;
      m++;
    }
    p0 += p1;
    p0 *= 0.5;
    c0 += c1;
    c0 *= 0.5;

    double range = Length(0, j);

    std::vector<double> trk_section;
    trk_section.push_back(c0.X());
    trk_section.push_back(c0.Y());
    trk_section.push_back(p0.X());
    trk_section.push_back(p0.Y());
    trk_section.push_back(p0.Z());
    trk_section.push_back(dEq);
    trk_section.push_back(dR);
    trk_section.push_back(range);

    for (auto const idx : indexes)
      dedx[idx] = trk_section;

    j = NextHit(j, view, inclDisabled);
  }

  return qSkipped;
}

pma::Track3D * pma::Track3D::GetRoot

(

)

Definition at line 1778 of file PmaTrack3D.cxx.

{
  pma::Track3D* trk = nullptr;

  if (fNodes.size()) {
    pma::Segment3D* seg = static_cast<pma::Segment3D*>(fNodes.front()->Prev());
    if (seg) trk = seg->Parent()->GetRoot();
    if (!trk) trk = this;
  }
  else
    throw cet::exception("pma::Track3D") << "Broken track, no nodes.";

  return trk;
}

double pma::Track3D::GetT0 ( ) const

inline

Definition at line 309 of file PmaTrack3D.h.

  {
    return fT0;
  }

ETag pma::Track3D::GetTag ( ) const

inlinenoexcept

Definition at line 67 of file PmaTrack3D.h.

  {
    return fTag;
  }

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

private

Calculate 3D position corresponding to 2D hit, return true if the 3D point is in the same TPC as the hit, false otherwise. Calculates also distance^2 between the hit and 2D projection of the track. NOTE: results are meaningful only if the function returns true.

Definition at line 2675 of file PmaTrack3D.cxx.

{
  TVector2 p2d = pma::WireDriftToCm(detProp,
                                    hit->WireID().Wire,
                                    hit->PeakTime(),
                                    hit->WireID().Plane,
                                    hit->WireID().TPC,
                                    hit->WireID().Cryostat);

  pma::Segment3D* seg = nullptr;
  double d2, min_d2 = 1.0e100;
  int tpc = (int)hit->WireID().TPC;
  int view = hit->WireID().Plane;
  for (size_t i = 0; i < fSegments.size(); i++) {
    if (fSegments[i]->TPC() != tpc) continue;

    d2 = fSegments[i]->GetDistance2To(p2d, view);
    if (d2 < min_d2) {
      min_d2 = d2;
      seg = fSegments[i];
    }
  }
  if (seg) {
    p3d = seg->GetUnconstrainedProj3D(p2d, view);
    dist2 = min_d2;

    pma::Node3D* prev = static_cast<pma::Node3D*>(seg->Prev());
    return prev->SameTPC(p3d); // 3D can be beyond the segment endpoints => in other TPC
  }

  return false;
}

bool pma::Track3D::HasRefPoint

(

TVector3 *

p

)

const

Definition at line 1830 of file PmaTrack3D.cxx.

{
  for (auto point : fAssignedPoints)
    if (point == p) return true;
  return false;
}

bool pma::Track3D::HasT0 ( ) const

inlinenoexcept

Check if the T0 has been set - enables us to distinguish between T0 set very close to zero or not set.

Definition at line 316 of file PmaTrack3D.h.

  {
    return fT0Flag;
  }

bool pma::Track3D::HasTagFlag ( ETag  value ) const

inlinenoexcept

Definition at line 72 of file PmaTrack3D.h.

  {
    return (fTag & value);
  }

bool pma::Track3D::HasTPC ( int  tpc ) const

inline

Definition at line 170 of file PmaTrack3D.h.

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

bool pma::Track3D::HasTwoViews

(

size_t

nmin = 1

)

const

Definition at line 443 of file PmaTrack3D.cxx.

{
  unsigned int nviews = 0;
  if (NHits(geo::kU) >= nmin) nviews++;
  if (NHits(geo::kV) >= nmin) nviews++;
  if (NHits(geo::kZ) >= nmin) nviews++;
  return nviews > 1;
}

double pma::Track3D::HitDxByView	(	size_t	index,
		unsigned int	view
	)		const

Length of the track part associated with index'th hit. Calculated as a half distance to the preceding hit plus half distance to the subsequent hit. In case of the first (last) hit - missing part is estimated as 1/4 of the distance to the next (previous) hit. NOTE: only hits from a given view are considered; other hits are accounted for segment lengths but overall dx is calculated between hits in given view.

Definition at line 1012 of file PmaTrack3D.cxx.

{
  if (index < size()) {
    return 0.5 * (HitDxByView(index, view, pma::Track3D::kForward) +
                  HitDxByView(index, view, pma::Track3D::kBackward));
  }
  else {
    mf::LogError("pma::Track3D") << "Hit index out of range.";
    return 0.0;
  }
}

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

private

Distance to the nearest subsequent (dir = Track3D::kForward) or preceeding (dir = Track3D::kBackward) hit in given view. In case of last (first) hit in this view the half-distance in opposite direction is returned. Parameter secondDir is only for internal protection - please leave the default value.

Definition at line 946 of file PmaTrack3D.cxx.

{
  pma::Hit3D* nexthit = nullptr;
  pma::Hit3D* hit = fHits[index];

  if (hit->View2D() != view) {
    mf::LogWarning("pma::Track3D") << "Function used with the hit not matching specified view.";
  }

  double dx = 0.0; // [cm]
  bool hitFound = false;
  int i = index;
  switch (dir) {
  case pma::Track3D::kForward:
    while (!hitFound && (++i < (int)size())) {
      nexthit = fHits[i];
      dx += sqrt(pma::Dist2(hit->Point3D(), nexthit->Point3D()));

      if (nexthit->View2D() == view)
        hitFound = true;
      else
        hitFound = false;

      hit = nexthit;
    }
    if (!hitFound) {
      if (!secondDir)
        dx = 0.5 * HitDxByView(index, view, pma::Track3D::kBackward, true);
      else {
        dx = Length();
        mf::LogWarning("pma::Track3D") << "Single hit in this view.";
      }
    }
    break;

  case pma::Track3D::kBackward:
    while (!hitFound && (--i >= 0)) {
      nexthit = fHits[i];
      dx += sqrt(pma::Dist2(hit->Point3D(), nexthit->Point3D()));

      if (nexthit->View2D() == view)
        hitFound = true;
      else
        hitFound = false;

      hit = nexthit;
    }
    if (!hitFound) {
      if (!secondDir)
        dx = 0.5 * HitDxByView(index, view, pma::Track3D::kForward, true);
      else {
        dx = Length();
        mf::LogWarning("pma::Track3D") << "Single hit in this view.";
      }
    }
    break;

  default: mf::LogError("pma::Track3D") << "Direction undefined."; break;
  }
  return dx;
}

int pma::Track3D::index_of

(

const pma::Hit3D *

hit

)

const

Definition at line 336 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < size(); i++)
    if (fHits[i] == hit) return (int)i;
  return -1;
}

int pma::Track3D::index_of

(

const pma::Node3D *

n

)

const

Definition at line 328 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < fNodes.size(); ++i)
    if (fNodes[i] == n) return (int)i;
  return -1;
}

bool pma::Track3D::InitFromHits ( detinfo::DetectorPropertiesData const &  detProp, int  tpc, int  cryo, float  initEndSegW = 0.05F  )

private

Definition at line 122 of file PmaTrack3D.cxx.

{
  art::ServiceHandle<geo::Geometry const> geom;

  float wtmp = fEndSegWeight;
  fEndSegWeight = initEndSegW;

  // endpoints for the first combination:
  TVector3 v3d_1(0., 0., 0.), v3d_2(0., 0., 0.);

  assert(!fHits.empty());

  pma::Hit3D* hit0_a = fHits.front();
  pma::Hit3D* hit1_a = fHits.front();

  float minX = detProp.ConvertTicksToX(hit0_a->PeakTime(), hit0_a->View2D(), tpc, cryo);
  float maxX = detProp.ConvertTicksToX(hit1_a->PeakTime(), hit1_a->View2D(), tpc, cryo);
  for (auto hit : fHits) {
    double const x = detProp.ConvertTicksToX(hit->PeakTime(), hit->View2D(), tpc, cryo);
    if (x < minX) {
      minX = x;
      hit0_a = hit;
    }
    if (x > maxX) {
      maxX = x;
      hit1_a = hit;
    }
  }

  pma::Hit3D* hit0_b = nullptr;
  pma::Hit3D* hit1_b = nullptr;

  float minDiff0 = 5000, minDiff1 = 5000;
  for (auto hit : fHits) {
    double const x = detProp.ConvertTicksToX(hit->PeakTime(), hit->View2D(), tpc, cryo);
    double diff = fabs(x - minX);
    if ((diff < minDiff0) && (hit->View2D() != hit0_a->View2D())) {
      minDiff0 = diff;
      hit0_b = hit;
    }
    diff = fabs(x - maxX);
    if ((diff < minDiff1) && (hit->View2D() != hit1_a->View2D())) {
      minDiff1 = diff;
      hit1_b = hit;
    }
  }

  if (hit0_a && hit0_b && hit1_a && hit1_b) {
    double x = 0.5 * (detProp.ConvertTicksToX(hit0_a->PeakTime(), hit0_a->View2D(), tpc, cryo) +
                      detProp.ConvertTicksToX(hit0_b->PeakTime(), hit0_b->View2D(), tpc, cryo));
    double y, z;
    geom->IntersectionPoint(
      hit0_a->Wire(), hit0_b->Wire(), hit0_a->View2D(), hit0_b->View2D(), cryo, tpc, y, z);
    v3d_1.SetXYZ(x, y, z);

    x = 0.5 * (detProp.ConvertTicksToX(hit1_a->PeakTime(), hit1_a->View2D(), tpc, cryo) +
               detProp.ConvertTicksToX(hit1_b->PeakTime(), hit1_b->View2D(), tpc, cryo));
    geom->IntersectionPoint(
      hit1_a->Wire(), hit1_b->Wire(), hit1_a->View2D(), hit1_b->View2D(), cryo, tpc, y, z);
    v3d_2.SetXYZ(x, y, z);

    ClearNodes();
    AddNode(detProp, v3d_1, tpc, cryo);
    AddNode(detProp, v3d_2, tpc, cryo);

    MakeProjection();
    UpdateHitsRadius();
    Optimize(detProp, 0, 0.01F, false, true, 100);
    SelectAllHits();
  }
  else {
    mf::LogVerbatim("pma::Track3D") << "Good hits not found.";
    fEndSegWeight = wtmp;
    return false;
  }

  if (pma::Dist2(fNodes.front()->Point3D(), fNodes.back()->Point3D()) < (0.3 * 0.3)) {
    mf::LogVerbatim("pma::Track3D") << "Short initial segment.";
    fEndSegWeight = wtmp;
    return false;
  }

  fEndSegWeight = wtmp;
  return true;
}

void pma::Track3D::InitFromMiddle ( detinfo::DetectorPropertiesData const &  detProp, int  tpc, int  cryo  )

private

Definition at line 300 of file PmaTrack3D.cxx.

{
  art::ServiceHandle<geo::Geometry const> geom;

  const auto& tpcGeo = geom->TPC(tpc, cryo);

  double minX = tpcGeo.MinX(), maxX = tpcGeo.MaxX();
  double minY = tpcGeo.MinY(), maxY = tpcGeo.MaxY();
  double minZ = tpcGeo.MinZ(), maxZ = tpcGeo.MaxZ();

  TVector3 v3d_1(0.5 * (minX + maxX), 0.5 * (minY + maxY), 0.5 * (minZ + maxZ));
  TVector3 v3d_2(v3d_1);

  TVector3 shift(5.0, 5.0, 5.0);
  v3d_1 += shift;
  v3d_2 -= shift;

  ClearNodes();
  AddNode(detProp, v3d_1, tpc, cryo);
  AddNode(detProp, v3d_2, tpc, cryo);

  MakeProjection();
  UpdateHitsRadius();

  Optimize(detProp, 0, 0.01F);
}

bool pma::Track3D::InitFromRefPoints ( detinfo::DetectorPropertiesData const &  detProp, int  tpc, int  cryo  )

private

Definition at line 212 of file PmaTrack3D.cxx.

{
  if (fAssignedPoints.size() < 2) return false;

  ClearNodes();

  TVector3 mean(0., 0., 0.), stdev(0., 0., 0.), p(0., 0., 0.);
  for (size_t i = 0; i < fAssignedPoints.size(); i++) {
    p = *(fAssignedPoints[i]);
    mean += p;
    p.SetXYZ(p.X() * p.X(), p.Y() * p.Y(), p.Z() * p.Z());
    stdev += p;
  }
  stdev *= 1.0 / fAssignedPoints.size();
  mean *= 1.0 / fAssignedPoints.size();
  p = mean;
  p.SetXYZ(p.X() * p.X(), p.Y() * p.Y(), p.Z() * p.Z());
  stdev -= p;

  double sx = stdev.X(), sy = stdev.Y(), sz = stdev.Z();
  if (sx >= 0.0)
    sx = sqrt(sx);
  else
    sx = 0.0;
  if (sy >= 0.0)
    sy = sqrt(sy);
  else
    sy = 0.0;
  if (sz >= 0.0)
    sz = sqrt(sz);
  else
    sz = 0.0;
  stdev.SetXYZ(sx, sy, sz);

  double scale = 2.0 * stdev.Mag();
  double iscale = 1.0 / scale;

  size_t max_index = 0;
  double norm2, max_norm2 = 0.0;
  std::vector<TVector3> data;
  for (size_t i = 0; i < fAssignedPoints.size(); i++) {
    p = *(fAssignedPoints[i]);
    p -= mean;
    p *= iscale;
    norm2 = p.Mag2();
    if (norm2 > max_norm2) {
      max_norm2 = norm2;
      max_index = i;
    }
    data.push_back(p);
  }

  double y = 0.0, kappa = 1.0, prev_kappa, kchg = 1.0;
  TVector3 w(data[max_index]);

  while (kchg > 0.0001)
    for (size_t i = 0; i < data.size(); i++) {
      y = (data[i] * w);
      w += (y / kappa) * (data[i] - y * w);

      prev_kappa = kappa;
      kappa += y * y;
      kchg = fabs((kappa - prev_kappa) / prev_kappa);
    }
  w *= 1.0 / w.Mag();

  TVector3 v1(w), v2(w);
  v1 *= scale;
  v1 += mean;
  v2 *= -scale;
  v2 += mean;
  std::sort(fAssignedPoints.begin(), fAssignedPoints.end(), pma::bSegmentProjLess(v1, v2));
  for (size_t i = 0; i < fAssignedPoints.size(); i++) {
    AddNode(detProp, *(fAssignedPoints[i]), tpc, cryo);
  }

  RebuildSegments();
  MakeProjection();

  if (size()) UpdateHitsRadius();

  Optimize(detProp, 0, 0.01F, false, true, 100);
  SelectAllHits();

  return true;
}

bool pma::Track3D::Initialize	(	detinfo::DetectorPropertiesData const &	detProp,
		float	initEndSegW = 0.05F
	)

Definition at line 78 of file PmaTrack3D.cxx.

{
  if (!HasTwoViews(2)) {
    mf::LogError("pma::Track3D") << "Need min. 2 hits per view, at least two views.";
    return false;
  }

  auto cryos = Cryos();
  if (cryos.size() > 1) {
    mf::LogError("pma::Track3D") << "Only one cryostat for now, please.";
    return false;
  }
  int cryo = cryos.front();

  auto tpcs = TPCs();
  if (tpcs.size() > 1) {
    mf::LogError("pma::Track3D") << "Only one TPC, please.";
    return false;
  }
  // single tpc, many tpc's are ok, but need to be handled from ProjectionMatchingAlg::buildMultiTPCTrack()
  int tpc = tpcs.front();

  if (InitFromRefPoints(detProp, tpc, cryo))
    mf::LogVerbatim("pma::Track3D") << "Track initialized with 3D reference points.";
  else if (InitFromHits(detProp, tpc, cryo, initEndSegW))
    mf::LogVerbatim("pma::Track3D") << "Track initialized with hit positions.";
  else {
    InitFromMiddle(detProp, tpc, cryo);
    mf::LogVerbatim("pma::Track3D") << "Track initialized in the module center.";
  }

  UpdateHitsRadius();
  return true;
}

void pma::Track3D::InsertNode	(	detinfo::DetectorPropertiesData const &	detProp,
		TVector3 const &	p3d,
		size_t	at_idx,
		unsigned int	tpc,
		unsigned int	cryo
	)

Definition at line 1426 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx =
    new pma::Node3D(detProp, p3d, tpc, cryo, false, fNodes[at_idx]->GetDriftShift());
  fNodes.insert(fNodes.begin() + at_idx, vtx);

  if (fNodes.size() > 1) RebuildSegments();
}

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

private

Definition at line 605 of file PmaTrack3D.cxx.

{
  for (size_t i = 0; i < fNodes.size() - 1; i++)
    if (fNodes[i]->NextCount() > 1) {
      for (size_t j = 0; j < fNodes[i]->NextCount(); j++) {
        pma::Segment3D* s = static_cast<pma::Segment3D*>(fNodes[i]->Next(j));
        if (s->Parent() != this) toSort.push_back(s->Parent());
      }
    }

  if (fNodes.back()->NextCount())
    for (size_t j = 0; j < fNodes.back()->NextCount(); j++) {
      pma::Segment3D* s = static_cast<pma::Segment3D*>(fNodes.back()->Next(j));
      toSort.push_back(s->Parent());
    }

  if (fNodes.front()->Prev()) {
    pma::Segment3D* s = static_cast<pma::Segment3D*>(fNodes.front()->Prev());
    toSort.push_back(s->Parent());
    s->Parent()->InternalFlip(toSort);
  }

  std::reverse(fNodes.begin(), fNodes.end());
  toSort.push_back(this);
  RebuildSegments();
}

bool pma::Track3D::IsAttachedTo

(

pma::Track3D const *

trk

)

const

Definition at line 1814 of file PmaTrack3D.cxx.

{
  if (trk == this) return true;

  std::vector<pma::Track3D const*> branchesThis, branchesTrk;

  if (!trk->GetBranches(branchesTrk)) return true; // has loop - better check the reason
  if (!GetBranches(branchesThis)) return true;     // has loop - better check the reason

  for (auto bThis : branchesThis)
    for (auto bTrk : branchesTrk)
      if (bThis == bTrk) return true;
  return false;
}

pma::Node3D* pma::Track3D::LastElement ( ) const

inline

Definition at line 348 of file PmaTrack3D.h.

  {
    return fNodes.back();
  }

double pma::Track3D::Length ( size_t  step = 1 ) const

inline

Definition at line 120 of file PmaTrack3D.h.

  {
    return Length(0, size() - 1, step);
  }

double pma::Track3D::Length	(	size_t	start,
		size_t	stop,
		size_t	step = 1
	)		const

Definition at line 878 of file PmaTrack3D.cxx.

{
  auto const nHits = size();
  if (nHits < 2) return 0.0;

  if (start > stop) {
    size_t tmp = stop;
    stop = start;
    start = tmp;
  }
  if (start >= nHits - 1) return 0.0;
  if (stop >= nHits) stop = nHits - 1;

  double result = 0.0;
  pma::Hit3D *h1 = nullptr, *h0 = fHits[start];

  if (!step) step = 1;
  size_t i = start + step;
  while (i <= stop) {
    h1 = fHits[i];
    result += sqrt(pma::Dist2(h1->Point3D(), h0->Point3D()));
    h0 = h1;
    i += step;
  }
  if (i - step < stop) // last step jumped beyond the stop point
  {                    // so need to add the last short segment
    result += sqrt(pma::Dist2(h0->Point3D(), back()->Point3D()));
  }
  return result;
}

void pma::Track3D::MakeFastProjection ( )

private

Definition at line 2990 of file PmaTrack3D.cxx.

{
  std::vector<std::pair<pma::Hit3D*, pma::Element3D*>> assignments;
  assignments.reserve(fHits.size());

  for (auto hi : fHits) {
    pma::Element3D* pe = nullptr;

    for (auto s : fSegments) {
      for (size_t j = 0; j < s->NHits(); ++j)
        if (hi == s->Hits()[j]) // look at next/prev vtx,seg,vtx
        {
          pe = s;
          double min_d2 = s->GetDistance2To(hi->Point2D(), hi->View2D());
          int const tpc = hi->TPC();

          pma::Node3D* nnext = static_cast<pma::Node3D*>(s->Next());
          if (nnext->TPC() == tpc) {
            double const d2 = nnext->GetDistance2To(hi->Point2D(), hi->View2D());
            if (d2 < min_d2) {
              min_d2 = d2;
              pe = nnext;
            }

            pma::Segment3D* snext = NextSegment(nnext);
            if (snext && (snext->TPC() == tpc)) {
              double const d2 = snext->GetDistance2To(hi->Point2D(), hi->View2D());
              if (d2 < min_d2) {
                min_d2 = d2;
                pe = snext;
              }

              nnext = static_cast<pma::Node3D*>(snext->Next());
              if (nnext->TPC() == tpc) {
                double const d2 = nnext->GetDistance2To(hi->Point2D(), hi->View2D());
                if (d2 < min_d2) {
                  min_d2 = d2;
                  pe = nnext;
                }
              }
            }
          }

          pma::Node3D* nprev = static_cast<pma::Node3D*>(s->Prev());
          if (nprev->TPC() == tpc) {
            double const d2 = nprev->GetDistance2To(hi->Point2D(), hi->View2D());
            if (d2 < min_d2) {
              min_d2 = d2;
              pe = nprev;
            }

            pma::Segment3D* sprev = PrevSegment(nprev);
            if (sprev && (sprev->TPC() == tpc)) {
              double const d2 = sprev->GetDistance2To(hi->Point2D(), hi->View2D());
              if (d2 < min_d2) {
                min_d2 = d2;
                pe = sprev;
              }

              nprev = static_cast<pma::Node3D*>(sprev->Prev());
              if (nprev->TPC() == tpc) {
                double const d2 = nprev->GetDistance2To(hi->Point2D(), hi->View2D());
                if (d2 < min_d2) {
                  min_d2 = d2;
                  pe = nprev;
                }
              }
            }
          }

          s->RemoveHitAt(j);
          break;
        }
      if (pe) break;
    }

    if (!pe)
      for (auto n : fNodes) {
        for (size_t j = 0; j < n->NHits(); ++j)
          if (hi == n->Hits()[j]) // look at next/prev seg,vtx,seg
          {
            pe = n;
            double d2, min_d2 = n->GetDistance2To(hi->Point2D(), hi->View2D());
            int tpc = hi->TPC();

            pma::Segment3D* snext = NextSegment(n);
            if (snext && (snext->TPC() == tpc)) {
              d2 = snext->GetDistance2To(hi->Point2D(), hi->View2D());
              if (d2 < min_d2) {
                min_d2 = d2;
                pe = snext;
              }

              pma::Node3D* nnext = static_cast<pma::Node3D*>(snext->Next());
              if (nnext->TPC() == tpc) {
                d2 = nnext->GetDistance2To(hi->Point2D(), hi->View2D());
                if (d2 < min_d2) {
                  min_d2 = d2;
                  pe = nnext;
                }

                snext = NextSegment(nnext);
                if (snext && (snext->TPC() == tpc)) {
                  d2 = snext->GetDistance2To(hi->Point2D(), hi->View2D());
                  if (d2 < min_d2) {
                    min_d2 = d2;
                    pe = snext;
                  }
                }
              }
            }

            pma::Segment3D* sprev = PrevSegment(n);
            if (sprev && (sprev->TPC() == tpc)) {
              d2 = sprev->GetDistance2To(hi->Point2D(), hi->View2D());
              if (d2 < min_d2) {
                min_d2 = d2;
                pe = sprev;
              }

              pma::Node3D* nprev = static_cast<pma::Node3D*>(sprev->Prev());
              if (nprev->TPC() == tpc) {
                d2 = nprev->GetDistance2To(hi->Point2D(), hi->View2D());
                if (d2 < min_d2) {
                  min_d2 = d2;
                  pe = nprev;
                }

                sprev = PrevSegment(nprev);
                if (sprev && (sprev->TPC() == tpc)) {
                  d2 = sprev->GetDistance2To(hi->Point2D(), hi->View2D());
                  if (d2 < min_d2) {
                    min_d2 = d2;
                    pe = sprev;
                  }
                }
              }
            }

            n->RemoveHitAt(j);
            break;
          }
        if (pe) break;
      }

    if (pe)
      assignments.emplace_back(hi, pe);
    else
      mf::LogWarning("pma::Track3D") << "Hit was not assigned to any element.";
  }

  for (auto const& a : assignments)
    a.second->AddHit(a.first);

  for (auto n : fNodes)
    n->UpdateHitParams();
  for (auto s : fSegments)
    s->UpdateHitParams();
}

void pma::Track3D::MakeProjection

(

)

Definition at line 2953 of file PmaTrack3D.cxx.

{
  for (auto n : fNodes)
    n->ClearAssigned(this);
  for (auto s : fSegments)
    s->ClearAssigned(this);

  pma::Element3D* pe = nullptr;

  bool skipFrontVtx = false, skipBackVtx = false;

  // skip outermost vertices if not branching
  if (!(fNodes.front()->IsFrozen()) && !(fNodes.front()->Prev()) &&
      (fNodes.front()->NextCount() == 1)) {
    skipFrontVtx = true;
  }
  if (!(fNodes.front()->IsFrozen()) && (fNodes.back()->NextCount() == 0)) { skipBackVtx = true; }

  for (auto h : fHits) // assign hits to nodes/segments
  {
    pe = GetNearestElement(h->Point2D(), h->View2D(), h->TPC(), skipFrontVtx, skipBackVtx);
    pe->AddHit(h);
  }

  for (auto p : fAssignedPoints) // assign ref points to nodes/segments
  {
    pe = GetNearestElement(*p);
    pe->AddPoint(p);
  }

  for (auto n : fNodes)
    n->UpdateHitParams();
  for (auto s : fSegments)
    s->UpdateHitParams();
}

void pma::Track3D::MakeProjectionInTree

(

bool

skipFirst = false

)

Definition at line 2200 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());
  }

  while (vtx) {
    auto segThis = NextSegment(vtx);

    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) seg->Parent()->MakeProjectionInTree(true);
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  MakeProjection();
}

unsigned int pma::Track3D::NEnabledHits

(

unsigned int

view = geo::kUnknown

)

const

Definition at line 433 of file PmaTrack3D.cxx.

{
  unsigned int n = 0;
  for (auto hit : fHits) {
    if (hit->IsEnabled() && ((view == geo::kUnknown) || (view == hit->View2D()))) n++;
  }
  return n;
}

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

Definition at line 910 of file PmaTrack3D.cxx.

{
  pma::Hit3D* hit = nullptr;
  if (index < -1) index = -1;
  while (++index < (int)size()) // look for the next index of hit from the view
  {
    hit = fHits[index];
    if (hit->View2D() == view) {
      if (inclDisabled)
        break;
      else if (hit->IsEnabled())
        break;
    }
  }
  return index;
}

pma::Segment3D * pma::Track3D::NextSegment

(

pma::Node3D *

vtx

)

const

Definition at line 1025 of file PmaTrack3D.cxx.

{
  pma::Segment3D* seg = nullptr;
  unsigned int nCount = vtx->NextCount();
  unsigned int k = 0;
  while (k < nCount) {
    seg = static_cast<pma::Segment3D*>(vtx->Next(k));
    if (seg && (seg->Parent() == this)) return seg;
    k++;
  }
  return 0;
}

unsigned int pma::Track3D::NHits

(

unsigned int

view

)

const

Definition at line 423 of file PmaTrack3D.cxx.

{
  unsigned int n = 0;
  for (auto hit : fHits) {
    if (hit->View2D() == view) n++;
  }
  return n;
}

std::vector<pma::Node3D*> const& pma::Track3D::Nodes ( ) const

inlinenoexcept

Definition at line 338 of file PmaTrack3D.h.

  {
    return fNodes;
  }

pma::Hit3D* pma::Track3D::operator[] ( size_t  index )

inline

Definition at line 98 of file PmaTrack3D.h.

{ return fHits[index]; }

pma::Hit3D const* pma::Track3D::operator[] ( size_t  index ) const

inline

Definition at line 99 of file PmaTrack3D.h.

{ return fHits[index]; }

double pma::Track3D::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
	)

Main optimization method.

Definition at line 1869 of file PmaTrack3D.cxx.

{
  if (!fNodes.size()) {
    mf::LogError("pma::Track3D") << "Track not initialized.";
    return 0.0;
  }

  if (!UpdateParams()) {
    mf::LogError("pma::Track3D") << "Track empty.";
    return 1.0e10;
  }
  double g0 = GetObjFunction(), g1 = 0.0;
  if (g0 == 0.0) return g0;

  // set branching flag only at the beginning, optimization is not changing
  // that and new nodes are not branching
  for (auto n : fNodes)
    n->SetVertexToBranching(setAllNodes);

  bool stop = false;
  fMinSegStop = fSegments.size();
  fMaxSegStop = (int)(size() / fMaxSegStopFactor) + 1;
  do {
    if (selSegHits || selVtxHits) SelectRndHits(selSegHits, selVtxHits);

    bool stepDone = true;
    unsigned int stepIter = 0;
    do {
      double gstep = 1.0;
      unsigned int iter = 0;
      while ((gstep > eps) && (iter < 1000)) {
        if ((fNodes.size() < 4) || (iter % 10 == 0))
          MakeProjection();
        else
          MakeFastProjection();

        if (!UpdateParams()) {
          mf::LogError("pma::Track3D") << "Track empty.";
          return 0.0;
        }

        for (auto n : fNodes)
          n->Optimize(fPenaltyValue, fEndSegWeight);

        g1 = g0;
        g0 = GetObjFunction();

        if (g0 == 0.0F) {
          MakeProjection();
          break;
        }
        gstep = fabs(g0 - g1) / g0;
        iter++;
      }

      stepIter++;
      if (fNodes.size() > 2) {
        stepDone = !(CheckEndSegment(pma::Track3D::kEnd) || CheckEndSegment(pma::Track3D::kBegin));
      }
    } while (!stepDone && (stepIter < 5));

    if (selAllHits) {
      selAllHits = false;
      selSegHits = 0;
      selVtxHits = 0;
      if (SelectAllHits()) continue;
    }

    switch (nNodes) {
    case 0: stop = true; break; // just optimize existing vertices

    case -1: // grow and optimize until automatic stop condition
      mf::LogVerbatim("pma::Track3D") << "optimized segments: " << fSegments.size();
      if ((fSegments.size() >= fSegStopValue) || (fSegments.size() >= fMaxSegStop)) { stop = true; }
      else {
        if (!AddNode(detProp)) stop = true;
      }
      break;

    default: // grow and optimize until fixed number of vertices is added
      if (nNodes > 12) {
        if (AddNode(detProp)) {
          MakeProjection();
          nNodes--;
        }
        else {
          mf::LogVerbatim("pma::Track3D") << "stop (3)";
          stop = true;
          break;
        }

        if (AddNode(detProp)) {
          MakeProjection();
          nNodes--;
          if (AddNode(detProp)) nNodes--;
        }
      }
      else if (nNodes > 4) {
        if (AddNode(detProp)) {
          MakeProjection();
          nNodes--;
        }
        else {
          mf::LogVerbatim("pma::Track3D") << "stop (2)";
          stop = true;
          break;
        }

        if (AddNode(detProp)) nNodes--;
      }
      else {
        if (AddNode(detProp)) { nNodes--; }
        else {
          mf::LogVerbatim("pma::Track3D") << "stop (1)";
          stop = true;
          break;
        }
      }
      break;
    }

  } while (!stop);

  MakeProjection();
  return GetObjFunction();
}

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

Definition at line 928 of file PmaTrack3D.cxx.

{
  pma::Hit3D* hit = nullptr;
  if (index > (int)size()) index = (int)size();
  while (--index >= 0) // look for the prev index of hit from the view
  {
    hit = fHits[index];
    if (hit->View2D() == view) {
      if (inclDisabled)
        break;
      else if (hit->IsEnabled())
        break;
    }
  }
  return index;
}

pma::Segment3D * pma::Track3D::PrevSegment

(

pma::Node3D *

vtx

)

const

Definition at line 1039 of file PmaTrack3D.cxx.

{
  if (vtx->Prev()) {
    auto seg = static_cast<pma::Segment3D*>(vtx->Prev());
    if (seg->Parent() == this) return seg;
  }
  return nullptr;
}

void pma::Track3D::push_back ( pma::Hit3D *  hit )

inline

Definition at line 90 of file PmaTrack3D.h.

  {
    hit->fParent = this;
    fHits.push_back(hit);
  }

bool pma::Track3D::push_back	(	const detinfo::DetectorPropertiesData &	detProp,
		const art::Ptr< recob::Hit > &	hit
	)

Definition at line 352 of file PmaTrack3D.cxx.

{
  for (auto const& trk_hit : fHits) {
    if (trk_hit->fHit == hit) return false;
  }
  pma::Hit3D* h3d = new pma::Hit3D(detProp, hit);
  h3d->fParent = this;
  fHits.push_back(h3d);
  return true;
}

void pma::Track3D::ReassignHitsInTree ( pma::Track3D *  plRoot = nullptr )

private

Definition at line 2150 of file PmaTrack3D.cxx.

{
  bool skipFirst;
  if (trkRoot)
    skipFirst = true;
  else {
    trkRoot = this;
    skipFirst = false;
  }

  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());
  }

  while (vtx) {
    auto segThis = NextSegment(vtx);

    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) seg->Parent()->ReassignHitsInTree(trkRoot);
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  double d0, dmin;
  pma::Hit3D* hit = nullptr;
  pma::Track3D* nearestTrk = nullptr;
  size_t i = 0;
  while (HasTwoViews(2) && (i < size())) {
    hit = fHits[i];
    d0 = hit->GetDistToProj();

    nearestTrk = GetNearestTrkInTree(hit->Point2D(), hit->View2D(), hit->TPC(), hit->Cryo(), dmin);
    if ((nearestTrk != this) && (dmin < 0.5 * d0)) {
      nearestTrk->push_back(release_at(i));
      mf::LogVerbatim("pma::Track3D") << "*** hit moved to another track ***";
    }
    else
      i++;
  }
  if (!size()) { mf::LogError("pma::Track3D") << "ALL hits moved to other tracks."; }
}

void pma::Track3D::RebuildSegments ( )

private

Definition at line 2427 of file PmaTrack3D.cxx.

{
  DeleteSegments();

  fSegments.reserve(fNodes.size() - 1);
  for (size_t i = 1; i < fNodes.size(); i++) {
    fSegments.push_back(new pma::Segment3D(this, fNodes[i - 1], fNodes[i]));
  }
}

pma::Hit3D * pma::Track3D::release_at

(

size_t

index

)

Definition at line 344 of file PmaTrack3D.cxx.

{
  pma::Hit3D* h3d = fHits[index];
  fHits.erase(fHits.begin() + index);
  return h3d;
}

void pma::Track3D::RemoveHits

(

const std::vector< art::Ptr< recob::Hit >> &

hits

)

Remove hits; removes also hit->node/seg assignments.

Definition at line 413 of file PmaTrack3D.cxx.

{
  unsigned int n = 0;
  for (auto const& hit : hits) {
    if (erase(hit)) n++;
  }
  if (n) MakeProjection();
}

bool pma::Track3D::RemoveNode

(

size_t

idx

)

Definition at line 1440 of file PmaTrack3D.cxx.

{
  if ((fNodes.size() > 1) && (idx < fNodes.size())) {
    pma::Node3D* vtx = fNodes[idx];
    fNodes.erase(fNodes.begin() + idx);
    RebuildSegments();

    if (!vtx->NextCount() && !vtx->Prev()) { delete vtx; }

    return true;
  }
  else
    return false;
}

std::vector<pma::Segment3D*> const& pma::Track3D::Segments ( ) const

inlinenoexcept

Definition at line 329 of file PmaTrack3D.h.

  {
    return fSegments;
  }

bool pma::Track3D::SelectAllHits

(

void

)

Definition at line 2942 of file PmaTrack3D.cxx.

{
  bool changed = false;
  for (auto h : fHits) {
    changed |= !(h->IsEnabled());
    h->SetEnabled(true);
  }
  return changed;
}

bool pma::Track3D::SelectHits

(

float

fraction = 1.0F

)

Definition at line 2889 of file PmaTrack3D.cxx.

{
  if (fraction < 0.0F) fraction = 0.0F;
  if (fraction > 1.0F) fraction = 1.0F;
  if (fraction < 1.0F) std::sort(fHits.begin(), fHits.end(), pma::bTrajectory3DDistLess());

  size_t nHitsColl = (size_t)(fraction * NHits(geo::kZ));
  size_t nHitsInd2 = (size_t)(fraction * NHits(geo::kV));
  size_t nHitsInd1 = (size_t)(fraction * NHits(geo::kU));
  size_t coll = 0, ind2 = 0, ind1 = 0;

  bool b, changed = false;
  for (auto h : fHits) {
    b = h->IsEnabled();
    if (fraction < 1.0F) {
      h->SetEnabled(false);
      switch (h->View2D()) {
      case geo::kZ:
        if (coll++ < nHitsColl) h->SetEnabled(true);
        break;
      case geo::kV:
        if (ind2++ < nHitsInd2) h->SetEnabled(true);
        break;
      case geo::kU:
        if (ind1++ < nHitsInd1) h->SetEnabled(true);
        break;
      }
    }
    else
      h->SetEnabled(true);

    changed |= (b != h->IsEnabled());
  }

  if (fraction < 1.0F) {
    FirstElement()->SelectAllHits();
    LastElement()->SelectAllHits();
  }
  return changed;
}

bool pma::Track3D::SelectRndHits	(	size_t	segmax,
		size_t	vtxmax
	)

Definition at line 2931 of file PmaTrack3D.cxx.

{
  bool changed = false;
  for (auto n : fNodes)
    changed |= n->SelectRndHits(vtxmax);
  for (auto s : fSegments)
    changed |= s->SelectRndHits(segmax);
  return changed;
}

void pma::Track3D::SetEndSegWeight ( float  value )

inlinenoexcept

Definition at line 401 of file PmaTrack3D.h.

  {
    fEndSegWeight = value;
  }

void pma::Track3D::SetMaxHitsPerSeg ( unsigned int  value )

inlinenoexcept

Definition at line 423 of file PmaTrack3D.h.

  {
    fMaxHitsPerSeg = value;
  }

void pma::Track3D::SetPenalty ( float  value )

inlinenoexcept

Definition at line 412 of file PmaTrack3D.h.

  {
    fPenaltyFactor = value;
  }

void pma::Track3D::SetT0FromDx	(	const detinfo::DetectorClocksData &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		double	dx
	)

Function to convert dx into dT0.

Definition at line 2377 of file PmaTrack3D.cxx.

{
  auto const* geom = lar::providerFrom<geo::Geometry>();
  const geo::TPCGeo& tpcGeo = geom->TPC(fNodes.front()->TPC(), fNodes.front()->Cryo());

  // GetXTicksCoefficient has a sign that we don't care about. We need to decide
  // the sign for ourselves based on the drift direction of the TPC
  double correctedSign = 0;
  if (tpcGeo.DetectDriftDirection() > 0) {
    if (dx > 0) { correctedSign = 1.0; }
    else {
      correctedSign = -1.0;
    }
  }
  else {
    if (dx > 0) { correctedSign = -1.0; }
    else {
      correctedSign = 1.0;
    }
  }

  // The magnitude of x in ticks is fine
  double dxInTicks = fabs(dx / detProp.GetXTicksCoefficient());
  // Now correct the sign
  dxInTicks = dxInTicks * correctedSign;
  // At this stage, we have xInTicks relative to the trigger time
  double dxInTime = dxInTicks * clockData.TPCClock().TickPeriod();
  // Reconstructed times are relative to the trigger (t=0), so this is our T0
  fT0 = dxInTime;

  mf::LogDebug("pma::Track3D") << dx << ", " << dxInTicks << ", " << correctedSign << ", " << fT0
                               << ", " << tpcGeo.DetectDriftDirection()
                               << " :: " << clockData.TriggerTime() << ", "
                               << clockData.TriggerOffsetTPC() << std::endl;

  // Mark this track as having a measured T0
  fT0Flag = true;
}

void pma::Track3D::SetTagFlag ( ETag  value )

inline

Definition at line 77 of file PmaTrack3D.h.

  {
    fTag = (ETag)(fTag | value);
  }

bool pma::Track3D::ShiftEndsToHits

(

)

Move the first/last Node3D to the first/last hit in the track; returns true if all OK, false if empty segments found.

Definition at line 2480 of file PmaTrack3D.cxx.

{
  pma::Element3D* el;
  pma::Node3D* vtx;

  if (!(fNodes.front()->Prev())) {
    el = GetNearestElement(front()->Point3D());
    vtx = dynamic_cast<pma::Node3D*>(el);
    if (vtx) {
      if (vtx == fNodes.front())
        fNodes.front()->SetPoint3D(front()->Point3D());
      else {
        mf::LogWarning("pma::Track3D") << "First hit is projected to inner node.";
        return false;
      }
    }
    else {
      pma::Segment3D* seg = dynamic_cast<pma::Segment3D*>(el);
      if (seg) {
        if (seg->Prev() == fNodes.front()) {
          double l0 = seg->Length();
          fNodes.front()->SetPoint3D(front()->Point3D());
          if ((seg->Length() < 0.2 * l0) && (fNodes.size() > 2)) {
            pma::Node3D* toRemove = fNodes[1];
            if (toRemove->NextCount() == 1) {
              mf::LogWarning("pma::Track3D")
                << "ShiftEndsToHits(): Short start segment, node removed.";
              fNodes.erase(fNodes.begin() + 1);

              RebuildSegments();
              MakeProjection();

              delete toRemove;
            }
            else {
              mf::LogWarning("pma::Track3D") << "Branching node, not removed.";
              return false;
            }
          }
        }
        else {
          mf::LogWarning("pma::Track3D") << "First hit is projected to inner segment.";
          return false;
        }
      }
    }
  }

  if (!(fNodes.back()->NextCount())) {
    el = GetNearestElement(back()->Point3D());
    vtx = dynamic_cast<pma::Node3D*>(el);
    if (vtx) {
      if (vtx == fNodes.back())
        fNodes.back()->SetPoint3D(back()->Point3D());
      else {
        mf::LogWarning("pma::Track3D") << "First hit is projected to inner node.";
        return false;
      }
    }
    else {
      pma::Segment3D* seg = dynamic_cast<pma::Segment3D*>(el);
      if (seg) {
        if (seg->Next() == fNodes.back()) {
          double l0 = seg->Length();
          fNodes.back()->SetPoint3D(back()->Point3D());
          if ((seg->Length() < 0.2 * l0) && (fNodes.size() > 2)) {
            size_t idx = fNodes.size() - 2;
            pma::Node3D* toRemove = fNodes[idx];
            if (toRemove->NextCount() == 1) {
              mf::LogWarning("pma::Track3D")
                << "ShiftEndsToHits(): Short end segment, node removed.";
              fNodes.erase(fNodes.begin() + idx);

              RebuildSegments();
              MakeProjection();

              delete toRemove;
            }
            else {
              mf::LogWarning("pma::Track3D") << "Branching node, not removed.";
              return false;
            }
          }
        }
        else {
          mf::LogWarning("pma::Track3D") << "First hit is projected to inner segment.";
          return false;
        }
      }
    }
  }

  return true;
}

size_t pma::Track3D::size ( void  ) const

inline

Definition at line 111 of file PmaTrack3D.h.

  {
    return fHits.size();
  }

void pma::Track3D::SortHits

(

void

)

Definition at line 2712 of file PmaTrack3D.cxx.

{
  std::vector<pma::Hit3D*> hits_tmp;
  hits_tmp.reserve(size());

  pma::Node3D* vtx = fNodes.front();
  pma::Segment3D* seg = NextSegment(vtx);
  while (vtx) {
    vtx->SortHits();
    for (size_t i = 0; i < vtx->NHits(); i++) {
      pma::Hit3D* h3d = &(vtx->Hit(i));
      if (h3d->fParent == this) hits_tmp.push_back(h3d);
    }

    if (seg) {
      seg->SortHits();
      for (size_t i = 0; i < seg->NHits(); i++) {
        pma::Hit3D* h3d = &(seg->Hit(i));
        if (h3d->fParent == this) hits_tmp.push_back(h3d);
      }

      vtx = static_cast<pma::Node3D*>(seg->Next());
      seg = NextSegment(vtx);
    }
    else
      break;
  }

  if (size() == hits_tmp.size()) {
    for (size_t i = 0; i < size(); i++) {
      fHits[i] = hits_tmp[i];
    }
  }
  else
    mf::LogError("pma::Track3D") << "Hit sorting problem.";
}

void pma::Track3D::SortHitsInTree

(

bool

skipFirst = false

)

Definition at line 2226 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());
  }

  while (vtx) {
    auto segThis = NextSegment(vtx);

    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) seg->Parent()->SortHitsInTree(true);
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  SortHits();
}

pma::Track3D * pma::Track3D::Split	(	detinfo::DetectorPropertiesData const &	detProp,
		size_t	idx,
		bool	try_start_at_idx = true
	)

Definition at line 1456 of file PmaTrack3D.cxx.

{
  if (!idx || (idx + 1 >= fNodes.size())) return 0;

  pma::Node3D* n = nullptr;
  pma::Track3D* t0 = new pma::Track3D();
  t0->fT0 = fT0;
  t0->fT0Flag = fT0Flag;
  t0->fTag = fTag;

  for (size_t i = 0; i < idx; ++i) {
    n = fNodes.front();
    n->ClearAssigned();

    pma::Segment3D* s = static_cast<pma::Segment3D*>(n->Prev());
    if (s && (s->Parent() == this)) s->RemoveNext(n);

    size_t k = 0;
    while (k < n->NextCount()) {
      s = static_cast<pma::Segment3D*>(n->Next(k));
      if (s->Parent() == this)
        n->RemoveNext(s);
      else
        k++;
    }

    fNodes.erase(fNodes.begin());
    t0->fNodes.push_back(n);
  }

  n = fNodes.front();
  t0->fNodes.push_back(
    new pma::Node3D(detProp, n->Point3D(), n->TPC(), n->Cryo(), false, n->GetDriftShift()));
  t0->RebuildSegments();
  RebuildSegments();

  size_t h = 0;
  while (h < size()) {
    pma::Hit3D* h3d = fHits[h];
    unsigned int view = h3d->View2D(), tpc = h3d->TPC(), cryo = h3d->Cryo();
    double dist2D_old = Dist2(h3d->Point2D(), view, tpc, cryo);
    double dist2D_new = t0->Dist2(h3d->Point2D(), view, tpc, cryo);

    if ((dist2D_new < dist2D_old) && t0->HasTPC(tpc))
      t0->push_back(release_at(h));
    else if (!HasTPC(tpc) && t0->HasTPC(tpc))
      t0->push_back(release_at(h));
    else
      h++;
  }

  bool passed = true;
  if (HasTwoViews() && t0->HasTwoViews()) {
    if (try_start_at_idx && t0->CanFlip()) {
      mf::LogVerbatim("pma::Track3D") << "  attach new t0 and this trk to a common start node";
      t0->Flip();
      passed = t0->AttachTo(fNodes.front());
    }
    else {
      mf::LogVerbatim("pma::Track3D") << "  attach this trk to the new t0 end node";
      passed = AttachTo(t0->fNodes.back());
    }
  }
  else {
    mf::LogVerbatim("pma::Track3D") << "  single-view track";
    passed = false;
  }

  if (!passed) {
    mf::LogVerbatim("pma::Track3D") << "  undo split";
    while (t0->size())
      push_back(t0->release_at(0));

    for (size_t i = 0; i < idx; ++i) {
      fNodes.insert(fNodes.begin() + i, t0->fNodes.front());
      t0->fNodes.erase(t0->fNodes.begin());
    }

    RebuildSegments();
    delete t0;
    t0 = 0;
  }

  return t0;
}

bool pma::Track3D::SwapVertices ( size_t  v0, size_t  v1  )

private

Definition at line 3210 of file PmaTrack3D.cxx.

{
  if (v0 == v1) return false;

  if (v0 > v1) {
    size_t vx = v0;
    v0 = v1;
    v1 = vx;
  }

  pma::Node3D* vtmp;
  if (v1 - v0 == 1) {
    pma::Segment3D* midSeg = NextSegment(fNodes[v0]);
    pma::Segment3D* prevSeg = PrevSegment(fNodes[v0]);
    pma::Segment3D* nextSeg = NextSegment(fNodes[v1]);

    fNodes[v1]->RemoveNext(nextSeg);
    midSeg->Disconnect();

    vtmp = fNodes[v0];
    fNodes[v0] = fNodes[v1];
    fNodes[v1] = vtmp;

    if (prevSeg) prevSeg->AddNext(fNodes[v0]);
    fNodes[v0]->AddNext(midSeg);
    midSeg->AddNext(fNodes[v1]);
    if (nextSeg) fNodes[v1]->AddNext(nextSeg);

    return false;
  }
  else {
    vtmp = fNodes[v0];
    fNodes[v0] = fNodes[v1];
    fNodes[v1] = vtmp;
    return true;
  }
}

unsigned int pma::Track3D::TestHits	(	detinfo::DetectorPropertiesData const &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	hits,
		double	dist = 0.4
	)		const

Count close 2D hits.

Definition at line 859 of file PmaTrack3D.cxx.

{
  if (hits.empty()) {
    mf::LogWarning("pma::Track3D") << "TestHits(): Empty cluster.";
    return 0;
  }

  TVector3 p3d;
  double tst, d2 = dist * dist;
  unsigned int nhits = 0;
  for (auto const& h : hits)
    if (GetUnconstrainedProj3D(detProp, h, p3d, tst) && tst < d2) ++nhits;

  return nhits;
}

double pma::Track3D::TestHitsMse	(	detinfo::DetectorPropertiesData const &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	hits,
		bool	normalized = true
	)		const

MSE of 2D hits.

Definition at line 833 of file PmaTrack3D.cxx.

{
  if (hits.empty()) {
    mf::LogWarning("pma::Track3D") << "TestHitsMse(): Empty cluster.";
    return -1.0;
  }

  double mse = 0.0;
  for (auto const& h : hits) {
    unsigned int cryo = h->WireID().Cryostat;
    unsigned int tpc = h->WireID().TPC;
    unsigned int view = h->WireID().Plane;
    unsigned int wire = h->WireID().Wire;
    float drift = h->PeakTime();

    mse += Dist2(pma::WireDriftToCm(detProp, wire, drift, view, tpc, cryo), view, tpc, cryo);
  }
  if (normalized)
    return mse / hits.size();
  else
    return mse;
}

std::vector< unsigned int > pma::Track3D::TPCs

(

)

const

Definition at line 453 of file PmaTrack3D.cxx.

{
  std::vector<unsigned int> tpc_idxs;
  for (auto hit : fHits) {
    unsigned int tpc = hit->TPC();

    bool found = false;
    for (unsigned int const tpc_idx : tpc_idxs)
      if (tpc_idx == tpc) {
        found = true;
        break;
      }

    if (!found) tpc_idxs.push_back(tpc);
  }
  return tpc_idxs;
}

double pma::Track3D::TuneFullTree	(	double	eps = 0.001,
		double	gmax = 50.0
	)

Definition at line 2279 of file PmaTrack3D.cxx.

{
  if (size_t depth = 1; !UpdateParamsInTree(false, depth)) {
    mf::LogError("pma::Track3D") << "TuneFullTree failed.";
    return -2; // negetive to tag destroyed tree
  }

  double g0 = GetObjFnInTree(), g1 = 0.0;
  if (!std::isfinite(g0)) {
    mf::LogError("pma::Track3D") << "Infinifte value of g.";
    return -2; // negetive to tag destroyed tree
  }
  if (g0 > gmax) {
    mf::LogWarning("pma::Track3D") << "Too high value of g: " << g0;
    return -1; // negetive to tag bad tree
  }
  if (g0 == 0.0) return g0;

  mf::LogVerbatim("pma::Track3D") << "Tune tree, g = " << g0;
  unsigned int stepIter = 0;
  do {
    double gstep = 1.0;
    unsigned int iter = 0;
    while ((gstep > eps) && (iter < 60)) {
      g1 = g0;
      g0 = TuneSinglePass();

      MakeProjectionInTree();

      if (size_t depth = 1; !UpdateParamsInTree(false, depth)) {
        g0 = -2;
        break;
      } // negetive to tag destroyed tree

      if (g0 == 0.0F) break;

      gstep = fabs(g0 - g1) / g0;
      iter++;
    }

    stepIter++;

  } while (stepIter < 5);

  MakeProjectionInTree();
  SortHitsInTree();

  if (g0 >= 0) { mf::LogVerbatim("pma::Track3D") << "  done, g = " << g0; }
  else {
    mf::LogError("pma::Track3D") << "TuneFullTree failed.";
  }
  return g0;
}

double pma::Track3D::TuneSinglePass

(

bool

skipFirst = false

)

Definition at line 2042 of file PmaTrack3D.cxx.

{
  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());
  }

  double g = 0.0;
  while (vtx) {
    vtx->Optimize(fPenaltyValue, fEndSegWeight);
    auto segThis = NextSegment(vtx);

    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) g += seg->Parent()->TuneSinglePass(true);
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  return g + GetObjFunction();
}

void pma::Track3D::UpdateHitsRadius ( )

private

Definition at line 3285 of file PmaTrack3D.cxx.

{
  std::vector<pma::Hit3D*> hitsColl, hitsInd1, hitsInd2;
  for (auto hit : fHits) {
    switch (hit->View2D()) {
    case geo::kZ: hitsColl.push_back(hit); break;
    case geo::kV: hitsInd2.push_back(hit); break;
    case geo::kU: hitsInd1.push_back(hit); break;
    }
  }
  fHitsRadius = std::max({pma::GetHitsRadius2D(hitsColl, true),
                          pma::GetHitsRadius2D(hitsInd2, true),
                          pma::GetHitsRadius2D(hitsInd1, true)});
}

bool pma::Track3D::UpdateParams ( )

private

Definition at line 3183 of file PmaTrack3D.cxx.

{
  size_t n = size();
  if (n == 0) {
    fPenaltyValue = 1;
    fSegStopValue = 1;
    return false;
  }

  float nCubeRoot = pow((double)n, 1.0 / 3.0);
  float avgDist2Root = sqrt(AverageDist2());
  if (avgDist2Root > 0) {
    fPenaltyValue = fPenaltyFactor * pow((double)fSegments.size(), 1.8) * avgDist2Root /
                    (fHitsRadius * nCubeRoot);

    fSegStopValue = (int)(fSegStopFactor * nCubeRoot * fHitsRadius / avgDist2Root);
    if (fSegStopValue < fMinSegStop) fSegStopValue = fMinSegStop;
    return true;
  }
  else {
    fPenaltyValue = 1;
    fSegStopValue = 1;
    return false;
  }
}

bool pma::Track3D::UpdateParamsInTree	(	bool	skipFirst,
		size_t &	depth
	)

Definition at line 2003 of file PmaTrack3D.cxx.

{
  constexpr size_t maxTreeDepth = 100; // really big tree...

  pma::Node3D* vtx = fNodes.front();

  if (skipFirst) {
    if (auto segThis = NextSegment(vtx)) vtx = static_cast<pma::Node3D*>(segThis->Next());

    if (++depth > maxTreeDepth) { mf::LogError("pma::Track3D") << "Tree depth above allowed max."; }
  }

  bool isOK = true;

  while (vtx) {
    auto segThis = NextSegment(vtx);

    for (size_t i = 0; i < vtx->NextCount(); i++) {
      auto seg = static_cast<pma::Segment3D*>(vtx->Next(i));
      if (seg != segThis) { isOK &= seg->Parent()->UpdateParamsInTree(true, depth); }
    }

    if (segThis)
      vtx = static_cast<pma::Node3D*>(segThis->Next());
    else
      break;
  }

  if (!UpdateParams()) {
    mf::LogError("pma::Track3D") << "Track empty.";
    isOK = false;
  }

  --depth;

  return isOK;
}

void pma::Track3D::UpdateProjection

(

void

)

Definition at line 3151 of file PmaTrack3D.cxx.

{
  for (auto n : fNodes)
    n->UpdateProjection();
  for (auto s : fSegments)
    s->UpdateProjection();
}

std::pair< TVector2, TVector2 > pma::Track3D::WireDriftRange	(	detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	view,
		unsigned int	tpc,
		unsigned int	cryo
	)		const

Rectangular region of the track 2D projection in view/tpc/cryo; first in the returned pair is (min_wire; min_drift), second is (max_wire; max_drift). Used for preselection of neighbouring hits in the track validation functions.

Definition at line 491 of file PmaTrack3D.cxx.

{
  std::pair<TVector2, TVector2> range(TVector2(0., 0.), TVector2(0., 0.));

  size_t n0 = 0;
  while ((n0 < fNodes.size()) && (fNodes[n0]->TPC() != (int)tpc))
    n0++;

  if (n0 < fNodes.size()) {
    size_t n1 = n0;
    while ((n1 < fNodes.size()) && (fNodes[n1]->TPC() == (int)tpc))
      n1++;

    if (n0 > 0) n0--;
    if (n1 == fNodes.size()) n1--;

    TVector2 p0 = fNodes[n0]->Projection2D(view);
    p0 = pma::CmToWireDrift(detProp, p0.X(), p0.Y(), view, tpc, cryo);

    TVector2 p1 = fNodes[n1]->Projection2D(view);
    p1 = pma::CmToWireDrift(detProp, p1.X(), p1.Y(), view, tpc, cryo);

    if (p0.X() > p1.X()) {
      double tmp = p0.X();
      p0.Set(p1.X(), p0.Y());
      p1.Set(tmp, p1.Y());
    }
    if (p0.Y() > p1.Y()) {
      double tmp = p0.Y();
      p0.Set(p0.X(), p1.Y());
      p1.Set(p1.X(), tmp);
    }

    range.first = p0;
    range.second = p1;
  }
  return range;
}

Member Data Documentation

std::vector<TVector3*> pma::Track3D::fAssignedPoints

private

Definition at line 487 of file PmaTrack3D.h.

float pma::Track3D::fEndSegWeight {0.05F}

private

Definition at line 509 of file PmaTrack3D.h.

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

private

Definition at line 485 of file PmaTrack3D.h.

float pma::Track3D::fHitsRadius {1.0F}

private

Definition at line 510 of file PmaTrack3D.h.

unsigned int pma::Track3D::fMaxHitsPerSeg {70}

private

Definition at line 499 of file PmaTrack3D.h.

unsigned int pma::Track3D::fMaxSegStop {2}

private

Definition at line 505 of file PmaTrack3D.h.

float pma::Track3D::fMaxSegStopFactor {8.0F}

private

Definition at line 501 of file PmaTrack3D.h.

unsigned int pma::Track3D::fMinSegStop {2}

private

Definition at line 504 of file PmaTrack3D.h.

std::vector<pma::Node3D*> pma::Track3D::fNodes

private

Definition at line 496 of file PmaTrack3D.h.

float pma::Track3D::fPenaltyFactor {1.0F}

private

Definition at line 500 of file PmaTrack3D.h.

float pma::Track3D::fPenaltyValue {0.1F}

private

Definition at line 508 of file PmaTrack3D.h.

std::vector<pma::Segment3D*> pma::Track3D::fSegments

private

Definition at line 497 of file PmaTrack3D.h.

float pma::Track3D::fSegStopFactor {0.2F}

private

Definition at line 507 of file PmaTrack3D.h.

unsigned int pma::Track3D::fSegStopValue {2}

private

Definition at line 503 of file PmaTrack3D.h.

double pma::Track3D::fT0 {}

private

Definition at line 512 of file PmaTrack3D.h.

bool pma::Track3D::fT0Flag {false}

private

Definition at line 513 of file PmaTrack3D.h.

ETag pma::Track3D::fTag {kNotTagged}

private

Definition at line 515 of file PmaTrack3D.h.

---

The documentation for this class was generated from the following files:

• PmaTrack3D.h
• PmaTrack3D.cxx