43 : fCalDataModuleLabel{pset.
get<std::string>(
"CalDataModuleLabel")}
85 unsigned int nPlanes = wireReadoutGeom.
Nplanes(tpcid);
105 const unsigned int nTimeTicks = wireVec.at(0).NSignal();
111 auto const i_plane = planeid.Plane;
113 std::stringstream ss_tmp_name, ss_tmp_title;
114 ss_tmp_name <<
"h_WireData_" << i_plane;
116 <<
";Wire Number;Time Tick";
118 auto const num_wires = wireReadoutGeom.
Nwires(planeid);
119 if (static_cast<unsigned int>(
WireData_histos[i_plane].GetNbinsX()) == num_wires) {
126 ss_tmp_title.str().c_str(),
136 for (
auto const& wire : wireVec) {
137 std::vector<geo::WireID> possible_wireIDs = wireReadoutGeom.
ChannelToWire(wire.Channel());
140 this_wireID = possible_wireIDs.at(0);
145 mf::LogError(
"CornerFinderAlg") <<
"Bail out! No Possible Wires!\n";
148 unsigned int i_plane = this_wireID.
Plane;
149 unsigned int i_wire = this_wireID.
Wire;
153 for (
unsigned int i_time = 0; i_time < nTimeTicks; i_time++) {
154 WireData_histos.at(i_plane).SetBinContent(i_wire, i_time, wire.Signal().at(i_time));
159 for (
unsigned int i_plane = 0; i_plane < wireReadoutGeom.
Nplanes(); i_plane++) {
187 for (
unsigned int tpc = 0; tpc < cryostat.NTPC(); ++tpc) {
194 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"Doing histogram " << histos <<
", of plane " << plane
195 <<
" with start points " << startx <<
" " << starty;
199 wireReadoutGeom.
Plane({cryostat.ID().Cryostat, tpc, plane}).View(),
204 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"Total feature points now is " << corner_vector.size();
214 std::vector<recob::EndPoint2D>& corner_vector,
226 struct compare_to_value {
227 compare_to_value(
int b) { this->b = b; }
228 bool operator()(
int i,
int j)
const {
return std::abs(b - i) <
std::abs(b - j); }
239 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"Working plane " <<
pid.Plane <<
".";
244 std::vector<int> cut_points_x{0};
245 std::vector<int> cut_points_y{0};
247 for (
int ix = 1; ix <= x_bins; ix++) {
255 if (jx == ix + fTrimming_buffer)
break;
262 for (
int iy = 1; iy <= y_bins; iy++) {
270 if (jy == iy + fTrimming_buffer)
break;
278 <<
"We have a total of " << cut_points_x.size() <<
" x cut points." 279 <<
"\nWe have a total of " << cut_points_y.size() <<
" y cut points.";
281 std::vector<int> x_low{1};
282 std::vector<int> x_high{x_bins};
283 std::vector<int> y_low{1};
284 std::vector<int> y_high{y_bins};
285 bool x_change =
true;
286 bool y_change =
true;
287 while (x_change || y_change) {
292 size_t current_size = x_low.size();
294 for (
size_t il = 0; il < current_size; il++) {
296 int comp_value = (x_high.at(il) + x_low.at(il)) / 2;
297 std::sort(cut_points_x.begin(), cut_points_x.end(), compare_to_value(comp_value));
299 if (cut_points_x.at(0) <= x_low.at(il) || cut_points_x.at(0) >= x_high.at(il))
continue;
302 x_low.at(il), cut_points_x.at(0), y_low.at(il), y_high.at(il));
304 cut_points_x.at(0), x_high.at(il), y_low.at(il), y_high.at(il));
306 x_low.push_back(cut_points_x.at(0));
307 x_high.push_back(x_high.at(il));
308 y_low.push_back(y_low.at(il));
309 y_high.push_back(y_high.at(il));
311 x_high[il] = cut_points_x.at(0);
316 x_high[il] = cut_points_x.at(0);
321 x_low[il] = cut_points_x.at(0);
326 current_size = x_low.size();
328 for (
size_t il = 0; il < current_size; il++) {
330 int comp_value = (y_high.at(il) - y_low.at(il)) / 2;
331 std::sort(cut_points_y.begin(), cut_points_y.end(), compare_to_value(comp_value));
333 if (cut_points_y.at(0) <= y_low.at(il) || cut_points_y.at(0) >= y_high.at(il))
continue;
336 x_low.at(il), x_high.at(il), y_low.at(il), cut_points_y.at(0));
338 x_low.at(il), x_high.at(il), cut_points_y.at(0), y_high.at(il));
340 y_low.push_back(cut_points_y.at(0));
341 y_high.push_back(y_high.at(il));
342 x_low.push_back(x_low.at(il));
343 x_high.push_back(x_high.at(il));
345 y_high[il] = cut_points_y.at(0);
350 y_high[il] = cut_points_y.at(0);
355 y_low[il] = cut_points_y.at(0);
361 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"First point in x is " << cut_points_x.at(0);
363 std::sort(cut_points_x.begin(), cut_points_x.end(), compare_to_value(x_bins / 2));
365 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"Now the first point in x is " << cut_points_x.at(0);
367 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"First point in y is " << cut_points_y.at(0);
369 std::sort(cut_points_y.begin(), cut_points_y.end(), compare_to_value(y_bins / 2));
371 MF_LOG_DEBUG(
"CornerFinderAlg") <<
"Now the first point in y is " << cut_points_y.at(0);
374 <<
"\nIntegral on the SW side is " 375 <<
WireData_histos.at(
pid.Plane).Integral(1, cut_points_x.at(0), 1, cut_points_y.at(0))
376 <<
"\nIntegral on the SE side is " 377 <<
WireData_histos.at(
pid.Plane).Integral(cut_points_x.at(0), x_bins, 1, cut_points_y.at(0))
378 <<
"\nIntegral on the NW side is " 379 <<
WireData_histos.at(
pid.Plane).Integral(1, cut_points_x.at(0), cut_points_y.at(0), y_bins)
380 <<
"\nIntegral on the NE side is " 382 cut_points_x.at(0), x_bins, cut_points_y.at(0), y_bins);
384 for (
size_t il = 0; il < x_low.size(); il++) {
386 std::stringstream h_name;
387 h_name <<
"h_" <<
pid.Cryostat <<
"_" <<
pid.TPC <<
"_" <<
pid.Plane <<
"_sub" << il;
388 TH2F h_tmp((h_name.str()).c_str(),
390 x_high.at(il) - x_low.at(il) + 1,
393 y_high.at(il) - y_low.at(il) + 1,
397 for (
int ix = 1; ix <= (x_high.at(il) - x_low.at(il) + 1); ix++) {
398 for (
int iy = 1; iy <= (y_high.at(il) - y_low.at(il) + 1); iy++) {
399 h_tmp.SetBinContent(ix,
402 y_low.at(il) + (iy - 1)));
407 std::make_tuple(
pid.Plane, h_tmp, x_low.at(il) - 1, y_low.at(il) - 1));
416 std::vector<geo::WireID>
const& wireIDs,
418 std::vector<recob::EndPoint2D>& corner_vector,
422 const int x_bins = h_wire_data.GetNbinsX();
423 const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
424 const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
426 const int y_bins = h_wire_data.GetNbinsY();
427 const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
428 const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
433 std::stringstream conversion_name;
435 std::stringstream dx_name;
437 std::stringstream dy_name;
439 std::stringstream cornerScore_name;
441 std::stringstream maxSuppress_name;
444 TH2F conversion_histo(conversion_name.str().c_str(),
445 "Image Conversion Histogram",
453 TH2F derivativeX_histo(dx_name.str().c_str(),
454 "Partial Derivatives (x)",
462 TH2F derivativeY_histo(dy_name.str().c_str(),
463 "Partial Derivatives (y)",
471 TH2D cornerScore_histo(cornerScore_name.str().c_str(),
480 TH2D maxSuppress_histo(maxSuppress_name.str().c_str(),
481 "Corner Points (Maximum Suppressed)",
493 cornerScore_histo, wireIDs, view, maxSuppress_histo, startx, starty);
499 TH2F
const& h_wire_data,
500 std::vector<geo::WireID>
const& wireIDs,
502 std::vector<recob::EndPoint2D>& corner_vector)
504 const int x_bins = h_wire_data.GetNbinsX();
505 const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
506 const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
508 const int y_bins = h_wire_data.GetNbinsY();
509 const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
510 const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
515 std::stringstream conversion_name;
517 std::stringstream dx_name;
519 std::stringstream dy_name;
521 std::stringstream cornerScore_name;
523 std::stringstream maxSuppress_name;
526 TH2F h_conversion(conversion_name.str().c_str(),
527 "Image Conversion Histogram",
534 TH2F h_derivative_x(dx_name.str().c_str(),
535 "Partial Derivatives (x)",
542 TH2F h_derivative_y(dy_name.str().c_str(),
543 "Partial Derivatives (y)",
550 TH2D h_cornerScore(cornerScore_name.str().c_str(),
551 "Feature Point Corner Score",
558 TH2D h_maxSuppress(maxSuppress_name.str().c_str(),
559 "Corner Points (Maximum Suppressed)",
573 std::stringstream LI_name;
575 TH2F h_lineIntegralScore(
576 LI_name.str().c_str(),
"Line Integral Score", x_bins,
x_min,
x_max, y_bins, y_min, y_max);
584 double temp_integral = 0;
586 const TF2 fConversion_TF2(
"fConversion_func",
fConversion_func.c_str(), -20, 20, -20, 20);
588 for (
int ix = 1; ix <= h_conversion.GetNbinsX(); ix++) {
589 for (
int iy = 1; iy <= h_conversion.GetNbinsY(); iy++) {
591 temp_integral = h_wire_data.Integral(ix, ix, iy, iy);
598 h_conversion.SetBinContent(ix, iy, temp_integral);
606 for (
int jy = iy - fConversion_func_neighborhood;
610 h_wire_data.GetBinContent(jx, jy) * fConversion_TF2.Eval(ix - jx, iy - jy);
613 h_conversion.SetBinContent(ix, iy, temp_integral);
618 if ((temp_integral > h_wire_data.GetBinContent(ix - 1, iy) &&
619 temp_integral > h_wire_data.GetBinContent(ix + 1, iy)) ||
620 (temp_integral > h_wire_data.GetBinContent(ix, iy - 1) &&
621 temp_integral > h_wire_data.GetBinContent(ix, iy + 1)))
622 h_conversion.SetBinContent(ix, iy, temp_integral);
628 if ((temp_integral > h_wire_data.GetBinContent(ix - 1, iy) &&
629 temp_integral > h_wire_data.GetBinContent(ix + 1, iy)) ||
630 (temp_integral > h_wire_data.GetBinContent(ix, iy - 1) &&
631 temp_integral > h_wire_data.GetBinContent(ix, iy + 1)))
637 h_conversion.SetBinContent(ix, iy, temp_integral);
650 TH2F& h_derivative_x,
651 TH2F& h_derivative_y)
653 const int x_bins = h_conversion.GetNbinsX();
654 const int y_bins = h_conversion.GetNbinsY();
662 h_derivative_x.SetBinContent(ix,
664 0.5 * (h_conversion.GetBinContent(ix + 1, iy) -
665 h_conversion.GetBinContent(ix - 1, iy)) +
666 0.25 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
667 h_conversion.GetBinContent(ix - 1, iy + 1)) +
668 0.25 * (h_conversion.GetBinContent(ix + 1, iy - 1) -
669 h_conversion.GetBinContent(ix - 1, iy - 1)));
670 h_derivative_y.SetBinContent(ix,
672 0.5 * (h_conversion.GetBinContent(ix, iy + 1) -
673 h_conversion.GetBinContent(ix, iy - 1)) +
674 0.25 * (h_conversion.GetBinContent(ix - 1, iy + 1) -
675 h_conversion.GetBinContent(ix - 1, iy - 1)) +
676 0.25 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
677 h_conversion.GetBinContent(ix + 1, iy - 1)));
680 h_derivative_x.SetBinContent(
683 12 * (h_conversion.GetBinContent(ix + 1, iy) - h_conversion.GetBinContent(ix - 1, iy)) +
684 8 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
685 h_conversion.GetBinContent(ix - 1, iy + 1)) +
686 8 * (h_conversion.GetBinContent(ix + 1, iy - 1) -
687 h_conversion.GetBinContent(ix - 1, iy - 1)) +
688 2 * (h_conversion.GetBinContent(ix + 1, iy + 2) -
689 h_conversion.GetBinContent(ix - 1, iy + 2)) +
690 2 * (h_conversion.GetBinContent(ix + 1, iy - 2) -
691 h_conversion.GetBinContent(ix - 1, iy - 2)) +
693 (h_conversion.GetBinContent(ix + 2, iy) - h_conversion.GetBinContent(ix - 2, iy)) +
694 4 * (h_conversion.GetBinContent(ix + 2, iy + 1) -
695 h_conversion.GetBinContent(ix - 2, iy + 1)) +
696 4 * (h_conversion.GetBinContent(ix + 2, iy - 1) -
697 h_conversion.GetBinContent(ix - 2, iy - 1)) +
698 1 * (h_conversion.GetBinContent(ix + 2, iy + 2) -
699 h_conversion.GetBinContent(ix - 2, iy + 2)) +
700 1 * (h_conversion.GetBinContent(ix + 2, iy - 2) -
701 h_conversion.GetBinContent(ix - 2, iy - 2)));
702 h_derivative_y.SetBinContent(
705 12 * (h_conversion.GetBinContent(ix, iy + 1) - h_conversion.GetBinContent(ix, iy - 1)) +
706 8 * (h_conversion.GetBinContent(ix - 1, iy + 1) -
707 h_conversion.GetBinContent(ix - 1, iy - 1)) +
708 8 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
709 h_conversion.GetBinContent(ix + 1, iy - 1)) +
710 2 * (h_conversion.GetBinContent(ix - 2, iy + 1) -
711 h_conversion.GetBinContent(ix - 2, iy - 1)) +
712 2 * (h_conversion.GetBinContent(ix + 2, iy + 1) -
713 h_conversion.GetBinContent(ix + 2, iy - 1)) +
715 (h_conversion.GetBinContent(ix, iy + 2) - h_conversion.GetBinContent(ix, iy - 2)) +
716 4 * (h_conversion.GetBinContent(ix - 1, iy + 2) -
717 h_conversion.GetBinContent(ix - 1, iy - 2)) +
718 4 * (h_conversion.GetBinContent(ix + 1, iy + 2) -
719 h_conversion.GetBinContent(ix + 1, iy - 2)) +
720 1 * (h_conversion.GetBinContent(ix - 2, iy + 2) -
721 h_conversion.GetBinContent(ix - 2, iy - 2)) +
722 1 * (h_conversion.GetBinContent(ix + 2, iy + 2) -
723 h_conversion.GetBinContent(ix + 2, iy - 2)));
727 <<
"Sobel derivative not supported for neighborhoods > 2.";
736 h_derivative_x.SetBinContent(
739 (h_conversion.GetBinContent(ix + 1, iy) - h_conversion.GetBinContent(ix - 1, iy)));
740 h_derivative_y.SetBinContent(
743 (h_conversion.GetBinContent(ix, iy + 1) - h_conversion.GetBinContent(ix, iy - 1)));
747 <<
"Local derivative not yet supported for neighborhoods > 1.";
760 float func_blur[11][11];
761 func_blur[0][0] = 0.000000;
762 func_blur[0][1] = 0.000000;
763 func_blur[0][2] = 0.000000;
764 func_blur[0][3] = 0.000001;
765 func_blur[0][4] = 0.000002;
766 func_blur[0][5] = 0.000004;
767 func_blur[0][6] = 0.000002;
768 func_blur[0][7] = 0.000001;
769 func_blur[0][8] = 0.000000;
770 func_blur[0][9] = 0.000000;
771 func_blur[0][10] = 0.000000;
772 func_blur[1][0] = 0.000000;
773 func_blur[1][1] = 0.000000;
774 func_blur[1][2] = 0.000004;
775 func_blur[1][3] = 0.000045;
776 func_blur[1][4] = 0.000203;
777 func_blur[1][5] = 0.000335;
778 func_blur[1][6] = 0.000203;
779 func_blur[1][7] = 0.000045;
780 func_blur[1][8] = 0.000004;
781 func_blur[1][9] = 0.000000;
782 func_blur[1][10] = 0.000000;
783 func_blur[2][0] = 0.000000;
784 func_blur[2][1] = 0.000004;
785 func_blur[2][2] = 0.000123;
786 func_blur[2][3] = 0.001503;
787 func_blur[2][4] = 0.006738;
788 func_blur[2][5] = 0.011109;
789 func_blur[2][6] = 0.006738;
790 func_blur[2][7] = 0.001503;
791 func_blur[2][8] = 0.000123;
792 func_blur[2][9] = 0.000004;
793 func_blur[2][10] = 0.000000;
794 func_blur[3][0] = 0.000001;
795 func_blur[3][1] = 0.000045;
796 func_blur[3][2] = 0.001503;
797 func_blur[3][3] = 0.018316;
798 func_blur[3][4] = 0.082085;
799 func_blur[3][5] = 0.135335;
800 func_blur[3][6] = 0.082085;
801 func_blur[3][7] = 0.018316;
802 func_blur[3][8] = 0.001503;
803 func_blur[3][9] = 0.000045;
804 func_blur[3][10] = 0.000001;
805 func_blur[4][0] = 0.000002;
806 func_blur[4][1] = 0.000203;
807 func_blur[4][2] = 0.006738;
808 func_blur[4][3] = 0.082085;
809 func_blur[4][4] = 0.367879;
810 func_blur[4][5] = 0.606531;
811 func_blur[4][6] = 0.367879;
812 func_blur[4][7] = 0.082085;
813 func_blur[4][8] = 0.006738;
814 func_blur[4][9] = 0.000203;
815 func_blur[4][10] = 0.000002;
816 func_blur[5][0] = 0.000004;
817 func_blur[5][1] = 0.000335;
818 func_blur[5][2] = 0.011109;
819 func_blur[5][3] = 0.135335;
820 func_blur[5][4] = 0.606531;
821 func_blur[5][5] = 1.000000;
822 func_blur[5][6] = 0.606531;
823 func_blur[5][7] = 0.135335;
824 func_blur[5][8] = 0.011109;
825 func_blur[5][9] = 0.000335;
826 func_blur[5][10] = 0.000004;
827 func_blur[6][0] = 0.000002;
828 func_blur[6][1] = 0.000203;
829 func_blur[6][2] = 0.006738;
830 func_blur[6][3] = 0.082085;
831 func_blur[6][4] = 0.367879;
832 func_blur[6][5] = 0.606531;
833 func_blur[6][6] = 0.367879;
834 func_blur[6][7] = 0.082085;
835 func_blur[6][8] = 0.006738;
836 func_blur[6][9] = 0.000203;
837 func_blur[6][10] = 0.000002;
838 func_blur[7][0] = 0.000001;
839 func_blur[7][1] = 0.000045;
840 func_blur[7][2] = 0.001503;
841 func_blur[7][3] = 0.018316;
842 func_blur[7][4] = 0.082085;
843 func_blur[7][5] = 0.135335;
844 func_blur[7][6] = 0.082085;
845 func_blur[7][7] = 0.018316;
846 func_blur[7][8] = 0.001503;
847 func_blur[7][9] = 0.000045;
848 func_blur[7][10] = 0.000001;
849 func_blur[8][0] = 0.000000;
850 func_blur[8][1] = 0.000004;
851 func_blur[8][2] = 0.000123;
852 func_blur[8][3] = 0.001503;
853 func_blur[8][4] = 0.006738;
854 func_blur[8][5] = 0.011109;
855 func_blur[8][6] = 0.006738;
856 func_blur[8][7] = 0.001503;
857 func_blur[8][8] = 0.000123;
858 func_blur[8][9] = 0.000004;
859 func_blur[8][10] = 0.000000;
860 func_blur[9][0] = 0.000000;
861 func_blur[9][1] = 0.000000;
862 func_blur[9][2] = 0.000004;
863 func_blur[9][3] = 0.000045;
864 func_blur[9][4] = 0.000203;
865 func_blur[9][5] = 0.000335;
866 func_blur[9][6] = 0.000203;
867 func_blur[9][7] = 0.000045;
868 func_blur[9][8] = 0.000004;
869 func_blur[9][9] = 0.000000;
870 func_blur[9][10] = 0.000000;
871 func_blur[10][0] = 0.000000;
872 func_blur[10][1] = 0.000000;
873 func_blur[10][2] = 0.000000;
874 func_blur[10][3] = 0.000001;
875 func_blur[10][4] = 0.000002;
876 func_blur[10][5] = 0.000004;
877 func_blur[10][6] = 0.000002;
878 func_blur[10][7] = 0.000001;
879 func_blur[10][8] = 0.000000;
880 func_blur[10][9] = 0.000000;
881 func_blur[10][10] = 0.000000;
883 double temp_integral_x = 0;
884 double temp_integral_y = 0;
890 <<
"WARNING...BlurNeighborhoods>10 not currently allowed. Shrinking to 10.";
894 TH2F* h_clone_derivative_x = (TH2F*)h_derivative_x.Clone(
"h_clone_derivative_x");
895 TH2F* h_clone_derivative_y = (TH2F*)h_derivative_y.Clone(
"h_clone_derivative_y");
900 for (
int ix = 1; ix <= h_derivative_x.GetNbinsX(); ix++) {
901 for (
int iy = 1; iy <= h_derivative_y.GetNbinsY(); iy++) {
911 h_clone_derivative_x->GetBinContent(jx, jy) * func_blur[(ix - jx) + 5][(iy - jy) + 5];
913 h_clone_derivative_y->GetBinContent(jx, jy) * func_blur[(ix - jx) + 5][(iy - jy) + 5];
916 h_derivative_x.SetBinContent(ix, iy, temp_integral_x);
917 h_derivative_y.SetBinContent(ix, iy, temp_integral_y);
921 delete h_clone_derivative_x;
922 delete h_clone_derivative_y;
931 TH2F
const& h_derivative_y,
934 const int x_bins = h_derivative_x.GetNbinsX();
935 const int y_bins = h_derivative_y.GetNbinsY();
938 double st_xx = 0., st_xy = 0., st_yy = 0.;
952 st_xx += h_derivative_x.GetBinContent(jx, jy) * h_derivative_x.GetBinContent(jx, jy);
953 st_yy += h_derivative_y.GetBinContent(jx, jy) * h_derivative_y.GetBinContent(jx, jy);
954 st_xy += h_derivative_x.GetBinContent(jx, jy) * h_derivative_y.GetBinContent(jx, jy);
963 st_xx -= h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
964 h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
965 st_xx += h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy) *
966 h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy);
968 st_yy -= h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
969 h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
970 st_yy += h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy) *
971 h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy);
973 st_xy -= h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
974 h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
975 st_xy += h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy) *
976 h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy);
981 h_cornerScore.SetBinContent(
985 h_cornerScore.SetBinContent(
988 (st_xx * st_yy - st_xy * st_xy) -
1003 TH2D
const& h_cornerScore,
1004 std::vector<geo::WireID> wireIDs,
1006 TH2D& h_maxSuppress,
1010 std::vector<recob::EndPoint2D> corner_vector;
1011 const int x_bins = h_cornerScore.GetNbinsX();
1012 const int y_bins = h_cornerScore.GetNbinsY();
1014 for (
int iy = 1; iy <= y_bins; iy++) {
1015 for (
int ix = 1; ix <= x_bins; ix++) {
1019 double temp_max = -1000;
1020 double temp_center_bin =
false;
1025 if (h_cornerScore.GetBinContent(jx, jy) > temp_max) {
1026 temp_max = h_cornerScore.GetBinContent(jx, jy);
1027 if (jx == ix && jy == iy)
1028 temp_center_bin =
true;
1030 temp_center_bin =
false;
1036 if (temp_center_bin) {
1043 time_tick, wireIDs[wire_number], h_cornerScore.GetBinContent(ix, iy), id, view, totalQ);
1044 corner_vector.push_back(
corner);
1046 h_maxSuppress.SetBinContent(ix, iy, h_cornerScore.GetBinContent(ix, iy));
1050 return corner_vector;
1059 float threshold)
const 1061 int x1 = hist.GetXaxis()->FindBin(begin_x);
1062 int y1 = hist.GetYaxis()->FindBin(begin_y);
1063 int x2 = hist.GetXaxis()->FindBin(end_x);
1064 int y2 = hist.GetYaxis()->FindBin(end_y);
1066 if (x1 == x2 &&
abs(y1 - y2) < 1
e-5)
return 0;
1074 int bin_counter = 0;
1078 float slope = (y2 -
y1) / ((
float)(x2 -
x1));
1080 for (
int ix = x1; ix <=
x2; ix++) {
1085 y_min = y1 + slope * (ix -
x1);
1086 y_max = y1 + slope * (ix + 1 -
x1);
1089 y_max = (y1 + 1) + slope * (ix - x1);
1090 y_min = (y1 + 1) + slope * (ix + 1 - x1);
1093 for (
int iy = y_min; iy <= y_max; iy++) {
1096 if (hist.GetBinContent(ix, iy) > threshold) fraction += 1.;
1103 for (
int iy = y_min; iy <= y_max; iy++) {
1105 if (hist.GetBinContent(x1, iy) > threshold) fraction += 1.;
1109 return fraction / bin_counter;
1115 TH2F
const& h_wire_data,
1116 std::vector<recob::EndPoint2D>
const& corner_vector,
1117 std::vector<recob::EndPoint2D>& corner_lineIntegralScore_vector,
1118 TH2F& h_lineIntegralScore)
const 1120 for (
auto const i_corner : corner_vector) {
1123 for (
auto const j_corner : corner_vector) {
1126 i_corner.WireID().Wire,
1127 i_corner.DriftTime(),
1128 j_corner.WireID().Wire,
1129 j_corner.DriftTime(),
1142 corner_lineIntegralScore_vector.push_back(
corner);
1144 h_lineIntegralScore.SetBinContent(h_wire_data.GetXaxis()->FindBin(i_corner.WireID().Wire),
1145 h_wire_data.GetYaxis()->FindBin(i_corner.DriftTime()),
std::string fDerivative_BlurFunc
std::vector< std::vector< geo::WireID > > WireData_IDs
void create_image_histo(TH2F const &h_wire_data, TH2F &h_conversion) const
std::string fConversion_func
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
void create_cornerScore_histogram(TH2F const &h_derivative_x, TH2F const &h_derivative_y, TH2D &h_cornerScore)
Encapsulate the construction of a single cyostat .
Float_t y1[n_points_granero]
std::vector< std::tuple< int, TH2F, int, int > > WireData_trimmed_histos
std::pair< float, float > minmax(const float a, const float b)
minmax
std::string fCornerScore_algorithm
std::vector< recob::EndPoint2D > perform_maximum_suppression(TH2D const &h_cornerScore, std::vector< geo::WireID > wireIDs, geo::View_t view, TH2D &h_maxSuppress, int startx=0, int starty=0) const
float fIntegral_fraction_threshold
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
int fDerivative_neighborhood
Float_t x1[n_points_granero]
void attach_feature_points(TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &, int startx=0, int starty=0)
int fCornerScore_neighborhood
The data type to uniquely identify a Plane.
std::string fCalDataModuleLabel
constexpr auto abs(T v)
Returns the absolute value of the argument.
WireID_t Wire
Index of the wire within its plane.
std::string fConversion_algorithm
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Geometry information for a single cryostat.
TH2F const & GetWireDataHist(unsigned int) const
void InitializeGeometry(geo::WireReadoutGeom const &)
float fIntegral_bin_threshold
void create_derivative_histograms(TH2F const &h_conversion, TH2F &h_derivative_x, TH2F &h_derivative_y)
Float_t y2[n_points_geant4]
void calculate_line_integral_score(TH2F const &h_wire_data, std::vector< recob::EndPoint2D > const &corner_vector, std::vector< recob::EndPoint2D > &corner_lineIntegralScore_vector, TH2F &h_lineIntegralScore) const
Access the description of the physical detector geometry.
View_t View() const
Which coordinate does this plane measure.
algorithm to find feature 2D points
Interface for a class providing readout channel mapping to geometry.
CornerFinderAlg(fhicl::ParameterSet const &pset)
float fCornerScore_Noble_epsilon
T get(std::string const &key) const
std::vector< TH1D > WireData_histos_ProjectionY
void get_feature_points(std::vector< recob::EndPoint2D > &, geo::WireReadoutGeom const &)
float line_integral(TH2F const &hist, int x1, float y1, int x2, float y2, float threshold) const
std::vector< TH2F > WireData_histos
void get_feature_points_LineIntegralScore(std::vector< recob::EndPoint2D > &, geo::WireReadoutGeom const &)
std::string fDerivative_method
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
Description of the physical geometry of one entire detector.
void GrabWires(std::vector< recob::Wire > const &wireVec, geo::WireReadoutGeom const &)
void create_smaller_histos(geo::WireReadoutGeom const &)
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
std::vector< TH1D > WireData_histos_ProjectionX
int fConversion_bins_per_input_x
float fCornerScore_Harris_kappa
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
Encapsulate the construction of a single detector plane .
unsigned int event_number
void get_feature_points_fast(std::vector< recob::EndPoint2D > &, geo::GeometryCore const &, geo::WireReadoutGeom const &)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
range_type< T > Iterate() const
float fConversion_threshold
int fDerivative_BlurNeighborhood
int fMaxSuppress_neighborhood
void attach_feature_points_LineIntegralScore(TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &)
double fTrimming_totalThreshold
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
int fMaxSuppress_threshold
Float_t x2[n_points_geant4]
int fConversion_bins_per_input_y
int fConversion_func_neighborhood
Interface to geometry for wire readouts .
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane .
float fTrimming_threshold