Gabriele Cosmo
138 lines
3.9 KiB
Python
138 lines
3.9 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import glob
|
|
|
|
def Normalize(f, x) :
|
|
dx = x[1:] - x[:-1]
|
|
integral = np.sum(f[:-1] * dx)
|
|
f_norm = f / integral
|
|
return f_norm
|
|
|
|
# When running in MT and outputting to cvs, each thread outputs to a different file.
|
|
# Therefore I need to cycle through each thread and append them all together.
|
|
|
|
energy = np.empty([0,1])
|
|
length = np.empty([0,1])
|
|
|
|
ntuple_name = "radioprotection_nt_102"
|
|
output_name = ntuple_name + "_t" + "*" + ".csv"
|
|
output_list = glob.glob(output_name)
|
|
|
|
for this_file in output_list :
|
|
energy_thread, length_thread = np.loadtxt(this_file, delimiter=',', unpack=True, usecols=(0,1))
|
|
|
|
energy = np.append(energy, energy_thread)
|
|
length = np.append(length, length_thread)
|
|
|
|
# Experimentally, the mean path length is calculated geometrically as a mean chord length.
|
|
# Here for convenience it's taken by averaging the effective path lengths.
|
|
mean_path_length = np.average(length)
|
|
|
|
# A conversion factor can be used to convert the target material to water or tissue equivalent.
|
|
# Its choice depends on the material and can be calculated in different ways.
|
|
# It is suggested that the user replace the following value with his own.
|
|
conversion_factor = 1.
|
|
|
|
# If the previous value is not overriden by the user, this script will attempt to read geometry.mac
|
|
# and provide a factor accordingly
|
|
if conversion_factor == 1. :
|
|
detector = np.array(["Diamond", "MicroDiamond", "Silicon", "SiliconBridge"])
|
|
factor = np.array([0.32, 0.32, 0.57, 0.57]) # conversion factor based on material stopping power
|
|
|
|
with open("geometry.mac") as search:
|
|
for line in search:
|
|
line = line.rstrip() # remove new line
|
|
|
|
for this in detector :
|
|
match = "/geometrySetup/selectDetector " + this
|
|
|
|
if match == line:
|
|
conversion_factor = factor[ detector == this ][0]
|
|
else:
|
|
conversion_factor = factor[ detector == "Diamond" ][0] # default detector type (no macro used)
|
|
|
|
|
|
y = energy * conversion_factor / mean_path_length
|
|
|
|
|
|
# The spectrum is now binned logarithmically, to avoid oscillations at higher energies
|
|
# (due to fewer counts) that wouldn't much meaning.
|
|
|
|
minimum = np.amin(y)
|
|
maximum = np.amax(y)
|
|
|
|
exp_start = np.floor(np.log10(minimum))
|
|
exp_end = np.ceil(np.log10(maximum))
|
|
|
|
n_decades = int(exp_end - exp_start)
|
|
|
|
# Number of logarithmic bins per decade:
|
|
# Higher values give better resolution, but lead to oscillations
|
|
# (especially at high energy) if your statistic has too few counts.
|
|
bins_per_dec = 60
|
|
|
|
n_bins = n_decades * bins_per_dec
|
|
|
|
y_bins = np.zeros(n_bins)
|
|
|
|
y_bins[0] = 10**exp_start
|
|
|
|
for i in range(1, n_bins) :
|
|
y_bins[i] = y_bins[i-1] * 10**( 1 / bins_per_dec )
|
|
|
|
|
|
# Create the histogram
|
|
|
|
# For now f is a number of counts...
|
|
f = np.histogram( y, bins=y_bins ) [0]
|
|
|
|
tot_counts = np.sum(f)
|
|
|
|
# ... so now I turn f into a density
|
|
bin_width = y_bins[1:] - y_bins[:-1]
|
|
f = f / bin_width
|
|
|
|
f = np.append(f, 0.) # give f and y_bins arrays the same size
|
|
|
|
# Normalize the spectra to unit area under the curve
|
|
f = Normalize(f, y_bins)
|
|
d = Normalize(y_bins*f, y_bins)
|
|
|
|
|
|
# Save to file
|
|
|
|
output_file = "analysed_spectra.csv"
|
|
header = "y[keV/um], f(y)[um/keV], d(y)[um/keV]"
|
|
|
|
np.savetxt( output_file, np.c_[ y_bins, f, d ], header=header, delimiter=',' )
|
|
|
|
|
|
# Plot
|
|
|
|
fig, ax1 = plt.subplots()
|
|
|
|
color = 'tab:blue'
|
|
|
|
ax1.semilogx( y_bins, y_bins * f, linewidth=0.5, color=color )
|
|
|
|
ax1.set_xlabel(r'$y \,\, [keV / \mu m]$')
|
|
ax1.set_ylabel(r'$y \cdot f(y) $', color=color)
|
|
ax1.tick_params(axis='y', labelcolor=color)
|
|
|
|
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
|
|
|
|
color = 'tab:red'
|
|
|
|
ax2.semilogx( y_bins, y_bins * d, linewidth=0.5, color=color )
|
|
|
|
ax2.set_ylabel(r'$y \cdot d(y) $', color=color)
|
|
ax2.tick_params(axis='y', labelcolor=color)
|
|
|
|
title = str(tot_counts) + " counts, " + str(bins_per_dec) + " bins per decade, " + str(conversion_factor) + " conversion factor"
|
|
fig.suptitle(title)
|
|
|
|
fig.tight_layout()
|
|
|
|
plt.subplots_adjust(top=0.92)
|
|
plt.show()
|