本文整理汇总了Python中ROOT.RooAddPdf.fitTo方法的典型用法代码示例。如果您正苦于以下问题:Python RooAddPdf.fitTo方法的具体用法?Python RooAddPdf.fitTo怎么用?Python RooAddPdf.fitTo使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类ROOT.RooAddPdf
的用法示例。
在下文中一共展示了RooAddPdf.fitTo方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: test_correlated_values
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
def test_correlated_values():
try:
import uncertainties
except ImportError:
raise SkipTest("uncertainties package is not installed")
from rootpy.stats.correlated_values import correlated_values
# construct pdf and toy data following example at
# http://root.cern.ch/drupal/content/roofit
# --- Observable ---
mes = RooRealVar("mes", "m_{ES} (GeV)", 5.20, 5.30)
# --- Parameters ---
sigmean = RooRealVar("sigmean", "B^{#pm} mass", 5.28, 5.20, 5.30)
sigwidth = RooRealVar("sigwidth", "B^{#pm} width", 0.0027, 0.001, 1.)
# --- Build Gaussian PDF ---
signal = RooGaussian("signal", "signal PDF", mes, sigmean, sigwidth)
# --- Build Argus background PDF ---
argpar = RooRealVar("argpar", "argus shape parameter", -20.0, -100., -1.)
background = RooArgusBG("background", "Argus PDF",
mes, RooFit.RooConst(5.291), argpar)
# --- Construct signal+background PDF ---
nsig = RooRealVar("nsig", "#signal events", 200, 0., 10000)
nbkg = RooRealVar("nbkg", "#background events", 800, 0., 10000)
model = RooAddPdf("model", "g+a",
RooArgList(signal,background),
RooArgList(nsig,nbkg))
# --- Generate a toyMC sample from composite PDF ---
data = model.generate(RooArgSet(mes), 2000)
# --- Perform extended ML fit of composite PDF to toy data ---
fitresult = model.fitTo(data, RooFit.Save(), RooFit.PrintLevel(-1))
nsig, nbkg = correlated_values(["nsig", "nbkg"], fitresult)
# Arbitrary math expression according to what the `uncertainties`
# package supports, automatically computes correct error propagation
sum_value = nsig + nbkg
value, error = sum_value.nominal_value, sum_value.std_dev
workspace = Workspace(name='workspace')
# import the data
assert_false(workspace(data))
with TemporaryFile():
workspace.Write()
示例2: get_num_sig_bkg
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [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]
示例3: fit
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [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
示例4: alpha
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
#.........这里部分代码省略.........
# 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)
a0Vjet = RooRealVar("a0Vjet", "width of the erf", -0.1, -5, 0)
a1Vjet = RooRealVar("a1Vjet", "width of the erf", 0.6, 0, 5)
a2Vjet = RooRealVar("a2Vjet", "width of the erf", -0.1, -1, 1)
# Define V+jets model
if fitFuncVjet == "ERFEXP": modelVjet = RooErfExpPdf("modelVjet", "error function for V+jets mass", J_mass, constVjet, offsetVjet, widthVjet)
elif fitFuncVjet == "EXP": modelVjet = RooExponential("modelVjet", "exp for V+jets mass", J_mass, constVjet)
elif fitFuncVjet == "POL": modelVjet = RooChebychev("modelVjet", "polynomial for V+jets mass", J_mass, RooArgList(a0Vjet, a1Vjet, a2Vjet))
elif fitFuncVjet == "POW": modelVjet = RooGenericPdf("modelVjet", "powerlaw for X mass", "@0^@1", RooArgList(J_mass, a0Vjet))
else:
print " ERROR! Pdf", fitFuncVjet, "is not implemented for Vjets"
exit()
# fit to main bkg in MC (whole range)
frVjet = modelVjet.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))
# integrals and number of events
iSBVjet = modelVjet.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange,HSBrange"))
iLSBVjet = modelVjet.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("LSBrange"))
iHSBVjet = modelVjet.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("HSBrange"))
iSRVjet = modelVjet.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("SRrange"))
iVRVjet = modelVjet.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg), RooFit.Range("VRrange"))
# Do not remove the following lines, integrals are computed here
iALVjet = modelVjet.createIntegral(jetMassArg, RooFit.NormSet(jetMassArg))
nSBVjet = iSBVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(SBcut)
nLSBVjet = iLSBVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(LSBcut)
nHSBVjet = iHSBVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(HSBcut)
nSRVjet = iSRVjet.getVal()/iALVjet.getVal()*setVjet.sumEntries(SRcut)
drawPlot("JetMass_Vjet", channel, J_mass, modelVjet, setVjet, binsJmass, frVjet)
if VERBOSE: print "********** Fit result [JET MASS Vjets] *"+"*"*40, "\n", frVjet.Print(), "\n", "*"*80
#*******************************************************#
# #
# VV, VH normalization #
# #
#*******************************************************#
# Variables for VV
# Error function and exponential to model the bulk
constVV = RooRealVar("constVV", "slope of the exp", -0.030, -0.1, 0.)
offsetVV = RooRealVar("offsetVV", "offset of the erf", 90., 1., 300.)
widthVV = RooRealVar("widthVV", "width of the erf", 50., 1., 100.)
erfrVV = RooErfExpPdf("baseVV", "error function for VV jet mass", J_mass, constVV, offsetVV, widthVV)
expoVV = RooExponential("baseVV", "error function for VV jet mass", J_mass, constVV)
# gaussian for the V mass peak
示例5: RooGenericPdf
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
bkg = RooGenericPdf('bkg','1/(exp(pow(@0/@1,@2))+1)',RooArgList(x,x0,p))
fsig= RooRealVar('fsig','fsig',0.5,0.,1.)
signal = RooAddPdf('signal','signal',sig,bkg,fsig)
# -----------------------------------------
# fit signal
canSname = 'can_Mjj'+str(mass)
canS = TCanvas(canSname,canSname,900,600)
gPad.SetLogy()
roohistSig = RooDataHist('roohist','roohist',RooArgList(x),hSig)
roohistSig.Print()
res_sig = signal.fitTo(roohistSig, RooFit.Save(ROOT.kTRUE))
res_sig.Print()
frame = x.frame()
roohistSig.plotOn(frame,RooFit.Binning(166))
signal.plotOn(frame)
signal.plotOn(frame,RooFit.Components('bkg'),RooFit.LineColor(ROOT.kRed),RooFit.LineWidth(2),RooFit.LineStyle(ROOT.kDashed))
#frame.GetXaxis().SetRangeUser(1118,6099)
frame.GetXaxis().SetRangeUser(1500,6000)
frame.GetXaxis().SetTitle('m_{jj} (GeV)')
frame.Draw()
parsSig = signal.getParameters(roohistSig)
parsSig.setAttribAll('Constant', True)
if histpdfSig:
示例6: range
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
g.GetYaxis().SetRangeUser(-0.12, 0)
from ROOT import kBlue
l.SetLineColor(kBlue)
l.Draw('same')
year_label.Draw()
plot_name = 'mean_res_st_ttrue_' + args[0][pos : -5] + '.pdf'
canvas.Print(os.path.join(plot_dir, plot_name), EmbedFonts = True)
else:
from ROOT import kGreen, kDashed
from P2VV.Utilities.Plotting import plot
from ROOT import TCanvas
for i in range(3):
result = model.fitTo(sdata, SumW2Error = False, **fitOpts)
if result.status() == 0 and abs(result.minNll()) < 5e5:
break
canvas = TCanvas("canvas", "canvas", 600, 400)
plot(canvas, t_diff_st, pdf = model, data = sdata, logy = True,
frameOpts = dict(Range = (-20, 20)),
yTitle = 'Candidates / (0.5)', dataOpts = dict(Binning = 80),
xTitle = '(t - t_{true}) / #sigma_{t}',
yScale = (1, 400000),
pdfOpts = dict(ProjWData = (RooArgSet(st), sdata, True)),
components = {'gexps' : dict(LineColor = kGreen, LineStyle = kDashed)})
year_label.Draw()
plot_name = 'tdiff_sigmat_' + args[0][pos : -5] + '.pdf'
canvas.Print(os.path.join(plot_dir, plot_name), EmbedFonts = True)
示例7: studyVqqResolution
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
#.........这里部分代码省略.........
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)
labels[ilab].SetFillStyle(0)
labels[ilab].SetTextFont(42)
labels[ilab].SetTextAlign(12)
示例8: RooGenericPdf
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
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))
rooDataHist = RooDataHist('rooDatahist','rooDathist',RooArgList(mjj),hData)
rooDataHist.Print()
if args.runFit:
res = model.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
if not args.decoBkg: res.Print()
# decorrelated background parameters for Bayesian limits
if args.decoBkg:
signal_norm.setConstant()
res = model.fitTo(rooDataHist, RooFit.Save(kTRUE), RooFit.Strategy(args.fitStrategy))
res.Print()
## temp workspace for the PDF diagonalizer
w_tmp = RooWorkspace("w_tmp")
deco = PdfDiagonalizer("deco",w_tmp,res)
# here diagonalizing only the shape parameters since the overall normalization is already decorrelated
background_deco = deco.diagonalize(background)
print "##################### workspace for decorrelation"
w_tmp.Print("v")
print "##################### original parameters"
示例9: rf501_simultaneouspdf
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
def rf501_simultaneouspdf():
signal_1, bkg_1, signal_2, bkg_2 = get_templates()
# C r e a t e m o d e l f o r p h y s i c s s a m p l e
# -------------------------------------------------------------
# Create observables
x = RooRealVar( "x", "x", 0, 200 )
x.setBins(n_bins)
nsig = RooRealVar( "nsig", "#signal events", N_signal_obs, 0., 2*N_data )
nbkg = RooRealVar( "nbkg", "#background events", N_bkg1_obs, 0., 2*N_data )
# Construct signal pdf
# mean = RooRealVar( "mean", "mean", mu4, 40, 200 )
# sigma = RooRealVar( "sigma", "sigma", sigma4, 0.1, 20 )
# gx = RooGaussian( "gx", "gx", x, mean, sigma )
roofit_signal_1 = RooDataHist( 'signal_1', 'signal_1', RooArgList(x), signal_1 )
signal_1_pdf = RooHistPdf ( 'signal_1_pdf' , 'signal_1_pdf', RooArgSet(x), roofit_signal_1)
# Construct background pdf
# mean_bkg = RooRealVar( "mean_bkg", "mean_bkg", mu3, 40, 200 )
# sigma_bkg = RooRealVar( "sigma_bkg", "sigma_bkg", sigma3, 0.1, 20 )
# px = RooGaussian( "px", "px", x, mean_bkg, sigma_bkg )
roofit_bkg_1 = RooDataHist( 'bkg_1', 'bkg_1', RooArgList(x), bkg_1 )
bkg_1_pdf = RooHistPdf ( 'bkg_1_pdf' , 'bkg_1_pdf', RooArgSet(x), roofit_bkg_1)
# Construct composite pdf
model = RooAddPdf( "model", "model", RooArgList( signal_1_pdf, bkg_1_pdf ), RooArgList( nsig, nbkg ) )
# C r e a t e m o d e l f o r c o n t r o l s a m p l e
# --------------------------------------------------------------
# Construct signal pdf.
# NOTE that sigma is shared with the signal sample model
y = RooRealVar( "y", "y", 0, 200 )
y.setBins(n_bins)
mean_ctl = RooRealVar( "mean_ctl", "mean_ctl", mu2, 0, 200 )
sigma_ctl = RooRealVar( "sigma", "sigma", sigma2, 0.1, 10 )
gx_ctl = RooGaussian( "gx_ctl", "gx_ctl", y, mean_ctl, sigma_ctl )
# Construct the background pdf
mean_bkg_ctl = RooRealVar( "mean_bkg_ctl", "mean_bkg_ctl", mu1, 0, 200 )
sigma_bkg_ctl = RooRealVar( "sigma_bkg_ctl", "sigma_bkg_ctl", sigma1, 0.1, 20 )
px_ctl = RooGaussian( "px_ctl", "px_ctl", y, mean_bkg_ctl, sigma_bkg_ctl )
# Construct the composite model
# f_ctl = RooRealVar( "f_ctl", "f_ctl", 0.5, 0., 20. )
model_ctl = RooAddPdf( "model_ctl", "model_ctl", RooArgList( gx_ctl, px_ctl ),
RooArgList( nsig, nbkg ) )
# G e t e v e n t s f o r b o t h s a m p l e s
# ---------------------------------------------------------------
real_data, real_data_ctl = get_data()
real_data_hist = RooDataHist( 'real_data_hist',
'real_data_hist',
RooArgList( x ),
real_data )
real_data_ctl_hist = RooDataHist( 'real_data_ctl_hist',
'real_data_ctl_hist',
RooArgList( y ),
real_data_ctl )
input_hists = MapStrRootPtr()
input_hists.insert( StrHist( "physics", real_data ) )
input_hists.insert( StrHist( "control", real_data_ctl ) )
# C r e a t e i n d e x c a t e g o r y a n d j o i n s a m p l e s
# ---------------------------------------------------------------------------
# Define category to distinguish physics and control samples events
sample = RooCategory( "sample", "sample" )
sample.defineType( "physics" )
sample.defineType( "control" )
# Construct combined dataset in (x,sample)
combData = RooDataHist( "combData", "combined data", RooArgList( x), sample ,
input_hists )
# C o n s t r u c t a s i m u l t a n e o u s p d f i n ( x , s a m p l e )
# -----------------------------------------------------------------------------------
# Construct a simultaneous pdf using category sample as index
simPdf = RooSimultaneous( "simPdf", "simultaneous pdf", sample )
# Associate model with the physics state and model_ctl with the control state
simPdf.addPdf( model, "physics" )
simPdf.addPdf( model_ctl, "control" )
#60093.048127 173.205689173 44.7112503776
# P e r f o r m a s i m u l t a n e o u s f i t
# ---------------------------------------------------
model.fitTo( real_data_hist,
RooFit.Minimizer( "Minuit2", "Migrad" ),
RooFit.NumCPU( 1 ),
# RooFit.Extended(),
# RooFit.Save(),
)
#.........这里部分代码省略.........
示例10: RooGaussian
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
# --- Build Gaussian PDF ---
signal = RooGaussian("signal", "signal PDF", mes, sigmean, sigwidth)
argpar = RooRealVar("argpar", "argus shape parameter", -20.0, -100., -1.)
background = RooArgusBG("background", "Argus PDF", mes, RooFit.RooConst(5.291), argpar)
# --- Construct signal+background PDF ---
nsig = RooRealVar("nsig", "#signal events", 200, 0., 10000)
nbkg = RooRealVar("nbkg", "#background events", 800, 0., 10000)
model = RooAddPdf("model", "g+a", RooArgList(signal, background), RooArgList(nsig, nbkg))
# --- Generate a toyMC sample from composite PDF ---
data = model.generate(RooArgSet(mes), 2000)
# --- Perform extended ML fit of composite PDF to toy data ---
model.fitTo(data)
# --- Plot toy data and composite PDF overlaid ---
mesframe = mes.frame()
data.plotOn(mesframe)
model.plotOn(mesframe)
model.plotOn(mesframe, RooFit.Components('background'), RooFit.LineStyle(kDashed))
mesframe.Draw()
print 'nsig:',nsig.getValV(), '+-', nsig.getError()
print 'nbkg:', nbkg.getValV(), '+-', nbkg.getError()
print 'mes:', mes.getValV(), '+-', mes.getError()
print 'mean:', sigmean.getValV(), '+-', sigmean.getError()
print 'width:', sigwidth.getValV(), '+-', sigwidth.getError()
print 'argpar:', argpar.getValV(), '+-', argpar.getError()
示例11: alpha
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
#.........这里部分代码省略.........
offsetVjet = RooRealVar("offsetVjet", "offset of the erf", 120.0, 80.0, 155.0)
if channel == "XWhenbb" or channel == "XZhmmb":
offsetVjet = RooRealVar("offsetVjet", "offset of the erf", 67.0, 50.0, 100.0)
if channel == "XWhmnb":
offsetVjet = RooRealVar("offsetVjet", "offset of the erf", 30.0, -50.0, 600.0)
if channel == "XZheeb":
offsetVjet.setMin(-400)
offsetVjet.setVal(0.0)
offsetVjet.setMax(1000)
widthVjet.setVal(1.0)
# Define V+jets model
if fitFuncVjet == "ERFEXP":
VjetMass = RooErfExpPdf("VjetMass", fitFuncVjet, J_mass, constVjet, offsetVjet, widthVjet)
elif fitFuncVjet == "EXP":
VjetMass = RooExponential("VjetMass", fitFuncVjet, J_mass, constVjet)
elif fitFuncVjet == "GAUS":
VjetMass = RooGaussian("VjetMass", fitFuncVjet, J_mass, offsetVjet, widthVjet)
elif fitFuncVjet == "POL":
VjetMass = RooChebychev("VjetMass", fitFuncVjet, J_mass, RooArgList(a0Vjet, a1Vjet, a2Vjet))
elif fitFuncVjet == "POW":
VjetMass = RooGenericPdf("VjetMass", fitFuncVjet, "@0^@1", RooArgList(J_mass, a0Vjet))
else:
print " ERROR! Pdf", fitFuncVjet, "is not implemented for Vjets"
exit()
if fitAltFuncVjet == "POL":
VjetMass2 = RooChebychev("VjetMass2", "polynomial for V+jets mass", J_mass, RooArgList(a0Vjet, a1Vjet, a2Vjet))
else:
print " ERROR! Pdf", fitAltFuncVjet, "is not implemented for Vjets"
exit()
# fit to main bkg in MC (whole range)
frVjet = VjetMass.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),
)
frVjet2 = VjetMass2.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),
)
if VERBOSE:
print "********** Fit result [JET MASS Vjets] *" + "*" * 40, "\n", frVjet.Print(), "\n", "*" * 80
# likelihoodScan(VjetMass, setVjet, [constVjet, offsetVjet, widthVjet])
# *******************************************************#
# #
# VV, VH normalization #
# #
# *******************************************************#
# Variables for VV
# Error function and exponential to model the bulk
constVV = RooRealVar("constVV", "slope of the exp", -0.030, -0.1, 0.0)
示例12: RooRealVar
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
p = RooRealVar('p','p',1,0,5)
x0 = RooRealVar('x0','x0',1000,100,5000)
bkg = RooGenericPdf('bkg','1/(exp(pow(@0/@1,@2))+1)',RooArgList(x,x0,p))
fsig= RooRealVar('fsig','fsig',0.5,0.,1.)
model = RooAddPdf('model','model',sig,bkg,fsig)
can = TCanvas('can_Mjj'+str(mass),'can_Mjj'+str(mass),900,600)
h.Draw()
gPad.SetLogy()
roohist = RooDataHist('roohist','roohist',RooArgList(x),h)
model.fitTo(roohist)
frame = x.frame()
roohist.plotOn(frame)
model.plotOn(frame)
model.plotOn(frame,RooFit.Components('bkg'),RooFit.LineColor(ROOT.kRed),RooFit.LineWidth(2),RooFit.LineStyle(ROOT.kDashed))
frame.Draw('same')
w = RooWorkspace('w','workspace')
getattr(w,'import')(model)
getattr(w,'import')(roohist)
w.Print()
w.writeToFile('RS'+str(mass)+'_workspace.root')
#----- keep the GUI alive ------------
if __name__ == '__main__':
rep = ''
示例13: fitChicSpectrum
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
def fitChicSpectrum(dataset,binname):
""" Fit chic spectrum"""
x = RooRealVar('Qvalue','Q',9.7,10.1)
x.setBins(80)
mean_1 = RooRealVar("mean_1","mean ChiB1",9.892,9,10,"GeV")
sigma_1 = RooRealVar("sigma_1","sigma ChiB1",0.0058,'GeV')
a1_1 = RooRealVar('#alpha1_1', '#alpha1_1', 0.748)
n1_1 = RooRealVar('n1_1', 'n1_1',2.8 )
a2_1 = RooRealVar('#alpha2_1', '#alpha2_1',1.739)
n2_1 = RooRealVar('n2_1', 'n2_1', 3.0)
deltam = RooRealVar('deltam','deltam',0.01943)
mean_2 = RooFormulaVar("mean_2","@[email protected]", RooArgList(mean_1,deltam))
sigma_2 = RooRealVar("sigma_2","sigma ChiB2",0.0059,'GeV')
a1_2 = RooRealVar('#alpha1_2', '#alpha1_2', 0.738)
n1_2 = RooRealVar('n1_2', 'n1_2', 2.8)
a2_2 = RooRealVar('#alpha2_2', '#alpha2_2', 1.699)
n2_2 = RooRealVar('n2_2', 'n2_2', 3.0)
parameters=RooArgSet()
parameters.add(RooArgSet(sigma_1, sigma_2))
parameters = RooArgSet(a1_1, a2_1, n1_1, n2_1)
parameters.add(RooArgSet( a1_2, a2_2, n1_2, n2_2))
chib1_pdf = My_double_CB('chib1', 'chib1', x, mean_1, sigma_1, a1_1, n1_1, a2_1, n2_1)
chib2_pdf = My_double_CB('chib2', 'chib2', x, mean_2, sigma_2, a1_2, n1_2, a2_2, n2_2)
#background
q01S_Start = 9.5
alpha = RooRealVar("#alpha","#alpha",1.5,-1,3.5)#0.2 anziche' 1
beta = RooRealVar("#beta","#beta",-2.5,-7.,0.)
q0 = RooRealVar("q0","q0",q01S_Start)#,9.5,9.7)
delta = RooFormulaVar("delta","TMath::Abs(@[email protected])",RooArgList(x,q0))
b1 = RooFormulaVar("b1","@0*(@[email protected])",RooArgList(beta,x,q0))
signum1 = RooFormulaVar( "signum1","( TMath::Sign( -1.,@[email protected] )+1 )/2.", RooArgList(x,q0) )
background = RooGenericPdf("background","Background", "signum1*pow(delta,#alpha)*exp(b1)", RooArgList(signum1,delta,alpha,b1) )
parameters.add(RooArgSet(alpha, beta, q0))
#together
chibs = RooArgList(chib1_pdf,chib2_pdf,background)
n_chib = RooRealVar("n_chib","n_chib",2075, 0, 100000)
ratio_21 = RooRealVar("ratio_21","ratio_21",0.5,0,1)
n_chib1 = RooFormulaVar("n_chib1","@0/([email protected])",RooArgList(n_chib, ratio_21))
n_chib2 = RooFormulaVar("n_chib2","@0/(1+1/@1)",RooArgList(n_chib, ratio_21))
n_background = RooRealVar('n_background','n_background',4550, 0, 50000)
ratio_list = RooArgList(n_chib1, n_chib2, n_background)
modelPdf = RooAddPdf('ModelPdf', 'ModelPdf', chibs, ratio_list)
frame = x.frame(RooFit.Title('m'))
range = x.setRange('range',9.7,10.1)
result = modelPdf.fitTo(dataset,RooFit.Save(),RooFit.Range('range'))
dataset.plotOn(frame,RooFit.MarkerSize(0.7))
modelPdf.plotOn(frame, RooFit.LineWidth(2) )
#plotting
canvas = TCanvas('fit', "", 1400, 700 )
canvas.Divide(1)
canvas.cd(1)
gPad.SetRightMargin(0.3)
gPad.SetFillColor(10)
modelPdf.paramOn(frame, RooFit.Layout(0.725,0.9875,0.9))
frame.Draw()
canvas.SaveAs( 'out-'+binname + '.png' )
示例14: rf501_simultaneouspdf
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
def rf501_simultaneouspdf():
# C r e a t e m o d e l f o r p h y s i c s s a m p l e
# -------------------------------------------------------------
# Create observables
x = RooRealVar("x", "x", 40, 200)
nsig = RooRealVar("nsig", "#signal events", 200, 0.0, 10000)
nbkg = RooRealVar("nbkg", "#background events", 800, 0.0, 200000)
# Construct signal pdf
mean = RooRealVar("mean", "mean", mu4, 40, 200)
sigma = RooRealVar("sigma", "sigma", sigma4, 0.1, 20)
gx = RooGaussian("gx", "gx", x, mean, sigma)
# Construct background pdf
mean_bkg = RooRealVar("mean_bkg", "mean_bkg", mu3, 40, 200)
sigma_bkg = RooRealVar("sigma_bkg", "sigma_bkg", sigma3, 0.1, 20)
px = RooGaussian("px", "px", x, mean_bkg, sigma_bkg)
# Construct composite pdf
model = RooAddPdf("model", "model", RooArgList(gx, px), RooArgList(nsig, nbkg))
# C r e a t e m o d e l f o r c o n t r o l s a m p l e
# --------------------------------------------------------------
# Construct signal pdf.
# NOTE that sigma is shared with the signal sample model
y = RooRealVar("y", "y", 40, 200)
mean_ctl = RooRealVar("mean_ctl", "mean_ctl", mu2, 40, 200)
sigma_ctl = RooRealVar("sigma", "sigma", sigma2, 0.1, 10)
gx_ctl = RooGaussian("gx_ctl", "gx_ctl", y, mean_ctl, sigma_ctl)
# Construct the background pdf
mean_bkg_ctl = RooRealVar("mean_bkg_ctl", "mean_bkg_ctl", mu1, 40, 200)
sigma_bkg_ctl = RooRealVar("sigma_bkg_ctl", "sigma_bkg_ctl", sigma1, 0.1, 20)
px_ctl = RooGaussian("px_ctl", "px_ctl", y, mean_bkg_ctl, sigma_bkg_ctl)
# Construct the composite model
# f_ctl = RooRealVar( "f_ctl", "f_ctl", 0.5, 0., 20. )
model_ctl = RooAddPdf("model_ctl", "model_ctl", RooArgList(gx_ctl, px_ctl), RooArgList(nsig, nbkg))
# G e t e v e n t s f o r b o t h s a m p l e s
# ---------------------------------------------------------------
real_data, real_data_ctl = get_data()
real_data_hist = RooDataHist("real_data_hist", "real_data_hist", RooArgList(x), real_data)
real_data_ctl_hist = RooDataHist("real_data_ctl_hist", "real_data_ctl_hist", RooArgList(y), real_data_ctl)
input_hists = MapStrRootPtr()
input_hists.insert(StrHist("physics", real_data))
input_hists.insert(StrHist("control", real_data_ctl))
# C r e a t e i n d e x c a t e g o r y a n d j o i n s a m p l e s
# ---------------------------------------------------------------------------
# Define category to distinguish physics and control samples events
sample = RooCategory("sample", "sample")
sample.defineType("physics")
sample.defineType("control")
# Construct combined dataset in (x,sample)
combData = RooDataHist("combData", "combined data", RooArgList(x), sample, input_hists)
# C o n s t r u c t a s i m u l t a n e o u s p d f i n ( x , s a m p l e )
# -----------------------------------------------------------------------------------
# Construct a simultaneous pdf using category sample as index
simPdf = RooSimultaneous("simPdf", "simultaneous pdf", sample)
# Associate model with the physics state and model_ctl with the control state
simPdf.addPdf(model, "physics")
simPdf.addPdf(model_ctl, "control")
# P e r f o r m a s i m u l t a n e o u s f i t
# ---------------------------------------------------
model.fitTo(real_data_hist)
summary = "fit in signal region\n"
summary += "nsig: " + str(nsig.getValV()) + " +- " + str(nsig.getError()) + "\n"
summary += "nbkg: " + str(nbkg.getValV()) + " +- " + str(nbkg.getError()) + "\n"
model_ctl.fitTo(real_data_ctl_hist)
summary += "fit in control region\n"
summary += "nsig: " + str(nsig.getValV()) + " +- " + str(nsig.getError()) + "\n"
summary += "nbkg: " + str(nbkg.getValV()) + " +- " + str(nbkg.getError()) + "\n"
# Perform simultaneous fit of model to data and model_ctl to data_ctl
simPdf.fitTo(combData)
summary += "Combined fit\n"
summary += "nsig: " + str(nsig.getValV()) + " +- " + str(nsig.getError()) + "\n"
summary += "nbkg: " + str(nbkg.getValV()) + " +- " + str(nbkg.getError()) + "\n"
# P l o t m o d e l s l i c e s o n d a t a s l i c e s
# ----------------------------------------------------------------
# Make a frame for the physics sample
frame1 = x.frame(RooFit.Bins(30), RooFit.Title("Physics sample"))
# Plot all data tagged as physics sample
combData.plotOn(frame1, RooFit.Cut("sample==sample::physics"))
# Plot "physics" slice of simultaneous pdf.
# NBL You _must_ project the sample index category with data using ProjWData
# as a RooSimultaneous makes no prediction on the shape in the index category
#.........这里部分代码省略.........
示例15: fit_gau2_che
# 需要导入模块: from ROOT import RooAddPdf [as 别名]
# 或者: from ROOT.RooAddPdf import fitTo [as 别名]
def fit_gau2_che(var, dataset, title='', print_pars=False, test=False,
mean_=None, sigma_=None, sigma1_=None, sigmaFraction_=None):
# define background
c0 = RooRealVar('c0', 'constant', 0.1, -1, 1)
c1 = RooRealVar('c1', 'linear', 0.6, -1, 1)
c2 = RooRealVar('c2', 'quadratic', 0.1, -1, 1)
c3 = RooRealVar('c3', 'c3', 0.1, -1, 1)
bkg = RooChebychev('bkg', 'background pdf', var,
RooArgList(c0, c1, c2, c3))
# define signal
val = 5.28
dmean = 0.05
valL = val - dmean
valR = val + dmean
if mean_ is None:
mean = RooRealVar("mean", "mean", val, valL, valR)
else:
mean = RooRealVar("mean", "mean", mean_)
val = 0.05
dmean = 0.02
valL = val - dmean
valR = val + dmean
if sigma_ is None:
sigma = RooRealVar('sigma', 'sigma', val, valL, valR)
else:
sigma = RooRealVar('sigma', 'sigma', sigma_)
if sigma1_ is None:
sigma1 = RooRealVar('sigma1', 'sigma1', val, valL, valR)
else:
sigma1 = RooRealVar('sigma1', 'sigma1', sigma1_)
peakGaus = RooGaussian("peakGaus", "peakGaus", var, mean, sigma)
peakGaus1 = RooGaussian("peakGaus1", "peakGaus1", var, mean, sigma1)
if sigmaFraction_ is None:
sigmaFraction = RooRealVar("sigmaFraction", "Sigma Fraction", 0.5, 0., 1.)
else:
sigmaFraction = RooRealVar("sigmaFraction", "Sigma Fraction", sigmaFraction_)
glist = RooArgList(peakGaus, peakGaus1)
peakG = RooAddPdf("peakG","peakG", glist, RooArgList(sigmaFraction))
listPeak = RooArgList("listPeak")
listPeak.add(peakG)
listPeak.add(bkg)
fbkg = 0.45
nEntries = dataset.numEntries()
val=(1-fbkg)* nEntries
listArea = RooArgList("listArea")
areaPeak = RooRealVar("areaPeak", "areaPeak", val, 0.,nEntries)
listArea.add(areaPeak)
nBkg = fbkg*nEntries
areaBkg = RooRealVar("areaBkg","areaBkg", nBkg, 0.,nEntries)
listArea.add(areaBkg)
model = RooAddPdf("model", "fit model", listPeak, listArea)
if not test:
fitres = model.fitTo(dataset, RooFit.Extended(kTRUE),
RooFit.Minos(kTRUE),RooFit.Save(kTRUE))
nbins = 35
frame = var.frame(nbins)
frame.GetXaxis().SetTitle("B^{0} mass (GeV/c^{2})")
frame.GetXaxis().CenterTitle()
frame.GetYaxis().CenterTitle()
frame.SetTitle(title)
mk_size = RooFit.MarkerSize(0.3)
mk_style = RooFit.MarkerStyle(kFullCircle)
dataset.plotOn(frame, mk_size, mk_style)
model.plotOn(frame)
as_bkg = RooArgSet(bkg)
cp_bkg = RooFit.Components(as_bkg)
line_style = RooFit.LineStyle(kDashed)
model.plotOn(frame, cp_bkg, line_style)
if print_pars:
fmt = RooFit.Format('NEU')
lyt = RooFit.Layout(0.65, 0.95, 0.92)
param = model.paramOn(frame, fmt, lyt)
param.getAttText().SetTextSize(0.02)
param.getAttText().SetTextFont(60)
frame.Draw()
#.........这里部分代码省略.........