SND@LHC Software
Loading...
Searching...
No Matches
ScifiCluster.py
Go to the documentation of this file.
1#!/usr/bin/python
2
3import ROOT#
4import numpy as np
5import math
6from math import exp, sqrt, log, ceil
7import random
8import matplotlib.pyplot as plt
9
10from ROOT.TMath import Landau, Gaus, Poisson
11from ShipGeoConfig import ConfigRegistry
12
13
16
17# Position of SiPMs in the detector (left/right and up/down)
18sipm_hor_pos = +1 # at X = +len_vert / 2
19sipm_vert_pos = -1 # at Y = -len_hor / 2
20
21# If reversed is true the channels counted in reverse way
22# (0->1535 => 1535<-0). It connected with how the modules are placed
23is_sipm_hor_reversed = True
24is_sipm_vert_reversed = True
25
26# SiPM geometry
27# hw - half width
28# pitch - one channel width
29# charr - array of 64 channels without gaps
30# edge - from left and right sides of the SiPM
31# gap - between two charr
32# array - 128 channels with the gap
33# biggap - gap between two arrays
34# ch_max_num - max channel in the SiPM (0-1535)
35# n_solid_sipms - number of "charr" arrays
36# Total SiPM size is 32.54 mm
37
38hw = 19.528 # Position of the first SiPM
39pitch = 0.025
40charr = 1.6
41edge = 0.016
42gap = 0.022 # Die gap
43array = gap + 2.*charr
44biggap = 0.038 # 0.006cm is the space between SiPMs, plus 2*edge (in CAD SiPMs are separated by 32.60 mm)
45ch_max_num = 128*12 - 1 # numbering from 0 to 1535
46n_solid_sipms = 12
47
48# SiPM geometry. If a point is in the solid part then the position will be counted.
49# Else false value returns.
50sipm_map = (
51 (-hw, -hw+edge, False),
52 (-hw+edge, -hw+edge+charr, (-0.5, 63.5)), # 1.1
53 (-hw+edge+charr, -hw+edge+charr+gap, False),
54 (-hw+edge+charr+gap, -hw+edge+array, (63.5, 127.5)), # 1.2
55 (-hw+edge+array, -hw+edge+array+biggap, False),
56 (-hw+edge+array+biggap, -hw+edge+array+biggap+charr, (127.5, 191.5)), # 2.1
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)), # 2.2
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)), # 3.1
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)), # 3.2
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)), # 4.1
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)), # 4.2
67
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)), # 5.1
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)), # 5.2
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)), # 6.1
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)), # 6.2
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)), # 7.1
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)), # 7.2
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)), # 8.1
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)), # 8.2
84
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)), # 9.1
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)), # 9.2
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)), # 10.1
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)), # 10.2
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)), # 11.1
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)), # 11.2
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)), #12.1
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)) #12.2
101
102)
103
104# If a point is in the dead space then the following channel value will be selected
105gaps_map = {
106 0.: -0.5,
107 1.: 63.5,
108 2.: 127.5,
109 3.: 191.5,
110 4.: 255.5,
111 5.: 319.5,
112 6.: 383.5,
113 7.: 447.5,
114
115 8.: 511.5,
116 9.: 575.5,
117 10.: 639.5,
118 11.: 703.5,
119 12.: 767.5,
120 13.: 831.5,
121 14.: 895.5,
122 15.: 959.5,
123
124 16.: 1023.5,
125 17.: 1087.5,
126 18.: 1151.5,
127 19.: 1215.5,
128 20.: 1279.5,
129 21.: 1343.5,
130 22.: 1407.5,
131 23.: 1471.5,
132 24.: 1535.5
133}
134
135# 1st mat gap is from -62.28mm -> -61.78mm ; 521.5 -> 523.5 channel
136# 2nd mat gap is from 71.22mm -> 71.72 mm ; 1046.5 -> 1048.5 channel
137# the gaps start and end inside these channels
138firstmatgap_startch = 521.5
139firstmatgap_endch = 523.5
140secmatgap_startch = 1046.5
141secmatgap_endch = 1048.5
142
143# Possible values of channel position and cluster width
144chpos_min = -0.5
145chpos_max = 1535.5
146cluster_width_min = 1
147cluster_width_max = 10
148
149#------------------ former values from SHiP's Target Tracker--------------------------------
150
151# Light yield attenuation A*exp(B*x) + C*exp(D*x)
152# Sigma depends on the distance from a SiPM and also the following parameters
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
157
158# The random width in the range is manually set as 1% of events
159# It covers cases where the width is 5-10
160random_width_percent = 0.01
161
162# Possible values of LY
163ly_min = 4.5
164ly_max = 104
165
166# Energy deposit converts in the light yield range linearly. Energy less then 0.18 MeV gives doesn't registered.
167# Energy more then 0.477 MeV converts randomly with approximation distribution.
168energy_range = 0.18, 0.477
169
170# The maximum value of the CDF integral
171CDF_integral = 0.0185640424
172
173# Parameters that are used for linear conversion of deposit energy to the range
174# equal to the light yield range
175ly_linear_params = 332.882, -58.7085
176
177# The coefficient between the maximums of CDF(E) and CDF(LY), a little differs
178k_cdfs_corr = 0.993076766938
179
180# The addition of a little randomness in the converting from the energy distibution to
181# the light yield distibution
182sigma_in_percent = 0.01
183
184# 4 sigma range includes 95% of events. It needs for creating a cluster
185sigma_from_width = 1 / 4.
186
187# The following parameters will be used in the converting from the energy distibution to the light yield distibution.
188# Dictionaries format is following - (x_min, x_max) : (parameters), approximating function.
189
190# ly_CDF_params - approximation of CDF of the light yield distribution.
191# If is doesn't used already.
192ly_CDF_params = {
193 (4.5, 13): (
194 (-13.2, 1.976),
195 lambda ly, *p: exp(p[0] + p[1]*sqrt(ly))
196 ),
197 (13, 104): (
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]
200 )
201}
202
203
204# Get a CDF value from LY (actually from an energy deposit which preliminarily linearly
205# converted to the range of the light yield (4.5 - 104 ph.e.)
206ly_CDF_landau_params = {
207 (4.5, 15): (
208 (0.001038, -0.000378, 3.53e-05),
209 lambda ly, *p: p[0] + p[1]*ly + p[2]*ly**2
210 ),
211 (15, 40): (
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
214 ),
215 (40, 104): (
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
218 )
219}
220
221# Get a light yield value from a CDF values
222CDF_ly_params = {
223 (0., 0.0006): (
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]
226 ),
227 (0.0006, 0.012): (
228 (158, 1.035, 0.24, 217),
229 lambda cdf, *p: p[0] * log(sqrt(cdf)+p[1]) + p[2]*exp(p[3]*cdf)
230 ),
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]
234 ),
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])
239 )
240}
241
242#------------------------------------------------------------------------------------
243
244
247
248snd_geo = ConfigRegistry.loadpy("$SNDSW_ROOT/geometry/shipLHC_geom_config.py")
249
250scifimat_width = snd_geo.Scifi.scifimat_width # 13.3 cm
251scifimat_length = snd_geo.Scifi.scifimat_length # 40 cm
252scifimat_gap = snd_geo.Scifi.scifimat_gap # 0.05 cm
253
254scifi_x = snd_geo.Scifi.xdim # 40 cm
255scifi_y = snd_geo.Scifi.ydim # 40 cm
256scifimat_z = snd_geo.Scifi.scifimat_z # 0.135 cm
257
258nmats = snd_geo.Scifi.nmats # 3
259n_planes = snd_geo.Scifi.nscifi # 5
260n_sipms = snd_geo.Scifi.nsipms # 12 per plane
261
262
267
268#from getGeoInformation.py
269xshift = -27.8
270yshift = +35.3
271
272
275
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):
278
279 """
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.
282 """
283
284 if reverse is True:
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)))
290 else:
291 mapindex = int(round(n_solid_sipms - (locpos + hw) / charr))
292 return gaps_map.get(mapindex, False)
293
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)
300 else:
301 mapindex2 = int(round(locpos + hw) / charr)
302 return gaps_map.get( mapindex2, False)
303
304def channel_to_cm(channelpos, sipm_map=sipm_map, reverse=False, pitch=pitch):
305
306 """
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.
309 """
310
311 if reverse is True:
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)
317 if reverse is False:
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
323
324def GetMatType(DetID):
325
326 """
327 It returns a type of a scifi module.
328 1 - vertical scifi assembly
329 0 - horizontal scifi assembly
330 """
331 return int((DetID % 100) // 10)
332
333
334def GetMatNum(DetID):
335
336 """
337 It returns an id (number) of a scifi module. In current version one plane has 3 vertical
338 and 3 horizontal scifi assemblies.
339 """
340
341 return int(DetID % 100 % 10)
342
344
345 """
346 It returns a length of a scifi mat. Identical for vertical and horizontal.
347 """
348
349 return scifimat_length
350
352
353 """
354 It returns a number of scifi mats in a plane. Identical for vertical and horizontal fiber planes.
355 """
356
357 return n_planes
358
360
361 """
362 It returns a number of SiPMs per plane. Identical for vertical and horizontal fiber planes.
363 """
364
365 return n_sipms
366
367def GetStationNum(DetID):
368
369 """
370 It returns an id of a plane. In current version the detector has 5 TT stations (numbering starts
371 at 1.
372 """
373 return int(DetID // 100)
374
375
378
379def global_to_local(DetID, globalpos):
380
381 """
382 It returns the local coordinates in one scifi assembly frame from global coordinates.
383 """
384
385 global xshift, yshift
386
387 if GetMatType(DetID) == 0:
388 return globalpos - yshift
389
390 if GetMatType(DetID) == 1:
391 return globalpos - xshift
392
393
394def local_to_global(DetID, localpos):
395
396 """
397 It returns the global coordinates from the local coordinates in one scifi assembly frame.
398 """
399
400 global xshift, yshift
401
402 if GetMatType(DetID) == 0:
403 return localpos + yshift
404
405 if GetMatType(DetID) == 1:
406 return localpos + xshift
407
408
409
410def ly_loss_mean(distance, params=ly_loss_params):
411
412 """
413 It returns the light yield depending on the distance to SiPMs
414 """
415
416 A1, k1, A2, k2 = params
417 return A1 * exp(k1 * distance / 100.) + A2 * exp(k2 * distance / 100.)
418
419def ly_attenuation(distance):
420
421 """
422 It return the light yield losses in percent depending on the distance to SiPMs
423 """
424
425 res_ly_loss = ly_loss_mean(distance) / ly_loss_mean(0.)
426 return res_ly_loss
427
428def cluster_width_mean(distance, params=cluster_width_mean_params):
429
430 """
431 It return a mean cluster width depending on the distance to SiPMs
432 """
433
434 A, B, C = params
435 return exp(A + B*distance) + C
436
437def cluster_width_sigma(distance, params=cluster_width_sigma_params):
438
439 """
440 It return a standard deviation of the mean cluster width depending on the distance to SiPMs
441 """
442
443 A, B = params
444 return A + B*distance
445
446def cluster_width_random(distance, ly, percent=random_width_percent,
447 cluster_width_min=cluster_width_min, cluster_width_max=cluster_width_max):
448
449 """
450 It generates a cluster. The cluster have 'ly' photoelectrons.
451 The cluster width depends on the distance to SiPM
452 """
453
454 mean = cluster_width_mean(distance)
455 sigma = cluster_width_sigma(distance)
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
459
460 # Generate again if the width < minimal width and the light yield < the width
461 while random_width < 1 or ly < random_width:
462 random_width = int(round(random.gauss(mean - 1, sigma))) + 1
463
464 return random_width
465
466def approx_function(var, approx_data):
467
468 """
469 This universal function substitutes the parameters to the function.
470 The parameters and the function are in the dictionary
471 """
472
473 for (left, right), (params, func) in approx_data.items():
474 if left <= var <= right:
475 return func(var, *params)
476 return False
477
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):
481
482 """
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
486
487 The linear converting range 0.18 MeV < dE < 0.477 MeV corresponds 4.5 < LY < 104 ph.e.
488
489 If energy more then 0.477 MeV the light yield calculated randomly (uniformly in the range)
490 according to the distribution
491
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)
494 """
495
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
500 cdf_raw = approx_function(ly_lin, ly_CDF_landau_params) * k_cdfs_corr
501 elif e_max <= energy:
502 cdf_raw = CDF_integral * np.random.uniform(0., 1.0)
503 else:
504 return 0.
505 cdf_random = random.gauss(cdf_raw, sigma_in_percent * CDF_integral)
506
507 # Generate again while the CDF value is not in the range
508 while cdf_random < 0 or cdf_random > CDF_integral:
509 cdf_random = random.gauss(cdf_raw, sigma_in_percent * CDF_integral)
510
511 return approx_function(cdf_random, CDF_ly_params)
512
513def cluster_generator(amplitude, width, wmp, cluster_width_max=cluster_width_max,
514 chpos_min=chpos_min, chpos_max=chpos_max):
515
516 """
517 It generates an event cluster with given weighted mean position in channels, width and
518 amplitude.
519
520 If right side of the cluster can be out of range, the maximum of the right side will be
521 right channel.
522
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.
527 """
528
529 if int(ceil(wmp + 0.5 + cluster_width_max)) < chpos_max:
530 max_fired_ch = int(ceil(wmp + 0.5 + cluster_width_max))
531 else:
532 max_fired_ch = int(ceil(chpos_max))
533 cluster = [0 for _ in range(max_fired_ch)]
534 mean = wmp
535 if width == 1:
536 if wmp == chpos_min:
537 fired_channel = 0
538 elif wmp == chpos_max:
539 fired_channel = int(chpos_max)
540 else:
541 fired_channel = int(wmp + 0.5)
542 cluster[fired_channel] += amplitude
543 else:
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
550
551 return cluster
552
553def is_realistic(cluster, width):
554
555 """
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.
558 """
559
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:
564 return False
565 else:
566 return True
567
568def create_cluster(amplitude, width, wmp):
569
570 """
571 The final function for creating a signal cluster
572
573 """
574 if wmp is False:
575 return False
576 if wmp is None:
577 return False
578
579 if not ((chpos_min < wmp) and (wmp< chpos_max)):
580 return False
581
582 if not ly_min <= amplitude >= width: return False
583 shifted_wmp = wmp + 0.5 # For right counting
584 cluster = cluster_generator(amplitude, width, shifted_wmp)
585
586 # Generate again if it doesn't look like real cluster
587 while is_realistic(cluster, width) is False:
588 cluster = cluster_generator(amplitude, width, shifted_wmp)
589
590 return cluster
591
592def fcheck_wall(cluster, wmp):
593
594 """
595 A function that prevents clusters from propagating to the gaps between the mats
596 """
597
598 if (wmp + 0.5)<= firstmatgap_startch: #1st mat
599 for chann in cluster:
600 if cluster[(chann + 0.5) > firstmatgap_startch] != 0.: cluster[chann] = 0.
601 return cluster
602 elif firstmatgap_endch <= (wmp + 0.5)<= secmatgap_startch: #2nd mat
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.
605 return cluster
606 elif (wmp + 0.5) >= secmatgap_endch: #3rd mat
607 for chann in cluster:
608 if cluster[(chann + 0.5)< secmatgap_endch] != 0.: cluster[chann] = 0.
609 return cluster
610 else:
611 return False #"False" is returned if wmp lies in the gaps
612
613def weigthed_mean_pos(cluster):
614
615 """
616 Calculate the weighted mean position of the cluster
617 """
618
619 if cluster is False:
620 return False
621 sumup = 0.
622 sumdown = 0.
623 wmp = 0.
624 for channel, value in enumerate(cluster):
625 if value > 0:
626 sumup += value * channel
627 sumdown += value
628 if sumdown != 0:
629 wmp = sumup / sumdown - 0.5
630 return wmp
631
632
635
637 # Constructor
638 def __init__(self, DetID, Xraw, Yraw, Edep, Xshift=0, Yshift=0):
639 self.DetID = DetID
640 self.Xraw = Xraw
641 self.Yraw = Yraw
642 self.Edep = Edep * 1000 # to MeV
643 self.ly = 0.
644 self.cluster = []
645 self.station = 0
646 self.matnum = 0
647 self.Xrec = None
648 self.Yrec = None
649 self.is_created = False
650 # Only to check the edge events
652 self.delta = -13
653 # The default values should be set
654 self.sipm_hor_pos = +1 # +1 at X = +len_vert / 2
655 self.sipm_vert_pos = +1 #-1 at Y = -len_hor / 2
656 self.is_sipm_hor_reversed = False # Channels and axis are co-directional or not
658 self.ly_min = 4.5
659 self.ly_max = 104
660 self.Xshift = Xshift
661 self.Yshift = Yshift
662
663 def SetLYRange(self, ly_min, ly_max):
664 self.ly_min = ly_min
665 self.ly_max = ly_max
666
667 def SetSipmPosition(self, sipm_hor_pos, sipm_vert_pos):
668 self.sipm_hor_pos = sipm_hor_pos
669 self.sipm_vert_pos = sipm_vert_pos
670
671 def SetSipmIsReversed(self, is_sipm_hor_reversed, is_sipm_vert_reversed):
672 self.is_sipm_hor_reversed = is_sipm_hor_reversed
673 self.is_sipm_vert_reversed = is_sipm_vert_reversed
674
675 # Cluster generator from the raw data. It returns False if it fails.
676 def ClusterGen(self, Print = False):
677
678 if GetMatType(self.DetID) is False:
679 return False
680 elif GetMatType(self.DetID) == 1: # vert
681 first_coord = self.Xraw
682 second_coord = self.Yraw
683 sipm_pos = self.sipm_vert_pos
684 is_sipm_reversed = self.is_sipm_hor_reversed
685 if Print: print("Detector ID -> " + str(self.DetID) + " -> Vertical")
686 elif GetMatType(self.DetID) == 0: # hor
687 first_coord = self.Yraw
688 second_coord = self.Xraw
689 sipm_pos = self.sipm_hor_pos
690 is_sipm_reversed = self.is_sipm_vert_reversed
691 if Print: print("Detector ID -> " + str(self.DetID) + " -> Horizontal")
692 self.matnum = GetMatNum(self.DetID)
693 nmats = GetMatQty()
694 matlen = GetMatLength()
695 self.station = GetStationNum(self.DetID)
696 if Print: print("First coordinate = " + str(first_coord))
697 localpos = global_to_local(self.DetID, first_coord)
698 if Print: print("Local Pos = " + str(localpos))
699
700 if sipm_pos == +1:
701 distance = GetMatLength()/2. - second_coord
702 elif sipm_pos == -1:
703 distance = GetMatLength()/2. + second_coord
704
705 channelpos = cm_to_channel(localpos, reverse=is_sipm_reversed)
706
707 self.ly = int(round(
708 edep_to_ly(self.Edep) * ly_attenuation(distance)
709 ))
710
711 if not self.ly_min < self.ly < self.ly_max:
712 return False
713
714 cluster_width = int(round(cluster_width_random(distance, ly=self.ly)))
715 if Print: print("Cluster width: " + str(cluster_width))
716 if Print: print("Channel position: " + str(channelpos))
717 cluster = create_cluster(self.ly, cluster_width, channelpos)
718
719 check_wall = fcheck_wall(cluster,channelpos)
720 #Cluster is returned if wmp is not in the 2 channels that are entirely in the mat gaps
721 if check_wall!=False:
722 self.is_created=True
723 cluster=check_wall
724 else:
725 self.is_created=False
726
727 wmp_of_cluster = weigthed_mean_pos(cluster)
728 recovery_localpos = channel_to_cm(wmp_of_cluster, reverse=is_sipm_reversed)
729
730 self.recovery_globalpos = local_to_global(self.DetID, recovery_localpos)
731 if Print: print("Position readout by the SiPMs => " + str(self.recovery_globalpos))
732 self.delta = first_coord - self.recovery_globalpos
733
734 # Some edge events may be reconstructed incorrectly
735 if abs(self.delta) > 1:
736 return False
737
738 if GetMatType(self.DetID) == 1: # vert
739 self.Xrec = self.recovery_globalpos
740 elif GetMatType(self.DetID) == 0: # hor
741 self.Yrec = self.recovery_globalpos
742
743 def GetType(self):
744
745 if GetMatType(self.DetID) == 0: #hor
746 return 0
747 else: #vert
748 return 1
749
750 def GetNSiPMs(self):
751
752 return GetSiPMs() #n_sipms
753
754 def GetXYZ(self):
755
756 if self.is_created:
757 return self.Xrec, self.Yrec, self.station
758 else:
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)
GetMatType(DetID)
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)
ly_attenuation(distance)
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)
GetStationNum(DetID)
cluster_width_sigma(distance, params=cluster_width_sigma_params)
approx_function(var, approx_data)