Gabriele Cosmo
1639 lines
60 KiB
Plaintext
1639 lines
60 KiB
Plaintext
{
|
|
"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)], y[(x>-fitrange_ang) & (x<fitrange_ang)]\n",
|
|
" try:\n",
|
|
" par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
|
|
" except:\n",
|
|
" par = (1, 0, np.inf)\n",
|
|
" xplot = np.linspace(-fitrange_ang, fitrange_ang, 200)\n",
|
|
" xplot2 = np.linspace(-200, -fitrange_ang, 2000)\n",
|
|
" xplot3 = np.linspace(fitrange_ang, 200, 2000)\n",
|
|
" label = (\"at crystal output\"+r\" ($\\it{\\gamma}$ only)\" if inputType==\"all\" else \\\n",
|
|
" \"at amorphous target input\"+r\" ($\\it{e}^-$ only)\")+\", Gaussian\\n\"+r\"fit (in $\\pm$ %.2f mrad)\" % \\\n",
|
|
" fitrange_ang+r\" $\\it{\\sigma}$ = %.2f mrad\" % abs(par[2])\n",
|
|
" if abs(par[2]) < np.inf:\n",
|
|
" ax[i].plot(xplot, gaussian(xplot, *par), c=\"C0\", alpha=0.5, label=label)\n",
|
|
" ax[i].plot(xplot2, gaussian(xplot2, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
|
|
" ax[i].plot(xplot3, gaussian(xplot3, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
|
|
" \n",
|
|
" divergenceIn[i] = abs(par[2])\n",
|
|
" angleStdIn[i] = series[i].std()\n",
|
|
"\n",
|
|
" ax[i].set_yscale(\"linear\")\n",
|
|
" ax[i].set_xlim(ang_plot_lim)\n",
|
|
" ax[i].set_title(\" \"+(\"horizontal\" if i==0 else \"vertical\"), loc=\"left\")\n",
|
|
" ax[i].legend(loc=\"upper right\", fontsize=\"xx-small\")\n",
|
|
" ax[i].grid(True)\n",
|
|
"\n",
|
|
"series = [data_out[(data_out.particle==\"e+\")].thx, data_out[(data_out.particle==\"e+\")].thy]\n",
|
|
"for i in (0, 1):\n",
|
|
" hist = ax[i].hist(series[i], bins=800, range=(-1000, 1000), histtype=\"step\", color=\"C1\")\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)], y[(x>-fitrange_ang) & (x<fitrange_ang)]\n",
|
|
" try:\n",
|
|
" par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
|
|
" except:\n",
|
|
" par = (1, 0, np.inf)\n",
|
|
" xplot = np.linspace(-fitrange_ang, fitrange_ang, 200)\n",
|
|
" xplot2 = np.linspace(-200, -fitrange_ang, 2000)\n",
|
|
" xplot3 = np.linspace(fitrange_ang, 200, 2000)\n",
|
|
" label = \"at amorphous target output\"+r\" ($\\it{e}^+$ only)\"+\", Gaussian\\n\"+r\"fit (in $\\pm$ %.2f mrad)\" % \\\n",
|
|
" fitrange_ang+r\" $\\it{\\sigma}$ = %.2f mrad\" % abs(par[2])\n",
|
|
" if abs(par[2]) < np.inf:\n",
|
|
" ax[i].plot(xplot, gaussian(xplot, *par), c=\"C1\", alpha=0.5, label=label)\n",
|
|
" ax[i].plot(xplot2, gaussian(xplot2, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
|
|
" ax[i].plot(xplot3, gaussian(xplot3, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
|
|
" divergence[i] = abs(par[2])\n",
|
|
" angleStd[i] = series[i].std()\n",
|
|
"\n",
|
|
" ax[i].set_yscale(\"linear\")\n",
|
|
" ax[i].set_xlim(ang_plot_lim)\n",
|
|
" ax[i].set_title(\" \"+(\"horizontal\" if i==0 else \"vertical\"), loc=\"left\")\n",
|
|
" ax[i].legend(loc=\"upper right\", fontsize=\"xx-small\")\n",
|
|
" ax[i].grid(True)\n",
|
|
"\n",
|
|
"fig.supylabel(\"Counts\", fontsize=\"small\")\n",
|
|
"ax[1].set_xlabel(r\"$\\it{\\theta}$ [mrad]\")\n",
|
|
"fig.tight_layout()\n",
|
|
"\n",
|
|
"if saveFigs:\n",
|
|
" plt.savefig(outpath + 'eBeam_theta_distrib_' + purename + '.jpg')\n",
|
|
"\n",
|
|
"print(\"e+ beam divergence along X and Y from Gauss fit in (-%.2f, %.2f) mrad range: [-%.2f, %.2f] mrad\" \\\n",
|
|
" % (fitrange_ang, fitrange_ang, divergence[0], divergence[1]))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "eeeef894-b2ab-4f60-90af-68faa33e6562",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Plot the spatial distribution"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "df9d997c-d378-40f2-83e8-558c5cf504b7",
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, ax = plt.subplots(ncols=2, nrows=3, figsize=(style['figure.figsize'][0], 14.0))\n",
|
|
"fig.subplots_adjust(wspace=-1)\n",
|
|
"\n",
|
|
"beamSizeX, beamSizeY = [0, 0, 0], [0, 0, 0]\n",
|
|
"\n",
|
|
"ax[0,0].set_title(r\"$\\it{e}^{+}$ only\")\n",
|
|
"ax[0,1].set_title(r\"$\\it{e}^{+}$ only\")\n",
|
|
"\n",
|
|
"ax[0,1].hist2d(data_out[(data_out.particle==\"e+\")].x, data_out[(data_out.particle==\"e+\")].y, \\\n",
|
|
" bins=nBins, range=(xRange, xRange), norm=LogNorm())\n",
|
|
"ax[0,1].grid(True)\n",
|
|
"ax[0,0].hist2d(data_out[(data_out.particle==\"e+\")].x, data_out[(data_out.particle==\"e+\")].y, \\\n",
|
|
" bins=nBins, range=(xRange1, xRange1), norm=LogNorm())\n",
|
|
"ax[0,0].grid(True)\n",
|
|
"\n",
|
|
"gs = ax[1,0].get_gridspec()\n",
|
|
"for ax0 in ax[1,:]:\n",
|
|
" ax0.remove()\n",
|
|
"axdown = fig.add_subplot(gs[1,:])\n",
|
|
"\n",
|
|
"if len(data_in) > 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)], y[(x>-fitrange) & (x<fitrange)]\n",
|
|
" try:\n",
|
|
" par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
|
|
" except:\n",
|
|
" par = (1, 0, np.inf)\n",
|
|
" xplot = np.linspace(-fitrange, fitrange, 200)\n",
|
|
" xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
|
|
" xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
|
|
" label = (\"at crystal output\"+r\" ($\\it{\\gamma}$ only)\" if inputType==\"all\" else \\\n",
|
|
" \"at amorphous target input\"+r\" ($\\it{e}^-$ only)\")+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
|
|
" if abs(par[2]) < np.inf:\n",
|
|
" axdown.plot(xplot, gaussian(xplot, *par), c=\"C0\", alpha=0.5, label=label)\n",
|
|
" axdown.plot(xplot2, gaussian(xplot2, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
|
|
" axdown.plot(xplot3, gaussian(xplot3, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
|
|
" beamSizeX[0] = par[2]\n",
|
|
"\n",
|
|
"hist = axdown.hist(data_out[(data_out.particle==\"e+\")].x, bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C1\")\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)], y[(x>-fitrange) & (x<fitrange)]\n",
|
|
"try:\n",
|
|
" par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
|
|
"except:\n",
|
|
" par = (1, 0, np.inf)\n",
|
|
"xplot = np.linspace(-fitrange, fitrange, 200)\n",
|
|
"xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
|
|
"xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
|
|
"label = \"at amorphous target output\"+r\" ($\\it{e}^{+}$ only)\"+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
|
|
"if abs(par[2]) < np.inf:\n",
|
|
" axdown.plot(xplot, gaussian(xplot, *par), c=\"C1\", alpha=0.5, label=label)\n",
|
|
" axdown.plot(xplot2, gaussian(xplot2, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
|
|
" axdown.plot(xplot3, gaussian(xplot3, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
|
|
"beamSizeX[1] = par[2]\n",
|
|
"\n",
|
|
"gs = ax[2,0].get_gridspec()\n",
|
|
"for ax0 in ax[2,:]:\n",
|
|
" ax0.remove()\n",
|
|
"axdown2 = fig.add_subplot(gs[2,:])\n",
|
|
"\n",
|
|
"if len(data_in) > 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)], y[(x>-fitrange) & (x<fitrange)]\n",
|
|
" try:\n",
|
|
" par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
|
|
" except:\n",
|
|
" par = (1, 0, np.inf)\n",
|
|
" xplot = np.linspace(-fitrange, fitrange, 200)\n",
|
|
" xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
|
|
" xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
|
|
" label = (\"at crystal output\"+r\" ($\\it{\\gamma}$ only)\" if inputType==\"all\" else \\\n",
|
|
" \"at amorphous target input\"+r\" ($\\it{e}^-$ only)\")+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
|
|
" if abs(par[2]) < np.inf:\n",
|
|
" axdown2.plot(xplot, gaussian(xplot, *par), c=\"C0\", alpha=0.5, label=label)\n",
|
|
" axdown2.plot(xplot2, gaussian(xplot2, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
|
|
" axdown2.plot(xplot3, gaussian(xplot3, *par), c=\"C0\", ls=\":\", alpha=0.5)\n",
|
|
" beamSizeY[0] = par[2]\n",
|
|
"\n",
|
|
"hist = axdown2.hist(data_out[(data_out.particle==\"e+\")].y, bins=nBins, histtype=\"step\", range=xRange, density=True, color=\"C1\")\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)], y[(x>-fitrange) & (x<fitrange)]\n",
|
|
"try:\n",
|
|
" par, _ = curve_fit(f=gaussian, xdata=x, ydata=y, bounds=(-np.inf, np.inf))\n",
|
|
"except:\n",
|
|
" par = (1, 0, np.inf)\n",
|
|
"xplot = np.linspace(-fitrange, fitrange, 200)\n",
|
|
"xplot2 = np.linspace(xRange[0], -fitrange, 2000)\n",
|
|
"xplot3 = np.linspace(fitrange, xRange[1], 2000)\n",
|
|
"label = \"at amorphous target output\"+r\" ($\\it{e}^{+}$ only)\"+\",\\nGaussian fit \"+r\"$\\it{\\sigma}$ = %.3f mm\" % abs(par[2]*10)\n",
|
|
"if abs(par[2]) < np.inf:\n",
|
|
" axdown2.plot(xplot, gaussian(xplot, *par), c=\"C1\", alpha=0.5, label=label)\n",
|
|
" axdown2.plot(xplot2, gaussian(xplot2, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
|
|
" axdown2.plot(xplot3, gaussian(xplot3, *par), c=\"C1\", ls=\":\", alpha=0.5)\n",
|
|
"beamSizeY[1] = par[2]\n",
|
|
"\n",
|
|
"ax[0,0].set_ylabel(r\"$\\it{y}$ at amorphous\"+\"\\ntarget output [cm]\")\n",
|
|
"ax[0,0].set_xlabel(r\" $\\it{x}$ at amorphous target output [cm]\", loc=\"left\")\n",
|
|
"axdown.set_xlabel(r\"$\\it{x}$ [cm]\")\n",
|
|
"axdown.set_ylabel(r\"$\\frac{\\mathrm{d}\\it{N}}{\\it{N}\\mathrm{d}\\it{x}}\\left[\\frac{1}{\\mathrm{cm}}\\right]$\")\n",
|
|
"axdown.set_xlim((xRange[0], xRange[1]))\n",
|
|
"axdown.set_ylim((0, yUpperLim))\n",
|
|
"axdown.legend(loc=\"upper left\", fontsize=16)\n",
|
|
"axdown.grid(True)\n",
|
|
"axdown2.set_xlabel(r\"$\\it{y}$ [cm]\")\n",
|
|
"axdown2.set_ylabel(r\"$\\frac{\\mathrm{d}\\it{N}}{\\it{N}\\mathrm{d}\\it{y}}\\left[\\frac{1}{\\mathrm{cm}}\\right]$\")\n",
|
|
"axdown2.set_xlim((xRange[0], xRange[1]))\n",
|
|
"axdown2.set_ylim((0, yUpperLim))\n",
|
|
"axdown2.legend(loc=\"upper left\", fontsize=16)\n",
|
|
"axdown2.grid(True)\n",
|
|
"fig.tight_layout()\n",
|
|
"\n",
|
|
"if saveFigs:\n",
|
|
" plt.savefig(outpath + 'eBeam_spatial_distrib_' + purename + '.jpg')\n",
|
|
"\n",
|
|
"beamSizeXStd = [data_in.x.std(), data_out[(data_out.particle==\"e+\")].x.std(), data_out[(data_out.particle==\"e+\")].x.std()]\n",
|
|
"beamSizeYStd = [data_in.y.std(), data_out[(data_out.particle==\"e+\")].y.std(), data_out[(data_out.particle==\"e+\")].y.std()]\n",
|
|
"\n",
|
|
"print(\"e+ beam size along X and Y from Gauss fit in (-%.2f, %.2f) mm range: [-%.2f, %.2f] mm\" % \\\n",
|
|
" (fitrange, fitrange, np.abs(beamSizeX[1]*10), np.abs(beamSizeY[1]*10)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "f33236b5-250d-4560-83b0-c9fc7574704f",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Plot the photon spectrum emitted by the Radiator crystal "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "30aabcff-7a56-41af-a3b4-1990b9f91dfa",
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"if len(data_in) > 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
|
|
}
|