Python RooWorkspace.set方法代码示例

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
fitfunc.SetParName(3,"Tail Z")
fitfunc.SetParName(4,"Tail length")


w = RooWorkspace("w")
w.factory("{0}[0,0.1]".format("uncM"))
w.factory("uncVZ[-100,100]")
w.factory("uncP[0,10]")
w.factory("cut[0,1]")

w.defineSet("myVars","{0},uncVZ".format("uncM"))

events = infile.Get("cut")
#eventsrad = radfile.Get("ntuple")
eventsprompt = promptfile.Get("cut")
dataset = RooDataSet("data","data",events,w.set("myVars"),"")

w.factory("Gaussian::vtx_model(uncVZ,mean[-50,50],sigma[0,50])")
gauss_pdf = w.pdf("vtx_model")
w.factory("EXPR::gaussExp('exp( ((@[email protected])<@3)*(-0.5*(@[email protected])^2/@2^2) + ((@[email protected])>[email protected])*(-0.5*@3^2/@2^2-(@[email protected]@3)/@4))',uncVZ,gauss_mean[-5,-20,20],gauss_sigma[5,1,50],exp_breakpoint[10,0,50],exp_length[3,0.5,20])")
gaussexp_pdf = w.pdf("gaussExp")
w.defineSet("obs_1d","uncVZ")
obs=w.set("obs_1d")
uncVZ = w.var("uncVZ")
uncVZ.setBins(200)
gauss_params = gauss_pdf.getParameters(obs)
gaussexp_params = gaussexp_pdf.getParameters(obs)

Medges = array.array('d')
Epsedges = array.array('d')
for i in range(0,nApMass+1):

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
class Wjj2DFitter:

    def __init__ (self, pars):
        self.pars = pars
        self.ws = RooWorkspace('wjj2dfitter')
        self.utils = Wjj2DFitterUtils(self.pars)
        self.useImportPars = False

        self.rangeString = None
        obs = []
        for v in self.pars.var:

            try:
                vName = self.pars.varNames[v]
            except AttributeError:
                vName = v

            obs.append(vName)
            var1 = self.ws.factory('%s[%f,%f]' % (vName, 
                                                  self.pars.varRanges[v][1], 
                                                  self.pars.varRanges[v][2])
                                   )
            var1.setUnit('GeV')
            try:
                var1.SetTitle(self.pars.varTitles[v])
            except AttributeError:
                var1.SetTitle('m_{jj}')
            var1.setPlotLabel(var1.GetTitle())
            if len(self.pars.varRanges[v][3]) > 1:
                vbinning = RooBinning(len(self.pars.varRanges[v][3]) - 1, 
                                   array('d', self.pars.varRanges[v][3]),
                                   '%sBinning' % vName)
                var1.setBinning(vbinning)
            else:
                var1.setBins(self.pars.varRanges[v][0])
            var1.Print()
            if v in self.pars.exclude:
                var1.setRange('signalRegion', self.pars.exclude[v][0],
                              self.pars.exclude[v][1])
                var1.setRange('lowSideband', var1.getMin(), 
                              self.pars.exclude[v][0])
                var1.setRange('highSideband', self.pars.exclude[v][1],
                              var1.getMax())
                self.rangeString = 'lowSideband,highSideband'

            if hasattr(self.pars, 'plotRanges'):
                var1.setRange('plotRange', self.pars.plotRanges[v][1],
                              self.pars.plotRanges[v][2])
                var1.setBins(self.pars.plotRanges[v][0], 'plotBins')
            else:
                var1.setRange('plotRange', var1.getMin(), var1.getMax())
                var1.setBins(var1.getBins(), 'plotBins')

        self.ws.defineSet('obsSet', ','.join(obs))

    def loadDataFromWorkspace(self, other, cut = None):
        #pull unbinned data from other workspace
        unbinnedData = other.data('data_unbinned')
        if not unbinnedData:
            unbinnedData = other.data('data_obs')

        if cut:
            unbinnedData = unbinnedData.reduce(cut)

        unbinnedData.Print()
        if self.pars.binData:
            #bin and import data
            unbinnedData.SetName('data_unbinned')
            getattr(self.ws, 'import')(unbinnedData)
            data = RooDataHist('data_obs', 'data_obs', other.set('obsSet'), 
                               unbinnedData)
            getattr(self.ws, 'import')(data)
        else:
            #just import data
            unbinnedData.SetName('data_obs')
            getattr(self.ws, 'import')(unbinnedData)

    def loadHistogramsFromWorkspace(self, other):
        #pull RooHist pdfs from other workspace
        pdfs = other.allPdfs()
        pdfIter = pdfs.createIterator()
        pdf = pdfIter.Next()
        while pdf:
            if pdf.IsA().InheritsFrom('RooHistPdf'):
                print 'importing',pdf.GetName(),'from old workspace'
                getattr(self.ws, 'import')(pdf)
            pdf = pdfIter.Next()

    def loadWorkspaceFromFile(self, filename, wsname = 'w', 
                              getFloatPars = True):
        print 'loading data workspace %s from file %s' % (wsname, filename)
        fin = TFile.Open(filename)
        if not fin:
            print 'failed to open the file',filename
            import os
            print 'cwd:',os.getcwd()
            print 'access of',filename,os.access(filename, os.R_OK)
            print 'list of root files in cwd'
            for f in os.listdir(os.getcwd()):
                if f[-5:] == '.root':
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
fcUpperHist=TH2D("fcUpperLimit","fcUpperLimit",n_massbins,xedges,n_epsbins,yedges)
plrPvalHist=TH2D("plrPval","plrPval",n_massbins,xedges,n_epsbins,yedges)
plrSigHist=TH2D("plrSig","plrSig",n_massbins,xedges,n_epsbins,yedges)
logplrHist=TH2D("logplr","logplr",n_massbins,xedges,n_epsbins,yedges)
candRescaledHist=TH1D("candidates_rescaled","candidates_rescaled",100,0,1.0)
candRescaled2DHist=TH2D("candidates_rescaled_2d","candidates_rescaled_2d",n_massbins,xedges,100,0,1.0)

w = RooWorkspace("w")
w.factory("{0}[0,0.1]".format(massVar))
w.factory("uncVZ[-100,100]")
w.factory("uncP[0,10]")
w.factory("cut[0,1]")

w.defineSet("myVars","{0},uncVZ".format(massVar))

dataset = RooDataSet("data","data",events,w.set("myVars"),"")

w.factory("Gaussian::vtx_model(uncVZ,mean[-50,50],sigma[0,50])")
gauss_pdf = w.pdf("vtx_model")
w.factory("EXPR::gaussExp('exp( ((@[email protected])<@3)*(-0.5*(@[email protected])^2/@2^2) + ((@[email protected])>[email protected])*(-0.5*@3^2/@2^2-(@[email protected]@3)/@4))',uncVZ,gauss_mean[-5,-20,20],gauss_sigma[5,1,50],exp_breakpoint[10,0,50],exp_length[3,0.5,20])")
gaussexp_pdf = w.pdf("gaussExp")
w.defineSet("obs_1d","uncVZ")
obs=w.set("obs_1d")
uncVZ = w.var("uncVZ")
uncVZ.setBins(200)
gauss_params = gauss_pdf.getParameters(obs)
gaussexp_params = gaussexp_pdf.getParameters(obs)

exppol4=TF1("exppol4","exp(pol4(0))",-5,100)
w.factory("EXPR::signal('(@0>{0})*exp(@1 + @2*@0 + @3*@0^2 + @4*@0^3 + @5*@0^4)',uncVZ,eff_p0[-1,1],eff_p1[-1,1],eff_p2[-1,1],eff_p3[-1,1],eff_p4[-1,1])".format(targetz))
w.factory("SUM::model(strength[0,1]*signal,gaussExp)")

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
def main(options,args):

    cfg = options.config
    workspaceName = cfg.get('Global','workspace')        
    
    ws = RooWorkspace(workspaceName)    

    #ws.Print("v")
    
    setupWorkspace(ws,options)        

    #create -log(likelihood)
    
    theNLL = ws.pdf('TopLevelPdf').createNLL(ws.data('allcountingdata'),
                                             RooFit.NumCPU(1),
                                             RooFit.ConditionalObservables(ws.set('condObs')),
                                             RooFit.Verbose(True))

    ws.saveSnapshot('standardmodel',ws.allVars())
    
    minuit = ROOT.RooMinuit(theNLL)
    minuit.setPrintLevel(1)
    minuit.setPrintEvalErrors(-1)
    minuit.setErrorLevel(.5)
    
    #find the values of the parameters that minimize the likelihood
    minuit.setStrategy(2)
    minuit.simplex()
    minuit.migrad()
    minuit.hesse()

    #ws.var('err_gl').setConstant(True)
    #ws.var('err_gs').setConstant(True)
    #ws.var('err_gb').setConstant(True)

    ws.defineSet('POI',
                 ROOT.RooArgSet(ws.var('%s_%s'%(cfg.get('Global','par1Name'),cfg.get('Global','couplingType'))),
                                ws.var('%s_%s'%(cfg.get('Global','par2Name'),cfg.get('Global','couplingType')))))

    ws.saveSnapshot('%s_fitresult'%cfg.get('Global','couplingType'),
                    ws.allVars())
        
    #create profile likelihood       
    level_68 = ROOT.TMath.ChisquareQuantile(.68,2)/2.0 # delta NLL for 68% confidence level for -log(LR)
    level_95 = ROOT.TMath.ChisquareQuantile(.95,2)/2.0 # delta NLL for 95% confidence level for -log(LR)

    print
    print '68% CL Delta-NLL 2 DOF=',level_68
    print '95% CL Delta-NLL 2 DOF=',level_95
    
    
    minuit.setPrintLevel(1)
    minuit.setPrintEvalErrors(-1)

    minuit.migrad()
    minuit.minos(ws.set('POI'))

    thePlot = minuit.contour(ws.var('%s_%s'%(cfg.get('Global','par1Name'),cfg.get('Global','couplingType'))),
                             ws.var('%s_%s'%(cfg.get('Global','par2Name'),cfg.get('Global','couplingType'))),
                             sqrt(2*level_95),sqrt(2*level_68)) # here the error is in sigmas 

    thePlot.SetName('%s_%s_%s_contour'%(cfg.get('Global','par1Name'),
                                        cfg.get('Global','par2Name'),
                                        cfg.get('Global','couplingType')))    
    
    thePlot.SetTitle('68% & 95% CL on the Best Fit Values of '+cfg.get('Global','par1Name')+' and '+cfg.get('Global','par2Name'))
    legend = ROOT.TLegend(2.01612903225806439e-01,7.86016949152542388e-01,
                          7.15725806451612989e-01,9.13135593220338992e-01)
    legend.SetNColumns(2)
    thePlot.addObject(legend)

    # 1-D Limits
    
    level_95 = ROOT.TMath.ChisquareQuantile(.95,1)/2.0 # delta NLL for -log(LR) with 1 dof
    print '95% CL Delta-NLL 1 DOF=',level_95
    minuit.setErrorLevel(level_95)

    #set 1-D limits on parameter 1 with parameter 2 == 0
    ws.var('%s_%s'%(cfg.get('Global','par2Name'),cfg.get('Global','couplingType'))).setVal(0.0)
    ws.var('%s_%s'%(cfg.get('Global','par2Name'),cfg.get('Global','couplingType'))).setConstant(True)
    minuit.minos(ws.set('POI'))

    parm1 = ws.var('%s_%s'%(cfg.get('Global','par1Name'),cfg.get('Global','couplingType')))

    print  'parameter 1 value: '+str(parm1.getVal())

    if not (0 < parm1.getVal()+parm1.getErrorHi() and 0 > parm1.getVal()+parm1.getErrorLo()):
        print '95% CL does not cover SM for parameter 1'
    else:
        print '95% CL covers SM for parameter 1'

    par1Line = ROOT.TLine(parm1.getVal()+parm1.getErrorLo(),0,
                          parm1.getVal()+parm1.getErrorHi(),0)
    par1Line.SetLineWidth(2)
    par1Line.SetLineColor(ROOT.kRed)
    
    thePlot.addObject(par1Line)

    #set 1-D limits on parameter 2 with parameter 1 == 0
    ws.var('%s_%s'%(cfg.get('Global','par2Name'),cfg.get('Global','couplingType'))).setConstant(False)
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
mc.SetPdf(simPdf)
mc.SetObservables(obs)
mc.SetParametersOfInterest(poi)
mc.SetNuisanceParameters(nuis)

plc = RooStats.ProfileLikelihoodCalculator(data, mc)
plc.SetConfidenceLevel(0.68)

NullSnapshotName = 'null_snapshot'
ws.defineSet('params', pars, False)
ws.saveSnapshot(NullSnapshotName, nullPars, True)
#ws.Print()

ws.loadSnapshot(NullSnapshotName)

nullSet = ws.set('params')
#nullSet.Print('v')
plc.SetNullParameters(nullSet)

htr = plc.GetHypoTest()
signif = htr.Significance()

pars.Print("v")
print "suppression significance:",signif

interval = plc.GetInterval()

# plPlot = RooStats.LikelihoodIntervalPlot(interval)
# plPlot.SetNPoints(20)

## print htrPlot

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
    compPdfs.append(combinedWS.factory('RooExtendPdf::%s_extended(%s, %s)' %\
                                           (comp, comp, norm.GetName())))

combinedWS.defineSet('obsSet', '%s,%s' % (other_var, var_name))

if opts.sigInject:
    combinedWS.var('r_signal').setVal(opts.sigInject)
combinedWS.var('r_signal').setError(0.1)
combinedWS.var('r_signal').setRange(-3., 9.)
combinedWS.var('r_signal').setConstant(False)

compNames = [ c.GetName() for c in compPdfs ]
compList = RooArgList(combinedWS.argSet(','.join(compNames)))
getattr(combinedWS, 'import')(RooAddPdf('total', 'total', compList))
combinedPdf = combinedWS.pdf('total')
genPars = combinedPdf.getParameters(combinedWS.set('obsSet'))
combinedWS.defineSet('params', genPars)

if opts.toy:
    #generate toy dataset
    print 'Generated parameters'
    combinedWS.set('params').Print('v')
    combinedWS.saveSnapshot("genPars", combinedWS.set('params'))

    data = combinedPdf.generate(combinedWS.set('obsSet'), 
                                RooFit.Name('data_obs'),
                                RooFit.Extended())
    data.Print('v')
    getattr(combinedWS, 'import')(data)

extraTag = ''

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import set [as 别名]
wspace . factory("SUM::model_pdf(s[0,15]*model_pdf_1,b[1,0,2]*model_pdf_2)")
wspace . factory("Lognormal::likelihood_b(b,1,3)")
wspace . factory("PROD::model_likelihood(model_pdf, likelihood_b)")
wspace . factory("Uniform::prior_pdf(s)")
  
# define observables
wspace . defineSet("observables","x")

# define parameters of interest
wspace . defineSet("poi","s")

# define nuisance parameters
wspace . defineSet("nuisance_parameters","b")

# load data
observables = RooArgList( wspace.set("observables") )
data = RooDataSet.read("counting_data_3.ascii", observables)
data . SetName("data")
getattr(wspace, 'import')(data)

# model config
modelConfig = RooStats.ModelConfig("counting_model_config")
modelConfig . SetWorkspace(wspace)
modelConfig . SetPdf(wspace.pdf("model_likelihood"))
modelConfig . SetPriorPdf(wspace.pdf("prior_pdf"))
modelConfig . SetParametersOfInterest(wspace.set("poi"))
modelConfig . SetNuisanceParameters(wspace.set("nuisance_parameters"))
getattr(wspace, 'import')(modelConfig, "counting_model_config")

# Bayesian Calculator
bc = RooStats.BayesianCalculator(data, modelConfig)