test for systematics

wiso · Dec 23, 2018 · 2769071 · 2769071
1 parent 25ed4cf
commit 2769071
Show file tree

Hide file tree

Showing 2 changed files with 48 additions and 1 deletion.
diff --git a/countingworkspace/countingworkspace.py b/countingworkspace/countingworkspace.py
@@ -239,9 +239,13 @@ def create_workspace(categories, processes,
     create_expected_number_of_signal_events(ws, categories, processes, expression_nsignal_gen=expression_nsignal_gen_with_sys)
 
     all_constrains = ROOT.RooArgSet()
+    all_globals = ROOT.RooArgSet()
     for sysname in sysnames:
-        _ = ws.factory('RooGaussian:constrain_sys{sysname}(global_{sysname}[0, -5, 5], theta_{sysname}, 1)'.format(sysname=sysname))
+        global_obs = ws.factory('global_{sysname}[0, -5, 5]'.format(sysname=sysname))
+        global_obs.setConstant()
+        _ = ws.factory('RooGaussian:constrain_sys{sysname}(global_{sysname}, theta_{sysname}, 1)'.format(sysname=sysname))
         all_constrains.add(_)
+        all_globals.add(global_obs)
     ws.defineSet('constrains', all_constrains)
     ws.factory('PROD:prod_constrains(%s)' % ','.join([v.GetName() for v in utils.iter_collection(all_constrains)]))
 

diff --git a/tests/test_all.py b/tests/test_all.py
@@ -242,6 +242,49 @@ def test_fit_asimov():
         np.testing.assert_allclose(ntrue, poi_fitted.getVal(), rtol=0.002)
 
 
+def test_fit_asimov_syst():
+    systematics_nsignal_gen = np.ones(NPROCESS) * 0.05
+    systematics_nsignal_gen[0] *= 2
+
+    ws = create_workspace(NCATEGORIES, NPROCESS, NTRUE, EFFICIENCIES, EXPECTED_BKG_CAT,
+                          systematics_nsignal_gen=[{'name': 'lumi', 'values': systematics_nsignal_gen}])
+
+    obs = ws.set('all_obs')
+    pdf = ws.obj('model')
+    assert obs
+    assert pdf
+    data_asimov = countingworkspace.utils.generate_asimov(ws)
+    pois = ws.obj('ModelConfig').GetParametersOfInterest()
+    assert pois
+
+    # not start the fit from the true values
+    for poi in iter_collection(pois):
+        poi.setVal(poi.getVal() * 1.1)
+
+    fr = pdf.fitTo(data_asimov, ROOT.RooFit.Save(), ROOT.RooFit.Hesse(True))
+    assert(fr.status() == 0)
+    pois_fitted = fr.floatParsFinal()
+    all_errors = []
+    for ntrue, poi_fitted in zip(NTRUE, iter_collection(pois_fitted)):
+        np.testing.assert_allclose(ntrue, poi_fitted.getVal(), rtol=0.002)
+        all_errors.append(poi_fitted.getError())
+
+    ws.loadSnapshot('initial')
+    for theta in iter_collection(ws.allVars().selectByName('theta*')):
+        theta.setVal(0)
+        theta.setConstant()
+
+    fr = pdf.fitTo(data_asimov, ROOT.RooFit.Save(), ROOT.RooFit.Hesse(True))
+    assert(fr.status() == 0)
+    pois_fitted = fr.floatParsFinal()
+    all_errors_stat = []
+    for ntrue, poi_fitted in zip(NTRUE, iter_collection(pois_fitted)):
+        np.testing.assert_allclose(ntrue, poi_fitted.getVal(), rtol=0.002)
+        all_errors_stat.append(poi_fitted.getError())    
+
+    sys_only_errors = (np.sqrt(np.array(all_errors)**2 - np.array(all_errors_stat)**2) / NTRUE)
+    np.testing.assert_allclose(sys_only_errors, systematics_nsignal_gen, rtol=2E-2)
+
 def test_create_workspace_raise():
     with pytest.raises(ValueError):
         create_workspace(