Python RooWorkspace.data方法代码示例

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import data [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 data [as 别名]
readData(ws, hidatafile, ppdatafile, cuts, mmin, mmax)

mass = ws.var('invariantMass')

ppBkgModel = 1
bkgModel = 0
if useKeys:
    bkgModel = 2

buildPdf(ws, False, ppBkgModel, True)
buildPdf(ws, True, bkgModel, True)
simPdf = buildSimPdf(ws, ws.cat('dataCat'))

#ws.Print()

data = ws.data('data').reduce('(QQsign==QQsign::PlusMinus)')
pars = simPdf.getParameters(data)
ws.var('nsig1_pp').setVal(90)
ws.var('nsig1_pp').setError(10)
ws.var('nsig1_hi').setVal(1200)
ws.var('nsig1_hi').setError(40)

ws.var('nbkg_pp').setVal(335*(mmax-mmin)/7.)
ws.var('nbkg_pp').setError(12)
ws.var('nbkg_hi').setVal(10000*(mmax-mmin)/7.)
ws.var('nbkg_hi').setError(100)

if len(simparamfile) > 0:
    pars.readFromFile(simparamfile)

ws.Print()

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import data [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 data [as 别名]
ws = RooWorkspace("ws","ws")

readData(ws, hidatafile, ppdatafile, cuts)

setBkgdPbPb = opts.bkgd_PbPb
setBkgdpp = opts.bkgd_pp
trkRotBkgd = opts.trkRot

buildPdf(ws, True, setBkgdPbPb, trkRotBkgd)  # pdf for PbPb
buildPdf(ws, False, setBkgdpp, False)   # pdf for pp
simPdf = buildSimPdf(ws, ws.cat('dataCat'))

mass = ws.var('invariantMass')
upsRap = ws.var('Rapidity')
##upsRap.setRange("RapidityBin",'%0.1f' % (opts.ymin),'%0.1f' % ( opts.ymax))
data = ws.data('data').reduce('(QQsign==QQsign::PlusMinus)&&(%s>%0.1f)&&(%s<%0.1f)' % (mass.GetName(), mmin, mass.GetName(), mmax))
##data = ws.data('data').reduce('(QQsign==QQsign::PlusMinus)')
data_ls = ws.data('data').reduce('(QQsign!=0)')
##data = ws.data('data').reduce('(QQsign==QQsign::PlusMinus)&&(%s>%0.1f)&&(%s<%0.1f)' % (upsRap, opts.ymin, upsRap, opts.ymax))
##data_ls = ws.data('data').reduce('(QQsign!=0)&&(%s>%0.1f)&&(%s<%0.1f)' % (upsRap, opts.ymin, upsRap, opts.ymax))

mass.setRange("fitRange",8.5,11.5)
mass.setRange(8.5,11.5)
##parsRap = pdf.getParameters(data)
pars = simPdf.getParameters(data)
#data_hi =  ws.data('data').reduce('(QQsign==QQsign::PlusMinus) && (dataCat == dataCat::hi)')
#data_hi.SetName('data_heavy')
#getattr(ws,'import')(data_hi)
ws.Print()
data.Print()

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import data [as 别名]
ppdatafile = 'data/dimuonTree_2011_pp.root'

cuts = '(muPlusPt > %0.1f) && (muMinusPt > %0.1f) && (abs(upsRapidity)<2.4) && (vProb > 0.05)' \
       % (opts.pt, opts.pt)
simparamfile = opts.paramfile
useKeys = opts.keys

ws = RooWorkspace("ws","ws")

readData(ws, hidatafile, ppdatafile, cuts)

buildPdf(ws, False, 1, True)
buildPdf(ws, True, 0, True)
simPdf = buildSimPdf(ws, ws.cat('dataCat'))

pars = simPdf.getParameters(ws.data('data'))
pars.readFromFile(simparamfile)

Npp_tot = Long(ws.var("nsig1_pp").getVal() + \
               ws.function("nsig2_pp").getVal() + \
               ws.function("nsig3_pp").getVal() + \
               ws.var("nbkg_pp").getVal() + 0.5)
Nhi_tot = Long(ws.var("nsig1_hi").getVal() + \
               ws.function("nsig2_hi").getVal() + \
               ws.function("nsig3_hi").getVal() + \
               ws.var("nbkg_hi").getVal() + 0.5)
dataCat = ws.cat('dataCat')
QQsign = ws.cat('QQsign')
reducedCols = RooArgSet(dataCat,QQsign)

dsToGen = ['os']

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import data [as 别名]
def plotStuff(plotList,plotstring, cutstring, plotfile, plotname, xlabel, ylabel, unitnorm):
    isFirst = True
    w = RooWorkspace("w")
    w.factory("x[-100,100]")
    for dataset in plotList:
        dataset[0].Draw(plotstring.format(dataset[1]),cutstring,"goff")
        hist = gDirectory.Get(dataset[1])
        data=RooDataHist(dataset[1],dataset[1],RooArgList(w.var("x")),hist)
        getattr(w,'import')(data)
    w.factory("HistPdf::triPdf(x,tri)")
    w.factory("HistPdf::wabPdf(x,wab)")
    w.factory("prod::triscale(a[0.3,0,10],{0})".format(w.data("tri").sum(False)))
    w.factory("prod::wabscale(b[0.1,0,10],{0})".format(w.data("wab").sum(False)))
    w.factory("SUM::sumModel(triscale*triPdf,wabscale*wabPdf)")

    w.pdf("sumModel").fitTo(w.data("data"),RooFit.SumW2Error(True),RooFit.Extended(True), RooFit.Verbose(False),RooFit.PrintLevel(-1))
    #w.pdf("sumModel").fitTo(w.data("data"),RooFit.Extended(True))
    #w.pdf("sumModel").fitTo(w.data("data"))
    frame=w.var("x").frame()
    w.data("data").plotOn(frame)
    #w.pdf("triPdf").plotOn(frame)
    #w.pdf("wabPdf").plotOn(frame)
    w.pdf("sumModel").plotOn(frame)
    w.pdf("sumModel").paramOn(frame)
    frame.SetTitle(gDirectory.Get("data").GetTitle())
    frame.Draw()
    c.Print(plotfile)
    c.Clear()

    dataHist = gDirectory.Get("data")
    triHist = gDirectory.Get("tri")
    wabHist = gDirectory.Get("wab")
    if legendright:
        leg = TLegend(0.7,0.75,0.9,0.9)
    else:
        leg = TLegend(0.1,0.75,0.3,0.9)
    hs = THStack("hs",plotname);
    for dataset in plotList:
        hist = gDirectory.Get(dataset[1])
        #hist.Sumw2()
        if unitnorm:
            hist.Scale(1.0/hist.Integral())
        else:
            hist.Scale(1.0/dataset[2])
            print "{0} {1} {2}".format(plotname,dataset[4],hist.Integral())
        hist.SetLineColor(dataset[3])
        leg.AddEntry(hist,dataset[4])
        hs.Add(hist)
        #hist.GetXaxis().SetTitle(xlabel)
        hist.GetYaxis().SetTitle(ylabel)

        #if isFirst:
            #hist.GetXaxis().SetTitle(xlabel)
            #hist.GetYaxis().SetTitle(ylabel)
            #hist.Draw()
        #else:
            #hist.Draw("same")
        isFirst = False
    sumHist = triHist.Clone("sum")
    sumHist.Add(wabHist)
    if unitnorm:
        sumHist.Scale(1.0/sumHist.Integral())
    sumHist.SetLineColor(6)
    leg.AddEntry(sumHist,"MC sum")
    hs.Add(sumHist)
    hs.Draw("nostack")
    hs.GetXaxis().SetTitle(xlabel)
    hs.GetYaxis().SetTitle(ylabel)
    leg.Draw()
    c.Print(plotfile)

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import data [as 别名]
  if tree.itype == -88:
    mc.add(getattr(ws,'set')('observables'))
  if tree.itype == 72 or tree.itype == 82:
    data.add(getattr(ws,'set')('observables'))

getattr(ws,'import')(mc)
getattr(ws,'import')(data)

ws.Print()

ws.factory("mean[5200,5300]")
ws.factory("CBShape:sig_cb1(mass,mean,sigma_1[10,200],alpha_1[0,1.],n_1[0.,50.])");
ws.factory("CBShape:sig_cb2(mass,mean,sigma_2[10,200],alpha_2[-1.,0.],n_2[0.,50.])");
ws.factory("SUM::mcsig( f[0.5,0,1]*sig_cb1, sig_cb2)");

ws.pdf("mcsig").fitTo(ws.data("mc"))

mcplot = ws.var("mass").frame()
ws.data("mc").plotOn(mcplot)
ws.pdf("mcsig").plotOn(mcplot)

# freeze obs
argset = ws.pdf("mcsig").getParameters(RooArgSet(ws.var("mass")))
ws.defineSet("signal_params",argset)
#RooArgSet(getattr(ws,'set')("signal_params")).setAttribAll("Constant")

ws.factory("Exponential:bkg(mass,exp_p1[-0.02,0.])")
ws.factory("SUM:pdf( sig_y[0,1000]*mcsig, bkg_y[0,2000]*bkg )")

ws.pdf("pdf").fitTo(ws.data("data"))

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

    pt_bins = [[12, 15],
               [15, 20]]
       
    eta_bins = [[0.0,1.4442]] # try to split material dependence
#                [1.560,2.5]]

    sieie = RooRealVar(options.showerShapeName,'#sigma_{i #eta i #eta}',.1,0,.15)
    sieie.setBins(200)
    pt    = RooRealVar(options.ptName,'Photon p_{T}',50,10,1000)
    eta   = RooRealVar(options.etaName,'Photon #eta',0.0,-3.0,3.0)
    chIso = RooRealVar(options.isoName, 'phoCHIso', 100, 0, 20)
    fakeMom = RooRealVar(options.momName, 'phohasFakemom', 0,1)
    
    #containor for the root variables
    vars = RooArgSet()
    vars.add(sieie)
    vars.add(pt)
    vars.add(eta)    
    vars.add(chIso)
    vars.add(fakeMom)

    allSig = None
    allBkg = None
    if options.weight is not None:
        weight = RooRealVar(options.weight,'',1)
        vars.add(weight)
        allSig = RooDataSet('allSig','All Signal Template Events',
                            vars,
                            RooFit.ImportFromFile(options.signalInput,
                                                  options.sigTreeName),
                            RooFit.WeightVar(options.weight)
                            )

        if options.useFakeMCBkg:
            fake = RooRealVar('Fake','',0.5,1.5)
            #RooRealVar isSidebands("isSidebands","isSidebands",0,1);
            #fake = RooRealVar(options.momName,"phohasFakemom",0,1)
            vars.add(fake)

        weight.setVal(1)

        ctauStates = RooCategory('ctauRegion','Cut Region in lifetime')
        ctauStates.defineType('phohasFakemom',0)
        #ctauStates.defineType('nonPrompt',1)
        
        allBkg = RooDataSet('allBkg','All Background Template Events',
                            RooArgSet(ctauStates,vars),
                            RooFit.ImportFromFile(options.backgroundInput,
                                                  options.bkgTreeName),
                            RooFit.WeightVar(options.weight)
                            )
    else:
        allSig = RooDataSet('allSig','All Signal Template Events',
                            vars,
                            RooFit.ImportFromFile(options.signalInput,
                                                  options.sigTreeName)
                            )
        
        allBkg = RooDataSet('allBkg','All Background Template Events',
                        vars,
                        RooFit.ImportFromFile(options.backgroundInput,
                                              options.bkgTreeName)
                        )
      
    output = TFile.Open(options.outputFile,'RECREATE')
    ws = RooWorkspace(options.outputFile[:options.outputFile.find('.root')],'Template Workspace')

    # put in the raw datasets, no cuts or manipulation
    getattr(ws,'import')(allSig,
                         RooFit.RenameVariable(options.showerShapeName,'Pho_SigmaIEtaIEta'),
                         RooFit.RenameVariable(options.ptName,'Pho_Pt'),
                         #RooFit.RenameVariable(options.isoName,'phoCHIso'),
                         RooFit.RenameVariable(options.etaName,'Pho_Eta'))
    getattr(ws,'import')(allBkg,
                         RooFit.RenameVariable(options.showerShapeName,'Pho_SigmaIEtaIEta'),
                         RooFit.RenameVariable(options.ptName,'Pho_Pt'),
                         #RooFit.RenameVariable(options.isoName,'phoCHIso'),
                         RooFit.RenameVariable(options.etaName,'Pho_Eta'))


    #loop through bins making all templates
    for etabin in eta_bins:
        for ptbin in pt_bins:
            
            if 'abs(Pho_Eta) < 1.4442':
                phoselsig = 'abs(Pho_Eta) > %f  && abs(Pho_Eta) < %f && Pho_Pt > %f && Pho_Pt < %f'%(etabin[0],
                                                                                              etabin[1],
                                                                                              ptbin[0],
                                                                                              ptbin[1])
#                if fake:
                phoselbkg = 'phoCHIso < 2  && abs(Pho_Eta) > %f  && abs(Pho_Eta) < %f && Pho_Pt > %f && Pho_Pt < %f'%(etabin[0],
                                                                                                         etabin[1],
                                                                                                         ptbin[0],
                                                                                                         ptbin[1])
                    
            postfix = '_Eta_%.4f_%.4f_Pt_%.2f_%.2f'%(etabin[0],etabin[1],ptbin[0],ptbin[1])
       
            binSig = ws.data('allSig').reduce(RooFit.Cut(phoselsig),
#.........这里部分代码省略.........