329 auto tracks = std::make_unique<std::vector<recob::Track>>();
330 auto allsp = std::make_unique<std::vector<recob::SpacePoint>>();
331 auto vtxs = std::make_unique<std::vector<recob::Vertex>>();
332 auto kinks = std::make_unique<
333 std::vector<recob::Vertex>>();
334 auto nodes = std::make_unique<std::vector<recob::Vertex>>();
335 auto t0s = std::make_unique<std::vector<anab::T0>>();
336 auto cosmicTags = std::make_unique<std::vector<anab::CosmicTag>>();
338 auto trk2hit_oldway = std::make_unique<
341 auto trk2hit = std::make_unique<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta>>();
343 auto trk2sp = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
344 auto trk2t0 = std::make_unique<art::Assns<recob::Track, anab::T0>>();
345 auto trk2ct = std::make_unique<art::Assns<recob::Track, anab::CosmicTag>>();
347 auto sp2hit = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
348 auto vtx2trk = std::make_unique<
352 std::make_unique<art::Assns<recob::Track, recob::Vertex>>();
354 auto pfps = std::make_unique<std::vector<recob::PFParticle>>();
356 auto pfp2clu = std::make_unique<art::Assns<recob::PFParticle, recob::Cluster>>();
357 auto pfp2vtx = std::make_unique<art::Assns<recob::PFParticle, recob::Vertex>>();
358 auto pfp2trk = std::make_unique<art::Assns<recob::PFParticle, recob::Track>>();
363 std::vector<art::Ptr<recob::Hit>> allhitlist;
377 size_t mvaLength = 0;
380 auto cluListHandle =
evt.getValidHandle<std::vector<recob::Cluster>>(
fCluModuleLabel);
381 auto splitCluHandle =
evt.getValidHandle<std::vector<recob::Cluster>>(
fEmModuleLabel);
385 pmalgTracker.init(hitsFromClusters, hitsFromEmParts);
390 if (init<4>(
evt, pmalgTracker)) { mvaLength = 4; }
391 else if (init<3>(
evt, pmalgTracker)) {
394 else if (init<2>(
evt, pmalgTracker)) {
398 throw cet::exception(
"PMAlgTrackMaker") <<
"No EM/track MVA data products." << std::endl;
403 auto cluListHandle =
evt.getValidHandle<std::vector<recob::Cluster>>(
fCluModuleLabel);
405 pmalgTracker.init(hitsFromClusters);
412 int retCode = pmalgTracker.build(clockData, detProp);
421 else if (retCode == 1)
429 auto const& result = pmalgTracker.result();
432 size_t spStart = 0, spEnd = 0;
433 double sp_pos[3], sp_err[6];
434 for (
size_t i = 0; i < 3; ++i)
436 for (
size_t i = 0; i < 6; ++i)
446 tracks->reserve(result.size());
447 for (
size_t trkIndex = 0; trkIndex < result.size(); ++trkIndex) {
476 if (trk->
size() < 2)
continue;
480 pdg = getPdgFromCnnOnHits<4>(
evt, *(result[trkIndex].Track()));
481 else if (mvaLength == 3)
482 pdg = getPdgFromCnnOnHits<3>(
evt, *(result[trkIndex].Track()));
483 else if (mvaLength == 2)
484 pdg = getPdgFromCnnOnHits<2>(
evt, *(result[trkIndex].Track()));
489 auto const trkPtr = make_trkptr(tracks->size() - 1);
493 t0s->push_back(
anab::T0(trk->
GetT0(), 0, 3, tracks->back().ID()));
495 auto const t0Ptr = make_t0ptr(t0s->size() - 1);
496 trk2t0->addSingle(trkPtr, t0Ptr);
502 std::vector<float> trkEnd0;
503 std::vector<float> trkEnd1;
506 int driftDir =
abs(tpc.DetectDriftDirection()) - 1;
508 for (
int i = 0; i < 3; ++i) {
512 if (i == driftDir) { shift = trk->
Nodes()[0]->GetDriftShift(); }
513 trkEnd0.push_back(trk->
Nodes()[0]->Point3D()[i] - shift);
514 trkEnd1.push_back(trk->
Nodes()[trk->
Nodes().size() - 1]->Point3D()[i] - shift);
519 for (
const auto t : tags) {
520 cosmicTags->emplace_back(trkEnd0, trkEnd1, 1, t);
521 auto const cosmicPtr = make_ctptr(cosmicTags->size() - 1);
522 trk2ct->addSingle(trkPtr, cosmicPtr);
527 for (
size_t h = 0, cnt = 0; h < trk->
size(); h++) {
532 trk2hit->addSingle(trkPtr, h3d->
Hit2DPtr(), metadata);
533 trk2hit_oldway->addSingle(
538 spStart = allsp->
size();
540 for (
size_t h = 0; h < trk->
size(); ++h) {
544 double hx = h3d->
Point3D().X();
545 double hy = h3d->
Point3D().Y();
546 double hz = h3d->
Point3D().Z();
548 if ((h == 0) || (std::fabs(sp_pos[0] - hx) > 1.0e-5) ||
549 (std::fabs(sp_pos[1] - hy) > 1.0e-5) || (std::fabs(sp_pos[2] - hz) > 1.0e-5)) {
567 spEnd = allsp->size();
569 if (spEnd > spStart)
util::CreateAssn(evt, *tracks, *allsp, *trk2sp, spStart, spEnd);
572 auto vsel = pmalgTracker.getVertices(
574 auto ksel = pmalgTracker.getKinks();
575 std::map<size_t, art::Ptr<recob::Vertex>> frontVtxs;
580 for (
auto const& v : vsel) {
581 xyz[0] = v.first.X();
582 xyz[1] = v.first.Y();
583 xyz[2] = v.first.Z();
584 mf::LogVerbatim(
"Summary") <<
" vtx:" << xyz[0] <<
":" << xyz[1] <<
":" << xyz[2]
585 <<
" (" << v.second.size() <<
" tracks)";
587 size_t vidx = vtxs->size();
588 vtxs->push_back(recob::Vertex(xyz, vidx));
590 auto const vptr = make_vtxptr(vidx);
591 if (!v.second.empty()) {
592 for (
const auto& vEntry : v.second) {
593 size_t tidx = vEntry.first;
594 bool isFront = vEntry.second;
596 if (isFront) frontVtxs[tidx] = vptr;
598 auto const tptr = make_trkptr(tidx);
599 vtx2trk->addSingle(vptr, tptr);
603 mf::LogWarning(
"PMAlgTrackMaker") <<
"No tracks found at this vertex.";
607 for (
auto const& k : ksel) {
608 xyz[0] = k.first.X();
609 xyz[1] = k.first.Y();
610 xyz[2] = k.first.Z();
611 mf::LogVerbatim(
"Summary") <<
" kink:" << xyz[0] <<
":" << xyz[1] <<
":" << xyz[2];
613 size_t kidx = kinks->size();
614 size_t tidx = k.second;
617 recob::Vertex(xyz, tidx));
619 auto const tptr = make_trkptr(tidx);
620 auto const kptr = make_kinkptr(kidx);
621 trk2kink->addSingle(tptr, kptr);
628 for (
size_t t = 0; t < result.size(); ++t) {
629 auto const& trk = *(result[t].Track());
630 for (
auto const* node : trk.
Nodes()) {
631 xyz[0] = node->Point3D().X();
632 xyz[1] = node->Point3D().Y();
633 xyz[2] = node->Point3D().Z();
639 for (
size_t t = 0; t < result.size(); ++t) {
641 if (result[t].Parent() >= 0) parentIdx = (
size_t)result[t].Parent();
643 std::vector<size_t> daughterIdxs;
644 for (
size_t idx : result[t].Daughters()) {
645 daughterIdxs.push_back(idx);
648 size_t pfpidx = pfps->size();
649 pfps->emplace_back((*tracks)[t].ParticleId(), pfpidx, parentIdx, daughterIdxs);
651 auto const pfpptr = make_pfpptr(pfpidx);
652 auto const tptr = make_trkptr(t);
653 pfp2trk->addSingle(pfpptr, tptr);
659 pfp2vtx->addSingle(pfpptr, vptr);
661 mf::LogWarning(
"PMAlgTrackMaker") <<
"Front vertex for PFParticle is missing.";
665 <<
" PFParticles created for reconstructed tracks.";
666 mf::LogVerbatim(
"Summary") <<
"Adding " << result.parents().size() <<
" primary PFParticles.";
667 for (
size_t t = 0; t < result.parents().size(); ++t) {
668 std::vector<size_t> daughterIdxs;
669 for (
size_t idx : result.parents()[t].Daughters()) {
670 daughterIdxs.push_back(idx);
673 size_t pfpidx = pfps->size();
675 pfps->emplace_back(0, pfpidx, parentIdx, daughterIdxs);
679 auto const pfpptr = make_pfpptr(pfpidx);
681 frontVtxs[daughterIdxs.front()];
683 pfp2vtx->addSingle(pfpptr, vptr);
685 mf::LogWarning(
"PMAlgTrackMaker") <<
"Front vertex for PFParticle is missing.";
688 mf::LogVerbatim(
"Summary") << pfps->size() <<
" PFParticles created in total.";
693 evt.put(std::move(tracks));
694 evt.put(std::move(allsp));
695 evt.put(std::move(vtxs));
698 evt.put(std::move(t0s));
699 evt.put(std::move(cosmicTags));
701 evt.put(std::move(trk2hit_oldway));
702 evt.put(std::move(trk2hit));
703 evt.put(std::move(trk2sp));
704 evt.put(std::move(trk2t0));
705 evt.put(std::move(trk2ct));
707 evt.put(std::move(sp2hit));
708 evt.put(std::move(vtx2trk));
711 evt.put(std::move(pfps));
712 evt.put(std::move(pfp2clu));
713 evt.put(std::move(pfp2vtx));
714 evt.put(std::move(pfp2trk));
bool SelectHits(float fraction=1.0F)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
static const std::string kKinksName
recob::Track convertFrom(const pma::Track3D &src, unsigned int tidx, int pdg=0)
static constexpr size_t kPFParticlePrimary
Define index to signify primary particle.
bool IsEnabled() const noexcept
Planes which measure X direction.
constexpr auto abs(T v)
Returns the absolute value of the argument.
TVector3 const & Point3D() const
Planes which measure Z direction.
fhicl::Atom< bool > RunVertexing
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Planes which measure Y direction.
std::vector< anab::CosmicTagID_t > getCosmicTag(const pma::Track3D::ETag pmaTag) const
geo::GeometryCore const * fGeom
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
pma::PMAlgStitching::Config fPmaStitchConfig
pma::ProjectionMatchingAlg::Config fPmaConfig
void push_back(Ptr< U > const &p)
bool HasTagFlag(ETag value) const noexcept
pma::PMAlgVertexing::Config fPmaVtxConfig
art::InputTag fCluModuleLabel
bool isNull() const noexcept
unsigned int FrontTPC() const
unsigned int FrontCryo() const
The data type to uniquely identify a TPC.
pma::PMAlgCosmicTagger::Config fPmaTaggingConfig
ETag GetTag() const noexcept
art::InputTag fHitModuleLabel
size_t CompleteMissingWires(detinfo::DetectorPropertiesData const &detProp, unsigned int view)
void push_back(pma::Hit3D *hit)
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
std::map< size_t, std::vector< double > > dedx_map
fhicl::Atom< float > TrackLikeThreshold
double GetRawdEdxSequence(std::map< size_t, std::vector< double >> &dedx, unsigned int view=geo::kZ, unsigned int skip=0, bool inclDisabled=false) const
double Dx() const noexcept
bool HasT0() const noexcept
art::InputTag fEmModuleLabel
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
art::InputTag fWireModuleLabel
art::Ptr< recob::Hit > const & Hit2DPtr() const
std::vector< TH1F * > fAdcInRejectedPoints
bool fSaveOnlyBranchingVtx
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
std::vector< pma::Node3D * > const & Nodes() const noexcept
pma::PMAlgTracker::Config fPmaTrackerConfig
2D representation of charge deposited in the TDC/wire plane
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
std::vector< TH1F * > fAdcInPassingPoints
cet::coded_exception< error, detail::translate > exception
static const std::string kNodesName