Python RooAddPdf.plotOn方法代码示例

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
def PlotSignalShapes(Selection):
    f__ = TFile.Open( "datacards/22June/2dPlots.root")
    signal_fname_1 = ("signals/22June/out_{sample:s}_syst.root", "cms_hgg_13TeV" )
    signal_fname_2 = ("signals/22June/out_ctcv_{sample:s}_syst.root" , "ctcv" )
    samples = {"thw":signal_fname_2, "thq":signal_fname_2, "tth":signal_fname_1 , "vh":signal_fname_1 }
    purity_h_name = "{sample:s}/"+Selection+"/h{sample:s}_"+Selection+"_purity_CtCv"
    purities = RooArgList()
    signalshapes = RooArgList()

    ctOverCvs = []

    mVar = None
    ctovercv_vals = None
    for sample in samples :
        purity = CtCvCpInfo("purity_" + sample)
        ctovercv_vals = sorted(purity.AllCtOverCVs.keys())
        purity.FillFrom2DHisto( f__.Get( purity_h_name.format( sample=sample ) ) )
        purity.GetCtOverCv()
        purities.add( purity.CtOverCvDataHistFunc )
        objsToKeep.append( purity )

        sFile = TFile.Open( samples[sample][0].format( sample=sample ) )
        ws = sFile.Get( samples[sample][1] )
        pdf = ws.pdf("RV{sample:s}_mh125".format( sample=sample) )
        objsToKeep.append(sFile)
        objsToKeep.append(ws)
        objsToKeep.append(pdf)
        signalshapes.add( pdf )

        ctOverCvs.append( ws.var( "CtOverCv" ) )
        mVar = ws.var("CMS_hgg_mass")
        
    ret = RooAddPdf("signal" , "signal" , signalshapes , purities )
    objsToKeep.append( ret )
    plot = mVar.frame()
    options = ""
    for ctovercv in ctovercv_vals :
        for var in ctOverCvs:
            var.setVal( ctovercv )
        name = "name%g" % ctovercv
        ret.plotOn( plot , RooFit.DrawOption(options) , RooFit.Name(name) )
        
        c = TCanvas()
        plot.Draw()
        c.SaveAs("a.gif+")

        if not "same" in options :
            options += " same"

    return c

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

#.........这里部分代码省略.........
    for i in ids:
        # RooFit definitions
        ## RooRealVars
        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))

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
    signal = RooAddPdf('signal','signal',sig,bkg,fsig)

    # -----------------------------------------
    # fit signal
    canSname = 'can_Mjj'+str(mass)
    if useSub:
      canSname = 'can_Sub_Mjj'+str(mass)
    canS = TCanvas(canSname,canSname,900,600)
    #gPad.SetLogy() 

    roohistSig = RooDataHist('roohist','roohist',RooArgList(x),hSig)

    signal.fitTo(roohistSig)
    frame = x.frame()
    roohistSig.plotOn(frame)
    signal.plotOn(frame)
    signal.plotOn(frame,RooFit.Components('bkg'),RooFit.LineColor(ROOT.kRed),RooFit.LineWidth(2),RooFit.LineStyle(ROOT.kDashed))
    frame.GetXaxis().SetRangeUser(900,4500)
    frame.GetXaxis().SetTitle('m_{jj} (GeV)')
    frame.Draw()

    parsSig = signal.getParameters(roohistSig)
    parsSig.setAttribAll('Constant', True)

if fitDat: 

    # -----------------------------------------
    # define parameters for background
    NBINS = 180
    p1 = RooRealVar('p1','p1',7,1,10)
    p2 = RooRealVar('p2','p2',5,1,10)

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
def dofit(roodataset, hname):
  
    mass_chib =  10.5103 # from PES uncorrected mass measurement

    deltaM    =  0.0105   #  MeV theoretical expectations
    ratio21   =  0.45     # same as chic2/chic1 and chib2/chib1
    
    # the following numbers are from an old 3P gun simulation
    # that needs to be re-done
    
    sigma1    =  0.003#0.0031 
    sigma2    =  0.003#0.0035
    
    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')
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
        dataInt3 = hData.Integral(hData.GetXaxis().FindBin(float(args.massMin)),hData.GetXaxis().FindBin(float(args.massMax)))
        background_norm3 = RooRealVar('background_norm3','background_norm3',dataInt3,0.,1e+08)

	background4 = RooGenericPdf('background4','(1/pow(@[email protected]*@0/%.1f+pow(@0/%.1f,2),@3))'%(sqrtS,sqrtS),RooArgList(mjj,p8,p9,p10))
        dataInt4 = hData.Integral(hData.GetXaxis().FindBin(float(args.massMin)),hData.GetXaxis().FindBin(float(args.massMax)))
        background_norm4 = RooRealVar('background_norm4','background_norm4',dataInt4,0.,1e+08)

	background5 = RooGenericPdf('background5','(pow(@0/%.1f,[email protected])*pow(1-pow(@0/%.1f,1/3),@2))'%(sqrtS,sqrtS),RooArgList(mjj,p11,p12))
        dataInt5 = hData.Integral(hData.GetXaxis().FindBin(float(args.massMin)),hData.GetXaxis().FindBin(float(args.massMax)))
        background_norm5 = RooRealVar('background_norm5','background_norm5',dataInt5,0.,1e+08)

	background6 = RooGenericPdf('background6','(pow(@0/%.1f,2)[email protected]*@0/%[email protected])'%(sqrtS,sqrtS),RooArgList(mjj,p13,p14))
        dataInt6 = hData.Integral(hData.GetXaxis().FindBin(float(args.massMin)),hData.GetXaxis().FindBin(float(args.massMax)))
        background_norm6 = RooRealVar('background_norm6','background_norm6',dataInt6,0.,1e+08)

	background7 = RooGenericPdf('background7','(([email protected]*@0/%.1f)*pow(@0/%.1f,@[email protected]*log(@0/%.1f)))'%(sqrtS,sqrtS,sqrtS),RooArgList(mjj,p15,p16,p17))
        dataInt7 = hData.Integral(hData.GetXaxis().FindBin(float(args.massMin)),hData.GetXaxis().FindBin(float(args.massMax)))
        background_norm7 = RooRealVar('background_norm7','background_norm7',dataInt7,0.,1e+08)

        # S+B model
        model = RooAddPdf("model","s+b",RooArgList(background,signal),RooArgList(background_norm,signal_norm))
	model2 = RooAddPdf("model2","s+b2",RooArgList(background2,signal),RooArgList(background_norm2,signal_norm))
	model3 = RooAddPdf("model3","s+b3",RooArgList(background3,signal),RooArgList(background_norm3,signal_norm))
	model4 = RooAddPdf("model4","s+b4",RooArgList(background4,signal),RooArgList(background_norm4,signal_norm))
	model5 = RooAddPdf("model5","s+b5",RooArgList(background5,signal),RooArgList(background_norm5,signal_norm))
	model6 = RooAddPdf("model6","s+b6",RooArgList(background6,signal),RooArgList(background_norm6,signal_norm))
	model7 = RooAddPdf("model7","s+b7",RooArgList(background7,signal),RooArgList(background_norm7,signal_norm))

        rooDataHist = RooDataHist('rooDatahist','rooDathist',RooArgList(mjj),hData)


        if args.runFit:
	    mframe = mjj.frame()
	    rooDataHist.plotOn(mframe, ROOT.RooFit.Name("setonedata"), ROOT.RooFit.Invisible())
	    res = model.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
	    model.plotOn(mframe, ROOT.RooFit.Name("model1"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kRed)) 
	    res2 = model2.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            model2.plotOn(mframe, ROOT.RooFit.Name("model2"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kOrange))
	    res3 = model3.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            model3.plotOn(mframe, ROOT.RooFit.Name("model3"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kGreen))
	    res4 = model4.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            model4.plotOn(mframe, ROOT.RooFit.Name("model4"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kBlue))
	    res5 = model5.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            model5.plotOn(mframe, ROOT.RooFit.Name("model5"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kViolet))
	    res6 = model6.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
#           model6.plotOn(mframe, ROOT.RooFit.Name("model6"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kPink))
	    res7 = model7.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
#           model7.plotOn(mframe, ROOT.RooFit.Name("model7"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kAzure))

	    rooDataHist2 = RooDataHist('rooDatahist2','rooDathist2',RooArgList(mjj),hData2)
	    rooDataHist2.plotOn(mframe, ROOT.RooFit.Name("data"))

	    canvas = TCanvas("cdouble", "cdouble", 800, 1000)

	    gStyle.SetOptStat(0);
            gStyle.SetOptTitle(0);
	    top = TPad("top", "top", 0., 0.5, 1., 1.)
	    top.SetBottomMargin(0.03)
	    top.Draw()
	    top.SetLogy()
            bottom = TPad("bottom", "bottom", 0., 0., 1., 0.5)
	    bottom.SetTopMargin(0.02)
	    bottom.SetBottomMargin(0.2)
	    bottom.Draw()

	    top.cd()

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

#pdf_ext_combine5.fitTo(data4, RooFit.Range("whole_range") )
pdf_ext_combine5.fitTo(data4_SB, RooFit.Range("left_side_band_region,right_side_band_region") )

print ""
print "after fit"
print "nGauss14: ", nGauss14.getVal()
print "nGauss15: ", nGauss15.getVal()
print "nGauss14 + nGauss15: ", nGauss14.getVal() + nGauss15.getVal()
print ""

# plot
#gauss12.plotOn(xframe8, RooFit.Normalization( n_generate * frac_combine4.getVal()   ,RooAbsReal.NumEvent),RooFit.LineColor(RooFit.kOrange))
#gauss13.plotOn(xframe8, RooFit.Normalization( n_generate*(1-frac_combine4.getVal() ) ,RooAbsReal.NumEvent),RooFit.LineColor(RooFit.kCyan))
pdf_combine4.plotOn(xframe8, RooFit.Normalization(n_generate ,RooAbsReal.NumEvent) ,RooFit.LineColor(RooFit.kBlue))

#data4.plotOn(xframe8)
data4_SB.plotOn(xframe8)

#pdf_ext_combine5.plotOn(xframe8, RooFit.Normalization(nGauss14.getVal() + nGauss15.getVal() ,RooAbsReal.NumEvent) ,RooFit.LineColor(RooFit.kRed))
pdf_ext_combine5.plotOn(xframe8, RooFit.Normalization(data4_SB.sumEntries() ,RooAbsReal.NumEvent) ,RooFit.LineColor(RooFit.kRed))
#pdf_ext_combine5.plotOn(xframe8, RooFit.Normalization(1.0,RooAbsReal.RelativeExpected) ,RooFit.LineColor(RooFit.kRed))

print ""
print "n_generate * frac_combine4.getVal(): ", n_generate * frac_combine4.getVal()
print "n_generate*(1-frac_combine4.getVal() ): ", n_generate*(1-frac_combine4.getVal() )
print "data4_SB: ", data4_SB.sumEntries()
print ""

# -------------------------------------------

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
                             RooFit.Normalization(sig_hist.GetEntries()/3,
                                                  RooAbsReal.Raw))
        fitter = lxgplus
    else:
        fitter = lxg
        fmip.setVal(1.0)
        fmip.setConstant(True)

    xf.Draw()
    # fitter.Print('v')
    # raise 'Stop here'
    fr = fitter.fitTo(sigDS, RooFit.Minos(False), RooFit.Save(True))
    fitter.plotOn(xf)
    if (fmip.getVal() < 0.9):
        lxgplus.plotOn(xf, RooFit.LineColor(kRed+2),
                       RooFit.LineStyle(kDashed),
                       RooFit.Components('ped*'))
        lxgplus.plotOn(xf, RooFit.LineColor(myBlue),
                       RooFit.LineStyle(kDashed),
                       RooFit.Components('lxg'))
        
    #lxgplus.paramOn(xf)
    
    ## c2 = TCanvas('c2', 'signal')
    xf.Draw()
    c2.Modified()
    c2.Update()
    ## fpar[0] = width.getVal()
    ## fpar[1] = mpv.getVal()
    ## fpar[2] = sig_hist.GetEntries()*(1-fped.getVal())
    ## fpar[3] = sigma.getVal()

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

#.........这里部分代码省略.........
    # 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)

    # Here I have as parameters N_12, ratio_12, N_background
    n_chib = RooRealVar("n_chib","n_chib",2075, 0, 100000)
    ratio_21 = RooRealVar("ratio_21","ratio_21",0.6, 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)
    parameters.add(RooArgSet(n_chib1, n_chib2, n_background))
    modelPdf = RooAddPdf('ModelPdf', 'ModelPdf', chibs, ratio_list)
    
    print 'Fitting to data'
    fit_region = x.setRange("fit_region",9.7,10.1)
    result=modelPdf.fitTo(data,RooFit.Save(), RooFit.Range("fit_region"))
    
        
    # define frame
    frame = x.frame()
    frame.SetNameTitle("fit_resonance","Fit Resonanace")
    frame.GetXaxis().SetTitle(x_axis_label )
    frame.GetYaxis().SetTitle( "Events/5 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.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)

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [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 plotOn [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 plotOn [as 别名]
print numEvts


# 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()

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
from ROOT import RooAddPdf
phiFrac = RooRealVar( 'phiFrac', 'phiFrac', 0.96, 0., 1. )
phiFrac.setError(0.003)
KKMassPdf = RooAddPdf( 'KKMassPdf', 'KKMassPdf', phiMassPdf, KKSWavePdf, phiFrac )

#Fit
fitResult = KKMassPdf.fitTo( sigData, SumW2Error = False, Save = True, Minos = False, **fitOpts )
fitResult.PrintSpecial( text = True )

#Plot
KKMassPlot = KKMassVar.frame(120)

if drawTotalDists :
    data.plotOn( KKMassPlot, MarkerStyle = kFullCircle, MarkerSize = 0.5, MarkerColor = kRed, LineWidth = 2, LineColor = kRed )
sigData.plotOn( KKMassPlot, MarkerStyle = kFullCircle, MarkerSize = 0.5, LineWidth = 2 )
KKMassPdf.plotOn( KKMassPlot, LineStyle = kSolid, LineWidth = 3, LineColor = kBlue                                   )
#KKMassPdf.plotOn( KKMassPlot, LineStyle = 7,      LineWidth = 3, LineColor = kRed,         Components = 'phiMassPdf' )
#KKMassPdf.plotOn( KKMassPlot, LineStyle = 5,      LineWidth = 3, LineColor = kMagenta + 3, Components = 'KKSWavePdf' )

binWidth = ( KKMassPlot.GetXaxis().GetXmax() - KKMassPlot.GetXaxis().GetXmin() ) / float( KKMassPlot.GetNbinsX() )
KKMassPlot.SetXTitle('m(K^{+}K^{-}) [MeV/c^{2}]')
KKMassPlot.SetYTitle('Candidates / (%.2g MeV/c^{2})' % binWidth )
KKMassPlot.SetMinimum(0.)
KKMassPlot.SetMaximum( 6200. if not drawTotalDists else 7000. )
KKMassPlot.SetTitleOffset( 1.10, 'x' )
KKMassPlot.SetTitleOffset( 1.15, 'y' )

KKMassPlot.SetMinimum(0.)

KKMassCanv = TCanvas('KKMassCanv')
KKMassCanv.SetLeftMargin(0.18)

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
 
# --- Construct signal+background PDF ---
nsig = RooRealVar("nsig", "#signal events", 200, 0., 10000)
nbkg = RooRealVar("nbkg", "#background events", 800, 0., 10000)
model = RooAddPdf("model", "g+a", RooArgList(signal, background), RooArgList(nsig, nbkg))

# --- Generate a toyMC sample from composite PDF ---
data = model.generate(RooArgSet(mes), 2000)
 
# --- Perform extended ML fit of composite PDF to toy data ---
model.fitTo(data)
 
# --- Plot toy data and composite PDF overlaid ---
mesframe = mes.frame()
data.plotOn(mesframe)
model.plotOn(mesframe)
model.plotOn(mesframe, RooFit.Components('background'), RooFit.LineStyle(kDashed))

mesframe.Draw()

print 'nsig:',nsig.getValV(), '+-', nsig.getError()
print 'nbkg:', nbkg.getValV(), '+-', nbkg.getError()
print 'mes:', mes.getValV(), '+-', mes.getError()
print 'mean:', sigmean.getValV(), '+-', sigmean.getError()
print 'width:', sigwidth.getValV(), '+-', sigwidth.getError()
print 'argpar:', argpar.getValV(), '+-', argpar.getError()


from time import sleep
sleep(10)

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
        rrv_number_Total_background_MC.setError(TMath.Sqrt(
                rrv_number_WJets.getError()* rrv_number_WJets.getError()+
                rrv_number_VV.getError()* rrv_number_VV.getError()+
                rrv_number_TTbar.getError()* rrv_number_TTbar.getError()+
                rrv_number_STop.getError() *rrv_number_STop.getError()
                ));

        model_Total_background_MC = RooAddPdf("model_Total_background_MC_xww","model_Total_background_MC_xww",RooArgList(old_workspace.pdf("WJets_xww_%s_%s"%(options.channel,options.category)), old_workspace.pdf("VV_xww_%s_%s"%(options.channel,options.category)),old_workspace.pdf("TTbar_xww_%s_%s"%(options.channel,options.category)),old_workspace.pdf("STop_xww_%s_%s"%(options.channel,options.category))),RooArgList(rrv_number_WJets,rrv_number_VV,rrv_number_TTbar,rrv_number_STop));

        rrv_number_signal.setVal(rrv_number_signal.getVal()*6.25);

        #### scale factor in order to scale MC to data in the final plot -> in order to avoid the normalization to data which is done by default in rooFit
        scale_number_Total_background_MC = rrv_number_Total_background_MC.getVal()/old_workspace.data(datasetname+"_xww_"+options.channel+"_"+options.category).sumEntries();
        scale_number_signal = rrv_number_signal.getVal()/old_workspace.data(datasetname+"_xww_"+options.channel+"_"+options.category).sumEntries();

        model_Total_background_MC.plotOn(mplot,RooFit.Normalization(scale_number_Total_background_MC),RooFit.Name("total_MC"),RooFit.Components("WJets_xww_%s_%s,VV_xww_%s_%s,TTbar_xww_%s_%s,STop_xww_%s_%s"%(options.channel,options.category,options.channel,options.category,options.channel,options.category,options.channel,options.category)),RooFit.DrawOption("L"), RooFit.LineColor(kRed), RooFit.VLines(),RooFit.LineWidth(2));

            
        model_signal_background_MC = RooAddPdf("model_signal_background_MC_xww","model_signal_background_MC_xww",RooArgList(model_pdf,model_Total_background_MC),RooArgList(rrv_number_signal,rrv_number_Total_background_MC));

        model_signal_background_MC.plotOn(mplot,RooFit.Normalization(scale_number_Total_background_MC+scale_number_signal),RooFit.Name("total_SpB_MC"),RooFit.Components("BulkWW_xww_%s_%s,model_Total_background_MC_xww"%(options.channel,options.category)),RooFit.DrawOption("L"), RooFit.LineColor(kBlue), RooFit.VLines(),RooFit.LineWidth(2),RooFit.LineStyle(7));

        model_pdf.plotOn(mplot,RooFit.Name("total_S_MC"),RooFit.Normalization(scale_number_signal),RooFit.DrawOption("L"), RooFit.LineColor(kGreen+2), RooFit.VLines(),RooFit.LineWidth(2),RooFit.LineStyle(kDashed));

        os.system("mkdir -p ./plots_signal_width");
        name = TString("check_%.3f_%.3f"%(mass[iMass],gammaVal));
        name.ReplaceAll("0.","0_");
        mplot.GetYaxis().SetRangeUser(1e-3,mplot.GetMaximum()*1.2);
        draw_canvas(mplot,"plots_signal_width/",name,0,1);
        
        new_workspace.writeToFile(new_file.GetName());

# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import plotOn [as 别名]
gauss5 = RooGaussian("gauss5","gaussian PDF",x,mean5,sigma5)


# add PDF gauss_combine1 = gauss4 + gauss5

frac_combine1   = RooRealVar("frac_combine1",   "fraction of gauss4 wrt gauss5", 0.7, 0.,   1.)

pdf_combine1 = RooAddPdf("pdf_combine1"," gauss4 + gauss5 ", RooArgList(gauss4 , gauss5 ), RooArgList( frac_combine1 ))


#pdf_combine1.plotOn(xframe3,RooFit.LineColor(RooFit.kBlue))

gauss4.plotOn(xframe3, RooFit.Normalization( frac_combine1.getVal()   ,RooAbsReal.Relative),RooFit.LineColor(RooFit.kOrange))
gauss5.plotOn(xframe3, RooFit.Normalization( 1-frac_combine1.getVal() ,RooAbsReal.Relative),RooFit.LineColor(RooFit.kCyan))
pdf_combine1.plotOn(xframe3, RooFit.Normalization(1.0,RooAbsReal.Relative) ,RooFit.LineColor(RooFit.kBlue))

# -------------------------------------------
# 4. use combine PDF to generate MC

xframe4 = x.frame(RooFit.Title("4. use combine PDF to generate MC"))

data2 = pdf_combine1.generate(RooArgSet(x),500)

#pdf_combine1.plotOn(xframe4, RooFit.Normalization(500, RooAbsReal.NumEvent) , RooFit.LineColor(RooFit.kOrange))
data2.plotOn(xframe4)
pdf_combine1.plotOn(xframe4,RooFit.LineColor(RooFit.kBlue))

# -------------------------------------------
# 5. use combine PDF to fit the toy MC