trkf::SpacePointAlg Class Reference

#include "SpacePointAlg.h"

Classes
struct	HitMCInfo

Public Member Functions
	SpacePointAlg (const fhicl::ParameterSet &pset)
	~SpacePointAlg ()
bool	filter () const
bool	merge () const
double	maxDT () const
double	maxS () const
int	minViews () const
bool	enableU () const
bool	enableV () const
bool	enableW () const
void	reconfigure (const fhicl::ParameterSet &pset)
void	update () const
double	correctedTime (const recob::Hit &hit) const
double	separation (const art::PtrVector< recob::Hit > &hits) const
bool	compatible (const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
void	fillSpacePoint (const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
void	fillSpacePoints (std::vector< recob::SpacePoint > &spts, std::multimap< double, KHitTrack > const &trackMap) const
	Fill a collection of space points. More...
void	fillComplexSpacePoint (const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
void	makeSpacePoints (const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
void	makeMCTruthSpacePoints (const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
const art::PtrVector< recob::Hit > &	getAssociatedHits (const recob::SpacePoint &spt) const
void	clearHitMap () const
int	numHitMap () const

Private Member Functions
void	makeSpacePoints (const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts, bool useMC) const

Private Attributes
double	fMaxDT
	Maximum time difference between planes. More...
double	fMaxS
	Maximum space separation between wires. More...
int	fMinViews
	Mininum number of views per space point. More...
bool	fEnableU
	Enable flag (U). More...
bool	fEnableV
	Enable flag (V). More...
bool	fEnableW
	Enable flag (W). More...
bool	fFilter
	Filter flag. More...
bool	fMerge
	Merge flag. More...
bool	fPreferColl
	Sort by collection wire. More...
std::map< const recob::Hit *, HitMCInfo >	fHitMCMap
std::map< int, art::PtrVector< recob::Hit > >	fSptHitMap

Detailed Description

Definition at line 77 of file SpacePointAlg.h.

Constructor & Destructor Documentation

trkf::SpacePointAlg::SpacePointAlg

(

const fhicl::ParameterSet &

pset

)

Definition at line 43 of file SpacePointAlg.cxx.

References reconfigure().

                                                              :
    fMaxDT(0.),
    fMaxS(0.),
    fMinViews(1000),
    fEnableU(false),
    fEnableV(false),
    fEnableW(false),
    fFilter(false),
    fMerge(false),
    fPreferColl(false)
    {
        reconfigure(pset);
    }

trkf::SpacePointAlg::~SpacePointAlg

(

)

Definition at line 60 of file SpacePointAlg.cxx.

    {
    }

Member Function Documentation

void trkf::SpacePointAlg::clearHitMap ( ) const

inline

Definition at line 145 of file SpacePointAlg.h.

Referenced by trkf::SeedFinderAlgorithm::FindSeeds(), cluster::ClusterMatchAlg::Match_SpacePoint(), trkf::Track3DKalmanHit::prepareForInput(), trkf::BezierTrackerModule::produce(), trkf::TrackKalmanCheater::produce(), and trkf::Track3DKalmanHit::produce().

{fSptHitMap.clear();}

bool trkf::SpacePointAlg::compatible	(	const art::PtrVector< recob::Hit > &	hits,
		bool	useMC = false
	)		const

Definition at line 308 of file SpacePointAlg.cxx.

References geo::CryostatID::Cryostat, fHitMCMap, fMaxDT, fMaxS, geo::WireGeo::GetCenter(), detinfo::DetectorProperties::GetXTicksOffset(), geo::WireGeo::HalfL(), trkf::SpacePointAlg::HitMCInfo::pchit, recob::Hit::PeakTime(), geo::PlaneID::Plane, s, art::PtrVector< T >::size(), t1, t2, geo::TPCID::TPC, trkf::SpacePointAlg::HitMCInfo::trackIDs, recob::Hit::View(), recob::Hit::WireID(), and geo::GeometryCore::WireIDToWireGeo().

Referenced by fillSpacePoints(), and makeSpacePoints().

    {
        // Get services.
        
        art::ServiceHandle<geo::Geometry> geom;
        const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
        
        int nhits = hits.size();
        
        // Fewer than two or more than three hits can never be compatible.
        
        bool result = nhits >= 2 && nhits <= 3;
        bool mc_ok = true;
        unsigned int tpc = 0;
        unsigned int cstat = 0;
        
        if(result) {
            
            // First do pairwise tests.
            // Do double loop over hits.
            
            for(int ihit1 = 0; result && ihit1 < nhits-1; ++ihit1) {
                const recob::Hit& hit1 = *(hits[ihit1]);
                geo::WireID hit1WireID = hit1.WireID();
                geo::View_t view1 = hit1.View();
                
                double t1 = hit1.PeakTime() - detprop->GetXTicksOffset(hit1WireID.Plane,hit1WireID.TPC,hit1WireID.Cryostat);
                
                // If using mc information, get a collection of track ids for hit 1.
                // If not using mc information, this section of code will trigger the
                // insertion of a single invalid HitMCInfo object into fHitMCMap.
                
                const HitMCInfo& mcinfo1 = fHitMCMap[(useMC ? &hit1 : 0)];
                const std::vector<int>& tid1 = mcinfo1.trackIDs;
                bool only_neg1 = tid1.size() > 0 && tid1.back() < 0;
                
                // Loop over second hit.
                
                for(int ihit2 = ihit1+1; result && ihit2 < nhits; ++ihit2) {
                    const recob::Hit& hit2 = *(hits[ihit2]);
                    geo::WireID hit2WireID = hit2.WireID();
                    geo::View_t view2 = hit2.View();
                    
                    // Test for same tpc and different views.
                    
                    result = result && hit1WireID.TPC == hit2WireID.TPC && view1 != view2 && hit1WireID.Cryostat == hit2WireID.Cryostat;
                    if(result) {
                        
                        // Remember which tpc and cryostat we are in.
                        
                        tpc = hit1WireID.TPC;
                        cstat = hit1WireID.Cryostat;
                        
                        double t2 = hit2.PeakTime() - detprop->GetXTicksOffset(hit2WireID.Plane,hit2WireID.TPC,hit2WireID.Cryostat);
                        
                        // Test maximum time difference.
                        
                        result = result && std::abs(t1-t2) <= fMaxDT;
                        
                        // Test mc truth.
                        
                        if(result && useMC) {
                            
                            // Test whether hits have a common parent track id.
                            
                            const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
                            std::vector<int> tid2 = mcinfo2.trackIDs;
                            bool only_neg2 = tid2.size() > 0 && tid2.back() < 0;
                            std::vector<int>::iterator it =
                            std::set_intersection(tid1.begin(), tid1.end(),
                                                  tid2.begin(), tid2.end(),
                                                  tid2.begin());
                            tid2.resize(it - tid2.begin());
                            
                            // Hits are compatible if they have parents in common.
                            // If the only parent id in common is negative (-999),
                            // then hits are compatible only if both hits have only
                            // negative parent tracks.
                            
                            bool only_neg3 = tid2.size() > 0 && tid2.back() < 0;
                            mc_ok = tid2.size() > 0 && (!only_neg3 || (only_neg1 && only_neg2));
                            result = result && mc_ok;
                            
                            // If we are still OK, check that either hit is
                            // the nearest neighbor of the other.
                            
                            if(result) {
                                result = mcinfo1.pchit[hit2WireID.Plane] == &hit2 ||
                                mcinfo2.pchit[hit1WireID.Plane] == &hit1;
                            }
                        }
                    }
                }
            }
            
            // If there are exactly three hits, and they pass pairwise tests, check
            // for spatial compatibility.
            
            if(result && nhits == 3) {
                
                // Loop over hits.
                
                double dist[3] = {0., 0., 0.};
                double sinth[3] = {0., 0., 0.};
                double costh[3] = {0., 0., 0.};
                
                for(int i=0; i<3; ++i) {
                    
                    // Get tpc, plane, wire.
                    
                    const recob::Hit& hit = *(hits[i]);
                    geo::WireID hitWireID = hit.WireID();
                    
                    const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
                    if ((hitWireID.TPC != tpc) || (hitWireID.Cryostat != cstat))
                        throw cet::exception("SpacePointAlg") << "compatible(): geometry mismatch\n";
                    
                    // Get angles and distance of wire.
                    
                    double hl = wgeom.HalfL();
                    double xyz[3];
                    double xyz1[3];
                    wgeom.GetCenter(xyz);
                    wgeom.GetCenter(xyz1, hl);
                    double s  = (xyz1[1] - xyz[1]) / hl;
                    double c = (xyz1[2] - xyz[2]) / hl;
                    sinth[hit.WireID().Plane] = s;
                    costh[hit.WireID().Plane] = c;
                    dist[hit.WireID().Plane] = xyz[2] * s - xyz[1] * c;
                }
                
                // Do space cut.
                
                double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0]
                            +(sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1]
                            +(sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
                
                result = result && std::abs(S) < fMaxS;
            }
        }
        
        // Done.
        
        return result;
    }

double trkf::SpacePointAlg::correctedTime

(

const recob::Hit &

hit

)

const

Definition at line 209 of file SpacePointAlg.cxx.

References geo::CryostatID::Cryostat, detinfo::DetectorProperties::GetXTicksOffset(), recob::Hit::PeakTime(), geo::PlaneID::Plane, geo::TPCID::TPC, and recob::Hit::WireID().

Referenced by trkf::SpacePointAna::analyze().

    {
        // Get services.
        
        art::ServiceHandle<geo::Geometry> geom;
        const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
        
        // Correct time for trigger offset and plane-dependent time offsets.
        
        double t = hit.PeakTime() - detprop->GetXTicksOffset(hit.WireID().Plane,hit.WireID().TPC,hit.WireID().Cryostat);
        
        return t;
    }

bool trkf::SpacePointAlg::enableU ( ) const

inline

Definition at line 93 of file SpacePointAlg.h.

Referenced by trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits(), trkf::SpacePointFinder::produce(), and trkf::SpacePointCheater::produce().

{return fEnableU;}

bool trkf::SpacePointAlg::enableV ( ) const

inline

Definition at line 94 of file SpacePointAlg.h.

Referenced by trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits(), trkf::SpacePointFinder::produce(), and trkf::SpacePointCheater::produce().

{return fEnableV;}

bool trkf::SpacePointAlg::enableW ( ) const

inline

Definition at line 95 of file SpacePointAlg.h.

References hits().

Referenced by trkf::SeedFinderAlgorithm::GetSeedsFromSortedHits(), trkf::SpacePointFinder::produce(), and trkf::SpacePointCheater::produce().

{return fEnableW;}

void trkf::SpacePointAlg::fillComplexSpacePoint	(	const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	sptv,
		int	sptid
	)		const

Definition at line 669 of file SpacePointAlg.cxx.

References art::PtrVector< T >::begin(), geo::CryostatID::Cryostat, geo::GeometryCore::Cryostat(), art::PtrVector< T >::end(), art::PtrVector< T >::front(), fSptHitMap, geo::WireGeo::GetCenter(), detinfo::DetectorProperties::GetXTicksCoefficient(), detinfo::DetectorProperties::GetXTicksOffset(), geo::WireGeo::HalfL(), hits(), geo::TPCGeo::Nplanes(), recob::Hit::PeakTime(), geo::PlaneID::Plane, s, recob::Hit::SigmaPeakTime(), art::PtrVector< T >::size(), sw, geo::CryostatGeo::TPC(), geo::TPCID::TPC, w, weight, recob::Hit::WireID(), geo::GeometryCore::WireIDToWireGeo(), and geo::GeometryCore::WirePitch().

Referenced by fillSpacePoints(), and makeSpacePoints().

    {
        // Get services.
        
        art::ServiceHandle<geo::Geometry> geom;
        const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
        
        // Calculate time pitch.
        
        double timePitch =    detprop->GetXTicksCoefficient();            // cm / tick
        
        // Figure out which tpc we are in.
        
        unsigned int tpc0 = 0;
        unsigned int cstat0 = 0;
        int nhits = hits.size();
        if(nhits > 0) {
            tpc0 = hits.front()->WireID().TPC;
            cstat0=hits.front()->WireID().Cryostat;
        }
        
        // Remember associated hits internally.
        
        if(fSptHitMap.count(sptid) != 0)
          throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): hit already present!\n";
        fSptHitMap[sptid] = hits;
        
        // Do a preliminary scan of hits.
        // Determine weight given to hits in each view.
        
        unsigned int nplanes = geom->Cryostat(cstat0).TPC(tpc0).Nplanes();
        std::vector<int> numhits(nplanes, 0);
        std::vector<double> weight(nplanes, 0.);
        
        for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
            ihit != hits.end(); ++ihit) {
            
            const recob::Hit& hit = **ihit;
            geo::WireID hitWireID = hit.WireID();
            /* // kept as assertions for performance reasons
             if (hitWireID.Cryostat != cstat0)
             throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): incompatible cryostat\n";
             if (hitWireID.TPC != tpc0);
             throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): incompatible TPC\n";
             if (hitWireID.Plane >= nplanes);
             throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): unknown plane\n";
             */
            assert(hitWireID.Cryostat == cstat0);
            assert(hitWireID.TPC == tpc0);
            assert(hitWireID.Plane < nplanes);
            ++numhits[hitWireID.Plane];
        }
        
        for(unsigned int plane = 0; plane < nplanes; ++plane) {
            double np = numhits[plane];
            if(np > 0.)
                weight[plane] = 1. / (np*np*np);
        }
        
        // Calculate position and error matrix.
        
        double xyz[3] = {0., 0., 0.};
        double errxyz[6] = {0.,
            0., 0.,
            0., 0., 0.};
        
        // Calculate x using drift times.
        // Loop over all hits and calculate the weighted average drift time.
        // Also calculate time variance and chisquare.
        
        double sumt2w = 0.;
        double sumtw = 0.;
        double sumw = 0.;
        
        for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
            ihit != hits.end(); ++ihit) {
            
            const recob::Hit& hit = **ihit;
            geo::WireID hitWireID = hit.WireID();
            
            // Correct time for trigger offset and view-dependent time offsets.
            
            double t0 = detprop->GetXTicksOffset(hitWireID.Plane,hitWireID.TPC,hitWireID.Cryostat);
            double t = hit.PeakTime() - t0;
            double et = hit.SigmaPeakTime();
            double w = weight[hitWireID.Plane]/(et*et);
            
            sumt2w += w*t*t;
            sumtw += w*t;
            sumw += w;
        }
        
        double drift_time = 0.;
        double var_time = 0.;
        double chisq = 0.;
        if(sumw != 0.) {
            drift_time = sumtw / sumw;
            var_time = sumt2w / sumw - drift_time * drift_time;
            if(var_time < 0.)
                var_time = 0.;
            chisq = sumt2w - sumtw * drift_time;
            if(chisq < 0.)
                chisq = 0.;
        }
        xyz[0] = drift_time * timePitch;
        errxyz[0] = var_time * timePitch * timePitch;
        
        // Calculate y, z using wires (need at least two hits).
        
        if(nhits >= 2) {
            
            // Calculate y and z by chisquare minimization of wire coordinates.
            
            double sw = 0.;    // sum w_i
            double sus = 0.;   // sum w_i u_i sin_th_i
            double suc = 0.;   // sum w_i u_i cos_th_i
            double sc2 = 0.;   // sum w_i cos2_th_i
            double ss2 = 0.;   // sum w_i sin2_th_i
            double ssc = 0.;   // sum w_i sin_th_i cos_th_i
            
            // Loop over points.
            
            for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
                ihit != hits.end(); ++ihit) {
                
                const recob::Hit& hit = **ihit;
                geo::WireID hitWireID = hit.WireID();
                const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
                
                // Calculate angle and wire coordinate in this view.
                
                double hl = wgeom.HalfL();
                double cen[3];
                double cen1[3];
                wgeom.GetCenter(cen);
                wgeom.GetCenter(cen1, hl);
                double s  = (cen1[1] - cen[1]) / hl;
                double c = (cen1[2] - cen[2]) / hl;
                double u = cen[2] * s - cen[1] * c;
                double eu = geom->WirePitch(hitWireID.Plane, hitWireID.TPC) / std::sqrt(12.);
                double w = weight[hitWireID.Plane] / (eu * eu);
                
                // Summations
                
                sw += w;
                sus += w*u*s;
                suc += w*u*c;
                sc2 += w*c*c;
                ss2 += w*s*s;
                ssc += w*s*c;
            }
            
            // Calculate y,z
            
            double denom = sc2 * ss2 - ssc * ssc;
            if(denom != 0.) {
                xyz[1] = (-suc * ss2 + sus * ssc) / denom;
                xyz[2] = (sus * sc2 - suc * ssc) / denom;
                errxyz[2] = ss2 / denom;
                errxyz[4] = ssc / denom;
                errxyz[5] = sc2 / denom;
            }
            
            // Make space point.
            
            recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
            sptv.push_back(spt);
        }
        return;
    }

void trkf::SpacePointAlg::fillSpacePoint	(	const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	sptv,
		int	sptid
	)		const

Definition at line 459 of file SpacePointAlg.cxx.

References art::PtrVector< T >::begin(), geo::CryostatID::Cryostat, art::PtrVector< T >::end(), fSptHitMap, geo::WireGeo::GetCenter(), detinfo::DetectorProperties::GetXTicksCoefficient(), detinfo::DetectorProperties::GetXTicksOffset(), geo::WireGeo::HalfL(), hits(), recob::Hit::PeakTime(), geo::PlaneID::Plane, s, recob::Hit::SigmaPeakTime(), art::PtrVector< T >::size(), sw, geo::TPCID::TPC, w, recob::Hit::WireID(), geo::GeometryCore::WireIDToWireGeo(), and geo::GeometryCore::WirePitch().

Referenced by fillSpacePoints(), and makeSpacePoints().

    {
        // Get services.
        
        art::ServiceHandle<geo::Geometry> geom;
        const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
        
        double timePitch=detprop->GetXTicksCoefficient();
        
        int nhits = hits.size();
        
        // Remember associated hits internally.
        
        if (fSptHitMap.find(sptid) != fSptHitMap.end())
            throw cet::exception("SpacePointAlg") << "fillSpacePoint(): hit already present!\n";
        fSptHitMap[sptid] = hits;
        
        // Calculate position and error matrix.
        
        double xyz[3] = {0., 0., 0.};
        double errxyz[6] = {0.,
            0., 0.,
            0., 0., 0.};
        
        // Calculate x using drift times.
        // Loop over all hits and calculate the weighted average drift time.
        // Also calculate time variance and chisquare.
        
        double sumt2w = 0.;
        double sumtw = 0.;
        double sumw = 0.;
        
        
        for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
            ihit != hits.end(); ++ihit) {
            
            const recob::Hit& hit = **ihit;
            geo::WireID hitWireID = hit.WireID();
            
            // Correct time for trigger offset and view-dependent time offsets.
            
            double t0 = detprop->GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
            double t = hit.PeakTime() - t0;
            double et = hit.SigmaPeakTime();
            double w = 1./(et*et);
            
            sumt2w += w*t*t;
            sumtw += w*t;
            sumw += w;
        }
        
        double drift_time = 0.;
        double var_time = 0.;
        double chisq = 0.;
        if(sumw != 0.) {
            drift_time = sumtw / sumw;
            //var_time = sumt2w / sumw - drift_time * drift_time;
            var_time = 1. / sumw;
            if(var_time < 0.)
                var_time = 0.;
            chisq = sumt2w - sumtw * drift_time;
            if(chisq < 0.)
                chisq = 0.;
        }
        xyz[0] = drift_time * timePitch;
        errxyz[0] = var_time * timePitch * timePitch;
        
        // Calculate y, z using wires (need at least two hits).
        
        if(nhits >= 2) {
            
            // Calculate y and z by chisquare minimization of wire coordinates.
            
            double sw = 0.;    // sum w_i
            double sus = 0.;   // sum w_i u_i sin_th_i
            double suc = 0.;   // sum w_i u_i cos_th_i
            double sc2 = 0.;   // sum w_i cos2_th_i
            double ss2 = 0.;   // sum w_i sin2_th_i
            double ssc = 0.;   // sum w_i sin_th_i cos_th_i
            
            // Loop over points.
            
            for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin();
                ihit != hits.end(); ++ihit) {
                
                const recob::Hit& hit = **ihit;
                geo::WireID hitWireID = hit.WireID();
                const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
                
                // Calculate angle and wire coordinate in this view.
                
                double hl = wgeom.HalfL();
                double cen[3];
                double cen1[3];
                wgeom.GetCenter(cen);
                wgeom.GetCenter(cen1, hl);
                double s  = (cen1[1] - cen[1]) / hl;
                double c = (cen1[2] - cen[2]) / hl;
                double u = cen[2] * s - cen[1] * c;
                double eu = geom->WirePitch(hitWireID.Plane, hitWireID.TPC) / std::sqrt(12.);
                double w = 1. / (eu * eu);
                
                // Summations
                
                sw += w;
                sus += w*u*s;
                suc += w*u*c;
                sc2 += w*c*c;
                ss2 += w*s*s;
                ssc += w*s*c;
            }
            
            // Calculate y,z
            
            double denom = sc2 * ss2 - ssc * ssc;
            if(denom != 0.) {
                xyz[1] = (-suc * ss2 + sus * ssc) / denom;
                xyz[2] = (sus * sc2 - suc * ssc) / denom;
                errxyz[2] = ss2 / denom;
                errxyz[4] = ssc / denom;
                errxyz[5] = sc2 / denom;
            }
            
            // Make space point.
            
            recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
            sptv.push_back(spt);
        }
        return;
    }

void trkf::SpacePointAlg::fillSpacePoints	(	std::vector< recob::SpacePoint > &	spts,
		std::multimap< double, KHitTrack > const &	trackMap
	)		const

Fill a collection of space points.

Fill a collection of space points.

Arguments:

spts - Collection of space points to fill. sptalg - Space point algorithm object.

This method uses the hits contained in this track to construct space points.

This method does not have any knowledge of what constitutes a good space point, except that Hits are required to be consecutive when sorted by path distance, and space points are required to pass compatibility tests used by the space point algorithm object. This method will make space points from either two or three Hits (even for three-plane detectors), if the space point algorithm is configured to allow it.

Definition at line 610 of file SpacePointAlg.cxx.

References art::PtrVector< T >::clear(), compatible(), fillComplexSpacePoint(), fillSpacePoint(), trkf::KHitWireX::getHit(), trkf::KHitTrack::getHit(), hits(), numHitMap(), art::PtrVector< T >::push_back(), art::PtrVector< T >::size(), and track.

Referenced by trkf::Track3DKalmanHit::createOutputs(), and trkf::TrackKalmanCheater::produce().

    {
        // Loop over KHitTracks.
        
        art::PtrVector<recob::Hit> hits;
        art::PtrVector<recob::Hit> compatible_hits;
        for(std::multimap<double, KHitTrack>::const_iterator it = trackMap.begin();
            it != trackMap.end(); ++it) {
            const KHitTrack& track = (*it).second;
            
            // Extrack Hit from track.
            
            const std::shared_ptr<const KHitBase>& hit = track.getHit();
            const KHitWireX* phit = dynamic_cast<const KHitWireX*>(&*hit);
            if(phit != 0) {
                const art::Ptr<recob::Hit> prhit = phit->getHit();
                
                // Test this hit for compatibility.
                
                hits.push_back(prhit);
                bool ok = this->compatible(hits);
                if(!ok) {
                    
                    // The new hit is not compatible.  Make a space point out of
                    // the last known compatible hits, provided there are at least
                    // two.
                    
                    if(compatible_hits.size() >= 2) {
                        this->fillSpacePoint(compatible_hits, spts, this->numHitMap());
                        compatible_hits.clear();
                    }
                    
                    // Forget about any previous hits.
                    
                    hits.clear();
                    hits.push_back(prhit);
                }
                
                // Update the list of known compatible hits.
                
                compatible_hits = hits;
            }
        }
        
        // Maybe make one final space point.
        
        if(compatible_hits.size() >= 2) {
            this->fillSpacePoint(compatible_hits, spts, this->numHitMap());
        }
    }

bool trkf::SpacePointAlg::filter ( ) const

inline

Definition at line 88 of file SpacePointAlg.h.

{return fFilter;}

const art::PtrVector< recob::Hit > & trkf::SpacePointAlg::getAssociatedHits

(

const recob::SpacePoint &

spt

)

const

Definition at line 1527 of file SpacePointAlg.cxx.

References fSptHitMap, and recob::SpacePoint::ID().

Referenced by trkf::SpacePointAna::analyze(), trkf::Track3DKalmanHit::createOutputs(), trkf::SeedFinderAlgorithm::FindSeedAtEnd(), trkf::SeedFinderAlgorithm::FindSeeds(), makeSpacePoints(), trkf::TCTrack::produce(), trkf::SpacePointCheater::produce(), trkf::SpacePointFinder::produce(), and trkf::TrackKalmanCheater::produce().

    {
        // It is an error if no hits are associated with this space point (throw exception).
        
        std::map<int, art::PtrVector<recob::Hit> >::const_iterator it =
        fSptHitMap.find(spt.ID());
        if(it == fSptHitMap.end())
        {
            mf::LogWarning("SpacePointAlg") <<"Looking for ID " << spt.ID()<< " from " << fSptHitMap.size()<<std::endl;
            throw cet::exception("SpacePointAlg") << "No Hits associated with space point.\n";
        }
        return (*it).second;
        
    }

void trkf::SpacePointAlg::makeMCTruthSpacePoints	(	const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	spts
	)		const

Definition at line 856 of file SpacePointAlg.cxx.

References makeSpacePoints().

Referenced by trkf::SpacePointAna::analyze(), and trkf::SpacePointCheater::produce().

    {
        makeSpacePoints(hits, spts, true);
    }

void trkf::SpacePointAlg::makeSpacePoints	(	const art::PtrVector< recob::Hit > &	hits,
		std::vector< recob::SpacePoint > &	spts
	)		const

Definition at line 846 of file SpacePointAlg.cxx.

Referenced by trkf::SpacePointAna::analyze(), trkf::SeedFinderAlgorithm::FindSeeds(), trkf::FeatureTracker::Get3DFeaturePoints(), makeMCTruthSpacePoints(), cluster::ClusterMatchAlg::Match_SpacePoint(), trkf::TCTrack::produce(), and trkf::SpacePointFinder::produce().

    {
        makeSpacePoints(hits, spts, false);
    }

void trkf::SpacePointAlg::makeSpacePoints ( const art::PtrVector< recob::Hit > &  hits, std::vector< recob::SpacePoint > &  spts, bool  useMC  ) const

private

Todo:

Why are we still checking on whether this is MC or not?

Todo:

Such checks should not be in reconstruction code.

Definition at line 866 of file SpacePointAlg.cxx.

References art::PtrVector< T >::begin(), evd::details::begin(), c1, c2, recob::SpacePoint::Chisq(), art::PtrVector< T >::clear(), compatible(), geo::CryostatID::Cryostat, geo::GeometryCore::Cryostat(), tca::debug, trkf::SpacePointAlg::HitMCInfo::dist2, art::PtrVector< T >::end(), art::PtrVector< T >::erase(), fEnableU, fEnableV, fEnableW, fFilter, fHitMCMap, fillComplexSpacePoint(), fillSpacePoint(), fMaxDT, fMaxS, fMerge, fMinViews, fPreferColl, fSptHitMap, getAssociatedHits(), geo::WireGeo::GetCenter(), detinfo::DetectorProperties::GetXTicksOffset(), geo::WireGeo::HalfL(), cheat::BackTrackerService::HitToAvgSimIDEs(), mf::isDebugEnabled(), geo::kCollection, geo::kU, geo::kV, geo::kZ, max, min, geo::GeometryCore::Ncryostats(), geo::TPCGeo::Nplanes(), geo::CryostatGeo::NTPC(), trkf::SpacePointAlg::HitMCInfo::pchit, recob::Hit::PeakTime(), geo::TPCGeo::Plane(), geo::PlaneID::Plane, art::PtrVector< T >::push_back(), art::PtrVector< T >::reserve(), geo::GeometryCore::SignalType(), cheat::BackTrackerService::SimIDEsToXYZ(), art::PtrVector< T >::size(), t1, t2, geo::CryostatGeo::TPC(), geo::TPCID::TPC, trkf::SpacePointAlg::HitMCInfo::trackIDs, update(), recob::Hit::View(), geo::PlaneGeo::Wire(), geo::WireID::Wire, recob::Hit::WireID(), geo::GeometryCore::WireIDToWireGeo(), geo::GeometryCore::WirePitch(), x, and trkf::SpacePointAlg::HitMCInfo::xyz.

    {
        // Get services.
        
        art::ServiceHandle<geo::Geometry> geom;
        const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
        
        // Clear space point to hit map.
        
        fSptHitMap.clear();
        
        // Print diagnostic information.
        
        update();
        
        // First make sure result vector is empty.
        
        spts.erase(spts.begin(), spts.end());
        
        // Statistics.
        
        int n2 = 0;  // Number of two-hit space points.
        int n3 = 0;  // Number of three-hit space points.
        int n2filt = 0;  // Number of two-hit space points after filtering/merging.
        int n3filt = 0;  // Number of three-hit space pointe after filtering/merging.
        
        // Sort hits into maps indexed by [cryostat][tpc][plane][wire].
        // If using mc information, also generate maps of sim::IDEs and mc
        // position indexed by hit.
        
        std::vector< std::vector<std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit> > > > > hitmap;
        fHitMCMap.clear();
        
        unsigned int ncstat = geom->Ncryostats();
        hitmap.resize(ncstat);
        for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
            unsigned int ntpc = geom->Cryostat(cstat).NTPC();
            hitmap[cstat].resize(ntpc);
            for(unsigned int tpc = 0; tpc < ntpc; ++tpc) {
                int nplane = geom->Cryostat(cstat).TPC(tpc).Nplanes();
                hitmap[cstat][tpc].resize(nplane);
            }
        }
        
        for(art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end(); ++ihit) {
            const art::Ptr<recob::Hit>& phit = *ihit;
            geo::View_t view = phit->View();
            if((view == geo::kU && fEnableU) ||
               (view == geo::kV && fEnableV) ||
               (view == geo::kZ && fEnableW)) {
                geo::WireID phitWireID = phit->WireID();
                hitmap[phitWireID.Cryostat][phitWireID.TPC][phitWireID.Plane].insert(std::make_pair(phitWireID.Wire, phit));
            }
        }
        
        // Fill mc information, including IDEs and closest neighbors
        // of each hit.
        if(useMC) {
            art::ServiceHandle<cheat::BackTrackerService> bt_serv;
            
            // First loop over hits and fill track ids and mc position.
            for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
                for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
                    int nplane = geom->Cryostat(cstat).TPC(tpc).Nplanes();
                    for(int plane = 0; plane < nplane; ++plane) {
                        for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator ihit = hitmap[cstat][tpc][plane].begin();
                            ihit != hitmap[cstat][tpc][plane].end(); ++ihit) {
                            const art::Ptr<recob::Hit>& phit = ihit->second;
                            const recob::Hit& hit = *phit;
                            HitMCInfo& mcinfo = fHitMCMap[&hit];   // Default HitMCInfo.
                            
                            // Fill default nearest neighbor information (i.e. none).
                            
                            mcinfo.pchit.resize(nplane, 0);
                            mcinfo.dist2.resize(nplane, 1.e20);
                            
                            // Get sim::IDEs for this hit.
                            
                            std::vector<sim::IDE> ides = bt_serv->HitToAvgSimIDEs(phit);
                            
                            // Get sorted track ids. for this hit.
                            
                            mcinfo.trackIDs.reserve(ides.size());
                            for(std::vector<sim::IDE>::const_iterator i = ides.begin();
                                i != ides.end(); ++i)
                                mcinfo.trackIDs.push_back(i->trackID);
                            sort(mcinfo.trackIDs.begin(), mcinfo.trackIDs.end());
                            
                            // Get position of ionization for this hit.
                            
                            try {
                                mcinfo.xyz = bt_serv->SimIDEsToXYZ(ides);
                            }
                            catch(cet::exception& x) {
                                mcinfo.xyz.clear();
                            }
                        } // end loop over ihit
                    }// end loop oer planes
                }// end loop over TPCs
            }// end loop over cryostats
            
            // Loop over hits again and fill nearest neighbor information for real.
            for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
                for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
                    int nplane = geom->Cryostat(cstat).TPC(tpc).Nplanes();
                    for(int plane = 0; plane < nplane; ++plane) {
                        for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator ihit = hitmap[cstat][tpc][plane].begin();
                            ihit != hitmap[cstat][tpc][plane].end(); ++ihit) {
                            const art::Ptr<recob::Hit>& phit = ihit->second;
                            const recob::Hit& hit = *phit;
                            HitMCInfo& mcinfo = fHitMCMap[&hit];
                            if(mcinfo.xyz.size() != 0) {
                                assert(mcinfo.xyz.size() == 3);
                                
                                // Fill nearest neighbor information for this hit.
                                
                                for(int plane2 = 0; plane2 < nplane; ++plane2) {
                                    for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator jhit = hitmap[cstat][tpc][plane2].begin();
                                        jhit != hitmap[cstat][tpc][plane2].end(); ++jhit) {
                                        const art::Ptr<recob::Hit>& phit2 = jhit->second;
                                        const recob::Hit& hit2 = *phit2;
                                        const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
                                        
                                        
                                        if(mcinfo2.xyz.size() != 0) {
                                            assert(mcinfo2.xyz.size() == 3);
                                            double dx = mcinfo.xyz[0] - mcinfo2.xyz[0];
                                            double dy = mcinfo.xyz[1] - mcinfo2.xyz[1];
                                            double dz = mcinfo.xyz[2] - mcinfo2.xyz[2];
                                            double dist2 = dx*dx + dy*dy + dz*dz;
                                            if(dist2 < mcinfo.dist2[plane2]) {
                                                mcinfo.dist2[plane2] = dist2;
                                                mcinfo.pchit[plane2] = &hit2;
                                            }
                                        }// end if mcinfo2.xyz valid
                                    }// end loop over jhit
                                }// end loop over plane2
                            }// end if mcinfo.xyz valid.
                        }// end loop over ihit
                    }// end loop over plane
                }// end loop over tpc
            }// end loop over cryostats
        }// end if MC
        
        // use mf::LogDebug instead of LOG_DEBUG because we reuse it in many lines
        // insertions are protected by mf::isDebugEnabled()
        mf::LogDebug debug("SpacePointAlg");
        if (mf::isDebugEnabled()) {
            debug << "Total hits = " << hits.size() << "\n\n";
            
            for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
                for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
                    int nplane = hitmap[cstat][tpc].size();
                    for(int plane = 0; plane < nplane; ++plane) {
                        debug << "TPC, Plane: " << tpc << ", " << plane
                        << ", hits = " << hitmap[cstat][tpc][plane].size() << "\n";
                    }
                }
            } // end loop over cryostats
        } // if debug
        
        // Make empty multimap from hit pointer on preferred
        // (most-populated or collection) plane to space points that
        // include that hit (used for sorting, filtering, and
        // merging).
        
        typedef const recob::Hit* sptkey_type;
        std::multimap<sptkey_type, recob::SpacePoint> sptmap;
        std::set<sptkey_type> sptkeys;              // Keys of multimap.
        
        // Loop over TPCs.
        for(unsigned int cstat = 0; cstat < ncstat; ++cstat){
            for(unsigned int tpc = 0; tpc < geom->Cryostat(cstat).NTPC(); ++tpc) {
                
                geo::TPCID tpcid(cstat, tpc);
                
                // Sort maps in increasing order of number of hits.
                // This is so that we can do the outer loops over hits
                // over the views with fewer hits.
                //
                // If config parameter PreferColl is true, treat the colleciton
                // plane as if it had the most hits, regardless of how many
                // hits it actually has.  This will force space points to be
                // filtered and merged with respect to the collection plane
                // wires.  It will also force space points to be sorted by
                // collection plane wire.
                
                int nplane = hitmap[cstat][tpc].size();
                std::vector<int> index(nplane);
                
                for(int i=0; i<nplane; ++i)
                    index[i] = i;
                
                for(int i=0; i<nplane-1; ++i) {
                    
                    for(int j=i+1; j<nplane; ++j) {
                        bool icoll = fPreferColl &&
                        geom->SignalType(geo::PlaneID(tpcid, index[i])) == geo::kCollection;
                        bool jcoll = fPreferColl &&
                        geom->SignalType(geo::PlaneID(tpcid, index[j])) == geo::kCollection;
                        if((hitmap[cstat][tpc][index[i]].size() > hitmap[cstat][tpc][index[j]].size() &&
                            !jcoll) || icoll) {
                            int temp = index[i];
                            index[i] = index[j];
                            index[j] = temp;
                        }
                    }
                }// end loop over i
                
                // how many views with hits?
                // This will allow for the special case where we might have only 2 planes of information and
                // still want space points even if a three plane TPC
                std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>& hitsByPlaneVec = hitmap[cstat][tpc];
                int nViewsWithHits(0);
                
                for(int i = 0; i < nplane; i++)
                {
                    if (hitsByPlaneVec[index[i]].size() > 0) nViewsWithHits++;
                }
                
                // If two-view space points are allowed, make a double loop
                // over hits and produce space points for compatible hit-pairs.
                
                if((nViewsWithHits == 2 || nplane == 2) && fMinViews <= 2) {
                    
                    // Loop over pairs of views.
                    for(int i=0; i<nplane-1; ++i) {
                        unsigned int plane1 = index[i];
                        
                        if (hitmap[cstat][tpc][plane1].empty()) continue;
                        
                        for(int j=i+1; j<nplane; ++j) {
                            unsigned int plane2 = index[j];
                            
                            if (hitmap[cstat][tpc][plane2].empty()) continue;
                            
                            // Get angle, pitch, and offset of plane2 wires.
                            const geo::WireGeo& wgeo2 = geom->Cryostat(cstat).TPC(tpc).Plane(plane2).Wire(0);
                            double hl2 = wgeo2.HalfL();
                            double xyz21[3];
                            double xyz22[3];
                            wgeo2.GetCenter(xyz21, -hl2);
                            wgeo2.GetCenter(xyz22, hl2);
                            double s2 = (xyz22[1] - xyz21[1]) / (2.*hl2);
                            double c2 = (xyz22[2] - xyz21[2]) / (2.*hl2);
                            double dist2 = -xyz21[1] * c2 + xyz21[2] * s2;
                            double pitch2 = geom->WirePitch(plane2, tpc, cstat);
                            
                            if(!fPreferColl && hitmap[cstat][tpc][plane1].size() > hitmap[cstat][tpc][plane2].size())
                                throw cet::exception("SpacePointAlg") << "makeSpacePoints(): hitmaps with incompatible size\n";
                            
                            
                            // Loop over pairs of hits.
                            
                            art::PtrVector<recob::Hit> hitvec;
                            hitvec.reserve(2);
                            
                            for(std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
                                ihit1 = hitmap[cstat][tpc][plane1].begin();
                                ihit1 != hitmap[cstat][tpc][plane1].end(); ++ihit1) {
                                
                                const art::Ptr<recob::Hit>& phit1 = ihit1->second;
                                geo::WireID phit1WireID = phit1->WireID();
                                const geo::WireGeo& wgeo = geom->WireIDToWireGeo(phit1WireID);
                                
                                // Get endpoint coordinates of this wire.
                                // (kept as assertions for performance reasons)
                                assert(phit1WireID.Cryostat == cstat);
                                assert(phit1WireID.TPC == tpc);
                                assert(phit1WireID.Plane == plane1);
                                double hl1 = wgeo.HalfL();
                                double xyz1[3];
                                double xyz2[3];
                                wgeo.GetCenter(xyz1, -hl1);
                                wgeo.GetCenter(xyz2, hl1);
                                
                                // Find the plane2 wire numbers corresponding to the endpoints.
                                
                                double wire21 = (-xyz1[1] * c2 + xyz1[2] * s2 - dist2) / pitch2;
                                double wire22 = (-xyz2[1] * c2 + xyz2[2] * s2 - dist2) / pitch2;
                                
                                int wmin = std::max(0., std::min(wire21, wire22));
                                int wmax = std::max(0., std::max(wire21, wire22) + 1.);
                                
                                std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
                                ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
                                ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
                                
                                for(; ihit2 != ihit2end; ++ihit2) {
                                    
                                    const art::Ptr<recob::Hit>& phit2 = ihit2->second;
                                    
                                    // Check current pair of hits for compatibility.
                                    // By construction, hits should always have compatible views
                                    // and times, but may not have compatible mc information.
                                    
                                    hitvec.clear();
                                    hitvec.push_back(phit1);
                                    hitvec.push_back(phit2);
                                    bool ok = compatible(hitvec,  useMC);
                                    if(ok) {
                                        
                                        // Add a space point.
                                        
                                        ++n2;
                                        
                                        // make a dummy vector of recob::SpacePoints
                                        // as we are filtering or merging and don't want to
                                        // add the created SpacePoint to the final collection just yet
                                        // This dummy vector will hold just one recob::SpacePoint,
                                        // which will go into the multimap and then the vector
                                        // will go out of scope.
                                        
                                        std::vector<recob::SpacePoint> sptv;
                                        fillSpacePoint(hitvec, sptv, sptmap.size());
                                        sptkey_type key = &*phit2;
                                        sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
                                        sptkeys.insert(key);
                                    }
                                }
                            }
                        }
                    }
                }// end if fMinViews <= 2
                
                // If three-view space points are allowed, make a triple loop
                // over hits and produce space points for compatible triplets.
                
                if(nplane >= 3 && fMinViews <= 3) {
                    
                    // Loop over triplets of hits.
                    
                    art::PtrVector<recob::Hit> hitvec;
                    hitvec.reserve(3);
                    
                    unsigned int plane1 = index[0];
                    unsigned int plane2 = index[1];
                    unsigned int plane3 = index[2];
                    
                    // Get angle, pitch, and offset of plane1 wires.
                    
                    const geo::WireGeo& wgeo1 = geom->Cryostat(cstat).TPC(tpc).Plane(plane1).Wire(0);
                    double hl1 = wgeo1.HalfL();
                    double xyz11[3];
                    double xyz12[3];
                    wgeo1.GetCenter(xyz11, -hl1);
                    wgeo1.GetCenter(xyz12, hl1);
                    double s1 = (xyz12[1] - xyz11[1]) / (2.*hl1);
                    double c1 = (xyz12[2] - xyz11[2]) / (2.*hl1);
                    double dist1 = -xyz11[1] * c1 + xyz11[2] * s1;
                    double pitch1 = geom->WirePitch(plane1, tpc, cstat);
                    const double TicksOffset1 = detprop->GetXTicksOffset(plane1,tpc,cstat);
                    
                    // Get angle, pitch, and offset of plane2 wires.
                    
                    const geo::WireGeo& wgeo2 = geom->Cryostat(cstat).TPC(tpc).Plane(plane2).Wire(0);
                    double hl2 = wgeo2.HalfL();
                    double xyz21[3];
                    double xyz22[3];
                    wgeo2.GetCenter(xyz21, -hl2);
                    wgeo2.GetCenter(xyz22, hl2);
                    double s2 = (xyz22[1] - xyz21[1]) / (2.*hl2);
                    double c2 = (xyz22[2] - xyz21[2]) / (2.*hl2);
                    double dist2 = -xyz21[1] * c2 + xyz21[2] * s2;
                    double pitch2 = geom->WirePitch(plane2, tpc, cstat);
                    const double TicksOffset2 = detprop->GetXTicksOffset(plane2,tpc,cstat);
                    
                    // Get angle, pitch, and offset of plane3 wires.
                    
                    const geo::WireGeo& wgeo3 = geom->Cryostat(cstat).TPC(tpc).Plane(plane3).Wire(0);
                    double hl3 = wgeo3.HalfL();
                    double xyz31[3];
                    double xyz32[3];
                    wgeo3.GetCenter(xyz31, -hl3);
                    wgeo3.GetCenter(xyz32, hl3);
                    double s3 = (xyz32[1] - xyz31[1]) / (2.*hl3);
                    double c3 = (xyz32[2] - xyz31[2]) / (2.*hl3);
                    double dist3 = -xyz31[1] * c3 + xyz31[2] * s3;
                    double pitch3 = geom->WirePitch(plane3, tpc, cstat);
                    const double TicksOffset3 = detprop->GetXTicksOffset(plane3,tpc,cstat);
                    
                    // Get sine of angle differences.
                    
                    double s12 = s1 * c2 - s2 * c1;   // sin(theta1 - theta2).
                    double s23 = s2 * c3 - s3 * c2;   // sin(theta2 - theta3).
                    double s31 = s3 * c1 - s1 * c3;   // sin(theta3 - theta1).
                    
                    // Loop over hits in plane1.
                    
                    std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
                    ihit1 = hitmap[cstat][tpc][plane1].begin(),
                    ihit1end = hitmap[cstat][tpc][plane1].end();
                    for(; ihit1 != ihit1end; ++ihit1) {
                        
                        unsigned int wire1 = ihit1->first;
                        const art::Ptr<recob::Hit>& phit1 = ihit1->second;
                        geo::WireID phit1WireID = phit1->WireID();
                        const geo::WireGeo& wgeo = geom->WireIDToWireGeo(phit1WireID);
                        
                        // Get endpoint coordinates of this wire from plane1.
                        // (kept as assertions for performance reasons)
                        assert(phit1WireID.Cryostat == cstat);
                        assert(phit1WireID.TPC == tpc);
                        assert(phit1WireID.Plane == plane1);
                        assert(phit1WireID.Wire == wire1);
                        double hl1 = wgeo.HalfL();
                        double xyz1[3];
                        double xyz2[3];
                        wgeo.GetCenter(xyz1, -hl1);
                        wgeo.GetCenter(xyz2, hl1);
                        
                        // Get corrected time and oblique coordinate of first hit.
                        
                        double t1 = phit1->PeakTime() - TicksOffset1;
                        double u1 = wire1 * pitch1 + dist1;
                        
                        // Find the plane2 wire numbers corresponding to the endpoints.
                        
                        double wire21 = (-xyz1[1] * c2 + xyz1[2] * s2 - dist2) / pitch2;
                        double wire22 = (-xyz2[1] * c2 + xyz2[2] * s2 - dist2) / pitch2;
                        
                        int wmin = std::max(0., std::min(wire21, wire22));
                        int wmax = std::max(0., std::max(wire21, wire22) + 1.);
                        
                        std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
                        ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
                        ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
                        
                        for(; ihit2 != ihit2end; ++ihit2) {
                            
                            int wire2 = ihit2->first;
                            const art::Ptr<recob::Hit>& phit2 = ihit2->second;
                            
                            // Get corrected time of second hit.
                            
                            double t2 = phit2->PeakTime() - TicksOffset2;
                            
                            // Check maximum time difference with first hit.
                            
                            bool dt12ok = std::abs(t1-t2) <= fMaxDT;
                            if(dt12ok) {
                                
                                // Test first two hits for compatibility before looping
                                // over third hit.
                                
                                hitvec.clear();
                                hitvec.push_back(phit1);
                                hitvec.push_back(phit2);
                                bool h12ok = compatible(hitvec, useMC);
                                if(h12ok) {
                                    
                                    // Get oblique coordinate of second hit.
                                    
                                    double u2 = wire2 * pitch2 + dist2;
                                    
                                    // Predict plane3 oblique coordinate and wire number.
                                    
                                    double u3pred = (-u1*s23 - u2*s31) / s12;
                                    double w3pred = (u3pred - dist3) / pitch3;
                                    double w3delta = std::abs(fMaxS / (s12 * pitch3));
                                    int w3min = std::max(0., std::ceil(w3pred - w3delta));
                                    int w3max = std::max(0., std::floor(w3pred + w3delta));
                                    
                                    std::map<unsigned int, art::Ptr<recob::Hit> >::const_iterator
                                    ihit3 = hitmap[cstat][tpc][plane3].lower_bound(w3min),
                                    ihit3end = hitmap[cstat][tpc][plane3].upper_bound(w3max);
                                    
                                    for(; ihit3 != ihit3end; ++ihit3) {
                                        
                                        int wire3 = ihit3->first;
                                        const art::Ptr<recob::Hit>& phit3 = ihit3->second;
                                        
                                        // Get corrected time of third hit.
                                        
                                        double t3 = phit3->PeakTime() - TicksOffset3;
                                        
                                        // Check time difference of third hit compared to first two hits.
                                        
                                        bool dt123ok = std::abs(t1-t3) <= fMaxDT && std::abs(t2-t3) <= fMaxDT;
                                        if(dt123ok) {
                                            
                                            // Get oblique coordinate of third hit and check spatial separation.
                                            
                                            double u3 = wire3 * pitch3 + dist3;
                                            double S = s23 * u1 + s31 * u2 + s12 * u3;
                                            bool sok = std::abs(S) <= fMaxS;
                                            if(sok) {
                                                
                                                // Test triplet for compatibility.
                                                
                                                hitvec.clear();
                                                hitvec.push_back(phit1);
                                                hitvec.push_back(phit2);
                                                hitvec.push_back(phit3);
                                                bool h123ok = compatible(hitvec,  useMC);
                                                if(h123ok) {
                                                    
                                                    // Add a space point.
                                                    
                                                    ++n3;
                                                    
                                                    // make a dummy vector of recob::SpacePoints
                                                    // as we are filtering or merging and don't want to
                                                    // add the created SpacePoint to the final collection just yet
                                                    // This dummy vector will hold just one recob::SpacePoint,
                                                    // which will go into the multimap and then the vector
                                                    // will go out of scope.
                                                    
                                                    std::vector<recob::SpacePoint> sptv;
                                                    fillSpacePoint(hitvec, sptv, sptmap.size()-1);
                                                    sptkey_type key = &*phit3;
                                                    sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
                                                    sptkeys.insert(key);
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }// end if fMinViews <= 3
                
                // Do Filtering.
                
                if(fFilter) {
                    
                    // Transfer (some) space points from sptmap to spts.
                    
                    spts.reserve(spts.size() + sptkeys.size());
                    
                    // Loop over keys of space point map.
                    // Space points that have the same key are candidates for filtering.
                    
                    for(std::set<sptkey_type>::const_iterator i = sptkeys.begin();
                        i != sptkeys.end(); ++i) {
                        sptkey_type key = *i;
                        
                        // Loop over space points corresponding to the current key.
                        // Choose the single best space point from among this group.
                        
                        double best_chisq = 0.;
                        const recob::SpacePoint* best_spt = 0;
                        
                        for(std::multimap<sptkey_type, recob::SpacePoint>::const_iterator j = sptmap.lower_bound(key);
                            j != sptmap.upper_bound(key); ++j) {
                            const recob::SpacePoint& spt = j->second;
                            if(best_spt == 0 || spt.Chisq() < best_chisq) {
                                best_spt = &spt;
                                best_chisq = spt.Chisq();
                            }
                        }
                        
                        // Transfer best filtered space point to result vector.
                        
                        if (!best_spt)
                            throw cet::exception("SpacePointAlg") << "makeSpacePoints(): no best point\n";
                        spts.push_back(*best_spt);
                        if(fMinViews <= 2)
                            ++n2filt;
                        else
                            ++n3filt;
                    }
                }// end if filtering
                
                // Do merging.
                
                else if(fMerge) {
                    
                    // Transfer merged space points from sptmap to spts.
                    
                    spts.reserve(spts.size() + sptkeys.size());
                    
                    // Loop over keys of space point map.
                    // Space points that have the same key are candidates for merging.
                    
                    for(std::set<sptkey_type>::const_iterator i = sptkeys.begin();
                        i != sptkeys.end(); ++i) {
                        sptkey_type key = *i;
                        
                        // Loop over space points corresponding to the current key.
                        // Make a collection of hits that is the union of the hits
                        // from each candidate space point.
                        
                        std::multimap<sptkey_type, recob::SpacePoint>::const_iterator
                        jSPT = sptmap.lower_bound(key), jSPTend = sptmap.upper_bound(key);
                        
                        art::PtrVector<recob::Hit> merged_hits;
                        for(; jSPT != jSPTend; ++jSPT) {
                            const recob::SpacePoint& spt = jSPT->second;
                            
                            // Loop over hits from this space points.
                            // Add each hit to the collection of all hits.
                            
                            const art::PtrVector<recob::Hit>& spt_hits = getAssociatedHits(spt);
                            merged_hits.reserve(merged_hits.size() + spt_hits.size()); // better than nothing, but not ideal
                            for(art::PtrVector<recob::Hit>::const_iterator k = spt_hits.begin();
                                k != spt_hits.end(); ++k) {
                                const art::Ptr<recob::Hit>& hit = *k;
                                merged_hits.push_back(hit);
                            }
                        }
                        
                        // Remove duplicates.
                        
                        std::sort(merged_hits.begin(), merged_hits.end());
                        art::PtrVector<recob::Hit>::iterator it = 
                        std::unique(merged_hits.begin(), merged_hits.end());
                        merged_hits.erase(it, merged_hits.end());
                        
                        // Construct a complex space points using merged hits.
                        
                        fillComplexSpacePoint(merged_hits, spts, sptmap.size() + spts.size()-1);
                        
                        if(fMinViews <= 2)
                            ++n2filt;
                        else
                            ++n3filt;
                    }
                }// end if merging
                
                // No filter, no merge.
                
                else {
                    
                    // Transfer all space points from sptmap to spts.
                    
                    spts.reserve(spts.size() + sptkeys.size());
                    
                    // Loop over space points.
                    
                    for(std::multimap<sptkey_type, recob::SpacePoint>::const_iterator j = sptmap.begin();
                        j != sptmap.end(); ++j) {
                        const recob::SpacePoint& spt = j->second;
                        spts.push_back(spt);
                    }
                    
                    // Update statistics.
                    
                    n2filt = n2;
                    n3filt = n3;
                }
            }// end loop over tpcs
        }// end loop over cryostats
        
        if (mf::isDebugEnabled()) {
            debug << "\n2-hit space points = " << n2 << "\n"
            << "3-hit space points = " << n3 << "\n"
            << "2-hit filtered/merged space points = " << n2filt << "\n"
            << "3-hit filtered/merged space points = " << n3filt;
        } // if debug
    }

double trkf::SpacePointAlg::maxDT ( ) const

inline

Definition at line 90 of file SpacePointAlg.h.

Referenced by trkf::BezierTrackerAlgorithm::MakeTracks(), and cluster::ClusterMatchAlg::Match_SpacePoint().

{return fMaxDT;}

double trkf::SpacePointAlg::maxS ( ) const

inline

Definition at line 91 of file SpacePointAlg.h.

{return fMaxS;}

bool trkf::SpacePointAlg::merge ( ) const

inline

Definition at line 89 of file SpacePointAlg.h.

Referenced by trkf::SpacePointAna::analyze(), and trkf::SpacePointAna::bookHistograms().

{return fMerge;}

int trkf::SpacePointAlg::minViews ( ) const

inline

Definition at line 92 of file SpacePointAlg.h.

Referenced by trkf::SpacePointCheater::produce(), and trkf::SpacePointFinder::produce().

{return fMinViews;}

int trkf::SpacePointAlg::numHitMap ( ) const

inline

Definition at line 148 of file SpacePointAlg.h.

Referenced by fillSpacePoints().

{return fSptHitMap.size();}

void trkf::SpacePointAlg::reconfigure

(

const fhicl::ParameterSet &

pset

)

Definition at line 67 of file SpacePointAlg.cxx.

References fEnableU, fEnableV, fEnableW, fFilter, fMaxDT, fMaxS, fMerge, fMinViews, fPreferColl, and fhicl::ParameterSet::get().

Referenced by trkf::TCTrack::reconfigure(), trkf::SpacePointCheater::reconfigure(), trkf::SpacePointFinder::reconfigure(), trkf::TrackKalmanCheater::reconfigure(), trkf::Track3DKalmanHit::reconfigure(), and SpacePointAlg().

    {
        // Get configuration parameters.
        
        fMaxDT = pset.get<double>("MaxDT");
        fMaxS = pset.get<double>("MaxS");
        
        fMinViews = pset.get<int>("MinViews");
        
        fEnableU = pset.get<bool>("EnableU");
        fEnableV = pset.get<bool>("EnableV");
        fEnableW = pset.get<bool>("EnableW");
        fFilter = pset.get<bool>("Filter");
        fMerge = pset.get<bool>("Merge");
        fPreferColl = pset.get<bool>("PreferColl");
        
        // Only allow one of fFilter and fMerge to be true.
        
        if(fFilter && fMerge)
            throw cet::exception("SpacePointAlg") << "Filter and Merge flags are both true.\n";
        
        // Report.
        
        mf::LogInfo("SpacePointAlg")
        << "SpacePointAlg configured with the following parameters:\n"
        << "  MaxDT = " << fMaxDT << "\n"
        << "  MaxS = " << fMaxS << "\n"
        << "  MinViews = " << fMinViews << "\n"
        << "  EnableU = " << fEnableU << "\n"
        << "  EnableV = " << fEnableV << "\n"
        << "  EnableW = " << fEnableW << "\n"
        << "  Filter = " << fFilter << "\n"
        << "  Merge = " << fMerge;
    }

double trkf::SpacePointAlg::separation

(

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

hits

)

const

Definition at line 225 of file SpacePointAlg.cxx.

References geo::CryostatID::Cryostat, geo::WireGeo::GetCenter(), geo::WireGeo::HalfL(), geo::PlaneID::Plane, s, art::PtrVector< T >::size(), geo::TPCID::TPC, recob::Hit::WireID(), and geo::GeometryCore::WireIDToWireGeo().

Referenced by trkf::SpacePointAna::analyze().

    {
        // Get geometry service.
        
        art::ServiceHandle<geo::Geometry> geom;
        
        // Trivial case - fewer than three hits.
        
        if(hits.size() < 3)
            return 0.;
        
        // Error case - more than three hits.
        
        if(hits.size() > 3) {
            mf::LogError("SpacePointAlg") << "Method separation called with more than three htis.";
            return 0.;
        }
        
        // Got exactly three hits.
        
        // Calculate angles and distance of each hit from origin.
        
        double dist[3] = {0., 0., 0.};
        double sinth[3] = {0., 0., 0.};
        double costh[3] = {0., 0., 0.};
        unsigned int cstats[3];
        unsigned int tpcs[3];
        unsigned int planes[3];
        
        for(int i=0; i<3; ++i) {
            
            // Get tpc, plane, wire.
            
            const recob::Hit& hit = *(hits[i]);
            const geo::WireGeo& wgeom = geom->WireIDToWireGeo(hit.WireID());
            cstats[i] = hit.WireID().Cryostat;
            tpcs[i] = hit.WireID().TPC;
            planes[i] = hit.WireID().Plane;
            
            // Check tpc and plane errors.
            
            for(int j=0; j<i; ++j) {
                
                if(cstats[j] != hit.WireID().Cryostat) {
                    mf::LogError("SpacePointAlg") << "Method separation called with hits from multiple cryostats..";
                    return 0.;
                }
                
                if(tpcs[j] != hit.WireID().TPC) {
                    mf::LogError("SpacePointAlg") << "Method separation called with hits from multiple tpcs..";
                    return 0.;
                }
                
                if(planes[j] == hit.WireID().Plane ) {
                    mf::LogError("SpacePointAlg") << "Method separation called with hits from the same plane..";
                    return 0.;
                }
            }
            
            // Get angles and distance of wire.
            
            double hl = wgeom.HalfL();
            double xyz[3];
            double xyz1[3];
            wgeom.GetCenter(xyz);
            wgeom.GetCenter(xyz1, hl);
            double s = (xyz1[1] - xyz[1]) / hl;
            double c = (xyz1[2] - xyz[2]) / hl;
            sinth[hit.WireID().Plane] = s;
            costh[hit.WireID().Plane] = c;
            dist[hit.WireID().Plane] = xyz[2] * s - xyz[1] * c;
        }
        
        double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0]
                    +(sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1]
                    +(sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
        return S;
    }

void trkf::SpacePointAlg::update

(

)

const

Definition at line 105 of file SpacePointAlg.cxx.

References geo::GeometryCore::Cryostat(), detinfo::DetectorProperties::GetXTicksOffset(), geo::kCollection, geo::kHorizontal, geo::kInduction, geo::kU, geo::kV, geo::kVertical, geo::kZ, geo::GeometryCore::Ncryostats(), geo::TPCGeo::Nplanes(), geo::CryostatGeo::NTPC(), geo::PlaneGeo::Orientation(), geo::TPCGeo::Plane(), geo::TPCGeo::Plane0Pitch(), geo::TPCGeo::PlaneLocation(), geo::GeometryCore::SignalType(), geo::WireGeo::ThetaZ(), geo::CryostatGeo::TPC(), geo::PlaneGeo::View(), geo::PlaneGeo::Wire(), and geo::TPCGeo::WirePitch().

Referenced by makeSpacePoints().

    {
        // Generate info report on first call only.
        
        static bool first = true;
        bool report = first;
        first = false;
        
        // Get services.
        
        art::ServiceHandle<geo::Geometry> geom;
        const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
        
        // Calculate and print geometry information.
        
        if(report)
            mf::LogInfo("SpacePointAlg") << "Updating geometry constants.\n";
        
        
        for(unsigned int cstat = 0; cstat < geom->Ncryostats(); ++cstat){
            
            // Loop over TPCs.
            
            unsigned int const ntpc = geom->Cryostat(cstat).NTPC();
            
            for(unsigned int tpc = 0; tpc < ntpc; ++tpc) {
                const geo::TPCGeo& tpcgeom = geom->Cryostat(cstat).TPC(tpc);
                
                // Loop over planes.
                
                unsigned int const nplane = tpcgeom.Nplanes();
                
                for(unsigned int plane = 0; plane < nplane; ++plane) {
                    geo::PlaneID planeid(cstat, tpc, plane);
                    const geo::PlaneGeo& pgeom = tpcgeom.Plane(planeid);
                    
                    // Fill view-dependent quantities.
                    
                    geo::View_t view = pgeom.View();
                    std::string viewname = "?";
                    if(view == geo::kU) {
                        viewname = "U";
                    }
                    else if(view == geo::kV) {
                        viewname = "V";
                    }
                    else if(view == geo::kZ) {
                        viewname = "Z";
                    }
                    else
                        throw cet::exception("SpacePointAlg") << "Bad view = "
                        << view << "\n";
                    
                    std::string sigtypename = "?";
                    geo::SigType_t sigtype = geom->SignalType(planeid);
                    if(sigtype == geo::kInduction)
                        sigtypename = "Induction";
                    else if(sigtype == geo::kCollection)
                        sigtypename = "Collection";
                    else
                        throw cet::exception("SpacePointAlg") << "Bad signal type = "
                        << sigtype << "\n";
                    
                    std::string orientname = "?";
                    geo::Orient_t orient = pgeom.Orientation();
                    if(orient == geo::kVertical)
                        orientname = "Vertical";
                    else if(orient == geo::kHorizontal)
                        orientname = "Horizontal";
                    else
                        throw cet::exception("SpacePointAlg") << "Bad orientation = "
                        << orient << "\n";
                    
                    if(report) {
                        const double* xyz = tpcgeom.PlaneLocation(plane);
                        mf::LogInfo("SpacePointAlg") << "\nCryostat, TPC, Plane: "
                        << cstat << ","
                        << tpc << ", "
                        << plane << "\n"
                        << "  View: " << viewname << "\n"
                        << "  SignalType: " << sigtypename << "\n"
                        << "  Orientation: " << orientname << "\n"
                        << "  Plane location: " << xyz[0] << "\n"
                        << "  Plane pitch: "
                        << tpcgeom.Plane0Pitch(plane) << "\n"
                        << "  Wire angle: "
                        << tpcgeom.Plane(plane).Wire(0).ThetaZ() << "\n"
                        << "  Wire pitch: " << tpcgeom.WirePitch() << "\n"
                        << "  Time offset: "
                        << detprop->GetXTicksOffset(plane,tpc,cstat) << "\n";
                    }
                    
                    if(orient != geo::kVertical)
                        throw cet::exception("SpacePointAlg")
                        << "Horizontal wire geometry not implemented.\n";
                }// end loop over planes
            }// end loop over tpcs
        }// end loop over cryostats
    }

Member Data Documentation

bool trkf::SpacePointAlg::fEnableU

private

Enable flag (U).

Definition at line 163 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

bool trkf::SpacePointAlg::fEnableV

private

Enable flag (V).

Definition at line 164 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

bool trkf::SpacePointAlg::fEnableW

private

Enable flag (W).

Definition at line 165 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

bool trkf::SpacePointAlg::fFilter

private

Filter flag.

Definition at line 166 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

std::map<const recob::Hit*, HitMCInfo> trkf::SpacePointAlg::fHitMCMap

mutableprivate

Definition at line 179 of file SpacePointAlg.h.

Referenced by compatible(), and makeSpacePoints().

double trkf::SpacePointAlg::fMaxDT

private

Maximum time difference between planes.

Definition at line 160 of file SpacePointAlg.h.

Referenced by compatible(), makeSpacePoints(), and reconfigure().

double trkf::SpacePointAlg::fMaxS

private

Maximum space separation between wires.

Definition at line 161 of file SpacePointAlg.h.

Referenced by compatible(), makeSpacePoints(), and reconfigure().

bool trkf::SpacePointAlg::fMerge

private

Merge flag.

Definition at line 167 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

int trkf::SpacePointAlg::fMinViews

private

Mininum number of views per space point.

Definition at line 162 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

bool trkf::SpacePointAlg::fPreferColl

private

Sort by collection wire.

Definition at line 168 of file SpacePointAlg.h.

Referenced by makeSpacePoints(), and reconfigure().

std::map<int, art::PtrVector<recob::Hit> > trkf::SpacePointAlg::fSptHitMap

mutableprivate

Definition at line 180 of file SpacePointAlg.h.

Referenced by fillComplexSpacePoint(), fillSpacePoint(), getAssociatedHits(), and makeSpacePoints().

---

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

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