Given a set of recob hits, run DBscan to form 3D clusters.
Recover the 2D hits from art and fill out the local data structures for the 3D clustering
230 cet::cpu_timer theClockMakeHits;
235 art::Handle<art::Assns<recob::Hit, recob::SpacePoint>> hitSpacePointAssnsHandle;
238 if (!hitSpacePointAssnsHandle.isValid())
return;
246 using SpacePointToWireIDMap = std::unordered_map<const recob::SpacePoint*, geo::WireID>;
248 SpacePointToWireIDMap spacePointToWireIDMap;
250 for (
const auto& assnPair : *hitSpacePointAssnsHandle) {
252 spacePointToWireIDMap[assnPair.second.get()] = assnPair.first->WireID();
257 using OldHitToNewHitMap = std::map<const recob::Hit*, const recob::Hit*>;
258 using SpacePointHitVecMap = std::map<const recob::SpacePoint*, std::vector<const recob::Hit*>>;
260 OldHitToNewHitMap oldHitToNewHitMap;
261 SpacePointHitVecMap spacePointHitVecMap;
264 std::unique_ptr<std::vector<recob::Hit>> newHitVecPtr(
new std::vector<recob::Hit>);
267 newHitVecPtr->reserve(3 * hitSpacePointAssnsHandle->size());
270 art::PtrMaker<recob::Hit> ptrMaker(
evt);
272 for (
auto& assnPair : *hitSpacePointAssnsHandle) {
273 const art::Ptr<recob::SpacePoint> spacePoint = assnPair.second;
274 const art::Ptr<recob::Hit>& recobHit = assnPair.first;
277 if (oldHitToNewHitMap.find(recobHit.get()) == oldHitToNewHitMap.end()) {
279 geo::WireID refWireID = spacePointToWireIDMap[spacePoint.get()];
283 const std::vector<geo::WireID>& wireIDs =
287 for (
const auto& wireID : wireIDs) {
288 if (wireID.TPC != refWireID.
TPC || wireID.Cryostat != refWireID.
Cryostat)
continue;
299 spacePointHitVecMap[spacePoint.get()].push_back(newHit);
301 recobHitToArtPtrMap[newHit] = ptrMaker(newHitVecPtr->size() - 1);
302 oldHitToNewHitMap[recobHit.get()] = newHit;
305 spacePointHitVecMap[spacePoint.get()].push_back(oldHitToNewHitMap[recobHit.get()]);
310 std::map<geo::PlaneID, double> planeOffsetMap;
312 auto const clock_data =
313 art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(
evt);
314 auto const det_prop =
315 art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(
evt, clock_data);
319 for (
size_t tpcIdx = 0; tpcIdx <
fGeometry->
NTPC(); tpcIdx++) {
324 det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 1)) -
325 det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 0));
327 det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 2)) -
328 det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 0));
331 <<
", plane 1: " << planeOffsetMap[
geo::PlaneID(cryoIdx, tpcIdx, 1)]
332 <<
", plane 2: " << planeOffsetMap[
geo::PlaneID(cryoIdx, tpcIdx, 2)] << std::endl;
334 << det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 0))
335 <<
", plane 1: " << det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 1))
336 <<
", plane 2: " << det_prop.GetXTicksOffset(
geo::PlaneID(cryoIdx, tpcIdx, 2))
342 using RecobHitTo2DHitMap = std::map<const recob::Hit*, const reco::ClusterHit2D*>;
344 RecobHitTo2DHitMap recobHitTo2DHitMap;
351 for (
auto& hitPair : oldHitToNewHitMap) {
357 planeOffsetMap.at(hitWireID.planeID()));
358 double xPosition(det_prop.ConvertTicksToX(
359 recobHit->
PeakTime(), hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat));
367 for (
auto& pointPair : spacePointHitVecMap) {
369 const std::vector<const recob::Hit*>& recobHitVec = pointPair.second;
371 if (recobHitVec.size() != 3) {
372 std::cout <<
"************>>>>>> do not have 3 hits associated to space point! "
373 << recobHitVec.size() <<
" ***************" << std::endl;
379 for (
const auto& recobHit : recobHitVec) {
386 unsigned int statusBits(0x7);
387 float avePeakTime(0.);
394 for (
size_t planeIdx = 0; planeIdx < 3; planeIdx++) {
397 wireIDVec[planeIdx] = hit2D->
WireID();
405 float weight = 1. / (hitRMS * hitRMS);
408 avePeakTime += peakTime * weight;
412 avePeakTime /= weightSum;
415 float hitChiSquare(0.);
416 float sigmaPeakTime(std::sqrt(1. / weightSum));
417 std::vector<float> hitDelTSigVec;
419 for (
const auto& hit2D : hitVector) {
421 float combRMS = std::sqrt(hitRMS * hitRMS - sigmaPeakTime * sigmaPeakTime);
423 float deltaTime = peakTime - avePeakTime;
424 float hitSig = deltaTime / combRMS;
426 hitChiSquare += hitSig * hitSig;
428 hitDelTSigVec.emplace_back(std::fabs(hitSig));
434 int lowMinIndex(std::numeric_limits<int>::max());
435 int lowMaxIndex(std::numeric_limits<int>::min());
436 int hiMinIndex(std::numeric_limits<int>::max());
437 int hiMaxIndex(std::numeric_limits<int>::min());
440 for (
const auto& hit2D : hitVector) {
444 lowMinIndex = std::min(hitStart, lowMinIndex);
445 lowMaxIndex = std::max(hitStart, lowMaxIndex);
446 hiMinIndex = std::min(hitStop + 1, hiMinIndex);
447 hiMaxIndex = std::max(hitStop + 1, hiMaxIndex);
451 if (hitChiSquare < m_maxHit3DChiSquare && hiMinIndex > lowMaxIndex) {
453 std::vector<float> chargeVec;
455 for (
const auto& hit2D : hitVector)
464 std::accumulate(chargeVec.begin(), chargeVec.end(), 0.) / float(chargeVec.size());
465 float overlapRange = float(hiMinIndex - lowMaxIndex);
466 float overlapFraction = overlapRange / float(hiMaxIndex - lowMinIndex);
469 std::vector<float> smallestChargeDiffVec;
470 std::vector<float> chargeAveVec;
471 float smallestDiff(std::numeric_limits<float>::max());
472 size_t chargeIndex(0);
474 for (
size_t idx = 0; idx < 3; idx++) {
475 size_t leftIdx = (idx + 2) % 3;
476 size_t rightIdx = (idx + 1) % 3;
478 smallestChargeDiffVec.push_back(
std::abs(chargeVec[leftIdx] - chargeVec[rightIdx]));
479 chargeAveVec.push_back(
float(0.5 * (chargeVec[leftIdx] + chargeVec[rightIdx])));
481 if (smallestChargeDiffVec.back() < smallestDiff) {
482 smallestDiff = smallestChargeDiffVec.back();
488 float deltaPeakTime =
489 hitVector[leftIdx]->getTimeTicks() - hitVector[rightIdx]->getTimeTicks();
495 float chargeAsymmetry = (chargeAveVec[chargeIndex] - chargeVec[chargeIndex]) /
496 (chargeAveVec[chargeIndex] + chargeVec[chargeIndex]);
499 if (chargeAsymmetry < -1. || chargeAsymmetry > 1.) {
502 std::cout <<
"============> Charge asymmetry out of range: " << chargeAsymmetry
503 <<
" <============" << std::endl;
505 <<
", Plane: " << hitWireID.
Plane <<
", Wire: " << hitWireID.
Wire << std::endl;
506 std::cout <<
" charge: " << chargeVec[0] <<
", " << chargeVec[1] <<
", "
507 << chargeVec[2] << std::endl;
508 std::cout <<
" index: " << chargeIndex <<
", smallest diff: " << smallestDiff
514 float deltaPeakTime = *std::max_element(hitDelTSigVec.begin(), hitDelTSigVec.end());
527 Eigen::Vector3f position(
528 float(spacePoint->
XYZ()[0]),
float(spacePoint->
XYZ()[1]),
float(spacePoint->
XYZ()[2]));
531 hitPairList.emplace_back(0,
552 hitRefiner.use_hits(std::move(newHitVecPtr));
555 hitRefiner.put_into();
563 theClockMakeHits.stop();
568 mf::LogDebug(
"Cluster3D") <<
">>>>> 3D hit building done, found " << hitPairList.size()
569 <<
" 3D Hits" << std::endl;
bool m_outputHistograms
Take the time to create and fill some histograms for diagnostics.
art::InputTag fSpacePointProducerLabel
Data members to follow.
float getTimeTicks() const
float chargeIntegral(float, float, float, float, int, int) const
Perform charge integration between limits.
geo::WireID WireID() const
float RMS() const
RMS of the hit shape, in tick units.
std::vector< float > m_chiSquare3DVec
std::vector< float > m_overlapRangeVec
The data type to uniquely identify a Plane.
const geo::WireID & WireID() const
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
CryostatID_t Cryostat
Index of cryostat.
WireID_t Wire
Index of the wire within its plane.
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
std::vector< float > m_qualityMetricVec
const recob::Hit * getHit() const
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
void clear()
clear the tuple vectors before processing next event
Class managing the creation of a new recob::Hit object.
std::vector< float > m_deltaTimeVec
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
const geo::Geometry * fGeometry
Hit2DVector m_clusterHit2DMasterVec
A class handling a collection of hits and its associations.
const Double32_t * XYZ() const
std::vector< int > m_smallIndexVec
std::vector< float > m_maxPullVec
PlaneID_t Plane
Index of the plane within its TPC.
std::vector< float > m_smallChargeDiffVec
art::InputTag fHitProducerLabel
std::vector< const reco::ClusterHit2D * > ClusterHit2DVec
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
float PeakTime() const
Time of the signal peak, in tick units.
std::vector< float > m_hitAsymmetryVec
int trigger_offset(DetectorClocksData const &data)
std::vector< float > m_overlapFractionVec
constexpr WireID()=default
Default constructor: an invalid TPC ID.
2D representation of charge deposited in the TDC/wire plane
TPCID_t TPC
Index of the TPC within its cryostat.
TTree * m_tupleTree
output analysis tree
BEGIN_PROLOG could also be cout
std::vector< float > fTimeVector
std::vector< float > m_spacePointChargeVec
Signal from collection planes.
unsigned getStatusBits() const
void setStatusBit(unsigned bits) const