643 {
644
645 std::multimap<int, int>
646 y_idx_quad_map;
647 std::multimap<int, int> x_idx_quad_map;
648
649 std::multimap<int, ldmx::TrigScintTrack> y_quad_map;
650 std::multimap<int, ldmx::TrigScintTrack> x_quad_map;
651
652
653 std::map<ldmx::TrigScintTrack, int> y_track_map;
654 std::map<ldmx::TrigScintTrack, int> x_track_map;
655
656 uint trk_idx = -1;
657 for (auto trk : tracks) {
658 trk_idx++;
659
660 if (trk.getCentroidX() == -1) {
661 if (verbose_)
662 ldmx_log(debug) << " -- In matchXYTracks found y track at "
663 << trk.getCentroidY() << "; mapping to quad "
664 << (int)trk.getCentroidY() / (n_bars_y_ / 2)
665 << " with trk index " << trk_idx;
666
667
668 y_quad_map.insert(
669 std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk));
670 y_track_map[trk] = trk_idx;
671 y_idx_quad_map.insert(
672 std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk_idx));
673
674 } else {
675
676 x_quad_map.insert(
677 std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk));
678 x_track_map[trk] = trk_idx;
679 x_idx_quad_map.insert(
680 std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk_idx));
681 if (verbose_)
682 ldmx_log(debug) << " -- In matchXYTracks found x track at (x,y) = ("
683 << trk.getCentroidX() << ", " << trk.getCentroidY()
684 << "); mapping to quad "
685 << (int)trk.getCentroidY() / (n_bars_y_ / 2)
686 << " with trk index " << trk_idx;
687 }
688 }
689
690
691
692
693
694
695
696
697
698 float x0 = 0;
699
700
701 float sx0 = fabs(x_start_);
702 float sx0_vert = fabs(bar_length_y_ / 2);
703
704
705
706 float sy0 = fabs(y_start_) / 4.;
707
708
709
710 for (auto yitr = y_quad_map.begin(); yitr != y_quad_map.end(); ++yitr) {
711 int n_yin_quad = y_quad_map.count((*yitr).first);
712 int n_xin_quad = x_quad_map.count((*yitr).first);
713 float y{-9999.}, sy{-9999.}, x{-9999.}, x1{-9999.}, x2{-9999.}, sx1{-9999.},
714 sx2{-9999.}, y1{-9999.}, y2{-9999.}, sy1{-9999.}, sy2{-9999.};
715
716 float y0 = (((*yitr).first * 8) * y_conv_factor_) + y_start_ + sy0;
717 float sx = 1. / 2 *
718 x_conv_factor_;
719
720
721
722
723 if (n_xin_quad == 0) {
724
725 x = x0;
726 sx = sx0_vert;
727 if (verbose_)
728 ldmx_log(debug) << "\t\t\t no x info in quad " << (*yitr).first
729 << "; will set x to middle of pad, pad half-width as "
730 "precision: set (x, sx)=("
731 << x << ", " << sx << ")";
732 }
733 else if (n_xin_quad ==
734 1) {
735
736
737
738 auto xitr = x_quad_map.find((*yitr).first);
739 x = ((*xitr).second).getCentroidX() * x_conv_factor_ + x_start_;
740
741 if (verbose_)
742 ldmx_log(debug) << "\t\t\t 1 x in quad " << (*yitr).first
743 << ", getting (x, sx)=(" << x << ", " << sx << ")";
744 }
745 else if (n_xin_quad == 2) {
746
747
748
749
750 auto xitr1 = x_quad_map.lower_bound((*yitr).first);
751 auto xitr2 = x_quad_map.upper_bound((*yitr).first);
752 xitr2--;
753
754 if (xitr1 != xitr2) {
755 x1 = ((*xitr1).second).getCentroidX() * x_conv_factor_ + x_start_;
756 x2 = ((*xitr2).second).getCentroidX() * x_conv_factor_ + x_start_;
757 sx1 = x_conv_factor_ / 2.;
758 sx2 = sx1;
759 x = (x1 + x2) / 2.;
760
761 sx = fabs(x1 - x2) / 2;
762 if (verbose_)
763 ldmx_log(debug) << "\t\t -- 2 x in quad: setting y track x "
764 "coordinate to midpoint";
765 }
766 }
767
768 if (n_xin_quad >= 3) {
769 x = x0;
770 sx = sx0;
771 if (verbose_)
772 ldmx_log(debug)
773 << "\t\t\t currently no x info assigned in ambiguous case of "
774 << n_xin_quad << "vertical bar track candidates in quad "
775 << (*yitr).first
776 << "; will set x to middle of pad, pad half-width as "
777 "precision: set (x, sx)=("
778 << x << ", " << sx << ")";
779 }
780
781
782
783 if (n_yin_quad == 1) {
784
785 y = ((*yitr).second).getCentroidY() * y_conv_factor_ + y_start_;
786 sy = ((*yitr).second).getResidual() * y_conv_factor_;
787
788
789 if (sy == 0) sy = 1. / 2 * y_conv_factor_;
790
791 if (n_xin_quad <= 1) {
792
793
794
795 if (n_xin_quad == 1) {
796 auto xidx = x_idx_quad_map.find((*yitr).first);
797 tracks.at((*xidx).second).setPosition(x, y);
798 tracks.at((*xidx).second).setSigmaXY(sx, sy);
799 }
800 if (verbose_)
801 ldmx_log(debug) << "\t\t\t in quad " << (*yitr).first
802 << ", set (x, y) = (" << x << ", " << y
803 << ") and (sx, sy) = " << sx << ", " << sy << ")";
804 auto yidx = y_idx_quad_map.find((*yitr).first);
805 tracks.at((*yidx).second).setPosition(x, y);
806 tracks.at((*yidx).second).setSigmaXY(sx, sy);
807 continue;
808 }
809 }
810
811 if (verbose_)
812 ldmx_log(debug) << "\t\t in quad " << (*yitr).first
813 << ", not single x,y tracks: " << n_xin_quad
814 << " of x and " << n_yin_quad << " of y";
815
816 if (n_yin_quad == 2) {
817
818
819 auto yitr1 = y_quad_map.lower_bound((*yitr).first);
820 auto yitr2 = y_quad_map.upper_bound((*yitr).first);
821 yitr2--;
822 y1 = ((*yitr1).second).getCentroidY() * y_conv_factor_ + y_start_;
823 y2 = ((*yitr2).second).getCentroidY() * y_conv_factor_ + y_start_;
824 sy1 = ((*yitr1).second).getResidual() * y_conv_factor_;
825 sy2 = ((*yitr2).second).getResidual() * y_conv_factor_;
826 if (sy1 == 0) sy1 = 1. / 2 * y_conv_factor_;
827 if (sy2 == 0) sy2 = 1. / 2 * y_conv_factor_;
828 y = (y1 + y2) / 2.;
829 sy = fabs(y1 - y2) / 2;
830 if (verbose_)
831 ldmx_log(debug)
832 << "\t\t -- 2 y in quad: setting x track y coordinate to midpoint";
833 }
834
835 if ((n_xin_quad == 0 || n_xin_quad >= 3) &&
836 (n_yin_quad == 2)) {
837 if (n_xin_quad == 0) {
838 if (verbose_)
839 ldmx_log(debug) << "\t\t -- No x tracks but 2 y tracks in quad: "
840 "unusual behaviour";
841 }
842 auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first);
843 auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first);
844 yidx2--;
845 tracks.at((*yidx1).second).setPosition(x, y1);
846 tracks.at((*yidx1).second).setSigmaXY(sx, sy1);
847 tracks.at((*yidx2).second).setPosition(x, y2);
848 tracks.at((*yidx2).second).setSigmaXY(sx, sy2);
849 continue;
850 }
851
852 if (n_yin_quad == 1 &&
853 n_xin_quad == 2) {
854
855
856
857
858 auto yidx = y_idx_quad_map.find((*yitr).first);
859 tracks.at((*yidx).second).setPosition(x, y);
860 tracks.at((*yidx).second).setSigmaXY(sx, sy);
861
862 int min_overlap_pe = 250;
863 if (((*yitr).second).getPE() < min_overlap_pe) {
864
865
866
867 y = y0;
868 sy = sy0;
869 if (verbose_)
870 ldmx_log(debug) << "\t\t -- Can't tell which x track should be "
871 "matched to single y track. Setting both x track "
872 "coordinates to y quadrant value:";
873 }
874 else if (verbose_)
875 ldmx_log(debug) << "\t\t -- Found large PE count ("
876 << ((*yitr).second).getPE() << " > " << min_overlap_pe
877 << "), suggesting overlap! Setting both x track "
878 "coordinates to y track value:";
879
880
881
882
883
884 if (verbose_)
885 ldmx_log(debug) << "\t\t -- (x1, x2, y) = (" << x1 << ", " << x2
886 << ", " << y << ") and (sx1, sx2, sy) = " << sx1 << ", "
887 << sx2 << ", " << sy << ")";
888
889
890 auto xidx1 = x_idx_quad_map.lower_bound((*yitr).first);
891 auto xidx2 = x_idx_quad_map.upper_bound((*yitr).first);
892 xidx2--;
893 tracks.at((*xidx1).second).setPosition(x1, y);
894 tracks.at((*xidx1).second).setSigmaXY(sx1, sy);
895 tracks.at((*xidx2).second).setPosition(x2, y);
896 tracks.at((*xidx2).second).setSigmaXY(sx2, sy);
897
898 }
899 else if (n_yin_quad == 2 && n_xin_quad == 1) {
900
901
902
903
904 auto xidx = x_idx_quad_map.find((*yitr).first);
905 tracks.at((*xidx).second).setPosition(x, y);
906 tracks.at((*xidx).second).setSigmaXY(sx, sy);
907
908 auto xitr = x_quad_map.lower_bound((*yitr).first);
909 int min_overlap_pe = 300;
910 if (((*xitr).second).getPE() < min_overlap_pe) {
911 if (verbose_)
912 ldmx_log(debug)
913 << "\t\t just 1 x track with not-unusual PE in the quad -- can't "
914 "match; setting mid-point values for x ";
915 x = x0;
916 sx = sx0;
917 }
918 else {
919
920
921
922
923 if (verbose_)
924 ldmx_log(debug) << "\t\t -- Found large PE count ("
925 << ((*xitr).second).getPE() << " > " << min_overlap_pe
926 << ") in x track, suggesting overlap! Setting both y "
927 "track coordinates to x track value:";
928 }
929 if (verbose_)
930 ldmx_log(debug) << "\t\t -- (x, y1, y2) = (" << x << ", " << y1 << ", "
931 << y2 << ") and (sx, sy1, sy2) = " << sx << ", " << sy1
932 << ", " << sy2 << ")";
933
934 auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first);
935 auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first);
936 yidx2--;
937 tracks.at((*yidx1).second).setPosition(x, y1);
938 tracks.at((*yidx1).second).setSigmaXY(sx, sy1);
939 tracks.at((*yidx2).second).setPosition(x, y2);
940 tracks.at((*yidx2).second).setSigmaXY(sx, sy2);
941
942 }
943 else if (n_yin_quad == 2 && n_xin_quad == 2) {
944
945 auto xidx1 = x_idx_quad_map.lower_bound((*yitr).first);
946 auto xidx2 = x_idx_quad_map.upper_bound((*yitr).first);
947 xidx2--;
948 auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first);
949 auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first);
950 yidx2--;
951
952 if (y_idx_quad_map.find((*yitr).first) == y_idx_quad_map.end())
953 ldmx_log(error) << "The two y tracks in the same quadrant at "
954 << (*yitr).first
955 << " appear to not be found in the y track map! "
956 "investigate. Note that yidx1.first = "
957 << (*yidx1).first
958 << " and yidx2.first = " << (*yidx2).first;
959 else {
960 tracks.at((*xidx1).second).setPosition(x1, y);
961 tracks.at((*xidx1).second).setSigmaXY(sx1, sy);
962 tracks.at((*xidx2).second).setPosition(x2, y);
963 tracks.at((*xidx2).second).setSigmaXY(sx2, sy);
964
965 tracks.at((*yidx1).second).setPosition(x, y1);
966 tracks.at((*yidx1).second).setSigmaXY(sx, sy1);
967 tracks.at((*yidx2).second).setPosition(x, y2);
968 tracks.at((*yidx2).second).setSigmaXY(sx, sy2);
969
970 if (verbose_)
971 ldmx_log(debug) << "\t\t -- in a 2 x 2 situaiton; midpoint y: " << y
972 << " for both x tracks, midpoint x: " << x
973 << " for both y tracks";
974 }
975 }
976
977 if (n_xin_quad > 2) {
978 if (verbose_)
979 ldmx_log(debug) << "\t\t -*-*-*- more than 2 x tracks in the same quad "
980 "-- nothing done about the x,y coordinates in this "
981 "situation -- implement if needed!!";
982 }
983 if (n_yin_quad > 2) {
984 if (verbose_)
985 ldmx_log(debug) << "\t\t -*-*-*- more than 2 y tracks in the same quad "
986 "-- nothing done about the x,y coordinates in this "
987 "situation -- implement if needed!!";
988 }
989
990 }
991
992 y_quad_map.clear();
993 x_quad_map.clear();
994
995
996}