Python RooRealVar.getVal方法代码示例

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
def get_num_sig_bkg(hist_DataTemplate,
                    hist_SignalTemplate,
                    hist_BackgdTemplate,
                    fit_range_min,
                    fit_range_max):
    '''Given 3 input histograms (TH1F), and a fit range, this function finds
    the amount of signal and background that sum up to the data histogram.
    It does histogram fits.'''
    # Find range of data template
    data_min = hist_DataTemplate.GetXaxis().GetXmin()
    data_max = hist_DataTemplate.GetXaxis().GetXmax()
    
    # Create basic variables
    x = RooRealVar("x","x",data_min,data_max)
    x.setBins(hist_DataTemplate.GetXaxis().GetNbins())  # Binned x values
    nsig = RooRealVar("nsig","number of signal events"    , 0, hist_DataTemplate.Integral())
    nbkg = RooRealVar("nbkg","number of background events", 0, hist_DataTemplate.Integral())
    
    # Create RooDataHists from input TH1Fs
    dh = RooDataHist("dh","dh",RooArgList(x),hist_DataTemplate)
    ds = RooDataHist("ds","ds",RooArgList(x),hist_SignalTemplate)
    db = RooDataHist("db","db",RooArgList(x),hist_BackgdTemplate)
    
    # Create Probability Distribution Functions from Monte Carlo
    sigPDF = RooHistPdf("sigPDF", "sigPDF", RooArgSet(x), ds)
    bkgPDF = RooHistPdf("bkgPDF", "bkgPDF", RooArgSet(x), db)
    
    model = RooAddPdf("model","(g1+g2)+a",RooArgList(bkgPDF,sigPDF),RooArgList(nbkg,nsig))
    
    # Find the edges of the bins that contain the fit range min/max
    data_min = hist_DataTemplate.GetXaxis().GetBinLowEdge(hist_DataTemplate.GetXaxis().FindFixBin(fit_range_min))
    data_max = hist_DataTemplate.GetXaxis().GetBinUpEdge(hist_DataTemplate.GetXaxis().FindFixBin(fit_range_max))
    
    r = model.fitTo(dh,RooFit.Save(),RooFit.Minos(0),RooFit.PrintEvalErrors(0),
                    RooFit.Extended(),RooFit.Range(data_min,data_max))
    r.Print("v")

    #print nsig.getVal(), nsig.getError(), nbkg.getVal(), nbkg.getError()
    #  Create pull distribution
    #mcstudy = RooMCStudy(model, RooArgSet(x), RooFit.Binned(1), RooFit.Silence(),
    #                     RooFit.Extended(),
    #                     RooFit.FitOptions(RooFit.Save(1),
    #                                       RooFit.PrintEvalErrors(0),
    #                                       RooFit.Minos(0))
    #                    )
    #mcstudy.generateAndFit(500)                          # Generate and fit toy MC
    #pull_dist = mcstudy.plotPull(nsig, -3.0, 3.0, 30, 1)  # make pull distribution
    pull_dist = None
    return [nsig.getVal(), nsig.getError(), nbkg.getVal(), nbkg.getError(), pull_dist]

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
def roofit_poisson_unbinned(data):
    """returns lambda, error of lambda"""

    x = RooRealVar('x', 'x', 0, max(data) * 10)
    lam = RooRealVar('lambda', 'lambda', 0.1, 0.000001, max(data))

    model = RooPoisson('model', 'model', x, lam)

    dataset = RooDataSet('data', 'data', RooArgSet(x))
    for val in data:
        x.setVal(val)
        dataset.add(RooArgSet(x))

    model.fitTo(dataset, RooFit.Save(), RooFit.PrintLevel(-1))
    return lam.getVal(), lam.getError()

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
def dump_simple_simul():
    # idea: take two data sets, both low counts, but same lambda
    # and fit said lambda simultanously with both sets

    r_x = RooRealVar('x', 'x', 0, bins)

    r_lam = RooRealVar('lam', 'lam', 0.5, 0.0, 1.5)
    r_lam1 = RooRealVar('lam1', 'lam1', 0.5, 0.0, 1.5)
    r_lam2 = RooRealVar('lam2', 'lam2', 0.5, 0.0, 1.5)

    model1 = RooPoisson('pois1', 'pois1', r_x, r_lam1)
    model2 = RooPoisson('pois2', 'pois2', r_x, r_lam2)

    r_index = RooCategory('index', 'index')
    r_index.defineType('1')
    r_index.defineType('2')

    simul_model = RooSimultaneous('model', 'model', r_index)
    simul_model.addPdf(model1, '1')
    simul_model.addPdf(model2, '2')

    data1 = RooDataSet('', '', RooArgSet(r_x))
    for val in unbinned_from_binned(xdata, ydata):
        r_x.setVal(val)
        data1.add(RooArgSet(r_x))
    r_index.setLabel('1')
    data1.addColumn(r_index)

    data2 = RooDataSet('', '', RooArgSet(r_x))
    for val in unbinned_from_binned(xdata, ydata2):
        r_x.setVal(val)
        data2.add(RooArgSet(r_x))
    r_index.setLabel('2')
    data2.addColumn(r_index)

    data1.append(data2)

    _result = simul_model.fitTo(data1)
    print(r_lam.getVal(), '+-', r_lam.getError(), lam, np.abs(r_lam.getVal() - lam) / lam)
    print(r_lam1.getVal(), '+-', r_lam1.getError(), lam, np.abs(r_lam1.getVal() - lam) / lam)
    print(r_lam2.getVal(), '+-', r_lam2.getError(), lam, np.abs(r_lam2.getVal() - lam) / lam)

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

#.........这里部分代码省略.........
        #    alpha.setConstant(True)
        N = RooRealVar("N","N",2.,0.,100.)
        #    N.setConstant(True)
        fsig = RooRealVar("fsig","fsig",0.9,0.,1.0)
        
        ## PDFs
        relBW = RooGenericPdf("relBW","relBW","@0/(pow(@0*@[email protected]*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",RooArgList(x,mean,width))
        CB = RooCBShape("CB","CB",x,meanCB,sigma,alpha,N)
        expo = RooExponential("expo","expo",x,slope)
        relBWTimesCB = RooFFTConvPdf("relBWTimesCB","relBWTimesCB",x,relBW,CB)
        relBWTimesCBPlusExp = RooAddPdf("relBWTimesCBPlusExp","relBWTimesCBPlusExp",relBWTimesCB,expo,fsig)

        # Fit
        frame = x.frame()
        h = histos[i]
        # h.Rebin(10)
        h.Sumw2()
        nbin = h.GetNbinsX()
        dh = RooDataHist("dh","dh",RooArgList(x),h)
        dh.plotOn(frame)
        relBWTimesCBPlusExp.fitTo(dh)
        relBWTimesCBPlusExp.plotOn(frame)
        relBWTimesCBPlusExp.paramOn(frame)

        # Plot
        c = TCanvas("c_"+i,"c_"+i) 
        c.SetFillColor(0)
        c.cd()
        frame.Draw()
        c.SaveAs(options.outDir+"/DimuonWithFit_"+i+".png")
        
        # Extract the result of the fit
        expParams = []
        expCoef = slope.getValV()
        fSig    = fsig.getValV()
        binLow = h.GetXaxis().FindBin(minM_fit)
        binHigh = h.GetXaxis().FindBin(maxM_fit)
        nEntries = h.Integral(binLow,binHigh)
        expParams = [expCoef,fSig,nEntries,binLow,binHigh]
        expPars[i] = expParams

        signalParams = [mean.getVal(),width.getVal(),meanCB.getVal(),sigma.getVal(),alpha.getVal(),N.getVal()]
        signalPars[i] = signalParams

        # Subtract the bkg from the histograms
        h_woBkg = h.Clone()
        h_bkg = TH1F("h_bkg_"+i,"Histogram of bkg events",nbin,xlow,xup)
        
        h_bkg.Sumw2()
        expNorm = (math.fabs(expCoef)*(1-fSig)*nEntries)/(math.exp(expCoef*minM_fit)-math.exp(expCoef*maxM_fit))
        for ibin in range(binLow,binHigh):
            w = integrateExp(expNorm,expCoef,h_bkg.GetBinLowEdge(ibin),h_bkg.GetBinLowEdge(ibin+1))
            h_bkg.SetBinContent(ibin,w)

        h_woBkg.Add(h_bkg,-1)
        nEvts_woBkg = h_woBkg.Integral(binLow,binHigh)
        h_woBkg.Scale(1/nEvts_woBkg)
        histosSubtr[i] = h_woBkg


        del expParams, c
        del relBWTimesCBPlusExp, relBW, CB, relBWTimesCB, expo
        del x, mean, width, sigma, fsig, N, alpha, slope, meanCB
        del frame, dh, h, h_woBkg, h_bkg

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
 
gauss3_ext = RooExtendPdf("gauss3_ext",  "extended p.d.f", gauss3,  nGauss3)

#gauss3_ext.plotOn(xframe2, RooFit.Normalization(1.0, RooAbsReal.RelativeExpected), RooFit.LineColor(RooFit.kGreen))

print "for gauss3 before fit"
print "mean3: ", mean3.Print()," sigma3: ",sigma3.Print()
print "nGauss3: ", nGauss3.Print()


# fit
gauss3_ext.fitTo(data1)

# plot

norm3 = nGauss3.getVal()

#gauss3_ext.plotOn(xframe2, RooFit.LineColor(RooFit.kRed))
#gauss3_ext.plotOn(xframe2, RooFit.Normalization(norm3, RooAbsReal.NumEvent), RooFit.LineColor(RooFit.kRed))
gauss3_ext.plotOn(xframe2, RooFit.Normalization(1.0, RooAbsReal.RelativeExpected), RooFit.LineColor(RooFit.kRed))

#gauss3.plotOn(xframe2, RooFit.Normalization(norm3, RooAbsReal.NumEvent), RooFit.LineColor(RooFit.kRed))

data1.plotOn(xframe2)
gauss1.plotOn(xframe2,RooFit.LineColor(RooFit.kBlue))

print "for data1"
print "number of events: ", data1.sumEntries()

print "for gauss3 after fit"
print "mean3: ", mean3.Print()," sigma3: ",sigma3.Print()

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

#.........这里部分代码省略.........
	    if args.widefit:	
		mframew = mjj.frame()
    	        rooDataHist2.plotOn(mframew, ROOT.RooFit.Name("data"))
                res6 = model.fitTo(rooDataHist2, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            	model.plotOn(mframew, ROOT.RooFit.Name("model1"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kRed))
            	res7 = model2.fitTo(rooDataHist2, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            	model2.plotOn(mframew, ROOT.RooFit.Name("model2"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kOrange))
            	res8 = model3.fitTo(rooDataHist2, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            	model3.plotOn(mframew, ROOT.RooFit.Name("model3"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kGreen))
            	res9 = model4.fitTo(rooDataHist2, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            	model4.plotOn(mframew, ROOT.RooFit.Name("model4"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kBlue))
            	res10 = model5.fitTo(rooDataHist2, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
            	model5.plotOn(mframew, ROOT.RooFit.Name("model5"), ROOT.RooFit.LineStyle(1), ROOT.RooFit.LineWidth(1), ROOT.RooFit.LineColor(ROOT.EColor.kViolet))

                if args.pyes:
                    c = TCanvas("c","c",800,800)
                    mframew.SetAxisRange(300.,1300.)
                    c.SetLogy()
#                   mframew.SetMaximum(10)
#                   mframew.SetMinimum(1)
                    mframew.Draw()
                    fitname = args.pdir+'/5funcfittowide_m'+str(mass)+fstr+'.pdf'
		    c.SaveAs(fitname)

                    cpull = TCanvas("cpull","cpull",800,800)
                    pulls = mframew.pullHist("data","model1")
                    pulls.Draw("ABX")
                    pullname = args.pdir+'/pullwidefit_m'+str(mass)+fstr+'.pdf'
                    cpull.SaveAs(pullname)


	    if args.chi2:
		    fullInt = model.createIntegral(RooArgSet(mjj))
		    norm = dataInt/fullInt.getVal()
		    chi1 = 0.
		    fullInt2 = model2.createIntegral(RooArgSet(mjj))
        	    norm2 = dataInt2/fullInt2.getVal()
	      	    chi2 = 0.
		    fullInt3 = model3.createIntegral(RooArgSet(mjj))
       		    norm3 = dataInt3/fullInt3.getVal()
	            chi3 = 0.
		    fullInt4 = model4.createIntegral(RooArgSet(mjj))
       		    norm4 = dataInt4/fullInt4.getVal()
         	    chi4 = 0.
		    fullInt5 = model5.createIntegral(RooArgSet(mjj))
	            norm5 = dataInt5/fullInt5.getVal()
     	            chi5 = 0.
		    for i in range(args.massMin, args.massMax):
        	        new = 0
			new2 = 0
			new3 = 0
			new4 = 0
			new5 = 0
			height = hData.GetBinContent(i)
	        	xLow = hData.GetXaxis().GetBinLowEdge(i)
			xUp = hData.GetXaxis().GetBinLowEdge(i+1)
			obs = height*(xUp-xLow)
			mjj.setRange("intrange",xLow,xUp)
			integ = model.createIntegral(RooArgSet(mjj),ROOT.RooFit.NormSet(RooArgSet(mjj)),ROOT.RooFit.Range("intrange"))
			exp = integ.getVal()*norm
			new = pow(exp-obs,2)/exp
                	chi1 = chi1 + new
			integ2 = model2.createIntegral(RooArgSet(mjj),ROOT.RooFit.NormSet(RooArgSet(mjj)),ROOT.RooFit.Range("intrange"))
                	exp2 = integ2.getVal()*norm2
                	new2 = pow(exp2-obs,2)/exp2
                	chi2 = chi2 + new2

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

#.........这里部分代码省略.........
    setDataSR = RooDataSet("setDataSR", "setDataSR", variables, RooFit.Cut(SRcut), RooFit.WeightVar(weight), RooFit.Import(treeData))
    setDataVR = RooDataSet("setDataVR", "setDataVR", variables, RooFit.Cut(VRcut), RooFit.WeightVar(weight), RooFit.Import(treeData)) # Observed in the VV mass, just for plotting purposes
    
    # same for the bkg datasets from MC, where we just apply the base selections (not blind)
    setVjet = RooDataSet("setVjet", "setVjet", variables, RooFit.Cut(baseCut), RooFit.WeightVar(weight), RooFit.Import(treeVjet))
    setVjetSB = RooDataSet("setVjetSB", "setVjetSB", variables, RooFit.Import(setVjet), RooFit.Cut(SBcut), RooFit.WeightVar(weight))
    setVjetSR = RooDataSet("setVjetSR", "setVjetSR", variables, RooFit.Import(setVjet), RooFit.Cut(SRcut), RooFit.WeightVar(weight))
    setVV = RooDataSet("setVV", "setVV", variables, RooFit.Cut(baseCut), RooFit.WeightVar(weight), RooFit.Import(treeVV))
    setVVSB = RooDataSet("setVVSB", "setVVSB", variables, RooFit.Import(setVV), RooFit.Cut(SBcut), RooFit.WeightVar(weight))
    setVVSR = RooDataSet("setVVSR", "setVVSR", variables, RooFit.Import(setVV), RooFit.Cut(SRcut), RooFit.WeightVar(weight))
    setTop = RooDataSet("setTop", "setTop", variables, RooFit.Cut(baseCut), RooFit.WeightVar(weight), RooFit.Import(treeTop))
    setTopSB = RooDataSet("setTopSB", "setTopSB", variables, RooFit.Import(setTop), RooFit.Cut(SBcut), RooFit.WeightVar(weight))
    setTopSR = RooDataSet("setTopSR", "setTopSR", variables, RooFit.Import(setTop), RooFit.Cut(SRcut), RooFit.WeightVar(weight))
    
    print "  Data events SB: %.2f" % setDataSB.sumEntries()
    print "  V+jets entries: %.2f" % setVjet.sumEntries()
    print "  VV, VH entries: %.2f" % setVV.sumEntries()
    print "  Top,ST entries: %.2f" % setTop.sumEntries()
    
    
    # the relative normalization of the varius bkg is taken from MC by counting all the events in the full fatJetMass range
    #coef = RooRealVar("coef", "coef", setVV.sumEntries()/setVjet.sumEntries(),0.,1.)
    coef_VV_Vjet = RooRealVar("coef2_1", "coef2_1", setVV.sumEntries()/setVjet.sumEntries(), 0., 1.)
    coef_Top_VVVjet = RooRealVar("coef3_21", "coef3_21", setTop.sumEntries()/(setVjet.sumEntries()+setVV.sumEntries()),0.,1.);
    coef_VV_Vjet.setConstant(True)
    coef_Top_VVVjet.setConstant(True)
    
    # Define entries
    entryVjet = RooRealVar("entryVjets",  "V+jets normalization", setVjet.sumEntries(), 0., 1.e6)
    entryVV = RooRealVar("entryVV",  "VV normalization", setVV.sumEntries(), 0., 1.e6)
    entryTop = RooRealVar("entryTop",  "Top normalization", setTop.sumEntries(), 0., 1.e6)
    
    entrySB = RooRealVar("entrySB",  "Data SB normalization", setDataSB.sumEntries(SBcut), 0., 1.e6)
    entrySB.setError(math.sqrt(entrySB.getVal()))
    
    entryLSB = RooRealVar("entryLSB",  "Data LSB normalization", setDataSB.sumEntries(LSBcut), 0., 1.e6)
    entryLSB.setError(math.sqrt(entryLSB.getVal()))

    entryHSB = RooRealVar("entryHSB",  "Data HSB normalization", setDataSB.sumEntries(HSBcut), 0., 1.e6)
    entryHSB.setError(math.sqrt(entryHSB.getVal()))
    
    #*******************************************************#
    #                                                       #
    #                    NORMALIZATION                      #
    #                                                       #
    #*******************************************************#
    
    # set reasonable ranges for J_mass and X_mass
    # these are used in the fit in order to avoid ROOFIT to look in regions very far away from where we are fitting 
    J_mass.setRange("h_reasonable_range", LOWMIN, HIGMAX)
    X_mass.setRange("X_reasonable_range", XBINMIN, XBINMAX)
    
    # Set RooArgSets once for all, see https://root.cern.ch/phpBB3/viewtopic.php?t=11758
    jetMassArg = RooArgSet(J_mass)
    
    #*******************************************************#
    #                                                       #
    #                 V+jets normalization                  #
    #                                                       #
    #*******************************************************#
    
    # Variables for V+jets
    constVjet   = RooRealVar("constVjet",   "slope of the exp",      -0.020, -1.,   0.)
    offsetVjet  = RooRealVar("offsetVjet",  "offset of the erf",     30.,   -50., 200.)
    widthVjet   = RooRealVar("widthVjet",   "width of the erf",     100.,     1., 200.)
    offsetVjet.setConstant(True)

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

#.........这里部分代码省略.........
        N_signal_before_fit = N_ttbar_before_fit + N_SingleTop_before_fit
        
        N_ttbar_error_before_fit = sum(histograms['TTJet'].errors())
        N_SingleTop_error_before_fit = sum(histograms['SingleTop'].errors())
        N_vjets_error_before_fit = sum(histograms['V+Jets'].errors())
        N_QCD_error_before_fit = sum(histograms['QCD'].errors())
        
        if (N_SingleTop_before_fit != 0):
            TTJet_SingleTop_ratio = N_ttbar_before_fit / N_SingleTop_before_fit
        else:
            print 'Bin ', variable_bin, ': ttbar/singleTop ratio undefined for %s channel! Setting to 0.' % channel
            TTJet_SingleTop_ratio = 0
            
        
        leptonAbsEta = RooRealVar("leptonAbsEta", "leptonAbsEta", 0., 2.4)
        # this has to move to dps.utils.Fitting.py
        vars = RooArgList()
        vars.add(leptonAbsEta)
        vars_set = RooArgSet()
        vars_set.add(leptonAbsEta)
        n_event_obs = histograms['data'].Integral()
        
        lowerBound = 0. 
        upperBound = n_event_obs + 10 * sqrt(n_event_obs)
        n_init = n_event_obs / 2.
        
        data = RooDataHist("data", "dataset with leptonAbsEta", vars, histograms['data'])
        rh_vj = RooDataHist("rh_vj", "vj", vars, histograms['V+Jets'])
        rh_qcd = RooDataHist("rh_qcd", "qcd", vars, histograms['QCD'])
        rh_signal = RooDataHist("rh_signal", "signal", vars, h_eta_signal)
        
        pdf_vj = RooHistPdf ("pdf_vj", "V+Jets pdf", vars_set, rh_vj, 0)
        pdf_qcd = RooHistPdf("pdf_qcd", "QCD pdf ", vars_set, rh_qcd, 0)
        pdf_signal = RooHistPdf("pdf_signal", "single top pdf", vars_set, rh_signal, 0)
        
        # RooRealVar(const char *name, const char *title, Double_t value, Double_t minValue, Double_t maxValue, const char *unit) :
        nSignal = RooRealVar("nSignal", "number of single top + ttbar events", N_signal_before_fit, lowerBound, upperBound, "event")
        nvj = RooRealVar  ("nvj", "number of V+Jets bgnd events", N_vjets_before_fit, lowerBound, upperBound, "event")
        nqcd = RooRealVar("nqcd", "number of QCD bgnd events", N_QCD_error_before_fit, lowerBound, upperBound, "event")
        
        model = RooAddPdf("model", "sig+vj+qcd",
                          RooArgList(pdf_signal, pdf_vj, pdf_qcd),
                          RooArgList(nSignal, nvj, nqcd)
                          )
        vj_constraint = RooGaussian("nvj_constraint", "nvj_constraint", nvj, RooFit.RooConst(N_vjets_before_fit), RooFit.RooConst(0.5 * N_vjets_before_fit))
        qcd_constraint = RooGaussian("nqcd_constraint", "nqcd_constraint", nqcd, RooFit.RooConst(N_qcd_before_fit), RooFit.RooConst(2 * N_qcd_before_fit))  
        model_with_constraints = RooProdPdf("model_with_constraints", "model with gaussian constraints",
                                            RooArgSet(model, vj_constraint, qcd_constraint), RooLinkedList())
        model_with_constraints.fitTo(data, RooFit.Minimizer("Minuit2", "Migrad"))  #WARNING: number of cores changes the results!!!
#         nll = model.createNLL(data, RooFit.NumCPU(2))
#         RooMinuit(nll).migrad()
#         frame1 = nSignal.frame(RooFit.Bins(100), RooFit.Range(lowerBound, n_event_obs), RooFit.Title("LL and profileLL in nSignal")) 
#         nll.plotOn(frame1, RooFit.ShiftToZero())
#         frame2 = nvj.frame(RooFit.Bins(100), RooFit.Range(lowerBound, n_event_obs), RooFit.Title("LL and profileLL in nvj"))
#         nll.plotOn(frame2, RooFit.ShiftToZero())
#         frame3 = nqcd.frame(RooFit.Bins(100), RooFit.Range(lowerBound, n_event_obs), RooFit.Title("LL and profileLL in nqcd"))
#         nll.plotOn(frame3, RooFit.ShiftToZero())  
#         
#         pll_nSignal = nll.createProfile(nSignal)
#         pll_nSignal.plotOn(frame1, RooFit.LineColor(2)) 
#         frame1.SetMinimum(0)
#         frame1.SetMaximum(3)
#         
#         pll_nvj = nll.createProfile(nvj)
#         pll_nvj.plotOn(frame2, RooFit.LineColor(2)) 
#         frame2.SetMinimum(0)
#         frame2.SetMaximum(3)
#         
#         pll_nqcd = nll.createProfile(nqcd)
#         pll_nqcd.plotOn(frame3, RooFit.LineColor(2)) 
#         frame3.SetMinimum(0)
#         frame3.SetMaximum(3)
#         c = TCanvas("profilell","profilell",1200, 400)
#         c.Divide(3)
#         c.cd(1)
#         frame1.Draw()
#         c.cd(2)
#         frame2.Draw()
#         c.cd(3)
#         frame3.Draw()
#         c.SaveAs('profileLL.png')
#         model.fitTo(data, RooFit.Minimizer("Minuit2", "Migrad"), RooFit.NumCPU(1))#WARNING: number of cores changes the results!!!
        fit_results = {}
        fit_results['signal'] = (nSignal.getVal(), nSignal.getError())
        fit_results['QCD'] = ufloat(nqcd.getVal(), nqcd.getError())
        fit_results['V+Jets'] = ufloat(nvj.getVal(), nvj.getError())
        
        N_ttbar, N_SingleTop = decombine_result(fit_results['signal'], TTJet_SingleTop_ratio)
        fit_results['signal'] = ufloat(nSignal.getVal(), nSignal.getError())
        fit_results['TTJet'] = ufloat(N_ttbar)
        fit_results['SingleTop'] = ufloat(N_SingleTop)
        
        if results == {}:  # empty
            for sample in fit_results.keys():
                results[sample] = [fit_results[sample]]
        else:
            for sample in fit_results.keys():
                results[sample].append(fit_results[sample])
        
    return results, None, None

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
iLeft = gauss11_ext.createIntegral(xArg, RooFit.NormSet(xArg), RooFit.Range("left_gaussian_range"))
iTotal = gauss11_ext.createIntegral(xArg, RooFit.NormSet(xArg), RooFit.Range("whole_range"))

print "with NormSet"
print "iLeft: ", iLeft.getVal()
print "iTotal: ", iTotal.getVal()

#iLeft = gauss11_ext.createIntegral(xArg, RooFit.Range("left_gaussian_range"))
#iTotal = gauss11_ext.createIntegral(xArg, RooFit.Range("whole_range"))

#print "without NormSet"
#print "iLeft: ", iLeft.getVal()
#print "iTotal: ", iTotal.getVal()

n_total_evaluate = nGauss11.getVal() / iLeft.getVal()
print "n_total_evaluate: ", n_total_evaluate

#gauss11_ext.plotOn(xframe7,RooFit.Normalization( n_total_evaluate , RooAbsReal.NumEvent), RooFit.Range("whole_range"), RooFit.LineColor(RooFit.kRed))
#gauss11_ext.plotOn(xframe7,RooFit.Normalization( nGauss11.getVal() , RooAbsReal.NumEvent), RooFit.Range("whole_range"), RooFit.LineColor(RooFit.kRed))
gauss11_ext.plotOn(xframe7,RooFit.Normalization( nGauss11.getVal() , RooAbsReal.NumEvent), RooFit.Range("left_gaussian_range"), RooFit.LineColor(RooFit.kRed))

#gauss11_ext.plotOn(xframe7,RooFit.Normalization( nGauss11.getVal() , RooAbsReal.NumEvent), RooFit.Range("whole_range"), RooFit.LineColor(RooFit.kRed),RooFit.DrawOption("F"), RooFit.FillColor(RooFit.kRed), RooFit.FillStyle(1001) )


#data3.plotOn(xframe7)
data3_half.plotOn(xframe7)

#gauss11.plotOn(xframe7,RooFit.LineColor(RooFit.kRed))

# -------------------------------------------

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

#.........这里部分代码省略.........
	# -------------------------------------------------------------------    
    	# draw the setVjet
    	Jmass_frame = J_mass.frame(RooFit.Title("test frame"))
    	setVjet.plotOn(Jmass_frame)    

        # ------------------------------------------------------------------- 
    	# use a PDF to fit the dataset
    	print "fitFuncVjet is", fitFuncVjet 

        constVjet_value_initial = -0.020
	offsetVjet_value_initial = 30.
	widthVjet_value_initial = 100.

    	constVjet_test   = RooRealVar("constVjet_test",   "slope of the exp", constVjet_value_initial , -1.,   0.)
    	offsetVjet_test  = RooRealVar("offsetVjet_test",  "offset of the erf", offsetVjet_value_initial,   -50., 200.)
    	widthVjet_test   = RooRealVar("widthVjet_test",   "width of the erf",  widthVjet_value_initial,     1., 200.)

        modelVjet_test = RooErfExpPdf("modelVjet_test", "error function for V+jets mass", J_mass, constVjet_test, offsetVjet_test, widthVjet_test)

#	constVjet_test.Print()
#	offsetVjet_test.Print()
#	widthVjet_test.Print()

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

        frVjet_test = modelVjet_test.fitTo(setVjet, RooFit.SumW2Error(True), RooFit.Range("h_reasonable_range"), RooFit.Strategy(2), RooFit.Minimizer("Minuit2"), RooFit.Save(1), RooFit.PrintLevel(1 if VERBOSE else -1))

#        constVjet_test.Print()
#        offsetVjet_test.Print()
#        widthVjet_test.Print()

	constVjet_value_fit_MC = constVjet_test.getVal()
	offsetVjet_value_fit_MC = offsetVjet_test.getVal()
	widthVjet_value_fit_MC = widthVjet_test.getVal()

	print "constVjet_value_fit_MC:", constVjet_value_fit_MC, "offsetVjet_value_fit_MC:",offsetVjet_value_fit_MC,"widthVjet_value_fit_MC:",widthVjet_value_fit_MC


        modelVjet_test.plotOn(Jmass_frame,RooFit.LineColor(4))

        # ------------------------------------------------------------------- 
	# use the shape of fit to generate the psudo-data

	Entries_pseudo_data = setVjet.sumEntries()
#        Entries_pseudo_data = 502
#        Entries_pseudo_data = 5021


        pseudo_data = modelVjet_test.generate(RooArgSet(J_mass),Entries_pseudo_data )
#	pseudo_data = modelVjet_test.generate(RooArgSet(J_mass),setVjet.sumEntries()) 
#        pseudo_data = modelVjet_test.generate(RooArgSet(J_mass),502 )
#        pseudo_data = modelVjet_test.generate(RooArgSet(J_mass),5021 )
#	pseudo_data.Print("v")

	Jmass_frame2 = J_mass.frame(RooFit.Title("test frame2"))
	pseudo_data.plotOn(Jmass_frame2)
        modelVjet_test.plotOn(Jmass_frame2,RooFit.LineColor(4))

        # -------------------------------------------------------------------
	# make another dataset that remove the signal region 
        pseudo_data_SB = RooDataSet("pseudo_data_SB", "pseudo_data_SB", RooArgSet(J_mass), RooFit.Import(pseudo_data), RooFit.Cut("fatjet1_prunedMassCorr<65 || fatjet1_prunedMassCorr>135") )
	pseudo_data_SB.plotOn(Jmass_frame2,RooFit.LineColor(2))

        # -------------------------------------------------------------------

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
#rh_qcd.plotOn(leptonAbsEtaframe3, RooFit.MarkerSize(1), RooFit.DataError(etype));
#pdf_tt.plotOn(leptonAbsEtaframe4);
#pdf_wj.plotOn(leptonAbsEtaframe5);
#pdf_qcd.plotOn(leptonAbsEtaframe6);
#data.plotOn(leptonAbsEtaframe8, RooFit.MarkerSize(1)); #plot pseudo-data
#model.plotOn(leptonAbsEtaframe8); #plot composite pdf (s+b model)
#model.plotOn(leptonAbsEtaframe8, RooFit.Components(RooArgSet(pdf_tt)), RooFit.LineStyle(kDashed), RooFit.LineColor(kBlue + 1));
#model.plotOn(leptonAbsEtaframe8, RooFit.Components(RooArgSet(pdf_wj)), RooFit.LineStyle(kDashed), RooFit.LineColor(kRed + 1));
#model.plotOn(leptonAbsEtaframe8, RooFit.Components(RooArgSet(pdf_qcd)), RooFit.LineStyle(kDashed), RooFit.LineColor(kOrange - 6));
      
# single top plots
#rh_stop.plotOn(leptonAbsEtaframe9, RooFit.MarkerSize(1), RooFit.DataError(etype));
#pdf_stop.plotOn(leptonAbsEtaframe10);

#ntt_fit = ntt.getVal();
nSignal_fit = nSignal.getVal();
nwj_fit = nwj.getVal();
nzj_fit = nzj.getVal();
nqcd_fit = nqcd.getVal();
#nstop_fit = nstop.getVal();


nSignal_fiterr = nSignal.getError();
#ntt_fiterr = ntt.getError();
nwj_fiterr = nwj.getError();
nzj_fiterr = nzj.getError();
nqcd_fiterr = nqcd.getError();
#nstop_fiterr = nstop.getError();

print "Total events in signal region: ", n_event_obs;
print 'DATA:', n_event_obs

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
                rrv_number_TTbar.getVal()+
                rrv_number_STop.getVal());

        rrv_number_Total_background_MC.setError(TMath.Sqrt(
                rrv_number_WJets.getError()* rrv_number_WJets.getError()+
                rrv_number_VV.getError()* rrv_number_VV.getError()+
                rrv_number_TTbar.getError()* rrv_number_TTbar.getError()+
                rrv_number_STop.getError() *rrv_number_STop.getError()
                ));

        model_Total_background_MC = RooAddPdf("model_Total_background_MC_xww","model_Total_background_MC_xww",RooArgList(old_workspace.pdf("WJets_xww_%s_%s"%(options.channel,options.category)), old_workspace.pdf("VV_xww_%s_%s"%(options.channel,options.category)),old_workspace.pdf("TTbar_xww_%s_%s"%(options.channel,options.category)),old_workspace.pdf("STop_xww_%s_%s"%(options.channel,options.category))),RooArgList(rrv_number_WJets,rrv_number_VV,rrv_number_TTbar,rrv_number_STop));

        rrv_number_signal.setVal(rrv_number_signal.getVal()*6.25);

        #### scale factor in order to scale MC to data in the final plot -> in order to avoid the normalization to data which is done by default in rooFit
        scale_number_Total_background_MC = rrv_number_Total_background_MC.getVal()/old_workspace.data(datasetname+"_xww_"+options.channel+"_"+options.category).sumEntries();
        scale_number_signal = rrv_number_signal.getVal()/old_workspace.data(datasetname+"_xww_"+options.channel+"_"+options.category).sumEntries();

        model_Total_background_MC.plotOn(mplot,RooFit.Normalization(scale_number_Total_background_MC),RooFit.Name("total_MC"),RooFit.Components("WJets_xww_%s_%s,VV_xww_%s_%s,TTbar_xww_%s_%s,STop_xww_%s_%s"%(options.channel,options.category,options.channel,options.category,options.channel,options.category,options.channel,options.category)),RooFit.DrawOption("L"), RooFit.LineColor(kRed), RooFit.VLines(),RooFit.LineWidth(2));

            
        model_signal_background_MC = RooAddPdf("model_signal_background_MC_xww","model_signal_background_MC_xww",RooArgList(model_pdf,model_Total_background_MC),RooArgList(rrv_number_signal,rrv_number_Total_background_MC));

        model_signal_background_MC.plotOn(mplot,RooFit.Normalization(scale_number_Total_background_MC+scale_number_signal),RooFit.Name("total_SpB_MC"),RooFit.Components("BulkWW_xww_%s_%s,model_Total_background_MC_xww"%(options.channel,options.category)),RooFit.DrawOption("L"), RooFit.LineColor(kBlue), RooFit.VLines(),RooFit.LineWidth(2),RooFit.LineStyle(7));

        model_pdf.plotOn(mplot,RooFit.Name("total_S_MC"),RooFit.Normalization(scale_number_signal),RooFit.DrawOption("L"), RooFit.LineColor(kGreen+2), RooFit.VLines(),RooFit.LineWidth(2),RooFit.LineStyle(kDashed));

        os.system("mkdir -p ./plots_signal_width");
        name = TString("check_%.3f_%.3f"%(mass[iMass],gammaVal));
        name.ReplaceAll("0.","0_");
        mplot.GetYaxis().SetRangeUser(1e-3,mplot.GetMaximum()*1.2);

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

#.........这里部分代码省略.........
    setVjetSR = RooDataSet(
        "setVjetSR", "setVjetSR", variables, RooFit.Import(setVjet), RooFit.Cut(SRcut), RooFit.WeightVar(weight)
    )
    setVV = RooDataSet(
        "setVV", "setVV", variables, RooFit.Cut(baseCut), RooFit.WeightVar(weight), RooFit.Import(treeVV)
    )
    setVVSB = RooDataSet(
        "setVVSB", "setVVSB", variables, RooFit.Import(setVV), RooFit.Cut(SBcut), RooFit.WeightVar(weight)
    )
    setVVSR = RooDataSet(
        "setVVSR", "setVVSR", variables, RooFit.Import(setVV), RooFit.Cut(SRcut), RooFit.WeightVar(weight)
    )
    setTop = RooDataSet(
        "setTop", "setTop", variables, RooFit.Cut(baseCut), RooFit.WeightVar(weight), RooFit.Import(treeTop)
    )
    setTopSB = RooDataSet(
        "setTopSB", "setTopSB", variables, RooFit.Import(setTop), RooFit.Cut(SBcut), RooFit.WeightVar(weight)
    )
    setTopSR = RooDataSet(
        "setTopSR", "setTopSR", variables, RooFit.Import(setTop), RooFit.Cut(SRcut), RooFit.WeightVar(weight)
    )

    print "  Data events SB: %.2f" % setDataSB.sumEntries()
    print "  V+jets entries: %.2f" % setVjet.sumEntries()
    print "  VV, VH entries: %.2f" % setVV.sumEntries()
    print "  Top,ST entries: %.2f" % setTop.sumEntries()

    nVV = RooRealVar("nVV", "VV normalization", setVV.sumEntries(SBcut), 0.0, 2 * setVV.sumEntries(SBcut))
    nTop = RooRealVar("nTop", "Top normalization", setTop.sumEntries(SBcut), 0.0, 2 * setTop.sumEntries(SBcut))
    nVjet = RooRealVar("nVjet", "Vjet normalization", setDataSB.sumEntries(), 0.0, 2 * setDataSB.sumEntries(SBcut))
    nVjet2 = RooRealVar("nVjet2", "Vjet2 normalization", setDataSB.sumEntries(), 0.0, 2 * setDataSB.sumEntries(SBcut))

    # Apply Top SF
    nTop.setVal(nTop.getVal() * topSF[nLept][nBtag])
    nTop.setError(nTop.getVal() * topSFErr[nLept][nBtag])

    # Define entries
    entryVjet = RooRealVar("entryVjets", "V+jets normalization", setVjet.sumEntries(), 0.0, 1.0e6)
    entryVV = RooRealVar("entryVV", "VV normalization", setVV.sumEntries(), 0.0, 1.0e6)
    entryTop = RooRealVar("entryTop", "Top normalization", setTop.sumEntries(), 0.0, 1.0e6)

    entrySB = RooRealVar("entrySB", "Data SB normalization", setDataSB.sumEntries(SBcut), 0.0, 1.0e6)
    entrySB.setError(math.sqrt(entrySB.getVal()))

    entryLSB = RooRealVar("entryLSB", "Data LSB normalization", setDataSB.sumEntries(LSBcut), 0.0, 1.0e6)
    entryLSB.setError(math.sqrt(entryLSB.getVal()))

    entryHSB = RooRealVar("entryHSB", "Data HSB normalization", setDataSB.sumEntries(HSBcut), 0.0, 1.0e6)
    entryHSB.setError(math.sqrt(entryHSB.getVal()))

    ###################################################################################
    #        _   _                                                                    #
    #       | \ | |                          | (_)         | | (_)                    #
    #       |  \| | ___  _ __ _ __ ___   __ _| |_ ___  __ _| |_ _  ___  _ __          #
    #       | . ` |/ _ \| '__| '_ ` _ \ / _` | | / __|/ _` | __| |/ _ \| '_ \         #
    #       | |\  | (_) | |  | | | | | | (_| | | \__ \ (_| | |_| | (_) | | | |        #
    #       |_| \_|\___/|_|  |_| |_| |_|\__,_|_|_|___/\__,_|\__|_|\___/|_| |_|        #
    #                                                                                 #
    ###################################################################################
    # fancy ASCII art thanks to, I guess, Jose

    # start by creating the fit models to get the normalization:
    # * MAIN and SECONDARY bkg are taken from MC by fitting the whole J_mass range
    # * The two PDFs are added together using the relative normalizations of the two bkg from MC
    # * DATA is then fit in the sidebands only using the combined bkg PDF
    # * The results of the fit are then estrapolated in the SR and the integral is evaluated.

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

#.........这里部分代码省略.........
    # a1 = RooRealVar('#alpha1', '#alpha1', a1_val)
    # a2 = RooRealVar('#alpha2', '#alpha2', a2_val)

    # n1 = RooRealVar('n1', 'n1', 1.7,1.0,5.)
    # n2 = RooRealVar('n2', 'n2', 1.7,0.,4.)
    # cb_pdf = My_double_CB('chib', 'chib', x, mean, sigma, a1, n1, a2, n2)
    # result = cb_pdf.fitTo(data, RooFit.Save())

    
    
    if returnOnlyFitResult:
        return result
    
    # define frame
    frame = x.frame()
    frame.SetNameTitle("fit_resonance","Fit Resonanace")
    frame.GetXaxis().SetTitle(x_axis_label )
    frame.GetYaxis().SetTitle( "Events/5 MeV " )
    frame.GetXaxis().SetTitleSize(0.04)
    frame.GetYaxis().SetTitleSize(0.04)
    frame.GetXaxis().SetTitleOffset(1.1)
    frame.GetXaxis().SetLabelSize(0.04)
    frame.GetYaxis().SetLabelSize(0.04)
    frame.SetLineWidth(1)
    frame.SetTitle(plotTitle)
    
    # plot things on frame
    data.plotOn(frame, RooFit.MarkerSize(0.7))
    cb_pdf.plotOn(frame, RooFit.LineWidth(2))
    

    # chiSquare legend
    chi2 = frame.chiSquare()
    probChi2 = TMath.Prob(chi2*ndof, ndof)
    chi2 = round(chi2,2)
    probChi2 = round(probChi2,2)
    leg = TLegend(0.3,0,.10,.10)
    leg.SetBorderSize(0)
    leg.SetFillStyle(0)
    cb_pdf.paramOn(frame, RooFit.Layout(0.17,0.56,0.93))
    leg.AddEntry(0,'#chi^{2} ='+str(chi2),'')
    leg.AddEntry(0,'Prob #chi^{2} = '+str(probChi2),'')
    leg.SetTextSize(0.04)
    frame.addObject(leg)

    if legendOnPlot:
        legend = TLegend(.08,.5,.55,.7)
        #legend.SetTextSize(0.04)
        legend.SetFillStyle(0)
        legend.SetBorderSize(0)
        #legend.AddEntry(0,'CMS','')
        legend.AddEntry(0,str(cuts.upsilon_pt_lcut)+' GeV < p_{T}(#Upsilon) < '+str(cuts.upsilon_pt_hcut)+' GeV','')
        #legend.AddEntry(0,'p_{T}(#Upsilon)<'+str(cuts.upsilon_pt_hcut),'')
        frame.addObject(legend)
        
        title = TLegend(.8,.75,.9,.93)
        title.SetTextSize(0.15)
        title.SetFillStyle(0)
        title.SetBorderSize(0)
        title.AddEntry(0,plotTitle,'')
        frame.addObject(title)
    
     # Canvas
    c1=TCanvas(output_name+output_suffix,output_name+output_suffix)
    frame.Draw()

    if drawPulls:
        #c1=TCanvas(output_name+output_suffix,output_name+output_suffix,700, 625)
        hpull = frame.pullHist()
        framePulls = x.frame()
        framePulls.SetTitle(';;Pulls')
        framePulls.GetYaxis().SetLabelSize(0.18)
        framePulls.GetYaxis().SetTitle('Pulls')
        framePulls.GetYaxis().SetTitleSize(0.18)
        framePulls.GetYaxis().SetTitleOffset(0.15)
        framePulls.GetYaxis().SetNdivisions(005)
        framePulls.GetXaxis().SetLabelSize(0.16)
        framePulls.GetXaxis().SetTitle('')
        line0 = TLine(9.8, 0, 9.96, 0)
        line0.SetLineColor(ROOT.kBlue)
        line0.SetLineWidth(2)
        framePulls.addObject(line0)
        framePulls.addPlotable(hpull,"P") 
        framePulls.SetMaximum(5)
        framePulls.SetMinimum(-5)
        pad1 = TPad("pad1", "The pad 80% of the height",0.0,0.2,1.0,1.0)
        pad1.cd()
        frame.Draw()
        pad2 = TPad("pad2", "The pad 20% of the height",0.0,0.01,1.0,0.2)
        pad2.cd()
        framePulls.Draw()
        c1.cd()
        pad1.Draw()
        pad2.Draw()
    #c1.SaveAs(output_name+output_suffix+'.png')
    print 'Chi2 = '+str(frame.chiSquare())

    pdf_parameters = CB_parameters(mean=mean.getVal(), sigma=sigma.getVal(), n1=n1.getVal(), n2=n2.getVal(), a1=a1.getVal(), a2=a2.getVal(), s_mean=mean.getError(), s_sigma=sigma.getError(), s_n1=n1.getError(), s_n2=n2.getError(), s_a1=a1.getError(), s_a2=a2.getError(), chiSquare = frame.chiSquare())
    #pdf_parameters.saveToFile("CB_"+output_name+output_suffix+".txt")
    return pdf_parameters, c1

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getVal [as 别名]
    for varList in condObs.itervalues() :
        valList = [ getObsVal( condObsVars[var], evObs ) for var in varList ]
        for comp in comps :
            pdfVals[comp] *= condObsDist( varList, valList, comp )

    # decide if event is signal or background
    totVal = sum( pdfVals.itervalues() )
    addVal = 0.
    genVal = randGen.Rndm()
    for comp in comps :
        addVal += pdfVals[comp]
        if addVal < genVal * totVal : continue
        if genMass :
            splitData[comp][stateName].add(dataObs)
        elif weightData :
            splitData[comp]['total'].add( dataObs, sigWeight.getVal() )
        else :
            splitData[comp]['total'].add(dataObs)
        break

if genMass :
    # generate J/psiKK mass
    genObs = RooArgSet( massPdfObs.find('mass') )
    dataObs.add(genObs)
    genData = RooDataSet( dataSet.GetName(), dataSet.GetTitle(), dataObs, WeightVar = 'wMC' )
    for comp in comps :
        splitData[comp]['total'] = RooDataSet( '%s_%s' % ( dataSet.GetName(), comp ), dataSet.GetTitle(), dataObs )
        for state in splitCat :
            stateName = state.GetName()
            splitData[comp][stateName] = pdfs[comp][stateName].generate( genObs, ProtoData = ( splitData[comp][stateName], False, False ) )
            splitData[comp]['total'].append( splitData[comp][stateName] )