fix: fix HEM

Author: andrzejnovak

README.md

@@ -58,7 +58,8 @@ bash build.sh
 # text2workspace.py -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel  --PO verbose --PO 'map=.*/*hcc*:r[1,-500,500]' model_combined.txt
 # text2workspace.py -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel  --PO verbose --PO 'map=.*/zcc:r[1,-5,5]' model_combined.txt
 
-combine -M FitDiagnostics --expectSignal 1 -d model_combined.root --cminDefaultMinimizerStrategy 0 --robustFit=1 -t -1 --toysFrequentist 
+combine -M FitDiagnostics -d model_combined.root --cminDefaultMinimizerStrategy 0 --robustFit=1 
+--setParameters z=1,r=1 -n "" -t -1 --toysFrequentist 
 combine -M FitDiagnostics --expectSignal 1 -d model_combined.root --cminDefaultMinimizerStrategy 0 --robustFit=1 --saveShapes --saveWithUncertainties -t -1 --toysFrequentist --setParameters z=1 -n ""
 combine -M Significance model_combined.root --expectSignal 1 --redefineSignalPOIs z -t -1 --toysFrequentist
 combineTool.py -M AsymptoticLimits -m 125 -d model_combined.root --expectSignal 1 --setParameters z=1 --redefineSignalPOIs r -t -1 --toysFrequentist 
@@ -83,6 +84,17 @@ combineTool.py -M Impacts -d model_combined.root -m 125 --redefineSignalPOIs z -
 plotImpacts.py -i impactsZ.json -o plots/impacts_out_Z
 ```
 
+Fitting Z unblinding
+```
+# Baseline
+combineTool.py -M Impacts -d model_combined.root -m 125 --doInitialFit --robustFit 1 --setParameterRanges r=-1,5 --cminDefaultMinimizerStrategy 0 --X-rtd FITTER_DYN_STEP --expectSignal 1 --redefineSignalPOIs z
+# Condor
+combineTool.py -M Impacts -d model_combined.root -m 125 --doFits --robustFit 1 --allPars --setParameterRanges r=-1,5  --redefineSignalPOIs z --cminDefaultMinimizerStrategy 0 --X-rtd MINIMIZER_analytic --job-mode condor --sub-opts='+JobFlavour = "workday"' --task-name ggHccZ --exclude 'rgx{qcdparams*}'
+# Collect
+combineTool.py -M Impacts -d model_combined.root -m 125 --redefineSignalPOIs z -o impactsZunbl.json
+plotImpacts.py -i impactsZunbl.json -o plots/impacts_out_Zunbl --blind
+```
+
 Fitting H
 ```
 # Baseline

new_Hxx.py

@@ -96,26 +96,35 @@ def shape_to_num(f, region, sName, ptbin, syst, mask, muon=False, bound=0.5):
     _nom_rate = np.sum(_nom[0] * mask)
     if _nom_rate < .1:
         return 1.0
+
+    _one_side = get_templ(f, region, sName, ptbin, syst=syst, muon=muon)
     _up = get_templ(f, region, sName, ptbin, syst=syst + "Up", muon=muon)
     _down = get_templ(f, region, sName, ptbin, syst=syst + "Down", muon=muon)
-    if _up is None and _down is None:
+    if _up is None and _down is None and _one_side is None:
         return None
     else:
-        if _down is not None:
-            _down_rate = np.sum(_down[0] * mask)
-        else:
-            _down_rate = _nom_rate
-        if _up is not None:
+        if _one_side is not None:
+            _up_rate = np.sum(_one_side[0] * mask)            
+            _diff = np.abs(_up_rate - _nom_rate)
+            magnitude = _diff / _nom_rate
+        elif _down is not None and _up is not None:
             _up_rate = np.sum(_up[0] * mask)
+            _down_rate = np.sum(_down[0] * mask)
+            _diff = np.abs(_up_rate - _nom_rate) + np.abs(_down_rate - _nom_rate)
+            magnitude = _diff / (2. * _nom_rate)
         else:
-            _up_rate = _nom_rate
-    _diff = np.abs(_up_rate - _nom_rate) + np.abs(_down_rate - _nom_rate)
-    magnitude = _diff / (2. * _nom_rate)
+            raise NotImplementedError
     if bound is not None:
         magnitude = min(magnitude, bound)
     return 1.0 + magnitude
 
 def get_templ(f, region, sample, ptbin, syst=None, muon=False):
+    if "16" in f.name:
+        year = 2016
+    elif "17" in f.name:
+        year = 2017
+    else:
+        year = 2018
     hist_name = '{}_{}'.format(sample, region)
     if syst is not None:
         hist_name += "_" + syst
@@ -126,8 +135,16 @@ def get_templ(f, region, sample, ptbin, syst=None, muon=False):
     try:
         f[hist_name]
     except:
-        print("{}Sample {}, {}, {}, {} not found.".format('(Muon) ' if muon else "",
-            sample, region, ptbin if not muon else "-", syst))
+        if syst is not None:
+            if "HEM" in syst and year in [2016, 2017]:  # always empty
+                pass
+            elif "HEM" in syst and ("Up" in syst or "Down" in syst) and year == 2018:  # always empty
+                pass
+            elif "L1Prefiring" in syst and year == 2018:  # always empty
+                pass
+            else:
+                print("{}Sample {}, {}, {}, {} not found.".format('(Muon) ' if muon else "",
+                    sample, region, ptbin if not muon else "-", syst))
         return None
     h_vals = f[hist_name].values
     h_edges = f[hist_name].edges
@@ -172,6 +189,9 @@ def dummy_rhalphabet(pseudo,
                      year=2017,
                      opts=None):
 
+    print(year)
+    assert year in ['2016', '2017', '2018']
+
     # Default lumi (needs at least one systematics for prefit)
     sys_lumi = rl.NuisanceParameter('CMS_lumi_13TeV_{}'.format(year), 'lnN')
     sys_lumi_correlated = rl.NuisanceParameter('CMS_lumi_13TeV_correlated', 'lnN')
@@ -389,7 +409,6 @@ def dummy_rhalphabet(pseudo,
 
     for ptbin in range(npt):
         for region in ['pass', 'fail']:
-            #ch = rl.Channel("ptbin%d%s" % (ptbin, region))
             ch = rl.Channel("ptbin{}{}{}".format(ptbin, region, year))
             model.addChannel(ch)
             if justZ:
@@ -953,13 +972,6 @@ def str2bool(v):
 
     parser.add_argument("--year", type=int, default=2017, help="Year")
 
-    parser.add_argument("--matched",
-                        type=str2bool,
-                        default='False',
-                        choices={True, False},
-                        help=("Use matched/unmatched templates"
-                              "(w/o there is some W/Z/H contamination from QCD)"))
-
     parser.add_argument("--mcstat",
                         type=str2bool,
                         default='True',