Walk Through Tutorial 2025 : Minimizing PROfit
This tutorial assumes that you are relatively familiar with the v1.0.x style XML configuration files. If not please read up on the basics at XML tutorial before starting to work though this document.
The purpose of this walk though is to show the commands and output for a variety of different combinations and studies currently possible in PROfit.
Note: This walk-through has been developed and tested on v1.0.3 PROfit release v1.0.3. We will update accordingly.
The primary XML we will be using is called PROfit_WalkThough_June2025.xml. It points to two simple SBND files and 2 simple ICARUS files for a simple numu disappearance only analysis. Take a look at the XML and see if you understand the systematics being used. Note the files are simple test files and not real physics for the most part!
The files, filelists and XML can be are found:
- SBND gpvms:
/exp/sbnd/data/users/markross/PROfit_Data/Tutorial2025 - uBooNE gpvms:
/exp/uboone/data/users/markross/PROfit_Data/Tutorial2025
Or ping me (Mark Ross-Lonergan) on slack for a tarball.
The standard PROfit executable command is of the form
PROfit -x my_xml_name.xml -t my_tag SUBCOMMAND
or longform
PROfit --xml my_xml_name.xml --tag my_tag SUBCOMMAND
The XML defines the entire analysis, MC files, channels, subchannels ..etc.. and the tag is just a user defined, human readable tag to keep track of which particular study/analysis that you are working on. The SUBCOMMAND is further instructions as to what to actually perform. The allowed subcommands can be found with PROfit --help along with other useful information, lets look at this now:
PROfit: a PROfessional, PROductive fitting and oscillation framework. Together let's minimize PROfit!
Usage: PROfit [OPTIONS] [SUBCOMMAND]
Options:
-h,--help Print this help message and exit
-x,--xml TEXT REQUIRED Input PROfit XML configuration file.
-v,--verbosity ENUM [3] Verbosity Level [1-4]->[Error,Warning,Info,Debug].
-t,--tag TEXT [PROfit] Analysis Tag used for output identification.
-o,--output TEXT [v1] Additional output filename quantifier
-n,--nthread UINT [1] Number of threads to parallelize over.
-m,--max UINT Max number of events to run over.
-c,--chi2 TEXT [PROchi] Which chi2 function to use. Options are PROchi or PROCNP
-d,--data TEXT Load from a seperate data xml/data file instead of signal injection. Only used with plot subcommand.
-i,--inject FLOAT ... Physics parameters to inject as true signal.
-s,--seed INT [-1] A global seed for PROseed rng. Default to -1 for hardware rng seed.
--inject-systs [TEXT,FLOAT] ...
Systematic shifts to inject. Map of name and shift value in sigmas. Only spline systs are supported right now.
--poisson-throw Do a Poisson stats throw of fake data.
--scale TEXT ... Scale detector POT by a given value.
--syst-list TEXT ... Excludes: --exclude-systs
Override list of systematics to use (note: all systs must be in the xml).
--exclude-systs TEXT ... Excludes: --syst-list
List of systematics to exclude.
--fit-options [TEXT,FLOAT] ...
Parameters for single, detailed global best fit LBFGSB.
--scan-fit-options [TEXT,FLOAT] ...
Parameters for simpier, multiple best fits in PROfile/surface LBFGSB.
-p,--preset TEXT Preset fitting params. Available `fast`, `good` and `overkill` .
-f,--rwfile TEXT File containing histograms for reweighting
-r,--mockrw TEXT ... Vector of reweights to use for mock data
--log TEXT File to save log to. Warning: Will overwrite this file.
--scale-by-width Scale histgrams by 1/(bin width).
--event-by-event Do you want to weight event-by-event?
--statonly Run a stats only surface instead of fitting systematics
--force Force loading binary data even if hash is incorrect (Be Careful!)
--no-xrootd Do not use XRootD, which is enabled by default
--shapeonly Excludes: --rateonly
Run a shape only analysis
--rateonly Excludes: --shapeonly
Run a rate only analysis
Subcommands:
process PROcess the MC and systematics in root files into binary data for future rapid loading.
surface Make a 2D surface scan of two physics parameters, profiling over all others.
profile Make a 1D profiled chi2 for each physics and nuisence parameter.
plot Make plots of CV, or injected point with error bars and covariance.
fc Run Feldman-Cousins for this injected signal
protest Testing ground for rapid quick tests.
Take a look at the subcommands at the end. This are all the individual options/things PROfit can do. I will now go though the SUBCOMMANDS one by one outlining their usage and outputs.
process has no additional options or flags.
process is the subcommand used to process all input MC and data root files as defined in the xml, and serialize into a small lightweight binary files. We don't want to read from root files more than necessary (once!). This is true for all CV montecarlo as well as systematics needed to form all covariance and splines for the given binning as defined in the xml.
Once process has been run once, and binary files produced, they are used instead of directly loading the root files in all future subcommands that share the same tag . This is indeed the point of a --tag, to indicate a single processing of the root files into serialized binary files.
process is unique in that it will run automatically if the binary files are not automatically found. So manually running process is not needed. When the binary files are created a hash of the XML is saved internally. If you then change the XML (e.g binning or MC root files) and attempt to load the incorrect binaries (e.g you keep the same tag) PROfit will complain and in this case you can manually reprocess with process command. Its good practice with PROfie to use --log mylog.log or pipe output to a logfile manually. The command is as follows
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG --log mylog.log process
For this exact test tutorial file, this takes about 40 seconds, but only has to be ran once. the expected output is two binary files
MYTAG_prop.bin MYTAG_syst.bin
If you modify the xml and try to run PROfit again with reprocessing you will get the following error:
3 main || Loading PROpeller and PROsysts from precalc binary input: MyTutorial_prop.bin
3 load || Serialization load of PROpeller from file MyTutorial_prop.bin took 0.0179727 seconds
3 main || Done loading. Config hash (1789188867) and binary loaded PROpeller (4220862266) are here.
1 main || ERROR config hash (1789188867) and binary loaded PROpeller (4220862266) not compatable!
If you get this error you can simply re-run process (or change tag, as this will create a new file also).
Subcommand plot has all global arguments and options above as well as some plot specific ones: PROfit plot --help :
Make plots of CV, or injected point with error bars and covariance.
Usage: PROfit plot [OPTIONS]
Options:
-h,--help Print this help message and exit
--with-splines Include graphs of splines in output.
--area-norm Make area normalized histograms.
First lets run a very simple plot of the CV
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG plot
Outputs are as follows
- MYTAG_v1_PROplot_CV.pdf (CV spectra for main reco variable)
- MYTAG_v1_other_0_PROplot_CV.pdf (CV spectra for all OTHER variables)
- MYTAG_v1_PROplot_ErrorBand.pdf (CV spectra + Error band for main reco variable)
- MYTAG_v1_PROplot_other_0_ErrorBand.pdf (CV spectra + Error band for all OTHER variables)
- MYTAG_v1_PROplot_Covar.pdf (All covariance Variables)
- MYTAG_v1_PROplot_CV.pdf (CV spectra for main reco variable)
- MYTAG_v1_PROplot.root (Root file containing above into)
The following is the output of MYTAG_v1_PROplot_Errorband.pdf
You will note there is "data", which here is just the expected CV (or asimov data) as we have not specified otherwise. Lets change that and inject a signal. We do this using --inject followed by the physics parameters values you want. E.g a 10eV^2 sterile with a sin^2thetamumu of 0.3. The physics parameters are defined in your model (e.g in PROmodel.h)
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG --inject 10. 0.3 plot
Now we get (a) an additional file MYTAG_v1_PROplot_Osc.pdf and also out ErrorBar plot now has "data" corresponding to the fake-data injected point.
You will note you overwrote the last PDF. This is not always wanted behavior. You could change the tag but that would mean remaking the binaries. This is where -o or --output is useful. This allows you to have a secondary tag used for printing out only and allows you to generate many different plots without overwriting.
You can also inject systematic shifts. E.g in this example we shift the RPA_CCQE up by 1 sigma and NormCCMEC down by 2 sigma. Note we use -o RPA1CCMEC2 to not overwrite the last plot.
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o RPA1CCMEC2 --inject 10. 0.3 --inject-systs RPA_CCQE 1.0 NormCCMEC -2.0 plot
This fake data now includes both a systematic shift as well as a injected signal:
We can also add the following for additional tweaks. This applies for both subcommand plot and those to follow:
-
--poisson-throw(Add a random Poisson Throw after injection and syst shifts) -
--scale(e.g--scale SBND 0.2Scale channels that match SBND by 0.2) -
--shapeonly(Shape only analysis) -
--rateonly(Rate only analysis, i.e 1 Bin norm)
Ok, so we can make spectra and inject signals. Lets perform a fit.
Perform a full global fit, as well as profiled results for all physics params and splines
Usage: PROfit profile [OPTIONS]
Options:
-h,--help Print this help message and exit
--syst-only Profile over nuisance parameters only
--mcmc-prefit Use MCMC to sample the systematic priors for the pre-fit error band.
Lets take that last injected signal as an example. We will use Poisson throws at the end, so to make it reproducible lets fix the random number seed with --seed 404. Lets also introduce the idea of multi-threading, as we are now going to be performing a full fit. Most fitting code is already threaded in PROfit, and you can use -n or --nthread INT to speed up your studies. Lets do that now, I will use 16 threads but tailor to your situation:
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o RPA1CCMEC2_rand --seed 404 -n 8 --inject 10. 0.3 --poisson-throw --inject-systs RPA_CCQE 1.0 NormCCMEC -2.0 --log fit.log profile
Take a look at all your outputs. First lets look at (a subcset) the MYTAG_RPA1CCMEC2_rand_PROfile.pdf
In all of these we show 1 parameter at a time (either physics or nuisance) and profile over remaining. The top row shows the physics parameters and delta chi^2 after profiling over all others. First 2 (blue) are over whole range, and next 2 (purple) are zoomed in to approx 1 sigma. For this example we injected a strong signal so we easily see the signal injected. Note the x-axis here is log scale.
The remaining show each nuisance/spline parameters. The dashed red is the prior error (by definition 1 sigma gaussian) and the black shows the postfit profiled delta chi^2. You can see the +1 sigma RPA_CCQE injection has a bit effect, but we fit well to the new injected value. You will see that injected -2.0 sigma NormCCMEC was partly covered but out BF here is closer to -1, but the SBNDNorm has also been shifted lower. This shows they are very correlated, and its easier to achieve the same good fit result by lowering both a little compared to NormCCMEC by 2 sigma.
Often we only care about the new +/- 1 sigma error bars and not the full deltachi^2 surface this can be shown in MYTAG_RPA1CCMEC2_rand_PROfile_1sigma.pdf
Here red star is injected point. Black is best fit point and blue is the initial seed. This gives you a quick glance at how things are being constrained as all nuisance variables are initially defined to have +/-1sigma error by construction pre-fit. Note this is strictly less info as the PROfile plot shown above.
To get an idea of how the whole fit is, we can take a look at MYTAG_RPA1CCMEC2_rand_PROfile_hists.pdf
which shows us a lot of interesting information. It contains the pre-fit error bars and pre-fit CV in blue/gray as well as the post fit CV and error bars in Red. The ratio below is the fake-data to post-fit showing you at a glance how good the resulting fit is.
If you want to get the exact values of the global best fit take a look inside your fit.log file. You should find a table outlining the details like below
278 3 main || ################################################
279 3 main || ########### Global Best Fit Results ############
280 3 main || ################################################
281 3 main || Global Best Fit chi^2: 15.116
282 3 main || at paramters:
283 3 main || #Deltam^{2} : 1.00376
284 3 main || sin^{2}2#theta_{#mu#mu} : -0.48838
285 3 main || GENIEReWeight_SBN_v1_multisigma_CoulombCCQE : -0.023345
286 3 main || GENIEReWeight_SBN_v1_multisigma_FrCEx_N : 0.272588
287 3 main || GENIEReWeight_SBN_v1_multisigma_FrCEx_pi : -0.152608
288 3 main || GENIEReWeight_SBN_v1_multisigma_MaCCRES : -0.483774
289 3 main || GENIEReWeight_SBN_v1_multisigma_MvCCRES : -0.388772
290 3 main || GENIEReWeight_SBN_v1_multisigma_NormCCMEC : -0.695212
291 3 main || GENIEReWeight_SBN_v1_multisigma_RPA_CCQE : 0.861459
292 3 main || GENIEReWeight_SBN_v1_multisigma_ZExpA1CCQE : 0.217941
293 3 main || GENIEReWeight_SBN_v1_multisigma_ZExpA2CCQE : 0.0516219
294 3 main || GENIEReWeight_SBN_v1_multisigma_ZExpA3CCQE : 0.0673006
295 3 main || GENIEReWeight_SBN_v1_multisigma_ZExpA4CCQE : -0.00862379
296 3 main || nu_ICARUS:0.02 : -0.50994
297 3 main || nu_SBND:0.02 : 0.366634
298 3 main || ################################################
Try injecting various signals and systematic yourself and seeing the resulting fit! If you don't inject any signal or spline shift, you will see the expected uncertainty band reduction for your expected sensitivity.
PROfit contains two separate levels of fitting parameters.
- Global fitting: Find the full global best fit. This is only done once and we want to make sure we have the BF correct. Uses better, but slower, configuration
- Local fitting: Using the global as seed, this is a quicker, but not as accurate, set of parameters for profiling/surface plotting.
All of these can be manually configured using
--fit-options PARAM VALUE or --scan-fit-options PARAM VALUE for global and local fitting parameters respectively.
The exact parameters and their meaning can be found Here for LBFGSB and below for internal PROfit variables
-
n_multistart: How many initial random points do we check, sampled by a latin-hypercube to guarantee random but no clumping. -
n_swarm_particles: How many of the above multistart points to release as an individual particle in the particle swarm (taking the best of the resulting chi^2 of above). -
n_swarm_iterations: How many iterations/steps to move each particle -
n_localfit: Take the best resultingn_localfitparticles from above, and perform a full L-BFGS-B fit. -
n_max_local_retries: If the fit fails, how many times to retry. Returns best obtained value after.
If you don't want to both with manually tweaking we do have three presets
-
--preset good(default good values) -
--preset fast(faster but not as accurate) -
--preset overkill(Much better but slower)
The defaults for these are given below, for both the local (scan) and global fits
Building a surface scan is very similar to the profile above in terms of the internal structure, fixing (now two) physics parameters one for X axis and one for Y axis, and at each point profiling over all other nuisance parameters.
Make a 2D surface scan of two physics parameters, profiling over all others.
Usage: PROfit surface [OPTIONS]
Options:
-h,--help Print this help message and exit
-g,--grid INT [40] Set grid size. If one dimension passed, grid assumed to be square, else rectangular
--xvar TEXT [sinsq2thmm] Name of variable to put on x-axis
--yvar TEXT [dmsq] Name of variable to put on x-axis
--xlims [FLOAT,FLOAT] Excludes: --xlo --xhi
Limits for x-axis
--ylims [FLOAT,FLOAT] Excludes: --ylo --yhi
Limits for y-axis
--xlo FLOAT [0.0001] Excludes: --xlims
Lower limit for x-axis
--xhi FLOAT [1] Excludes: --xlims
Upper limit for x-axis
--ylo FLOAT [0.01] Excludes: --ylims
Lower limit for y-axis
--yhi FLOAT [100] Excludes: --ylims
Upper limit for y-axis
--xlabel TEXT X-axis label
--ylabel TEXT Y-axis label
--logx,--linx{false} Specify if x-axis is logarithmic or linear (default log)
--logy,--liny{false} Specify if y-axis is logarithmic or linear (default log)
--brazil-band Run 1000 throws of stats+systs and draw 1 sigma and 2 sigma Brazil bands
--stat-throws Needs: --brazil-band
Only do stat throws for the Brazil band
--single-throw Needs: --brazil-band
Only run a single iteration of the Brazil band
--only-throw Needs: --brazil-band
Only run Brazil band throws and not the nominal surface
--from-many TEXT ... Needs: --brazil-band
Make Brazil band from many provided throws
For the most part these are intuitive enough and deal with the ranges and limits for the two physics parameters. Everything you have learnt above about injecting fake data and spline shifts applies also here. Lets first do a simple Asimov sensitivity study, i.e inject no signal and assume we observe the expected CV. The only new option here will be -g 60 which says I want to do a 60x60 grid. You can give rectangular grid values also. Note that here we are doing 3600 fits, each one with many particle swarms and L-BFGS-B minimizations over all nuisance parameters, so this can be slow. However, with multi-threading and a reasonable gridsize/limits you can iterate to a quick but sufficient processing time.
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o surfv1 --seed 404 -n 64 --log surf.log surface -g 60
We have a few outputs. Again in the surf.log you will get info on the global best fit in text form. But we also have a simple txt file output of the entire grid scanned MYTAG_surfv1_surface.txt. This will give you the delta chi^2 at each point scanned, as well as all the best-fit nuisance parameters at each point. E.g in this case, for the point (-0.8,1.6) the delta chi was 2.11783 and the 15 nuisence parameters had the following BF points:
xval yval chi2 p0 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p14
-0.8 1.6 2.11783 1.6 -0.8 0.0053304 0.118631 0.0811355 0.206839 0.149104 0.293632 0.0245438 -0.486896 -0.094708 0.115972 0.0472525 0.111157 0.0118355
This exact info is also saved in a root TTree inside MYTAG_surfv1_surf.root.
This file can be used to load the exact surface into your plotting routine of choice. We stress this: root is not great at making pretty contour plots. We intend for PROfit to give you the surface data, not to be the final plotting software. We do however give some low level plots, and some python matplotlib plotting routines as examples. One of these is inside MYTAG_surfv1_surface.pdf and shown below
Lets repeat this where we exclude an important systematic, say RPA.
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o surfv2_norpa --seed 404 --exclude-systs RPA_CCQE -n 64 --log surf.log surface -g 60
and one where we ONLY use RPA and expskin_Flux for example
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o surfv3_onlyrpa --seed 404 --syst-list expskin_Flux RPA_CCQE -n 64 --log surf.log surface -g 60
We now have three separate contour files. MYTAG_surfv1_surface.txt, MYTAG_surfv2_norpa_surface.txt and MYTAG_surfv3_onlyrpa_surface.txt. Here is a quick matplotlib example of making these into pretty contour plots! You will only need pip install numpy pandas matplotlib scipy.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator
from matplotlib.lines import Line2D
fig, ax = plt.subplots(figsize=(7, 7), dpi=150)
ax.set(xscale='log', yscale='log')
ax.set_xlabel(r'$\sin^2 2 \theta_{\mu \mu}$', fontsize=20)
ax.set_ylabel(r'$\Delta m_{41}^2$ [eV$^2$]', fontsize=20)
# === Define your inputs ===
file_list = [
"MYTAG_surfv1_surface.txt",
"MYTAG_surfv2_norpa_surface.txt",
"MYTAG_surfv3_onlyrpa_surface.txt"
]
labels = [
r"$90\%$ CL (SBN+ICARUS)",
r"No RPA CCQE",
r"Only RPA_CCQE and expskin_Flux"
]
colors = ['red','blue','lightgreen']
legend_proxies = []
# === Loop over files ===
for file, label, color in zip(file_list, labels, colors):
# Detect header
with open(file) as f:
for i, line in enumerate(f):
if line.strip().startswith("xval"):
data_start = i + 1
headers = line.strip().split()
break
# Read data
data = np.loadtxt(file, skiprows=data_start)
df = pd.DataFrame(data, columns=headers)
# Extract and prepare grid
x, y, z = df['xval'], df['yval'], df['chi2']
min_chi = z.min()
sig = df.pivot(index='xval', columns='yval', values='chi2').to_numpy()
xi = np.linspace(x.min(), x.max(), len(np.unique(x)))
yi = np.linspace(y.min(), y.max(), len(np.unique(y)))
interp = RegularGridInterpolator((xi, yi), sig - min_chi, method='pchip')
# Interpolation mesh
N = 125
xg = np.linspace(xi.min(), xi.max(), N)
yg = np.linspace(yi.min(), yi.max(), N)
xg_mesh, yg_mesh = np.meshgrid(xg, yg)
zg = interp((xg_mesh, yg_mesh))
# Plot contour at Δχ² = 1.6
level = [1.6]
ax.contour(10**xg, 10**yg, zg, levels=level, colors=[color], linewidths=2)
# Save proxy for legend
legend_proxies.append(Line2D([0], [0], color=color, lw=2))
# === Finalize plot ===
ax.legend(legend_proxies, labels, fontsize=16, loc='lower left').get_frame().set_linewidth(0.0)
ax.set_xlim(1e-2, 1)
ax.set_ylim(1e-2, 100)
plt.tight_layout()
plt.show()The resulting image is as follows.
Note that we do plan on integrating this more into PROfit directly. Especially if you use the python bindings this is trivial.
We can also do a signal injection study, in the same vein, but zoom in on the region of interest. For example a 10 eV^2 sterile with a few options. See if you can tell what these commands do, run some similar and modify the python script to plot successfully.
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o inject1 --inject 10 0.2 --seed 404 -n 16 --log surf1.log surface -g 60 --xlo 0.01 --ylo 1
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o inject2 --inject 10 0.2 --seed 404 --syst-list expskin_Flux RPA_CCQE SBNDNorm MaCCRES -n 16 --log surf2.log surface -g 60 --xlo 0.01 --ylo 1
PROfit -x PROfit_WalkThough_June2025.xml -t MYTAG -o inject3 --inject 10 0.2 --seed 404 --statonly -n 16 --log surf3.log surface -g 60 --xlo 0.01 --ylo 1
This results in something like this. Note we made quite a strong signal, so probably could have run a more zoomed in/denser grid!
