Python Fitter.getChisquareOverDoF方法代码示例

# 需要导入模块: from fitter import Fitter [as 别名]
# 或者: from fitter.Fitter import getChisquareOverDoF [as 别名]
def evalAngles():
    # load data and set errors
    datalist = loadCSVToList('../data/angles.txt')
    l = len(datalist)
    data = DataErrors().fromLists(list(zip(*datalist)[0]), list(zip(*datalist)[1]), ex=[0.5] * l, ey=[0] * l)
    data.setYErrorFunc(lambda x: np.sqrt(x))
    
    #draw graph
    c = TCanvas('c', '', 1280, 720)
    g = data.makeGraph('g', 'Winkel #alpha / #circ', 'Counts N')
    g.SetMarkerStyle(1)
    g.Draw('AP')
    
    #fit function
    fit = Fitter('f', '[0]+gaus(1)')
    fit.function.SetNpx(1000)  # for smoother curve
    fit.setParam(0, 'offset', 35)
    fit.setParam(1, 'ampl', 350)
    fit.setParam(2, 'theta', 180)
    fit.setParam(3, 'sigma', 20)
    fit.fit(g, 80, 280)
    fit.saveData('../calc/angles.txt', 'w')

    fit2 = Fitter('f', '[0]+1/2*[4]*(TMath::Erf(([1]+2*x-2*[2])/(2*sqrt(2)*[3]))+TMath::Erf(([1]-2*x+2*[2])/(2*sqrt(2)*[3])))')
    fit2.function.SetNpx(1000)  # for smoother curve
    fit2.function.SetLineColor(4)
    fit2.setParam(0, 'offset', 20)
    fit2.setParam(1, 'breite', 20)
    fit2.setParam(2, 'theta', 180)
    fit2.setParam(3, 'sigma', 5)
    fit2.setParam(4, 'amplitude', 350)
    fit2.fit(g, 80, 280, '+')
    fit2.saveData('../calc/angles_convolution.txt', 'w')
    
    l = TLegend(0.65, 0.55, 0.85, 0.85)
    l.AddEntry('g', 'Messwerte', 'p')
    l.AddEntry(fit.function, 'Fit mit Gaussverteilung', 'l')
    l.AddEntry(0, '#chi^2 / DoF = %.1f' % fit.getChisquareOverDoF(), '')
    l.AddEntry(0, '#alpha_{0,g} = (%.1f #pm %.1f) #circ' % (fit.params[2]['value'], fit.params[2]['error']), '')
    l.AddEntry(fit2.function, 'Fit mit Faltungsprodukt', 'l')
    l.AddEntry(0, '#chi^2 / DoF = %.1f' % fit2.getChisquareOverDoF(), '')
    l.AddEntry(0, '#alpha_{0,f} = (%.1f #pm %.1f) #circ' % (fit2.params[2]['value'], fit2.params[2]['error']), '')
    l.Draw()
    
    #print to file
    c.Update()
    c.Print('../img/angles.pdf')

# 需要导入模块: from fitter import Fitter [as 别名]
# 或者: from fitter.Fitter import getChisquareOverDoF [as 别名]
def evalCo():
    data = SzintData.fromPath('../data/co.TKA')
    data.prepare()

    c = TCanvas('c', '', 1280, 720)
    g = data.makeGraph('g', 'Kanalnummer', 'Z#ddot{a}hlrate / (1/s)')
    prepareGraph(g)
    g.SetMinimum(-0.01)
    g.GetXaxis().SetRangeUser(0, 8200)
    g.Draw('APX')

    c.Update()
    c.Print('../img/co_spectrum.pdf')

    fit = Fitter('f', 'pol1(0) + gaus(2) + gaus(5)')
    fit.function.SetLineWidth(2)
    fit.setParam(0, 'a', 0)
    fit.setParam(1, 'b', 0)
    fit.setParam(2, 'A1', 0.1)
    fit.setParam(3, 'c1', 2800)
    fit.setParam(4, 's1', 50)
    fit.setParam(5, 'A2', 0.1)
    fit.setParam(6, 'c2', 3250)
    fit.setParam(7, 's2', 75)
    fit.fit(g, 2650, 3500)
    fit.saveData('../calc/co_peaks.txt', 'w')

    with TxtFile('../calc/energy_co.txt', 'w') as f:
        f.writeline('\t', '1172', str(fit.params[3]['value']), str(fit.params[3]['error']))
        f.writeline('\t', '1332', str(fit.params[6]['value']), str(fit.params[6]['error']))

    l = TLegend(0.635, 0.62, 0.87, 0.87)
    l.AddEntry('g', 'Messung', 'p')
    l.AddEntry(fit.function, 'Fit mit y = a + b*x + gaus(x; A1, c1, s1)', 'l')
    l.AddEntry(fit.function, '                  + gaus(x; A2, c2, s2)', '')
    l.AddEntry(0, '#chi^{2} / DoF : %f' % fit.getChisquareOverDoF(), '')
    l.AddEntry(0, 'Peaks:', '')
    l.AddEntry(0, 'c1 = %.2f #pm %.2f' % (fit.params[3]['value'], fit.params[3]['error']), '')
    l.AddEntry(0, 'c2 = %.2f #pm %.2f' % (fit.params[6]['value'], fit.params[6]['error']), '')
    l.Draw()

    g.GetXaxis().SetRangeUser(2400, 3800)
    g.Draw('P')
    c.Update()
    c.Print('../img/co_peaks.pdf')

# 需要导入模块: from fitter import Fitter [as 别名]
# 或者: from fitter.Fitter import getChisquareOverDoF [as 别名]
def makeAreaFit():
    # calculate ares
    d_s = [1.000, 0.990, 0.990, 1.005, 1.000, 1.005]
    d_m = [1.700, 1.690, 1.695, 1.700, 1.705, 1.705]
    d_l = [2.880, 2.880, 2.875, 2.880, 2.880, 2.880]
    d = map(avgstd, [d_s, d_m, d_l])
    a = map(area, d)
    
    diaarea = DataErrors()
    for i in range(3):
        diaarea.addPoint(d[i][0], a[i][0], d[i][1], a[i][1])
    diaarea.saveDataToLaTeX(['Durchmesser $d$ / cm', '$s_d$ / cm', 'Fl\"ache $F / \\text{cm}^2$', '$s_F / \\text{cm}^2$'], 
                            ['%.4f', '%.4f', '%.4f', '%.4f'], 'Verschiedene Fl\"achen f\"ur die Samariummessung', 
                            'tab:data:samarium:area', '../src/data_samarium_areas.tex', 'w')
    
    with TxtFile('../fit/samarium.txt', 'w') as f:
        f.writeline('areas')
        f.writeline('=====')
        for b in a:
            f.writeline('\t', '%e' % b[0], TxtFile.PM, '%e' % b[1])
        f.writeline()
    
    
    #read data from files
    #file=[area, area error, path, time]
    files = []
    files.append([a[0][0], a[0][1], '../data/31_Sm_kl_1600_t3000.txt', 3000])
    files.append([a[1][0], a[1][1], '../data/34_Sm_m_1600_t2400.txt', 2400])
    files.append([a[2][0], a[2][1], '../data/08_Sm_ggrFl_1600-1600-0.txt', 1200])
    u = readSingleEntryFile('../data/09_Untergrund_1600-1600-0.txt') 
    tu = 3600
    d = DataErrors()
    for file in files:
        n = readSingleEntryFile(file[2])
        d.addPoint(file[0], n - u, file[1], np.sqrt(n / file[3] + u / tu))
        
    d.saveDataToLaTeX(['Fl\"ache $F / \\text{cm}^2$', '$s_F / \\text{cm}^2$', 'Z\"ahlrate $n / (1/\\text{s})$', '$s_n / (1/\\text{s})$'], 
                      ['%.4f', '%.4f', '%.3f', '%.3f'], 'Z\"ahlraten von \\samarium~f\"ur verschiedene Fl\"achen mit Fehlern', 
                      'tab:data:samarium', '../src/data_samarium.tex', 'w')
        
    c = TCanvas('c2', '', 800, 600)
    g = d.makeGraph('g', 'Fl#ddot{a}che F / cm^{2}', 'Z#ddot{a}hlrate n / (1/s)')
    g.Draw('AP')
    
    fit = Fitter('f', '[0]+[1]*x')
    fit.setParam(0, 'a', 0)
    fit.setParam(1, 'b', 0.05)
    fit.fit(g, 0, 30)
    fit.saveData('../fit/samarium.txt', 'a')
    
    a = fit.params[0]['value']
    sa = fit.params[0]['error']
    b = fit.params[1]['value']
    sb = fit.params[1]['error']

    l = TLegend(0.55, 0.15, 0.98, 0.5)
    l.AddEntry('g', '{}^{147} Samarium ohne Untergrund', 'p')  # TODO with error bar? (options +'e')
    l.AddEntry(fit.function, 'Fit mit n(F)=a+b*F', 'l')
    l.AddEntry(0, '#chi^{2} / DoF : %f' % fit.getChisquareOverDoF(), '')
    l.AddEntry(0, 'Paramter:', '') 
    l.AddEntry(0, 'a: %e #pm %e' % (a, sa), '')
    l.AddEntry(0, 'b: %e #pm %e' % (b, sb), '')   
    l.SetTextSize(0.03)
    l.Draw()
    
    c.Update()
    c.Print('../img/Samarium147-Flaechenabhaengigkeit.pdf', 'pdf')
    
    #calculation with fit parameters
    c  = 0.004025
    NA = 6.02214129e23
    m  = 2*150.36 + 3*15.999
    h = 0.1487
    t = (np.log(2) * c * NA * h) / (2 * m * b) / (3600 * 24 * 365.242)
    st = t * (sb / b)
    with TxtFile('../fit/samarium.txt', 'a') as f:
        f.writeline('calculation from fit')
        f.writeline('====================')
        f.writeline('\t', '%e' % t, TxtFile.PM, '%e' % st)
        f.writeline()
    
    # calculation from single data points
    sc = map(lambda p: calculateHalfLife(*p), d.points)
    with TxtFile('../fit/samarium.txt', 'a') as f:
        f.writeline('calculation from single data points')
        f.writeline('===================================')
        for s in sc:
            f.writeline('\t', '%e' % s[0], TxtFile.PM, '%e' % s[1])
        f.writeline()

# 需要导入模块: from fitter import Fitter [as 别名]
# 或者: from fitter.Fitter import getChisquareOverDoF [as 别名]
def makeMassFit():
    # config files
    # file = [mass, path]
    files = []
    files.append([2.0123, "../data/11_K_m9_3200_t420.txt", 420])
    files.append([2.0123, "../data/11b_K_m9_3200_t420.txt", 420])
    files.append([1.9047, "../data/13_K_m8_3200_t420.txt", 420])
    files.append([1.6812, "../data/15_K_m7_3200_t420.txt", 420])
    files.append([1.4827, "../data/17_K_m6_3200_t420.txt", 420])
    files.append([1.2952, "../data/19_K_m5_3200_t480.txt", 480])
    files.append([1.0993, "../data/21_K_m4_3200_t480.txt", 480])
    files.append([0.8086, "../data/23_K_m3_3200_t540.txt", 540])
    files.append([0.6954, "../data/25_K_m2_3200_t540.txt", 540])
    files.append([0.5007, "../data/27_K_m1_3200_t660.txt", 660])
    files.append([0.3030, "../data/29_K_m0_3200_t780.txt", 780])
    u = 0.760
    tu = 50

    d = DataErrors()

    for file in files:
        n = readSingleEntryFile(file[1])
        d.addPoint(file[0], n - u, 0.001, np.sqrt(n / file[2] + u / tu))

    d.saveDataToLaTeX(['Masse $m /$g', '$s_m /$g', 'Z\"ahlrate $n / (1/\\text{s})$', '$s_n / (1/\\text{s})$'],
                      ['%.3f', '%.3f', '%.3f', '%.3f'], 
                      'Z\"ahlraten von \\kalium\,f\"ur verschiedene Massen mit Fehlern', 'tab:data:kalium', '../src/data_kalium.tex', 'w')

    c = TCanvas('c2', '', 800, 600)
    g = d.makeGraph('g', 'Masse m / g', 'Z#ddot{a}hlrate n / (1/s)')
    g.SetMaximum(6)
    g.SetMinimum(2)
    g.Draw('AP')

    fit = Fitter('f', '[0]*(1-exp(-[1]*x))')
    fit.setParam(0, 'a')
    fit.setParam(1, 'b')
    fit.fit(g, 0.1, 2.5)
    fit.saveData('../fit/kalium.txt')
    
    a = fit.params[0]['value']
    sa = fit.params[0]['error']
    b = fit.params[1]['value']
    sb = fit.params[1]['error']

    l = TLegend(0.4, 0.15, 0.85, 0.5)
    l.AddEntry('g', '{}^{40} Kalium ohne Untergrund', 'p')  # TODO with error bar? (options +'e')
    l.AddEntry(fit.function, 'Fit mit n(m)=a(1-exp(-b*m))', 'l')
    l.AddEntry(0, '#chi^{2} / DoF : %f' % fit.getChisquareOverDoF(), '')
    l.AddEntry(0, 'Paramter:', '') 
    l.AddEntry(0, 'a: %e #pm %e' % (a, sa), '')
    l.AddEntry(0, 'b: %e #pm %e' % (b, sb), '')   
    l.SetTextSize(0.03)
    l.Draw()

    NA = 6.02214129e23
    hrel = 0.000118
    mrel = 39.0983 + 35.45
    f = 1.29
    rho = fit.getCorrMatrixElem(1, 0)
    thalf = (np.log(2) * NA * hrel * f) / (1.12 * mrel * 2 * a * b) / (3600 * 24 * 365.242)
    sthalf = thalf * np.sqrt((sa / a) ** 2 + (sb / b) ** 2 + 2 * rho * (sa / a) * (sb / b))

    with TxtFile.fromRelPath('../fit/kalium.txt', 'a') as f:
        f.writeline()
        f.writeline('calculations')
        f.writeline('============')
        f.writeline('\t', 'half-life of Kalium:', '%e' % thalf, TxtFile.PM, '%e' % sthalf)

    c.Update()
    c.Print('../img/Kalium40_Massenabhaengigkeit.pdf')

# 需要导入模块: from fitter import Fitter [as 别名]
# 或者: from fitter.Fitter import getChisquareOverDoF [as 别名]
def evalEnergyGauge():
    # get data
    datalist  = funcs.loadCSVToList('../calc/energy_na.txt')
    datalist += funcs.loadCSVToList('../calc/energy_co.txt')
    datalist += funcs.loadCSVToList('../calc/energy_eu.txt')
    
    #make latex table
    elems = ['Na', 'Na', 'Co', 'Co', 'Eu', 'Eu']
    with TxtFile('../src/energygauge.tex', 'w') as f:
        f.write2DArrayToLatexTable(zip(*([elems] + zip(*datalist))), ['Element', 'Literaturwert / keV', 'Kanal', 'Fehler auf Kanal'], 
                                   ['%s', '%.0f', '%.2f', '%.2f'], 'Referenzpeaks mit Literaturwerten', 'tab:energygauge')
    
    #convert do DataErrors
    datalist = zip(*datalist)
    data = DataErrors.fromLists(datalist[1], datalist[0], datalist[2], [0] * len(datalist[0]))
    
    c = TCanvas('c', '', 1280, 720)
    g = data.makeGraph('g', 'Kanalnummer', 'Energie / keV')
    g.GetXaxis().SetRangeUser(0, 3500)
    g.Draw('AP')
    
    fit = Fitter('f', 'pol1(0)')
    fit.function.SetNpx(1000)
    fit.setParam(0, 'a', 0)
    fit.setParam(1, 'b', 0.4)
    fit.fit(g, 0, 3500)
    fit.saveData('../calc/energy_gauge_lin.txt', 'w')
    
    l = TLegend(0.15, 0.6, 0.5, 0.85)
    l.AddEntry('g', 'Referenzpeaks', 'p')
    l.AddEntry(fit.function, 'Fit mit y = a + b*x', 'l')
    l.AddEntry(0, '#chi^{2} / DoF : %.0f' % fit.getChisquareOverDoF(), '')
    l.AddEntry(0, 'Parameter:', '')
    l.AddEntry(0, 'a = %.2f #pm %.2f' % (fit.params[0]['value'], fit.params[0]['error']), '')
    l.AddEntry(0, 'b = %.5f #pm %.5f' % (fit.params[1]['value'], fit.params[1]['error']), '')
    l.Draw()

    c.Update()
    c.Print('../img/energy_gauge_lin.pdf', 'pdf')
    
    fit2 = Fitter('f', 'pol2(0)')
    fit2.function.SetNpx(1000)
    fit2.setParam(0, 'a', 0)
    fit2.setParam(1, 'b', fit.params[1]['value'])
    fit2.setParam(2, 'c', 0)
    fit2.fit(g, 0, 3500)
    fit2.saveData('../calc/energy_gauge_lin.txt')
    
    l = TLegend(0.15, 0.575, 0.5, 0.85)
    l.AddEntry('g', 'Referenzpeaks', 'p')
    l.AddEntry(fit2.function, 'Fit mit y = a + b*x + c*x^2', 'l')
    l.AddEntry(0, '#chi^{2} / DoF : %.0f' % fit2.getChisquareOverDoF(), '')
    l.AddEntry(0, 'Parameter:', '')
    l.AddEntry(0, 'a = %.2f #pm %.2f' % (fit2.params[0]['value'], fit2.params[0]['error']), '')
    l.AddEntry(0, 'b = %.5f #pm %.5f' % (fit2.params[1]['value'], fit2.params[1]['error']), '')
    l.AddEntry(0, 'c = %.8f #pm %.8f' % (fit2.params[2]['value'], fit2.params[2]['error']), '')
    l.Draw()
    
    c.Update()
    c.Print('../img/energy_gauge_quad.pdf', 'pdf')
    
    #write raw data for reuse
    with TxtFile('../calc/energy_gauge_raw.txt', 'w') as f:
        f.writeline('\t', str(fit2.params[0]['value']), str(fit2.params[0]['error']))
        f.writeline('\t', str(fit2.params[1]['value']), str(fit2.params[1]['error']))
        f.writeline('\t', str(fit2.params[2]['value']), str(fit2.params[2]['error']))
        f.writeline(str(fit2.getCorrMatrixElem(0, 1))) #a, b
        f.writeline(str(fit2.getCorrMatrixElem(0, 2))) #a, c
        f.writeline(str(fit2.getCorrMatrixElem(1, 2))) #b, c