Python ROOT.RooWorkspace类代码示例

def main(options,args):
    gROOT.Reset()

    #load our super special Polarization PDF
    gROOT.ProcessLine('.L RooPolarizationPdf.cxx+')
    gROOT.ProcessLine('.L RooPolarizationConstraint.cxx+')
    
    #setup integration
    intConf = ROOT.RooAbsReal.defaultIntegratorConfig()
    #intConf.Print('v')
#    intConf.method1D().setLabel('RooAdaptiveGaussKronrodIntegrator1D')
    intConf.setEpsAbs(1e-13)
    intConf.setEpsRel(1e-13)
    print intConf.epsAbs()
    print intConf.epsRel()
#    intConf.method2D().setLabel('RooMCIntegrator')
#    intConf.methodND().setLabel('RooMCIntegrator')


    output = TFile.Open(options.workspaceName+'.root','RECREATE')

    print options.gen

	#BfracNow = jpsi.pTRange[1][1][1]

	#p_bin = 1

	#sigMaxMass = jpsi.polMassJpsi[p_bin]

    theWS = RooWorkspace(options.workspaceName,1)

    #save the polarization PDF code in the RooWorkspace
    theWS.importClassCode('RooPolarization*',True)

    buildDataAndCategories(theWS,options,args)    

    buildMassAndLifetimePDF(theWS)

#    if options.fitFrame is not None:
#        buildPolarizationPDF(theWS,options)

    #root is stupid
    output.cd()
    
    theWS.Print('v')


    ROOT.RooMsgService.instance().Print()

    doFit(theWS,options)
    
    theWS.Print('v')
    
    theWS.Write()
    output.Close()

    def update(self,wsname,fname,storelist=[]):
        """Stores the current model (and observable) into a ROOT.Workspace,
        persistified in a root file
        
        Parameters
        ----------
        wsname: str
             Name to be given to the ROOT.Workspace
        fname:  str
             Name to be given to the ROOT file
        storelist: list(ROOT.RooDataSet,ROOT.RooRealVar,...), optional
             List of RooFit objects to be stored in the same workspace
        """
        import os
        from ROOT import RooWorkspace, TFile

        # Create the ws with the model and data available
        w = RooWorkspace(wsname,'Workspace')
        wsImport = getattr(w,'import')
        # Put the models 
        for rawtuple in self.__models.values():
            model = rawtuple[0]
            anythingelse = rawtuple[1:]
            wsImport(model)
        # Put whatever the user want to store
        for _item in storelist:
            wsImport(_item)

        # Check if the fname already exist
        file_exist = os.path.isfile(fname)
        if file_exist:
            # Create a new auxiliar file, to provisionaly
            # store the ws 
            auxfname = '_auxroot.root'
            w.writeToFile(auxfname)
            # Previously delete the old Workspace, if any
            _rootfile = TFile(fname)
            _rootfile.Delete(wsname+';*')
            _rootfile.Delete("ProcessID*;*")
            _rootfile.Close()
            # Open the auxiliary file
            _aux = TFile(auxfname)
            # and copy the ws to the rootfile
            w = _aux.Get(wsname)
        # and copy the ws to the rootfile
        w.writeToFile(fname,False)

        if file_exist:
            # And closing and removing auxfile
            _aux.Close()
            os.remove(auxfname)

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

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

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

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

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

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

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

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

                del gaus_nll   
                del truth_pow_sigm
            
        getattr(ws,'import')(biasData)

# if being executed run bias study
if __name__ == '__main__':
    ntoys = int(sys.argv[1])    
    category = int(sys.argv[2])
    mass = float(sys.argv[3])
    channel = sys.argv[4]
    order = int(sys.argv[5])
    turnon = sys.argv[6] #fitted turn on type!!!
    truth = sys.argv[7] #truth model type!!!
    
    bs = RooWorkspace('bias_study')

    bs.factory("procWeight[0]")
    bs.factory("puWeight[0]")
    bs.factory("weight[0]")
    bs.factory("Mzg[100,180]")
    bs.var("Mzg").setRange("ROI",mass-1.5,mass+1.5)
    bs.var("Mzg").setBins(40000,"cache")
    bs.factory("Mz[0]")
    #bs.factory("dMzg[0,25]")
    #bs.factory("dMz[0,25]")
    bs.factory("r94cat[cat1=1,cat2=2,cat3=3,cat4=4]")
    bs.defineSet("observables",
                 "Mzg,Mz,r94cat,procWeight,puWeight")
    bs.defineSet("observables_weight",
                 "Mzg,Mz,r94cat,procWeight,puWeight,weight")

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

gROOT.LoadMacro("tdrstyle.C")
gROOT.ProcessLine("setTDRStyle()")
#TH1.SetDefaultSumw2()

gStyle.SetMarkerSize(0.7)
kBlue = 600

latex = TLatex()
latex.SetNDC()
latex.SetTextAlign(12)
latex.SetTextFont(42)
latex.SetTextSize(0.035)

ws = RooWorkspace("ws")
x = RooRealVar("x", "p_{T}^{reco}/p_{T}^{parton}", xlow,xup)
obs = RooArgList(x)
c1 = TCanvas("c1", "c1",600,600) ;
if logy:  c1.SetLogy()


# Loop over input files (currently disabled)
infile = TFile.Open(infilename)
for ih in xrange(ih_,ih_+1):
    h = getattr(infile, "hh1_%i" % ih)
    data = RooDataHist("data_%i" %ih, "", obs, h)
    getattr(ws, 'import')(data)

    hl = getattr(infile, "hl_%i" % ih)
    

setattr(RooWorkspace, "Import", getattr(RooWorkspace, "import"))

## Here starts the meat.

nentries = -1
## sTest = [-20, -10, -5, -2, -1, -0.5, 0, 0.5, 1, 2, 5, 10, 20]
## sTest = [-5, -2, -1, -0.5, 0, 0.5, 1, 2, 5]
sTest = [0]
phoPtRange = (12,15)

# chains = getChains('v11')
chains = pmvTrees.getChains('v15')
mcTree = chains['z']
dataTree = chains['data']

w = RooWorkspace('w')

## Define variables 
mmgMass = w.factory('mmgMass[40, 140]')
mmMass = w.factory('mmMass[10, 140]')
weight = w.factory('weight[1]')
phoScale = w.factory('phoScale[0,-50,50]')
weight.SetTitle('pileup.weight')

phoScale.setUnit('%')

## Photon scaling fraction, dlog(m_uuy)/dlog(E_y)
fPho = w.factory('fPho[0.15*91.2,0,100]')
fPhoFunc = w.factory('''FormulaVar::fPhoFunc(
    "mmgMass * (0.5 - 0.5 * mmMass^2 / mmgMass^2)",
    {mmMass, mmgMass}

    minPed = TMath.MinElement(dataTree.GetSelectedRows(), dataTree.GetV1())
    maxPed = TMath.MaxElement(dataTree.GetSelectedRows(), dataTree.GetV1())
    #print 'minPed:',minPed,'maxPed:',maxPed
    minPed = int(minPed) - 2.5
    maxPed = int(maxPed) + 2.5
    print 'minPed:',minPed,'maxPed:',maxPed
    while int((maxPed-minPed)/2.)*2. < (maxPed-minPed):
        maxPed += 1.
    print 'minPed:',minPed,'maxPed:',maxPed
    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)

RooAbsReal.defaultIntegratorConfig().getConfigSection(
    'RooAdaptiveGaussKronrodIntegrator1D').setRealValue('maxSeg', 1000)
RooAbsReal.defaultIntegratorConfig().method1D().setLabel(
    'RooAdaptiveGaussKronrodIntegrator1D')
RooAbsReal.defaultIntegratorConfig().method1DOpen().setLabel(
    'RooAdaptiveGaussKronrodIntegrator1D')

# seed the Random number generator
rndm = TRandom3(SEED + 1)
RooRandom.randomGenerator().SetSeed(int(rndm.Uniform(4294967295)))
del rndm

# start building the fit
from B2DXFitters.WS import WS

ws = RooWorkspace('ws')
one = WS(ws, RooConstVar('one', '1', 1.0))
zero = WS(ws, RooConstVar('zero', '0', 0.0))

# start by defining observables
time = WS(ws, RooRealVar('time', 'time [ps]', 0.2, 15.0))
qf = WS(ws, RooCategory('qf', 'final state charge'))
qf.defineType('h+', +1)
qf.defineType('h-', -1)
qt = WS(ws, RooCategory('qt', 'tagging decision'))
qt.defineType(      'B+', +1)
qt.defineType('Untagged',  0)
qt.defineType(      'B-', -1)

# now other settings
Gamma  = WS(ws, RooRealVar( 'Gamma',  'Gamma',  0.661)) # ps^-1

def setup_workspace(config):

  import ROOT
  from ROOT import RooWorkspace, gROOT, gStyle, RooAbsReal, RooMsgService, RooFit
  #from ROOT import RooFit, gROOT, gDirectory, gStyle, gPad, TTree, RooCmdArg,RooBinning
  #from ROOT import RooRealVar, RooMappedCategory, RooCategory, RooFormulaVar, RooAbsData
  #from ROOT import RooBMixDecay, RooMCStudy, RooAddModel, RooEffProd, RooMsgService
  #from ROOT import RooWorkspace, TCanvas, TFile, kFALSE, kTRUE, RooDataSet, TStopwatch
  #from ROOT import RooArgSet, RooArgList, RooRandom, RooMinuit, RooAbsReal, RooDataHist
  #from ROOT import TBrowser, TH2F, TF1, TH1F, RooGenericPdf, RooLinkedList
  from math import sqrt

  gROOT.SetStyle("Plain")
  gStyle.SetPalette(1)
  gStyle.SetOptStat(0)
  gStyle.SetOptFit(0)
  gStyle.SetOptStat(1111)
  gStyle.SetOptFit(10111)
  gStyle.SetOptTitle(1)

  #gROOT.ProcessLine(".L RooGaussianTrunk.cxx+")
  #gROOT.ProcessLine(".L RooCBShapeTrunk.cxx+")
  #gROOT.ProcessLine(".L RooChebychevTrunk.cxx+")
  #from ROOT import RooGaussianTrunk, RooChebychevTrunk, RooCBShapeTrunk

  #RooAbsReal.defaultIntegratorConfig().Print()
  RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-8)
  RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-8)
  #RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-6)
  #RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-6)
  RooAbsReal.defaultIntegratorConfig().Print()

  print "Numeric integration set up" #TODO: is the integration acceptable?

  ##This controls the logging output from RooFit
  #RooMsgService.instance().addStream(RooFit.DEBUG,RooFit.Topic(RooFit.Fitting))
  RooMsgService.instance().deleteStream(1)
  #RooMsgService.instance().addStream(RooFit.INFO,RooFit.Topic(RooFit.Generation + RooFit.Minization + RooFit.Plotting + RooFit.Fitting + RooFit.Integration + RooFit.LinkStateMgmt + RooFit.Eval + RooFit.Caching + RooFit.Optimization + RooFit.ObjectHandling + RooFit.InputArguments + RooFit.Tracing + RooFit.Contents + RooFit.DataHandling + RooFit.NumericIntegration))
  RooMsgService.instance().addStream(RooFit.INFO,RooFit.Topic(RooFit.LinkStateMgmt + RooFit.Caching + RooFit.ObjectHandling + RooFit.InputArguments + RooFit.Tracing))
  RooMsgService.instance().Print()


  print "Message service set up"

  w = RooWorkspace("w",False)

  w.factory("RAND[0,1]")
  
  if "norm" not in config["mode"]:
    D0_Mass = w.factory("D0_Mass[1815,1915]") 
  else:
    D0_Mass = w.factory("D0_Mass[1800,1930]")
  D0_Mass.setUnit("MeV")
  D0_Mass.setBins(60)
  
  Del_Mass = w.factory("Del_Mass[139,155]")
  Del_Mass.setUnit("MeV")
  Del_Mass.setBins(60)

  if "norm" not in config["mode"]:
    Dataset = w.factory("DataSet[BDT1,BDT2,BDT3]")
  else:
    Dataset = w.factory("DataSet[Norm]")

  w.factory("classID[Sig=0,Bkg=1]")
  w.factory("BDT_ada[-1,1]")
  w.factory("x1_PIDe[-2,20]")
  w.factory("x2_ProbNNmu[0,1]")

  #D0_Mass.setRange("blinded",1700.,1900.)
  
  if "norm" not in config["mode"]:
    dataCats = ["", "BDT1", "BDT2", "BDT3"]
  else:
    dataCats = ["", "Norm"]

  for data in dataCats:
    for dst_side in ["", "delsig", "delhigh", "dellow"]:
      for d_side in ["", "dsig", "dhigh", "dlow", "dhigh1", "dlow1", "dhigh2", "dlow2"]:
        name = data+dst_side+d_side

        if data == "BDT1":
          Dataset.setRange(name,"BDT1")
        elif data == "BDT2":
          Dataset.setRange(name,"BDT2")
        elif data == "BDT3":
          Dataset.setRange(name,"BDT3")
        elif data == "Norm":
          Dataset.setRange(name,"Norm")

        if dst_side == "delhigh":
          Del_Mass.setRange(name,148.,155.)
        elif dst_side == "delsig":
          Del_Mass.setRange(name,143.,148.)
        elif dst_side == "dellow":
          Del_Mass.setRange(name,139.,143.)

        if d_side == "dhigh2":
          D0_Mass.setRange(name,1910.,1930.)
        elif d_side == "dhigh1":
#.........这里部分代码省略.........