38 #include "art_root_io/TFileService.h" 79 Name(
"ProjectionMatchingAlg")};
90 Name(
"SaveOnlyBranchingVtx"),
91 Comment(
"use true to save only vertices interconnecting many tracks, otherwise vertex is " 92 "added to the front of each track")};
96 Comment(
"save track nodes (only for algorithm development purposes)")};
99 Name(
"HitModuleLabel"),
100 Comment(
"tag of unclustered hits, which were used to validate tracks")};
103 Comment(
"tag of recob::Wire producer.")};
106 Name(
"ClusterModuleLabel"),
107 Comment(
"tag of cluster collection, these clusters are used for track building")};
110 Name(
"EmClusterModuleLabel"),
111 Comment(
"EM-like clusters, will be excluded from tracking if provided")};
174 ,
fPmaConfig(config().ProjectionMatchingAlg())
183 produces<std::vector<recob::Track>>();
184 produces<std::vector<recob::SpacePoint>>();
185 produces<std::vector<recob::Vertex>>();
186 produces<std::vector<recob::Vertex>>(
kKinksName);
187 produces<std::vector<recob::Vertex>>(
kNodesName);
188 produces<std::vector<anab::T0>>();
189 produces<std::vector<anab::CosmicTag>>();
193 produces<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta>>();
195 produces<art::Assns<recob::Track, recob::SpacePoint>>();
196 produces<art::Assns<recob::SpacePoint, recob::Hit>>();
200 produces<art::Assns<recob::Track, recob::Vertex>>(
kKinksName);
201 produces<art::Assns<recob::Track, anab::T0>>();
202 produces<art::Assns<recob::Track, anab::CosmicTag>>();
204 produces<std::vector<recob::PFParticle>>();
205 produces<art::Assns<recob::PFParticle, recob::Cluster>>();
206 produces<art::Assns<recob::PFParticle, recob::Vertex>>();
207 produces<art::Assns<recob::PFParticle, recob::Track>>();
213 for (
size_t p = 0; p < wireReadoutGeom.MaxPlanes(); ++p) {
214 std::ostringstream ss1;
215 ss1 <<
"adc_plane_" << p;
217 (ss1.str() +
"_passing").c_str(),
"max adc around the point on track", 100., 0., 5.));
219 (ss1.str() +
"_rejected").c_str(),
"max adc around spurious point ", 100., 0., 5.));
232 int trkLikeIdx = hitResults->getIndex(
"track");
233 int emLikeIdx = hitResults->getIndex(
"em");
234 if ((trkLikeIdx < 0) || (emLikeIdx < 0)) {
236 <<
"No em/track labeled columns in MVA data products." << std::endl;
239 size_t nh[3] = {0, 0, 0};
240 for (
size_t hidx = 0; hidx < trk.
size(); ++hidx) {
241 ++nh[trk[hidx]->View2D()];
244 size_t best_view = 2;
245 if ((nh[0] >= nh[1]) && (nh[0] > 1.25 * nh[2])) best_view = 0;
246 if ((nh[1] >= nh[0]) && (nh[1] > 1.25 * nh[2])) best_view = 1;
248 std::vector<art::Ptr<recob::Hit>> trkHitPtrList;
249 trkHitPtrList.reserve(nh[best_view]);
250 for (
size_t hidx = 0; hidx < trk.
size(); ++hidx) {
251 if (trk[hidx]->View2D() == best_view) {
252 trkHitPtrList.emplace_back(trk[hidx]->Hit2DPtr());
255 auto vout = hitResults->getOutput(trkHitPtrList);
256 double trk_like = -1, trk_or_em = vout[trkLikeIdx] + vout[emLikeIdx];
258 trk_like = vout[trkLikeIdx] / trk_or_em;
272 if (!cluResults) {
return false; }
274 int trkLikeIdx = cluResults->getIndex(
"track");
275 int emLikeIdx = cluResults->getIndex(
"em");
276 if ((trkLikeIdx < 0) || (emLikeIdx < 0)) {
return false; }
279 cluResults->dataHandle(),
evt, cluResults->dataTag());
280 const auto& cnnOuts = cluResults->outputs();
281 std::vector<float> trackLike(cnnOuts.size());
282 for (
size_t i = 0; i < cnnOuts.size(); ++i) {
283 double trkOrEm = cnnOuts[i][trkLikeIdx] + cnnOuts[i][emLikeIdx];
284 if (trkOrEm > 0) { trackLike[i] = cnnOuts[i][trkLikeIdx] / trkOrEm; }
289 pmalgTracker.
init(hitsFromClusters, trackLike);
296 std::vector<anab::CosmicTagID_t> anabTags;
331 auto tracks = std::make_unique<std::vector<recob::Track>>();
332 auto allsp = std::make_unique<std::vector<recob::SpacePoint>>();
333 auto vtxs = std::make_unique<std::vector<recob::Vertex>>();
334 auto kinks = std::make_unique<
335 std::vector<recob::Vertex>>();
336 auto nodes = std::make_unique<std::vector<recob::Vertex>>();
337 auto t0s = std::make_unique<std::vector<anab::T0>>();
338 auto cosmicTags = std::make_unique<std::vector<anab::CosmicTag>>();
340 auto trk2hit_oldway = std::make_unique<
343 auto trk2hit = std::make_unique<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta>>();
345 auto trk2sp = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
346 auto trk2t0 = std::make_unique<art::Assns<recob::Track, anab::T0>>();
347 auto trk2ct = std::make_unique<art::Assns<recob::Track, anab::CosmicTag>>();
349 auto sp2hit = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
350 auto vtx2trk = std::make_unique<
354 std::make_unique<art::Assns<recob::Track, recob::Vertex>>();
356 auto pfps = std::make_unique<std::vector<recob::PFParticle>>();
358 auto pfp2clu = std::make_unique<art::Assns<recob::PFParticle, recob::Cluster>>();
359 auto pfp2vtx = std::make_unique<art::Assns<recob::PFParticle, recob::Vertex>>();
360 auto pfp2trk = std::make_unique<art::Assns<recob::PFParticle, recob::Track>>();
365 std::vector<art::Ptr<recob::Hit>> allhitlist;
379 size_t mvaLength = 0;
387 pmalgTracker.init(hitsFromClusters, hitsFromEmParts);
392 if (init<4>(evt, pmalgTracker)) { mvaLength = 4; }
393 else if (init<3>(evt, pmalgTracker)) {
396 else if (init<2>(evt, pmalgTracker)) {
400 throw cet::exception(
"PMAlgTrackMaker") <<
"No EM/track MVA data products." << std::endl;
407 pmalgTracker.init(hitsFromClusters);
414 int retCode = pmalgTracker.build(clockData, detProp);
423 else if (retCode == 1)
431 auto const& result = pmalgTracker.result();
434 size_t spStart = 0, spEnd = 0;
435 double sp_pos[3], sp_err[6];
436 for (
size_t i = 0; i < 3; ++i)
438 for (
size_t i = 0; i < 6; ++i)
448 tracks->reserve(result.size());
451 for (
size_t trkIndex = 0; trkIndex < result.size(); ++trkIndex) {
458 auto const& tpcid = tpcgeo.
ID();
459 auto has_plane = [&wireReadoutGeom, &tpcid](
auto const view) {
460 return wireReadoutGeom.HasPlane(
geo::PlaneID(tpcid, view));
484 if (trk->
size() < 2)
continue;
488 pdg = getPdgFromCnnOnHits<4>(
evt, *(result[trkIndex].Track()));
489 else if (mvaLength == 3)
490 pdg = getPdgFromCnnOnHits<3>(
evt, *(result[trkIndex].Track()));
491 else if (mvaLength == 2)
492 pdg = getPdgFromCnnOnHits<2>(
evt, *(result[trkIndex].Track()));
497 auto const trkPtr = make_trkptr(tracks->size() - 1);
501 t0s->push_back(
anab::T0(trk->
GetT0(), 0, 3, tracks->back().ID()));
503 auto const t0Ptr = make_t0ptr(t0s->size() - 1);
504 trk2t0->addSingle(trkPtr, t0Ptr);
510 std::vector<float> trkEnd0;
511 std::vector<float> trkEnd1;
516 for (
int i = 0; i < 3; ++i) {
520 if (i ==
to_int(axis)) { shift = trk->
Nodes()[0]->GetDriftShift(); }
521 trkEnd0.push_back(trk->
Nodes()[0]->Point3D()[i] - shift);
522 trkEnd1.push_back(trk->
Nodes()[trk->
Nodes().size() - 1]->Point3D()[i] - shift);
527 for (
const auto t : tags) {
528 cosmicTags->emplace_back(trkEnd0, trkEnd1, 1, t);
529 auto const cosmicPtr = make_ctptr(cosmicTags->size() - 1);
530 trk2ct->addSingle(trkPtr, cosmicPtr);
535 for (
size_t h = 0, cnt = 0; h < trk->
size(); h++) {
540 trk2hit->addSingle(trkPtr, h3d->
Hit2DPtr(), metadata);
541 trk2hit_oldway->addSingle(
546 spStart = allsp->
size();
548 for (
size_t h = 0; h < trk->
size(); ++h) {
552 double hx = h3d->
Point3D().X();
553 double hy = h3d->
Point3D().Y();
554 double hz = h3d->
Point3D().Z();
556 if ((h == 0) || (std::fabs(sp_pos[0] - hx) > 1.0
e-5) ||
557 (std::fabs(sp_pos[1] - hy) > 1.0
e-5) || (std::fabs(sp_pos[2] - hz) > 1.0
e-5)) {
575 spEnd = allsp->size();
577 if (spEnd > spStart)
util::CreateAssn(evt, *tracks, *allsp, *trk2sp, spStart, spEnd);
580 auto vsel = pmalgTracker.getVertices(
582 auto ksel = pmalgTracker.getKinks();
583 std::map<size_t, art::Ptr<recob::Vertex>> frontVtxs;
588 for (
auto const& v : vsel) {
589 xyz[0] = v.first.X();
590 xyz[1] = v.first.Y();
591 xyz[2] = v.first.Z();
592 mf::LogVerbatim(
"Summary") <<
" vtx:" << xyz[0] <<
":" << xyz[1] <<
":" << xyz[2]
593 <<
" (" << v.second.size() <<
" tracks)";
595 size_t vidx = vtxs->size();
596 vtxs->push_back(recob::Vertex(xyz, vidx));
598 auto const vptr = make_vtxptr(vidx);
599 if (!v.second.empty()) {
600 for (
const auto& vEntry : v.second) {
601 size_t tidx = vEntry.first;
602 bool isFront = vEntry.second;
604 if (isFront) frontVtxs[tidx] = vptr;
606 auto const tptr = make_trkptr(tidx);
607 vtx2trk->addSingle(vptr, tptr);
611 mf::LogWarning(
"PMAlgTrackMaker") <<
"No tracks found at this vertex.";
615 for (
auto const& k : ksel) {
616 xyz[0] = k.first.X();
617 xyz[1] = k.first.Y();
618 xyz[2] = k.first.Z();
619 mf::LogVerbatim(
"Summary") <<
" kink:" << xyz[0] <<
":" << xyz[1] <<
":" << xyz[2];
621 size_t kidx = kinks->size();
622 size_t tidx = k.second;
625 recob::Vertex(xyz, tidx));
627 auto const tptr = make_trkptr(tidx);
628 auto const kptr = make_kinkptr(kidx);
629 trk2kink->addSingle(tptr, kptr);
636 for (
size_t t = 0; t < result.size(); ++t) {
637 auto const&
trk = *(result[t].Track());
638 for (
auto const* node :
trk.Nodes()) {
639 xyz[0] = node->Point3D().X();
640 xyz[1] = node->Point3D().Y();
641 xyz[2] = node->Point3D().Z();
642 nodes->push_back(recob::Vertex(xyz, t));
647 for (
size_t t = 0; t < result.size(); ++t) {
649 if (result[t].Parent() >= 0) parentIdx = (
size_t)result[t].Parent();
651 std::vector<size_t> daughterIdxs;
652 for (
size_t idx : result[t].Daughters()) {
653 daughterIdxs.push_back(idx);
656 size_t pfpidx = pfps->size();
657 pfps->emplace_back((*tracks)[t].ParticleId(), pfpidx, parentIdx, daughterIdxs);
659 auto const pfpptr = make_pfpptr(pfpidx);
660 auto const tptr = make_trkptr(t);
661 pfp2trk->addSingle(pfpptr, tptr);
667 pfp2vtx->addSingle(pfpptr, vptr);
669 mf::LogWarning(
"PMAlgTrackMaker") <<
"Front vertex for PFParticle is missing.";
673 <<
" PFParticles created for reconstructed tracks.";
674 mf::LogVerbatim(
"Summary") <<
"Adding " << result.parents().size() <<
" primary PFParticles.";
675 for (
size_t t = 0; t < result.parents().size(); ++t) {
676 std::vector<size_t> daughterIdxs;
677 for (
size_t idx : result.parents()[t].Daughters()) {
678 daughterIdxs.push_back(idx);
681 size_t pfpidx = pfps->size();
683 pfps->emplace_back(0, pfpidx, parentIdx, daughterIdxs);
687 auto const pfpptr = make_pfpptr(pfpidx);
689 frontVtxs[daughterIdxs.front()];
691 pfp2vtx->addSingle(pfpptr, vptr);
693 mf::LogWarning(
"PMAlgTrackMaker") <<
"Front vertex for PFParticle is missing.";
696 mf::LogVerbatim(
"Summary") << pfps->size() <<
" PFParticles created in total.";
701 evt.
put(std::move(tracks));
702 evt.
put(std::move(allsp));
703 evt.
put(std::move(vtxs));
706 evt.
put(std::move(t0s));
707 evt.
put(std::move(cosmicTags));
709 evt.
put(std::move(trk2hit_oldway));
710 evt.
put(std::move(trk2hit));
711 evt.
put(std::move(trk2sp));
712 evt.
put(std::move(trk2t0));
713 evt.
put(std::move(trk2ct));
715 evt.
put(std::move(sp2hit));
716 evt.
put(std::move(vtx2trk));
719 evt.
put(std::move(pfps));
720 evt.
put(std::move(pfp2clu));
721 evt.
put(std::move(pfp2vtx));
722 evt.
put(std::move(pfp2trk));
PMAlgTrackMaker & operator=(PMAlgTrackMaker const &)=delete
bool SelectHits(float fraction=1.0F)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
static const std::string kKinksName
void produce(art::Event &e) override
recob::Track convertFrom(const pma::Track3D &src, unsigned int tidx, int pdg=0)
fhicl::Atom< art::InputTag > EmClusterModuleLabel
static constexpr size_t kPFParticlePrimary
Define index to signify primary particle.
Declaration of signal hit object.
fhicl::Atom< art::InputTag > WireModuleLabel
bool IsEnabled() const noexcept
fhicl::Atom< art::InputTag > HitModuleLabel
fhicl::Table< pma::PMAlgVertexing::Config > PMAlgVertexing
EDProducer(fhicl::ParameterSet const &pset)
Planes which measure X direction.
The data type to uniquely identify a Plane.
fhicl::Table< pma::PMAlgCosmicTagger::Config > PMAlgCosmicTagging
Geometry information for a single TPC.
fhicl::Atom< bool > SavePmaNodes
fhicl::Table< pma::ProjectionMatchingAlg::Config > ProjectionMatchingAlg
void produces(std::string const &instanceName={}, Persistable const persistable=Persistable::Yes)
fhicl::Atom< std::string > Validation
fhicl::Atom< art::InputTag > ClusterModuleLabel
TVector3 const & Point3D() const
Planes which measure Z direction.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
fhicl::Atom< bool > RunVertexing
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DriftAxis DriftAxisWithSign() const
Returns the expected drift direction based on geometry.
Planes which measure Y direction.
std::vector< anab::CosmicTagID_t > getCosmicTag(const pma::Track3D::ETag pmaTag) const
geo::GeometryCore const * fGeom
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
pma::PMAlgStitching::Config fPmaStitchConfig
#define DEFINE_ART_MODULE(klass)
pma::ProjectionMatchingAlg::Config fPmaConfig
void push_back(Ptr< U > const &p)
bool HasTagFlag(ETag value) const noexcept
pma::PMAlgVertexing::Config fPmaVtxConfig
int getPdgFromCnnOnHits(const art::Event &evt, const pma::Track3D &trk) const
art::InputTag fCluModuleLabel
bool isNull() const noexcept
constexpr int to_int(Coordinate const coord) noexcept
Enumerate the possible plane projections.
unsigned int FrontTPC() const
unsigned int FrontCryo() const
pma::PMAlgCosmicTagger::Config fPmaTaggingConfig
Description of the physical geometry of one entire detector.
Declaration of cluster object.
ETag GetTag() const noexcept
art::InputTag fHitModuleLabel
size_t CompleteMissingWires(detinfo::DetectorPropertiesData const &detProp, unsigned int view)
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
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
double GetRawdEdxSequence(std::map< size_t, std::vector< double >> &dedx, unsigned int view=geo::kZ, unsigned int skip=0, bool inclDisabled=false) const
Utility object to perform functions of association.
Provides recob::Track data product.
void init(const art::FindManyP< recob::Hit > &hitsFromClusters)
double Dx() const noexcept
bool HasT0() const noexcept
fhicl::Table< pma::PMAlgTracker::Config > PMAlgTracking
bool init(const art::Event &evt, pma::PMAlgTracker &pmalgTracker) const
art::InputTag fEmModuleLabel
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
fhicl::Table< pma::PMAlgStitching::Config > PMAlgStitching
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)
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
TPCID const & ID() const
Returns the identifier of this TPC.
PMAlgTrackMaker(Parameters const &config)
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
art framework interface to geometry description
std::vector< TH1F * > fAdcInPassingPoints
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane .
fhicl::Atom< bool > SaveOnlyBranchingVtx
static const std::string kNodesName