Python RooRealVar.getError方法代码示例

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

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
def dump_simple_simul():
    # idea: take two data sets, both low counts, but same lambda
    # and fit said lambda simultanously with both sets

    r_x = RooRealVar('x', 'x', 0, bins)

    r_lam = RooRealVar('lam', 'lam', 0.5, 0.0, 1.5)
    r_lam1 = RooRealVar('lam1', 'lam1', 0.5, 0.0, 1.5)
    r_lam2 = RooRealVar('lam2', 'lam2', 0.5, 0.0, 1.5)

    model1 = RooPoisson('pois1', 'pois1', r_x, r_lam1)
    model2 = RooPoisson('pois2', 'pois2', r_x, r_lam2)

    r_index = RooCategory('index', 'index')
    r_index.defineType('1')
    r_index.defineType('2')

    simul_model = RooSimultaneous('model', 'model', r_index)
    simul_model.addPdf(model1, '1')
    simul_model.addPdf(model2, '2')

    data1 = RooDataSet('', '', RooArgSet(r_x))
    for val in unbinned_from_binned(xdata, ydata):
        r_x.setVal(val)
        data1.add(RooArgSet(r_x))
    r_index.setLabel('1')
    data1.addColumn(r_index)

    data2 = RooDataSet('', '', RooArgSet(r_x))
    for val in unbinned_from_binned(xdata, ydata2):
        r_x.setVal(val)
        data2.add(RooArgSet(r_x))
    r_index.setLabel('2')
    data2.addColumn(r_index)

    data1.append(data2)

    _result = simul_model.fitTo(data1)
    print(r_lam.getVal(), '+-', r_lam.getError(), lam, np.abs(r_lam.getVal() - lam) / lam)
    print(r_lam1.getVal(), '+-', r_lam1.getError(), lam, np.abs(r_lam1.getVal() - lam) / lam)
    print(r_lam2.getVal(), '+-', r_lam2.getError(), lam, np.abs(r_lam2.getVal() - lam) / lam)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
def roofit_poisson_unbinned(data):
    """returns lambda, error of lambda"""

    x = RooRealVar('x', 'x', 0, max(data) * 10)
    lam = RooRealVar('lambda', 'lambda', 0.1, 0.000001, max(data))

    model = RooPoisson('model', 'model', x, lam)

    dataset = RooDataSet('data', 'data', RooArgSet(x))
    for val in data:
        x.setVal(val)
        dataset.add(RooArgSet(x))

    model.fitTo(dataset, RooFit.Save(), RooFit.PrintLevel(-1))
    return lam.getVal(), lam.getError()

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

#.........这里部分代码省略.........
    frTop = modelTop.fitTo(setTop, 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
    iSBTop = modelTop.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange,HSBrange"))
    iLSBTop = modelTop.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange"))
    iHSBTop = modelTop.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("HSBrange"))
    iSRTop = modelTop.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("SRrange"))
    iVRTop = modelTop.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("VRrange"))
    # Do not remove the following lines, integrals are computed here
    iALTop = modelTop.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg))
    nSBTop = iSBTop.getVal()/iALTop.getVal()*setTop.sumEntries(SBcut)
    nLSBTop = iLSBTop.getVal()/iALTop.getVal()*setTop.sumEntries(LSBcut)
    nHSBTop = iHSBTop.getVal()/iALTop.getVal()*setTop.sumEntries(HSBcut)
    nSRTop = iSRTop.getVal()/iALTop.getVal()*setTop.sumEntries(SRcut)
    
    drawPlot("JetMass_Top", channel, J_mass, modelTop, setTop, binsJmass, frTop)
    
    if VERBOSE: print "********** Fit result [JET MASS TOP] ***"+"*"*40, "\n", frTop.Print(), "\n", "*"*80
    
    #*******************************************************#
    #                                                       #
    #                 All bkg normalization                 #
    #                                                       #
    #*******************************************************#
    
    constVjet.setConstant(True)
    offsetVjet.setConstant(True)
    widthVjet.setConstant(True)
    a0Vjet.setConstant(True)
    a1Vjet.setConstant(True)
    a2Vjet.setConstant(True)
    
    constVV.setConstant(True)
    offsetVV.setConstant(True)
    widthVV.setConstant(True)
    meanVV.setConstant(True)
    sigmaVV.setConstant(True)
    fracVV.setConstant(True)
    meanVH.setConstant(True)
    sigmaVH.setConstant(True)
    fracVH.setConstant(True)
    
    constTop.setConstant(True)
    offsetTop.setConstant(True)
    widthTop.setConstant(True)
    meanW.setConstant(True)
    sigmaW.setConstant(True)
    fracW.setConstant(True)
    meanT.setConstant(True)
    sigmaT.setConstant(True)
    fracT.setConstant(True)
    
    
    # Final background model by adding the main+secondary pdfs (using 'coef': ratio of the secondary/main, from MC)
    model = RooAddPdf("model", "model", RooArgList(modelTop, modelVV, modelVjet), RooArgList(coef_Top_VVVjet, coef_VV_Vjet))#FIXME
    model.fixAddCoefRange("h_reasonable_range")
    
    # Extended fit model to data in SB
    # all the 3 sidebands (Low / High / the 2 combined) could be used
    # currently using the LOW+HIGH (the others are commented out)
    yieldLSB = RooRealVar("yieldLSB", "Lower SB normalization",  10, 0., 1.e6)
    yieldHSB = RooRealVar("yieldHSB", "Higher SB normalization", 10, 0., 1.e6)
    yieldSB  = RooRealVar("yieldSB",  "All SB normalization",    10, 0., 1.e6)
    #model_ext = RooExtendPdf("model_ext", "extended p.d.f",   model,  yieldLSB)
    #model_ext = RooExtendPdf("model_ext", "extended p.d.f",   model,  yieldHSB)
    model_ext = RooExtendPdf("model_ext", "extended p.d.f",   model,  yieldSB)
    #frMass = model_ext.fitTo(setDataSB, RooFit.ConditionalObservables(RooArgSet(J_mass)),RooFit.SumW2Error(True),RooFit.Extended(True),RooFit.Range("LSBrange"),RooFit.PrintLevel(-1))
    #frMass = model_ext.fitTo(setDataSB, RooFit.ConditionalObservables(RooArgSet(J_mass)),RooFit.SumW2Error(True),RooFit.Extended(True),RooFit.Range("HSBrange"),RooFit.PrintLevel(-1))
    #frMass = model_ext.fitTo(setDataSB, RooFit.ConditionalObservables(RooArgSet(J_mass)), RooFit.SumW2Error(True), RooFit.Extended(True), RooFit.Range("LSBrange,HSBrange"), RooFit.Strategy(2), RooFit.Minimizer("Minuit2"), RooFit.PrintLevel(1 if VERBOSE else -1))
    
    #print "********** Fit result [JET MASS DATA] **"+"*"*40
    #print frMass.Print()
    #print "*"*80
    
    # Calculate integral of the model obtained from the fit to data (fraction of PDF that is within a given region)
    #nSB = model_ext.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange,HSBrange"))
    #nSB = model_ext.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange"))
    #nSB = model_ext.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("HSBrange"))
    #nSR = model_ext.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("SRrange"))
    #nVR = model_ext.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("VRrange"))
    
    # scale the yieldSB from SB to SR using the ratio of the PDFs defined by the two integrals
    SRyield = RooFormulaVar("SRyield", "extrapolation to SR","(@[email protected]*@[email protected]*@4) * @5/@6 [email protected]*@[email protected]*@8", RooArgList(entrySB, entryVV, entryTop, iSBVV, iSBTop, iSRVjet, iSBVjet, iSRVV, iSRTop))
    VRyield = RooFormulaVar("VRyield", "extrapolation to VR","(@[email protected]*@[email protected]*@4) * @5/@6 [email protected]*@[email protected]*@8", RooArgList(entrySB, entryVV, entryTop, iSBVV, iSBTop, iVRVjet, iSBVjet, iVRVV, iVRTop))
    HSByield = RooFormulaVar("SRyield", "extrapolation to SR","(@[email protected]*@[email protected]*@4) * @5/@6 [email protected]*@[email protected]*@8", RooArgList(entryLSB, entryVV, entryTop, iLSBVV, iLSBTop, iHSBVjet, iLSBVjet, iHSBVV, iHSBTop))
    #   RooFormulaVar SRyield("SRyield","extrapolation to SR","(@0/@1)*@2",RooArgList(*nSR,*nSB,yieldLowerSB))
    #   RooFormulaVar SRyield("SRyield","extrapolation to SR","(@0/@1)*@2",RooArgList(*nSR,*nSB,yieldHigherSB))
    #SRyield = RooFormulaVar("SRyield", "extrapolation to SR","(@0/@1)*@2", RooArgList(nSR, nSB, entrySB))
    
    bkgYield            = SRyield.getVal()
    bkgYield_error      = math.sqrt(SRyield.getPropagatedError(frVjet)**2 + SRyield.getPropagatedError(frVV)**2 + SRyield.getPropagatedError(frTop)**2 + (entrySB.getError()*rSBSRVV)**2)
    bkgNorm             = entrySB.getVal() + SRyield.getVal() + VRyield.getVal()
    bkgYield_eig_norm   = RooRealVar("predSR_eig_norm", "expected yield in SR", bkgYield, 0., 1.e6)
    bkgYieldExt         = HSByield.getVal()
    
    drawPlot("JetMass", channel, J_mass, model, setDataSB, binsJmass, None, None, "", bkgNorm, True)

    
    print channel, "normalization = %.3f +/- %.3f, observed = %.0f" % (bkgYield, bkgYield_error, setDataSR.sumEntries() if not BLIND else -1)
    if VERBOSE: raw_input("Press Enter to continue...")

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
datacard.write('jmax 1\n')
datacard.write('kmax *\n')
datacard.write('---------------\n')
datacard.write('shapes * * '+wsFN+' w:$PROCESS\n')
datacard.write('---------------\n')
datacard.write('bin 1\n')    
datacard.write('observation -1\n')
datacard.write('------------------------------\n')
datacard.write('bin          1          1\n')          
datacard.write('process      signal     background\n')
datacard.write('process      0          1\n')          
datacard.write('rate         '+str(ExpectedSignalRate)+'      1\n')
datacard.write('------------------------------\n')      
#nuisance parameters --- gaussian prior
if bkgNuisance:
  datacard.write('p1  param    '+str(p1.getValV())+'   '+str(p1.getError())+'\n')
  datacard.write('p2  param    '+str(p2.getValV())+'   '+str(p2.getError())+'\n')
  datacard.write('p3  param    '+str(p3.getValV())+'   '+str(p3.getError())+'\n')
#flat parameters --- flat prior
if not (bkgNuisance or histpdfBkg): 
  datacard.write('background_norm  flatParam\n')
  datacard.write('p1  flatParam\n')
  datacard.write('p2  flatParam\n')
  datacard.write('p3  flatParam\n')


   
 
##----- keep the GUI alive ------------
#if __name__ == '__main__':
#  rep = ''

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

#.........这里部分代码省略.........
#		res5.Print()
#		res6.Print()
#		res7.Print()

            # decorrelated background parameters for Bayesian limits
            if args.decoBkg:
                signal_norm.setConstant()
                res = model.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
                res.Print()
                ## temp workspace for the PDF diagonalizer
                w_tmp = RooWorkspace("w_tmp")
                deco = PdfDiagonalizer("deco",w_tmp,res)
                # here diagonalizing only the shape parameters since the overall normalization is already decorrelated
                background_deco = deco.diagonalize(background)
                print "##################### workspace for decorrelation"
                w_tmp.Print("v")
                print "##################### original parameters"
                background.getParameters(rooDataHist).Print("v")
                print "##################### decorrelated parameters"
                # needed if want to evaluate limits without background systematics
                if args.fixBkg:
                    w_tmp.var("deco_eig1").setConstant()
                    w_tmp.var("deco_eig2").setConstant()
                    if not args.fixP3: w_tmp.var("deco_eig3").setConstant()
                background_deco.getParameters(rooDataHist).Print("v")
                print "##################### original pdf"
                background.Print()
                print "##################### decorrelated pdf"
                background_deco.Print()
                # release signal normalization
                signal_norm.setConstant(kFALSE)
                # set the background normalization range to +/- 5 sigma
                bkg_val = background_norm.getVal()
                bkg_error = background_norm.getError()
                background_norm.setMin(bkg_val-5*bkg_error)
                background_norm.setMax(bkg_val+5*bkg_error)
                background_norm.Print()
                # change background PDF names
                background.SetName("background_old")
                background_deco.SetName("background")

        # needed if want to evaluate limits without background systematics
        if args.fixBkg:
            background_norm.setConstant()
            p1.setConstant()
            p2.setConstant()
            p3.setConstant()

        # -----------------------------------------
        # dictionaries holding systematic variations of the signal shape
        hSig_Syst = {}
        hSig_Syst_DataHist = {}
        sigCDF = TGraph(hSig.GetNbinsX()+1)

        # JES and JER uncertainties
        if args.jesUnc != None or args.jerUnc != None:

            sigCDF.SetPoint(0,0.,0.)
            integral = 0.
            for i in range(1, hSig.GetNbinsX()+1):
                x = hSig.GetXaxis().GetBinLowEdge(i+1)
                integral = integral + hSig.GetBinContent(i)
                sigCDF.SetPoint(i,x,integral)

        if args.jesUnc != None:
            hSig_Syst['JESUp'] = copy.deepcopy(hSig)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
    sigmodel.fitTo(signal,RooFit.SumW2Error(True)) #,RooFit.Range("signal")
    chi2=RooChi2Var("chi2","chi2",sigmodel,signal,RooFit.DataError(RooAbsData.SumW2))
    nbins=data.numEntries()
    nfree=sigmodel.getParameters(data).selectByAttrib("Constant",False).getSize()
    s0.setConstant(True) 
    s1.setConstant(True) 

    if cut>=1:
      fullintegral=sumsighist.Integral()
    else:
      fullintegral=sighist.Integral()
    print "SIGNAL FRACTION",nsigref.getValV()/(nsigref.getValV()+nbkg.getValV())
    if nsigref.getValV()==0: continue

    if fit=="data":
      xframe=mass.frame(RooFit.Title("signal fraction in peak ="+str(int(nsigref.getValV()/(nsigref.getValV()+nbkg.getValV())*1000.)/1000.)+"+-"+str(int(nsigref.getError()/(nsigref.getValV()+nbkg.getValV())*1000.)/1000.)+", #chi^{2}/DOF = "+str(int((chi2.getVal()/(nbins-nfree))*10.)/10.)))
      signal.plotOn(xframe,RooFit.DataError(RooAbsData.SumW2))
      sigmodel.plotOn(xframe,RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      sigmodel.plotOn(xframe,RooFit.Components("sigbkg"+str(cut)),RooFit.LineStyle(kDashed),RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      sigmodel.plotOn(xframe,RooFit.Components("sig"),RooFit.LineStyle(kDotted),RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      canvas=TCanvas("c3","c3",0,0,600,600)
      xframe.Draw()
      canvas.SaveAs(prefix+"_"+plot[0]+str(cut)+"_sigfit.pdf")

    a0=RooRealVar("a0"+str(cut),"a0",100.,0.,1000.)
    a1=RooRealVar("a1"+str(cut),"a1",100.,0.,1000.)
    a2=RooRealVar("a2"+str(cut),"a2",50.,0.,1000.)
    a3=RooRealVar("a3"+str(cut),"a3",0.1,0.,1000.)
    b0=RooRealVar("b0"+str(cut),"b0",120.,0.,100.)
    b1=RooRealVar("b1"+str(cut),"b1",50.,0.,100.)
    b2=RooRealVar("b2"+str(cut),"b2",50.,0.,100.)

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

#.........这里部分代码省略.........
    # 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(), 
                        )
    summary = 'fit in signal region\n'
    summary += 'nsig: ' + str( nsig.getValV() ) + ' +- ' + str( nsig.getError() ) + '\n' 
    summary += 'nbkg: ' + str( nbkg.getValV() ) + ' +- ' + str( nbkg.getError() ) + '\n'
#     
#     model_ctl.fitTo( real_data_ctl_hist )
#     summary += 'fit in control region\n'
#     summary += 'nsig: ' + str( nsig.getValV() ) + ' +- ' + str( nsig.getError() ) + '\n' 
#     summary += 'nbkg: ' + str( nbkg.getValV() ) + ' +- ' + str( nbkg.getError() ) + '\n' 
# 
#     # Perform simultaneous fit of model to data and model_ctl to data_ctl
#     simPdf.fitTo( combData ) 
#     summary += 'Combined fit\n'
#     summary += 'nsig: ' + str( nsig.getValV() ) + ' +- ' + str( nsig.getError() ) + '\n' 
#     summary += 'nbkg: ' + str( nbkg.getValV() ) + ' +- ' + str( nbkg.getError() ) + '\n' 


    # 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" ),
                   RooFit.ProjWData( RooArgSet( sample ), combData ), ) 
    simPdf.plotOn( frame1, RooFit.Slice( sample, "physics" ),
                   RooFit.Components( "signal_1_pdf" ),
                   RooFit.ProjWData( RooArgSet( sample ), combData ),

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

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

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
                print "##################### decorrelated parameters"
                # needed if want to evaluate limits without background systematics
                if args.fixBkg:
                    w_tmp.var("deco_eig1").setConstant()
                    w_tmp.var("deco_eig2").setConstant()
                    if not args.fixP3: w_tmp.var("deco_eig3").setConstant()
                background_deco.getParameters(rooDataHist).Print("v")
                print "##################### original pdf"
                background.Print()
                print "##################### decorrelated pdf"
                background_deco.Print()
                # release signal normalization
                signal_norm.setConstant(kFALSE)
                # set the background normalization range to +/- 5 sigma
                bkg_val = background_norm.getVal()
                bkg_error = background_norm.getError()
                background_norm.setMin(bkg_val-5*bkg_error)
                background_norm.setMax(bkg_val+5*bkg_error)
                background_norm.Print()
                # change background PDF names
                background.SetName("background_old")
                background_deco.SetName("background")

        # needed if want to evaluate limits without background systematics
        if args.fixBkg:
            background_norm.setConstant()
            p1.setConstant()
            p2.setConstant()
            p3.setConstant()

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

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
fitResult = model.fitTo(data, RooFit.Minimizer("Minuit2", "Migrad"), 
                        RooFit.NumCPU(8), RooFit.Extended(True),
                        RooFit.SumW2Error(False), RooFit.Strategy(2), 
                        #verbosity
                        RooFit.PrintLevel(-1), RooFit.Warnings(False), RooFit.Verbose(False))

#ntt_fit = ntt.getVal();
nSignal_fit = nSignal.getVal();
nwj_fit = nwj.getVal();
#nzj_fit = nzj.getVal();
nVJ_fit = nVJ.getVal();
nqcd_fit = nqcd.getVal();
#nstop_fit = nstop.getVal();


nSignal_fiterr = nSignal.getError();
#ntt_fiterr = ntt.getError();
nwj_fiterr = nwj.getError();
nVJ_fiterr = nVJ.getError();
nqcd_fiterr = nqcd.getError();
#nstop_fiterr = nstop.getError();

print "Total events in signal region: ", n_event_obs;
print 'DATA:', n_event_obs
print 'N_tt:', n_ttbar
print 'N_signal', n_signal
print 'N_WJ:', n_wj
print 'N_VJ:', n_VJ
print 'N_sTop:', n_stop
print 'N_QCD:', n_qcd
print 'SumMC:', n_signal + n_VJ + n_qcd

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

#.........这里部分代码省略.........
        N_signal_before_fit = N_ttbar_before_fit + N_SingleTop_before_fit
        
        N_ttbar_error_before_fit = sum(histograms['TTJet'].errors())
        N_SingleTop_error_before_fit = sum(histograms['SingleTop'].errors())
        N_vjets_error_before_fit = sum(histograms['V+Jets'].errors())
        N_QCD_error_before_fit = sum(histograms['QCD'].errors())
        
        if (N_SingleTop_before_fit != 0):
            TTJet_SingleTop_ratio = N_ttbar_before_fit / N_SingleTop_before_fit
        else:
            print 'Bin ', variable_bin, ': ttbar/singleTop ratio undefined for %s channel! Setting to 0.' % channel
            TTJet_SingleTop_ratio = 0
            
        
        leptonAbsEta = RooRealVar("leptonAbsEta", "leptonAbsEta", 0., 2.4)
        # this has to move to dps.utils.Fitting.py
        vars = RooArgList()
        vars.add(leptonAbsEta)
        vars_set = RooArgSet()
        vars_set.add(leptonAbsEta)
        n_event_obs = histograms['data'].Integral()
        
        lowerBound = 0. 
        upperBound = n_event_obs + 10 * sqrt(n_event_obs)
        n_init = n_event_obs / 2.
        
        data = RooDataHist("data", "dataset with leptonAbsEta", vars, histograms['data'])
        rh_vj = RooDataHist("rh_vj", "vj", vars, histograms['V+Jets'])
        rh_qcd = RooDataHist("rh_qcd", "qcd", vars, histograms['QCD'])
        rh_signal = RooDataHist("rh_signal", "signal", vars, h_eta_signal)
        
        pdf_vj = RooHistPdf ("pdf_vj", "V+Jets pdf", vars_set, rh_vj, 0)
        pdf_qcd = RooHistPdf("pdf_qcd", "QCD pdf ", vars_set, rh_qcd, 0)
        pdf_signal = RooHistPdf("pdf_signal", "single top pdf", vars_set, rh_signal, 0)
        
        # RooRealVar(const char *name, const char *title, Double_t value, Double_t minValue, Double_t maxValue, const char *unit) :
        nSignal = RooRealVar("nSignal", "number of single top + ttbar events", N_signal_before_fit, lowerBound, upperBound, "event")
        nvj = RooRealVar  ("nvj", "number of V+Jets bgnd events", N_vjets_before_fit, lowerBound, upperBound, "event")
        nqcd = RooRealVar("nqcd", "number of QCD bgnd events", N_QCD_error_before_fit, lowerBound, upperBound, "event")
        
        model = RooAddPdf("model", "sig+vj+qcd",
                          RooArgList(pdf_signal, pdf_vj, pdf_qcd),
                          RooArgList(nSignal, nvj, nqcd)
                          )
        vj_constraint = RooGaussian("nvj_constraint", "nvj_constraint", nvj, RooFit.RooConst(N_vjets_before_fit), RooFit.RooConst(0.5 * N_vjets_before_fit))
        qcd_constraint = RooGaussian("nqcd_constraint", "nqcd_constraint", nqcd, RooFit.RooConst(N_qcd_before_fit), RooFit.RooConst(2 * N_qcd_before_fit))  
        model_with_constraints = RooProdPdf("model_with_constraints", "model with gaussian constraints",
                                            RooArgSet(model, vj_constraint, qcd_constraint), RooLinkedList())
        model_with_constraints.fitTo(data, RooFit.Minimizer("Minuit2", "Migrad"))  #WARNING: number of cores changes the results!!!
#         nll = model.createNLL(data, RooFit.NumCPU(2))
#         RooMinuit(nll).migrad()
#         frame1 = nSignal.frame(RooFit.Bins(100), RooFit.Range(lowerBound, n_event_obs), RooFit.Title("LL and profileLL in nSignal")) 
#         nll.plotOn(frame1, RooFit.ShiftToZero())
#         frame2 = nvj.frame(RooFit.Bins(100), RooFit.Range(lowerBound, n_event_obs), RooFit.Title("LL and profileLL in nvj"))
#         nll.plotOn(frame2, RooFit.ShiftToZero())
#         frame3 = nqcd.frame(RooFit.Bins(100), RooFit.Range(lowerBound, n_event_obs), RooFit.Title("LL and profileLL in nqcd"))
#         nll.plotOn(frame3, RooFit.ShiftToZero())  
#         
#         pll_nSignal = nll.createProfile(nSignal)
#         pll_nSignal.plotOn(frame1, RooFit.LineColor(2)) 
#         frame1.SetMinimum(0)
#         frame1.SetMaximum(3)
#         
#         pll_nvj = nll.createProfile(nvj)
#         pll_nvj.plotOn(frame2, RooFit.LineColor(2)) 
#         frame2.SetMinimum(0)
#         frame2.SetMaximum(3)
#         
#         pll_nqcd = nll.createProfile(nqcd)
#         pll_nqcd.plotOn(frame3, RooFit.LineColor(2)) 
#         frame3.SetMinimum(0)
#         frame3.SetMaximum(3)
#         c = TCanvas("profilell","profilell",1200, 400)
#         c.Divide(3)
#         c.cd(1)
#         frame1.Draw()
#         c.cd(2)
#         frame2.Draw()
#         c.cd(3)
#         frame3.Draw()
#         c.SaveAs('profileLL.png')
#         model.fitTo(data, RooFit.Minimizer("Minuit2", "Migrad"), RooFit.NumCPU(1))#WARNING: number of cores changes the results!!!
        fit_results = {}
        fit_results['signal'] = (nSignal.getVal(), nSignal.getError())
        fit_results['QCD'] = ufloat(nqcd.getVal(), nqcd.getError())
        fit_results['V+Jets'] = ufloat(nvj.getVal(), nvj.getError())
        
        N_ttbar, N_SingleTop = decombine_result(fit_results['signal'], TTJet_SingleTop_ratio)
        fit_results['signal'] = ufloat(nSignal.getVal(), nSignal.getError())
        fit_results['TTJet'] = ufloat(N_ttbar)
        fit_results['SingleTop'] = ufloat(N_SingleTop)
        
        if results == {}:  # empty
            for sample in fit_results.keys():
                results[sample] = [fit_results[sample]]
        else:
            for sample in fit_results.keys():
                results[sample].append(fit_results[sample])
        
    return results, None, None

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
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", 40, 200)
    nsig = RooRealVar("nsig", "#signal events", 200, 0.0, 10000)
    nbkg = RooRealVar("nbkg", "#background events", 800, 0.0, 200000)
    # 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)

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

    # Construct composite pdf
    model = RooAddPdf("model", "model", RooArgList(gx, px), 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", 40, 200)

    mean_ctl = RooRealVar("mean_ctl", "mean_ctl", mu2, 40, 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, 40, 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")

    # 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)
    summary = "fit in signal region\n"
    summary += "nsig: " + str(nsig.getValV()) + " +- " + str(nsig.getError()) + "\n"
    summary += "nbkg: " + str(nbkg.getValV()) + " +- " + str(nbkg.getError()) + "\n"

    model_ctl.fitTo(real_data_ctl_hist)
    summary += "fit in control region\n"
    summary += "nsig: " + str(nsig.getValV()) + " +- " + str(nsig.getError()) + "\n"
    summary += "nbkg: " + str(nbkg.getValV()) + " +- " + str(nbkg.getError()) + "\n"

    # Perform simultaneous fit of model to data and model_ctl to data_ctl
    simPdf.fitTo(combData)
    summary += "Combined fit\n"
    summary += "nsig: " + str(nsig.getValV()) + " +- " + str(nsig.getError()) + "\n"
    summary += "nbkg: " + str(nbkg.getValV()) + " +- " + str(nbkg.getError()) + "\n"

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

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getError [as 别名]
    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()
    ## fpar[4] = sig_hist.GetEntries()*fped.getVal()
    ## fpar[5] = ws.var('pedMean').getVal()
    ## fpar[6] = ws.var('pedWidth').getVal()
    
    ## maxx = Double(0.)
    ## fwhm = Double(0.)

    maxx = findMax(lxg, x) #, mpv.getVal(), width.getVal())
    maxxErr = maxx*mpv.getError()/mpv.getVal()

    if havePeds:
        maxx -= ws.var('pedMean').getVal()
        maxxErr = sqrt(ws.var('pedMean').getError()**2 + maxxErr**2)

    ## langaupro(fpar,maxx, fwhm)

    fr.Print('v')
    print 'for {0} MIP MPV: {1:0.2f} +/- {2:0.2f}'.format(HOTower,maxx,
                                                          maxxErr), \
          'S/N: {0:0.2f}'.format(maxx/savePedWidth),
          #'FWHM: {0:0.2f}'.format(fwhm),\
    if opts.sipm and havePeds:
        onePE = ws.var('peMean').getVal() # - ws.var('pedMean').getVal()
        peErr = ws.var('peMean').getError()