Python RooRealVar.setConstant方法代码示例

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
def buildPdf(ws, p):

    mass = RooRealVar("mass", "mass", p.minMass, p.maxMass)
    getattr(ws,'import')(mass)

    # Construct signal pdf
    mean = RooRealVar("mean", "mean", 90, 85, 95)
    width = RooRealVar("width", "width", 2.4952, 1, 3)
    width.setConstant(ROOT.kTRUE)
    sigma = RooRealVar("sigma", "sigma", 1.2, 0.2, 10)
    signalAll = RooVoigtian("signalAll", "signalAll", mass, mean, width, sigma)

    turnOnAll = RooRealVar("turnOnAll","turnOnAll", 80., 40., 150.)
    widthAll_bkg = RooRealVar("widthAll","widthAll", 2., 0., 50.)
    decayAll_bkg = RooRealVar("decayAll","decayAll", 80., 20., 150.)
    meanB = RooRealVar("meanB", "meanB", 90, 60, 130)
    sigmaB = RooRealVar("sigmaB", "sigmaB", 10, 1, 20)
    bkg_a1 = RooRealVar("bkg_a1", "bkg_a1", 0., -2., 2.)
    bkg_a2 = RooRealVar("bkg_a2", "bkg_a2", 0., -2., 2.)
    backgroundAll = RooGaussian("backgroundAll", "backgroundAll", mass, meanB, sigmaB)

    # Construct composite pdf
    sigAll = RooRealVar("sigAll", "sigAll", 2000, 0, 100000)
    bkgAll = RooRealVar("bkgAll", "bkgAll", 100, 0, 10000)
    modelAll = RooAddPdf("modelAll", "modelAll", RooArgList(signalAll, backgroundAll), RooArgList(sigAll, bkgAll))
    if p.NoBkgd:
        modelAll = RooAddPdf("modelAll", "modelAll", RooArgList(signalAll), RooArgList(sigAll))
    # Define pdf for all probes

    # Construct signal pdf.
    # NOTE that sigma is shared with the signal sample model
    signalPass = RooVoigtian("signalPass","signalPass",mass,mean,width,sigma)
    # Construct the background pdf
    backgroundPass = RooGaussian("backgroundPass", "backgroundPass", mass, meanB, sigmaB)

    # Construct the composite model
    efficiency = RooRealVar("efficiency","efficiency",0.9,0.3,1.)
    sigPass = RooFormulaVar("sigPass", "@0*@1", RooArgList(sigAll, efficiency))
    bkgPass = RooRealVar("bkgPass", "bkgPass", 100, 0, 10000)
    modelPass = RooAddPdf("modelPass", "modelPass", RooArgList(signalPass, backgroundPass), RooArgList(sigPass, bkgPass))
    if p.NoBkgd:
        modelPass = RooAddPdf("modelPass", "modelPass", RooArgList(signalPass), RooArgList(sigPass))

    frac = RooRealVar("frac", "frac", 0.8, 0., 1.)

    # Define combined pdf for simultaneous fit

    # Define category to distinguish physics and control samples events
    sample = RooCategory("sample","sample")
    sample.defineType("all")
    sample.defineType("pass")

    simPdf = RooSimultaneous("simPdf","simultaneous pdf",sample)

    # Associate model with the physics state and model_ctl with the control state
    simPdf.addPdf(modelAll,"all")
    simPdf.addPdf(modelPass,"pass")
    # ws.import(simPdf)
    getattr(ws,'import')(simPdf)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
def main():
    # Mjj0 of TT MC Bkg
    f1 = TFile("Merged_TT_TuneCUETP8M1_13TeV-powheg-pythia8-runallAnalysis.root")
    h_Mjj = f1.Get("histfacFatJet_ZLight/h_Mjj0")
    h_Mjj.GetXaxis().SetRangeUser(0,300)
    var_mean = h_Mjj.GetMean()

    # Build Gaussian PDF
    x     = RooRealVar(  'x',     'x',                0, 300 )
    mean  = RooRealVar(  'mean',  'mean of gaussian', var_mean )
    sigma = RooRealVar(  'sigma', 'width of gaussian', 5)
    gauss = RooGaussian( 'gauss', 'gaussian PDF', x, mean, sigma)
    data  = RooDataHist("data","Mjj dataset",RooArgList(x), h_Mjj);
    
    # Plot PDF
    xframe = x.frame(RooFit.Title("Gaussian p.d.f."))
    gauss.plotOn( xframe )
    gauss.plotOn(xframe,RooFit.LineColor(2)) 
    
    # Generate a toy MC set
    # data = gauss.generate( RooArgSet(x), 10000 )
    # Plot PDF and toy data overlaid
    xframe2 = x.frame(RooFit.Title("Gaussian p.d.f. with Mjj"))
    # data.plotOn( xframe2, RooLinkedList() )
    # data.plotOn( xframe2 )
    data.plotOn( xframe2 )
    gauss.plotOn( xframe2)
    # Fit PDF to toy
    mean.setConstant( kFALSE )
    sigma.setConstant( kFALSE )
    gauss.fitTo(data)
    
    c1 = TCanvas("c1","Example",800,400)
    c1.Divide(3)
    c1.cd(1)
    gPad.SetLeftMargin(0.15)
    xframe.GetYaxis().SetTitleOffset(1.6)
    xframe.Draw()
    
    c1.cd(2)
    gPad.SetLeftMargin(0.15)
    xframe2.GetYaxis().SetTitleOffset(1.6)
    xframe2.Draw() 

    c1.SaveAs('testMjj0.png')
    
    # # Print final value of parameters
    fout = TFile("output.root","recreate")
    c1.Write()
    fout.Close()

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

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
def fitPed(hist, ws, name='x'):
    maxBin = hist.GetMaximumBin()
    x = ws.var(name)
    #rds = ds.reduce('{1}<{0:0.2f}'.format(hist.GetBinLowEdge(maxBin+2),name))
    #rds.Print()
    x.setRange('ped_fit', x.getMin(), hist.GetBinLowEdge(maxBin+3))
    pedMean = RooRealVar('pedMean', 'mean_{ped}', hist.GetBinCenter(maxBin),
                         x.getMin(), x.getMax())
    pedMean.Print()
    pedWidth = RooRealVar('pedWidth', 'sigma_{ped}', 1., 0., 10.)
    pedWidth.Print()
    ped = RooGaussian('ped', 'ped', x, pedMean, pedWidth)

    pedMean.setConstant(False)
    ped.fitTo(ws.data('ds'), RooFit.Minos(False), RooFit.Range('ped_fit'),
              RooFit.PrintLevel(0))

    getattr(ws, 'import')(ped)

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

#.........这里部分代码省略.........
    masses.sort()

    # import ROOT stuff
    from ROOT import TFile, TH1F, TH1D, kTRUE, kFALSE
    from ROOT import RooRealVar, RooDataHist, RooArgList, RooArgSet, RooAddPdf, RooFit, RooGenericPdf, RooWorkspace, RooMsgService, RooHistPdf

    if not args.debug:
        RooMsgService.instance().setSilentMode(kTRUE)
        RooMsgService.instance().setStreamStatus(0,kFALSE)
        RooMsgService.instance().setStreamStatus(1,kFALSE)

    # input data file
    inputData = TFile(args.inputData)
    # input data histogram
    hData = inputData.Get(args.dataHistname)

    # input sig file
    inputSig = TFile(args.inputSig)

    sqrtS = args.sqrtS

    for mass in masses:

        print ">> Creating datacard and workspace for %s resonance with m = %i GeV..."%(args.final_state, int(mass))

        hSig = inputSig.Get( "h_" + args.final_state + "_" + str(int(mass)) )

        # calculate acceptance of the dijet mass cut
        sigAcc = hSig.Integral(hSig.GetXaxis().FindBin(args.massMin),hSig.GetXaxis().FindBin(args.massMax))/hSig.Integral(1,hSig.GetXaxis().FindBin(args.massMax))

        mjj = RooRealVar('mjj','mjj',args.massMin,args.massMax)

        rooSigHist = RooDataHist('rooSigHist','rooSigHist',RooArgList(mjj),hSig)
        rooSigHist.Print()
        signal = RooHistPdf('signal','signal',RooArgSet(mjj),rooSigHist)
        signal.Print()
        signal_norm = RooRealVar('signal_norm','signal_norm',0,-1e+04,1e+04)
        if args.fitBonly: signal_norm.setConstant()
        signal_norm.Print()

        p1 = RooRealVar('p1','p1',args.p1,0.,100.)
        p2 = RooRealVar('p2','p2',args.p2,0.,60.)
        p3 = RooRealVar('p3','p3',args.p3,-10.,10.)

        background = RooGenericPdf('background','(pow([email protected]/%.1f,@1)/pow(@0/%.1f,@[email protected]*log(@0/%.1f)))'%(sqrtS,sqrtS,sqrtS),RooArgList(mjj,p1,p2,p3))
        background.Print()
        dataInt = hData.Integral(hData.GetXaxis().FindBin(args.massMin),hData.GetXaxis().FindBin(args.massMax))
        background_norm = RooRealVar('background_norm','background_norm',dataInt,0.,1e+07)
        background_norm.Print()

        # S+B model
        model = RooAddPdf("model","s+b",RooArgList(background,signal),RooArgList(background_norm,signal_norm))

        rooDataHist = RooDataHist('rooDatahist','rooDathist',RooArgList(mjj),hData)
        rooDataHist.Print()

        if args.runFit:
            res = model.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            res.Print()

        dcName = 'datacard_' + args.final_state + '_m' + str(mass) + '.txt'
        wsName = 'workspace_' + args.final_state + '_m' + str(mass) + '.root'

        w = RooWorkspace('w','workspace')
        getattr(w,'import')(signal)
        getattr(w,'import')(background)
        getattr(w,'import')(background_norm)
        getattr(w,'import')(rooDataHist,RooFit.Rename("data_obs"))
        w.Print()
        w.writeToFile(os.path.join(args.output_path,wsName))

        # -----------------------------------------
        # write a datacard
        lumi = args.lumi
        signalCrossSection = 1. # set to 1. so that the limit on r can be interpreted as a limit on the signal cross section
        expectedSignalRate = signalCrossSection*lumi*sigAcc

        datacard = open(os.path.join(args.output_path,dcName),'w')
        datacard.write('imax 1\n')
        datacard.write('jmax 1\n')
        datacard.write('kmax *\n')
        datacard.write('---------------\n')
        datacard.write('shapes * * '+wsName+' 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')
        #flat parameters --- flat prior
        datacard.write('background_norm  flatParam\n')
        datacard.write('p1  flatParam\n')
        datacard.write('p2  flatParam\n')
        datacard.write('p3  flatParam\n')
        datacard.close()

    print '>> Datacards and workspaces created and stored in %s/.'%( os.path.join(os.getcwd(),args.output_path) )

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
def main():

    parser = OptionParser()
    parser.add_option("-d", "--dir", type="string", dest="outDir", metavar="DIR", default="./",
                      help="output directory for .png")
        

    (options, args) = parser.parse_args()

    if os.path.exists(options.outDir) and options.outDir!="./":
        print "Sorry, file ",options.outDir," already exist, choose another name\n"
        exit(1)
    else:
        os.system("mkdir -p "+options.outDir)


    """
    Set the style ...
    """

    myNewStyle = TStyle("myNewStyle","A better style for my plots")
    setStyle(myNewStyle)

    TH1F.SetDefaultSumw2(True)

    # Histogram range
    xlow = 70.
    xup = 110.
    
    # Mass range for fit
    minM_fit = 75.
    maxM_fit = 105.

    # Ratio plot range
    minR = 0.8
    maxR = 1.2

    # TLines for ratio plot
    line = TLine(minM_fit,1,maxM_fit,1)
    line.SetLineWidth(2)
    line.SetLineColor(2)

    # Canvas
    spectrum_height = 0.75
    ratio_correction = 1.4
    ratio_height = (1-spectrum_height)*ratio_correction
    xOffset = 0.08
    yOffset = 0.04

    cTest = TCanvas("cTest","cTest")
    
    c2 = TCanvas("c2","c2")
    c2.SetFillColor(0)
    c2.Divide(1,2)

    
    c3 = TCanvas("c3","c3")
    c3.SetFillColor(0)
    c3.Divide(1,2)

    # Files MuScleFit
    fDataBefore = TFile("h3_Z_data_beforeCorrection.root")
    fDataAfter  = TFile("h3_Z_data_afterCorrection.root")
    fMcBefore = TFile("h3_Z_mc_beforeCorrection.root")
    fMcAfter  = TFile("h3_Z_mc_afterCorrection.root")
    
                
    # Retrieve histograms and fit
    hDataBefore = fDataBefore.Get("h1_Z_data")
    hDataAfter  = fDataAfter.Get("h1_Z_data")
    hMcBefore   = fMcBefore.Get("h1_Z_mc")
    hMcAfter    = fMcAfter.Get("h1_Z_mc")

    # Identifiers
    ids = ["data_before","data_after","mc_before","mc_after"]

    # Create histograms dictionary
    histos = {}
    histos["data_before"] = hDataBefore
    histos["data_after"]  = hDataAfter
    histos["mc_before"]   = hMcBefore
    histos["mc_after"]    = hMcAfter

    histosSubtr = {}

    # Create parameters dictionary
    expPars = {}
    signalPars = {}

    for i in ids:
        # RooFit definitions
        ## RooRealVars
        x = RooRealVar("x","M_{#mu#mu} (GeV)",minM_fit,maxM_fit)
        mean = RooRealVar("mean","mean",91.19,87.,94.)
        meanCB = RooRealVar("meanCB","meanCB",0.,-10.,10.)
        meanCB.setConstant(True)
        width = RooRealVar("width","width",2.4952,2.3,2.6)
        width.setConstant(True)
        sigma = RooRealVar("sigma","sigma",1.3,0.001,3.)
        #    sigma.setConstant(True)
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
def mbc_dline_che(evtfile, mc, setMres, setGamma, setR, sp1, sp2, sp3, fa,
                  fb, setmd, setp, setxi, setN1, setN2, setNbkgd1, setNbkgd2,
                  title1, title2, epsfile, txtfile, ymin=0.5,
                  cuts=None, err_type='SYMM', test=False):
    
    from ROOT import (gROOT, RooRealVar, RooCategory, RooArgSet, RooDataSet,
                      RooFit, RooGaussian, RooArgList, RooAddPdf, RooSimultaneous,
                      RooArgusBG, RooFormulaVar, RooChebychev, RooAbsData,
                      RooDataHist, TCanvas, kRed, kBlue, kGreen, kMagenta,
                      TPaveText, RooDLineShape)
    set_root_style(stat=1, grid=0)

    mbc = RooRealVar('mbc', 'Beam constrained mass', 1.83, 1.89, 'GeV')
    ebeam = RooRealVar('ebeam','Ebeam', 1.8815, 1.892, 'GeV')

    dflav = RooCategory('dflav','D0 flavor')
    dflav.defineType('dflav',1)
    dflav.defineType('dbarflav',-1)

    if cuts != None:
        if 'kkmass' in cuts:
            kkmass = RooRealVar('kkmass', 'KK invariant mass', 0.97, 1.90, 'GeV')
            ras = RooArgSet(mbc, ebeam, kkmass, dflav)
            dataset = RooDataSet.read(evtfile, ras)
        else:
            raise NameError(cuts)

        sys.stdout.write('Using cuts: %s...' %cuts)
        dataset = dataset.reduce(cuts)
        sys.stdout.write(' selected %s events.\n' % dataset.numEntries())
    else:
        ras = RooArgSet(mbc, ebeam, dflav)
        dataset = RooDataSet.read(evtfile, ras)

    res = RooRealVar("datares", "datares", mc)
    mres = RooRealVar("mres","mres", setMres)
    gamma = RooRealVar('gamma', 'gamma', setGamma)

    r = RooRealVar('r', 'r', setR)
    sigmaE = RooRealVar("sigmaE","sigmaE", 0.0021)

    sigmap1 = RooRealVar("sigmap1","sigmap1", sp1, 0.002, 0.040)

    scalep2 = RooRealVar("scalep2","scalep2",2.00,1.500,5.500)
    scalep3 = RooRealVar("scalep3","scalep3",5.00,3.00,10.000)
   
    scalep2.setVal(sp2)
    scalep2.setConstant(1)
    scalep3.setVal(sp3)
    scalep3.setConstant(1)

    as12 = RooArgList(sigmap1,scalep2)
    sigmap2 = RooFormulaVar("sigmap2","sigma2","sigmap1*scalep2", as12)

    as123 = RooArgList(sigmap1,scalep2,scalep3)
    sigmap3 = RooFormulaVar("sigmap3","sigma3","sigmap1*scalep2*scalep3",
                            as123)
    
    md = RooRealVar("md","md", setmd,1.863,1.875)

    f2 = RooRealVar("f2","f2", fa)
    f3 = RooRealVar("f3","f3", fb)
    al23 = RooArgList(f2,f3)
    f1 = RooFormulaVar("f1","f1","1.0-f2-f3", al23)

    # Construct signal shape

    fcn1_1 = RooDLineShape("DLineshape1_1","DLineShape1_1",4,mbc,ebeam,
                           mres,gamma,r,sigmaE,sigmap1,md,res)
    fcn1_2 = RooDLineShape("DLineshape1_2","DLineShape1_2",4,mbc,ebeam,
                           mres,gamma,r,sigmaE,sigmap2,md,res)
    fcn1_3 = RooDLineShape("DLineshape1_3","DLineShape1_3",4,mbc,ebeam,
                           mres,gamma,r,sigmaE,sigmap3,md,res)

    fcn2_1 = RooDLineShape("DLineshape2_1","DLineShape2_1",4,mbc,ebeam,
                           mres,gamma,r,sigmaE,sigmap1,md,res)
    fcn2_2 = RooDLineShape("DLineshape2_2","DLineShape2_2",4,mbc,ebeam,
                           mres,gamma,r,sigmaE,sigmap2,md,res)
    fcn2_3 = RooDLineShape("DLineshape2_3","DLineShape2_3",4,mbc,ebeam,
                           mres,gamma,r,sigmaE,sigmap3,md,res)


    alf1_123 = RooArgList(fcn1_1,fcn1_2,fcn1_3)
    af12 = RooArgList(f1,f2) 
    sigpdf = RooAddPdf("signal1_3","signal1_3", alf1_123, af12)

    alf2_123 = RooArgList(fcn2_1,fcn2_2,fcn2_3)
    
    sigbarpdf = RooAddPdf("signal2_3","signal2_3", alf2_123, af12)

    con0 = RooRealVar('c0', 'constant', -1, 1)
    con1 = RooRealVar('c1', 'linear', -10, 10)
    con2 = RooRealVar('c2', 'quadratic', 1)

    bkgpdf = RooChebychev('bkgpdf', 'Background', mbc, RooArgList(con1, con2))
    
    bkgbarpdf = RooChebychev('bkgbarpdf', 'Background',
                             mbc, RooArgList(con1, con2))

    yld = RooRealVar('yld', 'D yield', 100, 0, 2000)
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
      nbkg2=RooRealVar("nbkg2"+name,"number of background events",qcdhist.Integral(),0,10*qcdhist.Integral())

      l0=RooRealVar("l0","l0",100.,0.,1000.)
      l1=RooRealVar("l1","l1",1.,0.,1000.)
      l2=RooRealVar("l2","l2",1.,0.,1000.)
      #sigbkg=RooPolynomial("sigbkg"+name,"bkg",mass,RooArgList(l0,l1,l2))
      #sigbkg=RooChebychev("sigbkg"+name,"bkg",mass,RooArgList(l0,l1))
      sigbkg=RooLogistics("sigbkg"+name,"bkg",mass,l0,l1)
      #sigbkg=RooExponential("sigbkg"+name,"sigbkg",mass,l0)
      nsigref=RooRealVar("nsigref"+name,"number of signal events",sighist.Integral(),0,10*sighist.Integral())
      nsigrefW=RooRealVar("nsigrefW"+name,"number of signal W events",sigWhist.Integral(),0,10*sigWhist.Integral())
      nsigrefZ=RooRealVar("nsigrefZ"+name,"number of signal Z events",sigZhist.Integral(),0,10*sigZhist.Integral())
      sigWfrac=RooRealVar("sigWfrac"+name,"fraction of W in signal",0.5,0,1.0)
      sig=RooAddPdf("sig"+name,"sig"+name,sigW,sigZ,sigWfrac) ;
      sigmodel=RooAddPdf("sigmodel"+name,"sig+sigbkg",RooArgList(sigbkg,sig),RooArgList(nsigbkg,nsigref))
      meanW.setConstant(False)
      meanWZ.setConstant(False) 
      widthWZ.setConstant(False) 
      sigWfrac.setConstant(False)

      sigmodelW=RooAddPdf("sigmodelW"+name,"sigW",RooArgList(sigbkg,sigW),RooArgList(nsigbkg,nsigrefW))
      sigmodelW.fitTo(signalW,RooFit.SumW2Error(True))
      sigmodelW.fitTo(signalW,RooFit.SumW2Error(True))
      sigmodelW.fitTo(signalW,RooFit.SumW2Error(True))
      xframe=mass.frame(RooFit.Title("             m="+str(int((meanW.getValV())*1000.)/1000.)+"#pm"+str(int(meanW.getError()*1000.)/1000.)+" GeV"))
      signalW.plotOn(xframe,RooFit.DataError(RooAbsData.SumW2))
      sigmodelW.plotOn(xframe,RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      canvas=TCanvas("c3","c3",0,0,600,600)
      xframe.GetYaxis().SetTitle("Events")
      xframe.Draw()
      canvas.SaveAs(prefix+"_"+plot[0]+name+"_sigWfit.pdf")

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

#.........这里部分代码省略.........
            c.Clear()
            c.SetWindowSize(1000,500)
            c.Divide(2,1)
            for i in [1,2] :
                c.cd(i)
                reg='barrel'
                if i==2: reg='endcap' 

                h=histos[r+k+reg]
                x=RooRealVar("x", h.GetXaxis().GetTitle(), h.GetXaxis().GetXmin(), h.GetXaxis().GetXmax())
                data=RooDataHist("data", "dataset with x", RooArgList(x), h)
                frame=x.frame()
                RooAbsData.plotOn( data, frame, RooFit.DataError(RooAbsData.SumW2) )

                mean1=RooRealVar("mean1","mean1",0,-0.5,0.5);
                sigma1=RooRealVar("sigma1","sigma1",0.1,0.01,1.0);
                gauss1=RooGaussian("g1","g",x,mean1,sigma1)
                
                if r=='dpt' or r=='den' :
                    mean2=RooRealVar("mean2","mean2",0,-0.5,0.5);
                    sigma2=RooRealVar("sigma2","sigma2",0.1,0.01,1.0);
                    alphacb=RooRealVar("alphacb","alphacb",1,0.1,3);
                    ncb=RooRealVar("ncb","ncb",4,1,100)
                    gauss2 = RooCBShape("cb2","cb",x,mean2,sigma2,alphacb,ncb);
                else:
                    mean1.setRange(0,0.5)
                    mean2=RooRealVar("mean2","mean",0,0,1);
                    sigma2=RooRealVar("sigma2","sigma",0.1,0.01,1.0);
                    gauss2=RooGaussian("g2","g",x,mean2,sigma2) ;

                frac = RooRealVar("frac","fraction",0.9,0.0,1.0)
                if data.sumEntries()<100 :
                    frac.setVal(1.0)
                    frac.setConstant(True)
                model = RooAddPdf("sum","g1+g2",RooArgList(gauss1,gauss2), RooArgList(frac))

                status=model.fitTo(data,RooFit.Save()).status()
                if status!=0 : continue

                model_cdf=model.createCdf(RooArgSet(x)) ;
                cl=0.90
                ul=0.5*(1.0+cl)
                closestToCL=1.0
                closestToUL=-1
                closestToMedianCL=1.0
                closestToMedian=-1
                for ibin in xrange(1,h.GetXaxis().GetNbins()*10):
                    xval=h.GetXaxis().GetXmin()+(ibin-1)*h.GetXaxis().GetBinWidth(ibin)/10.
                    x.setVal(xval)
                    cdfValToCL=math.fabs(model_cdf.getVal()-ul)
                    if cdfValToCL<closestToCL:
                        closestToCL=cdfValToCL
                        closestToUL=xval
                    cdfValToCL=math.fabs(model_cdf.getVal()-0.5)
                    if cdfValToCL<closestToMedianCL:
                        closestToMedianCL=cdfValToCL
                        closestToMedian=xval

                RooAbsPdf.plotOn(model,frame)
                frame.Draw()

                if i==1: drawHeader()
                labels.append( TPaveText(0.6,0.92,0.9,0.98,'brNDC') )
                ilab=len(labels)-1
                labels[ilab].SetName(r+k+'txt')
                labels[ilab].SetBorderSize(0)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
    nsig=RooRealVar("nsig"+str(cut),"number of signal events",100,0,1e10)
    nbkg=RooRealVar("nbkg"+str(cut),"number of background events",100,0,1e10)
    nbkg1=RooRealVar("nbkg1"+str(cut),"number of background events",100,0,1e10)
    nbkg2=RooRealVar("nbkg2"+str(cut),"number of background events",100,0,1e10)

    l0=RooRealVar("l0","l0",100.,0.,1000.)
    l1=RooRealVar("l1","l1",1.,0.,1000.)
    l2=RooRealVar("l2","l2",1.,0.,1000.)
    #sigbkg=RooPolynomial("sigbkg"+str(cut),"bkg",mass,RooArgList(l0,l1,l2))
    #sigbkg=RooChebychev("sigbkg"+str(cut),"bkg",mass,RooArgList(l0,l1))
    sigbkg=RooLogistics("sigbkg"+str(cut),"bkg",mass,l0,l1)
    #sigbkg=RooExponential("sigbkg"+str(cut),"sigbkg",mass,l0)
    nsigref=RooRealVar("nsigref","number of signal events",100,0,1e10)
    sigmodel=RooAddPdf("sigmodel"+str(cut),"sig+bkg",RooArgList(sigbkg,sig),RooArgList(nbkg,nsigref))
    s0.setConstant(False)
    s1.setConstant(False) 
    sigmodel.fitTo(signal,RooFit.SumW2Error(True)) #,RooFit.Range("signal")
    sigmodel.fitTo(signal,RooFit.SumW2Error(True)) #,RooFit.Range("signal")
    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())

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

#.........这里部分代码省略.........

        iGen_value = modelVjet_test.createIntegral(RooArgSet(J_mass), RooFit.Range("SRrange"))
        print "iGen_value:", iGen_value.getVal()
        
        iFit_value = modelVjet_test2.createIntegral(RooArgSet(J_mass), RooFit.Range("SRrange"))
        print "iFit_value:", iFit_value.getVal()


	Bias_value = ( iFit_value.getVal() - iGen_value.getVal() ) / iGen_value.getVal()
        print "Bias_value of SR:", Bias_value

#        iGen_value = modelVjet_test.createIntegral(RooArgSet(J_mass),RooFit.NormSet(RooArgSet(J_mass)), RooFit.Range("VRrange,SRrange"))
#	print "iGen_value:", iGen_value.getVal()

#        iGen_value = modelVjet_test.createIntegral(RooArgSet(J_mass),RooFit.NormSet(RooArgSet(J_mass)), RooFit.Range("VRrange,SRrange"))
#        print "iGen_value:", iGen_value.getVal()

#        iGen_value = modelVjet_test.createIntegral(RooArgSet(J_mass), RooFit.Range("h_reasonable_range"))
#        print "iGen_value:", iGen_value.getVal()

#        iGen_value = modelVjet_test.createIntegral(RooArgSet(J_mass))
#        print "iGen_value:", iGen_value.getVal()


        # ------------------------------------------------------------------- 
        # repeat thousand times to see bias distribution

        h_Bias = TH1D("h_Bias","h_Bias",80,-1,1);

        Jmass_frame4 = J_mass.frame(RooFit.Title("test frame4"))

	times_max = 50000

        constVjet_test.setConstant(True)
        offsetVjet_test.setConstant(True)
        widthVjet_test.setConstant(True)

        constVjet_test3   = RooRealVar("constVjet_test3",   "slope of the exp", constVjet_value_fit_MC , -1.,   0.)
        offsetVjet_test3  = RooRealVar("offsetVjet_test3",  "offset of the erf", offsetVjet_value_fit_MC ,   -50., 200.)
        widthVjet_test3   = RooRealVar("widthVjet_test3",   "width of the erf", widthVjet_value_fit_MC ,     1., 200.)


	for times in range(0,times_max):  

    		# inside loop
#		print "times:", times
		if times % 10 == 0 :
			print "Processing times:", times+1 ,"of", times_max 


		# generate pseudo-data
#		n_1_prime = gRandom->Poisson(n_1);
		Entries_pseudo_data_fluc = gRandom.Poisson( Entries_pseudo_data )
#		print "Entries_pseudo_data:", Entries_pseudo_data
#		print "Entries_pseudo_data_fluc:", Entries_pseudo_data_fluc

#		pseudo_data2 = modelVjet_test.generate(RooArgSet(J_mass),Entries_pseudo_data )
                pseudo_data2 = modelVjet_test.generate(RooArgSet(J_mass),Entries_pseudo_data_fluc )

#        	pseudo_data2.plotOn(Jmass_frame3,RooFit.LineColor(4),RooFit.Range("h_reasonable_range"))

                # take out VR+SR		
	        pseudo_data_SB2 = RooDataSet("pseudo_data_SB2", "pseudo_data_SB2", RooArgSet(J_mass), RooFit.Import(pseudo_data2), RooFit.Cut("fatjet1_prunedMassCorr<65 || fatjet1_prunedMassCorr>135") )

		# use other PDF to fit
#		print "constVjet_value_fit_MC:",constVjet_value_fit_MC

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

#.........这里部分代码省略.........

    # fit to secondary bkg in MC (whole range)
    frTop = TopMass.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),
    )

    if VERBOSE:
        print "********** Fit result [JET MASS TOP] ***" + "*" * 40, "\n", frTop.Print(), "\n", "*" * 80

    # likelihoodScan(TopMass, setTop, [offsetTop, widthTop])

    # *******************************************************#
    #                                                       #
    #                 All bkg normalization                 #
    #                                                       #
    # *******************************************************#

    #    nVjet.setConstant(False)
    #    nVjet2.setConstant(False)
    #
    #    constVjet.setConstant(False)
    #    offsetVjet.setConstant(False)
    #    widthVjet.setConstant(False)
    #    a0Vjet.setConstant(False)
    #    a1Vjet.setConstant(False)
    #    a2Vjet.setConstant(False)

    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)

    nVV.setConstant(True)
    nTop.setConstant(True)
    nVjet.setConstant(False)
    nVjet2.setConstant(False)

    # Final background model by adding the main+secondary pdfs (using 'coef': ratio of the secondary/main, from MC)
    TopMass_ext = RooExtendPdf("TopMass_ext", "extended p.d.f", TopMass, nTop)
    VVMass_ext = RooExtendPdf("VVMass_ext", "extended p.d.f", VVMass, nVV)
    VjetMass_ext = RooExtendPdf("VjetMass_ext", "extended p.d.f", VjetMass, nVjet)
    VjetMass2_ext = RooExtendPdf("VjetMass_ext", "extended p.d.f", VjetMass, nVjet2)
    BkgMass = RooAddPdf(
        "BkgMass", "BkgMass", RooArgList(TopMass_ext, VVMass_ext, VjetMass_ext), RooArgList(nTop, nVV, nVjet)
    )

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
           rrv_alpha1_CB = sig_param;
          elif TString(sig_param.GetName()).Contains("rrv_alpha2_CB_BulkG_WW"):
           rrv_alpha2_CB = sig_param;
          elif TString(sig_param.GetName()).Contains("rrv_n1_CB_BulkG_WW"):
           rrv_n1_CB = sig_param;
          elif TString(sig_param.GetName()).Contains("rrv_n2_CB_BulkG_WW"):
           rrv_n2_CB = sig_param;
          sig_param = sig_par.Next();

        ### copy all the datasets
        getattr(new_workspace,"import")(old_workspace.data(datasetname+"_xww_"+options.channel+"_"+options.category));

        ## Breit Wigner Core  --> mean and the resolution fixed
        rrv_width_BW = RooRealVar("rrv_width_BW_BulkG_WW_"+options.channel+"_"+options.category,"rrv_width_BW_BulkG_WW_"+options.channel+"_"+options.category,mass[iMass]*gammaVal);

        rrv_mean_BW.setConstant(kTRUE);
        rrv_width_BW.setConstant(kTRUE);

        ## Breit Wigner Core  --> mean sys are to be applied on the mean of BW
        rrv_mean_scale_p1 = RooRealVar("CMS_sig_p1_jes","CMS_sig_p1_jes",0);
        rrv_mean_scale_p1.setConstant(kTRUE);
        rrv_mean_scale_p1.setError(1);
         
        rrv_mean_scale_p2 = RooRealVar("CMS_sig_p1_scale_em","CMS_sig_p1_scale_em",0);
        rrv_mean_scale_p2.setConstant(kTRUE);
        rrv_mean_scale_p2.setError(1);
                                    
        rrv_mean_scale_X1 = RooRealVar("rrv_mean_shift_scale_lep_BW_"+options.channel+"_"+options.category,"rrv_mean_shift_scale_lep_BW_"+options.channel+"_"+options.category,float(0.001));
        rrv_mean_scale_X1.setConstant(kTRUE);
        rrv_mean_scale_X2 = RooRealVar("rrv_mean_shift_scale_jes_BW_"+options.channel+"_"+options.category,"rrv_mean_shift_scale_jes_BW_"+options.channel+"_"+options.category,float(0.033));
        rrv_mean_scale_X2.setConstant(kTRUE);

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
        print ">> Creating datacard and workspace for %s resonance with m = %i GeV..."%(args.signal_model, int(mass))

        # get signal shape
        hSig = signal_shapes_file.Get( "h_" + args.signal_model + "_" + str(int(mass)) )
        # normalize signal shape to the expected event yield (works even if input shapes are not normalized to unity)
        hSig.Scale(signalCrossSection*lumi/hSig.Integral()) # divide by a number that provides roughly an r value of 1-10
        print "[debug] hSig.Integral() = " + str(hSig.Integral())
        print "[debug] hSig.GetMean() = " + str(hSig.GetMean())
        print "[debug] hSig.GetRMS() = " + str(hSig.GetRMS())
        rooSigHist = RooDataHist('rooSigHist','rooSigHist',RooArgList(mjj),hSig)
        rooSigHist.Print()
        signal = RooHistPdf('signal','signal',RooArgSet(mjj),rooSigHist)
        signal.Print()
        signal_norm = RooRealVar('signal_norm','signal_norm',0,-1e+05,1e+05)
        if args.fitBonly: signal_norm.setConstant()
        signal_norm.Print()

        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.fixP3: p3.setConstant()

        background = RooGenericPdf('background','(pow([email protected]/%.1f,@1)/pow(@0/%.1f,@[email protected]*log(@0/%.1f)))'%(sqrtS,sqrtS,sqrtS),RooArgList(mjj,p1,p2,p3))
        background.Print()
        dataInt = hData.Integral(hData.GetXaxis().FindBin(float(args.massMin)),hData.GetXaxis().FindBin(float(args.massMax)))
        background_norm = RooRealVar('background_norm','background_norm',dataInt,0.,1e+08)
        background_norm.Print()

        # S+B model
        model = RooAddPdf("model","s+b",RooArgList(background,signal),RooArgList(background_norm,signal_norm))

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setConstant [as 别名]
    fmip = RooRealVar('fmip', 'f_{mip}', 0.95, 0., 1.)
    if havePeds:
        if opts.sipm:
            lxgplus = RooAddPdf('lxgplus', 'lxgplus', lxg,
                                ws.pdf('pedPlusOne'), fmip)
        else:
            lxgplus = RooAddPdf('lxgplus', 'lxgplus', lxg, ws.pdf('ped'), fmip)
        ws.pdf('ped').plotOn(xf, RooFit.LineColor(kRed),
                             RooFit.LineStyle(kDashed),
                             RooFit.Normalization(sig_hist.GetEntries()/3,
                                                  RooAbsReal.Raw))
        fitter = lxgplus
    else:
        fitter = lxg
        fmip.setVal(1.0)
        fmip.setConstant(True)

    xf.Draw()
    # fitter.Print('v')
    # raise 'Stop here'
    fr = fitter.fitTo(sigDS, RooFit.Minos(False), RooFit.Save(True))
    fitter.plotOn(xf)
    if (fmip.getVal() < 0.9):
        lxgplus.plotOn(xf, RooFit.LineColor(kRed+2),
                       RooFit.LineStyle(kDashed),
                       RooFit.Components('ped*'))
        lxgplus.plotOn(xf, RooFit.LineColor(myBlue),
                       RooFit.LineStyle(kDashed),
                       RooFit.Components('lxg'))
        
    #lxgplus.paramOn(xf)