Warning, /geant4/examples/extended/exoticphysics/channeling/ch5/analysis_ch5.ipynb is written in an unsupported language. File is not indexed.
0001 {
0002 "cells": [
0003 {
0004 "cell_type": "markdown",
0005 "id": "b1fed382-d8eb-4a41-a658-eb9e02ae77e6",
0006 "metadata": {},
0007 "source": [
0008 "# Read the output of a Geant4 **ch5** simulation"
0009 ]
0010 },
0011 {
0012 "cell_type": "markdown",
0013 "id": "29ffa6b6-2469-45d8-bce5-94e8b8d184cf",
0014 "metadata": {},
0015 "source": [
0016 "## Import the required libraries"
0017 ]
0018 },
0019 {
0020 "cell_type": "code",
0021 "execution_count": null,
0022 "id": "3764e2c2-8cb7-4b89-9d5c-1e1379a1207d",
0023 "metadata": {
0024 "tags": []
0025 },
0026 "outputs": [],
0027 "source": [
0028 "import numpy as np\n",
0029 "import pandas as pd\n",
0030 "import matplotlib.pyplot as plt\n",
0031 "#%matplotlib ipympl\n",
0032 "from matplotlib.ticker import AutoMinorLocator\n",
0033 "import os\n",
0034 "import uproot\n",
0035 "import math\n",
0036 "import re\n",
0037 "from scipy import interpolate\n",
0038 "\n",
0039 "# Set the number of digits to show in pandas dataframes\n",
0040 "pd.set_option('display.float_format', '{:.2f}'.format)\n",
0041 "\n",
0042 "import warnings\n",
0043 "warnings.simplefilter(\"ignore\", category=DeprecationWarning)\n",
0044 "pd.options.mode.chained_assignment = None"
0045 ]
0046 },
0047 {
0048 "cell_type": "markdown",
0049 "id": "cf2a5937-f6d1-4c85-ba11-9889be23b8cd",
0050 "metadata": {},
0051 "source": [
0052 "## Graphical settings for fancier plots"
0053 ]
0054 },
0055 {
0056 "cell_type": "code",
0057 "execution_count": null,
0058 "id": "9044a604-fb1b-4429-9b54-d399919a29b0",
0059 "metadata": {},
0060 "outputs": [],
0061 "source": [
0062 "import mplhep as hep\n",
0063 "\n",
0064 "style = hep.style.ROOT\n",
0065 "\n",
0066 "# figure size\n",
0067 "style['figure.figsize'] = (10.795, 7.000)\n",
0068 "\n",
0069 "# other sizes\n",
0070 "style['lines.markersize'] = 6\n",
0071 "style['lines.linewidth'] = 3\n",
0072 "\n",
0073 "# fonts\n",
0074 "style['font.family'] = \"serif\"\n",
0075 "style['axes.titlesize'] = \"x-small\"\n",
0076 "style['axes.labelsize'] = \"small\"\n",
0077 "style['xtick.labelsize'] = \"x-small\"\n",
0078 "style['ytick.labelsize'] = \"x-small\"\n",
0079 "style['legend.fontsize'] = \"x-small\"\n",
0080 "style['mathtext.fontset'] = \"dejavuserif\"\n",
0081 "\n",
0082 "# text locations\n",
0083 "style['xaxis.labellocation'] = 'center'\n",
0084 "style['yaxis.labellocation'] = 'center'\n",
0085 "\n",
0086 "# grid\n",
0087 "#style[\"axes.grid\"] = True\n",
0088 "\n",
0089 "# output plots\n",
0090 "out_format = \"pdf\"\n",
0091 "dpi = 800"
0092 ]
0093 },
0094 {
0095 "cell_type": "markdown",
0096 "id": "b12394a6-daae-4b8e-a458-28cae054a7a4",
0097 "metadata": {
0098 "tags": []
0099 },
0100 "source": [
0101 "## Set input path and base filename"
0102 ]
0103 },
0104 {
0105 "cell_type": "code",
0106 "execution_count": null,
0107 "id": "cbebf6ec-61ad-4e71-b856-d4afdb3c18d6",
0108 "metadata": {
0109 "tags": []
0110 },
0111 "outputs": [],
0112 "source": [
0113 "path = \"output/\"\n",
0114 "\n",
0115 "name = \"output_6GeV_W2.0mm_conv8.0mm_D50cm\" #name of the output root file (without extension)\n",
0116 "rad_th = 2.0 #[mm] thickness ot the radiator\n",
0117 "conv_th = 8.0 #[mm] thickness ot the converter"
0118 ]
0119 },
0120 {
0121 "cell_type": "markdown",
0122 "id": "69c5ad7f-f62b-4c95-a1ed-144dda88cec9",
0123 "metadata": {},
0124 "source": [
0125 "## Save and export settings"
0126 ]
0127 },
0128 {
0129 "cell_type": "code",
0130 "execution_count": null,
0131 "id": "9d7a6d4a-e1c4-4a57-a388-5e66a5c6242e",
0132 "metadata": {
0133 "tags": []
0134 },
0135 "outputs": [],
0136 "source": [
0137 "# Figure save option\n",
0138 "saveFigs = True\n",
0139 "\n",
0140 "# Option to convert file for RF-Track\n",
0141 "convertFileToRFTrack = True\n",
0142 "setGaussianTimeShape = True\n",
0143 "addID = False"
0144 ]
0145 },
0146 {
0147 "cell_type": "markdown",
0148 "id": "60605cbf-1850-4d71-a3f1-d3373e464913",
0149 "metadata": {},
0150 "source": [
0151 "## Set output path and desired variables"
0152 ]
0153 },
0154 {
0155 "cell_type": "code",
0156 "execution_count": null,
0157 "id": "57131664-c6b4-42d5-8c79-a12214c51ae6",
0158 "metadata": {
0159 "tags": []
0160 },
0161 "outputs": [],
0162 "source": [
0163 "# Stripped name\n",
0164 "purename = name.replace(\"output_\", \"\")\n",
0165 "\n",
0166 "# Define the output folder\n",
0167 "outpath = path + \"analysis_output/\"\n",
0168 "if not os.path.exists(outpath):\n",
0169 " os.makedirs(outpath)\n",
0170 "\n",
0171 "# Quantity to be stored as results\n",
0172 "case_list = []\n",
0173 "yield_pos_list = []\n",
0174 "yield_ele_list = []\n",
0175 "yield_ph_list = []\n",
0176 "yield_n_list = []\n",
0177 "pos_mean_E_list = []\n",
0178 "pos_spread_E_list = []\n",
0179 "pos_mean_div_fit_list = []\n",
0180 "pos_mean_size_fit_list = []\n",
0181 "Edep_rad_list = []\n",
0182 "Edep_conv_list = []\n",
0183 "PEDD_list = []"
0184 ]
0185 },
0186 {
0187 "cell_type": "markdown",
0188 "id": "8c779b26-1bc9-49c3-a7b9-7d1289460e12",
0189 "metadata": {
0190 "tags": []
0191 },
0192 "source": [
0193 "## Read the root file with the scoring ntuples for virtual detectors"
0194 ]
0195 },
0196 {
0197 "cell_type": "code",
0198 "execution_count": null,
0199 "id": "78e5a25f-0a7c-4af5-ad45-91eca1af2231",
0200 "metadata": {
0201 "tags": []
0202 },
0203 "outputs": [],
0204 "source": [
0205 "print(\"\\nopening rf:\", name ,\".root file ...\")\n",
0206 "\n",
0207 "rf = uproot.open(path + name + \".root\")\n",
0208 "rf_content = [item.split(';')[0] for item in rf.keys()]\n",
0209 "print('rf_content:', rf_content)\n",
0210 "\n",
0211 "#screenID: 2=after the Radiator (if used), 0=before the Converter, 1=after the Converter (or the single-crystal)\n",
0212 "#particle: Geant4 particle name\n",
0213 "#x,y are in [mm]\n",
0214 "#px,py,pz are in [MeV/c]\n",
0215 "#t is in [ns]\n",
0216 "#branches = [\"screenID\", \"particle\", \"x\", \"y\", \"px\", \"py\", \"pz\", \"t\", \"eventID\"]\n",
0217 "df = rf['scoring_ntuple'].arrays(library='pd')\n",
0218 "\n",
0219 "if 'eventID' in df:\n",
0220 " Nevents = len(np.array(df.eventID.unique()))\n",
0221 "else:\n",
0222 " Nevents = 1e4\n",
0223 "print(\"Nevents:\", Nevents) #number of simulated events (primary particles)\n",
0224 "\n",
0225 "df_rad_out = df[(df.screenID == 2)].copy().drop([\"screenID\"], axis=1)\n",
0226 "df_conv_in = df[(df.screenID == 0)].copy().drop([\"screenID\"], axis=1)\n",
0227 "df_conv_out = df[(df.screenID == 1)].copy().drop([\"screenID\"], axis=1)\n",
0228 "#df_wrong = df[(df.screenID == -1)].copy().drop([\"screenID\"], axis=1)\n",
0229 "del df"
0230 ]
0231 },
0232 {
0233 "cell_type": "markdown",
0234 "id": "9094f5fd-8bd2-4501-9757-561bfbe73130",
0235 "metadata": {},
0236 "source": [
0237 "### Dataframes with phase-spaces (of $\\gamma$/e-/e+/n) after the Radiator crystal\n",
0238 "note: if D = 0 -> this is empty!"
0239 ]
0240 },
0241 {
0242 "cell_type": "code",
0243 "execution_count": null,
0244 "id": "1043904e-ff55-4920-8349-099e6cff9a5c",
0245 "metadata": {
0246 "tags": []
0247 },
0248 "outputs": [],
0249 "source": [
0250 "print(np.array(df_rad_out.particle.unique()))\n",
0251 "print(df_rad_out.shape[0], '\\n')\n",
0252 "\n",
0253 "if df_rad_out.shape[0] > 0 :\n",
0254 " df_rad_out_n = df_rad_out[(df_rad_out.particle == \"neutron\")].copy()\n",
0255 " df_rad_out_ph = df_rad_out[(df_rad_out.particle == \"gamma\")].copy()\n",
0256 " df_rad_out_ele = df_rad_out[(df_rad_out.particle == \"e-\")].copy()\n",
0257 " df_rad_out_pos = df_rad_out[(df_rad_out.particle == \"e+\")].copy()\n",
0258 " \n",
0259 " print('neutrons:', df_rad_out_n.shape[0])\n",
0260 " print('photons:', df_rad_out_ph.shape[0])\n",
0261 " print('e-:', df_rad_out_ele.shape[0], ', e+:', df_rad_out_pos.shape[0], \\\n",
0262 " '-> e-e+-:', df_rad_out_ele.shape[0] + df_rad_out_pos.shape[0], '\\n')\n",
0263 "\n",
0264 " df_rad_out.head(5)"
0265 ]
0266 },
0267 {
0268 "cell_type": "markdown",
0269 "id": "f6a9b843-8841-4a80-b947-f458c97e989f",
0270 "metadata": {},
0271 "source": [
0272 "### Dataframes with phase-spaces (of $\\gamma$/e-/e+/n) before the Converter crystal"
0273 ]
0274 },
0275 {
0276 "cell_type": "code",
0277 "execution_count": null,
0278 "id": "d69706c0-e4d8-496d-850c-bd9e6de72cfc",
0279 "metadata": {
0280 "tags": []
0281 },
0282 "outputs": [],
0283 "source": [
0284 "print(np.array(df_conv_in.particle.unique()))\n",
0285 "print(df_conv_in.shape[0], '\\n')\n",
0286 "\n",
0287 "if df_conv_in.shape[0] > 0 :\n",
0288 " df_conv_in_n = df_conv_in[(df_conv_in.particle == \"neutron\")].copy()\n",
0289 " df_conv_in_ph = df_conv_in[(df_conv_in.particle == \"gamma\")].copy()\n",
0290 " df_conv_in_ele = df_conv_in[(df_conv_in.particle == \"e-\")].copy()\n",
0291 " df_conv_in_pos = df_conv_in[(df_conv_in.particle == \"e+\")].copy()\n",
0292 "\n",
0293 " print('neutrons:', df_conv_in_n.shape[0])\n",
0294 " print('photons:', df_conv_in_ph.shape[0])\n",
0295 " print('e-:', df_conv_in_ele.shape[0], ', e+:', df_conv_in_pos.shape[0], \\\n",
0296 " '-> e-/e+-:', df_conv_in_ele.shape[0] + df_conv_in_pos.shape[0], '\\n')\n",
0297 "\n",
0298 " df_conv_in.head(5)"
0299 ]
0300 },
0301 {
0302 "cell_type": "markdown",
0303 "id": "9dcc1d28-10fa-4f7d-9fa1-914ba08d96d4",
0304 "metadata": {},
0305 "source": [
0306 "### Dataframes with phase-spaces (of $\\gamma$/e-/e+/n) after the Converter crystal\n",
0307 "note: if the source is not hybrid -> this is empty!"
0308 ]
0309 },
0310 {
0311 "cell_type": "code",
0312 "execution_count": null,
0313 "id": "f5d9ce34-6ccf-4c4f-83e4-c2ef5cb4795d",
0314 "metadata": {
0315 "tags": []
0316 },
0317 "outputs": [],
0318 "source": [
0319 "print(np.array(df_conv_out.particle.unique()))\n",
0320 "print(df_conv_out.shape[0], '\\n')\n",
0321 "\n",
0322 "if df_conv_out.shape[0] > 0 :\n",
0323 " df_conv_out_n = df_conv_out[(df_conv_out.particle == \"neutron\")].copy()\n",
0324 " df_conv_out_ph = df_conv_out[(df_conv_out.particle == \"gamma\")].copy()\n",
0325 " df_conv_out_ele = df_conv_out[(df_conv_out.particle == \"e-\")].copy()\n",
0326 " df_conv_out_pos = df_conv_out[(df_conv_out.particle == \"e+\")].copy()\n",
0327 "\n",
0328 " print('neutrons:', df_conv_out_n.shape[0])\n",
0329 " print('photons:', df_conv_out_ph.shape[0])\n",
0330 " print('e-:', df_conv_out_ele.shape[0], ', e+:', df_conv_out_pos.shape[0], \\\n",
0331 " '-> e-e+-:', df_conv_out_ele.shape[0] + df_conv_out_pos.shape[0], '\\n')\n",
0332 "\n",
0333 " df_conv_out.head(5)"
0334 ]
0335 },
0336 {
0337 "cell_type": "markdown",
0338 "id": "d40fb008-776e-48d4-bcba-39b36ae7af58",
0339 "metadata": {},
0340 "source": [
0341 "## Plot the positron distribution at the Converter/Target exit [and correlate it with the photon distribution at the exit of the radiator crystal]"
0342 ]
0343 },
0344 {
0345 "cell_type": "markdown",
0346 "id": "3190fcce-76b8-48e5-92d1-93db9eec9fd1",
0347 "metadata": {},
0348 "source": [
0349 "### Import the required libraries"
0350 ]
0351 },
0352 {
0353 "cell_type": "code",
0354 "execution_count": null,
0355 "id": "7f1578f2-2913-4dc4-8c1a-7556fd84faeb",
0356 "metadata": {
0357 "tags": []
0358 },
0359 "outputs": [],
0360 "source": [
0361 "from scipy.optimize import curve_fit\n",
0362 "from matplotlib.colors import LogNorm\n",
0363 "\n",
0364 "def gaussian(x, a, b, c):\n",
0365 " import numpy as np\n",
0366 " return a*np.exp(-np.power(x - b, 2)/(2*np.power(c, 2)))\n",
0367 "\n",
0368 "### Set plot style ###\n",
0369 "style['figure.figsize'] = (style['figure.figsize'][0], 8.50) #manual tweak to figure height\n",
0370 "plt.style.use(style)"
0371 ]
0372 },
0373 {
0374 "cell_type": "markdown",
0375 "id": "ac6e8cb8-9888-4721-8639-71496b8c9208",
0376 "metadata": {},
0377 "source": [
0378 "### Options for this calculation"
0379 ]
0380 },
0381 {
0382 "cell_type": "code",
0383 "execution_count": null,
0384 "id": "6ad51514-f780-44f4-b72b-7f7f18fa1278",
0385 "metadata": {
0386 "tags": []
0387 },
0388 "outputs": [],
0389 "source": [
0390 "if 'conventional' in name:\n",
0391 " inputType = \"e-\" \n",
0392 "else:\n",
0393 " inputType = \"all\" #set a value different from \"all\" only if it is a conventional source\n",
0394 "\n",
0395 "bCutOnMomentum = False\n",
0396 "MomentumCut = 5 #MeV"
0397 ]
0398 },
0399 {
0400 "cell_type": "markdown",
0401 "id": "7ae3c2b7-cbfb-4183-953c-ecff0d9d8db0",
0402 "metadata": {},
0403 "source": [
0404 "### Define data_in and data_out dataframes and add more parameters to them"
0405 ]
0406 },
0407 {
0408 "cell_type": "code",
0409 "execution_count": null,
0410 "id": "f1b026f9-ad0c-4d2e-a0e5-df82c318f5bf",
0411 "metadata": {
0412 "tags": []
0413 },
0414 "outputs": [],
0415 "source": [
0416 "if df_rad_out.shape[0] > 0:\n",
0417 " data_in = df_rad_out.copy()\n",
0418 "else:\n",
0419 " if df_conv_out.shape[0] > 0:\n",
0420 " data_in = df_conv_in.copy()\n",
0421 " else:\n",
0422 " data_in = df_rad_out.copy() #set an empty dataframe, because this is a single-volume (crystal or conventional) source\n",
0423 " \n",
0424 "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",
0425 "data_in = data_in[data_in.pz >= 0] #selecting only events (that should be) from the input beam\n",
0426 "if bCutOnMomentum:\n",
0427 " data_in = data_in[data_in.P >= MomentumCut] #selecting only events with momentum > MomentumCut\n",
0428 "data_in.pz = data_in.pz/1000 #MeV -> GeV\n",
0429 "data_in.x = data_in.x/10 #mm -> cm\n",
0430 "data_in.y = data_in.y/10 #mm -> cm\n",
0431 "data_in[\"pt\"] = (data_in.px**2 + data_in.py**2)**0.5 #MeV\n",
0432 "data_in[\"thx\"] = np.arctan(data_in.px/1000 / data_in.pz)*1000 #mrad\n",
0433 "data_in[\"thy\"] = np.arctan(data_in.py/1000 / data_in.pz)*1000 #mrad\n",
0434 "data_in[\"tht\"] = np.arctan(data_in.pt/1000 / data_in.pz)*1000 #mrad\n",
0435 "\n",
0436 "if df_conv_out.shape[0] > 0:\n",
0437 " data_out = df_conv_out.copy() \n",
0438 "else: #this is the case of a single-volume (crystalline or conventional) source\n",
0439 " data_out = df_conv_in.copy() \n",
0440 "\n",
0441 "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",
0442 "if bCutOnMomentum:\n",
0443 " data_out = data_out[data_out.P >= MomentumCut] #selecting only events with momentum > MomentumCut\n",
0444 "data_out.pz = data_out.pz/1000 #MeV -> GeV\n",
0445 "data_out.x = data_out.x/10 #mm -> cm\n",
0446 "data_out.y = data_out.y/10 #mm -> cm\n",
0447 "data_out[\"pt\"] = (data_out.px**2 + data_out.py**2)**0.5 #MeV\n",
0448 "data_out[\"thx\"] = np.arctan(data_out.px/1000 / data_out.pz)*1000 #mrad\n",
0449 "data_out[\"thy\"] = np.arctan(data_out.py/1000 / data_out.pz)*1000 #mrad\n",
0450 "data_out[\"tht\"] = np.arctan(data_out.pt/1000 / data_out.pz)*1000 #mrad"
0451 ]
0452 },
0453 {
0454 "cell_type": "markdown",
0455 "id": "89f362e5-7855-4bdf-8edd-eea48d2c582b",
0456 "metadata": {},
0457 "source": [
0458 "### Plot the angular distribution"
0459 ]
0460 },
0461 {
0462 "cell_type": "code",
0463 "execution_count": null,
0464 "id": "af9627fe-7a0b-431b-a7a3-0fdb2bc67f08",
0465 "metadata": {
0466 "tags": []
0467 },
0468 "outputs": [],
0469 "source": [
0470 "# plot options\n",
0471 "fitrange_ang = 17.45 #mrad (for angular distribution) \n",
0472 "ang_plot_lim = (-50, 150) #mrad\n",
0473 "fitrange = 0.2 #cm (for spatial distribution)\n",
0474 "nBins = 200 #(for spatial distribution)\n",
0475 "xRange = (-0.45, 0.45) #cm\n",
0476 "xRange1 = (-10.5, 10.5) #cm\n",
0477 "yUpperLim = 14 #(for profiles)"
0478 ]
0479 },
0480 {
0481 "cell_type": "code",
0482 "execution_count": null,
0483 "id": "52dc1aa8-62b6-472d-bf28-dff2063a00df",
0484 "metadata": {
0485 "tags": []
0486 },
0487 "outputs": [],
0488 "source": [
0489 "fig, ax = plt.subplots(nrows=2, sharex=True)\n",
0490 "\n",
0491 "divergenceIn = [0, 0]\n",
0492 "angleStdIn = [0, 0]\n",
0493 "divergence = [0, 0]\n",
0494 "angleStd = [0, 0]\n",
0495 "\n",
0496 "if len(data_in) > 0:\n",
0497 " series = [data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].thx, \n",
0498 " data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].thy]\n",
0499 " for i in (0, 1):\n",
0500 " hist = ax[i].hist(series[i], bins=15000, range=(-200, 200), histtype=\"step\", color=\"C0\")\n",
0501 " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n",
0502 " x, y = x[(x>-fitrange_ang) & (x<fitrange_ang)], y[(x>-fitrange_ang) & (x<fitrange_ang)]\n",
0503 " try:\n",
0504 " par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
0505 " except:\n",
0506 " par = (1, 0, np.inf)\n",
0507 " xplot = np.linspace(-fitrange_ang, fitrange_ang, 200)\n",
0508 " xplot2 = np.linspace(-200, -fitrange_ang, 2000)\n",
0509 " xplot3 = np.linspace(fitrange_ang, 200, 2000)\n",
0510 " label = (\"at crystal output\"+r\" ($\\it{\\gamma}$ only)\" if inputType==\"all\" else \\\n",
0511 " \"at amorphous target input\"+r\" ($\\it{e}^-$ only)\")+\", Gaussian\\n\"+r\"fit (in $\\pm$ %.2f mrad)\" % \\\n",
0512 " fitrange_ang+r\" $\\it{\\sigma}$ = %.2f mrad\" % abs(par[2])\n",
0513 " if abs(par[2]) < np.inf:\n",
0514 " ax[i].plot(xplot, gaussian(xplot, *par), c=\"C0\", alpha=0.5, label=label)\n",
0515 " ax[i].plot(xplot2, gaussian(xplot2, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
0516 " ax[i].plot(xplot3, gaussian(xplot3, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
0517 " \n",
0518 " divergenceIn[i] = abs(par[2])\n",
0519 " angleStdIn[i] = series[i].std()\n",
0520 "\n",
0521 " ax[i].set_yscale(\"linear\")\n",
0522 " ax[i].set_xlim(ang_plot_lim)\n",
0523 " ax[i].set_title(\" \"+(\"horizontal\" if i==0 else \"vertical\"), loc=\"left\")\n",
0524 " ax[i].legend(loc=\"upper right\", fontsize=\"xx-small\")\n",
0525 " ax[i].grid(True)\n",
0526 "\n",
0527 "series = [data_out[(data_out.particle==\"e+\")].thx, data_out[(data_out.particle==\"e+\")].thy]\n",
0528 "for i in (0, 1):\n",
0529 " hist = ax[i].hist(series[i], bins=800, range=(-1000, 1000), histtype=\"step\", color=\"C1\")\n",
0530 " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n",
0531 " x, y = x[(x>-fitrange_ang) & (x<fitrange_ang)], y[(x>-fitrange_ang) & (x<fitrange_ang)]\n",
0532 " try:\n",
0533 " par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
0534 " except:\n",
0535 " par = (1, 0, np.inf)\n",
0536 " xplot = np.linspace(-fitrange_ang, fitrange_ang, 200)\n",
0537 " xplot2 = np.linspace(-200, -fitrange_ang, 2000)\n",
0538 " xplot3 = np.linspace(fitrange_ang, 200, 2000)\n",
0539 " label = \"at amorphous target output\"+r\" ($\\it{e}^+$ only)\"+\", Gaussian\\n\"+r\"fit (in $\\pm$ %.2f mrad)\" % \\\n",
0540 " fitrange_ang+r\" $\\it{\\sigma}$ = %.2f mrad\" % abs(par[2])\n",
0541 " if abs(par[2]) < np.inf:\n",
0542 " ax[i].plot(xplot, gaussian(xplot, *par), c=\"C1\", alpha=0.5, label=label)\n",
0543 " ax[i].plot(xplot2, gaussian(xplot2, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
0544 " ax[i].plot(xplot3, gaussian(xplot3, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
0545 " divergence[i] = abs(par[2])\n",
0546 " angleStd[i] = series[i].std()\n",
0547 "\n",
0548 " ax[i].set_yscale(\"linear\")\n",
0549 " ax[i].set_xlim(ang_plot_lim)\n",
0550 " ax[i].set_title(\" \"+(\"horizontal\" if i==0 else \"vertical\"), loc=\"left\")\n",
0551 " ax[i].legend(loc=\"upper right\", fontsize=\"xx-small\")\n",
0552 " ax[i].grid(True)\n",
0553 "\n",
0554 "fig.supylabel(\"Counts\", fontsize=\"small\")\n",
0555 "ax[1].set_xlabel(r\"$\\it{\\theta}$ [mrad]\")\n",
0556 "fig.tight_layout()\n",
0557 "\n",
0558 "if saveFigs:\n",
0559 " plt.savefig(outpath + 'eBeam_theta_distrib_' + purename + '.jpg')\n",
0560 "\n",
0561 "print(\"e+ beam divergence along X and Y from Gauss fit in (-%.2f, %.2f) mrad range: [-%.2f, %.2f] mrad\" \\\n",
0562 " % (fitrange_ang, fitrange_ang, divergence[0], divergence[1]))"
0563 ]
0564 },
0565 {
0566 "cell_type": "markdown",
0567 "id": "eeeef894-b2ab-4f60-90af-68faa33e6562",
0568 "metadata": {},
0569 "source": [
0570 "### Plot the spatial distribution"
0571 ]
0572 },
0573 {
0574 "cell_type": "code",
0575 "execution_count": null,
0576 "id": "df9d997c-d378-40f2-83e8-558c5cf504b7",
0577 "metadata": {
0578 "tags": []
0579 },
0580 "outputs": [],
0581 "source": [
0582 "fig, ax = plt.subplots(ncols=2, nrows=3, figsize=(style['figure.figsize'][0], 14.0))\n",
0583 "fig.subplots_adjust(wspace=-1)\n",
0584 "\n",
0585 "beamSizeX, beamSizeY = [0, 0, 0], [0, 0, 0]\n",
0586 "\n",
0587 "ax[0,0].set_title(r\"$\\it{e}^{+}$ only\")\n",
0588 "ax[0,1].set_title(r\"$\\it{e}^{+}$ only\")\n",
0589 "\n",
0590 "ax[0,1].hist2d(data_out[(data_out.particle==\"e+\")].x, data_out[(data_out.particle==\"e+\")].y, \\\n",
0591 " bins=nBins, range=(xRange, xRange), norm=LogNorm())\n",
0592 "ax[0,1].grid(True)\n",
0593 "ax[0,0].hist2d(data_out[(data_out.particle==\"e+\")].x, data_out[(data_out.particle==\"e+\")].y, \\\n",
0594 " bins=nBins, range=(xRange1, xRange1), norm=LogNorm())\n",
0595 "ax[0,0].grid(True)\n",
0596 "\n",
0597 "gs = ax[1,0].get_gridspec()\n",
0598 "for ax0 in ax[1,:]:\n",
0599 " ax0.remove()\n",
0600 "axdown = fig.add_subplot(gs[1,:])\n",
0601 "\n",
0602 "if len(data_in) > 0:\n",
0603 " hist = axdown.hist(data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].x, \\\n",
0604 " bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C0\") ;\n",
0605 " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n",
0606 " x, y = x[(x>-fitrange) & (x<fitrange)], y[(x>-fitrange) & (x<fitrange)]\n",
0607 " try:\n",
0608 " par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
0609 " except:\n",
0610 " par = (1, 0, np.inf)\n",
0611 " xplot = np.linspace(-fitrange, fitrange, 200)\n",
0612 " xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
0613 " xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
0614 " label = (\"at crystal output\"+r\" ($\\it{\\gamma}$ only)\" if inputType==\"all\" else \\\n",
0615 " \"at amorphous target input\"+r\" ($\\it{e}^-$ only)\")+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
0616 " if abs(par[2]) < np.inf:\n",
0617 " axdown.plot(xplot, gaussian(xplot, *par), c=\"C0\", alpha=0.5, label=label)\n",
0618 " axdown.plot(xplot2, gaussian(xplot2, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
0619 " axdown.plot(xplot3, gaussian(xplot3, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
0620 " beamSizeX[0] = par[2]\n",
0621 "\n",
0622 "hist = axdown.hist(data_out[(data_out.particle==\"e+\")].x, bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C1\")\n",
0623 "x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n",
0624 "x, y = x[(x>-fitrange) & (x<fitrange)], y[(x>-fitrange) & (x<fitrange)]\n",
0625 "try:\n",
0626 " par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
0627 "except:\n",
0628 " par = (1, 0, np.inf)\n",
0629 "xplot = np.linspace(-fitrange, fitrange, 200)\n",
0630 "xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
0631 "xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
0632 "label = \"at amorphous target output\"+r\" ($\\it{e}^{+}$ only)\"+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
0633 "if abs(par[2]) < np.inf:\n",
0634 " axdown.plot(xplot, gaussian(xplot, *par), c=\"C1\", alpha=0.5, label=label)\n",
0635 " axdown.plot(xplot2, gaussian(xplot2, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
0636 " axdown.plot(xplot3, gaussian(xplot3, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
0637 "beamSizeX[1] = par[2]\n",
0638 "\n",
0639 "gs = ax[2,0].get_gridspec()\n",
0640 "for ax0 in ax[2,:]:\n",
0641 " ax0.remove()\n",
0642 "axdown2 = fig.add_subplot(gs[2,:])\n",
0643 "\n",
0644 "if len(data_in) > 0:\n",
0645 " hist = axdown2.hist(data_in[(data_in.particle==\"gamma\" if inputType==\"all\" else data_in.particle==\"e-\")].y, \\\n",
0646 " bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C0\") ;\n",
0647 " x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n",
0648 " x, y = x[(x>-fitrange) & (x<fitrange)], y[(x>-fitrange) & (x<fitrange)]\n",
0649 " try:\n",
0650 " par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
0651 " except:\n",
0652 " par = (1, 0, np.inf)\n",
0653 " xplot = np.linspace(-fitrange, fitrange, 200)\n",
0654 " xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
0655 " xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
0656 " label = (\"at crystal output\"+r\" ($\\it{\\gamma}$ only)\" if inputType==\"all\" else \\\n",
0657 " \"at amorphous target input\"+r\" ($\\it{e}^-$ only)\")+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
0658 " if abs(par[2]) < np.inf:\n",
0659 " axdown2.plot(xplot, gaussian(xplot, *par), c=\"C0\", alpha=0.5, label=label)\n",
0660 " axdown2.plot(xplot2, gaussian(xplot2, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
0661 " axdown2.plot(xplot3, gaussian(xplot3, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
0662 " beamSizeY[0] = par[2]\n",
0663 "\n",
0664 "hist = axdown2.hist(data_out[(data_out.particle==\"e+\")].y, bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C1\")\n",
0665 "x, y = hist[1][:-1]+0.5*abs(hist[1][1]-hist[1][0]), hist[0]\n",
0666 "x, y = x[(x>-fitrange) & (x<fitrange)], y[(x>-fitrange) & (x<fitrange)]\n",
0667 "try:\n",
0668 " par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
0669 "except:\n",
0670 " par = (1, 0, np.inf)\n",
0671 "xplot = np.linspace(-fitrange, fitrange, 200)\n",
0672 "xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
0673 "xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
0674 "label = \"at amorphous target output\"+r\" ($\\it{e}^{+}$ only)\"+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
0675 "if abs(par[2]) < np.inf:\n",
0676 " axdown2.plot(xplot, gaussian(xplot, *par), c=\"C1\", alpha=0.5, label=label)\n",
0677 " axdown2.plot(xplot2, gaussian(xplot2, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
0678 " axdown2.plot(xplot3, gaussian(xplot3, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
0679 "beamSizeY[1] = par[2]\n",
0680 "\n",
0681 "ax[0,0].set_ylabel(r\"$\\it{y}$ at amorphous\"+\"\\ntarget output [cm]\")\n",
0682 "ax[0,0].set_xlabel(r\" $\\it{x}$ at amorphous target output [cm]\", loc=\"left\")\n",
0683 "axdown.set_xlabel(r\"$\\it{x}$ [cm]\")\n",
0684 "axdown.set_ylabel(r\"$\\frac{\\mathrm{d}\\it{N}}{\\it{N}\\mathrm{d}\\it{x}}\\left[\\frac{1}{\\mathrm{cm}}\\right]$\")\n",
0685 "axdown.set_xlim((xRange[0], xRange[1]))\n",
0686 "axdown.set_ylim((0, yUpperLim))\n",
0687 "axdown.legend(loc=\"upper left\", fontsize=16)\n",
0688 "axdown.grid(True)\n",
0689 "axdown2.set_xlabel(r\"$\\it{y}$ [cm]\")\n",
0690 "axdown2.set_ylabel(r\"$\\frac{\\mathrm{d}\\it{N}}{\\it{N}\\mathrm{d}\\it{y}}\\left[\\frac{1}{\\mathrm{cm}}\\right]$\")\n",
0691 "axdown2.set_xlim((xRange[0], xRange[1]))\n",
0692 "axdown2.set_ylim((0, yUpperLim))\n",
0693 "axdown2.legend(loc=\"upper left\", fontsize=16)\n",
0694 "axdown2.grid(True)\n",
0695 "fig.tight_layout()\n",
0696 "\n",
0697 "if saveFigs:\n",
0698 " plt.savefig(outpath + 'eBeam_spatial_distrib_' + purename + '.jpg')\n",
0699 "\n",
0700 "beamSizeXStd = [data_in.x.std(), data_out[(data_out.particle==\"e+\")].x.std(), data_out[(data_out.particle==\"e+\")].x.std()]\n",
0701 "beamSizeYStd = [data_in.y.std(), data_out[(data_out.particle==\"e+\")].y.std(), data_out[(data_out.particle==\"e+\")].y.std()]\n",
0702 "\n",
0703 "print(\"e+ beam size along X and Y from Gauss fit in (-%.2f, %.2f) mm range: [-%.2f, %.2f] mm\" % \\\n",
0704 " (fitrange, fitrange, np.abs(beamSizeX[1]*10), np.abs(beamSizeY[1]*10)))"
0705 ]
0706 },
0707 {
0708 "cell_type": "markdown",
0709 "id": "f33236b5-250d-4560-83b0-c9fc7574704f",
0710 "metadata": {},
0711 "source": [
0712 "## Plot the photon spectrum emitted by the Radiator crystal "
0713 ]
0714 },
0715 {
0716 "cell_type": "code",
0717 "execution_count": null,
0718 "id": "30aabcff-7a56-41af-a3b4-1990b9f91dfa",
0719 "metadata": {
0720 "tags": []
0721 },
0722 "outputs": [],
0723 "source": [
0724 "if len(data_in) > 0:\n",
0725 " df_ph_rad = data_in[(data_in.particle==\"gamma\")]\n",
0726 " nbin_gamma = 100\n",
0727 " range_gamma = (50, 500) #(50, 2450)\n",
0728 " IWantDensityInSpectrum = False\n",
0729 " IWantFluence = False\n",
0730 " df_ph_rad_sel = df_ph_rad[(df_ph_rad.P >= range_gamma[0]) & (df_ph_rad.P <= range_gamma[1])]\n",
0731 " ph_beam_area = 4*np.pi*np.std(df_ph_rad.x)*np.std(df_ph_rad.y)*100 \n",
0732 " spectrum_Eph_rad, edges_Eph_rad = np.histogram(df_ph_rad.P, density=IWantDensityInSpectrum, \\\n",
0733 " bins=nbin_gamma, range=range_gamma)\n",
0734 " bin_Eph_rad = edges_Eph_rad[:-1] + (edges_Eph_rad[1]-edges_Eph_rad[0])*0.5\n",
0735 " if not IWantDensityInSpectrum:\n",
0736 " spectrum_Eph_rad = spectrum_Eph_rad / Nevents\n",
0737 " if IWantFluence:\n",
0738 " spectrum_Eph_rad = spectrum_Eph_rad / ph_beam_area\n",
0739 " ylbl = 'ph/mm$^2/e$'\n",
0740 " else:\n",
0741 " ylbl = 'ph/e'\n",
0742 " else:\n",
0743 " ylbl = '1/N$\\\\times$dN/dE'\n",
0744 " spectral_int_Eph_rad = spectrum_Eph_rad * bin_Eph_rad\n",
0745 " plt.figure()\n",
0746 " plt.plot(bin_Eph_rad, spectrum_Eph_rad, alpha=0.75, label='')\n",
0747 " plt.xlabel('Energy [MeV]')\n",
0748 " plt.ylabel(ylbl)\n",
0749 " plt.title('spectrum of the photons emitted by the radiator crystal')\n",
0750 " plt.grid(True)\n",
0751 " plt.yscale('log')\n",
0752 " if saveFigs:\n",
0753 " plt.savefig(outpath + 'Eph_rad_' + purename + '.jpg')\n",
0754 " plt.show()"
0755 ]
0756 },
0757 {
0758 "cell_type": "markdown",
0759 "id": "18d93d3a-8365-4810-a21f-be495f5d264e",
0760 "metadata": {
0761 "tags": []
0762 },
0763 "source": [
0764 "## Histograms of Edep per event in the Radiator and the Converter\n",
0765 "note: pay attention in case of read_PS...these plots make not much sense!"
0766 ]
0767 },
0768 {
0769 "cell_type": "code",
0770 "execution_count": null,
0771 "id": "2df91e29-f46d-4dac-85c0-11dcbbc05946",
0772 "metadata": {
0773 "tags": []
0774 },
0775 "outputs": [],
0776 "source": [
0777 "NbinsEdep = 50\n",
0778 "fs = 20"
0779 ]
0780 },
0781 {
0782 "cell_type": "markdown",
0783 "id": "464d80e7-1622-4c82-b436-e415c5fc63b8",
0784 "metadata": {},
0785 "source": [
0786 "### Radiator"
0787 ]
0788 },
0789 {
0790 "cell_type": "code",
0791 "execution_count": null,
0792 "id": "7566fd8e-98a2-497c-8bd6-a659d9321c2e",
0793 "metadata": {
0794 "tags": []
0795 },
0796 "outputs": [],
0797 "source": [
0798 "if 'edep_rad' in rf_content: \n",
0799 " df_edep_rad = rf['edep_rad'].arrays(library='pd')\n",
0800 " edep_rad = np.array(df_edep_rad[\"edep\"])\n",
0801 " \n",
0802 " print(\"total Edep in the Radiator per event:\", round(np.sum(edep_rad)/Nevents, 2), \"MeV/e-\\n\")\n",
0803 " \n",
0804 " if len(df_edep_rad) > 0:\n",
0805 "\n",
0806 " edep_rad_cut = edep_rad[edep_rad > 1]\n",
0807 "\n",
0808 " plt.figure(figsize=(8, 6))\n",
0809 " h_rad = plt.hist(edep_rad_cut, bins=NbinsEdep,\n",
0810 " histtype=\"step\", label=\"\", density=False,\n",
0811 " edgecolor=\"blue\", linewidth=1.2, alpha=1.)\n",
0812 " plt.title('Edep per event in the radiator', fontsize=fs)\n",
0813 " plt.xlabel(\"Edep [MeV]\", fontsize=fs)\n",
0814 " plt.ylabel(\"Counts\", fontsize=fs)\n",
0815 " plt.xticks(fontsize=fs, rotation=0)\n",
0816 " plt.yticks(fontsize=fs, rotation=0)\n",
0817 " plt.grid(which=\"major\", color=\"gray\", linestyle=\"--\", linewidth=1)\n",
0818 " if saveFigs:\n",
0819 " plt.savefig(outpath + 'Edep_distrib_rad_' + purename + '.jpg')\n",
0820 " plt.show() \n",
0821 "else:\n",
0822 " edep_rad = []"
0823 ]
0824 },
0825 {
0826 "cell_type": "markdown",
0827 "id": "a5c0c89a-3e5c-4856-adf1-73cd6128916b",
0828 "metadata": {},
0829 "source": [
0830 "### Converter"
0831 ]
0832 },
0833 {
0834 "cell_type": "code",
0835 "execution_count": null,
0836 "id": "0b11b063-ffe4-4284-994a-54593e9b1daf",
0837 "metadata": {
0838 "tags": []
0839 },
0840 "outputs": [],
0841 "source": [
0842 "if 'edep_conv' in rf_content:\n",
0843 " df_edep_conv = rf['edep_conv'].arrays(library='pd')\n",
0844 " edep_conv = np.array(df_edep_conv[\"edep\"])\n",
0845 " \n",
0846 " print(\"total Edep in the Converter per event:\", round(np.sum(edep_conv)*1.e-3/Nevents, 2), \"GeV/e-\\n\")\n",
0847 " \n",
0848 " if len(df_edep_conv) > 0 and not \"_GT\" in name:\n",
0849 "\n",
0850 " edep_conv_cut = edep_conv[edep_conv > 0]\n",
0851 "\n",
0852 " plt.figure(figsize=(8, 6))\n",
0853 " h_conv = plt.hist(edep_conv_cut, bins=NbinsEdep,\n",
0854 " histtype=\"step\", label=\"\", density=False,\n",
0855 " edgecolor=\"blue\", linewidth=1.2, alpha=1.)\n",
0856 " plt.title('Edep per event in the converter', fontsize=fs)\n",
0857 " plt.xlabel(\"Edep [MeV]\", fontsize=fs)\n",
0858 " plt.ylabel(\"Counts\", fontsize=fs)\n",
0859 " plt.xticks(fontsize=fs, rotation=0)\n",
0860 " plt.yticks(fontsize=fs, rotation=0)\n",
0861 " plt.grid(which=\"major\", color=\"gray\", linestyle=\"--\", linewidth=1)\n",
0862 " if saveFigs:\n",
0863 " plt.savefig(outpath + 'Edep_distrib_conv_' + purename + '.jpg')\n",
0864 " plt.show()\n",
0865 "else:\n",
0866 " edep_conv = []"
0867 ]
0868 },
0869 {
0870 "cell_type": "markdown",
0871 "id": "6f646b9a-ec72-4500-adc1-8d998c6683aa",
0872 "metadata": {
0873 "tags": []
0874 },
0875 "source": [
0876 "## Settings to read the built-in BoxMesh scorer with Edep in the Absorber (Converter or Radiator)"
0877 ]
0878 },
0879 {
0880 "cell_type": "code",
0881 "execution_count": null,
0882 "id": "dcd7ef10-3916-4ce0-badc-3f140d9aff42",
0883 "metadata": {
0884 "tags": []
0885 },
0886 "outputs": [],
0887 "source": [
0888 "tranvsizeX = 100 #volume tranversal X size (mm)\n",
0889 "tranvsizeY = 100 #volume tranversal Y size (mm)\n",
0890 "\n",
0891 "print(\"tranverse sixe (X*Y) = %.1f x %.1f mm2\" % (tranvsizeX, tranvsizeY))\n",
0892 "\n",
0893 "if 'conv_th' in locals() and conv_th != '': \n",
0894 " tranvsizeZ = conv_th\n",
0895 "elif 'rad_th' in locals() and rad_th != '':\n",
0896 " tranvsizeZ = rad_th\n",
0897 "else:\n",
0898 " try:\n",
0899 " tranvsizeZ = float(name.split('W')[1].split('mm_')[0])\n",
0900 " except:\n",
0901 " print('Enter the thickness [mm] of the target (or the Radiator if it is not an hybrid source):')\n",
0902 " tranvsizeZ = float(input())\n",
0903 " \n",
0904 "print(\"tranvsizeZ = %.1f mm\\n\" % (tranvsizeZ))\n",
0905 " \n",
0906 "IWantPlotVoxel = False\n",
0907 "IWantPlotX0fract = False\n",
0908 "X0 = 3.504 #crystal radiation length (mm)\n",
0909 "\n",
0910 "normEvents = True; #normalize w.r.t. simulated events (set always to true to compare with previous results)\n",
0911 "\n",
0912 "axis_proj = 2 #0: axial, 1:sagittal, 2:coronal\n",
0913 "proj_type = 1 #1: mean, 2: sum\n",
0914 "p2m = 1 #voxels around which to mean\n",
0915 "\n",
0916 "plot_proj = True\n",
0917 "reader_verb = True\n",
0918 "\n",
0919 "titstr = ''\n",
0920 "\n",
0921 "cmap2D = 'viridis'\n",
0922 "Wfig = 12\n",
0923 "Hfig = 8\n",
0924 "fs = 20\n",
0925 "ms = 6"
0926 ]
0927 },
0928 {
0929 "cell_type": "markdown",
0930 "id": "26ce0632-da84-4982-a944-0680009b2746",
0931 "metadata": {},
0932 "source": [
0933 "### Function to read the result provided by the built-in BoxMesh scorer"
0934 ]
0935 },
0936 {
0937 "cell_type": "code",
0938 "execution_count": null,
0939 "id": "b511c1ad-a6eb-43f6-ac0e-21c88b1ca39a",
0940 "metadata": {
0941 "tags": []
0942 },
0943 "outputs": [],
0944 "source": [
0945 "def read_Edep_BoxMesh(filename, normEvents, Nevents, \n",
0946 " tranvsizeX, tranvsizeY, tranvsizeZ):\n",
0947 " \"\"\"\n",
0948 " Function to read the total Edep accumulated through a built-in BoxMesh scorer.\n",
0949 " It is general. It returns [data, (x,y,z)], where:\n",
0950 " data is a dictionary with 10 columns defined as follows\n",
0951 " {\"ind_x\", \"ind_y\", \"ind_z\", \"eDep\", \"x\", \"y\", \"z\", \"r\" ,\"eDep_err\", \"eDepDensity\", \"eDepDensity_err\"}\n",
0952 " and (x,y,z) are the coordinates [mm] of the voxel centers.\n",
0953 " \"\"\"\n",
0954 " \n",
0955 " def find(cond, N=1e7): \n",
0956 " \"\"\"\n",
0957 " function, equivalent to Matlab's one, that returns\n",
0958 " a list with the indices of an array/matrix\n",
0959 " that satisfy a given condition.\n",
0960 " \"\"\"\n",
0961 " cond = cond.reshape(np.size(cond), 1).ravel()\n",
0962 " index = list(np.argwhere(cond).ravel())\n",
0963 " return index[:N]\n",
0964 " \n",
0965 " # Retrieve data and put them into a dataframe\n",
0966 " data = pd.read_csv(\n",
0967 " filename, skiprows=3, names=[\"ind_x\", \"ind_y\", \"ind_z\", \"eDep\", \"eDep2\", \"Nentry\"]\n",
0968 " )\n",
0969 " \n",
0970 " # Calculate eDep uncertainty \n",
0971 " data[\"eDep_err\"] = (data[\"eDep2\"] - data[\"eDep\"]**2/data[\"Nentry\"])**0.5\n",
0972 " data.fillna(0, inplace=True)\n",
0973 " data = data.drop(columns=['eDep2', 'Nentry'])\n",
0974 "\n",
0975 " # Retrieve the number of voxels in each direction\n",
0976 " Nvoxel = len(data)\n",
0977 " zzzz = find(np.array(data[\"ind_z\"]) == 0, 2)\n",
0978 " if len(zzzz) > 1:\n",
0979 " nz = zzzz[1]\n",
0980 " else:\n",
0981 " nz = len(data[\"ind_y\"]) \n",
0982 " ny = np.max([round((np.max(find(np.array(data[\"ind_y\"])==0, nz+1))-1)/nz), 1])\n",
0983 " nx = max(np.array(data[\"ind_x\"]))+1\n",
0984 " \n",
0985 " # Set geometrical size\n",
0986 " Xmin = -tranvsizeX*0.5 \n",
0987 " Xmax = tranvsizeX*0.5 \n",
0988 " Ymin = -tranvsizeY*0.5 \n",
0989 " Ymax = tranvsizeY*0.5\n",
0990 " Zmin = 0\n",
0991 " Zmax = tranvsizeZ\n",
0992 "\n",
0993 " # Generate vectors of points linearly spaced (corrected in 24/12/2023)\n",
0994 " Xedges = np.linspace(Xmin, Xmax, nx+1) #[mm]\n",
0995 " Yedges = np.linspace(Ymin, Ymax, ny+1) #[mm]\n",
0996 " Zedges = np.linspace(Zmin, Zmax, nz+1) #[mm]\n",
0997 " dx = Xedges[1] - Xedges[0]\n",
0998 " dy = Yedges[1] - Yedges[0]\n",
0999 " dz = Zedges[1] - Zedges[0]\n",
1000 " #print(\"dx:\", dx, \"mm, dy:\", dy, \"mm, dz:\", dz, \"mm\")\n",
1001 " x = Xedges[:-1] + dx*0.5\n",
1002 " y = Yedges[:-1] + dy*0.5\n",
1003 " z = Zedges[:-1] + dz*0.5\n",
1004 "\n",
1005 " # Physical position inside the mesh data\n",
1006 " data[\"x\"] = dx * (data[\"ind_x\"] + 0.5) - tranvsizeX*0.5 #central x [mm]\n",
1007 " data[\"y\"] = dy * (data[\"ind_y\"] + 0.5) - tranvsizeY*0.5 #central y [mm]\n",
1008 " data[\"z\"] = dz * (data[\"ind_z\"] + 0.5) #central z [mm]\n",
1009 " data[\"r\"] = np.sqrt(data[\"x\"]**2 + data[\"y\"]**2) #central r [mm]\n",
1010 " \n",
1011 " # eDep Map\n",
1012 " if normEvents:\n",
1013 " data[\"eDep\"] = data[\"eDep\"] / Nevents\n",
1014 " data[\"eDep_err\"] = data[\"eDep_err\"] / Nevents\n",
1015 " data[\"eDepDensity\"] = data[\"eDep\"] / (dz * dy * dx) #[MeV/mm**3] or [MeV/(mm**3 event)]\n",
1016 " data[\"eDepDensity_err\"] = data[\"eDep_err\"] / (dz * dy * dx) #[MeV/mm**3] or [MeV/(mm**3 event)]\n",
1017 "\n",
1018 " return data, (x,y,z)"
1019 ]
1020 },
1021 {
1022 "cell_type": "markdown",
1023 "id": "1edbb447-5b72-4c6c-a5d7-a92cbe25d3bd",
1024 "metadata": {},
1025 "source": [
1026 "## Read the built-in BoxMesh scorer with Edep in the Absorber"
1027 ]
1028 },
1029 {
1030 "cell_type": "code",
1031 "execution_count": null,
1032 "id": "b6f22f16-1cdc-44b2-88bb-848fac3e52b3",
1033 "metadata": {
1034 "tags": []
1035 },
1036 "outputs": [],
1037 "source": [
1038 "# read BoxMesh scorer\n",
1039 "try:\n",
1040 " filename = path + name.replace(\"output\", \"Edep\") + \".txt\"\n",
1041 "except:\n",
1042 " print('Enter the name of the file with the BoxMesh scorer results:')\n",
1043 " filename = float(input())\n",
1044 "\n",
1045 "results_BoxMesh = read_Edep_BoxMesh(filename, normEvents, Nevents, \n",
1046 " tranvsizeX, tranvsizeY, tranvsizeZ)\n",
1047 "df_BoxMesh = results_BoxMesh[0]\n",
1048 "voxel_coord = results_BoxMesh[1]\n",
1049 "del results_BoxMesh\n",
1050 "\n",
1051 "# retrieve variables form BoxMesh scorer\n",
1052 "eDep = df_BoxMesh.eDep\n",
1053 "eDepDensity = df_BoxMesh.eDepDensity\n",
1054 "x = voxel_coord[0]\n",
1055 "y = voxel_coord[1]\n",
1056 "z = voxel_coord[2]\n",
1057 "if x.shape[0] > 1:\n",
1058 " dx = x[1] - x[0]\n",
1059 "else:\n",
1060 " dx = 0\n",
1061 "if y.shape[0] > 1:\n",
1062 " dy = y[1] - y[0]\n",
1063 "else:\n",
1064 " dy = 0\n",
1065 "if z.shape[0] > 1:\n",
1066 " dz = z[1] - z[0]\n",
1067 "else:\n",
1068 " dz = 0\n",
1069 "\n",
1070 "# variables for contour plot of Edep\n",
1071 "nx, ny , nz = [voxel_coord[i].shape[0] for i in range(3)]\n",
1072 "X, Z = np.meshgrid(x, z)\n",
1073 "icx = 0 if nx == 1 else int(X.shape[1]/2+1)\n",
1074 "eDepMap = np.transpose(np.array(eDep).reshape(nx,ny,nz)[icx,:,:])\n",
1075 "eDepDensityMap = np.transpose(np.array(eDepDensity).reshape(nx,ny,nz)[icx,:,:])\n",
1076 "\n",
1077 "# Zprofile of Edep\n",
1078 "ic = int(Z.shape[1]/2+1)\n",
1079 "profileZ_eDep = np.zeros(nz)\n",
1080 "profileZ_eDepDensity = np.zeros(nz)\n",
1081 "for j in range(nz):\n",
1082 " if proj_type == 1:\n",
1083 " profileZ_eDep[j] = np.mean(eDepMap[j][ic-p2m:ic+p2m])\n",
1084 " profileZ_eDepDensity[j] = np.mean(eDepDensityMap[j][ic-p2m:ic+p2m])\n",
1085 " else:\n",
1086 " profileZ_eDep[j] = np.sum(eDepMap[j][ic-p2m:ic+p2m])\n",
1087 " profileZ_eDepDensity[j] = np.sum(eDepDensityMap[j][ic-p2m:ic+p2m])\n",
1088 "\n",
1089 "PEDD = np.max(np.array(eDepDensity))\n",
1090 "print(\"PEDD:\", round(PEDD,2), \"MeV/(mm^3 e-)\")"
1091 ]
1092 },
1093 {
1094 "cell_type": "markdown",
1095 "id": "b4339b43-c89a-4f10-8639-f29416c5845c",
1096 "metadata": {},
1097 "source": [
1098 "### Contour and plofile plot of Edep (density) distribution"
1099 ]
1100 },
1101 {
1102 "cell_type": "code",
1103 "execution_count": null,
1104 "id": "4a89f66c-cf61-4d95-b4a4-e30769fa6f66",
1105 "metadata": {
1106 "tags": []
1107 },
1108 "outputs": [],
1109 "source": [
1110 "# set options\n",
1111 "if normEvents:\n",
1112 " clbl = \"Edep density per event [MeV/(mm$^3$ e-]\"\n",
1113 "else:\n",
1114 " clbl = \"Edep density [MeV/mm$^3$]\"\n",
1115 "\n",
1116 "# contour plot of Edep\n",
1117 "if eDepDensityMap.shape[0]>1 and eDepDensityMap.shape[1]>1: \n",
1118 " plt.figure(figsize=(Wfig, Hfig))\n",
1119 " Nlev = 50\n",
1120 " if IWantPlotVoxel:\n",
1121 " plt.contourf(eDepDensityMap, Nlev, cmap='viridis')\n",
1122 " xlabel = \"X [voxel]\"\n",
1123 " ylabel = \"Z [voxel]\" \n",
1124 " else:\n",
1125 " plt.contourf(X, Z, eDepDensityMap, Nlev, cmap='viridis')\n",
1126 " xlabel = \"X [mm]\"\n",
1127 " ylabel = \"Z [mm]\"\n",
1128 " #plt.gca().set_aspect('equal')\n",
1129 " plt.gca().invert_yaxis()\n",
1130 " cbar = plt.colorbar()\n",
1131 " cbar.ax.tick_params(labelsize=fs)\n",
1132 " cbar.set_label(clbl, fontsize=fs, rotation=90)\n",
1133 " plt.title(titstr, fontsize=fs)\n",
1134 " plt.xlabel(xlabel, fontsize=fs)\n",
1135 " plt.ylabel(ylabel, fontsize=fs, wrap=True)\n",
1136 " plt.xticks(fontsize=fs, rotation=0)\n",
1137 " plt.yticks(fontsize=fs, rotation=0)\n",
1138 " plt.gca().xaxis.set_minor_locator(AutoMinorLocator(5))\n",
1139 " plt.gca().yaxis.set_minor_locator(AutoMinorLocator(5))\n",
1140 " plt.gca().tick_params(axis=\"both\", which='major', direction='in', length=8)\n",
1141 " plt.gca().tick_params(axis=\"both\", which='minor', direction='in', length=4)\n",
1142 " plt.grid(False)\n",
1143 " if saveFigs:\n",
1144 " plt.savefig(outpath + 'Edep_proj_' + purename + '.jpg')\n",
1145 " plt.show()\n",
1146 "\n",
1147 "# plot Zprofile of Edep\n",
1148 "if IWantPlotVoxel:\n",
1149 " xplot = np.linspace(0, nz-1, nz)\n",
1150 " ylabel = \"Z [voxel]\" \n",
1151 "else:\n",
1152 " xplot = z\n",
1153 " ylabel = \"Z [mm]\"\n",
1154 "plt.figure(figsize=(Wfig, Hfig))\n",
1155 "plt.plot(xplot, profileZ_eDepDensity, label=\"\",\n",
1156 " color=\"blue\", linestyle='-', linewidth=2, alpha=0.75,\n",
1157 " marker='', markersize=ms, markerfacecolor=\"Blue\")\n",
1158 "plt.title(titstr, fontsize=fs)\n",
1159 "plt.xlabel(ylabel, fontsize=fs, wrap=True)\n",
1160 "plt.ylabel(clbl, fontsize=fs, wrap=True)\n",
1161 "plt.xticks(fontsize=fs, rotation=0)\n",
1162 "plt.yticks(fontsize=fs, rotation=0)\n",
1163 "plt.gca().xaxis.set_minor_locator(AutoMinorLocator(5))\n",
1164 "plt.gca().yaxis.set_minor_locator(AutoMinorLocator(5))\n",
1165 "plt.gca().tick_params(axis=\"both\", which='major', direction='in', length=8)\n",
1166 "plt.gca().tick_params(axis=\"both\", which='minor', direction='in', length=4)\n",
1167 "plt.grid(True)\n",
1168 "if saveFigs:\n",
1169 " plt.savefig(outpath + 'Edep_profileZ_' + purename + '.jpg')\n",
1170 "plt.show()"
1171 ]
1172 },
1173 {
1174 "cell_type": "markdown",
1175 "id": "96e92388-7a39-40b7-a1a5-25adfdab1e0c",
1176 "metadata": {},
1177 "source": [
1178 "### radius-vs-depth energy density"
1179 ]
1180 },
1181 {
1182 "cell_type": "code",
1183 "execution_count": null,
1184 "id": "4f6c664d-8140-4b16-aa7d-8cdd4a6d3d1c",
1185 "metadata": {
1186 "tags": []
1187 },
1188 "outputs": [],
1189 "source": [
1190 "# plot options\n",
1191 "rLim = 1.5 #mm\n",
1192 "bScaleWidth = False\n",
1193 "bLog = False\n",
1194 "upperLim = 20 #MeV/(mm^3 e-)\n",
1195 "ind_r_lsit = [0, 7, 15, 22] #6GeV: [0, 5, 10, 15]"
1196 ]
1197 },
1198 {
1199 "cell_type": "code",
1200 "execution_count": null,
1201 "id": "414084ec-d20e-4863-b513-d0221ad43e96",
1202 "metadata": {
1203 "tags": []
1204 },
1205 "outputs": [],
1206 "source": [
1207 "# radius-vs-depth energy density\n",
1208 "bData = df_BoxMesh[\"r\"] < rLim if rLim != None else df_BoxMesh[\"r\"] < 1e5\n",
1209 "Z, R = np.meshgrid(\n",
1210 " sorted(list(set(df_BoxMesh[bData][\"z\"]))),\n",
1211 " sorted(list(set(df_BoxMesh[bData][\"r\"]))),\n",
1212 ")\n",
1213 "ED = interpolate.griddata(\n",
1214 " (df_BoxMesh[bData][\"z\"], df_BoxMesh[bData][\"r\"]),\n",
1215 " df_BoxMesh[bData][\"eDepDensity\"], (Z, R), method='nearest'\n",
1216 ")\n",
1217 "\n",
1218 "if ED.shape[0]>1 and ED.shape[1]>1: \n",
1219 " width = 5*max(df_BoxMesh[\"z\"]) if bScaleWidth else 10\n",
1220 " fig, ax = plt.subplots(num=\"mesh\", figsize=[width, 8], nrows=2, sharex=True, constrained_layout=True)\n",
1221 "\n",
1222 " # 2d plot\n",
1223 " plotmesh = ax[0].pcolormesh(Z, R, ED, shading=\"nearest\", norm=LogNorm() if bLog else None)\n",
1224 " cbar = plt.colorbar(plotmesh, ax=ax[0])\n",
1225 " ax[0].grid()\n",
1226 " ax[0].set_ylim((0, max(df_BoxMesh[bData][\"r\"])))\n",
1227 "\n",
1228 " # 1d projection - energy density vs depth\n",
1229 " dr = (dx**2 + dy**2)**0.5\n",
1230 " for ind_r in ind_r_lsit:\n",
1231 " ind_r_eff = int(ind_r * 0.1 / dr)\n",
1232 " r_eff = dr * ind_r_eff\n",
1233 " bData = (df_BoxMesh[\"r\"] >= (dr * (ind_r_eff - 0.5))) & \\\n",
1234 " (df_BoxMesh[\"r\"] <= (dr * (ind_r_eff + 0.5))) \n",
1235 " temp = df_BoxMesh[bData][[\"z\", \"eDepDensity\"]].sort_values(by=[\"z\"])\n",
1236 " zlist = np.array(temp.z.unique())\n",
1237 " ymean = [np.mean(np.array(temp[temp[\"z\"] == zi][\"eDepDensity\"])) for zi in zlist]\n",
1238 " ax[1].plot(\n",
1239 " #df_BoxMesh[bData][\"z\"], df_BoxMesh[bData][\"eDepDensity\"],\n",
1240 " zlist, ymean, #corrected by gpaterno\n",
1241 " label = \"r = %.2f mm\" % r_eff,\n",
1242 " marker='', linestyle='-'\n",
1243 " )\n",
1244 " #ax[1].set_ylim((0, upperLim))\n",
1245 " ax[1].grid()\n",
1246 " ax[1].legend(loc=\"upper left\", fontsize=14)\n",
1247 " fs_local = 18\n",
1248 "\n",
1249 " ax[1].set_xlabel(\"z [mm]\", fontsize=fs_local)\n",
1250 " ax[0].set_ylabel(\"r [mm]\", fontsize=fs_local)\n",
1251 " ax[1].set_ylabel(r\"eDep [MeV/(mm$^3$ e-)]\", fontsize=fs_local)\n",
1252 " cbar.set_label(r\"eDep [MeV/(mm$^3$ e-)]\", rotation=90, labelpad=16, fontsize=fs_local)\n",
1253 " #fig.suptitle(\"e- beam features: %.2f GeV, %.1f mm\" % (Ein, sigma_in))\n",
1254 " if saveFigs:\n",
1255 " plt.savefig(outpath + 'Edep_profileR_' + purename + '.jpg')\n",
1256 " plt.show()"
1257 ]
1258 },
1259 {
1260 "cell_type": "markdown",
1261 "id": "4f319d8c-417b-4e96-a0f3-843939ed5c6b",
1262 "metadata": {},
1263 "source": [
1264 "## Get the positrons, analyze and export them"
1265 ]
1266 },
1267 {
1268 "cell_type": "markdown",
1269 "id": "26382273-291e-4998-ab0d-1eb30127bb50",
1270 "metadata": {},
1271 "source": [
1272 "### Get the positrons"
1273 ]
1274 },
1275 {
1276 "cell_type": "code",
1277 "execution_count": null,
1278 "id": "3a654a57-cdaf-4eb5-9660-bf3ead04891e",
1279 "metadata": {
1280 "tags": []
1281 },
1282 "outputs": [],
1283 "source": [
1284 "positrons = data_out[data_out.particle==\"e+\"].copy()\n",
1285 "print(\"yield_e+:\", round(positrons.shape[0] / Nevents, 2), \"\\n\")\n",
1286 "\n",
1287 "positrons"
1288 ]
1289 },
1290 {
1291 "cell_type": "markdown",
1292 "id": "b4e1a5ad-8e59-4d14-8d84-60789cbe9ad7",
1293 "metadata": {},
1294 "source": [
1295 "### Get positrons within a given energy (momentum) threshold"
1296 ]
1297 },
1298 {
1299 "cell_type": "code",
1300 "execution_count": null,
1301 "id": "5dc3eab2-a137-420d-a1c5-9b2a9b8c1670",
1302 "metadata": {
1303 "tags": []
1304 },
1305 "outputs": [],
1306 "source": [
1307 "p_th = 100 #MeV\n",
1308 "positrons[positrons[\"P\"] < 100]"
1309 ]
1310 },
1311 {
1312 "cell_type": "markdown",
1313 "id": "f1ae45cb-78ff-4633-ba38-e863e9e64b94",
1314 "metadata": {},
1315 "source": [
1316 "### Plot and Save the positron spectrum (taking into account angular and energy cuts)"
1317 ]
1318 },
1319 {
1320 "cell_type": "code",
1321 "execution_count": null,
1322 "id": "b455dff7-b783-4f28-b64a-c5a08aa55888",
1323 "metadata": {
1324 "tags": []
1325 },
1326 "outputs": [],
1327 "source": [
1328 "save_spectrum = False\n",
1329 "\n",
1330 "# Calculate positron beam energy distribution features\n",
1331 "E_pos = (positrons.P**2 + 0.511**2)**0.5 #MeV\n",
1332 "Emean_pos = np.mean(E_pos)\n",
1333 "Estdev = np.std(E_pos)\n",
1334 "Espread_pos = Estdev / Emean_pos\n",
1335 "print(\"Emean_pos:\", round(Emean_pos, 2), \"MeV\")\n",
1336 "print(\"Espread_pos:\", round(Espread_pos, 2), \"(sigma/mu)\")\n",
1337 "#angle_cut = np.mean(divergence) #mrad\n",
1338 "angle_cut = 140 #mrad\n",
1339 "positrons_cut = positrons[positrons[\"tht\"] < angle_cut]\n",
1340 "#mask = (positrons.x**2+positrons.y**2)**0.5 < np.abs(np.mean([beamSizeX[1], beamSizeY[1]]))\n",
1341 "#positrons_cut = positrons[mask] \n",
1342 "E_pos_cut = (positrons_cut.P**2 + 0.511**2)**0.5 #MeV\n",
1343 "Emean_pos_cut = np.mean(E_pos_cut)\n",
1344 "Estdev_cut = np.std(E_pos_cut)\n",
1345 "Espread_pos_cut = Estdev_cut / Emean_pos_cut\n",
1346 "print(\"Emean_pos_cut (tht < %.2f mrad):\" % angle_cut, round(Emean_pos_cut, 2), \"MeV\")\n",
1347 "print(\"Espread_pos_cut (tht < %.2f mrad):\" % angle_cut, round(Espread_pos_cut, 2), \"(sigma/mu)\")\n",
1348 "Eth = 100 # MeV\n",
1349 "NposEth = len(E_pos.values[E_pos.values < Eth]) / len(E_pos.values)\n",
1350 "NsigmaEth = Eth / (Emean_pos + Estdev)\n",
1351 "print(\"fraction of positrons with energy < %.0f MeV (%.2f*sigma): %.2f\\n\" % (Eth, NsigmaEth, NposEth))\n",
1352 "\n",
1353 "# Plot histogram of positrons energy\n",
1354 "nbin_pos = 100\n",
1355 "range_pos = (0, 100)\n",
1356 "IWantDensity = False\n",
1357 "plt.figure(figsize=(9, 6))\n",
1358 "h = plt.hist(E_pos, density=IWantDensity, bins=nbin_pos, range=range_pos, alpha=0.5, \\\n",
1359 " label='positron spectrum')\n",
1360 "h2 = plt.hist(E_pos_cut, density=IWantDensity, bins=nbin_pos, range=range_pos, alpha=0.5, \\\n",
1361 " label='positron spectrum within %.2f mrad' % (angle_cut))\n",
1362 "plt.legend(fontsize=14)\n",
1363 "plt.xlabel('Energy [MeV]')\n",
1364 "plt.ylabel('Counts [arb. units]')\n",
1365 "plt.title('')\n",
1366 "plt.grid(True)\n",
1367 "plt.yscale('log')\n",
1368 "if saveFigs:\n",
1369 " plt.savefig(outpath + 'eplusBeam_E_distrib_' + purename + '.jpg')\n",
1370 "plt.show()\n",
1371 "\n",
1372 "if not IWantDensity:\n",
1373 " print(\"sum of positron spectrum: %d\" % sum(h[0]))\n",
1374 " print(\"sum of positron spectrum within %.2f mrad: %d\\n\" % (angle_cut, sum(h2[0])))\n",
1375 "else:\n",
1376 " print(\"\\n\")\n",
1377 " \n",
1378 "# Write the positron spectrum to a file\n",
1379 "if save_spectrum:\n",
1380 " Eedges, spectrum = h[1], h[0]\n",
1381 " output_file = outpath + 'positron_spectrum_' + purename + '.txt'\n",
1382 " write_spectrum_to_G4file(Eedges, spectrum, output_file)\n",
1383 " _, spectrum_cut = h2[1], h2[0]\n",
1384 " output_file_cut = outpath + 'positron_spectrum_cut_%.1fmrad_' % angle_cut + purename + '.txt'\n",
1385 " write_spectrum_to_G4file(Eedges, spectrum_cut, output_file_cut)"
1386 ]
1387 },
1388 {
1389 "cell_type": "markdown",
1390 "id": "92ca628d-142a-4109-b6e7-bc202f1f7cbe",
1391 "metadata": {},
1392 "source": [
1393 "### Convert the positron beam in RF-Track format (for tracking in the Capture system...)"
1394 ]
1395 },
1396 {
1397 "cell_type": "code",
1398 "execution_count": null,
1399 "id": "a554a8e9-9640-4766-ae5a-6faa7e1619ef",
1400 "metadata": {
1401 "tags": []
1402 },
1403 "outputs": [],
1404 "source": [
1405 "convertFileToRFTrack = False\n",
1406 "\n",
1407 "# Prepare the file to feed the RF-Track code to tack the positrons inside the Capture System/Positron Linac with\n",
1408 "\"\"\"\n",
1409 "NOTE: In RF-Track x'=px/pz=tan(thx) and y'=py/pz=tan(thy) (see the manual at page 22).\n",
1410 " The approximation x'=thx and y'=thy holds for small-angle only,\n",
1411 " which is the case of particle beams in accelators. However, \n",
1412 " for example in a drift, we correctly have x2=x1+x'*(z2-z1) and y2=y1+y'*(z2-z1).\n",
1413 " In our case, positrons can exit the crystal also at large angles \n",
1414 " and, for some of them, px and/or py can be (much) greater than pz. Thus,\n",
1415 " x' and y' are not angles, but tangents (which can tend to infinity) and \n",
1416 " the fact that they are expressed in mrad means that they are multiplied by 1e-3.\n",
1417 " These particles would be lost in the capture section. However, if we want to\n",
1418 " calculate the tranverse phase space of the positron beam at the exit of the crystal,\n",
1419 " we must calculate properly (with arctan) thx, thy and set a proper angular cut.\n",
1420 "\"\"\"\n",
1421 "if convertFileToRFTrack:\n",
1422 " ## Add variables to the original positron dataframe\n",
1423 " positrons[\"xp[mrad]\"] = (positrons.px/(positrons.pz*1000))*1000\n",
1424 " positrons[\"yp[mrad]\"] = (positrons.py/(positrons.pz*1000))*1000\n",
1425 " positrons[\"p[MeV/c]\"] = positrons.P\n",
1426 " positrons[\"#x[mm]\"] = positrons[\"x\"]*10 #this is because they were previously converted form [mm] to [cm]\n",
1427 " positrons[\"y[mm]\"] = positrons[\"y\"]*10\n",
1428 " \n",
1429 " ## Select only a set of variables\n",
1430 " if addID:\n",
1431 " positrons[\"ID\"] = list(range(1, len(positrons)+1))\n",
1432 " positrons_out = positrons[[\"#x[mm]\", \"xp[mrad]\", \"y[mm]\", \"yp[mrad]\", \"p[MeV/c]\", \"ID\"]]\n",
1433 " else:\n",
1434 " positrons_out = positrons[[\"#x[mm]\", \"xp[mrad]\", \"y[mm]\", \"yp[mrad]\", \"p[MeV/c]\"]]\n",
1435 "\n",
1436 " ## Guassian sampling in time, to keep the RF phases of the Hybrid similar to the conventional scheme\n",
1437 " c = 299792458 #m/s\n",
1438 " if setGaussianTimeShape:\n",
1439 " mean_value = 17.762 #mm/c\n",
1440 " std_value = 1.1950 #mm/c\n",
1441 " num_sample = len(positrons_out)\n",
1442 " gaussian_samples = np.random.normal(mean_value, std_value, num_sample)\n",
1443 " t = gaussian_samples\n",
1444 " else: \n",
1445 " t = positrons.t * 1e-6 * c #ns -> mm/c\n",
1446 " positrons_out.insert(4, 't[mm/c]', t)\n",
1447 "\n",
1448 " ## Save the positron phase-space to a txt file\n",
1449 " name_out = name.replace('output', 'positrons')\n",
1450 " positrons_out.to_csv(outpath + name_out + \".dat\", index=False, sep=' ') \n",
1451 " print(\"positron beam phase-space converted to RF-Track format!\")"
1452 ]
1453 },
1454 {
1455 "cell_type": "markdown",
1456 "id": "4d7e77d2-b658-4e68-ab20-80c980242297",
1457 "metadata": {},
1458 "source": [
1459 "## Export the results of the simulation"
1460 ]
1461 },
1462 {
1463 "cell_type": "code",
1464 "execution_count": null,
1465 "id": "d7aa1f40-3ca6-4c74-9350-9ee4542010dc",
1466 "metadata": {
1467 "tags": []
1468 },
1469 "outputs": [],
1470 "source": [
1471 "# Fill the result list with those of the current simulation\n",
1472 "case_list.append(purename)\n",
1473 "yield_pos_list.append(round(data_out[(data_out.particle==\"e+\")].shape[0] / Nevents, 2))\n",
1474 "yield_ele_list.append(round(data_out[(data_out.particle==\"e-\")].shape[0] / Nevents, 2))\n",
1475 "yield_ph_list.append(round(data_out[(data_out.particle==\"gamma\")].shape[0] / Nevents, 2))\n",
1476 "yield_n_list.append(round(data_out[(data_out.particle==\"neutron\")].shape[0] / Nevents, 2))\n",
1477 "pos_mean_E_list.append(round(Emean_pos, 2))\n",
1478 "pos_spread_E_list.append(round(Espread_pos, 2)) \n",
1479 "pos_mean_div_fit_list.append(round(np.mean(divergence), 2))\n",
1480 "pos_mean_size_fit_list.append(round(np.mean([np.abs(beamSizeX[1]*10), np.abs(beamSizeY[1]*10)]), 2))\n",
1481 "Edep_rad_list.append(round(np.sum(edep_rad)/Nevents, 2))\n",
1482 "Edep_conv_list.append(round(np.sum(edep_conv)*1e-3/Nevents, 2))\n",
1483 "PEDD_list.append(round(PEDD, 2))"
1484 ]
1485 },
1486 {
1487 "cell_type": "code",
1488 "execution_count": null,
1489 "id": "f743be58-cfec-4993-bde2-c224f5bd9094",
1490 "metadata": {
1491 "tags": []
1492 },
1493 "outputs": [],
1494 "source": [
1495 "# Fill a dictionary with the results\n",
1496 "results = { \"case\": case_list,\n",
1497 " \"yield_e+\": yield_pos_list,\n",
1498 " \"yield_e-\": yield_ele_list,\n",
1499 " \"yield_ph\": yield_ph_list,\n",
1500 " \"yield_n\": yield_n_list,\n",
1501 " \"e+_mean_E[MeV]\": pos_mean_E_list,\n",
1502 " \"e+_spread_E[sigma/mu]\": pos_spread_E_list,\n",
1503 " \"e+_mean_div_fit[mrad]\": pos_mean_div_fit_list,\n",
1504 " \"e+_mean_size_fit[mm]\": pos_mean_size_fit_list,\n",
1505 " \"Edep_rad[MeV/e-]\": Edep_rad_list,\n",
1506 " \"Edep_conv[GeV/e-]\": Edep_conv_list,\n",
1507 " \"PEDD[MeV/(mm^3*e-)]\": PEDD_list\n",
1508 " }\n",
1509 "\n",
1510 "# Get a dataframe of the results\n",
1511 "df_results = pd.DataFrame.from_dict(results)\n",
1512 "\n",
1513 "# Export the results in a json file.\n",
1514 "# pickle works pretty much the same, \n",
1515 "# but it I cannot directly have a look at it.\n",
1516 "import json \n",
1517 "print(\"\\nWriting json file with these data:\")\n",
1518 "#for s in results:\n",
1519 "# print(\"\\t\", s, \":\", results[s])\n",
1520 "outname = outpath + name.replace(\"output\", \"results\") \n",
1521 "with open(outname + \".json\", \"w\") as file: \n",
1522 " json.dump(results, file)\n",
1523 "file.close()\n",
1524 "\n",
1525 "# Print results\n",
1526 "df_results"
1527 ]
1528 },
1529 {
1530 "cell_type": "markdown",
1531 "id": "e4e92d8d-24ee-44d3-b8b6-1e6ddc382119",
1532 "metadata": {},
1533 "source": [
1534 "## Read scoring_ntuple2: features of the particles that exit (towards any direction) the radiator crystal"
1535 ]
1536 },
1537 {
1538 "cell_type": "code",
1539 "execution_count": null,
1540 "id": "95150758-427a-4b05-b566-4fa9e282520b",
1541 "metadata": {
1542 "tags": []
1543 },
1544 "outputs": [],
1545 "source": [
1546 "analyze_scoring_ntuple2 = True\n",
1547 "\n",
1548 "if analyze_scoring_ntuple2:\n",
1549 " df2 = rf['scoring_ntuple2'].arrays(library='pd')\n",
1550 " df2['p'] = np.sqrt(df2['px']**2+df2['py']**2+df2['pz']**2)\n",
1551 " #df2.describe()"
1552 ]
1553 },
1554 {
1555 "cell_type": "code",
1556 "execution_count": null,
1557 "id": "f24917dd-72b6-49e8-8f50-95c15d645f8e",
1558 "metadata": {
1559 "tags": []
1560 },
1561 "outputs": [],
1562 "source": [
1563 "part_sel = 'gamma'\n",
1564 "p_low = 0 #MeV/c\n",
1565 "p_high = 100 #MeV/c\n",
1566 "\n",
1567 "if analyze_scoring_ntuple2:\n",
1568 " if len(df2) > 0:\n",
1569 " df2_sel = df2[(df2.particle == part_sel)].copy()\n",
1570 " df2_sel = df2_sel[(df2_sel.p >= p_low) & (df2_sel.p <= p_high)].copy()\n",
1571 " #df2_sel.describe()"
1572 ]
1573 },
1574 {
1575 "cell_type": "code",
1576 "execution_count": null,
1577 "id": "94034573-57ae-479a-a332-20fcf908b900",
1578 "metadata": {},
1579 "outputs": [],
1580 "source": [
1581 "if analyze_scoring_ntuple2 and len(df2) > 0:\n",
1582 " %matplotlib ipympl\n",
1583 "\n",
1584 " plt.figure(figsize=(10, 10))\n",
1585 " Npoint = 10000\n",
1586 " fs = 18\n",
1587 "\n",
1588 " ax3D = plt.axes(projection=\"3d\")\n",
1589 " scatt = ax3D.scatter3D(df2_sel.z[:Npoint], df_sel2.x[:Npoint], df2_sel.y[:Npoint], c=df2_sel.p[:Npoint], \\\n",
1590 " alpha=0.7, marker='.', cmap=plt.get_cmap('hsv'))\n",
1591 " cbar = fig.colorbar(scatt, ax=ax3D, shrink=0.5, aspect=5)\n",
1592 " cbar.set_label(\"p (MeV/c)\", rotation=90, labelpad=16, fontsize=fs)\n",
1593 " ax3D.set_xlabel(\"z (mm)\", fontsize=fs)\n",
1594 " ax3D.set_ylabel(\"x (mm)\", fontsize=fs)\n",
1595 " ax3D.set_zlabel(\"y (mm)\", fontsize=fs)\n",
1596 " ax3D.xaxis.set_inverted(True)\n",
1597 "\n",
1598 " elev = None\n",
1599 " azim = 120\n",
1600 " roll = None\n",
1601 " #ax3D.view_init(elev, azim, roll)\n",
1602 "\n",
1603 " if saveFigs:\n",
1604 " plt.savefig(outpath + '3Ddistrib_' + purename + '.jpg')\n",
1605 " plt.show() "
1606 ]
1607 },
1608 {
1609 "cell_type": "code",
1610 "execution_count": null,
1611 "id": "3c22a840-2fb2-4a95-aafe-dea4a9dce6d7",
1612 "metadata": {},
1613 "outputs": [],
1614 "source": []
1615 }
1616 ],
1617 "metadata": {
1618 "kernelspec": {
1619 "display_name": "Python 3 (ipykernel)",
1620 "language": "python",
1621 "name": "python3"
1622 },
1623 "language_info": {
1624 "codemirror_mode": {
1625 "name": "ipython",
1626 "version": 3
1627 },
1628 "file_extension": ".py",
1629 "mimetype": "text/x-python",
1630 "name": "python",
1631 "nbconvert_exporter": "python",
1632 "pygments_lexer": "ipython3",
1633 "version": "3.12.9"
1634 }
1635 },
1636 "nbformat": 4,
1637 "nbformat_minor": 5
1638 }
