Simultaneous fit to KA and PI production

[achiefa]

This PR implements the possibility of performing simultaneous fits of multiple hadrons.

Implementation detail

NNADparameterisation

We use the same neural network to parametrise multiple hadrons. This means that the outputs of the network needs to be split into flavours belonging to different hadrons. Hence, for each hadron we define a split matrix which is stored in the class attribute _SplitMatrices that is defined in the constructor of NNADparameterisation. _SplitMatrices is a std::map that maps a std::string (the name of the hadron) to the nnad:Matrix that encodes the filtering layer for the network outputs. Hence, given a neural network that parametrises KA and PI, we can select the flavours of the PI as follows

// Initialise NNAD ...
std::vector<double> nnx = NN.evaluate({0.5})
nnad::Matrix<double> SplitOutput = _SplitMatrices['PI'] * nnad::Matrix(nnx.size(), 1, nnx);

Then, one can apply flavour rotations to SplitOutput. Note that flavour rotations also need to be coupled with the hadron. Hence, we also include the map _Rotations, which is similar to _SplitMatrices but for flavour rotations.

Both _Rotations and _SplitMatrices are constructed in the constructor of NNADparameterisation. The reason is that both require runtime information taken from the runcard. For instance, the split matrices have to account for the number of flavours in each hadron and the total size of the network outputs, both specified in the runcard.
On the other hand, flavour rotations require the flavour map specified for each hadron (see "Runcard Example" to learn how to specify maps for different hadrons).

The data structure of _NNderivativeSets is also modified. Now it is

std::map<std::string, std::vector<apfel::Set<apfel::Distribution>>> _NNderivativeSets;

so that each hadron has its own apfel set. This requires changes in EvaluateOnGrid and DeriveOnGrid, which will loop over hadrons (see source code).

Chi2 and predictions

Since we have to handle predictions with different hadrons in the same fit, we need to ensure that convolutions are computed with the correct hadron apfel::Set. This is achieved by matching the hadron name from the DataHandler with the correct apfel::Set in NNADparameterisation when setting the FF in NangaParbat::ConvolutionTable. For instance, in AnalyticChiSquare we have

nalyticChiSquare::AnalyticChiSquare(std::vector<std::pair<NangaParbat::DataHandler *, NangaParbat::ConvolutionTable *>> DSVect,
                                       NangaParbat::Parameterisation * FFs,
                                       bool NUMERIC):
    NangaParbat::ChiSquare({}, FFs)
  {
    for (auto const &ds : DSVect)
      {
        try
          {
            ds.second->SetInputFFs(FFs->DistributionFunction(ds.first->GetHadron()));
          }
        catch (std::exception const &e)
          {
            std::cerr << e.what() << std::endl;
            std::cerr << "... " << ds.first->GetName() << std::endl;
            exit(1);
          }
      }

LHAPDFgrid

The script that creates the LHAPDF grids is also modified to handle multiple hadrons. Specifically, LHAPDFgrid will generate a grid for each hadron. The names of the sets are specified in the runcard for each hadron, under the key SetName (see runcard below).

Other improvements

Using multiple cores in ceres-solver (see Optimize.cc)

options.num_threads = numCores;

It is possible to specify the number of cores to use as a command-line argument

Optimize 2 MAPFF30_PI_KA_NNLO.yaml ../data fits --jobs 8

TODO

• Inter-hadron sum rules
• Compatibility with scripts in the Tests folder Many scripts in the Tests folder were already deprecated or outdated.
• Check that the implementation gives sensible results Test performed fitting KA and PI without sum rules. The two hadrons are compatible with previous separate fits (see ...)
• Backwards compatibility legacy version is tagged and the user is asked to use the old version of the code with legacy runcard

Runcard example

#################################
# Predictions
#################################
Predictions:
# Perturbative order, 0: LO, 1: NLO
perturbative order: 2

# Initial scale in GeV to be used for the DGLAP evolution of the FFs.
mu0: 5

# Quark thresholds
thresholds: [0, 0, 0, 1.51, 4.92]

# Strong coupling
alphas:
  aref: 0.118
  Qref: 91.1876

# Fine-structure constant
alphaem:
  aref: 0.00776578395589
  Qref: 91.1876

# APFEL++ grid
xgrid:
 - [100, 0.001, 3]
 - [33, 0.1023293, 3]
 - [48, 0.81283052, 3]

zgrid:
 - [100, 0.01, 3]
 - [34, 0.20892961, 3]
 - [48, 0.83176377, 3]

# PDF set
pdfset:
  name: NNPDF31_nnlo_pch_as_0118
  member: -1 #N>=0 for a specific member (0 for central) ; N<0 for all members to be used randomly(flat dist)
  
#################################
# Optimiser
#################################
# Parameters of the optimiser managed by ceres-solver
Optimizer:
max_num_iterations: 3000
chi2_tolerance: 3
#################################
# NNAD
#################################
NNAD:
# Initialisation seed
seed: 4
# Architecture
architecture: [1, 25, 14]
# The output function can be either the activation function of the
# NN (0), or linear (1), or quadratic (2). If this entry is absent a
# linear function is assumed.
output function: 2
# The flavour map gives the the specific combinations of
# distributions in the physical-basis (d, u, s, ..) to be fitted to
# the data. The number of combinations has to match the number of
# nodes of the output layer of the NN given in the
# architecture. When defining the flavour map, one should keep in
# mind that the code computes predictions using the QCD-evolution
# basis (Sigma, V, T3, ...). Therefore, in general,
# QCD-evolution-like combinations should be preferred. Moreover, the
# distributions are assumed to be for the positive charge of the a
# given hadronic species. The negative distributions are derived
# using charge conjugation (q->qbar).

flavour maps: 
  - {hadron: KA,
     map: [0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0, # - sb
           0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0, # - g
           0,  0,  0,  0,  1,  0,  0,  0,  0,  1,  0,  0,  0, # - s = ubar
           0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0, # - u
           0,  0,  0,  0,  0,  1,  0,  1,  0,  0,  0,  0,  0, # - d+
           0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0, # - c+
           0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0  # - b+
          #tb  bb  cb  sb  ub  db  g   d   u   s   c   b   t
            ],
     SetName: MAPFF30_KA_NNLO}
  - {hadron: PI,
     map: [0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0, # - db
           0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0, # - g
           0,  0,  0,  0,  1,  0,  0,  1,  0,  0,  0,  0,  0, # - d = ubar
           0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0, # - u
           0,  0,  0,  1,  0,  0,  0,  0,  0,  1,  0,  0,  0, # - s+
           0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0, # - c+
           0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0
          #tb  bb  cb  sb  ub  db  g   d   u   s   c   b   t
            ],
     SetName: MAPFF30_PI_NNLO}
#  combine: true
#################################
# Data
#################################
Data:
# Seed used for the replica generation and the splitting between
# training and validation (do not use a too large number here in
# order to avoid correlation in the replica generation).
seed: 0

# Datasets to be included in the fit along with specific cuts and
# traning fraction. Each single dataset can implement an arbitrary
# number of cuts determined by the name of the appropriate function
# (e.g. zcut) and the allowed range.
sets:
  - {name: "HERMES $K^-$ deuteron",       file: "HERMES_KA_MINUS_DEUTERON.yaml",         cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "HERMES $K^-$ proton",         file: "HERMES_KA_MINUS_PROTON.yaml",           cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "HERMES $K^+$ deuteron",       file: "HERMES_KA_PLUS_DEUTERON.yaml",          cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "HERMES $K^+$ proton",         file: "HERMES_KA_PLUS_PROTON.yaml",            cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "COMPASS $K^-$",               file: "COMPASS_KA_MINUS_DECORR.yaml",          cuts: [{name: Qcut, min: 2.0, max: 50}], training fraction: 0.5}
  - {name: "COMPASS $K^+$",               file: "COMPASS_KA_PLUS_DECORR.yaml",           cuts: [{name: Qcut, min: 2.0, max: 50}], training fraction: 0.5}
  - {name: "BELLE $K^\\pm$",              file: "BELLE_KA_PLUS_MINUS.yaml",              cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "BABAR conventional $K^\\pm$", file: "BABAR_KA_PLUS_MINUS_CONVENTIONAL.yaml", cuts: [{name: zcut, min: 0.2, max: 0.9}],   training fraction: 0.5}
  - {name: "TASSO 12 GeV $K^\\pm$",       file: "TASSO_12_KA_PLUS_MINUS.yaml",           cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 14 GeV $K^\\pm$",       file: "TASSO_14_KA_PLUS_MINUS.yaml",           cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 22 GeV $K^\\pm$",       file: "TASSO_22_KA_PLUS_MINUS.yaml",           cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TPC $K^\\pm$",                file: "TPC_KA_PLUS_MINUS.yaml",                cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 30 GeV $K^\\pm$",       file: "TASSO_30_KA_PLUS_MINUS.yaml",           cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 34 GeV $K^\\pm$",       file: "TASSO_34_KA_PLUS_MINUS.yaml",           cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 44 GeV $K^\\pm$",       file: "TASSO_44_KA_PLUS_MINUS.yaml",           cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TOPAZ $K^\\pm$",              file: "TOPAZ_KA_PLUS_MINUS.yaml",              cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "ALEPH $K^\\pm$",              file: "ALEPH_KA_PLUS_MINUS.yaml",              cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "DELPHI total $K^\\pm$",       file: "DELPHI_KA_PLUS_MINUS.yaml",             cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "DELPHI $uds$ $K^\\pm$",       file: "DELPHI_KA_PLUS_MINUS_UDS.yaml",         cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "DELPHI bottom $K^\\pm$",      file: "DELPHI_KA_PLUS_MINUS_B.yaml",           cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "OPAL $K^\\pm$",               file: "OPAL_KA_PLUS_MINUS.yaml",               cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "SLD total $K^\\pm$",          file: "SLD_KA_PLUS_MINUS.yaml",                cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "SLD $uds$ $K^\\pm$",          file: "SLD_KA_PLUS_MINUS_UDS.yaml",            cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}
  - {name: "SLD bottom $K^\\pm$",         file: "SLD_KA_PLUS_MINUS_B.yaml",              cuts: [{name: zcut, min: 0.02, max: 0.9}],  training fraction: 0.5}

  - {name: "HERMES $\\pi^-$ deuteron",  file: "HERMES_PI_MINUS_DEUTERON.yaml",   cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "HERMES $\\pi^-$ proton",    file: "HERMES_PI_MINUS_PROTON.yaml",     cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "HERMES $\\pi^+$ deuteron",  file: "HERMES_PI_PLUS_DEUTERON.yaml",    cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "HERMES $\\pi^+$ proton",    file: "HERMES_PI_PLUS_PROTON.yaml",      cuts: [{name: Qcut, min: 2.0, max: 50}, {name: zcut, min: 0.2, max: 0.8}], training fraction: 0.5}
  - {name: "COMPASS $\\pi^-$",          file: "COMPASS_PI_MINUS_DECORR.yaml",    cuts: [{name: Qcut, min: 2.0, max: 50}],    training fraction: 0.5}
  - {name: "COMPASS $\\pi^+$",          file: "COMPASS_PI_PLUS_DECORR.yaml",     cuts: [{name: Qcut, min: 2.0, max: 50}],    training fraction: 0.5}
  - {name: "BELLE $\\pi^\\pm$",         file: "BELLE_PI_PLUS_MINUS.yaml",        cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "BABAR prompt $\\pi^\\pm$",  file: "BABAR_PI_PLUS_MINUS_PROMPT.yaml", cuts: [{name: zcut, min: 0.076, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 12 GeV $\\pi^\\pm$",  file: "TASSO_12_PI_PLUS_MINUS.yaml",     cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 14 GeV $\\pi^\\pm$",  file: "TASSO_14_PI_PLUS_MINUS.yaml",     cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 22 GeV $\\pi^\\pm$",  file: "TASSO_22_PI_PLUS_MINUS.yaml",     cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TPC $\\pi^\\pm$",           file: "TPC_PI_PLUS_MINUS.yaml",          cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 30 GeV $\\pi^\\pm$",  file: "TASSO_30_PI_PLUS_MINUS.yaml",     cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 34 GeV $\\pi^\\pm$",  file: "TASSO_34_PI_PLUS_MINUS.yaml",     cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TASSO 44 GeV $\\pi^\\pm$",  file: "TASSO_44_PI_PLUS_MINUS.yaml",     cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "TOPAZ $\\pi^\\pm$",         file: "TOPAZ_PI_PLUS_MINUS.yaml",        cuts: [{name: zcut, min: 0.075, max: 0.9}], training fraction: 0.5}
  - {name: "ALEPH $\\pi^\\pm$",         file: "ALEPH_PI_PLUS_MINUS.yaml",        cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "DELPHI total $\\pi^\\pm$",  file: "DELPHI_PI_PLUS_MINUS.yaml",       cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "DELPHI $uds$ $\\pi^\\pm$",  file: "DELPHI_PI_PLUS_MINUS_UDS.yaml",   cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "DELPHI bottom $\\pi^\\pm$", file: "DELPHI_PI_PLUS_MINUS_B.yaml",     cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "OPAL $\\pi^\\pm$",          file: "OPAL_PI_PLUS_MINUS.yaml",         cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "SLD total $\\pi^\\pm$",     file: "SLD_PI_PLUS_MINUS.yaml",          cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "SLD $uds$ $\\pi^\\pm$",     file: "SLD_PI_PLUS_MINUS_UDS.yaml",      cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}
  - {name: "SLD bottom $\\pi^\\pm$",    file: "SLD_PI_PLUS_MINUS_B.yaml",        cuts: [{name: zcut, min: 0.02,  max: 0.9}], training fraction: 0.5}

[Copilot]

Copilot reviewed 12 out of 12 changed files in this pull request and generated 2 comments.

[Copilot]

Pull Request Overview

This PR extends the framework to support simultaneous fits of multiple hadrons by generalizing the neural network and LHAPDF interfaces. Key changes include:

• Updates to NNADparameterisation to handle multiple hadrons (via mapping output splits, rotations, and derivatives).
• Modifications to LHAPDFparameterisation and associated chi² routines to retrieve the appropriate distributions based on the hadron type.
• Changes in test and run utilities to remove redundant command‐line options and support hadron‐specific configurations.

Reviewed Changes

Copilot reviewed 15 out of 16 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
tests/Reweighting.cc	Added a warning message regarding unsupported multiple hadrons and adjusted LHAPDF initialization.
tests/PredsComparison.cc, ClosureTests.cc	Updated to call distribution functions with the hadron argument.
src/parameterisations/NNADparameterisation.cc	Refactored to split outputs and rotations by hadron; now uses maps for derivative sets and hadron output sizes.
src/parameterisations/LHAPDFparameterisation.cc	Adjusted constructors and distribution function definitions to use hadron-dependent maps.
run/* (Predictions.cc, Optimize.cc, LHAPDFGrid.cc, ComputeChi2s.cc)	Updated run utilities to work with the new multiple-hadron configuration and simplified usage.
inc/MontBlanc/*.h	Header updates reflecting new function signatures and internal data structures for multi-hadron support.

Files not reviewed (1)

• src/CMakeLists.txt: Language not supported

Comments suppressed due to low confidence (1)

tests/Reweighting.cc:38

• The warning message indicates 'undefined behaviour' when multiple hadrons are provided; consider rephrasing to clearly advise users on the limitations or required input format.

std::cout << "WARNING: this module does not handle multiple hadrons. If the more than one hadrons are "

[Copilot]

Indexing the unordered_map 'SetNamesMap' with the key 0 is likely incorrect since it is keyed by string; consider iterating over the map or using a valid string key for file path construction.