6from math
import exp, sqrt, log, ceil
8import matplotlib.pyplot
as plt
10from ROOT.TMath
import Landau, Gaus, Poisson
11from ShipGeoConfig
import ConfigRegistry
23is_sipm_hor_reversed =
True
24is_sipm_vert_reversed =
True
45ch_max_num = 128*12 - 1
51 (-hw, -hw+edge,
False),
52 (-hw+edge, -hw+edge+charr, (-0.5, 63.5)),
53 (-hw+edge+charr, -hw+edge+charr+gap,
False),
54 (-hw+edge+charr+gap, -hw+edge+array, (63.5, 127.5)),
55 (-hw+edge+array, -hw+edge+array+biggap,
False),
56 (-hw+edge+array+biggap, -hw+edge+array+biggap+charr, (127.5, 191.5)),
57 (-hw+edge+array+biggap+charr, -hw+edge+array+biggap+charr+gap,
False),
58 (-hw+edge+array+biggap+charr+gap, -hw+edge+array+biggap+array, (191.5, 255.5)),
59 (-hw+edge+2*array+biggap, -hw+edge+2*array+2*biggap,
False),
60 (-hw+edge+2*array+2*biggap, -hw+edge+2*array+2*biggap+charr, (255.5, 319.5)),
61 (-hw+edge+2*array+2*biggap+charr, -hw+edge+2*array+2*biggap+charr+gap,
False),
62 (-hw+edge+2*array+2*biggap+charr+gap, -hw+edge+2*array+2*biggap+array, (319.5, 383.5)),
63 (-hw+edge+3*array+2*biggap, -hw+edge+3*array+3*biggap,
False),
64 (-hw+edge+3*array+3*biggap, -hw+edge+3*array+3*biggap+charr, (383.5, 447.5)),
65 (-hw+edge+3*array+3*biggap+charr, -hw+edge+3*array+3*biggap+charr+gap,
False),
66 (-hw+edge+3*array+3*biggap+charr+gap, -hw+edge+3*array+3*biggap+array, (447.5, 511.5)),
68 (-hw+edge+4*array+3*biggap, -hw+edge+4*array+4*biggap,
False),
69 (-hw+edge+4*array+4*biggap, -hw+edge+4*array+4*biggap+charr, (511.5, 575.5)),
70 (-hw+edge+4*array+4*biggap+charr, -hw+edge+4*array+4*biggap+charr+gap,
False),
71 (-hw+edge+4*array+4*biggap+charr+gap, -hw+edge+4*array+4*biggap+array, (575.5, 639.5)),
72 (-hw+edge+5*array+4*biggap, -hw+edge+5*array+5*biggap,
False),
73 (-hw+edge+5*array+5*biggap, -hw+edge+5*array+5*biggap+charr, (639.5, 703.5)),
74 (-hw+edge+5*array+5*biggap+charr, -hw+edge+5*array+5*biggap+charr+gap,
False),
75 (-hw+edge+5*array+5*biggap+charr+gap, -hw+edge+5*array+5*biggap+array, (703.5, 767.5)),
76 (-hw+edge+6*array+5*biggap, -hw+edge+6*array+6*biggap,
False),
77 (-hw+edge+6*array+6*biggap, -hw+edge+6*array+6*biggap+charr, (767.5, 831.5)),
78 (-hw+edge+6*array+6*biggap+charr, -hw+edge+6*array+6*biggap+charr+gap,
False),
79 (-hw+edge+6*array+6*biggap+charr+gap, -hw+edge+6*array+6*biggap+array, (831.5, 895.5)),
80 (-hw+edge+7*array+6*biggap, -hw+edge+7*array+7*biggap,
False),
81 (-hw+edge+7*array+7*biggap, -hw+edge+7*array+7*biggap+charr, (895.5, 959.5)),
82 (-hw+edge+7*array+7*biggap+charr, -hw+edge+7*array+7*biggap+charr+gap,
False),
83 (-hw+edge+7*array+7*biggap+charr+gap, -hw+edge+7*array+7*biggap+array, (959.5, 1023.5)),
85 (-hw+edge+8*array+7*biggap, -hw+edge+8*array+8*biggap,
False),
86 (-hw+edge+8*array+8*biggap, -hw+edge+8*array+8*biggap+charr, (1023.5, 1087.5)),
87 (-hw+edge+8*array+8*biggap+charr, -hw+edge+8*array+8*biggap+charr+gap,
False),
88 (-hw+edge+8*array+8*biggap+charr+gap, -hw+edge+8*array+8*biggap+array, (1087.5, 1151.5)),
89 (-hw+edge+9*array+8*biggap, -hw+edge+9*array+9*biggap,
False),
90 (-hw+edge+9*array+9*biggap, -hw+edge+9*array+9*biggap+charr, (1151.5, 1215.5)),
91 (-hw+edge+9*array+9*biggap+charr, -hw+edge+9*array+9*biggap+charr+gap,
False),
92 (-hw+edge+9*array+9*biggap+charr+gap, -hw+edge+9*array+9*biggap+array, (1215.5, 1279.5)),
93 (-hw+edge+10*array+9*biggap, -hw+edge+10*array+10*biggap,
False),
94 (-hw+edge+10*array+10*biggap, -hw+edge+10*array+10*biggap+charr, (1279.5, 1343.5)),
95 (-hw+edge+10*array+10*biggap+charr, -hw+edge+10*array+10*biggap+charr+gap,
False),
96 (-hw+edge+10*array+10*biggap+charr+gap, -hw+edge+10*array+10*biggap+array, (1343.5, 1407.5)),
97 (-hw+edge+11*array+10*biggap, -hw+edge+11*array+11*biggap,
False),
98 (-hw+edge+11*array+11*biggap, -hw+edge+11*array+11*biggap+charr, (1407.5, 1471.5)),
99 (-hw+edge+11*array+11*biggap+charr, -hw+edge+11*array+11*biggap+charr+gap,
False),
100 (-hw+edge+11*array+11*biggap+charr+gap, -hw+edge+11*array+11*biggap+array, (1471.5, 1535.5))
138firstmatgap_startch = 521.5
139firstmatgap_endch = 523.5
140secmatgap_startch = 1046.5
141secmatgap_endch = 1048.5
147cluster_width_max = 10
153ly_loss_params = 20.78, -0.26, 7.89, -3.06
154ly_loss_sigma_params = 6.8482, -0.5757, 3.4458
155cluster_width_mean_params = 0.6661, -0.006955, 2.163
156cluster_width_sigma_params = 1.103, -0.0005751
160random_width_percent = 0.01
168energy_range = 0.18, 0.477
171CDF_integral = 0.0185640424
175ly_linear_params = 332.882, -58.7085
178k_cdfs_corr = 0.993076766938
182sigma_in_percent = 0.01
185sigma_from_width = 1 / 4.
195 lambda ly, *p: exp(p[0] + p[1]*sqrt(ly))
198 (0.0108183, -0.179752, -19.2, 0.00772965),
199 lambda ly, *p: p[0]*(1 - exp(p[1]*(ly+p[2]))) / (1 + exp(p[1]*(ly+p[2]))) + p[3]
206ly_CDF_landau_params = {
208 (0.001038, -0.000378, 3.53e-05),
209 lambda ly, *p: p[0] + p[1]*ly + p[2]*ly**2
212 (-0.001986, -0.0003014, 7.031e-05, -2.067e-06, 1.892e-08),
213 lambda ly, *p: p[0] + p[1]*ly + p[2]*ly**2 + p[3]*ly**3 + p[4]*ly**4
216 (-0.007149, 0.001056, -1.779e-05, 1.41e-07, -4.29e-10),
217 lambda ly, *p: p[0] + p[1]*ly + p[2]*ly**2 + p[3]*ly**3 + p[4]*ly**4
224 (89., 4.152, 0.0001574, -1.326e+04, 4.3),
225 lambda cdf, *p: p[0] * sqrt(p[1]*(cdf+p[2])) * (1-exp(p[3]*cdf)) + p[4]
228 (158, 1.035, 0.24, 217),
229 lambda cdf, *p: p[0] * log(sqrt(cdf)+p[1]) + p[2]*exp(p[3]*cdf)
231 (0.012, 0.018561405): (
232 (9.36, 335.984, -18100, -400, 15),
233 lambda cdf, *p: p[0] * log((p[1]-p[2]*cdf)/(p[1]+p[2]*cdf)) + p[3]*cdf + p[4]
235 (0.018561405, 0.0185640424): (
236 (9.36, 335.984, -18100, -400, 15),
237 lambda cdf, *p: (p[0] * log((p[1]-p[2]*0.018561405)/(p[1]+p[2]*0.018561405))
238 + p[3]*0.018561405 + p[4])
248snd_geo = ConfigRegistry.loadpy(
"$SNDSW_ROOT/geometry/shipLHC_geom_config.py")
250scifimat_width = snd_geo.Scifi.scifimat_width
251scifimat_length = snd_geo.Scifi.scifimat_length
252scifimat_gap = snd_geo.Scifi.scifimat_gap
254scifi_x = snd_geo.Scifi.xdim
255scifi_y = snd_geo.Scifi.ydim
256scifimat_z = snd_geo.Scifi.scifimat_z
258nmats = snd_geo.Scifi.nmats
259n_planes = snd_geo.Scifi.nscifi
260n_sipms = snd_geo.Scifi.nsipms
276def cm_to_channel(locpos, sipm_map=sipm_map, gaps_map=gaps_map, pitch=pitch, charr=charr,
277 reverse=False, ch_max_num=ch_max_num, n_solid_sipms=n_solid_sipms):
280 It converts a particle position (an event) measured in cm to a position measured
281 in channels. The SiPM map is used. The position is in the scifi module frame.
285 for left, right, ch_range
in sipm_map:
286 if left <= locpos <= right:
287 if not ch_range
is False:
288 ch_start = ch_range[0]
289 return int(round(ch_max_num - ((locpos-left) / pitch + ch_start)))
291 mapindex = int(round(n_solid_sipms - (locpos + hw) / charr))
292 return gaps_map.get(mapindex,
False)
294 elif reverse
is False:
295 for left, right, ch_range
in sipm_map:
296 if left <= locpos <= right:
297 if not ch_range
is False:
298 ch_start = ch_range[0]
299 return int(round(locpos-left) / pitch + ch_start)
301 mapindex2 = int(round(locpos + hw) / charr)
302 return gaps_map.get( mapindex2,
False)
304def channel_to_cm(channelpos, sipm_map=sipm_map, reverse=False, pitch=pitch):
307 It converts a particle position measured channels to a position measured
308 in cm. The SiPM map is used. The position is in the scifi module frame.
312 for left, _, ch_range
in sipm_map:
313 if not ch_range
is False:
314 ch_start, ch_end = ch_range
315 if ch_start <= channelpos <= ch_end:
316 return -(left + (channelpos - ch_start) * pitch)
318 for left, _, ch_range
in sipm_map:
319 if not ch_range
is False:
320 ch_start, ch_end = ch_range
321 if ch_start <= channelpos <= ch_end:
322 return left + (channelpos - ch_start) * pitch
327 It returns a type of a scifi module.
328 1 - vertical scifi assembly
329 0 - horizontal scifi assembly
331 return int((DetID % 100) // 10)
337 It returns an id (number) of a scifi module. In current version one plane has 3 vertical
338 and 3 horizontal scifi assemblies.
341 return int(DetID % 100 % 10)
346 It returns a length of a scifi mat. Identical for vertical and horizontal.
349 return scifimat_length
354 It returns a number of scifi mats in a plane. Identical for vertical and horizontal fiber planes.
362 It returns a number of SiPMs per plane. Identical for vertical and horizontal fiber planes.
370 It returns an id of a plane. In current version the detector has 5 TT stations (numbering starts
373 return int(DetID // 100)
382 It returns the local coordinates in one scifi assembly frame from global coordinates.
385 global xshift, yshift
388 return globalpos - yshift
391 return globalpos - xshift
397 It returns the global coordinates from the local coordinates in one scifi assembly frame.
400 global xshift, yshift
403 return localpos + yshift
406 return localpos + xshift
413 It returns the light yield depending on the distance to SiPMs
416 A1, k1, A2, k2 = params
417 return A1 * exp(k1 * distance / 100.) + A2 * exp(k2 * distance / 100.)
422 It return the light yield losses in percent depending on the distance to SiPMs
431 It return a mean cluster width depending on the distance to SiPMs
435 return exp(A + B*distance) + C
440 It return a standard deviation of the mean cluster width depending on the distance to SiPMs
444 return A + B*distance
447 cluster_width_min=cluster_width_min, cluster_width_max=cluster_width_max):
450 It generates a cluster. The cluster have 'ly' photoelectrons.
451 The cluster width depends on the distance to SiPM
456 if random.random() <= percent:
457 return random.randint(cluster_width_min, cluster_width_max)
458 random_width = int(round(random.gauss(mean - 1, sigma))) + 1
461 while random_width < 1
or ly < random_width:
462 random_width = int(round(random.gauss(mean - 1, sigma))) + 1
469 This universal function substitutes the parameters to the function.
470 The parameters and the function are in the dictionary
473 for (left, right), (params, func)
in approx_data.items():
474 if left <= var <= right:
475 return func(var, *params)
478def edep_to_ly(energy, CDF_integral=CDF_integral, energy_range=energy_range,
479 ly_linear_params=ly_linear_params, k_cdfs_corr=k_cdfs_corr, sigma_in_percent=sigma_in_percent,
480 ly_CDF_params=ly_CDF_params, CDF_ly_params=CDF_ly_params, ly_CDF_landau_params=ly_CDF_landau_params):
483 It returns the light yield calculated from the energy deposit. The calculations are based
484 on the equality of the cumulative distribution functions (CDF) :
485 energy => CDF(energy) => CDF(light yield) => ly
487 The linear converting range 0.18 MeV < dE < 0.477 MeV corresponds 4.5 < LY < 104 ph.e.
489 If energy more then 0.477 MeV the light yield calculated randomly (uniformly in the range)
490 according to the distribution
492 Also a little randomness is added to the CDF value with a normal distribution and
493 standard deviation with 'sigma_in_percent' (in percent of the whole range 0 - max CDF)
496 e_min, e_max = energy_range
497 A, C = ly_linear_params
498 if e_min < energy < e_max:
499 ly_lin = A * energy + C
501 elif e_max <= energy:
502 cdf_raw = CDF_integral * np.random.uniform(0., 1.0)
505 cdf_random = random.gauss(cdf_raw, sigma_in_percent * CDF_integral)
508 while cdf_random < 0
or cdf_random > CDF_integral:
509 cdf_random = random.gauss(cdf_raw, sigma_in_percent * CDF_integral)
514 chpos_min=chpos_min, chpos_max=chpos_max):
517 It generates an event cluster with given weighted mean position in channels, width and
520 If right side of the cluster can be out of range, the maximum of the right side will be
523 At first an array [0, 0, 0, ... ] is generated which corresponds to the channels.
524 Next the cluster generated in the array.
525 Final array will be like [0, 0, ..., 1, 2, 5, 1, 0, ...],
526 [0, 17, 0, ...] or etc.
529 if int(ceil(wmp + 0.5 + cluster_width_max)) < chpos_max:
530 max_fired_ch = int(ceil(wmp + 0.5 + cluster_width_max))
532 max_fired_ch = int(ceil(chpos_max))
533 cluster = [0
for _
in range(max_fired_ch)]
538 elif wmp == chpos_max:
539 fired_channel = int(chpos_max)
541 fired_channel = int(wmp + 0.5)
542 cluster[fired_channel] += amplitude
544 sigma = width * sigma_from_width
545 for _
in range(amplitude):
546 fired_channel = int(round(random.gauss(mean, sigma)))
547 while not 0 <= fired_channel < len(cluster):
548 fired_channel = int(round(random.gauss(mean, sigma)))
549 cluster[fired_channel] += 1
556 It returns TRUE if cluster is realistic: it doesn't have a gap between numders, like
557 [..., 0, 1, 2, 0, 0, 5, 6, ...], and it doens't have the light yield less then width.
560 cluster_only_values = [(channel, value)
for channel, value
in enumerate(cluster)
if value > 0]
561 first_channel, _ = cluster_only_values[0]
562 last_channel, _ = cluster_only_values[-1]
563 if len(cluster_only_values) != width
or (last_channel - first_channel + 1) != width:
571 The final function for creating a signal cluster
579 if not ((chpos_min < wmp)
and (wmp< chpos_max)):
582 if not ly_min <= amplitude >= width:
return False
583 shifted_wmp = wmp + 0.5
595 A function that prevents clusters from propagating to the gaps between the mats
598 if (wmp + 0.5)<= firstmatgap_startch:
599 for chann
in cluster:
600 if cluster[(chann + 0.5) > firstmatgap_startch] != 0.: cluster[chann] = 0.
602 elif firstmatgap_endch <= (wmp + 0.5)<= secmatgap_startch:
603 for chann
in cluster:
604 if (cluster[(chann + 0.5) < firstmatgap_endch] != 0.)
or (cluster[(chann + 0.5) > secmatgap_startch] != 0.): cluster[chann] = 0.
606 elif (wmp + 0.5) >= secmatgap_endch:
607 for chann
in cluster:
608 if cluster[(chann + 0.5)< secmatgap_endch] != 0.: cluster[chann] = 0.
616 Calculate the weighted mean position of the cluster
624 for channel, value
in enumerate(cluster):
626 sumup += value * channel
629 wmp = sumup / sumdown - 0.5
638 def __init__(self, DetID, Xraw, Yraw, Edep, Xshift=0, Yshift=0):
681 first_coord = self.
Xraw
682 second_coord = self.
Yraw
685 if Print: print(
"Detector ID -> " + str(self.
DetID) +
" -> Vertical")
687 first_coord = self.
Yraw
688 second_coord = self.
Xraw
691 if Print: print(
"Detector ID -> " + str(self.
DetID) +
" -> Horizontal")
696 if Print: print(
"First coordinate = " + str(first_coord))
698 if Print: print(
"Local Pos = " + str(localpos))
705 channelpos =
cm_to_channel(localpos, reverse=is_sipm_reversed)
715 if Print: print(
"Cluster width: " + str(cluster_width))
716 if Print: print(
"Channel position: " + str(channelpos))
721 if check_wall!=
False:
728 recovery_localpos =
channel_to_cm(wmp_of_cluster, reverse=is_sipm_reversed)
735 if abs(self.
delta) > 1:
759 return None,
None,
None
Class of SciFi Tracker ######################################################.
SetSipmIsReversed(self, is_sipm_hor_reversed, is_sipm_vert_reversed)
SetLYRange(self, ly_min, ly_max)
SetSipmPosition(self, sipm_hor_pos, sipm_vert_pos)
__init__(self, DetID, Xraw, Yraw, Edep, Xshift=0, Yshift=0)
ClusterGen(self, Print=False)
ly_loss_mean(distance, params=ly_loss_params)
channel_to_cm(channelpos, sipm_map=sipm_map, reverse=False, pitch=pitch)
cm_to_channel(locpos, sipm_map=sipm_map, gaps_map=gaps_map, pitch=pitch, charr=charr, reverse=False, ch_max_num=ch_max_num, n_solid_sipms=n_solid_sipms)
FUNCTIONS FOR SCIFICLASS ##################################################.
is_realistic(cluster, width)
cluster_width_mean(distance, params=cluster_width_mean_params)
cluster_width_random(distance, ly, percent=random_width_percent, cluster_width_min=cluster_width_min, cluster_width_max=cluster_width_max)
cluster_generator(amplitude, width, wmp, cluster_width_max=cluster_width_max, chpos_min=chpos_min, chpos_max=chpos_max)
create_cluster(amplitude, width, wmp)
local_to_global(DetID, localpos)
edep_to_ly(energy, CDF_integral=CDF_integral, energy_range=energy_range, ly_linear_params=ly_linear_params, k_cdfs_corr=k_cdfs_corr, sigma_in_percent=sigma_in_percent, ly_CDF_params=ly_CDF_params, CDF_ly_params=CDF_ly_params, ly_CDF_landau_params=ly_CDF_landau_params)
weigthed_mean_pos(cluster)
global_to_local(DetID, globalpos)
These functions assume that the detector is aligned with the beam ######.
fcheck_wall(cluster, wmp)
cluster_width_sigma(distance, params=cluster_width_sigma_params)
approx_function(var, approx_data)