BGO/overview_nmahepam.py
2026-04-28 11:48:03 +02:00

204 lines
No EOL
6 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import argparse
import sys
import numpy as np
import pandas as pd
import os.path
from scipy.optimize import curve_fit
from fit_functions import lan
# ---------------------------
# Argumente
# ---------------------------
parser = argparse.ArgumentParser()
parser.add_argument("filenumber", type=int)
parser.add_argument("-s", "--sumfit", action='store_true')
parser.add_argument("-e", "--export", action='store_true')
args = parser.parse_args()
# ---------------------------
# File mapping
# ---------------------------
if args.filenumber == 5:
file = "2026-03-10-ahbgoa-test-5"
elif args.filenumber == 7:
file = "2026-04-13-ahbgoa-test-7"
else:
print("Only filenumbers 5 and 7 are supported.")
sys.exit()
filename = "/net/asterix/home/asterix/ava/NMAHEPAM/hists/" + file + "_V"
# ---------------------------
# Layout & Positionen
# ---------------------------
pos = [3,2,1,6,5,4,11,10,9,8,7,14,13,12,17,16,15]
pos = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]
#alle Positionen verwenden
fpos = list(range(1,18))
# ---------------------------
# MAIN
# ---------------------------
def main():
fig = plt.figure(figsize=(19,11))
axes = fig.subplot_mosaic(
"""
..CCBBAA..
.FF.EE.DD.
KKJJIIHHGG
.NN.MM.LL.
..QQPPOO..
"""
)
axs = list(axes.values())
resV = 0.838214 / 2
energy = 'value'
a = 6*[40] + [10]
b = 6*[110] + [310]
startparams = 6*[[200,75,6,0.5,100]] + [[200,75,6,0.5,100]]
fitparams = {'P': [], 'mpv': [], 'err_mpv':[]}
for i in reversed(range(len(axs))):
axs[i].text(10, 1, pos[i], fontsize=14, color="grey")
axs[i].set_xlim([3,150])
axs[i].set_ylim([0.09,15])
axs[i].set_yscale('log')
axs[i].grid(True)
if pos[i] in fpos:
fitparams['P'].append(pos[i])
data = pd.read_csv(filename + "%02d" % pos[i] + ".Bhist", sep='\s+')
Itime = read_Itime(filename)
channels = list(data.columns)[1:]
for c in range(len(channels)):
y = data[channels[c]] * 3600 / resV / Itime
# ---------------------------
# FALL 1: kein sumfit → alles normal plotten
# ---------------------------
if not args.sumfit:
axs[i].step(data[energy], y, label=channels[c])
# ---------------------------
# FALL 2: sumfit → nur "sum" fitten
# ---------------------------
elif args.sumfit and channels[c] == 'SUM':
sig = np.sqrt(data[channels[c]] + 1)
par, cov = curve_fit(
lan,
data[energy][a[c]:b[c]],
data[channels[c]][a[c]:b[c]],
startparams[c],
sigma=sig[a[c]:b[c]],
absolute_sigma=True
)
print(str(channels[c])+': ',par)
fitparams['mpv'].append(par[1])
perr = np.sqrt(np.diag(cov))
fitparams['err_mpv'].append(perr[1])
axs[i].step(
data[energy],
y,
label=f"sum ; MPV={par[1]:.2f} ± {perr[1]:.2f}"
)
axs[i].plot(
data[energy][a[c]:b[c]],
lan(data[energy], *par)[a[c]:b[c]] * 3600 / resV / Itime,
'k--'
)
frame = pd.DataFrame(fitparams)
frame.to_csv('fitparameters/'+file+'_mpv_sums.txt', sep = ' ', index = False)
# ---------------------------
# FALL 3: sumfit aktiv → andere Kanäle nur plotten
# ---------------------------
else:
axs[i].step(data[energy], y)
axs[i].legend(fontsize=7)
# keine x-Ticks außer unten
no_x = [0,1,2,3,4,5,7,8,9,12]
for i in no_x:
plt.setp(axs[i].get_xticklabels(), visible=False)
# keine y-Ticks außer links
no_y = [1,2,4,5,7,8,9,10,12,13,15,16]
for i in no_y:
plt.setp(axs[i].get_yticklabels(), visible=False)
# ---------------------------
# Detektor-Labels (C1C6 im Uhrzeigersinn)
# ---------------------------
# oben
plt.figtext(0.51, 0.9, " B6 ", fontsize=18, ha='center', va='bottom',
bbox=dict(boxstyle="square", facecolor='grey', alpha=0.5))
# unten
plt.figtext(0.51, 0.07, " B3 ", fontsize=18, ha='center', va='top',
bbox=dict(boxstyle="square", facecolor='grey', alpha=0.5))
# links oben
plt.figtext(0.2, 0.75, " B5 ", fontsize=18, ha='center', va='bottom',
bbox=dict(boxstyle="square", facecolor='grey', alpha=0.5),
rotation=40)
# links unten
plt.figtext(0.2, 0.25, " B4 ", fontsize=18, ha='center', va='top',
bbox=dict(boxstyle="square", facecolor='grey', alpha=0.5),
rotation=-40)
# rechts oben
plt.figtext(0.82, 0.75, " B1 ", fontsize=18, ha='center', va='bottom',
bbox=dict(boxstyle="square", facecolor='grey', alpha=0.5),
rotation=-40)
# rechts unten
plt.figtext(0.82, 0.25, " B2 ", fontsize=18, ha='center', va='top',
bbox=dict(boxstyle="square", facecolor='grey', alpha=0.5),
rotation=40)
fig.suptitle(f"Overview position dependency '{file}'", fontsize=18)
fig.supxlabel("signal [mV]")
fig.supylabel("counts [1/h/mV]")
if args.export:
name = f"plots/overview-{args.filenumber}"
plt.savefig(name + ".png")
else:
plt.show()
# ---------------------------
# Itime
# ---------------------------
def read_Itime(file):
f = file + '.Itime'
while len(f) > 0:
if os.path.isfile(f + '.Itime'):
return int(pd.read_csv(f + '.Itime', sep='\s+', names=['time'])['time'][0])
f = f[:-1]
return 1
# ---------------------------
main()