Python RooRealVar.setRange方法代码示例

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
def NEvtsCalc( uPars, uEPars):
	Pars 			= {}
	Pars[0]			= RooRealVar('c1', 'Exponential constant', -1, -10, 0)
	Pars[1] 		= RooRealVar('ExpY', 'Background Yield', 100, 0, 10000000)
	Pars[2]			= RooRealVar('Mean', 'Voigtian Mean' , 90.0, 20, 180.0)
	Pars[3]			= RooRealVar('Width', 'Voigtian Width' , 5.0, 0.5, 40.0)
	Pars[4]			= RooRealVar('Sigma', 'Voigtian Sigma' , 5.0, 0.5, 40.0)
	Pars[5]			= RooRealVar('VoY', 'Signal Yield', 100, 0, 10000000)
	if len(Pars) != len(uPars):
		print 'The input array has a weird number of entries...'
		return 0
	if len(uPars) != len(uEPars):
		print 'The input arrays are of different sizes...'
		return 0
	for x in xrange(0, 6):
		Pars[x].setVal(uPars[x])
		Pars[x].setError(uEPars[x])
	v			= RooRealVar('v', 'Invariant Mass (GeV)', 60, 120)
	Voigt 			= RooVoigtian('Voigt', 'Voigtian - Signal', v, Pars[2], Pars[3], Pars[4])
## Calculate integral from -2sigma to 2sigma
	VStDev_r		= Voigt.sigma(v)
	VStDev			= VStDev_r.getVal()
	v.setRange("sobRange", Pars[2].getVal() - 2.*VStDev, Pars[2].getVal() + 2.*VStDev)
	integralSig     	= Voigt.createIntegral(RooArgSet(v), RooFit.NormSet(RooArgSet(v)), RooFit.Range("sobRange"))
	FinalNumber		= integralSig.getVal()*Pars[5].getVal()
	return FinalNumber

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
def doMCFit(dataSet, x_var, addTitlePlot=''):  
    
    cuts_str = ''
    data = dataSet.reduce( RooFit.Cut(cuts_str) )

    x=RooRealVar(x_var, 'm_{#tau}',1757,1797,'MeV')
    numBins = 100 # define here so that if I change it also the ndof change accordingly
    x.setBins(numBins)

    ######################################################
    # DEFINE PDF
    ######################################################
    
    w = RooWorkspace('w')
    getattr(w,'import')(x)
    w.factory('''RooDSCBShape::DSCB({0},
    #mu[1777, 1760,1790],
    #sigma[5,0,10],
    #alpha[1.2], n[50, 1, 150],
    #alpha, n
    )'''.format(x_var))
    #w.var('n').setConstant(False)
    signal = w.pdf('DSCB')
    # w.factory('''RooGaussian::GG({0},
    # #mu, #sigma
    # )'''.format(x_var))
    # signal = w.pdf('GG')
  
    # Fit
    fit_region = x.setRange('fit_region',1757,1797)
    result = signal.fitTo(dataSet, RooFit.Save(), RooFit.Range('fit_region'))

    # Frame
    frame = x.frame(RooFit.Title(' Combined mass KK#mu '+addTitlePlot))
    dataSet.plotOn(frame)
    signal.plotOn(frame, RooFit.LineWidth(2))

    # Legends
    signal.paramOn(frame, RooFit.Layout(0.6,0.9,0.9))
    chi2 = round(frame.chiSquare(),2)
    leg = TLegend(0.3,0,.10,.10)
    leg.SetBorderSize(0)
    leg.SetFillStyle(0)
    leg.AddEntry(0,'#chi^{2} ='+str(chi2),'')
    frame.addObject(leg)

    c1 = TCanvas('c1', 'c1')
    frame.Draw()
    c1.Update()

    for prm in ('#mu', '#sigma', '#alpha', 'n'): # TODO: automatize finding of variables from the function
        w.var(prm).setConstant()
    
    return w, c1

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
# gauss 10
mean10 = RooRealVar("mean10","mean of gaussian",3,-10,10)
sigma10 = RooRealVar("sigma10","width of gaussian",2,0.1,10)

gauss10 = RooGaussian("gauss10","gaussian PDF",x,mean10,sigma10)

# generate to MC 
data3 = gauss10.generate(RooArgSet(x),1000)

# split to half dataset
data3_half  = RooDataSet("data3_half", "data3 half", RooArgSet(x), RooFit.Import( data3 ), RooFit.Cut("x<3") )

# set range

x.setRange("whole_range",-10 ,10 )
x.setRange("right_gaussian_range",3 ,10 )
x.setRange("left_gaussian_range",-10 ,3 )

# another PDF gauss 11
mean11 = RooRealVar("mean11","mean of gaussian",2,-10,10)
sigma11 = RooRealVar("sigma11","width of gaussian",3,0.1,10)

gauss11 = RooGaussian("gauss11","gaussian PDF",x ,mean11 ,sigma11)

# fit partial range
#gauss11.fitTo(data3_half,RooFit.Range("left_gaussian_range"))

# convert to ext PDF then fit
nGauss11 = RooRealVar("nGauss11", "Gauss 11 ext pdf",750 , 0., 2000)
gauss11_ext = RooExtendPdf("gauss11_ext", "extended p.d.f", gauss11, nGauss11)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
class Fit2L:
    def __init__(self, file_names, br_name="mass"):
        self.file_names = file_names
        self.br_name = br_name
        min_x = 8
        max_x = 12
        self.obs = RooRealVar("obs", "m4l", min_x, max_x)
        self.nbins = int((max_x - min_x) * 20)
        self.obs.setBin(self.nbins)
        # self.ws = ROOT.RooWorkspace("combined", "combined")
        self.obs.setRange("fit", 8.5, 11.5)
        self.dummy_hists = []
        self.chi2_cut = 10000

    def build_model(self):
        # mean = RooRealVar("mean", "mass of 1S", 9.46, 9.2, 9.7)
        # sigma = RooRealVar("sigma", "sigma of gaussian", 0.14, 0.09, 0.3)
        mean = RooRealVar("mean", "mass of 1S", 9.48352)
        sigma = RooRealVar("sigma", "sigma of gaussian", 1.38574e-01)
        gaussian = ROOT.RooGaussian("gauss", "gauss", self.obs, mean, sigma)

        ## try Crystal Ball
        # alpha = RooRealVar("alpha", "alpha of CB", 5.9, 0, 100)
        # cb_n = RooRealVar("cb_n", "n of CB",  55.4, 0, 100)
        # gaussian = ROOT.RooCBShape("gauss", "gauss", self.obs, mean, sigma, alpha, cb_n)

        n_sig = RooRealVar("n_sig", "number of signal", 5000, 0, 100000)
        esig = ROOT.RooExtendPdf("esig", "esig", gaussian, n_sig)

        m2_shift = RooRealVar("m2_shift", "m2 shift", 0.56296)
        m2 = ROOT.RooFormulaVar("m2", "mass of 2S", "@0+ @1", RooArgList(mean, m2_shift))
        s2 = ROOT.RooFormulaVar("s2", "sigma*(unit+m2_shift)/mean", "@0*([email protected]/9.46)", RooArgList(sigma, m2_shift))
        g2 = ROOT.RooGaussian("g2", "gauss", self.obs, m2, s2)
        n2 = RooRealVar("n2", "number of 2S", 800, 100, 100000)
        # n2 = ROOT.RooFormulaVar("n2", "number of 2S" ,"n_sig*0.26", RooArgList(n_sig))
        esig2 = ROOT.RooExtendPdf("esig2", "esig2", g2, n2)

        m3_shift = RooRealVar("m3_shift", "m3 shift", 0.8949)
        m3 = ROOT.RooFormulaVar("m3", "mass of 3S", "mean + m3_shift", RooArgList(mean, m3_shift))
        s3 = ROOT.RooFormulaVar("s3", "sigma*(unit+m2_shift)/mean", "@0*([email protected]/9.46)", RooArgList(sigma, m3_shift))
        g3 = ROOT.RooGaussian("g3", "gauss", self.obs, m3, s3)
        # n3 = RooRealVar("n3", "number of 3S" , 50, 0, 1000)
        n3 = ROOT.RooFormulaVar("n3", "number of 3S", "n2*0.45", RooArgList(n2))
        esig3 = ROOT.RooExtendPdf("esig3", "esig3", g3, n3)

        n_bkg = RooRealVar("n_bkg", "number of bkg", 1000, 0, 1e6)
        # p0 = RooRealVar("p0", "p0", 5.8677e-02, -1E6, 1E6)
        # p1 = RooRealVar("p1", "p1", -5.086E-02, -1E6, 1E6)
        # p2 = RooRealVar("p2", "p2", -1.96e-02, -1E6, 1E6)
        # p3 = RooRealVar("p3", "p3", -1.08E-02, -1E6, 1E6)
        # p4 = RooRealVar("p4", "p4", -1.55E-02, -1E6, 1E6)

        # best fitted value using 2015+2016
        p0 = RooRealVar("p0", "p0", 5.86772e-02)
        p1 = RooRealVar("p1", "p1", -5.08695e-02)
        p2 = RooRealVar("p2", "p2", -1.95595e-02)
        p3 = RooRealVar("p3", "p3", -1.08105e-02)
        p4 = RooRealVar("p4", "p4", -1.54855e-02)
        bkg = ROOT.RooChebychev("bkg", "bkg", self.obs, RooArgList(p0, p1, p2, p3, p4))
        # bkg = ROOT.RooPolynomial("bkg", "bkg", self.obs, RooArgList(p0, p1, p2))
        ebkg = ROOT.RooExtendPdf("ebkg", "ebkg", bkg, n_bkg)
        model = ROOT.RooAddPdf("model", "model", RooArgList(esig, esig2, esig3, ebkg))
        getattr(self.ws, "import")(model)

    def get_data(self):
        """
        create RooDataSet from file_names
        """
        tree = ROOT.TChain("upsilon", "upsilon")
        if type(self.file_names) is list:
            for file_name in self.file_names:
                tree.Add(file_name)
        else:
            tree.Add(self.file_names)

        nentries = tree.GetEntries()
        print "total: ", nentries
        obs_set = RooArgSet(self.obs)
        data = ROOT.RooDataSet("data", "data", obs_set)
        for ientry in xrange(nentries):
            tree.GetEntry(ientry)
            for i, m4l in enumerate(getattr(tree, self.br_name)):
                # m4l = m4l/1000
                if m4l > self.obs.getMax() or m4l < self.obs.getMin():
                    continue
                if tree.chi2[i] > self.chi2_cut:
                    continue

                self.obs.setVal(m4l)
                data.add(obs_set)

        getattr(self.ws, "import")(data)

    def fit(self):
        if not hasattr(self, "ws"):
            self.ws = ROOT.RooWorkspace("combined", "combined")

        self.build_model()
        self.get_data()
        data = self.ws.obj("data")
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
                # apply onia pT cut
                pT = tree.pt[i]
                if not self.pass_onia_pt_cut(pT):
                    continue

                self.obs.setVal(m4l)
                data.add(obs_set)

        getattr(self.ws, "import")(data)
        print "selected events: ", data.sumEntries()
        fin.Close()

    def fit(self):
        print "my configuration:",self.get_ws_name()
        if not hasattr(self, "ws"):
            # try to look for the workspace and reused the data
            f1 = ROOT.TFile.Open(self.get_ws_name())
            self.ws = ROOT.RooWorkspace("combined", "combined")
            if f1:
                self.f1 = f1
                ws = f1.Get("combined")
                if ws:
                    getattr(self.ws, "import")(ws.obj('data'))
                if self.new_model:
                    print "building new models"
                    self.build_model()
                else:
                    getattr(self.ws, "import")(ws.obj(self.model_name))
            else:
                self.build_model()
                self.get_data()

        data = self.ws.obj("data")
        self.obs.setRange("fit_range", 8.2, 11.7)
        model = self.ws.obj(self.model_name)
        model.Print()
        print "total data:", data.sumEntries()
        nll = model.createNLL(data)

        #self.change_model()

        model.fitTo(data, ROOT.RooFit.Range("fit_range"))
        #model.fitTo(data)

        nll_uncondition = nll.getVal()
        self.ws.saveSnapshot("splusb", self.ws.allVars())

        if hasattr(self, "f1"):
            self.f1.Close()
        else:
            self.ws.writeToFile(self.get_ws_name())

    def plot(self, title):
        ## plot
        if not hasattr(self, "ws"):
            f1 = ROOT.TFile.Open(self.get_ws_name())
            self.ws = f1.Get("combined")

        self.ws.loadSnapshot("splusb")

        ## get number of signal and background in signal region
        obs = self.ws.var('obs')
        obs.setRange("signal", 9.2, 9.7)
        nsig = self.ws.obj("n_sig").getVal()
        nbkg = self.ws.obj("n_bkg").getVal()
        nall = self.ws.obj('model').expectedEvents(RooArgSet(obs))

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
            ws.pdf('ped').paramOn(xfped)

        c1 = TCanvas('c1', 'pedestal')
        xfped.Draw()

        makeHistPdf(ped_hist, ws, x)
        savePedWidth = ws.var('pedWidth').getVal()

    ws.Print()


dataTree.Draw('{0}>>sig_hist({1},{2:0.1f},{3:0.1f})'.format(HOTower,Nbins,
                                                            minMip,maxMip),
              sigCut, 'goff')
sig_hist = gDirectory.Get('sig_hist')
x.setRange(minMip, maxMip)

if (sig_hist.GetEntries > 0):
    ## fpar = array('d', [0.]*7)
    ## fparerr = array('d', [0.]*7)
    ## fpar[6] = pedRms

    ## chisqr = Double(0.)
    ## ndf = Long(0)

    ## if havePeds:
    ##     fit = langaupedfit(sig_hist, ped_hist, fpar, fparerr, chisqr, ndf)
    ##     #fit = preFitHisto(sig_hist, fpar, fparerr, chisqr, ndf)
    ## else:
    ##     fit = preFitHisto(sig_hist, fpar, fparerr, chisqr, ndf)

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

#.........这里部分代码省略.........
        a0 = RooRealVar('a0','a0',1.,-1.,1.) #,0.5,0.,1.)
        a1 = RooRealVar('a1','a1',0.1,-1.,1.) #-0.2,0.,1.)
        #a2 = RooRealVar('a2','a2',-0.1,1.,-1.)
        background = RooChebychev('background','Background',x,RooArgList(a0,a1))
        parameters.add(RooArgSet(a0, a1))
    else:
        parameters.add(RooArgSet(alpha, beta, q0))

    #together
    chibs = RooArgList(chib1_pdf,chib2_pdf,background)    
    
    # ndof
    floatPars = parameters.selectByAttrib("Constant",ROOT.kFALSE)
    ndof = numBins - floatPars.getSize() - 1

    # # Here I have as parameters N1, N2, and N_background
    # n_chib1 = RooRealVar("n_chib1","n_chib1",1250, 0, 50000)
    # n_chib2 =  RooRealVar("n_chib2","n_chib2",825, 0, 50000)
    # 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)

    # Here I have as parameters N_12, ratio_12, N_background
    n_chib = RooRealVar("n_chib","n_chib",2075, 0, 100000)
    ratio_21 = RooRealVar("ratio_21","ratio_21",0.6, 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)
    parameters.add(RooArgSet(n_chib1, n_chib2, n_background))
    modelPdf = RooAddPdf('ModelPdf', 'ModelPdf', chibs, ratio_list)
    
    print 'Fitting to data'
    fit_region = x.setRange("fit_region",9.7,10.1)
    result=modelPdf.fitTo(data,RooFit.Save(), RooFit.Range("fit_region"))
    
        
    # 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))
    chib1P_set = RooArgSet(chib1_pdf)
    modelPdf.plotOn(frame,RooFit.Components(chib1P_set), RooFit.LineColor(ROOT.kGreen+2), RooFit.LineStyle(2), RooFit.LineWidth(1))
    chib2P_set = RooArgSet(chib2_pdf)
    modelPdf.plotOn(frame, RooFit.Components(chib2P_set),RooFit.LineColor(ROOT.kRed), RooFit.LineStyle(2), RooFit.LineWidth(1))
    background_set =  RooArgSet(background)
    modelPdf.plotOn(frame,RooFit.Components(background_set), RooFit.LineColor(ROOT.kBlack), RooFit.LineStyle(2), RooFit.LineWidth(1))
    modelPdf.plotOn(frame, RooFit.LineWidth(2))
    frame.SetName("fit_resonance")  

    # Make numChib object
    numChib = NumChib(numChib=n_chib.getVal(), s_numChib=n_chib.getError(), ratio_21=ratio_21.getVal(), s_ratio_21=ratio_21.getError(), numBkg=n_background.getVal(), s_numBkg=n_background.getError(), corr_NB=result.correlation(n_chib, n_background),corr_NR=result.correlation(n_chib, ratio_21) , name='numChib'+output_suffix+ptBin_label,q0=q0.getVal(),s_q0=q0.getError(),alpha=alpha.getVal(),s_alpha=alpha.getError(), beta=beta.getVal(), s_beta=beta.getError(), chiSquare=frame.chiSquare())
    #numChib.saveToFile('numChib'+output_suffix+'.txt')

    if noPlots:

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
gSystem.ListLibraries()
fitresult = PDF.fitTo(dataset, RooFit.Optimize(0), RooFit.Strategy(2),
                      RooFit.Save(True), RooFit.NumCPU(1),
                      RooFit.SumW2Error(True), RooFit.Offset(True),
                      RooFit.Verbose(True))
fitresult.Print()
from helpers import FitStatus
print FitStatus(fitresult.status())


## Plot results
# # Use when debugging plots
# dataset = dataset.reduce(RooFit.EventRange(0, 100))

# reduce precision otherwise plotting doesn't work
time.setRange('fullrange', epsilon, 15.0)
RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-5)
RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-5)

dkcatset = RooArgSet(decaycat)
tframe1 = time.frame(RooFit.Range(0., time.getMax()), RooFit.Name('ptime'),
                     RooFit.Title('Ds#pi - powerlaw, DsK - powerlaw * %s ratio' %
                                  ratiofn))
dsetlist[0].plotOn(tframe1, RooFit.MarkerStyle(kOpenTriangleDown))
DsPi_Model.plotOn(tframe1, RooFit.LineColor(kBlue))
dsetlist[1].plotOn(tframe1, RooFit.MarkerStyle(kFullTriangleUp))
DsK_Model.plotOn(tframe1, RooFit.LineColor(kBlue+2))
# PDF.plotOn(tframe1, RooFit.Slice(decaycat, 'DsPi'),
#            RooFit.ProjWData(dkcatset, dataset, True),
#            RooFit.LineColor(kBlue))
# PDF.plotOn(tframe1, RooFit.Slice(decaycat, 'DsK'),

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

#.........这里部分代码省略.........
    print "  read tree:", treeName, "and trigger:", triName
    if ALTERNATIVE: print "  using ALTERNATIVE fit functions"
    print "-"*11*2
    
    # Silent RooFit
    RooMsgService.instance().setGlobalKillBelow(RooFit.FATAL)
    
    #*******************************************************#
    #                                                       #
    #              Variables and selections                 #
    #                                                       #
    #*******************************************************#
    
    # Define all the variables from the trees that will be used in the cuts and fits
    # this steps actually perform a "projection" of the entire tree on the variables in thei ranges, so be careful once setting the limits
    X_mass = RooRealVar(  massVar, "m_{X}" if nLept > 0 else "m_{T}^{X}", XBINMIN, XBINMAX, "GeV")
    J_mass = RooRealVar( "fatjet1_prunedMassCorr",       "corrected pruned mass", HBINMIN, HBINMAX, "GeV")
    CSV1 = RooRealVar(   "fatjet1_CSVR1",                           "",        -1.e99,   1.e4     )
    CSV2 = RooRealVar(   "fatjet1_CSVR2",                           "",        -1.e99,   1.e4     )
    nBtag = RooRealVar(  "fatjet1_nBtag",                           "",            0.,   4        )
    CSVTop = RooRealVar( "bjet1_CSVR",                              "",        -1.e99,   1.e4     )
    isZtoEE = RooRealVar("isZtoEE",                                 "",            0.,   2        )
    isZtoMM = RooRealVar("isZtoMM",                                 "",            0.,   2        )
    isWtoEN = RooRealVar("isWtoEN",                                 "",            0.,   2        )
    isWtoMN = RooRealVar("isWtoMN",                                 "",            0.,   2        )
    weight = RooRealVar( "eventWeightLumi",                         "",         -1.e9,   1.       )
    
    # Define the RooArgSet which will include all the variables defined before
    # there is a maximum of 9 variables in the declaration, so the others need to be added with 'add'
    variables = RooArgSet(X_mass, J_mass, CSV1, CSV2, nBtag, CSVTop)
    variables.add(RooArgSet(isZtoEE, isZtoMM, isWtoEN, isWtoMN, weight))
    
    # Define the ranges in fatJetMass - these will be used to define SB and SR
    J_mass.setRange("LSBrange", LOWMIN, LOWMAX)
    J_mass.setRange("HSBrange", HIGMIN, HIGMAX)
    J_mass.setRange("VRrange",  LOWMAX, SIGMIN)
    J_mass.setRange("SRrange",  SIGMIN, SIGMAX)
    J_mass.setBins(54)
    
    # Define the selection for the various categories (base + SR / LSBcut / HSBcut )
    baseCut = leptCut + " && " + btagCut + "&&" + topVeto
    massCut = massVar + ">%d" % XBINMIN
    baseCut += " && " + massCut
    
    # Cuts
    SRcut  = baseCut + " && %s>%d && %s<%d" % (J_mass.GetName(), SIGMIN, J_mass.GetName(), SIGMAX)
    LSBcut = baseCut + " && %s>%d && %s<%d" % (J_mass.GetName(), LOWMIN, J_mass.GetName(), LOWMAX)
    HSBcut = baseCut + " && %s>%d && %s<%d" % (J_mass.GetName(), HIGMIN, J_mass.GetName(), HIGMAX)
    SBcut  = baseCut + " && ((%s>%d && %s<%d) || (%s>%d && %s<%d))" % (J_mass.GetName(), LOWMIN, J_mass.GetName(), LOWMAX, J_mass.GetName(), HIGMIN, J_mass.GetName(), HIGMAX)
    VRcut  = baseCut + " && %s>%d && %s<%d" % (J_mass.GetName(), LOWMAX, J_mass.GetName(), SIGMIN)
    
    # Binning
    binsJmass = RooBinning(HBINMIN, HBINMAX)
    binsJmass.addUniform(HBINS, HBINMIN, HBINMAX)
    binsXmass = RooBinning(XBINMIN, XBINMAX)
    binsXmass.addUniform(binFact*XBINS, XBINMIN, XBINMAX)
    
    #*******************************************************#
    #                                                       #
    #                      Input files                      #
    #                                                       #
    #*******************************************************#
    
    # Import the files using TChains (separately for the bkg "classes" that we want to describe: here DY and VV+ST+TT)
    treeData = TChain(treeName)
    treeMC   = TChain(treeName)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
def doMCFit(inputfile_name, mass_chib, cuts, output_name='ChiB',  plotTitle = "#Chi_{b}", fittedVariable='qValue', returnOnlyFitResult = False, useOtherSignalParametrization = False, drawPulls = False, legendOnPlot=True): 

    mass_error = 0.15 
 
    print "Creating DataSet from file "+str(inputfile_name)
    dataSet = makeRooDataset(inputfile_name)

    if(fittedVariable == 'refittedMass'):
        x_var = 'rf1S_chib_mass'
        output_suffix = '_refit'
        x_axis_label= 'm_{#mu^{+} #mu^{-} #gamma} [GeV]'
    else:
        x_var = 'invm1S'
        output_suffix = '_qValue'
        x_axis_label = 'm_{#gamma #mu^{+} #mu^{-}} - m_{#mu^{+} #mu^{-}} + m^{PDG}_{#Upsilon}  [GeV]'

    cuts_str = str(cuts)
    #cuts_str = quality_cut +"photon_pt > 0.5 && abs(photon_eta) < 1.0 &&  abs(dimuon_rapidity) < 1.0 && dimuon_pt>5.0 && pi0_abs_mass > 0.025 &&  fabs(dz) < 0.1 "#&& numPrimaryVertices < 16"
    data = dataSet.reduce( RooFit.Cut(cuts_str) )
    
    print 'Creating pdf'
    x=RooRealVar(x_var,'m(#mu #mu #gamma) - m(#mu #mu) + m_{#Upsilon}',9.8,9.96)#9.7,10.1,'GeV')
    numBins = 32 # define here so that if I change it also the ndof change accordingly
    x.setBins(numBins)

    
    # Double sided Crystal Ball
    mean=RooRealVar("#mu","mean ChiB",mass_chib, mass_chib-mass_error,mass_chib+mass_error,"GeV")
    sigma=RooRealVar("#sigma","sigma ChiB",0.006, 0,0.2,'GeV')
    a1 = RooRealVar('#alpha1', '#alpha1', 0.75, 0, 3)
    a2 = RooRealVar('#alpha2', '#alpha2', 1.6, 0, 3)
    n1 = RooRealVar('n1', 'n1', 2.8)#, 1.0,4.0) # 2 per 2S
    n2 = RooRealVar('n2', 'n2', 3)#, 1.,4.0) # 1 per 2S
    parameters = RooArgSet(mean, sigma, a1, a2, n1, n2)
    
    cb_pdf = My_double_CB('chib', 'chib', x, mean, sigma, a1, n1, a2, n2)
    #cb_pdf = RooCBShape('chib', 'chib', x, mean, sigma, a1, n1)

    # ndof
    floatPars = parameters.selectByAttrib("Constant",ROOT.kFALSE)
    ndof = numBins - floatPars.getSize() - 1

    if useOtherSignalParametrization: # In this case I redefine cb_pdf
        n1 = RooRealVar('n1', 'n1', 2)
        cb = RooCBShape('cb1', 'cb1', x, mean, sigma, a1, n1)
        # I use a2 as the sigma of my second CB  
        a2 = RooRealVar('#alpha2', '#alpha2', 0.5, 0, 3)
        gauss = RooCBShape('cb2', 'cb2',x, mean, a2, a1, n1)
        # I use n2 as the ratio of cb1 with respect to cb2 
        n2 = RooRealVar('n2', 'n2',0.,1.)
        cb_pdf = RooAddPdf('chib','chib',RooArgList(cb, gauss),RooArgList(n2))
        #cb_pdf = cb
    
    print 'Fitting to data'
    fit_region = x.setRange("fit_region",9.8,9.96)
    result = cb_pdf.fitTo(data, RooFit.Save(), RooFit.Range("fit_region"))
    
    # a1_val = a1.getVal()
    # a2_val = a2.getVal()
    # 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),'')
#.........这里部分代码省略.........

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

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
linesigup.Draw()


# In[14]:


plotname = "plots/phi_" + histname
c.SaveAs(plotname  + ".png")
c.SaveAs(plotname  + ".eps")
c.SaveAs(plotname  + ".pdf")


# In[15]:


mass.setRange("signalrange",signallow,signalup)
mass.setRange("sideleftrange",leftlowside,leftupside)
mass.setRange("siderightrange",rightlowside,rightupside)
signalIntegralBkg = cheb.analyticalIntegral(cheb.getAnalyticalIntegral(RooArgSet(mass),RooArgSet(mass)),"signalrange")
leftsideIntegralBkg = cheb.analyticalIntegral(cheb.getAnalyticalIntegral(RooArgSet(mass),RooArgSet(mass)),"sideleftrange")
rightsideIntegralBkg = cheb.analyticalIntegral(cheb.getAnalyticalIntegral(RooArgSet(mass),RooArgSet(mass)),"siderightrange")


# In[16]:


totIntegralBkg= cheb.analyticalIntegral(cheb.getAnalyticalIntegral(RooArgSet(mass),RooArgSet(mass)))


# In[17]:

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
__all__ = ['model', 'variables']

from ROOT import RooRealVar, RooFormulaVar, RooExponential, RooArgList, RooArgSet, RooAddPdf, RooGaussian, RooProdPdf

mass_b = RooRealVar('B_M', 'B_M', 5200, 5450)
mass_b.setRange('leftband', 5200, 5279 - 50)
mass_b.setRange('signal', 5279 - 50, 5279 + 50)
mass_b.setRange('rightband', 5279 + 50, 5450)

mass_d = RooRealVar('D~0_M', 'D~0_M', 1750, 2000)

bkg_bmass_trans = RooFormulaVar('bkg_bmass_trans', 'bkg_bmass_trans', 'B_M - 5200', RooArgList(mass_b))
bkg_dmass_trans = RooFormulaVar('bkg_dmass_trans', 'bkg_dmass_trans', 'D~0_M - 5200', RooArgList(mass_d))
prompt_dmass_trans = RooFormulaVar('prompt_bmass_trans', 'prompt_bmass_trans', 'D~0_M - 5200', RooArgList(mass_d))

# Parameters
bmass_mu = RooRealVar('sig_bmass_mu', 'sig_bmass_mu', 5280, 5200, 5400)
sig_bmass_sigma = RooRealVar('sig_bmass_sigma', 'sig_bmass_sigma', 20, 10, 100)
sig_dmass_mu = RooRealVar('sig_dmass_mu', 'sig_dmass_mu', 1870, 1750, 2000)
sig_dmass_sigma = RooRealVar('sig_dmass_sigma', 'sig_dmass_sigma', 20, 10, 100)

prompt_bmass_sigma = RooRealVar('prompt_bmass_sigma', 'prompt_bmass_sigma', 20, 10, 100)
prompt_dmass_lamb = RooRealVar('prompt_bmass_lamb', 'prompt_bmass_lamb', -1.0e-03, -1, 0)

bkg_bmass_lamb = RooRealVar('bkg_bmass_lamb', 'bkg_bmass_lamb', -1.84819e-03, -1, 0)
bkg_dmass_lamb = RooRealVar('bkg_dmass_lamb', 'bkg_dmass_lamb', -1e-3, -1, 0)

# Model
sig_bmass = RooGaussian('sig_bmass', 'sig_bmass', mass_b, bmass_mu, sig_bmass_sigma)
sig_dmass = RooGaussian('sig_dmass', 'sig_dmass', mass_d, sig_dmass_mu, sig_dmass_sigma)

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
def fitLandauGaus(hist, full = False):

    ## c1 = ROOT.TCanvas()
    ## c1.Divide(2)
    ## c1.cd(1)
    ## hist.Draw()
    neg_landau = False
    if hist.GetMean() < 0.:
        neg_landau = True
    if neg_landau:
        hist = turnHisto(hist)

    hist.Rebin(2)
    hist.SetTitle('')
    hist.SetName('hSignal')
    ## c1.cd(2)
    ## hist.Draw('hist')
    ## c1.SaveAs('foobar.pdf')

    ### #if neg_landau:
    ### #    func = ROOT.TF1('my_landau','[0] * TMath::Landau(-x,[1],[2])', hist.GetXaxis().GetXmin(), hist.GetXaxis().GetXmax())
    ### #    func.SetParameters(1, hist.GetMean(), hist.GetRMS() )
    ### #else:
    ### func = ROOT.TF1('my_landau','[0] * TMath::Landau(x,[1],[2])', hist.GetXaxis().GetXmin(), hist.GetXaxis().GetXmax())
    ### func.SetParameters(1, hist.GetMean(), hist.GetRMS() )

    ### hist.Fit('my_landau','q')
    ### fit_res = []
    ### fit_res.append(func.GetParameter(0) if not neg_landau else     func.GetParameter(0))
    ### fit_res.append(func.GetParameter(1) if not neg_landau else -1.*func.GetParameter(1))
    ### fit_res.append(func.GetParameter(2) if not neg_landau else     func.GetParameter(2))
    ### return hist, fit_res

## ROOFIT VERSION

    xmin = hist.GetXaxis().GetXmin()
    xmax = hist.GetXaxis().GetXmax()
    mean = hist.GetMean()
    mp = hist.GetXaxis().GetBinCenter(hist.GetMaximumBin())
    rms = hist.GetRMS()
    flandau = ROOT.TF1('flandau','landau',mp-20,mp+40)
    flandau.SetLineWidth(1)
    flandau.SetLineColor(ROOT.kBlue)
    hist2 = hist.Clone(hist.GetName()+'_2')
    hist2.Scale(1./hist2.GetBinContent(hist2.GetMaximumBin()))
    hist2.Fit(flandau,'Q','',mp-20,mp+40)

    flandau2 = flandau.Clone('flandau2')
    flandau2.SetRange(0,500)
    flandau2.SetLineStyle(2)

    for i in range(flandau.GetNpar()):
        flandau2.SetParLimits(i,flandau.GetParameter(i),flandau.GetParameter(i))
    hist2.Fit(flandau2,'Q+')#,'same',mp-20,mp+40)
    for i in range(flandau.GetNpar()):
        hist2.GetFunction('flandau2').SetParameter(i,flandau.GetParameter(i))

    for i in range(flandau.GetNpar()):
        print flandau.GetParameter(i),flandau2.GetParameter(i)

    x   = RooRealVar('x', 'signal / adc', 0,500)
    x.setRange("signal",mp - 40, mp+90)
    x.setRange("draw",0,500)
    ral = RooArgList(x)
    dh  = RooDataHist('dh', 'dh', ral, RooFit.Import(hist))
    
    
    if full: 
        ml     = RooRealVar('ml', 'mean landau' , mp, mp-20., mp+30)
        sl     = RooRealVar('sl', 'sigma landau', 10, 1., 25.)
    else:
        ml     = RooRealVar('ml', 'mean landau' , mean, mean-40., mean)
        sl     = RooRealVar('sl', 'sigma landau', 10., 6., 14.)
    landau = RooLandau ('lx', 'lx', x, ml, sl)
    
    mean = 0
    if full: 
        mg     = RooRealVar ('mg', 'mean gaus' , 0,0,0)
        sg     = RooRealVar ('sg', 'sigma gaus', flandau.GetParameter(2), 0.1, 30.)
    else:
        mg     = RooRealVar ('mg', 'mean gaus' , 0,0,0) #mean, mean-30.,  mean+30.)
        sg     = RooRealVar ('sg', 'sigma gaus', 2., 0.1, 20.)
    gaus   = RooGaussian('gx', 'gx', x, mg, sg)
    
    x.setBins(1000,'cache')
    
    ## Construct landau (x) gauss
    lxg = RooFFTConvPdf('lxg','landau (x) gaus', x, landau, gaus)
    lxg.fitTo(dh,RooFit.Range("signal"))
    #,RooFit.Normalization(ROOT.RooAbsReal.NumEvent,1))
    a = lxg.getParameters(dh)

    print 'fit par0                                     %+6.1f'%flandau.GetParameter(0)
    print 'fit par1                                     %+6.1f'%flandau.GetParameter(1)
    print 'fit par2                                     %+6.1f'%flandau.GetParameter(2)
    print 'mp                                           %+6.1f'%mp
    print 'rms                                          %+6.1f'%rms
    print 'lxg.getParameters(dh).getRealValue(\'ml\'):  %+6.1f'% a.getRealValue('ml')
    print 'lxg.getParameters(dh).getRealValue(\'sl\'):  %+6.1f'% a.getRealValue('sl')
    print 'lxg.getParameters(dh).getRealValue(\'sg\'):  %+6.1f'% a.getRealValue('sg')
#.........这里部分代码省略.........

# 需要导入模块: from ROOT import RooRealVar [as 别名]
# 或者: from ROOT.RooRealVar import setRange [as 别名]
from ROOT import RooAbsPdf, RooGaussian, RooUniform
from ROOT import RooGenericPdf, RooEffProd, RooAddPdf, RooProdPdf, RooHistPdf
from ROOT import RooDataSet, RooDataHist
from ROOT import RooDecay, RooGaussModel

# my stuff
from factory import *

set_integrator_config()

# epsilon = sys.float_info.epsilon # python -> C++ doesn't like this
epsilon = 2E-4

# Observables
time = RooRealVar('time', 'B_{s} lifetime in ns', epsilon, 0.01+epsilon)
time.setRange('fullrange', epsilon, 0.01+epsilon)
# Limits determined from tree
dt = RooRealVar('dt', 'Error in lifetime measurement (ns)', 1E-5, 9E-5)
dt.setBins(100)

# Temporary RooArgSet to circumvent scoping issues for nested
# temporary objects.
timeargset = RooArgSet(time)
dtargset = RooArgSet(dt)

# Parameters
if accfn == 'powerlaw':
    turnon = RooRealVar('turnon', 'turnon', 1500., 500., 5000.)
    offset = RooRealVar('offset', 'offset', 0., -1E-3, 1E-3)
elif accfn == 'arctan':
    # turnon has a different range as it is in the denominator