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