LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TrajCluster_module.cc
Go to the documentation of this file.
1 
9 // C/C++ standard libraries
10 #include <string>
11 
12 // Framework libraries
16 #include "art_root_io/TFileService.h"
19 #include "fhiclcpp/ParameterSet.h"
20 
21 // LArSoft includes
37 
38 //root includes
39 #include "TTree.h"
40 
41 // ... more includes in the implementation section
42 
43 namespace cluster {
56  class TrajCluster : public art::EDProducer {
57  public:
58  explicit TrajCluster(fhicl::ParameterSet const& pset);
59 
60  private:
61  void produce(art::Event& evt) override;
62  void beginJob() override;
63  void endJob() override;
64 
65  tca::TrajClusterAlg fTCAlg; // define TrajClusterAlg object
66  TTree* showertree;
67  void GetHits(const std::vector<recob::Hit>& inputHits,
68  const geo::TPCID& tpcid,
69  std::vector<std::vector<unsigned int>>& tpcHits);
70  void GetHits(const geo::TPCID& tpcid,
71  const std::vector<recob::Slice>& inputSlices,
72  art::FindManyP<recob::Hit>& hitFromSlc,
73  std::vector<std::vector<unsigned int>>& tpcHits,
74  std::vector<int>& slcIDs);
75 
80 
81  unsigned int fMaxSliceHits;
85  }; // class TrajCluster
86 
87 } // namespace cluster
88 
89 //******************************************************************************
90 //*** implementation
91 //***
92 
93 // C/C++ standard libraries
94 #include <memory> // std::move()
95 
96 // Framework libraries
101 
102 //LArSoft includes
104 #include "lardata/ArtDataHelper/HitCreator.h" // recob::HitCollectionAssociator
108 #include "lardataobj/RecoBase/Hit.h"
112 
113 namespace cluster {
114 
115  struct HitLoc {
116  unsigned int index; // index of this entry in a sort vector
117  unsigned int ctp; // encoded Cryostat, TPC and Plane
118  unsigned int wire;
119  int tick; // hit StartTick using typedef int TDCtick_t in RawTypes.h
120  short localIndex; // defined in Hit.h
121  };
122 
123  //----------------------------------------------------------------------------
124  bool SortHits(HitLoc const& h1, HitLoc const& h2)
125  {
126  // sort by hit location (Cryostat, TPC, Plane, Wire, StartTick, hit LocalIndex)
127  if (h1.ctp != h2.ctp) return h1.ctp < h2.ctp;
128  if (h1.wire != h2.wire) return h1.wire < h2.wire;
129  if (h1.tick != h2.tick) return h1.tick < h2.tick;
130  return h1.localIndex < h2.localIndex;
131  } // SortHits
132 
133  //----------------------------------------------------------------------------
135  : EDProducer{pset}, fTCAlg{pset.get<fhicl::ParameterSet>("TrajClusterAlg")}
136  {
137  fHitModuleLabel = "NA";
138  if (pset.has_key("HitModuleLabel")) fHitModuleLabel = pset.get<art::InputTag>("HitModuleLabel");
139  fSliceModuleLabel = "NA";
140  if (pset.has_key("SliceModuleLabel"))
141  fSliceModuleLabel = pset.get<art::InputTag>("SliceModuleLabel");
142  fMaxSliceHits = UINT_MAX;
143  if (pset.has_key("MaxSliceHits")) fMaxSliceHits = pset.get<unsigned int>("MaxSliceHits");
144  fSpacePointModuleLabel = "NA";
145  if (pset.has_key("SpacePointModuleLabel"))
146  fSpacePointModuleLabel = pset.get<art::InputTag>("SpacePointModuleLabel");
148  if (pset.has_key("SpacePointHitAssnLabel"))
149  fSpacePointHitAssnLabel = pset.get<art::InputTag>("SpacePointHitAssnLabel");
150  fDoWireAssns = pset.get<bool>("DoWireAssns", true);
151  fDoRawDigitAssns = pset.get<bool>("DoRawDigitAssns", true);
152  fSaveAll2DVertices = false;
153  if (pset.has_key("SaveAll2DVertices")) fSaveAll2DVertices = pset.get<bool>("SaveAll2DVertices");
154 
155  // let HitCollectionAssociator declare that we are going to produce
156  // hits and associations with wires and raw digits
157  // (with no particular product label)
159  producesCollector(), "", fDoWireAssns, fDoRawDigitAssns);
160 
161  produces<std::vector<recob::Cluster>>();
162  produces<std::vector<recob::Vertex>>();
163  produces<std::vector<recob::EndPoint2D>>();
164  produces<std::vector<recob::Seed>>();
165  produces<std::vector<recob::Shower>>();
166  produces<art::Assns<recob::Cluster, recob::Hit>>();
167  produces<art::Assns<recob::Cluster, recob::EndPoint2D, unsigned short>>();
168  produces<art::Assns<recob::Cluster, recob::Vertex, unsigned short>>();
169  produces<art::Assns<recob::Shower, recob::Hit>>();
170 
171  produces<std::vector<recob::PFParticle>>();
172  produces<art::Assns<recob::PFParticle, recob::Cluster>>();
173  produces<art::Assns<recob::PFParticle, recob::Shower>>();
174  produces<art::Assns<recob::PFParticle, recob::Vertex>>();
175  produces<art::Assns<recob::PFParticle, recob::Seed>>();
176 
177  produces<art::Assns<recob::Slice, recob::Cluster>>();
178  produces<art::Assns<recob::Slice, recob::PFParticle>>();
179  produces<art::Assns<recob::Slice, recob::Hit>>();
180 
181  produces<std::vector<anab::CosmicTag>>();
182  produces<art::Assns<recob::PFParticle, anab::CosmicTag>>();
183 
184  // www: declear/create SpacePoint and association between SpacePoint and Hits from TrajCluster (Hit->SpacePoint)
185  produces<art::Assns<recob::SpacePoint, recob::Hit>>();
186  } // TrajCluster::TrajCluster()
187 
188  //----------------------------------------------------------------------------
190  {
192 
193  showertree = tfs->make<TTree>("showervarstree", "showerVarsTree");
195  }
196 
197  //----------------------------------------------------------------------------
199  {
200  std::vector<unsigned int> const& fAlgModCount = fTCAlg.GetAlgModCount();
201  std::vector<std::string> const& fAlgBitNames = fTCAlg.GetAlgBitNames();
202  if (fAlgBitNames.size() != fAlgModCount.size()) return;
203  mf::LogVerbatim myprt("TC");
204  myprt << "TrajCluster algorithm counts\n";
205  unsigned short icol = 0;
206  for (unsigned short ib = 0; ib < fAlgModCount.size(); ++ib) {
207  if (ib == tca::kKilled) continue;
208  myprt << std::left << std::setw(18) << fAlgBitNames[ib] << std::right << std::setw(10)
209  << fAlgModCount[ib] << " ";
210  ++icol;
211  if (icol == 4) {
212  myprt << "\n";
213  icol = 0;
214  }
215  } // ib
216  } // endJob
217 
218  //----------------------------------------------------------------------------
220  {
221  // Get a single hit collection from a HitsModuleLabel or multiple sets of "sliced" hits
222  // (aka clusters of hits that are close to each other in 3D) from a SliceModuleLabel.
223  // A pointer to the full hit collection is passed to TrajClusterAlg. The hits that are
224  // in each slice are reconstructed to find 2D trajectories (that become clusters),
225  // 2D vertices (EndPoint2D), 3D vertices, PFParticles and Showers. These data products
226  // are then collected and written to the event. Each slice is considered as an independent
227  // collection of hits with the additional requirement that all hits in a slice reside in
228  // one TPC
229 
230  // pointers to the slices in the event
231  std::vector<art::Ptr<recob::Slice>> slices;
232  std::vector<int> slcIDs;
233  unsigned int nInputHits = 0;
234 
235  // get a reference to the Hit collection
236  auto inputHits = art::Handle<std::vector<recob::Hit>>();
237  if (!evt.getByLabel(fHitModuleLabel, inputHits))
238  throw cet::exception("TrajClusterModule")
239  << "Failed to get a handle to hit collection '" << fHitModuleLabel.label() << "'\n";
240  nInputHits = (*inputHits).size();
241  if (!fTCAlg.SetInputHits(*inputHits, evt.run(), evt.event()))
242  throw cet::exception("TrajClusterModule")
243  << "Failed to process hits from '" << fHitModuleLabel.label() << "'\n";
244  // Try to determine the source of the hit collection using the assumption that it was
245  // derived from gaushit. If this is successful, pass the handle to TrajClusterAlg to
246  // recover hits that were incorrectly removed by disambiguation (DUNE)
247  if (fHitModuleLabel != "gaushit") {
248  auto sourceHits = art::Handle<std::vector<recob::Hit>>();
249  art::InputTag sourceModuleLabel("gaushit");
250  if (evt.getByLabel(sourceModuleLabel, sourceHits)) fTCAlg.SetSourceHits(*sourceHits);
251  } // look for gaushit collection
252 
253  // get an optional reference to the Slice collection
254  auto inputSlices = art::Handle<std::vector<recob::Slice>>();
255  if (fSliceModuleLabel != "NA") {
257  if (!evt.getByLabel(fSliceModuleLabel, inputSlices))
258  throw cet::exception("TrajClusterModule") << "Failed to get a inputSlices";
259  } // fSliceModuleLabel specified
260 
261  // get an optional reference to the SpacePoint collection
262  auto InputSpts = art::Handle<std::vector<recob::SpacePoint>>();
263  if (fSpacePointModuleLabel != "NA") {
264  if (!evt.getByLabel(fSpacePointModuleLabel, InputSpts))
265  throw cet::exception("TrajClusterModule") << "Failed to get a handle to SpacePoints\n";
266  tca::evt.sptHits.resize((*InputSpts).size(), {{UINT_MAX, UINT_MAX, UINT_MAX}});
267  art::FindManyP<recob::Hit> hitsFromSpt(InputSpts, evt, fSpacePointHitAssnLabel);
268  // TrajClusterAlg doesn't use the SpacePoint positions (only the assns to hits) but pass it
269  // anyway in case it is useful
270  fTCAlg.SetInputSpts(*InputSpts);
271  if (!hitsFromSpt.isValid())
272  throw cet::exception("TrajClusterModule")
273  << "Failed to get a handle to SpacePoint -> Hit assns\n";
274  // ensure that the assn is to the inputHit collection
275  auto& firstHit = hitsFromSpt.at(0)[0];
276  if (firstHit.id() != inputHits.id())
277  throw cet::exception("TrajClusterModule")
278  << "The SpacePoint -> Hit assn doesn't reference the input hit collection\n";
279  tca::evt.sptHits.resize((*InputSpts).size(), {{UINT_MAX, UINT_MAX, UINT_MAX}});
280  for (unsigned int isp = 0; isp < (*InputSpts).size(); ++isp) {
281  auto& hits = hitsFromSpt.at(isp);
282  for (unsigned short iht = 0; iht < hits.size(); ++iht) {
283  unsigned short plane = hits[iht]->WireID().Plane;
284  tca::evt.sptHits[isp][plane] = hits[iht].key();
285  } // iht
286  } // isp
287  } // fSpacePointModuleLabel specified
288 
289  if (nInputHits > 0) {
290  auto const clockData =
292  auto const detProp =
294  auto const* geom = lar::providerFrom<geo::Geometry>();
295  for (const auto& tpcid : geom->Iterate<geo::TPCID>()) {
296  // ignore protoDUNE dummy TPCs
297  if (geom->TPC(tpcid).DriftDistance() < 25.0) continue;
298  // a vector for the subset of hits in each slice in a TPC
299  // slice hits in this tpc
300  std::vector<std::vector<unsigned int>> sltpcHits;
301  if (inputSlices.isValid()) {
302  // get hits in this TPC and slice
303  art::FindManyP<recob::Hit> hitFromSlc(inputSlices, evt, fSliceModuleLabel);
304  GetHits(tpcid, *inputSlices, hitFromSlc, sltpcHits, slcIDs);
305  }
306  else {
307  // get hits in this TPC
308  // All hits are in one "fake" slice
309  GetHits(*inputHits, tpcid, sltpcHits);
310  slcIDs.resize(1);
311  slcIDs[0] = 1;
312  }
313  if (sltpcHits.empty()) continue;
314  for (unsigned short isl = 0; isl < sltpcHits.size(); ++isl) {
315  auto& tpcHits = sltpcHits[isl];
316  if (tpcHits.empty()) continue;
317  // only reconstruct slices with MC-matched hits?
318  // sort the slice hits by Cryostat, TPC, Wire, Plane, Start Tick and LocalIndex.
319  // This assumes that hits with larger LocalIndex are at larger Tick.
320  std::vector<HitLoc> sortVec(tpcHits.size());
321  for (unsigned int indx = 0; indx < tpcHits.size(); ++indx) {
322  auto& hit = (*inputHits)[tpcHits[indx]];
323  sortVec[indx].index = indx;
324  sortVec[indx].ctp = tca::EncodeCTP(hit.WireID());
325  sortVec[indx].wire = hit.WireID().Wire;
326  sortVec[indx].tick = hit.StartTick();
327  sortVec[indx].localIndex = hit.LocalIndex();
328  } // iht
329  std::sort(sortVec.begin(), sortVec.end(), SortHits);
330  std::vector tmp = tpcHits;
331  for (unsigned int ii = 0; ii < tpcHits.size(); ++ii)
332  tpcHits[ii] = tmp[sortVec[ii].index];
333  // clear the temp vector
334  tmp.resize(0);
335  sortVec.resize(0);
336  // look for a debug hit
337  if (tca::tcc.dbgStp) {
338  tca::debug.Hit = UINT_MAX;
339  for (unsigned short indx = 0; indx < tpcHits.size(); ++indx) {
340  auto& hit = (*inputHits)[tpcHits[indx]];
341  if ((int)hit.WireID().TPC == tca::debug.TPC &&
342  (int)hit.WireID().Plane == tca::debug.Plane &&
343  (int)hit.WireID().Wire == tca::debug.Wire &&
344  hit.PeakTime() > tca::debug.Tick - 10 && hit.PeakTime() < tca::debug.Tick + 10) {
345  std::cout << "Debug hit " << tpcHits[indx] << " found in slice ID " << slcIDs[isl];
346  std::cout << " RMS " << hit.RMS();
347  std::cout << " Multiplicity " << hit.Multiplicity();
348  std::cout << " GoodnessOfFit " << hit.GoodnessOfFit();
349  std::cout << "\n";
350  tca::debug.Hit = tpcHits[indx];
351  break;
352  } // Look for debug hit
353  } // iht
354  } // tca::tcc.dbgStp
355  fTCAlg.RunTrajClusterAlg(clockData, detProp, tpcHits, slcIDs[isl]);
356  } // isl
357  } // TPC
358  // stitch PFParticles between TPCs, create PFP start vertices, etc
360  if (tca::tcc.dbgSummary) tca::PrintAll(detProp, "TCM");
361  } // nInputHits > 0
362 
363  // Vectors to hold all data products that will go into the event
364  std::vector<recob::Hit> hitCol; // output hit collection
365  std::vector<recob::Cluster> clsCol;
366  std::vector<recob::PFParticle> pfpCol;
367  std::vector<recob::Vertex> vx3Col;
368  std::vector<recob::EndPoint2D> vx2Col;
369  std::vector<recob::Seed> sedCol;
370  std::vector<recob::Shower> shwCol;
371  std::vector<anab::CosmicTag> ctCol;
372  // a vector to correlate inputHits with output hits
373  std::vector<unsigned int> newIndex(nInputHits, UINT_MAX);
374 
375  // assns for those data products
376  // Cluster -> ...
377  std::unique_ptr<art::Assns<recob::Cluster, recob::Hit>> cls_hit_assn(
379  // unsigned short is the end to which a vertex is attached
380  std::unique_ptr<art::Assns<recob::Cluster, recob::EndPoint2D, unsigned short>> cls_vx2_assn(
382  std::unique_ptr<art::Assns<recob::Cluster, recob::Vertex, unsigned short>> cls_vx3_assn(
384  // Shower -> ...
385  std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> shwr_hit_assn(
387  // PFParticle -> ...
388  std::unique_ptr<art::Assns<recob::PFParticle, recob::Cluster>> pfp_cls_assn(
390  std::unique_ptr<art::Assns<recob::PFParticle, recob::Shower>> pfp_shwr_assn(
392  std::unique_ptr<art::Assns<recob::PFParticle, recob::Vertex>> pfp_vx3_assn(
394  std::unique_ptr<art::Assns<recob::PFParticle, anab::CosmicTag>> pfp_cos_assn(
396  std::unique_ptr<art::Assns<recob::PFParticle, recob::Seed>> pfp_sed_assn(
398  // Slice -> ...
399  std::unique_ptr<art::Assns<recob::Slice, recob::Cluster>> slc_cls_assn(
401  std::unique_ptr<art::Assns<recob::Slice, recob::PFParticle>> slc_pfp_assn(
403  std::unique_ptr<art::Assns<recob::Slice, recob::Hit>> slc_hit_assn(
405  // www: Hit -> SpacePoint
406  std::unique_ptr<art::Assns<recob::SpacePoint, recob::Hit>> sp_hit_assn(
408 
409  // temp struct to get the index of a 2D (or 3D vertex) into vx2Col (or vx3Col)
410  // given a slice index and a vertex ID (not UID)
411  struct slcVxStruct {
412  unsigned short slIndx;
413  int ID;
414  unsigned short vxColIndx;
415  };
416  std::vector<slcVxStruct> vx2StrList;
417  // vector to map 3V UID -> ID in each sub-slice
418  std::vector<slcVxStruct> vx3StrList;
419 
420  if (nInputHits > 0) {
421  unsigned short nSlices = fTCAlg.GetSlicesSize();
422  // define a hit collection begin index to pass to CreateAssn for each cluster
423  unsigned int hitColBeginIndex = 0;
424  for (unsigned short isl = 0; isl < nSlices; ++isl) {
425  unsigned short slcIndex = 0;
426  if (!slices.empty()) {
427  for (slcIndex = 0; slcIndex < slices.size(); ++slcIndex)
428  if (slices[slcIndex]->ID() == slcIDs[isl]) break;
429  if (slcIndex == slices.size()) continue;
430  }
431  auto& slc = fTCAlg.GetSlice(isl);
432  // See if there was a serious reconstruction failure that made the sub-slice invalid
433  if (!slc.isValid) continue;
434  // make EndPoint2Ds
435  for (auto& vx2 : slc.vtxs) {
436  if (vx2.ID <= 0) continue;
437  // skip complete 2D vertices?
438  if (!fSaveAll2DVertices && vx2.Vx3ID != 0) continue;
439  unsigned int vtxID = vx2.UID;
440  unsigned int wire = std::nearbyint(vx2.Pos[0]);
441  geo::PlaneID plID = tca::DecodeCTP(vx2.CTP);
442  geo::WireID wID = geo::WireID(plID.Cryostat, plID.TPC, plID.Plane, wire);
443  geo::View_t view = tca::tcc.geom->View(wID);
444  vx2Col.emplace_back((double)vx2.Pos[1] / tca::tcc.unitsPerTick, // Time
445  wID, // WireID
446  vx2.Score, // strength = score
447  vtxID, // ID
448  view, // View
449  0); // total charge - not relevant
450 
451  // fill the mapping struct
452  slcVxStruct tmp;
453  tmp.slIndx = isl;
454  tmp.ID = vx2.ID;
455  tmp.vxColIndx = vx2Col.size() - 1;
456  vx2StrList.push_back(tmp);
457 
458  } // vx2
459  // make Vertices
460  for (auto& vx3 : slc.vtx3s) {
461  if (vx3.ID <= 0) continue;
462  // ignore incomplete vertices
463  if (vx3.Wire >= 0) continue;
464  unsigned int vtxID = vx3.UID;
465  double xyz[3];
466  xyz[0] = vx3.X;
467  xyz[1] = vx3.Y;
468  xyz[2] = vx3.Z;
469  vx3Col.emplace_back(xyz, vtxID);
470  // fill the mapping struct
471  slcVxStruct tmp;
472  tmp.slIndx = isl;
473  tmp.ID = vx3.ID;
474  tmp.vxColIndx = vx3Col.size() - 1;
475  vx3StrList.push_back(tmp);
476  } // vx3
477  // Convert the tjs to clusters
478  for (auto& tj : slc.tjs) {
479  if (tj.AlgMod[tca::kKilled]) continue;
480  hitColBeginIndex = hitCol.size();
481  for (unsigned short ipt = tj.EndPt[0]; ipt <= tj.EndPt[1]; ++ipt) {
482  auto& tp = tj.Pts[ipt];
483  if (tp.Chg <= 0) continue;
484  // index of inputHits indices of hits used in one TP
485  std::vector<unsigned int> tpHits;
486  for (unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
487  if (!tp.UseHit[ii]) continue;
488  if (tp.Hits[ii] > slc.slHits.size() - 1) { break; } // bad slHits index
489  unsigned int allHitsIndex = slc.slHits[tp.Hits[ii]].allHitsIndex;
490  if (allHitsIndex > nInputHits - 1) { break; } // bad allHitsIndex
491  tpHits.push_back(allHitsIndex);
492  if (newIndex[allHitsIndex] != UINT_MAX) {
493  std::cout << "Bad Slice " << isl << " tp.Hits " << tp.Hits[ii] << " allHitsIndex "
494  << allHitsIndex;
495  std::cout << " old newIndex " << newIndex[allHitsIndex];
496  auto& oldhit = (*inputHits)[allHitsIndex];
497  std::cout << " old " << oldhit.WireID().Plane << ":" << oldhit.WireID().Wire << ":"
498  << (int)oldhit.PeakTime();
499  auto& newhit = hitCol[newIndex[allHitsIndex]];
500  std::cout << " new " << newhit.WireID().Plane << ":" << newhit.WireID().Wire << ":"
501  << (int)newhit.PeakTime();
502  std::cout << " hitCol size " << hitCol.size();
503  std::cout << "\n";
504  break;
505  }
506  } // ii
507  // Let the alg define the hit either by merging multiple hits or by a simple copy
508  // of a single hit from inputHits
509  // Merge hits in the TP that are on the same wire or create hits on multiple wires
510  // and update the old hits -> new hits assn (newIndex)
511  if (tj.AlgMod[tca::kHaloTj]) {
512  // dressed muon - don't merge hits
513  for (auto iht : tpHits) {
514  hitCol.push_back((*inputHits)[iht]);
515  newIndex[iht] = hitCol.size() - 1;
516  } // iht
517  }
518  else {
519  fTCAlg.MergeTPHits(tpHits, hitCol, newIndex);
520  }
521  } // tp
522  if (hitCol.empty()) continue;
523  // Sum the charge and make the associations
524  float sumChg = 0;
525  float sumADC = 0;
526  for (unsigned int indx = hitColBeginIndex; indx < hitCol.size(); ++indx) {
527  auto& hit = hitCol[indx];
528  sumChg += hit.Integral();
529  sumADC += hit.SummedADC();
530  if (!slices.empty() &&
531  !util::CreateAssn(evt, hitCol, slices[slcIndex], *slc_hit_assn, indx)) {
533  << "Failed to associate hits with Slice";
534  }
535  } // indx
536  geo::View_t view = hitCol[hitColBeginIndex].View();
537  auto& firstTP = tj.Pts[tj.EndPt[0]];
538  auto& lastTP = tj.Pts[tj.EndPt[1]];
539  int clsID = tj.UID;
540  if (tj.AlgMod[tca::kShowerLike]) clsID = -clsID;
541  // dressed muon - give the halo cluster the same ID as the parent
542  if (tj.AlgMod[tca::kHaloTj]) clsID = -tj.ParentID;
543  unsigned int nclhits = hitCol.size() - hitColBeginIndex + 1;
544  clsCol.emplace_back(firstTP.Pos[0], // Start wire
545  0, // sigma start wire
546  firstTP.Pos[1] / tca::tcc.unitsPerTick, // start tick
547  0, // sigma start tick
548  firstTP.AveChg, // start charge
549  firstTP.Ang, // start angle
550  0, // start opening angle (0 for line-like clusters)
551  lastTP.Pos[0], // end wire
552  0, // sigma end wire
553  lastTP.Pos[1] / tca::tcc.unitsPerTick, // end tick
554  0, // sigma end tick
555  lastTP.AveChg, // end charge
556  lastTP.Ang, // end angle
557  0, // end opening angle (0 for line-like clusters)
558  sumChg, // integral
559  0, // sigma integral
560  sumADC, // summed ADC
561  0, // sigma summed ADC
562  nclhits, // n hits
563  0, // wires over hits
564  0, // width (0 for line-like clusters)
565  clsID, // ID from TrajClusterAlg
566  view, // view
567  tca::DecodeCTP(tj.CTP), // planeID
568  recob::Cluster::Sentry // sentry
569  );
570  if (!util::CreateAssn(
571  evt, clsCol, hitCol, *cls_hit_assn, hitColBeginIndex, hitCol.size())) {
573  << "Failed to associate hits with cluster ID " << tj.UID;
574  } // exception
575  // make Slice -> cluster assn
576  if (!slices.empty()) {
577  if (!util::CreateAssn(evt, clsCol, slices[slcIndex], *slc_cls_assn)) {
579  << "Failed to associate slice with PFParticle";
580  } // exception
581  } // slices exist
582  // Make cluster -> 2V and cluster -> 3V assns
583  for (unsigned short end = 0; end < 2; ++end) {
584  if (tj.VtxID[end] <= 0) continue;
585  for (auto& vx2str : vx2StrList) {
586  if (vx2str.slIndx != isl) continue;
587  if (vx2str.ID != tj.VtxID[end]) continue;
588  if (!util::CreateAssnD(
589  evt, *cls_vx2_assn, clsCol.size() - 1, vx2str.vxColIndx, end)) {
591  << "Failed to associate cluster " << tj.UID << " with EndPoint2D";
592  } // exception
593  auto& vx2 = slc.vtxs[tj.VtxID[end] - 1];
594  if (vx2.Vx3ID > 0) {
595  for (auto vx3str : vx3StrList) {
596  if (vx3str.slIndx != isl) continue;
597  if (vx3str.ID != vx2.Vx3ID) continue;
598  if (!util::CreateAssnD(
599  evt, *cls_vx3_assn, clsCol.size() - 1, vx3str.vxColIndx, end)) {
601  << "Failed to associate cluster " << tj.UID << " with Vertex";
602  } // exception
603  break;
604  } // vx3str
605  } // vx2.Vx3ID > 0
606  break;
607  } // vx2str
608  } // end
609  } // tj (aka cluster)
610 
611  // make Showers
612  for (auto& ss3 : slc.showers) {
613  if (ss3.ID <= 0) continue;
615  shower.set_id(ss3.UID);
616  shower.set_total_energy(ss3.Energy);
617  shower.set_total_energy_err(ss3.EnergyErr);
618  shower.set_total_MIPenergy(ss3.MIPEnergy);
619  shower.set_total_MIPenergy_err(ss3.MIPEnergyErr);
620  shower.set_total_best_plane(ss3.BestPlane);
621  TVector3 dir = {ss3.Dir[0], ss3.Dir[1], ss3.Dir[2]};
622  shower.set_direction(dir);
623  TVector3 dirErr = {ss3.DirErr[0], ss3.DirErr[1], ss3.DirErr[2]};
624  shower.set_direction_err(dirErr);
625  TVector3 pos = {ss3.Start[0], ss3.Start[1], ss3.Start[2]};
626  shower.set_start_point(pos);
627  TVector3 posErr = {ss3.StartErr[0], ss3.StartErr[1], ss3.StartErr[2]};
628  shower.set_start_point_err(posErr);
629  shower.set_dedx(ss3.dEdx);
630  shower.set_dedx_err(ss3.dEdxErr);
631  shower.set_length(ss3.Len);
632  shower.set_open_angle(ss3.OpenAngle);
633  shwCol.push_back(shower);
634  // make the shower - hit association
635  std::vector<unsigned int> shwHits(ss3.Hits.size());
636  for (unsigned int iht = 0; iht < ss3.Hits.size(); ++iht)
637  shwHits[iht] = newIndex[ss3.Hits[iht]];
639  evt, *shwr_hit_assn, shwCol.size() - 1, shwHits.begin(), shwHits.end())) {
641  << "Failed to associate hits with Shower";
642  } // exception
643  } // ss3
644  } // slice isl
645 
646  // Add PFParticles now that clsCol is filled
647  for (unsigned short isl = 0; isl < nSlices; ++isl) {
648  unsigned short slcIndex = 0;
649  if (!slices.empty()) {
650  for (slcIndex = 0; slcIndex < slices.size(); ++slcIndex)
651  if (slices[slcIndex]->ID() == slcIDs[isl]) break;
652  if (slcIndex == slices.size()) continue;
653  }
654  auto& slc = fTCAlg.GetSlice(isl);
655  // See if there was a serious reconstruction failure that made the slice invalid
656  if (!slc.isValid) continue;
657  // make PFParticles
658  for (size_t ipfp = 0; ipfp < slc.pfps.size(); ++ipfp) {
659  auto& pfp = slc.pfps[ipfp];
660  if (pfp.ID <= 0) continue;
661  // parents and daughters are indexed within a slice so find the index offset in pfpCol
662  size_t self = pfpCol.size();
663  size_t offset = self - ipfp;
664  size_t parentIndex = UINT_MAX;
665  if (pfp.ParentUID > 0) parentIndex = pfp.ParentUID + offset - 1;
666  std::vector<size_t> dtrIndices(pfp.DtrUIDs.size());
667  for (unsigned short idtr = 0; idtr < pfp.DtrUIDs.size(); ++idtr)
668  dtrIndices[idtr] = pfp.DtrUIDs[idtr] + offset - 1;
669  pfpCol.emplace_back(pfp.PDGCode, self, parentIndex, dtrIndices);
670  auto pos = PosAtEnd(pfp, 0);
671  auto dir = DirAtEnd(pfp, 0);
672  double sp[] = {pos[0], pos[1], pos[2]};
673  double sd[] = {dir[0], dir[1], dir[2]};
674  double spe[] = {0., 0., 0.};
675  double sde[] = {0., 0., 0.};
676  sedCol.emplace_back(sp, sd, spe, sde);
677  // PFParticle -> clusters
678  std::vector<unsigned int> clsIndices;
679  for (auto tuid : pfp.TjUIDs) {
680  unsigned int clsIndex = 0;
681  for (clsIndex = 0; clsIndex < clsCol.size(); ++clsIndex)
682  if (abs(clsCol[clsIndex].ID()) == tuid) break;
683  if (clsIndex == clsCol.size()) continue;
684  clsIndices.push_back(clsIndex);
685  } // tjid
686  if (!util::CreateAssn(
687  evt, *pfp_cls_assn, pfpCol.size() - 1, clsIndices.begin(), clsIndices.end())) {
689  << "Failed to associate clusters with PFParticle";
690  } // exception
691  // PFParticle -> Vertex
692  if (pfp.Vx3ID[0] > 0) {
693  for (auto vx3str : vx3StrList) {
694  if (vx3str.slIndx != isl) continue;
695  if (vx3str.ID != pfp.Vx3ID[0]) continue;
696  std::vector<unsigned short> indx(1, vx3str.vxColIndx);
697  if (!util::CreateAssn(
698  evt, *pfp_vx3_assn, pfpCol.size() - 1, indx.begin(), indx.end())) {
700  << "Failed to associate PFParticle " << pfp.UID << " with Vertex";
701  } // exception
702  break;
703  } // vx3Index
704  } // start vertex exists
705  // PFParticle -> Seed
706  if (!sedCol.empty()) {
707  if (!util::CreateAssn(evt,
708  pfpCol,
709  sedCol,
710  *pfp_sed_assn,
711  sedCol.size() - 1,
712  sedCol.size(),
713  pfpCol.size() - 1)) {
715  << "Failed to associate seed with PFParticle";
716  } // exception
717  } // seeds exist
718  // PFParticle -> Slice
719  if (!slices.empty()) {
720  if (!util::CreateAssn(evt, pfpCol, slices[slcIndex], *slc_pfp_assn)) {
722  << "Failed to associate slice with PFParticle";
723  } // exception
724  } // slices exist
725  // PFParticle -> Shower
726  if (pfp.PDGCode == 1111) {
727  std::vector<unsigned short> shwIndex(1, 0);
728  for (auto& ss3 : slc.showers) {
729  if (ss3.ID <= 0) continue;
730  if (ss3.PFPIndex == ipfp) break;
731  ++shwIndex[0];
732  } // ss3
733  if (shwIndex[0] < shwCol.size()) {
734  if (!util::CreateAssn(
735  evt, *pfp_shwr_assn, pfpCol.size() - 1, shwIndex.begin(), shwIndex.end())) {
737  << "Failed to associate shower with PFParticle";
738  } // exception
739  } // valid shwIndex
740  } // pfp -> Shower
741  // PFParticle cosmic tag
742  if (tca::tcc.modes[tca::kTagCosmics]) {
743  std::vector<float> tempPt1, tempPt2;
744  tempPt1.push_back(-999);
745  tempPt1.push_back(-999);
746  tempPt1.push_back(-999);
747  tempPt2.push_back(-999);
748  tempPt2.push_back(-999);
749  tempPt2.push_back(-999);
750  ctCol.emplace_back(tempPt1, tempPt2, pfp.CosmicScore, anab::CosmicTagID_t::kNotTagged);
751  if (!util::CreateAssn(
752  evt, pfpCol, ctCol, *pfp_cos_assn, ctCol.size() - 1, ctCol.size())) {
754  << "Failed to associate CosmicTag with PFParticle";
755  }
756  } // cosmic tag
757  } // ipfp
758  } // isl
759 
760  // add the hits that weren't used in any slice to hitCol unless this is a
761  // special debugging mode and would be a waste of time
762  if (!slices.empty() && tca::tcc.recoSlice == 0) {
763  auto inputSlices = evt.getValidHandle<std::vector<recob::Slice>>(fSliceModuleLabel);
764  art::FindManyP<recob::Hit> hitFromSlc(inputSlices, evt, fSliceModuleLabel);
765  for (unsigned int allHitsIndex = 0; allHitsIndex < nInputHits; ++allHitsIndex) {
766  if (newIndex[allHitsIndex] != UINT_MAX) continue;
767  std::vector<unsigned int> oneHit(1, allHitsIndex);
768  fTCAlg.MergeTPHits(oneHit, hitCol, newIndex);
769  // find out which slice it is in
770  bool gotit = false;
771  for (size_t isl = 0; isl < slices.size(); ++isl) {
772  auto& hit_in_slc = hitFromSlc.at(isl);
773  for (auto& hit : hit_in_slc) {
774  if (hit.key() != allHitsIndex) continue;
775  gotit = true;
776  // Slice -> Hit assn
777  if (!util::CreateAssn(evt, hitCol, slices[isl], *slc_hit_assn)) {
779  << "Failed to associate old Hit with Slice";
780  } // exception
781  break;
782  } // hit
783  if (gotit) break;
784  } // isl
785  } // allHitsIndex
786  } // slices exist
787  else {
788  // no recob::Slices. Just copy the unused hits
789  for (unsigned int allHitsIndex = 0; allHitsIndex < nInputHits; ++allHitsIndex) {
790  if (newIndex[allHitsIndex] != UINT_MAX) continue;
791  std::vector<unsigned int> oneHit(1, allHitsIndex);
792  fTCAlg.MergeTPHits(oneHit, hitCol, newIndex);
793  } // allHitsIndex
794  } // recob::Slices
795  } // input hits exist
796 
797  // www: find spacepoint from hits (inputHits) through SpacePoint->Hit assns, then create association between spacepoint and trajcluster hits (here, hits in hitCol)
798  if (nInputHits > 0) {
799  // www: expecting to find spacepoint from hits (inputHits): SpacePoint->Hit assns
800  if (fSpacePointModuleLabel != "NA") {
802  // www: using sp from hit
803  for (unsigned int allHitsIndex = 0; allHitsIndex < nInputHits; ++allHitsIndex) {
804  if (newIndex[allHitsIndex] == UINT_MAX)
805  continue; // skip hits not used in slice (not TrajCluster hits)
806  auto& sp_from_hit = spFromHit.at(allHitsIndex);
807  for (auto& sp : sp_from_hit) {
808  // SpacePoint -> Hit assn
809  if (!util::CreateAssn(evt, hitCol, sp, *sp_hit_assn, newIndex[allHitsIndex])) {
811  << "Failed to associate new Hit with SpacePoint";
812  } // exception
813  } // sp
814  } // allHitsIndex
815  } // fSpacePointModuleLabel != "NA"
816  } // nInputHits > 0
817 
818  // clear the alg data structures
820 
821  // convert vectors to unique_ptrs
822  std::unique_ptr<std::vector<recob::Hit>> hcol(new std::vector<recob::Hit>(std::move(hitCol)));
823  std::unique_ptr<std::vector<recob::Cluster>> ccol(
824  new std::vector<recob::Cluster>(std::move(clsCol)));
825  std::unique_ptr<std::vector<recob::EndPoint2D>> v2col(
826  new std::vector<recob::EndPoint2D>(std::move(vx2Col)));
827  std::unique_ptr<std::vector<recob::Vertex>> v3col(
828  new std::vector<recob::Vertex>(std::move(vx3Col)));
829  std::unique_ptr<std::vector<recob::PFParticle>> pcol(
830  new std::vector<recob::PFParticle>(std::move(pfpCol)));
831  std::unique_ptr<std::vector<recob::Seed>> sdcol(
832  new std::vector<recob::Seed>(std::move(sedCol)));
833  std::unique_ptr<std::vector<recob::Shower>> scol(
834  new std::vector<recob::Shower>(std::move(shwCol)));
835  std::unique_ptr<std::vector<anab::CosmicTag>> ctgcol(
836  new std::vector<anab::CosmicTag>(std::move(ctCol)));
837 
838  // move the cluster collection and the associations into the event:
839  if (fHitModuleLabel != "NA") {
841  shcol.use_hits(std::move(hcol));
842  shcol.put_into(evt);
843  }
844  else {
846  shcol.use_hits(std::move(hcol));
847  shcol.put_into(evt);
848  }
849  evt.put(std::move(ccol));
850  evt.put(std::move(cls_hit_assn));
851  evt.put(std::move(v2col));
852  evt.put(std::move(v3col));
853  evt.put(std::move(scol));
854  evt.put(std::move(sdcol));
855  evt.put(std::move(shwr_hit_assn));
856  evt.put(std::move(cls_vx2_assn));
857  evt.put(std::move(cls_vx3_assn));
858  evt.put(std::move(pcol));
859  evt.put(std::move(pfp_cls_assn));
860  evt.put(std::move(pfp_shwr_assn));
861  evt.put(std::move(pfp_vx3_assn));
862  evt.put(std::move(pfp_sed_assn));
863  evt.put(std::move(slc_cls_assn));
864  evt.put(std::move(slc_pfp_assn));
865  evt.put(std::move(slc_hit_assn));
866  evt.put(std::move(ctgcol));
867  evt.put(std::move(pfp_cos_assn));
868  evt.put(std::move(sp_hit_assn)); // www: association between sp and hit (trjaclust)
869  } // TrajCluster::produce()
870 
872  void TrajCluster::GetHits(const std::vector<recob::Hit>& inputHits,
873  const geo::TPCID& tpcid,
874  std::vector<std::vector<unsigned int>>& tpcHits)
875  {
876  // Put hits in this TPC into a single "slice", unless a special debugging mode is specified to
877  // only reconstruct hits that are MC-matched
878  unsigned int tpc = tpcid.TPC;
879  tpcHits.resize(1);
880  for (size_t iht = 0; iht < inputHits.size(); ++iht) {
881  auto& hit = inputHits[iht];
882  if (hit.WireID().TPC == tpc) tpcHits[0].push_back(iht);
883  }
884  } // GetHits
885 
887  void TrajCluster::GetHits(const geo::TPCID& tpcid,
888  const std::vector<recob::Slice>& inputSlices,
889  art::FindManyP<recob::Hit>& hitFromSlc,
890  std::vector<std::vector<unsigned int>>& tpcHits,
891  std::vector<int>& slcIDs)
892  {
893  // Put the hits in all slices into tpcHits in this TPC
894  tpcHits.clear();
895  slcIDs.clear();
896  if (!hitFromSlc.isValid()) return;
897 
898  unsigned int tpc = tpcid.TPC;
899 
900  for (size_t isl = 0; isl < inputSlices.size(); ++isl) {
901  auto& hit_in_slc = hitFromSlc.at(isl);
902  if (hit_in_slc.size() < 3) continue;
903  int slcID = inputSlices[isl].ID();
904  for (auto& hit : hit_in_slc) {
905  if (hit->WireID().TPC != tpc) continue;
906  unsigned short indx = 0;
907  for (indx = 0; indx < slcIDs.size(); ++indx)
908  if (slcID == slcIDs[indx]) break;
909  if (indx == slcIDs.size()) {
910  slcIDs.push_back(slcID);
911  tpcHits.resize(tpcHits.size() + 1);
912  }
913  tpcHits[indx].push_back(hit.key());
914  } // hit
915  } // isl
916 
917  } // GetHits
918 
919  //----------------------------------------------------------------------------
921 
922 } // namespace cluster
void PrintAll(detinfo::DetectorPropertiesData const &detProp, std::string someText)
Definition: Utils.cxx:5367
void ClearResults()
Deletes all the results.
bool CreateAssnD(art::Event &evt, art::Assns< T, U, D > &assn, size_t first_index, size_t second_index, typename art::Assns< T, U, D >::data_t &&data)
Creates a single one-to-one association with associated data.
void set_start_point_err(const TVector3 &xyz_e)
Definition: Shower.h:137
void set_dedx_err(const std::vector< double > &q)
Definition: Shower.h:139
void set_direction_err(const TVector3 &dir_e)
Definition: Shower.h:135
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
Utilities related to art service access.
tca::TrajClusterAlg fTCAlg
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
TCConfig tcc
Definition: DataStructs.cxx:9
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
void set_total_energy(const std::vector< double > &q)
Definition: Shower.h:128
void RunTrajClusterAlg(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
art::InputTag fSliceModuleLabel
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
void set_total_MIPenergy_err(const std::vector< double > &q)
Definition: Shower.h:131
Float_t tmp
Definition: plot.C:35
Point3_t PosAtEnd(const PFPStruct &pfp, unsigned short end)
Definition: PFPUtils.cxx:3259
Cluster finding and building.
std::vector< std::string > const & GetAlgBitNames() const
bool SortHits(HitLoc const &h1, HitLoc const &h2)
std::vector< std::array< unsigned int, 3 > > sptHits
SpacePoint -> Hits assns by plane.
Definition: DataStructs.h:626
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248
void set_total_energy_err(const std::vector< double > &q)
Definition: Shower.h:129
static const SentryArgument_t Sentry
An instance of the sentry object.
Definition: Cluster.h:174
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
void set_id(const int id)
Definition: Shower.h:127
bool SetInputHits(std::vector< recob::Hit > const &inputHits, unsigned int run, unsigned int event)
void DefineShTree(TTree *t)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::string const & label() const noexcept
Definition: InputTag.cc:79
unsigned int Hit
set to the hit index in evt.allHits if a Plane:Wire:Tick match is found
Definition: DebugStruct.h:26
Helper functions to create a hit.
art::InputTag fHitModuleLabel
void hits()
Definition: readHits.C:15
void set_direction(const TVector3 &dir)
Definition: Shower.h:134
void use_hits(std::unique_ptr< std::vector< recob::Hit >> &&srchits)
Uses the specified collection as data product.
Definition: HitCreator.cxx:508
if(nlines<=0)
int Wire
Select hit Wire for debugging.
Definition: DebugStruct.h:24
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
float unitsPerTick
scale factor from Tick to WSE equivalent units
Definition: DataStructs.h:564
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
DebugStuff debug
Definition: DebugStruct.cxx:4
void put_into(art::Event &)
Moves the data into the event.
Definition: HitCreator.h:850
unsigned short GetSlicesSize() const
void set_length(const double &l)
Definition: Shower.h:140
int Plane
Select plane.
Definition: DebugStruct.h:22
void set_open_angle(const double &a)
Definition: Shower.h:141
A class handling a collection of hits and its associations.
Definition: HitCreator.h:787
EventNumber_t event() const
Definition: Event.cc:41
const geo::GeometryCore * geom
Definition: DataStructs.h:569
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Declaration of cluster object.
void set_total_best_plane(const int q)
Definition: Shower.h:132
Definition of data types for geometry description.
void set_total_MIPenergy(const std::vector< double > &q)
Definition: Shower.h:130
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
Detector simulation of raw signals on wires.
TCSlice const & GetSlice(unsigned short sliceIndex) const
art::InputTag fSpacePointModuleLabel
TH1F * h2
Definition: plot.C:44
ProducesCollector & producesCollector() noexcept
void GetHits(const std::vector< recob::Hit > &inputHits, const geo::TPCID &tpcid, std::vector< std::vector< unsigned int >> &tpcHits)
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.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
int Tick
Select hit PeakTime for debugging (< 0 for vertex finding)
Definition: DebugStruct.h:25
int TPC
Select TPC.
Definition: DebugStruct.h:21
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
TDirectory * dir
Definition: macro.C:5
geo::PlaneID DecodeCTP(CTP_t CTP)
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Produces clusters by the TrajCluster algorithm.
void produce(art::Event &evt) override
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:51
TH1F * h1
Definition: plot.C:41
void MergeTPHits(std::vector< unsigned int > &tpHits, std::vector< recob::Hit > &newHitCol, std::vector< unsigned int > &newHitAssns) const
void SetSourceHits(std::vector< recob::Hit > const &srcHits)
std::vector< PFPStruct > pfps
Definition: DataStructs.h:671
TCEvent evt
Definition: DataStructs.cxx:8
void set_start_point(const TVector3 &xyz)
Definition: Shower.h:136
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
void SetInputSpts(std::vector< recob::SpacePoint > const &sptHandle)
RunNumber_t run() const
Definition: Event.cc:29
Vector3_t DirAtEnd(const PFPStruct &pfp, unsigned short end)
Definition: PFPUtils.cxx:3251
short recoSlice
only reconstruct the slice with ID (0 = all)
Definition: DataStructs.h:581
tag cosmic rays
Definition: DataStructs.h:532
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
art::InputTag fSpacePointHitAssnLabel
void set_dedx(const std::vector< double > &q)
Definition: Shower.h:138
std::vector< unsigned int > const & GetAlgModCount() const
TrajCluster(fhicl::ParameterSet const &pset)