377 std::unique_ptr< std::vector<anab::TPCPurityInfo> > outputPtrVector(
new std::vector<anab::TPCPurityInfo>());
382 auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
389 art::Handle< std::vector<recob::Track>> trackHandle;
390 event.getByLabel(trackLabel, trackHandle);
392 if (!trackHandle.isValid())
continue;
397 using HitToSpacePointMap = std::unordered_map<const recob::Hit*,const recob::SpacePoint*>;
399 HitToSpacePointMap hitToSpacePointMap;
401 using PointCloud = std::vector<geo::Point_t>;
403 PointCloud pointCloud;
406 art::FindManyP<recob::Hit,recob::TrackHitMeta> trackHitAssns(trackHandle, event, trackLabel);
409 for(
size_t trackIdx = 0; trackIdx < trackHandle->size(); trackIdx++)
411 art::Ptr<recob::Track>
track(trackHandle,trackIdx);
413 const std::vector<art::Ptr<recob::Hit>>& trackHitsVec(trackHitAssns.at(
track.key()));
414 const std::vector<const recob::TrackHitMeta*> metaHitsVec = trackHitAssns.data(
track.key());
419 using TPCToHitMetaPairVecMap = std::unordered_map<unsigned int,HitMetaPairVec>;
421 TPCToHitMetaPairVecMap selectedHitMetaVecMap;
423 for(
size_t idx=0; idx<trackHitsVec.size(); idx++)
425 art::Ptr<recob::Hit>
hit(trackHitsVec.at(idx));
427 if (
hit->WireID().Plane ==
fSelectedPlane &&
hit->Multiplicity() == 1) selectedHitMetaVecMap[
hit->WireID().TPC].emplace_back(
hit,metaHitsVec.at(idx));
430 if (selectedHitMetaVecMap.empty())
continue;
434 TPCToHitMetaPairVecMap::iterator bestMapItr = selectedHitMetaVecMap.begin();
436 for(TPCToHitMetaPairVecMap::iterator mapItr = selectedHitMetaVecMap.begin(); mapItr != selectedHitMetaVecMap.end(); mapItr++)
438 if (mapItr->second.size() > bestMapItr->second.size()) bestMapItr = mapItr;
444 if (selectedHitMetaVec.size() <
fMinNumHits)
continue;
447 std::sort(selectedHitMetaVec.begin(),selectedHitMetaVec.end(),[](
const auto&
left,
const auto&
right){
return left.first->PeakTime() <
right.first->PeakTime();});
450 if (selectedHitMetaVec.back().first->PeakTime() - selectedHitMetaVec.front().first->PeakTime() <
fMinTickRange)
continue;
458 float firstHitTime(selectedHitMetaVec.front().first->PeakTime());
459 double maxDeltaX(1.5);
460 double wirePitch(0.3);
462 for(
const auto& hitMetaPair: selectedHitMetaVec)
464 unsigned int trkHitIndex = hitMetaPair.second->Index();
466 double cosTheta = -100.;
468 if (trkHitIndex != std::numeric_limits<unsigned int>::max() &&
track->HasValidPoint(trkHitIndex))
475 pointCloud.emplace_back(hitPos);
477 cosTheta =
std::abs(hitDir.Dot(wireDir));
481 deltaX = std::min(wirePitch / (1. - cosTheta), maxDeltaX);
483 else deltaX = maxDeltaX;
485 double charge =
fUseHitIntegral ? hitMetaPair.first->Integral() : hitMetaPair.first->SummedADC();
487 hitStatusChargePairVec.emplace_back(hitMetaPair,
StatusChargePair(
true,charge/deltaX));
488 hitPointDirTupleMap[hitMetaPair.first.get()] =
PointDirTuple(hitPos,hitDir,wireDir);
492 size_t numOrig = hitStatusChargePairVec.size();
497 if (lowCutIdx + 10 >= hiCutIdx)
499 std::cout <<
"*****>>>> lowCutIdx >= hiCutIdx: " << lowCutIdx <<
", " << hiCutIdx << std::endl;
512 PrincipalComponents2D pca;
523 if (!pca.getSvdOK())
continue;
526 PrincipalComponents3D pca3D;
532 double attenuation = eigenVectors.row(1)[1] / eigenVectors.row(1)[0];
536 unsigned minWire(100000);
538 std::set<unsigned> usedWiresSet;
540 for(
const auto& hitPair : hitStatusChargePairVec)
542 unsigned wire = hitPair.first.first->WireID().Wire;
544 usedWiresSet.insert(wire);
546 if (wire > maxWire) maxWire = wire;
547 if (wire < minWire) minWire = wire;
554 purityInfo.
Run =
event.run();
555 purityInfo.
Subrun =
event.subRun();
556 purityInfo.
Event =
event.event();
558 purityInfo.
TPC = wireID.
TPC;
559 purityInfo.
Wires = usedWiresSet.size();
560 purityInfo.
Ticks = hitStatusChargePairVec.back().first.first->PeakTime() - firstHitTime;
562 purityInfo.
FracError = std::sqrt(pca.getEigenValues()[0] / pca.getEigenValues()[1]);
564 outputPtrVector->emplace_back(purityInfo);
575 fWires = usedWiresSet.size();
576 fTicks = hitStatusChargePairVec.back().first.first->PeakTime() - firstHitTime;
578 fError = std::sqrt(pca.getEigenValues()[0] / pca.getEigenValues()[1]);
581 std::sort(hitStatusChargePairVec.begin(),hitStatusChargePairVec.end(),[](
const auto&
left,
const auto&
right){
return left.first.second->Index() <
right.first.second->Index();});
583 const geo::Point_t& trackStartPos =
track->LocationAtPoint(hitStatusChargePairVec.front().first.second->Index());
584 const geo::Vector_t trackStartDir =
track->DirectionAtPoint(hitStatusChargePairVec.front().first.second->Index());
593 const geo::Point_t& trackEndPos =
track->LocationAtPoint(hitStatusChargePairVec.back().first.second->Index());
594 const geo::Vector_t trackEndDir =
track->DirectionAtPoint(hitStatusChargePairVec.back().first.second->Index());
604 for(
size_t rowIdx = 0; rowIdx < 2; rowIdx++)
606 for(
size_t colIdx = 0; colIdx < 2; colIdx++)
fPCAAxes2D.emplace_back(eigenVectors.row(rowIdx)[colIdx]);
618 for(
size_t rowIdx = 0; rowIdx < 3; rowIdx++)
620 for(
size_t colIdx = 0; colIdx < 3; colIdx++)
fPCAAxes3D.emplace_back(eigenVectors3D.row(rowIdx)[colIdx]);
631 for(
const auto& hitPair : hitStatusChargePairVec)
633 fTickVec.emplace_back(hitPair.first.first->PeakTime());
634 fChargeVec.emplace_back(hitPair.second.second);
635 fDeltaXVec.emplace_back(hitPair.first.second->Dx());
638 fSnippetLengthVec.emplace_back(hitPair.first.first->EndTick() - hitPair.first.first->StartTick());
642 const geo::Vector_t& hitDir = std::get<1>(hitPointDirTupleMap[hitPair.first.first.get()]);
643 const geo::Vector_t& wireDir = std::get<2>(hitPointDirTupleMap[hitPair.first.first.get()]);
648 hitDirYZ /= std::sqrt(hitDirYZ.Mag2());
686 event.put(std::move(outputPtrVector));
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
int fWires
Number wires spanned.
std::vector< art::InputTag > fTrackLabelVec
labels for source of tracks
int fWireRange
Last - first wire number.
Eigen::Matrix3d EigenVectors
std::vector< double > fPCAAxes2D
Axes for PCA.
int fTPC
TPC for given track.
TTree * fDiagnosticTree
Pointer to our tree.
std::vector< double > fTrackEndDirYVec
Ending x direction of track.
std::vector< double > fTrackDirZVec
Starting x direction of track.
std::vector< double > fMeanPosition3D
Mean position used for PCA.
unsigned fMinNumHits
Minimum number of hits.
int fEventNumber
event number this event
std::vector< bool > fGoodHitVec
Hits were considered good.
std::vector< double > fGoodnessOfFitVec
Goodness of the hit's fit.
CryostatID_t Cryostat
Index of cryostat.
std::unordered_map< const recob::Hit *, PointDirTuple > HitPointDirTupleMap
std::vector< double > fEigenValues2D
Eigen values.
process_name use argoneut_mc_hitfinder track
std::tuple< geo::Point_t, geo::Vector_t, geo::Vector_t > PointDirTuple
std::vector< double > fTickVec
vector of ticks
double fAttenuation
Attenuation from calc.
std::pair< HitMetaPair, StatusChargePair > HitStatusChargePair
std::vector< double > fTrackDirXVec
Starting x direction of track.
void RejectOutliers(HitStatusChargePairVec &hitPairVector, const PrincipalComponents2D &pca) const
void GetPrincipalComponents2D(const HitStatusChargePairVec &hitPairVector, PrincipalComponents2D &pca) const
std::vector< double > fEigenValues3D
Eigen values 3D.
std::vector< double > fTrackEndXVec
Ending x position of track.
int fRunNumber
run number this event
std::vector< double > fTrackEndDirZVec
Ending x direction of track.
std::vector< double > fPCAAxes3D
Axes for PCA 3D.
std::vector< double > fDeltaXVec
Keep track of hits path length from track fit.
std::vector< double > fChargeVec
vector of hit charges
std::vector< int > fSnippetLengthVec
Lenght from start/end of hit.
std::vector< double > fTrackEndZVec
Ending z position of track.
bool fUseHitIntegral
Setting to false swaps to "SummedADC".
std::vector< double > fTrackDirYVec
Starting x direction of track.
int fTicks
Number ticks spanned.
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
std::vector< double > fTrackEndDirXVec
Ending x direction of track.
std::vector< HitMetaPair > HitMetaPairVec
std::pair< bool, double > StatusChargePair
std::vector< double > fCosThetaYZ
cos(thetaYZ) hit trajector to wire
Eigen::Matrix2d EigenVectors
void GetPrincipalComponents3D(const HitStatusChargePairVec &hitPairVector, HitPointDirTupleMap &, PrincipalComponents3D &pca) const
bool fDiagnosticTuple
Output the diagnostic tuple.
int fSubRunNumber
sub run number this event
float fMinRejectFraction
Reject this fraction of "low Q" hits.
int fTrackIdx
index of track
int fCryostat
Cryostat for given track.
unsigned fSelectedPlane
Select hits from this plane.
constexpr WireID()=default
Default constructor: an invalid TPC ID.
std::vector< HitStatusChargePair > HitStatusChargePairVec
std::vector< double > fTrackStartZVec
Starting z position of track.
std::vector< double > fTrackStartYVec
Starting y position of track.
TPCID_t TPC
Index of the TPC within its cryostat.
std::vector< double > fTrackStartXVec
Starting x position of track.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
float fMinTickRange
Require track to span some number of ticks.
float fMaxRejectFraction
Reject this fraction of "high Q" hits.
std::vector< int > fDegreesOfFreeVec
Degrees of freedom.
double fError
Error from calc.
std::vector< double > fMeanPosition2D
Mean position used for PCA.
BEGIN_PROLOG could also be cout
std::vector< double > fTrackEndYVec
Ending y position of track.
const geo::GeometryCore * fGeometry
float fSamplingRate
Recover the sampling rate from the clock data.