#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 37 of file PmaTrack3D.h.

Member Enumeration Documentation

enum pma::Track3D::EDirection

Enumerator
kForward
kBackward

Definition at line 40 of file PmaTrack3D.h.

40 { kForward = -1, kBackward = 1 };

pma::Track3D::kForward

Definition: PmaTrack3D.h:40

pma::Track3D::kBackward

Definition: PmaTrack3D.h:40

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 41 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 39 of file PmaTrack3D.h.

39 { kBegin = -1, kEnd = 1 };

pma::Track3D::kBegin

Definition: PmaTrack3D.h:39

pma::Track3D::kEnd

Definition: PmaTrack3D.h:39

Constructor & Destructor Documentation

pma::Track3D::Track3D ( )

default

pma::Track3D::Track3D ( const Track3D & src )

Definition at line 30 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 63 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 403 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 1264 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 352 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 1271 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 262 of file PmaTrack3D.h.

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

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

inline

Definition at line 267 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 2333 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 1655 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 1687 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 1702 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 1544 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 1584 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 1598 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 672 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 749 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 3159 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 103 of file PmaTrack3D.h.

   {
     return fHits.back();
   }

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

inline

Definition at line 161 of file PmaTrack3D.h.

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

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

inline

Definition at line 156 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 653 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 3248 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 2437 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 113 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 1192 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 471 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 2418 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 2749 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 2575 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 2603 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 1160 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 364 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 1745 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 340 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 533 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 632 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 98 of file PmaTrack3D.h.

   {
     return fHits.front();
   }

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

inline

Definition at line 150 of file PmaTrack3D.h.

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

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

inline

Definition at line 145 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 1793 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 378 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 393 of file PmaTrack3D.h.

   {
     return fEndSegWeight;
   }

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

inlinenoexcept

Definition at line 415 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 1837 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 2611 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 2652 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 2069 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 2251 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 1857 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 404 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 1048 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 1777 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 306 of file PmaTrack3D.h.

   {
     return fT0;
   }

ETag pma::Track3D::GetTag ( ) const

inlinenoexcept

Definition at line 64 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 2674 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 1829 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 313 of file PmaTrack3D.h.

   {
     return fT0Flag;
   }

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

inlinenoexcept

Definition at line 69 of file PmaTrack3D.h.

   {
     return (fTag & value);
   }

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

inline

Definition at line 167 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 442 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 1011 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 945 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 335 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 327 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 121 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 299 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 211 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 77 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 1425 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 604 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 1813 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 345 of file PmaTrack3D.h.

   {
     return fNodes.back();
   }

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

inline

Definition at line 117 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 877 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 2989 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 2952 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 2199 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 432 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 909 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 1024 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 422 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 335 of file PmaTrack3D.h.

   {
     return fNodes;
   }

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

inline

Definition at line 95 of file PmaTrack3D.h.

95 { return fHits[index]; }

ValidateOpDetReco.index

index

Definition: ValidateOpDetReco.py:379

pma::Track3D::fHits

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

Definition: PmaTrack3D.h:482

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

inline

Definition at line 96 of file PmaTrack3D.h.

96 { return fHits[index]; }

ValidateOpDetReco.index

index

Definition: ValidateOpDetReco.py:379

pma::Track3D::fHits

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

Definition: PmaTrack3D.h:482

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 1868 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 927 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 1038 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 87 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 351 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 2149 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 2426 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 343 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 412 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 1439 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 326 of file PmaTrack3D.h.

   {
     return fSegments;
   }

bool pma::Track3D::SelectAllHits ( void )

Definition at line 2941 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 2888 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 2930 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 398 of file PmaTrack3D.h.

   {
     fEndSegWeight = value;
   }

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

inlinenoexcept

Definition at line 420 of file PmaTrack3D.h.

   {
     fMaxHitsPerSeg = value;
   }

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

inlinenoexcept

Definition at line 409 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 2376 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 74 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 2479 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 108 of file PmaTrack3D.h.

   {
     return fHits.size();
   }

void pma::Track3D::SortHits ( void )

Definition at line 2711 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 2225 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 1455 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 3209 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 858 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 832 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 452 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 2278 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 2041 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 3284 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 3182 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 2002 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 3150 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 490 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 484 of file PmaTrack3D.h.

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

private

Definition at line 506 of file PmaTrack3D.h.

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

private

Definition at line 482 of file PmaTrack3D.h.

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

private

Definition at line 507 of file PmaTrack3D.h.

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

private

Definition at line 496 of file PmaTrack3D.h.

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

private

Definition at line 502 of file PmaTrack3D.h.

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

private

Definition at line 498 of file PmaTrack3D.h.

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

private

Definition at line 501 of file PmaTrack3D.h.

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

private

Definition at line 493 of file PmaTrack3D.h.

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

private

Definition at line 497 of file PmaTrack3D.h.

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

private

Definition at line 505 of file PmaTrack3D.h.

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

private

Definition at line 494 of file PmaTrack3D.h.

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

private

Definition at line 504 of file PmaTrack3D.h.

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

private

Definition at line 500 of file PmaTrack3D.h.

double pma::Track3D::fT0 {}

private

Definition at line 509 of file PmaTrack3D.h.

bool pma::Track3D::fT0Flag {false}

private

Definition at line 510 of file PmaTrack3D.h.

ETag pma::Track3D::fTag {kNotTagged}

private

Definition at line 512 of file PmaTrack3D.h.

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

larreco/larreco/RecoAlg/PMAlg/PmaTrack3D.h
larreco/larreco/RecoAlg/PMAlg/PmaTrack3D.cxx

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation