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