Python RooMsgService.instance方法代碼示例

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
    def model(self, ws, debug = 0):

        if self.sModel not in self.items:
            if debug>0:
                print self.legend, 'no Model section defined in the config file'
                print self.legend, 'no model will be created'
            return

        legend = '['+self.sModel+']:'

        # FIXME: find a better way to load custom PDFs
        #gROOT.ProcessLine(".L dijetQstarPdf.cxx+");
        #gROOT.ProcessLine(".L Qstar_qg.cxx+");
        #gROOT.ProcessLine(".L Qstar_qg_2.cxx+");
        #gROOT.ProcessLine(".L Jacobian_mt.cxx+");

        model_items = self.items[self.sModel]
        
        if (debug>0):
            print 'Loading model to the workspace', ws.GetName()

        # using RooStats::HLFactory to parse the model
        # make a temp file with a card file
        card_file_name = 'hlfCardFile.rs'
        with open(card_file_name, 'w') as cardfile:
            for item in model_items:
                print >> cardfile,  item[0], '=', item[1]

        hlfVerbosity = debug > 0
        #hlf = ROOT.RooStats.HLFactory('HLFactoryExost', card_file_name, True)
        hlf = RooStats.HLFactory('HLFactoryExost',
                                      ws,
                                      hlfVerbosity)

        self.modelAutoImportWarningExplanation(legend, debug)

        current_messaging_level = RooMsgService.instance().globalKillBelow() ;
        if debug <= 0:
            RooMsgService.instance().setGlobalKillBelow(RooFit.ERROR) ;

        hlf.ProcessCard(card_file_name)

        print legend, 'importing nonstandard class code into the workspace...'
        ws.importClassCode();
        #ws.importClassCode('Qstar_qg', true);
        print legend, '...code import done'

        # return messaging to the original level
        RooMsgService.instance().setGlobalKillBelow(current_messaging_level) ;

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
    def __init__(self, fit_data_collection, method="TMinuit"):
        RooMsgService.instance().setGlobalKillBelow(RooFit.FATAL)
        self.method = method
        self.logger = logging.getLogger("RooFit")
        self.constraints = fit_data_collection.constraints()
        self.constraint_type = ""
        self.saved_result = None

        self.fit_data_collection = fit_data_collection
        self.samples = fit_data_collection.mc_samples()
        self.normalisation = fit_data_collection.mc_normalisation()

        self.fit_data_1 = self.fit_data_collection.fit_data.items()[0][1]
        self.fit_boundaries = self.fit_data_1.fit_boundaries
        self.data_label = FitData.data_label
        self.histograms = self.fit_data_1.histograms_

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
def action(ws, action_name, settings = None):
    legend = '[Action]:'
    print legend, 'beginning the ', action_name, 'action'

    if action_name == 'bayes':

        if settings['status'] == 'fail':
            print legend, 'Bayesian calculator failed to configure'
            print legend, 'action stopped.'
            return

        legend = '[Bayesian calculator]:'

        #mconf = ws.obj('exostModelConfig')
        mconf = ws.obj(settings['model_config_name'])

        bCalc = ws.obj('exostBayes')

        # to suppress messages when pdf goes to zero
        RooMsgService.instance().setGlobalKillBelow(RooFit.FATAL)
  
        if mconf.GetParametersOfInterest().getSize() == 1:

            bInt = bCalc.GetInterval()


            cl = bCalc.ConfidenceLevel()
            print legend,  str(cl)+'% CL central interval: [', bInt.LowerLimit(), ' - ', bInt.UpperLimit(), ']'
            print legend, 'or', str(cl+(1.0-cl)/2), '% CL limits'

            # make a posterior plot
            if settings['do_posterior_plot']:

                plot = bCalc.GetPosteriorPlot()
                c1 = TCanvas('c1','Bayesian Calculator Result');
                c1.cd(1)
                plot.Draw()
                c1.SaveAs('bayesian_example_plot.'+settings['plot_format'])
  
        else:
            print legend, 'Error: Bayesian Calc. only supports one parameter of interest'

    print legend, 'the', action_name, 'action is done'

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
                  help='use keys pdf for background')
parser.add_option('-b', action='store_true', default=False, dest='b',
                  help='no x windows')
(opts, args) = parser.parse_args()

import pyroot_logon
from ROOT import gROOT
gROOT.ProcessLine('.L buildSimPdf.cc+')

from ROOT import readData,computeRatio,computeRatioError,buildPdf,\
     computeSumError,\
     RooWorkspace,RooFit,TCanvas,kRed,kGreen,kDashed,buildSimPdf,RooArgSet,\
     RooRealVar,RooMsgService, TMath, TFile, RooAbsReal
from math import sqrt

RooMsgService.instance().setGlobalKillBelow(RooFit.ERROR)

hidatafile = 'data/dimuonTree_150mub.root'
ppdatafile = 'data/dimuonTree_2011_pp.root'

mmin = 7.
mmax = 14.


cuts = '(muPlusPt > %0.1f) && (muMinusPt > %0.1f) && (abs(upsRapidity)<2.4) && (vProb > 0.05)' \
       % (opts.pt, opts.pt)
simparamfile = opts.paramfile
useKeys = opts.keys
## cuts = '(muPlusPt > 3.5) && (muMinusPt > 3.5) && (abs(upsRapidity)<2.4)'
## simparamfile = 'nom3.5SimFit.txt'

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
parser.add_option('-b', action='store_true', default=False, dest='b',
                  help='no x windows')
(opts, args) = parser.parse_args()

import pyroot_logon
from ROOT import gROOT
gROOT.ProcessLine('.L buildSimPdf.cc+')

from ROOT import readData,computeRatio,computeRatioError,buildPdf,\
     computeSumError,\
     RooWorkspace,RooFit,TCanvas,kRed,kGreen,kDashed,buildSimPdf,RooArgSet,\
     RooRealVar,RooMsgService, TMath, TFile, RooAbsReal, TGraph, \
     RooStats
from math import sqrt

RooMsgService.instance().setGlobalKillBelow(RooFit.FATAL)

hidatafile = 'data/dimuonTree_150mub.root'
ppdatafile = 'data/dimuonTree_2011_pp.root'

mmin = 7.
mmax = 14.

cuts = '(muPlusPt > %0.1f) && (muMinusPt > %0.1f) && (abs(upsRapidity)<2.4) && (vProb > 0.05)' \
       % (opts.pt, opts.pt)
simparamfile = opts.paramfile
useKeys = opts.keys
## cuts = '(muPlusPt > 3.5) && (muMinusPt > 3.5) && (abs(upsRapidity)<2.4)'
## simparamfile = 'nom3.5SimFit.txt'

ws = RooWorkspace("ws","ws")

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]

#.........這裏部分代碼省略.........
                            )
    mass_group.add_argument("--masslist",
                            help = "List containing mass information"
                            )

    args = parser.parse_args()

    # check if the output directory exists
    if not os.path.isdir( os.path.join(os.getcwd(),args.output_path) ):
        os.mkdir( os.path.join(os.getcwd(),args.output_path) )

    # mass points for which resonance shapes will be produced
    masses = []

    if args.massrange != None:
        MIN, MAX, STEP = args.massrange
        masses = range(MIN, MAX+STEP, STEP)
    elif args.masslist != None:
        # A mass list was provided
        print  "Will create mass list according to", args.masslist
        masslist = __import__(args.masslist.replace(".py",""))
        masses = masslist.masses
    else:
        masses = args.mass

    # sort masses
    masses.sort()

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

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

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

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

    sqrtS = args.sqrtS

    for mass in masses:

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

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

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

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

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

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
#import HWW2DConfig
config = __import__(opts.modeConfig)
import RooWjj2DFitter

from ROOT import TCanvas, RooFit, RooLinkedListIter, TMath, RooRandom, TFile, \
    RooDataHist, RooMsgService, TStopwatch, RooAbsPdf, RooAbsData, \
    RooWorkspace, RooArgList, RooAddPdf
import pulls

timer = TStopwatch()
timer.Start()

#RooAbsPdf.defaultIntegratorConfig().setEpsRel(1e-9)
#RooAbsPdf.defaultIntegratorConfig().setEpsAbs(1e-9)
if not opts.debug:
    RooMsgService.instance().setGlobalKillBelow(RooFit.WARNING)
    RooMsgService.instance().addStream(RooFit.ERROR,
                                       RooFit.Prefix(True), 
                                       RooFit.ClassName('RooExpPoly'),
                                       RooFit.OutputFile('/dev/null'))
    RooMsgService.instance().Print('v')

if hasattr(opts, "seed") and (opts.seed >= 0):
    print "random seed:", opts.seed
    RooRandom.randomGenerator().SetSeed(opts.seed)

mvaCutOverride = None
if hasattr(opts, "mvaCut"):
    mvaCutOverride = opts.mvaCut

mjjArgs = []

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
ntuple_file = None
if args[0] in input_data.keys():
    ntuple_file = input_data[args[0]]

from P2VV.Load import P2VVLibrary
from P2VV.Load import RooFitOutput
from P2VV.RooFitWrappers import *

from itertools import product
from ROOT import RooCBShape as CrystalBall
from P2VV.Parameterizations.GeneralUtils import valid_combinations
#from P2VV.Load import RooFitOutput

from ROOT import RooMsgService
RooMsgService.instance().addStream(RooFit.DEBUG,RooFit.Topic(RooFit.Generation))
## RooMsgService.instance().addStream(RooFit.DEBUG,RooFit.Topic(RooFit.Integration))

obj = RooObject( workspace = 'w')
w = obj.ws()

from math import pi
t = RealVar('time', Title = 'decay time', Unit='ps', Observable = True, MinMax=(0.3, 14))
m = RealVar('mass', Title = 'B mass', Unit = 'MeV', Observable = True, MinMax = (5250, 5550))
nPV = RealVar('nPV', Title = 'nPV', Observable = True, MinMax = (0, 15))
mpsi = RealVar('mdau1', Title = 'J/psi mass', Unit = 'MeV', Observable = True, MinMax = (3030, 3150))
st = RealVar('sigmat',Title = '#sigma(t)', Unit = 'ps', Observable = True, MinMax = (0.0001, 0.12))

# Categories
hlt1_biased = Category('hlt1_biased', States = {'biased' : 1, 'not_biased' : 0}, Observable = True)
hlt1_unbiased = Category('hlt1_unbiased_dec', States = {'unbiased' : 1, 'not_unbiased' : 0}, Observable = True)

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
# set constant what is supposed to be constant
setConstantIfSoConfigured(config, fitpdf)

# set up fitting options
fitopts = [ RooFit.Timer(), RooFit.Save(),
    RooFit.Strategy(config['FitConfig']['Strategy']),
    RooFit.Optimize(config['FitConfig']['Optimize']),
    RooFit.Offset(config['FitConfig']['Offset']),
    RooFit.NumCPU(config['FitConfig']['NumCPU']) ]

# set up blinding for data
fitopts.append(RooFit.Verbose(not (config['IsData'] and config['Blinding'])))
if config['IsData'] and config['Blinding']:
    from ROOT import RooMsgService
    RooMsgService.instance().setGlobalKillBelow(RooFit.WARNING)                                                                                                                             
    fitopts.append(RooFit.PrintLevel(-1))
fitOpts = RooLinkedList()
for o in fitopts: fitOpts.Add(o)

# fit
rawfitresult = fitpdf.fitTo(ds, fitOpts)

# pretty-print the result
from B2DXFitters.FitResult import getDsHBlindFitResult
result = getDsHBlindFitResult(config['IsData'], config['Blinding'],
    rawfitresult)
print result

# write raw fit result and workspace to separate ROOT files
from ROOT import TFile

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]

#.........這裏部分代碼省略.........
    mass_group.add_argument("--masslist",
                            help = "List containing mass information"
                            )

    args = parser.parse_args()

    fit_functions = args.fit_functions.split(",")

    # mass points for which resonance shapes will be produced
    masses = []

    if args.fitBonly:
        masses.append(750)
    else:
        if args.massrange != None:
            MIN, MAX, STEP = args.massrange
            masses = range(MIN, MAX+STEP, STEP)
        elif args.masslist != None:
            # A mass list was provided
            print  "Will create mass list according to", args.masslist
            masslist = __import__(args.masslist.replace(".py",""))
            masses = masslist.masses
        else:
            masses = args.mass

    # sort masses
    masses.sort()

    # import ROOT stuff
    from ROOT import gStyle, TFile, TH1F, TH1D, TGraph, kTRUE, kFALSE, TCanvas, TLegend, TPad, TLine
    from ROOT import RooHist, RooRealVar, RooDataHist, RooArgList, RooArgSet, RooAddPdf, RooProdPdf, RooEffProd, RooFit, RooGenericPdf, RooWorkspace, RooMsgService, RooHistPdf, RooExtendPdf

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

    # input data file
    #inputData = TFile(limit_config.get_data_input(args.analysis))
    # input data histogram
    #hData = inputData.Get(args.dataHistname)
    #hData.SetDirectory(0)
    data_file = TFile(analysis_config.get_b_histogram_filename(args.analysis, "BJetPlusX_2012"))
    hData = data_file.Get("BHistograms/h_pfjet_mjj")
    hData.SetDirectory(0)

    # input sig file
    if not args.fitSignal:
        print "[create_datacards] INFO : Opening resonance shapes file at " + limit_config.get_resonance_shapes(args.analysis, args.model)
        inputSig = TFile(limit_config.get_resonance_shapes(args.analysis, args.model), "READ")

    sqrtS = args.sqrtS

    # mass variable
    mjj = RooRealVar('mjj','mjj',float(args.massMin),float(args.massMax))

    # integrated luminosity and signal cross section
    lumi = args.lumi
    signalCrossSection = 1. # set to 1. so that the limit on r can be interpreted as a limit on the signal cross section

    if args.correctTrigger:
        trigger_efficiency_pdf = trigger_efficiency.get_pdf(args.analysis, mjj)
        trigger_efficiency_formula = trigger_efficiency.get_formula(args.analysis, mjj)
    else:
        trigger_efficiency_pdf = trigger_efficiency.get_trivial_pdf(mjj)
        trigger_efficiency_formula = trigger_efficiency.get_trivial_formula(mjj)

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
def alpha(channel):

    nElec = channel.count('e')
    nMuon = channel.count('m')
    nLept = nElec + nMuon
    nBtag = channel.count('b')
    
    # Channel-dependent settings
    # Background function. Semi-working options are: EXP, EXP2, EXPN, EXPTAIL
    if nLept == 0:
        treeName = 'SR'
        signName = 'XZh'
        colorVjet = sample['DYJetsToNuNu']['linecolor']
        triName = "HLT_PFMET"
        leptCut = "0==0"
        topVeto = selection["TopVetocut"]
        massVar = "X_cmass"
        binFact = 1
        #fitFunc = "EXP"
        #fitFunc = "EXP2"
        #fitFunc = "EXPN"
        #fitFunc = "EXPTAIL"
        fitFunc = "EXPN" if nBtag < 2 else "EXP"
        fitAltFunc = "EXPTAIL" if nBtag < 2 else "EXPTAIL"
        fitFuncVjet = "ERFEXP" if nBtag < 2 else "ERFEXP"
        fitFuncVV   = "EXPGAUS"
        fitFuncTop  = "GAUS2"
    elif nLept == 1:
        treeName = 'WCR'
        signName = 'XWh'
        colorVjet = sample['WJetsToLNu']['linecolor']
        triName = "HLT_Ele" if nElec > 0 else "HLT_Mu"
        leptCut = "isWtoEN" if nElec > 0 else "isWtoMN"
        topVeto = selection["TopVetocut"]
        massVar = "X_mass"
        binFact = 2
        if nElec > 0:
            fitFunc = "EXP" if nBtag < 2 else "EXP"
            fitAltFunc  = "EXPTAIL" if nBtag < 2 else "EXPTAIL"
        else:
            fitFunc = "EXPTAIL" if nBtag < 2 else "EXP"
            fitAltFunc  = "EXPN" if nBtag < 2 else "EXPTAIL"
        fitFuncVjet = "ERFEXP" if nBtag < 2 else "ERFEXP"
        fitFuncVV   = "EXPGAUS"
        fitFuncTop  = "GAUS3" if nBtag < 2 else "GAUS2"
    else:
        treeName = 'XZh'
        signName = 'XZh'
        colorVjet = sample['DYJetsToLL']['linecolor']
        triName = "HLT_Ele" if nElec > 0 else "HLT_Mu"
        leptCut = "isZtoEE" if nElec > 0 else "isZtoMM"
        topVeto = "0==0"
        massVar = "X_mass"
        binFact = 5
        if nElec > 0:
            fitFunc = "EXP" if nBtag < 2 else "EXP"
            fitAltFunc = "POW" if nBtag < 2 else "POW"
        else:
            fitFunc = "EXP" if nBtag < 2 else "EXP"
            fitAltFunc = "POW" if nBtag < 2 else "POW"
        fitFuncVjet = "ERFEXP" if nBtag < 2 else "EXP"
        fitFuncVV   = "EXPGAUS2"
        fitFuncTop  = "GAUS"
    
    btagCut = selection["2Btag"] if nBtag == 2 else selection["1Btag"]
    
    print "--- Channel", channel, "---"
    print "  number of electrons:", nElec, " muons:", nMuon, " b-tags:", nBtag
    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))
    
#.........這裏部分代碼省略.........

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
def alpha(channel):

    nElec = channel.count("e")
    nMuon = channel.count("m")
    nLept = nElec + nMuon
    nBtag = channel.count("b")

    # Channel-dependent settings
    # Background function. Semi-working options are: EXP, EXP2, EXPN, EXPTAIL
    if nLept == 0:
        treeName = "SR"
        signName = "XZh"
        colorVjet = sample["DYJetsToNuNu"]["linecolor"]
        triName = "HLT_PFMET"
        leptCut = "0==0"
        topVeto = selection["TopVetocut"]
        massVar = "X_cmass"
        binFact = 1
        fitFunc = "EXPN" if nBtag < 2 else "EXPN"
        fitAltFunc = "EXPTAIL" if nBtag < 2 else "EXPTAIL"
        fitFuncVjet = "ERFEXP" if nBtag < 2 else "EXP"
        fitAltFuncVjet = "POL" if nBtag < 2 else "POL"
        fitFuncVV = "EXPGAUS" if nBtag < 2 else "EXPGAUS"
        fitFuncTop = "GAUS2"
    elif nLept == 1:
        treeName = "WCR"
        signName = "XWh"
        colorVjet = sample["WJetsToLNu"]["linecolor"]
        triName = "HLT_Ele" if nElec > 0 else "HLT_Mu"
        leptCut = "isWtoEN" if nElec > 0 else "isWtoMN"
        topVeto = selection["TopVetocut"]
        massVar = "X_mass"
        binFact = 2
        if nElec > 0:
            fitFunc = "EXPTAIL" if nBtag < 2 else "EXPN"
            fitAltFunc = "EXPN" if nBtag < 2 else "POW"
        else:
            fitFunc = "EXPN" if nBtag < 2 else "EXPN"
            fitAltFunc = "EXPTAIL" if nBtag < 2 else "POW"
        fitFuncVjet = "ERFEXP" if nBtag < 2 else "EXP"
        fitAltFuncVjet = "POL" if nBtag < 2 else "POL"
        fitFuncVV = "EXPGAUS" if nBtag < 2 else "EXPGAUS"
        fitFuncTop = "GAUS3" if nBtag < 2 else "GAUS2"
    else:
        treeName = "XZh"
        signName = "XZh"
        colorVjet = sample["DYJetsToLL"]["linecolor"]
        triName = "HLT_Ele" if nElec > 0 else "HLT_Mu"
        leptCut = "isZtoEE" if nElec > 0 else "isZtoMM"
        topVeto = "X_dPhi>2.5"
        massVar = "X_mass"
        binFact = 2
        if nElec > 0:
            fitFunc = "EXPTAIL" if nBtag < 2 else "EXPTAIL"
            fitAltFunc = "POW" if nBtag < 2 else "POW"
        else:
            fitFunc = "EXPTAIL" if nBtag < 2 else "EXPTAIL"
            fitAltFunc = "POW" if nBtag < 2 else "POW"
        fitFuncVjet = "ERFEXP" if nBtag < 2 and nElec < 1 else "EXP"
        fitAltFuncVjet = "POL" if nBtag < 2 else "POL"
        fitFuncVV = "EXPGAUS2" if nBtag < 2 else "EXPGAUS2"
        fitFuncTop = "GAUS"

    btagCut = selection["2Btag"] if nBtag == 2 else selection["1Btag"]

    print "--- Channel", channel, "---"
    print "  number of electrons:", nElec, " muons:", nMuon, " b-tags:", nBtag
    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", "jet corrected pruned mass", HBINMIN, HBINMAX, "GeV")
    CSV1 = RooRealVar("fatjet1_CSVR1", "", -1.0e99, 1.0e4)
    CSV2 = RooRealVar("fatjet1_CSVR2", "", -1.0e99, 1.0e4)
    nB = RooRealVar("fatjet1_nBtag", "", 0.0, 4)
    CSVTop = RooRealVar("bjet1_CSVR", "", -1.0e99, 1.0e4)
    X_dPhi = RooRealVar("X_dPhi", "", 0.0, 3.15)
    isZtoEE = RooRealVar("isZtoEE", "", 0.0, 2)
    isZtoMM = RooRealVar("isZtoMM", "", 0.0, 2)
    isWtoEN = RooRealVar("isWtoEN", "", 0.0, 2)
    isWtoMN = RooRealVar("isWtoMN", "", 0.0, 2)
    weight = RooRealVar("eventWeightLumi", "", -1.0e9, 1.0)

    # 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, nB, CSVTop, X_dPhi)
    variables.add(RooArgSet(isZtoEE, isZtoMM, isWtoEN, isWtoMN, weight))
#.........這裏部分代碼省略.........

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
def setup_workspace():

  import ROOT
  from ROOT import RooWorkspace, gROOT, gStyle, RooAbsReal, RooMsgService, RooFit
  #from ROOT import RooFit, gROOT, gDirectory, gStyle, gPad, TTree, RooCmdArg,RooBinning
  #from ROOT import RooRealVar, RooMappedCategory, RooCategory, RooFormulaVar, RooAbsData
  #from ROOT import RooBMixDecay, RooMCStudy, RooAddModel, RooEffProd, RooMsgService
  #from ROOT import RooWorkspace, TCanvas, TFile, kFALSE, kTRUE, RooDataSet, TStopwatch
  #from ROOT import RooArgSet, RooArgList, RooRandom, RooMinuit, RooAbsReal, RooDataHist
  #from ROOT import TBrowser, TH2F, TF1, TH1F, RooGenericPdf, RooLinkedList

  gROOT.SetStyle("Plain")
  gStyle.SetPalette(1)
  gStyle.SetOptStat(0)
  gStyle.SetOptFit(0)
  gStyle.SetOptStat(1111)
  gStyle.SetOptFit(10111)
  gStyle.SetOptTitle(1)

  #RooAbsReal.defaultIntegratorConfig().Print()
  #RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-10)
  #RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-10)
  RooAbsReal.defaultIntegratorConfig().Print()

  print "Numeric integration set up" #TODO: is the integration acceptable?

  ##This controls the logging output from RooFit
  #RooMsgService.instance().addStream(RooFit.DEBUG,RooFit.Topic(RooFit.Fitting))
  RooMsgService.instance().deleteStream(1)
  #RooMsgService.instance().addStream(RooFit.INFO,RooFit.Topic(RooFit.Generation + RooFit.Minization + RooFit.Plotting + RooFit.Fitting + RooFit.Integration + RooFit.LinkStateMgmt + RooFit.Eval + RooFit.Caching + RooFit.Optimization + RooFit.ObjectHandling + RooFit.InputArguments + RooFit.Tracing + RooFit.Contents + RooFit.DataHandling + RooFit.NumericIntegration))
  RooMsgService.instance().addStream(RooFit.INFO,RooFit.Topic(RooFit.LinkStateMgmt + RooFit.Caching + RooFit.ObjectHandling + RooFit.InputArguments + RooFit.Tracing))
  RooMsgService.instance().Print()


  print "Message service set up"

  w = RooWorkspace("w",False)

  w.factory("RAND[0,1]")
  D0_M = w.factory("D0_M[1800,1930]") #TODO: define mass ranges slightly better
  Del_M = w.factory("Del_M[139,155]")
  D0_M.setUnit("MeV")
  Del_M.setUnit("MeV")

  for dst_side in ["", "dstsig", "dsthigh", "dstlow"]:
    for d_side in ["", "dsig", "dhigh", "dlow"]:
      name = dst_side+d_side
      if dst_side == "dsthigh":
        Del_M.setRange(name,148.,155.)
      elif dst_side == "dstsig":
        Del_M.setRange(name,143.,148.)
      elif dst_side == "dstlow":
        Del_M.setRange(name,139.,143.)
        
      if d_side == "dhigh":
        Del_M.setRange(name,1885.,1930.)
      elif d_side == "dsig":
        Del_M.setRange(name,1835.,1885.)
      elif d_side == "dlow":
        Del_M.setRange(name,1800.,1835.)




  w.defineSet("args","D0_M,Del_M")
  w.defineSet("argsPreCut","D0_M,Del_M,RAND")

  w.factory("RooGaussian::D0M_Sig_Gaus1(D0_M,D0M_Sig_Gaus_Mean[1865,1850,1880],D0M_Sig_Gaus_Sigma1[10,1,30])")
  w.factory("RooGaussian::D0M_Sig_Gaus2(D0_M,D0M_Sig_Gaus_Mean,D0M_Sig_Gaus_Sigma2[3,1,30])")
  w.factory("SUM::D0M_Sig_Gaus(D0M_Sig_Gaus1_Frac[0.8,0,1]*D0M_Sig_Gaus1,D0M_Sig_Gaus2)")
  
  
  w.factory("RooGaussian::D0M_MisId_Gaus1(D0_M,D0M_MisId_Gaus_Mean[1800,1740,1820],D0M_Sig_Gaus_Sigma1)")
  w.factory("RooGaussian::D0M_MisId_Gaus2(D0_M,D0M_MisId_Gaus_Mean,D0M_Sig_Gaus_Sigma2)")
  w.factory("SUM::D0M_MisId_Gaus(D0M_Sig_Gaus1_Frac*D0M_MisId_Gaus1,D0M_MisId_Gaus2)")
  
  
  w.factory("RooChebychev::D0M_Bkg_Poly(D0_M,{D0M_Bkg_Poly_a1[0,-1,1]})")
  


  w.factory("RooGaussian::DelM_Sig_Gaus1(Del_M,DelM_Sig_Gaus_Mean[145.5,143,148],DelM_Sig_Gaus_Sigma1[1,0,5] )")
  w.factory("RooGaussian::DelM_Sig_Gaus2(Del_M,DelM_Sig_Gaus_Mean,DelM_Sig_Gaus_Sigma2[.1,0,2] )")
  w.factory("SUM::DelM_Sig_Gaus(DelM_Sig_Gaus1_Frac[0.8,0,1]*DelM_Sig_Gaus1,DelM_Sig_Gaus2)")

  w.factory("RooDstD0BG::DelM_Bkg(Del_M,DelM_Bkg_m0[139.5,134,144],DelM_Bkg_c[80,0,1000],DelM_Bkg_a[-1,-100,10],DelM_Bkg_b[0.2,-0.2,10])")

  w.factory("PROD::Sig(DelM_Sig_Gaus,D0M_Sig_Gaus)")
  w.factory("PROD::Comb(DelM_Bkg,D0M_Bkg_Poly)")
  w.factory("PROD::MisId(DelM_Sig_Gaus,D0M_MisId_Gaus)")
  w.factory("PROD::Prompt(DelM_Bkg,D0M_Sig_Gaus)")
  

  w.factory("SUM::Final_PDF(N_Sig[10000,0,50000]*Sig,N_Prompt[5000,0,20000]*Prompt,N_Comb[10000,0,50000]*Comb,N_MisId[500,0,5000]*MisId)")


  return w

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
    def bayesian(self, ws, debug = 0):

        if self.sBayes not in self.items:
            if debug>0:
                print self.legend, 'no Bayesian Calculator section in the config file'
                print self.legend, 'cannot configure Bayesian calculator'
            return {'status':'fail'}

        legend = '['+self.sBayes+']:'

        _name = 'exostBayes'

        # getting items as a dictionary
        bayes_items = {}
        bayes_items.update(self.items[self.sBayes])
        
        print legend, 'Configuring Bayesian calculator... '

        # check if a CL is specified
        m_conf_name = self.check_value(self.sBayes, 'model_config')
        if m_conf_name == -1:
            print legend, 'Error: model config is not specified'
            print legend, 'Error:', self.sBayes, 'cannot be configured'
            return {'status':'fail'}
            

        # check if data is specified
        data_name = self.check_value(self.sBayes, 'data')

        data_valid = False
        if data_name != -1:
            data_valid = True
            if ws.data(data_name) != None:
                data = ws.data(data_name)
            else:
                data_valid = False
                print legend, 'Error: dataset', data_name, 'is not defined'

        if data_valid == False:
            print legend, 'Error:', self.sBayes, 'cannot be configured'
            return {'status':'fail'}


        # check if a CL is specified
        conf_level = self.check_value(self.sBayes, 'confidence_level')
        if conf_level == -1:
            print legend, 'Warning: desired confidence level is not specified, setting to 0.90'
            conf_level = 0.90


        # check if a posterior plot is requested
        post_plot = self.check_value(self.sBayes, 'posterior_plot')
        b_post_plot = False
        if post_plot == -1:
            print legend, 'will generate a posterior plot'
            b_post_plot = True
        elif post_plot == 'True':
            b_post_plot = True
        else:
            b_post_plot = False

        # plot format
        if post_plot:
            plot_format = self.check_value(self.sBayes, 'plot_format')
            if plot_format == -1:
                plot_format = 'png'


        # to suppress messgaes when pdf goes to zero
        RooMsgService.instance().setGlobalKillBelow(RooFit.FATAL)

        #mconf = ws.obj('exostModelConfig')
        mconf = ws.obj(m_conf_name)

        if mconf == None:
            if debug > 0:
                print self.legend, 'fail to get Model Config'
                print self.legend, 'unable to configure Bayesian calculator'
                print self.legend, 'calculator will not be imported into the workspace'
            return {'status':'fail'}
                    
        bCalc = RooStats.BayesianCalculator(data, mconf)
        bCalc.SetConfidenceLevel(float(conf_level))


        ############# FIXME: debug test
        #bInt = bCalc.GetInterval()
        #print legend, 'DEBUG3*******************************'
        # 
        #cl = bCalc.ConfidenceLevel()
        #print legend,  str(cl)+'% CL central interval: [', bInt.LowerLimit(), ' - ', bInt.UpperLimit(), ']'
        #print legend, 'or', str(cl+(1.0-cl)/2), '% CL limits'
        ##############################


        #size = 1.0 - float(conf_level)
        #bCalc.SetTestSize(size)
        
        # import the calculator into the Workspace
        getattr(ws, 'import')(bCalc, _name)
#.........這裏部分代碼省略.........

# 需要導入模塊: from ROOT import RooMsgService [as 別名]
# 或者: from ROOT.RooMsgService import instance [as 別名]
def setup_workspace(config):

  import ROOT
  from ROOT import RooWorkspace, gROOT, gStyle, RooAbsReal, RooMsgService, RooFit
  #from ROOT import RooFit, gROOT, gDirectory, gStyle, gPad, TTree, RooCmdArg,RooBinning
  #from ROOT import RooRealVar, RooMappedCategory, RooCategory, RooFormulaVar, RooAbsData
  #from ROOT import RooBMixDecay, RooMCStudy, RooAddModel, RooEffProd, RooMsgService
  #from ROOT import RooWorkspace, TCanvas, TFile, kFALSE, kTRUE, RooDataSet, TStopwatch
  #from ROOT import RooArgSet, RooArgList, RooRandom, RooMinuit, RooAbsReal, RooDataHist
  #from ROOT import TBrowser, TH2F, TF1, TH1F, RooGenericPdf, RooLinkedList
  from math import sqrt

  gROOT.SetStyle("Plain")
  gStyle.SetPalette(1)
  gStyle.SetOptStat(0)
  gStyle.SetOptFit(0)
  gStyle.SetOptStat(1111)
  gStyle.SetOptFit(10111)
  gStyle.SetOptTitle(1)

  #gROOT.ProcessLine(".L RooGaussianTrunk.cxx+")
  #gROOT.ProcessLine(".L RooCBShapeTrunk.cxx+")
  #gROOT.ProcessLine(".L RooChebychevTrunk.cxx+")
  #from ROOT import RooGaussianTrunk, RooChebychevTrunk, RooCBShapeTrunk

  #RooAbsReal.defaultIntegratorConfig().Print()
  RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-8)
  RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-8)
  #RooAbsReal.defaultIntegratorConfig().setEpsAbs(1e-6)
  #RooAbsReal.defaultIntegratorConfig().setEpsRel(1e-6)
  RooAbsReal.defaultIntegratorConfig().Print()

  print "Numeric integration set up" #TODO: is the integration acceptable?

  ##This controls the logging output from RooFit
  #RooMsgService.instance().addStream(RooFit.DEBUG,RooFit.Topic(RooFit.Fitting))
  RooMsgService.instance().deleteStream(1)
  #RooMsgService.instance().addStream(RooFit.INFO,RooFit.Topic(RooFit.Generation + RooFit.Minization + RooFit.Plotting + RooFit.Fitting + RooFit.Integration + RooFit.LinkStateMgmt + RooFit.Eval + RooFit.Caching + RooFit.Optimization + RooFit.ObjectHandling + RooFit.InputArguments + RooFit.Tracing + RooFit.Contents + RooFit.DataHandling + RooFit.NumericIntegration))
  RooMsgService.instance().addStream(RooFit.INFO,RooFit.Topic(RooFit.LinkStateMgmt + RooFit.Caching + RooFit.ObjectHandling + RooFit.InputArguments + RooFit.Tracing))
  RooMsgService.instance().Print()


  print "Message service set up"

  w = RooWorkspace("w",False)

  w.factory("RAND[0,1]")
  
  if "norm" not in config["mode"]:
    D0_Mass = w.factory("D0_Mass[1815,1915]") 
  else:
    D0_Mass = w.factory("D0_Mass[1800,1930]")
  D0_Mass.setUnit("MeV")
  D0_Mass.setBins(60)
  
  Del_Mass = w.factory("Del_Mass[139,155]")
  Del_Mass.setUnit("MeV")
  Del_Mass.setBins(60)

  if "norm" not in config["mode"]:
    Dataset = w.factory("DataSet[BDT1,BDT2,BDT3]")
  else:
    Dataset = w.factory("DataSet[Norm]")

  w.factory("classID[Sig=0,Bkg=1]")
  w.factory("BDT_ada[-1,1]")
  w.factory("x1_PIDe[-2,20]")
  w.factory("x2_ProbNNmu[0,1]")

  #D0_Mass.setRange("blinded",1700.,1900.)
  
  if "norm" not in config["mode"]:
    dataCats = ["", "BDT1", "BDT2", "BDT3"]
  else:
    dataCats = ["", "Norm"]

  for data in dataCats:
    for dst_side in ["", "delsig", "delhigh", "dellow"]:
      for d_side in ["", "dsig", "dhigh", "dlow", "dhigh1", "dlow1", "dhigh2", "dlow2"]:
        name = data+dst_side+d_side

        if data == "BDT1":
          Dataset.setRange(name,"BDT1")
        elif data == "BDT2":
          Dataset.setRange(name,"BDT2")
        elif data == "BDT3":
          Dataset.setRange(name,"BDT3")
        elif data == "Norm":
          Dataset.setRange(name,"Norm")

        if dst_side == "delhigh":
          Del_Mass.setRange(name,148.,155.)
        elif dst_side == "delsig":
          Del_Mass.setRange(name,143.,148.)
        elif dst_side == "dellow":
          Del_Mass.setRange(name,139.,143.)

        if d_side == "dhigh2":
          D0_Mass.setRange(name,1910.,1930.)
        elif d_side == "dhigh1":
#.........這裏部分代碼省略.........