Python RooRealVar.getValV方法代码示例

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


from time import sleep
sleep(10)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getValV [as 别名]
    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())
    if nsigref.getValV()==0: continue

    if fit=="data":
      xframe=mass.frame(RooFit.Title("signal fraction in peak ="+str(int(nsigref.getValV()/(nsigref.getValV()+nbkg.getValV())*1000.)/1000.)+"+-"+str(int(nsigref.getError()/(nsigref.getValV()+nbkg.getValV())*1000.)/1000.)+", #chi^{2}/DOF = "+str(int((chi2.getVal()/(nbins-nfree))*10.)/10.)))
      signal.plotOn(xframe,RooFit.DataError(RooAbsData.SumW2))
      sigmodel.plotOn(xframe,RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      sigmodel.plotOn(xframe,RooFit.Components("sigbkg"+str(cut)),RooFit.LineStyle(kDashed),RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      sigmodel.plotOn(xframe,RooFit.Components("sig"),RooFit.LineStyle(kDotted),RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))
      canvas=TCanvas("c3","c3",0,0,600,600)
      xframe.Draw()
      canvas.SaveAs(prefix+"_"+plot[0]+str(cut)+"_sigfit.pdf")

    a0=RooRealVar("a0"+str(cut),"a0",100.,0.,1000.)
    a1=RooRealVar("a1"+str(cut),"a1",100.,0.,1000.)
    a2=RooRealVar("a2"+str(cut),"a2",50.,0.,1000.)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getValV [as 别名]
datacard.write('jmax 1\n')
datacard.write('kmax *\n')
datacard.write('---------------\n')
datacard.write('shapes * * '+wsFN+' 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')      
#nuisance parameters --- gaussian prior
if bkgNuisance:
  datacard.write('p1  param    '+str(p1.getValV())+'   '+str(p1.getError())+'\n')
  datacard.write('p2  param    '+str(p2.getValV())+'   '+str(p2.getError())+'\n')
  datacard.write('p3  param    '+str(p3.getValV())+'   '+str(p3.getError())+'\n')
#flat parameters --- flat prior
if not (bkgNuisance or histpdfBkg): 
  datacard.write('background_norm  flatParam\n')
  datacard.write('p1  flatParam\n')
  datacard.write('p2  flatParam\n')
  datacard.write('p3  flatParam\n')


   
 
##----- keep the GUI alive ------------
#if __name__ == '__main__':
#  rep = ''

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

#.........这里部分代码省略.........
    # 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(), 
                        )
    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 
    # and can thus not be integrated
    simPdf.plotOn( frame1, RooFit.Slice( sample, "physics" ),
                   RooFit.ProjWData( RooArgSet( sample ), combData ), ) 
    simPdf.plotOn( frame1, RooFit.Slice( sample, "physics" ),
                   RooFit.Components( "signal_1_pdf" ),
                   RooFit.ProjWData( RooArgSet( sample ), combData ),

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

# In[10]:


tot.fitTo(dh)


# In[11]:


massFrame = mass.frame()
massFrame.SetTitle("Phi signal")
dh.plotOn(massFrame)
tot.plotOn(massFrame)
gauss.plotOn(massFrame,LineColor(kGreen),LineStyle(kDashed),Normalization((sFrac.getValV()*numEvts)/(numEvts)))
cheb.plotOn(massFrame,LineColor(kMagenta),LineStyle(kDotted),Normalization(((1.0-sFrac.getValV())*numEvts)/(numEvts)))
tot.paramOn(massFrame,Layout(0.60,0.99,0.75));
massFrame.Draw()


# In[12]:


plotmax = hist.GetMaximum()*1.05
sidesigma = sigma.getValV()
leftlowside = -7.*sidesigma + mean.getValV()
leftupside = -5.*sidesigma + mean.getValV()
rightlowside = +5.*sidesigma + mean.getValV()
rightupside = +7.*sidesigma + mean.getValV()

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getValV [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
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getValV [as 别名]
g_bp_pos_resolution_vs_p_5hit = TGraphErrors()

for bin in xrange(1, h_sp_residuals_bp_vs_p_5hit.GetXaxis().GetNbins()):

	projection = h_sp_residuals_bp_vs_p_5hit.ProjectionY("", bin, bin)
	
	if projection.GetEntries() < 100: continue

	sp_bp_residuals_data = RooDataHist("sp_bp_residuals_data", "Bend Plane Residuals Data", RooArgList(sp_bp_residuals), projection)
	sp_bp_residuals_plot = sp_bp_residuals.frame()
	sp_bp_residuals_data.plotOn(sp_bp_residuals_plot)
	
	resolution_model.fitTo(sp_bp_residuals_data)
	resolution_model.plotOn(sp_bp_residuals_plot)

	g_bp_pos_resolution_vs_p_5hit.SetPoint(bin-1, bin*.25, bp_core_sigma.getValV())
	g_bp_pos_resolution_vs_p_5hit.SetPointError(bin-1, 0, bp_core_sigma.getError())

	sp_bp_residuals_plot.Draw()
	canvas.Print("sp_residuals.pdf(")
	
	sp_bp_residuals_plot.remove()
	sp_bp_residuals_plot.remove()


g_bp_pos_resolution_vs_p_5hit.Draw("*Ae")
canvas.Print("sp_residuals.pdf(")

sp_nbp_residuals = RooRealVar("sp_residuals_nbp", "Non-Bend Plane Residuals (mm)", -20, 20)

nbp_core_mean = RooRealVar("nbp_core_mean", "nbp_core_mean", 0, -2, 2)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import getValV [as 别名]
      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")
      meanW.setConstant(True)
      widthWZ.setConstant(True) 
      sigmodelZ=RooAddPdf("sigmodelZ"+name,"sigZ",RooArgList(sigbkg,sigZ),RooArgList(nsigbkg,nsigrefZ))
      sigmodelZ.fitTo(signalZ,RooFit.SumW2Error(True))
      sigmodelZ.fitTo(signalZ,RooFit.SumW2Error(True))
      sigmodelZ.fitTo(signalZ,RooFit.SumW2Error(True))
      xframe=mass.frame(RooFit.Title("             #Delta m="+str(int((meanWZ.getValV())*1000.)/1000.)+"#pm"+str(int(meanWZ.getError()*1000.)/1000.)+" GeV"))
      signalZ.plotOn(xframe,RooFit.DataError(RooAbsData.SumW2))
      sigmodelZ.plotOn(xframe,RooFit.Normalization(1.0,RooAbsReal.RelativeExpected))