| Version 3 (modified by , 7 years ago) (diff) |
|---|
Page BUILD Built
GENIE NUISCOMP
PAGE BEING BUILT
The NUISANCE minimizer application can be used to tune MC models by directly comparing with published scattering data and using ROOT's minimizer libraries to find a best fit parameter set.
Wiki content
- Page BUILD Built
- GENIE NUISCOMP
- PAGE BEING BUILT
- Running NUISMIN
- Running a GENIE Minimization
Running NUISMIN
Author: Patrick Stowell
Date: June 2017
Versions: NUISANCE v2r0, GENIE 2.12.6
The following example details how to run NUISANCE and tune a simple model to MiniBooNE_CCQE_XSec_1DQ2_nu data, with additional examples on how to include penalty terms and perform fake data fits.
Each generator requires very slightly different ways to handle NUISANCE, therefore multiple versions of this tutorial have been provided. Please use the following links to choose what generator you would like to use.
Running a GENIE Minimization
In this example we will take pregenerated GENIE events as our starting model and run a tuning to find a best fit value for some of the CCQE Reweight dials in GENIE ReWeight?.
Generating GENIE Events
If we want to see how a given GENIE model behaves, first we need to generate events. This can be done using the standard gevgen application, using the appropriate target and flux for a given data sample.
If all you want to do is check your NUSIANCE is built correctly, you can skip this step by downloading MC files from our online storage area by following the steps found here:LinkToNUISANCEMCFiles
Running gevgen
The standard gevgen application options can be ran using
$ gevgen -h
Syntax:
gevgen [-h]
[-r run#]
-n nev
-e energy (or energy range)
-p neutrino_pdg
-t target_pdg
[-f flux_description]
[-o outfile_name]
[-w]
[--seed random_number_seed]
[--cross-sections xml_file]
[--event-generator-list list_name]
[--message-thresholds xml_file]
[--unphysical-event-mask mask]
[--event-record-print-level level]
[--mc-job-status-refresh-rate rate]
[--cache-file root_file]
We need to provide a flux and target list to GENIE when running gevgen.
We want to run comparisons to MiniBooNE muon-neutrino scattering data on a mineral oil target, therefore to generate events we pass it the default GENIE splines, the MiniBooNE flux in root format, and the required target and beam peg settings
$ source $GENIE_DIR/environment_setup.sh
$ export GXMLPATH=${GENIE_DIR}/genie_xsec/v2_12_0/NULL/DefaultPlusMECWithNC/data/
$ gevgen -n 2500000 \
-t 1000060120[0.85714],1000010010[0.14286] \
-p 14 --cross-sections $GXMLPATH/gxspl-FNALsmall.xml \
--event-generator-list Default -f MiniBooNE_numu_flux.root,numu_mb \
-e 0,10 -o gntp.2063030.ghep.root
For more information on how to generate events in GENIE please see: https://arxiv.org/abs/1510.05494
Once our events have been generated we can check that they have finished correctly by opening the output event file and checking it has a ‘gtree’ object and the checking flux spectrum file contains the correct histogram.
$ root gntp.2063030.ghep.root root [0] Attaching file gntp.2063030.ghep.root as _file0... Warning in <TClass::TClass>: no dictionary for class genie::NtpMCEventRecord is available Warning in <TClass::TClass>: no dictionary for class genie::NtpMCRecordI is available Warning in <TClass::TClass>: no dictionary for class genie::NtpMCRecHeader is available Warning in <TClass::TClass>: no dictionary for class genie::EventRecord is available Warning in <TClass::TClass>: no dictionary for class genie::GHepRecord is available Warning in <TClass::TClass>: no dictionary for class genie::Interaction is available Warning in <TClass::TClass>: no dictionary for class genie::InitialState is available Warning in <TClass::TClass>: no dictionary for class genie::Target is available Warning in <TClass::TClass>: no dictionary for class genie::ProcessInfo is available Warning in <TClass::TClass>: no dictionary for class genie::Kinematics is available Warning in <TClass::TClass>: no dictionary for class genie::XclsTag is available Warning in <TClass::TClass>: no dictionary for class genie::KPhaseSpace is available Warning in <TClass::TClass>: no dictionary for class genie::GHepParticle is available Warning in <TClass::TClass>: no dictionary for class pair<genie::EKineVar,double> is available root [1] _file0->ls(); TFile** gntp.2063030.ghep.root TFile* gntp.2063030.ghep.root KEY: genie::NtpMCTreeHeader header;1 GENIE output tree header KEY: TFolder gconfig;1 GENIE configs KEY: TFolder genv;1 GENIE user environment KEY: TTree gtree;1 GENIE MC Truth TTree, Format: [NtpMCEventRecord]
$ root input-flux.root root [0] Attaching file input-flux.root as _file0... root [1] _file0->ls(); TFile** input-flux.root TFile* input-flux.root KEY: TH1D spectrum;1 neutrino_flux
Now that the event samples have been generated correctly, we need to prepare them for use in NUISANCE.
Running PrepareGENIE
The standard gevgen application doesn’t save the total event rate predictions into the event file itself. NUISANCE needs these to correctly normalise predictions so before we can use these new events we need to prepare them.
The PrepareGENIE application is built when NUISANCE is built with GENIE support should be available after the NUISANCE environmental setup script is ran.
$ PrepareGENIE -h
PrepareGENIEEvents NUISANCE app.
Takes GHep Outputs and prepares events for NUISANCE.
PrepareGENIEEvents [-h,-help,--h,--help]
[-i inputfile1.root,inputfile2.root,inputfile3.root,...]
[-f flux_root_file.root,flux_hist_name]
[-t target1[frac1],target2[frac2],...]
Prepare Mode [Default] : Takes a single GHep file, reconstructs the original GENIE splines, and creates a duplicate file that
also contains the flux, event rate, and xsec predictions that NUISANCE needs.
Following options are required for Prepare Mode:
[ -i inputfile.root ] : Reads in a single GHep input file that needs the xsec calculation ran on it.
[ -f flux_file.root,hist_name ] : Path to root file containing the flux histogram the GHep records were generated with. A
simple method is to point this to the flux histogram genie generatrs '-f /path/to/events/input-flux.root,spectrum'.
[ -t target ] : Target that GHepRecords were generated with. Comma seperated list. E.g. for CH2
target=1000060120,1000010010,1000010010
The PrepareGENIE application, when ran, loops over all the events, extracts the cross-section as a function of energy for each discrete interaction mode and uses this information to reconstruct the cross-section splines for each target that were used to generate events.
These splines are then multiplied by specified flux and added according to the target definition provided to produce total flux and event rate predictions as a function of energy for the sample and saves them into the events file.
We want to prepare our MiniBooNE events so we pass in the event files, the input flux, and the CH2 target definition.
$ PrepareGENIE -i gntp.2063030.ghep.root -f input-flux.root,spectrum -t 1000060120,1000010010,1000010010
Now when we open our event file again, we should see the flux and event rate histograms are now saved into the file ready for NUISANCE to read them.
KEY: genie::NtpMCTreeHeader header;1 GENIE output tree header KEY: TFolder gconfig;1 GENIE configs KEY: TFolder genv;1 GENIE user environment KEY: TTree gtree;1 GENIE MC Truth TTree, Format: [NtpMCEventRecord] KEY: TDirectoryFile IndividualGENIESplines;1 IndividualGENIESplines KEY: TDirectoryFile TargetGENIESplines;1 TargetGENIESplines KEY: TH1F nuisance_xsec;1 KEY: TH1F nuisance_events;1 KEY: TH1F nuisance_flux;1
Choosing our samples
Now that we have an event sample we can load them load them into NUISANCE so that it can use them to form a joint likelihood by specifying them at run time.
Writing a card file
To specify samples we need to write a NUISANCE card file that lists all comparisons that should be made and the event files that should be used for each one.
We want to produce comparisons to MiniBooNE CCQE data, so first we should search the NUISANCE sample list.
The ‘nuissamples’ script is provided for easy access of the sample list. Running it without any arguments will return a full sample list of available data comparisons. Providing an additional argument will return only samples containing the provided substring.
We can list the MIniBooNE samples using
[stowell@hepgw1 ~]$ nuissamples MiniBooNE MiniBooNE_CCQE_XSec_1DQ2_nu MiniBooNE_CCQELike_XSec_1DQ2_nu MiniBooNE_CCQE_XSec_1DQ2_antinu MiniBooNE_CCQELike_XSec_1DQ2_antinu MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu MiniBooNE_CCQE_XSec_2DTcos_nu MiniBooNE_CCQELike_XSec_2DTcos_nu MiniBooNE_CCQE_XSec_2DTcos_antinu MiniBooNE_CCQELike_XSec_2DTcos_antinu MiniBooNE_CC1pip_XSec_1DEnu_nu MiniBooNE_CC1pip_XSec_1DQ2_nu MiniBooNE_CC1pip_XSec_1DTpi_nu MiniBooNE_CC1pip_XSec_1DTu_nu MiniBooNE_CC1pip_XSec_2DQ2Enu_nu MiniBooNE_CC1pip_XSec_2DTpiCospi_nu MiniBooNE_CC1pip_XSec_2DTpiEnu_nu MiniBooNE_CC1pip_XSec_2DTuCosmu_nu MiniBooNE_CC1pip_XSec_2DTuEnu_nu MiniBooNE_CC1pi0_XSec_1DEnu_nu MiniBooNE_CC1pi0_XSec_1DQ2_nu MiniBooNE_CC1pi0_XSec_1DTu_nu MiniBooNE_CC1pi0_XSec_1Dcosmu_nu MiniBooNE_CC1pi0_XSec_1Dcospi0_nu MiniBooNE_CC1pi0_XSec_1Dppi0_nu MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc MiniBooNE_NC1pi0_XSec_1Dcospi0_nu MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc MiniBooNE_NC1pi0_XSec_1Dppi0_antinu MiniBooNE_NC1pi0_XSec_1Dppi0_rhc MiniBooNE_NC1pi0_XSec_1Dppi0_nu MiniBooNE_NC1pi0_XSec_1Dppi0_fhc MiniBooNE_NCEL_XSec_Treco_nu
We only care about CC1pip data therefore the following samples are of interest
[stowell@hepgw1 ~]$ nuissamples MiniBooNE_CCQE_ MiniBooNE_CCQE_XSec_1DQ2_nu MiniBooNE_CCQE_XSec_1DQ2_antinu MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu MiniBooNE_CCQE_XSec_2DTcos_nu MiniBooNE_CCQE_XSec_2DTcos_antinu
In this example we will compare to the CCQE 1DQ2 distributions, but we could provide a large number of the samples seen in the lists if we wanted. When we specify multiple samples nuismin will automatically calculate the likelihood for each one and add them uncorrelated to form a joint likelihood.
We write our card file with these two datasets using the following sample object format:
sample NAME_OF_SAMPLE INPUT_TYPE:FILE_INPUT
So if our GENIE file is called 'MiniBooNE_FHC_numu_2.5M.root', our card file will be :
genie_tutorial.card
sample MiniBooNE_CCQE_XSec_1DQ2_nu GENIE:MiniBooNE_FHC_numu_2.5M.root sample MiniBooNE_CCQE_XSec_1DQ2_antinu GENIE:MiniBooNE_FHC_numu_2.5M.root
If you want to check this is valid you could run this card file with 'nuiscomp' as shown in the tutorial here HowToUseNUISCOMP-GENIE and check the nominal prediction looks okay.
Reweighting GENIE Comparisons
The comparisons application also allows different GENIE reweight parameters to be provided in our card file.
The format for fixed dial values is:
genie_parameter NAME VALUE STATE
or for free dial variations
genie_parameter NAME VALUE LOW HIGH STEP STATE
- NAME : Specifies the name of the reweight dial. These can be found in '$GENIE/src/ReWeight/GSyst.cxx'
- VALUE: Value of the reweight dial in units of 1-sigma variations (1-sigma defined by GENIE)
- LOW : Lower limit that the dial is allowed to float to.
- HIGH : Upper limit that the dial is allowed to float to.
- STEP : Step size passed to Minuit.
- STATE: State of this parameter dial, for most cases this should be left as 'FIX'.
For further examples on how to include these structures in card files please see Card File Examples.
In our example we shall generate another set of comparisons this time with two parameters shifted.
- Charged Current Quasielastic Axial Mass : Axial mass parameter used in the quasielastic form factor
- Charged Current Quasielastic Normalisation : Total normalisation of CCQE events.
First we look for the possible name in GENIE reweight:
$GENIE/src/ReWeight/GSyst.h
class GSyst {
public:
//......................................................................................
static string AsString(GSyst_t syst)
{
switch(syst) {
case ( kXSecTwkDial_MaNCEL ) : return "MaNCEL"; break;
...
case ( kXSecTwkDial_NormCCRES ) : return "NormCCRES"; break;
case ( kXSecTwkDial_MaCCRESshape ) : return "MaCCRESshape"; break;
case ( kXSecTwkDial_MvCCRESshape ) : return "MvCCRESshape"; break;
case ( kXSecTwkDial_MaCCRES ) : return "MaCCRES"; break;
case ( kXSecTwkDial_MvCCRES ) : return "MvCCRES"; break;
We can see the dials we are interested in this list, specified by the strings: 'MaCCQE' and 'NormCCQE' respectively.
Now we want to change these dials to +1 for MaCCES and -1 sigma for NormCCRES. Therefore we add the following lines to our card file:
genie_parameter MaCCRES +1.0 FIX genie_parameter NormCCRES -1.0 FIX
so that it now looks like the following
genie_tutorial.card
genie_parameter MaCCRES +1.0 FIX genie_parameter NormCCRES -1.0 FIX sample MiniBooNE_CC1pip_XSec_1DTpi_nu GENIE: gntp.2063030.ghep.root sample MiniBooNE_CC1pip_XSec_1DTu_nu GENIE: gntp.2063030.ghep.root
Setting up our reweight dials
Now that we have our samples, the next thing we want to do is choose some dials to vary. If we look inside
