39def plot_overlap(scifi, plot_outfile):
40 sGeo = ROOT.gGeoManager
41 scifi.SiPMOverlap()
42 scifi.SiPMmapping()
43
44 SiPMmap = scifi.GetSiPMmap()
45 FibresMap = scifi.GetFibresMap()
46 print("GetSiPMmap",len(SiPMmap))
47 print("GetFibresMap",len(FibresMap))
48
49
50
51
52 fiber_weight = getFiberWeight(FibresMap)
53
54 print("number of fibers in map", len(fiber_weight))
55
56
57
58
59
60
61
62
63 y_values = []
64 z_values = []
65
66 fig, ax = plt.subplots(figsize=(10, 10))
67 SiPMmapVol = sGeo.FindVolumeFast("SiPMmapVol")
68 for sipm in SiPMmapVol.GetNodes():
69 t0 = sipm.GetMatrix().GetTranslation()[0]
70 t1 = sipm.GetMatrix().GetTranslation()[1]
71 t2 = sipm.GetMatrix().GetTranslation()[2]
72
73 y0 = t1
74 z0 = t2
75
76 DX = sipm.GetVolume().GetShape().GetDX()
77 DY = sipm.GetVolume().GetShape().GetDY()
78 DZ = sipm.GetVolume().GetShape().GetDZ()
79
80
81
82 rect = Rectangle((y0 - DY, z0 - DZ), DY*2, DZ*2, edgecolor='blue', facecolor='none')
83 ax.add_patch(rect)
84
85 y_values.extend([y0 - DY, y0 + DY])
86 z_values.extend([z0 - DZ, z0 + DZ])
87
88
89
90 ScifiHorPlaneVol = sGeo.FindVolumeFast("ScifiHorPlaneVol1")
91 imat = 0
92 norm = colors.Normalize(vmin=0, vmax=1.4)
93 cmap = cm.get_cmap('RdYlGn')
94 for mat in ScifiHorPlaneVol.GetNodes():
95 mat_t0 = mat.GetMatrix().GetTranslation()[0]
96 mat_t1 = mat.GetMatrix().GetTranslation()[1]
97 mat_t2 = mat.GetMatrix().GetTranslation()[2]
98
99
100 for fiber in mat.GetVolume().GetNodes():
101 fiber_t0 = fiber.GetMatrix().GetTranslation()[0]
102 fiber_t1 = fiber.GetMatrix().GetTranslation()[1]
103 fiber_t2 = fiber.GetMatrix().GetTranslation()[2]
104
105
106 y0 = fiber_t1 + mat_t1
107 z0 = fiber_t2
108
109 fiber_radius = fiber.GetVolume().GetShape().GetDX()
110
111 fID = fiber.GetNumber()%100000 + imat*1e4
112 if fID in fiber_weight:
113 f_weight = fiber_weight[fID][0]
114 else:
115 f_weight = 0
116 color = cmap(norm(f_weight))
117 circ = Circle((y0, z0), radius=fiber_radius, edgecolor='red', facecolor=color, alpha=0.8)
118 ax.add_patch(circ)
119
120 imat+=1
121
122
123
124 print("SiPM Y axis position limit:", (min(y_values), max(y_values)))
125 print("SiPM Z axis position limit:", (min(z_values), max(z_values)))
126
127 xlim = (min(y_values)-0.1, max(y_values)+0.1)
128 ylim = (-0.09, 0.09)
129
130 fig_width = (xlim[1] - xlim[0]) * 10
131 fig_height = (ylim[1] - ylim[0]) * 10
132 fig.set_size_inches(fig_width, fig_height)
133
134 ax.set_xlim(*xlim)
135 ax.set_ylim(*ylim)
136 ax.set_aspect('equal')
137 ax.set_xlabel("Y axis (cm)")
138 ax.set_ylabel("Z axis (cm)")
139
140
141 yticks = list(ax.get_yticks())
142 yticks.append(0.081)
143 yticks = sorted(
set(yticks))
144 ax.set_yticks(yticks)
145
146 plt.grid(True)
147 plt.savefig(plot_outfile, dpi=300)
148 print(f"plot saved to {plot_outfile}")
149
150
set(INCLUDE_DIRECTORIES ${SYSTEM_INCLUDE_DIRECTORIES} ${VMC_INCLUDE_DIRS} ${CMAKE_SOURCE_DIR}/shipdata ${CMAKE_SOURCE_DIR}/shipLHC ${CMAKE_SOURCE_DIR}/analysis/cuts ${CMAKE_SOURCE_DIR}/analysis/tools ${FMT_INCLUDE_DIR}) include_directories($