【高能物理.ROOT】ROOT高斯拟合

ROOT拟合

如有错误，欢迎指正。
参考类似:假设检验
ROOT官方文档里面有很多关于ROOT的教程。其中就有关于拟合的示例代码，选取一两个看看就知道了。

简介

有时候有些模型的部分参数是未知的，需要通过拟合得到。一般的拟合是找到最大似然函数最大的位置作为参数的估计值。

• 最小二乘法（least squares）是一种常用的拟合方法，它通过最小化残差平方和来拟合数据。实际上就是假设了模型的误差项服从高斯分布，并通过极大似然估计的方法求解参数。没有对模型其他参数有限制。
• 卡方拟合(chi square fitting)跟最小二乘法类似,实际上是一种特殊的最大似然函数,只不过假设模型的统计误差服从高斯分布.

一般的高能物理是计数实验,计数的统计误差往往服从泊松分布,一般在计数小的时候，此时最大似然函数拟合中“正常项”就是泊松分布的样子.计数大的时候泊松分布趋近于高斯分布。用chi2拟合也没差了。

对于拟合的类型有很多种

• 一维拟合:一个数据只有一个值,这批数据符合某个分布,这个分布的参数是要估计的。
• 二维拟合:一个数据有两个个值。正常的直线拟合说的就是这种

拟合例子及代码

一维拟合高斯分布

拟合模型建立

    #模型建立
    m = ROOT.RooRealVar("m", "m", (变量范围[1]+变量范围[2])/2.0, 变量范围[1], 变量范围[2])
    s = ROOT.RooRealVar("s", "s", s0, 0, (变量范围[2]-变量范围[1])*2)
    # 变量=ROOT.RooRealVar(变量名, 变量名, 变量范围[1], 变量范围[2])
    model = ROOT.RooGaussian("Gauss", "Gauss", 变量, m, s)
    if type(fixm)!=type(None):
        m.setVal(fixm)
        m.setConstant(True)

拟合

    fitpara=model.fitTo(DATASet, Save=True)

结果输出

    frame = 变量.frame(Title=变量名, Bins=变量范围[0])
    DATASet.plotOn(frame)
    model.plotOn(frame,VisualizeError=(fitpara, 2),FillColor="kOrange")
    model.plotOn(frame,LineColor="g")

    if os.path.exists(保存位置):
        cforfit=ROOT.TCanvas(变量名, 变量名, 800, 800)
        ROOT.gPad.SetLeftMargin(0.15)
        frame.GetYaxis().SetTitleOffset(1.4)
        if 画图:frame.Draw()
        cforfit.SaveAs(os.path.join(保存位置,'fit_'+变量名+'.png'))

    print("minimized -log(L) value=",fitpara.minNll())
    print("estimated distance to minimum = ", fitpara.edm())
    result=[]
    for x in fitpara.floatParsFinal():
        print("拟合结果",x.getTitle(),"=",x.getValV(),"+-",x.getError())
        result.append([x.getTitle(),x.getValV(),x.getError()])
    print("cor:",fitpara.correlationMatrix())
    if type(fixm)==type(None):
        COR=fitpara.correlationMatrix()
        for i in range(2):
            for j in range(2):
                print(COR[i][j],end='\t')
            print()
        print("cov:",fitpara.covarianceMatrix())
        COV=fitpara.covarianceMatrix()
        for i in range(2):
            for j in range(2):
                print(COV[i][j],end='\t')
            print()

完整代码

import ROOT
import os

def 高斯拟合(DATASet,变量名,变量,变量范围,保存位置='',fixm=None,s0=1.0,画图=True):
    '''
    拟合其中的某个变量。拟合模型是高斯。
    DATASet:ROOT.RooDataSet
    变量名:'Mofll'
    变量范围:(bins,xmin,xmax)
    输出:
    '''
    #模型建立
    m = ROOT.RooRealVar("m", "m", (变量范围[1]+变量范围[2])/2.0, 变量范围[1], 变量范围[2])
    s = ROOT.RooRealVar("s", "s", s0, 0, (变量范围[2]-变量范围[1])*2)
    # 变量=ROOT.RooRealVar(变量名, 变量名, 变量范围[1], 变量范围[2])
    model = ROOT.RooGaussian("Gauss", "Gauss", 变量, m, s)
    if type(fixm)!=type(None):
        m.setVal(fixm)
        m.setConstant(True)

    fitpara=model.fitTo(DATASet, Save=True)

    frame = 变量.frame(Title=变量名, Bins=变量范围[0])
    DATASet.plotOn(frame)
    model.plotOn(frame,VisualizeError=(fitpara, 2),FillColor="kOrange")
    model.plotOn(frame,LineColor="g")

    if os.path.exists(保存位置):
        cforfit=ROOT.TCanvas(变量名, 变量名, 800, 800)
        ROOT.gPad.SetLeftMargin(0.15)
        frame.GetYaxis().SetTitleOffset(1.4)
        if 画图:frame.Draw()
        cforfit.SaveAs(os.path.join(保存位置,'fit_'+变量名+'.png'))

    print("minimized -log(L) value=",fitpara.minNll())
    print("estimated distance to minimum = ", fitpara.edm())
    result=[]
    for x in fitpara.floatParsFinal():
        print("拟合结果",x.getTitle(),"=",x.getValV(),"+-",x.getError())
        result.append([x.getTitle(),x.getValV(),x.getError()])
    print("cor:",fitpara.correlationMatrix())
    if type(fixm)==type(None):
        COR=fitpara.correlationMatrix()
        for i in range(2):
            for j in range(2):
                print(COR[i][j],end='\t')
            print()
        print("cov:",fitpara.covarianceMatrix())
        COV=fitpara.covarianceMatrix()
        for i in range(2):
            for j in range(2):
                print(COV[i][j],end='\t')
            print()
    return result

if __name__ == '__main__':
    # 产生数据
    x = ROOT.RooRealVar("x", "x", -10, 10)
    
    设定均值 = ROOT.RooRealVar("m", "m", 1.1)
    设定标准差 = ROOT.RooRealVar("s", "s", 2)
    设定模型 = ROOT.RooGaussian("model", "model", x, 设定均值, 设定标准差)
    
    DATASet = 设定模型.generate({x}, 1000)#没分bin数据
    x.setBins(25)
    
    # 进行拟合
    高斯拟合(DATASet,"x",x,(25,-10,10),
         保存位置='/home/csq/代码/vs_code/md文件/testcode/',
         画图=True)

结果

图片中黄色的是拟合参数的拟合误差$\pm 2\sigma$的模型变化行为。

输出结果：

minimized -log(L) value= 2149.229088005724
estimated distance to minimum =  1.912172983979974e-06
拟合结果 m = 1.153746032145044 +- 0.06565361253393576
拟合结果 s = 2.075960999790656 +- 0.04646477947196859
cor: Name: TMatrixTSym<double> Title: 
1.0     0.0007898421275319538
0.0007898421275319538   1.0
cov: Name: TMatrixTSym<double> Title: 
0.004310459607468626    2.4095056496145953e-06
2.4095056496145953e-06  0.0021589954668200222