Inheritance diagram for trkf::CCTrackMaker:

Classes
struct	clPar

struct	ClsChainPar

struct	MatchPars

struct	TrkPar

struct	vtxPar

Public Types
using	ModuleType = EDProducer

template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Member Functions
	CCTrackMaker (fhicl::ParameterSet const &pset)

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

void	fillProductDescriptions ()

void	registerProducts (ProductDescriptions &productsToRegister)

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

std::unique_ptr< Worker >	makeWorker (WorkerParams const &wp)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Protected Member Functions
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Private Member Functions
void	produce (art::Event &evt) override

void	PrintClusters (detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid) const

void	PrintTracks () const

void	MakeClusterChains (detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, geo::TPCID const &tpcid)

float	dXClTraj (art::FindManyP< recob::Hit > const &fmCluHits, unsigned short ipl, unsigned short icl1, unsigned short end1)

void	FillChgNear (detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)

void	FillWireHitRange (geo::TPCID const &tpcid)

void	FindMaybeVertices (geo::TPCID const &tpcid)

void	VtxMatch (detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, geo::TPCID const &tpcid)

void	PlnMatch (detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)

void	AngMatch (art::FindManyP< recob::Hit > const &fmCluHits)

void	MakeFamily (geo::TPCID const &tpcid)

void	TagCosmics (geo::TPCID const &tpcid)

void	FitVertices (detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)

void	FillEndMatch (detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid, MatchPars &match)

void	FillEndMatch2 (MatchPars &match)

float	ChargeAsym (std::array< float, 3 > &mChg)

bool	FindMissingCluster (unsigned short kpl, short &kcl, unsigned short &kend, float kWir, float kX, float okWir, float okX)

bool	DupMatch (MatchPars &match, unsigned short nplanes)

void	SortMatches (detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, unsigned short procCode, geo::TPCID const &tpcid)

void	FillTrkHits (art::FindManyP< recob::Hit > const &fmCluHits, unsigned short imat, unsigned short nplanes)

void	StoreTrack (detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, unsigned short imat, unsigned short procCode, geo::TPCID const &tpcid)

float	ChargeNear (geo::PlaneID const &planeid, unsigned short wire1, float time1, unsigned short wire2, float time2)

float	AngleFactor (float slope)

Private Attributes
std::string	fHitModuleLabel

std::string	fClusterModuleLabel

std::string	fVertexModuleLabel

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

geo::WireReadoutGeom const &	wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get()

TrackTrajectoryAlg	fTrackTrajectoryAlg

VertexFitAlg	fVertexFitAlg

unsigned short	algIndex

std::vector< short >	fMatchAlgs

std::vector< float >	fXMatchErr

std::vector< float >	fAngleMatchErr

std::vector< float >	fChgAsymFactor

std::vector< float >	fMatchMinLen

std::vector< bool >	fMakeAlgTracks

float	fMaxDAng

float	fChainMaxdX

float	fChainVtxAng

float	fMergeChgAsym

float	fMaxMergeError

float	fMergeErrorCut

float	fChgWindow {40}

float	fWirePitch

float	fFiducialCut

float	fDeltaRayCut

bool	fMakePFPs

unsigned short	fNVtxTrkHitsFit

float	fHitFitErrFac

bool	fuBCode

short	fDebugAlg

short	fDebugPlane

short	fDebugCluster

bool	fPrintAllClusters

bool	prt

std::array< unsigned int, 3 >	firstWire

std::array< unsigned int, 3 >	lastWire

std::array< unsigned int, 3 >	firstHit

std::array< unsigned int, 3 >	lastHit

std::array< std::vector< std::pair< int, int > >, 3 >	WireHitRange

std::vector< art::Ptr< recob::Hit > >	allhits

std::array< std::vector< clPar >, 3 >	cls

std::array< std::vector< ClsChainPar >, 3 >	clsChain

std::vector< vtxPar >	vtx

std::array< std::vector< unsigned short >, 3 >	vxCls

std::vector< TrkPar >	trk

std::array< std::vector< art::Ptr< recob::Hit > >, 3 >	trkHits

std::array< std::vector< art::Ptr< recob::Hit > >, 3 >	seedHits

std::array< float, 3 >	ChgNorm

std::vector< unsigned short >	pfpToTrkID

std::vector< MatchPars >	matcomb

Detailed Description

Definition at line 60 of file CCTrackMaker_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 26 of file Producer.h.

Constructor & Destructor Documentation

trkf::CCTrackMaker::CCTrackMaker ( fhicl::ParameterSet const & pset )

explicit

Definition at line 316 of file CCTrackMaker_module.cc.

References fAngleMatchErr, fChainMaxdX, fChainVtxAng, fChgAsymFactor, fClusterModuleLabel, fDebugAlg, fDebugCluster, fDebugPlane, fDeltaRayCut, fFiducialCut, fHitFitErrFac, fHitModuleLabel, fMakeAlgTracks, fMakePFPs, fMatchAlgs, fMatchMinLen, fMaxDAng, fMergeChgAsym, fNVtxTrkHitsFit, fPrintAllClusters, fuBCode, fVertexModuleLabel, and fXMatchErr.

                                                           : EDProducer{pset}
   {
     fHitModuleLabel = pset.get<std::string>("HitModuleLabel");
     fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
     fVertexModuleLabel = pset.get<std::string>("VertexModuleLabel");
     // track matching
     fMatchAlgs = pset.get<std::vector<short>>("MatchAlgs");
     fXMatchErr = pset.get<std::vector<float>>("XMatchErr");
     fAngleMatchErr = pset.get<std::vector<float>>("AngleMatchErr");
     fChgAsymFactor = pset.get<std::vector<float>>("ChgAsymFactor");
     fMatchMinLen = pset.get<std::vector<float>>("MatchMinLen");
     fMakeAlgTracks = pset.get<std::vector<bool>>("MakeAlgTracks");
     // Cluster merging
     fMaxDAng = pset.get<float>("MaxDAng");
     fChainMaxdX = pset.get<float>("ChainMaxdX");
     fChainVtxAng = pset.get<float>("ChainVtxAng");
     fMergeChgAsym = pset.get<float>("MergeChgAsym");
     // Cosmic ray tagging
     fFiducialCut = pset.get<float>("FiducialCut");
     fDeltaRayCut = pset.get<float>("DeltaRayCut");
     // make PFParticles
     fMakePFPs = pset.get<bool>("MakePFPs");
     // vertex fitting
     fNVtxTrkHitsFit = pset.get<unsigned short>("NVtxTrkHitsFit");
     fHitFitErrFac = pset.get<float>("HitFitErrFac");
     // uB code
     fuBCode = pset.get<bool>("uBCode");
     // debugging inputs
     fDebugAlg = pset.get<short>("DebugAlg");
     fDebugPlane = pset.get<short>("DebugPlane");
     fDebugCluster = pset.get<short>("DebugCluster");
     fPrintAllClusters = pset.get<bool>("PrintAllClusters");
 
     // Consistency check
     if (fMatchAlgs.size() > fXMatchErr.size() || fMatchAlgs.size() > fAngleMatchErr.size() ||
         fMatchAlgs.size() > fChgAsymFactor.size() || fMatchAlgs.size() > fMatchMinLen.size() ||
         fMatchAlgs.size() > fMakeAlgTracks.size()) {
       mf::LogError("CCTM") << "Incompatible fcl input vector sizes";
       return;
     }
     // Reality check
     for (unsigned short ii = 0; ii < fMatchAlgs.size(); ++ii) {
       if (fAngleMatchErr[ii] <= 0 || fXMatchErr[ii] <= 0) {
         mf::LogError("CCTM") << "Invalid matching parameters " << fAngleMatchErr[ii] << " "
                              << fXMatchErr[ii];
         return;
       }
     } // ii
 
     produces<std::vector<recob::PFParticle>>();
     produces<art::Assns<recob::PFParticle, recob::Track>>();
     produces<art::Assns<recob::PFParticle, recob::Cluster>>();
     produces<art::Assns<recob::PFParticle, recob::Seed>>();
     produces<art::Assns<recob::PFParticle, recob::Vertex>>();
     produces<std::vector<recob::Vertex>>();
     produces<std::vector<recob::Track>>();
     produces<art::Assns<recob::Track, recob::Hit>>();
     produces<std::vector<recob::Seed>>();
   }

Member Function Documentation

float trkf::CCTrackMaker::AngleFactor ( float slope )

private

Definition at line 3292 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), MakeClusterChains(), and PlnMatch().

   {
     float slp = fabs(slope);
     if (slp > 10.) slp = 30.;
     // return a value between 1 and 46
     return 1 + 0.05 * slp * slp;
   } // AngleFactor

void trkf::CCTrackMaker::AngMatch ( art::FindManyP< recob::Hit > const & fmCluHits )

private

float trkf::CCTrackMaker::ChargeAsym ( std::array< float, 3 > & mChg )

private

Definition at line 3038 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), and PlnMatch().

   {
     // find charge asymmetry between the cluster with the highest (lowest)
     // charge
     float big = 0, small = 1.E9;
     for (unsigned short ii = 0; ii < 3; ++ii) {
       if (mChg[ii] < small) small = mChg[ii];
       if (mChg[ii] > big) big = mChg[ii];
     }
     // chgAsym varies between 0 (perfect charge match) and 1 (dreadfull)
     return (big - small) / (big + small);
   } // CalculateChargeAsym

float trkf::CCTrackMaker::ChargeNear	(	geo::PlaneID const &	planeid,
		unsigned short	wire1,
		float	time1,
		unsigned short	wire2,
		float	time2
	)

private

Definition at line 3301 of file CCTrackMaker_module.cc.

References allhits, ChgNorm, geo::PlaneID::Plane, t1, and t2.

Referenced by FillEndMatch().

   {
     // returns the hit charge along a line between (wire1, time1) and
     // (wire2, time2)
 
     // put in increasing wire order (wire2 > wire1)
     unsigned short w1 = wire1;
     unsigned short w2 = wire2;
     double t1 = time1;
     double t2 = time2;
     double slp, prtime;
     if (w1 == w2) { slp = 0; }
     else {
       if (w1 > w2) {
         w1 = wire2;
         w2 = wire1;
         t1 = time2;
         t2 = time1;
       }
       slp = (t2 - t1) / (w2 - w1);
     }
 
     unsigned short wire;
 
     float chg = 0;
     for (unsigned short hit = 0; hit < allhits.size(); ++hit) {
       if (allhits[hit]->WireID().asPlaneID() != planeid) continue;
       wire = allhits[hit]->WireID().Wire;
       if (wire < w1) continue;
       if (wire > w2) continue;
       prtime = t1 + (wire - w1) * slp;
       if (prtime > allhits[hit]->PeakTimePlusRMS(3)) continue;
       if (prtime < allhits[hit]->PeakTimeMinusRMS(3)) continue;
       chg += ChgNorm[planeid.Plane] * allhits[hit]->Integral();
     } // hit
     return chg;
   } // ChargeNear

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::consumes ( InputTag const & tag )

protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

   {
     return collector_.consumes<T, BT>(tag);
   }

ConsumesCollector & art::ModuleBase::consumesCollector ( )

protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

   {
     return collector_;
   }

template<typename T , BranchType BT>

void art::ModuleBase::consumesMany ( )

protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

   {
     collector_.consumesMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::consumesView ( InputTag const & )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::consumesView ( InputTag const & tag )

inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

   {
     return collector_.consumesView<T, BT>(tag);
   }

void art::detail::Producer::doBeginJob ( SharedResources const & resources )

inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

   {
     setupQueues(resources);
     ProcessingFrame const frame{ScheduleID{}};
     beginJobWithFrame(frame);
   }

bool art::detail::Producer::doBeginRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

   {
     auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

bool art::detail::Producer::doBeginSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

   {
     auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

void art::detail::Producer::doEndJob ( )

inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     endJobWithFrame(frame);
   }

bool art::detail::Producer::doEndRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

   {
     auto r = rp.makeRun(mc, rp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

bool art::detail::Producer::doEndSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

   {
     auto sr = srp.makeSubRun(mc, srp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

bool art::detail::Producer::doEvent	(	EventPrincipal &	ep,
		ModuleContext const &	mc,
		std::atomic< std::size_t > &	counts_run,
		std::atomic< std::size_t > &	counts_passed,
		std::atomic< std::size_t > &	counts_failed
	)

inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

   {
     auto e = ep.makeEvent(mc);
     ++counts_run;
     ProcessingFrame const frame{mc.scheduleID()};
     produceWithFrame(e, frame);
     e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
     ++counts_passed;
     return true;
   }

void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const & fb )

inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const & fb )

inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseOutputFilesWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const & fb )

inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const & fb )

inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenOutputFilesWithFrame(fb, frame);
   }

bool trkf::CCTrackMaker::DupMatch	(	MatchPars &	match,
		unsigned short	nplanes
	)

private

Definition at line 2264 of file CCTrackMaker_module.cc.

References trkf::CCTrackMaker::MatchPars::Cls, trkf::CCTrackMaker::MatchPars::dAng, trkf::CCTrackMaker::MatchPars::dWir, trkf::CCTrackMaker::MatchPars::dX, trkf::CCTrackMaker::MatchPars::End, trkf::CCTrackMaker::MatchPars::Err, matcomb, trkf::CCTrackMaker::MatchPars::odAng, trkf::CCTrackMaker::MatchPars::odWir, trkf::CCTrackMaker::MatchPars::odX, trkf::CCTrackMaker::MatchPars::oErr, trkf::CCTrackMaker::MatchPars::oVtx, and trkf::CCTrackMaker::MatchPars::Vtx.

Referenced by PlnMatch(), and VtxMatch().

   {
     unsigned short nMatCl, nMiss;
     float toterr = match.Err + match.oErr;
     for (unsigned int imat = 0; imat < matcomb.size(); ++imat) {
       // check for exact matches
       if (match.Cls[0] == matcomb[imat].Cls[0] && match.Cls[1] == matcomb[imat].Cls[1] &&
           match.Cls[2] == matcomb[imat].Cls[2]) {
 
         // compare the error
         if (toterr < matcomb[imat].Err + matcomb[imat].oErr) {
           // keep the better one
           matcomb[imat].End[0] = match.End[0];
           matcomb[imat].End[1] = match.End[1];
           matcomb[imat].End[2] = match.End[2];
           matcomb[imat].Vtx = match.Vtx;
           matcomb[imat].dWir = match.dWir;
           matcomb[imat].dAng = match.dAng;
           matcomb[imat].dX = match.dX;
           matcomb[imat].Err = match.Err;
           matcomb[imat].oVtx = match.oVtx;
           matcomb[imat].odWir = match.odWir;
           matcomb[imat].odAng = match.odAng;
           matcomb[imat].odX = match.odX;
           matcomb[imat].oErr = match.oErr;
         }
         return true;
       } // test
       // check for a 3-plane match vs 2-plane match
       nMatCl = 0;
       nMiss = 0;
       for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
         if (match.Cls[ipl] >= 0) {
           if (match.Cls[ipl] == matcomb[imat].Cls[ipl] &&
               (match.End[0] == matcomb[imat].End[0] || match.End[1] == matcomb[imat].End[1]))
             ++nMatCl;
         }
         else {
           ++nMiss;
         }
       } // ipl
       if (nMatCl == 2 && nMiss == 1) return true;
     } // imat
     return false;
   } // DupMatch

float trkf::CCTrackMaker::dXClTraj	(	art::FindManyP< recob::Hit > const &	fmCluHits,
		unsigned short	ipl,
		unsigned short	icl1,
		unsigned short	end1
	)

private

Definition at line 1845 of file CCTrackMaker_module.cc.

References cls, and fWirePitch.

Referenced by MakeClusterChains().

   {
     // project cluster icl1 at end1 to find the best intersection with icl2
     float dw, dx, best = 999;
     std::vector<art::Ptr<recob::Hit>> clusterhits = fmCluHits.at(cls[ipl][icl1].EvtIndex);
     for (unsigned short hit = 0; hit < clusterhits.size(); ++hit) {
       dw = clusterhits[hit]->WireID().Wire - cls[ipl][icl1].Wire[end1];
       dx = fabs(cls[ipl][icl1].Time[end1] + dw * fWirePitch * cls[ipl][icl1].Slope[end1] -
                 clusterhits[hit]->PeakTime());
       if (dx < best) best = dx;
       if (dx < 0.01) break;
     } // hit
     return best;
   } // dXClTraj

void trkf::CCTrackMaker::FillChgNear	(	detinfo::DetectorPropertiesData const &	detProp,
		geo::TPCID const &	tpcid
	)

private

Definition at line 888 of file CCTrackMaker_module.cc.

References allhits, ChgNorm, cls, detinfo::DetectorPropertiesData::ConvertXToTicks(), dir, util::end(), fChgWindow, firstWire, fWirePitch, geo::Iterable< IterationPolicy, Transform >::Iterate(), w, WireHitRange, and wireReadoutGeom.

Referenced by produce().

   {
     // fills the CHgNear array
 
     unsigned short wire, nwires, indx;
     float dir, ctime, cx, chgWinLo, chgWinHi;
     float cnear;
 
     for (auto const& planeid : wireReadoutGeom.Iterate<geo::PlaneID>(tpcid)) {
       auto const ipl = planeid.Plane;
       for (unsigned short icl = 0; icl < cls[ipl].size(); ++icl) {
         // find the nearby charge at the US and DS ends
         nwires = cls[ipl][icl].Length / 2;
         if (nwires < 2) continue;
         // maximum of 30 wires for long clusters
         if (nwires > 30) nwires = 30;
         for (unsigned short end = 0; end < 2; ++end) {
           cnear = 0;
           // direction for adding/subtracting wire numbers
           dir = 1 - 2 * end;
           for (unsigned short w = 0; w < nwires; ++w) {
             wire = cls[ipl][icl].Wire[end] + dir * w;
             cx = cls[ipl][icl].X[end] + dir * w * cls[ipl][icl].Slope[end] * fWirePitch;
             ctime = detProp.ConvertXToTicks(cx, planeid);
             chgWinLo = ctime - fChgWindow;
             chgWinHi = ctime + fChgWindow;
             indx = wire - firstWire[ipl];
             if (WireHitRange[ipl][indx].first < 0) continue;
             unsigned int firhit = WireHitRange[ipl][indx].first;
             unsigned int lashit = WireHitRange[ipl][indx].second;
             for (unsigned int hit = firhit; hit < lashit; ++hit) {
               if (hit > allhits.size() - 1) {
                 mf::LogError("CCTM")
                   << "FillChgNear bad lashit " << lashit << " size " << allhits.size() << "\n";
                 continue;
               }
               if (allhits[hit]->PeakTime() < chgWinLo) continue;
               if (allhits[hit]->PeakTime() > chgWinHi) continue;
               cnear += allhits[hit]->Integral();
             } // hit
           }   // w
           cnear /= (float)(nwires - 1);
           if (cnear > cls[ipl][icl].Charge[end]) {
             cls[ipl][icl].ChgNear[end] = ChgNorm[ipl] * cnear / cls[ipl][icl].Charge[end];
           }
           else {
             cls[ipl][icl].ChgNear[end] = 0;
           }
         } // end
       }   // icl
     }     // ipl
 
   } // FillChgNear

void trkf::CCTrackMaker::FillEndMatch	(	detinfo::DetectorPropertiesData const &	detProp,
		geo::TPCID const &	tpcid,
		MatchPars &	match
	)

private

Definition at line 2579 of file CCTrackMaker_module.cc.

Referenced by PlnMatch(), and VtxMatch().

   {
     // fill the matching parameters for this cluster match. The calling routine
     // should set the match end vertex ID (if applicable) as well as the
     // cluster IDs and matching ends in each plane. Note that the matching variables
     // Note that dWir, dAng and dTim are not filled if there is a vertex (match.Vtx >= 0).
     // Likewise, odWir, odAng and odX are not filled if there is a vertex match
     // at the other end
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
 
     if (nplanes == 2) {
       FillEndMatch2(match);
       return;
     }
 
     std::array<short, 3> mVtx;
     std::array<short, 3> oVtx;
     std::array<float, 3> oWir;
     std::array<float, 3> oSlp;
     std::array<float, 3> oAng;
     std::array<float, 3> oX;
 
     std::array<float, 3> mChg;
 
     unsigned short ii, ipl, iend, jpl, jend, kpl, kend, oend;
     short icl, jcl, kcl;
 
     for (ipl = 0; ipl < 3; ++ipl) {
       mVtx[ipl] = -1;
       oVtx[ipl] = -1;
       oWir[ipl] = -66;
       oSlp[ipl] = -66;
       oAng[ipl] = -66;
       oX[ipl] = -66;
       mChg[ipl] = -1;
     } // ipl
 
     // initialize parameters that shouldn't have been set by the calling routine
     match.dWir = 0;
     match.dAng = 0;
     match.dX = 0;
     match.Err = 100;
     match.odWir = 0;
     match.odAng = 0;
     match.odX = 0;
     match.oErr = 100;
     match.oVtx = -1;
 
     if (prt) {
       mf::LogVerbatim myprt("CCTM");
       myprt << "FEM ";
       for (ipl = 0; ipl < nplanes; ++ipl) {
         myprt << " " << ipl << ":" << match.Cls[ipl] << ":" << match.End[ipl];
       }
     }
 
     short missingPlane = -1;
     unsigned short nClInPln = 0;
     // number of vertex matches at each end
     short aVtx = -1;
     unsigned short novxmat = 0;
     short aoVtx = -1;
     unsigned short nvxmat = 0;
     unsigned short nShortCl = 0;
     // fill the other end parameters in each plane
     for (ipl = 0; ipl < nplanes; ++ipl) {
       if (match.Cls[ipl] < 0) {
         missingPlane = ipl;
         continue;
       }
       ++nClInPln;
       icl = match.Cls[ipl];
       match.Chg[ipl] = clsChain[ipl][icl].TotChg;
       mChg[ipl] = clsChain[ipl][icl].TotChg;
       iend = match.End[ipl];
       mVtx[ipl] = clsChain[ipl][icl].VtxIndex[iend];
       if (clsChain[ipl][icl].Length < 6) ++nShortCl;
       if (mVtx[ipl] >= 0) {
         if (aVtx < 0) aVtx = mVtx[ipl];
         if (mVtx[ipl] == aVtx) ++nvxmat;
       }
       if (prt)
         mf::LogVerbatim("CCTM") << "FEM: m " << ipl << ":" << icl << ":" << iend << " Vtx "
                                 << mVtx[ipl] << "  Wir " << clsChain[ipl][icl].Wire[iend]
                                 << std::fixed << std::setprecision(3) << " Slp "
                                 << clsChain[ipl][icl].Slope[iend] << std::fixed
                                 << std::setprecision(1) << " X " << clsChain[ipl][icl].X[iend];
 
       oend = 1 - iend;
       oWir[ipl] = clsChain[ipl][icl].Wire[oend];
       oAng[ipl] = clsChain[ipl][icl].Angle[oend];
       oSlp[ipl] = clsChain[ipl][icl].Slope[oend];
       oX[ipl] = clsChain[ipl][icl].X[oend];
       oVtx[ipl] = clsChain[ipl][icl].VtxIndex[oend];
       if (oVtx[ipl] >= 0) {
         if (aoVtx < 0) aoVtx = oVtx[ipl];
         if (oVtx[ipl] == aoVtx) ++novxmat;
       }
 
       if (prt)
         mf::LogVerbatim("CCTM") << "     o " << ipl << ":" << icl << ":" << oend << " oVtx "
                                 << oVtx[ipl] << " oWir " << oWir[ipl] << std::fixed
                                 << std::setprecision(3) << " oSlp " << oSlp[ipl] << std::fixed
                                 << std::setprecision(1) << " oX " << oX[ipl] << " Chg "
                                 << (int)mChg[ipl];
 
     } // ipl
 
     bool isShort = (nShortCl > 1);
 
     if (nClInPln < 2) {
       mf::LogWarning("CCTM") << "Not enough matched planes supplied";
       return;
     }
 
     if (prt)
       mf::LogVerbatim("CCTM") << "FEM: Vtx m " << aVtx << " count " << nvxmat << " o " << aoVtx
                               << " count " << novxmat << " missingPlane " << missingPlane
                               << " nClInPln " << nClInPln;
 
     // perfect match
     if (nvxmat == 3 && novxmat == 3) {
       match.Vtx = aVtx;
       match.Err = 0;
       match.oVtx = aoVtx;
       match.oErr = 0;
       return;
     }
 
     // 2-plane vertex match?
     // factors applied to error = 1 (no vtx), 0.5 (2 pln vtx), 0.33 (3 pln vtx)
     float vxFactor = 1;
     float ovxFactor = 1;
     if (nClInPln == 3) {
       // a cluster in all 3 planes
       if (nvxmat == 3) {
         // and all vertex assignments agree at the match end
         match.Vtx = aVtx;
         vxFactor = 0.33;
       }
       if (novxmat == 3) {
         // and all vertex assignments agree at the other end
         match.oVtx = aoVtx;
         ovxFactor = 0.33;
       }
     }
     else {
       // a cluster in 2 planes
       if (nvxmat == 2) {
         match.Vtx = aVtx;
         vxFactor = 0.5;
       }
       if (novxmat == 2) {
         match.oVtx = aoVtx;
         ovxFactor = 0.5;
       }
     } // nClInPln
 
     // find wire, X and Time at both ends
 
     // Find the "other end" matching parameters with
     // two cases: a 3-plane match or a 2-plane match
     // and with/without an other end vertex
 
     float kWir, okWir;
     float kSlp, okSlp, kAng, okAng, okX, kX, kTim, okTim;
 
     if (nClInPln == 3) {
       ipl = 0;
       jpl = 1;
       kpl = 2;
     }
     else {
       // 2-plane match
       kpl = missingPlane;
       if (kpl == 0) {
         ipl = 1;
         jpl = 2;
       }
       else if (kpl == 1) {
         ipl = 2;
         jpl = 0;
       }
       else {
         ipl = 0;
         jpl = 1;
       } // kpl test
     }   // missing plane
     iend = match.End[ipl];
     jend = match.End[jpl];
     icl = match.Cls[ipl];
     jcl = match.Cls[jpl];
     if (nplanes > 2) {
       kcl = match.Cls[kpl];
       kend = match.End[kpl];
     }
 
     geo::PlaneID const iplaneid{tpcid, ipl};
     geo::PlaneID const jplaneid{tpcid, jpl};
     geo::PlaneID const kplaneid{tpcid, kpl};
     auto const& kplane = wireReadoutGeom.Plane(kplaneid);
     okSlp = wireReadoutGeom.ThirdPlaneSlope(ipl, oSlp[ipl], jpl, oSlp[jpl], tpcid);
     okAng = atan(okSlp);
     // handle the case near pi/2, where the errors on large slopes could result in
     // a wrong-sign kAng
     bool ignoreSign = (fabs(okSlp) > 10);
     if (ignoreSign) okAng = fabs(okAng);
     if (match.oVtx >= 0) {
       // a vertex exists at the other end
       okWir = kplane.WireCoordinate(geo::Point_t{0, vtx[match.oVtx].Y, vtx[match.oVtx].Z});
       okX = vtx[match.oVtx].X;
     }
     else {
       // no vertex at the other end
       auto const [ypos, zpos, _] =
         wireReadoutGeom
           .WireIDsIntersect(geo::WireID(iplaneid, oWir[ipl]), geo::WireID(jplaneid, oWir[jpl]))
           .value_or(geo::WireIDIntersection::invalid());
       okWir = 0.5 + kplane.WireCoordinate(geo::Point_t{0, ypos, zpos});
       okX = 0.5 * (oX[ipl] + oX[jpl]);
     }
     okTim = detProp.ConvertXToTicks(okX, kplaneid);
     if (prt)
       mf::LogVerbatim("CCTM") << "FEM: oEnd"
                               << " kpl " << kpl << " okSlp " << okSlp << " okAng " << okAng
                               << " okWir " << (int)okWir << " okX " << okX;
 
 
     kSlp = wireReadoutGeom.ThirdPlaneSlope(
       ipl, clsChain[ipl][icl].Slope[iend], jpl, clsChain[jpl][jcl].Slope[jend], kplaneid);
     kAng = atan(kSlp);
     if (ignoreSign) kAng = fabs(kAng);
     if (match.Vtx >= 0) {
       if (vtx.size() == 0 || (unsigned int)match.Vtx > vtx.size() - 1) {
         mf::LogError("CCTM") << "FEM: Bad match.Vtx " << match.Vtx << " vtx size " << vtx.size();
         return;
       }
       // a vertex exists at the match end
       kWir = kplane.WireCoordinate(geo::Point_t{0, vtx[match.Vtx].Y, vtx[match.Vtx].Z});
       kX = vtx[match.Vtx].X;
     }
     else {
       // no vertex at the match end
       auto const [ypos, zpos, _] =
         wireReadoutGeom
           .WireIDsIntersect(geo::WireID(iplaneid, clsChain[ipl][icl].Wire[iend]),
                             geo::WireID(jplaneid, clsChain[jpl][jcl].Wire[jend]))
           .value_or(geo::WireIDIntersection::invalid());
       kWir = 0.5 + kplane.WireCoordinate(geo::Point_t{0, ypos, zpos});
       kX = 0.5 * (clsChain[ipl][icl].X[iend] + clsChain[jpl][jcl].X[jend]);
     }
     kTim = detProp.ConvertXToTicks(kX, kplaneid);
     if (prt)
       mf::LogVerbatim("CCTM") << "FEM: mEnd"
                               << " kpl " << kpl << " kSlp " << kSlp << " kAng " << kAng << " kX "
                               << kX;
 
     // try to find a 3-plane match using this information
     if (nClInPln < 3 && FindMissingCluster(kpl, kcl, kend, kWir, kX, okWir, okX)) {
       nClInPln = 3;
       // update local variables
       match.Cls[kpl] = kcl;
       match.End[kpl] = kend;
       match.Chg[kpl] = clsChain[kpl][kcl].TotChg;
       mChg[kpl] = clsChain[kpl][kcl].TotChg;
       oend = 1 - kend;
       oWir[kpl] = clsChain[kpl][kcl].Wire[oend];
       oX[kpl] = clsChain[kpl][kcl].X[oend];
       oAng[kpl] = clsChain[kpl][kcl].Angle[oend];
       oSlp[kpl] = clsChain[kpl][kcl].Slope[oend];
     } // FindMissingCluster
 
     // decide whether to continue with a 2-plane match. The distance between match and other end should
     // be large enough to create a cluster
     if (nClInPln == 2 && fabs(okWir - kWir) > 3) return;
 
     // Calculate the cluster charge asymmetry. This factor will be applied
     // to the error of the end matches
     float chgAsym = 1;
     // Get the charge in the plane without a matching cluster
     if (nClInPln < 3 && mChg[missingPlane] <= 0) {
       if (missingPlane != kpl)
         mf::LogError("CCTM") << "FEM bad missingPlane " << missingPlane << " " << kpl << "\n";
       mChg[kpl] = ChargeNear(kplaneid, (unsigned short)kWir, kTim, (unsigned short)okWir, okTim);
       match.Chg[kpl] = mChg[kpl];
       if (prt)
         mf::LogVerbatim("CCTM") << "FEM: Missing cluster in " << kpl << " ChargeNear " << (int)kWir
                                 << ":" << (int)kTim << " " << (int)okWir << ":" << (int)okTim
                                 << " chg " << mChg[kpl];
       if (mChg[kpl] <= 0) return;
     }
 
     if (fChgAsymFactor[algIndex] > 0) {
       chgAsym = ChargeAsym(mChg);
       if (chgAsym > 0.5) return;
       chgAsym = 1 + fChgAsymFactor[algIndex] * chgAsym;
     }
 
     if (prt) mf::LogVerbatim("CCTM") << "FEM: charge asymmetry factor " << chgAsym;
     float sigmaX, sigmaA;
     float da, dx, dw;
 
     // check for vertex consistency at the match end
     aVtx = -1;
     bool allPlnVtxMatch = false;
     if (nClInPln == 3) {
       unsigned short nmvtx = 0;
       for (ii = 0; ii < nplanes; ++ii) {
         if (mVtx[ii] >= 0) {
           if (aVtx < 0) aVtx = mVtx[ii];
           ++nmvtx;
         }
       } // ii
       // same vertex in all planes
       if (nmvtx) allPlnVtxMatch = true;
     } // nClInPln
 
     // inflate the X match error to allow for missing one wire on the end of a cluster
     sigmaX = fXMatchErr[algIndex] + std::max(kSlp, (float)20);
     sigmaA = fAngleMatchErr[algIndex] * AngleFactor(kSlp);
     if (prt)
       mf::LogVerbatim("CCTM") << "bb " << algIndex << "  " << fXMatchErr[algIndex] << " "
                               << fAngleMatchErr[algIndex] << " kslp " << kSlp;
 
     if (nClInPln == 3) {
       kcl = match.Cls[kpl];
       kend = match.End[kpl];
       dw = kWir - clsChain[kpl][kcl].Wire[kend];
       match.dWir = dw;
       if (fabs(match.dWir) > 100) return;
       if (match.Vtx >= 0) { match.dX = kX - clsChain[kpl][kcl].X[kend]; }
       else {
         match.dX = std::abs(clsChain[ipl][icl].X[iend] - clsChain[jpl][jcl].X[jend]) +
                    std::abs(clsChain[ipl][icl].X[iend] - clsChain[kpl][kcl].X[kend]);
       }
       if (prt) mf::LogVerbatim("CCTM") << " dw " << dw << " dx " << match.dX;
       // TODO: Angle matching has problems with short clusters
       if (!isShort) {
         if (ignoreSign) { match.dAng = kAng - fabs(clsChain[kpl][kcl].Angle[kend]); }
         else {
           match.dAng = kAng - clsChain[kpl][kcl].Angle[kend];
         }
       } // !isShort
       da = fabs(match.dAng) / sigmaA;
       dx = fabs(match.dX) / sigmaX;
       if (allPlnVtxMatch) {
         // matched vertex. Use angle for match error
         match.Err = vxFactor * chgAsym * da / 3;
         if (prt) mf::LogVerbatim("CCTM") << " 3-pln w Vtx  match.Err " << match.Err;
       }
       else {
         dw /= 2;
         // divide by 9
         match.Err = vxFactor * chgAsym * sqrt(dx * dx + da * da + dw * dw) / 9;
         if (prt) mf::LogVerbatim("CCTM") << " 3-pln match.Err " << match.Err;
       }
     }
     else {
       // 2-plane match
       match.dWir = -1;
       match.dAng = -1;
       match.dX = clsChain[ipl][icl].X[iend] - clsChain[jpl][jcl].X[jend];
       // degrade error by 3 for 2-plane matches
       match.Err = 3 + vxFactor * chgAsym * fabs(match.dX) / sigmaX;
       if (prt) mf::LogVerbatim("CCTM") << " 2-pln Err " << match.Err;
     } // !(nClInPln == 3)
 
     if (nClInPln == 3) {
       // A cluster in all 3 planes
       dw = okWir - oWir[kpl];
       match.odWir = dw;
       if (match.oVtx >= 0) { match.odX = okX - oX[kpl]; }
       else {
         match.odX = std::abs(oX[ipl] - oX[jpl]) + std::abs(oX[ipl] - oX[kpl]);
       }
       if (prt)
         mf::LogVerbatim("CCTM") << " odw " << match.odWir << " odx " << match.odX << " sigmaX "
                                 << sigmaX;
       // TODO: CHECK FOR SHOWER-LIKE CLUSTER OTHER END MATCH
       if (!isShort) {
         if (ignoreSign) { match.odAng = okAng - fabs(oAng[kpl]); }
         else {
           match.odAng = okAng - oAng[kpl];
         }
       } // !isShort
       da = match.odAng / sigmaA;
       dx = fabs(match.odX) / sigmaX;
       // error for wire number match
       dw /= 2;
       // divide by number of planes with clusters * 3 for dx, da and dw
       match.oErr = ovxFactor * chgAsym * sqrt(dx * dx + da * da + dw * dw) / 9;
       if (prt) mf::LogVerbatim("CCTM") << " 3-pln match.oErr " << match.oErr;
     }
     else {
       // Only 2 clusters in 3 planes
       match.odX = (oX[ipl] - oX[jpl]) / sigmaX;
       match.oErr = 3 + ovxFactor * chgAsym * fabs(match.odX);
       if (prt) mf::LogVerbatim("CCTM") << " 2-pln match.oErr " << match.oErr;
     }
   } // FillEndMatch

void trkf::CCTrackMaker::FillEndMatch2 ( MatchPars & match )

private

Definition at line 2491 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), and VtxMatch().

   {
     // 2D version of FillEndMatch
 
     match.Err = 100;
     match.oErr = 100;
     match.Chg[2] = 0;
     match.dWir = 0;
     match.dAng = 0;
     match.odWir = 0;
     match.odAng = 0;
 
     unsigned short ipl = 0;
     unsigned short jpl = 1;
 
     if (match.Cls[0] < 0 || match.Cls[1] < 0) return;
 
     unsigned short icl = match.Cls[ipl];
     unsigned short iend = match.End[ipl];
     match.Chg[ipl] = clsChain[ipl][icl].TotChg;
     // cluster i match end
     float miX = clsChain[ipl][icl].X[iend];
     // cluster i other end
     unsigned short oiend = 1 - iend;
     float oiX = clsChain[ipl][icl].X[oiend];
 
     unsigned short jcl = match.Cls[jpl];
     unsigned short jend = match.End[jpl];
     match.Chg[jpl] = clsChain[jpl][jcl].TotChg;
     // cluster j match end
     float mjX = clsChain[jpl][jcl].X[jend];
     // cluster j other end
     unsigned short ojend = 1 - jend;
     float ojX = clsChain[jpl][jcl].X[ojend];
 
     // look for a match end vertex match
     match.Vtx = -1;
     if (clsChain[ipl][icl].VtxIndex[iend] >= 0 &&
         clsChain[ipl][icl].VtxIndex[iend] == clsChain[jpl][jcl].VtxIndex[jend]) {
       match.Vtx = clsChain[ipl][icl].VtxIndex[iend];
       miX = vtx[match.Vtx].X;
       mjX = vtx[match.Vtx].X;
     }
 
     // look for an other end vertex match
     match.oVtx = -1;
     if (clsChain[ipl][icl].VtxIndex[oiend] >= 0 &&
         clsChain[ipl][icl].VtxIndex[oiend] == clsChain[jpl][jcl].VtxIndex[ojend]) {
       match.oVtx = clsChain[ipl][icl].VtxIndex[oiend];
       oiX = vtx[match.oVtx].X;
       ojX = vtx[match.oVtx].X;
     }
 
     // find the charge asymmetry
     float chgAsym = 1;
     if (fChgAsymFactor[algIndex] > 0) {
       chgAsym = fabs(match.Chg[ipl] - match.Chg[jpl]) / (match.Chg[ipl] + match.Chg[jpl]);
       if (chgAsym > 0.5) return;
       chgAsym = 1 + fChgAsymFactor[algIndex] * chgAsym;
     }
 
     // find the error at the match end
     float maxSlp = fabs(clsChain[ipl][icl].Slope[iend]);
     if (fabs(clsChain[jpl][jcl].Slope[jend]) > maxSlp)
       maxSlp = fabs(clsChain[jpl][jcl].Slope[jend]);
     float sigmaX = fXMatchErr[algIndex] + std::max(maxSlp, (float)20);
     match.dX = fabs(miX - mjX);
     match.Err = chgAsym * match.dX / sigmaX;
 
     // find the error at the other end
     maxSlp = fabs(clsChain[ipl][icl].Slope[oiend]);
     if (fabs(clsChain[jpl][jcl].Slope[ojend]) > maxSlp)
       maxSlp = fabs(clsChain[jpl][jcl].Slope[ojend]);
     sigmaX = fXMatchErr[algIndex] + std::max(maxSlp, (float)20);
     match.odX = fabs(oiX - ojX);
     match.oErr = chgAsym * match.odX / sigmaX;
 
     if (prt)
       mf::LogVerbatim("CCTM") << "FEM2: m " << ipl << ":" << icl << ":" << iend << " miX " << miX
                               << " - " << jpl << ":" << jcl << ":" << jend << " mjX " << mjX
                               << " match.dX " << match.dX << " match.Err " << match.Err
                               << " chgAsym " << chgAsym << " o "
                               << " oiX " << oiX << " ojX " << ojX << " match.odX " << match.odX
                               << " match.oErr " << match.oErr << "\n";
 
   } // FillEndMatch2

void art::Modifier::fillProductDescriptions ( )

inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

   {
     ProductRegistryHelper::fillDescriptions(moduleDescription());
   }

void trkf::CCTrackMaker::FillTrkHits	(	art::FindManyP< recob::Hit > const &	fmCluHits,
		unsigned short	imat,
		unsigned short	nplanes
	)

private

Definition at line 3052 of file CCTrackMaker_module.cc.

References clear(), cls, clsChain, matcomb, prt, util::size(), and trkHits.

Referenced by SortMatches().

   {
     // Fills the trkHits vector using cluster hits associated with the match combo imat
 
     unsigned short ipl;
 
     for (ipl = 0; ipl < 3; ++ipl)
       trkHits[ipl].clear();
 
     if (imat > matcomb.size()) return;
 
     unsigned short indx;
     std::vector<art::Ptr<recob::Hit>> clusterhits;
     unsigned short icc, ccl, icl, ecl, iht, ii;
     short endOrder, fillOrder;
 
     if (prt)
       mf::LogVerbatim("CCTM") << "In FillTrkHits: matcomb " << imat << " cluster chains "
                               << matcomb[imat].Cls[0] << " " << matcomb[imat].Cls[1] << " "
                               << matcomb[imat].Cls[2];
 
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
       if (matcomb[imat].Cls[ipl] < 0) continue;
       // ccl is the cluster chain index
       ccl = matcomb[imat].Cls[ipl];
       // endOrder = 1 for normal order (hits added from US to DS), and -1 for reverse order
       endOrder = 1 - 2 * matcomb[imat].End[ipl];
       // re-order the sequence of cluster indices for reverse order
       if (endOrder < 0) {
         std::reverse(clsChain[ipl][ccl].ClsIndex.begin(), clsChain[ipl][ccl].ClsIndex.end());
         std::reverse(clsChain[ipl][ccl].Order.begin(), clsChain[ipl][ccl].Order.end());
         for (ii = 0; ii < clsChain[ipl][ccl].Order.size(); ++ii)
           clsChain[ipl][ccl].Order[ii] = 1 - clsChain[ipl][ccl].Order[ii];
       }
       if (ccl > clsChain[ipl].size() - 1) {
         mf::LogError("CCTM") << "Bad cluster chain index " << ccl << " in plane " << ipl;
         continue;
       }
       // loop over all the clusters in the chain
       for (icc = 0; icc < clsChain[ipl][ccl].ClsIndex.size(); ++icc) {
         icl = clsChain[ipl][ccl].ClsIndex[icc];
         if (icl > fmCluHits.size() - 1) {
           std::cout << "oops in FTH " << icl << " clsChain size " << clsChain[ipl].size() << "\n";
           exit(1);
         }
         ecl = cls[ipl][icl].EvtIndex;
         if (ecl > fmCluHits.size() - 1) {
           std::cout << "FTH bad EvtIndex " << ecl << " fmCluHits size " << fmCluHits.size() << "\n";
           continue;
         }
         clusterhits = fmCluHits.at(ecl);
         if (clusterhits.size() == 0) {
           std::cout << "FTH no cluster hits for EvtIndex " << cls[ipl][icl].EvtIndex << "\n";
           continue;
         }
         indx = trkHits[ipl].size();
         trkHits[ipl].resize(indx + clusterhits.size());
         // ensure the hit fill ordering is consistent
         fillOrder = 1 - 2 * clsChain[ipl][ccl].Order[icc];
         if (fillOrder == 1) {
           for (iht = 0; iht < clusterhits.size(); ++iht) {
             if (indx + iht > trkHits[ipl].size() - 1) std::cout << "FTH oops3\n";
             trkHits[ipl][indx + iht] = clusterhits[iht];
           }
         }
         else {
           for (ii = 0; ii < clusterhits.size(); ++ii) {
             iht = clusterhits.size() - 1 - ii;
             if (indx + ii > trkHits[ipl].size() - 1) std::cout << "FTH oops4\n";
             trkHits[ipl][indx + ii] = clusterhits[iht];
           } // ii
         }
       } // icc
       ii = trkHits[ipl].size() - 1;
       if (prt)
         mf::LogVerbatim("CCTM") << "plane " << ipl << " first p " << trkHits[ipl][0]->WireID().Plane
                                 << " w " << trkHits[ipl][0]->WireID().Wire << ":"
                                 << (int)trkHits[ipl][0]->PeakTime() << " last p "
                                 << trkHits[ipl][ii]->WireID().Plane << " w "
                                 << trkHits[ipl][ii]->WireID().Wire << ":"
                                 << (int)trkHits[ipl][ii]->PeakTime();
     } // ipl
 
     // TODO Check the ends of trkHits to see if there are missing hits that should have been included
     // in a cluster
 
   } // FillTrkHits

void trkf::CCTrackMaker::FillWireHitRange ( geo::TPCID const & tpcid )

private

Definition at line 3344 of file CCTrackMaker_module.cc.

References allhits, ChgNorm, DEFINE_ART_MODULE, firstHit, firstWire, lastHit, lastWire, geo::WireReadoutGeom::Nplanes(), geo::WireReadoutGeom::Nwires(), w, WireHitRange, and wireReadoutGeom.

Referenced by produce().

   {
     // fills the WireHitRange vector. Slightly modified version of the one in ClusterCrawlerAlg
 
     for (unsigned int ipl = 0; ipl < 3; ++ipl) {
       firstWire[ipl] = INT_MAX;
       lastWire[ipl] = 0;
       firstHit[ipl] = INT_MAX;
       lastHit[ipl] = 0;
       WireHitRange[ipl].clear();
       ChgNorm[ipl] = 0;
     }
 
     // find the first and last wire with a hit
     unsigned short oldipl = 0;
     for (unsigned int hit = 0; hit < allhits.size(); ++hit) {
       if (static_cast<geo::TPCID const&>(allhits[hit]->WireID()) != tpcid) continue;
       auto const ipl = allhits[hit]->WireID().Plane;
       if (allhits[hit]->WireID().Wire >
           wireReadoutGeom.Nwires(allhits[hit]->WireID().asPlaneID())) {
         if (lastWire[ipl] < firstWire[ipl]) {
           mf::LogError("CCTM") << "Invalid WireID().Wire " << allhits[hit]->WireID().Wire;
           return;
         }
       }
       //      ChgNorm[ipl] += allhits[hit]->Integral();
       if (ipl < oldipl) {
         mf::LogError("CCTM") << "Hits are not in increasing-plane order\n";
         return;
       }
       oldipl = ipl;
       if (firstHit[ipl] == INT_MAX) {
         firstHit[ipl] = hit;
         firstWire[ipl] = allhits[hit]->WireID().Wire;
       }
       lastHit[ipl] = hit;
       lastWire[ipl] = allhits[hit]->WireID().Wire;
     } // hit
 
     // xxx
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
       if (lastWire[ipl] < firstWire[ipl]) {
         mf::LogError("CCTM") << "Invalid first/last wire in plane " << ipl;
         return;
       }
     } // ipl
 
     // normalize charge in induction planes to the collection plane
     //    for(ipl = 0; ipl < nplanes; ++ipl) ChgNorm[ipl] = ChgNorm[nplanes - 1] / ChgNorm[ipl];
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl)
       ChgNorm[ipl] = 1;
 
     // now we can define the WireHitRange vector.
     int sflag, nwires, wire;
     unsigned int indx, thisWire, thisHit, lastFirstHit;
     //uint32_t chan;
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
       nwires = lastWire[ipl] - firstWire[ipl] + 1;
       WireHitRange[ipl].resize(nwires);
       // start by defining the "no hits on wire" condition
       sflag = -2;
       for (wire = 0; wire < nwires; ++wire)
         WireHitRange[ipl][wire] = std::make_pair(sflag, sflag);
       // overwrite with the "dead wires" condition
       sflag = -1;
       // next overwrite with the index of the first/last hit on each wire
       lastWire[ipl] = firstWire[ipl];
       thisHit = firstHit[ipl];
       lastFirstHit = firstHit[ipl];
       for (unsigned int hit = firstHit[ipl]; hit <= lastHit[ipl]; ++hit) {
         thisWire = allhits[hit]->WireID().Wire;
         if (thisWire > lastWire[ipl]) {
           indx = lastWire[ipl] - firstWire[ipl];
           int tmp1 = lastFirstHit;
           int tmp2 = thisHit;
           WireHitRange[ipl][indx] = std::make_pair(tmp1, tmp2);
           lastWire[ipl] = thisWire;
           lastFirstHit = thisHit;
         }
         else if (thisWire < lastWire[ipl]) {
           mf::LogError("CCTM") << "Hit not in proper order in plane " << ipl;
           exit(1);
         }
         ++thisHit;
       } // hit
       // define the last wire
       indx = lastWire[ipl] - firstWire[ipl];
       int tmp1 = lastFirstHit;
       ++lastHit[ipl];
       int tmp2 = lastHit[ipl];
       WireHitRange[ipl][indx] = std::make_pair(tmp1, tmp2);
       // add one to lastWire and lastHit for more standard indexing
       ++lastWire[ipl];
     } // ipl
 
     // error checking
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
       if (firstWire[ipl] < INT_MAX) continue;
       if (lastWire[ipl] > 0) continue;
       if (firstHit[ipl] < INT_MAX) continue;
       if (lastHit[ipl] > 0) continue;
       std::cout << "FWHR problem\n";
       exit(1); // FIXME: Should not exit from library code
     }          // ipl
 
     unsigned int nht = 0;
     std::vector<bool> hchk(allhits.size());
     for (unsigned int ii = 0; ii < hchk.size(); ++ii)
       hchk[ii] = false;
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
       for (unsigned int w = firstWire[ipl]; w < lastWire[ipl]; ++w) {
         indx = w - firstWire[ipl];
         if (indx > lastWire[ipl]) {
           std::cout << "FWHR bad index " << indx << "\n";
           exit(1);
         }
         // no hit on this wire
         if (WireHitRange[ipl][indx].first < 0) continue;
         unsigned int firhit = WireHitRange[ipl][indx].first;
         unsigned int lashit = WireHitRange[ipl][indx].second;
         for (unsigned int hit = firhit; hit < lashit; ++hit) {
           ++nht;
           if (hit > allhits.size() - 1) {
             std::cout << "FWHR: Bad hchk " << hit << " size " << allhits.size() << "\n";
             continue;
           }
           hchk[hit] = true;
           if (allhits[hit]->WireID().Plane != ipl || allhits[hit]->WireID().Wire != w) {
             std::cout << "FWHR bad plane " << allhits[hit]->WireID().Plane << " " << ipl
                       << " or wire " << allhits[hit]->WireID().Wire << " " << w << "\n";
             exit(1);
           }
         } // hit
       }   // w
     }     // ipl
 
     if (nht != allhits.size()) {
       std::cout << "FWHR hit count problem " << nht << " " << allhits.size() << "\n";
       for (unsigned int ii = 0; ii < hchk.size(); ++ii)
         if (!hchk[ii])
           std::cout << " " << ii << " " << allhits[ii]->WireID().Plane << " "
                     << allhits[ii]->WireID().Wire << " " << (int)allhits[ii]->PeakTime() << "\n";
       exit(1);
     }
 
   } // FillWireHitRange

void trkf::CCTrackMaker::FindMaybeVertices ( geo::TPCID const & tpcid )

private

Definition at line 1391 of file CCTrackMaker_module.cc.

References cls, util::end(), fWirePitch, geo::Iterable< IterationPolicy, Transform >::Iterate(), vtx, wireReadoutGeom, and X.

Referenced by produce().

   {
     // Project clusters to vertices and fill mVtxIndex. No requirement is
     // made that charge exists on the line between the Begin (End) of the
     // cluster and the vertex
 
     if (vtx.empty()) return;
 
     for (auto const& plane : wireReadoutGeom.Iterate<geo::PlaneGeo>(tpcid)) {
       auto const& planeID = plane.ID();
       unsigned int const ipl = planeID.Cryostat;
       for (unsigned int icl = 0; icl < cls[ipl].size(); ++icl) {
         for (unsigned int end = 0; end < 2u; ++end) {
           // ignore already attached clusters
           if (cls[ipl][icl].VtxIndex[end] >= 0) continue;
           short ibstvx = -1;
           float best = 1.;
           // index of the other end
           unsigned int oend = 1 - end;
           for (std::size_t ivx = 0; ivx < vtx.size(); ++ivx) {
             // ignore if the other end is attached to this vertex (which can happen with short clusters)
             if (cls[ipl][icl].VtxIndex[oend] == static_cast<short>(ivx)) continue;
             float const dWire = plane.WireCoordinate(geo::Point_t{0, vtx[ivx].Y, vtx[ivx].Z}) -
                                 cls[ipl][icl].Wire[end];
             if (end == 0) {
               if (dWire > 30 || dWire < -2) continue;
             }
             else {
               if (dWire < -30 || dWire > 2) continue;
             }
             // project the cluster to the vertex wire
             float const dX = fabs(cls[ipl][icl].X[end] +
                                   cls[ipl][icl].Slope[end] * fWirePitch * dWire - vtx[ivx].X);
             if (dX < best) {
               best = dX;
               ibstvx = ivx;
             }
           } // ivx
           if (ibstvx >= 0) {
             // attach
             cls[ipl][icl].VtxIndex[end] = ibstvx;
             cls[ipl][icl].mVtxIndex[end] = ibstvx;
           }
         } // end
       }   // icl
     }     // ipl
 
   } // FindMaybeVertices

bool trkf::CCTrackMaker::FindMissingCluster	(	unsigned short	kpl,
		short &	kcl,
		unsigned short &	kend,
		float	kWir,
		float	kX,
		float	okWir,
		float	okX
	)

private

Definition at line 2987 of file CCTrackMaker_module.cc.

References algIndex, clsChain, util::end(), fXMatchErr, and X.

Referenced by FillEndMatch().

   {
     // try to attach a missing cluster to the cluster chain kcl. kend is the "match end"
 
     unsigned short okend;
     float dxcut;
 
     if (kcl >= 0) return false;
 
     // Look for a missing cluster with loose cuts
     float kslp = fabs((okX - kX) / (okWir - kWir));
     if (kslp > 20) kslp = 20;
     // expand dX cut assuming there is a missing hit on the end of a cluster => 1 wire
     dxcut = 3 * fXMatchErr[algIndex] + kslp;
     unsigned short nfound = 0;
     unsigned short foundCl = 0, foundEnd = 0;
     for (unsigned short ccl = 0; ccl < clsChain[kpl].size(); ++ccl) {
       if (clsChain[kpl][ccl].InTrack >= 0) continue;
       // require a match at both ends
       for (unsigned short end = 0; end < 2; ++end) {
         okend = 1 - end;
         if (fabs(clsChain[kpl][ccl].Wire[end] - kWir) > 4) continue;
         if (fabs(clsChain[kpl][ccl].Wire[okend] - okWir) > 4) continue;
         // require at least one end to match
         if (fabs(clsChain[kpl][ccl].X[end] - kX) > dxcut &&
             fabs(clsChain[kpl][ccl].X[okend] - okX) > dxcut)
           continue;
         ++nfound;
         foundCl = ccl;
         foundEnd = end;
       } // end
     }   // ccl
     if (nfound == 0) return false;
     if (nfound > 1) {
       mf::LogVerbatim("CCTM") << "FindMissingCluster: Found too many matches. Write some code "
                               << nfound;
       return false;
     }
     kcl = foundCl;
     kend = foundEnd;
     return true;
 
   } // FindMissingCluster

void trkf::CCTrackMaker::FitVertices	(	detinfo::DetectorPropertiesData const &	detProp,
		geo::TPCID const &	tpcid
	)

private

Definition at line 799 of file CCTrackMaker_module.cc.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), util::end(), fHitFitErrFac, fNVtxTrkHitsFit, fVertexFitAlg, geo::WireReadoutGeom::Nplanes(), util::size(), trkf::VertexFitAlg::VertexFit(), vtx, and wireReadoutGeom.

Referenced by produce().

   {
     std::vector<std::vector<geo::WireID>> hitWID;
     std::vector<std::vector<double>> hitX;
     std::vector<std::vector<double>> hitXErr;
     TVector3 pos, posErr;
     std::vector<TVector3> trkDir;
     std::vector<TVector3> trkDirErr;
     float ChiDOF;
 
     if (fNVtxTrkHitsFit == 0) return;
 
     unsigned short indx, indx2, iht, nHitsFit;
 
     for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
       if (!vtx[ivx].Neutrino) continue;
       hitWID.clear();
       hitX.clear();
       hitXErr.clear();
       trkDir.clear();
       // find the tracks associated with this vertex
       unsigned int thePln, theTPC, theCst;
       for (unsigned short itk = 0; itk < trk.size(); ++itk) {
         for (unsigned short end = 0; end < 2; ++end) {
           if (trk[itk].VtxIndex[end] != ivx) continue;
           unsigned short itj = 0;
           if (end == 1) itj = trk[itk].TrjPos.size() - 1;
           // increase the size of the hit vectors by 1 for this track
           indx = hitX.size();
           hitWID.resize(indx + 1);
           hitX.resize(indx + 1);
           hitXErr.resize(indx + 1);
           trkDir.resize(indx + 1);
           trkDir[indx] = trk[itk].TrjDir[itj];
           auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
           for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
             if (trk[itk].TrkHits[ipl].size() < 2) continue;
             // make slots for the hits on this track in this plane
             nHitsFit = trk[itk].TrkHits[ipl].size();
             if (nHitsFit > fNVtxTrkHitsFit) nHitsFit = fNVtxTrkHitsFit;
             indx2 = hitWID[indx].size();
             hitWID[indx].resize(indx2 + nHitsFit);
             hitX[indx].resize(indx2 + nHitsFit);
             hitXErr[indx].resize(indx2 + nHitsFit);
             for (unsigned short ii = 0; ii < nHitsFit; ++ii) {
               if (end == 0) { iht = ii; }
               else {
                 iht = trk[itk].TrkHits[ipl].size() - ii - 1;
               }
               hitWID[indx][indx2 + ii] = trk[itk].TrkHits[ipl][iht]->WireID();
               thePln = trk[itk].TrkHits[ipl][iht]->WireID().Plane;
               theTPC = trk[itk].TrkHits[ipl][iht]->WireID().TPC;
               theCst = trk[itk].TrkHits[ipl][iht]->WireID().Cryostat;
               hitX[indx][indx2 + ii] = detProp.ConvertTicksToX(
                 trk[itk].TrkHits[ipl][iht]->PeakTime(), thePln, theTPC, theCst);
               hitXErr[indx][indx2 + ii] = fHitFitErrFac * trk[itk].TrkHits[ipl][iht]->RMS();
             } // ii
           }   // ipl
         }     // end
       }       // itk
       if (hitX.size() < 2) {
         mf::LogVerbatim("CCTM") << "Not enough hits to fit vtx " << ivx;
         continue;
       } // hitX.size() < 2
       pos(0) = vtx[ivx].X;
       pos(1) = vtx[ivx].Y;
       pos(2) = vtx[ivx].Z;
       fVertexFitAlg.VertexFit(hitWID, hitX, hitXErr, pos, posErr, trkDir, trkDirErr, ChiDOF);
       if (ChiDOF > 3000) continue;
       // update the vertex position
       vtx[ivx].X = pos(0);
       vtx[ivx].Y = pos(1);
       vtx[ivx].Z = pos(2);
       // and the track trajectory
       unsigned short fitTrk = 0;
       for (unsigned short itk = 0; itk < trk.size(); ++itk) {
         for (unsigned short end = 0; end < 2; ++end) {
           if (trk[itk].VtxIndex[end] != ivx) continue;
           unsigned short itj = 0;
           if (end == 1) itj = trk[itk].TrjPos.size() - 1;
           trk[itk].TrjDir[itj] = trkDir[fitTrk];
           ++fitTrk;
         } // end
       }   // itk
     }     // ivx
   }       // FitVertices

std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const

inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

   {
     return collector_.getConsumables();
   }

void trkf::CCTrackMaker::MakeClusterChains	(	detinfo::DetectorPropertiesData const &	detProp,
		art::FindManyP< recob::Hit > const &	fmCluHits,
		geo::TPCID const &	tpcid
	)

private

Definition at line 1441 of file CCTrackMaker_module.cc.

Referenced by produce().

   {
     unsigned short icl, icl1, icl2;
     float dw, dx, dWCut, dw1Max, dw2Max;
     float dA, dA2, dACut = fMaxDAng, chgAsymCut;
     float dXCut, chgasym, mrgErr;
     // long straight clusters
     bool ls1, ls2;
     bool gotprt = false;
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
       if (cls[ipl].size() > 1) {
         for (icl1 = 0; icl1 < cls[ipl].size() - 1; ++icl1) {
           prt = (fDebugAlg == 666 && ipl == fDebugPlane && icl1 == fDebugCluster);
           if (prt) gotprt = true;
           // maximum delta Wire overlap is 10% of the total length
           dw1Max = 0.6 * cls[ipl][icl1].Length;
           ls1 = (cls[ipl][icl1].Length > 100 &&
                  fabs(cls[ipl][icl1].Angle[0] - cls[ipl][icl1].Angle[1]) < 0.04);
           for (icl2 = icl1 + 1; icl2 < cls[ipl].size(); ++icl2) {
             ls2 = (cls[ipl][icl2].Length > 100 &&
                    fabs(cls[ipl][icl2].Angle[0] - cls[ipl][icl2].Angle[1]) < 0.04);
             dw2Max = 0.6 * cls[ipl][icl2].Length;
             // set overlap cut to be the shorter of the two
             dWCut = dw1Max;
             if (dw2Max < dWCut) dWCut = dw2Max;
             // but not exceeding 20 for very long clusters
             if (dWCut > 100) dWCut = 100;
             if (dWCut < 2) dWCut = 2;
             chgAsymCut = fMergeChgAsym;
             // Compare end 1 of icl1 with end 0 of icl2
 
             if (prt)
               mf::LogVerbatim("CCTM")
                 << "MCC P:C:W icl1 " << ipl << ":" << icl1 << ":" << cls[ipl][icl1].Wire[1]
                 << " vtx " << cls[ipl][icl1].VtxIndex[1] << " ls1 " << ls1 << " icl2 " << ipl << ":"
                 << icl2 << ":" << cls[ipl][icl2].Wire[0] << " vtx " << cls[ipl][icl2].VtxIndex[0]
                 << " ls2 " << ls2 << " dWCut " << dWCut;
             if (std::abs(cls[ipl][icl1].Wire[1] - cls[ipl][icl2].Wire[0]) > dWCut) continue;
             // ignore if the clusters begin/end on the same wire
             //            if(cls[ipl][icl1].Wire[1]  == cls[ipl][icl2].Wire[0]) continue;
             // or if the angle exceeds the cut
             float af = AngleFactor(cls[ipl][icl1].Slope[1]);
             dACut = fMaxDAng * af;
             dXCut = fChainMaxdX * 5 * af;
             dA = fabs(cls[ipl][icl1].Angle[1] - cls[ipl][icl2].Angle[0]);
             // compare the match angle at the opposite ends of the clusters.
             // May have a bad end/begin angle if there is a delta-ray in the middle
             dA2 = fabs(cls[ipl][icl1].Angle[0] - cls[ipl][icl2].Angle[1]);
 
             if (prt)
               mf::LogVerbatim("CCTM")
                 << " dA " << dA << " dA2 " << dA2 << " DACut " << dACut << " dXCut " << dXCut;
 
             if (dA2 < dA) dA = dA2;
             // ignore vertices that have only two associated clusters and
             // the angle is small
             if (dA < fChainVtxAng && cls[ipl][icl1].VtxIndex[1] >= 0) {
               dw = fWirePitch * (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1]);
               dx = cls[ipl][icl1].X[1] + cls[ipl][icl1].Slope[1] * dw * fWirePitch -
                    cls[ipl][icl2].X[0];
               unsigned short ivx = cls[ipl][icl1].VtxIndex[1];
               if (vtx[ivx].nClusInPln[ipl] == 2 && fabs(dx) < 1) {
                 cls[ipl][icl1].VtxIndex[1] = -2;
                 cls[ipl][icl2].VtxIndex[0] = -2;
                 vtx[ivx].nClusInPln[ipl] = 0;
                 if (prt) mf::LogVerbatim("CCTM") << " clobbered vertex " << ivx;
               } // vertices match
             }   // dA < 0.1 && ...
 
             // don't merge if a vertex exists at these ends
             if (cls[ipl][icl1].VtxIndex[1] >= 0) continue;
             if (cls[ipl][icl2].VtxIndex[0] >= 0) continue;
 
             // expand the angle match cut for clusters that appear to be stopping
             if (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1] < 3 &&
                 (cls[ipl][icl1].Length < 3 || cls[ipl][icl2].Length < 3)) {
               if (prt) mf::LogVerbatim("CCTM") << "Stopping cluster";
               dACut *= 1.5;
               chgAsymCut *= 1.5;
               dXCut *= 3;
             } // stopping US cluster
 
             // find the angle made by the endpoints of the clusters
             dw = fWirePitch * (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1]);
             if (dw != 0) {
               dx = cls[ipl][icl2].X[0] - cls[ipl][icl1].X[1];
               float dA3 = std::abs(atan(dx / dw) - cls[ipl][icl1].Angle[1]);
               if (prt) mf::LogVerbatim("CCTM") << " dA3 " << dA3;
               if (dA3 > dA) dA = dA3;
             }
 
             // angle matching
             if (dA > dACut) continue;
 
             if (prt)
               mf::LogVerbatim("CCTM")
                 << " rough dX " << fabs(cls[ipl][icl1].X[1] - cls[ipl][icl2].X[0]) << " cut = 20";
 
             // make a rough dX cut
             if (fabs(cls[ipl][icl1].X[1] - cls[ipl][icl2].X[0]) > 20) continue;
 
             // handle cosmic ray clusters that are broken at delta rays
             if (ls1 || ls2) {
               // tighter angle cuts but no charge cut
               if (dA > fChainVtxAng) continue;
             }
             else {
               chgasym = fabs(cls[ipl][icl1].Charge[1] - cls[ipl][icl2].Charge[0]);
               chgasym /= cls[ipl][icl1].Charge[1] + cls[ipl][icl2].Charge[0];
               if (prt) mf::LogVerbatim("CCTM") << " chgasym " << chgasym << " cut " << chgAsymCut;
               if (chgasym > chgAsymCut) continue;
             } // ls1 || ls2
             // project the longer cluster to the end of the shorter one
             if (cls[ipl][icl1].Length > cls[ipl][icl2].Length) {
               dw = fWirePitch * (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1]);
               dx = cls[ipl][icl1].X[1] + cls[ipl][icl1].Slope[1] * dw * fWirePitch -
                    cls[ipl][icl2].X[0];
             }
             else {
               dw = fWirePitch * (cls[ipl][icl1].Wire[1] - cls[ipl][icl2].Wire[0]);
               dx = cls[ipl][icl2].X[0] + cls[ipl][icl2].Slope[0] * dw * fWirePitch -
                    cls[ipl][icl1].X[1];
             }
 
             // handle overlapping clusters
             if (dA2 < 0.01 && abs(dx) > dXCut && dx < -1) {
               dx = dXClTraj(fmCluHits, ipl, icl1, 1);
               if (prt) mf::LogVerbatim("CCTM") << " new dx from dXClTraj " << dx;
             }
 
             if (prt)
               mf::LogVerbatim("CCTM")
                 << " X0 " << cls[ipl][icl1].X[1] << " slp " << cls[ipl][icl1].Slope[1] << " dw "
                 << dw << " oX " << cls[ipl][icl2].X[0] << " dx " << dx << " cut " << dXCut;
 
             if (fabs(dx) > dXCut) continue;
 
             // calculate a merge error that will be used to adjudicate between multiple merge attempts AngleFactor
             float xerr = dx / dXCut;
             float aerr = dA / dACut;
             mrgErr = xerr * xerr + aerr * aerr;
 
             if (prt)
               mf::LogVerbatim("CCTM")
                 << "icl1 mrgErr " << mrgErr << " MergeError " << cls[ipl][icl1].MergeError[1]
                 << " icl2 MergeError " << cls[ipl][icl2].MergeError[0];
 
             // this merge better than a previous one?
             if (mrgErr > cls[ipl][icl1].MergeError[1]) continue;
             if (mrgErr > cls[ipl][icl2].MergeError[0]) continue;
 
             // un-merge icl1 - this should always be true but check anyway
             if (cls[ipl][icl1].BrkIndex[1] >= 0) {
               unsigned short ocl = cls[ipl][icl1].BrkIndex[1];
               if (prt) mf::LogVerbatim("CCTM") << "clobber old icl1 BrkIndex " << ocl;
               if (cls[ipl][ocl].BrkIndex[0] == icl1) {
                 cls[ipl][ocl].BrkIndex[0] = -1;
                 cls[ipl][ocl].MergeError[0] = fMaxMergeError;
               }
               if (cls[ipl][ocl].BrkIndex[1] == icl1) {
                 cls[ipl][ocl].BrkIndex[1] = -1;
                 cls[ipl][ocl].MergeError[1] = fMaxMergeError;
               }
             } // cls[ipl][icl1].BrkIndex[1] >= 0
             cls[ipl][icl1].BrkIndex[1] = icl2;
             cls[ipl][icl1].MergeError[1] = mrgErr;
 
             // un-merge icl2
             if (cls[ipl][icl2].BrkIndex[0] >= 0) {
               unsigned short ocl = cls[ipl][icl2].BrkIndex[0];
               if (prt) mf::LogVerbatim("CCTM") << "clobber old icl2 BrkIndex " << ocl;
               if (cls[ipl][ocl].BrkIndex[0] == icl1) {
                 cls[ipl][ocl].BrkIndex[0] = -1;
                 cls[ipl][ocl].MergeError[0] = fMaxMergeError;
               }
               if (cls[ipl][ocl].BrkIndex[1] == icl1) {
                 cls[ipl][ocl].BrkIndex[1] = -1;
                 cls[ipl][ocl].MergeError[1] = fMaxMergeError;
               }
             } // cls[ipl][icl2].BrkIndex[0] >= 0
             cls[ipl][icl2].BrkIndex[0] = icl1;
             cls[ipl][icl2].MergeError[0] = mrgErr;
             if (prt) mf::LogVerbatim("CCTM") << " merge " << icl1 << " and " << icl2;
 
           } // icl2
         }   // icl1
 
         // look for broken clusters in which have a C shape similar to a Begin-Begin vertex. The clusters
         // will have large and opposite sign angles
         bool gotone;
         for (icl1 = 0; icl1 < cls[ipl].size() - 1; ++icl1) {
           gotone = false;
           for (icl2 = icl1 + 1; icl2 < cls[ipl].size(); ++icl2) {
             // check both ends
             for (unsigned short end = 0; end < 2; ++end) {
               // Ignore already identified broken clusters
               if (cls[ipl][icl1].BrkIndex[end] >= 0) continue;
               if (cls[ipl][icl2].BrkIndex[end] >= 0) continue;
               // require a large angle cluster
               if (fabs(cls[ipl][icl1].Angle[end]) < 1) continue;
               // and a second large angle cluster
               if (fabs(cls[ipl][icl2].Angle[end]) < 1) continue;
               if (prt)
                 mf::LogVerbatim("CCTM")
                   << "BrokenC: clusters " << cls[ipl][icl1].Wire[end] << ":"
                   << (int)cls[ipl][icl1].Time[end] << " " << cls[ipl][icl2].Wire[end] << ":"
                   << (int)cls[ipl][icl2].Time[end] << " angles " << cls[ipl][icl1].Angle[end] << " "
                   << cls[ipl][icl2].Angle[end] << " dWire "
                   << fabs(cls[ipl][icl1].Wire[end] - cls[ipl][icl2].Wire[end]);
               if (fabs(cls[ipl][icl1].Wire[end] - cls[ipl][icl2].Wire[end]) > 5) continue;
               // This is really crude but maybe OK
               // project 1 -> 2
               float dsl = cls[ipl][icl2].Slope[end] - cls[ipl][icl1].Slope[end];
               float fvw =
                 (cls[ipl][icl1].X[end] - cls[ipl][icl1].Wire[end] * cls[ipl][icl1].Slope[end] -
                  cls[ipl][icl2].X[end] + cls[ipl][icl2].Wire[end] * cls[ipl][icl2].Slope[end]) /
                 dsl;
               if (prt) mf::LogVerbatim("CCTM") << " fvw " << fvw;
               if (fabs(cls[ipl][icl1].Wire[end] - fvw) > 4) continue;
               if (fabs(cls[ipl][icl2].Wire[end] - fvw) > 4) continue;
               cls[ipl][icl1].BrkIndex[end] = icl2;
               // TODO This could use some improvement if necessary
               cls[ipl][icl1].MergeError[end] = 1;
               cls[ipl][icl2].BrkIndex[end] = icl1;
               cls[ipl][icl2].MergeError[end] = 1;
               gotone = true;
               dx = fabs(cls[ipl][icl1].X[end] - cls[ipl][icl2].X[end]);
               if (prt)
                 mf::LogVerbatim("CCTM")
                   << "BrokenC: icl1:W " << icl1 << ":" << cls[ipl][icl1].Wire[end] << " icl2:W "
                   << icl2 << ":" << cls[ipl][icl2].Wire[end] << " end " << end << " dx " << dx;
             } // end
             if (gotone) break;
           } // icl2
         }   // icl1
 
       } // cls[ipl].size() > 1
 
       // follow mother-daughter broken clusters and put them in the cluster chain array
       unsigned short mom, momBrkEnd, dtrBrkEnd, nit;
       short dtr;
 
       std::vector<bool> gotcl(cls[ipl].size());
       for (icl = 0; icl < cls[ipl].size(); ++icl)
         gotcl[icl] = false;
       if (prt)
         mf::LogVerbatim("CCTM") << "ipl " << ipl << " cls.size() " << cls[ipl].size() << "\n";
 
       std::vector<unsigned short> sCluster;
       std::vector<unsigned short> sOrder;
       for (icl = 0; icl < cls[ipl].size(); ++icl) {
         sCluster.clear();
         sOrder.clear();
         if (gotcl[icl]) continue;
         // don't start with a cluster broken at both ends
         if (cls[ipl][icl].BrkIndex[0] >= 0 && cls[ipl][icl].BrkIndex[1] >= 0) continue;
         for (unsigned int end = 0; end < 2; ++end) {
           if (cls[ipl][icl].BrkIndex[end] < 0) continue;
           if (cls[ipl][icl].MergeError[end] > fMergeErrorCut) continue;
           gotcl[icl] = true;
           mom = icl;
           // end where the mom is broken
           momBrkEnd = end;
           sCluster.push_back(mom);
           if (momBrkEnd == 1) {
             // typical case - broken at the DS end
             sOrder.push_back(0);
           }
           else {
             // broken at the US end
             sOrder.push_back(1);
           }
           dtr = cls[ipl][icl].BrkIndex[end];
           nit = 0;
           while (dtr >= 0 && dtr < (short)cls[ipl].size() && nit < cls[ipl].size()) {
             // determine which end of the dtr should be attached to mom
             for (dtrBrkEnd = 0; dtrBrkEnd < 2; ++dtrBrkEnd)
               if (cls[ipl][dtr].BrkIndex[dtrBrkEnd] == mom) break;
             if (dtrBrkEnd == 2) {
               gotcl[icl] = false;
               break;
             }
             // check for reasonable merge error
             if (cls[ipl][dtr].MergeError[dtrBrkEnd] < fMergeErrorCut) {
               sCluster.push_back(dtr);
               sOrder.push_back(dtrBrkEnd);
               gotcl[dtr] = true;
             }
             ++nit;
             // set up to check the new mom
             mom = dtr;
             // at the other end
             momBrkEnd = 1 - dtrBrkEnd;
             // with the new dtr
             dtr = cls[ipl][mom].BrkIndex[momBrkEnd];
           } // dtr >= 0 ...
           if (dtrBrkEnd == 2) continue;
         } // end
 
         if (!gotcl[icl]) {
           // a single unbroken cluster
           sCluster.push_back(icl);
           sOrder.push_back(0);
         }
 
         if (sCluster.size() == 0) {
           mf::LogError("CCTM") << "MakeClusterChains error in plane " << ipl << " cluster " << icl;
           return;
         }
 
         ClsChainPar ccp;
         // fill the struct parameters assuming this is a cluster chain containing only one cluster
         unsigned short jcl = sCluster[0];
         if (jcl > cls[ipl].size()) std::cout << "oops MCC\n";
         unsigned short oend;
         for (unsigned int end = 0; end < 2; ++end) {
           oend = end;
           if (sOrder[0] > 0) oend = 1 - end;
           ccp.Wire[end] = cls[ipl][jcl].Wire[oend];
           ccp.Time[end] = cls[ipl][jcl].Time[oend];
           ccp.X[end] = cls[ipl][jcl].X[oend];
           ccp.Slope[end] = cls[ipl][jcl].Slope[oend];
           ccp.Angle[end] = cls[ipl][jcl].Angle[oend];
           ccp.Dir[end] = cls[ipl][icl].Dir[oend];
           ccp.VtxIndex[end] = cls[ipl][jcl].VtxIndex[oend];
           ccp.ChgNear[end] = cls[ipl][jcl].ChgNear[oend];
           ccp.mBrkIndex[end] = cls[ipl][jcl].BrkIndex[oend];
         } // end
         ccp.Length = cls[ipl][icl].Length;
         ccp.TotChg = cls[ipl][icl].TotChg;
         ccp.InTrack = -1;
 
         for (unsigned short ii = 1; ii < sCluster.size(); ++ii) {
           jcl = sCluster[ii];
           if (jcl > cls[ipl].size()) std::cout << "oops MCC\n";
           // end is the end where the break is being mended
           unsigned int end = sOrder[ii];
           if (end > 1) std::cout << "oops2 MCC\n";
           oend = 1 - end;
           // update the parameters at the other end of the chain
           ccp.Wire[1] = cls[ipl][jcl].Wire[oend];
           ccp.Time[1] = cls[ipl][jcl].Time[oend];
           ccp.X[1] = cls[ipl][jcl].X[oend];
           ccp.Slope[1] = cls[ipl][jcl].Slope[oend];
           ccp.Angle[1] = cls[ipl][jcl].Angle[oend];
           ccp.Dir[1] = cls[ipl][jcl].Dir[oend];
           ccp.VtxIndex[1] = cls[ipl][jcl].VtxIndex[oend];
           ccp.ChgNear[1] = cls[ipl][jcl].ChgNear[oend];
           ccp.mBrkIndex[1] = cls[ipl][jcl].BrkIndex[oend];
           ccp.Length += cls[ipl][jcl].Length;
           ccp.TotChg += cls[ipl][jcl].TotChg;
         } // ii
         ccp.ClsIndex = sCluster;
         ccp.Order = sOrder;
         // redo the direction
         if (ccp.Time[1] > ccp.Time[0]) {
           ccp.Dir[0] = 1;
           ccp.Dir[1] = -1;
         }
         else {
           ccp.Dir[0] = -1;
           ccp.Dir[1] = 1;
         }
         clsChain[ipl].push_back(ccp);
 
       } // icl
 
       // re-index mBrkIndex to point to a cluster chain, not a cluster
       unsigned short brkCls;
       bool gotit;
       for (unsigned short ccl = 0; ccl < clsChain[ipl].size(); ++ccl) {
         for (unsigned short end = 0; end < 2; ++end) {
           if (clsChain[ipl][ccl].mBrkIndex[end] < 0) continue;
           brkCls = clsChain[ipl][ccl].mBrkIndex[end];
           gotit = false;
           // find this cluster index in a cluster chain
           for (unsigned short ccl2 = 0; ccl2 < clsChain[ipl].size(); ++ccl2) {
             if (ccl2 == ccl) continue;
             if (std::find(clsChain[ipl][ccl2].ClsIndex.begin(),
                           clsChain[ipl][ccl2].ClsIndex.end(),
                           brkCls) == clsChain[ipl][ccl2].ClsIndex.end())
               continue;
             // found it
             clsChain[ipl][ccl].mBrkIndex[end] = ccl2;
             gotit = true;
             break;
           } // ccl2
           if (!gotit)
             mf::LogError("CCTM") << "MCC: Cluster chain " << ccl << " end " << end
                                  << " Failed to find brkCls " << brkCls << " in plane " << ipl;
         } // end
       }   // ccl
 
     } // ipl
 
     if (gotprt) PrintClusters(detProp, tpcid);
     prt = false;
 
   } // MakeClusterChains

void trkf::CCTrackMaker::MakeFamily ( geo::TPCID const & tpcid )

private

Definition at line 944 of file CCTrackMaker_module.cc.

References tca::Length(), geo::WireReadoutGeom::Nplanes(), pfpToTrkID, util::size(), tmp, vtx, wireReadoutGeom, X, Y, and Z.

Referenced by produce().

   {
     // define the track/vertex and mother/daughter relationships
 
     unsigned short ivx, itr, ipl, ii, jtr;
     unsigned short nus, nds, nuhs, ndhs;
     float longUSTrk, longDSTrk, qual;
 
     // min distance^2 between a neutrino vertex candidate and a through
     // going muon
     float tgMuonCut2 = 9;
 
     // struct for neutrino vertex candidates
     struct NuVtx {
       unsigned short VtxIndex;
       unsigned short nUSTk;
       unsigned short nDSTk;
       unsigned short nUSHit;
       unsigned short nDSHit;
       float longUSTrk;
       float longDSTrk;
       float Qual;
     };
     std::vector<NuVtx> nuVtxCand;
 
     NuVtx aNuVtx;
 
     // analyze all of the vertices
     float best = 999, dx, dy, dz, dr;
     short imbest = -1;
     bool skipVtx;
     unsigned short itj;
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (ivx = 0; ivx < vtx.size(); ++ivx) {
       vtx[ivx].Neutrino = false;
       nus = 0;
       nds = 0;
       nuhs = 0;
       ndhs = 0;
       longUSTrk = 0;
       longDSTrk = 0;
       skipVtx = false;
       // skip vertices that are close to through-going muons
       for (itr = 0; itr < trk.size(); ++itr) {
         if (trk[itr].ID < 0) continue;
         if (trk[itr].PDGCode != 13) continue;
         for (itj = 0; itj < trk[itr].TrjPos.size(); ++itj) {
           dx = trk[itr].TrjPos[itj](0) - vtx[ivx].X;
           dy = trk[itr].TrjPos[itj](1) - vtx[ivx].Y;
           dz = trk[itr].TrjPos[itj](2) - vtx[ivx].Z;
           dr = dx * dx + dy * dy + dz * dz;
           if (dr < tgMuonCut2) {
             skipVtx = true;
             break;
           }
           if (skipVtx) break;
         } // itj
         if (skipVtx) break;
       } // itr
       if (skipVtx) continue;
       for (itr = 0; itr < trk.size(); ++itr) {
         if (trk[itr].ID < 0) continue;
         if (trk[itr].VtxIndex[0] == ivx) {
           // DS-going track
           ++nds;
           if (trk[itr].Length > longDSTrk) longDSTrk = trk[itr].Length;
           for (ipl = 0; ipl < nplanes; ++ipl)
             ndhs += trk[itr].TrkHits[ipl].size();
         } // DS-going track
         // Reject this vertex as a neutrino candidate if the track being
         // considered has a starting vertex
         if (trk[itr].VtxIndex[1] == ivx && trk[itr].VtxIndex[0] >= 0) {
           skipVtx = true;
           break;
         } // trk[itr].VtxIndex[0] > 0
         if (trk[itr].VtxIndex[1] == ivx && trk[itr].VtxIndex[0] < 0) {
           // US-going track w no US vertex
           ++nus;
           if (trk[itr].Length > longUSTrk) longUSTrk = trk[itr].Length;
           for (ipl = 0; ipl < nplanes; ++ipl)
             nuhs += trk[itr].TrkHits[ipl].size();
         } // US-going track
       }   // itr
       if (skipVtx) continue;
       if (nds == 0) continue;
       qual = 1 / (float)nds;
       qual /= (float)ndhs;
       if (nus > 0) qual *= (float)nuhs / (float)ndhs;
       if (qual < best) {
         best = qual;
         imbest = ivx;
       }
       if (nds > 0 && longDSTrk > 5) {
         // at least one longish track going DS
         aNuVtx.VtxIndex = ivx;
         aNuVtx.nUSTk = nus;
         aNuVtx.nDSTk = nds;
         aNuVtx.nUSHit = nuhs;
         aNuVtx.nDSHit = ndhs;
         aNuVtx.longUSTrk = longUSTrk;
         aNuVtx.longDSTrk = longDSTrk;
         aNuVtx.Qual = qual;
         nuVtxCand.push_back(aNuVtx);
       }
     } // ivx
     if (imbest < 0) return;
 
     // Found the neutrino interaction vertex
     ivx = imbest;
     vtx[ivx].Neutrino = true;
     // Put a 0 into the PFParticle -> track vector to identify a neutrino
     // track with no trajectory or hits
     pfpToTrkID.push_back(0);
 
     // list of DS-going current generation daughters so we can fix next
     // generation daughter assignments. This code section assumes that there
     // are no decays or 2ndry interactions with US-going primary tracks.
     std::vector<unsigned short> dtrGen;
     std::vector<unsigned short> dtrNextGen;
     for (itr = 0; itr < trk.size(); ++itr) {
       if (trk[itr].ID < 0) continue;
       if (trk[itr].VtxIndex[0] == ivx) {
         // DS-going primary track
         trk[itr].MomID = 0;
         // call every track coming from a neutrino vertex a proton
         trk[itr].PDGCode = 2212;
         pfpToTrkID.push_back(trk[itr].ID);
         dtrGen.push_back(itr);
       } // DS-going primary track
       if (trk[itr].VtxIndex[1] == ivx) {
         // US-going primary track
         trk[itr].MomID = 0;
         // call every track coming from a neutrino vertex a proton
         trk[itr].PDGCode = 2212;
         pfpToTrkID.push_back(trk[itr].ID);
         // reverse the track trajectory
         std::reverse(trk[itr].TrjPos.begin(), trk[itr].TrjPos.end());
         for (ii = 0; ii < trk[itr].TrjDir.size(); ++ii)
           trk[itr].TrjDir[ii] = -trk[itr].TrjDir[ii];
       } // DS-going primary track
     }   // itr
 
     if (dtrGen.size() == 0) return;
 
     unsigned short tmp, indx;
     unsigned short nit = 0;
 
     // follow daughters through all generations (< 10). Daughter tracks
     // may go US or DS
     while (nit < 10) {
       ++nit;
       dtrNextGen.clear();
       // Look for next generation daughters
       for (ii = 0; ii < dtrGen.size(); ++ii) {
         itr = dtrGen[ii];
         if (trk[itr].VtxIndex[1] >= 0) {
           // found a DS vertex
           ivx = trk[itr].VtxIndex[1];
           // look for a track associated with this vertex
           for (jtr = 0; jtr < trk.size(); ++jtr) {
             if (jtr == itr) continue;
             if (trk[jtr].VtxIndex[0] == ivx) {
               // DS-going track
               indx = trk[itr].DtrID.size();
               trk[itr].DtrID.resize(indx + 1);
               trk[itr].DtrID[indx] = jtr;
               trk[jtr].MomID = trk[itr].ID;
               // call all daughters pions
               trk[jtr].PDGCode = 211;
               dtrNextGen.push_back(jtr);
               pfpToTrkID.push_back(trk[jtr].ID);
             } // DS-going track
             if (trk[jtr].VtxIndex[1] == ivx) {
               // US-going track
               indx = trk[itr].DtrID.size();
               trk[itr].DtrID.resize(indx + 1);
               trk[itr].DtrID[indx] = jtr;
               trk[jtr].MomID = trk[itr].ID;
               // call all daughters pions
               trk[jtr].PDGCode = 211;
               dtrNextGen.push_back(jtr);
               pfpToTrkID.push_back(trk[jtr].ID);
               // reverse the track trajectory
               std::reverse(trk[jtr].TrjPos.begin(), trk[jtr].TrjPos.end());
               for (unsigned short jj = 0; jj < trk[jtr].TrjDir.size(); ++jj)
                 trk[jtr].TrjDir[jj] = -trk[jtr].TrjDir[jj];
               // interchange the trk - vtx assignments
               tmp = trk[jtr].VtxIndex[0];
               trk[jtr].VtxIndex[0] = trk[jtr].VtxIndex[1];
               trk[jtr].VtxIndex[1] = tmp;
             } // DS-going track
           }   // jtr
         }     // trk[itr].VtxIndex[0] >= 0
       }       // ii (itr)
       // break out if no next gen daughters found
       if (dtrNextGen.size() == 0) break;
       dtrGen = dtrNextGen;
     } // nit
 
   } // MakeFamily

std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const & wp )

inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

   {
     return doMakeWorker(wp);
   }

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::mayConsume ( InputTag const & tag )

protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

   {
     return collector_.mayConsume<T, BT>(tag);
   }

template<typename T , BranchType BT>

void art::ModuleBase::mayConsumeMany ( )

protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

   {
     collector_.mayConsumeMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const & )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const & tag )

inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

   {
     return collector_.mayConsumeView<T, BT>(tag);
   }

ModuleDescription const & art::ModuleBase::moduleDescription ( ) const

inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

   {
     if (md_.has_value()) {
       return *md_;
     }
 
     throw Exception{errors::LogicError,
                     "There was an error while calling moduleDescription().\n"}
       << "The moduleDescription() base-class member function cannot be called\n"
          "during module construction.  To determine which module is "
          "responsible\n"
          "for calling it, find the '<module type>:<module "
          "label>@Construction'\n"
          "tag in the message prefix above.  Please contact artists@fnal.gov\n"
          "for guidance.\n";
   }

void trkf::CCTrackMaker::PlnMatch	(	detinfo::DetectorPropertiesData const &	detProp,
		geo::TPCID const &	tpcid
	)

private

Definition at line 2074 of file CCTrackMaker_module.cc.

References util::abs(), algIndex, AngleFactor(), ChargeAsym(), trkf::CCTrackMaker::MatchPars::Chg, trkf::CCTrackMaker::MatchPars::Cls, clsChain, DupMatch(), trkf::CCTrackMaker::MatchPars::End, trkf::CCTrackMaker::MatchPars::Err, fChgAsymFactor, FillEndMatch(), fMatchMinLen, fuBCode, geom, geo::TPCGeo::HalfWidth(), geo::WireIDIntersection::invalid(), geo::kX, tca::Length(), matcomb, geo::WireReadoutGeom::Nplanes(), trkf::CCTrackMaker::MatchPars::oErr, trkf::CCTrackMaker::MatchPars::oVtx, geo::WireReadoutGeom::Plane(), prt, geo::WireReadoutGeom::ThirdPlaneSlope(), geo::GeometryCore::TPC(), trkf::CCTrackMaker::MatchPars::Vtx, geo::PlaneGeo::WireCoordinate(), geo::WireReadoutGeom::WireIDsIntersect(), wireReadoutGeom, and X.

Referenced by produce().

   {
     // Match clusters in all planes
     bool ignoreSign;
     float kSlp, kAng, kX, kWir, okWir;
     short idir, ioend, jdir, joend, kdir;
 
     double yp, zp;
     auto const& first_tpc = geom->TPC({0, 0});
     float tpcSizeY = first_tpc.HalfWidth();
     float tpcSizeZ = first_tpc.Length();
 
     float dxcut = 2;
     float dxkcut;
     float dwcut = 6;
     if (fuBCode) {
       dxcut = 20;
       dwcut = 60;
     }
 
     // temp array for making a rough charge asymmetry cut
     std::array<float, 3> mchg;
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
       geo::PlaneID const iplane_id{tpcid, ipl};
       for (unsigned short icl = 0; icl < clsChain[ipl].size(); ++icl) {
         if (clsChain[ipl][icl].InTrack >= 0) continue;
         // skip short clusters
         if (clsChain[ipl][icl].Length < fMatchMinLen[algIndex]) continue;
         unsigned short jpl = (ipl + 1) % nplanes;
         unsigned short kpl = (jpl + 1) % nplanes;
         geo::PlaneID const jplane_id{tpcid, jpl};
         for (unsigned short jcl = 0; jcl < clsChain[jpl].size(); ++jcl) {
           if (clsChain[jpl][jcl].InTrack >= 0) continue;
           // skip short clusters
           if (clsChain[jpl][jcl].Length < fMatchMinLen[algIndex]) continue;
           // make first charge asymmetry cut
           mchg[0] = clsChain[ipl][icl].TotChg;
           mchg[1] = clsChain[jpl][jcl].TotChg;
           mchg[2] = mchg[1];
           if (fChgAsymFactor[algIndex] > 0 && ChargeAsym(mchg) > 0.5) continue;
           for (unsigned short iend = 0; iend < 2; ++iend) {
             idir = clsChain[ipl][icl].Dir[iend];
             for (unsigned short jend = 0; jend < 2; ++jend) {
               jdir = clsChain[jpl][jcl].Dir[jend];
               if (idir != 0 && jdir != 0 && idir != jdir) continue;
               // make an X cut
               if (fabs(clsChain[ipl][icl].X[iend] - clsChain[jpl][jcl].X[jend]) > dxcut) continue;
               ioend = 1 - iend;
               joend = 1 - jend;
               // Find the expected third (k) plane parameters
               kSlp = wireReadoutGeom.ThirdPlaneSlope(
                 ipl, clsChain[ipl][icl].Slope[iend], jpl, clsChain[jpl][jcl].Slope[jend], tpcid);
               kAng = atan(kSlp);
 
               // Ensure the match end is within the TPC
               auto const intersection =
                 wireReadoutGeom
                   .WireIDsIntersect(
                     geo::WireID{iplane_id, (unsigned int)(0.5 + clsChain[ipl][icl].Wire[iend])},
                     geo::WireID{jplane_id, (unsigned int)(0.5 + clsChain[jpl][jcl].Wire[jend])})
                   .value_or(geo::WireIDIntersection::invalid());
               yp = intersection.y;
               zp = intersection.z;
               if (yp > tpcSizeY || yp < -tpcSizeY) continue;
               if (zp < 0 || zp > tpcSizeZ) continue;
 
               kX = 0.5 * (clsChain[ipl][icl].X[iend] + clsChain[jpl][jcl].X[jend]);
               auto const& kplane = wireReadoutGeom.Plane({tpcid, kpl});
               kWir = kplane.WireCoordinate(geo::Point_t{0, yp, zp});
 
               // now look at the other end
               auto const ointersection =
                 wireReadoutGeom
                   .WireIDsIntersect(
                     geo::WireID{iplane_id, (unsigned int)(0.5 + clsChain[ipl][icl].Wire[ioend])},
                     geo::WireID{jplane_id, (unsigned int)(0.5 + clsChain[jpl][jcl].Wire[joend])})
                   .value_or(geo::WireIDIntersection::invalid());
               yp = ointersection.y;
               zp = ointersection.z;
               if (yp > tpcSizeY || yp < -tpcSizeY) continue;
               if (zp < 0 || zp > tpcSizeZ) continue;
               okWir = kplane.WireCoordinate(geo::Point_t{0, yp, zp});
 
               if (prt)
                 mf::LogVerbatim("CCTM")
                   << "PlnMatch: chk i " << ipl << ":" << icl << ":" << iend << " idir " << idir
                   << " X " << clsChain[ipl][icl].X[iend] << " j " << jpl << ":" << jcl << ":"
                   << jend << " jdir " << jdir << " X " << clsChain[jpl][jcl].X[jend];
 
               if (prt)
                 mf::LogVerbatim("CCTM")
                   << "PlnMatch: chk j " << ipl << ":" << icl << ":" << iend << " " << jpl << ":"
                   << jcl << ":" << jend << " iSlp " << std::setprecision(2)
                   << clsChain[ipl][icl].Slope[iend] << " jSlp " << std::setprecision(2)
                   << clsChain[jpl][jcl].Slope[jend] << " kWir " << (int)kWir << " okWir "
                   << (int)okWir << " kSlp " << std::setprecision(2) << kSlp << " kAng "
                   << std::setprecision(2) << kAng << " kX " << std::setprecision(1) << kX;
 
               // handle the case near pi/2, where the errors on large slopes
               // could result in a wrong-sign kAng
               ignoreSign = (fabs(kSlp) > 1.5);
               if (ignoreSign) kAng = fabs(kAng);
               dxkcut = dxcut * AngleFactor(kSlp);
               bool gotkcl = false;
               for (unsigned short kcl = 0; kcl < clsChain[kpl].size(); ++kcl) {
                 if (clsChain[kpl][kcl].InTrack >= 0) continue;
                 // make second charge asymmetry cut
                 mchg[0] = clsChain[ipl][icl].TotChg;
                 mchg[1] = clsChain[jpl][jcl].TotChg;
                 mchg[2] = clsChain[kpl][kcl].TotChg;
                 if (fChgAsymFactor[algIndex] > 0 && ChargeAsym(mchg) > 0.5) continue;
                 for (unsigned short kend = 0; kend < 2; ++kend) {
                   kdir = clsChain[kpl][kcl].Dir[kend];
                   if (idir != 0 && kdir != 0 && idir != kdir) continue;
                   if (prt)
                     mf::LogVerbatim("CCTM")
                       << " kcl " << kcl << " kend " << kend << " dx "
                       << std::abs(clsChain[kpl][kcl].X[kend] - kX) << " dxkcut " << dxkcut;
                   if (std::abs(clsChain[kpl][kcl].X[kend] - kX) > dxkcut) continue;
                   // rough dWire cut
                   if (prt)
                     mf::LogVerbatim("CCTM")
                       << " kcl " << kcl << " kend " << kend << " dw "
                       << (clsChain[kpl][kcl].Wire[kend] - kWir) << " ignoreSign " << ignoreSign;
                   if (fabs(clsChain[kpl][kcl].Wire[kend] - kWir) > dwcut) continue;
                   if (prt) mf::LogVerbatim("CCTM") << " chk k " << kpl << ":" << kcl << ":" << kend;
                   MatchPars match;
                   match.Cls[ipl] = icl;
                   match.End[ipl] = iend;
                   match.Cls[jpl] = jcl;
                   match.End[jpl] = jend;
                   match.Cls[kpl] = kcl;
                   match.End[kpl] = kend;
                   match.Err = 100;
                   if (DupMatch(match, nplanes)) continue;
                   match.Chg[ipl] = 0;
                   match.Chg[jpl] = 0;
                   match.Chg[kpl] = 0;
                   match.Vtx = clsChain[ipl][icl].VtxIndex[iend];
                   match.oVtx = -1;
                   FillEndMatch(detProp, tpcid, match);
                   if (prt)
                     mf::LogVerbatim("CCTM") << " PlnMatch: match k " << kpl << ":" << match.Cls[kpl]
                                             << ":" << match.End[kpl] << " oChg " << match.Chg[kpl]
                                             << " mErr " << match.Err << " oErr " << match.oErr;
                   if (match.Chg[kpl] == 0) continue;
                   if (match.Err > 10 || match.oErr > 10) continue;
                   if (prt) mf::LogVerbatim("CCTM") << " dup? ";
                   if (DupMatch(match, nplanes)) continue;
                   matcomb.push_back(match);
                   gotkcl = true;
                 } // kend
               }   // kcl
               if (prt) mf::LogVerbatim("CCTM") << " PlnMatch: gotkcl " << gotkcl;
               if (!gotkcl) {
                 // make a 2-plane match and try again
                 MatchPars match;
                 match.Cls[ipl] = icl;
                 match.End[ipl] = iend;
                 match.Cls[jpl] = jcl;
                 match.End[jpl] = jend;
                 match.Cls[kpl] = -1;
                 match.End[kpl] = 0;
                 match.Err = 100;
                 if (DupMatch(match, nplanes)) continue;
                 match.Chg[ipl] = 0;
                 match.Chg[jpl] = 0;
                 match.Chg[kpl] = 0;
                 match.Vtx = clsChain[ipl][icl].VtxIndex[iend];
                 match.oVtx = -1;
                 FillEndMatch(detProp, tpcid, match);
                 if (prt)
                   mf::LogVerbatim("CCTM")
                     << " Tried 2-plane match"
                     << " k " << kpl << ":" << match.Cls[kpl] << ":" << match.End[kpl] << " Chg "
                     << match.Chg[kpl] << " Err " << match.Err << " match.oErr " << match.oErr;
                 if (match.Chg[kpl] <= 0) continue;
                 if (match.Err > 10 || match.oErr > 10) continue;
                 matcomb.push_back(match);
               } // !gotkcl
             }   // jend
           }     // iend
         }       // jcl
       }         // icl
     }           // ipl
   }             // PlnMatch

void trkf::CCTrackMaker::PrintClusters	(	detinfo::DetectorPropertiesData const &	detProp,
		geo::TPCID const &	tpcid
	)		const

private

Definition at line 3192 of file CCTrackMaker_module.cc.

References cls, clsChain, detinfo::DetectorPropertiesData::ConvertXToTicks(), util::end(), fPrintAllClusters, geo::Iterable< IterationPolicy, Transform >::Iterate(), geo::WireReadoutGeom::Nplanes(), art::right(), vtx, wireReadoutGeom, and X.

Referenced by MakeClusterChains(), and produce().

   {
 
     unsigned short iTime;
     mf::LogVerbatim myprt("CCTM");
     myprt << "******* PrintClusters *********  Num_Clusters_in             Wire:Time\n";
     myprt << "vtx  Index    X       Y       Z  Pln0 Pln1 Pln2          Pln0    Pln1    Pln2\n";
     for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
       myprt << std::right << std::setw(3) << ivx << std::setw(7) << ivx;
       myprt << std::fixed;
       myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].X;
       myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Y;
       myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Z;
       myprt << std::right << std::setw(5) << vtx[ivx].nClusInPln[0];
       myprt << std::right << std::setw(5) << vtx[ivx].nClusInPln[1];
       myprt << std::right << std::setw(5) << vtx[ivx].nClusInPln[2];
       myprt << "    ";
       for (auto const& plane : wireReadoutGeom.Iterate<geo::PlaneGeo>(tpcid)) {
         int time = 0.5 + detProp.ConvertXToTicks(vtx[ivx].X, plane.ID());
         int wire = plane.WireCoordinate(geo::Point_t{0, vtx[ivx].Y, vtx[ivx].Z});
         myprt << std::right << std::setw(7) << wire << ":" << time;
       }
 
       myprt << "\n";
     } // ivx
 
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
       myprt << ">>>>>>>>>> Cluster chains in Plane " << ipl << "\n";
       myprt << "ipl   ccl  Len    Chg    W0:T0     Ang0 Dir0  Vx0  mBk0   W1:T1     Ang1 Dir1  Vx1 "
                " mBk1  InTk    cls:Order \n";
       for (unsigned short ccl = 0; ccl < clsChain[ipl].size(); ++ccl) {
         myprt << std::right << std::setw(3) << ipl;
         myprt << std::right << std::setw(5) << ccl;
         myprt << std::right << std::setw(6) << clsChain[ipl][ccl].Length;
         myprt << std::right << std::setw(8) << (int)clsChain[ipl][ccl].TotChg;
         for (unsigned short end = 0; end < 2; ++end) {
           iTime = clsChain[ipl][ccl].Time[end];
           myprt << std::right << std::setw(5) << (int)clsChain[ipl][ccl].Wire[end] << ":"
                 << std::setprecision(1) << iTime;
           if (iTime < 10) { myprt << "   "; }
           else if (iTime < 100) {
             myprt << "  ";
           }
           else if (iTime < 1000)
             myprt << " ";
           myprt << std::right << std::setw(7) << std::setprecision(2)
                 << clsChain[ipl][ccl].Angle[end];
           myprt << std::right << std::setw(5) << clsChain[ipl][ccl].Dir[end];
           myprt << std::right << std::setw(5) << clsChain[ipl][ccl].VtxIndex[end];
           myprt << std::fixed << std::right << std::setw(6) << std::setprecision(1)
                 << clsChain[ipl][ccl].mBrkIndex[end];
         }
         myprt << std::right << std::setw(7) << clsChain[ipl][ccl].InTrack;
         myprt << "   ";
         for (unsigned short ii = 0; ii < clsChain[ipl][ccl].ClsIndex.size(); ++ii)
           myprt << " " << clsChain[ipl][ccl].ClsIndex[ii] << ":" << clsChain[ipl][ccl].Order[ii];
         myprt << "\n";
       } // ccl
       if (fPrintAllClusters) {
         myprt << ">>>>>>>>>> Clusters in Plane " << ipl << "\n";
         myprt << "ipl  icl  Evt   Len     Chg   W0:T0     Ang0 Dir0  Vx0  CN0   W1:T1      Ang1  "
                  "Dir1  Vx1  CN1 InTk Brk0 MrgEr0 Brk1 MrgEr1\n";
         for (unsigned short icl = 0; icl < cls[ipl].size(); ++icl) {
           myprt << std::right << std::setw(3) << ipl;
           myprt << std::right << std::setw(5) << icl;
           myprt << std::right << std::setw(5) << cls[ipl][icl].EvtIndex;
           myprt << std::right << std::setw(6) << cls[ipl][icl].Length;
           myprt << std::right << std::setw(8) << (int)cls[ipl][icl].TotChg;
           for (unsigned short end = 0; end < 2; ++end) {
             iTime = cls[ipl][icl].Time[end];
             myprt << std::right << std::setw(5) << (int)cls[ipl][icl].Wire[end] << ":" << iTime;
             if (iTime < 10) { myprt << "   "; }
             else if (iTime < 100) {
               myprt << "  ";
             }
             else if (iTime < 1000)
               myprt << " ";
             myprt << std::right << std::setw(7) << std::setprecision(2) << cls[ipl][icl].Angle[end];
             myprt << std::right << std::setw(5) << cls[ipl][icl].Dir[end];
             myprt << std::right << std::setw(5) << cls[ipl][icl].VtxIndex[end];
             myprt << std::fixed << std::right << std::setw(5) << std::setprecision(1)
                   << cls[ipl][icl].ChgNear[end];
           }
           myprt << std::fixed;
           myprt << std::right << std::setw(5) << cls[ipl][icl].InTrack;
           myprt << std::right << std::setw(5) << (int)cls[ipl][icl].BrkIndex[0];
           myprt << std::right << std::setw(7) << std::setprecision(1)
                 << cls[ipl][icl].MergeError[0];
           myprt << std::right << std::setw(5) << (int)cls[ipl][icl].BrkIndex[1];
           myprt << std::right << std::setw(7) << std::setprecision(1)
                 << cls[ipl][icl].MergeError[1];
           myprt << "\n";
         } // icl
       }   // fPrintAllClusters
     }     // ipl
   }       // PrintClusters

void trkf::CCTrackMaker::PrintTracks ( ) const

private

Definition at line 3143 of file CCTrackMaker_module.cc.

References art::right(), util::size(), and vtx.

Referenced by produce().

   {
     mf::LogVerbatim myprt("CCTM");
     myprt << "********* PrintTracks \n";
     myprt << "vtx  Index    X      Y      Z\n";
     for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
       myprt << std::right << std::setw(4) << ivx << std::setw(4) << vtx[ivx].EvtIndex;
       myprt << std::fixed;
       myprt << std::right << std::setw(10) << std::setprecision(1) << vtx[ivx].X;
       myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Y;
       myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Z;
       if (vtx[ivx].Neutrino) myprt << " Neutrino vertex";
       myprt << "\n";
     } // ivx
 
     myprt << ">>>>>>>>>> Tracks \n";
     myprt << "trk  ID  Proc nht nTrj  sX     sY     sZ     eX     eY     eZ  sVx eVx sGd eGd "
              "ChgOrd  dirZ Mom PDG     ClsIndices\n";
     for (unsigned short itr = 0; itr < trk.size(); ++itr) {
       myprt << std::right << std::setw(3) << itr << std::setw(4) << trk[itr].ID;
       myprt << std::right << std::setw(5) << std::setw(4) << trk[itr].Proc;
       unsigned short nht = 0;
       for (unsigned short ii = 0; ii < 3; ++ii)
         nht += trk[itr].TrkHits[ii].size();
       myprt << std::right << std::setw(5) << nht;
       myprt << std::setw(4) << trk[itr].TrjPos.size();
       myprt << std::fixed;
       myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[0](0);
       myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[0](1);
       myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[0](2);
       unsigned short itj = trk[itr].TrjPos.size() - 1;
       myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[itj](0);
       myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[itj](1);
       myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[itj](2);
       myprt << std::setw(4) << trk[itr].VtxIndex[0] << std::setw(4) << trk[itr].VtxIndex[1];
       myprt << std::setw(4) << trk[itr].GoodEnd[0];
       myprt << std::setw(4) << trk[itr].GoodEnd[1];
       myprt << std::setw(4) << trk[itr].ChgOrder;
       myprt << std::right << std::setw(10) << std::setprecision(3) << trk[itr].TrjDir[itj](2);
       myprt << std::right << std::setw(4) << trk[itr].MomID;
       myprt << std::right << std::setw(5) << trk[itr].PDGCode << "   ";
       for (unsigned short ii = 0; ii < trk[itr].ClsEvtIndices.size(); ++ii)
         myprt << " " << trk[itr].ClsEvtIndices[ii];
       myprt << "\n";
     } // itr
 
   } // PrintTracks

void trkf::CCTrackMaker::produce ( art::Event & evt )

overrideprivatevirtual

Implements art::EDProducer.

Definition at line 377 of file CCTrackMaker_module.cc.

   {
     fWirePitch = wireReadoutGeom.Plane({0, 0, 0}).WirePitch();
 
     auto tcol = std::make_unique<std::vector<recob::Track>>();
     auto thassn = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
     auto vcol = std::make_unique<std::vector<recob::Vertex>>();
     auto pcol = std::make_unique<std::vector<recob::PFParticle>>();
     auto ptassn = std::make_unique<art::Assns<recob::PFParticle, recob::Track>>();
     auto pcassn = std::make_unique<art::Assns<recob::PFParticle, recob::Cluster>>();
     auto psassn = std::make_unique<art::Assns<recob::PFParticle, recob::Seed>>();
     auto pvassn = std::make_unique<art::Assns<recob::PFParticle, recob::Vertex>>();
 
     // seed collection
     auto scol = std::make_unique<std::vector<recob::Seed>>();
     auto shassn = std::make_unique<art::Assns<recob::Seed, recob::Hit>>();
 
     // all hits
     art::Handle<std::vector<recob::Hit>> allhitsListHandle;
     // cluster list
     art::Handle<std::vector<recob::Cluster>> clusterListHandle;
     std::vector<art::Ptr<recob::Cluster>> clusterlist;
     // ClusterCrawler Vertices
     art::Handle<std::vector<recob::Vertex>> VtxListHandle;
     std::vector<art::Ptr<recob::Vertex>> vtxlist;
 
     // get Hits
     allhits.clear();
     if (evt.getByLabel(fHitModuleLabel, allhitsListHandle))
       art::fill_ptr_vector(allhits, allhitsListHandle);
 
     // get Clusters
     if (evt.getByLabel(fClusterModuleLabel, clusterListHandle))
       art::fill_ptr_vector(clusterlist, clusterListHandle);
     if (clusterlist.size() == 0) return;
     // get cluster - hit associations
     art::FindManyP<recob::Hit> fmCluHits(clusterListHandle, evt, fClusterModuleLabel);
 
     // get Vertices
     if (evt.getByLabel(fVertexModuleLabel, VtxListHandle))
       art::fill_ptr_vector(vtxlist, VtxListHandle);
     art::FindManyP<recob::Cluster, unsigned short> fmVtxCls(VtxListHandle, evt, fVertexModuleLabel);
 
     std::vector<CluLen> clulens;
 
     unsigned short itr, tID, tIndex;
 
     // maximum error (see MakeClusterChains) for considering clusters broken
     fMaxMergeError = 30;
     // Cut on the error for a definitely broken cluster. Clusters with fMergeErrorCut < MergeError < fMaxMergeError
     // are possibly broken clusters but we will consider these when matching between planes
     fMergeErrorCut = 10;
 
     std::vector<art::Ptr<recob::Hit>> tmpHits;
     std::vector<art::Ptr<recob::Cluster>> tmpCls;
     std::vector<art::Ptr<recob::Vertex>> tmpVtx;
 
     // vector for PFParticle constructor
     std::vector<size_t> dtrIndices;
 
     // some vectors for recob::Seed
     double sPos[3], sDir[3];
     double sErr[3] = {0, 0, 0};
 
     auto const detProp =
       art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
 
     // check consistency between clusters and associated hits
     std::vector<art::Ptr<recob::Hit>> clusterhits;
     for (std::size_t icl = 0; icl < clusterlist.size(); ++icl) {
       auto const ipl = clusterlist[icl]->Plane().Plane;
       clusterhits = fmCluHits.at(icl);
       if (clusterhits[0]->WireID().Wire != std::nearbyint(clusterlist[icl]->EndWire())) {
         std::cout << "CCTM Cluster-Hit End wire mis-match " << clusterhits[0]->WireID().Wire
                   << " vs " << std::nearbyint(clusterlist[icl]->EndWire()) << " Bail out! \n";
         return;
       }
       for (unsigned short iht = 0; iht < clusterhits.size(); ++iht) {
         if (clusterhits[iht]->WireID().Plane != ipl) {
           std::cout << "CCTM Cluster-Hit plane mis-match " << ipl << " vs "
                     << clusterhits[iht]->WireID().Plane << " on hit " << iht << " Bail out! \n";
           return;
         } // hit-cluster plane mis-match
       }   // iht
     }     // icl
     // end check consistency
 
     vtx.clear();
     trk.clear();
     for (auto const& tpcid : geom->Iterate<geo::TPCID>()) {
       unsigned int const nplanes = wireReadoutGeom.Nplanes(tpcid);
       if (nplanes > 3) continue;
       for (std::size_t ipl = 0; ipl < 3; ++ipl) {
         cls[ipl].clear();
         clsChain[ipl].clear();
         trkHits[ipl].clear();
       } // ipl
       // FillWireHitRange also calculates the charge in each plane
       FillWireHitRange(tpcid);
       for (auto const& planeid : wireReadoutGeom.Iterate<geo::PlaneID>(tpcid)) {
         clulens.clear();
         // sort clusters by increasing End wire number
         for (std::size_t icl = 0; icl < clusterlist.size(); ++icl) {
           if (clusterlist[icl]->Plane() != planeid) continue;
           CluLen clulen;
           clulen.index = icl;
           clulen.length = clusterlist[icl]->EndWire();
           clulens.push_back(clulen);
         }
         if (clulens.empty()) continue;
         // sort clusters
         std::sort(clulens.begin(), clulens.end(), lessThan);
         if (clulens.empty()) continue;
         for (unsigned short ii = 0; ii < clulens.size(); ++ii) {
           const unsigned short icl = clulens[ii].index;
           clPar clstr;
           clstr.EvtIndex = icl;
           recob::Cluster const& cluster = *(clusterlist[icl]);
           // Begin info -> end index 1 (DS)
           clstr.Wire[1] = cluster.StartWire();
           clstr.Time[1] = cluster.StartTick();
           clstr.X[1] = (float)detProp.ConvertTicksToX(cluster.StartTick(), planeid);
           clstr.Angle[1] = cluster.StartAngle();
           clstr.Slope[1] = std::tan(cluster.StartAngle());
           clstr.Dir[1] = 0;
           clstr.Charge[1] = ChgNorm[planeid.Plane] * cluster.StartCharge();
           // this will be filled later
           clstr.ChgNear[1] = 0;
           clstr.VtxIndex[1] = -1;
           clstr.mVtxIndex[1] = -1;
           clstr.BrkIndex[1] = -1;
           clstr.MergeError[1] = fMaxMergeError;
           // End info -> end index 0 (US)
           clstr.Wire[0] = cluster.EndWire();
           clstr.Time[0] = cluster.EndTick();
           clstr.X[0] = (float)detProp.ConvertTicksToX(cluster.EndTick(), planeid);
           clstr.Angle[0] = cluster.EndAngle();
           clstr.Slope[0] = std::tan(cluster.EndAngle());
           clstr.Dir[0] = 0;
           if (clstr.Time[1] > clstr.Time[0]) {
             clstr.Dir[0] = 1;
             clstr.Dir[1] = -1;
           }
           else {
             clstr.Dir[0] = -1;
             clstr.Dir[1] = 1;
           }
           clstr.Charge[0] = ChgNorm[planeid.Plane] * cluster.EndCharge();
           // this will be filled later
           clstr.ChgNear[1] = 0;
           clstr.VtxIndex[0] = -1;
           clstr.mVtxIndex[0] = -1;
           clstr.BrkIndex[0] = -1;
           clstr.MergeError[0] = fMaxMergeError;
           // other info
           clstr.InTrack = -1;
           clstr.Length = (unsigned short)(0.5 + clstr.Wire[1] - clstr.Wire[0]);
           clstr.TotChg = ChgNorm[planeid.Plane] * cluster.Integral();
           if (clstr.TotChg <= 0) clstr.TotChg = 1;
           clusterhits = fmCluHits.at(icl);
           if (clusterhits.size() == 0) {
             mf::LogError("CCTM") << "No associated cluster hits " << icl;
             continue;
           }
           // correct charge for missing cluster hits
           clstr.TotChg *= clstr.Length / (float)clusterhits.size();
           cls[planeid.Plane].push_back(clstr);
         } // ii (icl)
       }   // planeid
 
       FillChgNear(detProp, tpcid);
 
       // and finally the vertices
       double xyz[3];
       for (unsigned short ivx = 0; ivx < vtxlist.size(); ++ivx) {
         vtxPar aVtx;
         aVtx.ID = ivx + 1;
         aVtx.EvtIndex = ivx;
         vtxlist[ivx]->XYZ(xyz);
         aVtx.X = xyz[0];
         aVtx.Y = xyz[1];
         aVtx.Z = xyz[2];
         aVtx.nClusInPln[0] = 0;
         aVtx.nClusInPln[1] = 0;
         aVtx.nClusInPln[2] = 0;
         std::vector<art::Ptr<recob::Cluster>> const& vtxCls = fmVtxCls.at(ivx);
         std::vector<const unsigned short*> const& vtxClsEnd = fmVtxCls.data(ivx);
         for (unsigned short vcass = 0; vcass < vtxCls.size(); ++vcass) {
           auto const icl = vtxCls[vcass].key();
           // the cluster plane
           auto const ipl = vtxCls[vcass]->Plane().Plane;
           auto end = *vtxClsEnd[vcass];
           if (end > 1)
             throw cet::exception("CCTM")
               << "Invalid end data from vertex - cluster association" << end;
           bool gotit = false;
           // CCTM end is opposite of CC end
           end = 1 - end;
           for (unsigned short jcl = 0; jcl < cls[ipl].size(); ++jcl) {
             if (cls[ipl][jcl].EvtIndex == icl) {
               cls[ipl][jcl].VtxIndex[end] = ivx;
               ++aVtx.nClusInPln[ipl];
               gotit = true;
               break;
             } // index check
           }   // jcl
           if (!gotit)
             throw cet::exception("CCTM") << "Bad index from vertex - cluster association" << icl;
         } // icl
         vtx.push_back(aVtx);
       } // ivx
       // Find broken clusters
       MakeClusterChains(detProp, fmCluHits, tpcid);
       FindMaybeVertices(tpcid);
 
       // call algorithms in the specified order
       matcomb.clear();
       for (algIndex = 0; algIndex < fMatchAlgs.size(); ++algIndex) {
         if (fMatchAlgs[algIndex] == 1) {
           prt = (fDebugAlg == 1);
           VtxMatch(detProp, fmCluHits, tpcid);
           if (fMakeAlgTracks[algIndex]) SortMatches(detProp, fmCluHits, 1, tpcid);
         }
         else if (fMatchAlgs[algIndex] == 2) {
           prt = (fDebugAlg == 2);
           PlnMatch(detProp, tpcid);
           if (fMakeAlgTracks[algIndex]) SortMatches(detProp, fmCluHits, 2, tpcid);
         }
         if (prt) PrintClusters(detProp, tpcid);
       } // algIndex
       prt = false;
       pfpToTrkID.clear();
       // Determine the vertex/track hierarchy
       if (fMakePFPs) {
         TagCosmics(tpcid);
         MakeFamily(tpcid);
       }
       FitVertices(detProp, tpcid);
 
       // Make PFParticles -> tracks
       for (unsigned short ipf = 0; ipf < pfpToTrkID.size(); ++ipf) {
         // trackID of this PFP
         tID = pfpToTrkID[ipf];
         if (tID > trk.size() + 1) {
           mf::LogWarning("CCTM") << "Bad track ID";
           continue;
         }
         dtrIndices.clear();
         // load the daughter PFP indices
         mf::LogVerbatim("CCTM") << "PFParticle " << ipf << " tID " << tID;
         for (unsigned short jpf = 0; jpf < pfpToTrkID.size(); ++jpf) {
           itr = pfpToTrkID[jpf] - 1; // convert from track ID to track index
           if (trk[itr].MomID == tID) dtrIndices.push_back(jpf);
           if (trk[itr].MomID == tID)
             mf::LogVerbatim("CCTM") << " dtr jpf " << jpf << " itr " << itr;
         } // jpf
         unsigned short parentIndex = USHRT_MAX;
         if (tID == 0) {
           // make neutrino PFP USHRT_MAX == primary PFP
           recob::PFParticle pfp(14, ipf, parentIndex, dtrIndices);
           pcol->emplace_back(std::move(pfp));
           for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
             if (!vtx[ivx].Neutrino) continue;
             // make the vertex
             xyz[0] = vtx[ivx].X;
             xyz[1] = vtx[ivx].Y;
             xyz[2] = vtx[ivx].Z;
             size_t vStart = vcol->size();
             recob::Vertex vertex(xyz, vtx[ivx].ID);
             vcol->emplace_back(std::move(vertex));
             size_t vEnd = vcol->size();
             // PFParticle - vertex association
             util::CreateAssn(evt, *pcol, *vcol, *pvassn, vStart, vEnd);
             vtx[ivx].ID = -vtx[ivx].ID;
             break;
           } // ivx
         }
         else if (tID > 0) {
           // make non-neutrino PFP. Need to find the parent PFP index
           // Track index of this PFP
           tIndex = tID - 1;
           // trackID of this PFP parent
           for (unsigned short ii = 0; ii < pfpToTrkID.size(); ++ii) {
             if (pfpToTrkID[ii] == trk[tIndex].MomID) {
               parentIndex = ii;
               break;
             }
           } // ii
           // PFParticle
           mf::LogVerbatim("CCTM") << "daughters tID " << tID << " pdg " << trk[tIndex].PDGCode
                                   << " ipf " << ipf << " parentIndex " << parentIndex
                                   << " dtr size " << dtrIndices.size();
           recob::PFParticle pfp(trk[tIndex].PDGCode, ipf, parentIndex, dtrIndices);
           pcol->emplace_back(std::move(pfp));
           // track
           // make the track
           size_t tStart = tcol->size();
           recob::Track track(
             recob::TrackTrajectory(recob::tracking::convertCollToPoint(trk[tIndex].TrjPos),
                                    recob::tracking::convertCollToVector(trk[tIndex].TrjDir),
                                    recob::Track::Flags_t(trk[itr].TrjPos.size()),
                                    false),
             0,
             -1.,
             0,
             recob::tracking::SMatrixSym55(),
             recob::tracking::SMatrixSym55(),
             tID);
           tcol->emplace_back(std::move(track));
           size_t tEnd = tcol->size();
           // PFParticle - track association
           util::CreateAssn(evt, *pcol, *tcol, *ptassn, tStart, tEnd);
           // flag this track as already put in the event
           trk[tIndex].ID = -trk[tIndex].ID;
           // track -> hits association
           tmpHits.clear();
           for (unsigned int ipl = 0; ipl < nplanes; ++ipl)
             tmpHits.insert(
               tmpHits.end(), trk[tIndex].TrkHits[ipl].begin(), trk[tIndex].TrkHits[ipl].end());
           util::CreateAssn(evt, *tcol, tmpHits, *thassn);
           // Find seed hits and the end of the track that is best
           unsigned int end = 0;
           unsigned short itj = 0;
           if (end > 0) itj = trk[tIndex].TrjPos.size() - 1;
           for (unsigned short ii = 0; ii < 3; ++ii) {
             sPos[ii] = trk[tIndex].TrjPos[itj](ii);
             sDir[ii] = trk[tIndex].TrjDir[itj](ii);
           } // ii
           size_t sStart = scol->size();
           recob::Seed seed(sPos, sDir, sErr, sErr);
           scol->emplace_back(std::move(seed));
           size_t sEnd = scol->size();
           // PFP-seed association
           util::CreateAssn(evt, *pcol, *scol, *psassn, sStart, sEnd);
           // Seed-hit association
           tmpHits.clear();
           for (unsigned int ipl = 0; ipl < nplanes; ++ipl)
             tmpHits.insert(tmpHits.end(), seedHits[ipl].begin(), seedHits[ipl].end());
           util::CreateAssn(evt, *scol, tmpHits, *shassn);
           // cluster association
           // PFP-cluster association
           tmpCls.clear();
           for (unsigned short ii = 0; ii < trk[tIndex].ClsEvtIndices.size(); ++ii) {
             auto const icl = trk[tIndex].ClsEvtIndices[ii];
             tmpCls.push_back(clusterlist[icl]);
           } // ii
           util::CreateAssn(evt, *pcol, tmpCls, *pcassn);
         } // non-neutrino PFP
       }   // ipf
       // make non-PFP tracks
       for (unsigned short itr = 0; itr < trk.size(); ++itr) {
         // ignore already saved tracks
         if (trk[itr].ID < 0) continue;
         recob::Track track(
           recob::TrackTrajectory(recob::tracking::convertCollToPoint(trk[itr].TrjPos),
                                  recob::tracking::convertCollToVector(trk[itr].TrjDir),
                                  recob::Track::Flags_t(trk[itr].TrjPos.size()),
                                  false),
           0,
           -1.,
           0,
           recob::tracking::SMatrixSym55(),
           recob::tracking::SMatrixSym55(),
           trk[itr].ID);
         tcol->emplace_back(std::move(track));
         tmpHits.clear();
         for (unsigned int ipl = 0; ipl < nplanes; ++ipl)
           tmpHits.insert(tmpHits.end(), trk[itr].TrkHits[ipl].begin(), trk[itr].TrkHits[ipl].end());
         util::CreateAssn(evt, *tcol, tmpHits, *thassn);
       } // itr
       // make remnant vertices
       for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
         if (vtx[ivx].ID < 0) continue;
         xyz[0] = vtx[ivx].X;
         xyz[1] = vtx[ivx].Y;
         xyz[2] = vtx[ivx].Z;
         recob::Vertex vertex(xyz, vtx[ivx].ID);
         vcol->emplace_back(std::move(vertex));
       }
       if (fDebugAlg > 0) PrintTracks();
 
       double orphanLen = 0;
       for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
         for (std::size_t icl = 0; icl < cls[ipl].size(); ++icl) {
           if (cls[ipl][icl].Length > 40 && cls[ipl][icl].InTrack < 0) {
             orphanLen += cls[ipl][icl].Length;
             // unused cluster
             mf::LogVerbatim("CCTM")
               << "Orphan long cluster " << ipl << ":" << icl << ":" << cls[ipl][icl].Wire[0] << ":"
               << (int)cls[ipl][icl].Time[0] << " length " << cls[ipl][icl].Length;
           }
         } // icl
 
         cls[ipl].clear();
         clsChain[ipl].clear();
         std::cout << "Total orphan length " << orphanLen << "\n";
         trkHits[ipl].clear();
         seedHits[ipl].clear();
         vxCls[ipl].clear();
       } // ipl
     }   // tpcs
 
     evt.put(std::move(pcol));
     evt.put(std::move(ptassn));
     evt.put(std::move(pcassn));
     evt.put(std::move(pvassn));
     evt.put(std::move(psassn));
     evt.put(std::move(tcol));
     evt.put(std::move(thassn));
     evt.put(std::move(scol));
     evt.put(std::move(vcol));
 
     // final cleanup
     vtx.clear();
     trk.clear();
     allhits.clear();
     matcomb.clear();
     pfpToTrkID.clear();
 
   } // produce

void art::Modifier::registerProducts ( ProductDescriptions & productsToRegister )

inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

   {
     ProductRegistryHelper::registerProducts(productsToRegister,
                                             moduleDescription());
   }

void art::ModuleBase::setModuleDescription ( ModuleDescription const & md )

inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

   {
     md_ = md;
   }

void art::ModuleBase::sortConsumables ( std::string const & current_process_name )

inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

   {
     // Now that we know we have seen all the consumes declarations,
     // sort the results for fast lookup later.
     collector_.sortConsumables(current_process_name);
   }

void trkf::CCTrackMaker::SortMatches	(	detinfo::DetectorPropertiesData const &	detProp,
		art::FindManyP< recob::Hit > const &	fmCluHits,
		unsigned short	procCode,
		geo::TPCID const &	tpcid
	)

private

Definition at line 2311 of file CCTrackMaker_module.cc.

References util::begin(), clsChain, util::end(), trkf::fill(), FillTrkHits(), lessThan(), matcomb, geo::WireReadoutGeom::Nplanes(), prt, art::right(), util::size(), StoreTrack(), and wireReadoutGeom.

Referenced by produce(), and VtxMatch().

   {
     // sort cluster matches by increasing total match error. Find the minimum total error of all
     // cluster match combinations and make tracks from them
     CluLen merr;
     std::vector<CluLen> materr;
 
     if (matcomb.size() == 0) return;
 
     // sort by decreasing error
     for (std::size_t ii = 0; ii < matcomb.size(); ++ii) {
       merr.index = ii;
       merr.length = matcomb[ii].Err + matcomb[ii].oErr;
       materr.push_back(merr);
     } // ii
     std::sort(materr.begin(), materr.end(), lessThan);
 
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     if (prt) {
       mf::LogVerbatim myprt("CCTM");
       myprt << "SortMatches\n";
       myprt << "   ii    im  Vx   Err     dW     dA     dX  oVx   oErr    odW   odA    odX   Asym  "
                " icl   jcl   kcl \n";
       for (std::size_t ii = 0; ii < materr.size(); ++ii) {
         unsigned int im = materr[ii].index;
         float asym =
           fabs(matcomb[im].Chg[0] - matcomb[im].Chg[1]) / (matcomb[im].Chg[0] + matcomb[im].Chg[1]);
         asym *=
           fabs(matcomb[im].Chg[1] - matcomb[im].Chg[2]) / (matcomb[im].Chg[1] + matcomb[im].Chg[2]);
         myprt << std::fixed << std::right << std::setw(5) << ii << std::setw(5) << im
               << std::setw(4) << matcomb[im].Vtx << std::setw(7) << std::setprecision(2)
               << matcomb[im].Err << std::setw(7) << std::setprecision(1) << matcomb[im].dWir
               << std::setw(7) << std::setprecision(2) << matcomb[im].dAng << std::setw(7)
               << std::setprecision(2) << matcomb[im].dX << std::setw(4) << matcomb[im].oVtx
               << std::setw(7) << std::setprecision(2) << matcomb[im].oErr << std::setw(7)
               << std::setprecision(1) << matcomb[im].odWir << std::setw(7) << std::setprecision(2)
               << matcomb[im].odAng << std::setw(7) << std::setprecision(2) << matcomb[im].odX
               << std::setw(7) << std::setprecision(3) << asym << " 0:" << matcomb[im].Cls[0] << ":"
               << matcomb[im].End[0] << " 1:" << matcomb[im].Cls[1] << ":" << matcomb[im].End[1];
         if (nplanes > 2) myprt << " 2:" << matcomb[im].Cls[2] << ":" << matcomb[im].End[2];
         myprt << "\n";
       } // ii
     }   // prt
 
     // define an array to ensure clusters are only used once
     std::array<std::vector<bool>, 3> pclUsed;
     for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
       pclUsed[ipl].resize(clsChain[ipl].size());
     }
 
     // count the total number of clusters and length used in matcomb
     unsigned short matcombTotCl = 0;
     float matcombTotLen = 0;
     unsigned short icl;
     for (std::size_t ii = 0; ii < matcomb.size(); ++ii) {
       for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
         if (matcomb[ii].Cls[ipl] < 0) continue;
         icl = matcomb[ii].Cls[ipl];
         ++matcombTotCl;
         matcombTotLen += clsChain[ipl][icl].Length;
       }
     }
 
     if (prt)
       mf::LogVerbatim("CCTM") << "Number of clusters to match " << matcombTotCl << " total length "
                               << matcombTotLen;
 
     if (matcombTotLen <= 0) {
       mf::LogError("CCTM") << "SortMatches: bad matcomb total length " << matcombTotLen;
       return;
     }
 
     // vector of matcomb indices of unique cluster matches
     std::vector<unsigned short> matIndex;
     // vector of matcomb indices of unique cluster matches that have the best total error
     std::vector<unsigned short> bestMatIndex;
     float totLen, totErr, bestTotErr = 9999;
     // start with the best match
     unsigned short jj, jm, nused, jcl;
     // fraction of the length of all clustters in matcomb that are used in a match
     float fracLen;
 
     for (std::size_t ii = 0; ii < materr.size(); ++ii) {
       unsigned int im = materr[ii].index;
       matIndex.clear();
       // skip really bad matches
       if (matcomb[im].Err > bestTotErr) continue;
       totLen = 0;
       // initialize pclUsed and flag the clusters in this match
       for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
         // initialize to no clusters used
         std::fill(pclUsed[ipl].begin(), pclUsed[ipl].end(), false);
         // check for 2 plane match
         if (matcomb[im].Cls[ipl] < 0) continue;
         icl = matcomb[im].Cls[ipl];
         pclUsed[ipl][icl] = true;
         totLen += clsChain[ipl][icl].Length;
       } // ipl
       // Initialize the error sum
       totErr = matcomb[im].Err;
       // Save the index
       matIndex.push_back(im);
       // look for matches in the rest of the list that are not already matched.
       for (jj = 0; jj < materr.size(); ++jj) {
         if (jj == ii) continue;
         jm = materr[jj].index;
         // skip really bad matches
         if (matcomb[jm].Err > bestTotErr) continue;
         // check for non-unique cluster indices
         nused = 0;
         for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
           if (matcomb[jm].Cls[ipl] < 0) continue;
           jcl = matcomb[jm].Cls[ipl];
           if (pclUsed[ipl][jcl]) ++nused;
           // This cluster chain was used in a previous match
           if (nused > 0) break;
           totLen += clsChain[ipl][jcl].Length;
         } // ipl
         // at least one of the clusters in this match have been used
         if (nused != 0) continue;
         // found a match with an unmatched set of clusters. Update the total error and flag them
         totErr += matcomb[jm].Err;
         matIndex.push_back(jm);
         // Flag the clusters used and see if all of them are used
         for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
           if (matcomb[jm].Cls[ipl] < 0) continue;
           jcl = matcomb[jm].Cls[ipl];
           pclUsed[ipl][jcl] = true;
         } // ipl
       }   // jm
       if (totLen == 0) continue;
       nused = 0;
       for (unsigned int ipl = 0; ipl < nplanes; ++ipl) {
         for (unsigned short indx = 0; indx < pclUsed[ipl].size(); ++indx)
           if (pclUsed[ipl][indx]) ++nused;
       } // ipl
       if (totLen > matcombTotLen) std::cout << "Oops " << totLen << " " << matcombTotLen << "\n";
       // weight the total error by the total length of all clusters
       fracLen = totLen / matcombTotLen;
       totErr /= fracLen;
       if (prt) {
         mf::LogVerbatim myprt("CCTM");
         myprt << "match " << im << " totErr " << totErr << " nused " << nused << " fracLen "
               << fracLen << " totLen " << totLen << " mat: ";
         for (unsigned short indx = 0; indx < matIndex.size(); ++indx)
           myprt << " " << matIndex[indx];
       } // prt
       // check for more used clusters and a better total error
       if (totErr < bestTotErr) {
         bestTotErr = totErr;
         bestMatIndex = matIndex;
         if (nused == matcombTotCl) break;
         if (prt) {
           mf::LogVerbatim myprt("CCTM");
           myprt << "bestTotErr " << bestTotErr << " nused " << nused << " matcombTotCl "
                 << matcombTotCl << " mat: ";
           for (unsigned short indx = 0; indx < bestMatIndex.size(); ++indx)
             myprt << " " << bestMatIndex[indx];
         } // prt
         // stop looking if we have found everything
         if (fracLen > 0.999) break;
       } // totErr < bestTotErr
     }   // im
 
     if (bestTotErr > 9000) return;
 
     for (std::size_t ii = 0; ii < bestMatIndex.size(); ++ii) {
       unsigned int im = bestMatIndex[ii];
       FillTrkHits(fmCluHits, im, nplanes);
       // look for missing clusters?
       // store this track with processor code 1
       StoreTrack(detProp, fmCluHits, im, procCode, tpcid);
     } // ii
 
   } // SortMatches

void trkf::CCTrackMaker::StoreTrack	(	detinfo::DetectorPropertiesData const &	detProp,
		art::FindManyP< recob::Hit > const &	fmCluHits,
		unsigned short	imat,
		unsigned short	procCode,
		geo::TPCID const &	tpcid
	)

private

Definition at line 1864 of file CCTrackMaker_module.cc.

Referenced by SortMatches().

   {
     // store the current "under construction" track in the trk vector
 
     TrkPar newtrk;
 
     if (imat > matcomb.size() - 1) {
       mf::LogError("CCTM") << "Bad imat in StoreTrack";
       return;
     }
 
     // ensure there are at least 2 hits in at least 2 planes
     unsigned short nhitinpl = 0;
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (unsigned short ipl = 0; ipl < nplanes; ++ipl)
       if (trkHits[ipl].size() > 1) ++nhitinpl;
     if (nhitinpl < 2) {
       mf::LogError("CCTM") << "StoreTrack: Not enough hits in each plane\n";
       return;
     }
     if (prt)
       mf::LogVerbatim("CCTM") << "In StoreTrack: matcomb " << imat << " cluster chains "
                               << matcomb[imat].Cls[0] << " " << matcomb[imat].Cls[1] << " "
                               << matcomb[imat].Cls[2];
 
     // Track hit vectors for fitting the trajectory
     std::array<std::vector<geo::WireID>, 3> trkWID;
     std::array<std::vector<double>, 3> trkX;
     std::array<std::vector<double>, 3> trkXErr;
 
     // track trajectory for a track
     std::vector<TVector3> trkPos;
     std::vector<TVector3> trkDir;
 
     newtrk.ID = trk.size() + 1;
     newtrk.Proc = procCode;
     newtrk.TrkHits = trkHits;
     // BUG the double brace syntax is required to work around clang bug 21629
     // (https://bugs.llvm.org/show_bug.cgi?id=21629)
     newtrk.VtxIndex = {{-1, -1}};
     newtrk.ChgOrder = 0;
     newtrk.MomID = -1;
     // BUG the double brace syntax is required to work around clang bug 21629
     // (https://bugs.llvm.org/show_bug.cgi?id=21629)
     newtrk.EndInTPC = {{false, false}};
     newtrk.GoodEnd = {{false, false}};
     newtrk.DtrID = {0};
     newtrk.PDGCode = -1;
 
     unsigned short icl, iht;
 
     if (prt)
       mf::LogVerbatim("CCTM") << "CCTM: Make traj for track " << newtrk.ID << " procCode "
                               << procCode << " nhits in planes " << trkHits[0].size() << " "
                               << trkHits[1].size() << " " << trkHits[2].size();
     // make the track trajectory
     if (nplanes == 2) {
       trkWID[2].resize(0);
       trkX[2].resize(0);
       trkXErr[2].resize(0);
     }
     for (auto const& planeid : wireReadoutGeom.Iterate<geo::PlaneID>(tpcid)) {
       auto const ipl = planeid.Plane;
       trkWID[ipl].resize(trkHits[ipl].size());
       trkX[ipl].resize(trkHits[ipl].size());
       trkXErr[ipl].resize(trkHits[ipl].size());
       for (iht = 0; iht < trkHits[ipl].size(); ++iht) {
         trkWID[ipl][iht] = trkHits[ipl][iht]->WireID();
         trkX[ipl][iht] = detProp.ConvertTicksToX(trkHits[ipl][iht]->PeakTime(), planeid);
         trkXErr[ipl][iht] =
           fHitFitErrFac * trkHits[ipl][iht]->RMS() * trkHits[ipl][iht]->Multiplicity();
       } // iht
     }   // ipl
     fTrackTrajectoryAlg.TrackTrajectory(trkWID, trkX, trkXErr, trkPos, trkDir);
     if (trkPos.size() < 2) {
       mf::LogError("CCTM") << "StoreTrack: No trajectory points on failed track " << newtrk.ID
                            << " in StoreTrack: matcomb " << imat << " cluster chains "
                            << matcomb[imat].Cls[0] << " " << matcomb[imat].Cls[1] << " "
                            << matcomb[imat].Cls[2];
       // make a garbage trajectory
       trkPos.resize(2);
       trkPos[1](2) = 1;
       trkDir.resize(2);
       trkDir[1](2) = 1;
     }
     newtrk.TrjPos = trkPos;
     newtrk.TrjDir = trkDir;
 
     if (prt) mf::LogVerbatim("CCTM") << " number of traj points " << trkPos.size();
 
     // determine if each end is good in the sense that there are hits in each plane
     // that are consistent in time and are presumed to form a good 3D space point
     unsigned short nClose, indx, jndx;
     float xErr;
     for (unsigned int end = 0; end < 2; ++end) {
       nClose = 0;
       for (unsigned short ipl = 0; ipl < nplanes - 1; ++ipl) {
         if (trkX[ipl].size() == 0) continue;
         for (unsigned short jpl = ipl + 1; jpl < nplanes; ++jpl) {
           if (trkX[jpl].size() == 0) continue;
           if (end == 0) {
             indx = 0;
             jndx = 0;
           }
           else {
             indx = trkXErr[ipl].size() - 1;
             jndx = trkXErr[jpl].size() - 1;
           }
           xErr = 3 * (trkXErr[ipl][indx] + trkXErr[jpl][jndx]);
           if (std::abs(trkX[ipl][indx] - trkX[jpl][jndx]) < xErr) ++nClose;
         } // jpl
       }   // ipl
       if (nClose == nplanes) newtrk.GoodEnd[end] = true;
     } // end
 
     // set trajectory end points to a vertex if one exists
     unsigned short ivx, itj, ccl;
     float dx, dy, dz, dr0, dr1;
     unsigned short attachEnd;
     for (unsigned int end = 0; end < 2; ++end) {
       ivx = USHRT_MAX;
       if (end == 0 && matcomb[imat].Vtx >= 0) ivx = matcomb[imat].Vtx;
       if (end == 1 && matcomb[imat].oVtx >= 0) ivx = matcomb[imat].oVtx;
       if (ivx == USHRT_MAX) continue;
       // determine the proper end using the TrjPos order and brute force
       itj = 0;
       dx = vtx[ivx].X - newtrk.TrjPos[itj](0);
       dy = vtx[ivx].Y - newtrk.TrjPos[itj](1);
       dz = vtx[ivx].Z - newtrk.TrjPos[itj](2);
       dr0 = dx * dx + dy * dy + dz * dz;
       itj = newtrk.TrjPos.size() - 1;
       dx = vtx[ivx].X - newtrk.TrjPos[itj](0);
       dy = vtx[ivx].Y - newtrk.TrjPos[itj](1);
       dz = vtx[ivx].Z - newtrk.TrjPos[itj](2);
       dr1 = dx * dx + dy * dy + dz * dz;
       attachEnd = 1;
       if (dr0 < dr1) {
         itj = 0;
         attachEnd = 0;
         // a really bad match to the vertex
         if (dr0 > 5) return;
       }
       else {
         // a really bad match to the vertex
         if (dr1 > 5) return;
       }
       newtrk.TrjPos[itj](0) = vtx[ivx].X;
       newtrk.TrjPos[itj](1) = vtx[ivx].Y;
       newtrk.TrjPos[itj](2) = vtx[ivx].Z;
       newtrk.VtxIndex[attachEnd] = ivx;
       // correct the trajectory direction
       TVector3 dir;
       if (itj == 0) {
         dir = newtrk.TrjPos[1] - newtrk.TrjPos[0];
         newtrk.TrjDir[0] = dir.Unit();
       }
       else {
         dir = newtrk.TrjPos[itj - 1] - newtrk.TrjPos[itj];
         newtrk.TrjDir[itj] = dir.Unit();
       }
     } // end
 
     if (newtrk.VtxIndex[0] >= 0 && newtrk.VtxIndex[0] == newtrk.VtxIndex[1]) {
       mf::LogError("CCTM")
         << "StoreTrack: Trying to attach a vertex to both ends of a track. imat = " << imat;
       return;
     }
 
     // calculate the length
     newtrk.Length = 0;
     float norm;
     double X, Y, Z;
     for (unsigned short itj = 1; itj < newtrk.TrjPos.size(); ++itj) {
       X = newtrk.TrjPos[itj](0) - newtrk.TrjPos[itj - 1](0);
       Y = newtrk.TrjPos[itj](1) - newtrk.TrjPos[itj - 1](1);
       Z = newtrk.TrjPos[itj](2) - newtrk.TrjPos[itj - 1](2);
       norm = sqrt(X * X + Y * Y + Z * Z);
       newtrk.Length += norm;
     }
 
     // store the cluster -> track assignment
     newtrk.ClsEvtIndices.clear();
     for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
       if (matcomb[imat].Cls[ipl] < 0) continue;
       ccl = matcomb[imat].Cls[ipl];
       if (ccl > clsChain[ipl].size()) std::cout << "oops StoreTrack\n";
       clsChain[ipl][ccl].InTrack = newtrk.ID;
       for (unsigned short icc = 0; icc < clsChain[ipl][ccl].ClsIndex.size(); ++icc) {
         icl = clsChain[ipl][ccl].ClsIndex[icc];
         if (icl > cls[ipl].size()) std::cout << "oops StoreTrack\n";
         cls[ipl][icl].InTrack = newtrk.ID;
         if (cls[ipl][icl].EvtIndex > fmCluHits.size() - 1) {
           std::cout << "ooops2 store track EvtIndex " << cls[ipl][icl].EvtIndex << " size "
                     << fmCluHits.size() << " icl " << icl << "\n";
           continue;
         }
         newtrk.ClsEvtIndices.push_back(cls[ipl][icl].EvtIndex);
       } // icc
     }   // ipl
 
     if (prt) mf::LogVerbatim("CCTM") << " track ID " << newtrk.ID << " stored in StoreTrack";
 
     trk.push_back(newtrk);
   } // StoreTrack

void trkf::CCTrackMaker::TagCosmics ( geo::TPCID const & tpcid )

private

Definition at line 1146 of file CCTrackMaker_module.cc.

References fDeltaRayCut, fFiducialCut, geom, geo::TPCGeo::GetCenter(), geo::TPCGeo::HalfHeight(), geo::TPCGeo::HalfWidth(), geo::TPCGeo::Length(), tca::Length(), pfpToTrkID, and geo::GeometryCore::TPC().

Referenced by produce().

   {
     // Make cosmic ray PFParticles
     unsigned short ipf, itj;
     bool skipit = true;
 
     // Y,Z limits of the detector
     const geo::TPCGeo& thetpc = geom->TPC(tpcid);
     auto const world = thetpc.GetCenter();
 
     float XLo = world.X() - thetpc.HalfWidth() + fFiducialCut;
     float XHi = world.X() + thetpc.HalfWidth() - fFiducialCut;
     float YLo = world.Y() - thetpc.HalfHeight() + fFiducialCut;
     float YHi = world.Y() + thetpc.HalfHeight() - fFiducialCut;
     float ZLo = world.Z() - thetpc.Length() / 2 + fFiducialCut;
     float ZHi = world.Z() + thetpc.Length() / 2 - fFiducialCut;
 
     bool startsIn, endsIn;
 
     for (unsigned short itk = 0; itk < trk.size(); ++itk) {
       // ignore already used tracks
       if (trk[itk].ID < 0) continue;
       // ignore short tracks (< 10 cm)
       if (trk[itk].Length < 10) continue;
       // check for already identified PFPs
       skipit = false;
       for (ipf = 0; ipf < pfpToTrkID.size(); ++ipf) {
         if (pfpToTrkID[ipf] == trk[itk].ID) {
           skipit = true;
           break;
         }
       } // ipf
       if (skipit) continue;
       startsIn = true;
       if (trk[itk].TrjPos[0](0) < XLo || trk[itk].TrjPos[0](0) > XHi) startsIn = false;
       if (trk[itk].TrjPos[0](1) < YLo || trk[itk].TrjPos[0](1) > YHi) startsIn = false;
       if (trk[itk].TrjPos[0](2) < ZLo || trk[itk].TrjPos[0](2) > ZHi) startsIn = false;
       if (startsIn) continue;
       endsIn = true;
       itj = trk[itk].TrjPos.size() - 1;
       if (trk[itk].TrjPos[itj](0) < XLo || trk[itk].TrjPos[itj](0) > XHi) endsIn = false;
       if (trk[itk].TrjPos[itj](1) < YLo || trk[itk].TrjPos[itj](1) > YHi) endsIn = false;
       if (trk[itk].TrjPos[itj](2) < ZLo || trk[itk].TrjPos[itj](2) > ZHi) endsIn = false;
       if (endsIn) continue;
       // call it a cosmic muon
       trk[itk].PDGCode = 13;
       pfpToTrkID.push_back(trk[itk].ID);
     } // itk
 
     if (fDeltaRayCut <= 0) return;
 
     for (unsigned short itk = 0; itk < trk.size(); ++itk) {
       // find a tagged cosmic ray
       if (trk[itk].PDGCode != 13) continue;
 
     } // itk
 
   } // TagCosmics

void trkf::CCTrackMaker::VtxMatch	(	detinfo::DetectorPropertiesData const &	detProp,
		art::FindManyP< recob::Hit > const &	fmCluHits,
		geo::TPCID const &	tpcid
	)

private

Definition at line 1206 of file CCTrackMaker_module.cc.

References trkf::CCTrackMaker::MatchPars::Chg, clear(), trkf::CCTrackMaker::MatchPars::Cls, clsChain, DupMatch(), trkf::CCTrackMaker::MatchPars::End, trkf::CCTrackMaker::MatchPars::Err, FillEndMatch(), FillEndMatch2(), matcomb, geo::WireReadoutGeom::Nplanes(), trkf::CCTrackMaker::MatchPars::oErr, trkf::CCTrackMaker::MatchPars::oVtx, prt, SortMatches(), vtx, trkf::CCTrackMaker::MatchPars::Vtx, vxCls, wireReadoutGeom, and X.

Referenced by produce().

   {
     // Use vertex assignments to match clusters
     unsigned short ivx, ii, ipl, icl, jj, jpl, jcl, kk, kpl, kcl;
     short idir, iend, jdir, jend, kdir, kend, ioend;
 
     auto const nplanes = wireReadoutGeom.Nplanes(tpcid);
     for (ivx = 0; ivx < vtx.size(); ++ivx) {
 
       // list of cluster chains associated with this vertex in each plane
       for (ipl = 0; ipl < nplanes; ++ipl) {
         vxCls[ipl].clear();
         for (icl = 0; icl < clsChain[ipl].size(); ++icl) {
           if (clsChain[ipl][icl].InTrack >= 0) continue;
           for (iend = 0; iend < 2; ++iend) {
             if (clsChain[ipl][icl].VtxIndex[iend] == vtx[ivx].EvtIndex) vxCls[ipl].push_back(icl);
           } // end
         }   // icl
       }     // ipl
 
       if (prt) {
         mf::LogVerbatim myprt("CCTM");
         myprt << "VtxMatch: Vertex ID " << vtx[ivx].EvtIndex << "\n";
         for (ipl = 0; ipl < nplanes; ++ipl) {
           myprt << "ipl " << ipl << " cls";
           for (unsigned short ii = 0; ii < vxCls[ipl].size(); ++ii)
             myprt << " " << vxCls[ipl][ii];
           myprt << "\n";
         } // ipl
       }   // prt
       // match between planes, requiring clusters matched to this vertex
       iend = 0;
       jend = 0;
       bool gotkcl;
       float totErr;
       for (ipl = 0; ipl < nplanes; ++ipl) {
         if (nplanes == 2 && ipl > 0) continue;
         for (ii = 0; ii < vxCls[ipl].size(); ++ii) {
           icl = vxCls[ipl][ii];
           // ignore used clusters
           if (clsChain[ipl][icl].InTrack >= 0) continue;
           jpl = (ipl + 1) % nplanes;
           kpl = (jpl + 1) % nplanes;
           for (jj = 0; jj < vxCls[jpl].size(); ++jj) {
             jcl = vxCls[jpl][jj];
             if (clsChain[jpl][jcl].InTrack >= 0) continue;
             for (iend = 0; iend < 2; ++iend) {
               if (clsChain[ipl][icl].VtxIndex[iend] != vtx[ivx].EvtIndex) continue;
               ioend = 1 - iend;
               idir = clsChain[ipl][icl].Dir[iend];
               for (jend = 0; jend < 2; ++jend) {
                 if (clsChain[jpl][jcl].VtxIndex[jend] != vtx[ivx].EvtIndex) continue;
                 jdir = clsChain[jpl][jcl].Dir[jend];
                 if (idir != 0 && jdir != 0 && idir != jdir) continue;
                 // ignore outrageously bad other end X matches
                 if (fabs(clsChain[jpl][jcl].X[1 - jend] - clsChain[ipl][icl].X[ioend]) > 50)
                   continue;
                 MatchPars match;
                 match.Cls[ipl] = icl;
                 match.End[ipl] = iend;
                 match.Cls[jpl] = jcl;
                 match.End[jpl] = jend;
                 match.Vtx = ivx;
                 match.oVtx = -1;
                 // set large so that DupMatch doesn't get confused when called before FillEndMatch
                 match.Err = 1E6;
                 match.oErr = 1E6;
                 if (nplanes == 2) {
                   FillEndMatch2(match);
                   if (prt)
                     mf::LogVerbatim("CCTM")
                       << "FillEndMatch2: Err " << match.Err << " oErr " << match.oErr;
                   if (match.Err + match.oErr > 100) continue;
                   if (DupMatch(match, nplanes)) continue;
                   matcomb.push_back(match);
                   continue;
                 }
                 match.Cls[kpl] = -1;
                 match.End[kpl] = 0;
                 if (prt)
                   mf::LogVerbatim("CCTM")
                     << "VtxMatch: check " << ipl << ":" << icl << ":" << iend << " and " << jpl
                     << ":" << jcl << ":" << jend << " for cluster in kpl " << kpl;
                 gotkcl = false;
                 for (kk = 0; kk < vxCls[kpl].size(); ++kk) {
                   kcl = vxCls[kpl][kk];
                   if (clsChain[kpl][kcl].InTrack >= 0) continue;
                   for (kend = 0; kend < 2; ++kend) {
                     kdir = clsChain[kpl][kcl].Dir[kend];
                     if (idir != 0 && kdir != 0 && idir != kdir) continue;
                     if (clsChain[kpl][kcl].VtxIndex[kend] != vtx[ivx].EvtIndex) continue;
                     // rough check of other end match
                     // TODO: SHOWER-LIKE CLUSTER CHECK
                     match.Cls[kpl] = kcl;
                     match.End[kpl] = kend;
                     // first call to ignore redundant matches
                     if (DupMatch(match, nplanes)) continue;
                     FillEndMatch(detProp, tpcid, match);
                     // ignore if no signal at the other end
                     if (match.Chg[kpl] <= 0) continue;
                     if (match.Err + match.oErr > 100) continue;
                     // second call to keep matches with better error
                     if (DupMatch(match, nplanes)) continue;
                     matcomb.push_back(match);
                     gotkcl = true;
                     //                    break;
                   } // kend
                 }   // kk -> kcl
                 if (gotkcl) continue;
                 // look for a cluster that missed the vertex assignment
                 float best = 10;
                 short kbst = -1;
                 unsigned short kbend = 0;
                 if (prt)
                   mf::LogVerbatim("CCTM") << "VtxMatch: look for missed cluster chain in kpl";
                 for (kcl = 0; kcl < clsChain[kpl].size(); ++kcl) {
                   if (clsChain[kpl][kcl].InTrack >= 0) continue;
                   for (kend = 0; kend < 2; ++kend) {
                     kdir = clsChain[kpl][kcl].Dir[kend];
                     if (idir != 0 && kdir != 0 && idir != kdir) continue;
                     if (clsChain[kpl][kcl].VtxIndex[kend] >= 0) continue;
                     // make a rough dX cut at the match end
                     if (fabs(clsChain[kpl][kcl].X[kend] - vtx[ivx].X) > 5) continue;
                     // and at the other end
                     if (fabs(clsChain[kpl][kcl].X[1 - kend] - clsChain[ipl][icl].X[ioend]) > 50)
                       continue;
                     // check the error
                     match.Cls[kpl] = kcl;
                     match.End[kpl] = kend;
                     if (DupMatch(match, nplanes)) continue;
                     FillEndMatch(detProp, tpcid, match);
                     totErr = match.Err + match.oErr;
                     if (prt) {
                       mf::LogVerbatim myprt("CCTM");
                       myprt << "VtxMatch: Chk missing cluster match ";
                       for (unsigned short ii = 0; ii < nplanes; ++ii)
                         myprt << " " << ii << ":" << match.Cls[ii] << ":" << match.End[ii];
                       myprt << " Err " << match.Err << "\n";
                     }
                     if (totErr > 100) continue;
                     if (totErr < best) {
                       best = totErr;
                       kbst = kcl;
                       kbend = kend;
                     }
                   } // kend
                 }   // kcl
                 if (kbst >= 0) {
                   // found a decent match
                   match.Cls[kpl] = kbst;
                   match.End[kpl] = kbend;
                   FillEndMatch(detProp, tpcid, match);
                   matcomb.push_back(match);
                   // assign the vertex to this cluster
                   clsChain[kpl][kbst].VtxIndex[kbend] = ivx;
                   // and update vxCls
                   vxCls[kpl].push_back(kbst);
                 }
                 else {
                   // Try a 2 plane match if a 3 plane match didn't work
                   match.Cls[kpl] = -1;
                   match.End[kpl] = 0;
                   if (DupMatch(match, nplanes)) continue;
                   FillEndMatch(detProp, tpcid, match);
                   if (match.Err + match.oErr < 100) matcomb.push_back(match);
                 }
               } // jend
             }   // iend
           }     // jj
         }       // ii -> icl
       }         // ipl
 
       if (matcomb.size() == 0) continue;
       SortMatches(detProp, fmCluHits, 1, tpcid);
 
     } // ivx
 
     for (ipl = 0; ipl < 3; ++ipl)
       vxCls[ipl].clear();
 
   } // VtxMatch

Member Data Documentation

unsigned short trkf::CCTrackMaker::algIndex

private

Definition at line 79 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), FillEndMatch2(), FindMissingCluster(), PlnMatch(), and produce().

std::vector<art::Ptr<recob::Hit> > trkf::CCTrackMaker::allhits

private

Definition at line 124 of file CCTrackMaker_module.cc.

Referenced by ChargeNear(), FillChgNear(), FillWireHitRange(), and produce().

std::array<float, 3> trkf::CCTrackMaker::ChgNorm

private

Definition at line 207 of file CCTrackMaker_module.cc.

Referenced by ChargeNear(), FillChgNear(), FillWireHitRange(), and produce().

std::array<std::vector<clPar>, 3> trkf::CCTrackMaker::cls

private

Definition at line 146 of file CCTrackMaker_module.cc.

Referenced by dXClTraj(), FillChgNear(), FillTrkHits(), FindMaybeVertices(), MakeClusterChains(), PrintClusters(), produce(), and StoreTrack().

std::array<std::vector<ClsChainPar>, 3> trkf::CCTrackMaker::clsChain

private

Definition at line 166 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), FillEndMatch2(), FillTrkHits(), FindMissingCluster(), MakeClusterChains(), PlnMatch(), PrintClusters(), produce(), SortMatches(), StoreTrack(), and VtxMatch().

std::vector<float> trkf::CCTrackMaker::fAngleMatchErr

private

Definition at line 82 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and FillEndMatch().

float trkf::CCTrackMaker::fChainMaxdX

private

Definition at line 89 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and MakeClusterChains().

float trkf::CCTrackMaker::fChainVtxAng

private

Definition at line 90 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and MakeClusterChains().

std::vector<float> trkf::CCTrackMaker::fChgAsymFactor

private

Definition at line 83 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), FillEndMatch(), FillEndMatch2(), and PlnMatch().

float trkf::CCTrackMaker::fChgWindow {40}

private

Definition at line 95 of file CCTrackMaker_module.cc.

Referenced by FillChgNear().

std::string trkf::CCTrackMaker::fClusterModuleLabel

private

Definition at line 68 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and produce().

short trkf::CCTrackMaker::fDebugAlg

private

Definition at line 111 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), MakeClusterChains(), and produce().

short trkf::CCTrackMaker::fDebugCluster

private

Definition at line 113 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and MakeClusterChains().

short trkf::CCTrackMaker::fDebugPlane

private

Definition at line 112 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and MakeClusterChains().

float trkf::CCTrackMaker::fDeltaRayCut

private

Definition at line 99 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and TagCosmics().

float trkf::CCTrackMaker::fFiducialCut

private

Definition at line 98 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and TagCosmics().

float trkf::CCTrackMaker::fHitFitErrFac

private

Definition at line 105 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), FitVertices(), and StoreTrack().

std::string trkf::CCTrackMaker::fHitModuleLabel

private

Definition at line 67 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and produce().

std::array<unsigned int, 3> trkf::CCTrackMaker::firstHit

private

Definition at line 120 of file CCTrackMaker_module.cc.

Referenced by FillWireHitRange().

std::array<unsigned int, 3> trkf::CCTrackMaker::firstWire

private

Definition at line 118 of file CCTrackMaker_module.cc.

Referenced by FillChgNear(), and FillWireHitRange().

std::vector<bool> trkf::CCTrackMaker::fMakeAlgTracks

private

Definition at line 85 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and produce().

bool trkf::CCTrackMaker::fMakePFPs

private

Definition at line 101 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and produce().

std::vector<short> trkf::CCTrackMaker::fMatchAlgs

private

Definition at line 80 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and produce().

std::vector<float> trkf::CCTrackMaker::fMatchMinLen

private

Definition at line 84 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and PlnMatch().

float trkf::CCTrackMaker::fMaxDAng

private

Definition at line 88 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and MakeClusterChains().

float trkf::CCTrackMaker::fMaxMergeError

private

Definition at line 92 of file CCTrackMaker_module.cc.

Referenced by MakeClusterChains(), and produce().

float trkf::CCTrackMaker::fMergeChgAsym

private

Definition at line 91 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and MakeClusterChains().

float trkf::CCTrackMaker::fMergeErrorCut

private

Definition at line 93 of file CCTrackMaker_module.cc.

Referenced by MakeClusterChains(), and produce().

unsigned short trkf::CCTrackMaker::fNVtxTrkHitsFit

private

Definition at line 104 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and FitVertices().

bool trkf::CCTrackMaker::fPrintAllClusters

private

Definition at line 114 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and PrintClusters().

TrackTrajectoryAlg trkf::CCTrackMaker::fTrackTrajectoryAlg

private

Definition at line 75 of file CCTrackMaker_module.cc.

Referenced by StoreTrack().

bool trkf::CCTrackMaker::fuBCode

private

Definition at line 108 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and PlnMatch().

VertexFitAlg trkf::CCTrackMaker::fVertexFitAlg

private

Definition at line 76 of file CCTrackMaker_module.cc.

Referenced by FitVertices().

std::string trkf::CCTrackMaker::fVertexModuleLabel

private

Definition at line 69 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), and produce().

float trkf::CCTrackMaker::fWirePitch

private

Definition at line 96 of file CCTrackMaker_module.cc.

Referenced by dXClTraj(), FillChgNear(), FindMaybeVertices(), MakeClusterChains(), and produce().

std::vector<float> trkf::CCTrackMaker::fXMatchErr

private

Definition at line 81 of file CCTrackMaker_module.cc.

Referenced by CCTrackMaker(), FillEndMatch(), FillEndMatch2(), and FindMissingCluster().

art::ServiceHandle<geo::Geometry const> trkf::CCTrackMaker::geom

private

Definition at line 72 of file CCTrackMaker_module.cc.

Referenced by PlnMatch(), produce(), and TagCosmics().

std::array<unsigned int, 3> trkf::CCTrackMaker::lastHit

private

Definition at line 121 of file CCTrackMaker_module.cc.

Referenced by FillWireHitRange().

std::array<unsigned int, 3> trkf::CCTrackMaker::lastWire

private

Definition at line 119 of file CCTrackMaker_module.cc.

Referenced by FillWireHitRange().

std::vector<MatchPars> trkf::CCTrackMaker::matcomb

private

Definition at line 232 of file CCTrackMaker_module.cc.

Referenced by DupMatch(), FillTrkHits(), PlnMatch(), produce(), SortMatches(), StoreTrack(), and VtxMatch().

std::vector<unsigned short> trkf::CCTrackMaker::pfpToTrkID

private

Definition at line 212 of file CCTrackMaker_module.cc.

Referenced by MakeFamily(), produce(), and TagCosmics().

bool trkf::CCTrackMaker::prt

private

Definition at line 115 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), FillEndMatch2(), FillTrkHits(), MakeClusterChains(), PlnMatch(), produce(), SortMatches(), StoreTrack(), and VtxMatch().

std::array<std::vector<art::Ptr<recob::Hit> >, 3> trkf::CCTrackMaker::seedHits

private

Definition at line 205 of file CCTrackMaker_module.cc.

Referenced by produce().

std::vector<TrkPar> trkf::CCTrackMaker::trk

private

Definition at line 200 of file CCTrackMaker_module.cc.

std::array<std::vector<art::Ptr<recob::Hit> >, 3> trkf::CCTrackMaker::trkHits

private

Definition at line 203 of file CCTrackMaker_module.cc.

Referenced by FillTrkHits(), produce(), and StoreTrack().

std::vector<vtxPar> trkf::CCTrackMaker::vtx

private

Definition at line 179 of file CCTrackMaker_module.cc.

Referenced by FillEndMatch(), FillEndMatch2(), FindMaybeVertices(), FitVertices(), MakeClusterChains(), MakeFamily(), PrintClusters(), PrintTracks(), produce(), StoreTrack(), and VtxMatch().

std::array<std::vector<unsigned short>, 3> trkf::CCTrackMaker::vxCls

private

Definition at line 182 of file CCTrackMaker_module.cc.

Referenced by produce(), and VtxMatch().

std::array<std::vector<std::pair<int, int> >, 3> trkf::CCTrackMaker::WireHitRange

private

Definition at line 122 of file CCTrackMaker_module.cc.

Referenced by FillChgNear(), and FillWireHitRange().

geo::WireReadoutGeom const& trkf::CCTrackMaker::wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get()

private

Definition at line 73 of file CCTrackMaker_module.cc.

Referenced by FillChgNear(), FillEndMatch(), FillWireHitRange(), FindMaybeVertices(), FitVertices(), MakeClusterChains(), MakeFamily(), PlnMatch(), PrintClusters(), produce(), SortMatches(), StoreTrack(), and VtxMatch().

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

larreco/v10_01_12/source/larreco/TrackFinder/CCTrackMaker_module.cc

Classes

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation