diff --git a/validphys2/src/validphys/config.py b/validphys2/src/validphys/config.py index 67bc181326..ab0d06865a 100644 --- a/validphys2/src/validphys/config.py +++ b/validphys2/src/validphys/config.py @@ -1157,28 +1157,6 @@ def produce_loaded_user_covmat_path(self, user_covmat_path: str = ""): fileloc = l.check_vp_output_file(user_covmat_path) return fileloc - @configparser.explicit_node - def produce_covmat_custom(self, use_ht_uncertainties: bool = False): - if use_ht_uncertainties: - from validphys.theorycovariance.construction import thcov_ht - - return thcov_ht - else: - from validphys.theorycovariance.construction import covs_pt_prescrip - - return covs_pt_prescrip - - @configparser.explicit_node - def produce_combine_custom(self, use_ht_uncertainties: bool = False): - if use_ht_uncertainties: - from validphys.theorycovariance.construction import combine_by_type_ht - - return combine_by_type_ht - else: - from validphys.theorycovariance.construction import combine_by_type - - return combine_by_type - @configparser.explicit_node def produce_nnfit_theory_covmat( self, point_prescriptions: list = None, user_covmat_path: str = None @@ -1205,31 +1183,8 @@ def produce_nnfit_theory_covmat( from validphys.theorycovariance.construction import user_covmat_fitting f = user_covmat_fitting - elif use_ht_uncertainties: - # NOTE: this covmat is the same as for scale variations, which will result in a clash of - # table names if we wish to use them simultaneously - if use_user_uncertainties: - from validphys.theorycovariance.construction import total_theory_covmat_fitting - - f = total_theory_covmat_fitting - else: - from validphys.theorycovariance.construction import theory_covmat_custom_fitting - - f = theory_covmat_custom_fitting - @functools.wraps(f) - def res(*args, **kwargs): - return f(*args, **kwargs) - - # Set this to get the same filename regardless of the action. - res.__name__ = "theory_covmat" - return res - - @configparser.explicit_node - def produce_combine_by_type_custom(self, use_ht_uncertainties: bool = False): - if use_ht_uncertainties: - return validphys.theorycovariance.construction.combine_by_type_ht - return validphys.theorycovariance.construction.combine_by_type + return f def produce_fitthcovmat( self, use_thcovmat_if_present: bool = False, fit: (str, type(None)) = None @@ -1707,6 +1662,8 @@ def produce_theoryids(self, t0id, point_prescription): prescription. The options for the latter are defined in pointprescriptions.yaml. This hard codes the theories needed for each prescription to avoid user error.""" th = t0id.id + if point_prescription == 'power corrections': + return NSList([t0id], nskey="theoryid") lsv = yaml_safe.load(read_text(validphys.scalevariations, "scalevariationtheoryids.yaml")) @@ -1809,6 +1766,30 @@ def produce_total_phi_data(self, fitthcovmat): return validphys.results.total_phi_data_from_experiments return validphys.results.dataset_inputs_phi_data + # @configparser.explicit_node + def produce_power_corr_dict(self, pc_parameters=None): + """The parameters for the power corrections are given as a list. + This function converts this list into a dictionary with the keys + being the names of the types of power corrections (e.g. `H2p`, `H2d`,...). + """ + if pc_parameters is None: + return None + + pc_parameters_by_type = {} + # Loop over the parameterization for the power corrections in the runcard + for par in pc_parameters: + # Check that the length of shifts matches the length of nodes. + if len(par['yshift']) != len(par['nodes']): + raise ValueError( + f"The length of nodes does not match that of the list in {par['ht']}." + f"Check the runcard. Got {len(par['yshift'])} != {len(par['nodes'])}" + ) + + # Store parameters for each power correction + pc_parameters_by_type[par['ht']] = {'yshift': par['yshift'], 'nodes': par['nodes']} + + return pc_parameters_by_type + class Config(report.Config, CoreConfig): """The effective configuration parser class.""" diff --git a/validphys2/src/validphys/theorycovariance/construction.py b/validphys2/src/validphys/theorycovariance/construction.py index 6bd4db09f2..1b44c3ee1d 100644 --- a/validphys2/src/validphys/theorycovariance/construction.py +++ b/validphys2/src/validphys/theorycovariance/construction.py @@ -13,11 +13,9 @@ from reportengine.table import table pass -from validphys.results import results, results_central from validphys.core import PDF -from validphys.process_options import _Process from validphys.results import results, results_central -import validphys.theorycovariance.higher_twist_functions as ht_func +from validphys.theorycovariance.higher_twist_functions import compute_deltas_pc from validphys.theorycovariance.theorycovarianceutils import ( check_correct_theory_combination, check_fit_dataset_order_matches_grouped, @@ -51,10 +49,10 @@ def theory_covmat_dataset(results, results_central_bytheoryids, point_prescripti return thcovmat -ProcessInfo = namedtuple("ProcessInfo", ("preds", "namelist", "sizes", "data")) +ProcessInfo = namedtuple("ProcessInfo", ("preds", "namelist", "sizes", "data", "data_spec")) -def combine_by_type(each_dataset_results_central_bytheory): +def combine_by_type(each_dataset_results_central_bytheory, groups_data_by_process): """Groups the datasets by process and returns an instance of the ProcessInfo class Parameters @@ -72,6 +70,7 @@ def combine_by_type(each_dataset_results_central_bytheory): dataset_size = defaultdict(list) theories_by_process = defaultdict(list) ordered_names = defaultdict(list) + data_spec = defaultdict(list) for dataset in each_dataset_results_central_bytheory: name = dataset[0][0].name theory_centrals = [x[1].central_value for x in dataset] @@ -81,217 +80,23 @@ def combine_by_type(each_dataset_results_central_bytheory): theories_by_process[proc_type].append(theory_centrals) for key, item in theories_by_process.items(): theories_by_process[key] = np.concatenate(item, axis=1) - process_info = ProcessInfo( - preds=theories_by_process, namelist=ordered_names, sizes=dataset_size, data=None - ) - return process_info + # Store DataGroupSpecs instances + for group_proc in groups_data_by_process: + for exp_set in group_proc.datasets: + data_spec[exp_set.name] = exp_set -def combine_by_type_ht(each_dataset_results, groups_dataset_inputs_loaded_cd_with_cuts_byprocess): - """same as combine_by_type but now for a single theory and including commondata info""" - dataset_size = defaultdict(list) - theories_by_process = defaultdict(list) - cd_by_process = defaultdict(list) - ordered_names = defaultdict(list) - for dataset, cd in zip( - each_dataset_results, groups_dataset_inputs_loaded_cd_with_cuts_byprocess - ): - name = cd.setname - if name != dataset[0].name: - raise ValueError("The underlying datasets do not match!") - theory_centrals = [x.central_value for x in dataset] - dataset_size[name] = len(theory_centrals[0]) - proc_type = process_lookup(name) - ordered_names[proc_type].append(name) - cd_by_process[proc_type].append(cd.kinematics.values) - theories_by_process[proc_type].append(theory_centrals) - - for key in theories_by_process.keys(): - theories_by_process[key] = np.concatenate(theories_by_process[key], axis=1) - cd_by_process[key] = np.concatenate(cd_by_process[key], axis=0) process_info = ProcessInfo( - preds=theories_by_process, namelist=ordered_names, sizes=dataset_size, data=cd_by_process + preds=theories_by_process, + namelist=ordered_names, + sizes=dataset_size, + data=None, + data_spec=data_spec, ) return process_info -def thcov_shifts_ht(ht_parameters, ht_included_proc, ht_excluded_exp, groups_data_by_process, pdf): - """ - Same as `thcov_HT` but implementing theory covariance method for each node of the spline. - Note that 'groups_data_by_process' contains the same info as 'combine_by_type_ht'. At some - point we should use only one of them. - """ - groups_data = groups_data_by_process - start_proc_by_exp = defaultdict(list) - ndata_by_exp = defaultdict(list) - deltas = defaultdict(list) - running_index_tot = 0 - - HT = {} - HT_func = {} - - for par in ht_parameters: - if len(par['list']) != len(par['nodes']): - raise ValueError( - f"The length of nodes does not match that of the list in {par['ht']}. Check the runcard.\n \ - {len(par['list'])} vs. {len(par['nodes'])}" - ) - - HT[par['ht']] = {"list": par['list'], "nodes": par['nodes']} - - HT_func[par['ht']] = None - - # Initialise shifts according to 5pt - for idx_node, _ in enumerate(par['nodes']): - deltas[par['ht'] + f"({idx_node})"] = np.array([]) - - for idx_proc, group_proc in enumerate(groups_data): - for idx_exp, exp_set in enumerate(group_proc.datasets): - # For covmat construction - exp_ndata = exp_set.load_commondata().ndata - start_proc_by_exp[exp_set.name] = running_index_tot - running_index_tot += exp_ndata - ndata_by_exp[exp_set.name] = exp_ndata - - for ht in HT_func.keys(): - HT_func[ht] = ht_func.null_func(exp_ndata) - - if group_proc.name in ht_included_proc and exp_set.name not in ht_excluded_exp: - cd_table = exp_set.load_commondata().commondata_table - process_type = cd_table['process'].iloc[0] - - if isinstance(process_type, _Process): - process_type = process_type.name - - x = cd_table['kin1'].to_numpy() - q2 = cd_table['kin2'].to_numpy() - y = cd_table['kin3'].to_numpy() - - if x.size != exp_ndata: - raise ValueError("Problem with the number of data.") - - # NMC_NC_NOTFIXED_DW_EM-F2 - if process_type == "DIS_F2R": - HT_func['H2p'], HT_func['H2d'] = ht_func.DIS_F2R_ht( - exp_set, pdf, HT["H2p"], HT["H2d"], x, q2 - ) - - elif process_type == "DIS_F2P": - HT_func['H2p'] = ht_func.DIS_F2_ht(HT["H2p"], x, q2) - - elif process_type == "DIS_F2D": - HT_func['H2d'] = ht_func.DIS_F2_ht(HT["H2d"], x, q2) - - # EMC - elif process_type == "DIS_F2C": - HT_func['H2p'], HT_func['H2d'] = ht_func.DIS_F2C_ht(HT['H2p'], HT['H2d'], x, q2) - - # HERA NC - elif process_type in ["DIS_NCE", "DIS_NCP", "DIS_NCP_CH", "DIS_NCE_BT"]: - HT_func['H2p'], HT_func["HLp"] = ht_func.DIS_NC_ht( - HT['H2p'], HT['HLp'], x, q2, y - ) - - # CHORUS - elif process_type in ["DIS_SNU", "DIS_SNB"]: - # Lead target:: - A = 208.0 - Z = 82 - if process_type == "DIS_SNU": - l = 0 - elif process_type == "DIS_SNB": - l = 1 - - DIS_NU = ht_func.DIS_SNU(HT, (A, Z), (x, q2, y), Mh=0.938, Mw=80.398, lepton=l) - HT_func['H2p'] = DIS_NU.PC_2_p - HT_func['H2d'] = DIS_NU.PC_2_d - HT_func["HLp"] = DIS_NU.PC_L_p - HT_func["HLd"] = DIS_NU.PC_L_d - HT_func["H3p"] = DIS_NU.PC_3_p - HT_func["H3d"] = DIS_NU.PC_3_d - - # NuTeV - elif process_type in ["DIS_DM_NU", "DIS_DM_NB"]: - # Iron target - Z = 23.403 - A = 49.618 - if process_type == "DIS_SNU": - l = 0 - elif process_type == "DIS_SNB": - l = 1 - - DIS_NuTeV = ht_func.DIS_NUTEV( - HT, (A, Z), (x, q2, y), Mh=0.938, Mw=80.398, lepton=l - ) - HT_func['H2p'] = DIS_NuTeV.PC_2_p - HT_func['H2d'] = DIS_NuTeV.PC_2_d - HT_func["HLp"] = DIS_NuTeV.PC_L_p - HT_func["HLd"] = DIS_NuTeV.PC_L_d - HT_func["H3p"] = DIS_NuTeV.PC_3_p - HT_func["H3d"] = DIS_NuTeV.PC_3_d - - # HERA_CC - elif process_type in ["DIS_CCE", "DIS_CCP"]: - if process_type == "DIS_CCE": - l = 0 - elif process_type == "DIS_CCP": - l = 1 - DIS_CC_HERA = ht_func.DIS_HERA_CC(HT, (x, q2, y), Mh=0.938, Mw=80.398, lepton=l) - HT_func['H2p'] = DIS_CC_HERA.PC_2_p - HT_func['H2d'] = DIS_CC_HERA.PC_2_d - HT_func["HLp"] = DIS_CC_HERA.PC_L_p - HT_func["HLd"] = DIS_CC_HERA.PC_L_d - HT_func["H3p"] = DIS_CC_HERA.PC_3_p - HT_func["H3d"] = DIS_CC_HERA.PC_3_d - else: - raise Exception( - f"The process type `{process_type}` in `{exp_set.name} has not been implemented." - ) - - for ht in HT.keys(): - for idx_node in range(len(HT[ht]['nodes'])): - shifted_list = ht_func.beta_tilde_5pt(HT[ht]['list'], idx_node) - deltas[ht + f"({idx_node})"] = np.append( - deltas[ht + f"({idx_node})"], HT_func[ht](shifted_list) - ) - - # Construct the covariance matrix - covmats = defaultdict(list) - for exp_name_1, exp_idx_1 in start_proc_by_exp.items(): - for exp_name_2, exp_idx_2 in start_proc_by_exp.items(): - s = np.zeros(shape=(ndata_by_exp[exp_name_1], ndata_by_exp[exp_name_2])) - for shifts in deltas.keys(): - s += np.outer( - deltas[shifts][exp_idx_1 : exp_idx_1 + ndata_by_exp[exp_name_1]], - deltas[shifts][exp_idx_2 : exp_idx_2 + ndata_by_exp[exp_name_2]], - ) - - start_locs = (exp_idx_1, exp_idx_2) - covmats[start_locs] = s - - return covmats, deltas - - -def thcov_ht(thcov_shifts_ht, table_ht_deltas): - covmat, _ = thcov_shifts_ht - return covmat - - -@table -def table_ht_deltas(thcov_shifts_ht, procs_index, combine_by_type_custom): - _, deltas = thcov_shifts_ht - process_info = combine_by_type_custom - indexlist = [] - for procname in process_info.preds: - for datasetname in process_info.namelist[procname]: - slicer = procs_index.get_locs((procname, datasetname)) - indexlist += procs_index[slicer].to_list() - covmat_index = pd.MultiIndex.from_tuples(indexlist, names=procs_index.names) - df = pd.DataFrame(deltas, index=covmat_index, columns=deltas.keys()) - return df - - -def covmat_3fpt(name1, name2, deltas1, deltas2): +def covmat_3fpt(deltas1, deltas2): """Returns theory covariance sub-matrix for 3pt factorisation scale variation *only*, given two dataset names and collections of scale variation shifts""" @@ -414,6 +219,40 @@ def covmat_n3lo_ad(name1, name2, deltas1, deltas2): return 1 / norm * s +def covmat_power_corrections(deltas1, deltas2): + """Returns the the theory covariance sub-matrix for power + corrections. The two arguments `deltas1` and `deltas2` contain + the shifts for the firs and second experiment, respectively. + + The shifts are given in this form: + ``` + deltas1 = {shift1_label: array1_of_shifts1, + shift2_label: array1_of_shifts2, + shift3_label: array1_of_shifts3, + ...} + deltas2 = {shift1_label: array2_of_shifts1, + shift2_label: array2_of_shifts2, + shift3_label: array2_of_shifts3, + ...} + ``` + The sub-matrix is computed using the 5-point prescription, thus + + s = array1_of_shifts1 X array2_of_shifts1 + array1_of_shifts2 X array2_of_shifts2 + ... + + where `X` is the outer product. + """ + # Check that `deltas1` and `deltas2` have the same shifts + if deltas1.keys() != deltas2.keys(): + raise RuntimeError('The two dictionaries do not contain the same shifts.') + + size1 = next(iter(deltas1.values())).size + size2 = next(iter(deltas2.values())).size + s = np.zeros(shape=(size1, size2)) + for shift in deltas1.keys(): + s += np.outer(deltas1[shift], deltas2[shift]) + return s + + def compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2=None, deltas2=None): """Utility to compute the covariance matrix by prescription given the shifts with respect to the central value for a pair of processes. @@ -484,46 +323,82 @@ def compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2=None, del # alphas is correlated for all datapoints and the covmat construction is # therefore equivalent to that of the factorization scale variations s = covmat_3fpt(deltas1, deltas2) + elif point_prescription == 'power corrections': + # Shifts computed from power corrected predictions + s = covmat_power_corrections(deltas1, deltas2) return s @check_correct_theory_combination -def covs_pt_prescrip(combine_by_type, point_prescription): +def covs_pt_prescrip( + combine_by_type, + point_prescription, + pdf: PDF, + power_corr_dict, + pc_included_prosc, + pc_excluded_exps, +): """Produces the sub-matrices of the theory covariance matrix according to a point prescription which matches the number of input theories. - If 5 theories are provided, a scheme 'bar' or 'nobar' must be chosen in the runcard in order to specify the prescription. Sub-matrices correspond to applying the scale variation prescription to each pair of processes in turn, using a different procedure for the case where the processes are the same relative to when they are different.""" - process_info = combine_by_type + datagroup_spec = process_info.data_spec running_index = 0 - start_proc = defaultdict(list) - for name in process_info.preds: - size = len(process_info.preds[name][0]) - start_proc[name] = running_index - running_index += size covmats = defaultdict(list) - for name1 in process_info.preds: - for name2 in process_info.preds: - central1, *others1 = process_info.preds[name1] - deltas1 = list(other - central1 for other in others1) - central2, *others2 = process_info.preds[name2] - deltas2 = list(other - central2 for other in others2) - s = compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2, deltas2) - start_locs = (start_proc[name1], start_proc[name2]) - covmats[start_locs] = s + if point_prescription != 'power corrections': + start_proc = defaultdict(list) + for name in process_info.preds: + size = len(process_info.preds[name][0]) + start_proc[name] = running_index + running_index += size + + for name1 in process_info.preds: + for name2 in process_info.preds: + central1, *others1 = process_info.preds[name1] + deltas1 = list(other - central1 for other in others1) + central2, *others2 = process_info.preds[name2] + deltas2 = list(other - central2 for other in others2) + s = compute_covs_pt_prescrip(point_prescription, name1, deltas1, name2, deltas2) + start_locs = (start_proc[name1], start_proc[name2]) + covmats[start_locs] = s + + # For power corrections, the loops run over experimentes + else: + start_proc_by_exp = defaultdict(list) + for exp_name, data_spec in datagroup_spec.items(): + start_proc_by_exp[exp_name] = running_index + running_index += data_spec.load_commondata().ndata + + for exp_name1, data_spec1 in datagroup_spec.items(): + for exp_name2, data_spec2 in datagroup_spec.items(): + process_type1 = process_lookup(exp_name1) + process_type2 = process_lookup(exp_name2) + + is_excluded_exp = any(name in pc_excluded_exps for name in [exp_name1, exp_name2]) + is_included_proc = any( + proc not in pc_included_prosc for proc in [process_type1, process_type2] + ) + if not (is_excluded_exp or is_included_proc): + deltas1 = compute_deltas_pc(data_spec1, pdf, power_corr_dict) + deltas2 = compute_deltas_pc(data_spec2, pdf, power_corr_dict) + s = compute_covs_pt_prescrip( + point_prescription, exp_name1, deltas1, exp_name2, deltas2 + ) + start_locs = (start_proc_by_exp[exp_name1], start_proc_by_exp[exp_name2]) + covmats[start_locs] = s return covmats @table -def theory_covmat_custom(covmat_custom, procs_index, combine_by_type_custom): +def theory_covmat_custom_per_prescription(covs_pt_prescrip, procs_index, combine_by_type): """Takes the individual sub-covmats between each two processes and assembles them into a full covmat. Then reshuffles the order from ordering by process to ordering by experiment as listed in the runcard""" - process_info = combine_by_type_custom + process_info = combine_by_type # Construct a covmat_index based on the order of experiments as they are in combine_by_type # NOTE: maybe the ordering of covmat_index is always the same as that of procs_index? @@ -539,7 +414,7 @@ def theory_covmat_custom(covmat_custom, procs_index, combine_by_type_custom): # Put the covariance matrices between two process into a single covariance matrix total_datapoints = sum(process_info.sizes.values()) mat = np.zeros((total_datapoints, total_datapoints), dtype=np.float32) - for locs, cov in covmat_custom.items(): + for locs, cov in covs_pt_prescrip.items(): xsize, ysize = cov.shape mat[locs[0] : locs[0] + xsize, locs[1] : locs[1] + ysize] = cov df = pd.DataFrame(mat, index=covmat_index, columns=covmat_index) diff --git a/validphys2/src/validphys/theorycovariance/higher_twist_functions.py b/validphys2/src/validphys/theorycovariance/higher_twist_functions.py index bfdc4bbf2f..4a496330b8 100644 --- a/validphys2/src/validphys/theorycovariance/higher_twist_functions.py +++ b/validphys2/src/validphys/theorycovariance/higher_twist_functions.py @@ -1,21 +1,85 @@ +""" +This module contains the utilities for the computation of the shifts in +the theoretical predictions due to the effects of power corrections. Contrary +to scale variations, in the case of power corrections the shifts are not +computed using theories. The shifts are computed at "runtime" during vp-setupfit. + +The aim is that, after shifts being computed, the covmat can be constructed using +the same functions implemented for scale variations (e.g. `covs_pt_prescrip`). +The way shifts are computed depends also on the point prescription. In the case of +scale variations, the point prescription specifies the theories whereby shifts are +computed. In the case of power corrections, shifts and covmat constructions are +computed using a 5-point prescription extended to every parameter used to define +the power correction. + +This module comprehends a bunch of ``factory`` functions such as `DIS_F2_pc`. Each +of these functions returns another function that computes the shifts taking as arguments +the values of the parameters used to parametrise the power corrections. In +other words, these factory functions hard-code the dependence on the kinematic +and leave the dependence on the parameters free. In this way, the shifts +can be computed using different combinations of parameters (i.e. different prescriptions) +if needed. +""" + +from collections import defaultdict import operator import numpy as np import scipy.interpolate as scint from validphys.convolution import central_fk_predictions +from validphys.core import PDF, DataSetSpec +from validphys.process_options import _Process GEV_CM2_CONV = 3.893793e10 GF = 1.1663787e-05 # Fermi's constant [GeV^-2] - - -def beta_tilde_5pt(delta_h, idx): - shifted_list = [0 for _ in range(len(delta_h))] - shifted_list[idx] = delta_h[idx] - return shifted_list - - -def ht_parametrisation(delta_h: list, nodes: list, x: list, Q2: list): +Mh = 0.938 # Proton's mass in GeV/c^2 +MW = 80.398 # W boson mass in GeV/c^2 + +F2P_exps = ['SLAC_NC_NOTFIXED_P_EM-F2', 'BCDMS_NC_NOTFIXED_P_EM-F2'] +F2D_exps = ['SLAC_NC_NOTFIXED_D_EM-F2', 'BCDMS_NC_NOTFIXED_D_EM-F2'] +NC_SIGMARED_P_EM = ['NMC_NC_NOTFIXED_P_EM-SIGMARED', 'HERA_NC_318GEV_EM-SIGMARED'] +NC_SIGMARED_P_EP = [ + 'HERA_NC_225GEV_EP-SIGMARED', + 'HERA_NC_251GEV_EP-SIGMARED', + 'HERA_NC_300GEV_EP-SIGMARED', + 'HERA_NC_318GEV_EP-SIGMARED', +] +NC_SIGMARED_P_EAVG = ['HERA_NC_318GEV_EAVG_CHARM-SIGMARED', 'HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED'] + + +def dis_pc_func(delta_h: list, nodes: list, x: list, Q2: list): + """ + This function defines the parametrization used to model power corrections + for DIS-like processes. Currently, only the cubic spline is supported and + it is hard coded in this function. + + The initialization of the cubic spline requires a list of nodes, which contains + the array of the independent variables (e.g. x-Bjorken), + and a list of shifts that correspond to the dependent variables. Each pair (node, shift) + is a point in the plane. The ensemble of points will be interpolated according to the + cubic spline. + + Parameters + ---------- + delta_h: list + Shifts of the dependent variables at each node listed in `nodes`. These values correspond + to the `amplitude` of the power correction at each node. + nodes: list + List of nodes in the independent variables. For DIS-like processes, these are points + in the x-Bjorken. + x: list + List of x-Bjorken points where the power correction function is evaluated + Q2: list + List of scales where the power correction function is evaluated. Note that this list + is meant to be of the same length as `x`, and the two lists are meant to be considered + as pairs, e.g. (x1, Q2_1), (x2, Q2_2), ... . + + Returns + ------- + A list of power corrections for DIS-like processes where each point is evaluated at the + kinematic pair (x,Q2). + """ H = scint.CubicSpline(nodes, delta_h) H = np.vectorize(H) @@ -23,17 +87,95 @@ def ht_parametrisation(delta_h: list, nodes: list, x: list, Q2: list): return PC -def null_func(size): - """Auxiliary function used to return arrays of zeros for those datasets - that don't require higher twist.""" +# TODO Maybe we want to treat the function that parametrizes the PC +# as argument? +def DIS_F2_pc(pc2_nodes, x, q2): + """ + Returns the function that computes the shift for the ratio of structure + functions F2_d / F2_p. For this observable, power corrections are defined + such that + + F2 -> F2 + PC2, - def zeros(list): - return np.zeros(size) + and the shift is defined as - return zeros + Delta(F2) = (F2 + PC2) - F2 = PC2. + Note that, as in the case of `DIS_F2R_ht`, the shift still depends on the set + of parameters needed to define the parameterization of PC2. Also, this function + can be used to compute the shift for both proton and deuteron, provided that the + correct list of parameters is passed to the **curried** function. -def DIS_F2R_ht(experiment, pdf, H2p_dict, H2d_dict, x, q2): + The function used to parametrize the the power correction is `dis_pc_func` and + it is hard coded. + + Parameters + ---------- + + """ + + def PC_2(y_values): + result = dis_pc_func(y_values, pc2_nodes, x, q2) + return result + + return PC_2 + + +def DIS_F2R_pc(experiment, pdf, pc_2_p_nodes, pc_2_d_nodes, x, q2): + """ + Returns the function that computes the shift for the ratio of structure + functions F2_d / F2_p. For this observable, power corrections are defined + such that + + F2_d / F2_p -> (F2_d + PC2_d) / (F2_p + PC2_p) , + + and the shift is the defined as + + Delta(F2 ratio) = (F2_d + PC2_d) / (F2_p + PC2_p) - F2_d / F2_p . + + The shift is computed for a given set of kinematic variables specified + by the paris (x,Q2), but it still depends on the set of parameters need by + the power correction terms PC2_d and PC2_p. + + Note that this function does **not** return the power corrections for the + given kinematic points, but rather it returns another function where the + kinematic dependence has been **curried** (i.e. hard coded). This new function + takes as arguments the y-values of the nodes used to compute PC2_d and PC2_p + (see `delta_h` in `dis_pc_func`). Note that these y-values are not necessarily + the values listed in the runcard, as we can apply different point prescription. + For instance, we may want to pass in a set of y-values where the nodes are shifted + one at the time, leaving the others zero. The prescription is thus handled separately. + The returning function allows thus to compute Delta(F2 ratio)({...}_d, {...}_p), where + `{...}_d` and `{...}_p` are the sets of y-values for the parametrisation for the proton + and deuteron terms respectively. + + The function used to parametrize the the power correction is `dis_pc_func` and + it is hard coded. + + Parameters + ---------- + experiment: DataSetSpec + An instance of DataSetSpec used to extract information such as cuts + and fk tables. + pdf: PDF + An instance of the class PDF. This specifies the PDF to bo convoluted + with the FK tables. + pc_2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for the proton (see `dis_pc_func`). + pc_2_d_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for the deuteron (see `dis_pc_func`). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2_d and P2_p. + """ cuts = experiment.cuts fkspec_F2D, fkspec_F2P = experiment.fkspecs fk_F2D = fkspec_F2D.load_with_cuts(cuts) @@ -45,120 +187,620 @@ def DIS_F2R_ht(experiment, pdf, H2p_dict, H2d_dict, x, q2): F2P = np.concatenate(F2P.values) F2_ratio = operator.truediv(F2D, F2P) - def PC_2_p(list): - PC_p = ht_parametrisation(list, H2p_dict["nodes"], x, q2) - result = np.array(operator.truediv(F2D, np.sum([F2P, PC_p], axis=0)) - F2_ratio) + def func(y_values_d, y_values_p): + PC_d = dis_pc_func(y_values_d, pc_2_d_nodes, x, q2) + PC_p = dis_pc_func(y_values_p, pc_2_p_nodes, x, q2) + num = np.sum([F2D, PC_d], axis=0) + denom = np.sum([F2P, PC_p], axis=0) + result = np.array(operator.truediv(num, denom) - F2_ratio) return result - def PC_2_d(list): - PC_d = ht_parametrisation(list, H2d_dict["nodes"], x, q2) - result = np.array(operator.truediv(np.sum([F2D, PC_d], axis=0), F2P) - F2_ratio) - return result - - return PC_2_p, PC_2_d - + return func # PC_2_p, PC_2_d + + +def DIS_F2C_pc(pc2_p_nodes, pc2_d_nodes, x, q2): + """ + Builds the function used to compute the shifts for the charm + structure function measured by EMC. The process involved is + + mu^+ + Fe -> mu+^ + c cbar + X . + + This function works exactly as the previous functions used to + compute nuisance shifts. In this case, the constructed function + (`func` below) requires two lists of parameters for the proton + and the deuteron contribution. The reason being that in this process + the muon scatters off an iron target, and the power correction + contribution is a mixture of proton and deuteron nucleons. Hence, proton + and deuteron contribution are weighted by the appropriate atomic factor. + + Note that we are parametrising power corrections as proton and deuteron + targets. If we were to parametrize such contributions using, say, proton + and nucleon, than the weights would change. + + Parameters + ---------- + pc2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for the proton (see `dis_pc_func`). + pc2_d_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for the deuteron (see `dis_pc_func`). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2_d and P2_p. + """ + # Iron target + Z = 23.403 + A = 49.618 -def DIS_F2_ht(H2_dict, x, q2): - def PC_2(list): - result = ht_parametrisation(list, H2_dict['nodes'], x, q2) + def func(y_values_d, y_values_p): + PC2_d = dis_pc_func(y_values_d, pc2_d_nodes, x, q2) + PC2_p = dis_pc_func(y_values_p, pc2_p_nodes, x, q2) + result = 2 * (Z - A) / A * PC2_d + Z * PC2_p return result - return PC_2 + return func + + +def DIS_NC_XSEC_pc(pc2_nodes, pcL_nodes, pc3_nodes, lepton, x, q2, y): + """ + Builds the function used to compute the shifts for the DIS NC x-secs + delivered by HERA and NMC. The x-sec is reconstructed as calculated + in Yadism (see https://yadism.readthedocs.io/en/latest/theory/intro.html). + In particular, the x-sec is a linear combination of the structure functions + F_2, F_L, and F_3. The coefficients are also computed appropriately (see + link). The contribution of the power corrections is then + + Delta(x-sec) = x-sec_w_pc - x-sec_wo_pc = PC_2 + N_L * PC_L + N_3 * PC_3 + + where PC_i are the power corrections relative to the respective structure + functions and the N_i the respective coefficients (as defined in Yadism). + + This function works exactly as the previous functions used to + compute nuisance shifts. In addition, it requires the kinematic + invariant `y` to build the shift-function. + + Note that this function can be used for both proton and deuteron targets, + provided that the appropriate lists of nodes is given. + + Parameters + ---------- + pc2_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_2. + pcL_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_L. + pc3_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_3. + lepton: int + Whether the scattering particle is a lepton (0) or an anti-lepton(1). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + y: list[float] + Set of points in y where the power corrections will be evaluated. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2 and PL and P3. + """ + yp = 1 + np.power(1 - y, 2) + ym = 1 - np.power(1 - y, 2) + yL = np.power(y, 2) + N_L = -yL / yp # Coefficient for F_L + N_3 = np.power(-1, lepton) * ym / yp # Coefficient for F_3 + + def func(y_values_pc2, y_values_pcL, y_values_pc3): + PC_2 = dis_pc_func(y_values_pc2, pc2_nodes, x, q2) + PC_L = dis_pc_func(y_values_pcL, pcL_nodes, x, q2) + PC_3 = dis_pc_func(y_values_pc3, pc3_nodes, x, q2) + result = PC_2 + N_L * PC_L + N_3 * PC_3 + return result + return func + + +def DIS_CC_HERA_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, lepton, x, q2, y): + """ + Builds the function used to compute the shifts for the DIS CC x-secs + delivered by HERA. The x-sec is reconstructed as calculated + in Yadism (see https://yadism.readthedocs.io/en/latest/theory/intro.html). + In particular, the x-sec is a linear combination of the structure functions + F_2, F_L, and F_3. The coefficients are also computed appropriately (see + link). The contribution of the power corrections is then + + Delta(x-sec) = x-sec_w_pc - x-sec_wo_pc = N * (PC_2 + N_L * PC_L + N_3 * PC_3) + + where PC_i are the power corrections relative to the respective structure + functions and the N_i the respective coefficients (as defined in Yadism). + N is the overall normalization factor. + + For the HERA_CC_318GEV dataset, the target is always a proton. However, the + lepton may be either the electron (0) or the positron (1). This information + is needed in order to compute the coefficient N_3. + + Parameters + ---------- + pc2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_2. + pcL_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_L. + pc3_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_3. + lepton: int + Whether the scattering particle is a lepton (0) or an anti-lepton(1). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + y: list[float] + Set of points in y where the power corrections will be evaluated. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2 and PL and P3. + """ + yp = 1 + np.power(1 - y, 2) + ym = 1 - np.power(1 - y, 2) + yL = np.power(y, 2) + N = 1 / 4 * yp # Overall normalization + N_L = -yL / yp # Coefficient for F_L + N_3 = np.power(-1, lepton) * ym / yp # Coefficient for F_3 + + def func(y_values_pc2_p, y_values_pcL_p, y_values_pc3_p): + # Initialize power corrections for each structure function + PC2_p = dis_pc_func(y_values_pc2_p, pc2_p_nodes, x, q2) + PCL_p = dis_pc_func(y_values_pcL_p, pcL_p_nodes, x, q2) + PC3_p = dis_pc_func(y_values_pc3_p, pc3_p_nodes, x, q2) + + # Build the contribution to the x-sec of the power corrections + result = N * (PC2_p + N_L * PCL_p + N_3 * PC3_p) + return result -def DIS_F2C_ht(H2p_dict, H2d_dict, x, q2): + return func + + +def DIS_CC_NUTEV_pc( + pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, pc2_d_nodes, pcL_d_nodes, pc3_d_nodes, lepton, x, q2, y +): + """ + Builds the function used to compute the shifts for the DIS CC x-secs + delivered by NuTeV. The x-sec is reconstructed as calculated + in Yadism (see https://yadism.readthedocs.io/en/latest/theory/intro.html). + In particular, the x-sec is a linear combination of the structure functions + F_2, F_L, and F_3. The coefficients are also computed appropriately (see + link). Note that this experiment uses iron targets, and thus the coefficients + must take into account the nuclear mixture of porton and deuteron. The contribution + of the power corrections is then + + Delta(x-sec) = x-sec_w_pc - x-sec_wo_pc = N * (PC_2 + N_L * PC_L + N_3 * PC_3) + + where PC_i are the power corrections relative to the respective structure + functions (nuclear mixture implicit) and the N_i the respective coefficients (as defined in Yadism). + N is the overall normalization factor. + + For the NuTeV CC dataset, the target is always iron. However, the + lepton may be either the electron (0) or the positron (1). This + information is needed in order to compute the coefficient N_3. + + Parameters + ---------- + pc2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_2 of the proton. + pcL_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_L of the proton. + pc3_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_3 of the proton. + pc2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_2 of the deuteron. + pcL_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_L of the deuteron. + pc3_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_3 of the deuteron. + lepton: int + Whether the scattering particle is a lepton (0) or an anti-lepton(1). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + y: list[float] + Set of points in y where the power corrections will be evaluated. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2 and PL and P3 for proton and deuteron. + """ # Iron target Z = 23.403 A = 49.618 - - def PC_2_p(list): - result = ht_parametrisation(list, H2p_dict['nodes'], x, q2) + nuclear_factor = 2 * (Z - A) / A + yp = 1 + np.power(1 - y, 2) - 2 * np.power(x * y * Mh, 2) / q2 + ym = 1 - np.power(1 - y, 2) + yL = np.power(y, 2) + N_L = -yL / yp # Coefficient for F_L + N_3 = np.power(-1, lepton) * ym / yp # Coefficient for F_3 + + MW2 = np.power(MW, 2) + # Overall coefficient + # TODO: cross-check + N = 100 * yp / (2 * np.power(1 + q2 / MW2, 2)) + + def func( + y_values_pc2_p, + y_values_pcL_p, + y_values_pc3_p, + y_values_pc2_d, + y_values_pcL_d, + y_values_pc3_d, + ): + PC2_p = dis_pc_func(y_values_pc2_p, pc2_p_nodes, x, q2) + PCL_p = dis_pc_func(y_values_pcL_p, pcL_p_nodes, x, q2) + PC3_p = dis_pc_func(y_values_pc3_p, pc3_p_nodes, x, q2) + PC2_d = dis_pc_func(y_values_pc2_d, pc2_d_nodes, x, q2) + PCL_d = dis_pc_func(y_values_pcL_d, pcL_d_nodes, x, q2) + PC3_d = dis_pc_func(y_values_pc3_d, pc3_d_nodes, x, q2) + tmp_2 = Z * PC2_p + nuclear_factor * PC2_d + tmp_L = Z * PCL_p + nuclear_factor * PCL_d + tmp_3 = Z * PC3_p + nuclear_factor * PC3_d + result = N * (tmp_2 + N_L * tmp_L + N_3 * tmp_3) return result - def PC_2_d(list): - result = 2 * (Z - A) / A * ht_parametrisation(list, H2d_dict['nodes'], x, q2) + return func + + +# TODO This is function is really similar to the one +# defined for NUTEV CC. Can we reduce code repetitions? +def DIS_CC_CHORUS_pc( + pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, pc2_d_nodes, pcL_d_nodes, pc3_d_nodes, lepton, x, q2, y +): + """ + Same as DIS_CC_NUTEV_pc, but for CHORUS CC. + + Note that the difference here is in the definition of the overall + normalization N. + + Parameters + ---------- + pc2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_2 of the proton. + pcL_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_L of the proton. + pc3_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_3 of the proton. + pc2_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_2 of the deuteron. + pcL_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_L of the deuteron. + pc3_p_nodes: list[float] + The list of nodes in x-Bjorken used to define the parametrization of the + power correction for F_3 of the deuteron. + lepton: int + Whether the scattering particle is a lepton (0) or an anti-lepton(1). + x: list[float] + Set of points in x-Bjorken where the power corrections will be evaluated. + q2: list[float] + Set of points in Q2 where the power corrections will be evaluated. + y: list[float] + Set of points in y where the power corrections will be evaluated. + + Returns + ------- + The function the computes the shift for this observable. It depends on the + y-values for the parameterization of P2 and PL and P3 for proton and deuteron. + """ + # Lead target + A = 208.0 + Z = 82 + nuclear_factor = 2 * (Z - A) / A + yp = 1 + np.power(1 - y, 2) - 2 * np.power(x * y * Mh, 2) / q2 + ym = 1 - np.power(1 - y, 2) + yL = np.power(y, 2) + N_L = -yL / yp # Coefficient for F_L + N_3 = np.power(-1, lepton) * ym / yp # Coefficient for F_3 + + MW2 = np.power(MW, 2) + # Overall coefficient + # TODO: cross-check + N = GEV_CM2_CONV * (GF**2) * Mh / (2 * np.pi * np.power(1 + q2 / MW2, 2)) * yp + + def func( + y_values_pc2_p, + y_values_pcL_p, + y_values_pc3_p, + y_values_pc2_d, + y_values_pcL_d, + y_values_pc3_d, + ): + PC2_p = dis_pc_func(y_values_pc2_p, pc2_p_nodes, x, q2) + PCL_p = dis_pc_func(y_values_pcL_p, pcL_p_nodes, x, q2) + PC3_p = dis_pc_func(y_values_pc3_p, pc3_p_nodes, x, q2) + PC2_d = dis_pc_func(y_values_pc2_d, pc2_d_nodes, x, q2) + PCL_d = dis_pc_func(y_values_pcL_d, pcL_d_nodes, x, q2) + PC3_d = dis_pc_func(y_values_pc3_d, pc3_d_nodes, x, q2) + tmp_2 = Z * PC2_p + nuclear_factor * PC2_d + tmp_L = Z * PCL_p + nuclear_factor * PCL_d + tmp_3 = Z * PC3_p + nuclear_factor * PC3_d + result = N * (tmp_2 + N_L * tmp_L + N_3 * tmp_3) return result - return PC_2_p, PC_2_d - - -def DIS_NC_ht(H2_dict, HL_dict, x, q2, y): - yp = 1 + np.power(1 - y, 2) - yL = np.power(y, 2) - N_L = -yL / yp - - def PC_2(list): - return ht_parametrisation(list, H2_dict['nodes'], x, q2) - - def PC_L(list): - return N_L * ht_parametrisation(list, HL_dict['nodes'], x, q2) - - return PC_2, PC_L - - -class DIS_SNU: - def __init__(self, HT_dict, target_tuple, kin_tuple, Mh, Mw, lepton): # (A,Z) - - # Lead target - A = target_tuple[0] - Z = target_tuple[1] - x = kin_tuple[0] - q2 = kin_tuple[1] - y = kin_tuple[2] - self.nuclear_target = 2 * (Z - A) / A - self.Mh = Mh # 0.938 GeV - self.Mw2 = np.power(Mw, 2) # GeV^2 - self.yp = 1 + np.power(1 - y, 2) - 2 * np.power(x * y * Mh, 2) / q2 - self.yL = np.power(y, 2) - self.ym = 1 - np.power(1 - y, 2) - self.N = ( - GEV_CM2_CONV * (GF**2) * Mh / (2 * np.pi * np.power(1 + q2 / self.Mw2, 2)) * self.yp - ) - self.H2p_dict = HT_dict['H2p'] - self.H2d_dict = HT_dict['H2d'] - self.HLp_dict = HT_dict['HLp'] - self.HLd_dict = HT_dict['HLd'] - self.H3p_dict = HT_dict['H3p'] - self.H3d_dict = HT_dict['H3d'] - - self.x = x - self.q2 = q2 - self.l = lepton - - def PC_2_p(self, list): - norm = self.N - return norm * ht_parametrisation(list, self.H2p_dict['nodes'], self.x, self.q2) - - def PC_2_d(self, list): - norm = self.N * self.nuclear_target - return norm * ht_parametrisation(list, self.H2d_dict['nodes'], self.x, self.q2) - - def PC_L_p(self, list): - norm = -self.N * self.yL / self.yp - return norm * ht_parametrisation(list, self.HLp_dict['nodes'], self.x, self.q2) - - def PC_L_d(self, list): - norm = -self.N * self.yL / self.yp * self.nuclear_target - return norm * ht_parametrisation(list, self.HLd_dict['nodes'], self.x, self.q2) - - def PC_3_p(self, list): - norm = self.N * np.power(-1, self.l) * self.ym / self.yp * self.x - return norm * ht_parametrisation(list, self.H3p_dict['nodes'], self.x, self.q2) - - def PC_3_d(self, list): - norm = self.N * np.power(-1, self.l) * self.ym / self.yp * self.x * self.nuclear_target - return norm * ht_parametrisation(list, self.H3d_dict['nodes'], self.x, self.q2) - - -class DIS_NUTEV(DIS_SNU): - def __init__(self, HT_dict, target_tuple, kin_tuple, Mh, Mw, lepton): # (A,Z) - super().__init__(HT_dict, target_tuple, kin_tuple, Mh, Mw, lepton) - self.N = 50 * self.yp / np.power(1 + self.q2 / self.Mw2, 2) - - -class DIS_HERA_CC(DIS_SNU): - def __init__(self, HT_dict, kin_tuple, Mh, Mw, lepton): - super().__init__(HT_dict, (1, 1), kin_tuple, Mh, Mw, lepton) - y = kin_tuple[2] - self.yp = 1 + np.power(1 - y, 2) - N = 1 / 4 * self.yp + return func + + +def construct_pars_combs(parameters_dict): + """Construct the combination of parameters (the ones that parametrize the power + corrections) used to compute the shifts. + + Example + ------- + Given the following dictionary that specifies that power corrections + ``` + pc_dict = { + 'H1': {'list': [1,2], 'nodes': [0,1]} } + 'H2': {'list': [3,4,5], 'nodes': [0,1,2]} } + } + ``` + then this functions constructs a list as follows + ``` + pars_combs = [ + {'label': 'H1(1)', 'comb': {'H1': [1,0], 'H2': [0,0,0]}, + {'label': 'H1(2)', 'comb': {'H1': [0,1], 'H2': [0,0,0]}, + {'label': 'H2(1)', 'comb': {'H1': [0,0], 'H2': [3,0,0]}, + {'label': 'H2(2)', 'comb': {'H1': [0,0], 'H2': [0,4,0]}, + {'label': 'H2(3)', 'comb': {'H1': [0,0], 'H2': [0,0,5]}, + ] + ``` + + Parameters + ---------- + """ + combinations = [] + for key, values in parameters_dict.items(): + for i in range(len(values['yshift'])): + # Create a new dictionary with all keys and zeroed-out values + new_dict = {k: np.zeros_like(v['yshift']) for k, v in parameters_dict.items()} + # Set the specific value for the current index + label = key + f'({i})' + new_dict[key][i] = values['yshift'][i] + new_dict = {'label': label, 'comb': new_dict} + combinations.append(new_dict) + + return combinations + + +def compute_deltas_pc(dataset_sp: DataSetSpec, pdf: PDF, power_corr_dict: dict): + """ + Computes the shifts due to power corrections for a single dataset given + the set of parameters that model the power corrections. The result is + a dictionary containing as many arrays of shifts as the number of + combinations of the parameters. For instance, the final dictionary + may look like: + ``` + deltas1 = {comb1_label: array_of_shifts1, + comb2_label: array_of_shifts2, + comb3_label: array_of_shifts3, + ...} + ``` + Note that, as of now, we don't need to specify different prescriptions. + For that reason, the prescription adopted to construct the shifts is hard + coded in the function `construct_pars_combs`, and the prescription used to + compute the sub-matrix is hard-coded in `covmat_power_corrections`. + """ + + exp_name = dataset_sp.name + cd_table = dataset_sp.load_commondata().commondata_table + process_type = cd_table['process'].iloc[0] + if isinstance(process_type, _Process): + process_type = process_type.name + + pars_combs = construct_pars_combs(power_corr_dict) + deltas = defaultdict(list) + + pc_func = None + if process_type.startswith('DIS'): + pc2_p_nodes = power_corr_dict["H2p"]['nodes'] + pcL_p_nodes = power_corr_dict["HLp"]['nodes'] + pc3_p_nodes = power_corr_dict["H3p"]['nodes'] + pc2_d_nodes = power_corr_dict["H2d"]['nodes'] + pcL_d_nodes = power_corr_dict["HLd"]['nodes'] + pc3_d_nodes = power_corr_dict["H3d"]['nodes'] + + # TODO + # AFter the data re-implementation the name of the variables + # in the commondata table will change as indicated in the metadata. + # When this happens, this part must be updated. + x = cd_table['kin1'].to_numpy() + q2 = cd_table['kin2'].to_numpy() + y = cd_table['kin3'].to_numpy() + + # F2 ratio + if exp_name == "NMC_NC_NOTFIXED_EM-F2": + pc_func = DIS_F2R_pc(dataset_sp, pdf, pc2_p_nodes, pc2_d_nodes, x, q2) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['H2p'], pars_pc['comb']['H2d']) + + # F2 proton traget + elif exp_name in F2P_exps: + pc_func = DIS_F2_pc(pc2_p_nodes, x, q2) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['H2p']) + + # F2 deuteron traget + elif exp_name in F2D_exps: + pc_func = DIS_F2_pc(pc2_d_nodes, x, q2) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['H2d']) + + # EMC + elif exp_name.startswith('EMC_NC_250GEV'): + pc_func = DIS_F2C_pc(pc2_p_nodes, pc2_d_nodes, x, q2) + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func(pars_pc['comb']['H2p'], pars_pc['comb']['H2d']) + + # HERA and NMC SIGMARED NC + elif exp_name in np.concatenate([NC_SIGMARED_P_EM, NC_SIGMARED_P_EP, NC_SIGMARED_P_EAVG]): + # Electron + if exp_name in NC_SIGMARED_P_EM: + pc_func = DIS_NC_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, 0, x, q2, y) + # Positron + elif exp_name in NC_SIGMARED_P_EP: + pc_func = DIS_NC_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, 1, x, q2, y) + # Average positron and electron + # TODO + # Check if this is correct (ach) + elif NC_SIGMARED_P_EAVG: + + def average(y_values_pc2_p, y_values_pcL_p, y_values_pc3_p): + electron = DIS_NC_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, 0, x, q2, y) + positron = DIS_NC_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, 1, x, q2, y) + return ( + np.sum( + electron(y_values_pc2_p, y_values_pcL_p, y_values_pc3_p), + positron(y_values_pc2_p, y_values_pcL_p, y_values_pc3_p), + ) + / 2 + ) + + pc_func = average + else: + raise ValueError(f"{exp_name} not implemented.") + + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func( + pars_pc['comb']['H2p'], pars_pc['comb']['HLp'], pars_pc['comb']['H3p'] + ) + + # CHORUS + elif exp_name.startswith('CHORUS_CC'): + # Nu + if exp_name == 'CHORUS_CC_NOTFIXED_PB_NU-SIGMARED': + pc_func = DIS_CC_CHORUS_pc( + pc2_p_nodes, + pcL_p_nodes, + pc3_p_nodes, + pc2_d_nodes, + pcL_d_nodes, + pc3_d_nodes, + 0, + x, + q2, + y, + ) + # Nu bar + elif exp_name == 'CHORUS_CC_NOTFIXED_PB_NB-SIGMARED': + pc_func = DIS_CC_CHORUS_pc( + pc2_p_nodes, + pcL_p_nodes, + pc3_p_nodes, + pc2_d_nodes, + pcL_d_nodes, + pc3_d_nodes, + 1, + x, + q2, + y, + ) + else: + raise ValueError(f"{exp_name} not implemented.") + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func( + pars_pc['comb']['H2p'], + pars_pc['comb']['HLp'], + pars_pc['comb']['H3p'], + pars_pc['comb']['H2d'], + pars_pc['comb']['HLd'], + pars_pc['comb']['H3d'], + ) + + # NuTeV + elif exp_name.startswith('NUTEV_CC'): + # Nu + if exp_name == 'NUTEV_CC_NOTFIXED_FE_NU-SIGMARED': + pc_func = DIS_CC_NUTEV_pc( + pc2_p_nodes, + pcL_p_nodes, + pc3_p_nodes, + pc2_d_nodes, + pcL_d_nodes, + pc3_d_nodes, + 0, + x, + q2, + y, + ) + # Nu bar + elif exp_name == 'NUTEV_CC_NOTFIXED_FE_NB-SIGMARED': + pc_func = DIS_CC_NUTEV_pc( + pc2_p_nodes, + pcL_p_nodes, + pc3_p_nodes, + pc2_d_nodes, + pcL_d_nodes, + pc3_d_nodes, + 1, + x, + q2, + y, + ) + else: + raise ValueError(f"{exp_name} not implemented.") + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func( + pars_pc['comb']['H2p'], + pars_pc['comb']['HLp'], + pars_pc['comb']['H3p'], + pars_pc['comb']['H2d'], + pars_pc['comb']['HLd'], + pars_pc['comb']['H3d'], + ) + + # HERA_CC + elif exp_name.startswith('HERA_CC'): + # electron + if exp_name == 'HERA_CC_318GEV_EM-SIGMARED': + pc_func = DIS_CC_HERA_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, 0, x, q2, y) + # positron + elif exp_name == 'HERA_CC_318GEV_EP-SIGMARED': + pc_func = DIS_CC_HERA_XSEC_pc(pc2_p_nodes, pcL_p_nodes, pc3_p_nodes, 1, x, q2, y) + else: + raise ValueError(f"{exp_name} not implemented.") + + for pars_pc in pars_combs: + deltas[pars_pc['label']] = pc_func( + pars_pc['comb']['H2p'], pars_pc['comb']['HLp'], pars_pc['comb']['H3p'] + ) + + else: + raise ValueError( + f"The {process_type} observable for {exp_name} " "has not been implemented." + ) + + elif process_type == 'JET': + raise NotImplementedError("This part has not been implemented yet.") + + return deltas