LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
EMShower3D_module.cc
Go to the documentation of this file.
1 // Class: EMShower3D
3 // Module Type: producer
4 // File: EMShower3D_module.cc
5 // Authors: dorota.stefan@cern.ch robert.sulej@cern.ch
7 
13 #include "fhiclcpp/ParameterSet.h"
14 
24 
28 
29 #include <memory>
30 
33 
34 namespace {
35  struct IniSeg {
36  size_t idcl1;
37  size_t idcl2;
38  size_t idcl3;
39  size_t view1;
40  size_t view2;
41  size_t view3;
43  std::vector<art::Ptr<recob::Hit>> hits1;
44  std::vector<art::Ptr<recob::Hit>> hits2;
45  std::vector<art::Ptr<recob::Hit>> hits3;
46  };
47 }
48 
49 namespace ems {
50  class EMShower3D;
51 }
52 
54 public:
55  explicit EMShower3D(fhicl::ParameterSet const& p);
56 
57  EMShower3D(EMShower3D const&) = delete;
58  EMShower3D(EMShower3D&&) = delete;
59  EMShower3D& operator=(EMShower3D const&) = delete;
60  EMShower3D& operator=(EMShower3D&&) = delete;
61 
62 private:
63  void produce(art::Event& e) override;
64 
65  recob::Track ConvertFrom(detinfo::DetectorClocksData const& clock_data,
66  detinfo::DetectorPropertiesData const& det_prop,
67  pma::Track3D& src);
68  recob::Track ConvertFrom2(detinfo::DetectorClocksData const& clock_data,
69  detinfo::DetectorPropertiesData const& det_prop,
70  pma::Track3D& src);
71  recob::Cluster ConvertFrom(const std::vector<art::Ptr<recob::Hit>>& src);
72 
73  std::vector<ems::DirOfGamma*> CollectShower2D(detinfo::DetectorPropertiesData const& detProp,
74  art::Event const& e);
75 
76  void Link(detinfo::DetectorPropertiesData const& detProp, std::vector<ems::DirOfGamma*> input);
77 
78  // Remove tracks which are too close to each other
79  void Reoptimize(detinfo::DetectorPropertiesData const& detProp);
80 
81  void Make3DSeg(detinfo::DetectorPropertiesData const& detProp,
82  std::vector<ems::DirOfGamma*> pair);
83 
84  bool Validate(detinfo::DetectorPropertiesData const& detProp,
85  std::vector<ems::DirOfGamma*> input,
86  size_t id1,
87  size_t id2,
88  size_t c1,
89  size_t c2,
90  size_t plane3);
91 
92  void FilterOutSmallParts(detinfo::DetectorPropertiesData const& detProp,
93  double r2d,
94  const std::vector<art::Ptr<recob::Hit>>& hits_in,
96 
97  bool GetCloseHits(detinfo::DetectorPropertiesData const& detProp,
98  double r2d,
99  const std::vector<art::Ptr<recob::Hit>>& hits_in,
100  std::vector<size_t>& used,
101  std::vector<art::Ptr<recob::Hit>>& hits_out);
102 
103  bool Has(const std::vector<size_t>& v, size_t idx);
104 
105  size_t LinkCandidates(detinfo::DetectorPropertiesData const& detProp,
106  std::vector<ems::DirOfGamma*> input,
107  size_t id);
108 
109  std::vector<IniSeg> fInisegs;
110  std::vector<IniSeg> fSeltracks;
111  std::vector<IniSeg> fPMA3D;
112 
113  std::vector<std::vector<art::Ptr<recob::Hit>>> fClusters;
114 
115  std::vector<size_t> fClustersNotUsed;
116  std::vector<size_t> fTracksNotUsed;
117 
118  unsigned int fTrkIndex;
119  unsigned int fClIndex;
120  unsigned int fIniIndex;
121 
122  std::string fCluModuleLabel;
123  std::string fTrk3DModuleLabel;
124 
127 
129 };
130 
132  : EDProducer{p}
133  , fProjectionMatchingAlg(p.get<fhicl::ParameterSet>("ProjectionMatchingAlg"))
134  , fCalorimetryAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
135 {
136  fCluModuleLabel = p.get<std::string>("ClustersModuleLabel");
137  fTrk3DModuleLabel = p.get<std::string>("Trk3DModuleLabel");
138 
139  produces<std::vector<recob::Track>>();
140  produces<std::vector<recob::Vertex>>();
141  produces<std::vector<recob::Cluster>>();
142  produces<std::vector<recob::SpacePoint>>();
143  produces<art::Assns<recob::Track, recob::Hit>>();
144  produces<art::Assns<recob::Track, recob::Vertex>>();
145  produces<art::Assns<recob::Cluster, recob::Hit>>();
146  produces<art::Assns<recob::Track, recob::SpacePoint>>();
147  produces<art::Assns<recob::SpacePoint, recob::Hit>>();
148  produces<art::Assns<recob::Track, recob::Cluster>>();
149 }
150 
152 {
153  return recob::Cluster(0.0F,
154  0.0F,
155  0.0F,
156  0.0F,
157  0.0F,
158  0.0F,
159  0.0F,
160  0.0F,
161  0.0F,
162  0.0F,
163  0.0F,
164  0.0F,
165  0.0F,
166  0.0F,
167  0.0F,
168  0.0F,
169  0.0F,
170  0.0F,
171  src.size(),
172  0.0F,
173  0.0F,
174  fClIndex,
175  src[0]->View(),
176  src[0]->WireID().planeID());
177 }
178 
180  detinfo::DetectorPropertiesData const& detProp,
181  pma::Track3D& src)
182 {
183  auto const* geom = lar::providerFrom<geo::Geometry>();
184  double avdrift = (src.front()->Point3D().X() + src.back()->Point3D().X()) * 0.5;
185  unsigned int nplanes = geom->Nplanes({src.front()->Cryo(), src.front()->TPC()});
186  size_t nusedhitsmax = 0;
187  int bestplane = -1;
188  for (unsigned int p = 0; p < nplanes; ++p) {
189  unsigned int nusedP = 0;
190  fProjectionMatchingAlg.selectInitialHits(src, p, &nusedP);
191 
192  if (nusedP > nusedhitsmax) {
193  nusedhitsmax = nusedP;
194  bestplane = int(p);
195  }
196  }
197 
198  std::vector<std::vector<double>> vdedx;
199  std::vector<double> dedx;
200 
201  for (unsigned int p = 0; p < nplanes; ++p) {
202  unsigned int nusedP = 0;
203  double dqdxplane = fProjectionMatchingAlg.selectInitialHits(src, p, &nusedP);
204  double timeP = detProp.ConvertXToTicks(avdrift, p, src.front()->TPC(), src.front()->Cryo());
205  double dEdxplane = fCalorimetryAlg.dEdx_AREA(clock_data, detProp, dqdxplane, timeP, p);
206  dedx.push_back(dEdxplane);
207  if (int(p) == bestplane)
208  dedx.push_back(1);
209  else
210  dedx.push_back(0);
211  vdedx.push_back(dedx);
212  }
213 
214  std::vector<TVector3> xyz, dircos;
215 
216  for (size_t i = 0; i < src.size(); i++) {
217  xyz.push_back(src[i]->Point3D());
218 
219  if (i < src.size() - 1) {
220  TVector3 dc(src[i + 1]->Point3D());
221  dc -= src[i]->Point3D();
222  dc *= 1.0 / dc.Mag();
223  dircos.push_back(dc);
224  }
225  else
226  dircos.push_back(dircos.back());
227  }
228 
231  recob::Track::Flags_t(xyz.size()),
232  false),
233  0,
234  -1.,
235  0,
238  fIniIndex);
239 }
240 
242  detinfo::DetectorPropertiesData const& detProp,
243  pma::Track3D& src)
244 {
245  auto const* geom = lar::providerFrom<geo::Geometry>();
246 
247  double avdrift = (src.front()->Point3D().X() + src.back()->Point3D().X()) * 0.5;
248  unsigned int nplanes = geom->Nplanes({src.front()->Cryo(), src.front()->TPC()});
249  size_t nusedhitsmax = 0;
250  int bestplane = -1;
251  for (unsigned int p = 0; p < nplanes; ++p) {
252  unsigned int nusedP = 0;
253  fProjectionMatchingAlg.selectInitialHits(src, p, &nusedP);
254 
255  if (nusedP > nusedhitsmax) {
256  nusedhitsmax = nusedP;
257  bestplane = int(p);
258  }
259  }
260 
261  std::vector<std::vector<double>> vdedx;
262  std::vector<double> dedx;
263 
264  for (unsigned int p = 0; p < nplanes; ++p) {
265  unsigned int nusedP = 0;
266  double dqdxplane = fProjectionMatchingAlg.selectInitialHits(src, p, &nusedP);
267  double timeP = detProp.ConvertXToTicks(avdrift, p, src.front()->TPC(), src.front()->Cryo());
268  double dEdxplane = fCalorimetryAlg.dEdx_AREA(clockData, detProp, dqdxplane, timeP, p);
269  dedx.push_back(dEdxplane);
270  if (int(p) == bestplane)
271  dedx.push_back(1);
272  else
273  dedx.push_back(0);
274  vdedx.push_back(dedx);
275  }
276 
277  std::vector<TVector3> xyz, dircos;
278 
279  for (size_t i = 0; i < src.size(); i++) {
280  xyz.push_back(src[i]->Point3D());
281 
282  if (i < src.size() - 1) {
283  TVector3 dc(src[i + 1]->Point3D());
284  dc -= src[i]->Point3D();
285  dc *= 1.0 / dc.Mag();
286  dircos.push_back(dc);
287  }
288  else
289  dircos.push_back(dircos.back());
290  }
291 
294  recob::Track::Flags_t(xyz.size()),
295  false),
296  0,
297  -1.,
298  0,
301  fIniIndex);
302 }
303 
305 {
307  fSeltracks.clear();
308  fInisegs.clear();
309  fClusters.clear();
310  fPMA3D.clear();
311  fClustersNotUsed.clear();
312 
313  auto tracks = std::make_unique<std::vector<recob::Track>>();
314  auto vertices = std::make_unique<std::vector<recob::Vertex>>();
315  auto clusters = std::make_unique<std::vector<recob::Cluster>>();
316  auto allsp = std::make_unique<std::vector<recob::SpacePoint>>();
317 
318  auto trk2hit = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
319  auto trk2vtx = std::make_unique<art::Assns<recob::Track, recob::Vertex>>();
320  auto cl2hit = std::make_unique<art::Assns<recob::Cluster, recob::Hit>>();
321  auto trk2cl = std::make_unique<art::Assns<recob::Track, recob::Cluster>>();
322  auto trk2sp = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
323  auto sp2hit = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
324 
325  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
326  auto const detProp =
328 
330  fClustersNotUsed.clear();
331  fInisegs.clear();
333 
334  for (size_t id = 0; id < fCluListHandle->size(); id++) {
335  std::vector<art::Ptr<recob::Hit>> hitlist;
336  hitlist = fb.at(id);
337 
338  if (hitlist.size() > 5) fClustersNotUsed.push_back(id);
339  }
340 
341  std::vector<ems::DirOfGamma*> showernviews = CollectShower2D(detProp, e);
342 
343  Link(detProp, showernviews);
344 
345  while (fInisegs.size()) {
346  fSeltracks.push_back(fInisegs[0]);
347  fInisegs.erase(fInisegs.begin() + 0);
348  }
349 
350  Reoptimize(detProp);
351 
352  // conversion from pma track to recob::track
353 
354  size_t spStart = 0, spEnd = 0;
355  double sp_pos[3], sp_err[6], vtx_pos[3];
356  for (size_t i = 0; i < 6; i++)
357  sp_err[i] = 1.0;
358 
359  fTrkIndex = 0;
360 
361  for (auto const& trk : fSeltracks) {
362  tracks->push_back(ConvertFrom(clockData, detProp, *(trk.track)));
363 
364  vtx_pos[0] = trk.track->front()->Point3D().X();
365  vtx_pos[1] = trk.track->front()->Point3D().Y();
366  vtx_pos[2] = trk.track->front()->Point3D().Z();
367  vertices->emplace_back(vtx_pos, fTrkIndex);
368 
369  ++fTrkIndex;
370 
371  std::vector<art::Ptr<recob::Cluster>> cl2d;
372  cl2d.emplace_back(fCluListHandle, trk.idcl1);
373  cl2d.emplace_back(fCluListHandle, trk.idcl2);
374 
375  std::vector<art::Ptr<recob::Hit>> hits2d;
377 
378  spStart = allsp->size();
379  for (int h = trk.track->size() - 1; h >= 0; h--) {
380  pma::Hit3D* h3d = (*trk.track)[h];
381  hits2d.push_back(h3d->Hit2DPtr());
382 
383  if ((h == 0) || (sp_pos[0] != h3d->Point3D().X()) || (sp_pos[1] != h3d->Point3D().Y()) ||
384  (sp_pos[2] != h3d->Point3D().Z())) {
385  if (sp_hits.size()) // hits assigned to the previous sp
386  {
387  util::CreateAssn(e, *allsp, sp_hits, *sp2hit);
388  sp_hits.clear();
389  }
390  sp_pos[0] = h3d->Point3D().X();
391  sp_pos[1] = h3d->Point3D().Y();
392  sp_pos[2] = h3d->Point3D().Z();
393  allsp->push_back(recob::SpacePoint(sp_pos, sp_err, 1.0));
394  }
395  sp_hits.push_back(h3d->Hit2DPtr());
396  }
397  if (sp_hits.size()) // hits assigned to the last sp
398  {
399  util::CreateAssn(e, *allsp, sp_hits, *sp2hit);
400  }
401  spEnd = allsp->size();
402 
403  if (vertices->size()) {
404  size_t vtx_idx = (size_t)(vertices->size() - 1);
405  util::CreateAssn(e, *tracks, *vertices, *trk2vtx, vtx_idx, vtx_idx + 1);
406  }
407 
408  if (cl2d.size()) { util::CreateAssn(e, *tracks, cl2d, *trk2cl); }
409 
410  if (hits2d.size()) {
411  util::CreateAssn(e, *tracks, *allsp, *trk2sp, spStart, spEnd);
412  util::CreateAssn(e, *tracks, hits2d, *trk2hit);
413  }
414  }
415 
416  fIniIndex = fTrkIndex + 1;
417  for (auto const& trk : fPMA3D) {
418  tracks->push_back(ConvertFrom2(clockData, detProp, *(trk.track)));
419 
420  fIniIndex++;
421 
422  std::vector<art::Ptr<recob::Cluster>> cl2d;
423  cl2d.push_back(art::Ptr<recob::Cluster>(fCluListHandle, trk.idcl1));
424  cl2d.push_back(art::Ptr<recob::Cluster>(fCluListHandle, trk.idcl2));
425 
426  std::vector<art::Ptr<recob::Hit>> hits2d;
428 
429  spStart = allsp->size();
430  for (int h = trk.track->size() - 1; h >= 0; h--) {
431  pma::Hit3D* h3d = (*trk.track)[h];
432  hits2d.push_back(h3d->Hit2DPtr());
433 
434  if ((h == 0) || (sp_pos[0] != h3d->Point3D().X()) || (sp_pos[1] != h3d->Point3D().Y()) ||
435  (sp_pos[2] != h3d->Point3D().Z())) {
436  if (sp_hits.size()) // hits assigned to the previous sp
437  {
438  util::CreateAssn(e, *allsp, sp_hits, *sp2hit);
439  sp_hits.clear();
440  }
441  sp_pos[0] = h3d->Point3D().X();
442  sp_pos[1] = h3d->Point3D().Y();
443  sp_pos[2] = h3d->Point3D().Z();
444  allsp->push_back(recob::SpacePoint(sp_pos, sp_err, 1.0));
445  }
446  sp_hits.push_back(h3d->Hit2DPtr());
447  }
448  if (sp_hits.size()) // hits assigned to the last sp
449  {
450  util::CreateAssn(e, *allsp, sp_hits, *sp2hit);
451  }
452  spEnd = allsp->size();
453 
454  if (cl2d.size()) { util::CreateAssn(e, *tracks, cl2d, *trk2cl); }
455 
456  if (hits2d.size()) {
457  util::CreateAssn(e, *tracks, *allsp, *trk2sp, spStart, spEnd);
458  util::CreateAssn(e, *tracks, hits2d, *trk2hit);
459  }
460  }
461 
462  // create cluster from hits, which were an input to find initial part of the cascade.
463  fClIndex = 0;
464  for (auto const& cl : fClusters)
465  if (cl.size()) {
466  clusters->push_back(ConvertFrom(cl));
467  fClIndex++;
468 
469  util::CreateAssn(e, *clusters, cl, *cl2hit);
470  }
471 
472  for (unsigned int i = 0; i < showernviews.size(); i++)
473  delete showernviews[i];
474 
475  for (unsigned int i = 0; i < fSeltracks.size(); i++)
476  delete fSeltracks[i].track;
477 
478  for (unsigned int i = 0; i < fInisegs.size(); i++)
479  delete fInisegs[i].track;
480 
481  for (unsigned int i = 0; i < fPMA3D.size(); i++)
482  delete fPMA3D[i].track;
483  }
484 
485  e.put(std::move(tracks));
486  e.put(std::move(vertices));
487  e.put(std::move(clusters));
488  e.put(std::move(allsp));
489 
490  e.put(std::move(trk2hit));
491  e.put(std::move(trk2vtx));
492  e.put(std::move(cl2hit));
493  e.put(std::move(trk2cl));
494  e.put(std::move(trk2sp));
495  e.put(std::move(sp2hit));
496 }
497 
499 {
500  if (empty(fSeltracks)) return;
501  const float min_dist = 3.0F;
502  size_t ta = 0;
503  while (ta < (fSeltracks.size() - 1)) {
504  size_t tb = ta + 1;
505  bool found = false;
506  while (tb < fSeltracks.size()) {
507  if (ta == tb) {
508  tb++;
509  continue;
510  }
511 
512  TVector3 p1 = fSeltracks[ta].track->front()->Point3D();
513  TVector3 p2 = fSeltracks[tb].track->front()->Point3D();
514  float dist = std::sqrt(pma::Dist2(p1, p2));
515 
516  if (dist < min_dist)
517  if ((fSeltracks[ta].idcl1 == fSeltracks[tb].idcl1) ||
518  (fSeltracks[ta].idcl1 == fSeltracks[tb].idcl2) ||
519  (fSeltracks[ta].idcl2 == fSeltracks[tb].idcl1) ||
520  (fSeltracks[ta].idcl2 == fSeltracks[tb].idcl2)) {
521  found = true;
522  size_t view3 = fSeltracks[ta].view1;
523  size_t idcl3 = fSeltracks[ta].idcl1;
524  std::vector<art::Ptr<recob::Hit>> hits3 = fSeltracks[ta].hits1;
525  std::vector<art::Ptr<recob::Hit>> hits = fSeltracks[ta].hits1;
526  for (size_t h = 0; h < fSeltracks[ta].hits2.size(); ++h)
527  hits.push_back(fSeltracks[ta].hits2[h]);
528 
529  if ((fSeltracks[tb].view1 != fSeltracks[ta].view1) &&
530  (fSeltracks[tb].view1 != fSeltracks[ta].view2)) {
531  view3 = fSeltracks[tb].view1;
532  for (size_t h = 0; h < fSeltracks[tb].hits1.size(); ++h)
533  hits.push_back(fSeltracks[tb].hits1[h]);
534  }
535  if ((fSeltracks[tb].view2 != fSeltracks[ta].view1) &&
536  (fSeltracks[tb].view2 != fSeltracks[ta].view2)) {
537  view3 = fSeltracks[tb].view2;
538  for (size_t h = 0; h < fSeltracks[tb].hits2.size(); ++h)
539  hits.push_back(fSeltracks[tb].hits2[h]);
540  }
541 
542  if ((view3 == fSeltracks[ta].view1) || (view3 == fSeltracks[ta].view2)) {
543  delete fSeltracks[ta].track;
544  fSeltracks.erase(fSeltracks.begin() + ta);
545  }
546  else {
548 
549  if (pma::Dist2(track->back()->Point3D(), fSeltracks[ta].track->front()->Point3D()) <
550  pma::Dist2(track->front()->Point3D(), fSeltracks[ta].track->front()->Point3D()))
551  track->Flip();
552 
553  IniSeg initrack;
554  initrack.idcl1 = fSeltracks[ta].idcl1;
555  initrack.idcl3 = idcl3;
556  initrack.view1 = fSeltracks[ta].view1;
557  initrack.view3 = view3;
558  initrack.hits1 = fSeltracks[ta].hits1;
559  initrack.hits3 = hits3;
560  initrack.idcl2 = fSeltracks[ta].idcl2;
561  initrack.view2 = fSeltracks[ta].view2;
562  initrack.hits2 = fSeltracks[ta].hits2;
563  initrack.track = track;
564 
565  delete fSeltracks[tb].track;
566  delete fSeltracks[ta].track;
567  fSeltracks.erase(fSeltracks.begin() + tb);
568  fSeltracks.erase(fSeltracks.begin() + ta);
569  fSeltracks.push_back(initrack);
570  }
571  }
572 
573  if (found) break;
574  tb++;
575  }
576 
577  if (!found) ta++;
578  }
579 }
580 
581 std::vector<ems::DirOfGamma*> ems::EMShower3D::CollectShower2D(
582  detinfo::DetectorPropertiesData const& detProp,
583  art::Event const& e)
584 {
585  std::vector<ems::DirOfGamma*> input;
586 
589  for (unsigned int c = 0; c < fCluListHandle->size(); c++) {
590  std::vector<art::Ptr<recob::Hit>> hitlist;
591  hitlist = fb.at(c);
592 
593  if (hitlist.size() > 5) {
594  std::vector<art::Ptr<recob::Hit>> hits_out;
595  FilterOutSmallParts(detProp, 2.0, hitlist, hits_out);
596 
597  if (hits_out.size() > 5) {
598  fClusters.push_back(hits_out);
599 
600  ems::DirOfGamma* sh = new ems::DirOfGamma(detProp, hits_out, 14, c);
601 
602  if (sh->GetHits2D().size()) input.push_back(sh);
603  }
604  }
605  }
606  }
607 
608  return input;
609 }
610 
612  std::vector<ems::DirOfGamma*> input)
613 {
614  std::vector<std::vector<size_t>> saveids;
615  std::vector<size_t> saveidsnotusedcls;
616  size_t i = 0;
617 
618  while (i < input.size()) {
619  if (!input[i]->GetCandidates().size()) {
620  i++;
621  continue;
622  }
623 
624  double mindist = 1.0; // cm
625  std::vector<ems::DirOfGamma*> pairs;
626 
627  size_t startview = input[i]->GetFirstHit()->WireID().Plane;
628  size_t tpc = input[i]->GetFirstHit()->WireID().TPC;
629  size_t cryo = input[i]->GetFirstHit()->WireID().Cryostat;
630 
631  float t1 = detProp.ConvertTicksToX(input[i]->GetFirstHit()->PeakTime(), startview, tpc, cryo);
632 
633  unsigned int idsave = 0;
634  for (unsigned int j = 0; j < input.size(); j++) {
635  if (!input[j]->GetCandidates().size()) continue;
636 
637  size_t secondview = input[j]->GetFirstHit()->WireID().Plane;
638  size_t tpc_j = input[j]->GetFirstHit()->WireID().TPC;
639  size_t cryo_j = input[j]->GetFirstHit()->WireID().Cryostat;
640 
641  if ((i != j) && (secondview != startview) && (tpc == tpc_j) && (cryo == cryo_j)) {
642  float t2 =
643  detProp.ConvertTicksToX(input[j]->GetFirstHit()->PeakTime(), secondview, tpc_j, cryo_j);
644  float dist = fabs(t2 - t1);
645 
646  if (dist < mindist) {
647  mindist = dist;
648  pairs.clear();
649  pairs.push_back(input[i]);
650  pairs.push_back(input[j]);
651  idsave = j;
652  }
653  }
654  }
655 
656  bool exist = false;
657  for (unsigned int v = 0; v < saveids.size(); v++)
658  if ((saveids[v][0] == i) || (saveids[v][0] == idsave))
659  if ((saveids[v][1] == i) || (saveids[v][1] == idsave)) exist = true;
660 
661  if (pairs.size()) {
662  if (!exist) Make3DSeg(detProp, pairs);
663 
664  std::vector<size_t> ids;
665  ids.push_back(i);
666  ids.push_back(idsave);
667  saveids.push_back(ids);
668  }
669  else {
670  saveidsnotusedcls.push_back(i);
671  }
672 
673  i++;
674  }
675 
676  i = 0;
677  while (i < saveidsnotusedcls.size()) {
678  LinkCandidates(detProp, input, i);
679  i++;
680  }
681 }
682 
684  std::vector<ems::DirOfGamma*> input,
685  size_t id)
686 {
688 
689  size_t index = id;
690  bool found = false;
691 
692  if (input[id]->GetCandidates().size() < 2) { return index; }
693 
694  double mindist = 3.0; // cm
695  std::vector<ems::DirOfGamma*> pairs;
696 
697  size_t idcsave = 0;
698  size_t idcjsave = 0;
699  size_t c = 0;
700  size_t idsave = 0;
701  while (c < input[id]->GetCandidates().size()) {
702 
703  size_t startview = input[id]->GetCandidates()[c].GetPlane();
704  size_t tpc = input[id]->GetCandidates()[c].GetTPC();
705  size_t cryo = input[id]->GetCandidates()[c].GetCryo();
706 
707  float t1 = input[id]->GetCandidates()[c].GetPosition().Y(); // y --> drift in 2D space.
708 
709  // loop over 2D showers
710  for (size_t j = 0; j < input.size(); ++j) {
711  if (!input[j]->GetCandidates().size()) continue;
712  if (j == id) continue;
713 
714  // loop over candidates
715  for (size_t cj = 0; cj < input[j]->GetCandidates().size(); ++cj) {
716  size_t secondview = input[j]->GetCandidates()[cj].GetPlane();
717  size_t tpc_j = input[j]->GetCandidates()[cj].GetTPC();
718  size_t cryo_j = input[j]->GetCandidates()[cj].GetCryo();
719 
720  size_t thirdview = startview;
721 
722  const geo::CryostatGeo& cryostat = geom->Cryostat(geo::CryostatID(cryo));
723  for (size_t p = 0; p < cryostat.MaxPlanes(); p++)
724  if ((p == startview) || (p == secondview)) { continue; }
725  else {
726  thirdview = p;
727  break;
728  }
729 
730  if ((startview != secondview) && (tpc == tpc_j) && (cryo == cryo_j)) {
731  float t2 = input[j]->GetCandidates()[cj].GetPosition().Y();
732  float dist = fabs(t2 - t1);
733 
734  if ((dist < mindist) && Validate(detProp, input, id, j, c, cj, thirdview)) {
735  mindist = dist;
736  pairs.clear();
737  pairs.push_back(input[id]);
738  pairs.push_back(input[j]);
739  idsave = j;
740  index = j;
741  idcsave = c;
742  idcjsave = cj;
743  found = true;
744  }
745  }
746  }
747  }
748 
749  c++;
750  }
751 
752  if (found && pairs.size()) {
753  input[id]->SetIdCandidate(idcsave);
754  input[idsave]->SetIdCandidate(idcjsave);
755  Make3DSeg(detProp, pairs);
756  }
757 
758  return index;
759 }
760 
762  std::vector<ems::DirOfGamma*> pair)
763 {
764  if (pair.size() < 2) return;
765 
766  // to build a track correctly 2d hits must belong to the same tpc
767  size_t tpc1 = pair[0]->GetFirstHit()->WireID().TPC;
768  size_t tpc2 = pair[1]->GetFirstHit()->WireID().TPC;
769 
770  std::vector<art::Ptr<recob::Hit>> vec1 = pair[0]->GetIniHits();
771  std::vector<art::Ptr<recob::Hit>> vec2 = pair[1]->GetIniHits();
772 
773  if ((vec1.size() < 3) && (vec2.size() < 3)) return;
774 
775  std::vector<art::Ptr<recob::Hit>> hitscl1uniquetpc;
776  std::vector<art::Ptr<recob::Hit>> hitscl2uniquetpc;
777 
778  if (tpc1 == tpc2)
779  for (size_t i = 0; i < vec1.size(); ++i)
780  for (size_t j = 0; j < vec2.size(); ++j)
781  if ((vec1[i]->WireID().TPC == vec2[j]->WireID().TPC) && (tpc1 == vec2[j]->WireID().TPC)) {
782  hitscl1uniquetpc.push_back(vec1[i]);
783  hitscl2uniquetpc.push_back(vec2[j]);
784  }
785 
786  if ((hitscl1uniquetpc.size() > 2) && (hitscl2uniquetpc.size() > 2)) {
787  pma::Track3D* trk =
788  fProjectionMatchingAlg.buildSegment(detProp, hitscl1uniquetpc, hitscl2uniquetpc);
789 
790  // turn the track that front is at vertex - easier to handle associations.
791  if ((trk->back()->Hit2DPtr() == pair[0]->GetFirstHit()) ||
792  (trk->back()->Hit2DPtr() == pair[1]->GetFirstHit()))
793  trk->Flip();
794 
795  IniSeg initrack;
796  initrack.idcl1 = pair[0]->GetIdCl();
797  initrack.view1 = pair[0]->GetFirstHit()->WireID().Plane;
798  initrack.hits1 = hitscl1uniquetpc;
799  initrack.idcl2 = pair[1]->GetIdCl();
800  initrack.view2 = pair[1]->GetFirstHit()->WireID().Plane;
801  initrack.hits2 = hitscl2uniquetpc;
802  initrack.track = trk;
803 
804  fInisegs.push_back(initrack);
805  }
806 }
807 
809  std::vector<ems::DirOfGamma*> input,
810  size_t id1,
811  size_t id2,
812  size_t c1,
813  size_t c2,
814  size_t plane3)
815 {
816  bool result = false;
817  if (id1 == id2) return false;
818 
819  std::vector<art::Ptr<recob::Hit>> vec1 = input[id1]->GetCandidates()[c1].GetIniHits();
820  std::vector<art::Ptr<recob::Hit>> vec2 = input[id2]->GetCandidates()[c2].GetIniHits();
821 
822  if ((vec1.size() < 3) || (vec2.size() < 3)) return false;
823 
824  std::vector<art::Ptr<recob::Hit>> hitscl1uniquetpc;
825  std::vector<art::Ptr<recob::Hit>> hitscl2uniquetpc;
826 
827  size_t tpc = vec1[0]->WireID().TPC;
828  for (size_t i = 0; i < vec1.size(); ++i)
829  for (size_t j = 0; j < vec2.size(); ++j)
830  if ((vec1[i]->WireID().TPC == tpc) && (vec2[j]->WireID().TPC == tpc)) {
831  hitscl1uniquetpc.push_back(vec1[i]);
832  hitscl2uniquetpc.push_back(vec2[j]);
833  }
834 
835  if ((hitscl1uniquetpc.size() < 3) || (hitscl2uniquetpc.size() < 3)) return false;
836 
838  fProjectionMatchingAlg.buildSegment(detProp, hitscl1uniquetpc, hitscl2uniquetpc);
839  for (size_t i = 0; i < input.size(); ++i) {
840  std::vector<Hit2D*> hits2dcl = input[i]->GetHits2D();
841  for (size_t h = 0; h < hits2dcl.size(); ++h) {
842  TVector2 pfront = pma::GetProjectionToPlane(
843  track->front()->Point3D(), plane3, track->FrontTPC(), track->FrontCryo());
844  TVector2 pback = pma::GetProjectionToPlane(
845  track->back()->Point3D(), plane3, track->BackTPC(), track->BackCryo());
846  if ((pma::Dist2(hits2dcl[h]->GetPointCm(), pfront) < 1.0F) &&
847  (pma::Dist2(hits2dcl[h]->GetPointCm(), pback) < 1.0F)) {
848  result = true;
849  break;
850  }
851  }
852  }
853  delete track;
854  return result;
855 }
856 
857 bool ems::EMShower3D::Has(const std::vector<size_t>& v, size_t idx)
858 {
859  for (auto c : v)
860  if (c == idx) return true;
861  return false;
862 }
863 
865  double r2d,
866  const std::vector<art::Ptr<recob::Hit>>& hits_in,
867  std::vector<size_t>& used,
869 {
870 
871  hits_out.clear();
872 
873  const double gapMargin = 5.0; // can be changed to f(id_tpc1, id_tpc2)
874  size_t idx = 0;
875 
876  while ((idx < hits_in.size()) && Has(used, idx))
877  idx++;
878 
879  if (idx < hits_in.size()) {
880  hits_out.push_back(hits_in[idx]);
881  used.push_back(idx);
882 
883  double r2d2 = r2d * r2d;
884  double gapMargin2 = sqrt(2 * gapMargin * gapMargin);
885  gapMargin2 = (gapMargin2 + r2d) * (gapMargin2 + r2d);
886 
887  bool collect = true;
888  while (collect) {
889  collect = false;
890  for (size_t i = 0; i < hits_in.size(); i++)
891  if (!Has(used, i)) {
892  art::Ptr<recob::Hit> hi = hits_in[i];
893  TVector2 hi_cm = pma::WireDriftToCm(detProp,
894  hi->WireID().Wire,
895  hi->PeakTime(),
896  hi->WireID().Plane,
897  hi->WireID().TPC,
898  hi->WireID().Cryostat);
899 
900  bool accept = false;
901  // for (auto const& ho : hits_out)
902  for (size_t idx_o = 0; idx_o < hits_out.size(); idx_o++) {
903  art::Ptr<recob::Hit> ho = hits_out[idx_o];
904 
905  double d2 = pma::Dist2(hi_cm,
906  pma::WireDriftToCm(detProp,
907  ho->WireID().Wire,
908  ho->PeakTime(),
909  ho->WireID().Plane,
910  ho->WireID().TPC,
911  ho->WireID().Cryostat));
912 
913  if (hi->WireID().TPC == ho->WireID().TPC) {
914  if (d2 < r2d2) {
915  accept = true;
916  break;
917  }
918  }
919  else {
920  if (d2 < gapMargin2) {
921  accept = true;
922  break;
923  }
924  }
925  }
926  if (accept) {
927  collect = true;
928  hits_out.push_back(hi);
929  used.push_back(i);
930  }
931  }
932  }
933  return true;
934  }
935  else
936  return false;
937 }
938 
940  double r2d,
941  const std::vector<art::Ptr<recob::Hit>>& hits_in,
943 {
944  size_t min_size = hits_in.size() / 5;
945  if (min_size < 3) min_size = 3;
946 
947  std::vector<size_t> used;
948  std::vector<art::Ptr<recob::Hit>> close_hits;
949 
950  while (GetCloseHits(detProp, r2d, hits_in, used, close_hits)) {
951  if (close_hits.size() > min_size)
952  for (auto h : close_hits)
953  hits_out.push_back(h);
954  }
955 }
void Link(detinfo::DetectorPropertiesData const &detProp, std::vector< ems::DirOfGamma * > input)
Utilities related to art service access.
TTree * t1
Definition: plottest35.C:26
unsigned int MaxPlanes() const
Returns the largest number of planes among the TPCs in this cryostat.
Implementation of the Projection Matching Algorithm.
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:78
Declaration of signal hit object.
pma::Hit3D const * front() const
Definition: PmaTrack3D.h:87
std::vector< std::vector< art::Ptr< recob::Hit > > > fClusters
unsigned int Cryo() const noexcept
Definition: PmaHit3D.h:58
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:67
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
art::Handle< std::vector< recob::Cluster > > fCluListHandle
TVector3 const & Point3D() const
Definition: PmaHit3D.h:50
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:286
Set of hits with a 2D structure.
Definition: Cluster.h:69
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
CryostatGeo const & Cryostat(CryostatID const &cryoid=cryostat_zero) const
Returns the specified cryostat.
unsigned int BackTPC() const
Definition: PmaTrack3D.h:121
void FilterOutSmallParts(detinfo::DetectorPropertiesData const &detProp, double r2d, const std::vector< art::Ptr< recob::Hit >> &hits_in, std::vector< art::Ptr< recob::Hit >> &hits_out)
unsigned int fIniIndex
unsigned int fClIndex
std::vector< IniSeg > fInisegs
bool Validate(detinfo::DetectorPropertiesData const &detProp, std::vector< ems::DirOfGamma * > input, size_t id1, size_t id2, size_t c1, size_t c2, size_t plane3)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
unsigned int BackCryo() const
Definition: PmaTrack3D.h:122
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
void produce(art::Event &e) override
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::vector< Hit2D * > const & GetHits2D() const
Definition: DirOfGamma.h:133
double selectInitialHits(pma::Track3D &trk, unsigned int view=geo::kZ, unsigned int *nused=0) const
void hits()
Definition: readHits.C:15
TCanvas * c1
Definition: plotHisto.C:7
unsigned int fTrkIndex
EMShower3D(fhicl::ParameterSet const &p)
TCanvas * c2
Definition: plot_hist.C:75
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void Make3DSeg(detinfo::DetectorPropertiesData const &detProp, std::vector< ems::DirOfGamma * > pair)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
std::string fTrk3DModuleLabel
A trajectory in space reconstructed from hits.
double ConvertXToTicks(double X, int p, int t, int c) const
Provides recob::Track data product.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
unsigned int FrontTPC() const
Definition: PmaTrack3D.h:118
size_t LinkCandidates(detinfo::DetectorPropertiesData const &detProp, std::vector< ems::DirOfGamma * > input, size_t id)
unsigned int FrontCryo() const
Definition: PmaTrack3D.h:119
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Declaration of cluster object.
calo::CalorimetryAlg fCalorimetryAlg
std::vector< ems::DirOfGamma * > CollectShower2D(detinfo::DetectorPropertiesData const &detProp, art::Event const &e)
Definition: DirOfGamma.h:22
size_type size() const
Definition: PtrVector.h:302
TTree * t2
Definition: plottest35.C:36
TFile fb("Li6.root")
double ConvertTicksToX(double ticks, int p, int t, int c) const
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.
bool GetCloseHits(detinfo::DetectorPropertiesData const &detProp, double r2d, const std::vector< art::Ptr< recob::Hit >> &hits_in, std::vector< size_t > &used, std::vector< art::Ptr< recob::Hit >> &hits_out)
std::vector< IniSeg > fPMA3D
Implementation of the Projection Matching Algorithm.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
recob::Track ConvertFrom(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, pma::Track3D &src)
std::vector< size_t > fTracksNotUsed
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Contains all timing reference information for the detector.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
pma::ProjectionMatchingAlg fProjectionMatchingAlg
recob::Track ConvertFrom2(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, pma::Track3D &src)
std::vector< size_t > fClustersNotUsed
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:263
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
pma::Hit3D const * back() const
Definition: PmaTrack3D.h:88
art::Ptr< recob::Hit > const & Hit2DPtr() const
Definition: PmaHit3D.h:48
bool Has(const std::vector< size_t > &v, size_t idx)
std::vector< IniSeg > fSeltracks
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:992
bool Flip(const detinfo::DetectorPropertiesData &detProp, std::vector< pma::Track3D * > &allTracks)
Definition: PmaTrack3D.cxx:527
void Reoptimize(detinfo::DetectorPropertiesData const &detProp)
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
unsigned int TPC() const noexcept
Definition: PmaHit3D.h:59
Float_t e
Definition: plot.C:35
size_t size() const
Definition: PmaTrack3D.h:89
std::string fCluModuleLabel
Float_t track
Definition: plot.C:35
void clear()
Definition: PtrVector.h:533
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192