Comparing generators to data with nuiscomp

Written by Clarence on 22 April 2019.

Overview

A broad range of datasets have been implemented into the NUISANCE framework that support direct comparison between different generator Monte Carlo event samples and published data. The full list is updated for every git commit and available here.

nuiscomp is the simplest executable comparing data to generator predictions. The datasets and generator files are specified by providing different sample IDs and pointing to different files through a "card file".

You can also provide NUISANCE with systematic parameters and a variety of options outlined below to do closure tests, fitting systematics, and making error bands.

For the code-inclined, nuiscomp lives in nuisance/app/nuiscomp.cxx and is driven by nuisance/src/Routines/ComparisonRoutines.cxx. You can view it here.

This expands on the documentation provided by running

which is equivalent to

An example nuiscomp run

This example compares GENIE events to MiniBooNE data. The example is near identical for a different generator, you just have to Prepare appropriately and replace GENIE in the cardfile with your generator descriptor, e.g. NEUT, NUWRO, GIBUU.

Make a plain file mbcomp.card and fill it with:

You can now run nuiscomp with this card file by doing

Specifying a measurement

Changing the measurement is as easy as replacing the sample name field in the card file. The sample names can be found by running the nuissamples script in nuisance/scripts/nuissamples, making sure you have exported the NUISANCE environment variable to be your NUISANCE directory.

The output strings can then be used in the sample name field, e.g. sample name="MiniBooNE_CC1pip_XSec_1DQ2_nu" would specify MiniBooNE's CC1π⁺ cross section measurement in one dimension in Q² for neutrino.

The script essentially looks at the string comparison fields in nuisance/src/FCN/SampleList.cxx which is run by nuiscomp. Check it out if you're curious.

Options for a measurement

Samples have a normalisation option, which is applied through norm="scale", where scale is a number. It is only accepted if type="FREE" is specified, e.g.

The options given to type relate to how the sample is treated, either statistically (e.g. 𝜒² calculation without a covariance matrix) or normalisation scaling (e.g. doing a shape-only fit). These come in full effect in nuismin which is driven by a test-statistic which isn't the case for the simpler nuiscomp.

Use type="OPTION" with OPTION set to:

• FIX, FREE, SHAPE: Keep the sample normalisation fixed, free, or shape only (normalises the generator prediction to the data)
• FULL, DIAG: Use the full covariance or diagonal covariance matrix to calculate the test-statistic
• NORM: Add a normalisation penalty if one is specified in the sample, controlled by the private member fNormError, set in each sample's constructor.
  e.g. for MiniBooNE T_μcosθ_μ it is specified that the normalisation error was 10.1%, so is implemented here.

Sample dependent and debugging options

Some lesser used (sample dependent) types to provide to a sample are:

• Q2CORR: Apply a Q²-dependent correction if data supports it (currently just ANL, BNL, BEBC and FNAL bubble chamber CCQE measurements)
• MASK: Use provided masking histogram. This may be useful if you want to use a measurement that should mask out regions you don't want to use, e.g. low Q² in a bubble chamber

Some even lesser used options are related to how the sample is scaled. These options are always set in the sample implementation file and shouldn't have to be touched by a user for other than debugging purposes. They are:

• ENU1D: Treat the measurement as a total cross-section: that is scale toa bin-by-bin flux integrated cross-section rather than flux integrated over the entire range.
• NOWIDTH: No bin-width scaling is applied (so is e.g. cross-section/bin rather than cross-section/MeV)
• RAW: Measurement is raw events (not a cross-section), so does not scale MC to a differential cross-section. Instead it scales the MC integral to be the data integral for every iteration

Additionally, some samples have settings specific to each measurement. These are documented in the implementation. A good example is nuisance/src/ANL/ANL_CC1ppip_XSec_1DEnu_nu.cxx (found here), which has many data releases and corrections: W < 1.4 GeV, W < 1.6 GeV, and no W cut. Furthermore the W < 1.4 GeV and no W cut data has the ability to apply a E_𝝂 dependent correction (Phys. Rev. D 90, 112017 and Eur. Phys. J. C (2016) 76: 474. ).

Specifying a generator

To change the generator you first need to have generated events with said generator and linked NUISANCE to the generator. You can find notes on this here. You also need to have prepared generator files for NUISANCE use.

Once you have NUISANCE linked and prepared files you simply change the input field to match your generator (GENIE, NEUT, NUWRO, GIBUU or NUANCE), and new file location. For example input="NUWRO:PATH_TO_PREPARED_NUWRO_FILE.root".

Changing systematics via event reweighting in nuiscomp

An overview of the systematics present in NUISANCE is given here. Here we show how to include systematics in card file for nuiscomp.

Expanding on mbcomp.card with mbcomp_fit.card:

For nuiscomp we can additionally provide a state for each parameter. The valid state for parameters is just if we're specifying the nominal in terms of 1 sigma, fractional or absolute values.

• type="ABS" specifies it as an absolute value and uses the lookup in ${NUISANCE}/parameters/dial_conversion.card
• type="FRAC" specifies it as a fractional value
• Not specifying uses the default generator behaviour, often meaning relative the 1 sigma uncertainty (e.g. for GENIE these are found in ${GENIE}/src/ReWeight/GSystUncertainity.cxx).

N.B. nuismin expands on these options, as seen here.

Making fake-data via event reweighting

Fake-data can be provided to nuiscomp in two ways:

• Giving the -d parameter to nuiscomp (outlined later)
• Specifying a fakeparameter field

A card file for MiniBooNE CCQE in 1D Q² for neutrino, setting the fake-data to be GENIE with MaCCQE set to 1 sigma above the default, would be

Running nuiscomp

nuiscomp takes a number of command line arguments, one of which is required. It requires either

• -c mycard.card where mycard.card is the provided card file (see above)
• -i 'card file contents'

The -i option follows the same structure as a card file but is instead an argument. For example

would acheive the same effect as providing nuiscomp with the simple mbcomp.card earlier.

The complete list of arguments and their types is:

• -c string: the cardfile to read
• -o string: the output file to create (default: append ".root" to cardfile name)
• -n int: the number of events to run on (default: all)
• -d string: fakedata input, can be either (default: empty)
• "MC": for setting the data to the nominal MC prediction
• "fakedata.root": Direct path to root file containing a MC histogram called the measurement name with "_MC" at the end (the default format for NUISANCE MC histograms)
• -i string: card file-like XML commands (default: empty)
• -q string: config arguments, which overwrite nuisance/parameters/config.xml (default: empty)
• e +/- error level relative that in nuisance/parameters/config.xml
• v +/- verbosity level relative that in nuisance/parameters/config.xml

Running nuiscomp as

Inspecting nuiscomp logs

NUISANCE produces various output as it sets up, using different verbosity logging levels, defined in nuisance/parameters/config.xml. This section explains what lines indicate what.

In this example the home directory is /home/nuisance, built in /home/nuisance/build. The cardfile is called cardfile.xml, containing

Printing verbosity

The first lines set up the logger with verbosity levels, defined in /home/nuisance/parameters/config.xml. Here we've set it to VERBOSITY=4 and ERROR=2.

Reading the cardfile

NUISANCE then continues to process the cardfile and register the settings internally. It recognises our four GENIE parameters, their values and types. It then proceeds to find the sample that we specified, it's type, normalisation and file location

Loading the ReWeighter

The reweighting parameters are specified in the FitWeight class, which interfaces to GENIE and the reweighting. The GENIE parameters are each given a unique enum. The sample normalisation is also registered as a parameter and is given a unique enum.

Building the samples

We then progress to reading the input sample, its file location and loads up the GENIE event handler. We're provided with some information on the input file: the number of entries we're reading, the converted event rate and flux, and finally the Event/Flux which is equivalent to the flux averaged single bin cross-section.

We're then told where the data is being read from (/home/nuisance/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.root) and the histogram (q2qe_cross_section) which it contains. NUISANCE does a number of checks of the data and here notices that for some bins the provided covariance matrix does not agree with the error on the data histogram. In this case the mismatch is 0.1E-40, which isn't too concerning, although definitely is a problem!. Here we contacted the analyser (Dan Ruterbories, Rochester) and this is apparently due to rounding errors in making the covariance matrix. We suggest you do the same if noticing something similar!

The event loop

The main event loop starts once everything has been loaded up successfully. Events are extracted out of the input file, converted to the common format, and event characteristics are calculated (e.g. p_μ and cosθ_μ to construct Q²_QE). Each event passes through the signal definition defined in the sample, and the generator (MC) histograms are filled up.

The output prints it's running the Compare routine, and starts the first iteration.

NUISANCE reports the status of the event loop with 10 prints to screen per reconfigure. So running on 100,000 events will cause a print every 10,000 event. With every print, NUISANCE reports the interaction mode (M) and the weight (W) from the main loop. The interaction modes are documented here.

The bold underscored line with [M, W] = [1, 0.845794] is a CCQE event (M=1). Since we're running with a variation of the z-expansion parameters we'd expect a CCQE event to receive a weight different to 1, which is indeed the case (W=0.845794). This is useful for debugging purposes, and to see the progress during a fit with nuismin.

Finally the number of signal events (11083 out of 100000) is reported, and the time taken for a full reconfigure (12 s). Hence the looping and reweighting through GENIE ReWeight takes about 0.1 ms/event (N.B. without reweighting the time taken was 7 s).

The summary

The summary of the final state of the program is printed at the end of the routine.

We see a summary of the parameters and their values, and a converted value in real units (which uses nuisance/parameters/dial_conversion.card if available). As expected, we have four z-expansion parameters and one sample normalisation parameter. They are all set to the parameter values provided in the card file.

We also see a likelihood for the sample of -2logL=348.248/16, which means the χ2=348.248 and the number of bins N=16.

Inspecting nuiscomp output

Once the code has run nuiscomp saves the histograms to a root file (specified by the user giving the -o option).