Version 33 (modified by 8 years ago) (diff) | ,
---|
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.
Contact: Clarence Wret
You can also try: Patrick Stowell Callum Wilkinson, Luke Pickering, or Hayato-san.
NEUT howto
- Compiling NEUT and generating events
- Structure of card files
- Compiling NUISANCE against NEUT
- NEUT inputs for NUISANCE
- Using NEUT ReWeight with NUISANCE
- Joint Flux Inputs
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
, which lives in $NEUT_ROOT/src/neutsmpl/neutroot2
. Do
$ ./neutroot2 NEUT_CARD_FILE.card OUTPUT_ROOT_FILE.root
where NEUT_CARD_FILE.card
is a card file compliant with your chosen NEUT version
OUTPUT_ROOT_FILE.root
is the ROOT output file where the events are to be saved.
The output file contains all the event information and is the input for NUISANCE analyses.
NEUT 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 Winter2016 and at our 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 card-file, 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 incoming neutrino type, the input flux file, number of bound protons and neutrons and free protons, what interaction modes, and so on.
Here's an example card-file for generating all modes on T2K H2O. 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 above, you can link NUISANCE against NEUT.
Environment variables
To run NUISANCE with NEUT you'll need to provide environment variables $NEUT_ROOT
, $CERN
, $CERN_LEVEL
and $ROOTSYS
.
When building, set the environment variables and specify -DUSE_NEUT=1
when generating the makefiles with cmake
.
e.g. in my case
$ export NEUT_ROOT=/vols/software/NEUT_5.3.6 $ export CERN=/vols/software/CERNLIB $ export CERN_LEVEL=2005 $ source /vols/software/root/bin/thisroot.sh
where the last place could be replaced by
$ export ROOTSYS=/vols/software/root
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 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
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 GENIE and 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.
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
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 this? for more information on how this can be done.