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