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