100 std::vector<double>
const& Wire_2dvtx,
101 std::vector<double>
const& Time_2dvtx,
102 std::vector<double>
const& Plane_2dvtx);
109 std::vector<double>
const& merge_vtxY,
110 std::vector<double>
const& merge_vtxZ,
111 std::vector<double>
const& merge_vtxStgth);
137 std::vector<std::vector<int>>
Cls;
171 produces<std::vector<recob::Vertex>>();
172 produces<std::vector<recob::EndPoint2D>>();
173 produces<art::Assns<recob::EndPoint2D, recob::Hit>>();
177 produces<art::Assns<recob::Vertex, recob::Hit>>();
178 produces<art::Assns<recob::Vertex, recob::Shower>>();
179 produces<art::Assns<recob::Vertex, recob::Track>>();
182 Cls.resize(geom->Nplanes(), std::vector<int>());
198 if (tpc.Nplanes() > 2) {
200 std::cout <<
"yeah" << std::endl;
207 auto vcol = std::make_unique<std::vector<recob::Vertex>>();
208 auto epcol = std::make_unique<std::vector<recob::EndPoint2D>>();
209 auto assnep = std::make_unique<art::Assns<recob::EndPoint2D, recob::Hit>>();
210 auto assnsh = std::make_unique<art::Assns<recob::Vertex, recob::Shower>>();
211 auto assntr = std::make_unique<art::Assns<recob::Vertex, recob::Track>>();
212 auto assnh = std::make_unique<art::Assns<recob::Vertex, recob::Hit>>();
221 std::vector<art::Ptr<recob::EndPoint2D>> ccrawlerEndPoints;
228 mf::LogWarning(
"FeatureVertexFinder") <<
"Failed to get EndPoint2d's from Cluster Crawler";
237 std::vector<art::Ptr<recob::EndPoint2D>> cornerEndPoints;
244 mf::LogWarning(
"FeatureVertexFinder") <<
"Failed to get EndPoint2d's from Corner Finder";
259 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
271 mf::LogWarning(
"FeatureVertexFinder") <<
"Failed to get Cluster from default cluster module";
293 double tempxyz[3] = {
298 if (tempxyz[0] == 0 && tempxyz[1] == 0 && tempxyz[2] == 0) {
continue; }
300 vcol->push_back(the3Dvertex);
313 if (temp2dXYZ.X() == 0 && temp2dXYZ.Y() == 0 && temp2dXYZ.Z() == 0) {
continue; }
318 double EndPoint2d_TimeTick = detProp.ConvertXToTicks(temp2dXYZ.X(), planeID);
319 int EndPoint2d_Wire = 0;
320 int EndPoint2d_Channel = 0;
368 double tempxyz[3] = {
373 if (bail > 0) {
continue; }
374 if (tempxyz[0] == 0 && tempxyz[1] == 0 && tempxyz[2] == 0) {
continue; }
378 vcol->push_back(the3Dvertex);
395 double EndPoint2d_TimeTick = detProp.ConvertXToTicks(temp2dXYZ.X(), planeID);
396 int EndPoint2d_Wire = 0;
397 int EndPoint2d_Channel = 0;
431 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
433 for (
size_t i = 0; i < epcol->size(); ++i)
435 for (
size_t i = 0; i < vcol->size(); ++i)
438 evt.
put(std::move(epcol));
439 evt.
put(std::move(vcol));
440 evt.
put(std::move(assnep));
441 evt.
put(std::move(assntr));
442 evt.
put(std::move(assnsh));
443 evt.
put(std::move(assnh));
484 double y = 0.,
z = 0.;
485 double yy = 0.,
zz = 0.;
486 double yy2 = 0., zz2 = 0.;
487 double yy3 = 0., zz3 = 0.;
489 auto const numEndPoints =
size(EndPoints);
492 for (
size_t iendpt1 = 0; iendpt1 < numEndPoints; ++iendpt1) {
493 for (
size_t iendpt2 = iendpt1 + 1; iendpt2 < numEndPoints; ++iendpt2) {
495 auto const& endpt1 = *EndPoints[iendpt1];
496 auto const& endpt2 = *EndPoints[iendpt2];
502 if (endpt1_planeid.Plane == endpt2_planeid.Plane)
continue;
505 float tempXFeature1 = detProp.
ConvertTicksToX(endpt1.DriftTime(), endpt1_planeid);
506 float tempXFeature2 = detProp.
ConvertTicksToX(endpt2.DriftTime(), endpt2_planeid);
509 if (
std::abs(tempXFeature1 - tempXFeature2) >= 0.5)
continue;
535 for (
size_t iendpt3 = iendpt2 + 1; iendpt3 < numEndPoints; ++iendpt3) {
536 auto const& endpt3 = *EndPoints[iendpt3];
543 if (endpt3_planeid.Plane == endpt2_planeid.Plane ||
544 endpt3_planeid.Plane == endpt1_planeid.Plane)
549 float const tempXFeature3 = detProp.
ConvertTicksToX(endpt3.DriftTime(), endpt3_wireid);
551 if (
std::abs(tempXFeature3 - tempXFeature2) >= 1.0 ||
552 std::abs(tempXFeature3 - tempXFeature1) >= 1.0) {
568 detProp.
ConvertTicksToX(endpt1.DriftTime(), endpt1_wireid.asPlaneID()));
575 candidate_strength.push_back(endpt1.Strength() + endpt2.Strength() + endpt3.Strength());
580 if (iendpt1 < numEndPoints) { ++iendpt1; }
581 if (iendpt2 < numEndPoints) { ++iendpt2; }
582 if (iendpt3 < numEndPoints) { ++iendpt3; }
599 int nClustersFound = 0;
605 for (
size_t iclu = 0; iclu < RawClusters.
size(); ++iclu) {
607 std::vector<art::Ptr<recob::Hit>>
hit = fmhit.at(iclu);
611 if (hit.size() < 15) {
continue; }
616 if (view >= Cls.size()) {
617 std::cerr << __PRETTY_FUNCTION__ <<
"\033[93m Logic error in my code ... view " << view
618 <<
" not supported ! \033[00m" << std::endl;
622 Cls.at(RawClusters.
at(iclu)->View()).push_back(iclu);
626 std::vector<double> wires;
627 std::vector<double> times;
634 for (
size_t i = 0; i < hit.size(); ++i) {
635 wires.push_back(hit[i]->
WireID().Wire);
636 times.push_back(hit[i]->PeakTime());
647 TGraph* the2Dtrack =
new TGraph(n, &wires[0], ×[0]);
650 the2Dtrack->Fit(
"pol1",
"Q");
651 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
654 pol1->GetParameters(par);
655 parerror[1] = pol1->GetParError(1);
656 double FitChi2 = pol1->GetChisquare();
657 double FitNDF = pol1->GetNDF();
659 double fitGoodness = FitChi2 / FitNDF;
662 if (fitGoodness > 10) {
663 dtdwstart.push_back(std::tan(RawClusters[iclu]->StartAngle()));
676 <<
"Fitter failed, using the clusters default dTdW()";
678 dtdwstart.push_back(std::tan(RawClusters[iclu]->StartAngle()));
689 dtdwstart.push_back(std::tan(RawClusters[iclu]->StartAngle()));
699 for (
unsigned int i = 0; i < tpc.Nplanes(); ++i) {
703 if (Cls[i].
size() >= 1) {
707 for (
unsigned int j = 0; j < Cls[i].size(); ++j) {
709 Clu_Plane.push_back(RawClusters.
at(Cls.at(i).at(j))->View());
718 Clu_Length.push_back(std::sqrt(pow((RawClusters.
at(Cls.at(i).at(j))->StartWire() -
719 RawClusters.
at(Cls.at(i).at(j))->EndWire()) *
722 pow(RawClusters.
at(Cls.at(i).at(j))->StartTick() -
723 RawClusters.
at(Cls.at(i).at(j))->EndTick(),
730 RawClusters.
at(Cls.at(i).at(j))->StartTick() -
731 (
dtdwstart[Cls[i][j]] * RawClusters.
at(Cls.at(i).at(j))->StartWire()));
738 RawClusters.
at(Cls.at(i).at(j))->EndTick() -
739 (
dtdwstart[Cls[i][j]] * RawClusters.
at(Cls.at(i).at(j))->EndWire()));
762 for (
unsigned int n = nClustersFound;
n > 0;
n--) {
767 for (
unsigned int m = 0; m <
n; m++) {
776 float intersection_X =
783 float intersection_X2 =
791 if (intersection_X2 < 1) { intersection_X2 = -999; }
792 if (intersection_X2 > geom->
Nwires(planeid)) { intersection_X2 = -999; }
793 if (intersection_Y2 < 0) { intersection_Y2 = -999; }
795 if (intersection_X < 1) { intersection_X = -999; }
796 if (intersection_X > geom->
Nwires(planeid)) { intersection_X = -999; }
797 if (intersection_Y < 0) { intersection_Y = -999; }
804 std::vector<art::Ptr<recob::Hit>> hitClu1 = fmhit.at(m);
805 std::vector<art::Ptr<recob::Hit>> hitClu2 = fmhit.at(n);
812 intersection_X2 = -999;
813 intersection_Y2 = -999;
818 intersection_X = -999;
819 intersection_Y = -999;
828 intersection_X2 = -999;
829 intersection_Y2 = -999;
834 intersection_X = -999;
835 intersection_Y = -999;
839 mf::LogWarning(
"FeatureVertexFinder") <<
"FindManyHit Function faild";
840 intersection_X = -999;
841 intersection_Y = -999;
842 intersection_X2 = -999;
843 intersection_Y2 = -999;
850 if (intersection_X2 > 1 && intersection_Y2 > 0 &&
851 (intersection_X2 < geom->Nwires(planeid)) &&
862 if (intersection_X > 1 && intersection_Y > 0 &&
863 (intersection_X < geom->Nwires(planeid)) &&
881 std::vector<double>
const& Wire_2dvtx,
882 std::vector<double>
const& Time_2dvtx,
883 std::vector<double>
const& Plane_2dvtx)
886 std::vector<double> vtx_wire_merged;
887 std::vector<double> vtx_time_merged;
888 std::vector<double> vtx_plane_merged;
890 double y_coord = 0, z_coord = 0;
900 for (
size_t vtxloop1 = 0; vtxloop1 < Wire_2dvtx.size(); vtxloop1++) {
901 if (Wire_2dvtx[vtxloop1] < 0) {
continue; }
907 for (
size_t vtxloop2 = vtxloop1 + 1; vtxloop2 < Wire_2dvtx.size(); vtxloop2++) {
908 if (Wire_2dvtx[vtxloop2] < 0) {
continue; }
912 if (Plane_2dvtx[vtxloop1] == Plane_2dvtx[vtxloop2]) {
917 if (fabs(Wire_2dvtx[vtxloop1] - Wire_2dvtx[vtxloop2]) < 4) {
921 if (fabs(Time_2dvtx[vtxloop1] - Time_2dvtx[vtxloop2]) < 10) {
922 vtx_wire_merged.push_back(((Wire_2dvtx[vtxloop2] + Wire_2dvtx[vtxloop1]) / 2));
923 vtx_time_merged.push_back(((Time_2dvtx[vtxloop2] + Time_2dvtx[vtxloop1]) / 2));
924 vtx_plane_merged.push_back(Plane_2dvtx[vtxloop1]);
927 if (vtxloop2 < Wire_2dvtx.size()) { vtxloop2++; }
928 if (vtxloop1 < Wire_2dvtx.size()) { vtxloop1++; }
934 vtx_wire_merged.push_back(Wire_2dvtx[vtxloop1]);
935 vtx_time_merged.push_back(Time_2dvtx[vtxloop1]);
936 vtx_plane_merged.push_back(Plane_2dvtx[vtxloop1]);
942 uint32_t vtx1_channel = 0;
943 uint32_t vtx2_channel = 0;
950 for (
unsigned int vtx = vtx_wire_merged.size(); vtx > 0; --vtx) {
951 for (
unsigned int vtx1 = 0; vtx1 < vtx; ++vtx1) {
954 if (vtx_plane_merged[vtx1] == vtx_plane_merged[vtx])
continue;
960 geo::PlaneID const vtx1_planeid(tpcid, vtx_plane_merged[vtx1]);
961 geo::WireID const vtx1_wireid(vtx1_planeid, vtx_wire_merged[vtx1]);
966 mf::LogWarning(
"FeatureVertexFinder") <<
"PlaneWireToChannel Failed";
970 geo::PlaneID const vtx2_planeid(tpcid, vtx_plane_merged[vtx]);
971 geo::WireID const vtx2_wireid(vtx2_planeid, vtx_wire_merged[vtx]);
976 mf::LogWarning(
"FeatureVertexFinder") <<
"PlaneWireToChannel Failed";
982 mf::LogWarning(
"FeatureVertexFinder") <<
"match failed for some reason";
987 float tempXCluster1 = detProp.
ConvertTicksToX(vtx_time_merged[vtx1], vtx1_planeid);
988 float tempXCluster2 = detProp.
ConvertTicksToX(vtx_time_merged[vtx], vtx2_planeid);
1013 std::vector<double>
const& merge_vtxX,
1014 std::vector<double>
const& merge_vtxY,
1015 std::vector<double>
const& merge_vtxZ,
1016 std::vector<double>
const& merge_vtxStgth)
1019 std::vector<double> x_3dVertex_dupRemoved = {0.};
1020 std::vector<double> y_3dVertex_dupRemoved = {0.};
1021 std::vector<double> z_3dVertex_dupRemoved = {0.};
1022 std::vector<double> strength_dupRemoved = {0.};
1026 for (
size_t dup = 0; dup < merge_vtxX.size(); dup++) {
1028 float tempX_dup = merge_vtxX[dup];
1029 float tempY_dup = merge_vtxY[dup];
1030 float tempZ_dup = merge_vtxZ[dup];
1031 float tempStgth = merge_vtxStgth[dup];
1035 bool duplicate_found =
false;
1039 for (
size_t check = dup + 1; check < merge_vtxX.size(); check++) {
1044 if (
std::abs(merge_vtxX[check] - tempX_dup) < 0.1 &&
1045 std::abs(merge_vtxY[check] - tempY_dup) < 0.1 &&
1046 std::abs(merge_vtxZ[check] - tempZ_dup) < 0.1) {
1047 duplicate_found =
true;
1055 if (!duplicate_found && tempX_dup > 0) {
1056 x_3dVertex_dupRemoved.push_back(tempX_dup);
1057 y_3dVertex_dupRemoved.push_back(tempY_dup);
1058 z_3dVertex_dupRemoved.push_back(tempZ_dup);
1059 strength_dupRemoved.push_back(tempStgth);
1069 double tempX, tempY, tempZ, tempS;
1073 for (
size_t npri = 0; (npri < x_3dVertex_dupRemoved.size()) && flag; npri++) {
1075 for (
size_t mpri = 0; mpri < x_3dVertex_dupRemoved.size() - 1; mpri++) {
1077 if (strength_dupRemoved[mpri + 1] > strength_dupRemoved[mpri] ||
1078 (strength_dupRemoved[mpri + 1] == strength_dupRemoved[mpri] &&
1079 z_3dVertex_dupRemoved[mpri] > z_3dVertex_dupRemoved[mpri + 1])) {
1080 tempX = x_3dVertex_dupRemoved[mpri];
1081 x_3dVertex_dupRemoved[mpri] = x_3dVertex_dupRemoved[mpri + 1];
1082 x_3dVertex_dupRemoved[mpri + 1] = tempX;
1084 tempY = y_3dVertex_dupRemoved[mpri];
1085 y_3dVertex_dupRemoved[mpri] = y_3dVertex_dupRemoved[mpri + 1];
1086 y_3dVertex_dupRemoved[mpri + 1] = tempY;
1088 tempZ = z_3dVertex_dupRemoved[mpri];
1089 z_3dVertex_dupRemoved[mpri] = z_3dVertex_dupRemoved[mpri + 1];
1090 z_3dVertex_dupRemoved[mpri + 1] = tempZ;
1092 tempS = strength_dupRemoved[mpri];
1093 strength_dupRemoved[mpri] = strength_dupRemoved[mpri + 1];
1094 strength_dupRemoved[mpri + 1] = tempS;
1104 for (
size_t count = 0; count < x_3dVertex_dupRemoved.size(); count++) {
std::vector< double > Clu_Plane
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< double > TwoDvtx_wire
Encapsulate the construction of a single cyostat.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::vector< double > TwoDvtx_time
Declaration of signal hit object.
std::vector< double > Clu_EndPos_Wire
std::vector< double > candidate_y
EDProducer(fhicl::ParameterSet const &pset)
The data type to uniquely identify a Plane.
Geometry information for a single TPC.
void produce(art::Event &evt) override
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< double > candidate_strength
std::vector< double > Clu_Length
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
std::vector< double > Clu_Yintercept
WireID_t Wire
Index of the wire within its plane.
void MergeAndSort3dVtxCandidate(std::vector< double > const &merge_vtxX, std::vector< double > const &merge_vtxY, std::vector< double > const &merge_vtxZ, std::vector< double > const &merge_vtxStgth)
std::vector< double > Clu_EndPos_TimeTick
std::vector< double > MergeSort3dVtx_strength
std::vector< double > candidate_z
Definition of vertex object for LArSoft.
std::vector< double > Clu_StartPos_TimeTick
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
std::vector< std::vector< int > > Cls
void Find3dVtxFrom2dClusterVtxCand(detinfo::DetectorPropertiesData const &detProp, std::vector< double > const &Wire_2dvtx, std::vector< double > const &Time_2dvtx, std::vector< double > const &Plane_2dvtx)
#define DEFINE_ART_MODULE(klass)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
std::vector< double > Clu_StartPos_Wire
std::vector< double > MergeSort3dVtx_zpos
std::string fClusterModuleLabel
void Find2dClusterVertexCandidates(detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Cluster > RawClusters, art::FindManyP< recob::Hit > fmhit)
unsigned int NumberTimeSamples() const
Provides recob::Track data product.
std::string fCornerFinderModuleLabel
reference at(size_type n)
std::vector< double > dtdwstart
FeatureVertexFinder(fhicl::ParameterSet const &pset)
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
Declaration of cluster object.
Detector simulation of raw signals on wires.
std::vector< double > MergeSort3dVtx_ypos
double ConvertTicksToX(double ticks, int p, int t, int c) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
std::vector< double > MergeSort3dVtx_xpos
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
std::vector< double > Clu_Slope
std::vector< double > TwoDvtx_plane
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
raw::ChannelID_t NearestChannel(Point_t const &worldLoc, PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
std::vector< double > dtdwstartError
constexpr WireID()=default
Default constructor: an invalid TPC ID.
std::vector< double > candidate_x
std::string fCCrawlerEndPoint2dModuleLabel
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
std::string fHitModuleLabel
void Get3dVertexCandidates(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::EndPoint2D >> const &EndPoints, bool PlaneDet)
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane.
std::vector< double > Clu_Yintercept2