Inheritance diagram for reco3d::SpacePointSolver:

Public Types
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

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

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

virtual	~SpacePointSolver ()

void	produce (art::Event &evt)

void	beginJob ()

void	endJob ()

template<typename PROD , BranchType B = InEvent>
ProductID	getProductID (std::string const &instanceName={}) const

template<typename PROD , BranchType B>
ProductID	getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const

bool	modifiesEvent () const

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

template<typename T , art::BranchType BT>
art::ProductToken< T >	consumes (InputTag const &it)

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

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

template<typename T , art::BranchType BT>
art::ViewToken< T >	consumesView (InputTag const &it)

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

template<typename T , art::BranchType BT>
art::ProductToken< T >	mayConsume (InputTag const &it)

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

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

template<typename T , art::BranchType BT>
art::ViewToken< T >	mayConsumeView (InputTag const &it)

base_engine_t &	createEngine (seed_t seed)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)

seed_t	get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)

Static Public Member Functions
static cet::exempt_ptr< Consumer >	non_module_context ()

Protected Types
typedef std::map< const WireHit *, art::Ptr< recob::Hit > >	HitMap_t

Protected Member Functions
double	HitToXPos (const recob::Hit &hit, geo::TPCID tpc) const

bool	CloseDrift (double xa, double xb) const

bool	CloseSpace (geo::WireIDIntersection ra, geo::WireIDIntersection rb) const

void	FastForward (std::vector< InductionWireWithXPos >::iterator &it, double target, const std::vector< InductionWireWithXPos >::const_iterator &end) const

bool	ISect (int chanA, int chanB, geo::TPCID tpc, geo::WireIDIntersection &pt) const

bool	ISect (int chanA, int chanB, geo::TPCID tpc) const

void	AddNeighbours (const std::vector< SpaceCharge * > &spaceCharges) const

void	BuildSystemXUV (const std::vector< art::Ptr< recob::Hit >> &xhits, const std::vector< art::Ptr< recob::Hit >> &uhits, const std::vector< art::Ptr< recob::Hit >> &vhits, std::vector< CollectionWireHit * > &cwires, std::vector< InductionWireHit * > &iwires, bool incNei, HitMap_t &hitmap) const

void	BuildSystemXU (const std::vector< art::Ptr< recob::Hit >> &xhits, const std::vector< art::Ptr< recob::Hit >> &uhits, std::vector< CollectionWireHit * > &cwires, std::vector< InductionWireHit * > &iwires, bool incNei, HitMap_t &hitmap) const

bool	AddSpacePoint (const SpaceCharge &sc, int id, recob::ChargedSpacePointCollectionCreator &points) const
	return whether the point was inserted (only happens when it has charge) More...

void	FillSystemToSpacePoints (const std::vector< CollectionWireHit * > &cwires, recob::ChargedSpacePointCollectionCreator &pts) const

void	FillSystemToSpacePointsAndAssns (const std::vector< CollectionWireHit * > &cwires, const HitMap_t &hitmap, recob::ChargedSpacePointCollectionCreator &points, art::Assns< recob::SpacePoint, recob::Hit > &assn) const

CurrentProcessingContext const *	currentContext () const

void	validateConsumedProduct (BranchType const bt, ProductInfo const &pi)

void	prepareForJob (fhicl::ParameterSet const &pset)

void	showMissingConsumes () const

Protected Attributes
std::string	fHitLabel

bool	fFit

double	fAlpha

TH1 *	fDeltaX

const detinfo::DetectorProperties *	detprop

const geo::GeometryCore *	geom

Detailed Description

Definition at line 60 of file SpacePointSolver_module.cc.

Member Typedef Documentation

typedef std::map<const WireHit*, art::Ptr<recob::Hit> > reco3d::SpacePointSolver::HitMap_t

protected

Definition at line 90 of file SpacePointSolver_module.cc.

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>

inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

reco3d::SpacePointSolver::SpacePointSolver ( const fhicl::ParameterSet & pset )

explicit

Definition at line 135 of file SpacePointSolver_module.cc.

References recob::ChargedSpacePointCollectionCreator::produces().

   : fHitLabel(pset.get<std::string>("HitLabel")),
     fFit(pset.get<bool>("Fit")),
     fAlpha(pset.get<double>("Alpha"))
 {
   recob::ChargedSpacePointCollectionCreator::produces(*this, "pre");
   if(fFit){
     recob::ChargedSpacePointCollectionCreator::produces(*this);
     produces<art::Assns<recob::SpacePoint, recob::Hit>>();
     recob::ChargedSpacePointCollectionCreator::produces(*this, "noreg");
   }
 }

reco3d::SpacePointSolver::~SpacePointSolver ( )

virtual

Definition at line 149 of file SpacePointSolver_module.cc.

150 {

151 }

Member Function Documentation

void reco3d::SpacePointSolver::AddNeighbours ( const std::vector< SpaceCharge * > & spaceCharges ) const

protected

Definition at line 270 of file SpacePointSolver_module.cc.

References BuildSystemXUV(), Neighbour::fCoupling, SpaceCharge::fPred, Neighbour::fSC, SpaceCharge::fX, SpaceCharge::fY, SpaceCharge::fZ, reco3d::InductionWireWithXPos::operator<(), sqr(), x, y, and z.

Referenced by HitToXPos().

 {
   static const double kCritDist = 5;
 
   // Could use a QuadTree or VPTree etc, but seems like overkill
   class IntCoord
   {
   public:
     IntCoord(const SpaceCharge& sc)
       : fX(sc.fX/kCritDist),
         fY(sc.fY/kCritDist),
         fZ(sc.fZ/kCritDist)
     {
     }
 
     bool operator<(const IntCoord& i) const
     {
       return std::make_tuple(fX, fY, fZ) < std::make_tuple(i.fX, i.fY, i.fZ);
     }
 
     std::vector<IntCoord> Neighbours() const
     {
       std::vector<IntCoord> ret;
       for(int dx = -1; dx <= +1; ++dx){
         for(int dy = -1; dy <= +1; ++dy){
           for(int dz = -1; dz <= +1; ++dz){
             ret.push_back(IntCoord(fX+dx, fY+dy, fZ+dz));
           }
         }
       }
       return ret;
     }
   protected:
     IntCoord(int x, int y, int z) : fX(x), fY(y), fZ(z) {}
 
     int fX, fY, fZ;
   };
 
   std::map<IntCoord, std::vector<SpaceCharge*>> scMap;
   for(SpaceCharge* sc: spaceCharges){
     scMap[IntCoord(*sc)].push_back(sc);
   }
 
   std::cout << "Neighbour search..." << std::endl;
   std::cout << spaceCharges.size() << std::endl;
   // Now that we know all the space charges, can go through and assign neighbours
 
   int Ntests = 0;
   int Nnei = 0;
   for(SpaceCharge* sc1: spaceCharges){
     IntCoord ic(*sc1);
     for(IntCoord icn: ic.Neighbours()){
       for(SpaceCharge* sc2: scMap[icn]){
 
         ++Ntests;
 
         if(sc1 == sc2) continue;
         /*const*/ double dist2 = sqr(sc1->fX-sc2->fX) + sqr(sc1->fY-sc2->fY) + sqr(sc1->fZ-sc2->fZ);
 
         if(dist2 > sqr(kCritDist)) continue;
 
         if(dist2 == 0){
           std::cout << "ZERO DISTANCE SOMEHOW?" << std::endl;
           std::cout << sc1->fCWire << " " << sc1->fWire1 << " " << sc1->fWire2 << std::endl;
           std::cout << sc2->fCWire << " " << sc2->fWire1 << " " << sc2->fWire2 << std::endl;
           std::cout << dist2 << " " << sc1->fX << " " << sc2->fX << " " << sc1->fY << " " << sc2->fY << " " << sc1->fZ << " " << sc2->fZ << std::endl;
           continue;
           dist2 = sqr(kCritDist);
         }
 
         ++Nnei;
 
         // This is a pretty random guess
         const double coupling = exp(-sqrt(dist2)/2);
         sc1->fNeighbours.emplace_back(sc2, coupling);
 
         if(isnan(1/sqrt(dist2)) || isinf(1/sqrt(dist2))){
           std::cout << dist2 << " " << sc1->fX << " " << sc2->fX << " " << sc1->fY << " " << sc2->fY << " " << sc1->fZ << " " << sc2->fZ << std::endl;
           abort();
         }
       } // end for sc2
     } // end for icn
 
       // The neighbours lists use the most memory, so be careful to trim
     sc1->fNeighbours.shrink_to_fit();
   } // end for sc1
 
   for(SpaceCharge* sc: spaceCharges){
     for(Neighbour& nei: sc->fNeighbours){
       sc->fNeiPotential += nei.fCoupling * nei.fSC->fPred;
     }
   }
 
   std::cout << Ntests << " tests to find " << Nnei << std::endl;
 }

bool reco3d::SpacePointSolver::AddSpacePoint	(	const SpaceCharge &	sc,
		int	id,
		recob::ChargedSpacePointCollectionCreator &	points
	)		const

protected

return whether the point was inserted (only happens when it has charge)

Definition at line 640 of file SpacePointSolver_module.cc.

References recob::ChargedSpacePointCollectionCreator::add(), FillSystemToSpacePoints(), SpaceCharge::fPred, SpaceCharge::fX, SpaceCharge::fY, and SpaceCharge::fZ.

Referenced by BuildSystemXU().

 {
   static const double err[6] = {0,};
 
   const float charge = sc.fPred;
   if(charge == 0) return false;
 
   const double xyz[3] = {sc.fX, sc.fY, sc.fZ};
   points.add({ xyz, err, 0.0, id }, charge);
 
   return true;
 }

void reco3d::SpacePointSolver::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 154 of file SpacePointSolver_module.cc.

References art::TFileDirectory::make().

 {
   detprop = art::ServiceHandle<detinfo::DetectorPropertiesService>()->provider();
   geom = art::ServiceHandle<geo::Geometry>()->provider();
 
   art::ServiceHandle<art::TFileService> tfs;
   fDeltaX = tfs->make<TH1F>("deltax", ";#Deltax (cm)", 100, -2, +2);
 }

void reco3d::SpacePointSolver::BuildSystemXU	(	const std::vector< art::Ptr< recob::Hit >> &	xhits,
		const std::vector< art::Ptr< recob::Hit >> &	uhits,
		std::vector< CollectionWireHit * > &	cwires,
		std::vector< InductionWireHit * > &	iwires,
		bool	incNei,
		HitMap_t &	hitmap
	)		const

protected

Definition at line 547 of file SpacePointSolver_module.cc.

References AddSpacePoint(), evd::details::end(), InductionWireHit::fChannel, reco3d::InductionWireWithXPos::iwire, reco3d::InductionWireWithXPos::xpos, geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by BuildSystemXUV().

 {
   std::map<geo::TPCID, std::vector<art::Ptr<recob::Hit>>> xhits_by_tpc;
   for(auto& xhit: xhits){
     const std::vector<geo::TPCID> tpcs = geom->ROPtoTPCs(geom->ChannelToROP(xhit->Channel()));
     assert(tpcs.size() == 1);
     const geo::TPCID tpc = tpcs[0];
     xhits_by_tpc[tpc].push_back(xhit);
   }
 
   // Maps from TPC to the induction wires. Normally want to access them this
   // way.
   std::map<geo::TPCID, std::vector<InductionWireWithXPos>> uwires;
 
   for(auto& hit: uhits){
     const std::vector<geo::TPCID> tpcs = geom->ROPtoTPCs(geom->ChannelToROP(hit->Channel()));
 
     InductionWireHit* iwire = new InductionWireHit(hit->Channel(), hit->Integral());
     iwires.emplace_back(iwire);
     hitmap[iwire] = hit;
 
     for(geo::TPCID tpc: tpcs){
       if(xhits_by_tpc.count(tpc) == 0) continue;
 
       const double xpos = HitToXPos(*hit, tpc);
 
       uwires[tpc].emplace_back(iwire, xpos);
     } // end for tpc
   } // end for hit
 
   for(auto it = uwires.begin(); it != uwires.end(); ++it){
     std::sort(it->second.begin(), it->second.end());
   }
 
   for(auto it = xhits_by_tpc.begin(); it != xhits_by_tpc.end(); ++it){
     const geo::TPCID tpc = it->first;
     std::sort(it->second.begin(), it->second.end(),
               [this, tpc](art::Ptr<recob::Hit>& a, art::Ptr<recob::Hit>& b)
               {
                 return HitToXPos(*a, tpc) < HitToXPos(*b, tpc);
               });
   }
 
 
   std::cout << "Finding UV coincidences..." << std::endl;
   std::vector<SpaceCharge*> spaceCharges;
 
   for(auto it: xhits_by_tpc){
     const geo::TPCID tpc = it.first;
 
     auto uwires_begin = uwires[tpc].begin();
 
     for(auto& hit: it.second){
       const double xpos = HitToXPos(*hit, tpc);
 
       std::vector<SpaceCharge*> crossers;
 
       FastForward(uwires_begin, xpos, uwires[tpc].end());
 
       // Figure out which uwires intersect this xwire here.
       for(auto uit = uwires_begin; uit != uwires[tpc].end(); ++uit){
         InductionWireWithXPos uwire = *uit;
 
         if(uwire.xpos > xpos && !CloseDrift(uwire.xpos, xpos)) break;
 
         geo::WireIDIntersection ptXU;
         if(ISect(hit->Channel(), uwire.iwire->fChannel, tpc, ptXU)){
 
           // Don't have a cwire object yet, set it later
           SpaceCharge* sc = new SpaceCharge((xpos+uwire.xpos)/2, ptXU.y, ptXU.z,
                                             0, 0, uwire.iwire);
           spaceCharges.push_back(sc);
           crossers.push_back(sc);
         } // end for uwire
       } // end for uit
 
       CollectionWireHit* cwire = new CollectionWireHit(hit->Channel(), hit->Integral(), crossers);
       hitmap[cwire] = hit;
       cwires.push_back(cwire);
       for(SpaceCharge* sc: crossers) sc->fCWire = cwire;
     } // end for hit
   } // end for it (tpc)
 
   if(incNei) AddNeighbours(spaceCharges);
 }

void reco3d::SpacePointSolver::BuildSystemXUV	(	const std::vector< art::Ptr< recob::Hit >> &	xhits,
		const std::vector< art::Ptr< recob::Hit >> &	uhits,
		const std::vector< art::Ptr< recob::Hit >> &	vhits,
		std::vector< CollectionWireHit * > &	cwires,
		std::vector< InductionWireHit * > &	iwires,
		bool	incNei,
		HitMap_t &	hitmap
	)		const

protected

Definition at line 368 of file SpacePointSolver_module.cc.

References BuildSystemXU(), evd::details::end(), InductionWireHit::fChannel, reco3d::InductionWireWithXPos::iwire, geo::kU, geo::kV, reco3d::InductionWireWithXPos::operator<(), x, reco3d::InductionWireWithXPos::xpos, geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by AddNeighbours().

 {
   std::map<geo::TPCID, std::vector<art::Ptr<recob::Hit>>> xhits_by_tpc;
   for(auto& xhit: xhits){
     const std::vector<geo::TPCID> tpcs = geom->ROPtoTPCs(geom->ChannelToROP(xhit->Channel()));
     assert(tpcs.size() == 1);
     const geo::TPCID tpc = tpcs[0];
     xhits_by_tpc[tpc].push_back(xhit);
   }
 
   // Maps from TPC to the induction wires. Normally want to access them this
   // way.
   std::map<geo::TPCID, std::vector<InductionWireWithXPos>> uwires, vwires;
 
   for(auto& ihits: {uhits, vhits}){
       for(auto& hit: ihits){
       const std::vector<geo::TPCID> tpcs = geom->ROPtoTPCs(geom->ChannelToROP(hit->Channel()));
 
       // TODO: Empirically, total collection charge is about 5% high of total
       // induction charge, which might cause "spare" charge to go where it's
       // not wanted.
       InductionWireHit* iwire = new InductionWireHit(hit->Channel(), hit->Integral() * .95);
       iwires.emplace_back(iwire);
       hitmap[iwire] = hit;
 
       for(geo::TPCID tpc: tpcs){
         if(xhits_by_tpc.count(tpc) == 0) continue;
 
         const double xpos = HitToXPos(*hit, tpc);
 
         if(hit->View() == geo::kU) uwires[tpc].emplace_back(iwire, xpos);
         if(hit->View() == geo::kV) vwires[tpc].emplace_back(iwire, xpos);
       } // end for tpc
     } // end for hit
   } // end for U/V
 
   for(auto it = uwires.begin(); it != uwires.end(); ++it){
     std::sort(it->second.begin(), it->second.end());
   }
   for(auto it = vwires.begin(); it != vwires.end(); ++it){
     std::sort(it->second.begin(), it->second.end());
   }
 
   for(auto it = xhits_by_tpc.begin(); it != xhits_by_tpc.end(); ++it){
     const geo::TPCID tpc = it->first;
     std::sort(it->second.begin(), it->second.end(),
               [this, tpc](art::Ptr<recob::Hit>& a, art::Ptr<recob::Hit>& b)
               {
                 return HitToXPos(*a, tpc) < HitToXPos(*b, tpc);
               });
   }
 
 
   struct UVCrossing
   {
     geo::TPCID tpc;
     InductionWireHit *u, *v;
 
     bool operator<(const UVCrossing& x) const
     {
       return std::make_tuple(tpc, u, v) < std::make_tuple(x.tpc, x.u, x.v);
     }
   };
 
   // Build a table of UV crossers up front
   std::cout << "Building UV table..." << std::endl;
   std::map<UVCrossing, bool> isectUV;
   std::map<UVCrossing, geo::WireIDIntersection> ptsUV;
 
   for(auto it: uwires){
     const geo::TPCID tpc = it.first;
 
     auto vwires_begin = vwires[tpc].begin();
 
     for(InductionWireWithXPos uwire: uwires[tpc]){
 
       // Fast-forward up to the first vwire that could be relevant
       FastForward(vwires_begin, uwire.xpos, vwires[tpc].end());
 
       for(auto vit = vwires_begin; vit != vwires[tpc].end(); ++vit){
         const InductionWireWithXPos vwire = *vit;
 
         // No more vwires can be relevant, bail out
         if(vwire.xpos > uwire.xpos &&
            !CloseDrift(uwire.xpos, vwire.xpos)) break;
 
         const UVCrossing key = {tpc, uwire.iwire, vwire.iwire};
 
         isectUV[key] = ISect(uwire.iwire->fChannel, vwire.iwire->fChannel, tpc,
                              ptsUV[key]);
       } // end for vwire
     } // end for uwire
   } // end for tpc
 
   std::cout << "Finding XUV coincidences..." << std::endl;
   std::vector<SpaceCharge*> spaceCharges;
 
   for(auto it: xhits_by_tpc){
     const geo::TPCID tpc = it.first;
 
     auto uwires_begin = uwires[tpc].begin();
     auto vwires_begin = vwires[tpc].begin();
 
     for(auto& hit: it.second){
       const double xpos = HitToXPos(*hit, tpc);
 
       FastForward(uwires_begin, xpos, uwires[tpc].end());
       FastForward(vwires_begin, xpos, vwires[tpc].end());
 
       // Figure out which vwires intersect this xwire here so we don't do N^2
       // nesting inside the uwire loop below.
       std::vector<InductionWireWithXPos> vwires_cross;
       std::unordered_map<InductionWireHit*, geo::WireIDIntersection> ptsXV;
       vwires_cross.reserve(vwires[tpc].size()); // avoid reallocations
       for(auto vit = vwires_begin; vit != vwires[tpc].end(); ++vit){
         InductionWireWithXPos vwire = *vit;
 
         if(vwire.xpos > xpos && !CloseDrift(vwire.xpos, xpos)) break;
 
         if(ISect(hit->Channel(), vwire.iwire->fChannel, tpc, ptsXV[vwire.iwire]))
           vwires_cross.push_back(vwire);
       } // end for vwire
 
       std::vector<SpaceCharge*> crossers;
       for(auto uit = uwires_begin; uit != uwires[tpc].end(); ++uit){
         const InductionWireWithXPos uwire = *uit;
 
         if(uwire.xpos > xpos && !CloseDrift(uwire.xpos, xpos)) break;
 
         geo::WireIDIntersection ptXU;
         if(!ISect(hit->Channel(), uwire.iwire->fChannel, tpc, ptXU)) continue;
 
         for(const InductionWireWithXPos& vwire: vwires_cross){
 
           const geo::WireIDIntersection ptXV = ptsXV[vwire.iwire];
           if(!CloseSpace(ptXU, ptXV)) continue;
 
           if(!isectUV[{tpc, uwire.iwire, vwire.iwire}]) continue;
 
           const geo::WireIDIntersection ptUV = ptsUV[{tpc, uwire.iwire, vwire.iwire}];
 
           if(!CloseSpace(ptXU, ptUV) ||
              !CloseSpace(ptXV, ptUV)) continue;
 
           fDeltaX->Fill(xpos-uwire.xpos);
           fDeltaX->Fill(xpos-vwire.xpos);
           fDeltaX->Fill(uwire.xpos-vwire.xpos);
 
           // TODO exactly which 3D position to set for this point? This average
           // aleviates the problem with a single collection wire matching
           // multiple hits on an induction wire at different times.
 
           // Don't have a cwire object yet, set it later
           SpaceCharge* sc = new SpaceCharge((xpos+uwire.xpos+vwire.xpos)/3,
                                             ptXU.y, ptXU.z,
                                             0, uwire.iwire, vwire.iwire);
           spaceCharges.push_back(sc);
           crossers.push_back(sc);
         } // end for vwire
       } // end for uwire
 
       CollectionWireHit* cwire = new CollectionWireHit(hit->Channel(), hit->Integral(), crossers);
       hitmap[cwire] = hit;
       cwires.push_back(cwire);
       for(SpaceCharge* sc: crossers) sc->fCWire = cwire;
     } // end for hit
   } // end for it (tpc)
 
   if(incNei) AddNeighbours(spaceCharges);
 }

bool reco3d::SpacePointSolver::CloseDrift	(	double	xa,
		double	xb
	)		const

protected

Definition at line 220 of file SpacePointSolver_module.cc.

 {
   // Used to cut at 10 ticks (for basically empirical reasons). Reproduce that
   // in x.
   // Sampling rate is in ns/ticks
   // Drift velocity is in cm/us
   //  static const double k = 10*detprop->SamplingRate()*1e-3*detprop->DriftVelocity();
   const double k = 0.2;//0.4; // 1 sigma on deltax plot
 
   // TODO - figure out cut value
   return fabs(xa-xb) < k;
 }

bool reco3d::SpacePointSolver::CloseSpace	(	geo::WireIDIntersection	ra,
		geo::WireIDIntersection	rb
	)		const

protected

Definition at line 234 of file SpacePointSolver_module.cc.

References FastForward(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

 {
   TVector3 pa(ra.y, ra.z, 0);
   TVector3 pb(rb.y, rb.z, 0);
 
   // TODO - figure out cut value. Empirically .25 is a bit small
   //  return (pa-pb).Mag() < .5;
 
   return (pa-pb).Mag() < .35;
 }

template<typename T , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::consumes ( InputTag const & it )

inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::consumesMany ( )

inherited

Definition at line 162 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::consumesView ( InputTag const & it )

inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t seed )

inherited

Definition at line 26 of file EngineCreator.cc.

References art::EngineCreator::rng().

Referenced by evgen::CosmicsGen::CosmicsGen(), rndm::NuRandomService::createEngine(), cluster::fuzzyCluster::fuzzyCluster(), cluster::HoughLineFinder::HoughLineFinder(), art::MixFilter< T >::initEngine_(), larg4::LArG4::LArG4(), evgen::LightSource::LightSource(), evgen::NeutronOsc::NeutronOsc(), evgen::NucleonDecay::NucleonDecay(), opdet::OpMCDigi::OpMCDigi(), opdet::OptDetDigitizer::OptDetDigitizer(), phot::PhotonLibraryPropagation::PhotonLibraryPropagation(), detsim::SimDriftElectrons::SimDriftElectrons(), evgen::SingleGen::SingleGen(), evgen::SNNueAr40CCGen::SNNueAr40CCGen(), ToyOneShowerGen::ToyOneShowerGen(), and trkf::Track3DKalman::Track3DKalman().

 {
   return rng()->createEngine(placeholder_schedule_id(), seed);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make
	)

inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make,
		label_t const &	engine_label
	)

inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
 }

CurrentProcessingContext const * art::EDProducer::currentContext ( ) const

protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

   {
     return current_context_.get();
   }

void reco3d::SpacePointSolver::endJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 164 of file SpacePointSolver_module.cc.

165 {

166 }

void reco3d::SpacePointSolver::FastForward	(	std::vector< InductionWireWithXPos >::iterator &	it,
		double	target,
		const std::vector< InductionWireWithXPos >::const_iterator &	end
	)		const

protected

Definition at line 248 of file SpacePointSolver_module.cc.

Referenced by CloseSpace().

 {
   while(it != end && it->xpos < target && !CloseDrift(it->xpos, target)) ++it;
 }

void reco3d::SpacePointSolver::FillSystemToSpacePoints	(	const std::vector< CollectionWireHit * > &	cwires,
		recob::ChargedSpacePointCollectionCreator &	pts
	)		const

protected

Definition at line 657 of file SpacePointSolver_module.cc.

References FillSystemToSpacePointsAndAssns().

Referenced by AddSpacePoint().

 {
   int iPoint = 0;
   for(const CollectionWireHit* cwire: cwires){
     for(const SpaceCharge* sc: cwire->fCrossings){
       AddSpacePoint(*sc, iPoint++, points);
     } // for sc
   } // for cwire
 }

void reco3d::SpacePointSolver::FillSystemToSpacePointsAndAssns	(	const std::vector< CollectionWireHit * > &	cwires,
		const HitMap_t &	hitmap,
		recob::ChargedSpacePointCollectionCreator &	points,
		art::Assns< recob::SpacePoint, recob::Hit > &	assn
	)		const

protected

Definition at line 671 of file SpacePointSolver_module.cc.

References art::Assns< L, R, D >::addSingle(), SpaceCharge::fWire1, SpaceCharge::fWire2, and recob::ChargedSpacePointCollectionCreator::lastSpacePointPtr().

Referenced by FillSystemToSpacePoints().

 {
   int iPoint = 0;
   for(const CollectionWireHit* cwire: cwires){
     for(const SpaceCharge* sc: cwire->fCrossings){
       // fill the space point and reconstructed charge information;
       // if the point is filtered out, it's not inserted (no association either)
       if(!AddSpacePoint(*sc, iPoint++, points)) continue;
 
       // now fill the associations to the last added space point
       const auto& spsPtr = points.lastSpacePointPtr();
 
       const auto& hit = hitmap.at(cwire);
       assn.addSingle(spsPtr, hit);
 
       if(sc->fWire1){
         assn.addSingle(spsPtr, hitmap.at(sc->fWire1));
       }
       if(sc->fWire2){
         assn.addSingle(spsPtr, hitmap.at(sc->fWire2));
       }
     } // for sc
   } // for cwire
 }

EngineCreator::seed_t EngineCreator::get_seed_value	(	fhicl::ParameterSet const &	pset,
		char const	key[] = `"seed"`,
		seed_t const	implicit_seed = `-1`
	)

inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

 {
   auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
   return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
 }

template<typename PROD , BranchType B>

ProductID art::EDProducer::getProductID ( std::string const & instanceName = {} ) const

inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

   {
     return ProducerBase::getProductID<PROD, B>(moduleDescription_,
                                                instanceName);
   }

template<typename PROD , BranchType B>

ProductID art::ProducerBase::getProductID	(	ModuleDescription const &	moduleDescription,
		std::string const &	instanceName
	)		const

inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

   {
     auto const& pd =
       get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
     return pd.productID();
   }

double reco3d::SpacePointSolver::HitToXPos	(	const recob::Hit &	hit,
		geo::TPCID	tpc
	)		const

protected

Definition at line 256 of file SpacePointSolver_module.cc.

References AddNeighbours(), recob::Hit::Channel(), recob::Hit::PeakTime(), and w.

 {
   const std::vector<geo::WireID> ws = geom->ChannelToWire(hit.Channel());
   for(geo::WireID w: ws){
     if(geo::TPCID(w) == tpc)
       return detprop->ConvertTicksToX(hit.PeakTime(), w);
   }
 
   std::cout << "Wire does not exist on given TPC!" << std::endl;
   abort();
 }

bool reco3d::SpacePointSolver::ISect	(	int	chanA,
		int	chanB,
		geo::TPCID	tpc,
		geo::WireIDIntersection &	pt
	)		const

protected

Definition at line 170 of file SpacePointSolver_module.cc.

References pt.

 {
   // string has a default implementation for hash. Perhaps TPCID should
   // implement a hash function directly.
   typedef std::tuple<int, int, std::string/*geo::TPCID*/> Key_t;
   static std::unordered_map<Key_t, bool> bCache;
   static std::unordered_map<Key_t, geo::WireIDIntersection> pCache;
 
   // Prevent cache from growing without bound
   if(bCache.size() > 1e8){
     std::cout << "Clearing Intersection caches" << std::endl;
     bCache.clear();
     pCache.clear();
   }
 
   const Key_t key = std::make_tuple(chanA, chanB, std::string(tpc));
   if(bCache.count(key)){
     if(bCache[key]) pt = pCache[key];
     return bCache[key];
   }
 
   const std::vector<geo::WireID> awires = geom->ChannelToWire(chanA);
   const std::vector<geo::WireID> bwires = geom->ChannelToWire(chanB);
 
   for(geo::WireID awire: awires){
     if(geo::TPCID(awire) != tpc) continue;
     for(geo::WireID bwire: bwires){
       if(geo::TPCID(bwire) != tpc) continue;
 
       if(geom->WireIDsIntersect(awire, bwire, pt)){
         bCache[key] = true;
         pCache[key] = pt;
         return true;
       }
     }
   }
 
   bCache[key] = false;
   return false;
 }

bool reco3d::SpacePointSolver::ISect	(	int	chanA,
		int	chanB,
		geo::TPCID	tpc
	)		const

protected

Definition at line 213 of file SpacePointSolver_module.cc.

 {
   geo::WireIDIntersection junk;
   return ISect(chanA, chanB, tpc, junk);
 }

template<typename T , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::mayConsume ( InputTag const & it )

inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::mayConsumeMany ( )

inherited

Definition at line 205 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const & it )

inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

bool art::ProducerBase::modifiesEvent ( ) const

inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

     {
       return true;
     }

cet::exempt_ptr< art::Consumer > art::Consumer::non_module_context ( )

staticinherited

Definition at line 76 of file Consumer.cc.

 {
   static Consumer invalid{InvalidTag{}};
   return &invalid;
 }

void art::Consumer::prepareForJob ( fhicl::ParameterSet const & pset )

protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

 {
   if (!moduleContext_)
     return;
 
   pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
   for (auto& consumablesPerBranch : consumables_) {
     cet::sort_all(consumablesPerBranch);
   }
 }

void reco3d::SpacePointSolver::produce ( art::Event & evt )

virtual

Implements art::EDProducer.

Definition at line 701 of file SpacePointSolver_module.cc.

References e, art::fill_ptr_vector(), art::DataViewImpl::getByLabel(), hits(), Iterate(), geo::kCollection, geo::kU, geo::kV, geo::kY, geo::kZ, Metric(), art::Event::put(), and recob::ChargedSpacePointCollectionCreator::put().

 {
   art::Handle<std::vector<recob::Hit>> hits;
   std::vector<art::Ptr<recob::Hit> > hitlist;
   if(evt.getByLabel(fHitLabel, hits))
     art::fill_ptr_vector(hitlist, hits);
 
   recob::ChargedSpacePointCollectionCreator spcol_pre(evt, *this, "pre");
   recob::ChargedSpacePointCollectionCreator spcol_noreg(evt, *this, "noreg");
   recob::ChargedSpacePointCollectionCreator spcol(evt, *this);
   auto assns = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
 
   // Skip very small events
   if(hits->size() < 20){
     spcol_pre.put();
     if(fFit){
       spcol.put();
       evt.put(std::move(assns));
       spcol_noreg.put();
     }
     return;
   }
 
   art::ServiceHandle<geo::Geometry> geom;
 
   bool is2view = false;
   std::vector<art::Ptr<recob::Hit>> xhits, uhits, vhits;
   for(auto& hit: hitlist){
     if(isnan(hit->Integral()) || isinf(hit->Integral())){
       std::cout << "WARNING: bad recob::Hit::Integral() = "
                 << hit->Integral()
                 << ". Skipping." << std::endl;
       continue;
     }
 
     if(hit->SignalType() == geo::kCollection){
       // For DualPhase, both view are collection. Arbitrarily map V to the main
       // "X" view. For Argoneut and Lariat, collection=V is also the right
       // convention.
       if(hit->View() == geo::kZ){
         xhits.push_back(hit);
       }
       if(hit->View() == geo::kV){
         xhits.push_back(hit);
         is2view = true;
       }
       if(hit->View() == geo::kU || hit->View() == geo::kY){
         uhits.push_back(hit);
         is2view = true;
       }
     }
     else{
       if(hit->View() == geo::kU) uhits.push_back(hit);
       if(hit->View() == geo::kV) vhits.push_back(hit);
     }
   } // end for hit
 
   std::vector<CollectionWireHit*> cwires;
   // So we can find them all to free the memory
   std::vector<InductionWireHit*> iwires;
 
   HitMap_t hitmap;
   if(is2view)
     BuildSystemXU(xhits, uhits, cwires, iwires, fAlpha != 0, hitmap);
   else
     BuildSystemXUV(xhits, uhits, vhits, cwires, iwires, fAlpha != 0, hitmap);
 
   FillSystemToSpacePoints(cwires, spcol_pre);
   spcol_pre.put();
 
   if(fFit){
     std::cout << "Iterating..." << std::endl;
     double prevMetric = Metric(cwires, 0);//fAlpha);
     std::cout << "Begin: " << prevMetric << std::endl;
     for(int i = 0;; ++i){
       Iterate(cwires, 0);//fAlpha);
       const double metric = Metric(cwires, 0);//fAlpha);
       std::cout << i << " " << metric << std::endl;
       if(fabs(metric-prevMetric) < 1e-3*fabs(prevMetric)) break;
       if(i > 100) break;
       //    if(metric/prevMetric > .9999) break;
       prevMetric = metric;
     }
 
     FillSystemToSpacePoints(cwires, spcol_noreg);
     spcol_noreg.put();
 
     prevMetric = Metric(cwires, fAlpha);
     std::cout << "Begin: " << prevMetric << std::endl;
     for(int i = 0;; ++i){
       Iterate(cwires, fAlpha);
       const double metric = Metric(cwires, fAlpha);
       std::cout << i << " " << metric << std::endl;
       if(fabs(metric-prevMetric) < 1e-3*fabs(prevMetric)) break;
       if(i > 100) break;
       //    if(metric/prevMetric > .9999) break;
       prevMetric = metric;
     }
 
     FillSystemToSpacePointsAndAssns(cwires, hitmap, spcol, *assns);
     spcol.put();
     evt.put(std::move(assns));
   } // end if fFit
 
   for(InductionWireHit* i: iwires) delete i;
   for(CollectionWireHit* c: cwires) delete c;
 }

void art::Consumer::showMissingConsumes ( ) const

protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

 {
   if (!moduleContext_)
     return;
 
   // If none of the branches have missing consumes statements, exit early.
   if (std::all_of(cbegin(missingConsumes_),
                   cend(missingConsumes_),
                   [](auto const& perBranch) { return perBranch.empty(); }))
     return;
 
   constexpr cet::HorizontalRule rule{60};
   mf::LogPrint log{"MTdiagnostics"};
   log << '\n'
       << rule('=') << '\n'
       << "The following consumes (or mayConsume) statements are missing from\n"
       << module_context(moduleDescription_) << '\n'
       << rule('-') << '\n';
 
   cet::for_all_with_index(
     missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
       for (auto const& pi : perBranch) {
         log << "  "
             << assemble_consumes_statement(static_cast<BranchType>(i), pi)
             << '\n';
       }
     });
   log << rule('=');
 }

void art::Consumer::validateConsumedProduct	(	BranchType const	bt,
		ProductInfo const &	pi
	)

protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

 {
   // Early exits if consumes tracking has been disabled or if the
   // consumed product is an allowed consumable.
   if (!moduleContext_)
     return;
 
   if (cet::binary_search_all(consumables_[bt], pi))
     return;
 
   if (requireConsumes_) {
     throw Exception(errors::ProductRegistrationFailure,
                     "Consumer: an error occurred during validation of a "
                     "retrieved product\n\n")
       << "The following consumes (or mayConsume) statement is missing from\n"
       << module_context(moduleDescription_) << ":\n\n"
       << "  " << assemble_consumes_statement(bt, pi) << "\n\n";
   }
 
   missingConsumes_[bt].insert(pi);
 }