DarkNews is an event generator for new physics processes at accelerator neutrino experiments that simulates neutrino upscattering to heavy neutral leptons and their subsequent decays to single photons and di-lepton pairs.
DarkNews uses Vegas to generate weighted Monte Carlo samples of scattering and decay processes. Differential observables are implemented using analytical expressions with arbitrary interaction vertices, which are then specified at run-time based on the available models and the user's parameter choices. Processes involving heavy neutrinos N are calculated including contributions from the Standard Model Z and W bosons, as well as from a kinetically-mixed dark photon (Z'), and neutrino-N transition magnetic moments. DarkNews also computes the partial decay widths of heavy neutrinos for models with up to 3 heavy neutrinos.
Experiments as well as models are implemented on a case-by-case basis. The necessary ingredients to simulate upscattering or decay-in-flight rates are the active-neutrino flux, detector material, and geometry.
The full information of the event genration is saved to a pandas dataframes, but the user may also choose to print events to numpy ndarrays, as well as to HEPevt-formatted text files.
If you experience any problems or bugs, either open a new issue or contact [email protected].
Required dependencies:
The following dependencies (if missing) will be automatically installed during the main installation of the package:
- scipy
- pandas 1.0 or above
- pyarrow
- Cython
- vegas 5.1.1 or above
- Particle
- pyhepmc
- pyparsing
- dill
- matplotlib
DarkNews is available on PyPI so from your python3.7+ (virtual environment or otherwise), run
python3 -m pip install DarkNews
or if your pip version is already set to your preferred python version, simply pip install DarkNews
. This should install all dependencies for you.
troubleshooting: If you have any problems, try creating a brand new (conda or pyenv) environment, install the latest version of pip
, then pip install numpy first, and only then try to install pip install DarkNews.
If you experience any issues installing pyhepmc-ng, try installing numpy first, and then install pyhepmc-ng directly from Git using the following command: pip install git+https://[email protected]/scikit-hep/pyhepmc@master
.
Then pip install DarkNews.
If your installation is successful, you should be able to
- import the module from your python scripts or notebook with
import DarkNews
- run the command line tool
dn_gen
to generate events on the terminal.
To make changes to the package or contribute, you can clone the latest repository from GitHub, and from the main folder of the cloned directory, run:
python3 -m pip install -e .
The package will be installed locally in editable mode. Any changes to your local files in the repo will be automatically propagated to your DarkNews installation (except setup configurations). You may want to use Jupyter notebooks with auto reload
.
If you experience any problems with pip, you can resort to a local manual installation. After downloading the repository, from the main folder, you can run
python3 setup.py develop
to install it in developer mode (similar to editable mode above).
If you would like to output events to .parquet
files, you can install the following pip install DarkNews[parquet]
or pip install "DarkNews[parquet]"
.
A lot of information on the usage of the generator is provided by the example Jupyter notebooks in the repository.
You can download example notebooks from the repository's folder examples/
, or simply run
dn_get_examples
to download the examples folder in the current directory (requires git version >=2.19).
The main usage of DarkNews is covered in depth in the notebook Example_0_start_here.ipynb
.
We encourage you to read through the cells of the notebook.
It is possible to run the generator through the script bin/dn_gen
, passing the parameters as options.
dn_gen --mzprime=1.25 --m4=0.140 --neval=1000 --HNLtype=dirac --loglevel=INFO
Run dn_gen --help
to inspect the meaning of each parameter.
It is possible to run the generator by creating an instance of the DarkNews.GenLauncher.GenLauncher
class and calling its run
method.
from DarkNews import GenLauncher
gen_object = GenLauncher(mzprime=1.25, m4=0.140, neval=1000, HNLtype="dirac")
gen_object.run(loglevel="INFO")
The parameters are passed directly while instantiating the GenLauncher
object.
Some parameters (loglevel
, verbose
, logfile
, path
, overwrite_path
) related to the run itself can be passed also within the call of the run()
method.
If the argument pandas = True
(as it is by default), a dataframe containing the generated events is saved in the directory tree.
The dataframe is multi-indexed, where the second column is empty, then there is only the first index.
All dataframes contain the following process:
P_projectile
) + P_target
) P_decay_N_parent
) + P_recoil
)
which can be followed by a decay into di-leptons if decay_product=e+e-
or decay_product=mu+mu-
:
P_decay_N_parent
) P_decay_N_daughter
) + P_decay_ell_plus
) + P_decay_ell_minus
)
where decay_product='photon'
:
P_decay_N_parent
) P_decay_N_daughter
) + P_decay_photon
). Only one of the above decays is allowed per generation.
Here follows a complete description of the pandas dataframe:
Column | Subcolumn | type | description |
---|---|---|---|
P_projectile | 0, 1, 2, 3 | float |
4-momenta of beam neutrino |
P_decay_N_parent | 0, 1, 2, 3 | float |
4-momenta of HNL_parent |
P_target | 0, 1, 2, 3 | float |
4-momenta of nucleus |
P_recoil | 0, 1, 2, 3 | float |
4-momenta of recoiled nucleus |
P_decay_photon | 0, 1, 2, 3 | float |
4-momenta of photon (if it exists) |
P_decay_ell_minus | 0, 1, 2, 3 | float |
4-momenta of e- (if it exists) |
P_decay_ell_plus | 0, 1, 2, 3 | float |
4-momenta of e+ (if it exists) |
P_decay_N_daughter | 0, 1, 2, 3 | float |
4-momenta of HNL_daughter / nu_daughter |
pos_scatt | 0, 1, 2, 3 | float |
upscattering position |
pos_decay | 0, 1, 2, 3 | float |
decay position of primary particle (N_parent) -- no secondary decay position is saved. |
w_decay_rate_0 | float |
Weight of the decay rate of primary unstable particle: Σi wi = ΓN | |
w_decay_rate_1 | float |
Weight of the decay rate of secondary unstable particle: Σi wi = ΓX | |
w_event_rate | float |
Weight for the event rate: Σi wi = event rate | |
w_flux_avg_xsec | float |
Weight of the flux averaged cross section: Σi wi = int(sigma ⋅ flux) ⋅ exposure | |
target | string | Name of the target object, it will typically be a nucleus | |
target_pdgid | int |
PDG id of the target | |
scattering_regime | string | Regime can be coherent or p-elastic | |
helicity | string | Helicity process: can be flipping or conserving; flipping is suppressed | |
underlying_process | string | String of the underlying process, e.g, "nu(mu) + proton_in_C12 -> N4 + proton_in_C12 -> nu(mu) + e+ + e- + proton_in_C12" |
The pandas DataFrame also contains several useful information in attrs
. They can be accessed via
df.attrs['foo']
and are specific to the generation run. The attributes are detailed below:
Attrs | type | description |
---|---|---|
experiment | DarkNews.detector.Detector() |
The experiment class of DarkNews containing all information on the experimental setup, including neutrino fluxes, targets, exposure, and geometry (if implemented). |
model | DarkNews.model.HNLModel() |
The model class of DarkNews with all the couplings and vertices calculated from the user input. |
data_path | string | The path used to save the generation outputs. |
N{i}_ctau0 | float | The proper decay length of the i-th HNL in centimeters used in the generation of events, with i ={4,5,6}. |
This is an exhaustive list of the parameters that can be passed to the program.
They can be passed in the command line mode by prepending --
to the name.
This summary can also be accessed by running
dn_gen --help
The first column is the name of the parameter, the second is the type or the list of allowed values, the third is a brief explanation and the fourth is the default. Parameters marked as internal can not be specified as they are automatically computed on the basis of other parameters.
mzprime | float |
Mass of Z' | 1.25 |
m4 | float |
Mass of the fourth neutrino | 0.14 |
m5 | float |
Mass of the fifth neutrino | None |
m6 | float |
Mass of the sixth neutrino | None |
HNLtype | ["dirac", "majorana"] |
Dirac or majorana | "majorana" |
ue4 | float |
0.0 | |
ue5 | float |
0.0 | |
ue6 | float |
0.0 | |
umu4 | float |
0.0016202 | |
umu5 | float |
0.0033912 | |
umu6 | float |
0.0 | |
utau4 | float |
0.0 | |
utau5 | float |
0.0 | |
utau6 | float |
0.0 | |
ud4 | float |
1.0 | |
ud5 | float |
1.0 | |
ud6 | float |
1.0 |
gD | float |
U(1)d dark coupling gD | 1.0 |
alphaD | float |
U(1)d αdark = (gD2 / 4 π) | internal |
epsilon | float |
ε | 0.01 |
epsilon2 | float |
ε2 | internal |
alpha_epsilon2 | float |
αQED ⋅ ε2 | internal |
chi | float |
Kinetic mixing parameter | None |
mu_tr_e4 | float |
TMM mu_tr_e4 | 0.0 |
mu_tr_e5 | float |
TMM mu_tr_e5 | 0.0 |
mu_tr_e6 | float |
TMM mu_tr_e6 | 0.0 |
mu_tr_mu4 | float |
TMM mu_tr_mu4 | 0.0 |
mu_tr_mu5 | float |
TMM mu_tr_mu5 | 0.0 |
mu_tr_mu6 | float |
TMM mu_tr_mu6 | 0.0 |
mu_tr_tau4 | float |
TMM mu_tr_tau4 | 0.0 |
mu_tr_tau5 | float |
TMM mu_tr_tau5 | 0.0 |
mu_tr_tau6 | float |
TMM mu_tr_tau6 | 0.0 |
mu_tr_44 | float |
TMM mu_tr_44 | 0.0 |
mu_tr_45 | float |
TMM mu_tr_45 | 0.0 |
mu_tr_46 | float |
TMM mu_tr_46 | 0.0 |
mu_tr_55 | float |
TMM mu_tr_54 | 0.0 |
mu_tr_56 | float |
TMM mu_tr_55 | 0.0 |
mu_tr_56 | float |
TMM mu_tr_66 | 0.0 |
experiment | string |
The experiment to consider: see this section | "miniboone_fhc" |
nopelastic | bool |
Do not generate proton elastic events | False |
nocoh | bool |
Do not generate coherent events | False |
noHC | bool |
Do not include helicity conserving events | False |
noHF | bool |
Do not include helicity flipping events | False |
decay_products | ["e+e-","mu+mu-","photon"] |
Decay process of interest | "e+e-" |
enforce_prompt | bool |
Force particles to decay promptly | False |
nu_flavors | ["nu_e","nu_mu","nu_tau","nu_e_bar","nu_mu_bar","nu_tau_bar"] |
Projectile neutrino | ["nu_mu"] |
sample_geometry | sample_geometry |
Whether to sample the detector geometry using DarkNews. If False or a geometry function is not found, the upscattering position is assumed to be (0,0,0,0). | True |
loglevel | ["INFO", "WARNING", "ERROR", "DEBUG"] |
Logging level | "INFO" |
verbose | bool |
Verbose for logging | False |
logfile | string |
Path to log file; if not set, use std output | None |
neval | int |
Number of evaluations of integrand | 10000 |
nint | int |
Number of adaptive iterations | 20 |
neval_warmup | int |
Number of evaluations of integrand in warmup | 1000 |
nint_warmup | int |
Number of adaptive iterations in warmup | 10 |
seed | int |
numpy random number generator seed used in vegas | None |
pandas | bool |
Save pandas.DataFrame as .pckl file |
True |
parquet | bool |
Save pandas.DataFrame as .parquet file (engine=pyarrow) |
False |
numpy | bool |
Save events in .npy files |
False |
hepevt | bool |
If true, print events to HEPEVT-formatted text files (does not save event weights) | False |
hepevt_legacy | bool |
If true, print events to a legacy HEPEVT format (saving weights next to the number of particle in the event and without linebreaks in particle entries) | False |
hepmc2 | bool |
If true, prints events to HepMC2 format. | False |
hepmc3 | bool |
If true, prints events to HepMC3 format. | False |
hep_unweight | bool |
Unweigh events when printing in HEPEVT format (needs large statistics) | False |
unweighted_hep_events | int |
number of unweighted events to accept in any of the standard HEP formats. Has to be much smaller than neval for unweight procedure to work. | 100 |
sparse | int |
Specify the level of sparseness of the internal dataframe and output. Not supported for HEPevt. Allowed values are 0--3, where: 0 : keep all information, including event-by-event descriptions; 1 : keep all particle momenta, scattering and decay positions, and all weights; 2 : keep only neutrino energy (or 4-momentum in HEPMC/EVT), visible decay products and unstable particle momenta, scattering and decay positions, and all weights; 3 : keep only neutrino energy (or 4-momentum in HEPMC/EVT), visible decay products and unstable particle momenta, and all weights; 4 : keep only neutrino energy (or 4-momentum in HEPMC/EVT), unstable particle momenta, visible decay products momenta, and w_event_rate. Metadata is always kept if output is pickled. |
0 |
path | string |
Path where to save run's outputs | "./" |
make_summary_plots | bool |
if True, generates summary plots of kinematics in the path |
False |
It is possible to specify the parameters through a file, instead of writing each one as an option for the command line functionality or as an argument of GenLauncher
instance.
The parameter file should be provided as the argument param_file
of GenLauncher
or via the option --param-file
of the command line interface.
A template file template_parameters_file.txt
can be found in in the examples
directory.
In the following there are the main rules to specify the parameters:
- every line should be in the form of an assignment statement
<variable_name> = <variable_value>
-
comments can be specified with with
"#"
-
the different parameters can be specified with:
- string: by encapsulating each string with double or single quotes
"<string>"
or'<string>'
are equivalent, the escape character is\
(backslash). - booleans: by writing
True
orFalse
(it is case insensitive) - mathematical expression (which will results in
float
orint
numbers): see section below - lists: by encapsulating them into square brackets, separating each element with a comma; elements can be string, numbers, mathematical expressions or all of them together.
- string: by encapsulating each string with double or single quotes
-
When using mathematical expression, the following rules should be applied:
- numbers can be specified as usual:
5
is integer,5.0
is float (but5.
will result in an error),5e10
is the float number 5*10^10 +
,-
,*
,/
are the usual mathematical operators;^
is used to make powers (do not use**
);- it is possible to use round brackets
(
and)
e
(case-insensitive, isolated: not inside float numbers) is understood as pythonmath.e = 2.718281828
pi
(case-insensitive) is understood asmath.pi = 3.1415926535
sin(<expr>)
,cos(<expr>)
,tan(<expr>)
are the usual trigonometric functionsexp(<expr>)
is the usual exponentiationabs(<expr>)
is the absolute valuesgn(<expr>) = -1
if<expr> < -1e-100
,+1
if<expr> > 1e-100
,0
otherwisetrunc(<expr>)
returns the truncated float<expr>
to integerround(<expr>)
is the integer part of the float number<expr>
sum(<list>)
will sum each element of the list, returning a float number- any other string is intended as a variable which must have been previously defined (the file is scanned from top to bottom)
- in general it is possible to define any kind of variable, independently on those that will be actually used by the program, following the usual conventions for the variable name (use only letters, digits and underscores). Moreover, it's not possible to name variables after program-defined names, as for example the name of the functions.
- numbers can be specified as usual:
The following lines
hbar = 6.582119569e-25 # GeV s
c = 299792458.0 # m s^-1
will define two variables, named hbar
and c
with their values.
It is possible to write
a_certain_constant = hbar * c
to define a variable named a_certain_constant
with the value of the product between the pre-defined hbar
and c
variables from the example above.
It is possible to write any kind of possible expression, for example
a_variable = c^2 * 3.2e-4 / sin(PI/7) + 12 * exp( -2 * abs(hbar) )
obtaining a new variable a_variable
with the value of 66285419633555.3
The line
path = "my_directory/projects/this_project"
defines the path
variable, stored as the string "my_directory/projects/this_project"
.
The following lines are defining booleans (they are case insensitive), used to set the various switches:
pandas = True
numpy = false
The experiment to use can be specified in two ways through the experiment
argument (or --experiment
option accordingly):
- specifying a keyword for a pre-defined experiment among:
- DUNE FHC ND (
"dune_nd_fhc"
) - DUNE RHC ND (
"dune_nd_rhc"
) - SBND (
"sbnd"
) - MicroBooNE (
"microboone"
) - MINERVA FHC LE (
"minerva_le_fhc"
) - MINERVA FHC ME (
"minerva_me_fhc"
) - MiniBooNE FHC (
"miniboone_fhc"
) - NUMI FHC ME (
"minos_le_fhc"
) - NUMI FHC LE (
"minos_me_fhc"
) - ND280 FHC (
"nd280_fhc"
) - NOva FHC LE (
"nova_le_fhc"
) - FASERnu (
"fasernu"
)
- DUNE FHC ND (
- specifying the file path of an experiment file: every file should be specified using the same rules as for the parameters file, listed in the previous section.
A template file
template_custom_experiment.txt
can be found in in theexamples
directory. The following parameters must be present (in general it is possible to specify any number of parameters, but only the ones below would be relevant).
name | string |
Name of the experiment (your are free to use capital letters, when needed) |
fluxfile | string |
Path of the fluxes file with respect to the experiment file directory |
flux_norm | float |
Flux normalization factor: all fluxes should be normalized so that the units are nus/cm²/GeV/POT |
erange | list of float |
Neutrino energy range [<min>, <max>] in GeV |
nuclear_targets | list of string |
Detector materials in the form of "<element_name><mass_number>" (e.g. "Ar40" ) |
fiducial_mass | float |
Fiducial mass in tons |
fiducial_mass_per_target | list of float |
Fiducial mass for each target in the same order as the nuclear_targets parameter |
POTs | float |
Protons on target |
The only geometries currently implemented in DarkNews are:
- MiniBooNE -- a sphere with origin (0,0,0,0) and radius 574.6 cm.
- SBND -- a box with sides 4m x 4m x 5m.
- MicroBooNE -- the geometry of the cryostat. This is just a junction of a tube with two spherical caps.
Times of interactions are always set to 0, and any additional delay due to the N lifetime is taken in to account. All other experiments spit out events with spatial coordinates (0,0,0).
DarkNews relies on vegas to integrate and sample differential cross sections and decay rates.
The main point of contact with vegas is through the vegas.Integrator
class, which will receive the DarkNews integrands (e.g. DarkNews.integrands.UpscatteringHNLDecay()
), whose __call__()
method will compute the differential rates.
It is possible to set a seed for numpy's random number with the --seed
argument, which accepts integer values. This integer seed is then passed to numpy.seed().
By default, vegas uses numpy's random number generator, which is based on the Mersenne Twister pseudo-random number generator method.
The reproducibility of the generator output (i.e., same vegas samples) is only guaranteed for the same parameters and number of integrand evaluations, which effectively means that the user has to specify the same scope, model parameters, as well as the same number of neval, nint, neval_warmup and nint_warmup.