anaMamiRad.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **anaMamiRad**(iation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "import re\n",
    "import math\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SETTINGS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **run list**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataNames = {  # {physics_filenames : bkg_filenames} -- full statistics, select datasets to open with dataNamesToOpen\n",
    "\n",
    "# 2021 W 1mm 111 by Matek+IJCLab\n",
    "# 855 MeV only\n",
    "    \"data21/radiation/Amorph_W111_MATEK_IJC_06_11_2021_15h30_MCA1.txt\" : \"data21/radiation/background_06_11_2021_15h04_MCA1.txt\",  # after irradiation\n",
    "    \"data21/radiation/Amorph_W111_MATEK_IJC_06_11_2021_15h30_MCA6.txt\" : \"data21/radiation/background_06_11_2021_15h04_MCA6.txt\",  # after irradiation\n",
    "    \"data21/radiation/AXIS_W111_MATEK_IJC_06_11_2021_15h18_MCA1.txt\" : \"data21/radiation/background_06_11_2021_15h04_MCA1.txt\",  # after irradiation\n",
    "    \"data21/radiation/AXIS_W111_MATEK_IJC_06_11_2021_15h18_MCA6.txt\" : \"data21/radiation/background_06_11_2021_15h04_MCA6.txt\",  # after irradiation\n",
    "    \"data21/radiation/Amorph_W111_Matek_IJC_run1_mca1.txt\" : \"data21/radiation/background_run1_mca1.txt\",  # before irradiation\n",
    "    \"data21/radiation/Amorph_W111_Matek_IJC_run1_mca6.txt\" : \"data21/radiation/background_run1_mca6.txt\",  # before irradiation\n",
    "    \"data21/radiation/AXIS_W111_Matek_IJC_run1_mca1.txt\" : \"data21/radiation/background_run2_mca1.txt\",  # before irradiation\n",
    "    \"data21/radiation/AXIS_W111_Matek_IJC_run1_mca6.txt\" : \"data21/radiation/background_run2_mca6.txt\",  # before irradiation\n",
    "    \n",
    "# 2021 W 1mm 111 by Matek+Ferrara (same as 2020)\n",
    "# 855 MeV only\n",
    "    \"data21/radiation/Amorph_W111_Matek_Fe_run1_mca1.txt\" : \"data21/radiation/background_run2_mca1.txt\",\n",
    "    \"data21/radiation/Amorph_W111_Matek_Fe_run1_mca6.txt\" : \"data21/radiation/background_run2_mca6.txt\",\n",
    "    \"data21/radiation/AXIS_W111_Matek_Fe_run1_mca1.txt\" : \"data21/radiation/background_run3_mca1.txt\",\n",
    "    \"data21/radiation/AXIS_W111_Matek_Fe_run1_mca6.txt\" : \"data21/radiation/background_run3_mca6.txt\",\n",
    "    \n",
    "# 2020 W 1mm 111 by Matek\n",
    "# 855 MeV only\n",
    "    \"data20/radiation/W 1mm/W1mm_crad1.8197_phi-0.6903_70mm_MCA1.txt\" : \"data20/radiation/W 1mm/W1mm_background21_MCA1.txt\",  # axis\n",
    "    \"data20/radiation/W 1mm/W1mm_crad1.8197_phi-0.6903_70mm_MCA6.txt\" : \"data20/radiation/W 1mm/W1mm_background21_MCA6.txt\",\n",
    "    \"data20/radiation/W 1mm/W1mm_crad0.6197_phi4.3097_70mm_MCA1.txt\" : \"data20/radiation/W 1mm/W1mm_background21_MCA1.txt\",  # random\n",
    "    \"data20/radiation/W 1mm/W1mm_crad0.6197_phi4.3097_70mm_MCA6.txt\" : \"data20/radiation/W 1mm/W1mm_background21_MCA6.txt\",\n",
    "\n",
    "# 2019 YAP 0.5mm 110\n",
    "# 855 MeV only\n",
    "    # 855MeV, aperture 70mm, 4mm collimator in\n",
    "    \"data19/YAPaxis-mca1.txt\" : \"data19/Bkg8-09-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAPaxis-mca6.txt\" : \"data19/Bkg8-09-05-19-NB-mca6.txt\",\n",
    "    \"data19/YAPamorphous-mca1.txt\" : \"data19/Bkg8-09-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAPamorphous-mca6.txt\" : \"data19/Bkg8-09-05-19-NB-mca6.txt\",\n",
    "\n",
    "# 2019 PWO 0.5mm 021\n",
    "# 855 & 600 MeV\n",
    "    # 855MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/PWO-021-axis-855MeV-mca1.txt\" : \"data19/Bkg35-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWO-021-axis-855MeV-mca6.txt\" : \"data19/Bkg35-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/PWO-021-amorphous-855MeV-mca1.txt\" : \"data19/Bkg35-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWO-021-amorphous-855MeV-mca6.txt\" : \"data19/Bkg35-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 600MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/PWO-021-axis-600MeV-mca1.txt\" : \"data19/Bkg35-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWO-021-axis-600MeV-mca6.txt\" : \"data19/Bkg35-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/PWO-021-amorphous-600MeV-mca1.txt\" : \"data19/Bkg35-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWO-021-amorphous-600MeV-mca6.txt\" : \"data19/Bkg35-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "# 2019 PWO 0.5mm 110\n",
    "# 855 MeV only\n",
    "    # 855MeV, aperture 70mm, 4mm collimator in\n",
    "    \"data19/PWOaxis-mca1.txt\" : \"data19/Bkg8-09-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWOaxis-mca6.txt\" : \"data19/Bkg8-09-05-19-NB-mca6.txt\",\n",
    "    \"data19/PWOamorphous-mca1.txt\" : \"data19/Bkg8-09-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWOamorphous-mca6.txt\" : \"data19/Bkg8-09-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 855MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/PWOaxis-mca1bis.txt\" : \"data19/Bkg20-10-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWOaxis-mca6bis.txt\" : \"data19/Bkg20-10-05-19-NB-mca6.txt\",\n",
    "    \"data19/PWOamorphous-mca1bis.txt\" : \"data19/Bkg20-10-05-19-NB-mca1.txt\",\n",
    "    \"data19/PWOamorphous-mca6bis.txt\" : \"data19/Bkg20-10-05-19-NB-mca6.txt\",\n",
    "    \n",
    "# 2019 YAG:Ce 0.5mm 100 \n",
    "# 855 & 600 MeV\n",
    "    # 855MeV, aperture 70mm, 4mm collimator in\n",
    "    \"data19/YAGCEaxis-mca1.txt\" : \"data19/Bkg10-09-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAGCEaxis-mca6.txt\" : \"data19/Bkg10-09-05-19-NB-mca6.txt\",\n",
    "    \"data19/YAGCEamorphous-mca1.txt\" : \"data19/Bkg10-09-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAGCEamorphous-mca6.txt\" : \"data19/Bkg10-09-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 855MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/YAGCE2axis-mca1.txt\" : \"data19/Bkg18-10-05-19-NB-mca1.txt\",  \n",
    "    \"data19/YAGCE2axis-mca6.txt\" : \"data19/Bkg18-10-05-19-NB-mca6.txt\",\n",
    "    \"data19/YAGCE2amorphous-mca1.txt\" : \"data19/Bkg18-10-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAGCE2amorphous-mca6.txt\" : \"data19/Bkg18-10-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 600MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/YagCe-axis-600MeV-mca1.txt\" : \"data19/Bkg36-12-05-19-NB-mca1.txt\",  \n",
    "    \"data19/YagCe-axis-600MeV-mca6.txt\" : \"data19/Bkg36-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/YagCe-amorphous-600MeV-mca1.txt\" : \"data19/Bkg36-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/YagCe-amorphous-600MeV-mca6.txt\" : \"data19/Bkg36-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "# 2019 YAG 0.5mm 100 \n",
    "# 855 & 600 MeV\n",
    "    # 855MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/YAG-axis-855MeV-bis-mca1.txt\" : \"data19/Bkg44-12-05-19-NB-mca1.txt\",  \n",
    "    \"data19/YAG-axis-855MeV-bis-mca6.txt\" : \"data19/Bkg44-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/YAG-amorphous-855MeV-bis-mca1.txt\" : \"data19/Bkg44-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAG-amorphous-855MeV-bis-mca6.txt\" : \"data19/Bkg44-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 600MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/Yag-axis-600MeV-mca1.txt\" : \"data19/Bkg37-12-05-19-NB-mca1.txt\",  \n",
    "    \"data19/Yag-axis-600MeV-mca6.txt\" : \"data19/Bkg37-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/Yag-amorphous-600MeV-mca1.txt\" : \"data19/Bkg37-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/Yag-amorphous-600MeV-mca6.txt\" : \"data19/Bkg37-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "# 2019 YAG 0.1mm 100 \n",
    "# 855 & 600 MeV\n",
    "    # 855MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/YAGthin-axis-855MeV-mca1.txt\" : \"data19/Bkg43-12-05-19-NB-mca1.txt\",  \n",
    "    \"data19/YAGthin-axis-855MeV-mca6.txt\" : \"data19/Bkg43-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/YAGthin-amorphous-855MeV-mca1.txt\" : \"data19/Bkg43-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAGthin-amorphous-855MeV-mca6.txt\" : \"data19/Bkg43-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 600MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/YAGthin-axis-600MeV-mca1.txt\" : \"data19/Bkg27-11-05-19-NB-mca1.txt\",  \n",
    "    \"data19/YAGthin-axis-600MeV-mca6.txt\" : \"data19/Bkg27-11-05-19-NB-mca6.txt\",\n",
    "    \"data19/YAGthin-amorphous-600MeV-mca1.txt\" : \"data19/Bkg27-11-05-19-NB-mca1.txt\",\n",
    "    \"data19/YAGthin-amorphous-600MeV-mca6.txt\" : \"data19/Bkg27-11-05-19-NB-mca6.txt\",\n",
    "    \n",
    "# 2019 BGO 0.1mm 111 \n",
    "# 855 & 600 MeV\n",
    "    # 855MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/BGO111-axis-855MeV-mca1.txt\" : \"data19/Bkg46-12-05-19-NB-mca1.txt\",  \n",
    "    \"data19/BGO111-axis-855MeV-mca6.txt\" : \"data19/Bkg46-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/BGO111-amorphous-855MeV-mca1.txt\" : \"data19/Bkg46-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/BGO111-amorphous-855MeV-mca6.txt\" : \"data19/Bkg46-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "    # 600MeV, aperture 40mm, 4mm collimator out\n",
    "    \"data19/BGO111-axis-600MeV-mca1.txt\" : \"data19/Bkg39-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/BGO111-axis-600MeV-mca6.txt\" : \"data19/Bkg39-12-05-19-NB-mca6.txt\",\n",
    "    \"data19/BGO111-amorphous-600MeV-mca1.txt\" : \"data19/Bkg39-12-05-19-NB-mca1.txt\",\n",
    "    \"data19/BGO111-amorphous-600MeV-mca6.txt\" : \"data19/Bkg39-12-05-19-NB-mca6.txt\",\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **analysis settings**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsRatios = {  # list of run pairs whose ratio can be plotted -- {numerator_filenames : denominator_filenames}\n",
    "    \n",
    "# 2021 W 1mm 111 by Matek+IJCLab\n",
    "    \"data21/radiation/AXIS_W111_MATEK_IJC_06_11_2021_15h18_MCA1.txt\" : \"data21/radiation/Amorph_W111_MATEK_IJC_06_11_2021_15h30_MCA1.txt\",\n",
    "    \"data21/radiation/AXIS_W111_MATEK_IJC_06_11_2021_15h18_MCA6.txt\" : \"data21/radiation/Amorph_W111_MATEK_IJC_06_11_2021_15h30_MCA6.txt\",\n",
    "    \"data21/radiation/AXIS_W111_Matek_IJC_run1_mca1.txt\" : \"data21/radiation/Amorph_W111_Matek_IJC_run1_mca1.txt\",  \n",
    "    \"data21/radiation/AXIS_W111_Matek_IJC_run1_mca6.txt\" : \"data21/radiation/Amorph_W111_Matek_IJC_run1_mca6.txt\", \n",
    "    \n",
    "# 2021 W 1mm 111 by Matek+Ferrara (same as 2020)\n",
    "    \"data21/radiation/AXIS_W111_Matek_Fe_run1_mca1.txt\" : \"data21/radiation/Amorph_W111_Matek_Fe_run1_mca1.txt\",  \n",
    "    \"data21/radiation/AXIS_W111_Matek_Fe_run1_mca6.txt\" : \"data21/radiation/Amorph_W111_Matek_Fe_run1_mca6.txt\", \n",
    "    \n",
    "# 2020 W 1mm 111 by Matek\n",
    "    \"data20/radiation/W 1mm/W1mm_crad1.8197_phi-0.6903_70mm_MCA1.txt\" : \"data20/radiation/W 1mm/W1mm_crad0.6197_phi4.3097_70mm_MCA1.txt\",\n",
    "    \"data20/radiation/W 1mm/W1mm_crad1.8197_phi-0.6903_70mm_MCA6.txt\" : \"data20/radiation/W 1mm/W1mm_crad0.6197_phi4.3097_70mm_MCA6.txt\",\n",
    "    \n",
    "# 2019 YAP 0.5mm 110\n",
    "    \"data19/YAPaxis-mca1.txt\" : \"data19/YAPamorphous-mca1.txt\",\n",
    "    \"data19/YAPaxis-mca6.txt\" : \"data19/YAPamorphous-mca6.txt\",\n",
    "    \n",
    "# 2019 PWO 0.5mm 021\n",
    "    \"data19/PWO-021-axis-855MeV-mca1.txt\" : \"data19/PWO-021-amorphous-855MeV-mca1.txt\",\n",
    "    \"data19/PWO-021-axis-855MeV-mca6.txt\" : \"data19/PWO-021-amorphous-855MeV-mca6.txt\",\n",
    "    \"data19/PWO-021-axis-600MeV-mca1.txt\" : \"data19/PWO-021-amorphous-600MeV-mca1.txt\",\n",
    "    \"data19/PWO-021-axis-600MeV-mca6.txt\" : \"data19/PWO-021-amorphous-600MeV-mca6.txt\",\n",
    "    \n",
    "# 2019 PWO 0.5mm 110\n",
    "    \"data19/PWOaxis-mca1.txt\" : \"data19/PWOamorphous-mca1.txt\",\n",
    "    \"data19/PWOaxis-mca6.txt\" : \"data19/PWOamorphous-mca6.txt\",\n",
    "    \"data19/PWOaxis-mca1bis.txt\" : \"data19/PWOamorphous-mca1bis.txt\",\n",
    "    \"data19/PWOaxis-mca6bis.txt\" : \"data19/PWOamorphous-mca6bis.txt\",\n",
    "    \n",
    "# 2019 YAG:Ce 0.5mm 100 \n",
    "    \"data19/YAGCEaxis-mca1.txt\" : \"data19/YAGCEamorphous-mca1.txt\",\n",
    "    \"data19/YAGCEaxis-mca6.txt\" : \"data19/YAGCEamorphous-mca6.txt\",\n",
    "    \"data19/YAGCE2axis-mca1.txt\" : \"data19/YAGCE2amorphous-mca1.txt\",\n",
    "    \"data19/YAGCE2axis-mca6.txt\" : \"data19/YAGCE2amorphous-mca6.txt\",\n",
    "    \"data19/YagCe-axis-600MeV-mca1.txt\" : \"data19/YagCe-amorphous-600MeV-mca1.txt\",\n",
    "    \"data19/YagCe-axis-600MeV-mca6.txt\" : \"data19/YagCe-amorphous-600MeV-mca6.txt\",\n",
    "    \n",
    "# 2019 YAG 0.5mm 100 \n",
    "    \"data19/YAG-axis-855MeV-bis-mca1.txt\" : \"data19/YAG-amorphous-855MeV-bis-mca1.txt\",\n",
    "    \"data19/YAG-axis-855MeV-bis-mca6.txt\" : \"data19/YAG-amorphous-855MeV-bis-mca6.txt\",\n",
    "    \"data19/Yag-axis-600MeV-mca1.txt\" : \"data19/Yag-amorphous-600MeV-mca1.txt\",\n",
    "    \"data19/Yag-axis-600MeV-mca6.txt\" : \"data19/Yag-amorphous-600MeV-mca6.txt\",\n",
    "    \n",
    "# 2019 YAG 0.1mm 100 \n",
    "    \"data19/YAGthin-axis-855MeV-mca1.txt\" : \"data19/YAGthin-amorphous-855MeV-mca1.txt\",\n",
    "    \"data19/YAGthin-axis-855MeV-mca6.txt\" : \"data19/YAGthin-amorphous-855MeV-mca6.txt\",\n",
    "    \"data19/YAGthin-axis-600MeV-mca1.txt\" : \"data19/YAGthin-amorphous-600MeV-mca1.txt\",\n",
    "    \"data19/YAGthin-axis-600MeV-mca6.txt\" : \"data19/YAGthin-amorphous-600MeV-mca6.txt\",\n",
    "    \n",
    "# 2019 BGO 0.1mm 111 \n",
    "    \"data19/BGO111-axis-855MeV-mca1.txt\" : \"data19/BGO111-amorphous-855MeV-mca1.txt\",\n",
    "    \"data19/BGO111-axis-855MeV-mca6.txt\" : \"data19/BGO111-amorphous-855MeV-mca6.txt\",\n",
    "    \"data19/BGO111-axis-600MeV-mca1.txt\" : \"data19/BGO111-amorphous-600MeV-mca1.txt\",\n",
    "    \"data19/BGO111-axis-600MeV-mca6.txt\" : \"data19/BGO111-amorphous-600MeV-mca6.txt\",\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **data input & conditioning settings**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "dataPath = \"./\"  # files path (with / at the end)\n",
    "\n",
    "# selected datasets to be opened -- make sure each run has the correct background run referred to in dataNames\n",
    "dataNamesToOpen = [\n",
    "\"data21/radiation/AXIS_W111_Matek_IJC_run1_mca6.txt\",\n",
    "\"data21/radiation/AXIS_W111_Matek_IJC_run1_mca1.txt\",\n",
    "]\n",
    "\n",
    "linesHeader = 27  # nr. of header lines\n",
    "realTimeLine = 7  # header line with real time value\n",
    "\n",
    "bPrintExtInfo = False  # print extended information for each input file?\n",
    "\n",
    "nBin = 200  # number of data lines to be summed together into a single bin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DATA INPUT & CONDITIONING"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *also set detector calibration here below:*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "###### input data\n",
    "allData = set(dataNamesToOpen + [dataNames[s] for s in dataNames.keys() if (s in dataNamesToOpen)])\n",
    "dictData = {}\n",
    "\n",
    "if len(allData)==0:\n",
    "    print(\"no datasets to be opened!\")\n",
    "\n",
    "for i, s in enumerate(allData):\n",
    "    _, x, y, _ = np.loadtxt(dataPath + s, unpack=True, comments=\"#\")\n",
    "    \n",
    "    header = pd.read_csv(dataPath + s, sep=\"\\n\", squeeze=True, skiprows = lambda i : i > linesHeader)\n",
    "    t = float(re.findall(r\"[-+]?\\d*\\.\\d+|\\d+\", header[realTimeLine])[0])  # acquisition real time\n",
    "    \n",
    "    if bPrintExtInfo:\n",
    "        print(\"OPENED DATASET %d/%d -- %s:\" % (i+1, len(allData), s))\n",
    "        print(header.to_string())\n",
    "    else:\n",
    "        print(\"OPENED DATASET %d/%d -- %s -- acquisition real time t=%.2fs\" % (i+1, len(allData), s, t))\n",
    "    \n",
    "    x = x * 1e-3  # conversion into MeV\n",
    "    \n",
    "    # # # # # # # # # # # # # # # # # # # # # # # #\n",
    "    # recalibrations (after conversion into MeV)\n",
    "    \n",
    "    if (\"data21\" in s): # 2021\n",
    "        if ((\"mca1\" in s) or (\"MCA1\" in s)):\n",
    "            x = x * (1.460 / 0.919)\n",
    "        if ((\"mca6\" in s) or (\"MCA6\" in s)):\n",
    "            m = (2.615-1.460) / (3.876-2.169)\n",
    "            q = 2.615 - 3.876*m\n",
    "            x = m * x + q\n",
    "            \n",
    "    if (\"data20\" in s):  # 2020\n",
    "        if ((\"mca1\" in s) or (\"MCA1\" in s)):\n",
    "            m = (2.615-1.460) / (0.144-0.070)\n",
    "            q = 2.615 - 0.144*m\n",
    "            x = m * x + q\n",
    "        if ((\"mca6\" in s) or (\"MCA6\" in s)):\n",
    "            x = x * (1.460 / 43.23)\n",
    "            \n",
    "    # # # # # # # # # # # # # # # # # # # # # # # #\n",
    "\n",
    "    y = y / t  # spectrum normalisation\n",
    "    dictData[s] = [x, y, t]\n",
    "    print(\"---\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ANALYSIS & PLOTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def rebin(x0, y0, nBin):\n",
    "    x, y = np.empty(0), np.empty(0)\n",
    "    for i in range(0, math.floor(len(x0)/nBin)*nBin, nBin):\n",
    "        x = np.append(x, sum(x0[i:(i+nBin)])/nBin)\n",
    "        y = np.append(y, sum(y0[i:(i+nBin)]))\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **raw spectra** (normalised to the acquisition time only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "if (True):  # plot this?\n",
    "    \n",
    "    plt.close(\"spectra_raw\")\n",
    "    fig, ax = plt.subplots(num=\"spectra_raw\", figsize=[9, 7], nrows=2, sharex=True)\n",
    "\n",
    "    for s in dictData:\n",
    "        x, y = rebin(dictData[s][0], dictData[s][1], nBin)\n",
    "        ave = sum(x*y) / sum(y)\n",
    "        # print(\"%s RAW SPECTRUM \\n\\t--> ave=%f\" % (s, ave))  # uncomment to print average energy per spectrum\n",
    "        ax[0].plot(x, y, label=s)\n",
    "        ax[1].plot(x, y*x)\n",
    "\n",
    "    ax[0].set_yscale(\"log\")\n",
    "    ax[0].grid()\n",
    "    ax[0].set_ylabel(\"photon spectrum [1/s]\")\n",
    "    ax[1].grid()\n",
    "    ax[1].set_ylabel(\"spectral intensity [MeV/s]\")\n",
    "    ax[1].set_xlabel(\"energy [MeV]\")\n",
    "    fig.legend(bbox_to_anchor=(1, 0.5), loc='center left')\n",
    "    plt.tight_layout(pad=1)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **background-subtracted spectra** (normalised to the acquisition time and to the curve area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "if (True):  # plot this?\n",
    "    \n",
    "    plt.close(\"spectra_final\")\n",
    "    fig, ax = plt.subplots(num=\"spectra_final\", figsize=[9, 7], nrows=2, sharex=True)\n",
    "\n",
    "    for s in dataNamesToOpen:\n",
    "        xPhys, yPhys = rebin(dictData[s][0], dictData[s][1], nBin)\n",
    "        xBkg, yBkg = rebin(dictData[dataNames[s]][0], dictData[dataNames[s]][1], nBin)\n",
    "        x, y = xPhys, (yPhys-yBkg)/abs(xPhys[1]-xPhys[0])\n",
    "        ave = sum(x*y) / sum(y)\n",
    "        print(\"%s BKG-SUBTRACTED SPECTRUM \\n\\t--> ave=%fMeV\" % (s, ave))  # uncomment to print average energy per spectrum\n",
    "        ax[0].plot(x, y, label=s)\n",
    "        ax[1].plot(x, y*x)\n",
    "\n",
    "    ax[0].set_yscale(\"log\")\n",
    "    ax[0].grid()\n",
    "    ax[0].set_ylabel(r\"photon spectrum $[1 / (\\mathrm{MeV} \\cdot \\mathrm{s})]$\")\n",
    "    ax[1].grid()\n",
    "    ax[1].set_ylabel(\"spectral intensity [1/s]\")\n",
    "    ax[1].set_xlabel(\"energy [MeV]\")\n",
    "    fig.legend(bbox_to_anchor=(1, 0.5), loc='center left')\n",
    "    plt.tight_layout(pad=1)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **angle(s)-to-random ratios** (only if at least a pair of ratio partner runs is available)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "if (True):  # plot this?\n",
    "\n",
    "    plt.close(\"ratios\")\n",
    "\n",
    "    # case 0: no numerator-denominator pairs at all\n",
    "    if len([s for s in dataNamesToOpen if ((s in lsRatios) & ((lsRatios[s] in dataNamesToOpen) if s in lsRatios else False))])==0:\n",
    "        print(\"no ratio partner runs available --> no angle-to-random plot at all\")\n",
    "\n",
    "    # case 1: at least 1 numerator-denominator pair available\n",
    "    else:\n",
    "        fig, ax = plt.subplots(num=\"ratios\", figsize=[5, 4], sharex=True)\n",
    "        for s in [s for s in dataNamesToOpen if s in lsRatios]:  # loop on runs which have denominator partners\n",
    "            if not (lsRatios[s] in dataNamesToOpen):\n",
    "                print(\"%s (denominator for %s) not loaded --> load it to plot the ratio\" % (lsRatios[s], s))\n",
    "\n",
    "            else:\n",
    "                xPhysN, yPhysN = rebin(dictData[s][0], dictData[s][1], nBin)\n",
    "                xBkgN, yBkgN = rebin(dictData[dataNames[s]][0], dictData[dataNames[s]][1], nBin)\n",
    "                xPhysD, yPhysD = rebin(dictData[lsRatios[s]][0], dictData[lsRatios[s]][1], nBin)\n",
    "                xBkgD, yBkgD = rebin(dictData[dataNames[lsRatios[s]]][0], dictData[dataNames[lsRatios[s]]][1], nBin)\n",
    "                xN, yN, xD, yD = xPhysN, yPhysN-yBkgN, xPhysN, yPhysD-yBkgD\n",
    "                ax.plot(xD, yN / yD, label=\"%s / \\n%s\" % (s, lsRatios[s]))\n",
    "\n",
    "        ax.grid()\n",
    "        ax.set_ylabel(\"angle-to-random photon spectra ratio\")\n",
    "        ax.set_xlabel(\"energy [MeV]\")\n",
    "        fig.legend(bbox_to_anchor=(1, 0.5), loc='center left', fontsize=9.5)\n",
    "        plt.tight_layout()\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **fancy plots for papers**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *set boolean for the whole section and automatic labelling system entries in the first 2 cells below:*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bFancy = False  # plot all this fancy stuff?\n",
    "# note that single plots below are subject to individual booleans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####\n",
    "# also, automatic labelling system...\n",
    "\n",
    "if bFancy:\n",
    "    \n",
    "    dictLabels = {}\n",
    "    for s in dataNamesToOpen:\n",
    "        dictLabels.update({s: {}})\n",
    "        dictLabels[s].update({\"crys\": \"\", \"config\": \"\", \"energy\": \"\"})\n",
    "\n",
    "        # 2021 W 1mm 111 by Matek+IJCLab\n",
    "        if ((\"data21\" in s) & (\"IJC\" in s)):\n",
    "            if \"run\" in s:\n",
    "                dictLabels[s][\"crys\"] = r\"W $\\langle 111 \\rangle$ before irradiation\"  # could add thickness: $1~\\mathrm{mm}$ thick\n",
    "            if \"06_11_2021\" in s:\n",
    "                dictLabels[s][\"crys\"] = r\"W $\\langle 111 \\rangle$ after irradiation\"  # could add thickness: $1~\\mathrm{mm}$ thick\n",
    "            dictLabels[s][\"energy\"] = \"855 MeV\"\n",
    "            if \"AXIS\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"Amorph\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "\n",
    "        # 2019 BGO 0.1mm 111 \n",
    "        if \"data19/BGO111\" in s:\n",
    "            dictLabels[s][\"crys\"] = r\"BGO $\\langle 111 \\rangle$\"  # could add thickness: $0.1~\\mathrm{mm}$ thick\n",
    "            if \"axis\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"amorphous\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "            if \"855\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"855 MeV\"\n",
    "            if \"600\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"600 MeV\"\n",
    "\n",
    "        # 2019 YAG 0.5mm 100 \n",
    "        if (\"data19/YAG-\" in s) | (\"data19/Yag-\" in s):\n",
    "            dictLabels[s][\"crys\"] = r\"YAG undoped $\\langle 100 \\rangle$\"  # could add thickness: $0.5~\\mathrm{mm}$ thick\n",
    "            if \"axis\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"amorphous\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "            if \"855\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"855 MeV\"\n",
    "            if \"600\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"600 MeV\"\n",
    "\n",
    "        # 2019 YAG:Ce 0.5mm 100 \n",
    "        if (\"data19/YAGCE\" in s) | (\"data19/YagCe\" in s):\n",
    "            dictLabels[s][\"crys\"] = r\"YAG:Ce $\\langle 100 \\rangle$\"  # could add thickness: $0.5~\\mathrm{mm}$ thick\n",
    "            if \"axis\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"amorphous\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "            if \"600\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"600 MeV\"\n",
    "            else:\n",
    "                dictLabels[s][\"energy\"] = \"855 MeV\"\n",
    "\n",
    "        # 2019 YAG 0.1mm 100 \n",
    "        if \"data19/YAGthin\" in s:\n",
    "            dictLabels[s][\"crys\"] = r\"YAG undoped $\\langle 100 \\rangle$\"  # could add thickness: $0.1~\\mathrm{mm}$ thick\n",
    "            if \"axis\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"amorphous\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "            if \"855\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"855 MeV\"\n",
    "            if \"600\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"600 MeV\"\n",
    "\n",
    "        # 2019 PWO 0.5mm (but tilted) 021\n",
    "        if \"data19/PWO-021\" in s:\n",
    "            dictLabels[s][\"crys\"] = r\"PWO $\\langle 021 \\rangle$\"  # could add thickness: $0.5~\\mathrm{mm}$ thick\n",
    "            if \"axis\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"amorphous\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "            if \"855\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"855 MeV\"\n",
    "            if \"600\" in s:\n",
    "                dictLabels[s][\"energy\"] = \"600 MeV\"\n",
    "\n",
    "        # 2019 PWO 0.5mm (but tilted) 110\n",
    "        if \"data19/PWOa\" in s:\n",
    "            dictLabels[s][\"crys\"] = r\"PWO $\\langle 110 \\rangle$\"  # could add thickness: $0.5~\\mathrm{mm}$ thick\n",
    "            if \"axis\" in s:\n",
    "                dictLabels[s][\"config\"] = \"on-axis\"\n",
    "            if \"amorphous\" in s:\n",
    "                dictLabels[s][\"config\"] = \"random\"\n",
    "            dictLabels[s][\"energy\"] = \"855 MeV\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "if (bFancy & (True)):  # also, plot this in particular?\n",
    "\n",
    "    plt.close(\"spectrum_final_fancy\")\n",
    "    fig, ax = plt.subplots(num=\"spectrum_final_fancy\", figsize=[6.5, 4.5], sharex=True)\n",
    "\n",
    "    for s in dataNamesToOpen:\n",
    "        xPhys, yPhys = rebin(dictData[s][0], dictData[s][1], nBin)\n",
    "        xBkg, yBkg = rebin(dictData[dataNames[s]][0], dictData[dataNames[s]][1], nBin)\n",
    "        x, y = xPhys, (yPhys-yBkg)/abs(xPhys[1]-xPhys[0])\n",
    "        ave = sum(x*y) / sum(y)\n",
    "        label = None if dictLabels[s][\"crys\"]==\"\" else (dictLabels[s][\"config\"] if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))==1 else dictLabels[s][\"crys\"]+\", \"+dictLabels[s][\"energy\"]+\", \"+dictLabels[s][\"config\"])\n",
    "        ax.plot(x, y, label=label)\n",
    "                \n",
    "    ax.text(.5, .95,  dictLabels[s][\"crys\"]+\"\\nbeam energy \"+dictLabels[s][\"energy\"] if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))==1 else \"\", horizontalalignment='center', verticalalignment='top', transform=ax.transAxes, fontsize=12)\n",
    "    ax.set_yscale(\"log\")\n",
    "    ax.grid()\n",
    "    ax.set_ylabel(r\"$\\frac{\\mathrm{d}N}{T \\Delta  E \\mathrm{d}E}~\\left[ \\frac{1}{\\mathrm{MeV} \\cdot \\mathrm{s}} \\right]$\", fontsize=20)\n",
    "    ax.set_xlabel(r\"radiation energy $E$ [MeV]\", fontsize=12)\n",
    "    fig.legend(loc='upper right', fancybox=False, framealpha=1, edgecolor=\"k\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "if (bFancy & (True)):  # also, plot this in particular?\n",
    "\n",
    "    plt.close(\"spectralIntensity_final_fancy\")\n",
    "    fig, ax = plt.subplots(num=\"spectralIntensity_final_fancy\", figsize=[6.5, 4.5], sharex=True)\n",
    "\n",
    "    for s in dataNamesToOpen:\n",
    "        xPhys, yPhys = rebin(dictData[s][0], dictData[s][1], nBin)\n",
    "        xBkg, yBkg = rebin(dictData[dataNames[s]][0], dictData[dataNames[s]][1], nBin)\n",
    "        x, y = xPhys, (yPhys-yBkg)/abs(xPhys[1]-xPhys[0])\n",
    "        ave = sum(x*y) / sum(y)\n",
    "        label = None if dictLabels[s][\"crys\"]==\"\" else (dictLabels[s][\"config\"] if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))==1 else dictLabels[s][\"crys\"]+\", \"+dictLabels[s][\"energy\"]+\", \"+dictLabels[s][\"config\"])\n",
    "        ax.plot(x, y*x, label=label)\n",
    "        \n",
    "    ax.text(.5, .95,  dictLabels[s][\"crys\"]+\"\\nbeam energy \"+dictLabels[s][\"energy\"] if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))==1 else \"\", horizontalalignment='center', verticalalignment='top', transform=ax.transAxes, fontsize=12)\n",
    "    ax.grid()\n",
    "    ax.set_ylabel(r\"$\\frac{E \\mathrm{d}N}{T \\Delta  E \\mathrm{d}E}~\\left[ \\frac{1}{\\mathrm{s}} \\right]$\", fontsize=20)\n",
    "    ax.set_xlabel(r\"radiation energy $E$ [MeV]\", fontsize=12)\n",
    "    fig.legend(loc='upper right', fancybox=False, framealpha=1, edgecolor=\"k\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (bFancy & (True)):  # also, plot this in particular?\n",
    "    \n",
    "    plt.close(\"ratio_final_fancy\")\n",
    "    \n",
    "    # case 0: no numerator-denominator pairs at all\n",
    "    if len([s for s in dataNamesToOpen if ((s in lsRatios) & ((lsRatios[s] in dataNamesToOpen) if s in lsRatios else False))])==0:\n",
    "        print(\"no ratio partner runs available --> no angle-to-random plot at all\")\n",
    "\n",
    "    # case 1: at least 1 numerator-denominator pair available\n",
    "    else:\n",
    "        fig, ax = plt.subplots(num=\"ratio_final_fancy\", figsize=[6, 4.5], sharex=True)\n",
    "        for i, s in enumerate([s for s in dataNamesToOpen if s in lsRatios]):  # loop on runs which have denominator partners\n",
    "            if not (lsRatios[s] in dataNamesToOpen):\n",
    "                print(\"%s (denominator for %s) not loaded --> load it to plot the ratio\" % (lsRatios[s], s))\n",
    "\n",
    "            else:\n",
    "                xPhysN, yPhysN = rebin(dictData[s][0], dictData[s][1], nBin)\n",
    "                xBkgN, yBkgN = rebin(dictData[dataNames[s]][0], dictData[dataNames[s]][1], nBin)\n",
    "                xPhysD, yPhysD = rebin(dictData[lsRatios[s]][0], dictData[lsRatios[s]][1], nBin)\n",
    "                xBkgD, yBkgD = rebin(dictData[dataNames[lsRatios[s]]][0], dictData[dataNames[lsRatios[s]]][1], nBin)\n",
    "                xN, yN, xD, yD = xPhysN, yPhysN-yBkgN, xPhysN, yPhysD-yBkgD\n",
    "                ax.plot(xD, yN / yD, label=None if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))==1 else dictLabels[s][\"crys\"]+\", \"+dictLabels[s][\"energy\"])\n",
    "                if i==0:\n",
    "                    ax.grid()\n",
    "                    ax.set_ylabel(\"axis-to-random ratio\", fontsize=12)\n",
    "                    ax.set_xlabel(r\"radiation energy $E$ [MeV]\", fontsize=12)\n",
    "                if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))>1:\n",
    "                    plt.legend(fancybox=False, framealpha=1, edgecolor=\"k\")\n",
    "                ax.text(.5, .25,  dictLabels[s][\"crys\"]+\"\\nbeam energy \"+dictLabels[s][\"energy\"] if len(set([dictLabels[ss][\"crys\"]+dictLabels[ss][\"energy\"] for ss in dictLabels.keys() if not (dictLabels[ss][\"crys\"]==\"\")]))==1 else \"\", horizontalalignment='center', verticalalignment='top', transform=ax.transAxes, fontsize=12)\n",
    "                \n",
    "        plt.tight_layout()\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WHITEBOARD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}