Coverage for peeemtee/pmt_resp_func.py : 97%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# !/usr/bin/env python
3import numpy as np
4from scipy.stats.distributions import poisson
5from iminuit import Minuit
6from .tools import gaussian
9def fit_gaussian(x, y, print_level=1, calculate_hesse=False):
10 """
11 Fit a gaussian to data using iminuit migrad
13 Parameters
14 ----------
15 x: np.array
16 x values of the data
17 y: np.array
18 y values of the data
19 print_level: int, default: 1
20 0: quiet, 1: print fit details
22 Returns
23 -------
24 dict, dict:
25 optimal parameters {"mean": mean, "sigma": sigma, "A": A}
26 covariance matrix {("mean","mean"): cov_mean, ("mean", "sigma"): ...}
27 """
29 def make_quality_function(x, y):
30 def quality_function(mean, sigma, A):
31 return np.sum(((gaussian(x, mean, sigma, A) - y)) ** 2)
33 return quality_function
35 mean_start = x[y.argmax()]
36 above_half_max = x[y >= y.max() / 2]
37 if len(above_half_max) == 1:
38 sigma_start = np.abs(x[1] - x[0])
39 else:
40 sigma_start = (above_half_max.max() - above_half_max.min()) / 2.355
41 A_start = y.max() * np.sqrt(2 * np.pi) * sigma_start
42 qfunc = make_quality_function(x, y)
44 kwargs = {"mean": mean_start, "sigma": sigma_start, "A": A_start}
46 m = Minuit(
47 qfunc,
48 pedantic=False,
49 print_level=print_level,
50 **kwargs,
51 )
52 m.migrad()
53 if calculate_hesse:
54 m.hesse()
55 return m.values, m.covariance
56 return m.values, None
59class ChargeHistFitter(object):
60 """
61 Class that provides simple gaussian fit methods and
62 fit of pmt response function to pmt charge histogram.
64 Usage:
65 fitter = ChargeHistFitter()
66 # if fixed ped and spe is required for higher nphe measurements
67 fitter.fix_ped_spe(ped_mean, ped_sigma, spe_charge, spe_sigma)
68 fitter.pre_fit(x, y)
69 fitter.fit_pmt_resp_func(x, y)
71 Plotting the results:
72 plt.plot(x, fitter.pmt_resp_func(x, **fitter.popt_prf))
74 """
76 def __init__(self):
77 self.fixed_ped_spe = False
79 def _calculate_single_gaussian_shape(self, n):
80 popt = self.popt_prf
81 A = poisson.pmf(int(n), popt["nphe"]) * popt["entries"]
82 sigma = np.sqrt(n * popt["spe_sigma"] ** 2 + popt["ped_sigma"] ** 2)
83 mean = n * popt["spe_charge"] + popt["ped_mean"]
84 return mean, sigma, A
86 def single_gaussian_values(self, x, n):
87 """
88 Calculates single gaussian y-values corresponding to the input x-values
90 Parameters
91 ----------
92 x: np.array
93 1D array of x-values the gaussian y-values should be calculated for
94 n: int
95 multiplicity of the gaussian the y-values should be calculated for
97 Returns
98 -------
99 y: np.array
100 y-values of the gaussian
102 """
103 mean, sigma, A = self._calculate_single_gaussian_shape(n)
104 return gaussian(x, mean, sigma, A)
106 def pmt_resp_func(
107 self, x, nphe, ped_mean, ped_sigma, spe_charge, spe_sigma, entries
108 ):
109 func = 0.0
110 for i in range(self.n_gaussians):
111 pois = poisson.pmf(int(i), nphe)
112 sigma = np.sqrt(i * spe_sigma ** 2 + ped_sigma ** 2)
113 arg = (x - (i * spe_charge + ped_mean)) / sigma
114 func += pois / sigma * np.exp(-0.5 * arg ** 2)
115 func = entries * func / np.sqrt(2 * np.pi)
116 return func
118 def pmt_resp_func_uap(
119 self,
120 x,
121 nphe,
122 ped_mean,
123 ped_sigma,
124 spe_charge,
125 spe_sigma,
126 entries,
127 uap_mean,
128 uap_sigma,
129 uap_A,
130 ):
131 func = 0.0
132 for i in range(self.n_gaussians):
133 pois = poisson.pmf(int(i), nphe)
134 sigma = np.sqrt(i * spe_sigma ** 2 + ped_sigma ** 2)
135 arg = (x - (i * spe_charge + ped_mean)) / sigma
136 func += pois / sigma * np.exp(-0.5 * arg ** 2)
137 func = entries * func / np.sqrt(2 * np.pi)
138 func += gaussian(x, uap_mean, uap_sigma, uap_A)
139 return func
141 def fix_ped_spe(self, ped_mean, ped_sigma, spe_charge, spe_sigma):
142 """
143 Fixes ped and spe in fit_pmt_resp_func and sets fixed parameters
145 Parameters
146 ----------
147 ped_mean: float
148 mean of gaussian fit of pedestal
149 ped_sigma: float
150 sigma of gaussian fit of pedestal
151 spe_charge: float
152 charge of spe (spe_mean - spe_charge)
153 spe_sigma: float
154 sigma of gaussian fit of spe peak
156 """
157 self.fixed_ped_spe = True
158 self.popt_ped = {"mean": ped_mean, "sigma": ped_sigma}
159 self.popt_spe = {"sigma": spe_sigma}
160 self.spe_charge = spe_charge
162 def pre_fit(
163 self,
164 x,
165 y,
166 valley=None,
167 spe_upper_bound=None,
168 n_sigma=3,
169 print_level=1,
170 calculate_hesse=False,
171 ):
172 """
173 Performs single gaussian fits to pedestal and single p.e. peaks
175 Parameters
176 ----------
177 x: np.array
178 bin centres of the charge histogram
179 y: np.array
180 bin counts of the charge histogram
181 valley: float (default=None)
182 user set x position to split hisogram in pedestal and spe
183 spe_upper_bound: float (default=None)
184 user set upper bound for gaussian fit of spe peak
185 n_sigma: float
186 spe fit range starts at:
187 mean of pedestal + n_sigma * sigma of pedestal
188 print_level: int, default: 1
189 0: quiet, 1: print fit details
190 """
192 if self.fixed_ped_spe:
193 popt, pcov = fit_gaussian(
194 x,
195 y,
196 print_level=print_level,
197 calculate_hesse=calculate_hesse,
198 )
199 self.popt_gauss = popt
200 self.nphe = popt["mean"] / self.spe_charge
201 if self.nphe < 0:
202 self.nphe = 0
203 self.entries = np.max(y)
204 self.n_gaussians = int(self.nphe * 2)
205 if print_level > 0:
206 print(self.entries)
208 else:
210 if valley is None:
211 x_ped, y_ped = x, y
212 else:
213 cond = x < valley
214 x_ped, y_ped = x[cond], y[cond]
216 popt_ped, pcov_ped = fit_gaussian(
217 x_ped,
218 y_ped,
219 print_level=print_level,
220 calculate_hesse=calculate_hesse,
221 )
223 if valley is None:
224 if spe_upper_bound is None:
225 cond = x > (popt_ped["mean"] + n_sigma * popt_ped["sigma"])
226 else:
227 cond = (
228 x > (popt_ped["mean"] + n_sigma * popt_ped["sigma"])
229 ) & (x < spe_upper_bound)
230 else:
231 if spe_upper_bound is None:
232 cond = x > valley
233 else:
234 cond = (x > valley) & (x < spe_upper_bound)
235 x_spe, y_spe = x[cond], y[cond]
237 popt_spe, pcov_spe = fit_gaussian(
238 x_spe,
239 y_spe,
240 print_level=print_level,
241 calculate_hesse=calculate_hesse,
242 )
244 self.popt_ped = popt_ped
245 self.pcov_ped = pcov_ped
246 self.popt_spe = popt_spe
247 self.pcov_spe = pcov_spe
248 self.opt_ped_values = gaussian(x, **popt_ped)
249 self.opt_spe_values = gaussian(x, **popt_spe)
251 self.spe_charge = popt_spe["mean"] - popt_ped["mean"]
252 self.nphe = -np.log(popt_ped["A"] / (popt_ped["A"] + popt_spe["A"]))
253 if self.nphe < 0:
254 self.nphe = 0
255 self.n_gaussians = 10
257 def fit_pmt_resp_func(
258 self,
259 x,
260 y,
261 n_gaussians=None,
262 print_level=1,
263 mod=False,
264 strong_limits=True,
265 fixed_parameters=["ped_mean", "ped_sigma"],
266 ):
267 """
268 Performs fit of pmt response function to charge histogram
270 Parameters
271 ----------
272 x: np.array
273 bin centres of the charge histogram
274 y: np.array
275 bin counts of the charge histogram
276 n_gaussians: int
277 number of gaussians to be fitted
278 mod: string
279 if False: no modification
280 "uap": fits an additional underamplified pulse gaussian
281 print_level: int, default: 1
282 0: quiet, 1: print fit details
283 """
284 if n_gaussians:
285 self.n_gaussians = n_gaussians
286 if not mod:
287 func = self.pmt_resp_func
289 if mod == "uap":
290 func = self.pmt_resp_func_uap
292 self.used_fit_function = func
294 def make_quality_function(x, y, mod):
295 if not mod:
297 def quality_function(
298 nphe, ped_mean, ped_sigma, spe_charge, spe_sigma, entries
299 ):
300 model = func(
301 x,
302 nphe,
303 ped_mean,
304 ped_sigma,
305 spe_charge,
306 spe_sigma,
307 entries,
308 )
309 mask = y != 0
310 return np.sum(((model[mask] - y[mask])) ** 2 / y[mask])
312 if mod == "uap":
314 def quality_function(
315 nphe,
316 spe_charge,
317 ped_mean,
318 ped_sigma,
319 spe_sigma,
320 entries,
321 uap_mean,
322 uap_sigma,
323 uap_A,
324 ):
325 model = func(
326 x,
327 nphe,
328 ped_mean,
329 ped_sigma,
330 spe_charge,
331 spe_sigma,
332 entries,
333 uap_mean,
334 uap_sigma,
335 uap_A,
336 )
337 mask = y != 0
338 return np.sum(((model[mask] - y[mask])) ** 2 / y[mask])
340 return quality_function
342 qfunc = make_quality_function(x, y, mod=mod)
343 self.qfunc = qfunc
344 if self.fixed_ped_spe:
345 entries_start = self.entries
346 else:
347 entries_start = self.popt_ped["A"] + self.popt_spe["A"]
349 kwargs = {
350 "nphe": self.nphe,
351 "spe_charge": self.spe_charge,
352 "spe_sigma": self.popt_spe["sigma"],
353 "entries": entries_start,
354 "ped_mean": self.popt_ped["mean"],
355 "ped_sigma": self.popt_ped["sigma"],
356 }
357 if strong_limits:
358 kwargs["limit_nphe"] = (0, self.nphe * 2)
359 kwargs["limit_spe_charge"] = (
360 self.spe_charge / 2,
361 self.spe_charge * 2,
362 )
363 # kwargs["limit_spe_sigma"] = (0, self.popt_spe["sigma"] * 2)
364 kwargs["limit_spe_sigma"] = (0, self.spe_charge * 0.8)
365 kwargs["limit_entries"] = (0, entries_start * 2)
367 for parameter in fixed_parameters:
368 kwargs[f"fix_{parameter}"] = True
369 if self.fixed_ped_spe:
370 kwargs["fix_spe_charge"] = True
371 kwargs["fix_spe_sigma"] = True
372 kwargs["fix_ped_mean"] = True
373 kwargs["fix_ped_sigma"] = True
375 if mod == "uap":
376 kwargs["uap_mean"] = self.popt_spe["mean"] / 5
377 kwargs["uap_sigma"] = self.popt_spe["sigma"] / 5
378 kwargs["uap_A"] = entries_start / 200
379 kwargs["limit_uap_mean"] = (
380 self.popt_ped["mean"],
381 self.popt_spe["mean"] / 3,
382 )
383 kwargs["limit_uap_A"] = (0, entries_start / 100)
384 kwargs["limit_uap_sigma"] = (0, self.popt_spe["sigma"] / 3)
386 self.m = Minuit(
387 qfunc,
388 pedantic=False,
389 print_level=print_level,
390 throw_nan=True,
391 **kwargs,
392 )
393 try:
394 self.m.migrad()
395 except RuntimeError:
396 self.success = False
397 else:
398 self.success = True
399 # self.m.hesse()
400 self.popt_prf = self.m.values
401 self.opt_prf_values = func(x, **self.m.values)
402 self.pcov_prf = self.m.covariance