dna/jetcounter/plot.py

0001 #!/usr/bin/env python3
0002 # coding: utf-8
0003 # *********************************************************************
0004 # To execute this script simply type at your machine's prompt:
0005 #   python plot.py
0006 # OR
0007 #   python3 plot.py
0008 # This script needs:
0009 # 1. A bunch of CShistory_t* text files created during simulation.
0010 # 2. (optionally) energy.spectrum file provided with the example.
0011 # *********************************************************************
0012
0013 import os
0014 import numpy as np
0015 import matplotlib
0016 matplotlib.use('Agg')  # to work also in sessions without display
0017 import matplotlib.pyplot as plt
0018
0019
0020 fs = []
0021 for filename in os.listdir():
0022     if "CShistory_t" in filename:
0023         fs.append(np.loadtxt(filename))
0024 events_data = np.concatenate(fs)
0025
0026 nevents = events_data.shape[0]
0027
0028 # ###
0029 # ICSD plot
0030
0031 cluster_sizes = events_data[:, 0]
0032
0033 M1 = cluster_sizes.mean()
0034
0035 bin_edges = np.arange(cluster_sizes.max() + 2)
0036
0037 y, bin_edges = np.histogram(cluster_sizes, bins=bin_edges, density=True)
0038
0039 icsd_fig = plt.figure(figsize=(10, 5))
0040 plt.errorbar(bin_edges[:-1], y, yerr=np.sqrt(y / nevents),
0041              fmt='d', label=f"M₁ = {M1:.3g}")
0042
0043 plt.ylim(.3/cluster_sizes.size, 1)
0044 plt.yscale("log")
0045
0046 plt.title("ICSD")
0047 plt.xlabel("cluster size")
0048 plt.ylabel("probability")
0049 plt.legend()
0050
0051 icsd_fig.savefig("ICSD_py.png")
0052
0053
0054 # ###
0055 # energy spectra
0056 initial_energies = events_data[:, 1]
0057 interaction_energies = events_data[:, 2]
0058 final_energies = events_data[:, 3]
0059
0060 labels = ["initial", "interaction", "final"]
0061
0062 spec_fig = plt.figure(figsize=(10, 5))
0063 for i, energies in enumerate((initial_energies, interaction_energies, final_energies)):
0064     y, bin_edges = np.histogram(energies, bins=int(np.sqrt(nevents)), density=True)
0065     bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2
0066     plt.errorbar(bin_centers, y, yerr=np.sqrt(y / nevents), drawstyle='steps-mid',
0067                  label="{} <E> = {:.3g} MeV".format(labels[i], energies.mean()))
0068
0069
0070 try:
0071     source_spectrum_data = np.loadtxt("energy.spectrum")
0072     length, gain, offset = source_spectrum_data[:3]
0073     cumulative_counts = source_spectrum_data[3:]
0074     counts = np.diff(cumulative_counts)
0075     bins = np.arange(length-1)*gain+offset
0076
0077     normalization_factor = counts.sum()*gain
0078     density = counts/normalization_factor
0079
0080     mean_energy = np.sum(bins*density)*gain
0081
0082     plt.plot(bins, density, drawstyle='steps-mid', label="input <E> = {:.3g} MeV".format(mean_energy))
0083 except FileNotFoundError:
0084     print("Input spectrum file 'energy.spectrum' not found!")
0085     pass
0086
0087
0088 plt.title("Energy spectra")
0089 plt.xlabel("energy [MeV]")
0090 plt.ylabel("probability density [1/MeV]")
0091 plt.legend()
0092
0093 spec_fig.savefig("energy_spectra.png")