{ "cells": [ { "cell_type": "markdown", "id": "b1fed382-d8eb-4a41-a658-eb9e02ae77e6", "metadata": {}, "source": [ "# Read the output of a Geant4 **ch5** simulation" ] }, { "cell_type": "markdown", "id": "29ffa6b6-2469-45d8-bce5-94e8b8d184cf", "metadata": {}, "source": [ "## Import the required libraries" ] }, { "cell_type": "code", "execution_count": null, "id": "3764e2c2-8cb7-4b89-9d5c-1e1379a1207d", "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "#%matplotlib ipympl\n", "from matplotlib.ticker import AutoMinorLocator\n", "import os\n", "import uproot\n", "import math\n", "import re\n", "from scipy import interpolate\n", "\n", "# Set the number of digits to show in pandas dataframes\n", "pd.set_option('display.float_format', '{:.2f}'.format)\n", "\n", "import warnings\n", "warnings.simplefilter(\"ignore\", category=DeprecationWarning)\n", "pd.options.mode.chained_assignment = None" ] }, { "cell_type": "markdown", "id": "cf2a5937-f6d1-4c85-ba11-9889be23b8cd", "metadata": {}, "source": [ "## Graphical settings for fancier plots" ] }, { "cell_type": "code", "execution_count": null, "id": "9044a604-fb1b-4429-9b54-d399919a29b0", "metadata": {}, "outputs": [], "source": [ "import mplhep as hep\n", "\n", "style = hep.style.ROOT\n", "\n", "# figure size\n", "style['figure.figsize'] = (10.795, 7.000)\n", "\n", "# other sizes\n", "style['lines.markersize'] = 6\n", "style['lines.linewidth'] = 3\n", "\n", "# fonts\n", "style['font.family'] = \"serif\"\n", "style['axes.titlesize'] = \"x-small\"\n", "style['axes.labelsize'] = \"small\"\n", "style['xtick.labelsize'] = \"x-small\"\n", "style['ytick.labelsize'] = \"x-small\"\n", "style['legend.fontsize'] = \"x-small\"\n", "style['mathtext.fontset'] = \"dejavuserif\"\n", "\n", "# text locations\n", "style['xaxis.labellocation'] = 'center'\n", "style['yaxis.labellocation'] = 'center'\n", "\n", "# grid\n", "#style[\"axes.grid\"] = True\n", "\n", "# output plots\n", "out_format = \"pdf\"\n", "dpi = 800" ] }, { "cell_type": "markdown", "id": "b12394a6-daae-4b8e-a458-28cae054a7a4", "metadata": { "tags": [] }, "source": [ "## Set input path and base filename" ] }, { "cell_type": "code", "execution_count": null, "id": "cbebf6ec-61ad-4e71-b856-d4afdb3c18d6", "metadata": { "tags": [] }, "outputs": [], "source": [ "path = \"output/\"\n", "\n", "name = \"output_6GeV_W2.0mm_conv8.0mm_D50cm\" #name of the output root file (without extension)\n", "rad_th = 2.0 #[mm] thickness ot the radiator\n", "conv_th = 8.0 #[mm] thickness ot the converter" ] }, { "cell_type": "markdown", "id": "69c5ad7f-f62b-4c95-a1ed-144dda88cec9", "metadata": {}, "source": [ "## Save and export settings" ] }, { "cell_type": "code", "execution_count": null, "id": "9d7a6d4a-e1c4-4a57-a388-5e66a5c6242e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Figure save option\n", "saveFigs = True\n", "\n", "# Option to convert file for RF-Track\n", "convertFileToRFTrack = True\n", "setGaussianTimeShape = True\n", "addID = False" ] }, { "cell_type": "markdown", "id": "60605cbf-1850-4d71-a3f1-d3373e464913", "metadata": {}, "source": [ "## Set output path and desired variables" ] }, { "cell_type": "code", "execution_count": null, "id": "57131664-c6b4-42d5-8c79-a12214c51ae6", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Stripped name\n", "purename = name.replace(\"output_\", \"\")\n", "\n", "# Define the output folder\n", "outpath = path + \"analysis_output/\"\n", "if not os.path.exists(outpath):\n", " os.makedirs(outpath)\n", "\n", "# Quantity to be stored as results\n", "case_list = []\n", "yield_pos_list = []\n", "yield_ele_list = []\n", "yield_ph_list = []\n", "yield_n_list = []\n", "pos_mean_E_list = []\n", "pos_spread_E_list = []\n", "pos_mean_div_fit_list = []\n", "pos_mean_size_fit_list = []\n", "Edep_rad_list = []\n", "Edep_conv_list = []\n", "PEDD_list = []" ] }, { "cell_type": "markdown", "id": "8c779b26-1bc9-49c3-a7b9-7d1289460e12", "metadata": { "tags": [] }, "source": [ "## Read the root file with the scoring ntuples for virtual detectors" ] }, { "cell_type": "code", "execution_count": null, "id": "78e5a25f-0a7c-4af5-ad45-91eca1af2231", "metadata": { "tags": [] }, "outputs": [], "source": [ "print(\"\\nopening rf:\", name ,\".root file ...\")\n", "\n", "rf = uproot.open(path + name + \".root\")\n", "rf_content = [item.split(';')[0] for item in rf.keys()]\n", "print('rf_content:', rf_content)\n", "\n", "#screenID: 2=after the Radiator (if used), 0=before the Converter, 1=after the Converter (or the single-crystal)\n", "#particle: Geant4 particle name\n", "#x,y are in [mm]\n", "#px,py,pz are in [MeV/c]\n", "#t is in [ns]\n", "#branches = [\"screenID\", \"particle\", \"x\", \"y\", \"px\", \"py\", \"pz\", \"t\", \"eventID\"]\n", "df = rf['scoring_ntuple'].arrays(library='pd')\n", "\n", "if 'eventID' in df:\n", " Nevents = len(np.array(df.eventID.unique()))\n", "else:\n", " Nevents = 1e4\n", "print(\"Nevents:\", Nevents) #number of simulated events (primary particles)\n", "\n", "df_rad_out = df[(df.screenID == 2)].copy().drop([\"screenID\"], axis=1)\n", "df_conv_in = df[(df.screenID == 0)].copy().drop([\"screenID\"], axis=1)\n", "df_conv_out = df[(df.screenID == 1)].copy().drop([\"screenID\"], axis=1)\n", "#df_wrong = df[(df.screenID == -1)].copy().drop([\"screenID\"], axis=1)\n", "del df" ] }, { "cell_type": "markdown", "id": "9094f5fd-8bd2-4501-9757-561bfbe73130", "metadata": {}, "source": [ "### Dataframes with phase-spaces (of $\\gamma$/e-/e+/n) after the Radiator crystal\n", "note: if D = 0 -> this is empty!" ] }, { "cell_type": "code", "execution_count": null, "id": "1043904e-ff55-4920-8349-099e6cff9a5c", "metadata": { "tags": [] }, "outputs": [], "source": [ "print(np.array(df_rad_out.particle.unique()))\n", "print(df_rad_out.shape[0], '\\n')\n", "\n", "if df_rad_out.shape[0] > 0 :\n", " df_rad_out_n = df_rad_out[(df_rad_out.particle == \"neutron\")].copy()\n", " df_rad_out_ph = df_rad_out[(df_rad_out.particle == \"gamma\")].copy()\n", " df_rad_out_ele = df_rad_out[(df_rad_out.particle == \"e-\")].copy()\n", " df_rad_out_pos = df_rad_out[(df_rad_out.particle == \"e+\")].copy()\n", " \n", " print('neutrons:', df_rad_out_n.shape[0])\n", " print('photons:', df_rad_out_ph.shape[0])\n", " print('e-:', df_rad_out_ele.shape[0], ', e+:', df_rad_out_pos.shape[0], \\\n", " '-> e-e+-:', df_rad_out_ele.shape[0] + df_rad_out_pos.shape[0], '\\n')\n", "\n", " df_rad_out.head(5)" ] }, { "cell_type": "markdown", "id": "f6a9b843-8841-4a80-b947-f458c97e989f", "metadata": {}, "source": [ "### Dataframes with phase-spaces (of $\\gamma$/e-/e+/n) before the Converter crystal" ] }, { "cell_type": "code", "execution_count": null, "id": "d69706c0-e4d8-496d-850c-bd9e6de72cfc", "metadata": { "tags": [] }, "outputs": [], "source": [ "print(np.array(df_conv_in.particle.unique()))\n", "print(df_conv_in.shape[0], '\\n')\n", "\n", "if df_conv_in.shape[0] > 0 :\n", " df_conv_in_n = df_conv_in[(df_conv_in.particle == \"neutron\")].copy()\n", " df_conv_in_ph = df_conv_in[(df_conv_in.particle == \"gamma\")].copy()\n", " df_conv_in_ele = df_conv_in[(df_conv_in.particle == \"e-\")].copy()\n", " df_conv_in_pos = df_conv_in[(df_conv_in.particle == \"e+\")].copy()\n", "\n", " print('neutrons:', df_conv_in_n.shape[0])\n", " print('photons:', df_conv_in_ph.shape[0])\n", " print('e-:', df_conv_in_ele.shape[0], ', e+:', df_conv_in_pos.shape[0], \\\n", " '-> e-/e+-:', df_conv_in_ele.shape[0] + df_conv_in_pos.shape[0], '\\n')\n", "\n", " df_conv_in.head(5)" ] }, { "cell_type": "markdown", "id": "9dcc1d28-10fa-4f7d-9fa1-914ba08d96d4", "metadata": {}, "source": [ "### Dataframes with phase-spaces (of $\\gamma$/e-/e+/n) after the Converter crystal\n", "note: if the source is not hybrid -> this is empty!" ] }, { "cell_type": "code", "execution_count": null, "id": "f5d9ce34-6ccf-4c4f-83e4-c2ef5cb4795d", "metadata": { "tags": [] }, "outputs": [], "source": [ "print(np.array(df_conv_out.particle.unique()))\n", "print(df_conv_out.shape[0], '\\n')\n", "\n", "if df_conv_out.shape[0] > 0 :\n", " df_conv_out_n = df_conv_out[(df_conv_out.particle == \"neutron\")].copy()\n", " df_conv_out_ph = df_conv_out[(df_conv_out.particle == \"gamma\")].copy()\n", " df_conv_out_ele = df_conv_out[(df_conv_out.particle == \"e-\")].copy()\n", " df_conv_out_pos = df_conv_out[(df_conv_out.particle == \"e+\")].copy()\n", "\n", " print('neutrons:', df_conv_out_n.shape[0])\n", " print('photons:', df_conv_out_ph.shape[0])\n", " print('e-:', df_conv_out_ele.shape[0], ', e+:', df_conv_out_pos.shape[0], \\\n", " '-> e-e+-:', df_conv_out_ele.shape[0] + df_conv_out_pos.shape[0], '\\n')\n", "\n", " df_conv_out.head(5)" ] }, { "cell_type": "markdown", "id": "d40fb008-776e-48d4-bcba-39b36ae7af58", "metadata": {}, "source": [ "## Plot the positron distribution at the Converter/Target exit [and correlate it with the photon distribution at the exit of the radiator crystal]" ] }, { "cell_type": "markdown", "id": "3190fcce-76b8-48e5-92d1-93db9eec9fd1", "metadata": {}, "source": [ "### Import the required libraries" ] }, { "cell_type": "code", "execution_count": null, "id": "7f1578f2-2913-4dc4-8c1a-7556fd84faeb", "metadata": { "tags": [] }, "outputs": [], "source": [ "from scipy.optimize import curve_fit\n", "from matplotlib.colors import LogNorm\n", "\n", "def gaussian(x, a, b, c):\n", " import numpy as np\n", " return a*np.exp(-np.power(x - b, 2)/(2*np.power(c, 2)))\n", "\n", "### Set plot style ###\n", "style['figure.figsize'] = (style['figure.figsize'][0], 8.50) #manual tweak to figure height\n", "plt.style.use(style)" ] }, { "cell_type": "markdown", "id": "ac6e8cb8-9888-4721-8639-71496b8c9208", "metadata": {}, "source": [ "### Options for this calculation" ] }, { "cell_type": "code", "execution_count": null, "id": "6ad51514-f780-44f4-b72b-7f7f18fa1278", "metadata": { "tags": [] }, "outputs": [], "source": [ "if 'conventional' in name:\n", " inputType = \"e-\" \n", "else:\n", " inputType = \"all\" #set a value different from \"all\" only if it is a conventional source\n", "\n", "bCutOnMomentum = False\n", "MomentumCut = 5 #MeV" ] }, { "cell_type": "markdown", "id": "7ae3c2b7-cbfb-4183-953c-ecff0d9d8db0", "metadata": {}, "source": [ "### Define data_in and data_out dataframes and add more parameters to them" ] }, { "cell_type": "code", "execution_count": null, "id": "f1b026f9-ad0c-4d2e-a0e5-df82c318f5bf", "metadata": { "tags": [] }, "outputs": [], "source": [ "if df_rad_out.shape[0] > 0:\n", " data_in = df_rad_out.copy()\n", "else:\n", " if df_conv_out.shape[0] > 0:\n", " data_in = df_conv_in.copy()\n", " else:\n", " data_in = df_rad_out.copy() #set an empty dataframe, because this is a single-volume (crystal or conventional) source\n", " \n", "data_in[\"P\"] = (data_in.px*data_in.px + data_in.py*data_in.py + data_in.pz*data_in.pz)**0.5 #MeV \n", "data_in = data_in[data_in.pz >= 0] #selecting only events (that should be) from the input beam\n", "if bCutOnMomentum:\n", " data_in = data_in[data_in.P >= MomentumCut] #selecting only events with momentum > MomentumCut\n", "data_in.pz = data_in.pz/1000 #MeV -> GeV\n", "data_in.x = data_in.x/10 #mm -> cm\n", "data_in.y = data_in.y/10 #mm -> cm\n", "data_in[\"pt\"] = (data_in.px**2 + data_in.py**2)**0.5 #MeV\n", "data_in[\"thx\"] = np.arctan(data_in.px/1000 / data_in.pz)*1000 #mrad\n", "data_in[\"thy\"] = np.arctan(data_in.py/1000 / data_in.pz)*1000 #mrad\n", "data_in[\"tht\"] = np.arctan(data_in.pt/1000 / data_in.pz)*1000 #mrad\n", "\n", "if df_conv_out.shape[0] > 0:\n", " data_out = df_conv_out.copy() \n", "else: #this is the case of a single-volume (crystalline or conventional) source\n", " data_out = df_conv_in.copy() \n", "\n", "data_out[\"P\"] = (data_out.px*data_out.px+data_out.py*data_out.py+data_out.pz*data_out.pz)**0.5 #MeV\n", "if bCutOnMomentum:\n", " data_out = data_out[data_out.P >= MomentumCut] #selecting only events with momentum > MomentumCut\n", "data_out.pz = data_out.pz/1000 #MeV -> GeV\n", "data_out.x = data_out.x/10 #mm -> cm\n", "data_out.y = data_out.y/10 #mm -> cm\n", "data_out[\"pt\"] = (data_out.px**2 + data_out.py**2)**0.5 #MeV\n", "data_out[\"thx\"] = np.arctan(data_out.px/1000 / data_out.pz)*1000 #mrad\n", "data_out[\"thy\"] = np.arctan(data_out.py/1000 / data_out.pz)*1000 #mrad\n", "data_out[\"tht\"] = np.arctan(data_out.pt/1000 / data_out.pz)*1000 #mrad" ] }, { "cell_type": "markdown", "id": "89f362e5-7855-4bdf-8edd-eea48d2c582b", "metadata": {}, "source": [ "### Plot the angular distribution" ] }, { "cell_type": "code", "execution_count": null, "id": "af9627fe-7a0b-431b-a7a3-0fdb2bc67f08", "metadata": { "tags": [] }, "outputs": [], "source": [ "# plot options\n", "fitrange_ang = 17.45 #mrad (for angular distribution) \n", "ang_plot_lim = (-50, 150) #mrad\n", "fitrange = 0.2 #cm (for spatial distribution)\n", "nBins = 200 #(for spatial distribution)\n", "xRange = (-0.45, 0.45) #cm\n", "xRange1 = (-10.5, 10.5) #cm\n", "yUpperLim = 14 #(for profiles)" ] }, { "cell_type": "code", "execution_count": null, "id": "52dc1aa8-62b6-472d-bf28-dff2063a00df", "metadata": { "tags": [] }, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows=2, sharex=True)\n", "\n", "divergenceIn = [0, 0]\n", "angleStdIn = [0, 0]\n", "divergence = [0, 0]\n", "angleStd = [0, 0]\n", "\n", "if len(data_in) > 0:\n", " series = [data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].thx, \n", " data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].thy]\n", " for i in (0, 1):\n", " hist = ax[i].hist(series[i], bins=15000, range=(-200, 200), histtype=\"step\", color=\"C0\")\n", " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n", " x, y = x[(x>-fitrange_ang) & (x-fitrange_ang) & (x-fitrange_ang) & (x-fitrange_ang) & (x 0:\n", " hist = axdown.hist(data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].x, \\\n", " bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C0\") ;\n", " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n", " x, y = x[(x>-fitrange) & (x-fitrange) & (x-fitrange) & (x-fitrange) & (x 0:\n", " hist = axdown2.hist(data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].y, \\\n", " bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C0\") ;\n", " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n", " x, y = x[(x>-fitrange) & (x-fitrange) & (x-fitrange) & (x-fitrange) & (x 0:\n", " df_ph_rad = data_in[(data_in.particle==\"gamma\")]\n", " nbin_gamma = 100\n", " range_gamma = (50, 500) #(50, 2450)\n", " IWantDensityInSpectrum = False\n", " IWantFluence = False\n", " df_ph_rad_sel = df_ph_rad[(df_ph_rad.P >= range_gamma[0]) & (df_ph_rad.P <= range_gamma[1])]\n", " ph_beam_area = 4*np.pi*np.std(df_ph_rad.x)*np.std(df_ph_rad.y)*100 \n", " spectrum_Eph_rad, edges_Eph_rad = np.histogram(df_ph_rad.P, density=IWantDensityInSpectrum, \\\n", " bins=nbin_gamma, range=range_gamma)\n", " bin_Eph_rad = edges_Eph_rad[:-1] + (edges_Eph_rad[1]-edges_Eph_rad[0])*0.5\n", " if not IWantDensityInSpectrum:\n", " spectrum_Eph_rad = spectrum_Eph_rad / Nevents\n", " if IWantFluence:\n", " spectrum_Eph_rad = spectrum_Eph_rad / ph_beam_area\n", " ylbl = 'ph/mm$^2/e$'\n", " else:\n", " ylbl = 'ph/e'\n", " else:\n", " ylbl = '1/N$\\\\times$dN/dE'\n", " spectral_int_Eph_rad = spectrum_Eph_rad * bin_Eph_rad\n", " plt.figure()\n", " plt.plot(bin_Eph_rad, spectrum_Eph_rad, alpha=0.75, label='')\n", " plt.xlabel('Energy [MeV]')\n", " plt.ylabel(ylbl)\n", " plt.title('spectrum of the photons emitted by the radiator crystal')\n", " plt.grid(True)\n", " plt.yscale('log')\n", " if saveFigs:\n", " plt.savefig(outpath + 'Eph_rad_' + purename + '.jpg')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "18d93d3a-8365-4810-a21f-be495f5d264e", "metadata": { "tags": [] }, "source": [ "## Histograms of Edep per event in the Radiator and the Converter\n", "note: pay attention in case of read_PS...these plots make not much sense!" ] }, { "cell_type": "code", "execution_count": null, "id": "2df91e29-f46d-4dac-85c0-11dcbbc05946", "metadata": { "tags": [] }, "outputs": [], "source": [ "NbinsEdep = 50\n", "fs = 20" ] }, { "cell_type": "markdown", "id": "464d80e7-1622-4c82-b436-e415c5fc63b8", "metadata": {}, "source": [ "### Radiator" ] }, { "cell_type": "code", "execution_count": null, "id": "7566fd8e-98a2-497c-8bd6-a659d9321c2e", "metadata": { "tags": [] }, "outputs": [], "source": [ "if 'edep_rad' in rf_content: \n", " df_edep_rad = rf['edep_rad'].arrays(library='pd')\n", " edep_rad = np.array(df_edep_rad[\"edep\"])\n", " \n", " print(\"total Edep in the Radiator per event:\", round(np.sum(edep_rad)/Nevents, 2), \"MeV/e-\\n\")\n", " \n", " if len(df_edep_rad) > 0:\n", "\n", " edep_rad_cut = edep_rad[edep_rad > 1]\n", "\n", " plt.figure(figsize=(8, 6))\n", " h_rad = plt.hist(edep_rad_cut, bins=NbinsEdep,\n", " histtype=\"step\", label=\"\", density=False,\n", " edgecolor=\"blue\", linewidth=1.2, alpha=1.)\n", " plt.title('Edep per event in the radiator', fontsize=fs)\n", " plt.xlabel(\"Edep [MeV]\", fontsize=fs)\n", " plt.ylabel(\"Counts\", fontsize=fs)\n", " plt.xticks(fontsize=fs, rotation=0)\n", " plt.yticks(fontsize=fs, rotation=0)\n", " plt.grid(which=\"major\", color=\"gray\", linestyle=\"--\", linewidth=1)\n", " if saveFigs:\n", " plt.savefig(outpath + 'Edep_distrib_rad_' + purename + '.jpg')\n", " plt.show() \n", "else:\n", " edep_rad = []" ] }, { "cell_type": "markdown", "id": "a5c0c89a-3e5c-4856-adf1-73cd6128916b", "metadata": {}, "source": [ "### Converter" ] }, { "cell_type": "code", "execution_count": null, "id": "0b11b063-ffe4-4284-994a-54593e9b1daf", "metadata": { "tags": [] }, "outputs": [], "source": [ "if 'edep_conv' in rf_content:\n", " df_edep_conv = rf['edep_conv'].arrays(library='pd')\n", " edep_conv = np.array(df_edep_conv[\"edep\"])\n", " \n", " print(\"total Edep in the Converter per event:\", round(np.sum(edep_conv)*1.e-3/Nevents, 2), \"GeV/e-\\n\")\n", " \n", " if len(df_edep_conv) > 0 and not \"_GT\" in name:\n", "\n", " edep_conv_cut = edep_conv[edep_conv > 0]\n", "\n", " plt.figure(figsize=(8, 6))\n", " h_conv = plt.hist(edep_conv_cut, bins=NbinsEdep,\n", " histtype=\"step\", label=\"\", density=False,\n", " edgecolor=\"blue\", linewidth=1.2, alpha=1.)\n", " plt.title('Edep per event in the converter', fontsize=fs)\n", " plt.xlabel(\"Edep [MeV]\", fontsize=fs)\n", " plt.ylabel(\"Counts\", fontsize=fs)\n", " plt.xticks(fontsize=fs, rotation=0)\n", " plt.yticks(fontsize=fs, rotation=0)\n", " plt.grid(which=\"major\", color=\"gray\", linestyle=\"--\", linewidth=1)\n", " if saveFigs:\n", " plt.savefig(outpath + 'Edep_distrib_conv_' + purename + '.jpg')\n", " plt.show()\n", "else:\n", " edep_conv = []" ] }, { "cell_type": "markdown", "id": "6f646b9a-ec72-4500-adc1-8d998c6683aa", "metadata": { "tags": [] }, "source": [ "## Settings to read the built-in BoxMesh scorer with Edep in the Absorber (Converter or Radiator)" ] }, { "cell_type": "code", "execution_count": null, "id": "dcd7ef10-3916-4ce0-badc-3f140d9aff42", "metadata": { "tags": [] }, "outputs": [], "source": [ "tranvsizeX = 100 #volume tranversal X size (mm)\n", "tranvsizeY = 100 #volume tranversal Y size (mm)\n", "\n", "print(\"tranverse sixe (X*Y) = %.1f x %.1f mm2\" % (tranvsizeX, tranvsizeY))\n", "\n", "if 'conv_th' in locals() and conv_th != '': \n", " tranvsizeZ = conv_th\n", "elif 'rad_th' in locals() and rad_th != '':\n", " tranvsizeZ = rad_th\n", "else:\n", " try:\n", " tranvsizeZ = float(name.split('W')[1].split('mm_')[0])\n", " except:\n", " print('Enter the thickness [mm] of the target (or the Radiator if it is not an hybrid source):')\n", " tranvsizeZ = float(input())\n", " \n", "print(\"tranvsizeZ = %.1f mm\\n\" % (tranvsizeZ))\n", " \n", "IWantPlotVoxel = False\n", "IWantPlotX0fract = False\n", "X0 = 3.504 #crystal radiation length (mm)\n", "\n", "normEvents = True; #normalize w.r.t. simulated events (set always to true to compare with previous results)\n", "\n", "axis_proj = 2 #0: axial, 1:sagittal, 2:coronal\n", "proj_type = 1 #1: mean, 2: sum\n", "p2m = 1 #voxels around which to mean\n", "\n", "plot_proj = True\n", "reader_verb = True\n", "\n", "titstr = ''\n", "\n", "cmap2D = 'viridis'\n", "Wfig = 12\n", "Hfig = 8\n", "fs = 20\n", "ms = 6" ] }, { "cell_type": "markdown", "id": "26ce0632-da84-4982-a944-0680009b2746", "metadata": {}, "source": [ "### Function to read the result provided by the built-in BoxMesh scorer" ] }, { "cell_type": "code", "execution_count": null, "id": "b511c1ad-a6eb-43f6-ac0e-21c88b1ca39a", "metadata": { "tags": [] }, "outputs": [], "source": [ "def read_Edep_BoxMesh(filename, normEvents, Nevents, \n", " tranvsizeX, tranvsizeY, tranvsizeZ):\n", " \"\"\"\n", " Function to read the total Edep accumulated through a built-in BoxMesh scorer.\n", " It is general. It returns [data, (x,y,z)], where:\n", " data is a dictionary with 10 columns defined as follows\n", " {\"ind_x\", \"ind_y\", \"ind_z\", \"eDep\", \"x\", \"y\", \"z\", \"r\" ,\"eDep_err\", \"eDepDensity\", \"eDepDensity_err\"}\n", " and (x,y,z) are the coordinates [mm] of the voxel centers.\n", " \"\"\"\n", " \n", " def find(cond, N=1e7): \n", " \"\"\"\n", " function, equivalent to Matlab's one, that returns\n", " a list with the indices of an array/matrix\n", " that satisfy a given condition.\n", " \"\"\"\n", " cond = cond.reshape(np.size(cond), 1).ravel()\n", " index = list(np.argwhere(cond).ravel())\n", " return index[:N]\n", " \n", " # Retrieve data and put them into a dataframe\n", " data = pd.read_csv(\n", " filename, skiprows=3, names=[\"ind_x\", \"ind_y\", \"ind_z\", \"eDep\", \"eDep2\", \"Nentry\"]\n", " )\n", " \n", " # Calculate eDep uncertainty \n", " data[\"eDep_err\"] = (data[\"eDep2\"] - data[\"eDep\"]**2/data[\"Nentry\"])**0.5\n", " data.fillna(0, inplace=True)\n", " data = data.drop(columns=['eDep2', 'Nentry'])\n", "\n", " # Retrieve the number of voxels in each direction\n", " Nvoxel = len(data)\n", " zzzz = find(np.array(data[\"ind_z\"]) == 0, 2)\n", " if len(zzzz) > 1:\n", " nz = zzzz[1]\n", " else:\n", " nz = len(data[\"ind_y\"]) \n", " ny = np.max([round((np.max(find(np.array(data[\"ind_y\"])==0, nz+1))-1)/nz), 1])\n", " nx = max(np.array(data[\"ind_x\"]))+1\n", " \n", " # Set geometrical size\n", " Xmin = -tranvsizeX*0.5 \n", " Xmax = tranvsizeX*0.5 \n", " Ymin = -tranvsizeY*0.5 \n", " Ymax = tranvsizeY*0.5\n", " Zmin = 0\n", " Zmax = tranvsizeZ\n", "\n", " # Generate vectors of points linearly spaced (corrected in 24/12/2023)\n", " Xedges = np.linspace(Xmin, Xmax, nx+1) #[mm]\n", " Yedges = np.linspace(Ymin, Ymax, ny+1) #[mm]\n", " Zedges = np.linspace(Zmin, Zmax, nz+1) #[mm]\n", " dx = Xedges[1] - Xedges[0]\n", " dy = Yedges[1] - Yedges[0]\n", " dz = Zedges[1] - Zedges[0]\n", " #print(\"dx:\", dx, \"mm, dy:\", dy, \"mm, dz:\", dz, \"mm\")\n", " x = Xedges[:-1] + dx*0.5\n", " y = Yedges[:-1] + dy*0.5\n", " z = Zedges[:-1] + dz*0.5\n", "\n", " # Physical position inside the mesh data\n", " data[\"x\"] = dx * (data[\"ind_x\"] + 0.5) - tranvsizeX*0.5 #central x [mm]\n", " data[\"y\"] = dy * (data[\"ind_y\"] + 0.5) - tranvsizeY*0.5 #central y [mm]\n", " data[\"z\"] = dz * (data[\"ind_z\"] + 0.5) #central z [mm]\n", " data[\"r\"] = np.sqrt(data[\"x\"]**2 + data[\"y\"]**2) #central r [mm]\n", " \n", " # eDep Map\n", " if normEvents:\n", " data[\"eDep\"] = data[\"eDep\"] / Nevents\n", " data[\"eDep_err\"] = data[\"eDep_err\"] / Nevents\n", " data[\"eDepDensity\"] = data[\"eDep\"] / (dz * dy * dx) #[MeV/mm**3] or [MeV/(mm**3 event)]\n", " data[\"eDepDensity_err\"] = data[\"eDep_err\"] / (dz * dy * dx) #[MeV/mm**3] or [MeV/(mm**3 event)]\n", "\n", " return data, (x,y,z)" ] }, { "cell_type": "markdown", "id": "1edbb447-5b72-4c6c-a5d7-a92cbe25d3bd", "metadata": {}, "source": [ "## Read the built-in BoxMesh scorer with Edep in the Absorber" ] }, { "cell_type": "code", "execution_count": null, "id": "b6f22f16-1cdc-44b2-88bb-848fac3e52b3", "metadata": { "tags": [] }, "outputs": [], "source": [ "# read BoxMesh scorer\n", "try:\n", " filename = path + name.replace(\"output\", \"Edep\") + \".txt\"\n", "except:\n", " print('Enter the name of the file with the BoxMesh scorer results:')\n", " filename = float(input())\n", "\n", "results_BoxMesh = read_Edep_BoxMesh(filename, normEvents, Nevents, \n", " tranvsizeX, tranvsizeY, tranvsizeZ)\n", "df_BoxMesh = results_BoxMesh[0]\n", "voxel_coord = results_BoxMesh[1]\n", "del results_BoxMesh\n", "\n", "# retrieve variables form BoxMesh scorer\n", "eDep = df_BoxMesh.eDep\n", "eDepDensity = df_BoxMesh.eDepDensity\n", "x = voxel_coord[0]\n", "y = voxel_coord[1]\n", "z = voxel_coord[2]\n", "if x.shape[0] > 1:\n", " dx = x[1] - x[0]\n", "else:\n", " dx = 0\n", "if y.shape[0] > 1:\n", " dy = y[1] - y[0]\n", "else:\n", " dy = 0\n", "if z.shape[0] > 1:\n", " dz = z[1] - z[0]\n", "else:\n", " dz = 0\n", "\n", "# variables for contour plot of Edep\n", "nx, ny , nz = [voxel_coord[i].shape[0] for i in range(3)]\n", "X, Z = np.meshgrid(x, z)\n", "icx = 0 if nx == 1 else int(X.shape[1]/2+1)\n", "eDepMap = np.transpose(np.array(eDep).reshape(nx,ny,nz)[icx,:,:])\n", "eDepDensityMap = np.transpose(np.array(eDepDensity).reshape(nx,ny,nz)[icx,:,:])\n", "\n", "# Zprofile of Edep\n", "ic = int(Z.shape[1]/2+1)\n", "profileZ_eDep = np.zeros(nz)\n", "profileZ_eDepDensity = np.zeros(nz)\n", "for j in range(nz):\n", " if proj_type == 1:\n", " profileZ_eDep[j] = np.mean(eDepMap[j][ic-p2m:ic+p2m])\n", " profileZ_eDepDensity[j] = np.mean(eDepDensityMap[j][ic-p2m:ic+p2m])\n", " else:\n", " profileZ_eDep[j] = np.sum(eDepMap[j][ic-p2m:ic+p2m])\n", " profileZ_eDepDensity[j] = np.sum(eDepDensityMap[j][ic-p2m:ic+p2m])\n", "\n", "PEDD = np.max(np.array(eDepDensity))\n", "print(\"PEDD:\", round(PEDD,2), \"MeV/(mm^3 e-)\")" ] }, { "cell_type": "markdown", "id": "b4339b43-c89a-4f10-8639-f29416c5845c", "metadata": {}, "source": [ "### Contour and plofile plot of Edep (density) distribution" ] }, { "cell_type": "code", "execution_count": null, "id": "4a89f66c-cf61-4d95-b4a4-e30769fa6f66", "metadata": { "tags": [] }, "outputs": [], "source": [ "# set options\n", "if normEvents:\n", " clbl = \"Edep density per event [MeV/(mm$^3$ e-]\"\n", "else:\n", " clbl = \"Edep density [MeV/mm$^3$]\"\n", "\n", "# contour plot of Edep\n", "if eDepDensityMap.shape[0]>1 and eDepDensityMap.shape[1]>1: \n", " plt.figure(figsize=(Wfig, Hfig))\n", " Nlev = 50\n", " if IWantPlotVoxel:\n", " plt.contourf(eDepDensityMap, Nlev, cmap='viridis')\n", " xlabel = \"X [voxel]\"\n", " ylabel = \"Z [voxel]\" \n", " else:\n", " plt.contourf(X, Z, eDepDensityMap, Nlev, cmap='viridis')\n", " xlabel = \"X [mm]\"\n", " ylabel = \"Z [mm]\"\n", " #plt.gca().set_aspect('equal')\n", " plt.gca().invert_yaxis()\n", " cbar = plt.colorbar()\n", " cbar.ax.tick_params(labelsize=fs)\n", " cbar.set_label(clbl, fontsize=fs, rotation=90)\n", " plt.title(titstr, fontsize=fs)\n", " plt.xlabel(xlabel, fontsize=fs)\n", " plt.ylabel(ylabel, fontsize=fs, wrap=True)\n", " plt.xticks(fontsize=fs, rotation=0)\n", " plt.yticks(fontsize=fs, rotation=0)\n", " plt.gca().xaxis.set_minor_locator(AutoMinorLocator(5))\n", " plt.gca().yaxis.set_minor_locator(AutoMinorLocator(5))\n", " plt.gca().tick_params(axis=\"both\", which='major', direction='in', length=8)\n", " plt.gca().tick_params(axis=\"both\", which='minor', direction='in', length=4)\n", " plt.grid(False)\n", " if saveFigs:\n", " plt.savefig(outpath + 'Edep_proj_' + purename + '.jpg')\n", " plt.show()\n", "\n", "# plot Zprofile of Edep\n", "if IWantPlotVoxel:\n", " xplot = np.linspace(0, nz-1, nz)\n", " ylabel = \"Z [voxel]\" \n", "else:\n", " xplot = z\n", " ylabel = \"Z [mm]\"\n", "plt.figure(figsize=(Wfig, Hfig))\n", "plt.plot(xplot, profileZ_eDepDensity, label=\"\",\n", " color=\"blue\", linestyle='-', linewidth=2, alpha=0.75,\n", " marker='', markersize=ms, markerfacecolor=\"Blue\")\n", "plt.title(titstr, fontsize=fs)\n", "plt.xlabel(ylabel, fontsize=fs, wrap=True)\n", "plt.ylabel(clbl, fontsize=fs, wrap=True)\n", "plt.xticks(fontsize=fs, rotation=0)\n", "plt.yticks(fontsize=fs, rotation=0)\n", "plt.gca().xaxis.set_minor_locator(AutoMinorLocator(5))\n", "plt.gca().yaxis.set_minor_locator(AutoMinorLocator(5))\n", "plt.gca().tick_params(axis=\"both\", which='major', direction='in', length=8)\n", "plt.gca().tick_params(axis=\"both\", which='minor', direction='in', length=4)\n", "plt.grid(True)\n", "if saveFigs:\n", " plt.savefig(outpath + 'Edep_profileZ_' + purename + '.jpg')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "96e92388-7a39-40b7-a1a5-25adfdab1e0c", "metadata": {}, "source": [ "### radius-vs-depth energy density" ] }, { "cell_type": "code", "execution_count": null, "id": "4f6c664d-8140-4b16-aa7d-8cdd4a6d3d1c", "metadata": { "tags": [] }, "outputs": [], "source": [ "# plot options\n", "rLim = 1.5 #mm\n", "bScaleWidth = False\n", "bLog = False\n", "upperLim = 20 #MeV/(mm^3 e-)\n", "ind_r_lsit = [0, 7, 15, 22] #6GeV: [0, 5, 10, 15]" ] }, { "cell_type": "code", "execution_count": null, "id": "414084ec-d20e-4863-b513-d0221ad43e96", "metadata": { "tags": [] }, "outputs": [], "source": [ "# radius-vs-depth energy density\n", "bData = df_BoxMesh[\"r\"] < rLim if rLim != None else df_BoxMesh[\"r\"] < 1e5\n", "Z, R = np.meshgrid(\n", " sorted(list(set(df_BoxMesh[bData][\"z\"]))),\n", " sorted(list(set(df_BoxMesh[bData][\"r\"]))),\n", ")\n", "ED = interpolate.griddata(\n", " (df_BoxMesh[bData][\"z\"], df_BoxMesh[bData][\"r\"]),\n", " df_BoxMesh[bData][\"eDepDensity\"], (Z, R), method='nearest'\n", ")\n", "\n", "if ED.shape[0]>1 and ED.shape[1]>1: \n", " width = 5*max(df_BoxMesh[\"z\"]) if bScaleWidth else 10\n", " fig, ax = plt.subplots(num=\"mesh\", figsize=[width, 8], nrows=2, sharex=True, constrained_layout=True)\n", "\n", " # 2d plot\n", " plotmesh = ax[0].pcolormesh(Z, R, ED, shading=\"nearest\", norm=LogNorm() if bLog else None)\n", " cbar = plt.colorbar(plotmesh, ax=ax[0])\n", " ax[0].grid()\n", " ax[0].set_ylim((0, max(df_BoxMesh[bData][\"r\"])))\n", "\n", " # 1d projection - energy density vs depth\n", " dr = (dx**2 + dy**2)**0.5\n", " for ind_r in ind_r_lsit:\n", " ind_r_eff = int(ind_r * 0.1 / dr)\n", " r_eff = dr * ind_r_eff\n", " bData = (df_BoxMesh[\"r\"] >= (dr * (ind_r_eff - 0.5))) & \\\n", " (df_BoxMesh[\"r\"] <= (dr * (ind_r_eff + 0.5))) \n", " temp = df_BoxMesh[bData][[\"z\", \"eDepDensity\"]].sort_values(by=[\"z\"])\n", " zlist = np.array(temp.z.unique())\n", " ymean = [np.mean(np.array(temp[temp[\"z\"] == zi][\"eDepDensity\"])) for zi in zlist]\n", " ax[1].plot(\n", " #df_BoxMesh[bData][\"z\"], df_BoxMesh[bData][\"eDepDensity\"],\n", " zlist, ymean, #corrected by gpaterno\n", " label = \"r = %.2f mm\" % r_eff,\n", " marker='', linestyle='-'\n", " )\n", " #ax[1].set_ylim((0, upperLim))\n", " ax[1].grid()\n", " ax[1].legend(loc=\"upper left\", fontsize=14)\n", " fs_local = 18\n", "\n", " ax[1].set_xlabel(\"z [mm]\", fontsize=fs_local)\n", " ax[0].set_ylabel(\"r [mm]\", fontsize=fs_local)\n", " ax[1].set_ylabel(r\"eDep [MeV/(mm$^3$ e-)]\", fontsize=fs_local)\n", " cbar.set_label(r\"eDep [MeV/(mm$^3$ e-)]\", rotation=90, labelpad=16, fontsize=fs_local)\n", " #fig.suptitle(\"e- beam features: %.2f GeV, %.1f mm\" % (Ein, sigma_in))\n", " if saveFigs:\n", " plt.savefig(outpath + 'Edep_profileR_' + purename + '.jpg')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "4f319d8c-417b-4e96-a0f3-843939ed5c6b", "metadata": {}, "source": [ "## Get the positrons, analyze and export them" ] }, { "cell_type": "markdown", "id": "26382273-291e-4998-ab0d-1eb30127bb50", "metadata": {}, "source": [ "### Get the positrons" ] }, { "cell_type": "code", "execution_count": null, "id": "3a654a57-cdaf-4eb5-9660-bf3ead04891e", "metadata": { "tags": [] }, "outputs": [], "source": [ "positrons = data_out[data_out.particle==\"e+\"].copy()\n", "print(\"yield_e+:\", round(positrons.shape[0] / Nevents, 2), \"\\n\")\n", "\n", "positrons" ] }, { "cell_type": "markdown", "id": "b4e1a5ad-8e59-4d14-8d84-60789cbe9ad7", "metadata": {}, "source": [ "### Get positrons within a given energy (momentum) threshold" ] }, { "cell_type": "code", "execution_count": null, "id": "5dc3eab2-a137-420d-a1c5-9b2a9b8c1670", "metadata": { "tags": [] }, "outputs": [], "source": [ "p_th = 100 #MeV\n", "positrons[positrons[\"P\"] < 100]" ] }, { "cell_type": "markdown", "id": "f1ae45cb-78ff-4633-ba38-e863e9e64b94", "metadata": {}, "source": [ "### Plot and Save the positron spectrum (taking into account angular and energy cuts)" ] }, { "cell_type": "code", "execution_count": null, "id": "b455dff7-b783-4f28-b64a-c5a08aa55888", "metadata": { "tags": [] }, "outputs": [], "source": [ "save_spectrum = False\n", "\n", "# Calculate positron beam energy distribution features\n", "E_pos = (positrons.P**2 + 0.511**2)**0.5 #MeV\n", "Emean_pos = np.mean(E_pos)\n", "Estdev = np.std(E_pos)\n", "Espread_pos = Estdev / Emean_pos\n", "print(\"Emean_pos:\", round(Emean_pos, 2), \"MeV\")\n", "print(\"Espread_pos:\", round(Espread_pos, 2), \"(sigma/mu)\")\n", "#angle_cut = np.mean(divergence) #mrad\n", "angle_cut = 140 #mrad\n", "positrons_cut = positrons[positrons[\"tht\"] < angle_cut]\n", "#mask = (positrons.x**2+positrons.y**2)**0.5 < np.abs(np.mean([beamSizeX[1], beamSizeY[1]]))\n", "#positrons_cut = positrons[mask] \n", "E_pos_cut = (positrons_cut.P**2 + 0.511**2)**0.5 #MeV\n", "Emean_pos_cut = np.mean(E_pos_cut)\n", "Estdev_cut = np.std(E_pos_cut)\n", "Espread_pos_cut = Estdev_cut / Emean_pos_cut\n", "print(\"Emean_pos_cut (tht < %.2f mrad):\" % angle_cut, round(Emean_pos_cut, 2), \"MeV\")\n", "print(\"Espread_pos_cut (tht < %.2f mrad):\" % angle_cut, round(Espread_pos_cut, 2), \"(sigma/mu)\")\n", "Eth = 100 # MeV\n", "NposEth = len(E_pos.values[E_pos.values < Eth]) / len(E_pos.values)\n", "NsigmaEth = Eth / (Emean_pos + Estdev)\n", "print(\"fraction of positrons with energy < %.0f MeV (%.2f*sigma): %.2f\\n\" % (Eth, NsigmaEth, NposEth))\n", "\n", "# Plot histogram of positrons energy\n", "nbin_pos = 100\n", "range_pos = (0, 100)\n", "IWantDensity = False\n", "plt.figure(figsize=(9, 6))\n", "h = plt.hist(E_pos, density=IWantDensity, bins=nbin_pos, range=range_pos, alpha=0.5, \\\n", " label='positron spectrum')\n", "h2 = plt.hist(E_pos_cut, density=IWantDensity, bins=nbin_pos, range=range_pos, alpha=0.5, \\\n", " label='positron spectrum within %.2f mrad' % (angle_cut))\n", "plt.legend(fontsize=14)\n", "plt.xlabel('Energy [MeV]')\n", "plt.ylabel('Counts [arb. units]')\n", "plt.title('')\n", "plt.grid(True)\n", "plt.yscale('log')\n", "if saveFigs:\n", " plt.savefig(outpath + 'eplusBeam_E_distrib_' + purename + '.jpg')\n", "plt.show()\n", "\n", "if not IWantDensity:\n", " print(\"sum of positron spectrum: %d\" % sum(h[0]))\n", " print(\"sum of positron spectrum within %.2f mrad: %d\\n\" % (angle_cut, sum(h2[0])))\n", "else:\n", " print(\"\\n\")\n", " \n", "# Write the positron spectrum to a file\n", "if save_spectrum:\n", " Eedges, spectrum = h[1], h[0]\n", " output_file = outpath + 'positron_spectrum_' + purename + '.txt'\n", " write_spectrum_to_G4file(Eedges, spectrum, output_file)\n", " _, spectrum_cut = h2[1], h2[0]\n", " output_file_cut = outpath + 'positron_spectrum_cut_%.1fmrad_' % angle_cut + purename + '.txt'\n", " write_spectrum_to_G4file(Eedges, spectrum_cut, output_file_cut)" ] }, { "cell_type": "markdown", "id": "92ca628d-142a-4109-b6e7-bc202f1f7cbe", "metadata": {}, "source": [ "### Convert the positron beam in RF-Track format (for tracking in the Capture system...)" ] }, { "cell_type": "code", "execution_count": null, "id": "a554a8e9-9640-4766-ae5a-6faa7e1619ef", "metadata": { "tags": [] }, "outputs": [], "source": [ "convertFileToRFTrack = False\n", "\n", "# Prepare the file to feed the RF-Track code to tack the positrons inside the Capture System/Positron Linac with\n", "\"\"\"\n", "NOTE: In RF-Track x'=px/pz=tan(thx) and y'=py/pz=tan(thy) (see the manual at page 22).\n", " The approximation x'=thx and y'=thy holds for small-angle only,\n", " which is the case of particle beams in accelators. However, \n", " for example in a drift, we correctly have x2=x1+x'*(z2-z1) and y2=y1+y'*(z2-z1).\n", " In our case, positrons can exit the crystal also at large angles \n", " and, for some of them, px and/or py can be (much) greater than pz. Thus,\n", " x' and y' are not angles, but tangents (which can tend to infinity) and \n", " the fact that they are expressed in mrad means that they are multiplied by 1e-3.\n", " These particles would be lost in the capture section. However, if we want to\n", " calculate the tranverse phase space of the positron beam at the exit of the crystal,\n", " we must calculate properly (with arctan) thx, thy and set a proper angular cut.\n", "\"\"\"\n", "if convertFileToRFTrack:\n", " ## Add variables to the original positron dataframe\n", " positrons[\"xp[mrad]\"] = (positrons.px/(positrons.pz*1000))*1000\n", " positrons[\"yp[mrad]\"] = (positrons.py/(positrons.pz*1000))*1000\n", " positrons[\"p[MeV/c]\"] = positrons.P\n", " positrons[\"#x[mm]\"] = positrons[\"x\"]*10 #this is because they were previously converted form [mm] to [cm]\n", " positrons[\"y[mm]\"] = positrons[\"y\"]*10\n", " \n", " ## Select only a set of variables\n", " if addID:\n", " positrons[\"ID\"] = list(range(1, len(positrons)+1))\n", " positrons_out = positrons[[\"#x[mm]\", \"xp[mrad]\", \"y[mm]\", \"yp[mrad]\", \"p[MeV/c]\", \"ID\"]]\n", " else:\n", " positrons_out = positrons[[\"#x[mm]\", \"xp[mrad]\", \"y[mm]\", \"yp[mrad]\", \"p[MeV/c]\"]]\n", "\n", " ## Guassian sampling in time, to keep the RF phases of the Hybrid similar to the conventional scheme\n", " c = 299792458 #m/s\n", " if setGaussianTimeShape:\n", " mean_value = 17.762 #mm/c\n", " std_value = 1.1950 #mm/c\n", " num_sample = len(positrons_out)\n", " gaussian_samples = np.random.normal(mean_value, std_value, num_sample)\n", " t = gaussian_samples\n", " else: \n", " t = positrons.t * 1e-6 * c #ns -> mm/c\n", " positrons_out.insert(4, 't[mm/c]', t)\n", "\n", " ## Save the positron phase-space to a txt file\n", " name_out = name.replace('output', 'positrons')\n", " positrons_out.to_csv(outpath + name_out + \".dat\", index=False, sep=' ') \n", " print(\"positron beam phase-space converted to RF-Track format!\")" ] }, { "cell_type": "markdown", "id": "4d7e77d2-b658-4e68-ab20-80c980242297", "metadata": {}, "source": [ "## Export the results of the simulation" ] }, { "cell_type": "code", "execution_count": null, "id": "d7aa1f40-3ca6-4c74-9350-9ee4542010dc", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Fill the result list with those of the current simulation\n", "case_list.append(purename)\n", "yield_pos_list.append(round(data_out[(data_out.particle==\"e+\")].shape[0] / Nevents, 2))\n", "yield_ele_list.append(round(data_out[(data_out.particle==\"e-\")].shape[0] / Nevents, 2))\n", "yield_ph_list.append(round(data_out[(data_out.particle==\"gamma\")].shape[0] / Nevents, 2))\n", "yield_n_list.append(round(data_out[(data_out.particle==\"neutron\")].shape[0] / Nevents, 2))\n", "pos_mean_E_list.append(round(Emean_pos, 2))\n", "pos_spread_E_list.append(round(Espread_pos, 2)) \n", "pos_mean_div_fit_list.append(round(np.mean(divergence), 2))\n", "pos_mean_size_fit_list.append(round(np.mean([np.abs(beamSizeX[1]*10), np.abs(beamSizeY[1]*10)]), 2))\n", "Edep_rad_list.append(round(np.sum(edep_rad)/Nevents, 2))\n", "Edep_conv_list.append(round(np.sum(edep_conv)*1e-3/Nevents, 2))\n", "PEDD_list.append(round(PEDD, 2))" ] }, { "cell_type": "code", "execution_count": null, "id": "f743be58-cfec-4993-bde2-c224f5bd9094", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Fill a dictionary with the results\n", "results = { \"case\": case_list,\n", " \"yield_e+\": yield_pos_list,\n", " \"yield_e-\": yield_ele_list,\n", " \"yield_ph\": yield_ph_list,\n", " \"yield_n\": yield_n_list,\n", " \"e+_mean_E[MeV]\": pos_mean_E_list,\n", " \"e+_spread_E[sigma/mu]\": pos_spread_E_list,\n", " \"e+_mean_div_fit[mrad]\": pos_mean_div_fit_list,\n", " \"e+_mean_size_fit[mm]\": pos_mean_size_fit_list,\n", " \"Edep_rad[MeV/e-]\": Edep_rad_list,\n", " \"Edep_conv[GeV/e-]\": Edep_conv_list,\n", " \"PEDD[MeV/(mm^3*e-)]\": PEDD_list\n", " }\n", "\n", "# Get a dataframe of the results\n", "df_results = pd.DataFrame.from_dict(results)\n", "\n", "# Export the results in a json file.\n", "# pickle works pretty much the same, \n", "# but it I cannot directly have a look at it.\n", "import json \n", "print(\"\\nWriting json file with these data:\")\n", "#for s in results:\n", "# print(\"\\t\", s, \":\", results[s])\n", "outname = outpath + name.replace(\"output\", \"results\") \n", "with open(outname + \".json\", \"w\") as file: \n", " json.dump(results, file)\n", "file.close()\n", "\n", "# Print results\n", "df_results" ] }, { "cell_type": "markdown", "id": "e4e92d8d-24ee-44d3-b8b6-1e6ddc382119", "metadata": {}, "source": [ "## Read scoring_ntuple2: features of the particles that exit (towards any direction) the radiator crystal" ] }, { "cell_type": "code", "execution_count": null, "id": "95150758-427a-4b05-b566-4fa9e282520b", "metadata": { "tags": [] }, "outputs": [], "source": [ "analyze_scoring_ntuple2 = True\n", "\n", "if analyze_scoring_ntuple2:\n", " df2 = rf['scoring_ntuple2'].arrays(library='pd')\n", " df2['p'] = np.sqrt(df2['px']**2+df2['py']**2+df2['pz']**2)\n", " #df2.describe()" ] }, { "cell_type": "code", "execution_count": null, "id": "f24917dd-72b6-49e8-8f50-95c15d645f8e", "metadata": { "tags": [] }, "outputs": [], "source": [ "part_sel = 'gamma'\n", "p_low = 0 #MeV/c\n", "p_high = 100 #MeV/c\n", "\n", "if analyze_scoring_ntuple2:\n", " if len(df2) > 0:\n", " df2_sel = df2[(df2.particle == part_sel)].copy()\n", " df2_sel = df2_sel[(df2_sel.p >= p_low) & (df2_sel.p <= p_high)].copy()\n", " #df2_sel.describe()" ] }, { "cell_type": "code", "execution_count": null, "id": "94034573-57ae-479a-a332-20fcf908b900", "metadata": {}, "outputs": [], "source": [ "if analyze_scoring_ntuple2 and len(df2) > 0:\n", " %matplotlib ipympl\n", "\n", " plt.figure(figsize=(10, 10))\n", " Npoint = 10000\n", " fs = 18\n", "\n", " ax3D = plt.axes(projection=\"3d\")\n", " scatt = ax3D.scatter3D(df2_sel.z[:Npoint], df_sel2.x[:Npoint], df2_sel.y[:Npoint], c=df2_sel.p[:Npoint], \\\n", " alpha=0.7, marker='.', cmap=plt.get_cmap('hsv'))\n", " cbar = fig.colorbar(scatt, ax=ax3D, shrink=0.5, aspect=5)\n", " cbar.set_label(\"p (MeV/c)\", rotation=90, labelpad=16, fontsize=fs)\n", " ax3D.set_xlabel(\"z (mm)\", fontsize=fs)\n", " ax3D.set_ylabel(\"x (mm)\", fontsize=fs)\n", " ax3D.set_zlabel(\"y (mm)\", fontsize=fs)\n", " ax3D.xaxis.set_inverted(True)\n", "\n", " elev = None\n", " azim = 120\n", " roll = None\n", " #ax3D.view_init(elev, azim, roll)\n", "\n", " if saveFigs:\n", " plt.savefig(outpath + '3Ddistrib_' + purename + '.jpg')\n", " plt.show() " ] }, { "cell_type": "code", "execution_count": null, "id": "3c22a840-2fb2-4a95-aafe-dea4a9dce6d7", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }