Python ROOT.RooAddPdf类代码示例

def test_correlated_values():

    try:
        import uncertainties
    except ImportError:
        raise SkipTest("uncertainties package is not installed")
    from rootpy.stats.correlated_values import correlated_values

    # construct pdf and toy data following example at
    # http://root.cern.ch/drupal/content/roofit

    # --- Observable ---
    mes = RooRealVar("mes", "m_{ES} (GeV)", 5.20, 5.30)

    # --- Parameters ---
    sigmean = RooRealVar("sigmean", "B^{#pm} mass", 5.28, 5.20, 5.30)
    sigwidth = RooRealVar("sigwidth", "B^{#pm} width", 0.0027, 0.001, 1.)

    # --- Build Gaussian PDF ---
    signal = RooGaussian("signal", "signal PDF", mes, sigmean, sigwidth)

    # --- Build Argus background PDF ---
    argpar = RooRealVar("argpar", "argus shape parameter", -20.0, -100., -1.)
    background = RooArgusBG("background", "Argus PDF",
                            mes, RooFit.RooConst(5.291), argpar)

    # --- 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 ---
    fitresult = model.fitTo(data, RooFit.Save(), RooFit.PrintLevel(-1))

    nsig, nbkg = correlated_values(["nsig", "nbkg"], fitresult)

    # Arbitrary math expression according to what the `uncertainties`
    # package supports, automatically computes correct error propagation
    sum_value = nsig + nbkg
    value, error = sum_value.nominal_value, sum_value.std_dev

    workspace = Workspace(name='workspace')
    # import the data
    assert_false(workspace(data))
    with TemporaryFile():
        workspace.Write()

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

def get_num_sig_bkg(hist_DataTemplate,
                    hist_SignalTemplate,
                    hist_BackgdTemplate,
                    fit_range_min,
                    fit_range_max):
    '''Given 3 input histograms (TH1F), and a fit range, this function finds
    the amount of signal and background that sum up to the data histogram.
    It does histogram fits.'''
    # Find range of data template
    data_min = hist_DataTemplate.GetXaxis().GetXmin()
    data_max = hist_DataTemplate.GetXaxis().GetXmax()
    
    # Create basic variables
    x = RooRealVar("x","x",data_min,data_max)
    x.setBins(hist_DataTemplate.GetXaxis().GetNbins())  # Binned x values
    nsig = RooRealVar("nsig","number of signal events"    , 0, hist_DataTemplate.Integral())
    nbkg = RooRealVar("nbkg","number of background events", 0, hist_DataTemplate.Integral())
    
    # Create RooDataHists from input TH1Fs
    dh = RooDataHist("dh","dh",RooArgList(x),hist_DataTemplate)
    ds = RooDataHist("ds","ds",RooArgList(x),hist_SignalTemplate)
    db = RooDataHist("db","db",RooArgList(x),hist_BackgdTemplate)
    
    # Create Probability Distribution Functions from Monte Carlo
    sigPDF = RooHistPdf("sigPDF", "sigPDF", RooArgSet(x), ds)
    bkgPDF = RooHistPdf("bkgPDF", "bkgPDF", RooArgSet(x), db)
    
    model = RooAddPdf("model","(g1+g2)+a",RooArgList(bkgPDF,sigPDF),RooArgList(nbkg,nsig))
    
    # Find the edges of the bins that contain the fit range min/max
    data_min = hist_DataTemplate.GetXaxis().GetBinLowEdge(hist_DataTemplate.GetXaxis().FindFixBin(fit_range_min))
    data_max = hist_DataTemplate.GetXaxis().GetBinUpEdge(hist_DataTemplate.GetXaxis().FindFixBin(fit_range_max))
    
    r = model.fitTo(dh,RooFit.Save(),RooFit.Minos(0),RooFit.PrintEvalErrors(0),
                    RooFit.Extended(),RooFit.Range(data_min,data_max))
    r.Print("v")

    #print nsig.getVal(), nsig.getError(), nbkg.getVal(), nbkg.getError()
    #  Create pull distribution
    #mcstudy = RooMCStudy(model, RooArgSet(x), RooFit.Binned(1), RooFit.Silence(),
    #                     RooFit.Extended(),
    #                     RooFit.FitOptions(RooFit.Save(1),
    #                                       RooFit.PrintEvalErrors(0),
    #                                       RooFit.Minos(0))
    #                    )
    #mcstudy.generateAndFit(500)                          # Generate and fit toy MC
    #pull_dist = mcstudy.plotPull(nsig, -3.0, 3.0, 30, 1)  # make pull distribution
    pull_dist = None
    return [nsig.getVal(), nsig.getError(), nbkg.getVal(), nbkg.getError(), pull_dist]

#.........这里部分代码省略.........
            c.SetWindowSize(1000,500)
            c.Divide(2,1)
            for i in [1,2] :
                c.cd(i)
                reg='barrel'
                if i==2: reg='endcap' 

                h=histos[r+k+reg]
                x=RooRealVar("x", h.GetXaxis().GetTitle(), h.GetXaxis().GetXmin(), h.GetXaxis().GetXmax())
                data=RooDataHist("data", "dataset with x", RooArgList(x), h)
                frame=x.frame()
                RooAbsData.plotOn( data, frame, RooFit.DataError(RooAbsData.SumW2) )

                mean1=RooRealVar("mean1","mean1",0,-0.5,0.5);
                sigma1=RooRealVar("sigma1","sigma1",0.1,0.01,1.0);
                gauss1=RooGaussian("g1","g",x,mean1,sigma1)
                
                if r=='dpt' or r=='den' :
                    mean2=RooRealVar("mean2","mean2",0,-0.5,0.5);
                    sigma2=RooRealVar("sigma2","sigma2",0.1,0.01,1.0);
                    alphacb=RooRealVar("alphacb","alphacb",1,0.1,3);
                    ncb=RooRealVar("ncb","ncb",4,1,100)
                    gauss2 = RooCBShape("cb2","cb",x,mean2,sigma2,alphacb,ncb);
                else:
                    mean1.setRange(0,0.5)
                    mean2=RooRealVar("mean2","mean",0,0,1);
                    sigma2=RooRealVar("sigma2","sigma",0.1,0.01,1.0);
                    gauss2=RooGaussian("g2","g",x,mean2,sigma2) ;

                frac = RooRealVar("frac","fraction",0.9,0.0,1.0)
                if data.sumEntries()<100 :
                    frac.setVal(1.0)
                    frac.setConstant(True)
                model = RooAddPdf("sum","g1+g2",RooArgList(gauss1,gauss2), RooArgList(frac))

                status=model.fitTo(data,RooFit.Save()).status()
                if status!=0 : continue

                model_cdf=model.createCdf(RooArgSet(x)) ;
                cl=0.90
                ul=0.5*(1.0+cl)
                closestToCL=1.0
                closestToUL=-1
                closestToMedianCL=1.0
                closestToMedian=-1
                for ibin in xrange(1,h.GetXaxis().GetNbins()*10):
                    xval=h.GetXaxis().GetXmin()+(ibin-1)*h.GetXaxis().GetBinWidth(ibin)/10.
                    x.setVal(xval)
                    cdfValToCL=math.fabs(model_cdf.getVal()-ul)
                    if cdfValToCL<closestToCL:
                        closestToCL=cdfValToCL
                        closestToUL=xval
                    cdfValToCL=math.fabs(model_cdf.getVal()-0.5)
                    if cdfValToCL<closestToMedianCL:
                        closestToMedianCL=cdfValToCL
                        closestToMedian=xval

                RooAbsPdf.plotOn(model,frame)
                frame.Draw()

                if i==1: drawHeader()
                labels.append( TPaveText(0.6,0.92,0.9,0.98,'brNDC') )
                ilab=len(labels)-1
                labels[ilab].SetName(r+k+'txt')
                labels[ilab].SetBorderSize(0)
                labels[ilab].SetFillStyle(0)

	def fit(self, save_to, signal_name=None, fix_p3=False, fit_range=[300., 1200.], fit_strategy=1):
		# Run a RooFit fit

		# Create background PDF
		p1 = RooRealVar('p1','p1',args.p1,0.,100.)
		p2 = RooRealVar('p2','p2',args.p2,0.,60.)
		p3 = RooRealVar('p3','p3',args.p3,-10.,10.)
		if args.fix_p3:
			p3.setConstant()
		background_pdf = RooGenericPdf('background_pdf','(pow([email protected]/%.1f,@1)/pow(@0/%.1f,@[email protected]*log(@0/%.1f)))'%(self.collision_energy,self.collision_energy,self.collision_energy),RooArgList(self.mjj_,p1,p2,p3))
		background_pdf.Print()
		data_integral = data_histogram.Integral(data_histogram.GetXaxis().FindBin(float(fit_range[0])),data_histogram.GetXaxis().FindBin(float(fit_range[1])))
		background_norm = RooRealVar('background_norm','background_norm',data_integral,0.,1e+08)
		background_norm.Print()

		# Create signal PDF and fit model
		if signal_name:
			signal_pdf = RooHistPdf('signal_pdf', 'signal_pdf', RooArgSet(self.mjj_), self.signal_roohistograms_[signal_name])
			signal_pdf.Print()
			signal_norm = RooRealVar('signal_norm','signal_norm',0,-1e+05,1e+05)
			signal_norm.Print()
			model = RooAddPdf("model","s+b",RooArgList(background_pdf,signal_pdf),RooArgList(background_norm,signal_norm))
		else:
			model = RooAddPdf("model","b",RooArgList(background_pdf),RooArgList(background_norm))

		# Run fit
		res = model.fitTo(data_, RooFit.Save(kTRUE), RooFit.Strategy(fit_strategy))

		# Save to workspace
		self.workspace_ = RooWorkspace('w','workspace')
		#getattr(w,'import')(background,ROOT.RooCmdArg())
		getattr(self.workspace_,'import')(background_pdf,RooFit.Rename("background"))
		getattr(self.workspace_,'import')(background_norm,ROOT.RooCmdArg())
		getattr(self.workspace_,'import')(self.data_roohistogram_,RooFit.Rename("data_obs"))
		getattr(self.workspace_, 'import')(model, RooFit.Rename("model"))
		if signal_name:
			getattr(self.workspace_,'import')(signal_roohistogram,RooFit.Rename("signal"))
			getattr(self.workspace_,'import')(signal_pdf,RooFit.Rename("signal_pdf"))
			getattr(self.workspace_,'import')(signal_norm,ROOT.RooCmdArg())
	
		self.workspace_.Print()
		self.workspace_.writeToFile(save_to)
		if signal_name:
			roofit_results[signal_name] = save_to
		else:
			roofit_results["background"] = save_to

if fitSig: 

    # define parameters for signal fit
    m = RooRealVar('mean','mean',float(mass),float(mass)-200,float(mass)+200)
    s = RooRealVar('sigma','sigma',0.1*float(mass),0,10000)
    a = RooRealVar('alpha','alpha',1,-10,10)
    n = RooRealVar('n','n',1,0,100)
    sig = RooCBShape('sig','sig',x,m,s,a,n)        

    p  = RooRealVar('p','p',1,0,5)
    x0 = RooRealVar('x0','x0',1000,100,5000)

    bkg = RooGenericPdf('bkg','1/(exp(pow(@0/@1,@2))+1)',RooArgList(x,x0,p))

    fsig= RooRealVar('fsig','fsig',0.5,0.,1.)
    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)

#.........这里部分代码省略.........
    nLSBVjet = iLSBVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(LSBcut)
    nHSBVjet = iHSBVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(HSBcut)
    nSRVjet = iSRVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(SRcut)
    
    drawPlot("JetMass_Vjet", channel, J_mass, modelVjet, setVjet, binsJmass, frVjet)

    if VERBOSE: print "********** Fit result [JET MASS Vjets] *"+"*"*40, "\n", frVjet.Print(), "\n", "*"*80
    
    #*******************************************************#
    #                                                       #
    #                 VV, VH normalization                  #
    #                                                       #
    #*******************************************************#
    
    # Variables for VV
    # Error function and exponential to model the bulk
    constVV  = RooRealVar("constVV",  "slope of the exp",  -0.030, -0.1,   0.)
    offsetVV = RooRealVar("offsetVV", "offset of the erf", 90.,     1., 300.)
    widthVV  = RooRealVar("widthVV",  "width of the erf",  50.,     1., 100.)
    erfrVV   = RooErfExpPdf("baseVV", "error function for VV jet mass", J_mass, constVV, offsetVV, widthVV)
    expoVV   = RooExponential("baseVV", "error function for VV jet mass", J_mass, constVV)
    # gaussian for the V mass peak
    meanVV   = RooRealVar("meanVV",   "mean of the gaussian",           90.,    60., 100.)
    sigmaVV  = RooRealVar("sigmaVV",  "sigma of the gaussian",          10.,     6.,  30.)
    fracVV   = RooRealVar("fracVV",   "fraction of gaussian wrt erfexp", 3.2e-1, 0.,   1.)
    gausVV   = RooGaussian("gausVV",  "gaus for VV jet mass", J_mass, meanVV, sigmaVV)
    # gaussian for the H mass peak
    meanVH   = RooRealVar("meanVH",   "mean of the gaussian",           125.,   100., 150.)
    sigmaVH  = RooRealVar("sigmaVH",  "sigma of the gaussian",           30.,     5.,  40.)
    fracVH   = RooRealVar("fracVH",   "fraction of gaussian wrt erfexp",  1.5e-2, 0.,   1.)
    gausVH   = RooGaussian("gausVH",  "gaus for VH jet mass", J_mass, meanVH, sigmaVH)
    
    # Define VV model
    if fitFuncVV == "ERFEXPGAUS": modelVV  = RooAddPdf("modelVV",   "error function + gaus for VV jet mass", RooArgList(gausVV, erfrVV), RooArgList(fracVV))
    elif fitFuncVV == "ERFEXPGAUS2": modelVV  = RooAddPdf("modelVV",   "error function + gaus + gaus for VV jet mass", RooArgList(gausVH, gausVV, erfrVV), RooArgList(fracVH, fracVV))
    elif fitFuncVV == "EXPGAUS": modelVV  = RooAddPdf("modelVV",   "error function + gaus for VV jet mass", RooArgList(gausVV, expoVV), RooArgList(fracVV))
    elif fitFuncVV == "EXPGAUS2": modelVV  = RooAddPdf("modelVV",   "error function + gaus + gaus for VV jet mass", RooArgList(gausVH, gausVV, expoVV), RooArgList(fracVH, fracVV))
    else:
        print "  ERROR! Pdf", fitFuncVV, "is not implemented for VV"
        exit()
    
    # fit to secondary bkg in MC (whole range)
    frVV = modelVV.fitTo(setVV, RooFit.SumW2Error(True), RooFit.Range("h_reasonable_range"), RooFit.Strategy(2), RooFit.Minimizer("Minuit2"), RooFit.Save(1), RooFit.PrintLevel(1 if VERBOSE else -1))
    
    # integrals and number of events
    iSBVV = modelVV.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange,HSBrange"))
    iLSBVV = modelVV.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange"))
    iHSBVV = modelVV.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("HSBrange"))
    iSRVV = modelVV.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("SRrange"))
    iVRVV = modelVV.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("VRrange"))
    # Do not remove the following lines, integrals are computed here
    iALVV = modelVV.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg))
    nSBVV = iSBVV.getVal()/iALVV.getVal()*setVV.sumEntries(SBcut)
    nLSBVV = iLSBVV.getVal()/iALVV.getVal()*setVV.sumEntries(LSBcut)
    nHSBVV = iHSBVV.getVal()/iALVV.getVal()*setVV.sumEntries(HSBcut)
    nSRVV = iSRVV.getVal()/iALVV.getVal()*setVV.sumEntries(SRcut)
    rSBSRVV = nSRVV/nSBVV
    
    drawPlot("JetMass_VV", channel, J_mass, modelVV, setVV, binsJmass, frVV)
    
    if VERBOSE: print "********** Fit result [JET MASS VV] ****"+"*"*40, "\n", frVV.Print(), "\n", "*"*80
    
    #*******************************************************#
    #                                                       #
    #                 Top, ST normalization                 #
    #                                                       #

def findOnePe(hist, ws, name='x', Npe = 1):
    fitPed(hist, ws, name)
    x = ws.var(name)

    ped = ws.pdf('ped')
    pedWidth = ws.var('pedWidth')

    pdfs = RooArgList(ped)
    pdfList = []

    fped = RooRealVar('fped', 'f_{ped}', 0.8, 0., 1.)
    fractions = RooArgList(fped)
    fList = []
    peList = []

    peMean = RooRealVar('peMean', 'mean_{pe}', 6., 0., 20.)
    peWidth = RooRealVar('peWidth', 'width_{pe}', pedWidth.getVal(), 0., 10.)

    for i in range(0, Npe):
        pem = RooFormulaVar('pem{0}'.format(i+1), '@0+{0}*@1'.format(i+1),
                            RooArgList(ws.var('pedMean'), peMean))
        peList.append(pem)
        npepdf = RooGaussian('pe{0}pdf'.format(i+1), 'pe{0}pdf'.format(i+1),
                             x, pem, pedWidth)
        pdfs.add(npepdf)
        pdfList.append(npepdf)

        fnpe = RooRealVar('fpe{0}'.format(i+1), 'fpe{0}'.format(i+1),
                          0.5, -0.1, 1.0)
        fractions.add(fnpe)
        fList.append(fnpe)

    #bgMean = RooRealVar("bgMean", "bgMean", 6.0, x.getMin(), x.getMax())
    bgScale = RooRealVar("bgScale", "bgScale", 0.5, -1.0, Npe + 1.0)
    bgMean = RooFormulaVar("bgMean", "@[email protected]*@2",
                           RooArgList(peMean, ws.var('pedMean'), bgScale))
    bgWidthL = RooRealVar("bgWidthL", "bgWidthL", pedWidth.getVal()*2,
                          0., 25.)
    bgWidthR = RooRealVar("bgWidthR", "bgWidthR", pedWidth.getVal()*7,
                          0., 25.)

    bgGauss = RooBifurGauss("bgGauss", "bgGauss", x, bgMean,
                            bgWidthR, bgWidthR)

    if (Npe > 1):
        pdfs.add(bgGauss)
    else:
        fractions.remove(fractions.at(fractions.getSize()-1))

##     pem = RooFormulaVar('pem', '@[email protected]', RooArgList(peMean, ws.var('pedMean')))
##     firstPe = RooGaussian('firstPe', 'firstPe', x, pem, peWidth)

##     pdfs.Print("v")
##     fractions.Print("v")
    pedPlusOne = RooAddPdf('pedPlusOne', 'pedPlusOne', pdfs, fractions, True)

    ## pedWidth = ped.GetParameter(2)
    ## pedMean = ped.GetParameter(1)
    ## pedA = ped.GetParameter(0)
    
    secondMax = hist.GetMaximumBin() + 1
    goingDown = True
    maxVal = hist.GetBinContent(secondMax)
    foundMax = False
    while (not foundMax) and (secondMax < hist.GetNbinsX()):
        tmpVal = hist.GetBinContent(secondMax+1)
        if (tmpVal < maxVal):
            if not goingDown:
                foundMax = True
            else:
                goingDown = True
                maxVal = tmpVal
                secondMax += 1
        elif (tmpVal > maxVal):
            goingDown = False
            maxVal = tmpVal
            secondMax += 1
        else:
            maxVal = tmpVal
            secondMax += 1

    secondMaxx = hist.GetBinCenter(secondMax)
    print 'found 2nd maximum in bin',secondMax,'value',secondMaxx

##     peMean.setVal(secondMaxx)
##     bgMean.setVal(secondMaxx*0.6)
    x.setRange('pedPlus_fit', x.getMin(), ws.var('pedMean').getVal()+pedWidth.getVal()*6.*(Npe+0))

    pedPlusOne.fitTo(ws.data('ds'), RooFit.Minos(False),
                     RooFit.Range('pedPlus_fit'),
                     RooFit.PrintLevel(1))

    getattr(ws, 'import')(pedPlusOne)

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

def rf501_simultaneouspdf():
    # C r e a t e   m o d e l   f o r   p h y s i c s   s a m p l e
    # -------------------------------------------------------------

    # Create observables
    x = RooRealVar( "x", "x", -8, 8 ) 

    # Construct signal pdf
    mean = RooRealVar( "mean", "mean", 0, -8, 8 ) 
    sigma = RooRealVar( "sigma", "sigma", 0.3, 0.1, 10 ) 
    gx = RooGaussian( "gx", "gx", x, mean, sigma ) 

    # Construct background pdf
    a0 = RooRealVar( "a0", "a0", -0.1, -1, 1 ) 
    a1 = RooRealVar( "a1", "a1", 0.004, -1, 1 ) 
    px = RooChebychev( "px", "px", x, RooArgList( a0, a1 ) ) 

    # Construct composite pdf
    f = RooRealVar( "f", "f", 0.2, 0., 1. ) 
    model = RooAddPdf( "model", "model", RooArgList( gx, px ), RooArgList( f ) ) 



    # C r e a t e   m o d e l   f o r   c o n t r o l   s a m p l e
    # --------------------------------------------------------------

    # Construct signal pdf. 
    # NOTE that sigma is shared with the signal sample model
    mean_ctl = RooRealVar( "mean_ctl", "mean_ctl", -3, -8, 8 ) 
    gx_ctl = RooGaussian( "gx_ctl", "gx_ctl", x, mean_ctl, sigma ) 

    # Construct the background pdf
    a0_ctl = RooRealVar( "a0_ctl", "a0_ctl", -0.1, -1, 1 ) 
    a1_ctl = RooRealVar( "a1_ctl", "a1_ctl", 0.5, -0.1, 1 ) 
    px_ctl = RooChebychev( "px_ctl", "px_ctl", x, RooArgList( a0_ctl, a1_ctl ) ) 

    # Construct the composite model
    f_ctl = RooRealVar( "f_ctl", "f_ctl", 0.5, 0., 1. ) 
    model_ctl = RooAddPdf( "model_ctl", "model_ctl", RooArgList( gx_ctl, px_ctl ),
                           RooArgList( f_ctl ) ) 
    


    # G e n e r a t e   e v e n t s   f o r   b o t h   s a m p l e s 
    # ---------------------------------------------------------------

    # Generate 1000 events in x and y from model
    data = model.generate( RooArgSet( x ), 100 ) 
    data_ctl = model_ctl.generate( RooArgSet( x ), 2000 ) 



    # C r e a t e   i n d e x   c a t e g o r y   a n d   j o i n   s a m p l e s 
    # ---------------------------------------------------------------------------
    # Define category to distinguish physics and control samples events
    sample = RooCategory( "sample", "sample" ) 
    sample.defineType( "physics" ) 
    sample.defineType( "control" ) 

    # Construct combined dataset in (x,sample)
    combData = RooDataSet( "combData", "combined data", RooArgSet(x), RooFit.Index( sample ),
                          RooFit.Import( "physics", data ),
                          RooFit.Import( "control", data_ctl ) ) 



    # C o n s t r u c t   a   s i m u l t a n e o u s   p d f   i n   ( x , s a m p l e )
    # -----------------------------------------------------------------------------------

    # Construct a simultaneous pdf using category sample as index
    simPdf = RooSimultaneous( "simPdf", "simultaneous pdf", sample ) 

    # Associate model with the physics state and model_ctl with the control state
    simPdf.addPdf( model, "physics" ) 
    simPdf.addPdf( model_ctl, "control" ) 



    # P e r f o r m   a   s i m u l t a n e o u s   f i t
    # ---------------------------------------------------

    # Perform simultaneous fit of model to data and model_ctl to data_ctl
    simPdf.fitTo( combData ) 



    # P l o t   m o d e l   s l i c e s   o n   d a t a    s l i c e s 
    # ----------------------------------------------------------------

    # Make a frame for the physics sample
    frame1 = x.frame( RooFit.Bins( 30 ), RooFit.Title( "Physics sample" ) ) 

    # Plot all data tagged as physics sample
    combData.plotOn( frame1, RooFit.Cut( "sample==sample::physics" ) ) 

    # Plot "physics" slice of simultaneous pdf. 
    # NBL You _must_ project the sample index category with data using ProjWData 
    # as a RooSimultaneous makes no prediction on the shape in the index category 
    # and can thus not be integrated
    simPdf.plotOn( frame1, RooFit.Slice( sample, "physics" ),
#.........这里部分代码省略.........

def rf501_simultaneouspdf():
    signal_1, bkg_1, signal_2, bkg_2 = get_templates()
    # C r e a t e   m o d e l   f o r   p h y s i c s   s a m p l e
    # -------------------------------------------------------------

    # Create observables
    x = RooRealVar( "x", "x", 0, 200 ) 
    x.setBins(n_bins)
    nsig = RooRealVar( "nsig", "#signal events", N_signal_obs, 0., 2*N_data )
    nbkg = RooRealVar( "nbkg", "#background events", N_bkg1_obs, 0., 2*N_data )

    # Construct signal pdf
#     mean = RooRealVar( "mean", "mean", mu4, 40, 200 ) 
#     sigma = RooRealVar( "sigma", "sigma", sigma4, 0.1, 20 )
#     gx = RooGaussian( "gx", "gx", x, mean, sigma ) 
    roofit_signal_1 = RooDataHist( 'signal_1', 'signal_1', RooArgList(x), signal_1 )
    signal_1_pdf = RooHistPdf ( 'signal_1_pdf' , 'signal_1_pdf', RooArgSet(x), roofit_signal_1) 

    # Construct background pdf
#     mean_bkg = RooRealVar( "mean_bkg", "mean_bkg", mu3, 40, 200 ) 
#     sigma_bkg = RooRealVar( "sigma_bkg", "sigma_bkg", sigma3, 0.1, 20 ) 
#     px = RooGaussian( "px", "px", x, mean_bkg, sigma_bkg ) 
    roofit_bkg_1 = RooDataHist( 'bkg_1', 'bkg_1', RooArgList(x), bkg_1 )
    bkg_1_pdf = RooHistPdf ( 'bkg_1_pdf' , 'bkg_1_pdf', RooArgSet(x), roofit_bkg_1) 

    # Construct composite pdf
    model = RooAddPdf( "model", "model", RooArgList( signal_1_pdf, bkg_1_pdf ), RooArgList( nsig, nbkg ) ) 



    # C r e a t e   m o d e l   f o r   c o n t r o l   s a m p l e
    # --------------------------------------------------------------

    # Construct signal pdf. 
    # NOTE that sigma is shared with the signal sample model
    y = RooRealVar( "y", "y", 0, 200 )
    y.setBins(n_bins)
    mean_ctl = RooRealVar( "mean_ctl", "mean_ctl", mu2, 0, 200 ) 
    sigma_ctl = RooRealVar( "sigma", "sigma", sigma2, 0.1, 10 ) 
    gx_ctl = RooGaussian( "gx_ctl", "gx_ctl", y, mean_ctl, sigma_ctl ) 

    # Construct the background pdf
    mean_bkg_ctl = RooRealVar( "mean_bkg_ctl", "mean_bkg_ctl", mu1, 0, 200 ) 
    sigma_bkg_ctl = RooRealVar( "sigma_bkg_ctl", "sigma_bkg_ctl", sigma1, 0.1, 20 ) 
    px_ctl = RooGaussian( "px_ctl", "px_ctl", y, mean_bkg_ctl, sigma_bkg_ctl ) 

    # Construct the composite model
#     f_ctl = RooRealVar( "f_ctl", "f_ctl", 0.5, 0., 20. ) 
    model_ctl = RooAddPdf( "model_ctl", "model_ctl", RooArgList( gx_ctl, px_ctl ),
                           RooArgList( nsig, nbkg ) ) 
    


    # G e t   e v e n t s   f o r   b o t h   s a m p l e s 
    # ---------------------------------------------------------------
    real_data, real_data_ctl = get_data()
    real_data_hist = RooDataHist( 'real_data_hist',
                                 'real_data_hist',
                                 RooArgList( x ),
                                 real_data )
    real_data_ctl_hist = RooDataHist( 'real_data_ctl_hist',
                                     'real_data_ctl_hist',
                                     RooArgList( y ),
                                     real_data_ctl )
    input_hists = MapStrRootPtr()
    input_hists.insert( StrHist( "physics", real_data ) )
    input_hists.insert( StrHist( "control", real_data_ctl ) )

    # C r e a t e   i n d e x   c a t e g o r y   a n d   j o i n   s a m p l e s 
    # ---------------------------------------------------------------------------
    # Define category to distinguish physics and control samples events
    sample = RooCategory( "sample", "sample" ) 
    sample.defineType( "physics" ) 
    sample.defineType( "control" ) 

    # Construct combined dataset in (x,sample)
    combData = RooDataHist( "combData", "combined data", RooArgList( x), sample ,
                           input_hists )


    # C o n s t r u c t   a   s i m u l t a n e o u s   p d f   i n   ( x , s a m p l e )
    # -----------------------------------------------------------------------------------

    # Construct a simultaneous pdf using category sample as index
    simPdf = RooSimultaneous( "simPdf", "simultaneous pdf", sample ) 

    # Associate model with the physics state and model_ctl with the control state
    simPdf.addPdf( model, "physics" ) 
    simPdf.addPdf( model_ctl, "control" ) 

#60093.048127    173.205689173    44.7112503776

    # P e r f o r m   a   s i m u l t a n e o u s   f i t
    # ---------------------------------------------------
    model.fitTo( real_data_hist,
                RooFit.Minimizer( "Minuit2", "Migrad" ),
                        RooFit.NumCPU( 1 ),
#                         RooFit.Extended(),
#                         RooFit.Save(), 
                        )
#.........这里部分代码省略.........

h.Rebin(20)
x = RooRealVar('mjj','mjj',900,4500)

m = RooRealVar('mean','mean',float(mass),float(mass)-200,float(mass)+200)
s = RooRealVar('sigma','sigma',0.1*float(mass),0,10000)
a = RooRealVar('alpha','alpha',1,-10,10)
n = RooRealVar('n','n',1,0,100)
sig = RooCBShape('sig','sig',x,m,s,a,n)        

p  = RooRealVar('p','p',1,0,5)
x0 = RooRealVar('x0','x0',1000,100,5000)

bkg = RooGenericPdf('bkg','1/(exp(pow(@0/@1,@2))+1)',RooArgList(x,x0,p))

fsig= RooRealVar('fsig','fsig',0.5,0.,1.)
model = RooAddPdf('model','model',sig,bkg,fsig)

can = TCanvas('can_Mjj'+str(mass),'can_Mjj'+str(mass),900,600)
h.Draw()
gPad.SetLogy() 

roohist = RooDataHist('roohist','roohist',RooArgList(x),h)


model.fitTo(roohist)
frame = x.frame()
roohist.plotOn(frame)
model.plotOn(frame)
model.plotOn(frame,RooFit.Components('bkg'),RooFit.LineColor(ROOT.kRed),RooFit.LineWidth(2),RooFit.LineStyle(ROOT.kDashed))
frame.Draw('same')

#.........这里部分代码省略.........
    )

    if VERBOSE:
        print "********** Fit result [JET MASS Vjets] *" + "*" * 40, "\n", frVjet.Print(), "\n", "*" * 80

    # likelihoodScan(VjetMass, setVjet, [constVjet, offsetVjet, widthVjet])

    # *******************************************************#
    #                                                       #
    #                 VV, VH normalization                  #
    #                                                       #
    # *******************************************************#

    # Variables for VV
    # Error function and exponential to model the bulk
    constVV = RooRealVar("constVV", "slope of the exp", -0.030, -0.1, 0.0)
    offsetVV = RooRealVar("offsetVV", "offset of the erf", 90.0, 1.0, 300.0)
    widthVV = RooRealVar("widthVV", "width of the erf", 50.0, 1.0, 100.0)
    erfrVV = RooErfExpPdf("baseVV", "error function for VV jet mass", J_mass, constVV, offsetVV, widthVV)
    expoVV = RooExponential("baseVV", "error function for VV jet mass", J_mass, constVV)
    # gaussian for the V mass peak
    meanVV = RooRealVar("meanVV", "mean of the gaussian", 90.0, 60.0, 100.0)
    sigmaVV = RooRealVar("sigmaVV", "sigma of the gaussian", 10.0, 6.0, 30.0)
    fracVV = RooRealVar("fracVV", "fraction of gaussian wrt erfexp", 3.2e-1, 0.0, 1.0)
    gausVV = RooGaussian("gausVV", "gaus for VV jet mass", J_mass, meanVV, sigmaVV)
    # gaussian for the H mass peak
    meanVH = RooRealVar("meanVH", "mean of the gaussian", 125.0, 100.0, 150.0)
    sigmaVH = RooRealVar("sigmaVH", "sigma of the gaussian", 10.0, 5.0, 50.0)
    fracVH = RooRealVar("fracVH", "fraction of gaussian wrt erfexp", 1.5e-2, 0.0, 1.0)
    gausVH = RooGaussian("gausVH", "gaus for VH jet mass", J_mass, meanVH, sigmaVH)

    # Define VV model
    if fitFuncVV == "ERFEXPGAUS":
        VVMass = RooAddPdf("VVMass", fitFuncVV, RooArgList(gausVV, erfrVV), RooArgList(fracVV))
    elif fitFuncVV == "ERFEXPGAUS2":
        VVMass = RooAddPdf("VVMass", fitFuncVV, RooArgList(gausVH, gausVV, erfrVV), RooArgList(fracVH, fracVV))
    elif fitFuncVV == "EXPGAUS":
        VVMass = RooAddPdf("VVMass", fitFuncVV, RooArgList(gausVV, expoVV), RooArgList(fracVV))
    elif fitFuncVV == "EXPGAUS2":
        VVMass = RooAddPdf("VVMass", fitFuncVV, RooArgList(gausVH, gausVV, expoVV), RooArgList(fracVH, fracVV))
    else:
        print "  ERROR! Pdf", fitFuncVV, "is not implemented for VV"
        exit()

    # fit to secondary bkg in MC (whole range)
    frVV = VVMass.fitTo(
        setVV,
        RooFit.SumW2Error(True),
        RooFit.Range("h_reasonable_range"),
        RooFit.Strategy(2),
        RooFit.Minimizer("Minuit2"),
        RooFit.Save(1),
        RooFit.PrintLevel(1 if VERBOSE else -1),
    )

    if VERBOSE:
        print "********** Fit result [JET MASS VV] ****" + "*" * 40, "\n", frVV.Print(), "\n", "*" * 80

    # *******************************************************#
    #                                                       #
    #                 Top, ST normalization                 #
    #                                                       #
    # *******************************************************#

    # Variables for Top
    # Error Function * Exponential to model the bulk

        rrv_number_signal = old_workspace.var("rate_BulkWW_xww_for_unbin");
        
        rrv_number_Total_background_MC = RooRealVar("rrv_number_Total_background_MC_xww","rrv_number_Total_background_MC_xww",
                rrv_number_WJets.getVal()+
                rrv_number_VV.getVal()+
                rrv_number_TTbar.getVal()+
                rrv_number_STop.getVal());

        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));

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' )