Python RooAddPdf.paramOn方法代码示例

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]

#.........这里部分代码省略.........
        x = RooRealVar("x","M_{#mu#mu} (GeV)",minM_fit,maxM_fit)
        mean = RooRealVar("mean","mean",91.19,87.,94.)
        meanCB = RooRealVar("meanCB","meanCB",0.,-10.,10.)
        meanCB.setConstant(True)
        width = RooRealVar("width","width",2.4952,2.3,2.6)
        width.setConstant(True)
        sigma = RooRealVar("sigma","sigma",1.3,0.001,3.)
        #    sigma.setConstant(True)
        slope = RooRealVar("slope","slope",-0.1,-1.0,0.)
        #    slope.setConstant(True)
        alpha = RooRealVar("alpha","alpha",1.,0.,30.)
        #    alpha.setConstant(True)
        N = RooRealVar("N","N",2.,0.,100.)
        #    N.setConstant(True)
        fsig = RooRealVar("fsig","fsig",0.9,0.,1.0)
        
        ## PDFs
        relBW = RooGenericPdf("relBW","relBW","@0/(pow(@0*@[email protected]*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",RooArgList(x,mean,width))
        CB = RooCBShape("CB","CB",x,meanCB,sigma,alpha,N)
        expo = RooExponential("expo","expo",x,slope)
        relBWTimesCB = RooFFTConvPdf("relBWTimesCB","relBWTimesCB",x,relBW,CB)
        relBWTimesCBPlusExp = RooAddPdf("relBWTimesCBPlusExp","relBWTimesCBPlusExp",relBWTimesCB,expo,fsig)

        # Fit
        frame = x.frame()
        h = histos[i]
        # h.Rebin(10)
        h.Sumw2()
        nbin = h.GetNbinsX()
        dh = RooDataHist("dh","dh",RooArgList(x),h)
        dh.plotOn(frame)
        relBWTimesCBPlusExp.fitTo(dh)
        relBWTimesCBPlusExp.plotOn(frame)
        relBWTimesCBPlusExp.paramOn(frame)

        # Plot
        c = TCanvas("c_"+i,"c_"+i) 
        c.SetFillColor(0)
        c.cd()
        frame.Draw()
        c.SaveAs(options.outDir+"/DimuonWithFit_"+i+".png")
        
        # Extract the result of the fit
        expParams = []
        expCoef = slope.getValV()
        fSig    = fsig.getValV()
        binLow = h.GetXaxis().FindBin(minM_fit)
        binHigh = h.GetXaxis().FindBin(maxM_fit)
        nEntries = h.Integral(binLow,binHigh)
        expParams = [expCoef,fSig,nEntries,binLow,binHigh]
        expPars[i] = expParams

        signalParams = [mean.getVal(),width.getVal(),meanCB.getVal(),sigma.getVal(),alpha.getVal(),N.getVal()]
        signalPars[i] = signalParams

        # Subtract the bkg from the histograms
        h_woBkg = h.Clone()
        h_bkg = TH1F("h_bkg_"+i,"Histogram of bkg events",nbin,xlow,xup)
        
        h_bkg.Sumw2()
        expNorm = (math.fabs(expCoef)*(1-fSig)*nEntries)/(math.exp(expCoef*minM_fit)-math.exp(expCoef*maxM_fit))
        for ibin in range(binLow,binHigh):
            w = integrateExp(expNorm,expCoef,h_bkg.GetBinLowEdge(ibin),h_bkg.GetBinLowEdge(ibin+1))
            h_bkg.SetBinContent(ibin,w)

        h_woBkg.Add(h_bkg,-1)

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]

#.........这里部分代码省略.........
    frame.GetYaxis().SetLabelSize(0.04)
    frame.SetLineWidth(1)
    frame.SetTitle(plotTitle) 
    
    # plot things on frame
    data.plotOn(frame, RooFit.MarkerSize(0.7))
    chib1P_set = RooArgSet(chib1_pdf)
    modelPdf.plotOn(frame,RooFit.Components(chib1P_set), RooFit.LineColor(ROOT.kGreen+2), RooFit.LineStyle(2), RooFit.LineWidth(1))
    chib2P_set = RooArgSet(chib2_pdf)
    modelPdf.plotOn(frame, RooFit.Components(chib2P_set),RooFit.LineColor(ROOT.kRed), RooFit.LineStyle(2), RooFit.LineWidth(1))
    background_set =  RooArgSet(background)
    modelPdf.plotOn(frame,RooFit.Components(background_set), RooFit.LineColor(ROOT.kBlack), RooFit.LineStyle(2), RooFit.LineWidth(1))
    modelPdf.plotOn(frame, RooFit.LineWidth(2))
    frame.SetName("fit_resonance")  

    # Make numChib object
    numChib = NumChib(numChib=n_chib.getVal(), s_numChib=n_chib.getError(), ratio_21=ratio_21.getVal(), s_ratio_21=ratio_21.getError(), numBkg=n_background.getVal(), s_numBkg=n_background.getError(), corr_NB=result.correlation(n_chib, n_background),corr_NR=result.correlation(n_chib, ratio_21) , name='numChib'+output_suffix+ptBin_label,q0=q0.getVal(),s_q0=q0.getError(),alpha=alpha.getVal(),s_alpha=alpha.getError(), beta=beta.getVal(), s_beta=beta.getError(), chiSquare=frame.chiSquare())
    #numChib.saveToFile('numChib'+output_suffix+'.txt')

    if noPlots:
        chi2 = frame.chiSquare()
        del frame
        return numChib, chi2
    
    # Legend
    parameters_on_legend = RooArgSet()#n_chib, ratio_21, n_background)
    if massFreeToChange:
        #parameters_on_legend.add(scale_mean)
        parameters_on_legend.add(mean_1)
        #parameters_on_legend.add(mean_2)
    if sigmaFreeToChange:
        parameters_on_legend.add(scale_sigma)
    if massFreeToChange or sigmaFreeToChange:
        modelPdf.paramOn(frame, RooFit.Layout(0.1,0.6,0.2),RooFit.Parameters(parameters_on_legend))
    
    if printLegend: #chiquadro, prob, numchib...
        leg = TLegend(0.48,0.75,0.97,0.95)
        leg.SetBorderSize(0)
        #leg.SetTextSize(0.04)
        leg.SetFillStyle(0)
        chi2 = frame.chiSquare()
        probChi2 = TMath.Prob(chi2*ndof, ndof)
        chi2 = round(chi2,2)
        probChi2 = round(probChi2,2)
        leg.AddEntry(0,'#chi^{2} = '+str(chi2),'')
        leg.AddEntry(0,'Prob #chi^{2} = '+str(probChi2),'')
        N_bkg, s_N_bkg = roundPair(numChib.numBkg, numChib.s_numBkg)
        leg.AddEntry(0,'N_{bkg} = '+str(N_bkg)+' #pm '+str(s_N_bkg),'')
        N_chib12, s_N_chib12 = roundPair(numChib.numChib, numChib.s_numChib)
        leg.AddEntry(0,'N_{#chi_{b}} = '+str(N_chib12)+' #pm '+str(s_N_chib12),'')
        Ratio = numChib.calcRatio()
        s_Ratio = numChib.calcRatioError()
        Ratio, s_Ratio = roundPair(Ratio, s_Ratio)
        leg.AddEntry(0,'N_{2}/N_{1} = '+str(Ratio)+' #pm '+str(s_Ratio),'')

        if printSigReso: # Add Significance
            Sig =  numChib.calcSignificance()
            s_Sig = numChib.calcSignificanceError()
            Sig, s_Sig = roundPair(Sig, s_Sig)
            leg.AddEntry(0,'Sig = '+str(Sig)+' #pm '+str(s_Sig),'')
            if(Chib1_parameters.sigma>Chib2_parameters.sigma):
                Reso = Chib1_parameters.sigma * 1000 # So it's in MeV and not in GeV
                s_Reso = Chib1_parameters.s_sigma * 1000 # So it's in MeV and not in GeV
            else:
                Reso = Chib2_parameters.sigma * 1000 # So it's in MeV and not in GeV
                s_Reso = Chib2_parameters.s_sigma * 1000 # So it's in MeV and not in GeV

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]

#.........这里部分代码省略.........
    alpha1    =  0.95
    alpha2    =  1.12

    n         =  2.5  


    mass1_v   = RooRealVar('mchi1','m_{#chi1}',mass_chib)
    deltaM_v  = RooRealVar('deltaM','#Delta_{m}',deltaM,0.005,0.015)
    mass2_v   = RooFormulaVar('mchi2','@[email protected]',RooArgList(mass1_v,deltaM_v))
    sigma1_v  = RooRealVar('sigma1','#sigma_1',sigma1)
    sigma2_v  = RooRealVar('sigma2','#sigma_2',sigma2)

    alpha1_v  = RooRealVar('alpha1','#alpha_1',alpha1)
    alpha2_v  = RooRealVar('alpha2','#alpha_2',alpha2)

    n_v       = RooRealVar('n','n',n)

    ratio21_v = RooRealVar('ratio21','r_{21}',ratio21)

    
    x = RooRealVar("invm3S","#chi_{b} Data",10.4,10.7)

    # choose here binning of mass plot
    x.setBins(150)


    #signal pdf
    chib1 = RooCBShape('chib1','chib1',x,mass1_v,sigma1_v,alpha1_v,n_v)
    chib2 = RooCBShape('chib2','chib2',x,mass2_v,sigma2_v,alpha2_v,n_v)
 
    
    # define background
    q01S_Start = 10.4
    alpha =    RooRealVar("#alpha","#alpha",1.5,0.2,3.5)
    beta =     RooRealVar("#beta","#beta",-2.5,-7.,0.)
    #q0   =      RooRealVar("q0","q0",q01S_Start,q01S_Start-0.05,q01S_Start+0.05)
    q0   =      RooRealVar("q0","q0",q01S_Start)
    delta =     RooFormulaVar("delta","TMath::Abs(@[email protected])",RooArgList(x,q0))
    b1 =        RooFormulaVar("b1","@0*(@[email protected])",RooArgList(beta,x,q0))
    signum1 =   RooFormulaVar( "signum1","( TMath::Sign( -1.,@[email protected] )+1 )/2.", RooArgList(x,q0) )

    background = RooGenericPdf("background","Background", "signum1*pow(delta,#alpha)*exp(b1)", RooArgList(signum1,delta,alpha,b1) )


 
 
    n_evts_1 = RooRealVar('N_{3P_{1}}','N_{3P_{1}}',50,30,1000)
    n_evts_2 = RooFormulaVar('N_{3P_{2}}','@0*@1',RooArgList(n_evts_1,ratio21_v))
    n_bck    = RooRealVar('nbkg','n_{bkg}',500,0,100000)


    #build final pdf
    modelPdf = RooAddPdf('ModelPdf', 'ModelPdf', RooArgList(chib1,chib2,background),RooArgList(n_evts_1,n_evts_2,n_bck))
    
    # fit
    low_cut = x.setRange("low_cut",10.4,10.7)
    result = modelPdf.fitTo(roodataset, RooFit.Save(), RooFit.Range("low_cut") )
   
    frame = x.frame(RooFit.Title("m(#chi_{b}(3P))"))
    roodataset.plotOn(frame, RooFit.MarkerSize(0.7))
    modelPdf.plotOn(frame, RooFit.LineWidth(1))


    modelPdf.plotOn(frame, RooFit.LineWidth(2) )

    frame.GetXaxis().SetTitle('m_{#gamma #mu^{+} #mu^{-}} - m_{#mu^{+} #mu^{-}} + m^{PDG}_{#Upsilon(3S)}  [GeV/c^{2}]' )
    #frame.GetYaxis().SetTitle( "Events/15.0 MeV " )
    frame.GetXaxis().SetTitleSize(0.04)
    frame.GetYaxis().SetTitleSize(0.04)
    frame.GetXaxis().SetTitleOffset(1.1)
    frame.GetXaxis().SetLabelSize(0.04)
    frame.GetYaxis().SetLabelSize(0.04)

    frame.SetLineWidth(1)
    frame.SetName("fit_resonance")

    chi2 = frame.chiSquare()
    chi2 = round(chi2,2)
    leg=TLegend(0.50,0.7,0.60,0.8)
    leg.AddEntry(0,'#chi^{2} ='+str(chi2),'')
    leg.SetBorderSize(0)
    leg.SetFillColor(0)
    leg.SetTextSize(0.06)

    gROOT.SetStyle("Plain")

    frame.SaveAs(str(hname) + '.root')

#   param_set = RooArgSet(n_evts_Roo4, m_chib[1][3],alpha, beta, q0)

    canvas = TCanvas('fit', "", 1400, 700 )
    canvas.Divide(1)
    canvas.cd(1)
    gPad.SetRightMargin(0.3)
    gPad.SetFillColor(10)
#   modelPdf.paramOn(frame, RooFit.Layout(0.725,0.9875,0.9), RooFit.Parameters(param_set))
    modelPdf.paramOn(frame, RooFit.Layout(0.725,0.9875,0.9))
    frame.Draw()
    leg.Draw("same")
    canvas.SaveAs( str(hname) + '.png' )

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]
def fitChicSpectrum(dataset,binname):
    """ Fit chic spectrum"""


    x = RooRealVar('Qvalue','Q',9.7,10.1)
    x.setBins(80)




    mean_1 = RooRealVar("mean_1","mean ChiB1",9.892,9,10,"GeV")
    sigma_1 = RooRealVar("sigma_1","sigma ChiB1",0.0058,'GeV')
    a1_1 = RooRealVar('#alpha1_1', '#alpha1_1', 0.748)
    n1_1 = RooRealVar('n1_1', 'n1_1',2.8 )
    a2_1 = RooRealVar('#alpha2_1', '#alpha2_1',1.739)
    n2_1 = RooRealVar('n2_1', 'n2_1', 3.0)


    deltam = RooRealVar('deltam','deltam',0.01943)
    
    mean_2 = RooFormulaVar("mean_2","@[email protected]", RooArgList(mean_1,deltam))
    sigma_2 = RooRealVar("sigma_2","sigma ChiB2",0.0059,'GeV')
    a1_2 = RooRealVar('#alpha1_2', '#alpha1_2', 0.738)
    n1_2 = RooRealVar('n1_2', 'n1_2', 2.8)
    a2_2 = RooRealVar('#alpha2_2', '#alpha2_2', 1.699)
    n2_2 = RooRealVar('n2_2', 'n2_2', 3.0)

    
    parameters=RooArgSet()
    
    parameters.add(RooArgSet(sigma_1, sigma_2))
    parameters = RooArgSet(a1_1, a2_1, n1_1, n2_1)
    parameters.add(RooArgSet( a1_2, a2_2, n1_2, n2_2))
 
    chib1_pdf = My_double_CB('chib1', 'chib1', x, mean_1, sigma_1, a1_1, n1_1, a2_1, n2_1)
    chib2_pdf = My_double_CB('chib2', 'chib2', x, mean_2, sigma_2, a1_2, n1_2, a2_2, n2_2)

    
    #background
    q01S_Start = 9.5
    alpha   =   RooRealVar("#alpha","#alpha",1.5,-1,3.5)#0.2 anziche' 1
    beta    =   RooRealVar("#beta","#beta",-2.5,-7.,0.)
    q0      =   RooRealVar("q0","q0",q01S_Start)#,9.5,9.7)
    delta   =   RooFormulaVar("delta","TMath::Abs(@[email protected])",RooArgList(x,q0))
    b1      =   RooFormulaVar("b1","@0*(@[email protected])",RooArgList(beta,x,q0))
    signum1 =   RooFormulaVar( "signum1","( TMath::Sign( -1.,@[email protected] )+1 )/2.", RooArgList(x,q0) )
    
    
    background = RooGenericPdf("background","Background", "signum1*pow(delta,#alpha)*exp(b1)", RooArgList(signum1,delta,alpha,b1) )

    parameters.add(RooArgSet(alpha, beta, q0))

    #together
    chibs = RooArgList(chib1_pdf,chib2_pdf,background)    

    

    n_chib = RooRealVar("n_chib","n_chib",2075, 0, 100000)
    ratio_21 = RooRealVar("ratio_21","ratio_21",0.5,0,1)
    n_chib1 = RooFormulaVar("n_chib1","@0/([email protected])",RooArgList(n_chib, ratio_21))
    n_chib2 = RooFormulaVar("n_chib2","@0/(1+1/@1)",RooArgList(n_chib, ratio_21))
    n_background = RooRealVar('n_background','n_background',4550, 0, 50000)
    ratio_list = RooArgList(n_chib1, n_chib2, n_background)


    modelPdf = RooAddPdf('ModelPdf', 'ModelPdf', chibs, ratio_list)


    frame = x.frame(RooFit.Title('m'))
    range = x.setRange('range',9.7,10.1)
    result = modelPdf.fitTo(dataset,RooFit.Save(),RooFit.Range('range'))
    dataset.plotOn(frame,RooFit.MarkerSize(0.7))

    modelPdf.plotOn(frame, RooFit.LineWidth(2) )

    
    #plotting
    canvas = TCanvas('fit', "", 1400, 700 )
    canvas.Divide(1)
    canvas.cd(1)
    gPad.SetRightMargin(0.3)
    gPad.SetFillColor(10)
    modelPdf.paramOn(frame, RooFit.Layout(0.725,0.9875,0.9))
    frame.Draw()
    canvas.SaveAs( 'out-'+binname + '.png' )

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]
def fit_gau2_che(var, dataset, title='', print_pars=False, test=False,
                 mean_=None, sigma_=None, sigma1_=None, sigmaFraction_=None):
    # define background

    c0 = RooRealVar('c0', 'constant', 0.1, -1, 1)
    c1 = RooRealVar('c1', 'linear', 0.6, -1, 1)
    c2 = RooRealVar('c2', 'quadratic', 0.1, -1, 1)
    c3 = RooRealVar('c3', 'c3', 0.1, -1, 1)

    bkg = RooChebychev('bkg', 'background pdf', var,
                       RooArgList(c0, c1, c2, c3))
    
    # define signal
    val = 5.28
    dmean = 0.05 
    valL = val - dmean
    valR = val + dmean

    if mean_ is None:
        mean = RooRealVar("mean", "mean", val, valL, valR)
    else:
        mean = RooRealVar("mean", "mean", mean_)

    val = 0.05
    dmean = 0.02
    valL = val - dmean
    valR = val + dmean

    if sigma_ is None:
        sigma = RooRealVar('sigma', 'sigma', val, valL, valR)
    else:
        sigma = RooRealVar('sigma', 'sigma', sigma_)

    if sigma1_ is None:
        sigma1 = RooRealVar('sigma1', 'sigma1', val, valL, valR)
    else:
        sigma1 = RooRealVar('sigma1', 'sigma1', sigma1_)

    peakGaus = RooGaussian("peakGaus", "peakGaus", var, mean, sigma)
    peakGaus1 = RooGaussian("peakGaus1", "peakGaus1", var, mean, sigma1)    
    
    if sigmaFraction_ is None:
        sigmaFraction = RooRealVar("sigmaFraction", "Sigma Fraction", 0.5, 0., 1.)
    else:
        sigmaFraction = RooRealVar("sigmaFraction", "Sigma Fraction", sigmaFraction_)

    glist = RooArgList(peakGaus, peakGaus1)
    peakG = RooAddPdf("peakG","peakG", glist, RooArgList(sigmaFraction))
    
    listPeak = RooArgList("listPeak")
    
    listPeak.add(peakG)
    listPeak.add(bkg)
    
    fbkg = 0.45
    nEntries = dataset.numEntries()

    val=(1-fbkg)* nEntries
    listArea = RooArgList("listArea")
    
    areaPeak = RooRealVar("areaPeak", "areaPeak", val, 0.,nEntries)
    listArea.add(areaPeak)

    nBkg = fbkg*nEntries
    areaBkg = RooRealVar("areaBkg","areaBkg", nBkg, 0.,nEntries)
    
    listArea.add(areaBkg)
    model = RooAddPdf("model", "fit model", listPeak, listArea)

    if not test:
        fitres = model.fitTo(dataset, RooFit.Extended(kTRUE),
                             RooFit.Minos(kTRUE),RooFit.Save(kTRUE))

    nbins = 35
    frame = var.frame(nbins)

    frame.GetXaxis().SetTitle("B^{0} mass (GeV/c^{2})")
    frame.GetXaxis().CenterTitle()
    frame.GetYaxis().CenterTitle()
    frame.SetTitle(title)

    mk_size = RooFit.MarkerSize(0.3)
    mk_style = RooFit.MarkerStyle(kFullCircle)
    dataset.plotOn(frame, mk_size, mk_style)

    model.plotOn(frame)

    as_bkg = RooArgSet(bkg)
    cp_bkg = RooFit.Components(as_bkg)
    line_style = RooFit.LineStyle(kDashed)
    model.plotOn(frame, cp_bkg, line_style)

    if print_pars:
        fmt = RooFit.Format('NEU')
        lyt = RooFit.Layout(0.65, 0.95, 0.92)
        param = model.paramOn(frame, fmt, lyt)
        param.getAttText().SetTextSize(0.02)
        param.getAttText().SetTextFont(60)
    
    frame.Draw()
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]
# In[10]:


tot.fitTo(dh)


# In[11]:


massFrame = mass.frame()
massFrame.SetTitle("Phi signal")
dh.plotOn(massFrame)
tot.plotOn(massFrame)
gauss.plotOn(massFrame,LineColor(kGreen),LineStyle(kDashed),Normalization((sFrac.getValV()*numEvts)/(numEvts)))
cheb.plotOn(massFrame,LineColor(kMagenta),LineStyle(kDotted),Normalization(((1.0-sFrac.getValV())*numEvts)/(numEvts)))
tot.paramOn(massFrame,Layout(0.60,0.99,0.75));
massFrame.Draw()


# In[12]:


plotmax = hist.GetMaximum()*1.05
sidesigma = sigma.getValV()
leftlowside = -7.*sidesigma + mean.getValV()
leftupside = -5.*sidesigma + mean.getValV()
rightlowside = +5.*sidesigma + mean.getValV()
rightupside = +7.*sidesigma + mean.getValV()

signallow = -3.*sidesigma + mean.getValV()
signalup = +3.*sidesigma + mean.getValV()

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import paramOn [as 别名]
#nk1 = model3.createNLL(data)
#RooMinuit(nk1).migrad() 
#Rk1 = k2.frame()

#nk1.plotOn(Rk1,RooFit.ShiftToZero()) 
#cRcc = TCanvas("Rcc", "Rcc", 500, 500)

#Rk1.SetMaximum(4.);Rk1.SetMinimum(0)
#Rk1.GetXaxis().SetTitle("nttbar/nmc")
#Rk1.SetTitle("")
#Rk1.Draw()

cR11 = TCanvas("R11", "R", 500, 500)
xframe = x.frame()
data.plotOn(xframe, RooFit.DataError(RooAbsData.SumW2) ) 
model3.paramOn(xframe, RooFit.Layout(0.65,0.9,0.9) )
model3.plotOn(xframe)
chi2 = xframe.chiSquare(2)
ndof = xframe.GetNbinsX()
print "chi2 = "+ str(chi2)
print "ndof = "+ str(ndof)
xframe.Draw()



print "k1:"+str(k1.getVal())+", err:"+str(k1.getError())+", init:"+str(rttbar)
print "k2:"+str(k2.getVal())+", err:"+str(k2.getError())
n_mctotal = n_ttbar+n_background

print "####################"
print "ttbar = " + str(n_ttbar)