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