pma::ProjectionMatchingAlg Class Reference

#include "ProjectionMatchingAlg.h"

Classes
struct	Config

Public Member Functions
	ProjectionMatchingAlg (const Config &config)
	ProjectionMatchingAlg (const fhicl::ParameterSet &pset)
double	validate_on_adc (const pma::Track3D &trk, const img::DataProviderAlg &adcImage, float thr) const
double	validate_on_adc_test (const pma::Track3D &trk, const img::DataProviderAlg &adcImage, const std::vector< art::Ptr< recob::Hit > > &hits, TH1F *histoPassing, TH1F *histoRejected) const
double	validate (const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit > > &hits) const
double	validate (const TVector3 &p0, const TVector3 &p1, const std::vector< art::Ptr< recob::Hit > > &hits, unsigned int testView, unsigned int tpc, unsigned int cryo) const
double	twoViewFraction (pma::Track3D &trk) const
unsigned int	testHits (const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit > > &hits, double eps=1.0) const
	Count the number of hits that are closer than eps * fHitTestingDist2D to the track 2D projection. More...
bool	isContained (const pma::Track3D &trk, float margin=0.0F) const
pma::Track3D *	buildTrack (const std::vector< art::Ptr< recob::Hit > > &hits_1, const std::vector< art::Ptr< recob::Hit > > &hits_2=std::vector< art::Ptr< recob::Hit > >()) const
	wire planes; number of segments used to create the track depends on the number of hits. More...
pma::Track3D *	buildMultiTPCTrack (const std::vector< art::Ptr< recob::Hit > > &hits) const
	as far as hits origin from at least two wire planes. More...
pma::Track3D *	buildShowerSeg (const std::vector< art::Ptr< recob::Hit > > &hits, const pma::Vector3D &vtx) const
	Build a shower segment from sets of hits and attached to the provided vertex. More...
pma::Track3D *	buildSegment (const std::vector< art::Ptr< recob::Hit > > &hits_1, const std::vector< art::Ptr< recob::Hit > > &hits_2=std::vector< art::Ptr< recob::Hit > >()) const
pma::Track3D *	buildSegment (const std::vector< art::Ptr< recob::Hit > > &hits_1, const std::vector< art::Ptr< recob::Hit > > &hits_2, const TVector3 &point) const
pma::Track3D *	buildSegment (const std::vector< art::Ptr< recob::Hit > > &hits, const TVector3 &point) const
void	FilterOutSmallParts (double r2d, const std::vector< art::Ptr< recob::Hit > > &hits_in, std::vector< art::Ptr< recob::Hit > > &hits_out, const TVector2 &vtx2d) const
void	RemoveNotEnabledHits (pma::Track3D &trk) const
pma::Track3D *	extendTrack (const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit > > &hits, bool add_nodes) const
	Add more hits to an existing track, reoptimize, optionally add more nodes. More...
void	guideEndpoints (pma::Track3D &trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > &hits) const
	Add 3D reference points to clean endpoints of a track (both need to be in the same TPC). More...
void	guideEndpoints (pma::Track3D &trk, pma::Track3D::ETrackEnd endpoint, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > &hits) const
	Add 3D reference points to clean endpoint of a track. More...
std::vector< pma::Hit3D * >	trimTrackToVolume (pma::Track3D &trk, TVector3 p0, TVector3 p1) const
bool	alignTracks (pma::Track3D &first, pma::Track3D &second) const
void	mergeTracks (pma::Track3D &dst, pma::Track3D &src, bool reopt) const
void	autoFlip (pma::Track3D &trk, pma::Track3D::EDirection dir=Track3D::kForward, double thr=0.0, unsigned int n=0) const
double	selectInitialHits (pma::Track3D &trk, unsigned int view=geo::kZ, unsigned int *nused=0) const

Private Member Functions
bool	chkEndpointHits (int wire, int wdir, double drift_x, int view, unsigned int tpc, unsigned int cryo, const pma::Track3D &trk, const std::vector< art::Ptr< recob::Hit > > &hits) const
bool	addEndpointRef (pma::Track3D &trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > &hits, std::pair< int, int > const *wires, double const *xPos, unsigned int tpc, unsigned int cryo) const
bool	GetCloseHits (double r2d, const std::vector< art::Ptr< recob::Hit > > &hits_in, std::vector< size_t > &used, std::vector< art::Ptr< recob::Hit > > &hits_out) const
bool	Has (const std::vector< size_t > &v, size_t idx) const
void	ShortenSeg (pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const
bool	TestTrk (pma::Track3D &trk, const geo::TPCGeo &tpcgeom) const

Static Private Member Functions
static size_t	getSegCount (size_t trk_size)

Private Attributes
double	fOptimizationEps
double	fFineTuningEps
double	fTrkValidationDist2D
double	fHitTestingDist2D
double	fMinTwoViewFraction
geo::GeometryCore const *	fGeom
detinfo::DetectorProperties const *	fDetProp

Detailed Description

Definition at line 52 of file ProjectionMatchingAlg.h.

Constructor & Destructor Documentation

pma::ProjectionMatchingAlg::ProjectionMatchingAlg

(

const Config &

config

)

Definition at line 17 of file ProjectionMatchingAlg.cxx.

References fDetProp, fFineTuningEps, fHitTestingDist2D, pma::ProjectionMatchingAlg::Config::FineTuningEps, fMinTwoViewFraction, fOptimizationEps, fTrkValidationDist2D, pma::ProjectionMatchingAlg::Config::HitTestingDist2D, pma::ProjectionMatchingAlg::Config::HitWeightU, pma::ProjectionMatchingAlg::Config::HitWeightV, pma::ProjectionMatchingAlg::Config::HitWeightZ, geo::kU, geo::kV, geo::kZ, pma::ProjectionMatchingAlg::Config::MinTwoViewFraction, pma::ProjectionMatchingAlg::Config::NodeMargin3D, pma::ProjectionMatchingAlg::Config::OptimizationEps, pma::Node3D::SetMargin(), pma::Element3D::SetOptFactor(), and pma::ProjectionMatchingAlg::Config::TrkValidationDist2D.

                                                                                            :
    fGeom( &*(art::ServiceHandle<geo::Geometry>()) ),
        fDetProp(lar::providerFrom<detinfo::DetectorPropertiesService>())
{
        fOptimizationEps = config.OptimizationEps();
        fFineTuningEps = config.FineTuningEps();

        fTrkValidationDist2D = config.TrkValidationDist2D();
        fHitTestingDist2D = config.HitTestingDist2D();

        fMinTwoViewFraction = config.MinTwoViewFraction();

        pma::Node3D::SetMargin(config.NodeMargin3D());

        pma::Element3D::SetOptFactor(geo::kU, config.HitWeightU());
        pma::Element3D::SetOptFactor(geo::kV, config.HitWeightV());
        pma::Element3D::SetOptFactor(geo::kZ, config.HitWeightZ());

        fDetProp = lar::providerFrom<detinfo::DetectorPropertiesService>();
}

pma::ProjectionMatchingAlg::ProjectionMatchingAlg ( const fhicl::ParameterSet &  pset )

inline

Definition at line 108 of file ProjectionMatchingAlg.h.

References hits(), twoViewFraction(), validate(), validate_on_adc(), validate_on_adc_test(), and lar::dump::vector().

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

Member Function Documentation

bool pma::ProjectionMatchingAlg::addEndpointRef ( pma::Track3D &  trk, const std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > &  hits, std::pair< int, int > const *  wires, double const *  xPos, unsigned int  tpc, unsigned int  cryo  ) const

private

Definition at line 996 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddRefPoint(), chkEndpointHits(), fGeom, hits(), geo::GeometryCore::IntersectionPoint(), pma::Track3D::size(), x, y, and z.

Referenced by autoFlip(), and guideEndpoints().

{
        double x = 0.0, y = 0.0, z = 0.0;
        std::vector< std::pair<int, unsigned int> > wire_view;
        for (unsigned int i = 0; i < 3; i++)
                if (wires[i].first >= 0)
        {
                const auto hiter = hits.find(i);
                if (hiter != hits.end())
                {
                        if (chkEndpointHits(wires[i].first, wires[i].second, xPos[i], i, tpc, cryo, trk, hiter->second))
                        {
                                x += xPos[i];
                                wire_view.push_back(std::pair<int, unsigned int>(wires[i].first, i));
                        }
                }
        }
        if (wire_view.size() > 1)
        {
                x /= wire_view.size();
                fGeom->IntersectionPoint(
                        wire_view[0].first, wire_view[1].first,
                        wire_view[0].second, wire_view[1].second,
                        cryo, tpc, y, z);
                trk.AddRefPoint(x, y, z);
                mf::LogVerbatim("ProjectionMatchingAlg") << "trk tpc:" << tpc << " size:" << trk.size()
                        << " add ref.point (" << x << "; " << y << "; " << z << ")";
                return true;
        }
        else
        {
                mf::LogVerbatim("ProjectionMatchingAlg") << "trk tpc:" << tpc << " size:" << trk.size()
                        << " wire-plane-parallel track, but need two clean views of endpoint";
                return false;
        }
}

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

Flip tracks to get second as a continuation of first; returns false if not possible (tracks in reversed order).

Definition at line 1265 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), pma::Dist2(), pma::Track3D::Flip(), pma::Track3D::front(), and pma::Hit3D::Point3D().

Referenced by isContained(), and mergeTracks().

{
        unsigned int k = 0;
        double distFF = pma::Dist2(first.front()->Point3D(), second.front()->Point3D());
        double dist = distFF;

        double distFB = pma::Dist2(first.front()->Point3D(), second.back()->Point3D());
        if (distFB < dist) { k = 1; dist = distFB; }

        double distBF = pma::Dist2(first.back()->Point3D(), second.front()->Point3D());
        if (distBF < dist) { k = 2; dist = distBF; }

        double distBB = pma::Dist2(first.back()->Point3D(), second.back()->Point3D());
        if (distBB < dist) { k = 3; dist = distBB; }

        switch (k) // flip to get dst end before src start, do not merge if track's order reversed
        {
                case 0: first.Flip(); break;
                case 1: mf::LogError("PMAlgTrackMaker") << "Tracks in reversed order."; return false;
                case 2: break;
                case 3: second.Flip(); break;
                default: throw cet::exception("pma::ProjectionMatchingAlg") << "alignTracks: should never happen." << std::endl;
        }
        return true;
}

void pma::ProjectionMatchingAlg::autoFlip ( pma::Track3D &  trk, pma::Track3D::EDirection  dir = Track3D::kForward, double  thr = 0.0, unsigned int  n = 0  ) const

inline

Try to correct track direction of the stopping particle: dir: kForward - particle stop is at the end of the track; kBackward - particle stop is at the beginning of the track; dQ/dx difference has to be above thr to actually flip the track; compares dQ/dx of n hits at each end of the track (default is based on the track length).

Definition at line 235 of file ProjectionMatchingAlg.h.

References addEndpointRef(), pma::Track3D::AutoFlip(), chkEndpointHits(), dir, GetCloseHits(), getSegCount(), Has(), hits(), geo::kZ, n, selectInitialHits(), ShortenSeg(), TestTrk(), and lar::dump::vector().

                                                            { trk.AutoFlip(dir, thr, n); };

pma::Track3D * pma::ProjectionMatchingAlg::buildMultiTPCTrack

(

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

hits

)

const

as far as hits origin from at least two wire planes.

Build a track from sets of hits, multiple TPCs are OK (like taken from PFParticles),

Definition at line 481 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), buildSegment(), d, pma::Dist2(), fFineTuningEps, fOptimizationEps, pma::Track3D::front(), getSegCount(), hits(), mergeTracks(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), pma::Track3D::Segments(), pma::Track3D::SelectAllHits(), pma::Track3D::size(), and pma::Track3D::SortHits().

Referenced by pma::PMAlgFitter::buildTracks(), and isContained().

{
        std::map< unsigned int, std::vector< art::Ptr<recob::Hit> > > hits_by_tpc;
        for (auto const & h : hits)
        {
                hits_by_tpc[h->WireID().TPC].push_back(h);
        }

        std::vector< pma::Track3D* > tracks;
        for(auto const & hsel : hits_by_tpc)
        {
                //pma::Track3D* trk = buildTrack(hsel.second);
                pma::Track3D* trk = buildSegment(hsel.second);
                if (trk) tracks.push_back(trk);
        }

        bool need_reopt = false;
        while (tracks.size() > 1)
        {
                need_reopt = true;

                pma::Track3D* first = tracks.front();
                pma::Track3D* best = 0;
                double d, dmin = 1.0e12;
                size_t t_best = 0, cfg = 0;
                for (size_t t = 1; t < tracks.size(); ++t)
                {
                        pma::Track3D* second = tracks[t];

                        d = pma::Dist2(first->front()->Point3D(), second->front()->Point3D());
                        if (d < dmin) { dmin = d; best = second; t_best = t; cfg = 0; }

                        d = pma::Dist2(first->front()->Point3D(), second->back()->Point3D());
                        if (d < dmin) { dmin = d; best = second; t_best = t; cfg = 1; }

                        d = pma::Dist2(first->back()->Point3D(), second->front()->Point3D());
                        if (d < dmin) { dmin = d; best = second; t_best = t; cfg = 2; }

                        d = pma::Dist2(first->back()->Point3D(), second->back()->Point3D());
                        if (d < dmin) { dmin = d; best = second; t_best = t; cfg = 3; }
                }
                if (best)
                {
                        switch (cfg)
                        {
                                default:
                                case 0:
                                case 1:
                                        mergeTracks(*best, *first, false);
                                        tracks[0] = best; delete first; break;

                                case 2:
                                case 3:
                                        mergeTracks(*first, *best, false);
                                        delete best; break;
                        }
                        tracks.erase(tracks.begin() + t_best);
                }
                else break; // should not happen
        }

        pma::Track3D* trk = 0;
        if (!tracks.empty())
        {
                trk = tracks.front();
                if (need_reopt)
                {
                        double g = trk->Optimize(0, fOptimizationEps);
                        mf::LogVerbatim("ProjectionMatchingAlg") << "  reopt after merging initial tpc segments: done, g = " << g;
                }

                int nSegments = getSegCount(trk->size()) - trk->Segments().size() + 1;
                if (nSegments > 0) // need to add segments
                {
                        double g = 0.0;
                        size_t nNodes = (size_t)( nSegments - 1 ); // n nodes to add
                        if (nNodes)
                        {
                                mf::LogVerbatim("ProjectionMatchingAlg") << "  optimize trk (add " << nSegments << " seg)";

                                g = trk->Optimize(nNodes, fOptimizationEps, false, true, 25, 10);   // build nodes
                                mf::LogVerbatim("ProjectionMatchingAlg") << "  nodes done, g = " << g;
                                trk->SelectAllHits();
                        }
                        g = trk->Optimize(0, fFineTuningEps);      // final tuning
                        mf::LogVerbatim("ProjectionMatchingAlg") << "  tune done, g = " << g;
                }
                trk->SortHits();
        }
        return trk;
}

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

Build a straight segment from two sets of hits (they should origin from two wire planes); method is intendet for short tracks or shower initial parts, where only a few hits per plane are available and there is no chance to see a curvature or any other features.

Definition at line 854 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddHits(), fFineTuningEps, pma::Track3D::HasTwoViews(), pma::Track3D::Initialize(), pma::Track3D::Optimize(), pma::Track3D::SetEndSegWeight(), pma::Track3D::SortHits(), and pma::Track3D::TPCs().

Referenced by buildMultiTPCTrack(), buildSegment(), buildShowerSeg(), shower::EMShowerAlg::ConstructTrack(), isContained(), ems::EMShower3D::Make3DSeg(), shower::EMShowerAlg::MakeShower(), ems::EMShower3D::Reoptimize(), and ems::EMShower3D::Validate().

{
        pma::Track3D* trk = new pma::Track3D();
        trk->SetEndSegWeight(0.001F);
        trk->AddHits(hits_1);
        trk->AddHits(hits_2);

        if (trk->HasTwoViews() &&
           (trk->TPCs().size() == 1)) // now only in single tpc
        {
                if (!trk->Initialize(0.001F))
                {
                        mf::LogWarning("ProjectionMatchingAlg") << "initialization failed, delete segment";
                        delete trk;
                        return 0;
                }
                trk->Optimize(0, fFineTuningEps);

                trk->SortHits();
                return trk;
        }
        else
        {
                mf::LogWarning("ProjectionMatchingAlg") << "need at least two views in single tpc";
                delete trk;
                return 0;
        }
}

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

Build a straight segment from two sets of hits (they should origin from two wire planes), starting from a given point (like vertex known from another algorithm); method is intendet for short tracks or shower initial parts, where only a few hits per plane are available and there is no chance to see a curvature or any other features.

Definition at line 886 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), buildSegment(), pma::Dist2(), fFineTuningEps, pma::Track3D::Flip(), pma::Track3D::front(), pma::Track3D::Nodes(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), and pma::Track3D::SortHits().

{
        pma::Track3D* trk = buildSegment(hits_1, hits_2);

        if (trk)
        {
                double dfront = pma::Dist2(trk->front()->Point3D(), point);
                double dback = pma::Dist2(trk->back()->Point3D(), point);
                if (dfront > dback) trk->Flip();

                trk->Nodes().front()->SetPoint3D(point);
                trk->Nodes().front()->SetFrozen(true);
                trk->Optimize(0, fFineTuningEps);

                trk->SortHits();
        }
        return trk;
}

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

Build a straight segment from set of hits (they should origin from two wire planes at least), starting from a given point.

Definition at line 909 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::back(), buildSegment(), pma::Dist2(), fFineTuningEps, pma::Track3D::Flip(), pma::Track3D::front(), hits(), pma::Track3D::Nodes(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), and pma::Track3D::SortHits().

{
        pma::Track3D* trk = buildSegment(hits);

        if (trk)
        {
                double dfront = pma::Dist2(trk->front()->Point3D(), point);
                double dback = pma::Dist2(trk->back()->Point3D(), point);
                if (dfront > dback) trk->Flip();

                trk->Nodes().front()->SetPoint3D(point);
                trk->Nodes().front()->SetFrozen(true);
                trk->Optimize(0, fFineTuningEps);

                trk->SortHits();
        }
        return trk;
}

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

Build a shower segment from sets of hits and attached to the provided vertex.

Definition at line 576 of file ProjectionMatchingAlg.cxx.

References buildSegment(), geo::GeometryCore::Cryostat(), fGeom, FilterOutSmallParts(), geo::GeometryCore::FindCryostatAtPosition(), geo::GeometryCore::FindTPCAtPosition(), pma::GetProjectionToPlane(), geo::TPCGeo::HasPlane(), geo::GeometryCore::HasTPC(), hits(), ShortenSeg(), pma::Track3D::size(), geo::CryostatGeo::TPC(), and geo::TPCID::TPC.

Referenced by pma::PMAlgFitter::buildShowers(), and isContained().

{
        double vtxarray[3] {vtx.X(), vtx.Y(), vtx.Z()};

        if (!fGeom->HasTPC(fGeom->FindTPCAtPosition(vtxarray))) return 0;

        TVector3 vtxv3(vtx.X(), vtx.Y(), vtx.Z());
        
        const size_t tpc = fGeom->FindTPCAtPosition(vtxarray).TPC;
        const size_t cryo = fGeom->FindCryostatAtPosition(vtxarray);
        const geo::TPCGeo& tpcgeom = fGeom->Cryostat(cryo).TPC(tpc);

        // use only hits from tpc where the vtx is
        std::vector<art::Ptr<recob::Hit> > hitstpc;
        for (size_t h = 0; h < hits.size(); ++h)
                if (hits[h]->WireID().TPC == tpc)
                        hitstpc.push_back(hits[h]);

        if (!hitstpc.size()) return 0;
        
        std::vector<art::Ptr<recob::Hit> > hitstrk;
        size_t view = 0; size_t countviews = 0;
        while (view < 3)
        {
                mf::LogVerbatim("ProjectionMatchinAlg") << "collecting hits from view: " << view;
                if (!tpcgeom.HasPlane(view)) {++view; continue;}

                // select hits only for a single view
                std::vector<art::Ptr<recob::Hit> > hitsview;
                for (size_t h = 0; h < hitstpc.size(); ++h)
                        if (hitstpc[h]->WireID().Plane == view)
                                hitsview.push_back(hitstpc[h]);
                if (!hitsview.size()) {++view; continue;}

                // filter our small groups of hits, far from main shower                
                std::vector<art::Ptr<recob::Hit> > hitsfilter;
                TVector2 proj_pr = pma::GetProjectionToPlane(vtxv3, view, tpc, cryo);
        
                mf::LogVerbatim("ProjectionMatchinAlg") << "start filter out: ";
                FilterOutSmallParts(2.0, hitsview, hitsfilter, proj_pr);
                mf::LogVerbatim("ProjectionMatchingAlg") << "after filter out";

                for (size_t h = 0; h < hitsfilter.size(); ++h)
                        hitstrk.push_back(hitsfilter[h]);

                if (hitsfilter.size() >= 5) countviews++;

                ++view;
        }

        if (!hitstrk.size() || (countviews < 2)) 
        {
                mf::LogWarning("ProjectionMatchinAlg") << "too few hits, segment not built";
                return 0;
        }

        // hits are prepared, finally segment is built

        pma::Track3D* trk = new pma::Track3D();
        trk = buildSegment(hitstrk, vtxv3);
        
        // make shorter segment to estimate direction more precise
        ShortenSeg(*trk, tpcgeom);

        if (trk->size() < 10)
        {
                mf::LogWarning("ProjectionMatchingAlg") << "too short segment, delete segment";
                delete trk;
                return 0;
        }

        return trk;
}

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

wire planes; number of segments used to create the track depends on the number of hits.

Build a track from two sets of hits from single TPC, hits should origin from at least two

Definition at line 424 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddHits(), f, fFineTuningEps, fMinTwoViewFraction, fOptimizationEps, getSegCount(), pma::Track3D::Initialize(), geo::kU, geo::kV, geo::kZ, pma::Track3D::NHits(), pma::Track3D::Optimize(), pma::Track3D::SelectAllHits(), pma::Track3D::size(), pma::Track3D::SortHits(), pma::Track3D::TPCs(), and twoViewFraction().

Referenced by isContained(), and pma::PMAlgTracker::matchCluster().

{
        pma::Track3D* trk = new pma::Track3D(); // track candidate
        trk->AddHits(hits_1);
        trk->AddHits(hits_2);

        mf::LogVerbatim("ProjectionMatchingAlg") << "track size: " << trk->size();
        std::vector< unsigned int > tpcs = trk->TPCs();
        for (size_t t = 0; t < tpcs.size(); ++t)
        {
                mf::LogVerbatim("ProjectionMatchingAlg") << "  tpc:" << tpcs[t];
        }
        mf::LogVerbatim("ProjectionMatchingAlg")
                << "  #coll:" << trk->NHits(geo::kZ)
                << "  #ind2:" << trk->NHits(geo::kV)
                << "  #ind1:" << trk->NHits(geo::kU);

        size_t nSegments = getSegCount(trk->size());
        size_t nNodes = (size_t)( nSegments - 1 ); // n nodes to add

        mf::LogVerbatim("ProjectionMatchingAlg") << "  initialize trk";
        if (!trk->Initialize())
        {
                mf::LogWarning("ProjectionMatchingAlg") << "  initialization failed, delete trk";
                delete trk;
                return 0;
        }

        double f = twoViewFraction(*trk);
        if (f > fMinTwoViewFraction)
        {
                double g = 0.0;
                mf::LogVerbatim("ProjectionMatchingAlg") << "  optimize trk (" << nSegments << " seg)";
                if (nNodes)
                {
                        g = trk->Optimize(nNodes, fOptimizationEps, false, true, 25, 10);   // build nodes
                        mf::LogVerbatim("ProjectionMatchingAlg") << "  nodes done, g = " << g;
                        trk->SelectAllHits();
                }
                g = trk->Optimize(0, fFineTuningEps);              // final tuning
                mf::LogVerbatim("ProjectionMatchingAlg") << "  tune done, g = " << g;

                trk->SortHits();
                // trk->ShiftEndsToHits(); // not sure if useful already here
                return trk;
        }
        else
        {
                mf::LogVerbatim("ProjectionMatchingAlg") << "  clusters do not match, f = " << f;
                delete trk;
                return 0;
        }
}

bool pma::ProjectionMatchingAlg::chkEndpointHits ( int  wire, int  wdir, double  drift_x, int  view, unsigned int  tpc, unsigned int  cryo, const pma::Track3D &  trk, const std::vector< art::Ptr< recob::Hit > > &  hits  ) const

private

Definition at line 963 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorProperties::ConvertTicksToX(), fDetProp, hits(), pma::Track3D::size(), and x.

Referenced by addEndpointRef(), and autoFlip().

{
        size_t nCloseHits = 0;
        int forwWires = 3, backWires = -1;
        double xMargin = 0.4;
        for (auto h : hits)
                if ((view == (int)h->WireID().Plane) &&
                    (tpc == h->WireID().TPC) &&
                        (cryo == h->WireID().Cryostat))
        {
                bool found = false;
                for (size_t ht = 0; ht < trk.size(); ht++)
                        if (trk[ht]->Hit2DPtr().key() == h.key())
                {
                        found = true; break;
                }
                if (found) continue;

                int dw = wdir * (h->WireID().Wire - wire);
                if ((dw <= forwWires) && (dw >= backWires))
                {
                        double x = fDetProp->ConvertTicksToX(h->PeakTime(), view, tpc, cryo);
                        if (fabs(x - drift_x) < xMargin) nCloseHits++;
                }
        }
        if (nCloseHits > 1) return false;
        else return true;
}

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

Add more hits to an existing track, reoptimize, optionally add more nodes.

Definition at line 931 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddHits(), fFineTuningEps, fOptimizationEps, getSegCount(), hits(), geo::kU, geo::kV, geo::kZ, pma::Track3D::NHits(), pma::Track3D::Nodes(), pma::Track3D::Optimize(), and pma::Track3D::size().

Referenced by pma::PMAlgTracker::extendTrack(), isContained(), and pma::PMAlgTracker::reassignHits_1().

{
        pma::Track3D* copy = new pma::Track3D(trk);
        copy->AddHits(hits);

        mf::LogVerbatim("ProjectionMatchingAlg") << "ext. track size: " << copy->size()
                << "  #coll:" << copy->NHits(geo::kZ)
                << "  #ind2:" << copy->NHits(geo::kV)
                << "  #ind1:" << copy->NHits(geo::kU);

        if (add_nodes)
        {
                size_t nSegments = getSegCount(copy->size());
                int nNodes = nSegments - copy->Nodes().size() + 1; // n nodes to add
                if (nNodes < 0) nNodes = 0;

                if (nNodes)
                {
                        mf::LogVerbatim("ProjectionMatchingAlg") << "  add " << nNodes << " nodes";
                        copy->Optimize(nNodes, fOptimizationEps);
                }
        }
        double g = copy->Optimize(0, fFineTuningEps);
        mf::LogVerbatim("ProjectionMatchingAlg") << "  reopt done, g = " << g;

        return copy;
}

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

Get rid of small groups of hits around cascades; used to calculate cascade starting direction using the compact core cluster.

Definition at line 654 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), GetCloseHits(), and pma::WireDriftToCm().

Referenced by buildShowerSeg(), and isContained().

{
        size_t min_size = hits_in.size() / 5;
        if (min_size < 3) min_size = 3;

        std::vector<size_t> used;
        std::vector< std::vector< art::Ptr<recob::Hit> > > sub_groups;
        std::vector< art::Ptr<recob::Hit> > close_hits;

        float mindist2 = 99999.99; size_t id_sub_groups_save = 0; size_t id = 0;
        while (GetCloseHits(r2d, hits_in, used, close_hits)) 
        {
                
                sub_groups.emplace_back(close_hits);

                for (size_t h = 0; h < close_hits.size(); ++h)
                {
                        TVector2 hi_cm = pma::WireDriftToCm(close_hits[h]->WireID().Wire, 
                        close_hits[h]->PeakTime(), 
                        close_hits[h]->WireID().Plane, 
                        close_hits[h]->WireID().TPC, 
                        close_hits[h]->WireID().Cryostat);

                        float dist2 = pma::Dist2(hi_cm, vtx2d);
                        if (dist2 < mindist2)
                        {
                                id_sub_groups_save = id;
                                mindist2 = dist2;
                        }
                }
                
                id++;
        }

        for (size_t i = 0; i < sub_groups.size(); ++i)
        {
                if (i == id_sub_groups_save)
                {
                        for (auto h : sub_groups[i]) hits_out.push_back(h);
                }
        }

        for (size_t i = 0; i < sub_groups.size(); ++i)
        {
                if ((i != id_sub_groups_save) && (hits_out.size() < 10) && (sub_groups[i].size() > min_size))
                        for (auto h : sub_groups[i]) hits_out.push_back(h);
        }
}

bool pma::ProjectionMatchingAlg::GetCloseHits ( double  r2d, const std::vector< art::Ptr< recob::Hit > > &  hits_in, std::vector< size_t > &  used, std::vector< art::Ptr< recob::Hit > > &  hits_out  ) const

private

Definition at line 709 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), fDetProp, fGeom, detinfo::DetectorProperties::GetXTicksCoefficient(), Has(), recob::Hit::PeakTime(), geo::TPCGeo::Plane(), geo::GeometryCore::TPC(), geo::WireID::Wire, recob::Hit::WireID(), and geo::PlaneGeo::WirePitch().

Referenced by autoFlip(), and FilterOutSmallParts().

{
        hits_out.clear();

        const double gapMargin = 5.0; // can be changed to f(id_tpc1, id_tpc2)
        size_t idx = 0;

        while ((idx < hits_in.size()) && Has(used, idx)) idx++;

        if (idx < hits_in.size())
        {               
                hits_out.push_back(hits_in[idx]);
                used.push_back(idx);

                unsigned int tpc = hits_in[idx]->WireID().TPC;
                unsigned int cryo = hits_in[idx]->WireID().Cryostat;
                unsigned int view = hits_in[idx]->WireID().Plane;
                double wirePitch = fGeom->TPC(tpc, cryo).Plane(view).WirePitch();
                double driftPitch = fDetProp->GetXTicksCoefficient(tpc, cryo);

                double r2d2 = r2d*r2d;
                double gapMargin2 = sqrt(2 * gapMargin*gapMargin);
                gapMargin2 = (gapMargin2 + r2d)*(gapMargin2 + r2d);

                bool collect = true;
                while (collect)
                {
                        collect = false;
                        for (size_t i = 0; i < hits_in.size(); i++)
                                if (!Has(used, i))
                                {
                                        art::Ptr<recob::Hit> const & hi = hits_in[i];
                                        TVector2 hi_cm(wirePitch * hi->WireID().Wire, driftPitch * hi->PeakTime());

                                        bool accept = false;
                
                                        for (auto const & ho : hits_out)                                        
                                        {

                                                TVector2 ho_cm(wirePitch * ho->WireID().Wire, driftPitch * ho->PeakTime());
                                                double d2 = pma::Dist2(hi_cm, ho_cm);
                                                
                                                if (d2 < r2d2) { accept = true; break; }
                                        }
                                        if (accept)
                                        {
                                                collect = true;
                                                hits_out.push_back(hi);
                                                used.push_back(i);
                                        }
                                }
                }
                return true;
        }
        else return false;
}

size_t pma::ProjectionMatchingAlg::getSegCount ( size_t  trk_size )

staticprivate

Definition at line 415 of file ProjectionMatchingAlg.cxx.

Referenced by autoFlip(), buildMultiTPCTrack(), buildTrack(), and extendTrack().

{
        int nSegments = (int)( 0.8 * trk_size / sqrt(trk_size) );

        if (nSegments > 1) return (size_t)nSegments;
        else return 1;
}

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

Add 3D reference points to clean endpoints of a track (both need to be in the same TPC).

Definition at line 1036 of file ProjectionMatchingAlg.cxx.

References addEndpointRef(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), detinfo::DetectorProperties::ConvertTicksToX(), fDetProp, fFineTuningEps, fGeom, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), pma::Segment3D::GetDirection3D(), geo::TPCGeo::HasPlane(), hits(), pma::Element3D::Length(), pma::Track3D::NextHit(), pma::Track3D::Optimize(), pma::Track3D::PrevHit(), pma::Track3D::Segments(), pma::Track3D::size(), and geo::GeometryCore::TPC().

Referenced by pma::PMAlgTrackingBase::guideEndpoints(), and isContained().

{
        unsigned int tpc = trk.FrontTPC(), cryo = trk.FrontCryo();
        if ((tpc != trk.BackTPC()) || (cryo != trk.BackCryo()))
        {
                mf::LogWarning("ProjectionMatchingAlg") << "Please, apply before TPC stitching.";
                return;
        }

        const double maxCosXZ = 0.992546; // 7 deg

        pma::Segment3D* segFront = trk.Segments().front();
        if (trk.Segments().size() > 2)
        {
                pma::Segment3D* segFront1 = trk.Segments()[1];
                if ((segFront->Length() < 0.8) && (segFront1->Length() > 5.0))
                        segFront = segFront1;
        }
        pma::Vector3D dirFront = segFront->GetDirection3D();
        pma::Vector3D dirFrontXZ(dirFront.X(), 0., dirFront.Z());
        dirFrontXZ *= 1.0 / dirFrontXZ.R();

        pma::Segment3D* segBack = trk.Segments().back();
        if (trk.Segments().size() > 2)
        {
                pma::Segment3D* segBack1 = trk.Segments()[trk.Segments().size() - 2];
                if ((segBack->Length() < 0.8) && (segBack1->Length() > 5.0))
                        segBack = segBack1;
        }
        pma::Vector3D dirBack = segBack->GetDirection3D();
        pma::Vector3D dirBackXZ(dirBack.X(), 0., dirBack.Z());
        dirBackXZ *= 1.0 / dirBackXZ.R();

        if ((fabs(dirFrontXZ.Z()) < maxCosXZ) && (fabs(dirBackXZ.Z()) < maxCosXZ))
        {
                return; // front & back are not parallel to wire planes => exit
        }

        unsigned int nPlanesFront = 0, nPlanesBack = 0;
        std::pair<int, int> wiresFront[3], wiresBack[3]; // wire index; index direction
        double xFront[3], xBack[3];

        for (unsigned int i = 0; i < 3; i++)
        {
                bool frontPresent = false, backPresent = false;
                if (fGeom->TPC(tpc, cryo).HasPlane(i))
                {
                        int idxFront0 = trk.NextHit(-1, i);
                        int idxBack0 = trk.PrevHit(trk.size(), i);
                        if ((idxFront0 >= 0) && (idxFront0 < (int)trk.size()) &&
                            (idxBack0 >= 0) && (idxBack0 < (int)trk.size()))
                        {
                                int idxFront1 = trk.NextHit(idxFront0, i);
                                int idxBack1 = trk.PrevHit(idxBack0, i);
                                if ((idxFront1 >= 0) && (idxFront1 < (int)trk.size()) &&
                                    (idxBack1 >= 0) && (idxBack1 < (int)trk.size()))
                                {
                                        int wFront0 = trk[idxFront0]->Wire();
                                        int wBack0 = trk[idxBack0]->Wire();

                                        int wFront1 = trk[idxFront1]->Wire();
                                        int wBack1 = trk[idxBack1]->Wire();

                                        wiresFront[i].first = wFront0;
                                        wiresFront[i].second = wFront0 - wFront1;
                                        xFront[i] = fDetProp->ConvertTicksToX(trk[idxFront0]->PeakTime(), i, tpc, cryo);

                                        wiresBack[i].first = wBack0;
                                        wiresBack[i].second = wBack0 - wBack1;
                                        xBack[i] = fDetProp->ConvertTicksToX(trk[idxBack0]->PeakTime(), i, tpc, cryo);

                                        if (wiresFront[i].second)
                                        {
                                                if (wiresFront[i].second > 0) wiresFront[i].second = 1;
                                                else wiresFront[i].second = -1;

                                                frontPresent = true;
                                                nPlanesFront++;
                                        }

                                        if (wiresBack[i].second)
                                        {
                                                if (wiresBack[i].second > 0) wiresBack[i].second = 1;
                                                else wiresBack[i].second = -1;

                                                backPresent = true;
                                                nPlanesBack++;
                                        }
                                }
                        }
                }
                if (!frontPresent) { wiresFront[i].first = -1; }
                if (!backPresent) { wiresBack[i].first = -1; }
        }

        bool refAdded = false;
        if ((nPlanesFront > 1) && (fabs(dirFrontXZ.Z()) >= maxCosXZ))
        {
                refAdded |= addEndpointRef(trk, hits, wiresFront, xFront, tpc, cryo);
        }

        if ((nPlanesBack > 1) && (fabs(dirBackXZ.Z()) >= maxCosXZ))
        {
                refAdded |= addEndpointRef(trk, hits, wiresBack, xBack, tpc, cryo);
        }
        if (refAdded)
        {
                mf::LogVerbatim("ProjectionMatchingAlg") << "guide wire-plane-parallel track endpoints";
                double g = trk.Optimize(0, 0.1 * fFineTuningEps);
                mf::LogVerbatim("ProjectionMatchingAlg") << "  done, g = " << g;
        }
}

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

Add 3D reference points to clean endpoint of a track.

Definition at line 1151 of file ProjectionMatchingAlg.cxx.

References addEndpointRef(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), detinfo::DetectorProperties::ConvertTicksToX(), fDetProp, fFineTuningEps, fGeom, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), pma::Segment3D::GetDirection3D(), geo::TPCGeo::HasPlane(), hits(), pma::Track3D::kBegin, pma::Element3D::Length(), pma::Track3D::NextHit(), pma::Track3D::Optimize(), art::productstatus::present(), pma::Track3D::PrevHit(), pma::Track3D::Segments(), pma::Track3D::size(), and geo::GeometryCore::TPC().

{
        const double maxCosXZ = 0.992546; // 7 deg

        unsigned int tpc, cryo;
        pma::Segment3D* seg0 = 0;
        pma::Segment3D* seg1 = 0;

        if (endpoint == pma::Track3D::kBegin)
        {
                tpc = trk.FrontTPC(), cryo = trk.FrontCryo();
                seg0 = trk.Segments().front();
                if (trk.Segments().size() > 2)
                {
                        seg1 = trk.Segments()[1];
                }
        }
        else
        {
                tpc = trk.BackTPC(), cryo = trk.BackCryo();
                seg0 = trk.Segments().back();
                if (trk.Segments().size() > 2)
                {
                        seg1 = trk.Segments()[trk.Segments().size() - 2];
                }
        }
        if (seg1 && (seg0->Length() < 0.8) && (seg1->Length() > 5.0))
        {
                seg0 = seg1;
        }
        pma::Vector3D dir0 = seg0->GetDirection3D();
        pma::Vector3D dir0XZ(dir0.X(), 0., dir0.Z());
        dir0XZ *= 1.0 / dir0XZ.R();

        if (fabs(dir0XZ.Z()) < maxCosXZ) { return; } // not parallel to wire planes => exit

        unsigned int nPlanes = 0;
        std::pair<int, int> wires[3]; // wire index; index direction
        double x0[3];

        for (unsigned int i = 0; i < 3; i++)
        {
                bool present = false;
                if (fGeom->TPC(tpc, cryo).HasPlane(i))
                {
                        int idx0 = -1, idx1 = -1;
                        if (endpoint == pma::Track3D::kBegin)
                        {
                                idx0 = trk.NextHit(-1, i);
                        }
                        else
                        {
                                idx0 = trk.PrevHit(trk.size(), i);
                        }

                        if ((idx0 >= 0) && (idx0 < (int)trk.size()) &&
                            (trk[idx0]->TPC() == tpc) && (trk[idx0]->Cryo() == cryo))
                        {
                                if (endpoint == pma::Track3D::kBegin)
                                {
                                        idx1 = trk.NextHit(idx0, i);
                                }
                                else
                                {
                                        idx1 = trk.PrevHit(idx0, i);
                                }

                                if ((idx1 >= 0) && (idx1 < (int)trk.size()) &&
                                    (trk[idx1]->TPC() == tpc) && (trk[idx1]->Cryo() == cryo))
                                {
                                        int w0 = trk[idx0]->Wire();
                                        int w1 = trk[idx1]->Wire();

                                        wires[i].first = w0;
                                        wires[i].second = w0 - w1;
                                        x0[i] = fDetProp->ConvertTicksToX(trk[idx0]->PeakTime(), i, tpc, cryo);

                                        if (wires[i].second)
                                        {
                                                if (wires[i].second > 0) wires[i].second = 1;
                                                else wires[i].second = -1;

                                                present = true;
                                                nPlanes++;
                                        }
                                }
                        }
                }
                if (!present) { wires[i].first = -1; }
        }

        if ((nPlanes > 1) && (fabs(dir0XZ.Z()) >= maxCosXZ) &&
             addEndpointRef(trk, hits, wires, x0, tpc, cryo))
        {
                mf::LogVerbatim("ProjectionMatchingAlg") << "guide wire-plane-parallel track endpoint";
                double g = trk.Optimize(0, 0.1 * fFineTuningEps);
                mf::LogVerbatim("ProjectionMatchingAlg") << "  done, g = " << g;
        }
}

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

private

Definition at line 834 of file ProjectionMatchingAlg.cxx.

Referenced by autoFlip(), and GetCloseHits().

{
    for (auto c : v) if (c == idx) return true;
    return false;
}

bool pma::ProjectionMatchingAlg::isContained ( const pma::Track3D &  trk, float  margin = 0.0F  ) const

inline

Test if hits at the track endpoinds do not stick out of TPC which they belong to. Here one can implement some configurable margin if needed for real data imeprfections.

Definition at line 153 of file ProjectionMatchingAlg.h.

References alignTracks(), pma::Track3D::back(), buildMultiTPCTrack(), buildSegment(), buildShowerSeg(), buildTrack(), extendTrack(), FilterOutSmallParts(), pma::Track3D::FirstElement(), pma::Track3D::front(), guideEndpoints(), hits(), pma::Track3D::LastElement(), mergeTracks(), pma::Hit3D::Point3D(), RemoveNotEnabledHits(), pma::Node3D::SameTPC(), trimTrackToVolume(), and lar::dump::vector().

Referenced by pma::PMAlgTracker::matchCluster(), and pma::PMAlgTracker::reassignHits_1().

        {
                return (trk.FirstElement()->SameTPC(trk.front()->Point3D(), margin) &&
                        trk.LastElement()->SameTPC(trk.back()->Point3D(), margin));
        }

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

Add src to dst as it was its continuation; nodes of src are added to dst after its own nodes, hits of src are added to hits of dst, then dst is reoptimized.

Definition at line 1291 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::AddNode(), alignTracks(), pma::Track3D::BackCryo(), pma::Track3D::BackTPC(), pma::Element3D::Cryo(), pma::Dist2(), fFineTuningEps, pma::Track3D::FrontCryo(), pma::Track3D::FrontTPC(), pma::Element3D::IsFrozen(), pma::Track3D::Length(), pma::Track3D::MakeProjection(), n, pma::Track3D::Nodes(), pma::Track3D::Optimize(), pma::Track3D::push_back(), pma::Track3D::ShiftEndsToHits(), pma::Track3D::size(), pma::Track3D::SortHits(), and pma::Element3D::TPC().

Referenced by buildMultiTPCTrack(), isContained(), and pma::PMAlgTracker::mergeCoLinear().

{
        if (!alignTracks(dst, src)) return;

        unsigned int tpc = src.FrontTPC();
        unsigned int cryo = src.FrontCryo();
        double lmean = dst.Length() / (dst.Nodes().size() - 1);
        if ((pma::Dist2(
                        dst.Nodes().back()->Point3D(),
                        src.Nodes().front()->Point3D()) > 0.5 * lmean) ||
            (tpc != dst.BackTPC()) || (cryo != dst.BackCryo()))
        {
                dst.AddNode(src.Nodes().front()->Point3D(), tpc, cryo);
                if (src.Nodes().front()->IsFrozen())
                        dst.Nodes().back()->SetFrozen(true);
        }
        for (size_t n = 1; n < src.Nodes().size(); n++)
        {
                pma::Node3D* node = src.Nodes()[n];

                dst.AddNode(src.Nodes()[n]->Point3D(),
                            node->TPC(), node->Cryo());

                if (node->IsFrozen())
                        dst.Nodes().back()->SetFrozen(true);
        }
        for (size_t h = 0; h < src.size(); h++)
        {
                dst.push_back(src[h]->Hit2DPtr());
        }
        if (reopt)
        {
                double g = dst.Optimize(0, fFineTuningEps);
                mf::LogVerbatim("ProjectionMatchingAlg") << "  reopt after merging done, g = " << g;
        }
        else
        {
                dst.MakeProjection();
        }

        dst.SortHits();
        dst.ShiftEndsToHits();

        dst.MakeProjection();
        dst.SortHits();
}

void pma::ProjectionMatchingAlg::RemoveNotEnabledHits

(

pma::Track3D &

trk

)

const

Definition at line 842 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::release_at(), and pma::Track3D::size().

Referenced by isContained(), and ShortenSeg().

{
        size_t h = 0;
        while (h < trk.size())
        {
                if (trk[h]->IsEnabled()) ++h;
                else (trk.release_at(h));
        }
}

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

Intendet to calculate dQ/dx in the initial part of EM cascade; collection view is used by default, but it works also with other projections.

Definition at line 1339 of file ProjectionMatchingAlg.cxx.

References util::absDiff(), pma::Hit3D::GetDistToProj(), pma::Hit3D::GetSegFraction(), pma::Track3D::HitDxByView(), pma::Track3D::NextHit(), pma::Track3D::size(), pma::Hit3D::SummedADC(), pma::Hit3D::TagOutlier(), pma::Hit3D::View2D(), and pma::Hit3D::Wire().

Referenced by autoFlip(), ems::EMShower3D::ConvertFrom(), and ems::EMShower3D::ConvertFrom2().

{
        for (size_t i = 0; i < trk.size(); i++)
        {
                pma::Hit3D* hit = trk[i];
                if (hit->View2D() == view)
                {
                        if ((hit->GetDistToProj() > 0.5) || // more than 0.5cm away away from the segment
                            (hit->GetSegFraction() < -1.0)) // projects before segment start (to check!!!)
                                hit->TagOutlier(true);
                        else hit->TagOutlier(false);
                }
        }

        unsigned int nhits = 0;
        double last_x, dx = 0.0, last_q, dq = 0.0, dqdx = 0.0;
        int ih = trk.NextHit(-1, view);

        pma::Hit3D* hit = trk[ih];
        pma::Hit3D* lastHit = hit;

        if ((ih >= 0) && (ih < (int)trk.size()))
        {
                hit->TagOutlier(true);

                ih = trk.NextHit(ih, view);
                while ((dx < 2.5) && (ih >= 0) && (ih < (int)trk.size()))
                {
                        hit = trk[ih];

                        if (util::absDiff(hit->Wire(), lastHit->Wire()) > 2)
                                break; // break on gap in wire direction

                        last_x = trk.HitDxByView(ih, view);
                        last_q = hit->SummedADC();
                        if (dx + last_x < 3.0)
                        {
                                dq += last_q;
                                dx += last_x;
                                nhits++;
                        }
                        else break;

                        lastHit = hit;
                        ih = trk.NextHit(ih, view);
                }
                while ((ih >= 0) && (ih < (int)trk.size()))
                {
                        hit = trk[ih];
                        hit->TagOutlier(true);
                        ih = trk.NextHit(ih, view);
                }
        }
        else { mf::LogError("ProjectionMatchingAlg") << "Initial part selection failed."; }

        if (!nhits) { mf::LogError("ProjectionMatchingAlg") << "Initial part too short to select useful hits."; }

        if (dx > 0.0) dqdx = dq / dx;

  if (nused) (*nused) = nhits;

        //std::cout << "nhits=" << nhits << ", dq=" << dq << ", dx=" << dx << std::endl;
        return dqdx;
}

void pma::ProjectionMatchingAlg::ShortenSeg ( pma::Track3D &  trk, const geo::TPCGeo &  tpcgeom  ) const

private

Definition at line 772 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), pma::Track3D::front(), pma::Track3D::GetMse(), pma::Track3D::Length(), pma::Track3D::Optimize(), pma::Hit3D::Point3D(), RemoveNotEnabledHits(), pma::Track3D::size(), pma::Track3D::SortHits(), and TestTrk().

Referenced by autoFlip(), and buildShowerSeg().

{
        double mse = trk.GetMse();
        mf::LogWarning("ProjectionMatchingAlg") << "initial value of mse: " << mse;
        while ((mse > 0.5) && TestTrk(trk, tpcgeom))
        {
                mse = trk.GetMse();
                for (size_t h = 0; h < trk.size(); ++h)
                        if (std::sqrt(pma::Dist2(trk[h]->Point3D(), trk.front()->Point3D())) 
                                        > 0.8*trk.Length()) trk[h]->SetEnabled(false);

                RemoveNotEnabledHits(trk);

                // trk.Optimize(0.0001, false); // BUG: first argument missing; tentatively:
                trk.Optimize(0, 0.0001, false);
                trk.SortHits();

                mf::LogWarning("ProjectionMatchingAlg") << " mse: " << mse;
                if (mse == trk.GetMse()) break;
                else mse = trk.GetMse();
        }

        RemoveNotEnabledHits(trk);
}

unsigned int pma::ProjectionMatchingAlg::testHits ( const pma::Track3D &  trk, const std::vector< art::Ptr< recob::Hit > > &  hits, double  eps = 1.0  ) const

inline

Count the number of hits that are closer than eps * fHitTestingDist2D to the track 2D projection.

Definition at line 146 of file ProjectionMatchingAlg.h.

References fHitTestingDist2D, hits(), and pma::Track3D::TestHits().

Referenced by pma::PMAlgTracker::matchCluster(), pma::PMAlgTracker::matchTrack(), pma::PMAlgTracker::reassignHits_1(), and pma::PMAlgTracker::reassignHits_2().

        { return trk.TestHits(hits, eps * fHitTestingDist2D); }

bool pma::ProjectionMatchingAlg::TestTrk ( pma::Track3D &  trk, const geo::TPCGeo &  tpcgeom  ) const

private

Definition at line 799 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), pma::Track3D::front(), geo::TPCGeo::HasPlane(), pma::Track3D::NEnabledHits(), pma::Hit3D::Point3D(), pma::Track3D::size(), and fhicl::has_insertion_operator_impl::test().

Referenced by autoFlip(), and ShortenSeg().

{
        bool test = false; 

        if (tpcgeom.HasPlane(0) && tpcgeom.HasPlane(1))
        {
                if ((trk.NEnabledHits(0) > 5) && (trk.NEnabledHits(1) > 5)) 
                        test = true; 
        }
        else if (tpcgeom.HasPlane(0) && tpcgeom.HasPlane(2))
        {
                if ((trk.NEnabledHits(0) > 5) && (trk.NEnabledHits(2) > 5)) 
                        test = true; 
        }
        else if (tpcgeom.HasPlane(1) && tpcgeom.HasPlane(2))
        {
                if ((trk.NEnabledHits(1) > 5) && (trk.NEnabledHits(2) > 5)) 
                        test = true;
        }

        double length = 0.0;
        for (size_t h = 0; h < trk.size(); ++h)
        {
                if (!trk[h]->IsEnabled()) break;
                length = std::sqrt(pma::Dist2(trk.front()->Point3D(), trk[h]->Point3D()));
        }

        mf::LogWarning("ProjectionMatchingAlg") << "length of segment: " << length;
        if (length < 3.0) test = false; // cm
        
        return test;
}

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

Definition at line 1254 of file ProjectionMatchingAlg.cxx.

Referenced by isContained().

{
        std::vector< pma::Hit3D* > trimmedHits;



        return trimmedHits;
}

double pma::ProjectionMatchingAlg::twoViewFraction

(

pma::Track3D &

trk

)

const

Calculate the fraction of trajectory seen by two 2D projections at least; even a prfect track starts/stops with the hit from one 2D view, then hits from other views come, which results with the fraction value high, but always < 1.0; wrong cluster matchings or incomplete tracks give significantly lower values.

Definition at line 396 of file ProjectionMatchingAlg.cxx.

References pma::Track3D::DisableSingleViewEnds(), pma::Track3D::Length(), pma::Track3D::SelectHits(), and pma::Track3D::size().

Referenced by buildTrack(), pma::PMAlgTracker::matchCluster(), and ProjectionMatchingAlg().

{
        trk.SelectHits();
        trk.DisableSingleViewEnds();

        size_t idx = 0;
        while ((idx < trk.size() - 1) && !trk[idx]->IsEnabled()) idx++;
        double l0 = trk.Length(0, idx + 1);

        idx = trk.size() - 1;
        while ((idx > 1) && !trk[idx]->IsEnabled()) idx--;
        double l1 = trk.Length(idx - 1, trk.size() - 1);

        trk.SelectHits();

        return 1.0 - (l0 + l1) / trk.Length();
}

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

Calculate the fraction of the track that is closer than fTrkValidationDist2D to any hit from hits in their plane (a plane that was not used to build the track). Hits should be preselected, so all belong to the same plane.

Definition at line 231 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorProperties::ConvertTicksToX(), pma::Track3D::Cryos(), pma::Dist2(), pma::Track3D::Dist2(), f, fDetProp, fGeom, pma::Track3D::front(), fTrkValidationDist2D, pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), hits(), geo::TPCGeo::Plane(), geo::GeometryCore::PlaneWireToChannel(), pma::Hit3D::Point3D(), pma::Node3D::SameTPC(), pma::Track3D::Segments(), pma::Element3D::TPC(), geo::GeometryCore::TPC(), pma::Track3D::TPCs(), geo::GeometryCore::WireCoordinate(), pma::Track3D::WireDriftRange(), and geo::PlaneGeo::WirePitch().

Referenced by ProjectionMatchingAlg(), ems::EMShower3D::Validate(), and pma::PMAlgTracker::validate().

{
    if (hits.empty()) { return 0; }

        double max_d = fTrkValidationDist2D;
        double d2, max_d2 = max_d * max_d;
        unsigned int nAll = 0, nPassed = 0;
        unsigned int testPlane = hits.front()->WireID().Plane;

        std::vector< unsigned int > trkTPCs = trk.TPCs();
        std::vector< unsigned int > trkCryos = trk.Cryos();
        std::map< std::pair< unsigned int, unsigned int >, std::pair< TVector2, TVector2 > > ranges;
        std::map< std::pair< unsigned int, unsigned int >, double > wirePitch;
        for (auto c : trkCryos)
                for (auto t : trkTPCs)
                {
                        ranges[std::pair< unsigned int, unsigned int >(t, c)] = trk.WireDriftRange(testPlane, t, c);
                        wirePitch[std::pair< unsigned int, unsigned int >(t, c)] = fGeom->TPC(t, c).Plane(testPlane).WirePitch();
                }

        unsigned int tpc, cryo;
        std::map< std::pair< unsigned int, unsigned int >, std::vector< pma::Vector2D > > all_close_points;

        for (const auto h : hits)
                if (h->WireID().Plane == testPlane)
        {
                tpc = h->WireID().TPC;
                cryo = h->WireID().Cryostat;
                std::pair< unsigned int, unsigned int > tpc_cryo(tpc, cryo);
                std::pair< TVector2, TVector2 > rect = ranges[tpc_cryo];

                if ((h->WireID().Wire > rect.first.X() - 10) &&  // chceck only hits in the rectangle around
                    (h->WireID().Wire < rect.second.X() + 10) && // the track projection, it is faster than
                    (h->PeakTime() > rect.first.Y() - 100) &&    // calculation of trk.Dist2(p2d, testPlane)
                    (h->PeakTime() < rect.second.Y() + 100))
                {
                        TVector2 p2d(wirePitch[tpc_cryo] * h->WireID().Wire, fDetProp->ConvertTicksToX(h->PeakTime(), testPlane, tpc, cryo));

                        d2 = trk.Dist2(p2d, testPlane, tpc, cryo);

                        if (d2 < max_d2) all_close_points[tpc_cryo].emplace_back(p2d.X(), p2d.Y());
                }
        }

        auto const & channelStatus = art::ServiceHandle< lariov::ChannelStatusService >()->GetProvider();

        double step = 0.3;
        // then check how points-close-to-the-track-projection are distributed along the
        // track, namely: are there track sections crossing empty spaces, except dead wires?
        pma::Vector3D p(trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
        for (auto const * seg : trk.Segments())
        {
            if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
            {
                p = seg->End(); continue;
            }
            pma::Vector3D p0 = seg->Start();
            pma::Vector3D p1 = seg->End();

            pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());

                tpc = seg->TPC(); cryo = seg->Cryo();

                pma::Vector3D dc = step * seg->GetDirection3D();

                auto const & points = all_close_points[std::pair< unsigned int, unsigned int >(tpc, cryo)];

            double f = pma::GetSegmentProjVector(p, p0, p1);

                double wirepitch = fGeom->TPC(tpc, cryo).Plane(testPlane).WirePitch();
                while ((f < 1.0) && node->SameTPC(p))
                {
                        pma::Vector2D p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
                        geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
                        if (fGeom->HasWire(wireID))
                        {
                          raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
                          if (channelStatus.IsGood(ch))
                          {
                              if (points.size())
                                  {
                                          p2d.SetX(wirepitch * p2d.X());
                                          for (const auto & h : points)
                                          {
                                                  d2 = pma::Dist2(p2d, h);
                                                  if (d2 < max_d2) { nPassed++; break; }
                                          }
                                  }
                                nAll++;
                          }
                          //else mf::LogVerbatim("ProjectionMatchingAlg")
                          //    << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
                        }

                        p += dc; f = pma::GetSegmentProjVector(p, p0, p1);
                }
                p = seg->End();
        }

        if (nAll > 0)
        {
                double v = nPassed / (double)nAll;
                mf::LogVerbatim("ProjectionMatchingAlg") << "  trk fraction ok: " << v;
                return v;
        }
        else return 1.0;
}

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

Calculate the fraction of the 3D segment that is closer than fTrkValidationDist2D to any hit from hits in the testPlane of TPC/Cryo. Hits from the testPlane are preselected by this function among all provided (so a bit slower than fn above).

Definition at line 341 of file ProjectionMatchingAlg.cxx.

References pma::Dist2(), f, fGeom, fTrkValidationDist2D, pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), hits(), geo::TPCGeo::Plane(), geo::GeometryCore::PlaneWireToChannel(), geo::GeometryCore::TPC(), geo::GeometryCore::WireCoordinate(), pma::WireDriftToCm(), and geo::PlaneGeo::WirePitch().

{
        double step = 0.3;
        double max_d = fTrkValidationDist2D;
        double d2, max_d2 = max_d * max_d;
        unsigned int nAll = 0, nPassed = 0;

        TVector3 p(p0);
        TVector3 dc(p1); dc -= p;
        dc *= step / dc.Mag();

        auto const & channelStatus = art::ServiceHandle< lariov::ChannelStatusService >()->GetProvider();

        double f = pma::GetSegmentProjVector(p, p0, p1);
        double wirepitch = fGeom->TPC(tpc, cryo).Plane(testPlane).WirePitch();
        while (f < 1.0)
        {
                TVector2 p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
                geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
                if (fGeom->HasWire(wireID))
                {
                        raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
                        if (channelStatus.IsGood(ch))
                        {
                                p2d.Set(wirepitch * p2d.X(), p2d.Y());
                                for (const auto & h : hits)
                                        if ((h->WireID().Plane == testPlane) &&
                                            (h->WireID().TPC == tpc) &&
                                            (h->WireID().Cryostat == cryo)) 
                                {
                                        d2 = pma::Dist2(p2d, pma::WireDriftToCm(h->WireID().Wire, h->PeakTime(), testPlane, tpc, cryo));
                                        if (d2 < max_d2) { nPassed++; break; }
                                }
                                nAll++;
                        }
                        //else mf::LogVerbatim("ProjectionMatchingAlg")
                        //      << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
                }

                p += dc; f = pma::GetSegmentProjVector(p, p0, p1);
        }

        if (nAll > 3) // validate actually only if 2D projection in testPlane has some minimum length
        {
                double v = nPassed / (double)nAll;
                mf::LogVerbatim("ProjectionMatchingAlg") << "  segment fraction ok: " << v;
                return v;
        }
        else return 1.0;
}

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

Calculate the fraction of the track that is close to non-empty pixel (above thr value) in the ADC image of the testView (a view that was not used to build the track).

Definition at line 39 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorProperties::ConvertXToTicks(), pma::Track3D::Cryos(), f, fDetProp, fGeom, pma::Track3D::front(), pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), img::DataProviderAlg::Plane(), geo::GeometryCore::PlaneWireToChannel(), pma::Hit3D::Point3D(), img::DataProviderAlg::poolMax(), pma::Node3D::SameTPC(), pma::Track3D::Segments(), pma::Element3D::TPC(), pma::Track3D::TPCs(), and geo::GeometryCore::WireCoordinate().

Referenced by ProjectionMatchingAlg(), and pma::PMAlgTracker::validate().

{
        unsigned int nAll = 0, nPassed = 0;
        unsigned int testPlane = adcImage.Plane();

        std::vector< unsigned int > trkTPCs = trk.TPCs();
        std::vector< unsigned int > trkCryos = trk.Cryos();

        unsigned int tpc, cryo;

        auto const & channelStatus = art::ServiceHandle< lariov::ChannelStatusService >()->GetProvider();

        double step = 0.3;
        // check how pixels with a high signal are distributed along the track
        // namely: are there track sections crossing empty spaces, except dead wires?
        pma::Vector3D p(trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
        for (auto const * seg : trk.Segments())
        {
            if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
            {
                p = seg->End(); continue;
            }
            pma::Vector3D p0 = seg->Start();
            pma::Vector3D p1 = seg->End();

            pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());

                tpc = seg->TPC(); cryo = seg->Cryo();

                pma::Vector3D dc = step * seg->GetDirection3D();

                double f = pma::GetSegmentProjVector(p, p0, p1);
                while ((f < 1.0) && node->SameTPC(p))
                {
                        pma::Vector2D p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
                        geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());

                        int widx = (int)p2d.X();
                        int didx = (int)fDetProp->ConvertXToTicks(p2d.Y(), testPlane, tpc, cryo);

                        if (fGeom->HasWire(wireID))
                        {
                          raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
                          if (channelStatus.IsGood(ch))
                          {
                              float max_adc = adcImage.poolMax(widx, didx, 2); // +/- 2 wires, can be parameterized
                              if (max_adc > thr) nPassed++;

                                  nAll++;
                          }
                          //else mf::LogVerbatim("ProjectionMatchingAlg")
                          //    << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
                        }

                        p += dc; f = pma::GetSegmentProjVector(p, p0, p1);
                }

                p = seg->End(); // need to have it at the end due to the p in the first iter set to the hit position, not segment start
        }

        if (nAll > 0)
        {
                double v = nPassed / (double)nAll;
                mf::LogVerbatim("ProjectionMatchingAlg") << "  trk fraction ok: " << v;
                return v;
        }
        else return 1.0;
}

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

Calculate the fraction of the track that is closer than fTrkValidationDist2D to any hit from hits in the testView (a view that was not used to build the track). Creates also histograms of values in pixels for the passing and rejected points on the track, so the threshold value for the ADC-based calibration can be estimated.

Definition at line 110 of file ProjectionMatchingAlg.cxx.

References detinfo::DetectorProperties::ConvertTicksToX(), detinfo::DetectorProperties::ConvertXToTicks(), pma::Track3D::Cryos(), pma::Dist2(), pma::Track3D::Dist2(), f, fDetProp, fGeom, pma::Track3D::front(), fTrkValidationDist2D, pma::GetSegmentProjVector(), geo::GeometryCore::HasWire(), hits(), img::DataProviderAlg::Plane(), geo::TPCGeo::Plane(), geo::GeometryCore::PlaneWireToChannel(), pma::Hit3D::Point3D(), img::DataProviderAlg::poolMax(), pma::Node3D::SameTPC(), pma::Track3D::Segments(), pma::Element3D::TPC(), geo::GeometryCore::TPC(), pma::Track3D::TPCs(), geo::GeometryCore::WireCoordinate(), pma::Track3D::WireDriftRange(), and geo::PlaneGeo::WirePitch().

Referenced by ProjectionMatchingAlg(), and pma::PMAlgTracker::validate().

{
        double max_d = fTrkValidationDist2D;
        double d2, max_d2 = max_d * max_d;
        unsigned int nAll = 0, nPassed = 0;
        unsigned int testPlane = adcImage.Plane();

        std::vector< unsigned int > trkTPCs = trk.TPCs();
        std::vector< unsigned int > trkCryos = trk.Cryos();
        std::map< std::pair< unsigned int, unsigned int >, std::pair< TVector2, TVector2 > > ranges;
        std::map< std::pair< unsigned int, unsigned int >, double > wirePitch;
        for (auto c : trkCryos)
                for (auto t : trkTPCs)
                {
                        ranges[std::pair< unsigned int, unsigned int >(t, c)] = trk.WireDriftRange(testPlane, t, c);
                        wirePitch[std::pair< unsigned int, unsigned int >(t, c)] = fGeom->TPC(t, c).Plane(testPlane).WirePitch();
                }

        unsigned int tpc, cryo;
        std::map< std::pair< unsigned int, unsigned int >, std::vector< pma::Vector2D > > all_close_points;

        for (const auto h : hits)
                if (h->WireID().Plane == testPlane)
        {
                tpc = h->WireID().TPC;
                cryo = h->WireID().Cryostat;
                std::pair< unsigned int, unsigned int > tpc_cryo(tpc, cryo);
                std::pair< TVector2, TVector2 > rect = ranges[tpc_cryo];

                if ((h->WireID().Wire > rect.first.X() - 10) &&  // chceck only hits in the rectangle around
                    (h->WireID().Wire < rect.second.X() + 10) && // the track projection, it is faster than
                    (h->PeakTime() > rect.first.Y() - 100) &&    // calculation of trk.Dist2(p2d, testPlane)
                    (h->PeakTime() < rect.second.Y() + 100))
                {
                        TVector2 p2d(wirePitch[tpc_cryo] * h->WireID().Wire, fDetProp->ConvertTicksToX(h->PeakTime(), testPlane, tpc, cryo));

                        d2 = trk.Dist2(p2d, testPlane, tpc, cryo);

                        if (d2 < max_d2) { all_close_points[tpc_cryo].emplace_back(p2d.X(), p2d.Y()); }
                }
        }

        auto const & channelStatus = art::ServiceHandle< lariov::ChannelStatusService >()->GetProvider();

        double step = 0.3;
        // then check how points-close-to-the-track-projection are distributed along the
        // track, namely: are there track sections crossing empty spaces, except dead wires?
        pma::Vector3D p(trk.front()->Point3D().X(), trk.front()->Point3D().Y(), trk.front()->Point3D().Z());
        for (auto const * seg : trk.Segments())
        {
            if (seg->TPC() < 0) // skip segments between tpc's, look only at those contained in tpc
            {
                p = seg->End(); continue;
            }
            pma::Vector3D p0 = seg->Start();
            pma::Vector3D p1 = seg->End();

            pma::Node3D* node = static_cast<pma::Node3D*>(seg->Prev());

                tpc = seg->TPC(); cryo = seg->Cryo();

                pma::Vector3D dc = step * seg->GetDirection3D();

                auto const & points = all_close_points[std::pair< unsigned int, unsigned int >(tpc, cryo)];

            double f = pma::GetSegmentProjVector(p, p0, p1);

                double wirepitch = fGeom->TPC(tpc, cryo).Plane(testPlane).WirePitch();
                while ((f < 1.0) && node->SameTPC(p))
                {
                        pma::Vector2D p2d(fGeom->WireCoordinate(p.Y(), p.Z(), testPlane, tpc, cryo), p.X());
                        geo::WireID wireID(cryo, tpc, testPlane, (int)p2d.X());
                                
                        int widx = (int)p2d.X();
                        int didx = (int)fDetProp->ConvertXToTicks(p2d.Y(), testPlane, tpc, cryo);

                        if (fGeom->HasWire(wireID))
                        {
                          raw::ChannelID_t ch = fGeom->PlaneWireToChannel(wireID);
                          if (channelStatus.IsGood(ch))
                          {
                              bool is_close = false;
                              float max_adc = adcImage.poolMax(widx, didx, 2);

                              if (points.size())
                                  {
                                          p2d.SetX(wirepitch * p2d.X());
                                          for (const auto & h : points)
                                          {
                                                  d2 = pma::Dist2(p2d, h);
                                                  if (d2 < max_d2) { is_close = true; nPassed++; break; }
                                          }
                                  }
                                  nAll++;
                                          
                                  // now fill the calibration histograms
                                  if (is_close) { if (histoPassing) histoPassing->Fill(max_adc); }
                                  else { if (histoRejected) histoRejected->Fill(max_adc); }
                          }
                          //else mf::LogVerbatim("ProjectionMatchingAlg")
                          //    << "crossing BAD CHANNEL (wire #" << (int)p2d.X() << ")" << std::endl;
                        }

                        p += dc; f = pma::GetSegmentProjVector(p, p0, p1);
                }
                p = seg->End();
        }

        if (nAll > 0)
        {
                double v = nPassed / (double)nAll;
                mf::LogVerbatim("ProjectionMatchingAlg") << "  trk fraction ok: " << v;
                return v;
        }
        else return 1.0;
}

Member Data Documentation

detinfo::DetectorProperties const* pma::ProjectionMatchingAlg::fDetProp

private

Definition at line 288 of file ProjectionMatchingAlg.h.

Referenced by chkEndpointHits(), GetCloseHits(), guideEndpoints(), ProjectionMatchingAlg(), validate(), validate_on_adc(), and validate_on_adc_test().

double pma::ProjectionMatchingAlg::fFineTuningEps

private

Definition at line 279 of file ProjectionMatchingAlg.h.

Referenced by buildMultiTPCTrack(), buildSegment(), buildTrack(), extendTrack(), guideEndpoints(), mergeTracks(), and ProjectionMatchingAlg().

geo::GeometryCore const* pma::ProjectionMatchingAlg::fGeom

private

Definition at line 287 of file ProjectionMatchingAlg.h.

Referenced by addEndpointRef(), buildShowerSeg(), GetCloseHits(), guideEndpoints(), validate(), validate_on_adc(), and validate_on_adc_test().

double pma::ProjectionMatchingAlg::fHitTestingDist2D

private

Definition at line 282 of file ProjectionMatchingAlg.h.

Referenced by ProjectionMatchingAlg(), and testHits().

double pma::ProjectionMatchingAlg::fMinTwoViewFraction

private

Definition at line 284 of file ProjectionMatchingAlg.h.

Referenced by buildTrack(), and ProjectionMatchingAlg().

double pma::ProjectionMatchingAlg::fOptimizationEps

private

Definition at line 276 of file ProjectionMatchingAlg.h.

Referenced by buildMultiTPCTrack(), buildTrack(), extendTrack(), and ProjectionMatchingAlg().

double pma::ProjectionMatchingAlg::fTrkValidationDist2D

private

Definition at line 281 of file ProjectionMatchingAlg.h.

Referenced by ProjectionMatchingAlg(), validate(), and validate_on_adc_test().

---

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

• larreco/v06_64_02/source/larreco/RecoAlg/ProjectionMatchingAlg.h
• larreco/v06_64_02/source/larreco/RecoAlg/ProjectionMatchingAlg.cxx