Skip to content

Commit

Permalink
fix?
Browse files Browse the repository at this point in the history
  • Loading branch information
VanyaBelyaev committed Dec 26, 2024
1 parent e28002a commit a350207
Show file tree
Hide file tree
Showing 3 changed files with 24 additions and 79 deletions.
58 changes: 7 additions & 51 deletions ostap/fitting/roostats.py
Original file line number Diff line number Diff line change
Expand Up @@ -64,42 +64,7 @@
if (3,0) <= sys.version_info : from itertools import zip_longest
else : from itertools import izip_longest as zip_longest
# =============================================================================
## Helper function to add constraints to RooSimultaneous
def simfit_with_constraints ( pdf , *constraints ) :
""" Helper function to add constraints to RooSimultaneous
"""

assert constraints and all ( v and isinstance ( v , ROOT.RooAbsPdf ) for v in constraints ), \
"Invalid constraintss %s" % str ( constraints )

assert pdf and isinstance ( pdf , ROOT.RooSimultaneous ), \
"PDF is not RooSimultaneous!"

name = '%s_with_constraints' % pdf.name
title = '%s with constraints' % pdf.title
cat = pdf.indexCat()

newpdf = ROOT.RooSimultaneous ( name , title , cat )
keep = []

labels = cat.labels()
for label, index in cat.items() :
print ( 'LABEL/INDEX', label, index )
cpdf = pdf.getPdf ( label )
lst = ROOT.RooArgList()
lst.add ( cpdf )
for c in constraints : lst.add ( c )
nname = '%s_with_constraints' % cpdf.name
ntitle = '%s with constraints' % cpdf.title
nc = ROOT.RooProdPdf ( nname , ntitle , lst )
newpdf.addPdf ( nc , str ( label ) )
keep.append ( ( nc , lst ) )

keep = tuple ( keep )
return newpdf, keep

# =============================================================================
# @class ModelConfig
## @class ModelConfig
# Helper class to create `RooStats::ModelConfig`
# @see RooStats::ModelConfig
class ModelConfig(object):
Expand Down Expand Up @@ -145,8 +110,7 @@ def __init__ ( self ,
mc = ROOT.RooStats.ModelConfig ( mcname , mctitle , ws )
self.__mc = mc

print ( 'MODEL-POI/1', poi )


self.__pdf = pdf
self.__poi = poi
self.__input_observables = observables
Expand Down Expand Up @@ -179,18 +143,12 @@ def __init__ ( self ,
if not cp in lgobs : lgobs.add ( cp )
cnts.add ( c )

if isinstance ( self.__raw_pdf , ROOT.RooSimultaneous ) :
logger.warning ("specific treatment of Roosimultaneous!")
self.__final_pdf, aux = simfit_with_constraints ( self.__raw_pdf , *constraints )
self.__cs = constraints, aux

else :
clst = ROOT.RooArgList ()
clst.add ( self.__raw_pdf )
for c in cnts : clst.add ( c )
self.__final_pdf = ROOT.RooProdPdf ( '%s_constrained' % pdf.name ,
clst = ROOT.RooArgList ()
clst.add ( self.__raw_pdf )
for c in cnts : clst.add ( c )
self.__final_pdf = ROOT.RooProdPdf ( '%s_constrained' % pdf.name ,
"PDF with %d constraints" % len ( cnts ) , clst )
self.__cs = constraints, cnts, clst
self.__cs = constraints, cnts, clst

## attention, redefine/update
constraints = cnts
Expand Down Expand Up @@ -224,8 +182,6 @@ def __init__ ( self ,
mc.SetParametersOfInterest ( pois )
self.__ps = params , pars, pois

print ( 'MODEL-POI/2', pois )

## (6) Nuisance parameters <---------- ATTENTION!! move them to global observables!
pars = [ v for v in self.pdf_params ( dataset ) if not v in observables and not v in pois ] ## and not v in lgobs ]
nuis = ROOT.RooArgSet ()
Expand Down
3 changes: 0 additions & 3 deletions ostap/fitting/simfit.py
Original file line number Diff line number Diff line change
Expand Up @@ -1470,9 +1470,6 @@ def fit_options ( self ) :
def fit_options ( self , value ) :
self.pdf.fit_options = value


#

# =============================================================================
if '__main__' == __name__ :

Expand Down
42 changes: 17 additions & 25 deletions ostap/fitting/tests/test_fitting_roostats2.py
Original file line number Diff line number Diff line change
Expand Up @@ -85,29 +85,23 @@ def test_point_limit_ac1() :

## create ModelConfig for 'S+B' model
model_sb = ModelConfig ( pdf = the_model ,
poi = rS.S , ## parameter of interest
poi = the_model.S , ## parameter of interest
dataset = data ,
name = 'S+B' ,
snapshot = rS ) ## ATTENTION!

print ( 'POI-1', the_model.S , model_sb.poi )

with FIXVAR ( the_model.S ) :
the_model.S = 0
rB , _ = the_model.fitTo ( data )
print ( 'POI-2', the_model.S , rB.S )

# # create ModelConfig for 'B-only' model
model_b = ModelConfig ( pdf = the_model ,
poi = rB.S , ## parameter of interest
dataset = data ,
workspace = model_sb.workspace ,
name = 'B-only' ,
snapshot = rB )
print ( 'POI-3', the_model.S , model_b.poi )

print ( 'POI-4', the_model.S , model_sb.poi , model_b.poi )


## create ModelConfig for 'B-only' model
model_b = ModelConfig ( pdf = the_model ,
poi = the_model.S , ## parameter of interest
dataset = data ,
workspace = model_sb.workspace ,
name = 'B-only' ,
snapshot = rB )

logger.info ( 'Model config %s\n%s' % ( model_sb.name , model_sb.table ( prefix = '# ' ) ) )
logger.info ( 'Model config %s\n%s' % ( model_b.name , model_b .table ( prefix = '# ' ) ) )

Expand Down Expand Up @@ -175,13 +169,13 @@ def test_point_limit_ac2() :
the_model.S = 0
rB , _ = the_model.fitTo ( data )

## create ModelConfig for 'B-only' model
model_b = ModelConfig ( pdf = the_model ,
poi = the_model.S , ## parameter of interest
dataset = data ,
workspace = model_sb.workspace ,
name = 'B-only' ,
snapshot = rB )
## create ModelConfig for 'B-only' model
model_b = ModelConfig ( pdf = the_model ,
poi = the_model.S , ## parameter of interest
dataset = data ,
workspace = model_sb.workspace ,
name = 'B-only' ,
snapshot = rB )

logger.info ( 'Model config %s\n%s' % ( model_sb.name , model_sb.table ( prefix = '# ' ) ) )
logger.info ( 'Model config %s\n%s' % ( model_b.name , model_b .table ( prefix = '# ' ) ) )
Expand Down Expand Up @@ -625,15 +619,13 @@ def test_point_limit3 () :
test_point_limit_ac1 ()
test_point_limit_ac2 ()

"""
test_point_limit_fc ()
test_point_limit_hc ()

test_point_limit_pl ()

test_point_limit2 ()
test_point_limit3 ()
"""

import ostap.logger.table as T
title = 'Summary of 90%CL Upper Limits'
Expand Down

0 comments on commit a350207

Please sign in to comment.