...
 
Commits (5)
# HDTV Change Log
## 18.04
- Use uncertainties library (@opapst)
- Large overhaul of compiled libraries (@nima) (Closes #21)
- Fix efficiency fit and nuclide data (@jmayer)
- Other small fixes
## 18.01
- Changed building process in setup.py to allow for package building
......
......@@ -9,6 +9,7 @@ try:
except ImportError:
import urllib
import hdtv.ui
from uncertainties import ufloat
from hdtv.database.common import *
......@@ -16,13 +17,8 @@ def SearchNuclide(nuclide):
"""
Opens table of nuclides with peak energies, gives back the peak energies of the nuclide and its intensities.
"""
Energies = [] # the energies of the right nuclide will be saved in this list
EnergiesError = []
Intensities = []
IntensitiesError = []
Halflive = 0
HalfliveError = 0
source = "DDEP"
out = {'nuclide': nuclide, 'transitions': []}
if nuclide == "Ra-226":
nuclide = "Ra-226D"
......@@ -38,33 +34,33 @@ def SearchNuclide(nuclide):
except:
raise hdtv.cmdline.HDTVCommandError("Error looking up nuclide {}".format(nuclide))
for line in data.split("\r\n"):
sep = line.split(" ; ")
if str(sep[0]) == "Half-life (s)":
Halflive = float(sep[1])
HalfliveError = float(sep[2])
out['halflife'] = ufloat(float(sep[1]), float(sep[2]))
if str(sep[0]) == "Reference":
source = str(sep[1])
source = source.replace(source[-1], "")
source = source.replace(source[-1], "")
out['reference'] = str(sep[1])
try:
if str(sep[4]) == "g":
Energies.append(float(sep[0]))
energy = ufloat(float(sep[0]), 0)
try:
EnergiesError.append(float(sep[1]))
energy.std_dev = float(sep[1])
except BaseException:
EnergiesError.append(0)
# because it is given in %
Intensities.append(float(sep[2]) / 100)
pass
intensity = ufloat(float(sep[2]) / 100, 0)
try:
IntensitiesError.append(
float(sep[3]) / 100) # because it is given in %
intensity.std_dev = float(sep[3]) / 100
except BaseException:
IntensitiesError.append(0)
except BaseException:
pass
out['transitions'].append({'energy': energy, 'intensity' : intensity})
except:
pass
return Energies, EnergiesError, Intensities, IntensitiesError, Halflive, HalfliveError, source
return out
......@@ -4,44 +4,28 @@
import json
import hdtv.util
import hdtv.ui
from uncertainties import ufloat
from hdtv.database.common import *
# TODO: integrate this in a class?
def SearchNuclide(nuclide):
"""
Opens table of nuclides with peak energies, gives back the peak energies of the nuclide and its intensities.
"""
# opens json File with a list of nuclieds and their energies
tableOfEnergiesFile = open(os.path.join(hdtv.datadir, "IAEA.json"))
tableOfEnergiesStr = tableOfEnergiesFile.read() # Sting oft the json file
Energies = [] # the energies of the right nuclide will be saved in this list
EnergiesError = []
Intensities = []
IntensitiesError = []
Halflive = 0
HalfliveError = 0
source = "IAEA"
# compares the name of every nuclide with the given one and saves the
# energies in the list
for tableOfEnergiesData in json.loads(tableOfEnergiesStr):
if tableOfEnergiesData['nuclide'] == nuclide:
Halflive = float(tableOfEnergiesData['halflife'])
HalfliveError = float(tableOfEnergiesData['halflife_uncertainty'])
Energies = [float(data['energy'])
for data in tableOfEnergiesData['transitions']]
EnergiesError = [float(data['energy_uncertainty'])
for data in tableOfEnergiesData['transitions']]
Intensities = [float(data['intensity'])
for data in tableOfEnergiesData['transitions']]
IntensitiesError = [float(data['intensity_uncertainty'])
for data in tableOfEnergiesData['transitions']]
return Energies, EnergiesError, Intensities, IntensitiesError, Halflive, HalfliveError, source
break
# Get all data as dict from json file
alldata = json.loads(open(os.path.join(hdtv.datadir, "IAEA.json")).read())
# error message if nuclide is not in the table
if Energies == []:
# Select nuclide
data = next(x for x in alldata if x['nuclide'] == nuclide)
if not data:
errorText = "There is no nuclide called " + nuclide + " in the table."
raise hdtv.cmdline.HDTVCommandError(errorText)
# Use ufloat to represent values with uncertainties
transitions = [{'energy': ufloat(t['energy'], t['energy_uncertainty']),
'intensity': ufloat(t['intensity'], t['intensity_uncertainty'])}
for t in data['transitions']]
data['transitions'] = transitions
data['halflife'] = ufloat(data['halflife'], data['halflife_uncertainty'])
del data['halflife_uncertainty']
return data
......@@ -146,9 +146,9 @@ class _Efficiency(object):
hasXerrors = False
# Convert energies to array needed by ROOT
try:
list(map(lambda x: E.append(x[0].value.nominal_value), self._fitInput))
list(map(lambda x: E.append(x[0].nominal_value), self._fitInput))
list(map(lambda x: delta_E.append(
x[0].value.std_dev), self._fitInput))
x[0].std_dev), self._fitInput))
list(map(lambda x: EN.append(0.0), self._fitInput))
hasXerrors = True
except AttributeError: # energies does not seem to be ufloat list
......@@ -157,9 +157,9 @@ class _Efficiency(object):
# Convert efficiencies to array needed by ROOT
try:
list(map(lambda x: eff.append(x[1].value.nominal_value), self._fitInput))
list(map(lambda x: eff.append(x[1].nominal_value), self._fitInput))
list(map(lambda x: delta_eff.append(
x[1].value.std_dev), self._fitInput))
x[1].std_dev), self._fitInput))
list(map(lambda x: effN.append(0.0), self._fitInput))
except AttributeError: # energies does not seem to be ufloat list
list(map(lambda x: eff.append(x[1]), self._fitInput))
......
......@@ -48,38 +48,23 @@ def SearchNuclide(nuclide, database):
return data
def TabelOfNuclide(nuclide, Energies, EnergiesError, Intensities,
IntensitiesError, Halflive, HalfliveError, source):
def TableOfNuclide(data):
"""
Creates a table of the given data.
"""
tabledata = list()
Data = [[], []]
data = ufloat(Halflive, HalfliveError)
for j in range(0, len(Energies)):
Data[0].append(ufloat(Energies[j], EnergiesError[j]))
Data[1].append(ufloat(
Intensities[j], IntensitiesError[j]))
# for table option values are saved
tableline = dict()
tableline["Energy"] = Data[0][j]
tableline["Intensity"] = Data[1][j]
tabledata.append(tableline)
result_header = "Data of the nuclide " + \
str(nuclide) + " of the data source " + \
str(source) + "." + "\n" + "Halflife: " + str(data)
print()
result_header = ' '.join(["\n",
"Nuclide:", data['nuclide'], "\n",
"Halflife:", str(data['halflife']), "\n",
"Reference:", data['reference'], "\n",
])
table = hdtv.util.Table(
data=tabledata,
keys=[
"Energy",
"Intensity"],
data=data['transitions'],
keys=["energy", "intensity"],
extra_header=result_header,
sortBy=None,
ignoreEmptyCols=False)
hdtv.ui.msg(str(table))
......@@ -136,20 +121,13 @@ def MatchPeaksAndEnergies(peaks, energies, sigma):
return(accordance)
# Peak is the Energy of the fitted Peak, Vol its volume
def MatchPeaksAndIntensities(
peaks, peak_ids, energies, intensities, intensities_error, sigma=0.5):
# and Intensity and Energy the data from the chart
def MatchFitsAndTransitions(fits, transitions, sigma=0.5):
"""
Combines peaks with the right intensities from the table (with searchEnergie).
Combines peaks with the right intensities.
"""
return list(zip([
[
[
peak,
ufloat(intensity, intensity_error),
peak_id,
] for peak, peak_id in zip(peaks, peak_ids)
if abs(energy - peak) <= sigma
] for energy, intensity, intensity_error in zip(
energies, intensities, intensities_error)]))
return [
{'fit': fit, 'transition': transition}
for fit in fits
for transition in transitions
if abs(transition['energy'] - fit.ExtractParams()[0][0]['pos']) <= sigma
]
This diff is collapsed.
......@@ -66,9 +66,10 @@ void DisplayObj::Remove(View1D *view) { Remove(view->GetDisplayStack()); }
//! Add the object to the display stack
void DisplayObj::Draw(DisplayStack *stack) {
auto &objects = stack->fObjects;
auto zindex = GetZIndex();
auto pos = std::find_if(
objects.begin(),
objects.end(), [zindex = GetZIndex()](const DisplayObj *obj) {
objects.end(), [&zindex](const DisplayObj *obj) {
return obj->GetZIndex() > zindex;
});
objects.insert(pos, this);
......@@ -101,9 +102,10 @@ void DisplayObj::ToBottom(View1D *view) { ToBottom(view->GetDisplayStack()); }
//! display stack.
void DisplayObj::ToTop(DisplayStack *stack) {
auto &objects = stack->fObjects;
auto zindex = GetZIndex();
auto pos = std::find_if(
objects.begin(),
objects.end(), [zindex = GetZIndex()](const DisplayObj *obj) {
objects.end(), [&zindex](const DisplayObj *obj) {
return obj->GetZIndex() > zindex;
});
if (pos != objects.end() && *pos == this) {
......@@ -128,9 +130,10 @@ void DisplayObj::ToTop() {
void DisplayObj::ToBottom(DisplayStack *stack) {
auto &objects = stack->fObjects;
auto zindex = GetZIndex();
auto pos = std::find_if(
objects.begin(),
objects.end(), [zindex = GetZIndex()](const DisplayObj *obj) {
objects.end(), [&zindex](const DisplayObj *obj) {
return !(obj->GetZIndex() < zindex);
});
if (pos != objects.end() && *pos == this) {
......
......@@ -53,8 +53,8 @@ void Painter::DrawFunction(DisplayFunc *dFunc, int x1, int x2) {
double norm = fUseNorm ? dFunc->GetNorm() : 1.0;
// Do x axis clipping
x1 = std::min(x1, EtoX(dFunc->GetMinE()));
x2 = std::max(x2, EtoX(dFunc->GetMaxE()));
x1 = std::max(x1, EtoX(dFunc->GetMinE()));
x2 = std::min(x2, EtoX(dFunc->GetMaxE()));
int ly, cy;
ch = dFunc->E2Ch(XtoE(x1 - 0.5));
......@@ -93,8 +93,8 @@ void Painter::DrawSpectrum(DisplaySpec *dSpec, int x1, int x2) {
int lClip = fYBase;
// Do x axis clipping
x1 = std::min(x1, EtoX(dSpec->GetMinE()));
x2 = std::max(x2, EtoX(dSpec->GetMaxE()));
x1 = std::max(x1, EtoX(dSpec->GetMinE()));
x2 = std::min(x2, EtoX(dSpec->GetMaxE()));
switch (fViewMode) {
case kVMSolid:
......
......@@ -86,19 +86,32 @@ def RebuildLibraries(libdir):
for name in modules:
BuildLibrary(name, libdir)
def PrepareBuild(libdir):
# Create base directory
if not os.path.exists(libdir):
os.makedirs(libdir)
utildir = os.path.join(libdir, "util")
if os.path.exists(utildir):
shutil.rmtree(utildir)
srcdir = os.path.dirname(__file__)
shutil.copytree(os.path.join(srcdir, "util"), utildir)
shutil.copy(os.path.join(srcdir, "Makefile.def"), libdir)
shutil.copy(os.path.join(srcdir, "Makefile.body"), libdir)
def BuildLibrary(name, libdir):
PrepareBuild(libdir)
dir = os.path.join(libdir, name)
hdtv.ui.info("Rebuild library %s in %s" % ((libfmt % name), dir))
# Create base directory
if not os.path.exists(libdir):
os.makedirs(libdir)
# Copy everything
# Remove existing plugin
if os.path.exists(dir):
shutil.rmtree(dir)
# Copy everything
shutil.copytree(os.path.join(os.path.dirname(__file__), name), dir)
shutil.copyfile(os.path.join(os.path.dirname(__file__), name, "../Makefile.def"), os.path.join(dir, "../Makefile.def"))
shutil.copyfile(os.path.join(os.path.dirname(__file__), name, "../Makefile.body"), os.path.join(dir, "../Makefile.body"))
# Make library
subprocess.check_call(['make', 'clean', '-j', '--silent'], cwd=dir)
subprocess.check_call(['make', '-j', '--silent'], cwd=dir)
......
......@@ -284,8 +284,8 @@ double EEFitter::EvalBg(const double *x, const double *p) const {
auto last = first - fIntBgDeg - 1;
return sum + std::accumulate(std::reverse_iterator<const double *>(first),
std::reverse_iterator<const double *>(last),
0.0, [x = *x](double bg, double param) {
return std::fma(bg, x, param);
0.0, [&x](double bg, double param) {
return std::fma(bg, *x, param);
});
}
......
......@@ -233,7 +233,7 @@ double TheuerkaufFitter::Eval(const double *x, const double *p) const {
sum += std::accumulate(
std::reverse_iterator<const double *>(p + fNumParams),
std::reverse_iterator<const double *>(p + fNumParams - fIntBgDeg - 1),
0.0, [x = *x](double bg, double param) { return bg * x + param; });
0.0, [&x](double bg, double param) { return bg * *x + param; });
// Evaluate peaks
return std::accumulate(fPeaks.begin(), fPeaks.end(), sum,
......@@ -252,7 +252,7 @@ double TheuerkaufFitter::EvalBg(const double *x, const double *p) const {
sum += std::accumulate(
std::reverse_iterator<const double *>(p + fNumParams),
std::reverse_iterator<const double *>(p + fNumParams - fIntBgDeg - 1),
0.0, [x = *x](double bg, double param) { return bg * x + param; });
0.0, [&x](double bg, double param) { return bg * *x + param; });
// Evaluate steps in peaks
return std::accumulate(fPeaks.begin(), fPeaks.end(), sum,
......
# -*- coding: utf-8 -*-
VERSION = "18.01"
VERSION = "18.04"
......@@ -21,6 +21,8 @@
import sys
import os
import tempfile
import shutil
sys.path.append(os.path.join(os.path.dirname(__file__), 'helpers'))
......@@ -28,3 +30,13 @@ def pytest_configure():
print('Update Root Include Path ...')
import hdtv.rootext
hdtv.rootext.UpdateRootIncludePath()
print('Force Library Rebuild ...')
os.environ["XDG_CACHE_HOME"] = tempfile.mkdtemp()
import hdtv.rootext.dlmgr
hdtv.rootext.dlmgr.RebuildLibraries(hdtv.rootext.dlmgr.usrlibdir)
def pytest_sessionfinish(session, exitstatus):
tmpdir = os.getenv("XDG_CACHE_HOME")
if tmpdir != "" and os.path.exists(tmpdir):
# fingers crossed
shutil.rmtree(tmpdir)
......@@ -166,9 +166,9 @@ def test_cmd_cal_pos_read(calfile):
def test_cmd_nuclide(nuclide):
f, ferr = hdtvcmd("nuclide {}".format(nuclide))
assert ferr == ""
assert "Data of the nuclide" in f
assert "Energy" in f
assert "Intensity" in f
assert "Nuclide:" in f
assert "energy" in f
assert "intensity" in f
@pytest.mark.parametrize("parameters, function", [
("0.1,0.2,0.3,0.4,0.5", "pow"),
......