43 std::vector<art::Ptr<recob::Hit>> hits1;
44 std::vector<art::Ptr<recob::Hit>> hits2;
45 std::vector<art::Ptr<recob::Hit>> hits3;
82 std::vector<ems::DirOfGamma*> pair);
85 std::vector<ems::DirOfGamma*> input,
100 std::vector<size_t>& used,
103 bool Has(
const std::vector<size_t>& v,
size_t idx);
106 std::vector<ems::DirOfGamma*> input,
113 std::vector<std::vector<art::Ptr<recob::Hit>>>
fClusters;
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>>();
176 src[0]->WireID().planeID());
183 auto const* geom = lar::providerFrom<geo::Geometry>();
186 size_t nusedhitsmax = 0;
188 for (
unsigned int p = 0; p < nplanes; ++p) {
189 unsigned int nusedP = 0;
192 if (nusedP > nusedhitsmax) {
193 nusedhitsmax = nusedP;
198 std::vector<std::vector<double>> vdedx;
199 std::vector<double> dedx;
201 for (
unsigned int p = 0; p < nplanes; ++p) {
202 unsigned int nusedP = 0;
206 dedx.push_back(dEdxplane);
207 if (
int(p) == bestplane)
211 vdedx.push_back(dedx);
214 std::vector<TVector3> xyz, dircos;
216 for (
size_t i = 0; i < src.
size(); i++) {
217 xyz.push_back(src[i]->
Point3D());
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);
226 dircos.push_back(dircos.back());
245 auto const* geom = lar::providerFrom<geo::Geometry>();
249 size_t nusedhitsmax = 0;
251 for (
unsigned int p = 0; p < nplanes; ++p) {
252 unsigned int nusedP = 0;
255 if (nusedP > nusedhitsmax) {
256 nusedhitsmax = nusedP;
261 std::vector<std::vector<double>> vdedx;
262 std::vector<double> dedx;
264 for (
unsigned int p = 0; p < nplanes; ++p) {
265 unsigned int nusedP = 0;
269 dedx.push_back(dEdxplane);
270 if (
int(p) == bestplane)
274 vdedx.push_back(dedx);
277 std::vector<TVector3> xyz, dircos;
279 for (
size_t i = 0; i < src.
size(); i++) {
280 xyz.push_back(src[i]->
Point3D());
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);
289 dircos.push_back(dircos.back());
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>>();
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>>();
335 std::vector<art::Ptr<recob::Hit>> hitlist;
341 std::vector<ems::DirOfGamma*> showernviews =
CollectShower2D(detProp, e);
343 Link(detProp, showernviews);
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++)
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);
371 std::vector<art::Ptr<recob::Cluster>> cl2d;
375 std::vector<art::Ptr<recob::Hit>> hits2d;
378 spStart = allsp->
size();
379 for (
int h =
trk.track->size() - 1; h >= 0; h--) {
383 if ((h == 0) || (sp_pos[0] != h3d->
Point3D().X()) || (sp_pos[1] != h3d->
Point3D().Y()) ||
384 (sp_pos[2] != h3d->
Point3D().Z())) {
390 sp_pos[0] = h3d->
Point3D().X();
391 sp_pos[1] = h3d->
Point3D().Y();
392 sp_pos[2] = h3d->
Point3D().Z();
401 spEnd = allsp->size();
403 if (vertices->size()) {
404 size_t vtx_idx = (size_t)(vertices->size() - 1);
422 std::vector<art::Ptr<recob::Cluster>> cl2d;
426 std::vector<art::Ptr<recob::Hit>> hits2d;
429 spStart = allsp->
size();
430 for (
int h =
trk.track->size() - 1; h >= 0; h--) {
434 if ((h == 0) || (sp_pos[0] != h3d->
Point3D().X()) || (sp_pos[1] != h3d->
Point3D().Y()) ||
435 (sp_pos[2] != h3d->
Point3D().Z())) {
441 sp_pos[0] = h3d->
Point3D().X();
442 sp_pos[1] = h3d->
Point3D().Y();
443 sp_pos[2] = h3d->
Point3D().Z();
452 spEnd = allsp->size();
472 for (
unsigned int i = 0; i < showernviews.size(); i++)
473 delete showernviews[i];
475 for (
unsigned int i = 0; i < fSeltracks.size(); i++)
476 delete fSeltracks[i].
track;
478 for (
unsigned int i = 0; i <
fInisegs.size(); i++)
481 for (
unsigned int i = 0; i < fPMA3D.size(); i++)
482 delete fPMA3D[i].
track;
485 e.
put(std::move(tracks));
486 e.
put(std::move(vertices));
487 e.
put(std::move(clusters));
488 e.
put(std::move(allsp));
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));
501 const float min_dist = 3.0F;
512 TVector3 p1 =
fSeltracks[ta].track->front()->Point3D();
513 TVector3 p2 =
fSeltracks[tb].track->front()->Point3D();
524 std::vector<art::Ptr<recob::Hit>> hits3 =
fSeltracks[ta].hits1;
526 for (
size_t h = 0; h <
fSeltracks[ta].hits2.size(); ++h)
532 for (
size_t h = 0; h <
fSeltracks[tb].hits1.size(); ++h)
538 for (
size_t h = 0; h <
fSeltracks[tb].hits2.size(); ++h)
555 initrack.idcl3 = idcl3;
557 initrack.view3 = view3;
559 initrack.hits3 = hits3;
563 initrack.track =
track;
585 std::vector<ems::DirOfGamma*> input;
590 std::vector<art::Ptr<recob::Hit>> hitlist;
593 if (hitlist.size() > 5) {
594 std::vector<art::Ptr<recob::Hit>> hits_out;
597 if (hits_out.size() > 5) {
602 if (sh->
GetHits2D().size()) input.push_back(sh);
612 std::vector<ems::DirOfGamma*> input)
614 std::vector<std::vector<size_t>> saveids;
615 std::vector<size_t> saveidsnotusedcls;
618 while (i < input.size()) {
619 if (!input[i]->GetCandidates().size()) {
624 double mindist = 1.0;
625 std::vector<ems::DirOfGamma*> pairs;
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;
631 float t1 = detProp.
ConvertTicksToX(input[i]->GetFirstHit()->PeakTime(), startview, tpc, cryo);
633 unsigned int idsave = 0;
634 for (
unsigned int j = 0; j < input.size(); j++) {
635 if (!input[j]->GetCandidates().size())
continue;
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;
641 if ((i != j) && (secondview != startview) && (tpc == tpc_j) && (cryo == cryo_j)) {
643 detProp.
ConvertTicksToX(input[j]->GetFirstHit()->PeakTime(), secondview, tpc_j, cryo_j);
644 float dist = fabs(t2 - t1);
646 if (dist < mindist) {
649 pairs.push_back(input[i]);
650 pairs.push_back(input[j]);
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;
664 std::vector<size_t> ids;
666 ids.push_back(idsave);
667 saveids.push_back(ids);
670 saveidsnotusedcls.push_back(i);
677 while (i < saveidsnotusedcls.size()) {
684 std::vector<ems::DirOfGamma*> input,
692 if (input[
id]->GetCandidates().
size() < 2) {
return index; }
694 double mindist = 3.0;
695 std::vector<ems::DirOfGamma*> pairs;
701 while (c < input[
id]->GetCandidates().
size()) {
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();
707 float t1 = input[id]->GetCandidates()[c].GetPosition().Y();
710 for (
size_t j = 0; j < input.size(); ++j) {
711 if (!input[j]->GetCandidates().size())
continue;
712 if (j ==
id)
continue;
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();
720 size_t thirdview = startview;
723 for (
size_t p = 0; p < cryostat.
MaxPlanes(); p++)
724 if ((p == startview) || (p == secondview)) {
continue; }
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);
734 if ((dist < mindist) &&
Validate(detProp, input,
id, j, c, cj, thirdview)) {
737 pairs.push_back(input[
id]);
738 pairs.push_back(input[j]);
752 if (found && pairs.size()) {
753 input[id]->SetIdCandidate(idcsave);
754 input[idsave]->SetIdCandidate(idcjsave);
762 std::vector<ems::DirOfGamma*> pair)
764 if (pair.size() < 2)
return;
767 size_t tpc1 = pair[0]->GetFirstHit()->WireID().TPC;
768 size_t tpc2 = pair[1]->GetFirstHit()->WireID().TPC;
770 std::vector<art::Ptr<recob::Hit>> vec1 = pair[0]->GetIniHits();
771 std::vector<art::Ptr<recob::Hit>> vec2 = pair[1]->GetIniHits();
773 if ((vec1.size() < 3) && (vec2.size() < 3))
return;
775 std::vector<art::Ptr<recob::Hit>> hitscl1uniquetpc;
776 std::vector<art::Ptr<recob::Hit>> hitscl2uniquetpc;
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]);
786 if ((hitscl1uniquetpc.size() > 2) && (hitscl2uniquetpc.size() > 2)) {
791 if ((trk->
back()->
Hit2DPtr() == pair[0]->GetFirstHit()) ||
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;
809 std::vector<ems::DirOfGamma*> input,
817 if (id1 == id2)
return false;
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();
822 if ((vec1.size() < 3) || (vec2.size() < 3))
return false;
824 std::vector<art::Ptr<recob::Hit>> hitscl1uniquetpc;
825 std::vector<art::Ptr<recob::Hit>> hitscl2uniquetpc;
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]);
835 if ((hitscl1uniquetpc.size() < 3) || (hitscl2uniquetpc.size() < 3))
return false;
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) {
846 if ((
pma::Dist2(hits2dcl[h]->GetPointCm(), pfront) < 1.0F) &&
847 (
pma::Dist2(hits2dcl[h]->GetPointCm(), pback) < 1.0F)) {
860 if (c == idx)
return true;
867 std::vector<size_t>& used,
873 const double gapMargin = 5.0;
876 while ((idx < hits_in.size()) &&
Has(used, idx))
879 if (idx < hits_in.size()) {
880 hits_out.push_back(hits_in[idx]);
883 double r2d2 = r2d * r2d;
884 double gapMargin2 = sqrt(2 * gapMargin * gapMargin);
885 gapMargin2 = (gapMargin2 + r2d) * (gapMargin2 + r2d);
890 for (
size_t i = 0; i < hits_in.size(); i++)
902 for (
size_t idx_o = 0; idx_o < hits_out.size(); idx_o++) {
920 if (d2 < gapMargin2) {
928 hits_out.push_back(hi);
944 size_t min_size = hits_in.size() / 5;
945 if (min_size < 3) min_size = 3;
947 std::vector<size_t> used;
948 std::vector<art::Ptr<recob::Hit>> close_hits;
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);
void Link(detinfo::DetectorPropertiesData const &detProp, std::vector< ems::DirOfGamma * > input)
Utilities related to art service access.
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)
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Declaration of signal hit object.
pma::Hit3D const * front() const
std::vector< std::vector< art::Ptr< recob::Hit > > > fClusters
unsigned int Cryo() const noexcept
TrackTrajectory::Flags_t Flags_t
::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
CryostatID_t Cryostat
Index of cryostat.
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)
Set of hits with a 2D structure.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
Geometry information for a single cryostat.
CryostatGeo const & Cryostat(CryostatID const &cryoid=cryostat_zero) const
Returns the specified cryostat.
unsigned int BackTPC() const
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)
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={})
unsigned int BackCryo() const
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
void produce(art::Event &e) override
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
std::vector< Hit2D * > const & GetHits2D() const
double selectInitialHits(pma::Track3D &trk, unsigned int view=geo::kZ, unsigned int *nused=0) const
EMShower3D(fhicl::ParameterSet const &p)
#define DEFINE_ART_MODULE(klass)
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)
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)
unsigned int FrontTPC() const
size_t LinkCandidates(detinfo::DetectorPropertiesData const &detProp, std::vector< ems::DirOfGamma * > input, size_t id)
unsigned int FrontCryo() const
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.
Declaration of cluster object.
calo::CalorimetryAlg fCalorimetryAlg
std::vector< ems::DirOfGamma * > CollectShower2D(detinfo::DetectorPropertiesData const &detProp, art::Event const &e)
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.
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)
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
art::Ptr< recob::Hit > const & Hit2DPtr() const
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)
void Reoptimize(detinfo::DetectorPropertiesData const &detProp)
TPCID_t TPC
Index of the TPC within its cryostat.
unsigned int TPC() const noexcept
std::string fCluModuleLabel
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
The data type to uniquely identify a cryostat.