Is there a more efficient way (or obvious way using TF1 that I missed) to build TF1s from currying with C-style formula strings and information from a histogram (using PyROOT) then how I have done in this example (also attached TF1_from_functions.py (4.0 KB))? The inefficient part seems to be having to parse the formula content for C-style math functions and convert them to ROOT library functions (e.g., pow to TMath::Power) and parameter arrays ([0] to p[0]).
Is there a smarter or better way to do this?
import re
import ROOT
from ROOT import gROOT
from ROOT import TCanvas, TF1, TH1F
def build_full_model(model_function, hists):
"""
Currying function to build a compilable function for TF1's constructor
https://root.cern.ch/doc/master/classTF1.html#aa8905d28455ed7be02019f20b9cc5827
Args:
model_function (`str`): The C-style formula expression
hists (`list` of `TH1`s): The histograms to add into the model
Returns:
`TF1`: A model of the function and the histograms
`int`: The number of free parameters in the model
"""
# Search for number inside of [] (i.e., [0])
n_function_parameters = max(
[int(item) for item in re.findall(r'\[([0-9_]+)\]', model_function)])
for C_op, ROOT_op in zip(['exp', 'log', 'pow'],
['Exp', 'Log', 'Power']):
model_function = model_function.replace(C_op, 'ROOT.TMath.' + ROOT_op)
for idx in range(n_function_parameters + 1):
model_function = model_function.replace(
'[{}]'.format(idx), 'p[{}]'.format(idx))
model_function = model_function.replace(' ', '')
for operator in ['+', '-', '*']:
model_function = model_function.replace('{}x'.format(operator),
'{}x[0]'.format(operator))
model_function = model_function.replace('x/', 'x[0]/')
def func(x, p):
hist_value = [
hist.GetBinContent(hist.FindBin(x[0])) for hist in hists]
hist_terms = ' + '.join(
['p[{}] * {}'.format(n_function_parameters + idx + 1, hist_value[idx])
for idx in range(len(hists))])
full_model = ' + '.join([model_function, hist_terms])
return eval(full_model)
n_parameters = n_function_parameters + len(hists) + 1
return func, n_parameters
if __name__ == '__main__':
gROOT.SetBatch(ROOT.kTRUE)
# Set error ignore level above kInfo
ROOT.gErrorIgnoreLevel = ROOT.kWarning
x_range = [0, 10]
falling_exponential = TF1(
'falling_exponential', 'exp([0]-[1]*x)', x_range[0], x_range[1])
# Expect to have parameters live in a JSON config file
# so this dictionary is mocking a possible structure
parameters = {
'p0': {
'value': 5,
'name': 'constant'
},
'p1': {
'value': 0.3,
'name': 'decay rate'
},
}
falling_exponential.SetParameters(*[parameters['p{}'.format(param)]['value']
for param in range(len(parameters))])
falling_exponential.SetParNames(*[parameters['p{}'.format(param)]['name']
for param in range(len(parameters))])
falling_exponential.SetLineColor(ROOT.kBlack)
# falling_exponential.Print()
gaussian = TF1(
'gaussian', '[0]*TMath::Gaus(x,[1],[2],0)', x_range[0], x_range[1])
gaussian.SetParameters(7, 3, 2)
gaussian_hist = gaussian.GetHistogram()
gauss_2 = TF1('gauss_2', 'TMath::Gaus(x,5,1)', x_range[0], x_range[1])
rand_hist = TH1F('rand_hist', '',
int((x_range[1] - x_range[0]) / 0.2), x_range[0], x_range[1])
rand_hist.FillRandom('gauss_2', 10000)
rand_hist.Scale(200 / rand_hist.Integral())
rand_hist.SetLineColor(ROOT.kViolet)
model_func, n_model_params = build_full_model(
falling_exponential.GetFormula().GetTitle(), [gaussian_hist, rand_hist])
model = TF1('model', model_func, x_range[0], x_range[1], n_model_params)
model.SetParameters(
falling_exponential.GetParameter(0), falling_exponential.GetParameter(1), 1, 1)
model.SetParNames('c0', 'c1', 'mu')
model.Print()
model.SetLineColor(ROOT.kBlue)
canvas = TCanvas()
canvas.cd()
gaussian_hist.SetMinimum(0)
gaussian_hist.SetMaximum(1.4 * model.GetMaximum(x_range[0], x_range[1]))
gaussian_hist.SetTitle('Results from currying')
gaussian_hist.Draw()
rand_hist.Draw('SAME')
falling_exponential.Draw('SAME')
model.Draw('SAME')
canvas.SaveAs('currying_example.pdf')
ROOT Version: 6.11