= How to generate events with NEUT and feed them to NUISANCE =

This page outlines how to use NEUT to create a sample of interactions on your favourite target. We also outline how to feed this generated sample through the NUISANCE framework and get good results.

{{{
#!div class="important"
Contact: [mailto:c.wret14@imperial.ac.uk Clarence Wret]

You can also try: [mailto:p.stowell@sheffield.ac.uk Patrick Stowell]
[mailto:callum.wilkinson@lhep.unibe.ch Callum Wilkinson], [mailto:luke.pickering08@imperial.ac.uk Luke Pickering],
or [mailto:hayato@icrr.u-tokyo.ac.jp Hayato-san].
}}}


[[PageOutline(2-3,NEUT howto,inline,numbered)]]


[=#point_compile]
== Compiling NEUT and generating events ==

This guide assumes you've got your hand on a copy of NEUT and CERNLIB. To request a copy, please kindly inquire to Hayato-san at `hayato@icrr.u-tokyo.ac.jp`.

The main neutrino-nucleus interaction generator lives in `$NEUT_ROOT/src/neutsmpl/neutroot2`, which is the compiled program with source `$NEUT_ROOT/src/neutsmpl/neutroot.F`. 

=== Dependencies ===
NEUT has two dependencies: `ROOT` and `CERNLIB`. 

You need to have `$ROOTSYS`, `$CERN` and `$CERN_LEVEL` environment variables set-up to compile NEUT. These should all be specified by the user in `$NEUT_ROOT/src/neutsmpl/EnvMakeneutsmpl.csh`, or available to the system by

{{{
$ export ROOTSYS=/LOCATION/OF/ROOT
$ export CERNLIB=/LOCATION/OF/CERNLIB
$ export CERN_LEVEL=/CERNLIB/VERSION
}}}


`$ROOTSYS` should be set to your ROOT installation directory(e.g. sourcing `YOUR_DIR/root/bin/thisroot.sh`). 

`$CERNLIB` is your CERNLIB installation folder. 

`$CERN_LEVEL` is the folder inside `$CERNLIB` which you've installed: in my case (2017) the folder is called `2005`, so I set `export CERN_LEVEL=2005`.

=== Compiling ===
Once the environment variables are set, compilation is done by running 
{{{
$ ./Makeneutsmpl
}}}

in `$NEUT_ROOT/src/neutsmpl/`. This sources the above `EnvMakeneutsmpl.csh`, runs through a host of makefiles, sets up symlinks and builds the project.

=== Running the executable ===
You can now run the main executable (`neutroot2`) by:

{{{
$ ./neutroot2 NEUT_CARD_FILE.card OUTPUT_ROOT_FILE.root
}}}
 
where `NEUT_CARD_FILE.card` is a card file compliant with your chosen NEUT version, and `OUTPUT_ROOT_FILE.root` is the ROOT output file where the events are to be saved. 

This output file contains all the event information and is the input for NUISANCE analyses.

=== Card-file examples ===
Examples of NEUT card files can be found in the `$NEUT_ROOT/src/neutsmpl/Cards` directory.

Alternatively, have a look at what team NUISANCE used for 2016 analyses at 
[wiki:Winter2016] and at our [https://imperialcollegelondon.app.box.com/files/0/f/13226969736/neut_cards box].


== Structure of card files ==

NEUT has a load of defaults set on run-time if the user does not supply them, which means you're running with default NEUT.

Selecting different models is done in the cardfile, and all the inputs that a user can specify is found in `src/neutcore/necard.F`.

The user can specify how many events to generate, the input flux in a `TH1D`, the incoming neutrino type, the flux file, number of bound protons and neutrons and free protons, what interaction modes, and so on.

Here's a list for generating all modes on T2K H,,2,,O. The "C" at the start of each line is just a comment and will not be parsed by NEUT.


{{{
C Number of events ; EVCT-NEVT
EVCT-NEVT  100000
C
C Particle Code ; EVCT-IDPT
EVCT-IDPT 14
EVCT-MPOS 2
EVCT-RAD  100.
EVCT-MDIR 1
EVCT-DIR 0. 0. 1.
C
EVCT-MPV 3
EVCT-FILENM  'nuisance_fluxes/t2kflux_2016_plus250kA.root'
EVCT-HISTNM  'enu_nd280_numu'
EVCT-INMEV 0
C
C **** TARGET INFORMATION ****
NEUT-NUMBNDN 8
NEUT-NUMBNDP 8
NEUT-NUMFREP 2
NEUT-NUMATOM 16
C
C **** WHAT INTERACTION MODES DO WE WANT? (0 = ALL) ****
NEUT-MODE   0 
C **** WHAT CCQE MODEL ARE WE USING (2 = RFG WITH BBBA05) ***
NEUT-MDLQE 2
C
C **** RANDOM NUMBER AND IF WE WANT NEUT TO BE QUIET ****
NEUT-RAND 1
NEUT-QUIET 2
}}}

The `EVCT-MPV 3` says the flux is coming from a file, `EVCT-FILENM` specifies the file location of the flux, `EVCT-HISTNM` is the name of the TH1D inside `EVCT-FILENM`, and `EVCT-INMEV` specifies if the input histogram is in MeV or not.


== Compiling NUISANCE against NEUT ==
Once you have a working NEUT installation and generated events (see [#point_compile above], you can link NUISANCE against NEUT.

=== Environment variables ===
To run NUISANCE with NEUT you'll need to provide environment variables `$NEUT_ROOT`, `$CERNLIB` and `$CERN_LEVEL`.

When building, set the environment variables and specify `-DUSE_NEUT=1` when generating the makefiles with `cmake`.

=== Compiling ===

{{{
cd YOUR_NUISANCE_PATH
mkdir -pv build
cd build
cmake ../ -DUSE_NEUT=1
source $uname/setup.sh
make
make install
make docs
}}}

For multiple generator support (e.g. NEUT and GENIE enabled), set the environment variables for the other generators, and compile with

{{{
cmake ../ -DUSE_NEUT=1 -DUSE_GENIE=1
}}}

instead.


== NEUT inputs for NUISANCE ==

=== NEUT output required by NUISANCE ===
The output of `neutroot2`, contains a ROOT file containing `nevents` `neutvect` objects which has all the information about each event in it. This is read in by NUISANCE to generate predictions.

`neutroot2` also saves the predicted flux (`flux_numu`) and event rate (`evtrt_numu`) predictions into the output file, which NUISANCE also uses to convert event rates into cross-sections.

=== Using NEUT events with NUISANCE ===
To read NEUT inputs the generator needs to be specified in the NUISANCE input card file when specifying a sample. e.g. for MiniBooNE CCQE

{{{
sample MiniBooNE_CCQE_XSec_1DQ2_nu  NEUT:/path/to/neut/input/file.root
}}}


== Using NEUT !ReWeight with NUISANCE ==
NUISANCE has full support of running NEUT with its reweighting library. 

Users can generate error-bands, perform generator tunings to specific data-sets, and so on all through NUISANCE.

=== Compiling the reweight classes in NEUT ===
To compile NEUT's reweighting engine, make sure you've got the main executables compiled (see [#point_compile above]), set your `$ROOTSYS` environment varaible and do

{{{
$ cd $NEUT_ROOT/src/reweight
$ make clean
$ make all
}}}

to compile the reweighting engine on top of NEUT itself.

=== Running NUISANCE with NEUT !ReWeight ===
{{{
#!div class="important"
We recommend talking to Hayato-san or us (or other NEUT experts) before fitting parameters to make sure that they are not deprecated.
}}}

The syntax for running NEUT with reweighting in NUISANCE is very similar to how it's done in [wiki:HowToGenie GENIE] and [wiki:HowToNuWro NuWro].

The user specifies the parameter(s) and the range which should be varied in the card file as:

{{{
neut_parameter  MaCCQE  0.0  -4.0  4.0  1.0  FREE
}}}

where `MaCCQE` is the parameter to be varied, `0.0` is the starting values, `-4.0` is the minimum value, `4.0` is the maximum value, `1.0` is the step-size to take, and `FREE` means the parameter is free to vary in the minimisation. 

We can also specify `FIX` for the parameters if we want to reweight a systematic parameter to some value other than nominal and keep it there. 

{{{
#!div class="important"
'''The parameters:'''

The whole list of tunable parameters can be found in `$NEUT_ROOT/src/reweight/NSyst.h`. The user specifies the `return`ed `std::string`s (e.g. `MaCCQE` and `CA5RES`) in the card file.

The starting, minimum, maximum and step values are all in units of 1-sigma uncertainties. So `-4.0` for the minimum means -4*1-sigma ranges, and `0.0` for the starting value means the nominal value nominal.

The hard-coded 1-sigma uncertainty defaults for NEUT can be found in `$NEUT_ROOT/src/reweight/NSystUncertainity.cc` and the `NSystUncertainty::SetDefaults(void)` function. These are in the physical units (e.g. MaCCQE has +/- 0.165 GeV uncertainty in NEUT for 1 sigma)

The hard-coded nominal defaults for NEUT can be found in `$NEUT_ROOT/src/reweight/NFortFns.cc`. These are also the physical units of the parameter (e.g. MaCCQE default is 1.21 GeV in NEUT).

The actual physical parameter value is found by 
{{{
#!span class="code"
`par_value = nominal * (1 + par_var * 1-sigma)`
}}}

The output of the fit is in the units of `par_var`, not `par_value`. So an output for `MaCCQE` of `0.7` would equate to `1.21*(1+1*0.165) = 1.35` GeV.
}}}


== Joint Flux Inputs ==

`neutroot2` does not natively support combined beams or targets. 

If you want to use measurement classes that require this (e.g. MINERvA CC0pi nue+nuebar) you will have to setup a JOINT input. See [wiki:HowToJointInput this] for more information on how this can be done.