Skip to content
Snippets Groups Projects
Commit 23a9e2e0 authored by Greg Henning's avatar Greg Henning
Browse files

fix guess and return in adaptive fit

parent 3821a72a
Branches
No related merge requests found
......@@ -24,7 +24,7 @@ import matplotlib.pyplot as plt
from pyhisto.histogram import Histogram
from pyhisto.tools import GaussRealFit
def adaptive_fit(
histo_file: str,
......@@ -68,7 +68,11 @@ def adaptive_fit(
# if no parasites,
if parasites is None or len(parasites) == 0:
# fall back to realfit
the_fit = GaussRealFit(the_histo)
the_fit = GaussRealFit(
the_histo,
_xguess=xpeak,
_stdevguess=guessstdev,
)
if plot is not None:
the_fit.plot(plot)
return yaml.dump(
......@@ -87,14 +91,19 @@ def adaptive_fit(
variables_limits = []
variables_parameters = []
# .. Background is either linear or constant
bg_guess = sum(the_histo)/(1.+ len(the_histo))
# print(f"# {bg_guess=}")
if bg_type == "lin":
functions_to_fit.append(lambda x, B, SL: B + SL * (x - xpeak))
variables_guess.extend([the_histo[0].count, 0.0])
variables_guess.extend([
bg_guess,
0.0
])
variables_limits.extend([(0.0, np.inf), (-np.inf, np.inf)])
variables_parameters.append(slice(0, 2))
else:
functions_to_fit.append(lambda x, B: B)
variables_guess.extend([the_histo[0].count])
variables_guess.extend([bg_guess])
variables_limits.extend([(0.0, np.inf)])
variables_parameters.append(slice(0, 1))
# .. Main peak
......@@ -103,9 +112,19 @@ def adaptive_fit(
/ np.sqrt(2.0 * np.pi * R * R)
* np.exp(-((x - M) ** 2.0) / 2.0 / R / R)
)
sum_guess = max(
1,
(
sum(the_histo) - (
the_histo[0].count +
the_histo[-1].count
)* len(the_histo) / 2.0
)
)
# print(f"# {sum_guess=}")
variables_guess.extend(
[
max(1, sum(the_histo) - the_histo[0].count * len(the_histo)),
sum_guess,
xpeak,
guessstdev,
]
......@@ -129,7 +148,7 @@ def adaptive_fit(
)
variables_guess.extend(
[
(max(1, sum(the_histo) - the_histo[0].count * len(the_histo)) / 2.0),
sum_guess / 2.0,
extra_peak,
guessstdev,
]
......@@ -268,14 +287,31 @@ if __name__ == "__main__":
try:
print(adaptive_fit(**vars(the_arguments)))
except:
# trying to get a guess of the integral
the_histo = Histogram(
fromfile=the_arguments.histo_file
).slice(
the_arguments.xmin,
the_arguments.xmax
)
sum_guess = max(
1,
(
sum(the_histo) - (
the_histo[0].count +
the_histo[-1].count
) * len(the_histo) / 2.0
)
)
# done
print(
"""
integral: 0.0
integral_u: 0.0
mean: 0.0
mean_u: 0.0
stdev: 0.0
stdev_u: 0.0
f"""
integral: {sum_guess}
integral_u: {2.35 * (sum_guess ** 0.5)}
mean: {the_arguments.xpeak}
mean_u: 1.0e30
stdev: {the_arguments.guessstdev}
stdev_u: 1.0e30
error: True
"""
)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment