EIC code displayed by LXR master/​geant4/​examples/​advanced/​exp_microdosimetry/​1_plot_distributions.py	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-30 09:19:57

 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()

This page was automatically generated by the 2.3.7 LXR engine. The LXR team