Warning, /geant4/examples/extended/exoticphysics/channeling/ch2/analysis_ch2.ipynb is written in an unsupported language. File is not indexed.
0001 {
0002 "cells": [
0003 {
0004 "cell_type": "code",
0005 "execution_count": 1,
0006 "id": "d52f88aa",
0007 "metadata": {
0008 "tags": []
0009 },
0010 "outputs": [],
0011 "source": [
0012 "# Read and plot the simulation results of particle interactions in Oriented Crystals\n",
0013 "# obatined through example ch2, which is baed on G4ChannelingFastSimModel."
0014 ]
0015 },
0016 {
0017 "cell_type": "code",
0018 "execution_count": 2,
0019 "id": "c010d890",
0020 "metadata": {
0021 "tags": []
0022 },
0023 "outputs": [],
0024 "source": [
0025 "import numpy as np\n",
0026 "import pandas as pd\n",
0027 "import matplotlib.pyplot as plt\n",
0028 "import os\n",
0029 "import uproot\n",
0030 "from matplotlib.colors import LogNorm # optional, for log color scaling"
0031 ]
0032 },
0033 {
0034 "cell_type": "code",
0035 "execution_count": 3,
0036 "id": "bf061146",
0037 "metadata": {
0038 "tags": []
0039 },
0040 "outputs": [],
0041 "source": [
0042 "#################################### INPUT FILE #########################################\n",
0043 "# Set path and filename of the simulation file\n",
0044 "G4_sim_path = \"\"\n",
0045 "root_file = \"results\"\n",
0046 "\n",
0047 "# Set whether to save plots (without displaying them) or just display them\n",
0048 "save_fig = True\n",
0049 "fig_path = G4_sim_path\n",
0050 "#########################################################################################"
0051 ]
0052 },
0053 {
0054 "cell_type": "code",
0055 "execution_count": 4,
0056 "id": "3b298eab",
0057 "metadata": {
0058 "tags": []
0059 },
0060 "outputs": [
0061 {
0062 "name": "stdout",
0063 "output_type": "stream",
0064 "text": [
0065 "rf_content: ['crystal', 'detector_primaries', 'detector_photons', 'detector_secondaries', 'missed_crystal'] \n",
0066 "\n"
0067 ]
0068 }
0069 ],
0070 "source": [
0071 "# Create directory where to strore the figures if it does not exist\n",
0072 "if fig_path != '' and not os.path.exists(fig_path):\n",
0073 " os.makedirs(fig_path)\n",
0074 " print('created fig_path:', fig_path)\n",
0075 " \n",
0076 "# Open the simulation output root file \n",
0077 "rf = uproot.open(G4_sim_path + root_file + '.root')\n",
0078 "rf_content = [item.split(';')[0] for item in rf.keys()]\n",
0079 "print('rf_content:', rf_content, '\\n')"
0080 ]
0081 },
0082 {
0083 "cell_type": "code",
0084 "execution_count": 5,
0085 "id": "599eb756",
0086 "metadata": {
0087 "tags": []
0088 },
0089 "outputs": [],
0090 "source": [
0091 "# Import the scoring ntuples and convert them into pandas dataframes\n",
0092 "branches = [\"eventID\", \"volume\", \"x\", \"y\", \"angle_x\", \"angle_y\", \\\n",
0093 " \"Ekin\" , \"particle\", \"particleID\", \"parentID\"]\n",
0094 "branchesprimary = branches + [\"incoming_angle_x\", \"deflection_angle_x\", \\\n",
0095 " \"incoming_angle_y\", \"deflection_angle_y\"]\n",
0096 "\n",
0097 "df_in = rf['crystal'].arrays(branches, library='pd')\n",
0098 "df_prim = rf['detector_primaries'].arrays(branchesprimary, library='pd')\n",
0099 "df_ph = rf['detector_photons'].arrays(branches, library='pd')\n",
0100 "df_sec = rf['detector_secondaries'].arrays(branches, library='pd')\n",
0101 "df_missed = rf['missed_crystal'].arrays(branches, library='pd')"
0102 ]
0103 },
0104 {
0105 "cell_type": "code",
0106 "execution_count": 6,
0107 "id": "07c70dd2",
0108 "metadata": {
0109 "tags": []
0110 },
0111 "outputs": [],
0112 "source": [
0113 "#########################################################################################\n",
0114 "# Plot angle_x distribution of primaries at the detector after interaction with a crystal\n",
0115 "\n",
0116 "############# INPUT #############\n",
0117 "# Feel free to modify according to your needs\n",
0118 "Nmax = 100000000 #max number of events to elaborate\n",
0119 "\n",
0120 "# Feel free to replace df_prim by df_in, df_ph, df_sec or df_missed\n",
0121 "# Feel free to replace \"angle_x\" by other ntuples from branches and\n",
0122 "# from branchesprimary (for df_prim) \n",
0123 "# ONLY NUMERIC VALUES\n",
0124 "datax = df_prim[\"angle_x\"][:Nmax]*1.e3 #mrad <= rad (feel free to modify the coefficient)\n",
0125 "\n",
0126 "# Feel free to modify the number of bins and the plot range\n",
0127 "NbinTheta = 100\n",
0128 "rangeTheta = [-1, 2] #mrad\n",
0129 "\n",
0130 "# Set whether to use linear o log scale\n",
0131 "use_log_y = False # set True for LogNorm color scale\n",
0132 "\n",
0133 "# Feel free to modify the names of axes\n",
0134 "plt_xlabel = '$\\\\theta_x$ [mrad]'\n",
0135 "plt_ylabel = 'PDF: 1/N dN/d$\\\\theta_x$ [mrad]$^{-1}$'\n",
0136 "\n",
0137 "# Feel free to modify the filename to save the plot\n",
0138 "filename = 'thetaXdistribution.pdf'\n",
0139 "\n",
0140 "#some plt parameters\n",
0141 "fs = 16\n",
0142 "lw = 2\n",
0143 "#################################\n",
0144 "\n",
0145 "# Create 1D histogram\n",
0146 "thetaXdistrib, thetaEdges = np.histogram(datax.values, \\\n",
0147 " bins=NbinTheta, range=rangeTheta, density=True)\n",
0148 "thetabin = thetaEdges[:-1] + (thetaEdges[1]-thetaEdges[0])*0.5\n",
0149 "plt.figure(figsize=(9, 6))\n",
0150 "plt.grid()\n",
0151 "plt.plot(thetabin, thetaXdistrib, linewidth=lw, alpha=1, label='')\n",
0152 "plt.xlabel(plt_xlabel, fontsize=fs)\n",
0153 "plt.ylabel(plt_ylabel, fontsize=fs)\n",
0154 "\n",
0155 "# Set log scale\n",
0156 "if use_log_y:\n",
0157 " plt.yscale('log',base=2) \n",
0158 "\n",
0159 "# Save the plot or just show it\n",
0160 "if save_fig:\n",
0161 " plt.savefig(fig_path + filename)\n",
0162 " plt.close() "
0163 ]
0164 },
0165 {
0166 "cell_type": "code",
0167 "execution_count": 7,
0168 "id": "c9f9c18a",
0169 "metadata": {
0170 "tags": []
0171 },
0172 "outputs": [],
0173 "source": [
0174 "#########################################################################################\n",
0175 "# angle_x_in - angle_x_defl distribution of primaries at the detector after interaction with a crystal\n",
0176 "\n",
0177 "############# INPUT #############\n",
0178 "# Feel free to modify according to your needs\n",
0179 "Nmax = 100000000 #max number of events to elaborate\n",
0180 "\n",
0181 "# Example data (replace these with your real arrays)\n",
0182 "# datatetaxin and datatetadeflx must be the same length\n",
0183 "datatetaxin = df_prim[\"incoming_angle_x\"][:Nmax]*1.e3 #mrad <= rad (feel free to modify the coefficient)\n",
0184 "datatetadeflx = df_prim[\"deflection_angle_x\"][:Nmax]*1.e3 #mrad <= rad (feel free to modify the coefficient)\n",
0185 "\n",
0186 "# Feel free to modify the number of bins and the plot range\n",
0187 "NbinTheta = 50\n",
0188 "\n",
0189 "# Feel free to modify the plot range\n",
0190 "xrange = (-0.1, 0.1)\n",
0191 "yrange = (-1, 2)\n",
0192 "\n",
0193 "# Set whether to use linear o log scale\n",
0194 "use_log_color = True # set True for LogNorm color scale\n",
0195 "\n",
0196 "# Feel free to modify the names of axes\n",
0197 "plt_xlabel2 = '$\\\\theta_{x in}$ [mrad]'\n",
0198 "plt_ylabel2 = '$\\\\theta_{x defl}$ [mrad]'\n",
0199 "\n",
0200 "# Feel free to modify the filename to save the plot\n",
0201 "filename2 = 'thetaXin_thetaXdefl.pdf'\n",
0202 "\n",
0203 "#some plt parameters\n",
0204 "fs = 16\n",
0205 "lw = 2\n",
0206 "#################################\n",
0207 "\n",
0208 "# Create 2D histogram\n",
0209 "plt.figure(figsize=(8, 6))\n",
0210 "hist = plt.hist2d(\n",
0211 " datatetaxin,\n",
0212 " datatetadeflx,\n",
0213 " bins=NbinTheta,\n",
0214 " density=True,\n",
0215 " range=[xrange, yrange],\n",
0216 " norm=LogNorm() if use_log_color else None,\n",
0217 " cmap='jet'\n",
0218 ")\n",
0219 "\n",
0220 "# Add colorbar (PDF scale)\n",
0221 "cbar = plt.colorbar()\n",
0222 "cbar.set_label('PDF', fontsize=fs)\n",
0223 "\n",
0224 "# Labels and title\n",
0225 "plt.xlabel(plt_xlabel2, fontsize=fs)\n",
0226 "plt.ylabel(plt_ylabel2, fontsize=fs)\n",
0227 "\n",
0228 "# Save the plot or just show it\n",
0229 "if save_fig:\n",
0230 " plt.savefig(fig_path + filename2)\n",
0231 " plt.close() "
0232 ]
0233 },
0234 {
0235 "cell_type": "code",
0236 "execution_count": null,
0237 "id": "3098a395",
0238 "metadata": {
0239 "tags": []
0240 },
0241 "outputs": [],
0242 "source": [
0243 "################################################################################################\n",
0244 "# Plot spectrum\n",
0245 "\n",
0246 "############# INPUT #############\n",
0247 "# Feel free to modify the collimator parameters\n",
0248 "# !!! related only to real secondary photons from results.root, not from Spectrum.dat\n",
0249 "# For Spectrum.dat, see the options in the simulation macro.\n",
0250 "apply_collimation = True\n",
0251 "coll_angle = 2.3183 #mrad\n",
0252 "\n",
0253 "# Feel free to modify\n",
0254 "NbinE = 20\n",
0255 "rangeE = [0, 10] #MeV\n",
0256 "\n",
0257 "# path of the spectrum file obtained using all the Baier-Katkov integration photons\n",
0258 "BK_spectrum_file = \"Spectrum.dat\"\n",
0259 "#################################\n",
0260 "\n",
0261 "# Array with photon energies and angles\n",
0262 "Eph = df_ph['Ekin'].values #MeV \n",
0263 "\n",
0264 "Nph = len(Eph)\n",
0265 "print(\"number of emitted photons:\", Nph)\n",
0266 "thetaX_ph = df_ph['angle_x'].values*1e3 #rad -> mrad\n",
0267 "thetaY_ph = df_ph['angle_y'].values*1e3 #rad -> mrad\n",
0268 "\n",
0269 "# Take only the photons inside the collimator acceptance\n",
0270 "theta_ph = np.sqrt(thetaX_ph**2 + thetaY_ph**2) \n",
0271 "if apply_collimation: \n",
0272 " thetaX_ph = thetaX_ph[theta_ph <= coll_angle]\n",
0273 " thetaY_ph = thetaY_ph[theta_ph <= coll_angle]\n",
0274 " Eph = Eph[theta_ph <= coll_angle]\n",
0275 " theta_ph = theta_ph[theta_ph <= coll_angle]\n",
0276 "\n",
0277 "# Calculate the scored photon energy spectrum\n",
0278 "spectrum0, EbinEdges = np.histogram(Eph, bins=NbinE, range=rangeE, density=False)\n",
0279 "Ebin = EbinEdges[:-1] + (EbinEdges[1]-EbinEdges[0])*0.5\n",
0280 "stepx = Ebin[1]-Ebin[0]\n",
0281 "Nprimaries = df_in[\"Ekin\"].size\n",
0282 "\n",
0283 "spectrum = spectrum0 / (Nprimaries*stepx)\n",
0284 "spectral_intensity = Ebin * spectrum\n",
0285 "\n",
0286 "# Statistical uncertainties: sqrt(N)\n",
0287 "spectrum_err = np.sqrt(spectrum0) / (Nprimaries * stepx)\n",
0288 "spectral_intensity_err = Ebin * spectrum_err\n",
0289 "\n",
0290 "# Read the spectrum file obtained using all the Bair-Katkov integration photons\n",
0291 "BK_spectrum = np.loadtxt(G4_sim_path+BK_spectrum_file, dtype='float', comments='#', \\\n",
0292 " delimiter=' ', skiprows=1, unpack=True)\n",
0293 "E_ext = BK_spectrum[0]\n",
0294 "S_ext = BK_spectrum[1]\n",
0295 "\n",
0296 "# Plot the photon energy spectrum\n",
0297 "fig = plt.figure(figsize=(13, 6))\n",
0298 "fs = 16\n",
0299 "lw = 2\n",
0300 "bw = 0.6\n",
0301 "color0 = '#B1B3FB'\n",
0302 "\n",
0303 "#!!! The spectrum is normalized on the total radiation probability W_rad \n",
0304 "#(from MinPhotonEnergy to the energy of the primary particle) which is equivalent to the radiation yield.\n",
0305 "#The integral is equal to W_rad, not to 1!\n",
0306 "plt.subplot(1,2,1)\n",
0307 "plt.bar(Ebin, spectrum, width=bw, color=color0, linewidth=lw, alpha=1, label='secondary photons')\n",
0308 "plt.errorbar(Ebin, spectrum, yerr=spectrum_err, fmt='o', color='k', capsize=3)\n",
0309 "plt.xlim(rangeE)\n",
0310 "plt.plot(E_ext, S_ext, 'r-', lw=2.5, label='from '+BK_spectrum_file)\n",
0311 "plt.title('Emitted photon spectrum')\n",
0312 "plt.xlabel('$E$ [MeV]', fontsize=fs)\n",
0313 "plt.ylabel('$dW_{rad}/dE$ [MeV$^{-1}$]', fontsize=fs)\n",
0314 "plt.legend()\n",
0315 "#plt.yscale('log')\n",
0316 "\n",
0317 "#The spectral intensity is the spectrum above multiplied by the energy,\n",
0318 "#so the spectral intensity of bremsstrahlung is nearly constant.\n",
0319 "plt.subplot(1,2,2)\n",
0320 "plt.bar(Ebin, spectral_intensity, width=bw, color=color0, linewidth=lw, alpha=1, label='secondary photons')\n",
0321 "plt.errorbar(Ebin, spectral_intensity, yerr=spectral_intensity_err, fmt='o', color='k', capsize=3)\n",
0322 "plt.plot(E_ext, E_ext * S_ext, 'r-', lw=2.5, label='from '+BK_spectrum_file)\n",
0323 "plt.title('Emitted photon spectral intensity')\n",
0324 "plt.xlabel('$E$ [MeV]', fontsize=fs)\n",
0325 "plt.ylabel('$E dW_{rad}/dE$', fontsize=fs)\n",
0326 "plt.xlim(rangeE)\n",
0327 "plt.legend()\n",
0328 "#plt.yscale('log')\n",
0329 "if save_fig:\n",
0330 " plt.savefig(fig_path + 'spectrum.pdf')\n",
0331 " plt.close() "
0332 ]
0333 },
0334 {
0335 "cell_type": "code",
0336 "execution_count": null,
0337 "id": "7ddc1e1f",
0338 "metadata": {},
0339 "outputs": [],
0340 "source": []
0341 }
0342 ],
0343 "metadata": {
0344 "kernelspec": {
0345 "display_name": "Python 3 (ipykernel)",
0346 "language": "python",
0347 "name": "python3"
0348 },
0349 "language_info": {
0350 "codemirror_mode": {
0351 "name": "ipython",
0352 "version": 3
0353 },
0354 "file_extension": ".py",
0355 "mimetype": "text/x-python",
0356 "name": "python",
0357 "nbconvert_exporter": "python",
0358 "pygments_lexer": "ipython3",
0359 "version": "3.12.9"
0360 }
0361 },
0362 "nbformat": 4,
0363 "nbformat_minor": 5
0364 }
