LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
PMAlgTrajFitter_module.cc
Go to the documentation of this file.
1 // Class: PMAlgTrajFitter
3 // Module Type: producer
4 // File: PMAlgTrajFitter_module.cc
5 // Author: D.Stefan (Dorota.Stefan@ncbj.gov.pl) and R.Sulej (Robert.Sulej@cern.ch), May 2015
6 //
7 // Creates 3D tracks using Projection Matching Algorithm, please see
8 // RecoAlg/ProjectionMatchingAlg.h for basics of the PMA algorithm and its settings.
9 //
10 // Progress:
11 // July 2016: fit-only module separated from PMAlgTrackMaker_module, common tracking
12 // functionality moved to algorithm class;
13 // note, that some part of vertexing can be still reasonable in this module:
14 // - use hierarchy of input PFP's to create vertices, join tracks, reoptimize
15 // - detect vertices/kinks on fitted trajectories
16 //
18 
23 #include "fhiclcpp/types/Atom.h"
24 #include "fhiclcpp/types/Table.h"
26 
27 // LArSoft includes
41 
43 
44 #include <memory>
45 
46 namespace trkf {
47 
49  public:
50  struct Config {
51  using Name = fhicl::Name;
53 
55  Name("ProjectionMatchingAlg")};
56 
58 
60 
62  Name("HitModuleLabel"),
63  Comment("tag of unclustered hits, which were used to produce PFPs and clusters")};
64 
66  Name("PfpModuleLabel"),
67  Comment("tag of the input PFParticles and associated clusters")};
68 
70  Name("SaveOnlyBranchingVtx"),
71  Comment("use true to save only vertices interconnecting many tracks, otherwise vertex is "
72  "added to the front of each track")};
73 
75  Name("SavePmaNodes"),
76  Comment("save track nodes (only for algorithm development purposes)")};
77  };
79 
80  explicit PMAlgTrajFitter(Parameters const& config);
81 
82  PMAlgTrajFitter(PMAlgTrajFitter const&) = delete;
83  PMAlgTrajFitter(PMAlgTrajFitter&&) = delete;
84  PMAlgTrajFitter& operator=(PMAlgTrajFitter const&) = delete;
86 
87  private:
88  void produce(art::Event& e) override;
89 
90  // ******************** fcl parameters ***********************
91  art::InputTag fHitModuleLabel; // tag for unclustered hit collection
92  art::InputTag fPfpModuleLabel; // tag for PFParticle and cluster collections
93 
97 
98  bool fSaveOnlyBranchingVtx; // for debugging, save only vertices which connect many tracks
99  bool fSavePmaNodes; // for debugging, save only track nodes
100 
101  // ********** instance names (collections, assns) ************
102  static const std::string kKinksName; // kinks on tracks
103  static const std::string kNodesName; // pma nodes
104 
105  // *********************** geometry **************************
107  };
108  // -------------------------------------------------------------
109  const std::string PMAlgTrajFitter::kKinksName = "kink";
110  const std::string PMAlgTrajFitter::kNodesName = "node";
111  // -------------------------------------------------------------
112 
114  : EDProducer{config}
115  , fHitModuleLabel(config().HitModuleLabel())
116  , fPfpModuleLabel(config().PfpModuleLabel())
117  , fPmaConfig(config().ProjectionMatchingAlg())
118  , fPmaFitterConfig(config().PMAlgFitting())
119  , fPmaVtxConfig(config().PMAlgVertexing())
120  , fSaveOnlyBranchingVtx(config().SaveOnlyBranchingVtx())
121  , fSavePmaNodes(config().SavePmaNodes())
122  {
123  produces<std::vector<recob::Track>>();
124  produces<std::vector<recob::SpacePoint>>();
125  produces<std::vector<recob::Vertex>>(); // no instance name for interaction vertices
126  produces<std::vector<recob::Vertex>>(kKinksName); // collection of kinks on tracks
127  produces<std::vector<recob::Vertex>>(kNodesName); // collection of pma nodes
128 
130  recob::Hit>>(); // ****** REMEMBER to remove when FindMany improved ******
131  produces<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta>>();
132 
133  produces<art::Assns<recob::Track, recob::SpacePoint>>();
134  produces<art::Assns<recob::SpacePoint, recob::Hit>>();
135  produces<
137  recob::Track>>(); // no instance name for assns of tracks to interaction vertices
138  produces<art::Assns<recob::Track, recob::Vertex>>(kKinksName); // assns of kinks to tracks
139 
140  produces<art::Assns<recob::PFParticle, recob::Track>>();
141  }
142  // ------------------------------------------------------
143 
145  {
146  // ---------------- Create data products ------------------
147  auto tracks = std::make_unique<std::vector<recob::Track>>();
148  auto allsp = std::make_unique<std::vector<recob::SpacePoint>>();
149  auto vtxs = std::make_unique<std::vector<recob::Vertex>>(); // interaction vertices
150  auto kinks = std::make_unique<
151  std::vector<recob::Vertex>>(); // kinks on tracks (no new particles start in kinks)
152  auto nodes = std::make_unique<std::vector<recob::Vertex>>(); // pma nodes
153 
154  auto trk2hit_oldway = std::make_unique<
156  recob::Hit>>(); // ****** REMEMBER to remove when FindMany improved ******
157  auto trk2hit = std::make_unique<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta>>();
158 
159  auto trk2sp = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
160 
161  auto sp2hit = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
162  auto vtx2trk = std::make_unique<
164  recob::Track>>(); // one or more tracks (particles) start in the vertex
165  auto trk2kink =
166  std::make_unique<art::Assns<recob::Track, recob::Vertex>>(); // one or more kinks on the track
167 
168  auto pfp2trk = std::make_unique<art::Assns<recob::PFParticle, recob::Track>>();
169 
170  // ------------------- Collect inputs ---------------------
171  art::Handle<std::vector<recob::Hit>> allHitListHandle;
174  std::vector<art::Ptr<recob::Hit>> allhitlist;
175  if (!(evt.getByLabel(fHitModuleLabel,
176  allHitListHandle) && // all hits used to make clusters and PFParticles
177  evt.getByLabel(fPfpModuleLabel, cluListHandle) && // clusters associated to PFParticles
178  evt.getByLabel(fPfpModuleLabel, pfparticleHandle))) // and finally PFParticles
179  {
180  throw cet::exception("PMAlgTrajFitter")
181  << "Not all required data products found in the event." << std::endl;
182  }
183 
184  art::fill_ptr_vector(allhitlist, allHitListHandle);
185 
186  art::FindManyP<recob::Hit> hitsFromClusters(cluListHandle, evt, fPfpModuleLabel);
187  art::FindManyP<recob::Cluster> clustersFromPfps(pfparticleHandle, evt, fPfpModuleLabel);
188  art::FindManyP<recob::Vertex> vtxFromPfps(pfparticleHandle, evt, fPfpModuleLabel);
189 
190  auto const detProp =
192 
193  // -------------- PMA Fitter for this event ---------------
194  auto pmalgFitter = pma::PMAlgFitter(allhitlist,
195  *cluListHandle,
196  *pfparticleHandle,
197  hitsFromClusters,
198  clustersFromPfps,
199  vtxFromPfps,
200  fPmaConfig,
202  fPmaVtxConfig);
203 
204  // ------------------ Do the job here: --------------------
205  int retCode = pmalgFitter.build(detProp);
206  // --------------------------------------------------------
207  switch (retCode) {
208  case -2: mf::LogError("Summary") << "problem"; break;
209  case -1: mf::LogWarning("Summary") << "no input"; break;
210  case 0: mf::LogVerbatim("Summary") << "no tracks done"; break;
211  default:
212  if (retCode < 0)
213  mf::LogVerbatim("Summary") << "unknown result";
214  else if (retCode == 1)
215  mf::LogVerbatim("Summary") << retCode << " track ready";
216  else
217  mf::LogVerbatim("Summary") << retCode << " tracks ready";
218  break;
219  }
220 
221  // ---------- Translate output to data products: ----------
222  auto const& result = pmalgFitter.result();
223  if (!result.empty()) // ok, there is something to save
224  {
225  size_t spStart = 0, spEnd = 0;
226  double sp_pos[3], sp_err[6];
227  for (size_t i = 0; i < 3; i++)
228  sp_pos[i] = 1.0;
229  for (size_t i = 0; i < 6; i++)
230  sp_err[i] = 1.0;
231 
232  // use the following to create PFParticle <--> Track associations;
233  std::map<size_t, std::vector<art::Ptr<recob::Track>>> pfPartToTrackVecMap;
234 
235  //auto const make_trkptr = art::PtrMaker<recob::Track>(evt, *this); // PtrMaker Step #1
236 
237  tracks->reserve(result.size());
238  for (size_t trkIndex = 0; trkIndex < result.size(); ++trkIndex) {
239  pma::Track3D* trk = result[trkIndex].Track();
240 
241  trk->SelectHits(); // just in case, set all to enabled
242  unsigned int itpc = trk->FrontTPC(), icryo = trk->FrontCryo();
243  geo::TPCID const tpcid{icryo, itpc};
244  auto has_plane = [this, &tpcid](auto const view) {
245  return wireReadoutGeom.HasPlane(geo::PlaneID(tpcid, view));
246  };
247 
248  if (has_plane(geo::kU)) trk->CompleteMissingWires(detProp, geo::kU);
249  if (has_plane(geo::kV)) trk->CompleteMissingWires(detProp, geo::kV);
250  if (has_plane(geo::kZ)) trk->CompleteMissingWires(detProp, geo::kZ);
251 
252  // gc: make sure no tracks are created with less than 2 points
253  if (trk->size() < 2) continue;
254 
255  tracks->push_back(pma::convertFrom(*trk, trkIndex));
256 
257  //auto const trkPtr = make_trkptr(tracks->size() - 1); // PtrMaker Step #2
258 
259  size_t trkIdx = tracks->size() - 1; // stuff for assns:
260  art::ProductID trkId = evt.getProductID<std::vector<recob::Track>>();
261  art::Ptr<recob::Track> trkPtr(trkId, trkIdx, evt.productGetter(trkId));
262 
263  //gc: save associated hits in the same order as trajectory points
264  for (size_t h = 0, cnt = 0; h < trk->size(); h++) {
265  pma::Hit3D* h3d = (*trk)[h];
266  if (!h3d->IsEnabled()) continue;
267 
268  recob::TrackHitMeta metadata(cnt++, h3d->Dx());
269  trk2hit->addSingle(trkPtr, h3d->Hit2DPtr(), metadata);
270  trk2hit_oldway->addSingle(
271  trkPtr, h3d->Hit2DPtr()); // ****** REMEMBER to remove when FindMany improved ******
272  }
273 
275  spStart = allsp->size();
276  for (size_t h = 0; h < trk->size(); ++h) {
277  pma::Hit3D* h3d = (*trk)[h];
278  if (!h3d->IsEnabled()) continue;
279 
280  double hx = h3d->Point3D().X();
281  double hy = h3d->Point3D().Y();
282  double hz = h3d->Point3D().Z();
283 
284  if ((h == 0) || (std::fabs(sp_pos[0] - hx) > 1.0e-5) ||
285  (std::fabs(sp_pos[1] - hy) > 1.0e-5) || (std::fabs(sp_pos[2] - hz) > 1.0e-5)) {
286  if (sp_hits.size()) // hits assigned to the previous sp
287  {
288  util::CreateAssn(evt, *allsp, sp_hits, *sp2hit);
289  sp_hits.clear();
290  }
291  sp_pos[0] = hx;
292  sp_pos[1] = hy;
293  sp_pos[2] = hz;
294  allsp->push_back(recob::SpacePoint(sp_pos, sp_err, 1.0));
295  }
296  sp_hits.push_back(h3d->Hit2DPtr());
297  }
298 
299  if (sp_hits.size()) // hits assigned to the last sp
300  {
301  util::CreateAssn(evt, *allsp, sp_hits, *sp2hit);
302  }
303  spEnd = allsp->size();
304 
305  if (spEnd > spStart) util::CreateAssn(evt, *tracks, *allsp, *trk2sp, spStart, spEnd);
306 
307  // if there is a PFParticle collection then recover PFParticle and add info to map
308  if (result[trkIndex].Key() > -1) {
309  size_t trackIdx = tracks->size() - 1;
310  art::ProductID trackId = evt.getProductID<std::vector<recob::Track>>();
311  art::Ptr<recob::Track> trackPtr(trackId, trackIdx, evt.productGetter(trackId));
312  pfPartToTrackVecMap[result[trkIndex].Key()].push_back(trackPtr);
313  }
314  }
315 
316  auto vid = evt.getProductID<std::vector<recob::Vertex>>();
317  auto kid = evt.getProductID<std::vector<recob::Vertex>>(kKinksName);
318  auto const* kinkGetter = evt.productGetter(kid);
319 
320  auto tid = evt.getProductID<std::vector<recob::Track>>();
321  auto const* trkGetter = evt.productGetter(tid);
322 
323  auto vsel = pmalgFitter.getVertices(
324  fSaveOnlyBranchingVtx); // vtx pos's with vector of connected track idxs
325  auto ksel = pmalgFitter.getKinks(); // pairs of kink position - associated track idx
326  std::map<size_t, art::Ptr<recob::Vertex>> frontVtxs; // front vertex ptr for each track index
327 
328  if (fPmaFitterConfig.RunVertexing()) // save vertices and vtx-trk assns
329  {
330  double xyz[3];
331  for (auto const& v : vsel) {
332  xyz[0] = v.first.X();
333  xyz[1] = v.first.Y();
334  xyz[2] = v.first.Z();
335  mf::LogVerbatim("Summary") << " vtx:" << xyz[0] << ":" << xyz[1] << ":" << xyz[2]
336  << " (" << v.second.size() << " tracks)";
337 
338  size_t vidx = vtxs->size();
339  vtxs->push_back(recob::Vertex(xyz, vidx));
340 
341  art::Ptr<recob::Vertex> vptr(vid, vidx, evt.productGetter(vid));
342  if (vptr.isNull()) mf::LogWarning("PMAlgTrajFitter") << "Vertex ptr is null.";
343  if (!v.second.empty()) {
344  for (const auto& vEntry : v.second) {
345  size_t tidx = vEntry.first;
346  bool isFront = vEntry.second;
347 
348  if (isFront) frontVtxs[tidx] = vptr; // keep ptr of the front vtx
349 
350  art::Ptr<recob::Track> tptr(tid, tidx, trkGetter);
351  vtx2trk->addSingle(vptr, tptr);
352  }
353  }
354  else
355  mf::LogWarning("PMAlgTrajFitter") << "No tracks found at this vertex.";
356  }
357  mf::LogVerbatim("Summary") << vtxs->size() << " vertices ready";
358 
359  for (auto const& k : ksel) {
360  xyz[0] = k.first.X();
361  xyz[1] = k.first.Y();
362  xyz[2] = k.first.Z();
363  mf::LogVerbatim("Summary") << " kink:" << xyz[0] << ":" << xyz[1] << ":" << xyz[2];
364 
365  size_t kidx = kinks->size();
366  size_t tidx = k.second; // track idx on which this kink was found
367 
368  kinks->push_back(
369  recob::Vertex(xyz, tidx)); // save index of track (will have color of trk in evd)
370 
371  art::Ptr<recob::Track> tptr(tid, tidx, trkGetter);
372  art::Ptr<recob::Vertex> kptr(kid, kidx, kinkGetter);
373  trk2kink->addSingle(tptr, kptr);
374  }
375  mf::LogVerbatim("Summary") << ksel.size() << " kinks ready";
376  }
377 
378  if (fSavePmaNodes) {
379  double xyz[3];
380  for (size_t t = 0; t < result.size(); ++t) {
381  auto const& trk = *(result[t].Track());
382  for (auto const* node : trk.Nodes()) {
383  xyz[0] = node->Point3D().X();
384  xyz[1] = node->Point3D().Y();
385  xyz[2] = node->Point3D().Z();
386  nodes->push_back(recob::Vertex(xyz, t));
387  }
388  }
389  }
390 
391  for (const auto& pfParticleItr : pfPartToTrackVecMap) {
392  art::Ptr<recob::PFParticle> pfParticle(pfparticleHandle, pfParticleItr.first);
393  mf::LogVerbatim("PMAlgTrajFitter")
394  << "PFParticle key: " << pfParticle.key() << ", self: " << pfParticle->Self()
395  << ", #tracks: " << pfParticleItr.second.size();
396 
397  if (!pfParticle.isNull())
398  util::CreateAssn(evt, pfParticle, pfParticleItr.second, *pfp2trk);
399  else
400  mf::LogError("PMAlgTrajFitter")
401  << "Error in PFParticle lookup, pfparticle index: " << pfParticleItr.first
402  << ", key: " << pfParticle.key();
403  }
404  }
405 
406  evt.put(std::move(tracks));
407  evt.put(std::move(allsp));
408  evt.put(std::move(vtxs));
409  evt.put(std::move(kinks), kKinksName);
410  evt.put(std::move(nodes), kNodesName);
411 
412  evt.put(std::move(trk2hit_oldway)); // ****** REMEMBER to remove when FindMany improved ******
413  evt.put(std::move(trk2hit));
414  evt.put(std::move(trk2sp));
415 
416  evt.put(std::move(sp2hit));
417  evt.put(std::move(vtx2trk));
418  evt.put(std::move(trk2kink), kKinksName);
419 
420  evt.put(std::move(pfp2trk));
421  }
422 
424 
425 } // namespace trkf
fhicl::Table< pma::ProjectionMatchingAlg::Config > ProjectionMatchingAlg
bool SelectHits(float fraction=1.0F)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
size_t Self() const
Returns the index of this particle.
Definition: PFParticle.h:88
recob::Track convertFrom(const pma::Track3D &src, unsigned int tidx, int pdg=0)
ProductID getProductID(std::string const &instance_name="") const
Planes which measure V.
Definition: geo_types.h:132
Declaration of signal hit object.
bool IsEnabled() const noexcept
Definition: PmaHit3D.h:86
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
void produces(std::string const &instanceName={}, Persistable const persistable=Persistable::Yes)
Class to keep data related to recob::Hit associated with recob::Track.
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
Planes which measure Z direction.
Definition: geo_types.h:134
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
Data related to recob::Hit associated with recob::Track.The purpose is to collect several variables t...
Definition: TrackHitMeta.h:43
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool HasPlane(PlaneID const &planeid) const
Returns whether we have the specified plane.
void produce(art::Event &e) override
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
Planes which measure U.
Definition: geo_types.h:131
EDProductGetter const * productGetter(ProductID const pid) const
fhicl::Table< pma::PMAlgFitter::Config > PMAlgFitting
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
geo::WireReadoutGeom const & wireReadoutGeom
Interface for a class providing readout channel mapping to geometry.
fhicl::Atom< bool > RunVertexing
key_type key() const noexcept
Definition: Ptr.h:166
bool isNull() const noexcept
Definition: Ptr.h:211
fhicl::Table< pma::PMAlgVertexing::Config > PMAlgVertexing
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:118
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:119
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
Declaration of cluster object.
size_type size() const
Definition: PtrVector.h:302
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.
pma::PMAlgVertexing::Config fPmaVtxConfig
static const std::string kNodesName
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Provides recob::Track data product.
double Dx() const noexcept
Definition: PmaHit3D.h:69
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
pma::ProjectionMatchingAlg::Config fPmaConfig
art::Ptr< recob::Hit > const & Hit2DPtr() const
Definition: PmaHit3D.h:48
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:992
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
size_t size() const
Definition: PmaTrack3D.h:89
PMAlgTrajFitter(Parameters const &config)
static const std::string kKinksName
fhicl::Atom< art::InputTag > HitModuleLabel
void clear()
Definition: PtrVector.h:533
PMAlgTrajFitter & operator=(PMAlgTrajFitter const &)=delete
art framework interface to geometry description
pma::PMAlgFitter::Config fPmaFitterConfig
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane .
fhicl::Atom< art::InputTag > PfpModuleLabel