491 def Exec(self, opt) :
492 self.kalman_tracks.Clear('C')
493
494
495 if self.scale>1 and self.standalone:
496 if ROOT.gRandom.Rndm() > 1.0/self.scale: return
497
498 self.events_run += 1
499 hit_collection = {"pos" : [[], [], []],
500 "d" : [[], [], []],
501 "vert" : [],
502 "index" : [],
503 "system" : [],
504 "detectorID" : [],
505 "B" : [[], [], []],
506 "time": [],
507 "mask": []}
508
509 if ("us" in self.hits_to_fit) or ("ds" in self.hits_to_fit) or ("ve" in self.hits_to_fit) :
510
511 for i_hit, muFilterHit in enumerate(self.MuFilterHits) :
512
513 if muFilterHit.GetSystem() == 1 :
514 if "ve" not in self.hits_to_fit :
515 continue
516 elif muFilterHit.GetSystem() == 2 :
517 if "us" not in self.hits_to_fit :
518 continue
519 elif muFilterHit.GetSystem() == 3 :
520 if "ds" not in self.hits_to_fit :
521 continue
522 else :
523 if self.logger.IsLogNeeded(ROOT.fair.Severity.warn):
524 print("WARNING! Unknown MuFilter system!!")
525
526 self.mufiDet.GetPosition(muFilterHit.GetDetectorID(), self.a, self.b)
527
528 hit_collection["pos"][0].append(self.a.X())
529 hit_collection["pos"][1].append(self.a.Y())
530 hit_collection["pos"][2].append(self.a.Z())
531
532 hit_collection["B"][0].append(self.b.X())
533 hit_collection["B"][1].append(self.b.Y())
534 hit_collection["B"][2].append(self.b.Z())
535
536 hit_collection["vert"].append(muFilterHit.isVertical())
537 hit_collection["system"].append(muFilterHit.GetSystem())
538
539 hit_collection["d"][0].append(self.MuFilter_ds_dx)
540 hit_collection["d"][2].append(self.MuFilter_ds_dz)
541
542 hit_collection["index"].append(i_hit)
543
544 hit_collection["detectorID"].append(muFilterHit.GetDetectorID())
545 hit_collection["mask"].append(False)
546
547 times = []
548
549 if muFilterHit.GetSystem() == 3 :
550 hit_collection["d"][1].append(self.MuFilter_ds_dx)
551 for ch in range(self.MuFilter_ds_nSiPMs_hor):
552 if muFilterHit.isVertical() and ch==self.MuFilter_ds_nSiPMs_vert: break
553 if self.isMC:
554 times.append(muFilterHit.GetAllTimes()[ch])
555 else: times.append(muFilterHit.GetAllTimes()[ch]*6.25)
556
557 else :
558 hit_collection["d"][1].append(self.MuFilter_us_dy)
559 for ch in range(self.MuFilter_us_nSiPMs):
560 if self.isMC:
561 times.append(muFilterHit.GetAllTimes()[ch])
562 else: times.append(muFilterHit.GetAllTimes()[ch]*6.25)
563 hit_collection["time"].append(times)
564
565 if "sf" in self.hits_to_fit :
566 if self.Scifi_meas:
567
568 self.clusScifi.Clear()
569 self.scifiCluster()
570
571
572 for i_clust, scifiCl in enumerate(self.clusScifi) :
573 scifiCl.GetPosition(self.a,self.b)
574
575 hit_collection["pos"][0].append(self.a.X())
576 hit_collection["pos"][1].append(self.a.Y())
577 hit_collection["pos"][2].append(self.a.Z())
578
579 hit_collection["B"][0].append(self.b.X())
580 hit_collection["B"][1].append(self.b.Y())
581 hit_collection["B"][2].append(self.b.Z())
582
583
584 hit_collection["d"][0].append(scifiCl.GetN()*self.Scifi_dx)
585 hit_collection["d"][1].append(scifiCl.GetN()*self.Scifi_dy)
586 hit_collection["d"][2].append(self.Scifi_dz)
587
588 if int(scifiCl.GetFirst()/100000)%10==1:
589 hit_collection["vert"].append(True)
590 else: hit_collection["vert"].append(False)
591 hit_collection["index"].append(i_clust)
592
593 hit_collection["system"].append(0)
594 hit_collection["detectorID"].append(scifiCl.GetFirst())
595 hit_collection["mask"].append(False)
596
597 times = []
598 if self.isMC : times.append(scifiCl.GetTime()/6.25)
599 else: times.append(scifiCl.GetTime())
600 hit_collection["time"].append(times)
601
602 else:
603 if self.hits_for_triplet == 'sf' and self.hits_to_fit == 'sf':
604
605 N_plane_ZY = {1:0, 2:0, 3:0, 4:0, 5:0}
606 N_plane_ZX = {1:0, 2:0, 3:0, 4:0, 5:0}
607 for scifiHit in self.ScifiHits:
608 if not scifiHit.isValid(): continue
609 if scifiHit.isVertical():
610 N_plane_ZX[scifiHit.GetStation()] += 1
611 else:
612 N_plane_ZY[scifiHit.GetStation()] += 1
613 if self.mask_plane:
614 mask_plane_ZY = []
615 mask_plane_ZX = []
616
617 N_plane_ZY = dict(sorted(N_plane_ZY.items(), key=lambda item: item[1], reverse = True))
618 N_plane_ZX = dict(sorted(N_plane_ZX.items(), key=lambda item: item[1], reverse = True))
619
620 n_zx = self.Scifi_nPlanes - list(N_plane_ZX.values()).count(0)
621 n_zy = self.Scifi_nPlanes - list(N_plane_ZY.values()).count(0)
622
623 if n_zx < self.min_planes_hit or n_zy < self.min_planes_hit: return
624
625 for ii in range(n_zx-self.min_planes_hit):
626 if list(N_plane_ZX.values())[ii] > self.Nhits_per_plane:
627 mask_plane_ZX.append(list(N_plane_ZX.keys())[ii])
628 for ii in range(n_zy-self.min_planes_hit):
629 if list(N_plane_ZY.values())[ii] > self.Nhits_per_plane:
630 mask_plane_ZY.append(list(N_plane_ZY.keys())[ii])
631
632
633 for i_hit, scifiHit in enumerate(self.ScifiHits) :
634 if not scifiHit.isValid(): continue
635 self.scifiDet.GetSiPMPosition(scifiHit.GetDetectorID(), self.a, self.b)
636 hit_collection["pos"][0].append(self.a.X())
637 hit_collection["pos"][1].append(self.a.Y())
638 hit_collection["pos"][2].append(self.a.Z())
639
640 hit_collection["B"][0].append(self.b.X())
641 hit_collection["B"][1].append(self.b.Y())
642 hit_collection["B"][2].append(self.b.Z())
643
644 hit_collection["d"][0].append(self.Scifi_dx)
645 hit_collection["d"][1].append(self.Scifi_dy)
646 hit_collection["d"][2].append(self.Scifi_dz)
647
648 hit_collection["vert"].append(scifiHit.isVertical())
649 hit_collection["index"].append(i_hit)
650
651 hit_collection["system"].append(0)
652
653 hit_collection["detectorID"].append(scifiHit.GetDetectorID())
654
655 if self.hits_for_triplet == 'sf' and self.hits_to_fit == 'sf' and self.mask_plane:
656 if (scifiHit.isVertical()==0 and scifiHit.GetStation() in mask_plane_ZY) or (scifiHit.isVertical() and scifiHit.GetStation() in mask_plane_ZX):
657 hit_collection["mask"].append(True)
658 else: hit_collection["mask"].append(False)
659 else:
660 hit_collection["mask"].append(False)
661
662 times = []
663
664 if self.isMC : times.append(scifiHit.GetTime())
665 else: times.append(scifiHit.GetTime()*6.25)
666 hit_collection["time"].append(times)
667
668
669 if len(hit_collection['pos'][0])==0: return
670
671
672 for key, item in hit_collection.items() :
673 if key == 'vert' :
674 this_dtype = np.bool_
675 elif key == 'mask' :
676 this_dtype = np.bool_
677 elif key == "index" or key == "system" or key == "detectorID" :
678 this_dtype = np.int32
679 elif key != 'time' :
680 this_dtype = np.float32
681 if key== 'time':
682 length = max(map(len, item))
683 hit_collection[key] = np.stack(np.array([xi+[None]*(length-len(xi)) for xi in item]), axis = 1)
684 else: hit_collection[key] = np.array(item, dtype = this_dtype)
685
686
687 triplet_condition_system = []
688 if "sf" in self.hits_for_triplet :
689 triplet_condition_system.append(0)
690 if "ve" in self.hits_for_triplet :
691 triplet_condition_system.append(1)
692 if "us" in self.hits_for_triplet :
693 triplet_condition_system.append(2)
694 if "ds" in self.hits_for_triplet :
695 triplet_condition_system.append(3)
696
697
698 for i_muon in range(self.max_reco_muons) :
699
700 triplet_hits_horizontal = np.logical_and( ~hit_collection["vert"],
701 np.isin(hit_collection["system"], triplet_condition_system) )
702 triplet_hits_vertical = np.logical_and( hit_collection["vert"],
703 np.isin(hit_collection["system"], triplet_condition_system) )
704
705 n_planes_ZY = numPlanesHit(hit_collection["system"][triplet_hits_horizontal],
706 hit_collection["detectorID"][triplet_hits_horizontal])
707 if n_planes_ZY < self.min_planes_hit :
708 break
709
710 n_planes_ZX = numPlanesHit(hit_collection["system"][triplet_hits_vertical],
711 hit_collection["detectorID"][triplet_hits_vertical])
712 if n_planes_ZX < self.min_planes_hit :
713 break
714
715
716 muon_hits_horizontal = np.logical_and( np.logical_and( ~hit_collection["vert"], ~hit_collection["mask"]),
717 np.isin(hit_collection["system"], [1, 2, 3]))
718 muon_hits_vertical = np.logical_and( np.logical_and( hit_collection["vert"], ~hit_collection["mask"]),
719 np.isin(hit_collection["system"], [1, 2, 3]))
720 scifi_hits_horizontal = np.logical_and( np.logical_and( ~hit_collection["vert"], ~hit_collection["mask"]),
721 np.isin(hit_collection["system"], [0]))
722 scifi_hits_vertical = np.logical_and( np.logical_and( hit_collection["vert"], ~hit_collection["mask"]),
723 np.isin(hit_collection["system"], [0]))
724
725
726 ZY = np.dstack([np.concatenate([np.tile(hit_collection["pos"][2][muon_hits_horizontal], self.muon_weight),
727 hit_collection["pos"][2][scifi_hits_horizontal]]),
728 np.concatenate([np.tile(hit_collection["pos"][1][muon_hits_horizontal], self.muon_weight),
729 hit_collection["pos"][1][scifi_hits_horizontal]])])[0]
730
731 d_ZY = np.dstack([np.concatenate([np.tile(hit_collection["d"][2][muon_hits_horizontal], self.muon_weight),
732 hit_collection["d"][2][scifi_hits_horizontal]]),
733 np.concatenate([np.tile(hit_collection["d"][1][muon_hits_horizontal], self.muon_weight),
734 hit_collection["d"][1][scifi_hits_horizontal]])])[0]
735
736 ZX = np.dstack([np.concatenate([np.tile(hit_collection["pos"][2][muon_hits_vertical], self.muon_weight),
737 hit_collection["pos"][2][scifi_hits_vertical]]),
738 np.concatenate([np.tile(hit_collection["pos"][0][muon_hits_vertical], self.muon_weight),
739 hit_collection["pos"][0][scifi_hits_vertical]])])[0]
740
741 d_ZX = np.dstack([np.concatenate([np.tile(hit_collection["d"][2][muon_hits_vertical], self.muon_weight),
742 hit_collection["d"][2][scifi_hits_vertical]]),
743 np.concatenate([np.tile(hit_collection["d"][0][muon_hits_vertical], self.muon_weight),
744 hit_collection["d"][0][scifi_hits_vertical]])])[0]
745
746 is_scaled = False
747 ZY_hough = self.h_ZY.fit_randomize(ZY, d_ZY, self.n_random, is_scaled, self.draw)
748 ZX_hough = self.h_ZX.fit_randomize(ZX, d_ZX, self.n_random, is_scaled, self.draw)
749
750 tol = self.tolerance
751
752
753 if len(hit_collection["detectorID"]) <= self.max_n_Scifi_hits and self.hits_for_triplet == 'sf' and self.hits_to_fit == 'sf' :
754
755 if max(N_plane_ZX.values()) <= self.max_n_hits_plane and max(N_plane_ZY.values()) <= self.max_n_hits_plane :
756 tol = 5*self.tolerance
757
758 if self.hits_for_triplet == 'ds' and self.hits_to_fit == 'ds' :
759
760 N_plane_ZY = {0:0, 1:0, 2:0, 3:0}
761 N_plane_ZX = {0:0, 1:0, 2:0, 3:0}
762 for item in range(len(hit_collection["detectorID"])):
763 if hit_collection["vert"][item]: N_plane_ZX[(hit_collection["detectorID"][item]%10000)//1000] += 1
764 else: N_plane_ZY[(hit_collection["detectorID"][item]%10000)//1000] += 1
765 if max(N_plane_ZX.values()) <= self.max_n_hits_plane and max(N_plane_ZY.values()) <= self.max_n_hits_plane :
766 tol = 3*self.tolerance
767
768
769 track_hits_for_triplet_ZY = hit_finder(ZY_hough[0], ZY_hough[1],
770 np.dstack([hit_collection["pos"][2][triplet_hits_horizontal],
771 hit_collection["pos"][1][triplet_hits_horizontal]]),
772 np.dstack([hit_collection["d"][2][triplet_hits_horizontal],
773 hit_collection["d"][1][triplet_hits_horizontal]]), tol)
774
775 track_hits_for_triplet_ZX = hit_finder(ZX_hough[0], ZX_hough[1],
776 np.dstack([hit_collection["pos"][2][triplet_hits_vertical],
777 hit_collection["pos"][0][triplet_hits_vertical]]),
778 np.dstack([hit_collection["d"][2][triplet_hits_vertical],
779 hit_collection["d"][0][triplet_hits_vertical]]), tol)
780
781 n_planes_hit_ZY = numPlanesHit(hit_collection["system"][triplet_hits_horizontal][track_hits_for_triplet_ZY],
782 hit_collection["detectorID"][triplet_hits_horizontal][track_hits_for_triplet_ZY])
783
784 n_planes_hit_ZX = numPlanesHit(hit_collection["system"][triplet_hits_vertical][track_hits_for_triplet_ZX],
785 hit_collection["detectorID"][triplet_hits_vertical][track_hits_for_triplet_ZX])
786
787
788 if (self.hits_to_fit == 'sf' and len(hit_collection["detectorID"]) <= self.max_n_Scifi_hits and \
789 (n_planes_hit_ZY < self.min_planes_hit or n_planes_hit_ZX < self.min_planes_hit)):
790
791 is_scaled = True
792 ZY_hough = self.h_ZY.fit_randomize(ZY, d_ZY, self.n_random, is_scaled, self.draw)
793 ZX_hough = self.h_ZX.fit_randomize(ZX, d_ZX, self.n_random, is_scaled, self.draw)
794
795
796 track_hits_for_triplet_ZY = hit_finder(ZY_hough[0], ZY_hough[1],
797 np.dstack([hit_collection["pos"][2][triplet_hits_horizontal],
798 hit_collection["pos"][1][triplet_hits_horizontal]]),
799 np.dstack([hit_collection["d"][2][triplet_hits_horizontal],
800 hit_collection["d"][1][triplet_hits_horizontal]]), tol)
801
802 n_planes_hit_ZY = numPlanesHit(hit_collection["system"][triplet_hits_horizontal][track_hits_for_triplet_ZY],
803 hit_collection["detectorID"][triplet_hits_horizontal][track_hits_for_triplet_ZY])
804 if n_planes_hit_ZY < self.min_planes_hit:
805 break
806
807 track_hits_for_triplet_ZX = hit_finder(ZX_hough[0], ZX_hough[1],
808 np.dstack([hit_collection["pos"][2][triplet_hits_vertical],
809 hit_collection["pos"][0][triplet_hits_vertical]]),
810 np.dstack([hit_collection["d"][2][triplet_hits_vertical],
811 hit_collection["d"][0][triplet_hits_vertical]]), tol)
812 n_planes_hit_ZX = numPlanesHit(hit_collection["system"][triplet_hits_vertical][track_hits_for_triplet_ZX],
813 hit_collection["detectorID"][triplet_hits_vertical][track_hits_for_triplet_ZX])
814
815 if n_planes_hit_ZY < self.min_planes_hit or n_planes_hit_ZX < self.min_planes_hit:
816 break
817
818
819
820
821
822 track_hits_ZY = hit_finder(ZY_hough[0], ZY_hough[1],
823 np.dstack([hit_collection["pos"][2][~hit_collection["vert"]],
824 hit_collection["pos"][1][~hit_collection["vert"]]]),
825 np.dstack([hit_collection["d"][2][~hit_collection["vert"]],
826 hit_collection["d"][1][~hit_collection["vert"]]]), tol)
827
828 track_hits_ZX = hit_finder(ZX_hough[0], ZX_hough[1],
829 np.dstack([hit_collection["pos"][2][hit_collection["vert"]],
830 hit_collection["pos"][0][hit_collection["vert"]]]),
831 np.dstack([hit_collection["d"][2][hit_collection["vert"]],
832 hit_collection["d"][0][hit_collection["vert"]]]), tol)
833
834 posM = ROOT.TVector3(0, 0, 0.)
835 momM = ROOT.TVector3(0,0,100.)
836
837
838 covM = ROOT.TMatrixDSym(6)
839 if self.hits_to_fit.find('sf') >= 0:
840 res = self.kalman_sigmaScifi_spatial
841 if self.hits_to_fit == 'ds':
842 res = self.kalman_sigmaMufiDS_spatial
843 for i in range(3): covM[i][i] = res*res
844 for i in range(3,6): covM[i][i] = ROOT.TMath.Power(res / (4.*2.) / ROOT.TMath.Sqrt(3), 2)
845 rep = ROOT.genfit.RKTrackRep(13)
846
847
848 state = ROOT.genfit.MeasuredStateOnPlane(rep)
849 rep.setPosMomCov(state, posM, momM, covM)
850
851
852 seedState = ROOT.TVectorD(6)
853 seedCov = ROOT.TMatrixDSym(6)
854 rep.get6DStateCov(state, seedState, seedCov)
855
856 theTrack = ROOT.genfit.Track(rep, seedState, seedCov)
857
858
859 hit_z = np.concatenate([hit_collection["pos"][2][hit_collection["vert"]][track_hits_ZX],
860 hit_collection["pos"][2][~hit_collection["vert"]][track_hits_ZY]])
861
862 hit_A0 = np.concatenate([hit_collection["pos"][0][hit_collection["vert"]][track_hits_ZX],
863 hit_collection["pos"][0][~hit_collection["vert"]][track_hits_ZY]])
864
865 hit_A1 = np.concatenate([hit_collection["pos"][1][hit_collection["vert"]][track_hits_ZX],
866 hit_collection["pos"][1][~hit_collection["vert"]][track_hits_ZY]])
867
868 hit_B0 = np.concatenate([hit_collection["B"][0][hit_collection["vert"]][track_hits_ZX],
869 hit_collection["B"][0][~hit_collection["vert"]][track_hits_ZY]])
870
871 hit_B1 = np.concatenate([hit_collection["B"][1][hit_collection["vert"]][track_hits_ZX],
872 hit_collection["B"][1][~hit_collection["vert"]][track_hits_ZY]])
873
874 hit_B2 = np.concatenate([hit_collection["B"][2][hit_collection["vert"]][track_hits_ZX],
875 hit_collection["B"][2][~hit_collection["vert"]][track_hits_ZY]])
876
877 hit_detid = np.concatenate([hit_collection["detectorID"][hit_collection["vert"]][track_hits_ZX],
878 hit_collection["detectorID"][~hit_collection["vert"]][track_hits_ZY]])
879
880 kalman_spatial_sigma = np.concatenate([hit_collection["d"][0][hit_collection["vert"]][track_hits_ZX] / 12**0.5,
881 hit_collection["d"][1][~hit_collection["vert"]][track_hits_ZY] / 12**0.5])
882
883
884 kalman_max_dis = np.concatenate([((hit_collection["d"][0][hit_collection["vert"]][track_hits_ZX]/2.)**2 +
885 (hit_collection["d"][2][hit_collection["vert"]][track_hits_ZX]/2.)**2)**0.5,
886 ((hit_collection["d"][1][~hit_collection["vert"]][track_hits_ZY]/2.)**2 +
887 (hit_collection["d"][2][~hit_collection["vert"]][track_hits_ZY]/2.)**2)**0.5])
888
889 hitID = 0
890
891 hit_time = {}
892 for ch in range(hit_collection["time"].shape[0]):
893 hit_time[ch] = np.concatenate([hit_collection["time"][ch][hit_collection["vert"]][track_hits_ZX],
894 hit_collection["time"][ch][~hit_collection["vert"]][track_hits_ZY]])
895
896 for i_z_sorted in hit_z.argsort() :
897 tp = ROOT.genfit.TrackPoint()
898 hitCov = ROOT.TMatrixDSym(7)
899 hitCov[6][6] = kalman_spatial_sigma[i_z_sorted]**2
900
901 measurement = ROOT.genfit.WireMeasurement(ROOT.TVectorD(7, array('d', [hit_A0[i_z_sorted],
902 hit_A1[i_z_sorted],
903 hit_z[i_z_sorted],
904 hit_B0[i_z_sorted],
905 hit_B1[i_z_sorted],
906 hit_B2[i_z_sorted],
907 0.])),
908 hitCov,
909 1,
910 6,
911 tp)
912
913 measurement.setMaxDistance(kalman_max_dis[i_z_sorted])
914 measurement.setDetId(int(hit_detid[i_z_sorted]))
915 measurement.setHitId(int(hitID))
916 hitID += 1
917 tp.addRawMeasurement(measurement)
918 theTrack.insertPoint(tp)
919
920 if not theTrack.checkConsistency():
921 theTrack.Delete()
922 raise RuntimeError("Kalman fitter track consistency check failed.")
923
924
925 self.kalman_fitter.processTrack(theTrack)
926
927 fitStatus = theTrack.getFitStatus()
928 if not fitStatus.isFitConverged() and 0>1:
929 theTrack.Delete()
930 raise RuntimeError("Kalman fit did not converge.")
931
932
933 theTrack.SetUniqueID(self.track_type)
934 if fitStatus.isFitConverged():
935 if self.genfitTrack: self.kalman_tracks.Add(theTrack)
936 else :
937
938 this_track = ROOT.sndRecoTrack(theTrack)
939 pointTimes = ROOT.std.vector(ROOT.std.vector('float'))()
940 for n, i_z_sorted in enumerate(hit_z.argsort()):
941 t_per_hit = []
942 for ch in range(len(hit_time)):
943 if hit_time[ch][i_z_sorted] != None:
944 t_per_hit.append(hit_time[ch][i_z_sorted])
945 pointTimes.push_back(t_per_hit)
946 this_track.setRawMeasTimes(pointTimes)
947 this_track.setTrackType(self.track_type)
948
949 self.kalman_tracks[i_muon] = this_track
950
951 theTrack.Delete()
952
953
954
955 index_to_remove_ZX = np.where(np.in1d(hit_collection["detectorID"], hit_collection["detectorID"][hit_collection["vert"]][track_hits_ZX]))[0]
956 index_to_remove_ZY = np.where(np.in1d(hit_collection["detectorID"], hit_collection["detectorID"][~hit_collection["vert"]][track_hits_ZY]))[0]
957
958 index_to_remove = np.concatenate([index_to_remove_ZX, index_to_remove_ZY])
959
960
961 for key in hit_collection.keys() :
962 if len(hit_collection[key].shape) == 1 :
963 hit_collection[key] = np.delete(hit_collection[key], index_to_remove)
964 elif len(hit_collection[key].shape) == 2 :
965 hit_collection[key] = np.delete(hit_collection[key], index_to_remove, axis = 1)
966 else :
967 raise Exception("Wrong number of dimensions found when deleting hits in iterative muon identification algorithm.")
968