159 auto tracks = std::make_unique< std::vector<recob::Track> >();
160 auto allsp = std::make_unique< std::vector<recob::SpacePoint> >();
161 auto vtxs = std::make_unique< std::vector<recob::Vertex> >();
162 auto kinks = std::make_unique< std::vector<recob::Vertex> >();
163 auto nodes = std::make_unique< std::vector<recob::Vertex> >();
165 auto trk2hit_oldway = std::make_unique< art::Assns<recob::Track, recob::Hit> >();
166 auto trk2hit = std::make_unique< art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> >();
168 auto trk2sp = std::make_unique< art::Assns<recob::Track, recob::SpacePoint> >();
170 auto sp2hit = std::make_unique< art::Assns<recob::SpacePoint, recob::Hit> >();
171 auto vtx2trk = std::make_unique< art::Assns<recob::Vertex, recob::Track> >();
172 auto trk2kink = std::make_unique< art::Assns<recob::Track, recob::Vertex> >();
174 auto pfp2trk = std::make_unique< art::Assns< recob::PFParticle, recob::Track> >();
180 std::vector< art::Ptr<recob::Hit> > allhitlist;
185 throw cet::exception(
"PMAlgTrajFitter") <<
"Not all required data products found in the event." << std::endl;
196 *cluListHandle, *pfparticleHandle,
197 hitsFromClusters, clustersFromPfps, vtxFromPfps,
201 int retCode = pmalgFitter.build();
210 else if (retCode == 1)
mf::LogVerbatim(
"Summary") << retCode <<
" track ready";
216 auto const & result = pmalgFitter.result();
219 size_t spStart = 0, spEnd = 0;
220 double sp_pos[3], sp_err[6];
221 for (
size_t i = 0; i < 3; i++) sp_pos[i] = 1.0;
222 for (
size_t i = 0; i < 6; i++) sp_err[i] = 1.0;
225 std::map< size_t, std::vector< art::Ptr<recob::Track> > > pfPartToTrackVecMap;
229 tracks->reserve(result.size());
230 for (
size_t trkIndex = 0; trkIndex < result.size(); ++trkIndex)
241 if (trk->
size()<2)
continue;
247 size_t trkIdx = tracks->size() - 1;
248 art::ProductID trkId = getProductID< std::vector<recob::Track> >();
252 for (
size_t h = 0, cnt = 0; h < trk->
size(); h++)
258 trk2hit->addSingle(trkPtr, h3d->
Hit2DPtr(), metadata);
259 trk2hit_oldway->addSingle(trkPtr, h3d->
Hit2DPtr());
263 spStart = allsp->
size();
264 for (
size_t h = 0; h < trk->
size(); ++h)
269 double hx = h3d->
Point3D().X();
270 double hy = h3d->
Point3D().Y();
271 double hz = h3d->
Point3D().Z();
274 (std::fabs(sp_pos[0] - hx) > 1.0e-5) ||
275 (std::fabs(sp_pos[1] - hy) > 1.0e-5) ||
276 (std::fabs(sp_pos[2] - hz) > 1.0e-5))
283 sp_pos[0] = hx; sp_pos[1] = hy; sp_pos[2] = hz;
293 spEnd = allsp->size();
295 if (spEnd > spStart)
util::CreateAssn(*
this, evt, *tracks, *allsp, *trk2sp, spStart, spEnd);
298 if (result[trkIndex].Key() > -1)
300 size_t trackIdx = tracks->size() - 1;
301 art::ProductID trackId = getProductID< std::vector<recob::Track> >();
303 pfPartToTrackVecMap[result[trkIndex].Key()].push_back(trackPtr);
307 auto vid = getProductID< std::vector<recob::Vertex> >();
308 auto kid = getProductID< std::vector<recob::Vertex> >(
kKinksName);
309 auto const* kinkGetter = evt.productGetter(kid);
311 auto tid = getProductID< std::vector<recob::Track> >();
312 auto const* trkGetter = evt.productGetter(tid);
315 auto ksel = pmalgFitter.getKinks();
316 std::map< size_t, art::Ptr<recob::Vertex> > frontVtxs;
321 for (
auto const & v : vsel)
323 xyz[0] = v.first.X(); xyz[1] = v.first.Y(); xyz[2] = v.first.Z();
325 <<
" vtx:" << xyz[0] <<
":" << xyz[1] <<
":" << xyz[2]
326 <<
" (" << v.second.size() <<
" tracks)";
328 size_t vidx = vtxs->size();
332 if (vptr.isNull())
mf::LogWarning(
"PMAlgTrajFitter") <<
"Vertex ptr is null.";
333 if (!v.second.empty())
335 for (
const auto & vEntry : v.second)
337 size_t tidx = vEntry.first;
338 bool isFront = vEntry.second;
340 if (isFront) frontVtxs[tidx] = vptr;
343 vtx2trk->addSingle(vptr, tptr);
346 else mf::LogWarning(
"PMAlgTrajFitter") <<
"No tracks found at this vertex.";
350 for (
auto const & k : ksel)
352 xyz[0] = k.first.X(); xyz[1] = k.first.Y(); xyz[2] = k.first.Z();
353 mf::LogVerbatim(
"Summary") <<
" kink:" << xyz[0] <<
":" << xyz[1] <<
":" << xyz[2];
355 size_t kidx = kinks->size();
356 size_t tidx = k.second;
362 trk2kink->addSingle(tptr, kptr);
370 for (
size_t t = 0; t < result.size(); ++t)
372 auto const & trk = *(result[t].Track());
373 for (
auto const * node : trk.
Nodes())
375 xyz[0] = node->Point3D().X(); xyz[1] = node->Point3D().Y(); xyz[2] = node->Point3D().Z();
381 for (
const auto & pfParticleItr : pfPartToTrackVecMap)
384 mf::LogVerbatim(
"PMAlgTrajFitter") <<
"PFParticle key: " << pfParticle.key()
385 <<
", self: " << pfParticle->Self() <<
", #tracks: " << pfParticleItr.second.size();
387 if (!pfParticle.isNull())
util::CreateAssn(*
this, evt, pfParticle, pfParticleItr.second, *pfp2trk);
388 else mf::LogError(
"PMAlgTrajFitter") <<
"Error in PFParticle lookup, pfparticle index: " 389 << pfParticleItr.first <<
", key: " << pfParticle.key();
393 evt.put(std::move(tracks));
394 evt.put(std::move(allsp));
395 evt.put(std::move(vtxs));
399 evt.put(std::move(trk2hit_oldway));
400 evt.put(std::move(trk2hit));
401 evt.put(std::move(trk2sp));
403 evt.put(std::move(sp2hit));
404 evt.put(std::move(vtx2trk));
407 evt.put(std::move(pfp2trk));
bool SelectHits(float fraction=1.0F)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool HasPlane(unsigned int iplane) const
Returns whether a plane with index iplane is present in this TPC.
recob::Track convertFrom(const pma::Track3D &src, unsigned int tidx, int pdg=0)
unsigned int FrontTPC(void) const
bool fSaveOnlyBranchingVtx
bool IsEnabled(void) const
std::vector< pma::Node3D * > const & Nodes(void) const
Planes which measure Z direction.
art::ServiceHandle< geo::Geometry > fGeom
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
art::Ptr< recob::Hit > const & Hit2DPtr(void) const
Definition of vertex object for LArSoft.
void push_back(Ptr< U > const &p)
art::InputTag fPfpModuleLabel
fhicl::Atom< bool > RunVertexing
size_t CompleteMissingWires(unsigned int view)
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
unsigned int FrontCryo(void) const
void push_back(pma::Hit3D *hit)
TVector3 const & Point3D(void) const
pma::PMAlgVertexing::Config fPmaVtxConfig
static const std::string kNodesName
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
pma::ProjectionMatchingAlg::Config fPmaConfig
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
static const std::string kKinksName
pma::PMAlgFitter::Config fPmaFitterConfig
cet::coded_exception< error, detail::translate > exception
art::InputTag fHitModuleLabel