GEANT4 Parameter Tuning Using Professor

The Geant4 toolkit is used extensively in high energy physics to simulate the passage of particles through matter and to predict effects such as detector efficiencies and smearing. Geant4 uses many underlying models to predict particle interaction kinematics, and uncertainty in these models leads to uncertainty in high energy physics measurements. The Geant4 collaboration recently made free parameters in some models accessible through partnership with Geant4 developers. We present a study of the impact of varying parameters in three Geant4 hadronic physics models on agreement with thin target datasets and describe fits to these datasets using the Professor model tuning framework. We find that varying parameters produces substantially better agreement with some datasets, but that more degrees of freedom are required for full agreement. This work is a first step towards a common framework for propagating uncertainties in Geant4 models to high energy physics measurements, and we outline future work required to complete that goal.


Introduction
G 4 [3,5,6] is a toolkit for simulating the passage of particle through matter that is used in a variety of fields. In high energy physics, it is routinely used to construct detailed simulations of particle detectors that are used to estimate quantities such as detector smearing, acceptance and efficiency. In neutrino physics, it is additionally used to simulate neutrino beamlines and to predict the neutrino flux through detectors. G 4 typically accepts an input list of particles (with initial 4-positions and 4-momenta) and a geometry from the user. Each particle is stepped through the geometry, with probabilities of a variety of different kinds of interactions to happen at each step. Once an interaction occurs, the final state particles and their momenta are determined from a model of that process based on theory and tuned to data where possible. Various options for both the interaction cross sections and the final state models are available.
Many of these models, and in particular those describing interactions of hadrons on nuclei, carry significant uncertainties. When G 4-based simulations are used to correct observed data distributions and to extract measurements of fundamental parameters, these models introduce uncertainty in the measurements. These uncertainties have historically been difficult to quantify, with experimenters typically guessing at what aspects of the G 4 simulation have the biggest impact on a measurement, and then attempting to quantify the uncertainty in those effects.
-1 -With the release of G 4 version 10.4, the G 4 collaboration made it possible to access some underlying parameters in several models, and to vary these parameters. This opens the possibility of fitting these parameters to datasets, to extract optimal values of the parameters and uncertainties in those parameters that could be propagated to physics measurements.
One significant hurdle to this goal is the fact that the available G 4 parameters are not reweightable. In order to produce a simulation with a varied set of parameters, one must execute a completely new simulation; it is not currently possible to vary the models by applying weights to events. With simulations of typical experimental setup taking hours, fits of many parameters to many datasets are computationally challenging. Professor is a framework for parameter fitting that has been successfully used to tune non-reweightable parameters of models in hadron collider event generators such as Pythia [21,22].
This article describes a first attempt to apply the Professor framework to G 4 hadronic model parameter tuning. Section 2 describes the models and parameters considered; thin target datasets considered in this study are listed in Section 3; Section 4 compares scans of these parameters to available datasets; Section 5 describes the results of Professor fits. While this represents the first step towards a common framework for assessing G 4 model parameter uncertainty in high energy physics measurement, a broad program of work would be needed to fully achieve that goal. The major components of that program are outlined in Section 6

Available Geant4 Parameters
As of release series 10.4, G 4 provides a configuration interface for the hadronic models Bertini Cascade, PreCompound, and Fritioff (FTF). A configuration interface is also available to control and/or tweak modeling of the electromagnetic (EM) processes, but those are not considered here. The configuration interfaces differ from model to model, and guidance for varying the models is available in Geant4 user documentation [17,18]. The models also include a number of switches and parameters with discrete value options; the data comparisons and fits described here use the default value of the switches and discrete parameters.
In Geant4 release 10.4, FTF offers 16 configurable parameters governing baryon projectiles only. These parameters, their default values in Version 10.4, and the limits set by model developers are listed in Table 1. Several of the parameters are associated with modeling nuclear destruction in baryon-nucleus interactions. In these interactions, the probability of the nucleons to be involved in a reggeon cascade is given by where s i and s j are projections of the radii of i-th and j-th nucleons on the impact parameter plane and the coefficient C nd is Other parameters involve modeling of momentum distributions of the nucleons involved in the cascade, which is described in greater detail in [17]; one key distribution is the average transverse momentum square of an ejected nucleon, which can be parameterized as where y lab is the rapidity of the projectile nucleus in the rest frame of the target nucleus.

Datasets Considered
There is a large body of data that could be used to tune hadronic models in G 4. Upon consultation with the model developers, we chose thin target datasets published by the HARP, IAEA, NA49, NA61, and ITEP771 collaborations, as detailed in table 2.

Sensitivity to Parameters in Modeling Thin Target Data
Some of the available Geant4 parameters have a greater impact on predictions of the thin target datasets listed in Table 3 than others. Because the computing resources required for parameter fits grow with the number of parameters considered, it is useful to understand which parameters create the greatest changes in model predictions, and are therefore most important for inclusion in fits.
To estimate the sensitivity of each parameter to thin target predictions, we produced simulations of several thin-target measurements where all other parameters are fixed to their default values and the parameter in question was varied within the limits specified by the model developers (or between 50% and 200% of the default value if limits are not specified). An example of such a parameter scan is shown in Figure 1 for several quantities measured by NA49 for 158 GeV proton interactions on carbon and for the FTF parameter BARYON_NONDIFF_M_PROJ. Example parameter scans for parameters in the Bertini and PreCompound are shown in Figures 2 and 3, respectively.
For each distribution considered in the parameter scan, we calculate the asymmetry for a particular bin as |y max − y min |/|y max + y min | where y max (y min ) is the maximum value of the predicted distributions in bin in question. This is then averaged over all bins and all distributions to produce a total asymmetry, shown in Table 1. Based on these sensitivities and other details of the parameter scans, the parameters noted with an * in Table 1 were chosen for inclusion in fits.

Parameter Fits using Professor
Generating a model prediction for all of the thin target datasets for a given set of parameters requires multiple CPU hours. Traditional fitting methods wherein many possible parameter values are scanned, simulations are performed, and chi-square differences between models and data are minimized, is not computationally feasible. Instead, the Professor [15] fitting framework was used. Professor parameterizes the prediction for a given bin of an observable distribution using the ndegree polynomial of parameters. The coefficients of the polynomial are analytically evaluated by the singular value decomposition using a simulated data ensemble in which the underlying parameters are thrown randomly within their a priori probability distributions (which are taken to be flat distributions within parameter limits for the Geant4 parameters discussed here). The resulting polynomials can then be used to construct predictions of each dataset given any set of parameters, substantially reducing the time required to produce dataset predictions from hours to a fraction of a second, and allowing traditional chi-square minimization.   Predictions for 158 GeV proton interactions on Carbon in default Geant4 (red) and with the FTF parameter BARYON_NONDIFF_M_PROJ varied within the its allowed limits (green); NA49 measurements [7,13] of the same quantities are also shown (black).

Parameter Fits
Three "global fits" were performed, one for each of the three models considered here. In each case, the parameters of the model in question were fit to all the datasets listed in Table 2 simultaneously. Example fit results are shown in Figures 4-6 for selected datasets. Comparisons of additional datasets with these global fit results are available in Appendix A. The appendix includes results compared to all fitted datasets except those considered in the Bertini fit, where only ITEP 7.5 GeV proton data are shown for the sake of brevity. In each figure, data points from the relevant datasets are compared with G 4 predictions with both the default and best fit parameters. The an error band on the fit result is taken from parameter errors returned by the fit, accounting for correlations. However, because the best fit chi-square per degree of freedom is much more than one, these error bands are not complete and do not include uncertainty associated with the models' inability to produce good agreement with data [20]. A summary of the default and best-fit chi-squares is Figure 2. Predictions for 5 GeV π − interactions on Carbon using default Geant4 (red) and with the Bertini parameter RadiusScale varied within the its allowed limits (green). HARP measurements [7,13] of the same quantities are also shown (black). available in Tables 5-3. Fits were also performed individually for several datasets; the chi-squares and best fit parameters for those fits are also available in Tables 5-3. While improvement in fit quality is obtained by fitting each dataset separately, fit quality is still poor for most datasets, and the best fit parameters vary significantly between these fits.
Geant4 predictions using the best fit parameters do improve agreement with the data in some areas. For example, Bertini agreement with ITEP proton scattering data is much improved for many (but not all) nuclei. In other cases, the fit causes data Monte Carlo agreement to be substantially worse, e.g. in comparisons with NA61 31 GeV pC → π − X data using FTF. In other cases, e.g. in predictions of IAEA 1.5 GeV pC → nX data, the best fit is close to the default prediction.

A Dependence of Fit Results
It is clear from the fits described above that expanded Geant4 parameter space will be required to achieve good agreement with many datasets. One natural expansion is to allow different parameter values for different nuclear targets. To explore this possibility, separate Bertini model parameter fits were performed to IAEA 1.5 GeV proton datasets on several nuclear targets. The best fit parameters versus nuclear mass number (A) is shown in Figure 7. In some cases, there are clear trends versus nuclear mass. The preferred value of the FermiScale parameter decreases as nuclear mass increases, whereas RadiusScale is nearly flat (perhaps because the fitter pushes that parameter to the lower GeV proton interactions on Carbon using default Geant4 (red) and with the PreCompound parameter LevelDensity varied within the its allowed limits (green). IAEA measurements [19] of the same quantities are also shown (black).    limit of the allowed range for that parameter). It appears that better parameter fits could be obtained if Geant4 offered the possibility of setting different parameters for different nuclear targets.

Comparison of Professor Fit Results with G 4 Simulations
As discussed above, the Professor fits do not actually run Geant4 simulations as part of the fitting algorithm. Rather, interpolation between many simulations with randomly chosen parameters is used to estimate the result one would obtain with a true Geant4 simulation.
To study the efficacy of the interpolation, we compared the best-fit predictions returned by Professor with full Geant4 simulations run with the best-fit parameters. In general, Professor does a reasonably good job of predicting the full simulation, but breaks down in some areas of phase space, as shown in Figure 8. We find this to be particularly likely when parameters are pushed to the edges of the allowed limits. Increasing the order of the polynomial fit function used by Professor may improve this behavior.

Future Work
This work illustrates that recently available variable parameters in the Bertini, PreCompound and FTF models in Geant4 provide a powerful mechanism for tuning Geant4 predictions to data. The Figure 6. Results of the global FTF parameter fit, compared to NA49 31 GeV pC → K 0 S X and pC → K − X data. Data points are shown in black; default Geant4 in red and Geant4 with best fit parameters in blue; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. work was motivated by a desire for a framework for propagating Geant4 model uncertainties to measurements, similar to that provided by the GENIE neutrino event generator [8,9]. Many steps remain before that could be a reality, including • Additional datasets and degrees of freedom in fits: While the fits described in Section 5 do improve data/simulation agreement, it is clear that the parameters currently available in Geant4 are not sufficient to bring predictions into agreement with the data. Fits to individual datasets indicate that simple extensions of parameters, e.g. allowing different parameters for different projectiles or target nuclei would improve agreement, but that these simple extensions would need to be supplemented by additional new parameters, and/or significantly expanded ranges of existing parameters. Also, there are clear limitations of the current parameters; FTF provides parameters that alter models only for baryon projectiles, for example. There is also additional thin target data that could be included in additional fits, that may require yet more parameters beyond what is indicated here.
• Inclusion of cross sections All of the parameters discussed here modify the final states simulated by Geant4 given that some type of interaction has occurred. They do not modify the probability that an interaction or particular type of interaction will occur. Any complete assessment of Geant4 uncertainties would need to include uncertainties in the cross sections -11 - themselves, not just final state models. In fact, the interaction cross sections are among the most important sources of uncertainty in some high energy physics measurements, and experiments have developed event weighting mechanisms to assess their impact [1,4].
• Expansion to more G 4 models The models discussed here are those used by the Geant FTFP_BERT physics list to simulate the passage through matter of hadronic particles with energy between 0 and 100 TeV. They were chosen because they are important models for a variety of planned and future high energy physics experiments. But there are a number of interaction types that are not covered by the models or parameters discussed here. Electromagnetic models, for example, although fairly well understood, can lead to uncertainties in modeling off quantities such as shower shape in detectors. Other Geant4 hadronic models that could be included in future fits are INCL++ (Inter-nuclear Cascade), QGS (Quark Gluon String) and BIC (Binary Light Ion Cascade).
• Treatment of correlated errors in datasets All of the fits described above assume that there are no bin-to-bin correlations in the thin-target data. This approach was taken because the experiments reporting thin target data have generally not provided covariance matrices. However, it is likely that some of the uncertainties in these datasets are correlated. For example, the NA49 data included a 3.8% systematic uncertainty which may be correlated across bins. Future efforts at fitting G 4 parameters should consider the possibility of bin-to-bin correlations in the datasets, which may substantially alter the fit results.
• Reweightable parameters A major hurdle to using the parameter framework currently available in Geant4 for propagation of model-related systematics to physics measurements is that the parameters are not reweightable. That is, when a parameter is changed, the user must run an entirely new simulation with the parameter change. Thus, to propagate G 4 systematics to measurements, experiments would have to produce many different varied simulations. Currently, the computational resources required to generate a single simulation make this prohibitively expensive for most high energy physics collaborations. Although modification of models and assessment of systematic uncertainties via event weighting has drawbacks (e.g. weighting can never create events in phase space that was not generated in the first place), it is the only way many experiments can realistically propagate model uncertainties at present. Implementation of a reweighting engine for the Geant4 parameters would therefore dramatically improve the feasibility of using these parameters for error propagation. One method of reweighting that effectively modifies Geant4 models is reweighting of double or triple differential cross sections [1,4]. The advent of high performance computing may make generation of alternative samples more feasible in the future. Figure 8. Results of the best-fit Professor prediction (blue line) from the Bertini Global fit, compared to IAEA 1.5 GeV pPb → nX data (black points) and to a Geant4 simulation run with the best fit parameters extracted from the Professor fit. In general, the professor prediction agrees with a full simulation, but the Professor prediction fails in some areas of phase space.

Conclusion
We have studied the hadronic model parameters recently made available by developers of the Bertini,PreCompound, and FTF models in Geant4. These parameters facilitate variation of the final state content of hadronic interactions in detector and beam simulations. We have varied each of the parameters within the ranges set by developers and compared the resulting simulations to an array of thin target datasets. We have identified the parameters that create the largest variations in predictions, and have tuned those parameters to the data using the Professor fitting framework. Although agreement with data is improved, model predictions are still quite far from data in many areas of phase space. Steps for further work that would develop this parameter infrastructure into a framework for error propagation have also been outlined. The toolkit and fitting framework used for these studies is available for use through correspondence with the authors.                Figure 28. Results of the global FTF parameter fit, compared to NA49 158 GeV pC → π − X, pC → π + X, and pC → pX, pC →pX and pC → nX data. Data points are shown in black; default Geant4 is red and Geant4 with best fit parameters in blue; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. Figure 29. Results of the global FTF parameter fit, compared to NA61 31 GeV pC → K − X, pC → K + X and pC → ΛX data in bins of final state hadron angle. Data points are shown in black; default Geant4 is red and Geant4 with best fit parameters in black; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. Figure 30. Results of the global FTF parameter fit, compared to NA61 31 GeV pC → ΛX, pC → π − X, and pC → π + X data in bins of final state hadron angle. Data points are shown in black; default Geant4 is red and Geant4 with best fit parameters in black; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. Figure 31. Results of the global FTF parameter fit, compared to NA61 31 GeV pC → π + X and pC → pX data in bins of final state hadron angle. Data points are shown in black; default Geant4 is red and Geant4 with best fit parameters in black; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. Figure 32. Results of the global FTF parameter fit, compared to HARP 5 GeV pC → π + X and pC → π − X data in bins of final state hadron angle. Data points are shown in black; default Geant4 is red and the global fit result in blue; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. Figure 33. Results of the global FTF parameter fit, compared to HARP 5 GeV pC → π − X and pPb → π + X data in bins of final state hadron angle. Data points are shown in black; default Geant4 is red and the global fit result in blue; the green band shows uncertainties propagated from parameter uncertainties returned by the fit. Figure 34. Results of the global FTF parameter fit, compared to HARP 5 GeV pPb → π + X and pPb → π − X data in bins of final state hadron angle. Data points are shown in black; default Geant4 is red and the global fit result in blue; the green band shows uncertainties propagated from parameter uncertainties returned by the fit.  -44 -