New ht thcovmat

[achiefa]

This PR allows for the inclusion of theory uncertainties due to the effect of power corrections. The theory covariance matrix is constructed by computing the shifts for the theoretical predictions as done for the MHOUs. The shift is computed at the level of the structure functions. Then, the shifts for the structure functions are combined to reconstruct the shift for the xsec. For this reason, the calculation of the shift depends on the dataset. Currently, only 1-JET and DIS (NC and CC) are supported.

At the level of the runcard, power corrections are specified as follows

theorycovmatconfig:
  point_prescriptions: ["power corrections"]
  pc_parameters:
  - {ht: H2p, yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
  - {ht: H2d, yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
  - {ht: HLp, yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
  - {ht: HLd, yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
  - {ht: H3p, yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
  - {ht: H3d, yshift: [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0], nodes: [0.0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1]}
  - {ht: Hj, yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]}
  func_type: "linear" 
  pc_included_procs: ["JETS", "DIS NC", "DIS CC"]
  pc_excluded_exps: [HERA_NC_318GEV_EAVG_CHARM-SIGMARED,
                     HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED]
  pdf: 240701-02-rs-nnpdf40-baseline
  use_thcovmat_in_fitting: true
  use_thcovmat_in_sampling: true

For each process/sf implemented, we need to specify a series of information:

1. The name of the power correction, under the key ht. This is useful to identify the type of power correction, in particular when computing the posterior.
2. The values in yshift are the magnitudes of the prior.
3. nodes contains the points where the prior is shifted

The array pc_included_procs specifies the processes for which the shifts are computed. I've also implemented the possibility to exclude some particular datasets within the selected processes, and this can be done by specifying the names in pc_excluded_exps.

The key func_type is temporary and will be deleted once we decide which function to use to construct the prior.

All the relevant details for this PR are contained in the module higher_twist_functions.py. For each observable for which the shift has to be computed, I implemented a factory that constructs a function which will then compute the shift. I thought it was kind of necessary to make the shift dependent on the shift parameters (yshift and nodes) and on the prescription according to which we vary the parameters (which for now is fixed), and not on the kinematics (I think this is known as currying in computer science).

TO DO

• Add documentation
• Remove func_type
• What happens if power_corrections is specified but the parameters are not given?
• Change the name ht to name (?)
• Something else?

[github-actions]

Greetings from your nice fit 🤖 !
I have good news for you, I just finished my tasks:

• Fit Name: NNBOT-2ad60a7b9-2024-07-15
• Fit Report wrt master: https://vp.nnpdf.science/f0TNPeX3SlO5uwbJkr7I8Q==
• Fit Report wrt latest stable reference: https://vp.nnpdf.science/X6jx-wHfT5iJuz14nPatnQ==
• Fit Data: https://data.nnpdf.science/fits/NNBOT-2ad60a7b9-2024-07-15.tar.gz

Check the report carefully, and please buy me a ☕ , or better, a GPU 😉!

[RoyStegeman]

keyword arguments (kwargs) you don't use at this level, you don't have to specify. That's what **kwargs is for.

[RoyStegeman]

Maybe reference Eq. 6 of https://arxiv.org/pdf/1701.05838

[RoyStegeman]

Where is this function? Also I think you can do N2=NL=1, but not NL,N2=1

[achiefa]

The DIS.py is not in sync with the new way shifts are computed. Hence, that function that adds the shift in the theory predictions can be removed. I'll work on that later on, as we discussed. I should have deleted them.

[RoyStegeman]

Why do we need this here?

[achiefa]

I don't remember why I implemented that...I think I can remove it

[RoyStegeman]

Ah yes, we may want to start thinking of a long term solution as well. A tree of if-statments that grows exponentially for each source of theory uncertainty (scales, alphas, HT) is not super nice

[achiefa]

Hi @roy, I've updated the code so that now the covmat for power corrections can be constructed as in the case of scale variations. Please, look and let me know if something is unclear. I'd like if you could double-check the functions that compute the shifts, in particular when normalisation factors and conversion factors are used.

There are few things that I still don't like here and there, but you're free to propose modifications.

[achiefa]

I had to add additional dependencies here for the power correction case. I don't know if there is a clever solution, maybe extracting these quantities from the context?

[achiefa]

I think the tests fail because there is not default for pc_included_procs (typo) and pc_excluded_exps. What do you suggest to do? Should I add two functions in config.py that default these two values when not specified in the runcard?

[github-actions]

Greetings from your nice fit 🤖 !
I have good news for you, I just finished my tasks:

• Fit Name: NNBOT-701fd3dd4-2025-01-14
• Fit Report wrt master: https://vp.nnpdf.science/LhL3COHERwuMwHEEeGvymg==
• Fit Report wrt latest stable reference: https://vp.nnpdf.science/-y7-UpCPSFGPxzP_41VXLA==
• Fit Data: https://data.nnpdf.science/fits/NNBOT-701fd3dd4-2025-01-14.tar.gz

Check the report carefully, and please buy me a ☕ , or better, a GPU 😉!

[scarlehoff]

Hi @achiefa thanks for this.

I've left a few comments on the code. The 1000line higher_twist_functions I haven't look into in much detail, but one thing that I'd like to mention is please avoid using commondata_table. This is basically an object containing all commondata information like in the old days which was necessary to add for some of validphys functions but the idea is to instead get only the information you want.

For instance, the process type is part of the metadata now (no need to load the data and uncertainties!)
Same for the kinematics (although by the time you need the kinematics you often have already read the whole commondata). You can get the kinematics directly with cd.metadata.load_kinematics() and then you already have x, Q, and y for DIS as god intended!

[scarlehoff]

Please remove the comments (and modify the comments at the top of the file to explain what this runcard is).

The meta information is also missing.

[scarlehoff]

This information should come from the theory I believe?

[achiefa]

Yes, I agree. But I'm afraid that fetching this information from the theory would overcomplicate the code. But maybe you know a better way to do so (?)

[scarlehoff]

from validphys.loader import Loader
tid = Loader().check_theoryID(40_000_000)
mw = tid.get_description()["MW"]

You can also load the theories directly without checking them (which will download the full theory) but I guess you have already some theory locally to use these functions?