Python RooWorkspace.var方法代码示例

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
 def generate_bw(self, peak, width, low, high, nbins, nev):
     ws = RooWorkspace('ws', 'ws')
     ws.factory('BreitWigner::bw(mass[200, 600], peak[350, 200, 700], width[50])')
     ws.var('mass').setRange(low,high)
     ws.var('peak').setVal(peak)
     ws.var('width').setVal(width)
     _bw = ws.pdf('bw')
     ws.var('mass').setBins(nbins)
     _ds = _bw.generateBinned(RooArgSet(ws.var('mass')), nev)
     _h = _ds.createHistogram('hBw', ws.var('mass'))
     return _h

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
def doMCFit(dataSet, x_var, addTitlePlot=''):  
    
    cuts_str = ''
    data = dataSet.reduce( RooFit.Cut(cuts_str) )

    x=RooRealVar(x_var, 'm_{#tau}',1757,1797,'MeV')
    numBins = 100 # define here so that if I change it also the ndof change accordingly
    x.setBins(numBins)

    ######################################################
    # DEFINE PDF
    ######################################################
    
    w = RooWorkspace('w')
    getattr(w,'import')(x)
    w.factory('''RooDSCBShape::DSCB({0},
    #mu[1777, 1760,1790],
    #sigma[5,0,10],
    #alpha[1.2], n[50, 1, 150],
    #alpha, n
    )'''.format(x_var))
    #w.var('n').setConstant(False)
    signal = w.pdf('DSCB')
    # w.factory('''RooGaussian::GG({0},
    # #mu, #sigma
    # )'''.format(x_var))
    # signal = w.pdf('GG')
  
    # Fit
    fit_region = x.setRange('fit_region',1757,1797)
    result = signal.fitTo(dataSet, RooFit.Save(), RooFit.Range('fit_region'))

    # Frame
    frame = x.frame(RooFit.Title(' Combined mass KK#mu '+addTitlePlot))
    dataSet.plotOn(frame)
    signal.plotOn(frame, RooFit.LineWidth(2))

    # Legends
    signal.paramOn(frame, RooFit.Layout(0.6,0.9,0.9))
    chi2 = round(frame.chiSquare(),2)
    leg = TLegend(0.3,0,.10,.10)
    leg.SetBorderSize(0)
    leg.SetFillStyle(0)
    leg.AddEntry(0,'#chi^{2} ='+str(chi2),'')
    frame.addObject(leg)

    c1 = TCanvas('c1', 'c1')
    frame.Draw()
    c1.Update()

    for prm in ('#mu', '#sigma', '#alpha', 'n'): # TODO: automatize finding of variables from the function
        w.var(prm).setConstant()
    
    return w, c1

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
 def generate_poly(self, b, low, high, nbins, nev):
     ws = RooWorkspace('ws', 'ws')
     ws.factory('Polynomial::bgPol(mass[200, 500], b[-0.001])')
     ws.var('mass').setRange(low,high)
     ws.var('b').setVal(b)
     _pdf = ws.pdf('bgPol')
     ws.var('mass').setBins(nbins)
     _ds = _pdf.generateBinned(RooArgSet(ws.var('mass')), nev)
     _h = _ds.createHistogram('hBgpol', ws.var('mass'))
     return _h

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
def test():
    workspace = RooWorkspace("testworkspace")

    model = buildModel(workspace)
    mass = workspace.var('mass')
    # data = model.generate(RooArgSet(mass), 1000)
    data = getData(workspace, nevents = 10000,
                   sigma1 = 0.015, mean1=1.,
                   sigma2 = 0.015, mean2=1.)


    model.fitTo(data, NumCPU(3))

    mframe = mass.frame(Range(50,130))
    data.plotOn(mframe)
    model.plotOn(mframe) #, Range(50,130)) #, LineColor(ROOT.kRed), LineStyle(ROOT.kDashed))
    model.paramOn(mframe,
                  Format('NEU', AutoPrecision(2) ),
                  Layout(.55, 0.92, 0.92) )

    mframe.Draw()
    return workspace

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
compPdfs = []
for comp in components:
    compPdf = fitter_mWW.ws.pdf(comp)
    getattr(combinedWS, 'import')(compPdf)
    newPdf = combinedWS.pdf(comp)
    newPdf.SetName('%s_%s' % (comp, var_name))
    norm = combinedWS.function('f_%s_norm' % comp)
    combinedWS.factory('PROD::%s(%s_%s, %s)' % (comp, comp, other_var,
                                                newPdf.GetName()))
    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')

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

#.........这里部分代码省略.........
            process_lines['rate']   += ( delimiter + str(float(self.d_input[p_name]['rate']) * self.lumi_scaling) )
            process_lines['sys']     = "#systematics line: not implemented yet!!!"


        return process_lines


    def _get_observation(self):
        """
        Read the data from trees and applies a cut.
        So far, we only get rate directly as a number.
        """
        return self.d_input['observation']['rate']

    def _get_systematics_lines(self):
        """
        Find systematics and construct a table/dict
        """
        systematics_lines_list = []
        sys_dict = self.d_input['systematics']
        #loop on keys, i.e. sys names and append value if process found, otherwise, append '-'
        for sys_id in sys_dict.keys():
            values = []
            for sig_id in self.signal_process_list:
                try:
                    value = sys_dict[sys_id][sig_id]
                except KeyError:
                    value = '-'
                values.append(str(value))

            for bkg_id in self.bkg_process_list:
                try:
                    value = sys_dict[sys_id][bkg_id]
                except KeyError:
                    value = '-'
                values.append(str(value))

            if sys_dict[sys_id]['type'].startswith('param'): values=[]
            systematics_lines_list.append('{0} {1} {2}'.format(sys_id, sys_dict[sys_id]['type'],string.join(values,' ') ))

            self.log.debug('Systematic line: {0} '.format(systematics_lines_list[-1]))        #show the last one


        systematics_lines = ''
        n_systematics = 0
        for line in systematics_lines_list:
            systematics_lines += line
            systematics_lines += '\n'
            n_systematics += 1
        return (n_systematics, systematics_lines)


    def make_workspace(self):
        """Make RooWorkspace and dump to a file"""

        gSystem.AddIncludePath("-I$CMSSW_BASE/src/ ");
        gSystem.Load("$CMSSW_BASE/lib/slc5_amd64_gcc472/libHiggsAnalysisCombinedLimit.so");
        gSystem.AddIncludePath("-I$ROOFITSYS/include");



        self.w = RooWorkspace('w')
        #run all functions_and_definitions:
        for factory_statement in self.d_input['functions_and_definitions']:
            self.w.factory(factory_statement)

        for p in self.process_list:
            self.log.debug('Checking for shape in {0}/{1}'.format(self.datacard_name, p))
            try:
                self.d_input[p]['shape']
            except KeyError:
                pass
            else:
                if self.d_input[p]['shape']:
                    self.shapes_exist = True
                    self.w.factory(self.d_input[p]['shape'])
        self.log.debug('Printing workspace...')

        self.data_obs = self.w.pdf('ggH').generate(RooArgSet(self.w.var('mass4l')), self._get_observation())

        self.data_obs.SetNameTitle('data_obs','data_obs')
        getattr(self.w,'import')(self.data_obs)
        if self.DEBUG:
            print 20*"----"
            self.w.Print()
            print 20*"----"
        self.w.writeToFile(self.shapes_output_file)
        self.log.debug('Datacard workspace written to {0}'.format(self.shapes_output_file))


    def scale_lumi_by(self, lumi_scaling):
        """
        Scales luminosity in datacards by a fixed factor. This can be
        used to get exclusion limits projections with higher luminosities.
        """
        self.lumi_scaling = lumi_scaling
        if self.lumi_scaling != 1.0:
            self.card_header+='Rates in datacard are scaled by a factor of {0}'.format(self.lumi_scaling)

        self.log.debug('Rates in datacards will be scaled by a factor of {0}'.format(self.lumi_scaling))

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
def main():
    # usage description
    usage = "Example: ./scripts/createDatacards.py --inputData inputs/rawhistV7_Run2015D_scoutingPFHT_UNBLINDED_649_838_JEC_HLTplusV7_Mjj_cor_smooth.root --dataHistname mjj_mjjcor_gev --inputSig inputs/ResonanceShapes_gg_13TeV_Scouting_Spring15.root -f gg -o datacards -l 1866 --lumiUnc 0.027 --massrange 1000 1500 50 --runFit --p1 5 --p2 7 --p3 0.4 --massMin 838 --massMax 2037 --fitStrategy 2"

    # input parameters
    parser = ArgumentParser(description='Script that creates combine datacards and corresponding RooFit workspaces',epilog=usage)

    parser.add_argument("--inputData", dest="inputData", required=True,
                        help="Input data spectrum",
                        metavar="INPUT_DATA")

    parser.add_argument("--dataHistname", dest="dataHistname", required=True,
                        help="Data histogram name",
                        metavar="DATA_HISTNAME")

    parser.add_argument("--inputSig", dest="inputSig", required=True,
                        help="Input signal shapes",
                        metavar="INPUT_SIGNAL")

    parser.add_argument("-f", "--final_state", dest="final_state", required=True,
                        help="Final state (e.g. qq, qg, gg)",
                        metavar="FINAL_STATE")

    parser.add_argument("-f2", "--type", dest="atype", required=True, help="Type (e.g. hG, lG, hR, lR)")

    parser.add_argument("-o", "--output_path", dest="output_path", required=True,
                        help="Output path where datacards and workspaces will be stored",
                        metavar="OUTPUT_PATH")

    parser.add_argument("-l", "--lumi", dest="lumi", required=True,
                        default=1000., type=float,
                        help="Integrated luminosity in pb-1 (default: %(default).1f)",
                        metavar="LUMI")

    parser.add_argument("--massMin", dest="massMin",
                        default=500, type=int,
                        help="Lower bound of the mass range used for fitting (default: %(default)s)",
                        metavar="MASS_MIN")

    parser.add_argument("--massMax", dest="massMax",
                        default=1200, type=int,
                        help="Upper bound of the mass range used for fitting (default: %(default)s)",
                        metavar="MASS_MAX")

    parser.add_argument("--p1", dest="p1",
                        default=5.0000e-03, type=float,
                        help="Fit function p1 parameter (default: %(default)e)",
                        metavar="P1")

    parser.add_argument("--p2", dest="p2",
                        default=9.1000e+00, type=float,
                        help="Fit function p2 parameter (default: %(default)e)",
                        metavar="P2")

    parser.add_argument("--p3", dest="p3",
                        default=5.0000e-01, type=float,
                        help="Fit function p3 parameter (default: %(default)e)",
                        metavar="P3")

    parser.add_argument("--lumiUnc", dest="lumiUnc",
                        required=True, type=float,
                        help="Relative uncertainty in the integrated luminosity",
                        metavar="LUMI_UNC")

    parser.add_argument("--jesUnc", dest="jesUnc",
                        type=float,
                        help="Relative uncertainty in the jet energy scale",
                        metavar="JES_UNC")

    parser.add_argument("--jerUnc", dest="jerUnc",
                        type=float,
                        help="Relative uncertainty in the jet energy resolution",
                        metavar="JER_UNC")

    parser.add_argument("--sqrtS", dest="sqrtS",
                        default=13000., type=float,
                        help="Collision center-of-mass energy (default: %(default).1f)",
                        metavar="SQRTS")

    parser.add_argument("--fixP3", dest="fixP3", default=False, action="store_true", help="Fix the fit function p3 parameter")

    parser.add_argument("--runFit", dest="runFit", default=False, action="store_true", help="Run the fit")

    parser.add_argument("--fitBonly", dest="fitBonly", default=False, action="store_true", help="Run B-only fit")

    parser.add_argument("--fixBkg", dest="fixBkg", default=False, action="store_true", help="Fix all background parameters")

    parser.add_argument("--decoBkg", dest="decoBkg", default=False, action="store_true", help="Decorrelate background parameters")

    parser.add_argument("--fitStrategy", dest="fitStrategy", type=int, default=1, help="Fit strategy (default: %(default).1f)")

    parser.add_argument("--debug", dest="debug", default=False, action="store_true", help="Debug printout")

    parser.add_argument("--postfix", dest="postfix", default='', help="Postfix for the output file names (default: %(default)s)")

    parser.add_argument("--pyes", dest="pyes", default=False, action="store_true", help="Make files for plots")

    parser.add_argument("--jyes", dest="jyes", default=False, action="store_true", help="Make files for JES/JER plots")

    parser.add_argument("--pdir", dest="pdir", default='testarea', help="Name a directory for the plots (default: %(default)s)")
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
## cuts = '(muPlusPt > 3.5) && (muMinusPt > 3.5) && (abs(upsRapidity)<2.4)'
## simparamfile = 'nom3.5SimFit.txt'

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)

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
mmin = 7.
mmax = 14.


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

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

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)

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
    dataTree.Draw('{0}>>pedhist({1},{2:0.1f},{3:0.1f}'.format(HOTower,
                                                              int(maxPed-minPed),
                                                              minPed,maxPed),
                  pedCut, 'goff')

    pedhist = gDirectory.Get('pedhist')
    ws = RooWorkspace('ws')
    x = RooRealVar('x', 'energy', minPed, maxPed, 'fC')
    x.Print()

    ds = fillDataSet(dataTree.GetV1(), x, dataTree.GetSelectedRows())
    getattr(ws, 'import')(ds)
    findOnePe(pedhist, ws, Npe=opts.npe)

    pedPlusOne = ws.pdf('pedPlusOne')
    peMean = ws.var('peMean')
    pedMean = ws.var('pedMean')
    ## pedhist.Draw()
    ## onePeF.Draw('same')

    xf = x.frame(x.getMin(), x.getMax(), int(x.getMax()-x.getMin()))
    ds.plotOn(xf)
    pedPlusOne.plotOn(xf) #,
                      #RooFit.Range('pedPlus_fit'))
    pedPlusOne.plotOn(xf, #RooFit.Range('pedPlus_fit'),
                      RooFit.Components("bg*"),
                      RooFit.LineColor(kRed),
                      RooFit.LineStyle(kDashed))
    xf.Draw()
    gPad.Update()

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
w.factory("{0}[0,0.18]".format(massVar))
w.factory("{0}[-100,100]".format(vtxVar))
w.factory("uncP[0,10]")
w.factory("cut[0,1]")

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

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

w.factory("Gaussian::vtx_model({0},mean[-50,50],sigma[0,50])".format(vtxVar))
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))',{0},gauss_mean[0,-20,20],gauss_sigma[5,1,50],exp_breakpoint[10,0,50],exp_length[3,0.5,20])".format(vtxVar))
gaussexp_pdf = w.pdf("gaussExp")
w.defineSet("obs_1d",vtxVar)
obs=w.set("obs_1d")
vtxZRefitUnc = w.var(vtxVar)
vtxZRefitUnc.setBins(200)
gauss_params = gauss_pdf.getParameters(obs)
gaussexp_params = gaussexp_pdf.getParameters(obs)

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

w.defineSet("poi","strength")
w.factory("nbkg[0,1e6]")
nbins_histpdf = 200
bkgHist = TH1D("bkgHist","bkgHist",nbins_histpdf,-100,100)
sigHist = TH1D("sigHist","sigHist",nbins_histpdf,-100,100)

pdf = w.pdf("model")

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [as 别名]
dataset = RooDataSet("data","data",events,w.set("myVars"),"cut==1")

#frame=w.var("uncVZ").frame()
#c.SetLogy()
#dataset.plotOn(frame)
#frame.Draw()
#c.SaveAs("test.png")
#binning=3
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(1000)
gauss_params = gauss_pdf.getParameters(obs)
gaussexp_params = gaussexp_pdf.getParameters(obs)

for i in range(1,n_massbins):
#for i in range(0,totalH.GetXaxis().GetNbins()-binning+2,binning):
    print i
    centermass=yieldhist.GetXaxis().GetBinCenter(i)
    massrange=2.5*(massres_a*centermass + massres_b)
    lowedge = centermass-massrange/2.0
    highedge = centermass+massrange/2.0
    name="Radiative vertex Z, mass [{}, {}] GeV".format(lowedge,highedge)
    #dataset = RooDataSet("data","data",events,w.set("myVars"),"abs(uncM-{0})<{1}".format(centermass,massrange/2))
    dataInRange = dataset.reduce(obs,"abs(uncM-{0})<{1}".format(centermass,massrange/2))
    if dataInRange.sumEntries()<100:

# 需要导入模块: from ROOT import RooWorkspace [as 别名]
# 或者: from ROOT.RooWorkspace import var [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 var [as 别名]
#!/usr/bin/env python

from ROOT import TFile, TCanvas
from ROOT import RooDataSet, RooWorkspace, RooArgSet

tf = TFile.Open('AnalysisOut.root')
tree = tf.Get('AnalysisTree')

ws = RooWorkspace("w","w")
observables = RooArgSet()
ws.defineSet("observables",observables)

ws.factory("mass[5050,6000]")
getattr(ws,'set')("observables").add(ws.var("mass"))
ws.factory("gamgams_pt[0,40e3]")
getattr(ws,'set')("observables").add(ws.var("gamgams_pt"))


mc = RooDataSet('mc','',getattr(ws,'set')('observables'))
data = RooDataSet('data','',getattr(ws,'set')('observables'))

for ev in range(tree.GetEntries()):
  tree.GetEntry(ev)
  if tree.itype != -88 and tree.itype!=72 and tree.itype!=82: continue
  if tree.bdtoutput<0.2: continue
  if tree.B0_MM < ws.var("mass").getMin() or tree.B0_MM > ws.var("mass").getMax(): continue

  ws.var("mass").setVal(tree.B0_MM)
  ws.var("gamgams_pt").setVal(tree.gamgams_PT)

  if tree.itype == -88:

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