Commit 7f978101 authored by Erik Asbjørn Mikkelsen Jensen's avatar Erik Asbjørn Mikkelsen Jensen
Browse files

added chart of nuclides example

parent addfc762
......@@ -4,8 +4,8 @@ LDLIBS = $$(root-config --glibs)
DATATABLES = ame16/mass16.txt ame16/rct1-16.txt ame16/rct2-16.txt nubase16/nubase2016.txt
EXAMPLE_DIR = examples
EXAMPLE_FIGURES = qec-all.pdf qec-no-estimates.pdf
EXAMPLE_DATAFILES = qec.dat
EXAMPLE_FIGURES = qec-all.pdf qec-no-estimates.pdf beta-delayed.pdf
EXAMPLE_DATAFILES = qec.dat beta-delayed.dat
EXPANSION_DIR = expansion
EXPANSION_EXS = expand-nuchart-QEC expand-nuchart-QECp expand-nuchart-QEC2p expand-nuchart-QECa expand-nuchart-QB2n expand-nuchart-QBa
EXS = dataprinter
......@@ -16,7 +16,7 @@ EXS = dataprinter
default: nuchart.root dataprinter
all: default example1
all: default example1 example2
nuchart.root: treemaker.py $(addprefix $(EXPANSION_DIR)/, $(EXPANSION_EXS))
ifeq ($(ROOTSYS), ) # source thisroot.sh if ROOTSYS is not defined
......@@ -25,7 +25,6 @@ else
python3 treemaker.py
endif
$(foreach exec, $(addprefix $(EXPANSION_DIR)/, $(EXPANSION_EXS)), ./$(exec); )
$(RM) -r __pycache__
example1: $(EXAMPLE_DIR)/graph-example-qec.py # PHONEY target avoids multiple executions with two real targets ('qec-all.pdf', 'qec-no-estimates.pdf'), but results in 'example1' always being executed if make is run on 'all' or 'example1'
python3 $<
......@@ -33,7 +32,15 @@ example1: $(EXAMPLE_DIR)/graph-example-qec.py # PHONEY target avoids multiple ex
$(EXAMPLE_DIR)/graph-example-qec.py: $(EXAMPLE_DIR)/qec.dat
$(EXAMPLE_DIR)/qec.dat: dataprinter nuchart.root Makefile
./$< nuchart.root a A Z QEC QEC_est "A <= 70" "QEC >= 0" > $@
./$< nuchart.root a A Z QEC QEC_est "A <= 70" "QEC > 0" > $@
example2: $(EXAMPLE_DIR)/graph-example-beta-delayed.py
python3 $<
$(EXAMPLE_DIR)/graph-example-beta-delayed.py: $(EXAMPLE_DIR)/beta-delayed.dat
$(EXAMPLE_DIR)/beta-delayed.dat: dataprinter nuchart.root Makefile
./$< nuchart.root a Z N QECp QBn QEC2p QB2n QECa QBa QECp_est QBn_est QEC2p_est QB2n_est QECa_est QBa_est half_life_stbl "Z <= 50" "N <= 50" > $@
dataprinter: dataprinter.cxx
$(CXX) $(CXXFLAGS) -o $@ $< $(LDLIBS)
......
# TODO
# Nuchart
This project extracts mass excesses, various Q-values and separation energies,
half-lives, spins and parities of all nuclides documented in the *AME* and *NUBASE*
data tables of 2016.
The extracted values are stored in a *ROOT* tree called ***nuchart.root***.
* Make example with N:Z:E colour plot (see Uffe Bergmann) using matshow or pcolor or something else in matplotlib
* Make QBa:A figure
* Do TODO in top of *datagetter.py*
It is possible to expand the tree with more values derived from the values already
present in the tree.
For example, one could calculate Q-values of beta-delayed two-neutron emission by
using the Q-values of beta decay and subsequent two-neutron separation energies
stored in the tree.
This specific expansion, among others, is carried out in the C++ files in the
***expansion*** directory.
Those files serve as general templates for doing further user-specified expansions.
As an added bonus, so to say, the ***dataprinter*** executable made in this project
can take any ROOT file containing a tree, and, given various parameters and
conditions on the parameters, produce output in easily human- and machine-interpretable
ASCII format containing just the parameters fulfilling the given conditions.
This enables further manipulation/analysis of the data outside of ROOT.
# DONE
***dataprinter*** is utilised in the ***examples*** directory (via the project
***Makefile***), in which some figures are produced (outside of ROOT, in python)
as examples of the values one can extract from *Nuchart*.
Please see ***dataprinter.cxx*** or perhaps the ***Makefile*** for further information
on how to use it.
### Prerequisites
The prerequisites for producing ***nuchart.root*** and ***dataprinter*** are
* GNU Make
* A fairly recent C++ compiler
* Python 3.0 or above
* ROOT with the PyROOT module installed
In order to produce the example figures, the Python libraries NumPy and Matplotlib
are also required.
### Building
In order to build ***nuchart.root*** and ***dataprinter*** run
```shell script
make
```
in the project's main directory.
In order to also produce the example figures run
```shell script
make all
```
in the project's main directory.
### Project status
While the project is mostly ready to be used, there are still a few details that
need fixing or tweaking.
These are listed below.
#### To do
* Do todo in top of *datagetter.py* and *graph-example-beta-delayed.py*
#### Done
* Everything prior
* Make generic script which writes desired values from root file to a text file
* Combine subprojects into one big project
* Make QEC example
* Tidied root folder up
* Make chart of nuclides example with beta-delayed particle emission (mostly done)
* Everything prior ✔
* Make generic script which writes desired values from root file to a text file ✔
* Combine subprojects into one big project ✔
* Make QEC example ✔
* Tidied root folder up ✔
\ No newline at end of file
import numpy as np
import matplotlib.pyplot as plt
# TODO: some points on the final plot are missing, and some should probably be
# removed based on beta stability considerations.
# this example still serves as a basic draft for making parameter-specific
# chart of nuclides figures, however.
n, Z, N, QECp, QBn, QEC2p, QB2n, QECa, QBa, QECp_est, QBn_est, QEC2p_est, \
QB2n_est, QECa_est, QBa_est, half_life_stbl = np.loadtxt('examples/beta-delayed.dat', unpack=True)
filters = [(half_life_stbl == 1), ((QECp > 0) | (QBn > 0)) & ((QEC2p < 0) | (QB2n < 0)),
((QEC2p > 0) | (QB2n > 0)), ((QECa > 0) | (QBa > 0))]
labels = ['Stable', '$\mathrm{EC\\alpha \ / \ \\beta^{-}\\alpha}$',
'$\mathrm{ECp \ / \ \\beta^{-}n}$',
'$\mathrm{EC2p \ / \ \\beta^{-}2n\ (or\ ECp \ / \ \\beta^{-}n)}$']
rectcenters = [0.5, 0.5, 0.3, 0.3]
rectsizes = [1.0, 1.0, 0.6, 0.6]
linewidths = [0.0, 0.0, 0.0, 0.0]
edgecolors = ['k', 'tomato', 'limegreen', 'dodgerblue', ]
facecolors = ['k', 'tomato', 'limegreen', 'dodgerblue', ]
xmin = np.min(N)
xmax = np.max(N)
ymin = np.min(Z)
ymax = np.max(Z)
ax = plt.gca()
ax.set_xticks(np.arange(xmin, xmax + 1, 10))
ax.set_yticks(np.arange(xmin, xmax + 1, 10))
ax.set_xticks(np.arange(xmin, xmax + 1, 5), minor=True)
ax.set_yticks(np.arange(xmin, xmax + 1, 5), minor=True)
ax.set_axisbelow(True)
ax.grid(which='both', linestyle='--', alpha=0.5)
for i in range(len(filters)):
plt.plot(-100, -100, marker='s', markeredgecolor=edgecolors[i],
markerfacecolor=facecolors[i], linestyle='', linewidth=linewidths[i], label=labels[i])
for x, y in zip(N[filters[i]], Z[filters[i]]):
ax.add_patch(plt.Rectangle((x - rectcenters[i], y - rectcenters[i]),
rectsizes[i], rectsizes[i], edgecolor=edgecolors[i],
facecolor=facecolors[i], linewidth=linewidths[i]))
ax.set_aspect('equal')
plt.xlim(xmin - 1, xmax + 1)
plt.ylim(xmin - 1, xmax + 1)
plt.xlabel('$N$')
plt.ylabel('$Z$')
handles, labels = plt.gca().get_legend_handles_labels() # tak til Ian Hincks: https://stackoverflow.com/questions/22263807/how-is-order-of-items-in-matplotlib-legend-determined\
order = [0, 2, 3, 1]
plt.title('Note: Incomplete figure. See comment in python code.')
plt.legend([handles[idx] for idx in order], [labels[idx] for idx in order])
plt.savefig('examples/beta-delayed.pdf')
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment