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