diff --git a/.gitignore b/.gitignore index 80b79d73..4ce8dcca 100644 --- a/.gitignore +++ b/.gitignore @@ -101,3 +101,5 @@ ENV/ .mypy_cache/ *.root *.png +*.pdf +*.DS_Store diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000..5394f134 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,31 @@ +# Read the Docs configuration file for Sphinx projects +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the OS, Python version and other tools you might need +build: + os: ubuntu-22.04 + tools: + python: "3.11" + +# Build documentation in the "docs/" directory with Sphinx +sphinx: + configuration: documentation/conf.py + # You can configure Sphinx to use a different builder, for instance use the dirhtml builder for simpler URLs + # builder: "dirhtml" + # Fail on all warnings to avoid broken references + # fail_on_warning: true + +# Optionally build your docs in additional formats such as PDF and ePub +# formats: +# - pdf +# - epub + +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +python: + install: + - requirements: documentation/requirements.txt \ No newline at end of file diff --git a/Dockerfile b/Dockerfile index 2c7fb38e..91af0e82 100644 --- a/Dockerfile +++ b/Dockerfile @@ -24,6 +24,7 @@ RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\ RUN source $COMMON_BUILD_PREFIX/setup.sh &&\ pip install iminuit &&\ + pip install numericalunits &&\ /bin/true ######################## diff --git a/README.md b/README.md index cac850c4..d4af8611 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ Mermithid is an extension of [morpho](https://github.com/morphoorg/morpho) that ## Requirements -You will need to install via a package manager (such as apt-get): +If you are not using a container with pre-installed dependencies, you will need to install via a package manager (such as apt-get): - python (3.x; 2.7.x support not guaranteed) - python-pip @@ -19,6 +19,9 @@ These are two possible ways of installing and working with mermithid. ### Virtual environment installation +Before installing, clone subdirecories recursively: ``git submodule update --init --recursive``. +
Then install mermithid in your environment: + 1. Cicada and Phylloxera need to be installed in a sub directory: ```bash @@ -45,16 +48,24 @@ These are two possible ways of installing and working with mermithid. ### Docker installation -Docker provides a uniform test bed for development and bug testing. Please use this environment to testing/resolving bugs. +Create a docker image and start a mermithid container 1. Install Docker (Desktop version): -2. Clone and pull the latest master version of mermithid -3. Inside the mermithid folder, execute `docker-compose run mermithid`. The container prompter should appear at the end of the installation. A directory (`mermithid_share`) should be created in your home and mounted under the `/host` folder: you can modify this by editing the docker-compose file. -4. When reinstalling, you can remove the image using `docker rmi mermithid_mermithid` +3. Clone and pull the latest main version of mermithid or the feature branch you want to work with +4. Go to the cloned directory: ``cd mermithid`` +5. Pull the submodules: ``git submodule update --init --recursive`` +6. Build docker image: ``docker build --no-cache -t mermithid: .`` +7. To start the container and mount a directory for data sharing with your host into the container do: +
```docker run --rm -it -v ~/mermithid_share:/host mermithid: /bin/bash``` + +An alternative to steps 6 and 7 is to use docker-compose by executing: ``docker-compose run mermithid``. + +
In both cases files saved in ```/host``` are shared with the host in ```~/mermithid_share```. + ### Running mermithid -In both cases, you need to set the paths right for using these software. For example in the docker container: +For running mermithid you need to set the paths right for using these software. For example in the docker container: ```bash source $MERMITHID_BUILD_PREFIX/setup.sh @@ -64,4 +75,11 @@ source $MERMITHID_BUILD_PREFIX/bin/this_phylloxera.sh ## Quick start and examples -Mermithid works a-la morpho, where the operations on data are defined using processors. Each processor should be defined with a name, then should have its attributes configured using a dictionary before being run. Examples of how to use mermithid can be found in the "tests" folder. +Mermithid works a-la morpho, where the operations on data are defined using processors. Each processor should be defined with a name, then should have its attributes configured using a dictionary before being run. Examples of how to use mermithid can be found in the "tests" and the "test_analysis" folders. + + +## Easy development + +To develop mermithid without having to rebuild the container, share the repository on the host with the container by starting it with: ```docker run --rm -it -v ~/mermithid_share:/host -v ~/repos/mermithid:/mermithid mermithid: /bin/bash```. This assumes that mermithid was cloned to ``~/repos``. + +After sourcing the setup scripts, modify the PYTHONPATH: ```export PYTHONPATH=/mermithid:$PYTHONPATH```. Now changes made on the host will directly be used by the container. diff --git a/docker-compose.yaml b/docker-compose.yaml index 3b701007..cb13bac2 100644 --- a/docker-compose.yaml +++ b/docker-compose.yaml @@ -1,9 +1,20 @@ -version: "3" - +#version: "3.9" +#Check running containers docker ps (detached) -d +#Attach to container docker attach [name] services: mermithid: - build: . - command: "/bin/bash" + image: mermithid:sensitivity_branch + container_name: mermithid_sensitivity + stdin_open: true # Interactive search + tty: true # Allocate TTY volumes: - # share a subdirectory from the host to /host in the docker (can be edited) - - ~/mermithid_share:/host + - ~/repos/mermithid:/mermithid + - ~/repos/mermithid_share:/host_data + - ~/repos/termite:/termite + command: > + bash -c " + source /usr/local/p8/mermithid/v1.2.3/setup.sh && + source /usr/local/p8/mermithid/v1.2.3/bin/this_phylloxera.sh && + source /usr/local/p8/mermithid/v1.2.3/bin/this_cicada.sh && + exec bash + " diff --git a/documentation/Makefile b/documentation/Makefile index 8f0bb615..3fd5dcea 100644 --- a/documentation/Makefile +++ b/documentation/Makefile @@ -8,8 +8,6 @@ SPHINXPROJ = mermithid SOURCEDIR = . BUILDDIR = build -PY_CMD = import better_apidoc; better_apidoc.main(['better_apidoc','-t', '_templates','--force','--separate','-o','better_apidoc_out','../$(SPHINXPROJ)']) # overwrite existing files -APIDOC_CMD = python -c "$(PY_CMD)" # Put it first so that "make" without argument is like "make help". help: @@ -22,5 +20,3 @@ help: %: Makefile @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) -betterapi: - @$(APIDOC_CMD) \ No newline at end of file diff --git a/documentation/conf.py b/documentation/conf.py index 1209ea2b..f801a83b 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -1,311 +1,31 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- +# Configuration file for the Sphinx documentation builder. # -# morpho documentation build configuration file, created by -# sphinx-quickstart on Mon Aug 25 15:26:40 2014. -# -# This file is execfile()d with the current directory set to its -# containing dir. -# -# Note that not all possible configuration values are present in this -# autogenerated file. -# -# All configuration values have a default; values that are commented out -# serve to show the default. - -import sys -import os - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# sys.path.insert(0, os.path.abspath('.')) - - -# -- General configuration ------------------------------------------------ - -# If your documentation needs a minimal Sphinx version, state it here. -# needs_sphinx = '1.0' - -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom -# ones. -extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.todo', - 'sphinx.ext.viewcode', - # 'sphinxcontrib.programoutput', - 'sphinx.ext.napoleon', - # 'sphinxarg.ext', -] - -# better_apidoc.main([ -# 'better_apidoc', -# '-t', '_templates', # path to jinja templates for generated package and module rst files -# '--force', # overwrite existing files -# '--separate', # split the modules into their own files -# # output location, be sure to update index.rst if you change this -# '-o', 'better_apidoc_out', -# '../morpho', # path to the package containing modules to document ##TODO update this with your path -# ]) - - -def run_apidoc(_): - """Generage API documentation""" - import better_apidoc - better_apidoc.main( - ['better-apidoc', '-t', './_templates', '--force', '--no-toc', - '--separate', '-o', 'better_apidoc_out', '../mermithid']) - - -def setup(app): - app.connect('builder-inited', run_apidoc) - -todo_include_todos = True - -autoclass_content = "both" - -# Add any paths that contain templates here, relative to this directory. -# templates_path = ['_templates'] - -# The suffix of source filenames. -source_suffix = '.rst' +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html -# The encoding of source files. -# source_encoding = 'utf-8-sig' +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information -# The master toctree document. -master_doc = 'index' - -# General information about the project. -# TODO change this line to match the package name project = 'mermithid' -copyright = '2018, The Project 8 Collaboration' - -# The version info for the project you're documenting, acts as replacement for -# |version| and |release|, also used in various other places throughout the -# built documents. -# -# The short X.Y version. -# version = #'0.0.1' -# import pkg_resources -# release: The full version, including alpha/beta/rc tags. -try: - import subprocess - release = subprocess.check_output( - ['git', 'describe', '--long']).decode('utf-8').strip() - version = subprocess.check_output( - ['git', 'describe', '--abbrev=0']).decode('utf-8').strip() -except Exception as e: - print("error message is:\n{}".format(e.message)) - version = "unknown" - release = "unknown" -print('version/release are: {}/{}'.format(version, release)) - -# The language for content autogenerated by Sphinx. Refer to documentation -# for a list of supported languages. -# language = None - -# There are two options for replacing |today|: either, you set today to some -# non-false value, then it is used: -# today = '' -# Else, today_fmt is used as the format for a strftime call. -# today_fmt = '%B %d, %Y' - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -exclude_patterns = ['_build'] - -# The reST default role (used for this markup: `text`) to use for all -# documents. -# default_role = None - -# If true, '()' will be appended to :func: etc. cross-reference text. -# add_function_parentheses = True - -# If true, the current module name will be prepended to all description -# unit titles (such as .. function::). -# add_module_names = True - -# If true, sectionauthor and moduleauthor directives will be shown in the -# output. They are ignored by default. -# show_authors = False - -# The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' - -# A list of ignored prefixes for module index sorting. -# modindex_common_prefix = [] - -# If true, keep warnings as "system message" paragraphs in the built documents. -# keep_warnings = False - - -# -- Options for HTML output ---------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# try: -import sphinx_rtd_theme -html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] -html_theme = 'sphinx_rtd_theme' -# except ImportError: -# html_theme = 'haiku' - -# Theme options are theme-specific and customize the look and feel of a theme -# further. For a list of options available for each theme, see the -# documentation. -# html_theme_options = {} - -# Add any paths that contain custom themes here, relative to this directory. -# html_theme_path = [] - -# The name for this set of Sphinx documents. If None, it defaults to -# " v documentation". -# html_title = None - -# A shorter title for the navigation bar. Default is the same as html_title. -# html_short_title = None - -# The name of an image file (relative to this directory) to place at the top -# of the sidebar. -# html_logo = None - -# The name of an image file (within the static path) to use as favicon of the -# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 -# pixels large. - -# TODO put your favicon here -# html_favicon = 'morpho_favicon.ico' - -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -# html_static_path = ['_static'] - -# Add any extra paths that contain custom files (such as robots.txt or -# .htaccess) here, relative to this directory. These files are copied -# directly to the root of the documentation. -# html_extra_path = [] - -# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, -# using the given strftime format. -# html_last_updated_fmt = '%b %d, %Y' - -# If true, SmartyPants will be used to convert quotes and dashes to -# typographically correct entities. -# html_use_smartypants = True - -# Custom sidebar templates, maps document names to template names. -# html_sidebars = {} - -# Additional templates that should be rendered to pages, maps page names to -# template names. -# html_additional_pages = {} - -# If false, no module index is generated. -# html_domain_indices = True - -# If false, no index is generated. -# html_use_index = True - -# If true, the index is split into individual pages for each letter. -# html_split_index = False - -# If true, links to the reST sources are added to the pages. -# html_show_sourcelink = True - -# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. -# html_show_sphinx = True - -# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. -# html_show_copyright = True - -# If true, an OpenSearch description file will be output, and all pages will -# contain a tag referring to it. The value of this option must be the -# base URL from which the finished HTML is served. -# html_use_opensearch = '' - -# This is the file name suffix for HTML files (e.g. ".xhtml"). -# html_file_suffix = None - -# Output file base name for HTML help builder. -htmlhelp_basename = 'mermithiddoc' - - -# -- Options for LaTeX output --------------------------------------------- - -latex_elements = { - # The paper size ('letterpaper' or 'a4paper'). - # 'papersize': 'letterpaper', - - # The font size ('10pt', '11pt' or '12pt'). - # 'pointsize': '10pt', - - # Additional stuff for the LaTeX preamble. - # 'preamble': '', -} - -# Grouping the document tree into LaTeX files. List of tuples -# (source start file, target name, title, # author, documentclass [howto, manual, or own class]). -latex_documents = [ - ('index', 'mermithid.tex', 'mermithid Documentation', - 'The Project 8 Collaboration', 'manual'), -] - -# The name of an image file (relative to this directory) to place at the top of -# the title page. -# latex_logo = None - -# For "manual" documents, if this is true, then toplevel headings are parts, -# not chapters. -# latex_use_parts = False - -# If true, show page references after internal links. -# latex_show_pagerefs = False - -# If true, show URL addresses after external links. -# latex_show_urls = False - -# Documents to append as an appendix to all manuals. -# latex_appendices = [] - -# If false, no module index is generated. -# latex_domain_indices = True - - -# -- Options for manual page output --------------------------------------- +copyright = '2024, The Project 8 Collaboration' +author = 'The Project 8 Collaboration' +release = '2018' -# One entry per manual page. List of tuples -# (source start file, name, description, authors, manual section). -man_pages = [ - ('index', 'mermithid', 'mermithid Documentation', - ['The Project 8 Collaboration'], 1) -] +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration -# If true, show URL addresses after external links. -# man_show_urls = False +extensions = [] +templates_path = ['_templates'] +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store', '.md'] -# -- Options for Texinfo output ------------------------------------------- -# Grouping the document tree into Texinfo files. List of tuples -# (source start file, target name, title, author, -# dir menu entry, description, category) -texinfo_documents = [ - ('index', 'mermithid', 'mermithid Documentation', - 'The Project 8 Collaboration', 'mermithid', 'One line description of project.', - 'Miscellaneous'), -] -# Documents to append as an appendix to all manuals. -# texinfo_appendices = [] +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output -# If false, no module index is generated. -# texinfo_domain_indices = True +html_theme = 'sphinx_rtd_theme' #'alabaster' +html_static_path = ['_static'] -# How to display URL addresses: 'footnote', 'no', or 'inline'. -# texinfo_show_urls = 'footnote' -# If true, do not generate a @detailmenu in the "Top" node's menu. -# texinfo_no_detailmenu = False +extensions = ['sphinxcontrib.contentui', 'sphinx_rtd_theme'] \ No newline at end of file diff --git a/documentation/contribute.rst b/documentation/contribute.rst index 2085fd76..dd6c9c53 100644 --- a/documentation/contribute.rst +++ b/documentation/contribute.rst @@ -34,7 +34,8 @@ Other Conventions In mermithid 1, __init__.py files are set up such that :: - from package import * + + from package import * will import all functions from all subpackages and modules into the namespace. If a package contains the subpackages "subpackage1" and "subpackage2", and the modules "module1" and "module2", then the __init__.py file should include imports of the form: :: @@ -54,5 +55,5 @@ will import all modules into the namespace, but it will not directly import the __all__ = ["module1", "module2"] -In this case, functions would be called via module1.function_name(). If one wants all of the functions from module1 in the namespace, then they can include "from package.module1 import *" at the top of their code. This change to more explicit imports should prevent any issues with function names clashing as mermithid grows. +In this case, functions would be called via module1.function_name(). If one wants all of the functions from module1 in the namespace, then they can include "from package.module1 import \*" at the top of their code. This change to more explicit imports should prevent any issues with function names clashing as mermithid grows. diff --git a/documentation/index.rst b/documentation/index.rst index c9ae02e5..d4667d0e 100644 --- a/documentation/index.rst +++ b/documentation/index.rst @@ -1,13 +1,28 @@ -Welcome to mermithid's documentation! -==================================== +.. mermithid documentation master file, created by + sphinx-quickstart on Mon May 20 13:35:50 2024. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. -Contents: +Welcome to mermithid's documentation! +===================================== .. toctree:: - :maxdepth: 2 + :maxdepth: 2 + :caption: Contents: + + intro + howmermithidworks + install + contribute + sensitivity + sensitivity_configurable_parameters + validation_log + + + +Indices and tables +================== - intro - install - contribute - validation_log - better_apidoc_out/modules +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/documentation/requirements.txt b/documentation/requirements.txt new file mode 100644 index 00000000..0f433de6 --- /dev/null +++ b/documentation/requirements.txt @@ -0,0 +1,2 @@ +sphinxcontrib-contentui +sphinx_rtd_theme \ No newline at end of file diff --git a/documentation/sensitivity.rst b/documentation/sensitivity.rst new file mode 100644 index 00000000..5d54b159 --- /dev/null +++ b/documentation/sensitivity.rst @@ -0,0 +1,84 @@ +-------------------------------------- +Neutrino Mass Sensitivity Calculation +-------------------------------------- + + +Scripts used in sensitivity calculations +------------------------------------------- + +The mermithid sensitivity calculations are written in python and a python script can be used to configure the calculation and make sensitivity plots. +Mermithid has several processors for this purpose, in `mermithid/mermithid/processors/Sensitivity/`_. For cavity sensitivity calculations, use the `CavitySensitivityCurveProcessor`_. + +.. _mermithid/mermithid/processors/Sensitivity/: https://github.com/project8/mermithid/tree/feature/sensitivity_curve/mermithid/processors/Sensitivity +.. _CavitySensitivityCurveProcessor: https://github.com/project8/mermithid/blob/feature/sensitivity_curve/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py + +Mermithid processors are all designed to be used in the same way: + +1. Define a dictionary with the processor's configurable parameters; +2. Instantiate a processor and pass the configuration dictionary to its ``Configure()`` method; +3. Call the processor's ``Run()`` method to make it perform its task. + +In `mermithid/tests`_, the ``Sensitivity_test.py`` script contains several examples of how to perform sensitivity calculations using mermithid processors. In ``test_SensitivityCurveProcessor``, the ``CavitySensitivityCurveProcessor`` is used as described above to calculate and plot sensitivity to the neutrino mass as a function of gas density. +Other working examples for creating sensivitiy plots vs. frequency or exposure can be found in `mermithid/test_analyses/Cavity_Sensitivity_analysis.py`_. + +.. _mermithid/tests: https://github.com/project8/mermithid/blob/feature/sensitivity_curve/tests +.. _mermithid/test_analyses/Cavity_Sensitivity_analysis.py: https://github.com/project8/mermithid/blob/feature/sensitivity_curve/test_analysis/Cavity_Sensitivity_analysis.py + +The dictionary to configure the processor (see #1, above) can read in a separate configuration file (``.cfg``) with sensitivity-specific parameters. The documentation page `sensitivity_configurable_parameters.rst`_ describes all parameters in such a ``.cfg``. + +.. _sensitivity_configurable_parameters.rst: https://github.com/project8/mermithid/blob/feature/sensitivity_curve/documentation/sensitivity_configurable_parameters.rst + + +Analytic sensitivity formula +----------------------------------- +Test + + +Systematic uncertainty contributions +------------------------------------- + +The following contributions to energy broadening of the beta spectrum are included: + +1. ``sigma_trans``: Translational Doppler broadening due to thermal motion of tritium atoms or molecules. +2. ``sigma_f``: Energy broadening due to mis-reconstruction of the start frequencies of events from the range of electron pitch angles (after axial frequency corrections), and from the process of fitting chirp tracks that have finite length and SNR to find where the tracks start. +3. ``sigma_B``: Energy broadening due to radial, azimuthal, and temporal varation of the magnetic field. This is the remaining broadening after any event-by-event corrections for the field have been performed. +4. ``sigma_Miss``: Energy broadening due to the detection of events after one or more tracks have been missed. +5. ``sigma_Plasma``: Energy broadening due to plasma effects of the charged particles in the source volume of an atomic tritium experiment. + +In the list above, each variable that starts with ``sigma`` is a standard deviation of a distribution of energies. The measured spectrum is the convolution of the underlying beta spectrum with such distributions of energies. Asymmetries in these distributions are not accounted for in this approximate model. + +Each ``sigma`` has an associated ``delta``, which is the uncertainty on ``sigma`` from calibration, theory, and/or simulation (see the previous section). For example, ``sigma_trans`` has an associated variable ``delta_sigma_trans``. + +In the script ``mermithid/mermithid/misc/SensitivityCavityFormulas.py``, lists of these energy broadening standard deviations (sigmas) and uncertainties on them (deltas) are returned by the ``get_systematics`` method of the ``CavitySensitivity`` class. + +Contributions 4 and 5 are simply inputted in the sensitivity configuration file; we do not yet have a way to calculate these in mermithid. Contributions 1, 2, and 3 are calculated in mermithid, as described below. + + +Translational Doppler broadening (``sigma_trans``) +======================================================== +The thermal translational motion of tritium atoms causes a Doppler broadening of the :math:`{\beta}` energy spectrum. ``sigma_trans`` is the standard deviation of this broadening distribution. There are two options for how to include translational Doppler broadening in your sensitivity calculations, in mermithid: + +1. Manually input values for ``sigma_trans`` and its uncertainty ``delta_trans``, calculated outside of mermithid. +This is done in the ``DopplerBroadening`` section of the configuration file, by setting ``UseFixedValue`` to ``True`` and providing values for ``Default_Systematic_Smearing`` and ``Default_Systematic_Uncertainty``. + +2. Have mermithid calculate ``sigma_trans``, for you (default option). + - For an atomic tritium experiment, mermithid will also calculate ``delta_trans``, conservatively assuming that the gas temperature uncertainty will be determined simply by our knowledge that the gas is not so hot that it escapes the atom trap. These calculations for an atomic experiment are described further, below. + - For a molecular tritium experiment, you need to input a number ``fraction_uncertainty_on_doppler_broadening`` (which equals ``delta_trans``/``sigma_trans``) in the ``DopplerBroadening`` section of the configuration file. + +Calculation of ``sigma_trans`` for option 2: +For a thermalized source gas, the translational Doppler broadening is described by a Gaussian with standard deviation +:math:`{\sigma_{\text{trans}} = \sqrt{\frac{p_{\text{rec}}^2}{2m_T}2 k_B T}}`, +where :math:`{m_T}` is the mass of tritium and :math:`{T}` is the gas temperature. + + +Calculation of ``delta_trans`` for option 2, with atomic T: + + + +Track start frequency determination and pitch angle correction (``sigma_f``) +==================================================================================== + + +Radial, azimuthal, and temporal field broadening (``sigma_B``) +==================================================================================== + diff --git a/documentation/sensitivity_configurable_parameters.rst b/documentation/sensitivity_configurable_parameters.rst new file mode 100644 index 00000000..ab42affd --- /dev/null +++ b/documentation/sensitivity_configurable_parameters.rst @@ -0,0 +1,115 @@ +---------------------------------------------------- +Configuring mermithid sensitivity calculation +---------------------------------------------------- + +The sensitivity calculation is configured using a configuration file. The configuration files for Project 8 live in the termite repository: https://github.com/project8/termite/tree/feature/sensitivity_config_files/sensitivity_config_files + +Our main goal is the calcualtion of sensitivity in a cavity experiment. The configurations below are to be used for the CavitySensitivity class in https://github.com/project8/mermithid/blob/feature/sensitivity_curve/mermithid/misc/SensitivityCavityFormulas.py +This class is used by the CavitySensitivityCurveProcessor and the SensitivityParameterScanProcessor. + +Structure of a config file +-------------------------- + +Configuration files have several sections: + + +* Experiment +* Efficiency +* FrequencyExtraction +* DopplerBroadening +* MagneticField +* FinalStates +* MissingTracks +* PlasmaEffects + +Each section has a number of parameters that are used for the calcualtion of the sensitivty contribution of the respective section. + + +Parameters +---------- + +Below is a list of the parameters with a short description of what role they play in the calculation. + +**Experiment** + +* ``L_over_D``: The ratio of the length of the cavity to the diameter of the cavity. Together with the frequency of the cavity's TE011 mode, this parameter determines the length and volume of a single cavity. It also impacts the cavity Q-factor which is determined by the bandwidth needed to observe the axial frequency of the trapped electrons. +* ``livetime``: Run duration of the experiment. Together with the total volume, efficiency, and the gas density, this determines the statistical power of the experiment. +* ``n_cavities``: Number of identical cavities in the experimennt. This parameter multiplies the single cavity volume to give the total experimental volume. +* ``background_rate_per_ev``: Background rate per electronvolt for the entire experiment. It is not automultiplied by the number of channels (cavities) in the experiment. +* ``number_density``: This is the mean tritium gas density between the electron trapping coils. The density make be different at axial positions past the trap coils. The number_density together with the total volume, livetime, and the efficiency determines the statistical power of the experiment. Gas density also determines the track length and therefore the frequency resolution. The sensitivity curve processor can optimize this parameter to maximize the sensitivity. In that case this number is overwritten in the calculation. +* ``sri_factor``: The statistical rate increase factor articifially increases the number of observed events (it multiplies the total efficiency). It is highly recommended to set it to 1. +* ``atomic``: If true, the calculation is done for atomic tritium. If false, moecular tritium is assumed. This affects the number of decays per gas molecule/atom (2 for molecular 1 for atomic), the track length in a given gas density (via electron scattering cross section), and the width of the final ground state. + + +**Efficiency** + +* ``usefixedvalue``: If true, fixed efficiency is used. If false, the efficiency is the product of radial, detection, and trapping efficienc. The trapping efficiency is calculated from the minimum pitch angle. +* ``fixed_efficiency``: For example, set to roughly 2% for a 88deg minimum trapped pitch angle, assuming 100% detection efficiency of the trapped angles. +* ``radial_efficiency``: Typically set to 0.67 from a calcualtion done for a 325MHz cavity with Halbach bite and radial cut on power of > 0.5 * maximum power. +* ``detection_efficiency``: Fraction of events that is not detected. + +**FrequencyExtraction** + +We use the CRLB for calculating the frequency resolution. The CRLB is calculated from the signal to noise ratio (SNR) and track length. The sensitivity calculations therefore include a calculation of the noise and signal power. + +* ``usefixedvalue``: If true all parameters below are ignored but ``default_systematic_smearing`` and ``default_systematic_uncertainty``. +* ``default_systematic_smearing``: Used if ``usefixedvalue`` is true. Units must be eV. +* ``default_systematic_uncertainty``: Used if ``usefixedvalue`` is true. Units must be eV +* ``usefixeduncertainty``: If true, the uncertainty on the frequency extraction is fixed to the value of ``fixed_relativ_uncertainty``. False will result in an error because no calculation is currently implemented. +* ``fixed_relativ_uncertainty``: If ``usefixeduncertainty`` is true, this relative value is used as the uncertainty on the frequency extraction. +* ``crlb_on_sidebands``: If true, the total resolution from frequency extraction includes the correction from axial field variation based on the CRLB resolution of sidebands. +* ``sideband_power_fraction``: Fraction of power in the sidebands. This is used to calculate the CRLB resolution of the sidebands. +* ``sideband_order``: The order of the sidebands used to calculate the axial field correction. 2 gives better precision but would require lower Q which in turn affects SNR. However, the sideband order is currently not included in the Q calculation. +* ``amplifier_temperature``: Used to calculate noise temperature. +* ``quantum_amp_efficiency``: Used to calculate noise temperature. +* ``att_cir_db``: Used to calculate noise temperature. +* ``att_line_db``: Used to calculate noise temperature. +* ``cavity_temperature``: Used to calculate noise temperature. Has to be compatible with the gas species (molecular tritium freezes below 30K). +* ``unloaded_q``: The unloaded Q of hte cavity. From the axial frequency of the minimum trapped angle (depending on ``L_over_D`` and the cavity frequency) and the bandwidth needed to observe it, the loaded Q-factor is calculated. The unloaded Q-factor is used to calculate the coupling and therefore impacts the SNR. +* ``minimum_angle_in_bandwidth``: Minimum pitch angle of which the axial frequency is contained in the bandwidth. Impacts the loaded Q and therefore the SNR. In the future it will be linked to the efficiency above. +* ``crlb_scaling_factor``: Arbitrary factor to scale the resolution. +* ``magnetic_field_smearing``: A magnetic field smearing in eV can be added here. + +**DopplerBroadening** + +* ``usefixedvalue``: If True ``default_systematic_smearing`` and ``default_systematic_uncertainty`` are used. +* ``default_systematic_smearing``: Default systematic broadening for this category. Units must be eV. +* ``default_systematic_uncertainty``: Default systematic uncertainty for this category. Units must be eV. +* ``gas_temperature``: Temperature of the source gas. This should only be different from the cavity temperature if the gas is not in thermal equilibrium with the cavity. The gas temperature is used to calculate the Doppler broadening. +* ``gas_temperature_uncertainty``: Absolute uncertainty of the gas temperature. +* ``fraction_uncertainty_on_doppler_broadening``: Fractional uncertainty on the Doppler broadening. + + +**MagneticField** + +* ``usefixedvalue``: If True ``default_systematic_smearing`` and ``default_systematic_uncertainty`` are used. +* ``default_systematic_smearing``: Default systematic broadening for this category. Units must be eV. +* ``default_systematic_uncertainty``: Default systematic uncertainty for this category. Units must be eV. +* ``nominal_field``: Determines the CRES and cavity TE011 mode frequency. The cavity dimensions are derived from this and ``L_over_D`` +* ``useinhomogeneity``: True +* ``fraction_uncertainty_on_field_broadening``: Fractional uncertainty on field inhomogeneity. Applies to all parameters below +* ``sigma_meanb``: Fixed input in eV. Magnetic field instability (which is not fully corrected using live calibration) and unknown wiggles in the z-field profile, relative to a smooth trap shape. +* ``sigmae_r``: Fixed input in eV. Energy broadening from radial field inhomogeneity that remains after radial reconstruction. Accounts for both the uncertainty on each electron's radius and the uncertainty on the radial field profile. +* ``sigmae_theta``: Fixed input in eV. Energy broadening remaining after theta reconstruction, from electrons with lower pitch angles exploring high fields. Accounts for both the uncertainty on theta and uncertainties on the trap depth/boxiness. +* ``sigmae_phi``: Fixed input in eV. Energy broadening from phi field inhomogeneity that remains after phi reconstruction. + +**FinalStates** + +* ``ground_state_width_uncertainty_fraction``: Uncertainty on the ground state width. Recommended to use 0.001. + + +The sections below have so far not been used and are assumed to be negligible. + +**MissingTracks** + +* ``usefixedvalue``: If True ``default_systematic_smearing`` and ``default_systematic_uncertainty`` are used. +* ``default_systematic_smearing``: Default systematic broadening for this category. Units must be eV. +* ``default_systematic_uncertainty``: Default systematic uncertainty for this category. Units must be eV. + +**PlasmaEffects** + +* ``usefixedvalue``: If True ``default_systematic_smearing`` and ``default_systematic_uncertainty`` are used. +* ``default_systematic_smearing``: Default systematic broadening for this category. Units must be eV. +* ``default_systematic_uncertainty``: Default systematic uncertainty for this category. Units must be eV. + + diff --git a/documentation/validation_log.rst b/documentation/validation_log.rst index e148dc91..f408e467 100644 --- a/documentation/validation_log.rst +++ b/documentation/validation_log.rst @@ -8,7 +8,7 @@ Version: v1.2.3 ~~~~~~~~~~~~~~~~ Release Date: Tues July 20 2021 -'''''''''''''''''''''''''''''' +'''''''''''''''''''''''''''''''''' Fixes: ''''''''''''' @@ -112,7 +112,7 @@ Version: v1.1.8 ~~~~~~~~~~~~~~~ Release Date: Thur Apr 18 2019 -''''''''''''''''''''''''''''' +''''''''''''''''''''''''''''''''' New features: ''''''''''''' diff --git a/mermithid/cavity/HannekeFunctions.py b/mermithid/cavity/HannekeFunctions.py new file mode 100644 index 00000000..be185df4 --- /dev/null +++ b/mermithid/cavity/HannekeFunctions.py @@ -0,0 +1,139 @@ +import numpy as np +from scipy import special +from numericalunits import T, kB, hbar, e, me, eV, c0, eps0, m, Hz, MHz + +# Physics constants +endpoint = 18563.251*eV # 30472.604*eV # Krypton +bessel_derivative_zero = special.jnp_zeros(0, 1)[0] + +def magneticfield_from_frequency(cyclotron_frequency, kin_energy=endpoint): + # magnetic field for a given cyclotron frequency + return 2*np.pi*me*cyclotron_frequency/e*gamma(kin_energy) + +def gamma(kin_energy): + return kin_energy/(me*c0**2) + 1 + +def beta_factor(kin_energy=endpoint): + # electron speed at kin_energy + return np.sqrt(kin_energy**2+2*kin_energy*me*c0**2)/(kin_energy+me*c0**2) + +def larmor_radius(magnetic_field, kin_energy=endpoint, pitch=np.pi/2): + transverse_speed = beta_factor(kin_energy)*c0*np.sin(pitch) + # Larmor radius, exact, can also be apporximated by beta*Xprime01*cavity_radius. + return gamma(kin_energy)*me*transverse_speed/(e*magnetic_field) + +# Hanneke factor for TE011 mode, for an electron at fixed location (r_position, z_position). +# Note: We don't "halve the power," given this from Rick: https://3.basecamp.com/3700981/buckets/3107037/uploads/8101664058. +def hanneke_factor_TE011(r_position, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, mode_frequency): + global bessel_derivative_zero + # Calculate the lambda_mnp_squared factor + mode_p = 1 # TE011 mode + z_L = l_cav/2 + # Calculate the lambda_mnp_squared factor + classical_electron_radius_c_squared = e**2 / (4 * np.pi * eps0 * me) + lambda_bessel_part = (1 / (special.jvp(0, bessel_derivative_zero, n=2) * special.jv(0, bessel_derivative_zero))) + constant_factor = - 2 * lambda_bessel_part * classical_electron_radius_c_squared + #constant_factor_unitless = constant_factor/m**3/Hz**2 + #print("Hanneke prefactor [units /m**3/Hz**2]", constant_factor_unitless) + + lambda_mnp_squared = constant_factor / (z_L * r_cav**2) # This should be in units of Hz^2 + angular_part = special.jvp(0, bessel_derivative_zero * r_position/r_cav)**2 + axial_part = np.sin(mode_p * np.pi/2 * (z_position/z_L + 1))**2 + + # Efficienctly deal with 2D scans. + if hasattr(axial_part, "__len__") and hasattr(angular_part, "__len__"): + lambda_mnp_squared *= np.outer(angular_part,axial_part) + else: + lambda_mnp_squared *= (angular_part*axial_part) + # Calculate the damping factor + delta = 2*loaded_Q*( cyclotron_frequency/mode_frequency-1) + return lambda_mnp_squared, delta + +# Calculate the radiated power for an electron at fixed location (r_position, z_position) with energy tranverse_kinetic_energy. +def hanneke_radiated_power(r_position, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, tranverse_kinetic_energy, mode_frequency=None): + if mode_frequency is None: + # Assume that the center of the mode and the cyclotron frequency are identical + mode_frequency = cyclotron_frequency + + lambda_mnp_squared, delta = hanneke_factor_TE011(r_position, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, mode_frequency=mode_frequency) + return tranverse_kinetic_energy*(2*loaded_Q/(1+delta**2))*lambda_mnp_squared/(mode_frequency*np.pi*2) + +# Calculate the radiated power for an electron at fixed location (r_position, z_position) with energy tranverse_kinetic_energy. +# Averaging over the larmor power included. +def larmor_orbit_averaged_hanneke_power(r_position, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, + kinetic_energy=endpoint, pitch=np.pi/2, mode_frequency=None, n_points=100): + if mode_frequency is None: + # Assume that the center of the mode and the cyclotron frequency are identical + mode_frequency = cyclotron_frequency + tranverse_kinetic_energy = kinetic_energy*np.sin(pitch)**2 + + magnetic_field = magneticfield_from_frequency(cyclotron_frequency, kinetic_energy) + electron_orbit_r = larmor_radius(magnetic_field, kin_energy=kinetic_energy, pitch=pitch) + + random_angles = (np.linspace(0,1,n_points)+np.random.rand())*2*np.pi # equally spaced ponts on a circle with random offset + x_random = electron_orbit_r*np.cos(random_angles) + y_random = electron_orbit_r*np.sin(random_angles) + + r_scalar = False + if not hasattr(r_position, "__len__"): + r_scalar = True + r_position = np.array([r_position]) + if hasattr(z_position, "__len__"): + hanneke_powers = np.empty((len(r_position),len(z_position))) + else: + hanneke_powers = np.empty(len(r_position)) + + for i,r_pos_center in enumerate(r_position): + r_pos_orbit = np.sqrt((x_random+r_pos_center)**2 + y_random**2) + # Check for points outside the cavity + if np.any(np.abs(r_pos_orbit) > r_cav): + if hasattr(z_position, "__len__"): + hanneke_powers[i] = np.zeros(len(z_position)) + else: + hanneke_powers[i] = 0 + else: + hanneke_power = hanneke_radiated_power(r_pos_orbit, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, tranverse_kinetic_energy, mode_frequency=mode_frequency) + hanneke_powers[i] = np.mean(hanneke_power.T, axis=-1) + if r_scalar: + hanneke_powers = hanneke_powers[0] + + return hanneke_powers +""" +def larmor_orbit_averaged_hanneke_power(r_position, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, + kinetic_energy=endpoint, pitch=np.pi/2, mode_frequency=None, n_points=100): + if mode_frequency is None: + # Assume that the center of the mode and the cyclotron frequency are identical + mode_frequency = cyclotron_frequency + tranverse_kinetic_energy = kinetic_energy*np.sin(pitch)**2 + + magnetic_field = magneticfield_from_frequency(cyclotron_frequency, kinetic_energy) + electron_orbit_r = larmor_radius(magnetic_field, kin_energy=kinetic_energy, pitch=pitch) + + random_angles = (np.linspace(0,1,n_points)+np.random.rand())*2*np.pi # equally spaced ponts on a circle with random offset + x_random = electron_orbit_r*np.cos(random_angles) + y_random = electron_orbit_r*np.sin(random_angles) + if hasattr(z_position, "__len__"): + hanneke_powers = np.empty((len(r_position),len(z_position))) + else: + hanneke_powers = np.empty(len(r_position)) + + for i,r_pos_center in enumerate(r_position): + r_pos_orbit = np.sqrt((x_random+r_pos_center)**2 + y_random**2) + # Check for points outside the cavity + if np.any(np.abs(r_pos_orbit) > r_cav): + if hasattr(z_position, "__len__"): + hanneke_powers[i] = np.zeros(len(z_position)) + else: + hanneke_powers[i] = 0 + else: + hanneke_power = hanneke_radiated_power(r_pos_orbit, z_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, tranverse_kinetic_energy, mode_frequency=mode_frequency) + hanneke_powers[i] = np.mean(hanneke_power.T, axis=-1) + return hanneke_powers +""" + +# Calculate the average radiated power for an electron with radius r_position in a box trap: +def larmor_orbit_averaged_hanneke_power_box(r_position, loaded_Q, l_cav, r_cav, cyclotron_frequency, + kinetic_energy=endpoint, pitch=np.pi/2, mode_frequency=None, n_points=100): + + return larmor_orbit_averaged_hanneke_power(r_position, 0, loaded_Q, l_cav, r_cav, cyclotron_frequency, + kinetic_energy=kinetic_energy, pitch=pitch, mode_frequency=mode_frequency, n_points=n_points)/2 diff --git a/mermithid/cavity/__init__.py b/mermithid/cavity/__init__.py new file mode 100644 index 00000000..78b46a30 --- /dev/null +++ b/mermithid/cavity/__init__.py @@ -0,0 +1,6 @@ +''' +''' + +from __future__ import absolute_import + +from . import HannekeFunctions \ No newline at end of file diff --git a/mermithid/misc/CRESFunctions_numericalunits.py b/mermithid/misc/CRESFunctions_numericalunits.py new file mode 100644 index 00000000..5f65d3e4 --- /dev/null +++ b/mermithid/misc/CRESFunctions_numericalunits.py @@ -0,0 +1,58 @@ +''' +Miscellaneous functions for CRES conversions +Author: C. Claessens +Date: 05/02/2024 + +Note: The functions here use the numericalunits package. They are used by the Sensitivity calculations. +For CRES functions not using numericalunits import from ConversionFunctions.py. +''' + +from __future__ import absolute_import + +import numpy as np + +from morpho.utilities import morphologging +logger = morphologging.getLogger(__name__) + +from mermithid.misc.Constants_numericalunits import * + + +def gamma(kin_energy): + return kin_energy/(me*c0**2) + 1 + +def beta(kin_energy): + # electron speed at kin_energy + return np.sqrt(kin_energy**2+2*kin_energy*me*c0**2)/(kin_energy+me*c0**2) + +def frequency(kin_energy, magnetic_field): + # cyclotron frequency + return e/(2*np.pi*me)/gamma(kin_energy)*magnetic_field + +def cyclotron_radius(freq, kin_energy): + v = beta(kin_energy)*c0 + return v/freq/(2*np.pi) + +def wavelength(kin_energy, magnetic_field): + return c0/frequency(kin_energy, magnetic_field) + +def kin_energy(freq, magnetic_field): + return (e*c0**2/(2*np.pi*freq)*magnetic_field - me*c0**2) + +def rad_power(kin_energy, pitch, magnetic_field): + # electron radiation power + f = frequency(kin_energy, magnetic_field) + b = beta(kin_energy) + Pe = 2*np.pi*(e*f*b*np.sin(pitch/rad))**2/(3*eps0*c0*(1-b**2)) + return Pe + +def track_length(rho, kin_energy=None, molecular=True): + if kin_energy is None: + kin_energy = tritium_endpoint_molecular if molecular else tritium_endpoint_atomic + crosssect = tritium_electron_crosssection_molecular if molecular else tritium_electron_crosssection_atomic + return 1 / (rho * crosssect * beta(kin_energy) * c0) + +def sin2theta_sq_to_Ue4_sq(sin2theta_sq): + return 0.5*(1-np.sqrt(1-sin2theta_sq)) + +def Ue4_sq_to_sin2theta_sq(Ue4_sq): + return 4*Ue4_sq*(1-Ue4_sq) diff --git a/mermithid/misc/Constants.py b/mermithid/misc/Constants.py index a967fae4..bf60d290 100644 --- a/mermithid/misc/Constants.py +++ b/mermithid/misc/Constants.py @@ -33,7 +33,7 @@ def Vud(): return 0.97425 #CKM element #Beta decay-specific physical constants def QT(): return 18563.251 #For atomic tritium (eV), from Bodine et al. (2015) -def QT2(): return 18573.24 #For molecular tritium (eV), Bodine et al. (2015) +def QT2(): return 18574.01 #For molecular tritium (eV), Bodine et al. (2015) and Meyers et al. Discussion: https://projecteight.slack.com/archives/CG5TY2UE7/p1649963449399179 def Rn(): return 2.8840*10**(-3) #Helium-3 nuclear radius in units of me, from Kleesiek et al. (2018): https://arxiv.org/pdf/1806.00369.pdf def M_3He_in_me(): return 5497.885 #Helium-3 mass in units of me, Kleesiek et al. (2018) def atomic_num(): return 2. #For helium-3 diff --git a/mermithid/misc/Constants_numericalunits.py b/mermithid/misc/Constants_numericalunits.py new file mode 100644 index 00000000..0f8d2ac1 --- /dev/null +++ b/mermithid/misc/Constants_numericalunits.py @@ -0,0 +1,46 @@ +''' +Some constants useful for various things... +The constants here use the numericalunits package. For constants not using this package import form Constants.py +''' + +import numpy as np + +from numericalunits import e, me, c0, eps0, kB, hbar +from numericalunits import meV, eV, keV, MeV, mm, cm, m, ns, s, Hz, kHz, MHz, GHz, amu, nJ +from numericalunits import nT, uT, mT, T, mK, K, C, F, g, W +from numericalunits import hour, year, day, s, ms +from numericalunits import mu0, NA, kB, hbar, me, c0, e, eps0, hPlanck + + +T0 = -273.15*K + +tritium_livetime = 5.605e8*s +tritium_mass_atomic = 3.016* amu *c0**2 +tritium_electron_crosssection_atomic = 9.e-23*m**2 #Hamish extrapolated to 18.6keV using Shah et al. (1987): https://iopscience.iop.org/article/10.1088/0022-3700/20/14/022 +tritium_endpoint_atomic = 18563.251*eV +last_1ev_fraction_atomic = 2.06e-13/eV**3 #From https://arxiv.org/pdf/2012.14341 #Version from Rene; unsure of source: 2.067914e-13/eV**3 + +tritium_mass_molecular = 6.032099 * amu *c0**2 +tritium_electron_crosssection_molecular = 3.67*1e-22*m**2 #[Inelastic from Aseev (2000) for T2] + [Elastic from Liu (1987) for H2, extrapolated by Elise to 18.6keV] +tritium_endpoint_molecular = 18574.01*eV +last_1ev_fraction_molecular = 1.69e-13/eV**3 #From https://arxiv.org/pdf/2012.14341 #Version from Rene; unsure of source: 1.67364e-13/eV**3 + +ground_state_width = 0.436 * eV +ground_state_width_uncertainty = 0.001*0.436*eV + +gyro_mag_ratio_proton = 42.577*MHz/T + +# units that do not show up in numericalunits +# missing pre-factors +fW = W*1e-15 + +# unitless units, relative fractions +pc = 0.01 +ppm = 1e-6 +ppb = 1e-9 +ppt = 1e-12 +ppq = 1e-15 + +# radian and degree which are also not really units +rad = 1 +deg = np.pi/180 \ No newline at end of file diff --git a/mermithid/misc/ConversionFunctions.py b/mermithid/misc/ConversionFunctions.py index c1a59a85..88590af6 100644 --- a/mermithid/misc/ConversionFunctions.py +++ b/mermithid/misc/ConversionFunctions.py @@ -33,3 +33,5 @@ def Energy(F, B=1): else: gamma = (e()*B)/(2.0*np.pi*emass_kg) * 1/(F) return (gamma-1)*m_electron() + + diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py index 5e36fe6e..4bd26fe2 100644 --- a/mermithid/misc/FakeTritiumDataFunctions.py +++ b/mermithid/misc/FakeTritiumDataFunctions.py @@ -140,7 +140,10 @@ def ephasespace(K, Q): Tritium beta spectrum definition """ -#Beta spectrum with a lower energy bound Kmin +#spectral_rate_in_window: Beta spectrum with a lower energy bound Kmin +#This function only is exact for the atomic case, or for a molecular case in the last +# ~18 eV of the spectrum (where molecular final states can be included later as a small +# gaussian broadening). def spectral_rate_in_window(K, Q, mnu, Kmin): if Q-mnu > K > Kmin: return GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K, Q)*(Q - K)*np.sqrt((Q - K)**2 - (mnu)**2) diff --git a/mermithid/misc/__init__.py b/mermithid/misc/__init__.py index afb4ebbe..fa5aaa83 100644 --- a/mermithid/misc/__init__.py +++ b/mermithid/misc/__init__.py @@ -8,3 +8,5 @@ from . import TritiumFormFactor from . import FakeTritiumDataFunctions from . import ConversionFunctions +from . import CRESFunctions_numericalunits +from . import Constants_numericalunits \ No newline at end of file diff --git a/mermithid/processors/Sensitivity/AnalyticSensitivityEstimation.py b/mermithid/processors/Sensitivity/AnalyticSensitivityEstimation.py new file mode 100644 index 00000000..8e129e2f --- /dev/null +++ b/mermithid/processors/Sensitivity/AnalyticSensitivityEstimation.py @@ -0,0 +1,127 @@ +''' +Calculate analytic sensitivity +function. +Author: C. Claessens +Date:12/16/2021 + +More description +''' + +from __future__ import absolute_import + + +import numpy as np + + +# Numericalunits is a package to handle units and some natural constants +# natural constants + +from numericalunits import meV, eV, T + + +# morpho imports +from morpho.utilities import morphologging, reader +from morpho.processors import BaseProcessor +from mermithid.sensitivity.SensitivityFormulas import Sensitivity + + +logger = morphologging.getLogger(__name__) + + + +__all__ = [] +__all__.append(__name__) + +class AnalyticSensitivityEstimation(BaseProcessor, Sensitivity): + ''' + Description + Args: + + Inputs: + + Output: + + ''' + # first do the BaseProcessor __init__ + def __init__(self, name, *args, **kwargs): + BaseProcessor.__init__(self, name, *args, **kwargs) + + def InternalConfigure(self, params): + ''' + Configure + ''' + # file paths + self.config_file_path = reader.read_param(params, 'config_file_path', "required") + + # Configuration is done in __init__ of SensitivityClass + Sensitivity.__init__(self, self.config_file_path) + + return True + + + + def InternalRun(self): + + self.results = {'CL90_limit': self.CL90()/eV, 'm_beta_squared_uncertainty': self.sensitivity()/(eV**2)} + + return True + + + # Sensitivity formulas: + # These are the functions that have so far been in the SensitivityClass but would not be used by a fake data experiment + # I did not include DeltaEWidth here becuase I consider it more general information about an experiment that could be used by the fake data studies too + + def StatSens(self): + """Pure statistic sensitivity assuming Poisson count experiment in a single bin""" + sig_rate = self.SignalRate() + DeltaE = self.DeltaEWidth() + sens = 2/(3*sig_rate*self.Experiment.LiveTime)*np.sqrt(sig_rate*self.Experiment.LiveTime*DeltaE + +self.BackgroundRate()*self.Experiment.LiveTime/DeltaE) + return sens + + def SystSens(self): + """Pure systematic componenet to sensitivity""" + labels, sigmas, deltas = self.get_systematics() + sens = 4*np.sqrt(np.sum((sigmas*deltas)**2)) + return sens + + def sensitivity(self, **kwargs): + """Combined statisical and systematic uncertainty. + Using kwargs settings in namespaces can be changed. + Example how to change number density which lives in namespace Experiment: + self.sensitivity(Experiment={"number_density": rho}) + """ + for sect, options in kwargs.items(): + for opt, val in options.items(): + self.__dict__[sect].__dict__[opt] = val + + StatSens = self.StatSens() + SystSens = self.SystSens() + + # Standard deviation on a measurement of m_beta**2 + sigma_m_beta_2 = np.sqrt(StatSens**2 + SystSens**2) + return sigma_m_beta_2 + + def CL90(self, **kwargs): + """ Gives 90% CL upper limit on neutrino mass.""" + # If integrate gaussian -infinity to 1.28*sigma, the result is 90%. + # https://www.wolframalpha.com/input?i=1%2F2+%281+%2B+erf%281.281551%2Fsqrt%282%29%29%29 + # So we assume that 1.28*sigma_{m_beta^2} is the 90% upper limit on m_beta^2. + # Note: 1.64 update, to prevent gaining from a sensitivity from a negative mass result + ### We now use 1.64 = 1.28**2 + return np.sqrt(1.281551**2*self.sensitivity(**kwargs)) + + def sterial_m2_limit(self, Ue4_sq): + return np.sqrt(1.281551**2*np.sqrt((self.StatSens()/Ue4_sq)**2 + self.SystSens()**2)) + + + # print functions + def print_statistics(self): + print("Statistic", " "*18, "%.2f"%(np.sqrt(self.StatSens())/meV), "meV") + + def print_systematics(self): + labels, sigmas, deltas = self.get_systematics() + + print() + for label, sigma, delta in zip(labels, sigmas, deltas): + print(label, " "*(np.max([len(l) for l in labels])-len(label)), "%8.2f"%(sigma/meV), "+/-", "%8.2f"%(delta/meV), "meV") diff --git a/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py b/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py new file mode 100644 index 00000000..091d6958 --- /dev/null +++ b/mermithid/processors/Sensitivity/CavitySensitivityCurveProcessor.py @@ -0,0 +1,1123 @@ +''' +Calculate sensitivity curve and plot vs. number density, exposure, livetime, or frequency. + +Author: C. Claessens, T. E. Weiss +Date: 06/07/2023 +Updated: April 2025 +''' + +from __future__ import absolute_import + + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from copy import deepcopy + + +# Numericalunits is a package to handle units and some natural constants +# natural constants +from numericalunits import e, me, c0, eps0, kB, hbar +from numericalunits import meV, eV, keV, MeV, cm, m, mm +from numericalunits import nT, uT, mT, T, mK, K, C, F, g, W +from numericalunits import hour, year, day, ms, ns, s, Hz, kHz, MHz, GHz +ppm = 1e-6 +ppb = 1e-9 + +# morpho imports +from morpho.utilities import morphologging, reader +from morpho.processors import BaseProcessor +from mermithid.sensitivity.SensitivityCavityFormulas import CavitySensitivity + + +logger = morphologging.getLogger(__name__) + + + +__all__ = [] +__all__.append(__name__) + +class CavitySensitivityCurveProcessor(BaseProcessor): + ''' + Description + Args: + + Inputs: + + Output: + + ''' + def InternalConfigure(self, params): + ''' + Configure + ''' + # file paths + self.config_file_path = reader.read_param(params, 'config_file_path', "required") + self.comparison_config_file_path = reader.read_param(params, 'comparison_config_file_path', '') + self.plot_path = reader.read_param(params, 'plot_path', "required") + self.PhaseII_path = reader.read_param(params, 'PhaseII_config_path', '') + + # labels + self.main_curve_upper_label = reader.read_param(params, 'main_curve_upper_label', r"molecular"+"\n"+r"$V_\mathrm{eff} = 2\, \mathrm{cm}^3$"+"\n"+r"$\sigma_B = 7\,\mathrm{ppm}$") + self.main_curve_lower_label = reader.read_param(params, 'main_curve_lower_label', r"$\sigma_B = 1\,\mathrm{ppm}$") + self.comparison_curve_label = reader.read_param(params, 'comparison_curve_label', r"atomic"+"\n"+r"$V_\mathrm{eff} = 5\, \mathrm{m}^3$"+"\n"+r"$\sigma_B = 0.13\,\mathrm{ppm}$") + self.main_curve_color = reader.read_param(params, 'main_curve_color', "darkblue") + self.comparison_curve_colors = reader.read_param(params,'comparison_curve_colors', ["blue", "darkred", "red"]) + self.main_curve_linestyle = reader.read_param(params, 'main_curve_linestyle', "solid") + self.comparison_curve_linestyles = reader.read_param(params,'comparison_curve_linestyles', ["solid", "solid", "solid"]) + self.main_curve_marker = reader.read_param(params, 'main_curve_marker', "d") + self.comparison_curve_markers = reader.read_param(params,'comparison_curve_markers', ["d", "d", "d"]) + + # options + self.optimize_main_density = reader.read_param(params, 'optimize_main_density', True) + self.optimize_comparison_density = reader.read_param(params, 'optimize_comparison_density', True) + self.verbose = reader.read_param(params, 'verbose', True) + self.comparison_curve = reader.read_param(params, 'comparison_curve', False) + self.B_error = reader.read_param(params, 'B_inhomogeneity', 7e-6) + self.B_error_uncertainty = reader.read_param(params, 'B_inhom_uncertainty', 0.05) + self.sigmae_theta_r = reader.read_param(params, 'sigmae_theta_r', 0.159) #eV + self.configure_sigma_theta_r = reader.read_param(params, "configure_sigma_theta_r", False) + + # main plot configurations + self.figsize = reader.read_param(params, 'figsize', (6,6)) + self.legend_location = reader.read_param(params, 'legend_location', 'upper left') + self.legend_bbox_to_anchor = reader.read_param(params, 'legend_bbox_to_anchor', (0.15,0,1.1,0.87)) + self.fontsize = reader.read_param(params, 'fontsize', 12) + + self.density_axis = reader.read_param(params, "density_axis", True) + self.frequency_axis = reader.read_param(params, "frequency_axis", False) + self.exposure_axis = reader.read_param(params, "exposure_axis", False) + self.livetime_axis = reader.read_param(params, "livetime_axis", False) + self.ncavities_livetime_axis = reader.read_param(params, "ncavities_livetime_axis", False) + self.ncav_eff_time_axis = reader.read_param(params, "ncav_eff_time_axis", False) + self.track_length_axis = reader.read_param(params, 'track_length_axis', True) + self.atomic_axis = reader.read_param(params, 'atomic_axis', False) + self.molecular_axis = reader.read_param(params, 'molecular_axis', False) + self.magnetic_field_axis = reader.read_param(params, 'magnetic_field_axis', False) + + self.density_range = reader.read_param(params, 'density_range', [1e14,1e21]) + self.year_range = reader.read_param(params, "year_range", [0.1, 20]) + self.exposure_range = reader.read_param(params, "exposure_range", [1e-10, 1e3]) + self.frequency_range = reader.read_param(params, "frequency_range", [1e6, 1e9]) + self.det_thresh_range = reader.read_param(params, 'det_thresh_range', [1, 160]) + + self.ylim = reader.read_param(params, 'y_limits', [1e-2, 1e2]) + + self.label_x_position = reader.read_param(params, 'label_x_position', 5e19) + self.upper_label_y_position = reader.read_param(params, 'upper_label_y_position', 5) + self.lower_label_y_position = reader.read_param(params, 'lower_label_y_position', 2.3) + self.goal_x_pos = reader.read_param(params, "goals_x_position", {"LFA threshold (0.7 eV)": 0.11, "Phase IV (0.04 eV)": 0.11}) + self.goals_y_rel_position = reader.read_param(params, "goals_y_rel_position", {"LFA threshold (0.7 eV)": 0.855, "Phase IV (0.04 eV)": 0.855}) + + self.comparison_label_y_position = reader.read_param(params, 'comparison_label_y_position', 0.044) + self.comparison_label_x_position = reader.read_param(params, 'comparison_label_x_position', 5e16) + if self.comparison_label_x_position == None: + self.comparison_label_x_position = reader.read_param(params, 'label_x_position', 5e16) + + self.add_PhaseII = reader.read_param(params, "add_PhaseII", False) + self.add_point_at_configured_density = reader.read_param(params, "add_point_at_configured_density", True) + self.add_1year_1cav_point_to_last_ref = reader.read_param(params, "add_1year_1cav_point_to_last_ref", False) + + # key parameter plots + self.make_key_parameter_plots = reader.read_param(params, 'plot_key_parameters', False) + + if self.density_axis: + self.add_sens_line = self.add_density_sens_line + logger.info("Plotting sensitivity vs. density") + elif self.frequency_axis: + self.add_sens_line = self.add_frequency_sens_line + logger.info("Plotting sensitivity vs. frequency") + if self.optimize_main_density==False or self.optimize_comparison_density==False: + logger.warning("You told me not to optimize the density, but I am doing so anyway! We need to fix this.") + elif self.exposure_axis or self.livetime_axis or self.ncavities_livetime_axis or self.ncav_eff_time_axis: + self.add_sens_line = self.add_exposure_sens_line + logger.info("Plotting sensitivity vs. some exposure (e.g., livetime or ncavities*efficiency*livetime)") + else: raise ValueError("No axis specified") + + + # goals + self.goals = reader.read_param(params, "goals", {}) + + + # setup sensitivities + if self.add_PhaseII: + self.sens_PhaseII = CavitySensitivity(self.PhaseII_path, verbose=self.verbose) + + + self.sens_main = CavitySensitivity(self.config_file_path, verbose=self.verbose) + self.sens_main_is_atomic = self.sens_main.Experiment.atomic + + if self.comparison_curve: + ref = [] + for file in self.comparison_config_file_path: + ref.append(CavitySensitivity(file, verbose=self.verbose)) + + self.sens_ref = ref + is_atomic = [] + for i in range(len(self.sens_ref)): + is_atomic.append(self.sens_ref[i].Experiment.atomic) + self.sens_ref_is_atomic = is_atomic + else: + self.sens_ref_is_atomic = [False] + + # check atomic and molecular + if self.molecular_axis: + if not self.sens_main_is_atomic: + #self.molecular_sens = self.sens_main + logger.info("Main curve is molecular") + elif False not in self.sens_ref_is_atomic: + logger.warn("No experiment is configured to be molecular") + #self.molecular_sens = self.sens_ref + #logger.info("Comparison curve is molecular") + #else: + + # Setting natoms_per_particle for atomic and molecular cases, to use later + # when calculating and printing event rates + if self.sens_main_is_atomic: + self.sens_main_natoms_per_particle = 1 + else: + self.sens_main_natoms_per_particle = 2 + if self.comparison_curve: + for i in range(len(self.sens_ref)): + if self.sens_ref_is_atomic[i]: + self.sens_ref[i].natoms_per_particle = 1 + else: + self.sens_ref[i].natoms_per_particle = 2 + + if self.atomic_axis: + if self.sens_main_is_atomic: + #self.atomic_sens = self.sens_main + logger.info("Main curve is atomic") + elif True in self.sens_ref_is_atomic: + #self.atomic_sens = self.sens_ref + logger.info("A comparison curve is atomic") + else: + logger.warn("No experiment is configured to be atomic") + + # densities, exposures, runtimes + self.rhos = np.logspace(np.log10(self.density_range[0]), np.log10(self.density_range[1]), 80)/m**3 + self.exposures = np.logspace(np.log10(self.exposure_range[0]), np.log10(self.exposure_range[1]), 100)*m**3*year + self.years = np.logspace(np.log10(self.year_range[0]), np.log10(self.year_range[1]), 100)*year + self.frequencies = np.logspace(np.log10(self.frequency_range[0]), np.log10(self.frequency_range[1]), 20)*Hz + self.thresholds = np.linspace(self.det_thresh_range[0], self.det_thresh_range[1], 40) + + ncav_temp = self.sens_main.Experiment.n_cavities + self.ncavities_livetime_range = [ncav_temp*self.year_range[0], ncav_temp*self.year_range[1]] + self.ncav_eff_time_range = [ncav_temp*self.year_range[0]*0.01, ncav_temp*self.year_range[1]*0.1] + self.ncav_eff_times = np.logspace(np.log10(self.ncav_eff_time_range[0]), np.log10(self.ncav_eff_time_range[1]), 100)*year + + return True + + + + def InternalRun(self): + + #Calculating the sensitivity for parameters in the main config file + self.sens_with_configured_density_and_thresh = self.sens_main.CL90() + + #Optimizing the detection threshold for the main config file + #Before the density optimization + self.track_duration_before_opt = self.sens_main.track_length(self.sens_main.Experiment.number_density) + #threshs_before_opt = np.linspace(0.5*SNR_track_duration, 1.0*SNR_track_duration, 10) + thresh_limits = [self.sens_main.CL90(Threshold={"detection_threshold": th}) for th in self.thresholds] + index = np.argmin(thresh_limits) + self.thresh_opt_before_n_opt = self.thresholds[index] + self.sens_with_configured_density_and_opt_thresh = thresh_limits[index] + self.sens_main.Threshold.detection_threshold = self.thresh_opt_before_n_opt + + logger.info("Sensitivity before density or threshold optimization: {} eV".format(self.sens_with_configured_density_and_thresh/eV)) + logger.info("Sensitivity with optimized threshold, before density opt: {} eV".format(self.sens_with_configured_density_and_opt_thresh/eV)) + logger.info("Optimum threshold: {}".format(self.thresh_opt_before_n_opt)) + logger.info("Scanned thresholds: {}".format(self.thresholds)) + logger.info("Systematics before density optimization and after threshold optimization:") + self.sens_main.print_systematics() + logger.info("Number density before optimization: {} \m^3".format(self.sens_main.Experiment.number_density*m**3)) + logger.info("Corresponding track length before density optimization: {} s".format(self.track_duration_before_opt/s)) + logger.info("***Efficiencies and bkgd before density optimization, after threshold opt:***") + self.sens_main.print_Efficiencies() + self.sens_main.assign_background_rate_from_threshold() + self.sens_main.BackgroundRate() + logger.info('RF background: {}/eV/s'.format(self.sens_main.RF_background_rate_per_eV*eV*s)) + logger.info('Total background: {}/eV/s'.format(self.sens_main.background_rate*eV*s)) + logger.info("***Done printing pre-optimization***") + + + #Optimizing the detection threshold for the comparison config files + #Before the density optimization + if self.comparison_curve: + for i in range(len(self.sens_ref)): + thresh_limits = [self.sens_ref[i].CL90(Threshold={"detection_threshold": th}) for th in self.thresholds] + self.sens_ref[i].sens_with_configured_density_and_opt_thresh = np.min(thresh_limits) + + + # create main plot + self.create_plot() + + # optionally add Phase II curve and point to exposure plot + if self.exposure_axis and self.add_PhaseII: + self.add_Phase_II_exposure_sens_line(self.sens_PhaseII) + + # add second and third x axis for track lengths + if self.density_axis and self.track_length_axis: + self.add_track_length_axis() + + # add magnetic field axis if frequency axis + if self.frequency_axis and self.magnetic_field_axis: + self.add_magnetic_field_axis() + + # add goals + for key, value in self.goals.items(): + logger.info('Adding goal: {} = {}'.format(key, value)) + self.add_goal(value, key, self.goal_x_pos[key], self.goals_y_rel_position[key]) + + #Add point at configured density + if self.density_axis and self.add_point_at_configured_density: + self.ax.scatter([self.sens_main.Experiment.number_density*m**3], [self.sens_with_configured_density_and_opt_thresh/eV], marker="s", s=25, color=self.main_curve_color, label="Operating density", zorder=4) #label="Density: {:.{}f}".format(self.Experiment.number_density*m**3, 1) + if self.comparison_curve: + for i in range(len(self.sens_ref)): + self.ax.scatter([self.sens_ref[i].Experiment.number_density*m**3], [self.sens_ref[i].sens_with_configured_density_and_opt_thresh/eV], marker="s", s=25, color=self.comparison_curve_colors[i], zorder=4) + + # optimize density + if self.optimize_main_density: + thresh_opt_main = [] + limits = [] + for rho in self.rhos: + thresh_limits = [] + for thresh in self.thresholds: + self.sens_main.Threshold.detection_threshold = thresh + thresh_limits.append(self.sens_main.CL90(Experiment={"number_density": rho})) + index = np.argmin(thresh_limits) + thresh_opt_main.append(self.thresholds[index]) + limits.append(thresh_limits[index]) + #limit = [self.sens_main.CL90(Experiment={"number_density": rho})/eV for rho in self.rhos] + opt_index = np.argmin(limits) + rho_opt = self.rhos[opt_index] + self.sens_main.Experiment.number_density = rho_opt + self.sens_main.Threshold.detection_threshold = thresh_opt_main[opt_index] + logger.info("Optimized thresholds at each density: {}".format(thresh_opt_main)) + logger.info("Optimized threshold at optimum n: {}".format(thresh_opt_main[opt_index])) + + else: + #Use the optimized detection threshold even if the density isn't optimized + self.sens_main.Threshold.detection_threshold = self.thresh_opt_before_n_opt + + + # if B is list plot line for each B + if self.configure_sigma_theta_r and (isinstance(self.sigmae_theta_r, list) or isinstance(self.sigmae_theta_r, np.ndarray)): + N = len(self.sigmae_theta_r) + for a, color in self.range(0, N): + #sig = self.sens_main.BToKeErr(self.sens_main.MagneticField.nominal_field*self.B_error[a], self.sens_main.MagneticField.nominal_field) + #self.sens_main.MagneticField.usefixedvalue = True + #self.sens_main.MagneticField.default_systematic_smearing = sig + #self.sens_main.MagneticField.default_systematic_uncertainty = 0.05*sig + self.sens_main.MagneticField.sigmae_r = self.sigmae_theta_r[a] * eV + self.sens_main.MagneticField.sigmae_theta = 0 * eV + if self.livetime_axis: + self.add_sens_line(self.sens_main, exposure_type="livetime", color=color) + elif self.ncavities_livetime_axis: + self.add_sens_line(self.sens_main, exposure_type="ncavities_livetime", color=color) + elif self.ncav_eff_time_axis: + self.add_sens_line(self.sens_main, exposure_type="ncav_eff_time", color=color) + else: + self.add_sens_line(self.sens_main, color=color) + #print("sigmae_theta_r:", self.sens_main.MagneticField.sigmae_r/eV) + self.sens_main.print_systematics() + self.add_text(self.label_x_position, self.upper_label_y_position, self.main_curve_upper_label, color="darkblue") + self.add_text(self.label_x_position, self.lower_label_y_position, self.main_curve_lower_label, color="darkred") + + else: + if self.configure_sigma_theta_r: + self.sens_main.MagneticField.sigmaer = self.sigmae_theta_r * eV + self.sens_main.MagneticField.sigmae_theta = 0 * eV + if self.livetime_axis: + self.add_sens_line(self.sens_main, exposure_type="livetime", color=self.main_curve_color, marker=self.main_curve_marker, linestyle=self.main_curve_linestyle, label=self.main_curve_upper_label) + elif self.ncavities_livetime_axis: + self.add_sens_line(self.sens_main, exposure_type="ncavities_livetime", color=self.main_curve_color, marker=self.main_curve_marker, linestyle=self.main_curve_linestyle, label=self.main_curve_upper_label) + elif self.ncav_eff_time_axis: + self.add_sens_line(self.sens_main, exposure_type="ncav_eff_time", color=self.main_curve_color, marker=self.main_curve_marker, linestyle=self.main_curve_linestyle, label=self.main_curve_upper_label) + else: + self.add_sens_line(self.sens_main, color=self.main_curve_color, linestyle=self.main_curve_linestyle, label=self.main_curve_upper_label, zorder=3) + + # PRINT OPTIMUM RESULTS + + if self.optimize_main_density: + #This is done above - should confirm whether it actually needs to be done again, or not + self.sens_main.Experiment.number_density = rho_opt + self.sens_main.Threshold.detection_threshold = thresh_opt_main[opt_index] + self.sens_main.EffectiveVolume() + #MAY NEED TO DO THE ABOVE FOR THE SET DENSITY IN THE CONFIG FILE, WHEN NOT OPTIMIZING MAIN DENSITY + + # if the magnetic field uncertainties were configured above, set them back to the first value in the list + if self.configure_sigma_theta_r and (isinstance(self.sigmae_theta_r, list) or isinstance(self.sigmae_theta_r, np.ndarray)): + self.sens_main.MagneticField.sigmae_r = self.sigmae_theta_r[0] * eV + self.sens_main.MagneticField.sigmae_theta = 0 * eV + elif self.configure_sigma_theta_r: + self.sens_main.MagneticField.sigmaer = self.sigmae_theta_r * eV + self.sens_main.MagneticField.sigmae_theta = 0 * eV + + if self.optimize_main_density: + logger.info('MAIN CURVE at optimum density:') + rho = rho_opt + else: + logger.info('MAIN CURVE at configured density (not optimum):') + rho = self.sens_main.Experiment.number_density + logger.info('veff = {} m**3, rho = {} /m**3:'.format(self.sens_main.effective_volume/(m**3), rho*(m**3))) + logger.info("Loaded Q: {}".format(self.sens_main.loaded_q)) + logger.info("Axial frequency for minimum detectable angle: {} MHz".format(self.sens_main.required_bw_axialfrequency/MHz)) + logger.info("Total bandwidth: {} MHz".format(self.sens_main.required_bw/MHz)) + logger.info('Larmor power = {} W, Hanneke power = {} W'.format(self.sens_main.larmor_power/W, self.sens_main.signal_power/W)) + logger.info('Hanneke / Larmor power = {}'.format(self.sens_main.signal_power/self.sens_main.larmor_power)) + logger.info("eta = slope*track_length/(2*omega_c) = {}".format(self.sens_main.eta)) + + if self.sens_main.FrequencyExtraction.crlb_on_sidebands: + logger.info("Trap p: {}".format(self.sens_main.p)) + logger.info("Trap q: {}".format(self.sens_main.q)) + #print(self.sens_main.q_array) + logger.info("Uncertainty from determination of f_carrier and f_lsb, due to noise: {} eV".format(self.sens_main.sigma_K_noise/eV)) + + self.sens_main.print_SNRs(rho) + self.sens_main.print_Efficiencies() + + if self.exposure_axis or self.livetime_axis or self.ncavities_livetime_axis or self.ncav_eff_time_axis: + logger.info("NUMBERS BELOW ARE FOR THE HIGHEST-EXPOSURE POINT ON THE CURVE:") + logger.info('CL90 limit: {} eV'.format(self.sens_main.CL90(Experiment={"number_density": rho})/eV)) + logger.info('Tritium in Veff: {}'.format(rho*self.sens_main.effective_volume)) + self.sens_main.assign_background_rate_from_threshold() + self.sens_main.BackgroundRate() + logger.info('RF background: {}/eV/s'.format(self.sens_main.RF_background_rate_per_eV*eV*s)) + logger.info('Total background: {}/eV/s'.format(self.sens_main.background_rate*eV*s)) + logger.info('Total signal: {}'.format(rho*self.sens_main.effective_volume* + self.sens_main.Experiment.LiveTime/ + self.sens_main.tau_tritium*self.sens_main_natoms_per_particle)) + logger.info('Signal in last eV: {}'.format(self.sens_main.last_1ev_fraction*eV**3* + rho*self.sens_main.effective_volume* + self.sens_main.Experiment.LiveTime/ + self.sens_main.tau_tritium*self.sens_main_natoms_per_particle)) + + self.sens_main.print_statistics() + self.sens_main.print_systematics() + + if self.make_key_parameter_plots: + if not self.optimize_main_density: + logger.info("Set optimize_main_density to True to make key parameter plots") + else: + # First key parameter plot: Stat and Syst vs. density + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = [], [], [] + for i in range(len(self.rhos)): + temp_rho = deepcopy(self.sens_main.Experiment.number_density) + self.sens_main.Experiment.number_density = self.rhos[i] + self.sens_main.Threshold.detection_threshold = thresh_opt_main[i] + labels, sigmas, deltas = self.sens_main.get_systematics() + sigma_startf.append(sigmas[1]) + stat_on_mbeta2.append(self.sens_main.StatSens()) + syst_on_mbeta2.append(self.sens_main.SystSens()) + self.sens_main.Experiment.number_density = temp_rho + + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = np.array(sigma_startf), np.array(stat_on_mbeta2), np.array(syst_on_mbeta2) + fig = plt.figure() + plt.loglog(self.rhos*m**3, stat_on_mbeta2/eV**2, label='Statistical uncertainty') + plt.loglog(self.rhos*m**3, syst_on_mbeta2/eV**2, label='Systematic uncertainty') + plt.xlabel(r"Number density $n\, \, (\mathrm{m}^{-3})$") + plt.ylabel(r"Standard deviation in $m_\beta^2$ (eV$^2$)") + plt.legend() + plt.tight_layout() + plt.savefig("stat_and_syst_vs_density_LFA_threshold.png", dpi=300) + + fig = plt.figure() + plt.loglog(self.rhos*m**3, sigma_startf/eV) + plt.xlabel(r"Number density $n\, \, (\mathrm{m}^{-3})$") + plt.ylabel(r"Resolution from $f$ reconstruction, axial field (eV)") + plt.tight_layout() + plt.savefig("resolution_from_CRLB_vs_density_LFA_threshold.png", dpi=300) + """ + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = [], [], [] + for i in range(len(self.rhos)): + temp_rho = deepcopy(self.sens_ref[1].Experiment.number_density) + self.sens_ref[1].Experiment.number_density = self.rhos[i] + self.sens_ref[1].Threshold.detection_threshold = thresh_opt_main[i] + labels, sigmas, deltas = self.sens_ref[1].get_systematics() + sigma_startf.append(sigmas[1]) + stat_on_mbeta2.append(self.sens_ref[1].StatSens()) + syst_on_mbeta2.append(self.sens_ref[1].SystSens()) + self.sens_ref[1].Experiment.number_density = temp_rho + + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = np.array(sigma_startf), np.array(stat_on_mbeta2), np.array(syst_on_mbeta2) + fig = plt.figure() + plt.loglog(self.rhos*m**3, stat_on_mbeta2/eV**2, label='Statistical uncertainty') + plt.loglog(self.rhos*m**3, syst_on_mbeta2/eV**2, label='Systematic uncertainty') + plt.xlabel(r"Number density $n\, \, (\mathrm{m}^{-3})$") + plt.ylabel(r"Standard deviation in $m_\beta^2$ (eV$^2$)") + plt.legend() + plt.tight_layout() + plt.savefig("stat_and_syst_vs_density_module1ofPhaseIV.png") + + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = [], [], [] + for i in range(len(self.rhos)): + temp_rho = deepcopy(self.sens_ref[2].Experiment.number_density) + self.sens_ref[2].Experiment.number_density = self.rhos[i] + self.sens_ref[2].Threshold.detection_threshold = thresh_opt_main[i] + labels, sigmas, deltas = self.sens_ref[2].get_systematics() + sigma_startf.append(sigmas[1]) + stat_on_mbeta2.append(self.sens_ref[2].StatSens()) + syst_on_mbeta2.append(self.sens_ref[2].SystSens()) + self.sens_ref[2].Experiment.number_density = temp_rho + + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = np.array(sigma_startf), np.array(stat_on_mbeta2), np.array(syst_on_mbeta2) + fig = plt.figure() + plt.loglog(self.rhos*m**3, stat_on_mbeta2/eV**2, label='Statistical uncertainty') + plt.loglog(self.rhos*m**3, syst_on_mbeta2/eV**2, label='Systematic uncertainty') + plt.xlabel(r"Number density $n\, \, (\mathrm{m}^{-3})$") + plt.ylabel(r"Standard deviation in $m_\beta^2$ (eV$^2$)") + plt.legend() + plt.tight_layout() + plt.savefig("stat_and_syst_vs_density_10PhaseIVcavities.png") + """ + + # Optimize comparison curves over density + if self.comparison_curve: + for i in range(len(self.sens_ref)): + if self.optimize_comparison_density: + thresh_opt_comparison = [] + limits_ref = [] + for rho in self.rhos: + thresh_limits = [] + for thresh in self.thresholds: + self.sens_ref[i].Threshold.detection_threshold = thresh + thresh_limits.append(self.sens_ref[i].CL90(Experiment={"number_density": rho})) + index = np.argmin(thresh_limits) + thresh_opt_comparison.append(self.thresholds[index]) + limits_ref.append(thresh_limits[index]) + #limit_ref = [self.sens_ref[i].CL90(Experiment={"number_density": rho})/eV for rho in self.rhos] + limit2_index = np.argmin(limits_ref) + rho_opt_ref = self.rhos[limit2_index] + limit2 = limits_ref[limit2_index] #self.sens_ref[i].CL90(Experiment={"number_density": rho_opt_ref}) + logger.info("Optimized thresholds at each density: {}".format(thresh_opt_comparison)) + logger.info("Optimized threshold at optimum n: {}".format(thresh_opt_comparison[limit2_index])) + else: + thresh_limits = [self.sens_ref[i].CL90(Threshold={"detection_threshold": th}) for th in self.thresholds] + index2 = np.argmin(thresh_limits) + self.sens_ref[i].Threshold.detection_threshold = self.thresholds[index2] + limit2 = thresh_limits[index2] + + if self.optimize_comparison_density: + logger.info('COMPARISON CURVE at optimum density:') + rho2 = rho_opt_ref + self.sens_ref[i].Experiment.number_density = rho2 + self.sens_ref[i].Threshold.detection_threshold = thresh_opt_comparison[limit2_index] + self.sens_ref[i].EffectiveVolume() + else: + logger.info('COMPARISON CURVE at configured density (not optimum):') + rho2 = self.sens_ref[i].Experiment.number_density + logger.info("Optimized threshold at configured density: {}".format(self.thresholds[index2])) + + logger.info('veff = {} m**3, rho = {} /m**3:'.format(self.sens_ref[i].effective_volume/(m**3), rho2*(m**3))) + logger.info("Loaded Q: {}".format(self.sens_ref[i].loaded_q)) + logger.info('Larmor power = {} W, Hanneke power = {} W'.format(self.sens_ref[i].larmor_power/W, self.sens_ref[i].signal_power/W)) + logger.info('Hanneke / Larmor power = {}'.format(self.sens_ref[i].signal_power/self.sens_ref[i].larmor_power)) + + if self.sens_ref[i].FrequencyExtraction.crlb_on_sidebands: + logger.info("Uncertainty from determination of f_carrier and f_lsb, due to noise: {} eV".format(self.sens_ref[i].sigma_K_noise/eV)) + + self.sens_ref[i].print_SNRs(rho2) + self.sens_ref[i].print_Efficiencies() + #if self.exposure_axis or self.livetime_axis or self.ncavities_livetime_axis or self.ncav_eff_time_axis: + # logger.info("NUMBERS BELOW ARE FOR THE HIGHEST-EXPOSURE POINT ON THE CURVE:") + logger.info('CL90 limit: {} eV'.format(limit2/eV)) + logger.info('Tritium in Veff: {}'.format(rho2*self.sens_ref[i].effective_volume)) + self.sens_ref[i].assign_background_rate_from_threshold() + self.sens_ref[i].BackgroundRate() + logger.info('RF background: {}/eV/s'.format(self.sens_ref[i].RF_background_rate_per_eV*eV*s)) + logger.info('Total background: {}/eV/s'.format(self.sens_ref[i].background_rate*eV*s)) + logger.info('Total signal: {}'.format(rho2*self.sens_ref[i].effective_volume* + self.sens_ref[i].Experiment.LiveTime/ + self.sens_ref[i].tau_tritium*self.sens_ref[i].natoms_per_particle)) + logger.info('Signal in last eV: {}'.format(self.sens_ref[i].last_1ev_fraction*eV**3* + rho2*self.sens_ref[i].effective_volume* + self.sens_ref[i].Experiment.LiveTime/ + self.sens_ref[i].tau_tritium*self.sens_ref[i].natoms_per_particle)) + + self.sens_ref[i].print_statistics() + self.sens_ref[i].print_systematics() + + #Add comparison curves + label=self.comparison_curve_label + if self.livetime_axis: + limits = self.add_sens_line(self.sens_ref[i], exposure_type="livetime", plot_key_params=False, color=self.comparison_curve_colors[i], marker=self.comparison_curve_markers[i], linestyle=self.comparison_curve_linestyles[i], label=label[i]) + elif self.ncavities_livetime_axis: + limits = self.add_sens_line(self.sens_ref[i], exposure_type="ncavities_livetime", plot_key_params=False, color=self.comparison_curve_colors[i], marker=self.comparison_curve_markers[i], linestyle=self.comparison_curve_linestyles[i], label=label[i]) + elif self.ncav_eff_time_axis: + limits = self.add_sens_line(self.sens_ref[i], exposure_type="ncav_eff_time", plot_key_params=False, color=self.comparison_curve_colors[i], marker=self.comparison_curve_markers[i], linestyle=self.comparison_curve_linestyles[i], label=label[i]) + else: + limits = self.add_sens_line(self.sens_ref[i], plot_key_params=False, color=self.comparison_curve_colors[i], linestyle=self.comparison_curve_linestyles[i], label=label[i], zorder=5) + #self.ax.text(self.comparison_label_x_position[a], self.comparison_label_y_position[a], label[a], color=colors[a], fontsize=9.5) + + if self.exposure_axis or self.livetime_axis or self.ncavities_livetime_axis or self.ncav_eff_time_axis: + if self.add_1year_1cav_point_to_last_ref: + logger.info("Going to add single exposure point") + self.add_single_exposure_point(self.sens_ref[-1], color=self.comparison_curve_colors[-1]) + + + #self.add_comparison_curve(label=self.comparison_curve_label) + #self.add_arrow(self.sens_main) + + # save plot + self.save(self.plot_path) + + if self.verbose: + self.print_disclaimers() + + return True + + + def print_disclaimers(self): + logger.info("Disclaimers / assumptions:") + logger.info("1. Often, in practice, 'sigmae_r' is used to stand-in for combined radial, \ + azimuthal, and temporal magnetic field spectral broadening effects. Check \ + your experiment config file *and* your processor config dictionary to see \ + if that is the case.") + logger.info("2. Trap design is not yet linked to cavity L/D in the sensitivity model. So, \ + the model does *not* capture how reducing L/D worsens the resolution.") + logger.info("3. In reality, the frequency resolution could be worse or somewhat better \ + than predicted by the chirp CRLB calculation used here. See work by Florian.") + logger.info("4. The analytic sensitivity formula accounts for energy resolution contributions \ + that are *normally distributed*. (Energy resolution = std of the response fn \ + that broadens the spectrum.) To account for asymmetric contributions, generate \ + spectra with MC sampling and then analyze them. This can be done in mermithid.") + logger.info("5. The best-fit mbeta is assumed to be zero when converting to a 90\% limit.") + if self.density_axis: + logger.info("6. This sensitivity formula does not work for very small numbers of counts, \ + because the analytic formula assumes Gaussian statistics. In typical Phase IV \ + scenarios, if the minimum allowed density is 1e-20 atoms/m^3, the optimization \ + over density still works.") + logger.info("Once you have read these disclaimers and are familiar with them, you can set \ + verbose==False in your config dictionary to stop seeing them.") + + + def create_plot(self): + # setup axis + plt.rcParams.update({'font.size': self.fontsize}) + self.fig, self.ax = plt.subplots(figsize=self.figsize) + ax = self.ax + ax.set_xscale("log") + ax.set_yscale("log") + if self.density_axis: + logger.info("Adding density axis") + ax.set_xlim(self.rhos[0]*m**3, self.rhos[-1]*m**3) + + if self.atomic_axis and self.molecular_axis: + axis_label = r"(Atomic / molecular) number density $n\, \, (\mathrm{m}^{-3})$" + elif self.atomic_axis: + axis_label = r"Average atom number density in sensitive volume$\, \, (\mathrm{m}^{-3})$" + elif self.molecular_axis: + axis_label = r"Molecular number density $n\, \, (\mathrm{m}^{-3})$" + else: + axis_label = r"Number density $n\, \, (\mathrm{m}^{-3})$" + + ax.set_xlabel(axis_label) + ax.set_ylim(self.ylim) + ax.set_ylabel(r"90% CL on $m_\beta$ (eV)") + + + elif self.frequency_axis: + logger.info("Adding frequency axis") + ax.set_xlim(self.frequencies[0]/Hz, self.frequencies[-1]/Hz) + ax.set_xlabel("TE011 frequency (Hz)") + #ax.set_ylim((np.array(self.ylim)**2/np.sqrt(1.64))) + ax.set_ylim((np.array(self.ylim)**2/1.64)) + ax.set_ylabel(r"Standard deviation in $m_\beta^2$ (eV$^2$)") + + self.ax2 = ax.twinx() + self.ax2.set_ylabel(r"90% CL on $m_\beta$ (eV)") + self.ax2.set_yscale("log") + self.ax2.set_ylim(self.ylim) + + elif self.exposure_axis or self.livetime_axis or self.ncavities_livetime_axis or self.ncav_eff_time_axis: + if self.exposure_axis: + logger.info("Adding exposure axis") + ax.set_xlim(self.exposure_range[0], self.exposure_range[-1]) + axis_label = r"Effective Volume $\times$ Time (m$^3$y)" #r"Efficiency $\times$ Volume $\times$ Time (m$^3$y)" + elif self.livetime_axis: + logger.info("Adding livetime axis") + ax.set_xlim(self.year_range[0], self.year_range[-1]) + axis_label = r"Livetime (years)" + #The annotations below are temporary, for pre-CDR plots + #ax.annotate("1 year", xy=[1, 0.26],xytext=[0.6, 0.16], fontsize=13, arrowprops=dict(arrowstyle="->"), bbox=dict(pad=0, facecolor="none", edgecolor="none")) + #ax.annotate("", xy=[1, 0.11],xytext=[0.8, 0.155], fontsize=13, arrowprops=dict(arrowstyle="->")) + #ax.annotate("8 years", xy=[8, 0.041],xytext=[4, 0.053], fontsize=13, arrowprops=dict(arrowstyle="->"), bbox=dict(pad=0, facecolor="none", edgecolor="none")) + elif self.ncavities_livetime_axis: + logger.info("Adding ncavities*livetime axis") + ax.set_xlim(self.ncavities_livetime_range[0], self.ncavities_livetime_range[-1]) + axis_label = r"Number of Cavities $\times$ Livetime (years)" + elif self.ncav_eff_time_axis: + logger.info("Adding ncavities*efficiency*livetime axis") + ax.set_xlim(self.ncav_eff_time_range[0], self.ncav_eff_time_range[-1]) + axis_label = r"Number of Cavities $\times$ Efficiency $\times$ Livetime (years)" + + ax.set_xlabel(axis_label) + ax.set_ylabel(r"90% CL on $m_\beta$ (eV)") + ax.set_ylim(self.ylim) + + self.ax2 = ax.twinx() + self.ax2.set_ylabel(r"Standard deviation in $m_\beta^2$ (eV$^2$)") + self.ax2.set_yscale("log") + self.ax2.set_ylim((np.array(self.ylim)**2/1.64)) + + + if self.make_key_parameter_plots: + + if self.density_axis: + + self.kp_fig, self.kp_ax = plt.subplots(1,2, figsize=(10,6)) + self.kp_fig.tight_layout() + + for ax in self.kp_ax: + ax.set_xlim(self.rhos[0]*m**3, self.rhos[-1]*m**3) + ax.set_xscale("log") + ax.set_yscale("log") + axis_label = r"Number density $n\, \, (\mathrm{m}^{-3})$" + ax.set_xlabel(axis_label) + + self.kp_ax[0].set_ylabel('Resolution (meV)') + self.kp_ax[1].set_ylabel('Track analysis duration (ms)') + + elif self.frequency_axis: + + self.kp_fig, self.kp_ax = plt.subplots(2,2, figsize=(10,10)) + + + self.kp_ax = self.kp_ax.flatten() + for ax in self.kp_ax: + ax.set_xlim(self.frequencies[0]/Hz, self.frequencies[-1]/Hz) + ax.set_xscale("log") + ax.set_yscale("log") + axis_label = "TE011 frequency (Hz)" + ax.set_xlabel(axis_label) + ax.axvline(325e6, linestyle="--", color="grey") + self.kp_ax[0].set_ylabel('Resolution (meV)') + self.kp_ax[1].set_ylabel(r'Optimum desnity (1/m$^3$)') + self.kp_ax[2].set_ylabel(r'Total and effective (dashed) Volume (m$^3$)') + self.kp_ax[3].set_ylabel('Noise temperature (K)') + + self.kp_fig.tight_layout() + + def add_track_length_axis(self): + + if self.atomic_axis: + ax2 = self.ax.twiny() + ax2.set_xscale("log") + #ax2.set_xlabel("(Atomic) track duration (s)") + ax2.set_xlabel("Track duration (s)") + + if self.sens_main_is_atomic: + ax2.set_xlim(self.sens_main.track_length(self.rhos[0])/s, + self.sens_main.track_length(self.rhos[-1])/s) + else: + if self.comparison_curve: + for sens in self.sens_ref: + if sens.Experiment.atomic: + ax2.set_xlim(sens.track_length(self.rhos[0])/s, + sens.track_length(self.rhos[-1])/s) + + if self.molecular_axis: + ax3 = self.ax.twiny() + ax3.set_xscale("log") + ax3.set_xlabel("(Molecular) track duration (s)") + + if self.atomic_axis: + ax3.spines["top"].set_position(("axes", 1.2)) + ax3.set_frame_on(True) + ax3.patch.set_visible(False) + for sp in ax3.spines.values(): + sp.set_visible(False) + ax3.spines["top"].set_visible(True) + + if not self.sens_main_is_atomic: + ax3.set_xlim(self.sens_main.track_length(self.rhos[0])/s, + self.sens_main.track_length(self.rhos[-1])/s) + else: + if self.comparison_curve: + for sens in self.sens_ref: + if not sens.Experiment.atomic: + ax3.set_xlim(sens.track_length(self.rhos[0])/s, + sens.track_length(self.rhos[-1])/s) + + if not self.molecular_axis and not self.atomic_axis: + logger.warning("No track length axis added since neither atomic nor molecular was requested") + + self.fig.tight_layout() + + def add_magnetic_field_axis(self): + ax3 = self.ax.twiny() + ax3.set_xscale("log") + ax3.set_xlabel("Magnetic field strength (T)") + + gamma = self.sens_main.T_endpoint/(me*c0**2) + 1 + ax3.set_xlim(self.frequencies[0]/(e/(2*np.pi*me)/gamma)/T, self.frequencies[-1]/(e/(2*np.pi*me)/gamma)/T) + + #def add_comparison_curve(self, label, color='k'): + + def add_arrow(self, sens): + #I am not sure why the two line below are here + if not hasattr(self, "opt_ref"): + self.add_comparison_curve() + + def get_relative(val, axis): + xmin, xmax = self.ax.get_xlim() if axis == "x" else self.ax.get_ylim() + return (np.log10(val)-np.log10(xmin))/(np.log10(xmax)-np.log10(xmin)) + + """ + rho_IV = self.rhos[self.opt_ref] + track_length_IV = self.sens_ref.track_length(rho_IV) + track_length_III = self.sens_main.track_length(rho_IV) + rho_III = rho_IV*track_length_III/track_length_IV + limit_III = sens.CL90(Experiment={"number_density": rho_III}) + + x_start = get_relative(rho_IV*m**3, "x") + y_start = get_relative(2*limit_III/eV, "y") + x_stop = get_relative(rho_III*m**3, "x") + y_stop = get_relative(limit_III/eV, "y") + plt.arrow(x_start, y_start, x_stop-x_start, y_stop-y_start, + transform = self.ax.transAxes, + facecolor = 'black', + edgecolor='k', + length_includes_head=True, + head_width=0.02, + head_length=0.03, + ) + """ + + def add_goal(self, value, label, x_pos, y_rel_position): + self.ax.axhline(value, color="darkgray", ls="solid", linewidth="1", zorder=1) + self.ax.text(x_pos, y_rel_position*value, label, fontsize=self.fontsize - 3) + + def add_density_sens_line(self, sens, plot_key_params=False, **kwargs): + limits = [] + resolutions = [] + crlb_window = [] + crlb_max_window = [] + #crlb_slope_zero_window = [] + #det_effs = [] + + temp_rho = deepcopy(sens.Experiment.number_density) + thresh_opt = [] + for rho in self.rhos: + thresh_limits = [] + #thresholds = np.linspace(10, 100, 10) #0.3*SNR_track_duration, 0.8*SNR_track_duration, 10) + for thresh in self.thresholds: + sens.Threshold.detection_threshold = thresh + thresh_limits.append(sens.CL90(Experiment={"number_density": rho})/eV) + index = np.argmin(thresh_limits) + thresh_opt.append(self.thresholds[index]) + limits.append(thresh_limits[index]) #(sens.CL90(Experiment={"number_density": rho})/eV) + resolutions.append(sens.sigma_K_noise/meV) + crlb_window.append(sens.mean_track_duration/ms) + crlb_max_window.append(sens.time_window/ms) + #crlb_slope_zero_window.append(sens.time_window_slope_zero/ms) + #det_effs.append(self.sens_main.detection_efficiency) + + #print(det_effs) + sens.Experiment.number_density = temp_rho + self.ax.plot(self.rhos*m**3, limits, **kwargs) + #logger.info('Minimum limit at {}: {}'.format(self.rhos[np.argmin(limits)]*m**3, np.min(limits))) + + if self.make_key_parameter_plots and plot_key_params: + self.kp_ax[0].plot(self.rhos*m**3, resolutions, **kwargs) + + self.kp_ax[1].plot(self.rhos*m**3, crlb_max_window, color='red', marker='.') + #self.kp_ax[1].plot(self.rhos*m**3, crlb_slope_zero_window, color='green', marker='.') + self.kp_ax[1].plot(self.rhos*m**3, crlb_window, linestyle="--", marker='.', **kwargs) + return limits + + def add_single_exposure_point(self, sens, livetime=1*year, n_cavities=1, color="red"): + logger.info("Adding exposure point") + + if livetime/year % 1 == 0: + livetime_for_label = int(round(livetime/year)) + else: + livetime_for_label = round(livetime/year, 1) + + if sens.Experiment.atomic: + #label="Atomic, reaching pilot-T target".format(n_cavities, livetime_for_label) + label="{} Phase IV cavity, {} year of livetime".format(n_cavities, livetime_for_label) + else: + label="Molecular, {} cavity, {} year".format(n_cavities, livetime_for_label) + + #exposure_fraction = n_cavities*livetime/(sens.Experiment.n_cavities*sens.Experiment.livetime) + + sens.Experiment.n_cavities = n_cavities + sens.Experiment.livetime = livetime + sens.CavityVolume() + sens.EffectiveVolume() + sens.CavityPower() + + #limit = [sens.sensitivity(Experiment={"number_density": rho})/eV**2 for rho in self.rhos] + #opt = np.argmin(limit) + #rho_opt = self.rhos[opt] + #sens.Experiment.number_density = rho_opt + + #sigma_mbeta2 = sens.sensitivity()/eV**2 + + standard_exposure = sens.EffectiveVolume()*sens.Experiment.livetime/m**3/year #*exposure_fraction + #lt = standard_exposure/sens.EffectiveVolume() + #sens.Experiment.livetime = lt + limit = sens.CL90()/eV #DETECTION THRESHOLD MAY NOT BE OPTIMIZED - NEED TO CHECK + + if self.exposure_axis: + self.ax.scatter([standard_exposure], [limit], marker="s", s=25, color=color, label=label, zorder=10) + logger.info("Exposure and mass limit for single point: {}, {}".format(standard_exposure, limit)) + elif self.livetime_axis: + self.ax.scatter([sens.Experiment.livetime/year], [limit], marker="s", s=25, color=color, label=label, zorder=10) + logger.info("Livetime and mass limit for single point: {}, {}".format(sens.Experiment.livetime/year, limit)) + elif self.ncavities_livetime_axis: + ncav_time = sens.Experiment.n_cavities*sens.Experiment.livetime/year + self.ax.scatter([ncav_time], [limit], marker="d", s=25, color='k', label=label, zorder=10) + logger.info("Ncavities*livetime and mass limit for single point: {}, {}".format(ncav_time, limit)) + elif self.ncav_eff_time_axis: + ncav_eff_time = sens.EffectiveVolume()/sens.TrapVolume()*sens.Experiment.n_cavities*sens.Experiment.livetime/year + self.ax.scatter([ncav_eff_time], [limit], marker="d", s=25, color='k', label=label, zorder=10) + logger.info("Ncavities*efficiency*livetime and mass limit for single point: {}, {}".format(ncav_eff_time, limit)) + + sens.print_statistics() + sens.print_systematics() + + def add_exposure_sens_line(self, sens, exposure_type=None, plot_key_params=False, **kwargs): + + #Optimization happens before calling this function, so I'm commenting out the below + """ + sigma_mbeta = [sens.sensitivity(Experiment={"number_density": rho})/eV**2 for rho in self.rhos] + opt = np.argmin(sigma_mbeta) + rho_opt = self.rhos[opt] + sens.Experiment.number_density = rho_opt + """ + #logger.info("Years: {}".format(sens.Experiment.livetime/year)) + + #Calcuating sensitivity for the exposure in the config file + limit_inputted_exposure = sens.CL90()/eV + limits = [] + if exposure_type=="livetime": + #self.ax.scatter([sens.Experiment.livetime/year], [np.min(sigma_mbeta)], s=40, marker="d", zorder=20, **kwargs) + self.ax.scatter([sens.Experiment.livetime/year], [limit_inputted_exposure], s=40, zorder=20, marker=kwargs["marker"], color=kwargs["color"], label=kwargs["label"]) + for lt in self.years: + sens.Experiment.livetime = lt + #limits.append(sens.sensitivity()/eV**2) + limits.append(sens.CL90()/eV) + self.ax.plot(self.years/year, limits, color=kwargs["color"], linestyle=kwargs["linestyle"]) + elif exposure_type=="ncavities_livetime": + self.ax.scatter([sens.Experiment.n_cavities*sens.Experiment.livetime/year], [limit_inputted_exposure], s=40, zorder=20, marker=kwargs["marker"], color=kwargs["color"], label=kwargs["label"]) + ncavities_livetime = [] + year_list = self.years/sens.Experiment.n_cavities + for lt in year_list: + sens.Experiment.livetime = lt + ncavities_livetime.append(sens.Experiment.n_cavities*lt/year) + limits.append(sens.CL90()/eV) + self.ax.plot(ncavities_livetime, limits, color=kwargs["color"], linestyle=kwargs["linestyle"]) + elif exposure_type=="ncav_eff_time": + ncav_eff = sens.EffectiveVolume()/sens.TrapVolume()*sens.Experiment.n_cavities + self.ax.scatter([ncav_eff*sens.Experiment.livetime/year], [limit_inputted_exposure], s=40, zorder=20, marker=kwargs["marker"], color=kwargs["color"], label=kwargs["label"]) + ncav_eff_time = [] + year_list = self.ncav_eff_times/ncav_eff + for lt in year_list: + sens.Experiment.livetime = lt + ncav_eff_time.append(ncav_eff*lt/year) + limits.append(sens.CL90()/eV) + self.ax.plot(ncav_eff_time, limits, color=kwargs["color"], linestyle=kwargs["linestyle"]) + else: + years = [] + standard_exposure = sens.EffectiveVolume()*sens.Experiment.livetime/m**3/year + #self.ax.scatter([standard_exposure], [np.min(sigma_mbeta)], s=40, marker="d", zorder=20, **kwargs) + self.ax.scatter([standard_exposure], [limits], s=40, zorder=20, marker=kwargs["marker"], color=kwargs["color"], label=kwargs["label"]) + for ex in self.exposures: + lt = ex/sens.EffectiveVolume() + years.append(lt/year) + sens.Experiment.livetime = lt + #sigma_mbetas.append(sens.sensitivity()/eV**2) + limits.append(sens.CL90()/eV) + self.ax.plot(self.exposures/m**3/year, limits, color=kwargs["color"], linestyle=kwargs["linestyle"]) #label="{} density = {:.1e} {}".format(gas, rho_opt*m**3, unit)) + + + + def add_Phase_II_exposure_sens_line(self, sens): + logger.warning("Adding Phase II sensitivity") + sens.Experiment.number_density = 2.09e17/m**3 + sens.effective_volume = 1.2*mm**3 + sens.Experiment.sri_factor = 1 #0.389*0.918*0.32 + sens.Experiment.livetime = 7185228*s + sens.CRLB_constant = 180 + + standard_exposure = sens.effective_volume*sens.Experiment.livetime/m**3/year + #sens.sensitivity() + #sens.CL90(Experiment={"number_density": 2.09e17/m**3})/eV + + sens.print_systematics() + sens.print_statistics() + sens.print_SNRs() + sens.print_Efficiencies() + + logger.info("Phase II sensitivity for exposure {} calculated: {}".format(standard_exposure, sens.sensitivity()/eV**2)) + + # Phase II experimental results from frequentist analysis + phaseIIsens = 9822 + phaseIIsense_error = 1520 + exposure_error = np.sqrt((standard_exposure*0.008)**2 + (standard_exposure*0.09)**2) + + self.ax.errorbar([standard_exposure], [phaseIIsens], yerr=[phaseIIsense_error], + marker="None", color='k', linestyle="None", elinewidth=3, + label="Phase II (measured)") + + + sigma_mbetas = [] + years = [] + for ex in self.exposures: + lt = ex/sens.effective_volume + years.append(lt/year) + sens.Experiment.livetime = lt + sigma_mbetas.append(sens.sensitivity()/eV**2) + #exposures.append(sens.EffectiveVolume()/m**3*sens.Experiment.livetime/year) + + + # unit = r"m$^{-3}$" + # gas = r"T$_2$" + self.ax.plot(self.exposures/m**3/year, sigma_mbetas, color='k', linestyle=':')#, label="{} density = {:.1e} {}".format(gas, 7.5e16, unit)) + + def get_relative(val, axis): + xmin, xmax = self.ax.get_xlim() if axis == "x" else self.ax.get_ylim() + return (np.log10(val)-np.log10(xmin))/(np.log10(xmax)-np.log10(xmin)) + + + """x_start = get_relative(1e-3, "x") + y_start = get_relative(6e1, "y") + x_stop = get_relative(1e-4, "x") + y_stop = get_relative(3e1, "y")""" + + x_start = get_relative(1e-7, "x") + y_start = get_relative(7e3, "y") + x_stop = get_relative(2e-8, "x") + y_stop = get_relative(1.3e3, "y") + self.ax.arrow(x_start, y_start, x_stop-x_start, y_stop-y_start, + transform = self.ax.transAxes, + facecolor = 'black', + edgecolor='k', + length_includes_head=True, + head_width=0.01, + head_length=0.01, + ) + + self.ax.annotate("Phase II T$_2$ density \nand resolution", xy=[x_start*0.9, y_start*1.01],textcoords="axes fraction", fontsize=13) + + def add_frequency_sens_line(self, sens, plot_key_params=True, **kwargs): + limits = [] + resolutions = [] + crlb_window = [] + crlb_max_window = [] + #crlb_slope_zero_window = [] + total_volumes = [] + effective_volumes = [] + opt_rho = [] + noise_power = [] + + configured_magnetic_field = sens.MagneticField.nominal_field + + gamma = sens.T_endpoint/(me*c0**2) + 1 + temp_rho = deepcopy(sens.Experiment.number_density) + temp_field = deepcopy(sens.MagneticField.nominal_field) + for freq in self.frequencies: + magnetic_field = freq/(e/(2*np.pi*me)/gamma) + sens.MagneticField.nominal_field = magnetic_field + + logger.info("Frequency: {:2e} Hz, Magnetic field: {:2e} T".format(freq/Hz, magnetic_field/T)) + + # Calculate new cavity properties + sens.CavityRadius() + total_volumes.append(sens.CavityVolume()/m**3) + effective_volumes.append(sens.EffectiveVolume()/m**3) + sens.CavityPower() + + # Optimization happens before calling this function, so I'm commenting out the below + # Nevermind. But the optimization should be optional, as with other variables on the x-axis. + # Note that the detection threshold is not optimized here. + rho_limits = [sens.CL90(Experiment={"number_density": rho})/eV for rho in self.rhos] + limits.append(np.min(rho_limits)) + this_optimum_rho = self.rhos[np.argmin(rho_limits)] + if this_optimum_rho == self.rhos[0] or this_optimum_rho == self.rhos[-1]: + raise ValueError("Cannot optimize density. Ideal value {:2e} at edge of range".format(this_optimum_rho*m**3)) + opt_rho.append(this_optimum_rho*m**3) + sens.Experiment.number_density = this_optimum_rho + + # other quantities + resolutions.append(sens.syst_frequency_extraction()[0]/meV) + crlb_window.append(sens.mean_track_duration/ms) + crlb_max_window.append(sens.time_window/ms) + #crlb_slope_zero_window.append(sens.time_window_slope_zero/ms) + noise_power.append(sens.noise_temp/K) + + # set rho and f back + sens.Experiment.number_density = temp_rho + sens.MagneticField.nominal_field = temp_field + sens.CavityRadius() + sens.CavityPower() + + + self.ax2.plot(self.frequencies/Hz, limits, **kwargs) + #logger.info('Minimum limit at {}: {}'.format(self.rhos[np.argmin(limits)]*m**3, np.min(limits))) + + # change cavity back to config file + """sens.MagneticField.nominal_field = configured_magnetic_field + sens.CavityRadius() + sens.CavityVolume() + sens.EffectiveVolume() + sens.CavityPower()""" + + if self.make_key_parameter_plots and plot_key_params: + self.kp_ax[0].plot(self.frequencies/Hz, resolutions, **kwargs) + + self.kp_ax[1].plot(self.frequencies/Hz, opt_rho, **kwargs) + self.kp_ax[2].plot(self.frequencies/Hz, total_volumes, **kwargs) + self.kp_ax[2].plot(self.frequencies/Hz, effective_volumes, linestyle="--", **kwargs) + self.kp_ax[3].plot(self.frequencies/Hz, noise_power, linestyle="-", **kwargs) + return limits + + def add_text(self, x, y, text, color="k"): #, fontsize=9.5 + self.ax.text(x, y, text, color=color) + + def range(self, start, stop): + cmap = matplotlib.cm.get_cmap('Spectral') + norm = matplotlib.colors.Normalize(vmin=start, vmax=stop) + return [(idx, cmap(norm(idx))) for idx in range(start, stop)] + + def save(self, savepath, **kwargs): + logger.info("Saving") + legend_fontsize = self.fontsize - 2 + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=self.legend_bbox_to_anchor, fontsize=legend_fontsize) + #handles, labels = self.fig.gca().get_legend_handles_labels() + #print(labels) + #order = [0,3,1,2] + #[handles[idx] for idx in order],[labels[idx] for idx in order], + """ + if self.density_axis: + if self.track_length_axis: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(0.15,0,1.1,0.87), fontsize=legend_fontsize) #(0.15,0,1,0.85) + else: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(0.15,0,1.,0.765), fontsize=legend_fontsize) + elif self.frequency_axis: + if self.magnetic_field_axis: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(0.14,0,1,0.85), fontsize=legend_fontsize) + else: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(0.14,0,1,0.95), fontsize=legend_fontsize) + else: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(-0.,0,0.86,0.955), fontsize=legend_fontsize) #(-0.,0,0.88,0.955) #(-0.,0,0.89,0.97) + """ + self.fig.tight_layout() + #keywords = ", ".join(["%s=%s"%(key, value) for key, value in kwargs.items()]) + metadata = {"Author": "p8/mermithid", + "Title": "Neutrino mass sensitivity", + "Subject":"90% CL upper limit on neutrino mass assuming true mass is zero." + } + #"Keywords": keywords} + if savepath is not None: + self.fig.savefig(savepath.replace(".pdf", ".png"), dpi=300, metadata=metadata) + self.fig.savefig(savepath.replace(".png", ".pdf"), bbox_inches="tight", metadata=metadata) + + if self.make_key_parameter_plots: + self.kp_fig.savefig('key_parameters.pdf') + + diff --git a/mermithid/processors/Sensitivity/ConstantSensitivityParameterPlots.py b/mermithid/processors/Sensitivity/ConstantSensitivityParameterPlots.py new file mode 100644 index 00000000..cbd26cef --- /dev/null +++ b/mermithid/processors/Sensitivity/ConstantSensitivityParameterPlots.py @@ -0,0 +1,275 @@ +''' +Calculate analytic sensitivity +function. +Author: C. Claessens +Date:12/16/2021 + +More description +''' + +from __future__ import absolute_import + +import matplotlib.pyplot as plt +import numpy as np +from scipy.optimize import minimize +import time +import warnings +warnings.filterwarnings("ignore") + +# Numericalunits is a package to handle units and some natural constants +# natural constants + +from numericalunits import meV, eV, m, T + + +# morpho imports +from morpho.utilities import morphologging, reader +from morpho.processors import BaseProcessor +from mermithid.processors.Sensitivity import AnalyticSensitivityEstimation +from mermithid.sensitivity.SensitivityFormulas import Sensitivity + + +logger = morphologging.getLogger(__name__) + + + +__all__ = [] +__all__.append(__name__) + +class ConstantSensitivityParameterPlots(AnalyticSensitivityEstimation): + ''' + Description + Args: + + Inputs: + + Output: + + ''' + # first do the BaseProcessor __init__ + def __init__(self, name, *args, **kwargs): + BaseProcessor.__init__(self, name, *args, **kwargs) + + def InternalConfigure(self, params): + ''' + Configure + ''' + # file paths + self.config_file_path = reader.read_param(params, 'config_file_path', "required") + self.sensitivity_target = reader.read_param(params, 'sensitivity_target', [0.4**2/1.64, 0.7**2/1.64, 1**2/1.64] ) + self.initialInhomogeneity = reader.read_param(params, 'initial_Inhomogeneity', 1e-8) + + self.veff_range = reader.read_param(params, 'veff_range', [0.001, 1]) + self.density_range = reader.read_param(params, 'density_range', [5.e15,1.e19]) + self.BError_range = reader.read_param(params, 'BError_range', [1e-8,1e-3]) + + # Configuration is done in __init__ of SensitivityClass + Sensitivity.__init__(self, self.config_file_path) + + self.rhos = np.linspace(self.density_range[0], self.density_range[1], 1000)/(m**3) + self.veffs = np.logspace(np.log10(self.veff_range[0]), np.log10(self.veff_range[1]), 25)*m**3 + self.Berrors = np.logspace(np.log10(self.BError_range[0]), np.log10(self.BError_range[1]), 2000) + + freq_res = np.empty(np.shape(self.rhos)) + freq_res_delta = np.empty(np.shape(self.rhos)) + + for i in range(len(self.rhos)): + self.Experiment.number_density = self.rhos[i] + freq_res[i], freq_res_delta[i] = self.syst_frequency_extraction() + + + """plt.figure() + plt.subplot(121) + plt.plot(self.rhos*m**3, freq_res/eV) + plt.xlabel('Density (/m³') + plt.ylabel('Energy resolution from frequency uncertainty (eV)') + plt.xscale('log') + + plt.subplot(122) + plt.plot(self.rhos*m**3, freq_res_delta/eV) + plt.xlabel('Density (/m³)') + plt.ylabel('Uncertainty on energy resolution from frequency uncertainty (eV)') + plt.xscale('log') + plt.tight_layout() + + plt.savefig('Frequency_resolution_vs_rho.pdf') + plt.show()""" + + return True + + + + def InternalRun(self): + + + + #self.results = {'CL90_limit': self.CL90()/eV, 'm_beta_squared_uncertainty': self.sensitivity()/(eV**2)} + #print(self.results) + + self.needed_Bs = np.empty((len(self.sensitivity_target), len(self.veffs))) + self.needed_res = np.empty((len(self.sensitivity_target), len(self.veffs))) + self.needed_res_sigma = np.empty((len(self.sensitivity_target), len(self.veffs))) + self.needed_freq_res = np.empty((len(self.sensitivity_target), len(self.veffs))) + self.needed_freq_res_sigma = np.empty((len(self.sensitivity_target), len(self.veffs))) + self.rho_opts = np.empty((len(self.sensitivity_target), len(self.veffs))) + self.CLs = np.empty((len(self.sensitivity_target), len(self.veffs))) + index = [] + + for j in range(len(self.sensitivity_target)): + for i, veff in enumerate(self.veffs): + logger.info('\nVeff = {} m³'.format(veff/(m**3))) + self.Experiment.v_eff = veff + self.MagneticField.inhomogenarity = self.initialInhomogeneity + self.rho_opts[j, i]=self.FindOptimumPressure() + self.BHomogeneityAndResNeeded(self.sensitivity_target[j]) + n = 0 + drho = self.density_range[1]/(m**3)*1.5 + + while np.abs(drho)> self.rho_opts[j, i] * 0.001 and n<10: + logger.info('Iteration: {}'.format(n)) + n+=1 + + new_rho = self.FindOptimumPressure() + drho = self.rho_opts[j, i]-new_rho + self.rho_opts[j, i] = new_rho + self.needed_Bs[j, i], self.needed_res[j, i], self.needed_res_sigma[j, i], self.needed_freq_res[j, i], self.needed_freq_res_sigma[j, i] = self.BHomogeneityAndResNeeded(self.sensitivity_target[j]) + + self.CLs[j, i] = self.CL90() + + + + index.append(np.where((self.CLs[j]/eVself.Berrors[1]))) + logger.info('Achieved 90CL limits: {}'.format(self.CLs[j]/eV)) + time.sleep(1) + + + plt.figure(figsize=(10, 5)) + + plt.subplot(121) + + for j in range(len(self.sensitivity_target)): + plt.plot(self.veffs[index[j]]/(m**3), self.needed_Bs[j][index[j]], label='90% CL = {} eV'.format(np.round(np.sqrt(1.64*self.sensitivity_target[j]),1))) + #plt.scatter(self.veffs/(m**3), self.needed_Bs, c=self.CLs/eV, marker='.') + #plt.colorbar() + plt.xlabel('Effective Volume (m³)') + plt.ylabel('Required field homogeneity') + plt.xscale('log') + plt.yscale('log') + plt.legend() + plt.tight_layout() + #plt.savefig('B_vs_veff.pdf') + + plt.subplot(122) + for j in range(len(self.sensitivity_target)): + plt.plot(self.veffs[index[j]]/(m**3), self.rho_opts[j][index[j]]*m**3, label='90% CL = {} eV'.format(np.round(np.sqrt(1.64*self.sensitivity_target[j]),1))) + plt.xlabel('Effective Volume (m³)') + plt.ylabel('Optimum number density (1/m³)') + plt.xscale('log') + #plt.yscale('log') + plt.legend() + plt.tight_layout() + + plt.savefig('B_rho_vs_veff.pdf') + plt.savefig('B_rho_vs_veff.png', dpi=300) + + plt.figure(figsize=(10, 5)) + + plt.subplot(121) + for j in range(len(self.sensitivity_target)): + plt.plot(self.veffs[index[j]]/(m**3), self.needed_res[j][index[j]]/eV, label='90% CL = {} eV'.format(np.round(np.sqrt(1.64*self.sensitivity_target[j]),1))) + plt.xlabel('Effective Volume (m³)') + plt.ylabel('Total energy resolution (eV)') + plt.xscale('log') + #plt.yscale('log') + plt.legend() + plt.tight_layout() + + plt.subplot(122) + for j in range(len(self.sensitivity_target)): + plt.plot(self.veffs[index[j]]/(m**3), self.needed_res_sigma[j][index[j]]/eV, label='90% CL = {} eV'.format(np.round(np.sqrt(1.64*self.sensitivity_target[j]),1))) + plt.xlabel('Effective Volume (m³)') + plt.ylabel('Uncertainty on total energy resolution (eV)') + plt.xscale('log') + #plt.yscale('log') + plt.legend() + plt.tight_layout() + plt.savefig('res_vs_veff.pdf') + plt.savefig('res_vs_veff.png', dpi=300) + + plt.show() + return True + + def SensVSrho(self, rho): + self.Experiment.number_density = rho + return self.CL90(Experiment={"number_density": rho}) + + def DistanceToTargetVSBError(self, BError, sensitivity_target): + sig = self.BToKeErr(self.MagneticField.nominal_field*BError, self.MagneticField.nominal_field) + self.MagneticField.usefixedvalue = True + self.MagneticField.default_systematic_smearing = sig + self.MagneticField.default_systematic_uncertainty = 0.05*sig + #print(self.sensitivity()/(eV**2), self.sensitivity_target) + return np.abs(self.sensitivity()/(eV**2)-sensitivity_target) + + def FindOptimumPressure(self): + limit = [self.SensVSrho(rho)/eV for rho in self.rhos] + opt_rho_index = np.argmin(limit) + + rho_opt = self.rhos[opt_rho_index] + if rho_opt == self.rhos[0] or rho_opt == self.rhos[-1]: + raise ValueError('Optimum rho {} is on edge or range'.format(rho_opt*m**3)) + + result = minimize(self.SensVSrho, rho_opt, method='Nelder-Mead') + if result.success: + logger.info('\tReplacing numerical value by actual optimiuation result.') + rho_opt = result.x[0] + + limit = self.CL90(Experiment={"number_density": rho_opt})/eV + logger.info('\tOptimum density: {}'.format(rho_opt*m**3)) + return rho_opt + + def FindMaxAllowedBerror(self, sensitivity_target): + distance_to_target = [self.DistanceToTargetVSBError(Berror, sensitivity_target) for Berror in self.Berrors] + optBerror_index = np.argmin(distance_to_target) + Berror_opt = self.Berrors[optBerror_index] + + + + result = minimize(self.DistanceToTargetVSBError, Berror_opt, args=(sensitivity_target), method='Nelder-Mead') + if result.success and result.x[0]>0: + Berror_opt = result.x[0] + + sig = self.BToKeErr(self.MagneticField.nominal_field*Berror_opt, self.MagneticField.nominal_field) + self.MagneticField.usefixedvalue = True + self.MagneticField.default_systematic_smearing = sig + self.MagneticField.default_systematic_uncertainty = 0.05*sig + + return Berror_opt + + + def BHomogeneityAndResNeeded(self, sensitivity_target): + neededBres = self.FindMaxAllowedBerror(sensitivity_target) + labels, sigmas, deltas = self.get_systematics() + # #b_sigma_tmp = (sigma_sys/4)**2 + # needed_total_sigma = 0 + # for i,l in enumerate(labels): + # #print('\t', l, sigmas[i]/eV, deltas[i]) + # if l != 'Magnetic Field': + # #b_sigma_tmp -= (sigmas[i]*deltas[i])**2 + # needed_total_sigma += sigmas[i]**2 + # else: + # sigma_b_old = sigmas[i] + # delta_b_old = deltas[i] + + # needed_sigma_b = sigma_b_old #np.sqrt(b_sigma_tmp)/delta_b_old + needed_sigma_b = sigmas[np.where(labels=='Magnetic Field')] + needed_b_error = self.KeToBerr(needed_sigma_b, self.MagneticField.nominal_field) + # needed_total_sigma = np.sqrt(needed_sigma_b**2 + needed_total_sigma) + needed_total_sigma = np.sqrt(np.sum(sigmas**2)) + needed_total_delta = np.sqrt(np.sum(deltas**2)) + + needed_freq_sigma = sigmas[np.where(labels=='Start Frequency Resolution')] + needed_freq_delta = deltas[np.where(labels=='Start Frequency Resolution')] + return needed_b_error/self.MagneticField.nominal_field, needed_total_sigma, needed_total_delta, needed_freq_sigma, needed_freq_delta + + diff --git a/mermithid/processors/Sensitivity/SensitivityCurveProcessor.py b/mermithid/processors/Sensitivity/SensitivityCurveProcessor.py new file mode 100644 index 00000000..ddb2f129 --- /dev/null +++ b/mermithid/processors/Sensitivity/SensitivityCurveProcessor.py @@ -0,0 +1,313 @@ +''' +Calculate sensitivity curve and plot vs. pressure +function. +Author: R. Reimann, C. Claessens +Date:11/17/2020 + +More description +''' + +from __future__ import absolute_import + + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np + + +# Numericalunits is a package to handle units and some natural constants +# natural constants +from numericalunits import e, me, c0, eps0, kB, hbar +from numericalunits import meV, eV, keV, MeV, cm, m, ns, s, Hz, kHz, MHz, GHz +from numericalunits import nT, uT, mT, T, mK, K, C, F, g, W +from numericalunits import hour, year, day +ppm = 1e-6 +ppb = 1e-9 + +# morpho imports +from morpho.utilities import morphologging, reader +from morpho.processors import BaseProcessor +from mermithid.sensitivity.SensitivityFormulas import Sensitivity + + +logger = morphologging.getLogger(__name__) + + + +__all__ = [] +__all__.append(__name__) + +class SensitivityCurveProcessor(BaseProcessor): + ''' + Description + Args: + + Inputs: + + Output: + + ''' + def InternalConfigure(self, params): + ''' + Configure + ''' + # file paths + self.config_file_path = reader.read_param(params, 'config_file_path', "required") + self.comparison_config_file_path = reader.read_param(params, 'comparison_config_file_path', '') + self.plot_path = reader.read_param(params, 'plot_path', "required") + + + # labels + self.main_curve_upper_label = reader.read_param(params, 'main_curve_upper_label', r"molecular"+"\n"+r"$V_\mathrm{eff} = 2\, \mathrm{cm}^3$"+"\n"+r"$\sigma_B = 7\,\mathrm{ppm}$") + self.main_curve_lower_label = reader.read_param(params, 'main_curve_lower_label', r"$\sigma_B = 1\,\mathrm{ppm}$") + self.comparison_curve_label = reader.read_param(params, 'comparison_curve_label', r"atomic"+"\n"+r"$V_\mathrm{eff} = 5\, \mathrm{m}^3$"+"\n"+r"$\sigma_B = 0.13\,\mathrm{ppm}$") + + + # options + self.comparison_curve = reader.read_param(params, 'comparison_curve', False) + self.B_error = reader.read_param(params, 'B_inhomogeneity', 7e-6) + self.B_error_uncertainty = reader.read_param(params, 'B_inhom_uncertainty', 0.05) + + + # plot configurations + self.figsize = reader.read_param(params, 'figsize', (6,6)) + self.density_range = reader.read_param(params, 'density_range', [1e14,1e21]) + self.ylim = reader.read_param(params, 'y_limits', [1e-2, 1e2]) + self.track_length_axis = reader.read_param(params, 'track_length_axis', True) + self.atomic_axis = reader.read_param(params, 'atomic_axis', False) + self.molecular_axis = reader.read_param(params, 'molecular_axis', False) + self.label_x_position = reader.read_param(params, 'label_x_position', 5e19) + self.upper_label_y_position = reader.read_param(params, 'upper_label_y_position', 5) + self.lower_label_y_position = reader.read_param(params, 'lower_label_y_position', 2.3) + + + # goals + self.goals = reader.read_param(params, "goals", {}) + + + # setup sensitivities + self.sens_main = Sensitivity(self.config_file_path) + self.sens_main_is_atomic = self.sens_main.Experiment.atomic + + if self.comparison_curve: + self.sens_ref = Sensitivity(self.comparison_config_file_path) + self.sens_ref_is_atomic = self.sens_ref.Experiment.atomic + + # check atomic and molecular + if self.molecular_axis: + if not self.sens_main_is_atomic: + self.molecular_sens = self.sens_main + logger.info("Main curve is molecular") + elif not self.sens_ref_is_atomic: + self.molecular_sens = self.sens_ref + logger.info("Comparison curve is molecular") + else: + raise ValueError("No experiment is configured to be molecular") + + if self.atomic_axis: + if self.sens_main_is_atomic: + self.atomic_sens = self.sens_main + logger.info("Main curve is atomic") + elif self.sens_ref_is_atomic: + self.atomic_sens = self.sens_ref + logger.info("Comparison curve is atomic") + else: + raise ValueError("No experiment is configured to be atomic") + + # densities + self.rhos = np.logspace(np.log10(self.density_range[0]), np.log10(self.density_range[1]), 100)/m**3 + + return True + + + + def InternalRun(self): + + self.create_plot() + + # add second and third x axis for track lengths + if self.track_length_axis: + self.add_track_length_axis() + + # add line for comparison using second config + if self.comparison_curve: + self.add_comparison_curve(label=self.comparison_curve_label) + #self.add_arrow(self.sens_main) + + for key, value in self.goals.items(): + logger.info('Adding goal: {}'.format(key)) + self.add_goal(value*eV, key) + + # if B is list plot line for each B + if isinstance(self.B_error, list) or isinstance(self.B_error, np.ndarray): + N = len(self.B_error) + for a, color in self.range(0, N): + sig = self.sens_main.BToKeErr(self.sens_main.MagneticField.nominal_field*self.B_error[a], self.sens_main.MagneticField.nominal_field) + self.sens_main.MagneticField.usefixedvalue = True + self.sens_main.MagneticField.default_systematic_smearing = sig + self.sens_main.MagneticField.default_systematic_uncertainty = 0.05*sig + self.add_sens_line(self.sens_main, color=color) + self.add_text(self.label_x_position, self.upper_label_y_position, self.main_curve_upper_label) + self.add_text(self.label_x_position, self.lower_label_y_position, self.main_curve_lower_label) + + else: + sig = self.sens_main.BToKeErr(self.sens_main.MagneticField.nominal_field*self.B_error[a], self.sens_main.MagneticField.nominal_field) + self.sens_main.MagneticField.usefixedvalue = True + self.sens_main.MagneticField.default_systematic_smearing = sig + self.sens_main.MagneticField.default_systematic_uncertainty = 0.05*sig + self.add_sens_line(self.sens_main, color='blue') + self.add_text(self.label_x_position, self.upper_label_y_position, self.main_curve_upper_label) + + + # save plot + self.save(self.plot_path) + + # print number of events + limit = [self.sens_main.CL90(Experiment={"number_density": rho})/eV for rho in self.rhos] + self.opt_ref = np.argmin(limit) + + rho_opt = self.rhos[self.opt_ref] + + + logger.info('Main curve (veff = {} cm**3, rho = {} /m**3):'.format(self.sens_main.Experiment.v_eff/(cm**3), rho_opt*(m**3))) + logger.info('CL90 limit: {}'.format(self.sens_main.CL90(Experiment={"number_density": rho_opt})/eV)) + logger.info('T2 in Veff: {}'.format(rho_opt*self.sens_main.Experiment.v_eff)) + logger.info('Total signal: {}'.format(rho_opt*self.sens_main.Experiment.v_eff* + self.sens_main.Experiment.LiveTime/ + self.sens_main.tau_tritium*2)) + logger.info('Signal in last eV: {}'.format(self.sens_main.last_1ev_fraction*eV**3* + rho_opt*self.sens_main.Experiment.v_eff* + self.sens_main.Experiment.LiveTime/ + self.sens_main.tau_tritium*2)) + + self.sens_main.print_statistics() + self.sens_main.print_systematics() + + return True + + + def create_plot(self): + # setup axis + self.fig, self.ax = plt.subplots(figsize=self.figsize) + ax = self.ax + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlim(self.rhos[0]*m**3, self.rhos[-1]*m**3) + ax.set_ylim(self.ylim) + #ax.grid() + + if self.atomic_axis and self.molecular_axis: + axis_label = r"(atomic / molecular) number density $\rho\, /\, \mathrm{m}^{-3}$" + elif self.atomic_axis: + axis_label = r"(atomic) number density $\rho\, /\, \mathrm{m}^{-3}$" + elif self.molecular_axis: + axis_label = r"(molecular) number density $\rho\, /\, \mathrm{m}^{-3}$" + else: + axis_label = r"number density $\rho\, /\, \mathrm{m}^{-3}$" + + ax.set_xlabel(axis_label) + ax.set_ylabel(r"90% CL $m_\beta$ / eV") + + def add_track_length_axis(self): + if self.atomic_axis: + ax2 = self.ax.twiny() + ax2.set_xscale("log") + ax2.set_xlabel("(atomic) track length / s") + ax2.set_xlim(self.atomic_sens.track_length(self.rhos[0])/s, + self.atomic_sens.track_length(self.rhos[-1])/s) + + if self.molecular_axis: + ax3 = self.ax.twiny() + + if self.atomic_axis: + ax3.spines["top"].set_position(("axes", 1.2)) + ax3.set_frame_on(True) + ax3.patch.set_visible(False) + for sp in ax3.spines.values(): + sp.set_visible(False) + ax3.spines["top"].set_visible(True) + + ax3.set_xscale("log") + ax3.set_xlabel("(molecular) track length / s") + ax3.set_xlim(self.molecular_sens.track_length(self.rhos[0])/s, + self.molecular_sens.track_length(self.rhos[-1])/s) + + else: + logger.warning("No track length axis added since neither atomic nor molecular was requested") + self.fig.tight_layout() + + def add_comparison_curve(self, label, color='k'): + limit = [self.sens_ref.CL90(Experiment={"number_density": rho})/eV for rho in self.rhos] + self.opt_ref = np.argmin(limit) + + rho_opt = self.rhos[self.opt_ref] + logger.info('Ref. curve (veff = {} m**3):'.format(self.sens_ref.Experiment.v_eff/(m**3))) + logger.info('T in Veff: {}'.format(rho_opt*self.sens_ref.Experiment.v_eff)) + logger.info('Ref. total signal: {}'.format(rho_opt*self.sens_ref.Experiment.v_eff* + self.sens_ref.Experiment.LiveTime/ + self.sens_ref.tau_tritium)) + + self.ax.plot(self.rhos*m**3, limit, color=color) + self.ax.axvline(self.rhos[self.opt_ref]*m**3, ls=":", color="gray", alpha=0.4) + + #self.ax.axhline(0.04, color="gray", ls="--") + #self.ax.text(3e14, 0.044, "Phase IV (40 meV)") + self.ax.text(1.5e14, 0.11, label) + + def add_arrow(self, sens): + if not hasattr(self, "opt_ref"): + self.add_comparison_curve() + + def get_relative(val, axis): + xmin, xmax = self.ax.get_xlim() if axis == "x" else self.ax.get_ylim() + return (np.log10(val)-np.log10(xmin))/(np.log10(xmax)-np.log10(xmin)) + + rho_IV = self.rhos[self.opt_ref] + track_length_IV = self.sens_ref.track_length(rho_IV) + track_length_III = self.sens_main.track_length(rho_IV) + rho_III = rho_IV*track_length_III/track_length_IV + limit_III = sens.CL90(Experiment={"number_density": rho_III}) + + x_start = get_relative(rho_IV*m**3, "x") + y_start = get_relative(2*limit_III/eV, "y") + x_stop = get_relative(rho_III*m**3, "x") + y_stop = get_relative(limit_III/eV, "y") + plt.arrow(x_start, y_start, x_stop-x_start, y_stop-y_start, + transform = self.ax.transAxes, + facecolor = 'black', + edgecolor='k', + length_includes_head=True, + head_width=0.02, + head_length=0.03, + ) + + def add_goal(self, value, label): + self.ax.axhline(value/eV, color="gray", ls="--") + self.ax.text(3e18, 0.8*value/eV, label) + + def add_sens_line(self, sens, **kwargs): + limits = [sens.CL90(Experiment={"number_density": rho})/eV for rho in self.rhos] + self.ax.plot(self.rhos*m**3, limits, **kwargs) + logger.info('Minimum limit at {}: {}'.format(self.rhos[np.argmin(limits)]*m**3, np.min(limits))) + + def add_text(self, x, y, text): + self.ax.text(x, y, text) + + def range(self, start, stop): + cmap = matplotlib.cm.get_cmap('Spectral') + norm = matplotlib.colors.Normalize(vmin=start, vmax=stop-1) + return [(idx, cmap(norm(idx))) for idx in range(start, stop)] + + def save(self, savepath, **kwargs): + self.fig.tight_layout() + #keywords = ", ".join(["%s=%s"%(key, value) for key, value in kwargs.items()]) + metadata = {"Author": "p8/mermithid", + "Title": "Neutrino mass sensitivity", + "Subject":"90% CL upper limit on neutrino mass assuming true mass is zero." + } + #"Keywords": keywords} + if savepath is not None: + self.fig.savefig(savepath.replace(".pdf", ".png"), dpi=300, metadata=metadata) + self.fig.savefig(savepath.replace(".png", ".pdf"), bbox_inches="tight", metadata=metadata) + + diff --git a/mermithid/processors/Sensitivity/SensitivityParameterScanProcessor.py b/mermithid/processors/Sensitivity/SensitivityParameterScanProcessor.py new file mode 100644 index 00000000..9be815bd --- /dev/null +++ b/mermithid/processors/Sensitivity/SensitivityParameterScanProcessor.py @@ -0,0 +1,473 @@ +''' +Scan a parameter and calculate the sensitivity curve for each value of the parameter. +Author: C. Claessens +Date: 03/14/2024 + +More description +''' + +from __future__ import absolute_import + + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import os +import concurrent.futures + + + +# Numericalunits is a package to handle units and some natural constants +# natural constants +from numericalunits import e, me, c0, eps0, kB, hbar +from numericalunits import meV, eV, keV, MeV, cm, m, mm +from numericalunits import nT, uT, mT, T, mK, K, C, F, g, W +from numericalunits import hour, year, day, ms, ns, s, Hz, kHz, MHz, GHz +ppm = 1e-6 +ppb = 1e-9 +deg = np.pi/180 + +# morpho imports +from morpho.utilities import morphologging, reader +from morpho.processors import BaseProcessor +from mermithid.sensitivity.SensitivityFormulas import Sensitivity +from mermithid.sensitivity.SensitivityCavityFormulas import CavitySensitivity + + +logger = morphologging.getLogger(__name__) + + + +__all__ = [] +__all__.append(__name__) + + +def return_var_name(variable): + for name in globals(): + if eval(name) == variable: + return name + + + +class SensitivityParameterScanProcessor(BaseProcessor): + ''' + Description + Args: + + Inputs: + + Output: + + ''' + def InternalConfigure(self, params): + ''' + Configure + ''' + # file paths + self.config_file_path = reader.read_param(params, 'config_file_path', "required") + self.plot_path = reader.read_param(params, 'plot_path', "required") + + + # options + self.scan_parameter_name = reader.read_param(params, 'scan_parameter_name', 'MagneticField.sigmae_r') + self.scan_parameter_range = reader.read_param(params, "scan_parameter_range", [0.1, 2.1]) + self.scan_parameter_steps = reader.read_param(params, "scan_parameter_steps", 3) + self.scan_parameter_scale = reader.read_param(params, "scan_parameter_scale", "lin") + scan_parameter_unit = reader.read_param(params, "scan_parameter_unit", eV) + + + self.scan_parameter_unit_string = return_var_name(scan_parameter_unit) + print("Unit string: ", self.scan_parameter_unit_string) + self.scan_parameter_unit = scan_parameter_unit + + # main plot configurations + self.figsize = reader.read_param(params, 'figsize', (6,6)) + self.legend_location = reader.read_param(params, 'legend_location', 'upper left') + self.fontsize = reader.read_param(params, 'fontsize', 12) + self.plot_sensitivity_scan_on_log_scale = reader.read_param(params, 'plot_sensitivity_scan_on_log_scale', True) + + self.density_axis = reader.read_param(params, "density_axis", True) + self.track_length_axis = reader.read_param(params, 'track_length_axis', True) + self.atomic_axis = reader.read_param(params, 'atomic_axis', False) + self.molecular_axis = reader.read_param(params, 'molecular_axis', False) + + self.density_range = reader.read_param(params, 'density_range', [1e14,1e21]) + self.ylim = reader.read_param(params, 'y_limits', [1e-2, 1e2]) + + + + # key parameter plots + self.make_key_parameter_plots = reader.read_param(params, 'plot_key_parameters', False) + + if self.density_axis: + self.add_sens_line = self.add_density_sens_line + logger.info("Doing density lines") + + + + # goals + self.goals = reader.read_param(params, "goals", {}) + self.goal_x_pos = reader.read_param(params, "goals_x_position", 1e14) + self.goals_y_rel_position = reader.read_param(params, "goals_y_rel_position", 0.75) + + + # setup sensitivities + + self.cavity = reader.read_param(params, 'cavity', True) + + if self.cavity: + self.sens_main = CavitySensitivity(self.config_file_path) + else: + self.sens_main = Sensitivity(self.config_file_path) + self.sens_main_is_atomic = self.sens_main.Experiment.atomic + + + # Setting natoms_per_particle for atomic and molecular cases, to use later + # when calculating and printing event rates + if self.sens_main_is_atomic: + self.sens_main_natoms_per_particle = 1 + else: + self.sens_main_natoms_per_particle = 2 + + + # check atomic and molecular + if self.molecular_axis: + if not self.sens_main_is_atomic: + #self.molecular_sens = self.sens_main + logger.info("Main curve is molecular") + + + if self.atomic_axis: + if self.sens_main_is_atomic: + #self.atomic_sens = self.sens_main + logger.info("Main curve is atomic") + else: + logger.warn("No experiment is configured to be atomic") + + # densities, exposures, runtimes + self.rhos = np.logspace(np.log10(self.density_range[0]), np.log10(self.density_range[1]), 100)/m**3 + if(self.scan_parameter_scale == "lin"): + self.scan_parameter_values = np.linspace(self.scan_parameter_range[0], self.scan_parameter_range[1], self.scan_parameter_steps)*self.scan_parameter_unit + elif(self.scan_parameter_scale == "log"): + self.scan_parameter_values = np.logspace(np.log10(self.scan_parameter_range[0]), np.log10(self.scan_parameter_range[1]), self.scan_parameter_steps)*self.scan_parameter_unit + else: + logger.warn("Unexpected parameter scale, assuming linear") + self.scan_parameter_values = np.linspace(self.scan_parameter_range[0], self.scan_parameter_range[1], self.scan_parameter_steps)*self.scan_parameter_unit + return True + + + + def InternalRun(self): + + self.create_plot(self.scan_parameter_values/self.scan_parameter_unit) + # add second and third x axis for track lengths + if self.density_axis and self.track_length_axis: + self.add_track_length_axis() + + # add goals to density plot + for key, value in self.goals.items(): + logger.info('Adding goal: {} = {}'.format(key, value)) + self.add_goal(value, key) + + # arrays for results values for output + self.optimum_limits = [] + self.optimum_rhos = [] + + self.noise_temp = [] + self.SNR = [] + self.track_duration = [] + self.total_sigma = [] + self.sys_lim = [] + + for i, color in self.range(self.scan_parameter_values/self.scan_parameter_unit): + parameter_value = self.scan_parameter_values[i] + + category, param = self.scan_parameter_name.split(".") + + # read current value of param + try: + current_value = self.sens_main.__dict__[category].__dict__[param] + logger.info(f"Current value of {param}: {current_value/self.scan_parameter_unit}") + except KeyError as e: + logger.error(f"Parameter {param} not found in {category}") + raise e + + + # Set to scan param value + self.sens_main.__dict__[category].__dict__[param] = parameter_value + # Re-calc the cavity init with the new param + self.sens_main.CalcDefaults(overwrite=True) + # Ensure param scan value unchanged + read_back = self.sens_main.__dict__[category].__dict__[param] + #setattr(self.sens_main, self.scan_parameter_name, parameter_value) + #read_back = getattr(self.sens_main, self.scan_parameter_name) + logger.info(f"Setting {self.scan_parameter_name} to {parameter_value/self.scan_parameter_unit} and reading back: {read_back/ self.scan_parameter_unit}") + + # pitch angle set equal + if(param == "min_pitch_used_in_analysis"): + self.sens_main.__dict__["FrequencyExtraction"].__dict__["minimum_angle_in_bandwidth"] = parameter_value + + # If the scanned param isn't trap length, calc trap length for cavity L/D + if (param != "trap_length"): + self.sens_main.TrapLength() + + logger.info("Calculating cavity experiment radius, volume, effective volume, power") + self.sens_main.CavityRadius() + self.sens_main.CavityVolume() + self.sens_main.EffectiveVolume() + self.sens_main.CavityPower() + + # optimize density + logger.info("Optimizing density") + limit = [self.sens_main.CL90(Experiment={"number_density": rho}) for rho in self.rhos] + self.optimum_limits.append(np.min(limit)) + opt = np.argmin(limit) + rho_opt = self.rhos[opt] + self.optimum_rhos.append(rho_opt) + + # add main curve + logger.info("Drawing main curve") + label = f"{param} = {parameter_value/self.scan_parameter_unit:.2f} {self.scan_parameter_unit_string}" + self.add_sens_line(self.sens_main, label=label, color=color) + + if self.make_key_parameter_plots: + logger.info("Making key parameter plots") + # First key parameter plot: Stat and Syst vs. density + + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = [], [], [] + + for n in self.rhos: + self.sens_main.CL90(Experiment={"number_density": n}) + labels, sigmas, deltas = self.sens_main.get_systematics() + sigma_startf.append(sigmas[1]) + stat_on_mbeta2.append(self.sens_main.StatSens()) + syst_on_mbeta2.append(self.sens_main.SystSens()) + + sigma_startf, stat_on_mbeta2, syst_on_mbeta2 = np.array(sigma_startf), np.array(stat_on_mbeta2), np.array(syst_on_mbeta2) + fig = plt.figure() + plt.loglog(self.rhos*m**3, stat_on_mbeta2/eV**2, label='Statistical uncertainty') + plt.loglog(self.rhos*m**3, syst_on_mbeta2/eV**2, label='Systematic uncertainty') + plt.xlabel(r"Number density $n\, \, (\mathrm{m}^{-3})$") + plt.ylabel(r"Standard deviation in $m_\beta^2$ (eV$^2$)") + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(self.plot_path, f"{param}_{parameter_value/self.scan_parameter_unit}_stat_and_syst_vs_density.pdf")) + + + + logger.info('Experiment info:') + # set optimum density back + self.sens_main.CL90(Experiment={"number_density": rho_opt}) + logger.info('veff = {} m**3, rho = {} /m**3:'.format(self.sens_main.effective_volume/(m**3), rho_opt*(m**3))) + logger.info("Loaded Q: {}".format(self.sens_main.loaded_q)) + logger.info("Axial frequency for minimum detectable angle: {} MHz".format(self.sens_main.required_bw_axialfrequency/MHz)) + logger.info("Total bandwidth: {} MHz".format(self.sens_main.required_bw/MHz)) + logger.info('Larmor power = {} W, Hanneke power = {} W'.format(self.sens_main.larmor_power/W, self.sens_main.signal_power/W)) + logger.info('Hanneke / Larmor power = {}'.format(self.sens_main.signal_power/self.sens_main.larmor_power)) + + if self.sens_main.FrequencyExtraction.crlb_on_sidebands: + logger.info("Uncertainty from determination of f_carrier and f_lsb, due to noise: {} eV".format(self.sens_main.sigma_K_noise/eV)) + + noise_temp, SNR, track_duration = self.sens_main.print_SNRs(rho_opt) + # Store relevant values + self.noise_temp.append(noise_temp/K) + self.SNR.append(SNR) + self.track_duration.append(track_duration/ms) + + logger.info('CL90 limit: {}'.format(self.sens_main.CL90(Experiment={"number_density": rho_opt})/eV)) + logger.info('T2 in Veff: {}'.format(rho_opt*self.sens_main.effective_volume)) + logger.info('Total signal: {}'.format(rho_opt*self.sens_main.effective_volume* + self.sens_main.Experiment.LiveTime/ + self.sens_main.tau_tritium*self.sens_main_natoms_per_particle)) + logger.info('Signal in last eV: {}'.format(self.sens_main.last_1ev_fraction*eV**3* + rho_opt*self.sens_main.effective_volume* + self.sens_main.Experiment.LiveTime/ + self.sens_main.tau_tritium*self.sens_main_natoms_per_particle)) + + self.sens_main.print_Efficiencies() + self.sens_main.print_statistics() + systematic_limit, total_sigma = self.sens_main.print_systematics() + self.sys_lim.append(systematic_limit) + self.total_sigma.append(total_sigma) + + self.save("sensitivity_vs_density_for_{}_scan.pdf".format(param)) + + + # plot and print best limits + self.results = {"scan_parameter": self.scan_parameter_name, "scan parameter_unit": self.scan_parameter_unit_string, + "scan_parameter_values": self.scan_parameter_values, "optimum_limits_eV": np.array(self.optimum_limits)/eV, + "optimum_densities/m3": np.array(self.optimum_rhos)*(m**3), + "Noise Temperatures/K": np.array(self.noise_temp), + "SNRs 1eV from temperature": np.array(self.SNR), "track durations": np.array(self.track_duration), + "Systematic limits": np.array(self.sys_lim), "Total Sigmas": np.array(self.total_sigma)} + + results_array = [np.array(self.scan_parameter_values/self.scan_parameter_unit),np.array(self.noise_temp),np.array(self.SNR), + np.array(self.optimum_rhos)*(m**3),np.array(self.track_duration),np.array(self.total_sigma), + 1000*np.array(self.optimum_limits)/eV,np.array(self.sys_lim)] + fmt_array = '%.2g','%.3f','%.1f','%.2E','%.2f','%.1f','%.1f','%.1f' + header_string = 'Param Value, Noise temperature /K, SNR, Gas Density /m3, Track Length /ms, Resolution, Sensitivity /meV, Systematic Limit /meV' + np.savetxt("results_array_{}.csv".format(param),np.array(results_array).T,delimiter=',',fmt=fmt_array,header=header_string) + logger.info("Scan parameter: {}".format(self.scan_parameter_name)) + logger.info("Tested parameter values: {}".format(self.scan_parameter_values/self.scan_parameter_unit)) + logger.info("Best limits: {}".format(np.array(self.optimum_limits)/eV)) + + plt.figure(figsize=self.figsize) + #plt.title("Sensitivity vs. {}".format(self.scan_parameter_name)) + plt.plot(self.scan_parameter_values/self.scan_parameter_unit, np.array(self.optimum_limits)/eV, marker=".", label="Density optimized scenarios") + plt.xlabel(f"{param} ({self.scan_parameter_unit_string})", fontsize=self.fontsize) + plt.ylabel(r"90% CL on $m_\beta$ (eV)", fontsize=self.fontsize) + if self.plot_sensitivity_scan_on_log_scale: + plt.yscale("log") + # TODO log x here + if(self.scan_parameter_scale == "log"): + plt.xscale("log") + for key, value in self.goals.items(): + logger.info('Adding goal: {} = {}'.format(key, value)) + plt.axhline(value, label=key, color="grey", linestyle="--") + plt.legend(fontsize=self.fontsize) + plt.tight_layout() + plt.savefig(os.path.join(self.plot_path, f"{param}_scan_optimum_limits.pdf")) + plt.show() + + + + + + + return True + + + def create_plot(self, param_range=[]): + # setup axis + plt.rcParams.update({'font.size': self.fontsize}) + self.fig, self.ax = plt.subplots(figsize=self.figsize) + ax = self.ax + ax.set_xscale("log") + ax.set_yscale("log") + if self.density_axis: + logger.info("Adding density axis") + ax.set_xlim(self.rhos[0]*m**3, self.rhos[-1]*m**3) + + if self.atomic_axis and self.molecular_axis: + axis_label = r"(Atomic / molecular) number density $n\, \, (\mathrm{m}^{-3})$" + elif self.atomic_axis: + axis_label = r"(Atomic) number density $n\, \, (\mathrm{m}^{-3})$" + elif self.molecular_axis: + axis_label = r"(Molecular) number density $n\, \, (\mathrm{m}^{-3})$" + else: + axis_label = r"Number density $n\, \, (\mathrm{m}^{-3})$" + + ax.set_xlabel(axis_label) + ax.set_ylim(self.ylim) + ax.set_ylabel(r"90% CL on $m_\beta$ (eV)") + + if len(param_range)>4: + # add colorbar with colors from self.range + cmap = matplotlib.cm.get_cmap('Spectral') + norm = matplotlib.colors.Normalize(vmin=np.min(param_range), vmax=np.max(param_range)) + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + self.fig.colorbar(sm, ticks=np.round(param_range, 2), label=f"{self.scan_parameter_name} ({self.scan_parameter_unit_string})") + + + + + def add_track_length_axis(self): + + if self.atomic_axis: + + ax2 = self.ax.twiny() + ax2.set_xscale("log") + ax2.set_xlabel("(Atomic) track length (s)") + ax2.set_xlim(self.sens_main.track_length(self.rhos[0])/s, + self.sens_main.track_length(self.rhos[-1])/s) + + if self.molecular_axis: + ax3 = self.ax.twiny() + + if self.atomic_axis: + ax3.spines["top"].set_position(("axes", 1.2)) + ax3.set_frame_on(True) + ax3.patch.set_visible(False) + for sp in ax3.spines.values(): + sp.set_visible(False) + ax3.spines["top"].set_visible(True) + + ax3.set_xscale("log") + ax3.set_xlabel("(Molecular) track length (s)") + ax3.set_xlim(self.sens_main.track_length(self.rhos[0])/s, + self.sens_main.track_length(self.rhos[-1])/s) + + if not self.atomic_axis and not self.molecular_axis: + logger.warning("No track length axis added since neither atomic nor molecular was requested") + self.fig.tight_layout() + + + + def add_goal(self, value, label): + self.ax.axhline(value, color="gray", ls="--") + self.ax.text(self.goal_x_pos, self.goals_y_rel_position*value, label) + + def add_density_sens_line(self, sens, plot_key_params=False, **kwargs): + limits = [] + resolutions = [] + crlb_window = [] + crlb_max_window = [] + crlb_slope_zero_window = [] + + for rho in self.rhos: + limits.append(sens.CL90(Experiment={"number_density": rho})/eV) + resolutions.append(sens.sigma_K_noise/meV) + crlb_window.append(sens.best_time_window/ms) + crlb_max_window.append(sens.time_window/ms) + crlb_slope_zero_window.append(sens.time_window_slope_zero/ms) + + + self.ax.plot(self.rhos*m**3, limits, **kwargs) + rho_opt = self.rhos[np.argmin(limits)] + # set experiment to optimum density + sens.CL90(Experiment={"number_density": rho_opt}) + logger.info('Minimum limit at {}: {}'.format(rho_opt*m**3, np.min(limits))) + + if self.make_key_parameter_plots and plot_key_params: + self.kp_ax[0].plot(self.rhos*m**3, resolutions, **kwargs) + + self.kp_ax[1].plot(self.rhos*m**3, crlb_max_window, color='red', marker='.') + self.kp_ax[1].plot(self.rhos*m**3, crlb_slope_zero_window, color='green', marker='.') + self.kp_ax[1].plot(self.rhos*m**3, crlb_window, linestyle="--", marker='.', **kwargs) + return limits + + + + def add_text(self, x, y, text, color="k"): #, fontsize=9.5 + self.ax.text(x, y, text, color=color) + + def range(self, param_range): + cmap = matplotlib.cm.get_cmap('Spectral') + norm = matplotlib.colors.Normalize(vmin=0, vmax=len(param_range)-1) + return [(idx, cmap(norm(idx))) for idx, _ in enumerate(param_range)] + + def save(self, filename, **kwargs): + + if self.density_axis: + if self.scan_parameter_steps < 5: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(0.15,0,1,0.765)) + + else: + legend=self.fig.legend(loc=self.legend_location, framealpha=0.95, bbox_to_anchor=(-0.,0,0.89,0.97)) + + + + #keywords = ", ".join(["%s=%s"%(key, value) for key, value in kwargs.items()]) + metadata = {"Author": "p8/mermithid", + "Title": "Neutrino mass sensitivity", + "Subject":"90% CL upper limit on neutrino mass assuming true mass is zero." + } + #"Keywords": keywords} + + self.fig.tight_layout() + self.fig.savefig(os.path.join(self.plot_path, filename), bbox_inches="tight", metadata=metadata) + self.fig.savefig(os.path.join(self.plot_path, filename.replace(".pdf", ".png")), bbox_inches="tight", metadata=metadata) + + + diff --git a/mermithid/processors/Sensitivity/__init__.py b/mermithid/processors/Sensitivity/__init__.py new file mode 100644 index 00000000..918dfef2 --- /dev/null +++ b/mermithid/processors/Sensitivity/__init__.py @@ -0,0 +1,10 @@ +''' +''' + +from __future__ import absolute_import + +from .SensitivityCurveProcessor import SensitivityCurveProcessor +from .CavitySensitivityCurveProcessor import CavitySensitivityCurveProcessor +from .AnalyticSensitivityEstimation import AnalyticSensitivityEstimation +from .ConstantSensitivityParameterPlots import ConstantSensitivityParameterPlots +from .SensitivityParameterScanProcessor import SensitivityParameterScanProcessor diff --git a/mermithid/processors/__init__.py b/mermithid/processors/__init__.py index bc325afd..d7083508 100644 --- a/mermithid/processors/__init__.py +++ b/mermithid/processors/__init__.py @@ -8,3 +8,4 @@ from . import plots from . import misc from . import Fitters +from . import Sensitivity diff --git a/mermithid/processors/misc/FrequencyShifter.py b/mermithid/processors/misc/FrequencyShifter.py index dfc900aa..d65d46ef 100644 --- a/mermithid/processors/misc/FrequencyShifter.py +++ b/mermithid/processors/misc/FrequencyShifter.py @@ -32,7 +32,7 @@ def InternalRun(self): self.results = {} return True if self.frequency_shift == 0: - looger.warning("Zero frequency shift!") + logger.warning("Zero frequency shift!") self.results = self.frequencies return True for frequency in self.frequencies: diff --git a/mermithid/sensitivity/SensitivityCavityFormulas.py b/mermithid/sensitivity/SensitivityCavityFormulas.py new file mode 100755 index 00000000..f64e4a9f --- /dev/null +++ b/mermithid/sensitivity/SensitivityCavityFormulas.py @@ -0,0 +1,1001 @@ +''' +Class calculating neutrino mass sensitivities based on analytic formulas from CDR. +Author: R. Reimann, C. Claessens, T. E. Weiss, W. Van De Pontseele +Date: 06/07/2023 +Updated: December 2024 + +The statistical method and formulas are described in +CDR (CRES design report, Section 1.3) https://www.overleaf.com/project/5b9314afc673d862fa923d53. +''' +import numpy as np +from scipy.stats import ncx2, chi2 +from scipy.special import roots_laguerre +import matplotlib.pyplot as plt + +from mermithid.misc.Constants_numericalunits import * +from mermithid.misc.CRESFunctions_numericalunits import * +from mermithid.cavity.HannekeFunctions import * +from mermithid.sensitivity.SensitivityFormulas import * + + + +try: + from morpho.utilities import morphologging + logger = morphologging.getLogger(__name__) +except: + print("Run without morpho!") + + + +# Wouters functinos +def db_to_pwr_ratio(q_db): + return 10**(q_db/10) + +def axial_motion(magnetic_field, pitch, trap_length, minimum_trapped_pitch, kin_energy, flat_fraction=0.5, trajectory = None): + # returns the axial motion frequency and a trajectory of point along the axial motion + # also return the average magnetic field seen by the electron + # from z=0 to z=cavity_length/2 with npoints set by the trajectory variable + # See LUCKEY write-up for a little more on Talia's "flat fraction" trap model + + # Input parameters: + # pitch and minimum_trapped_pitch are in radians + + # Axial motion: + z_w = trap_length/2 + speed = beta(kin_energy)*c0 + transverse_speed = speed*np.cos(pitch) + tan_min = np.tan(minimum_trapped_pitch) + # Axial frequency + time_flat = z_w*flat_fraction/transverse_speed + time_harmonic = np.pi*z_w*(1-flat_fraction)*tan_min/(2*speed*np.sin(pitch)) + axial_frequency = 1/4/(time_flat+time_harmonic) + + #Average magnetic field: + magnetic_field_avg_harm = magnetic_field/2*(1+1/np.sin(pitch)**2) + magnetic_field_avg = (magnetic_field_avg_harm*time_harmonic + magnetic_field*time_flat)/(time_harmonic+time_flat) + + # Trajectory: + if trajectory is None: + z_t = None + else: + omega_harm = speed*np.sin(pitch)/z_w/tan_min/(1-flat_fraction) + time = np.linspace(0, time_flat+time_harmonic, trajectory) + z_t = np.heaviside(time_flat-time, 0.5)*time*transverse_speed +\ + np.heaviside(time-time_flat, 0.5)*(z_w*flat_fraction + z_w*(1-flat_fraction)*tan_min/np.tan(pitch)*np.sin(omega_harm*(time-time_flat))) + + return axial_frequency, magnetic_field_avg, z_t + +def magnetic_field_flat_harmonic(z, magnetic_field, trap_length, minimum_trapped_pitch, flat_fraction=0.5): + z_w = trap_length/2 + a = z_w*(1-flat_fraction)*np.tan(minimum_trapped_pitch) + return magnetic_field*(1+np.heaviside(np.abs(z)-z_w*flat_fraction, 0.5)*(np.abs(z)-z_w*flat_fraction)**2/a**2) + + +def axial_frequency_box(length, kin_energy, max_pitch_angle=86*np.pi/180): + #pitch_max = max_pitch_angle/180*np.pi + return (beta(kin_energy)*c0*np.cos(max_pitch_angle)) / (2*length) + +def mean_field_frequency_variation(cyclotron_frequency, length_diameter_ratio, max_pitch_angle=86*np.pi/180, q=0.16): + # Because of the different electron trajectories in the trap, + # An electron will see a slightly different magnetic field + # depending on its position in the trap, especially the pitch angle. + # This is a rough estimate of the mean field variation, inspired by calculation performed by Rene. + #y = (90-max_pitch_angle)/4 + phi_rad = (np.pi/2-max_pitch_angle) + return q*phi_rad**2*cyclotron_frequency*(10/length_diameter_ratio) + #return 0.002*y**2*cyclotron_frequency*(10/length_diameter_ratio) + +# Noise power entering the amplifier, inclding the transmitted noise from the cavity and the reflected noise from the circulator. +# Insertion loss is included. +def Pn_dut_entrance(t_cavity, + t_amplifier, + att_line_db, + att_cir_db, + coupling, + freq, + bandwidth, + loaded_Q): + att_cir = db_to_pwr_ratio(att_cir_db) + att_line = db_to_pwr_ratio(att_line_db) + assert( (np.all(att_cir<=1)) & (np.all(att_line<=1)) ) + + # Calculate the noise at the cavity + Pn_cav = Pn_cavity(t_cavity, coupling, loaded_Q, bandwidth, freq) + Pn_circulator = t_effective(t_amplifier, freq)*kB*bandwidth + Pn_circulator_after_reflection = Pn_reflected(Pn_f(Pn_circulator,t_amplifier, t_cavity,att_line, bandwidth), coupling, loaded_Q, bandwidth, freq) + # Propagate the noise over the line towards the circulator + Pn_entrance = Pn_f(Pn_circulator_after_reflection+Pn_cav,t_cavity,t_amplifier,att_line,bandwidth) + # Apply the effect of the circulator + return Pn_f(Pn_entrance,t_amplifier,t_amplifier,att_cir,bandwidth) + +# Noise power genrated in the cavity integrated over the bandwidth that couples into the readout line. +def Pn_cavity(t_cavity, coupling, loaded_Q, bandwidth, freq): + return kB*t_effective(t_cavity, freq)*4*coupling/(1+coupling)**2*freq/loaded_Q*np.arctan(loaded_Q*bandwidth/freq) + +# Noise power reflecting of the cavity +def Pn_reflected(Pn_incident, coupling, loaded_Q, bandwidth, freq): + + reflection_coefficient = 1-freq/loaded_Q/bandwidth*np.arctan(loaded_Q*bandwidth/freq)*4*coupling/(1+coupling)**2 + return Pn_incident*reflection_coefficient + +# Power at the end of a lossy line with temperature gradient +def Pn_f(Pn_i,t_i,t_f,a,bandwidth): # eq 10 + if hasattr(a, "__len__") or a!=1: + return Pn_i+ kB*bandwidth*(t_f-t_i)*(1+ (1-a)/np.log(a))+ (t_i*kB*bandwidth-Pn_i)*(1-a) + else: + return Pn_i*np.ones_like(t_f) + +# Effective temperature taking the quantum photon population into account. +def t_effective(t_physical, cyclotron_frequency): + quantum = 2*np.pi*hbar*cyclotron_frequency/kB + #for numerical stability + if np.all(quantum/t_physical < 1e-2): + return t_physical + else: + return quantum*(1/2+1/(np.exp(quantum/t_physical)-1)) + + +# Calculate threshold z to trap electrons born at some pitch angle theta_start +# Electrons are trapped if they start at z values less than this threshold +# (Considering one axial side of the trap) +def max_z_to_trap_vs_theta_start(theta_start, trap_length, minimum_trapped_pitch, flat_fraction=0.5): + z_w = trap_length/2 + sec_min = 1/np.cos(minimum_trapped_pitch) + sin2_min = np.sin(minimum_trapped_pitch)**2 + sin2_start = np.sin(theta_start)**2 + return z_w*flat_fraction + z_w*(1-flat_fraction)*sec_min*np.sqrt(sin2_start - sin2_min) + +def dist_of_theta_start_after_trapping(theta_start, trap_length, minimum_trapped_pitch, flat_fraction=0.5): + # Distribution of theta_start for electrons born uniformly along z + # Multiplied by sin(theta_start), to account for the birth pitch angles - is this + # correct? Is another normalization needed after multiplying by sin(theta_start)? + z_threshold = max_z_to_trap_vs_theta_start(theta_start, trap_length, minimum_trapped_pitch, flat_fraction) + return z_threshold/(trap_length/2)*np.sin(theta_start) + +def theta_bottom_from_theta_start(theta_start, B_min, B_start): + return np.arcsin(np.sin(theta_start)*np.sqrt(B_min/B_start)) + +def dist_of_theta_bottom_after_trapping(B_min, theta_start_array, trap_length, minimum_trapped_pitch, flat_fraction=0.5, n_z_start=1000, n_theta_bottom=10): + z_start_array = np.linspace(0, trap_length/2, n_z_start) + B_start_array = magnetic_field_flat_harmonic(z_start_array, B_min, trap_length, minimum_trapped_pitch, flat_fraction) + theta_bottoms = [] + for theta_start in theta_start_array: + theta_bottoms.append(theta_bottom_from_theta_start(theta_start, B_min, B_start_array)) + theta_bottoms = np.array(theta_bottoms) + theta_bottoms_bin_centers = np.linspace(minimum_trapped_pitch, np.pi/2, n_theta_bottom) + bin_size = (np.pi/2 - minimum_trapped_pitch)/n_theta_bottom + prob_theta_bottom = np.zeros(len(theta_bottoms_bin_centers)) + for i in range(len(theta_bottoms)): + for j in range(len(theta_bottoms[0])): + for k in range(len(theta_bottoms_bin_centers)): + if (theta_bottoms[i][j] >= theta_bottoms_bin_centers[k]-bin_size/2) and (theta_bottoms[i][j] < theta_bottoms_bin_centers[k]+bin_size/2): + prob_theta_bottom[k] += dist_of_theta_start_after_trapping(theta_start_array[i], trap_length, minimum_trapped_pitch, flat_fraction) + normalization = np.sum(prob_theta_bottom) + prob_theta_bottom = prob_theta_bottom/normalization #Is this the correct approach? + return theta_bottoms_bin_centers, prob_theta_bottom + +""" +figure = plt.figure() +theta_start_array = np.linspace(87*deg, np.pi/2, 1000) +prob_theta_start = dist_of_theta_start_after_trapping(theta_start_array, 4.05*m, 87*deg, flat_fraction=0.75) +plt.scatter(theta_start_array/deg, prob_theta_start) +plt.xlabel("Starting pitch angle $\\theta_{start}$ ($\degree$)", fontsize=14) +plt.ylabel("Probability (arb. units)", fontsize=14) +plt.savefig("test_theta_start_dist.png", dpi=300) +plt.show() +""" + + +# Trapping efficiency from axial field variation. +def trapping_efficiency(z_range, bg_magnetic_field, min_pitch_angle, trap_flat_fraction = 0.5): + + """ + Calculate the trapping efficiency for a given trap length and flat fraction. + + The trapping efficiency is computed using the formula: + epsilon(z) = sqrt(1 - B(z)/B_max(z)) + where B(z) is the magnetic field at position z, and B_max(z) is the maximum magnetic field along the z axis. + + Parameters + ---------- + z_range : float + The axial range (in z-direction, from trap center) over which electron trapping happens. + bg_magnetic_field : float + The background magnetic field. + min_pitch_angle : float + Minimum pitch angle to be trapped. + trap_flat_fraction : float, optional + Flat fraction of the trap. Default is 0.5. + + Returns + ------- + mean_efficiency : float + The mean trapping efficiency across the trap z-range. + + Notes + ----- + The magnetic field profile is computed using the `magnetic_field_flat_harmonic` function, currently it only produces z-profile of the trap without radial variation. + No radial field variation was assumed for this calculation. + The mean trapping efficiency is averaged over the region where the trapping field exists. + """ + + zs = np.linspace(-z_range, z_range, 500) + + profiles = [] + #Collect z profile of the magnetic field + for z in zs: + profiles.append(magnetic_field_flat_harmonic(z, bg_magnetic_field, z_range*2, min_pitch_angle, trap_flat_fraction)) + + #Calculate maximum trapping field along z (Bz_max) + maximum_Bz = max(profiles) + + #Calculate mean trapping efficiency using mean of epsilon(z) = sqrt(1-B(z)/B_max(z)) at z = 0 + mean_efficiency = np.mean(np.array([np.sqrt(1-b_at_z/maximum_Bz) for b_at_z in profiles])) + + return mean_efficiency + + +def track_duration_distribution(track_duration, mean_track_duration): + # The "1/mean_track_duration" factor normalizes the distribution + return (1/mean_track_duration)*np.exp(-track_duration/mean_track_duration) + + + +############################################################################### +class CavitySensitivity(Sensitivity): + """ + Documentation: + * Phase IV sensitivity document: https://www.overleaf.com/project/5de3e02edd267500011b8cc4 + * Talia's sensitivity script: https://3.basecamp.com/3700981/buckets/3107037/documents/2388170839 + * Nick's CRLB for frequency resolution: https://3.basecamp.com/3700981/buckets/3107037/uploads/2009854398 + * Molecular contamination in atomic tritium: https://3.basecamp.com/3700981/buckets/3107037/documents/3151077016 + """ + def __init__(self, config_path, verbose=True): + Sensitivity.__init__(self, config_path, verbose=verbose) + + # Calc non-config parameters outside of init function: + ## Allows re-calcing params if config values changed later, e.g. param scans + self.CalcDefaults(overwrite=False) + + # Add any additional initialization to this function, NOT __INIT__!! + def CalcDefaults(self, overwrite=False): + ### + #Initialization related to the effective volume: + ### + self.Jprime_0 = 3.8317 + self.cavity_freq = frequency(self.T_endpoint, self.MagneticField.nominal_field) + self.CavityRadius() + + #Get trap length from cavity length if not specified + if ((not hasattr(self.Experiment, 'trap_length')) or overwrite): + self.Experiment.trap_length = 0.8 * 2 * self.cavity_radius * self.Experiment.cavity_L_over_D + logger.info("Calc'd trap length: {} m".format(round(self.Experiment.trap_length/m, 3), 2)) + + self.Efficiency = NameSpace({opt: eval(self.cfg.get('Efficiency', opt)) for opt in self.cfg.options('Efficiency')}) + self.CavityVolume() + self.CavityPower() + + #Calculate position dependent trapping efficiency + self.pos_dependent_trapping_efficiency = trapping_efficiency( z_range = self.Experiment.trap_length /2, + bg_magnetic_field = self.MagneticField.nominal_field, + min_pitch_angle = self.FrequencyExtraction.minimum_angle_in_bandwidth, + trap_flat_fraction = self.MagneticField.trap_flat_fraction + ) + + #We may decide to remove the "Threshold" section and move the threshold-related parameters to the "Efficiency" section. + if not self.Efficiency.usefixedvalue: + self.Threshold = NameSpace({opt: eval(self.cfg.get('Threshold', opt)) for opt in self.cfg.options('Threshold')}) + + #Cyclotron radius is sometimes used in the effective volume calculation + self.cyc_rad = cyclotron_radius(self.cavity_freq, self.T_endpoint) + + #Assigning the background constant if it's not in the config file + if hasattr(self.Experiment, "bkgd_constant"): + self.bkgd_constant = self.Experiment.bkgd_constant + logger.info("Using background rate constant of {}/eV/s".format(self.bkgd_constant)) + else: + self.bkgd_constant = 1 + logger.info("Using background rate constant of 1/eV/s") + + # Need to get power fractions before calculating effective volume (given impact on detection efficiency) + # Power fractions are relative to the power of a 90° carrier electron + # If average_power_fractions==True, use carrier and sideband pitch power fractions averaged over the usable pitch angle range. + # If average_power_fractions==False, instead read in a csv file with power fractions vs. pitch angle from simulations. + # This file should have three columns: pitch angle (in degrees), carrier power, sideband power. + if hasattr(self.FrequencyExtraction, "use_average_power_fractions"): + if self.FrequencyExtraction.use_average_power_fractions: + logger.info("Using average carrier and sideband power fractions") + else: + logger.info("Using carrier and sideband (power fractions vs. pitch angle) from file") + # Read powers from the file and then scale them by power of maximum + # pitch angle to get power fractions. + theta_array, carrier_power_array, sideband_power_array = [], [], [] + power_file = open(self.FrequencyExtraction.powers_vs_theta_file, 'r') + for i in power_file.readlines()[1:]: # Skip header line + line = i.strip() + theta_array.append(float(line.split(",")[0])) # In degrees + carrier_power_array.append(float(line.split(",")[1])) + sideband_power_array.append(float(line.split(",")[2])) + power_file.close() + self.theta_array = np.array(theta_array)*deg # The "*deg" multiplies by np.pi/180 + carrier_power_array = np.array(carrier_power_array) + sideband_power_array = np.array(sideband_power_array) + # The calculation below assumes that the file contains a pitch angle very close to 90 degrees: + max_theta_index = np.argmax(self.theta_array) + self.carrier_power_fraction_array = carrier_power_array / carrier_power_array[max_theta_index] + self.sideband_power_fraction_array = sideband_power_array / carrier_power_array[max_theta_index] + + # Calculating distribution of pitch angles at the bottom of the trap, after trapping + theta_start_array = np.linspace(self.FrequencyExtraction.minimum_angle_in_bandwidth, np.pi/2, self.Efficiency.n_theta_start_for_trapped_pitch_dist) + self.theta_bottoms_bin_centers, self.prob_theta_bottom = dist_of_theta_bottom_after_trapping(self.MagneticField.nominal_field, theta_start_array, self.Experiment.trap_length, self.FrequencyExtraction.minimum_angle_in_bandwidth, flat_fraction=self.MagneticField.trap_flat_fraction, n_z_start=self.Efficiency.n_z_for_trapped_pitch_dist, n_theta_bottom=self.Efficiency.n_theta_bottom_for_trapped_pitch_dist) + + # Determining which probability corresponds to each pitch angle in self.theta_array + self.prob_theta_array = np.interp(self.theta_array, self.theta_bottoms_bin_centers, self.prob_theta_bottom) + + # Plotting pitch angle distribution + figure = plt.figure() + plt.scatter(self.theta_bottoms_bin_centers/deg, self.prob_theta_bottom, s=3, label="Binned distribution", color='red') + plt.scatter(self.theta_array/deg, self.prob_theta_array, s=2, marker='v', label="Interpolated to $\\theta_{bottom}$ values in simulations", color='blue') + plt.xlabel("Pitch angle at bottom of trap $\\theta_{bottom}$ ($\degree$)", fontsize=14) + plt.ylabel("Probability density", fontsize=14) + plt.legend(fontsize=12, loc='lower center') + plt.tight_layout() + plt.savefig("theta_bottom_dist_interpolated_{}.png".format(self.Experiment.exp_label), dpi=300) + + #### + #Initialization related to the energy resolution: + #### + #No longer using this CRLB_constant. If this change sticks, will remove it. + self.CRLB_constant = 6 + if hasattr(self.FrequencyExtraction, "crlb_constant"): + self.CRLB_constant = self.FrequencyExtraction.crlb_constant + logger.info("Using configured CRLB constant") + + # Number of steps in pitch angle between min_pitch and pi/2 for the frequency noise uncertainty calculation + self.pitch_steps = 100 + if hasattr(self.FrequencyExtraction, "pitch_steps"): + self.pitch_steps = self.FrequencyExtraction.pitch_steps + logger.info("Using configured pitch_steps value") + + + #Calculate the resolution contribution from noise (affecting frequency reconstruction), + #accounting for a magnetic field correction based on sideband frequencies (if that option is chosen). + #This needs to be done before calculating the effective volume, since the reconstruction efficiency + #is a factor in the effective volume. + self.syst_frequency_extraction() + + ### + #Back to effective volume: + ### + #Calculate the effective volume and print out related quantities + self.EffectiveVolume() + logger.info("Trap radius: {} cm".format(round(self.cavity_radius/cm, 3), 2)) + logger.info("Total trap volume: {} m^3".format(self.total_trap_volume/m**3)) + logger.info("Cyclotron radius: {}m".format(self.cyc_rad/m)) + if self.use_cyc_rad: + logger.info("Using cyclotron radius as unusable distance from wall, for radial efficiency calculation") + + #Just calculated for comparison + self.larmor_power = rad_power(self.T_endpoint, np.pi/2, self.MagneticField.nominal_field) # currently not used + + # Determining whether to use a fixed detection efficiency or calculate it from the detection threshold + if not self.Efficiency.usefixedvalue: + if self.Threshold.use_detection_threshold: + logger.info("Overriding any detection eff and RF background in the config file; calculating these from the detection_threshold.") + else: + logger.info("Using the detection eff and RF background rate from the config file.") + + + # CAVITY + def CavityRadius(self): + axial_mode_index = 1 + self.cavity_radius = c0/(2*np.pi*self.cavity_freq)*np.sqrt(self.Jprime_0**2+axial_mode_index**2*np.pi**2/(4*self.Experiment.cavity_L_over_D**2)) + return self.cavity_radius + + def CavityVolume(self): + #radius = 0.5*wavelength(self.T_endpoint, self.MagneticField.nominal_field) + self.total_cavity_volume = 2*self.cavity_radius*self.Experiment.cavity_L_over_D*np.pi*(self.cavity_radius)**2*self.Experiment.n_cavities + + logger.info("Frequency: {} MHz".format(round(self.cavity_freq/MHz, 3))) + logger.info("Wavelength: {} cm".format(round(wavelength(self.T_endpoint, self.MagneticField.nominal_field)/cm, 3))) + logger.info("Cavity radius: {} cm".format(round(self.cavity_radius/cm, 3))) + logger.info("Cavity length: {} cm".format(round(2*self.cavity_radius*self.Experiment.cavity_L_over_D/cm, 3))) + logger.info("Total cavity volume: {} m^3".format(round(self.total_cavity_volume/m**3, 3)))\ + + return self.total_cavity_volume + + + # ELECTRON TRAP + def TrapVolume(self): + # Total volume of the electron traps in all cavities + self.total_trap_volume = self.Experiment.trap_length*np.pi*(self.cavity_radius)**2*self.Experiment.n_cavities + return self.total_trap_volume + + + + def EffectiveVolume(self): + self.total_trap_volume = self.TrapVolume() + + if self.Efficiency.usefixedvalue: + self.effective_volume = self.total_trap_volume * self.Efficiency.fixed_efficiency + self.use_cyc_rad = False + self.RF_background_rate_per_eV = self.Experiment.RF_background_rate_per_eV + else: + #Detection efficiency + if self.Threshold.use_detection_threshold: + #Calculating the detection efficiency given the SNR of data and the threshold. + #If the config file contains a detection effciency or RF background rate, they are overridden. + self.assign_background_rate_from_threshold() + self.assign_detection_efficiency_from_threshold() + else: + #Using the inputted detection efficiency and RF background rate from the config file. + self.detection_efficiency = self.Efficiency.detection_efficiency + self.RF_background_rate_per_eV = self.Experiment.RF_background_rate_per_eV + + + #Radial efficiency + if self.Efficiency.unusable_dist_from_wall >= self.cyc_rad: + self.radial_efficiency = (self.cavity_radius - self.Efficiency.unusable_dist_from_wall)**2/self.cavity_radius**2 + self.use_cyc_rad = False + else: + self.radial_efficiency = (self.cavity_radius - self.cyc_rad)**2/self.cavity_radius**2 + self.use_cyc_rad = True + + #Efficiency from a cut during analysis on the axial frequency + self.fa_cut_efficiency = trapping_efficiency(z_range = self.Experiment.trap_length /2, + bg_magnetic_field = self.MagneticField.nominal_field, + min_pitch_angle = self.Efficiency.min_pitch_used_in_analysis, + trap_flat_fraction = self.MagneticField.trap_flat_fraction + )/self.pos_dependent_trapping_efficiency + + #The effective volume includes the three efficiency factors above, as well as the trapping efficiency + self.effective_volume = self.total_trap_volume*self.radial_efficiency*self.detection_efficiency*self.fa_cut_efficiency*self.pos_dependent_trapping_efficiency*self.recon_efficiency + + # The "signal rate improvement" factor can be toggled to test the increase in statistics required to reach some sensitivity + self.effective_volume*=self.Experiment.sri_factor + return self.effective_volume + + + def BoxTrappingEfficiency(self): + self.box_trapping_efficiency = np.cos(self.FrequencyExtraction.minimum_angle_in_bandwidth) + return self.box_trapping_efficiency + + def TrapLength(self): + self.Experiment.trap_length = 0.8 * 2 * self.cavity_radius * self.Experiment.cavity_L_over_D + logger.info("Calc'd trap length: {} m".format(round(self.Experiment.trap_length/m, 3), 2)) + + def CavityPower(self): + #Jprime_0 = 3.8317 + max_ax_freq, mean_field, z_t = axial_motion(self.MagneticField.nominal_field, + self.FrequencyExtraction.minimum_angle_in_bandwidth, + self.Experiment.trap_length, + self.FrequencyExtraction.minimum_angle_in_bandwidth, + self.T_endpoint, flat_fraction=self.MagneticField.trap_flat_fraction, trajectory = 1000) #1000 + + #The np.random.triangluar function weights the radii, accounting for the fact that there are more electrons at large radii than small ones + r_sample_size = 50 + if((not self.Efficiency.calculate_det_eff_for_sampled_radii) or (self.Efficiency.usefixedvalue)): r_sample_size = 1000 + power_vs_r_with_zeros = np.mean(larmor_orbit_averaged_hanneke_power(np.random.triangular(0, self.cavity_radius, self.cavity_radius, size=r_sample_size), + z_t, self.CavityLoadedQ(), + 2*self.Experiment.cavity_L_over_D*self.cavity_radius, + self.cavity_radius, + self.cavity_freq), axis=1) + #Remove zeros, since these represent electrons that hit the cavity wall and are not detected + self.signal_power_vs_r = power_vs_r_with_zeros[power_vs_r_with_zeros != 0] + + self.signal_power = np.mean(self.signal_power_vs_r) + return self.signal_power + + + def CavityLoadedQ(self): + # Using Wouter's calculation: + # Total required bandwidth is the sum of the endpoint region and the axial frequency. + # I will assume the bandwidth is dominated by the sidebands and not by the energy ROI + + #self.loaded_q =1/(0.22800*((90-self.FrequencyExtraction.minimum_angle_in_bandwidth)*np.pi/180)**2+2**2*0.01076**2/(4*0.22800)) + + endpoint_frequency = self.cavity_freq + #required_bw_axialfrequency = axial_frequency(self.Experiment.cavity_L_over_D*self.CavityRadius()*2, + # self.T_endpoint, + # self.FrequencyExtraction.minimum_angle_in_bandwidth/deg) + max_ax_freq, mean_field, _ = axial_motion(self.MagneticField.nominal_field, + self.FrequencyExtraction.minimum_angle_in_bandwidth, + self.Experiment.trap_length, + self.FrequencyExtraction.minimum_angle_in_bandwidth, + self.T_endpoint, flat_fraction=self.MagneticField.trap_flat_fraction) + required_bw_axialfrequency = max_ax_freq*self.FrequencyExtraction.sideband_order + self.required_bw_axialfrequency = required_bw_axialfrequency + required_bw_meanfield = required_bw_meanfield = np.abs(frequency(self.T_endpoint, mean_field) - endpoint_frequency) + required_bw = np.add(required_bw_axialfrequency,required_bw_meanfield) # Broadcasting + self.required_bw = required_bw + + # Cavity coupling + self.loaded_q = endpoint_frequency/required_bw # FWHM + return self.loaded_q + + # SENSITIVITY + # see parent class in SensitivityFormulas.py + + + # SYSTEMATICS + # Generic systematics are implemented in the parent class in SensitivityFormulas.py + + def calculate_tau_snr(self, time_window, power_fraction=1, tau_snr_array_for_radii=False): + """ + power_fraction may be used as a carrier or a sideband power fraction, + relative to the power of a 90 degree carrier. + """ + endpoint_frequency = self.cavity_freq + + # Cavity coupling + self.CavityLoadedQ() + coupling = self.FrequencyExtraction.unloaded_q/self.loaded_q-1 + + # Attenuation frequency dependence at hoc method + #att_cir_db = -0.3 + #att_line_db = -0.05 + att_cir_db_freq = self.FrequencyExtraction.att_cir_db*(1+endpoint_frequency/(10*GHz)) + att_line_db_freq = self.FrequencyExtraction.att_line_db*(1+endpoint_frequency/(10*GHz)) + + # Noise power for bandwidth set by density/track length + fft_bandwidth = 3/time_window #(delta f) is the frequency bandwidth of interest. We have a main carrier and 2 axial side bands, so 3*(FFT bin width) + self.fft_bandwidth = fft_bandwidth + tn_fft = Pn_dut_entrance(self.FrequencyExtraction.cavity_temperature, + self.FrequencyExtraction.amplifier_temperature, + att_line_db_freq,att_cir_db_freq, + coupling, + endpoint_frequency, + fft_bandwidth,self.loaded_q)/kB/fft_bandwidth + + # Noise temperature of amplifier + tn_amplifier = endpoint_frequency*hbar*2*np.pi/kB/self.FrequencyExtraction.quantum_amp_efficiency + tn_system_fft = tn_amplifier+tn_fft + self.noise_temp = tn_system_fft + + # Pe = rad_power(self.T_endpoint, self.FrequencyExtraction.pitch_angle, self.MagneticField.nominal_field) + # logger.info("Power: {}".format(Pe/W)) + if tau_snr_array_for_radii: + Pe = self.signal_power_vs_r * power_fraction + else: + Pe = self.signal_power * power_fraction + + P_signal_received = Pe*db_to_pwr_ratio(att_cir_db_freq+att_line_db_freq) + self.received_power = P_signal_received + tau_snr = kB*tn_system_fft/P_signal_received + self.noise_energy = kB*tn_system_fft + + # end of Wouter's calculation + return tau_snr + + """ + def print_SNRs(self, rho_opt): + tau_snr = self.calculate_tau_snr(self.time_window, sideband_power_fraction=1) + logger.info("tau_SNR: {}s".format(tau_snr/s)) + eV_bandwidth = np.abs(frequency(self.T_endpoint, self.MagneticField.nominal_field) - frequency(self.T_endpoint + 1*eV, self.MagneticField.nominal_field)) + SNR_1eV = 1/eV_bandwidth/tau_snr + track_duration = track_length(rho_opt, self.T_endpoint, molecular=(not self.Experiment.atomic)) + SNR_track_duration = track_duration/tau_snr + SNR_1ms = 0.001*s/tau_snr + logger.info("SNR for 1eV bandwidth: {}".format(SNR_1eV)) + logger.info("SNR 1 eV from temperatures:{}".format(self.received_power/(self.noise_energy*eV_bandwidth))) + logger.info("Track duration: {}ms".format(track_duration/ms)) + logger.info("Sampling duration for 1eV: {}ms".format(1/eV_bandwidth/ms)) + logger.info("SNR for track duration: {}".format(SNR_track_duration)) + logger.info("SNR for 1 ms: {}".format(SNR_1ms)) + logger.info("Received power: {}W".format(self.received_power/W)) + logger.info("Noise power in 1eV: {}W".format(self.noise_energy*eV_bandwidth/W)) + logger.info("Noise temperature: {}K".format(self.noise_temp/K)) + logger.info("Opimtum energy window: {} eV".format(self.DeltaEWidth()/eV)) + """ + + def frequency_variance_from_CRLB(self, tau_SNR, track_duration): + self.eta = self.slope*track_duration/(4*self.cavity_freq*np.pi) + if self.eta < 1e-6: + # This is for the case where the track is flat (almost no slope), and where we + # treat it as a pure sinusoid (don't fit the slope when extracting the frequency). + # Applies for a complex signal. + return self.FrequencyExtraction.CRLB_scaling_factor*(6*tau_SNR/track_duration**3)/(2*np.pi)**2 + else: + # Non-zero, fitted slope. + # Doesn't assume that alpha*T/2 << omega_c, since it includes the 5*eta/(1-eta) term in Eq. 25 of Joe's write-up: https://3.basecamp.com/3700981/buckets/3107037/documents/6331876030. + # CODE IMPLEMENTATION NEEDS TO BE DOUBLE-CHECKED BY CONSIDERING AN EXPERIMENT WITH LARGE-ISH ETA. + # The first term relies on the relation delta_t_start = sqrt(20)*tau_snr. This is from Equation 6.40 of Nick's thesis, + # derived in Appendix A and verified with an MC study. + # Using a factor of 23 instead of 20, from re-calculating Nick's integrals (though this derivation is approximate). + # Nick's derivation uses an expression for P_fa assuming the phase is known. + # The phase won't be known, but it's more difficult to determine the unknown-phase expression. + # Working on that. + return self.FrequencyExtraction.CRLB_scaling_factor*(23*(self.slope*tau_SNR)**2 + tau_SNR/track_duration**3*(96 - 6*5*self.eta/(1+self.eta)))/(2*np.pi)**2 + + + def syst_frequency_extraction(self): + # cite{https://3.basecamp.com/3700981/buckets/3107037/uploads/2009854398} (Section 1.2, p 7-9) + # Are we double counting the antenna collection efficiency? We use it here. Does it also impact the effective volume, v_eff ? + + if self.FrequencyExtraction.UseFixedValue: + sigma = self.FrequencyExtraction.Default_Systematic_Smearing + delta = self.FrequencyExtraction.Default_Systematic_Uncertainty + return sigma, delta + + #Sample track duration values from exponential distribution + self.mean_track_duration = track_length(self.Experiment.number_density, self.T_endpoint, molecular=(not self.Experiment.atomic)) + if hasattr(self.FrequencyExtraction, "use_full_track_duration_distribution") and self.FrequencyExtraction.use_full_track_duration_distribution: + n_track_duration_samples = 5000 + self.time_window = np.random.exponential(self.mean_track_duration/s, n_track_duration_samples)*s + #Order self.time_window from longest to shortest track durations. That will later allow us to + #break a loop when none of the resolution from noise values (for different thetas) will be below + #the threshold anymore + self.time_window = np.sort(self.time_window)[::-1] + else: + self.time_window = self.mean_track_duration*np.ones(1) + + endpoint_frequency = self.cavity_freq + + # Converting the max allowed energy reconstruction stdev to a max allowed frequency variance + self.max_allowed_f_var = (self.FrequencyExtraction.max_allowed_energy_std_recon*(2*np.pi*endpoint_frequency**2)/e/self.MagneticField.nominal_field/c0**2)**2 + + # using Pe and alpha (aka slope) from above + Pe = self.signal_power #/self.FrequencyExtraction.mode_coupling_efficiency + + self.slope = endpoint_frequency * 2 * np.pi * Pe/me/c0**2 # track slope + + #self.time_window_slope_zero = abs(self.cavity_freq-frequency(self.T_endpoint+20*meV, self.MagneticField.nominal_field))/self.slope + + #tau_snr_full_length = [] + self.var_f_c_CRLB = [] + if self.FrequencyExtraction.use_average_power_fractions: + for window in self.time_window: + tau_snr_full_length = self.calculate_tau_snr(window, self.FrequencyExtraction.carrier_power_fraction) + #Calculate the frequency variance from the CRLB + self.var_f_c_CRLB.append(self.frequency_variance_from_CRLB(tau_snr_full_length, window)) + else: + for window in self.time_window: + tau_snr_full_length = self.calculate_tau_snr(window, self.carrier_power_fraction_array)[:len(self.theta_array)-1] #Cut out theta=pi/2, since sideband power is 0 there, resulting in infinite tau_snr. + #Calculate the frequency variance from the CRLB + self.var_f_c_CRLB.append(self.frequency_variance_from_CRLB(tau_snr_full_length, window)) + + # sigma_f from pitch angle reconstruction + if self.FrequencyExtraction.crlb_on_sidebands: + #Calculate noise contribution to uncertainty, including energy correction for pitch angle. + #This comes from section 6.1.9 of the CDR. + + tau_snr_full_length_sideband = [] + if self.FrequencyExtraction.use_average_power_fractions: + for window in self.time_window: + tau_snr_full_length_sideband.append(self.calculate_tau_snr(window, self.FrequencyExtraction.sideband_power_fraction)) + else: + for window in self.time_window: + tau_snr_full_length_sideband_temp = self.calculate_tau_snr(window, self.sideband_power_fraction_array) + tau_snr_full_length_sideband.append(tau_snr_full_length_sideband_temp[:len(self.theta_array)-1]) #Cut out theta=pi/2, since sideband power is 0 there, resulting in infinite tau_snr. + + # Defining variables trhat are the same regardless of track duration: + + m = self.FrequencyExtraction.sideband_order + + # Defining array of pitch angle complements (pi/2 - theta) used when calculating + # the parameters describing the track shape (p and q) + thetas_for_p_and_q_calc = np.linspace(self.FrequencyExtraction.minimum_angle_in_bandwidth, 90*deg, self.pitch_steps) + pitch_comps_for_p_and_q_calc = np.pi/2 - thetas_for_p_and_q_calc + + # Defining array of pitch angle complement values over which we calculate the + # resolution contribution from noise. + if self.FrequencyExtraction.use_average_power_fractions: + pitch_comps = pitch_comps_for_p_and_q_calc + else: + pitch_comps = np.pi/2 - self.theta_array[:len(self.theta_array)-1] #Cut out theta=pi/2 since sideband power is 0 there, resulting in infinite tau_snr. + + #Define the trap parameter p based on the relation between the trap length and the cavity mode + #This p is for a box trap + self.p_box = np.pi*beta(self.T_endpoint)*self.cavity_radius/self.Jprime_0/self.Experiment.trap_length + + #Now find p for the actual trap that we have + #Using the average p across the pitch angle range + ax_freq_array, mean_field_array, z_t = axial_motion(self.MagneticField.nominal_field, + thetas_for_p_and_q_calc, self.Experiment.trap_length, + self.FrequencyExtraction.minimum_angle_in_bandwidth, + self.T_endpoint, flat_fraction=self.MagneticField.trap_flat_fraction) + fc0_endpoint = self.cavity_freq + p_array = ax_freq_array/fc0_endpoint/pitch_comps_for_p_and_q_calc #An array + if self.FrequencyExtraction.use_average_power_fractions: + p_array = p_array[:1] #Cut out theta=pi/2 (ill defined there) + self.p = np.mean(p_array) + + # Now calculating q for the trap that we have + # Using the q for the minimum trapped pitch angle + fc_endpoint_min_theta = frequency(self.T_endpoint, mean_field_array[0]) + self.q = (fc_endpoint_min_theta/fc0_endpoint - 1)/(pitch_comps_for_p_and_q_calc[0])**2 + + # Derivative of f_c0 (frequency corrected to B-field at bottom of the trap) with respect to f_c + dfc0_dfc_array = 0.5*(1 - (1 - 4*self.q*pitch_comps/m/self.p + self.q*pitch_comps**2)/(1 - self.q*pitch_comps**2)) + + # Derivative of f_c0 with respect to f_lsb (lower sideband frequency) + dfc0_dlsb_array = 0.5 - 2*self.q*pitch_comps/m/self.p/(1 - self.q*pitch_comps**2) + + #Cut out theta=pi/2 since sideband power is 0 there, resulting in infinite tau_snr. + prob_theta_array_without_pi_over_2 = self.prob_theta_array[:len(self.theta_array)-1] + + #Now loop over track durations + recon_eff, sigma_f_noise = [], [] + for i in range(len(self.time_window)): + # (sigmaf_lsb)^2: + var_f_sideband_crlb = self.frequency_variance_from_CRLB(tau_snr_full_length_sideband[i], window) + + # Noise variance term from the carrier frequency uncertainty + var_noise_from_fc_array = dfc0_dfc_array**2*self.var_f_c_CRLB[i] + + # Noise variance term from the lower sideband frequency uncertainty + var_noise_from_flsb_array = dfc0_dlsb_array**2*var_f_sideband_crlb + + # Total uncertainty for each pitch angle + var_f_noise_array = var_noise_from_fc_array + var_noise_from_flsb_array + + #Cut out values of var_f_noise_array that are larger than max_allowed_f_var + valid_var_indices = var_f_noise_array <= self.max_allowed_f_var + + #Break the loop if no valid var_f_noise_array values remain + #Since we sorted the track durations from longest to shortest, + #no valid values will remian for later (shorter) track durations. + if np.sum(valid_var_indices) == 0: + break + + valid_vars = var_f_noise_array[valid_var_indices] + # Weight the reconstruction efficiency by the pitch-angle probability density + total_prob = np.sum(prob_theta_array_without_pi_over_2) + recon_eff.append(np.sum(prob_theta_array_without_pi_over_2[valid_var_indices]) / total_prob) + + #Choose the values of prob_theta_array_without_pi_over_2 that correspond to the valid_vars + prob_theta_array_temp = prob_theta_array_without_pi_over_2[valid_var_indices] + + # Next, we average over sigma_noise values. + # This is a quadrature sum average weighted by the pitch angle distribution, + # reflecting that the detector response function could be constructed by sampling + # from many normal distributions with different standard deviations (sigma_noise_array), + # then finding the standard deviation of the full group of sampled values. + # IS THE BELOW CORRECT? + sigma_f_noise.append(np.sqrt(np.sum(valid_vars*prob_theta_array_temp)/np.sum(prob_theta_array_temp))) + + #Average over track durations + self.recon_efficiency = np.mean(recon_eff) + self.sigma_f_noise = np.mean(sigma_f_noise) + + else: + #Cut out values of self.var_f_c_CRLB that are larger than self.max_allowed_f_var. + #Average over the remaining values. + #Save the percentage of values that were not cut out as self.recon_efficiency + valid_var_indices = self.var_f_c_CRLB <= self.max_allowed_f_var + valid_vars = self.var_f_c_CRLB[valid_var_indices] + self.recon_efficiency = len(valid_vars)/len(self.var_f_c_CRLB) + self.sigma_f_noise = np.sqrt(np.mean(valid_vars)) + + # Convert uncertainty from frequency to energy + self.sigma_K_noise = e*self.MagneticField.nominal_field/(2*np.pi*endpoint_frequency**2)*self.sigma_f_noise*c0**2 + + + # combined sigma_f in eV + sigma_f = np.sqrt(self.sigma_K_noise**2 + self.FrequencyExtraction.magnetic_field_smearing**2) + # delta_sigma_f = np.sqrt((delta_sigma_K_f_CRLB**2 + self.FrequencyExtraction.magnetic_field_uncertainty**2)/2) + if self.FrequencyExtraction.usefixeduncertainty: + return sigma_f, self.FrequencyExtraction.fixed_relativ_uncertainty*sigma_f + else: + raise NotImplementedError("Uncertainty on CRLB for cavity noise calculation is not implemented.") + + def syst_magnetic_field(self): + """ + Magnetic field uncertanty is in principle generic but its impact on efficiency depends on reconstruction and therefore on detector technology. + """ + # magnetic field uncertainties can be decomposed in several part + # * true magnetic field inhomogeneity + # (would be there also without a trap) + # * magnetic field calibration has uncertainties + # (would be there also without a trap) + # * position / pitch angle reconstruction has uncertainties + # (this can even be the degenerancy we see for harmonic traps) + # (depends on trap shape) + + if self.MagneticField.UseFixedValue: + sigma = self.MagneticField.Default_Systematic_Smearing + delta = self.MagneticField.Default_Systematic_Uncertainty + return sigma, delta + + B = self.MagneticField.nominal_field + if self.MagneticField.useinhomogeneity: + frac_uncertainty = self.MagneticField.fraction_uncertainty_on_field_broadening + sigma_meanB = self.MagneticField.sigma_meanb + sigmaE_meanB = self.BToKeErr(sigma_meanB*B, B) + sigmaE_r = self.MagneticField.sigmae_r + sigmaE_theta = self.MagneticField.sigmae_theta + sigmaE_phi = self.MagneticField.sigmae_phi + sigma = np.sqrt(sigmaE_meanB**2 + sigmaE_r**2 + sigmaE_theta**2 + sigmaE_phi**2) + return sigma, frac_uncertainty*sigma + else: + return 0, 0 + + def det_efficiency_track_duration(self): + """ + Detection efficiency implemented based on René's slides, with faster and stable implementation using Gauss-Laguerre quadrature (G-L method): + https://3.basecamp.com/3700981/buckets/3107037/documents/8013439062 + Gauss-Laguerre Quadrature: https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature + + The following changes were made to the original integral to fit the G-L method: + Original integrand: ∫[0 to \inf] ncx2(df=2, nc=t/τ).sf(thres) * (1/μ) * exp(-t/μ) dt + + Where. t = track_duration, μ (\mu) = mean_track_duration, τ (\tau) = tau_snr_ex_carrier, thres = detection_threshold + + We do the change of variable, x = t / μ. So, t = x μ, or, dt = μ dx + + Substituting into the original integral: + ∫[0 to \inf] ncx2(df=2, nc=xμ/τ).sf(thres) * (1/μ) * exp(-x) μ dx + The μ's cancel out, and the integral takes the form: + ∫[0 to \inf] f(x) * exp(-x) dx + where, f(x) = ncx2(df=2, nc=xμ/τ).sf(thres) + + Parameters: None + + Returns: avg_efficiency (float): SNR and threshold dependent detection efficieny. + + Notes: Also check the antenna paper for more details. Especially the section on the signal detection with matched filtering. + """ + # Calculate the mean track duration + # TO-DO: Only do the lines below ones for a given density; don't repeat for each threshold being scanned ... + #mean_track_duration = track_length(self.Experiment.number_density, self.T_endpoint, molecular=(not self.Experiment.atomic)) + if self.FrequencyExtraction.use_average_power_fractions: + tau_snr_ex_total = self.calculate_tau_snr(self.mean_track_duration, self.FrequencyExtraction.carrier_power_fraction + self.FrequencyExtraction.sideband_power_fraction, tau_snr_array_for_radii=self.Efficiency.calculate_det_eff_for_sampled_radii) + else: + tau_snr_ex_total = self.calculate_tau_snr(self.mean_track_duration, self.carrier_power_fraction_array + self.sideband_power_fraction_array, tau_snr_array_for_radii=self.Efficiency.calculate_det_eff_for_sampled_radii) + if isinstance(tau_snr_ex_total, float): + tau_snr_ex_total = [tau_snr_ex_total] + + # Roots and weights for the Laguerre polynomial + x, w = roots_laguerre(100) #n=100 is the number of quadrature points + + # Scale the track duration to match the form of Gauss-Laguerre quadrature + scaled_x = x * self.mean_track_duration # scaled_x = xμ + + # Evaluate the non-central chi-squared dist values at the scaled quadrature points + sf_values = np.array([ncx2(df=2, nc=2 * scaled_x / tau_snr).sf(self.Threshold.detection_threshold) for tau_snr in tau_snr_ex_total]) + + # Calculate and return the integration result from weighted sum + eff_for_each_r_and_theta = np.sum(w * sf_values, axis=1) + + # Average efficiencies over the sampled electron radii and pitch angles. + # Calculation below accounts for trapped pitch angle distribution (self.prob_theta_array). + # Weighting for radial distribution is accounted for in sampling, earlier. + if self.FrequencyExtraction.use_average_power_fractions: + avg_efficiency = np.mean(eff_for_each_r_and_theta) + else: + if not self.Efficiency.calculate_det_eff_for_sampled_radii: + avg_efficiency = np.sum(self.prob_theta_array * eff_for_each_r_and_theta)/sum(self.prob_theta_array) + else: + #Sum over radii with equal weights, and sum over pitch angles with probability weights + #I'm not sure if I get the axes right, below. + avg_efficiency = np.sum(eff_for_each_r_and_theta, axis=0)/len(self.signal_power_vs_r) + avg_efficiency = np.sum(self.prob_theta_array * avg_efficiency)/sum(self.prob_theta_array) + return avg_efficiency + + def assign_detection_efficiency_from_threshold(self): + self.detection_efficiency = self.det_efficiency_track_duration() + return self.detection_efficiency + + def rf_background_rate_cavity(self): + # Detection efficiency implemented based on René's slides + # https://3.basecamp.com/3700981/buckets/3107037/documents/8013439062 + # Also check the antenna paper for more details, especially the section + # on the signal detection with matched filtering. + # The background constant will need to be determined from Monte Carlo simulations. + return chi2(df=2).sf(self.Threshold.detection_threshold)*self.bkgd_constant/(eV*s) + + def assign_background_rate_from_threshold(self): + self.RF_background_rate_per_eV = self.rf_background_rate_cavity() + return self.RF_background_rate_per_eV + + + + + # PRINTS + def print_SNRs(self, rho=None): + #logger.warning("Deprecation warning: This function does not modify the number density in the Experiment namespace. Values printed are for pre-set number density.") + + logger.info("**SNR parameters**:") + if rho == None: + mean_track_duration = self.mean_track_duration + logger.info("SNR-related parameters are printed for pre-set number density.") + else: + mean_track_duration = track_length(rho, self.T_endpoint, molecular=(not self.Experiment.atomic)) + + tau_snr_90deg = self.calculate_tau_snr(mean_track_duration, power_fraction=1) + #For an example carrier: + if self.FrequencyExtraction.use_average_power_fractions: + tau_snr_ex_carrier = self.calculate_tau_snr(mean_track_duration, self.FrequencyExtraction.carrier_power_fraction) + else: + tau_snr_ex_carrier = np.mean(self.calculate_tau_snr(mean_track_duration, self.carrier_power_fraction_array)) + + eV_bandwidth = np.abs(self.cavity_freq - frequency(self.T_endpoint + 1*eV, self.MagneticField.nominal_field)) + SNR_1eV_90deg = 1/eV_bandwidth/tau_snr_90deg + SNR_track_duration_90deg = mean_track_duration/tau_snr_90deg + SNR_1ms_90deg = 0.001*s/tau_snr_90deg + + SNR_1eV_ex_carrier = 1/eV_bandwidth/tau_snr_ex_carrier + SNR_track_duration_ex_carrier = mean_track_duration/tau_snr_ex_carrier + SNR_1ms_ex_carrier = 0.001*s/tau_snr_ex_carrier + + logger.info("Number density: {} m^-3".format(self.Experiment.number_density*m**3)) + logger.info("Track duration: {}ms".format(mean_track_duration/ms)) + logger.info("tau_SNR for 90° carrier: {}s".format(tau_snr_90deg/s)) + logger.info("tau_SNR for carrier used in calculation (see config file): {}s".format(tau_snr_ex_carrier/s)) + logger.info("Sampling duration for 1eV: {}ms".format(1/eV_bandwidth/ms)) + + logger.info("Received power for 90° carrier: {}W".format(self.received_power/W)) + logger.info("Noise temperature: {}K".format(self.noise_temp/K)) + logger.info("Noise power in 1eV: {}W".format(self.noise_energy*eV_bandwidth/W)) + logger.info("SNRs of carriers (90°, used in calc) for 1eV bandwidth: {}, {}".format(SNR_1eV_90deg, SNR_1eV_ex_carrier)) + #logger.info("SNR 1 eV from temperatures:{}".format(self.received_power/(self.noise_energy*eV_bandwidth))) + logger.info("SNRs of carriers (90°, used in calc) for track duration at optimum density: {}, {}".format(SNR_track_duration_90deg, SNR_track_duration_ex_carrier)) + logger.info("SNR of carriers (90°, used in calc) for 1 ms: {}, {}".format(SNR_1ms_90deg, SNR_1ms_ex_carrier)) + + + logger.info("Optimum energy window: {} eV".format(self.DeltaEWidth()/eV)) + + #logger.info("CRLB if slope is nonzero and needs to be fitted: {} Hz".format(np.sqrt(self.var_f_CRLB_slope_fitted)/Hz)) + #logger.info("CRLB constant: {}".format(self.CRLB_constant)) + logger.info("**Done printing SNR parameters.**") + + return self.noise_temp, SNR_1eV_90deg, mean_track_duration + + + def print_Efficiencies(self): + + logger.info("Effective volume: {} mm^3".format(round(self.effective_volume/mm**3, 3))) + logger.info("Total efficiency: {}".format(self.effective_volume/self.total_trap_volume)) + + if not self.Efficiency.usefixedvalue: + # radial and detection efficiency are configured in the config file + logger.info("Reconstruction efficiency: {}".format(self.recon_efficiency)) + logger.info("Radial efficiency: {}".format(self.radial_efficiency)) + logger.info("Detection efficiency: {}".format(self.detection_efficiency)) + #logger.info("Detection efficiency integration error: {}".format(self.abs_err)) + logger.info("Trapping efficiency: {}".format(self.pos_dependent_trapping_efficiency)) + logger.info("Efficiency from axial frequency cut: {}".format(self.fa_cut_efficiency)) + logger.info("SRI factor: {}".format(self.Experiment.sri_factor)) + + + + +""" # Cramer-Rao lower bound / how much worse are we than the lower bound +ScalingFactorCRLB = self.FrequencyExtraction.CRLB_scaling_factor +ts = self.FrequencyExtraction.track_timestep +# "This is apparent in the case of resonant patch antennas and cavities, in which the time scale of the signal onset is set by the Q-factor of the resonant structure." +# You can get it from the finite impulse response of the antennas from HFSS +Gdot = self.FrequencyExtraction.track_onset_rate + +fEndpoint = frequency(self.T_endpoint, self.MagneticField.nominal_field) +betae = beta(self.T_endpoint) +Pe = rad_power(self.T_endpoint, self.FrequencyExtraction.pitch_angle, self.MagneticField.nominal_field) +alpha_approx = fEndpoint * 2 * np.pi * Pe/me/c0**2 # track slope +# quantum limited noise +sigNoise = np.sqrt((2*pi*fEndpoint*hbar*self.FrequencyExtraction.amplifier_noise_scaling+kB*self.FrequencyExtraction.antenna_noise_temperature)/ts) # noise level +Amplitude = np.sqrt(self.FrequencyExtraction.epsilon_collection*Pe) +Nsteps = 1 / (self.Experiment.number_density * self.Te_crosssection*betae*c0*ts) # Number of timesteps of length ts + +# sigma_f from Cramer-Rao lower bound in Hz +sigma_f_CRLB = (ScalingFactorCRLB /(2*np.pi) * sigNoise/Amplitude * np.sqrt(alpha_approx**2/(2*Gdot) ++ 96.*Nsteps/(ts**2*(Nsteps**4-5*Nsteps**2+4))))""" + + +""" +CRLB_constant = 6 +sigma_CRLB_slope_zero = np.sqrt((CRLB_constant*tau_snr_part_length/self.time_window_slope_zero**3)/(2*np.pi)**2)*self.FrequencyExtraction.CRLB_scaling_factor + +sigma_f_CRLB = np.min([sigma_CRLB_slope_zero, sigma_f_CRLB_slope_fitted]) + +# logger.info("CRLB options are: {} , {}".format(sigma_CRLB_slope_zero/Hz, sigma_f_CRLB_slope_fitted/Hz)) +self.best_time_window=[self.time_window_slope_zero, self.time_window][np.argmin([sigma_CRLB_slope_zero, sigma_f_CRLB_slope_fitted])]""" + +"""# uncertainty in alpha +delta_alpha = 6*sigNoise/(Amplitude*ts**2) * np.sqrt(10/(Nsteps*(Nsteps**4-5*Nsteps**2+4))) +# uncetainty in sigma_f in Hz due to uncertainty in alpha +delta_sigma_f_CRLB = delta_alpha * alpha_approx *sigNoise**2/(8*np.pi**2*Amplitude**2*Gdot*sigma_f_CRLB*ScalingFactorCRLB**2)""" + + +"""fc_endpoint_array = frequency(self.T_endpoint, mean_field_array) +self.q_array = 1/phis**2*(fc_endpoint_array/fc0_endpoint - 1) +self.q = np.mean(self.q_array[1:])""" diff --git a/mermithid/sensitivity/SensitivityFormulas.py b/mermithid/sensitivity/SensitivityFormulas.py new file mode 100755 index 00000000..7747cac0 --- /dev/null +++ b/mermithid/sensitivity/SensitivityFormulas.py @@ -0,0 +1,433 @@ +''' +Class calculating neutrino mass sensitivities based on analytic formulas from CDR. +Author: R. Reimann, C. Claessens, T. E. Weiss +Date: Nov. 17, 2020 +Updated: December 2024 + +The statistical method and formulas are described in +CDR (CRES design report, Section 1.3) https://www.overleaf.com/project/5b9314afc673d862fa923d53. +''' +import numpy as np +import configparser + +from mermithid.misc.Constants_numericalunits import * +from mermithid.misc.CRESFunctions_numericalunits import * + +try: + from morpho.utilities import morphologging + logger = morphologging.getLogger(__name__) +except: + print("Run without morpho!") + +class NameSpace(object): + def __init__(self, iteritems): + if type(iteritems) == dict: + iteritems = iteritems.items() + for k, v in iteritems: + setattr(self, k.lower(), v) + def __getattribute__(self, item): + return object.__getattribute__(self, item.lower()) + +############################################################################### +class Sensitivity(object): + """ + Documentation: + * Phase IV sensitivity document: https://www.overleaf.com/project/5de3e02edd267500011b8cc4 + * Talias sensitivity script: https://3.basecamp.com/3700981/buckets/3107037/documents/2388170839 + * Nicks CRLB for frequency resolution: https://3.basecamp.com/3700981/buckets/3107037/uploads/2009854398 + * Molecular contamination in atomic tritium: https://3.basecamp.com/3700981/buckets/3107037/documents/3151077016 + """ + def __init__(self, config_path, verbose=True): + self.cfg = configparser.ConfigParser() + with open(config_path, 'r') as configfile: + self.cfg.read_file(configfile) + + # display configuration + try: + if verbose: + logger.info("Config file content:") + for sect in self.cfg.sections(): + logger.info(' Section: {}'.format(sect)) + for k,v in self.cfg.items(sect): + logger.info(' {} = {}'.format(k,v)) + except: + pass + + self.Experiment = NameSpace({opt: eval(self.cfg.get('Experiment', opt)) for opt in self.cfg.options('Experiment')}) + + # seetings fro molecular or atomic tritium + self.tau_tritium = tritium_livetime + if self.Experiment.atomic: + self.T_mass = tritium_mass_atomic + self.Te_crosssection = tritium_electron_crosssection_atomic + self.T_endpoint = tritium_endpoint_atomic + self.last_1ev_fraction = last_1ev_fraction_atomic + else: + self.T_mass = tritium_mass_molecular + self.Te_crosssection = tritium_electron_crosssection_molecular + self.T_endpoint = tritium_endpoint_molecular + self.last_1ev_fraction = last_1ev_fraction_molecular + + # effective volume if configured + if hasattr(self.Experiment, "v_eff"): + self.effective_volume = self.Experiment.v_eff # v_eff can be configured in Experiment section + + + self.DopplerBroadening = NameSpace({opt: eval(self.cfg.get('DopplerBroadening', opt)) for opt in self.cfg.options('DopplerBroadening')}) + self.FrequencyExtraction = NameSpace({opt: eval(self.cfg.get('FrequencyExtraction', opt)) for opt in self.cfg.options('FrequencyExtraction')}) + self.MagneticField = NameSpace({opt: eval(self.cfg.get('MagneticField', opt)) for opt in self.cfg.options('MagneticField')}) + self.MissingTracks = NameSpace({opt: eval(self.cfg.get('MissingTracks', opt)) for opt in self.cfg.options('MissingTracks')}) + self.PlasmaEffects = NameSpace({opt: eval(self.cfg.get('PlasmaEffects', opt)) for opt in self.cfg.options('PlasmaEffects')}) + + if self.cfg.has_section('FinalStates'): + self.FinalStates = NameSpace({opt: eval(self.cfg.get('FinalStates', opt)) for opt in self.cfg.options('FinalStates')}) + + if not self.Experiment.atomic and not self.cfg.has_section('FinalStates'): + logger.warning(f"No configuration of ground state width uncertainty. Using default value {ground_state_width_uncertainty/eV} eV") + + + # SENSITIVITY + def SignalRate(self): + """signal events in the energy interval before the endpoint, scale with DeltaE**3""" + self.EffectiveVolume() + signal_rate = self.Experiment.number_density*self.effective_volume*self.last_1ev_fraction/self.tau_tritium + if not self.Experiment.atomic: + if hasattr(self.Experiment, 'gas_fractions'): + avg_n_T_atoms = self.AvgNumTAtomsPerParticle_MolecularExperiment(self.Experiment.gas_fractions, self.Experiment.H2_type_gas_fractions) + signal_rate *= avg_n_T_atoms + else: + signal_rate *= 2 + if hasattr(self.Experiment, 'active_gas_fraction'): + signal_rate *= self.Experiment.active_gas_fraction + return signal_rate + + def BackgroundRate(self): + """background rate, can be calculated from multiple components. + Currently, RF noise and cosmic ray backgrounds are included. + Assumes that background rate is constant over considered energy / frequency range.""" + self.cosmic_ray_background = self.Experiment.cosmic_ray_bkgd_per_tritium_particle*self.Experiment.number_density*self.effective_volume + self.background_rate = self.RF_background_rate_per_eV + self.cosmic_ray_background + return self.background_rate + + def SignalEvents(self): + """Number of signal events.""" + return self.SignalRate()*self.Experiment.LiveTime*self.DeltaEWidth()**3 + + def BackgroundEvents(self): + """Number of background events.""" + return self.BackgroundRate()*self.Experiment.LiveTime*self.DeltaEWidth() + + def DeltaEWidth(self): + """optimal energy bin width""" + labels, sigmas, deltas = self.get_systematics() + return np.sqrt(self.BackgroundRate()/self.SignalRate() + + 8*np.log(2)*(np.sum(sigmas**2))) + + def StatSens(self): + """Pure statistic sensitivity assuming Poisson count experiment in a single bin + As defined, it needs to be squared before being added to the systematic component""" + sig_rate = self.SignalRate() + DeltaE = self.DeltaEWidth() + sens = 2/(3*sig_rate*self.Experiment.LiveTime)*np.sqrt(sig_rate*self.Experiment.LiveTime*DeltaE + +self.BackgroundRate()*self.Experiment.LiveTime/DeltaE) + return sens + + def SystSens(self): + """Pure systematic component to sensitivity + As defined, it needs to be squared before being added to the statistical component""" + labels, sigmas, deltas = self.get_systematics() + sens = 4*np.sqrt(np.sum((sigmas*deltas)**2)) + return sens + + def sensitivity(self, **kwargs): + """Combined statisical and systematic uncertainty. + Using kwargs settings in namespaces can be changed. + Example how to change number density which lives in namespace Experiment: + self.sensitivity(Experiment={"number_density": rho}) + """ + for sect, options in kwargs.items(): + for opt, val in options.items(): + self.__dict__[sect].__dict__[opt] = val + + StatSens = self.StatSens() + SystSens = self.SystSens() + + # Standard deviation on a measurement of m_beta**2 assuming a mass of zero + sigma_m_beta_2 = np.sqrt(StatSens**2 + SystSens**2) + return sigma_m_beta_2 + + def CL90(self, **kwargs): + """ Gives 90% CL upper limit on neutrino mass.""" + # 90% of gaussian are contained in +-1.64 sigma region + #return np.sqrt(np.sqrt(1.64)*self.sensitivity(**kwargs)) + return np.sqrt(1.64*self.sensitivity(**kwargs)) + + def sterial_m2_limit(self, Ue4_sq): + return np.sqrt(1.64*np.sqrt((self.StatSens()/Ue4_sq)**2 + self.SystSens()**2)) + + # PHYSICS Functions + + def AvgNumTAtomsPerParticle_MolecularExperiment(self, gas_fractions, H2_type_gas_fractions): + """ + Given gas composition info (H2 vs. other gases, and how much of each H2-type isotopolog), returns an average number of tritium atoms per gas particle. + + Inputs: + - gas_fractions: dict of composition fractions of each gas (different from scatter fractions!); all H2 isotopologs are combined under key 'H2' + - H2_type_gas_fractions: dict with fraction of each isotopolog, out of total amount of H2 + """ + H2_iso_avg_num = 0 + for (key, val) in H2_type_gas_fractions.items(): + if key=='T2': + H2_iso_avg_num += 2*val + elif key=='HT' or key=='DT': + H2_iso_avg_num += val + elif key=='H2' or key=='HD' or key=='D2': + pass + return gas_fractions['H2']*H2_iso_avg_num + + def BToKeErr(self, BErr, B): + return e*BErr/(2*np.pi*frequency(self.T_endpoint, B)/c0**2) + + def KeToBerr(self, KeErr, B): + return KeErr/e*(2*np.pi*frequency(self.T_endpoint, B)/c0**2) + + def track_length(self, rho): + return track_length(rho, self.T_endpoint, not self.Experiment.atomic) + + # SYSTEMATICS + + def get_systematics(self): + """ Returns list of energy broadenings (sigmas) and + uncertainties on these energy broadenings (deltas) + for all considered systematics. We need to make sure + that we do not include effects twice or miss any + important effect. + + Returns: + * list of labels + * list of energy broadenings + * list of energy broadening uncertainties + """ + + # Different types of uncertainty contributions + sigma_trans, delta_sigma_trans = self.syst_doppler_broadening() + sigma_f, delta_sigma_f = self.syst_frequency_extraction() + sigma_B, delta_sigma_B = self.syst_magnetic_field() + sigma_Miss, delta_sigma_Miss = self.syst_missing_tracks() + sigma_Plasma, delta_sigma_Plasma = self.syst_plasma_effects() + + labels = ["Thermal Doppler Broadening", "Noise (f_carrier and f_lsb)", "Magnetic Field", "Missing Tracks", "Plasma Effects"] + sigmas = [sigma_trans, sigma_f, sigma_B, sigma_Miss, sigma_Plasma] + deltas = [delta_sigma_trans, delta_sigma_f, delta_sigma_B, delta_sigma_Miss, delta_sigma_Plasma] + + if not self.Experiment.atomic: + labels.append("Molecular final state") + sigmas.append(ground_state_width) + if self.cfg.has_section('FinalStates') and hasattr(self.FinalStates, "ground_state_width_uncertainty_fraction"): + deltas.append(self.FinalStates.ground_state_width_uncertainty_fraction*ground_state_width) + else: + deltas.append(ground_state_width_uncertainty) + + return np.array(labels), np.array(sigmas), np.array(deltas) + + def print_statistics(self): + print("Contribution to sigma_(m_beta^2)", " "*18, "%.2f"%(self.StatSens()/meV**2), "meV^2 ->", "%.2f"%(np.sqrt(self.StatSens())/meV), "meV") + print("Statistical mass limit", " "*18, "%.2f"%(np.sqrt(1.64*self.StatSens())/meV), "meV") + + def print_systematics(self): + labels, sigmas, deltas = self.get_systematics() + + print() + sigma_squared = 0 + for label, sigma, delta in zip(labels, sigmas, deltas): + print(label, " "*(np.max([len(l) for l in labels])-len(label)), "%8.2f"%(sigma/meV), "+/-", "%8.2f"%(delta/meV), "meV") + sigma_squared += sigma**2 + sigma_total = np.sqrt(sigma_squared) + print("Total sigma", " "*(np.max([len(l) for l in labels])-len("Total sigma")), "%8.2f"%(sigma_total/meV),) + #try: + # print("(Contribution from axial variation: ", "%8.2f"%(self.sigma_K_reconstruction/meV)," meV)") + #except AttributeError: + # pass + print("Contribution to sigma_(m_beta^2)", " "*18, "%.2f"%(self.SystSens()/meV**2), "meV^2 ->", "%.2f"%(np.sqrt(self.SystSens())/meV), "meV") + print("Systematic mass limit", " "*18, "%.2f"%(np.sqrt(1.64*self.SystSens())/meV), "meV") + logger.info("Max allowed cyclotron frequency uncertainty: {} Hz".format(np.sqrt(self.max_allowed_f_var)/Hz)) + logger.info("Corresponding reconstruction efficiency: {}".format(self.recon_efficiency)) + return np.sqrt(1.64*self.SystSens())/meV, np.sqrt(np.sum(sigmas**2))/meV + + def syst_doppler_broadening(self): + # estimated standard deviation of Doppler broadening distribution from + # translational motion of tritium atoms / molecules + # Predicted uncertainty on standard deviation, based on expected precision + # of temperature knowledge + if self.DopplerBroadening.UseFixedValue: + sigma = self.DopplerBroadening.Default_Systematic_Smearing + delta = self.DopplerBroadening.Default_Systematic_Uncertainty + return sigma, delta + + # termal doppler broardening + gasTemp = self.DopplerBroadening.gas_temperature + mass_T = self.T_mass + endpoint = self.T_endpoint + + # these factors are mainly negligible in the recoil equation below + E_rec = 3.409 * eV # maximal value # same for molecular tritium? + mbeta = 0*eV # term neglidible + betanu = 1 # neutrinos are fast + # electron-neutrino correlation term: 1 + 0.105(6)*betae*cosThetanu + # => expectation value of cosThetanu = 0.014 + cosThetaenu = 0.014 + + Ke = endpoint + betae = np.sqrt(Ke**2+2*Ke*me*c0**2)/(Ke+me*c0**2) ## electron speed at energy Ke + Emax = endpoint + me*c0**2 + Ee = endpoint + me*c0**2 + p_rec = np.sqrt( Emax**2-me**2*c0**4 + (Emax - Ee - E_rec)**2 - mbeta**2 + 2*Ee*(Emax - Ee - E_rec)*betae*betanu*cosThetaenu ) + sigma_trans = np.sqrt(p_rec**2/(2*mass_T)*2*kB*gasTemp) + + if self.Experiment.atomic == True: + delta_trans = np.sqrt(p_rec**2/(2*mass_T)*kB/gasTemp*self.DopplerBroadening.gas_temperature_uncertainty**2) + else: + delta_trans = sigma_trans*self.DopplerBroadening.fraction_uncertainty_on_doppler_broadening + return sigma_trans, delta_trans + + """ + def syst_frequency_extraction(self): + # cite{https://3.basecamp.com/3700981/buckets/3107037/uploads/2009854398} (Section 1.2, p 7-9) + # Are we double counting the antenna collection efficiency? We use it here. Does it also impact the effective volume, v_eff ? + # Are we double counting the effect of magnetic field uncertainty here? + # Is 'sigma_f_Bfield' the same as 'rRecoErr', 'delta_rRecoErr', 'rRecoPhiErr', 'delta_rRecoPhiErr'? + + if self.FrequencyExtraction.UseFixedValue: + sigma = self.FrequencyExtraction.Default_Systematic_Smearing + delta = self.FrequencyExtraction.Default_Systematic_Uncertainty + return sigma, delta + + # Cramer-Rao lower bound / how much worse are we than the lower bound + ScalingFactorCRLB = self.FrequencyExtraction.CRLB_scaling_factor + ts = self.FrequencyExtraction.track_timestep + # "This is apparent in the case of resonant patch antennas and cavities, in which the time scale of the signal onset is set by the Q-factor of the resonant structure." + # You can get it from the finite impulse response of the antennas from HFSS + Gdot = self.FrequencyExtraction.track_onset_rate + + fEndpoint = frequency(self.T_endpoint, self.MagneticField.nominal_field) + betae = beta(self.T_endpoint) + Pe = rad_power(self.T_endpoint, self.FrequencyExtraction.pitch_angle, self.MagneticField.nominal_field) + alpha_approx = fEndpoint * 2 * np.pi * Pe/me/c0**2 # track slope + # quantum limited noise + sigNoise = np.sqrt((2*np.pi*fEndpoint*hbar*self.FrequencyExtraction.amplifier_noise_scaling+kB*self.FrequencyExtraction.antenna_noise_temperature)/ts) # noise level + Amplitude = np.sqrt(self.FrequencyExtraction.epsilon_collection*Pe) + Nsteps = 1 / (self.Experiment.number_density * self.Te_crosssection*betae*c0*ts) # Number of timesteps of length ts + + # sigma_f from Cramer-Rao lower bound in Hz + sigma_f_CRLB = (ScalingFactorCRLB /(2*np.pi) * sigNoise/Amplitude * np.sqrt(alpha_approx**2/(2*Gdot) + + 96.*Nsteps/(ts**2*(Nsteps**4-5*Nsteps**2+4)))) + # uncertainty in alpha + delta_alpha = 6*sigNoise/(Amplitude*ts**2) * np.sqrt(10/(Nsteps*(Nsteps**4-5*Nsteps**2+4))) + # uncetainty in sigma_f in Hz due to uncertainty in alpha + delta_sigma_f_CRLB = delta_alpha * alpha_approx *sigNoise**2/(8*np.pi**2*Amplitude**2*Gdot*sigma_f_CRLB*ScalingFactorCRLB**2) + + # sigma_f from Cramer-Rao lower bound in eV + sigma_K_f_CRLB = e*self.MagneticField.nominal_field/(2*np.pi*fEndpoint**2)*sigma_f_CRLB*c0**2 + delta_sigma_K_f_CRLB = e*self.MagneticField.nominal_field/(2*np.pi*fEndpoint**2)*delta_sigma_f_CRLB*c0**2 + + # combined sigma_f in eV + sigma_f = np.sqrt(sigma_K_f_CRLB**2 + self.FrequencyExtraction.magnetic_field_smearing**2) + delta_sigma_f = np.sqrt((delta_sigma_K_f_CRLB**2 + self.FrequencyExtraction.magnetic_field_uncertainty**2)/2) + + # the magnetic_field_smearing and uncertainty added here consider the following effect: + # thinking in terms of a phase II track, there is some smearing of the track / width of the track which influences the frequency extraction + # this does not account for any effect comming from converting frequency to energy + # the reason behind the track width / smearing is the change in B field that the electron sees within one axial oscillation. + # Depending on the trap shape this smearing may be different. + + return sigma_f, delta_sigma_f + """ + + def syst_magnetic_field(self): + + # magnetic field uncertainties can be decomposed in several part + # * true magnetic field inhomogeneity + # (would be there also without a trap) + # * magnetic field calibration has uncertainties + # (would be there also without a trap) + # * position / pitch angle reconstruction has uncertainties + # (this can even be the degenerancy we see for harmonic traps) + # (depends on trap shape) + + if self.MagneticField.UseFixedValue: + sigma = self.MagneticField.Default_Systematic_Smearing + delta = self.MagneticField.Default_Systematic_Uncertainty + return sigma, delta + + B = self.MagneticField.nominal_field + if self.MagneticField.useinhomogenarity: + inhomogenarity = self.MagneticField.inhomogenarity + sigma = self.BToKeErr(inhomogenarity*B, B) + return sigma, 0.05*sigma + + BMapErr = self.MagneticField.probe_repeatability # Probe Repeatability + delta_BMapErr = self.MagneticField.probe_resolution # Probe resolution + + BFlatErr = self.MagneticField.BFlatErr # averaging over the flat part of the field + delta_BFlatErr = self.MagneticField.relative_uncertainty_BFlatErr*BFlatErr # UPDATE ? + + Delta_t_since_calib = self.MagneticField.time_since_calibration + shiftBdot = self.MagneticField.shift_Bdot + smearBdot = self.MagneticField.smear_Bdot + delta_shiftBdot = self.MagneticField.uncertainty_shift_Bdot + delta_smearBdot = self.MagneticField.uncertainty_smearBdot + BdotErr = Delta_t_since_calib * np.sqrt(shiftBdot**2 + smearBdot**2) + delta_BdotErr = Delta_t_since_calib**2/BdotErr * np.sqrt(shiftBdot**2 * delta_shiftBdot**2 + smearBdot**2 * delta_smearBdot**2) + + # position uncertainty is linear in wavelength + # position uncertainty is nearly constant w.r.t. radial position + # based on https://3.basecamp.com/3700981/buckets/3107037/uploads/3442593126 + rRecoErr = self.MagneticField.rRecoErr + delta_rRecoErr = self.MagneticField.relative_Uncertainty_rRecoErr * rRecoErr + + rRecoPhiErr = self.MagneticField.rRecoPhiErr + delta_rRecoPhiErr = self.MagneticField.relative_uncertainty_rRecoPhiErr * rRecoPhiErr + + rProbeErr = self.MagneticField.rProbeErr + delta_rProbeErr = self.MagneticField.relative_uncertainty_rProbeErr * rProbeErr + + rProbePhiErr = self.MagneticField.rProbePhiErr + delta_rProbePhiErr = self.MagneticField.relative_uncertainty_rProbePhiErr * rProbePhiErr + + Berr = np.sqrt(BMapErr**2 + + BFlatErr**2 + + BdotErr**2 + + rRecoErr**2 + + rRecoPhiErr**2 + + rProbeErr**2 + + rProbePhiErr**2) + + delta_Berr = 1/Berr * np.sqrt(BMapErr**2 * delta_BMapErr**2 + + BFlatErr**2 * delta_BFlatErr**2 + + BdotErr**2 * delta_BdotErr**2 + + rRecoErr**2 * delta_rRecoErr**2 + + rRecoPhiErr**2 * delta_rRecoPhiErr**2 + + rProbeErr**2 * delta_rProbeErr**2 + + rProbePhiErr**2 * delta_rProbePhiErr**2) + + return self.BToKeErr(Berr, B), self.BToKeErr(delta_Berr, B) + + def syst_missing_tracks(self): + # this systematic should describe the energy broadening due to the line shape. + # Line shape is caused because you miss the first n tracks but then detect the n+1 + # track and you assume that this is the start frequency. + # This depends on the gas composition, density and cross-section. + if self.MissingTracks.UseFixedValue: + sigma = self.MissingTracks.Default_Systematic_Smearing + delta = self.MissingTracks.Default_Systematic_Uncertainty + return sigma, delta + else: + raise NotImplementedError("Missing track systematic is not implemented.") + + def syst_plasma_effects(self): + if self.PlasmaEffects.UseFixedValue: + sigma = self.PlasmaEffects.Default_Systematic_Smearing + delta = self.PlasmaEffects.Default_Systematic_Uncertainty + return sigma, delta + else: + raise NotImplementedError("Plasma effect sysstematic is not implemented.") diff --git a/mermithid/sensitivity/__init__.py b/mermithid/sensitivity/__init__.py new file mode 100644 index 00000000..f5cebfc9 --- /dev/null +++ b/mermithid/sensitivity/__init__.py @@ -0,0 +1,7 @@ +''' +''' + +from __future__ import absolute_import + +from . import SensitivityFormulas +from . import SensitivityCavityFormulas \ No newline at end of file diff --git a/test_analysis/CCA_Sensitivity.py b/test_analysis/CCA_Sensitivity.py new file mode 100644 index 00000000..670c2470 --- /dev/null +++ b/test_analysis/CCA_Sensitivity.py @@ -0,0 +1,96 @@ +""" +Script to make Sensitivty plots for cavity experiments +Author: C. Claessens, T. Weiss +Date: October 6, 2023 +""" + + +from morpho.utilities import morphologging, parser +logger = morphologging.getLogger(__name__) + +from mermithid.processors.Sensitivity import CavitySensitivityCurveProcessor + +import numpy as np + +# Configuration for CCA Sensitivity vs. density plot +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_CCA_Experiment.cfg", + #"config_file_path": "/host_repos/sensitivity_branches/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + "plot_path": "./cca_sensitivity_vs_density_curve.pdf", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": True, + "atomic_axis": False, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 20], + "density_range": [1e13,1e21], + "efficiency_range": [0.0001, 1], + #"density_range": [1e8, 1e12], + "main_curve_upper_label": r"CCA", #r"Molecular"+"\n"+"Reaching target", + "comparison_curve_label": [#"Molecular, reaching target", + "Atomic, alternative scenario 1"], + #"Atomic, reaching target"], + # #["Molecular"+"\n"+"Conservative", "Atomic"+"\n"+"Conservative", r"Atomic"+"\n"+"Reaching target"], + "goals": {"Pilot T goal (0.1 eV)": 0.1, "Phase IV (0.04 eV)": 0.04}, + "comparison_curve": False, + "comparison_curve_colors": ["red"], + "comparison_config_file_path": [#"/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + "/termite/sensitivity_config_files/Config_atomic_find_factor_22_Experiment_conservative.cfg"], #, + #"/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + "comparison_label_y_position": [2, 0.105, 0.046], #[2, 0.105, 0.046], + "comparison_label_x_position": [4.5e15, 7e14, 7e14], #[4.5e15, 2.2e16, 1e15], + #"sigmae_theta_r": 0.159, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 4e14, #4e14, #0.02, #1e14, + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": True + } +sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +sens_curve.Configure(sens_config_dict) +sens_curve.Run() + + +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_CCA_Experiment.cfg", #Config_PhaseIII_325MHz_Experiment_conservative.cfg", + "plot_path": "./sensitivity_vs_exposure_curve_CCA.pdf", + # optional + "figsize": (10,6), + "fontsize": 15, + "legend_location": "upper right", + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "density_axis": False, + "exposure_axis": True, + "cavity": True, + "add_PhaseII": False, + #"add_1year_1cav_point_to_last_ref": True, + "y_limits": [10e-3, 500], + "density_range": [1e12, 1e19], + "exposure_range": [1e-11, 1e10], + "main_curve_upper_label": r"CCA", #"Molecular, conservative", + #"main_curve_lower_label": r"$\sigma^\bar{B}_\mathrm{reco} = 0.07\,\mathrm{eV}$", + "comparison_curve": False, + "upper_label_y_position": 0.7, + "label_x_position": 0.015, #1.5e15, #0.02, #1e14, + "goals_x_position": 0.2e-10, #2e12 #0.0002 + "goals_y_rel_position": 0.4 + } +sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +sens_curve.Configure(sens_config_dict) +sens_curve.Run() + + + + + + + + + + diff --git a/test_analysis/Cavity_Sensitivity_analysis.py b/test_analysis/Cavity_Sensitivity_analysis.py new file mode 100644 index 00000000..f003c160 --- /dev/null +++ b/test_analysis/Cavity_Sensitivity_analysis.py @@ -0,0 +1,474 @@ +""" +Script to make sensitivity plots for cavity experiments +Author: C. Claessens, T. E. Weiss +Date: October 6, 2023 +""" + + +from morpho.utilities import morphologging, parser +logger = morphologging.getLogger(__name__) + +from mermithid.processors.Sensitivity import CavitySensitivityCurveProcessor + +import numpy as np + +# Configuration for Sensitivity vs. exposure plot +# Phase III and IV at 325 MHz +# Phase II point and curve for comparison + +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", #Config_PhaseIII_325MHz_Experiment_conservative.cfg", + "plot_path": "./sensitivity_vs_exposure_curve_PhaseIV_87deg_min_pitch.pdf", + # optional + "figsize": (10,6), + "fontsize": 15, + "legend_location": "upper right", + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "density_axis": False, + "exposure_axis": True, + "cavity": True, + "add_PhaseII": True, + "add_1year_1cav_point_to_last_ref": True, + "PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg", + "y_limits": [10e-3, 500], + "density_range": [0, 1e19], + "exposure_range": [1e-11, 1e4], + "main_curve_upper_label": r"Phase IV scenario", #"Molecular, conservative", + #"main_curve_lower_label": r"$\sigma^\bar{B}_\mathrm{reco} = 0.07\,\mathrm{eV}$", + "goals": {"Phase IV (0.04 eV)": 0.04}, + "comparison_curve": False, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment_conservative.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg" + ], + "comparison_curve_label": [r"Molecular, reaching PIII target", "Atomic, conservative", "Atomic, reaching PIV target"], + "comparison_curve_colors": ["blue", "darkred", "red"], + #"B_inhomogeneity": np.linspace(0.1, 2.1, 10)*1e-6, + #"B_inhom_uncertainty": 0.01, + #"sigmae_theta_r": 0.159, #in eV, energy broadening from theta and r reconstruction + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 0.015, #1.5e15, #0.02, #1e14, + "goals_x_position": 0.2e-10, #2e12 #0.0002 + "goals_y_rel_position": 0.4 + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + +# Configuration for Sensitivity vs. frequency plot +sens_config_dict2 = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "plot_path": "./sensitivity_vs_frequency_Oct-22-2024.pdf", + # optional + "figsize": (9,6), + "fontsize": 15, + "legend_location": "upper left", + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "density_axis": False, + "frequency_axis": True, + "magnetic_field_axis": True, + "cavity": True, + "add_PhaseII": False, + "add_1year_1cav_point_to_last_ref": False, + "PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg", + "y_limits": [10e-3, 10], + "frequency_range": [1e7, 20e9], + #"efficiency_range": [0.0001, 1], + "density_range": [1e7, 1e20], + "main_curve_upper_label": r"Molecular, conservative", + "main_curve_lower_label": r"$\sigma^\bar{B}_\mathrm{reco} = 0.07\,\mathrm{eV}$", + "goals": {"Phase IV (0.04 eV)": (0.04**2/1.64)}, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment_conservative.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + "comparison_curve_label": [r"Molecular, reaching PIII target", "Atomic, conservative", "Atomic, reaching PIV target"], + "comparison_curve_colors": ["blue", "darkred", "red"], + #"config_file_path": "/host_repos/sensitivity_branches/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg", + #"B_inhomogeneity": np.linspace(0.1, 2.1, 10)*1e-6, + #"B_inhom_uncertainty": 0.01, + "sigmae_theta_r": 0.159, #in eV, energy broadening from theta and r reconstruction + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 1e8, #1.5e15, #0.02, #1e14, + "goals_x_position": 1e9, #2e12 #0.0002 + "goals_y_rel_position": 0.4, + "plot_key_parameters": False + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict2) +#sens_curve.Run() + +# Configuration for Sensitivity vs. density plot for different B sigmas +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment_conservative.cfg", + "plot_path": "./sensitivity_vs_density_curve_with_sigmaBcorrs.pdf", + # optional + "track_length_axis": True, + "molecular_axis": True, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,1e19], + "efficiency_range": [0.0001, 1], + #"density_range": [1e8, 1e12], + "main_curve_upper_label": r"Molecular"+"\n"+"2 years"+"\n"+r"$\sigma^B_\mathrm{corr} = 1\,\mathrm{eV}$", + "main_curve_lower_label": r"$\sigma^B_\mathrm{corr} = 0.16\,\mathrm{eV}$", + "comparison_curve_label": [r"Atomic"+"\n"+r"10 $\times$ 3 years"+"\n"+r"$\sigma^B_\mathrm{corr} = 0.16\,\mathrm{eV}$"], + "goals": {"Phase III (0.2 eV)": 0.2, "Phase IV (0.04 eV)": 0.04}, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment_conservative.cfg"], + #"comparison_config_file_path": "/host_repos/sensitivity_branches/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg", + #"B_inhomogeneity": np.linspace(0.1, 2.1, 10)*1e-6, + #"B_inhom_uncertainty": 0.01, + "configure_sigma_theta_r": True, + "sigmae_theta_r": np.linspace(0.16, 1., 10), #in eV, energy broadening from theta and r reconstruction + "comparison_label_y_position": [0.044], + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 1.5e15, #0.02, #1e14, + "goals_x_position": 2e12, #0.0002 + "plot_key_parameters": False + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + +# Configuration for Sensitivity vs. density plot +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam_threshold.cfg", #"/termite/sensitivity_config_files/Config_LFA_Experiment_1GHz.cfg", #Config_atomic_325MHz_Experiment_conservative.cfg", + "plot_path": "./LFA_and_PhaseIV_sensitivity_vs_density_target_and_threshold_Mar-24-2026.pdf", + # optional + "figsize": (7.5,6.4), + "fontsize": 15, + "track_length_axis": False, + "legend_location": "upper left", + "legend_bbox_to_anchor": (0.155,-0.01,1.12,0.885), #(0.17,-0.01,1.1,0.87), + "molecular_axis": False, + "atomic_axis": True, + "density_axis": True, + "track_length_axis": True, + "effs_for_sampled_radii": True, + "y_limits": [2e-2, 6.5], + "density_range": [3e14,3e18], #5e13 + "det_thresh_range": [5, 115], + "add_point_at_configured_density": True, + "main_curve_upper_label": r"LFA, threshold scenario", #560 MHz #Phase III scenario: 1 GHz", + "goals": {"LFA threshold (0.7 eV)": 0.7, "Phase IV (0.04 eV)": 0.04}, #"Pilot T goal (0.1 eV)": 0.1, + "goals_x_position": {"LFA threshold (0.7 eV)": 2.6e16, "Phase IV (0.04 eV)": 4e14}, #6e14, #3.3e14, #5.5e13, + "goals_y_rel_position": {"LFA threshold (0.7 eV)": 1.1, "Phase IV (0.04 eV)": 0.79}, #0.755 + "comparison_curve": False, + "verbose": False, + "main_curve_color": "blue", + "comparison_curve_colors": ["blue", "darkred", "black"], + "main_curve_linestyle": "dashed", + "comparison_curve_linestyles": ["solid", "dotted", "dashdot"], + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam.cfg", "/termite/sensitivity_config_files/Config_PIVmodule1_150MHz_minpitch_87deg.cfg", "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg"], + "comparison_curve_label": [r"LFA, target scenario", r"One full-size module", r"Phase IV: Ten full-size modules"], #: 150 MHz + "comparison_label_y_position": [2, 0.105, 0.046], #[2, 0.105, 0.046], + "comparison_label_x_position": [4.5e15, 7e14, 7e14], #[4.5e15, 2.2e16, 1e15], + #"sigmae_theta_r": 0.159, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 4e14, #4e14, #0.02, #1e14, + "plot_key_parameters": False, + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam_threshold.cfg", #"/termite/sensitivity_config_files/Config_LFA_Experiment.cfg", + "plot_path": "./LFA_and_PhaseIV_sensitivity_vs_livetime_curve_target_and_threshold_Apr-14-2026.pdf", #ncav-eff-time + "exposure_axis": True, + # optional + "figsize": (8.3, 6.3), #(10,6), + "fontsize": 15, + "legend_location": "upper right", + "legend_bbox_to_anchor": (-0.,0,0.86,0.971), #last entry: 0.955 + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "exposure_axis": False, + "density_axis": False, + "ncav_eff_time_axis": False, + "ncavities_livetime_axis": False, + "livetime_axis": True, + "cavity": True, + "add_PhaseII": False, + "PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg", + "add_1year_1cav_point_to_last_ref": False, + "y_limits": [2e-2, 3], + #"density_range": [1e12,1e19], + "year_range": [0.1,35], + "main_curve_upper_label": r"LFA, threshold: $1.7\,$m$^3$, 1 yr", # 560 MHz, $V = 1.7\,$m$^3$ + "goals": {"LFA threshold (0.7 eV)": 0.7, "Phase IV (0.04 eV)": 0.04}, + "goals_x_position": {"LFA threshold (0.7 eV)": 4.5, "Phase IV (0.04 eV)": 0.108}, #6e14, #3.3e14, #5.5e13, + "goals_y_rel_position": {"LFA threshold (0.7 eV)": 0.83, "Phase IV (0.04 eV)": 0.83}, #6e14, #3.3e14, #5.5e13, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam.cfg", "/termite/sensitivity_config_files/Config_PIVmodule1_150MHz_minpitch_87deg.cfg", "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg"], #"/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam_threshold.cfg", + "comparison_curve_label": [r"LFA, target: $1.7\,$m$^3$, 1 yr", r'One full-size module: $100\,$m$^3$, 1 yr', r"Phase IV$-$Ten full-size modules: $1000\,$m$^3$, 8 yrs"], #150 MHz, $V = 94\,$m$^3$ + "main_curve_color": "blue", + "comparison_curve_colors": ["blue", "darkred", "black"], + "main_curve_linestyle": "dashed", + "comparison_curve_linestyles": ["solid", "dotted", "dashdot"], + "main_curve_marker": "d", + "comparison_curve_markers": ["o", "^", "X"], + "optimize_main_density": False, + "optimize_comparison_density": False, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 0.115, + #"goals_x_position": 0.12, #4e-2, #<-- Number for ncav*eff*time #0.11, <-- Number for ncavities*livetime + #"goals_y_rel_position": 0.86, #0.84, <-- Number for ncav*eff*time #0.81, <-- Number for ncavities*livetime + } +sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +sens_curve.Configure(sens_config_dict) +sens_curve.Run() + + +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", #"/termite/sensitivity_config_files/Config_LFA_Experiment.cfg", + "plot_path": "./PhaseIV_scenario_sensitivity_vs_livetime_curve_March-4-2025.pdf", #ncav-eff-time + "exposure_axis": True, + # optional + "figsize": (8.3, 6.3), #(10,6), + "fontsize": 15, + "legend_location": "upper right", + "legend_bbox_to_anchor": (-0.,0,0.86,0.955), + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "exposure_axis": False, + "density_axis": False, + "ncav_eff_time_axis": False, + "ncavities_livetime_axis": False, + "livetime_axis": True, + "cavity": True, + "add_PhaseII": False, + "PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg", + "add_1year_1cav_point_to_last_ref": False, + "y_limits": [1e-2, 1.5e-1], + #"density_range": [1e12,1e19], + "year_range": [0.1,5e2], + "main_curve_upper_label": r"Default Phase IV scenario", + "goals": {"Phase IV (0.04 eV)": 0.04}, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/P4_scenario_USR_sigmaeR.cfg", "/termite/sensitivity_config_files/P4_scenario_USR_nonexposure.cfg"], + "comparison_curve_label": [r"Less field broadening, fewer events", r"Cavity magic"], #, "One Phase IV cavity"], + "main_curve_color": "red", + "comparison_curve_colors": ["darkorange", "purple"], + "optimize_main_density": False, + "optimize_comparison_density": True, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 0.115, + "goals_x_position": 1.1e-1, #<-- Number for ncav*eff*time #0.11, <-- Number for ncavities*livetime + "goals_y_rel_position": 0.89, #0.84, <-- Number for ncav*eff*time #0.81, <-- Number for ncavities*livetime + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + + +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_LFA_Experiment_1GHz.cfg", #"/termite/sensitivity_config_files/Config_LFA_Experiment.cfg", + "plot_path": "./lfa_and_PhaseIV_sensitivity_vs_exposure_curve_Oct-22-2024.pdf", + "exposure_axis": True, + # optional + "figsize": (10,6), + "fontsize": 15, + "legend_location": "upper right", + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "exposure_axis": True, + "density_axis": False, + "cavity": True, + "add_PhaseII": True, + "PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg", + "add_1year_1cav_point_to_last_ref": False, + "y_limits": [10e-3, 500], + #"density_range": [1e12,1e19], + "exposure_range": [1e-11, 1e4], + #"main_curve_upper_label": r"LFA (atomic T)$\,-\,$1 GHz", + "main_curve_upper_label": r"Phase III scenario: 1 GHz", + #"main_curve_lower_label": r"$\sigma^\bar{B}_\mathrm{reco} = 0.16\,\mathrm{eV}$", + "goals": {"Phase IV (0.04 eV)": (0.04)}, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam.cfg", + "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg"], + "comparison_curve_label": [r"Phase III scenario: 560 MHz", r"Phase IV scenario: 150 MHz"], + "main_curve_color": "blue", + "comparison_curve_colors": ["red", "black"], + "optimize_main_density": False, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 0.015, #1.5e15, #0.02, #1e14, + "goals_x_position": 0.2e-10, #2e12 #0.0002 + "goals_y_rel_position": 0.4, + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + +# Configuration for Sensitivity vs. density plot for best possible molecular scenario +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment_best_case.cfg", + "plot_path": "./sensitivity_vs_density_T2_best_case_curve.pdf", + # optional + "figsize": (6.7,6), + "track_length_axis": False, + "molecular_axis": True, + "atomic_axis": False, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,3e18], + "efficiency_range": [0.0001, 1], + "main_curve_upper_label": r"Molecular, best-case scenario", + "goals": {"Phase III (0.2 eV)": 0.2, "Phase IV (0.04 eV)": 0.04}, + "comparison_curve": False, + "sigmae_theta_r": 0.159, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 4e14, + "goals_x_position": 1.2e12, + "plot_key_parameters": False + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + +# Configuration for Sensitivity vs. densty plot for molecular best case and Phase III atomic scenarios +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg", + "plot_path": "./sensitivity_vs_density_T2_best_case_curve_comparison.pdf", + # optional + "figsize": (6.7,6), + "legend_location": "upper left", + "track_length_axis": True, + "molecular_axis": True, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [3e-2, 0.5], + "density_range": [1e14, 7e17],#[1e12,3e18], + "efficiency_range": [0.0001, 1], + "main_curve_upper_label": r"Atomic, $\Delta B_{r, \phi, t}=0$, rate 'boosted' $\times 2$", + "goals": {"Phase IV (0.04 eV)": 0.04}, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment_best_case.cfg"], + "comparison_curve_label": [r"Molecular, same conditions"], + "comparison_label_y_position": [2, 0.105, 0.046], + "comparison_label_x_position": [4.5e15, 7e14, 7e14], + "sigmae_theta_r": 0.0, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 4e14, + "goals_x_position": 1.2e14, + "goals_y_rel_position": 1.1, + "plot_key_parameters": False + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + +# Configuration for Sensitivity vs. density plot with 0 positional field uncertainty +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg", + "plot_path": "./sensitivity_vs_density_T2_no-DeltaB-r-phi-t.pdf", + # optional + "figsize": (7.5,5), + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [3e-2, 0.5], #[2e-2, 4], + "density_range": [1e14, 7e17], #[1e12,3e18], + "efficiency_range": [0.0001, 1], + "main_curve_upper_label": r" ", #r"Atomic, $\Delta B_{r, \phi, t}=0$, rate 'boosted' $\times 2$", + "goals": {"Phase IV (0.04 eV)": 0.04}, + "comparison_curve": False, + "sigmae_theta_r": 0.0, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 4e14, + "goals_x_position": 1.2e14, + "goals_y_rel_position": 1.1, + #"legend_location": "upper left", + "plot_key_parameters": True + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + +# Configuration for CCA Sensitivity vs. density plot +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_CCA_pitchAngUncertVerification.cfg", + "plot_path": "./cca_sensitivity_vs_density_curve_new_config.pdf", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": True, + "atomic_axis": False, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 20], + "density_range": [1e13,1e21], + "efficiency_range": [0.0001, 1], + #"density_range": [1e8, 1e12], + "main_curve_upper_label": r"CCA", #r"Molecular"+"\n"+"Reaching target", + "comparison_curve_label": [#"Molecular, reaching target", + "Atomic, alternative scenario 1"], + #"Atomic, reaching target"], + # #["Molecular"+"\n"+"Conservative", "Atomic"+"\n"+"Conservative", r"Atomic"+"\n"+"Reaching target"], + "goals": {"Pilot T goal (0.1 eV)": 0.1, "Phase IV (0.04 eV)": 0.04}, + "comparison_curve": False, + "comparison_curve_colors": ["red"], + "comparison_config_file_path": [#"/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + "/termite/sensitivity_config_files/Config_atomic_find_factor_22_Experiment_conservative.cfg"], #, + #"/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + "comparison_label_y_position": [2, 0.105, 0.046], #[2, 0.105, 0.046], + "comparison_label_x_position": [4.5e15, 7e14, 7e14], #[4.5e15, 2.2e16, 1e15], + #"sigmae_theta_r": 0.159, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 4e14, #4e14, #0.02, #1e14, + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": True + } + +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + + + + + + diff --git a/test_analysis/LFA_Sensitivity.py b/test_analysis/LFA_Sensitivity.py new file mode 100644 index 00000000..1e752c52 --- /dev/null +++ b/test_analysis/LFA_Sensitivity.py @@ -0,0 +1,112 @@ +""" +Script to make Sensitivty plots for cavity experiments +Author: C. Claessens, T. Weiss +Date: October 6, 2023 +""" + + +from morpho.utilities import morphologging, parser +logger = morphologging.getLogger(__name__) + +from mermithid.processors.Sensitivity import CavitySensitivityCurveProcessor + +import numpy as np + +# Configuration for CCA Sensitivity vs. density plot +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/config_LFA_Experiment_1GHz.cfg", #Config_LFA_Experiment_1GHz.cfg", + #"config_file_path": "/host_repos/sensitivity_branches/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + "plot_path": "./lfa_sensitivity_vs_density_curve.pdf", + # optional + "figsize": (7.0,5), #Was (7, 6) + "track_length_axis": True, + "molecular_axis": False, + "atomic_axis": False, + "density_axis": True, + "optimize_main_density": True, + "cavity": True, + #"y_limits": [2e-2, 20], + "y_limits": [2e-1, 20], + "density_range": [1e13,1e20], + #"efficiency_range": [0.0001, 1], + #"density_range": [1e8, 1e12], + "configure_sigma_theta_r": True, + "sigmae_theta_r": np.linspace(0.1, 0.5, 10), #in eV, energy broadening from theta and r reconstruction + "main_curve_upper_label": r"$\sigma^B_\mathrm{corr} = 0.5\,\mathrm{eV}$", #r"Molecular"+"\n"+"Reaching target", + "main_curve_lower_label": r"$\sigma^B_\mathrm{corr} = 0.1\,\mathrm{eV}$", #r"Molecular"+"\n"+"Reaching target", + #"comparison_curve_label": [#"Molecular, reaching target", + # "Atomic, alternative scenario 1"], + # #"Atomic, reaching target"], + # # #["Molecular"+"\n"+"Conservative", "Atomic"+"\n"+"Conservative", r"Atomic"+"\n"+"Reaching target"], + #"goals": {"Pilot T goal (0.1 eV)": 0.1, "Phase IV (0.04 eV)": 0.04}, + "comparison_curve": False, + "comparison_curve_colors": ["red"], + #"comparison_config_file_path": [#"/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment.cfg", + # "/termite/sensitivity_config_files/Config_atomic_find_factor_22_Experiment_conservative.cfg"], #, + # #"/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + "comparison_label_y_position": [2, 0.105, 0.046], #[2, 0.105, 0.046], + "comparison_label_x_position": [4.5e15, 7e14, 7e14], #[4.5e15, 2.2e16, 1e15], + #"sigmae_theta_r": 0.159, + "lower_label_y_position": 0.4, + "upper_label_y_position": 4, + "label_x_position": 4e14, #4e14, #0.02, #1e14, + "goals_x_position": 1.2e14, #0.0002 + "plot_key_parameters": True + } +#sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +#sens_curve.Configure(sens_config_dict) +#sens_curve.Run() + + +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam.cfg", #"/termite/sensitivity_config_files/Config_LFA_Experiment.cfg", + "plot_path": "./lfa_with_BNL_constraints_sensitivity_vs_exposure_curve_PhaseIV.pdf", + "exposure_axis": True, + # optional + "figsize": (10,6), + "fontsize": 15, + "legend_location": "upper right", + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "density_axis": False, + "cavity": True, + "add_1year_1cav_point_to_last_ref": False, + "y_limits": [10e-3, 500], + "density_range": [1e12,1e19], + "exposure_range": [1e-11, 1e4], + #"main_curve_upper_label": r"LFA (atomic T)$\,-\,$1 GHz", + "main_curve_upper_label": r"Phase III scenario: 1 GHz, L/D=9", + #"main_curve_lower_label": r"$\sigma^\bar{B}_\mathrm{reco} = 0.16\,\mathrm{eV}$", + "goals": {"Phase III (0.4 eV)": (0.4**2/1.64), "Phase IV (0.04 eV)": (0.04**2/1.64)}, + "comparison_curve": False, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam.cfg", + "/termite/sensitivity_config_files/Config_LFA_Experiment_max_BNL_diam.cfg"], + # "/termite/sensitivity_config_files/Config_atomic_find_factor_22_Experiment_conservative.cfg"], + #"comparison_config_file_path": ["/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment_conservative.cfg", + # "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment_conservative.cfg", + # "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + #"comparison_curve_label": [r"Molecular, conservative", "Atomic, conservative", "Atomic, reaching PIV target"], + "comparison_curve_label": [r"Phase III scenario: 500 MHz, L/D=5", r"Phase III scenario: 560 MHz, L/D=8"], #r"Phase IV scenario: 150 MHz, L/D=5"], + #"comparison_curve_colors": ["blue", "darkred", "red"], + "comparison_curve_colors": ["blue", "darkviolet"], + "optimize_main_density": False, + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 0.015, #1.5e15, #0.02, #1e14, + "goals_x_position": 0.2e-10, #2e12 #0.0002 + "goals_y_rel_position": 0.4, + "add_PhaseII": False + #"PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg" + } +sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") +sens_curve.Configure(sens_config_dict) +sens_curve.Run() + + + + + + diff --git a/test_analysis/Sensitivity_parameter_scan.py b/test_analysis/Sensitivity_parameter_scan.py new file mode 100644 index 00000000..17d9eeff --- /dev/null +++ b/test_analysis/Sensitivity_parameter_scan.py @@ -0,0 +1,163 @@ +""" +Script to make Sensitivty plots for cavity experiments +Author: C. Claessens, T. Weiss +Date: October 6, 2023 +""" + + +from morpho.utilities import morphologging, parser +logger = morphologging.getLogger(__name__) + +from mermithid.processors.Sensitivity import SensitivityParameterScanProcessor + +import numpy as np + +# import all the scanned parameter units +from numericalunits import e, me, c0, eps0, kB, hbar +from numericalunits import meV, eV, keV, MeV, cm, m, mm +from numericalunits import nT, uT, mT, T, mK, K, C, F, g, W +from numericalunits import hour, year, day, ms, ns, s, Hz, kHz, MHz, GHz +deg = np.pi/180 + + + + +# Configuration for magnetic homogeneity scan +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_find_factor_22_Experiment_conservative.cfg", + "plot_path": ".", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": False, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,3e18], + #"density_range": [1e8, 1e12], + "goals": {"Phase IV (0.04 eV)": 0.04}, + #"scan_parameter_name": "MagneticField.sigmae_r", + #"scan_parameter_range": [0.1, 5], + #"scan_parameter_steps": 4, + #"scan_parameter_unit": eV, + "scan_parameter_name": "Experiment.l_over_d", + "scan_parameter_range": [5, 15], + "scan_parameter_steps": 5, + "scan_parameter_unit": 1, + "plot_sensitivity_scan_on_log_scale": False, + "label_x_position": 4e14, #4e14, #0.02, #1e14, + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": False + } +#sens_scan = SensitivityParameterScanProcessor("sensitivity_curve_processor") +#sens_scan.Configure(sens_config_dict) +#sens_scan.Run() + + +# Configuration for minimum pitch angle scan +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "plot_path": ".", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": False, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,3e18], + #"density_range": [1e8, 1e12], + "goals": {"Phase IV (0.04 eV)": 0.04}, + #"scan_parameter_name": "FrequencyExtraction.minimum_angle_in_bandwidth", + #"scan_parameter_range": [83, 89], + #"scan_parameter_steps": 10, + #"scan_parameter_unit": deg, + "scan_parameter_name": "Experiment.l_over_d", + "scan_parameter_range": [5, 15], + "scan_parameter_steps": 5, + "scan_parameter_unit": 1, + "plot_sensitivity_scan_on_log_scale": False, + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": False + } +#sens_scan = SensitivityParameterScanProcessor("sensitivity_curve_processor") +#sens_scan.Configure(sens_config_dict) +#sens_scan.Run() + + +# Configuration for L over D scan +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "plot_path": ".", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": False, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,3e18], + #"density_range": [1e8, 1e12], + "goals": {"Phase IV (0.04 eV)": 0.04}, + #"scan_parameter_name": "Experiment.l_over_d", + #"scan_parameter_range": [5, 8], + #"scan_parameter_steps": 3, + #"scan_parameter_unit": 1, + "scan_parameter_name": "FrequencyExtraction.sideband_power_fraction", + "scan_parameter_range": [0, .4], + "scan_parameter_steps": 9, # This is the one that currently has text output + "scan_parameter_unit": 1, + "plot_sensitivity_scan_on_log_scale": False, + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": False + } +#sens_scan = SensitivityParameterScanProcessor("sensitivity_curve_processor") +#sens_scan.Configure(sens_config_dict) +#sens_scan.Run() + +# Configuration for efficiency scan +sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + #"config_file_path": "/termite/sensitivity_config_files/Config_atomic_find_factor_22_Experiment_conservative.cfg", + "plot_path": ".", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": False, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,3e18], + #"density_range": [1e8, 1e12], + "goals": {"Phase IV (0.04 eV)": 0.04}, + "scan_parameter_name": "MagneticField.sigmae_r", + "scan_parameter_range": [.02,0.18], + "scan_parameter_steps": 5, + "scan_parameter_scale": "lin", + "scan_parameter_unit": eV, + "plot_sensitivity_scan_on_log_scale": False, + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": False + } +sens_scan = SensitivityParameterScanProcessor("sensitivity_curve_processor") +sens_scan.Configure(sens_config_dict) +sens_scan.Run() + +print(sens_scan.results) + + + + + + + + + diff --git a/tests/Misc_test.py b/tests/Misc_test.py index d1f05019..33bd7c73 100644 --- a/tests/Misc_test.py +++ b/tests/Misc_test.py @@ -11,8 +11,6 @@ from morpho.utilities import morphologging, parser logger = morphologging.getLogger(__name__) -import matplotlib.pyplot as plt -import numpy as np class MiscTest(unittest.TestCase): @@ -39,7 +37,6 @@ def test_FreqConversionTest(self): - if __name__ == '__main__': args = parser.parse_args(False) @@ -54,4 +51,4 @@ def test_FreqConversionTest(self): stderr_lb=args.stderr_verbosity, propagate=False) - unittest.main() \ No newline at end of file + unittest.main() diff --git a/tests/Sensitivity_test.py b/tests/Sensitivity_test.py new file mode 100644 index 00000000..d31cd348 --- /dev/null +++ b/tests/Sensitivity_test.py @@ -0,0 +1,155 @@ +""" +Script to test the Sensitivty processors +Author: C. Claessens, T. Weiss +Date: December 12, 2021 +""" + +import unittest + +from morpho.utilities import morphologging, parser +logger = morphologging.getLogger(__name__) + +import numpy as np + +class SensitivityTest(unittest.TestCase): + + def test_SensitivityCurveProcessor(self): + from mermithid.processors.Sensitivity import CavitySensitivityCurveProcessor + + # Configuration for Sensitivity vs. exposure plot + # Phase III and IV at 325 MHz + # Phase II point and curve for comparison + + sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "plot_path": "./sensitivity_vs_exposure_curve.pdf", + # optional + "figsize": (10,6), + "fontsize": 15, + "legend_location": "upper right", + "track_length_axis": False, + "molecular_axis": False, + "atomic_axis": False, + "density_axis": False, + "cavity": True, + "add_PhaseII": True, + "add_1year_1cav_point_to_last_ref": True, + "PhaseII_config_path": "/termite/sensitivity_config_files/Config_PhaseII_Experiment.cfg", + "y_limits": [10e-3, 500], + "density_range": [1e12,1e19], + "exposure_range": [1e-11, 1e4], + "main_curve_upper_label": r"Molecular, conservative", + "main_curve_lower_label": r"$\sigma^\bar{B}_\mathrm{reco} = 0.07\,\mathrm{eV}$", + "goals": {"Phase III (0.2 eV)": (0.2**2/1.64), "Phase IV (0.04 eV)": (0.04**2/1.64)}, + "comparison_curve": True, + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment_conservative.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + "comparison_curve_label": [r"150Hz, PIII 87 degrees", "Atomic, conservative", "Atomic, reaching PIV target"], + "comparison_curve_colors": ["blue", "darkred", "red"], + #"B_inhomogeneity": np.linspace(0.1, 2.1, 10)*1e-6, + #"B_inhom_uncertainty": 0.01, + "sigmae_theta_r": 0.159, #in eV, energy broadening from theta and r reconstruction + "lower_label_y_position": 0.17, + "upper_label_y_position": 0.7, + "label_x_position": 0.015, #1.5e15, #0.02, #1e14, + "goals_x_position": 0.2e-10, #2e12 #0.0002 + "goals_y_rel_position": 0.4 + } + #sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") + #sens_curve.Configure(sens_config_dict) + #sens_curve.Run() + + + # Configuration for Sensitivity vs. density plot + # Find optimum density for different scenarios + sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_PhaseIII_325MHz_Experiment_conservative.cfg", + "plot_path": "./sensitivity_vs_density_curve.pdf", + # optional + "figsize": (7.0,6), + "track_length_axis": True, + "molecular_axis": True, + "atomic_axis": True, + "density_axis": True, + "cavity": True, + "y_limits": [2e-2, 4], + "density_range": [1e12,3e18], + "main_curve_upper_label": r"Molecular, conservative", + "comparison_curve_label": ["Molecular, reaching target", + "Atomic, conservative", + "Atomic, reaching target"], + "goals": {"Phase III (0.2 eV)": 0.2, "Phase IV (0.04 eV)": 0.04}, + "comparison_curve": True, + #"comparison_curve_colors": ["red"], + "comparison_config_file_path": ["/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment_conservative.cfg", + "/termite/sensitivity_config_files/Config_atomic_325MHz_Experiment.cfg"], + "goals_x_position": 1.2e12, #0.0002 + "plot_key_parameters": False + } + sens_curve = CavitySensitivityCurveProcessor("sensitivity_curve_processor") + sens_curve.Configure(sens_config_dict) + sens_curve.Run() + + + + def test_SensitivityProcessor(self): + from mermithid.processors.Sensitivity import AnalyticSensitivityEstimation + + # Calculates sensitivity for given configuration (no parameter optimization or scanning) + # Uses the non-cavity specific sensitivity calculations + + sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg" + } + # sens = AnalyticSensitivityEstimation("sensitivity_processor") + # sens.Configure(sens_config_dict) + # sens.Run() + + # sens.print_statistics() + # sens.print_systematics() + + # results = sens.results + # logger.info(results) + + def test_ConstantSensitivityCurvesProcessor(self): + from mermithid.processors.Sensitivity import ConstantSensitivityParameterPlots + + # Makes plots of key parameters vs. configured inputs with fixed target sensitivity + # Uses the non-cavity specific sensitivity calculations + + sens_config_dict = { + # required + "config_file_path": "/termite/sensitivity_config_files/Config_atomic_150MHz_minpitch_87deg.cfg", + "sensitivity_target": [0.4**2/1.64]#, 0.7**2/1.64, 1**2/1.64] + } + #sens = ConstantSensitivityParameterPlots("sensitivity_processor") + #sens.Configure(sens_config_dict) + #sens.Run() + + #sens.print_statistics() + #sens.print_systematics() + + + + + +if __name__ == '__main__': + + args = parser.parse_args(False) + + + logger = morphologging.getLogger('morpho', + level=args.verbosity, + stderr_lb=args.stderr_verbosity, + propagate=False) + logger = morphologging.getLogger(__name__, + level=args.verbosity, + stderr_lb=args.stderr_verbosity, + propagate=False) + + unittest.main()