User-defined function fit

[MITSUKI000]

Hi, I was trying to do a DSCB fit and the code below failed.

import math
def double_sided_crystal_ball(x, p):
    # parameters: p[0]: constants, p[1]: mean, p[2]: std 
    #             p[3]: a1(low),  p[4]: n1(low),  p[5]: a2(high)  p[6]: n2(high)
    s = (x[0] - p[1]) / p[2]
    a1, n1, a2, n2 = p[3], p[4], p[5], p[6]
    A1 = math.pow(n1 / abs(a1), n1) * math.exp(-a1**2 / 2)
    A2 = math.pow(n2 / abs(a2), n2) * math.exp(-a2**2 / 2)
    B1 = math.pow((n1 / abs(a1) - abs(a1) - s), -n1)
    B2 = math.pow((n2 / abs(a2) - abs(a2) - s), -n2)
    if s < -a1:
        result = p[0] * A1 * B1
    elif s > a2:
        result = p[0] * A2 * B2
    else:
        result = p[0] * math.exp(-s**2 / 2)
        
    return result

double_sided_crystal_ball = r.TF1("dscb", double_sided_crystal_ball, 100e3, 150e3, 7)
double_sided_crystal_ball.SetParameters(1, 125e3, 1e3, 1, 10, 1, 10)
invMass.Fit("dscb")

And I got this Error:

TypeError: none of the 2 overloaded methods succeeded. Full details:
  TFitResultPtr TH1::Fit(TF1* f1, const char* option = "", const char* goption = "", double xmin = 0, double xmax = 0) =>
    TypeError: could not convert argument 1
  TFitResultPtr TH1::Fit(const char* formula, const char* option = "", const char* goption = "", double xmin = 0, double xmax = 0) =>
    ValueError: math domain error

I also tried

invMass.Fit(double_sided_crystal_ball)

and got the same error.

[MITSUKI000]

The math domain error is due to negative roots so I modified my code as:

def double_sided_crystal_ball(x, p):
    # parameters: p[0]: constants, p[1]: mean, p[2]: std 
    #             p[3]: a1(low),  p[4]: n1(low),  p[5]: a2(high)  p[6]: n2(high)
    s = (x[0] - p[1]) / p[2]
    a1, n1, a2, n2 = p[3], p[4], p[5], p[6]
    
    print(p[4])
    
    result = 1
        
    if (n1 / abs(a1) - abs(a1) - s) > 0 and (n2 / abs(a2) - abs(a2) + s) > 0:

        A1 = math.pow((n1 / abs(a1)), n1) * math.exp(-a1**2 / 2)
        A2 = math.pow((n2 / abs(a2)), n2) * math.exp(-a2**2 / 2)
        B1 = math.pow((n1 / abs(a1) - abs(a1) - s), -n1)
        B2 = math.pow((n2 / abs(a2) - abs(a2) + s), -n2)

        if s < -a1:
            result = p[0] * A1 * B1
        elif s > a2:
            result = p[0] * A2 * B2
        else:
            result = p[0] * math.exp(-s**2 / 2)

    return result

fit1 = r.TF1("fit1", double_sided_crystal_ball, 120e3, 130e3, 3 )
fit1.SetParameters(1, 125e3, 1e3, 1, 10, 1, 10)
fit1.FixParameter(4, 1)
fit1.FixParameter(5, 1)
fit1.SetParLimits(6, 1, 2)
fit1.SetParLimits(7, 0.1, 50)
invMass_MC.Fit("fit1", "B")
invMass_MC.Draw("E")

However, the methods SetParameters and SetParLimits are not working, the printed values of p[4] is always the following values

[Wile_E_Coyote]

fit1 = r.TF1("fit1", double_sided_crystal_ball, 120e3, 130e3, 7)

BTW. Parameters are numbered from 0 not from 1.

[MITSUKI000]

Thanks it works now