SpacePointHit3DBuilder class definiton. More...

Inheritance diagram for lar_cluster3d::SpacePointHit3DBuilder:

Public Types
enum	TimeValues { COLLECTARTHITS = 0, BUILDTHREEDHITS = 1, NUMTIMEVALUES }
	enumerate the possible values for time checking if monitoring timing More...

using	RecobHitToPtrMap = std::map< const recob::Hit *, art::Ptr< recob::Hit >>
	Defines a structure mapping art representation to internal. More...

Public Member Functions
	SpacePointHit3DBuilder (fhicl::ParameterSet const &pset)
	Constructor. More...

	~SpacePointHit3DBuilder ()
	Destructor. More...

void	configure (const fhicl::ParameterSet &) override
	Interface for configuring the particular algorithm tool. More...

void	Hit3DBuilder (const art::Event &evt, reco::HitPairList &hitPairList, RecobHitToPtrMap &) const override
	Given a set of recob hits, run DBscan to form 3D clusters. More...

float	getTimeToExecute (IHit3DBuilder::TimeValues index) const override
	If monitoring, recover the time to execute a particular function. More...

Private Types
using	Hit2DVector = std::vector< reco::ClusterHit2D >

Private Attributes
art::InputTag	m_spacePointTag
	Data members to follow. More...

bool	m_enableMonitoring

std::vector< float >	m_timeVector

Hit2DVector	m_clusterHit2DMasterVec

geo::Geometry *	m_geometry

const detinfo::DetectorProperties *	m_detector

Detailed Description

SpacePointHit3DBuilder class definiton.

Definition at line 39 of file SpacePointHit3DBuilder_tool.cc.

Member Typedef Documentation

using lar_cluster3d::SpacePointHit3DBuilder::Hit2DVector = std::vector<reco::ClusterHit2D>

private

Definition at line 71 of file SpacePointHit3DBuilder_tool.cc.

using lar_cluster3d::IHit3DBuilder::RecobHitToPtrMap = std::map<const recob::Hit*, art::Ptr<recob::Hit>>

inherited

Defines a structure mapping art representation to internal.

Definition at line 44 of file IHit3DBuilder.h.

Member Enumeration Documentation

enum lar_cluster3d::IHit3DBuilder::TimeValues

inherited

enumerate the possible values for time checking if monitoring timing

Enumerator
COLLECTARTHITS
BUILDTHREEDHITS
NUMTIMEVALUES

Definition at line 57 of file IHit3DBuilder.h.

                     {COLLECTARTHITS   = 0,
                      BUILDTHREEDHITS  = 1,
                      NUMTIMEVALUES
     };

Constructor & Destructor Documentation

lar_cluster3d::SpacePointHit3DBuilder::SpacePointHit3DBuilder ( fhicl::ParameterSet const & pset )

explicit

Constructor.

Parameters

pset

Definition at line 88 of file SpacePointHit3DBuilder_tool.cc.

References configure().

 {
     this->configure(pset);
 }

lar_cluster3d::SpacePointHit3DBuilder::~SpacePointHit3DBuilder ( )

Destructor.

Definition at line 95 of file SpacePointHit3DBuilder_tool.cc.

96 {

97 }

Member Function Documentation

void lar_cluster3d::SpacePointHit3DBuilder::configure ( const fhicl::ParameterSet & )

overridevirtual

Interface for configuring the particular algorithm tool.

Parameters

ParameterSet The input set of parameters for configuration

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 101 of file SpacePointHit3DBuilder_tool.cc.

References fhicl::ParameterSet::get(), m_detector, m_enableMonitoring, m_geometry, and m_spacePointTag.

Referenced by SpacePointHit3DBuilder().

 {
     m_spacePointTag    = pset.get<art::InputTag>("SpacePointTag");
     m_enableMonitoring = pset.get<bool>         ("EnableMonitoring",    true);
     
     art::ServiceHandle<geo::Geometry> geometry;
     
     m_geometry = &*geometry;
     m_detector = lar::providerFrom<detinfo::DetectorPropertiesService>();
 }

float lar_cluster3d::SpacePointHit3DBuilder::getTimeToExecute ( IHit3DBuilder::TimeValues index ) const

inlineoverridevirtual

If monitoring, recover the time to execute a particular function.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 67 of file SpacePointHit3DBuilder_tool.cc.

References m_timeVector.

67 {return m_timeVector.at(index);}

lar_cluster3d::SpacePointHit3DBuilder::m_timeVector

std::vector< float > m_timeVector

Definition: SpacePointHit3DBuilder_tool.cc:79

void lar_cluster3d::SpacePointHit3DBuilder::Hit3DBuilder	(	const art::Event &	evt,
		reco::HitPairList &	hitPairList,
		RecobHitToPtrMap &	recobHitToArtPtrMap
	)		const

overridevirtual

Given a set of recob hits, run DBscan to form 3D clusters.

Parameters

hitPairList	The input list of 3D hits to run clustering on
clusterParametersList	A list of cluster objects (parameters from associated hits)

Recover the 2D hits from art and fill out the local data structures for the 3D clustering

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 113 of file SpacePointHit3DBuilder_tool.cc.

References lar_cluster3d::IHit3DBuilder::BUILDTHREEDHITS, detinfo::DetectorProperties::ConvertTicksToX(), DEFINE_ART_CLASS_TOOL, art::Ptr< T >::get(), art::DataViewImpl::getByLabel(), reco::ClusterHit2D::getHit(), detinfo::DetectorProperties::GetXTicksOffset(), art::Handle< T >::isValid(), geo::kU, geo::kV, geo::kW, m_clusterHit2DMasterVec, m_detector, m_enableMonitoring, m_geometry, m_spacePointTag, m_timeVector, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), lar_cluster3d::IHit3DBuilder::NUMTIMEVALUES, geo::PlaneID::Plane, reco::ClusterHit2D::SHAREDINTRIPLET, detinfo::DetectorProperties::TriggerOffset(), reco::ClusterHit2D::USEDINTRIPLET, weight, recob::Hit::WireID(), and recob::SpacePoint::XYZ().

 {
     m_timeVector.resize(NUMTIMEVALUES, 0.);
     
     cet::cpu_timer theClockMakeHits;
     
     if (m_enableMonitoring) theClockMakeHits.start();
 
     // Start by recovering the associations between space points and hits
     art::Handle< art::Assns<recob::Hit, recob::SpacePoint> > hitSpacePointAssnsHandle;
     evt.getByLabel(m_spacePointTag, hitSpacePointAssnsHandle);
     
     if (!hitSpacePointAssnsHandle.isValid()) return;
     
     // First step is to loop through and get a mapping between space points and associated hits
     // and, importantly, a list of unique hits (and mapping between art ptr and hit)
     using SpacePointHitVecMap = std::map<const recob::SpacePoint*, std::vector<const recob::Hit*>>;
     using RecobHitSet         = std::set<const recob::Hit*>;
     
     SpacePointHitVecMap spacePointHitVecMap;
     RecobHitSet         recobHitSet;
 
     for(auto& assnPair : *hitSpacePointAssnsHandle)
     {
         const art::Ptr<recob::SpacePoint> spacePoint = assnPair.second;
         const art::Ptr<recob::Hit>&       recobHit   = assnPair.first;
 
         spacePointHitVecMap[spacePoint.get()].push_back(recobHit.get());
         recobHitSet.insert(recobHit.get());
         recobHitToArtPtrMap[recobHit.get()] = recobHit;
     }
     
     // We'll want to correct the hit times for the plane offsets
     // (note this is already taken care of when converting to position)
     std::map<geo::PlaneID, double> planeOffsetMap;
     
     // Initialize the plane to hit vector map
     for(size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++)
     {
         for(size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++)
         {
             // What we want here are the relative offsets between the planes
             // Note that plane 0 is assumed the "first" plane and is the reference
             planeOffsetMap[geo::PlaneID(cryoIdx,tpcIdx,0)] = 0.;
             planeOffsetMap[geo::PlaneID(cryoIdx,tpcIdx,1)] = m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,1))
                                                            - m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0));
             planeOffsetMap[geo::PlaneID(cryoIdx,tpcIdx,2)] = m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,2))
                                                            - m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0));
             
             std::cout << "***> plane 0 offset: " << planeOffsetMap[geo::PlaneID(cryoIdx,tpcIdx,0)] << ", plane 1: " << planeOffsetMap[geo::PlaneID(cryoIdx,tpcIdx,1)] << ", plane 2: " << planeOffsetMap[geo::PlaneID(cryoIdx,tpcIdx,2)] << std::endl;
             std::cout << "     Det prop plane 0: " << m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0)) << ", plane 1: "  << m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,1)) << ", plane 2: " << m_detector->GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,2)) << ", Trig: " << m_detector->TriggerOffset() << std::endl;
         }
     }
 
     // We need temporary mapping from recob::Hit's to our 2D hits
     using RecobHitTo2DHitMap = std::map<const recob::Hit*,const reco::ClusterHit2D*>;
     
     RecobHitTo2DHitMap recobHitTo2DHitMap;
     
     // Set the size of the container for our hits
     m_clusterHit2DMasterVec.clear();
     m_clusterHit2DMasterVec.reserve(recobHitSet.size());
 
     // Now go throught the list of unique hits and create the 2D hits we'll use
     for(auto& recobHit : recobHitSet)
     {
         const geo::WireID& hitWireID(recobHit->WireID());
         
         double hitPeakTime(recobHit->PeakTime() - planeOffsetMap[hitWireID.planeID()]);
         double xPosition(m_detector->ConvertTicksToX(recobHit->PeakTime(), hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat));
         
         m_clusterHit2DMasterVec.emplace_back(0, 0., 0., xPosition, hitPeakTime, *recobHit);
         
         recobHitTo2DHitMap[recobHit]  = &m_clusterHit2DMasterVec.back();
     }
     
     // Now we can go through the space points and build our 3D hits
     for(auto& pointPair : spacePointHitVecMap)
     {
         const recob::SpacePoint*              spacePoint  = pointPair.first;
         const std::vector<const recob::Hit*>& recobHitVec = pointPair.second;
         
         if (recobHitVec.size() != 3)
         {
             std::cout << "************>>>>>> do not have 3 hits associated to space point! " << recobHitVec.size() << " ***************" << std::endl;
             continue;
         }
         
         std::vector<const reco::ClusterHit2D*> hit2DVec(recobHitVec.size());
         
         for(const auto& recobHit : recobHitVec)
         {
             const reco::ClusterHit2D* hit2D = recobHitTo2DHitMap.at(recobHit);
             
             hit2DVec[hit2D->getHit().WireID().Plane] = hit2D;
         }
         
         // Weighted average, delta, sigmas, chisquare, kitchen sink, refrigerator for beer, etc.
         float avePeakTime(0.);
         float weightSum(0.);
 
         for(const auto& hit2D : hit2DVec)
         {
             float hitSigma = hit2D->getHit().RMS();
             float weight   = 1. / (hitSigma * hitSigma);
 
             avePeakTime   += weight * hit2D->getTimeTicks();
             weightSum     += weight;
         }
         
         avePeakTime   /= weightSum;
 
         // Armed with the average peak time, now get hitChiSquare and the sig vec
         float              hitChiSquare(0.);
         float              sigmaPeakTime(std::sqrt(1./weightSum));
         
         for(const auto& hit2D : hit2DVec)
         {
             float hitRMS    = hit2D->getHit().RMS();
             float combRMS   = std::sqrt(hitRMS*hitRMS - sigmaPeakTime*sigmaPeakTime);
             float peakTime  = hit2D->getTimeTicks();
             float deltaTime = peakTime - avePeakTime;
             float hitSig    = deltaTime / combRMS;
             
             hitChiSquare += hitSig * hitSig;
         }
 
         // The x position is a weighted sum but the y-z position is simply the average
         float position[]  = { float(spacePoint->XYZ()[0]), float(spacePoint->XYZ()[1]), float(spacePoint->XYZ()[2])};
         float totalCharge = hit2DVec[0]->getHit().Integral() + hit2DVec[1]->getHit().Integral() + hit2DVec[2]->getHit().Integral();
             
         reco::ClusterHit2DVec hitVector;
         
         hitVector.resize(3,NULL);
             
         // Make sure we have the hits
         hitVector.at(hit2DVec[0]->getHit().WireID().Plane) = hit2DVec[0];
         hitVector.at(hit2DVec[1]->getHit().WireID().Plane) = hit2DVec[1];
         hitVector.at(hit2DVec[2]->getHit().WireID().Plane) = hit2DVec[2];
             
         // And get the wire IDs
         std::vector<geo::WireID> wireIDVec = {geo::WireID(0,0,geo::kU,0), geo::WireID(0,0,geo::kV,0), geo::WireID(0,0,geo::kW,0)};
             
         for(const auto& hit : hitVector)
         {
             wireIDVec[hit->getHit().WireID().Plane] = hit->getHit().WireID();
                 
             if (hit->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET) hit->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
                 
             hit->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
         }
             
         unsigned int statusBits(0x7);
             
         // For compiling at the moment
         std::vector<float> hitDelTSigVec = {0.,0.,0.};
 
         hitDelTSigVec[0] = std::fabs(hitVector[0]->getTimeTicks() - 0.5 * (hitVector[1]->getTimeTicks() + hitVector[2]->getTimeTicks()));
         hitDelTSigVec[1] = std::fabs(hitVector[1]->getTimeTicks() - 0.5 * (hitVector[2]->getTimeTicks() + hitVector[0]->getTimeTicks()));
         hitDelTSigVec[2] = std::fabs(hitVector[2]->getTimeTicks() - 0.5 * (hitVector[0]->getTimeTicks() + hitVector[1]->getTimeTicks()));
         
         float deltaPeakTime = *std::min_element(hitDelTSigVec.begin(),hitDelTSigVec.end());
 
         // Create the 3D cluster hit
         hitPairList.emplace_back(new reco::ClusterHit3D(0,
                                                         statusBits,
                                                         position,
                                                         totalCharge,
                                                         avePeakTime,
                                                         deltaPeakTime,
                                                         sigmaPeakTime,
                                                         hitChiSquare,
                                                         0.,
                                                         0.,
                                                         hitVector,
                                                         hitDelTSigVec,
                                                         wireIDVec));
     }
     
     if (m_enableMonitoring)
     {
         theClockMakeHits.stop();
         
         m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time();
     }
     
     mf::LogDebug("Cluster3D") << ">>>>> 3D hit building done, found " << hitPairList.size() << " 3D Hits" << std::endl;
 
     return;
 }