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