channeling/ch2/analysis_ch2.py

0001 #!/usr/bin/env python
0002 # coding: utf-8
0003
0004 # Read and plot the simulation results of particle interactions in Oriented Crystals
0005 # obatined through example ch2, which is baed on G4ChannelingFastSimModel.
0006
0007 import numpy as np
0008 import pandas as pd
0009 import matplotlib.pyplot as plt
0010 import os
0011 import uproot
0012
0013 ################################### INPUT ############################################
0014 # Set path and filename of the simulation file
0015 G4_sim_path = ""
0016 root_file = "results"
0017
0018 Nmax = 1e5 #max number of events to elaborate
0019
0020 save_fig = True
0021 fig_path = G4_sim_path
0022
0023 apply_collimation = False
0024 coll_angle = 20/8627 #rad
0025 NbinE = 25
0026 rangeE = [0, 10] #MeV
0027
0028 NbinTheta = 100
0029 rangeTheta = [-1, 2] #mrad
0030 ######################################################################################
0031
0032 # Create figure directory if it does not exist
0033 if fig_path != '' and not os.path.exists(fig_path):
0034     os.makedirs(fig_path)
0035     print('created fig_path:', fig_path)
0036
0037 # Open the simulation output root file
0038 rf = uproot.open(G4_sim_path + root_file + '.root')
0039 rf_content = [item.split(';')[0] for item in rf.keys()]
0040 print('rf_content:', rf_content, '\n')
0041
0042 # Import the scoring ntuples and convert them into pandas dataframes
0043 branches = ["eventID", "volume", "x", "y", "angle_x", "angle_y",
0044             "Ekin" , "particle", "particleID", "parentID"]
0045 df_in = rf['crystal'].arrays(branches, library='pd')
0046 df_out = rf['detector'].arrays(branches, library='pd')
0047 df_ph = rf['detector_photons'].arrays(branches, library='pd')
0048
0049 # Define in and out dataframes
0050 df_in_all_primary = df_in[df_in.parentID == 0]
0051 df_out_all_primary = df_out[df_out.parentID == 0]
0052 Nmax = min([int(Nmax), len(df_out_all_primary)])
0053 df_in_primary = df_in_all_primary[:Nmax]
0054 df_out_primary = df_out_all_primary[:Nmax]
0055
0056 # Select only the columns useful for deflection
0057 df_in_primary_sel = df_in_primary[["eventID", "angle_x", "angle_y"]]
0058 df_out_primary_sel = df_out_primary[["eventID", "angle_x", "angle_y"]]
0059 del df_in_primary, df_out_primary
0060
0061 # Array with photon energies and angles
0062 Eph = df_ph['Ekin'].values #MeV
0063 Nph = len(Eph)
0064 print("number of emitted photons:", Nph)
0065 thetaX_ph = df_ph['angle_x'].values*1e3 #rad -> mrad
0066 thetaY_ph = df_ph['angle_y'].values*1e3 #rad -> mrad
0067
0068 # Take only the photons inside the collimator acceptance
0069 theta_ph = np.sqrt(thetaX_ph**2 + thetaY_ph**2)
0070 if apply_collimation:
0071     thetaX_ph = thetaX_ph[theta_ph <= coll_angle]
0072     thetaY_ph = thetaY_ph[theta_ph <= coll_angle]
0073     Eph = Eph[theta_ph <= coll_angle]
0074     theta_ph = theta_ph[theta_ph <= coll_angle]
0075
0076 # Calculate the scored photon energy spectrum
0077 spectrum, EbinEdges = np.histogram(Eph, bins=NbinE, range=rangeE, density=True)
0078 Ebin = EbinEdges[:-1] + (EbinEdges[1]-EbinEdges[0])*0.5
0079 spectral_intensity = Ebin * spectrum
0080
0081 # Plot the photon energy spectrum
0082 fig = plt.figure(figsize=(13, 6))
0083 fs = 16
0084 lw = 2
0085 bw = 0.6
0086 plt.subplot(1,2,1)
0087 plt.bar(Ebin, spectrum, width=bw, linewidth=lw, alpha=1, label='')
0088 plt.title('Emitted photon spectrum')
0089 plt.xlabel('E [MeV]', fontsize=fs)
0090 plt.ylabel('1/N$\\times$dN/dE', fontsize=fs)
0091 plt.yscale('log')
0092 plt.subplot(1,2,2)
0093 plt.bar(Ebin, spectral_intensity, width=bw, linewidth=lw, alpha=1, label='')
0094 plt.title('Emitted photon spectral intensity')
0095 plt.xlabel('E [MeV]', fontsize=fs)
0096 plt.ylabel('1/N$\\times$dW/dE', fontsize=fs)
0097 plt.yscale('log')
0098 if save_fig:
0099     plt.savefig(fig_path + 'spectrum.jpg')
0100 plt.close()
0101
0102 # Plot angle_x distribution at the detector
0103 thetaXdistrib, thetaEdges = np.histogram(df_out_primary_sel["angle_x"].values*1e3, \
0104                                          bins=NbinTheta, range=rangeTheta, density=True)
0105 thetabin = thetaEdges[:-1] + (thetaEdges[1]-thetaEdges[0])*0.5
0106 plt.figure(figsize=(9, 6))
0107 plt.plot(thetabin, thetaXdistrib, linewidth=lw, alpha=1, label='')
0108 plt.xlabel('$\\theta_X$ [mrad]', fontsize=fs)
0109 plt.ylabel('1/N$\\times$dN/d$\\theta_X$', fontsize=fs)
0110 if save_fig:
0111     plt.savefig(fig_path + 'thetaXdistribution.jpg')
0112 plt.close()
0113