Merge pull request #1927 from NNPDF/W_hel_rap_200804174

W-> l v pseudo rapidity (200804174)

validphys2/src/validphys/datafiles/new_commondata/CMS_WPWM_13TEV_ETA/data_WM.yaml

@@ -0,0 +1,19 @@
+data_central:
+- 1637.0
+- 1647.0
+- 1655.0
+- 1646.0
+- 1633.0
+- 1635.0
+- 1630.0
+- 1621.0
+- 1628.0
+- 1605.0
+- 1602.0
+- 1574.0
+- 1561.0
+- 1533.0
+- 1479.0
+- 1423.0
+- 1384.0
+- 1278.0

validphys2/src/validphys/datafiles/new_commondata/CMS_WPWM_13TEV_ETA/data_WP.yaml

@@ -0,0 +1,19 @@
+data_central:
+- 1956.0
+- 1968.0
+- 1973.0
+- 1976.0
+- 1970.0
+- 1978.0
+- 1986.0
+- 1988.0
+- 2013.0
+- 2012.0
+- 2032.0
+- 2025.0
+- 2032.0
+- 2038.0
+- 2025.0
+- 2026.0
+- 2051.0
+- 1976.0

validphys2/src/validphys/datafiles/new_commondata/CMS_WPWM_13TEV_ETA/filter.py

@@ -0,0 +1,95 @@
+import yaml
+
+from filter_utils import get_kinematics, get_data_values, get_systematics
+
+
+def filter_CMS_W_13TEV_data_kinetic(figure):
+    """
+    writes data central values and kinematics
+    to respective .yaml file
+    """
+    with open("metadata.yaml", "r") as file:
+        metadata = yaml.safe_load(file)
+
+    version = metadata["hepdata"]["version"]
+    tables = metadata["implemented_observables"][0]["tables"]
+
+    kin = get_kinematics(version, figure)
+    central_values = get_data_values(version, figure)
+
+    data_central_yaml = {"data_central": central_values}
+    kinematics_yaml = {"bins": kin}
+
+    # write central values and kinematics to yaml file
+    if figure == "17a":
+        with open("data_WP.yaml", "w") as file:
+            yaml.dump(data_central_yaml, file, sort_keys=False)
+
+        with open("kinematics_WP.yaml", "w") as file:
+            yaml.dump(kinematics_yaml, file, sort_keys=False)
+
+    elif figure == "17b":
+        with open("data_WM.yaml", "w") as file:
+            yaml.dump(data_central_yaml, file, sort_keys=False)
+
+        with open("kinematics_WM.yaml", "w") as file:
+            yaml.dump(kinematics_yaml, file, sort_keys=False)
+
+
+def filter_CMS_W_13TEV_uncertainties(observable, figure):
+    """
+    writes uncertainties to respective .yaml file
+
+    Parameters
+    ----------
+    observable : str
+        eg. "W+" or "W-"
+    """
+
+    with open("metadata.yaml", "r") as file:
+        metadata = yaml.safe_load(file)
+    version = metadata["hepdata"]["version"]
+
+    systematics = get_systematics(observable, version, figure)
+
+    # error definition
+    error_definitions = {}
+    errors = []
+
+    for sys in systematics:
+
+        error_definitions[sys[0]['name']] = {
+            "description": f"{sys[0]['name']}",
+            "treatment": "ADD",
+            "type": "CORR",
+        }
+
+    #
+    for i in range(metadata['implemented_observables'][0]['ndata']):
+        error_value = {}
+
+        for sys in systematics:
+            error_value[sys[0]['name']] = float(sys[0]['values'][i])
+
+        errors.append(error_value)
+
+    uncertainties_yaml = {"definitions": error_definitions, "bins": errors}
+
+    # write uncertainties
+    if observable == "W+":
+        with open(f"uncertainties_WP.yaml", 'w') as file:
+            yaml.dump(uncertainties_yaml, file, sort_keys=False)
+    elif observable == "W-":
+        with open(f"uncertainties_WM.yaml", 'w') as file:
+            yaml.dump(uncertainties_yaml, file, sort_keys=False)
+
+
+if __name__ == "__main__":
+    # WP data
+    filter_CMS_W_13TEV_data_kinetic(figure="17a")
+    filter_CMS_W_13TEV_uncertainties(observable="W+", figure="17a")
+
+    # WM data
+    filter_CMS_W_13TEV_data_kinetic(figure="17b")
+    filter_CMS_W_13TEV_uncertainties(observable="W-", figure="17b")
+

validphys2/src/validphys/datafiles/new_commondata/CMS_WPWM_13TEV_ETA/filter_utils.py

@@ -0,0 +1,170 @@
+import yaml
+import uproot
+import numpy as np
+
+
+def get_kinematics(version, figure):
+    """
+    returns the relevant kinematics values.
+
+    Parameters
+    ----------
+    version : int
+            integer read from metadata.yaml that
+            indicated the version of the hepdata
+            tables
+
+    figure : str
+        eg. 17a or 17b
+
+    Returns
+    -------
+    list
+        list containing the kinematic values for all
+        hepdata tables
+    """
+    kin = []
+
+    hepdata_table = f"rawdata/HEPData-ins1810913-v{version}-Figure_{figure}.yaml"
+
+    with open(hepdata_table, 'r') as file:
+        input = yaml.safe_load(file)
+
+    for eta in input["independent_variables"][0]['values']:
+        kin_value = {
+            'eta': {'min': eta['low'], 'mid': 0.5 * (eta['low'] + eta['high']), 'max': eta['high']},
+            'm_W2': {'min': None, 'mid': 6460.5, 'max': None},
+            'sqrts': {'min': None, 'mid': 13000.0, 'max': None},
+        }
+
+        kin.append(kin_value)
+
+    return kin
+
+
+def get_data_values(version, figure):
+    """
+    returns the central data.
+
+    Parameters
+    ----------
+    version : int
+            integer read from metadata.yaml that
+            indicated the version of the hepdata
+            tables
+
+    figure : str
+        eg. 17a or 17b
+
+    Returns
+    -------
+    list
+        list containing the central values for all
+        hepdata tables
+
+    """
+
+    data_central = []
+
+    hepdata_table = f"rawdata/HEPData-ins1810913-v{version}-Figure_{figure}.yaml"
+
+    with open(hepdata_table, 'r') as file:
+        input = yaml.safe_load(file)
+
+    values = input['dependent_variables'][0]['values']
+
+    for value in values:
+        data_central.append(value['value'])
+
+    return data_central
+
+
+def decompose_covmat(covmat):
+    """Given a covmat it return an array sys with shape (ndat,ndat)
+    giving ndat correlated systematics for each of the ndat point.
+    The original covmat is obtained by doing sys@sys.T"""
+
+    lamb, mat = np.linalg.eig(covmat)
+    sys = np.multiply(np.sqrt(lamb), mat)
+    return sys
+
+
+def get_systematics(observable, version, figure):
+    """
+    Following the CMS advice we take the covariance matrix from
+    https://cms-results.web.cern.ch/cms-results/public-results/publications/SMP-18-012/index.html
+
+    The root file sumpois 2dxsec contains the needed covariance matrix.
+
+    Parameters
+    ----------
+    observable : str
+        The observable for which the covariance matrix is needed
+        can be either "W+" or "W-"
+
+    """
+
+    # Open the ROOT file
+    file = uproot.open("rawdata/covariance_coefficients_sumpois_2dxsec.root")
+
+    # access the TH2D histogram
+    histogram = file["sub_cov_matrix"]
+
+    # Get the labels along x-axis
+    x_labels = histogram.axis(0).labels()
+
+    # Select rows whose label starts with "POI" and contains the observable
+    if observable == "W+":
+        poi_indices = [
+            i for i, label in enumerate(x_labels) if "POI, $W^{+}$, $|\\eta^{l}|" in label
+        ]
+
+    elif observable == "W-":
+        poi_indices = [
+            i for i, label in enumerate(x_labels) if "POI, $W^{-}$, $|\\eta^{l}|" in label
+        ]
+
+    # Extract the submatrix
+    submatrix = histogram.values()[poi_indices][:, poi_indices]
+
+    # Convert submatrix to numpy array
+    submatrix_array = np.array(submatrix)
+
+    # Get Luminosity covariance matrix
+    if observable=="W+":
+        with open("rawdata/HEPData-ins1810913-v1-Impacts_Figure_A23a.yaml", "r") as file:
+            impacts = yaml.safe_load(file)
+    elif observable=="W-":
+        with open("rawdata/HEPData-ins1810913-v1-Impacts_Figure_A23b.yaml", "r") as file:
+            impacts = yaml.safe_load(file)
+
+
+
+    lumi_unc = np.array([val['value'] for val in impacts['dependent_variables'][2]['values']])
+
+    hepdata_table = f"rawdata/HEPData-ins1810913-v{version}-Figure_{figure}.yaml"
+
+    with open(hepdata_table, 'r') as file:
+        input = yaml.safe_load(file)
+
+    values = input['dependent_variables'][0]['values']
+    values = np.array([val['value'] for val in values])
+
+    lumi_unc *= values / 100
+    lumi_covmat = lumi_unc[:, np.newaxis] @ lumi_unc[:, np.newaxis].T
+
+    artificial_uncertainties = np.real(decompose_covmat(lumi_covmat+submatrix_array))
+
+    uncertainties = []
+
+    for i, unc in enumerate(artificial_uncertainties.T):
+
+        name = f"artificial_uncertainty_{i}"
+        values = [unc[i] for i in range(len(unc))]
+        uncertainties.append([{"name": name, "values": values}])
+
+    return uncertainties
+
+
+if __name__ == "__main__":
+    get_systematics(observable="W+", version=1, figure='17a')