diff --git a/.gitignore b/.gitignore index 9044f01..31f0f32 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,8 @@ data/* .ipynb_checkpoints/ notebooks/root_tests visualization/SiphraReader/build +__pycache__/ +*.pyc # Do not ignore the sample_data. Track only .dat and .json files !data/test_data data/test_data/* diff --git a/Output.root b/Output.root deleted file mode 100644 index c0b8213..0000000 Binary files a/Output.root and /dev/null differ diff --git a/analysis/__init__.py b/analysis/__init__.py index b5e352a..f516f04 100644 --- a/analysis/__init__.py +++ b/analysis/__init__.py @@ -1,3 +1,4 @@ from .fit import * +from .calibration import * -__all__ = ['gauspeak', 'bg_exp', 'peak_and_bg', 'fit_peak_expbg'] \ No newline at end of file +__all__ = ['gauspeak', 'bg_exp', 'peak_and_bg', 'fit_peak_expbg', 'calibration_fit', 'calibrated_histogram', 'calibrated_acquisition' ] \ No newline at end of file diff --git a/analysis/__pycache__/__init__.cpython-314.pyc b/analysis/__pycache__/__init__.cpython-314.pyc deleted file mode 100644 index 8555aba..0000000 Binary files a/analysis/__pycache__/__init__.cpython-314.pyc and /dev/null differ diff --git a/analysis/__pycache__/fit.cpython-314.pyc b/analysis/__pycache__/fit.cpython-314.pyc deleted file mode 100644 index f865946..0000000 Binary files a/analysis/__pycache__/fit.cpython-314.pyc and /dev/null differ diff --git a/analysis/calibration.py b/analysis/calibration.py new file mode 100644 index 0000000..c680b9a --- /dev/null +++ b/analysis/calibration.py @@ -0,0 +1,66 @@ +import ROOT +import numpy as np + +def calibration_fit(histogram, energy_ranges, energies): + """Function to create a linear calibration fit based on a histogram. Returns slope and constant of linear fit. + Inputs: Histogram to base calibration on. + Ranges within which the peaks of the histogram are located, input as a list if tuples. AT LEAST 2 points. + Known energies of these peaks, in MeV.""" + + channels = [] + for i in range(len(energy_ranges)): + cal_fit=ROOT.TF1("cal_fit_" + str(i), "gaus", energy_ranges[i][0], energy_ranges[i][1]) + histogram.Fit(cal_fit, "R+S") + channels.append(cal_fit.GetParameter('Mean')) + + channels = np.array(channels, dtype='float64') + energies = np.array(energies, dtype='float64') + + graph = ROOT.TGraph(len(channels), channels, energies) + fit = ROOT.TF1("calib", "pol1", min(channels), max(channels)) + graph.Fit(fit) + + a = fit.GetParameter(1) + b = fit.GetParameter(0) + + return [a, b] + + +def calibrated_histogram(linear_fit, acquisition, n_of_bins): + a = linear_fit[0] + b = linear_fit[1] + data_cal = a * (acquisition['s']/len(acquisition.active_chs)) + b + + emax = a * n_of_bins + b + emin = b + + hist_cal = ROOT.TH1F("h_cal", "Calibrated Spectrum", n_of_bins, emin, emax) + hist_cal.Fill(data_cal) + + return(hist_cal) + + +def calibrated_acquisition(linear_fit, acquisition, n_of_bins): + a = linear_fit[0] + b = linear_fit[1] + data_cal = a * (acquisition['s']/len(acquisition.active_chs)) + b + + emax = a * n_of_bins + b + emin = b + + return(data_cal, emin, emax) + + +def energy_resolution(hist, peak_ranges, peak_energies): + """Calculates energy resolution. Input peak_ranges as a list of tuples. + Outputs a list of resolutions for the different energies.""" + resolutions = [] + for i in range (len(peak_ranges)): + resolution_fit = ROOT.TF1("res_fit_" + str(i), "gaus", peak_ranges[i][0], peak_ranges[i][1]) + hist.Fit(resolution_fit, "R+S") + + std_dev = resolution_fit.GetParameter('Sigma') + resolution = (2.355 * std_dev)/peak_energies[i] + resolutions.append(resolution) + + return resolutions \ No newline at end of file diff --git a/file_converters/__pycache__/d2a_decoder.cpython-314.pyc b/file_converters/__pycache__/d2a_decoder.cpython-314.pyc deleted file mode 100644 index 2772b04..0000000 Binary files a/file_converters/__pycache__/d2a_decoder.cpython-314.pyc and /dev/null differ diff --git a/noteboks/baseline_subtraction.ipynb b/noteboks/baseline_subtraction.ipynb index ced9272..53e5e9f 100644 --- a/noteboks/baseline_subtraction.ipynb +++ b/noteboks/baseline_subtraction.ipynb @@ -14,7 +14,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 5, "id": "initial_id", "metadata": { "ExecuteTime": { @@ -47,7 +47,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 6, "id": "7aa83758-b8f5-4a20-96c3-c8381ea71536", "metadata": {}, "outputs": [], @@ -64,7 +64,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 7, "id": "68a9cd08dc972f65", "metadata": { "ExecuteTime": { @@ -85,7 +85,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 8, "id": "cae21a3ee684d69f", "metadata": { "ExecuteTime": { @@ -97,26 +97,46 @@ } }, "outputs": [ + { + "ename": "TypeError", + "evalue": "none of the 2 overloaded methods succeeded. Full details:\n void TH1::FillN(Int_t, const Double_t*, const Double_t*, const Double_t*, Int_t) =>\n TypeError: takes at least 5 arguments (3 given)\n could not convert argument to buffer or nullptr", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mTypeError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[8]\u001b[39m\u001b[32m, line 24\u001b[39m\n\u001b[32m 22\u001b[39m hist_singlech = ROOT.TH1F(\u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mCh. \u001b[39m\u001b[38;5;132;01m{\u001b[39;00macq.active_chs[\u001b[32m0\u001b[39m]\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mname\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m\"\u001b[39m, \u001b[33m\"\u001b[39m\u001b[33m\"\u001b[39m, BITS9, \u001b[32m0\u001b[39m, BITS12)\n\u001b[32m 23\u001b[39m hist.Fill(acq[ch]/\u001b[38;5;28mlen\u001b[39m(acq.active_chs))\n\u001b[32m---> \u001b[39m\u001b[32m24\u001b[39m \u001b[43mhist_singlech\u001b[49m\u001b[43m.\u001b[49m\u001b[43mFill\u001b[49m\u001b[43m(\u001b[49m\u001b[43macq\u001b[49m\u001b[43m[\u001b[49m\u001b[43macq\u001b[49m\u001b[43m.\u001b[49m\u001b[43mactive_chs\u001b[49m\u001b[43m[\u001b[49m\u001b[32;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 25\u001b[39m hist.Scale(\u001b[32m1\u001b[39m/acq.exposure)\n\u001b[32m 26\u001b[39m hist_singlech.Scale(\u001b[32m1\u001b[39m/acq.exposure)\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/COMCUBE/lib/python3.14/site-packages/ROOT/_pythonization/_th1.py:223\u001b[39m, in \u001b[36m_FillWithNumpyArray\u001b[39m\u001b[34m(self, *args)\u001b[39m\n\u001b[32m 219\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(weights) != \u001b[38;5;28mlen\u001b[39m(data):\n\u001b[32m 220\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[32m 221\u001b[39m \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mLength mismatch: data length (\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mlen\u001b[39m(data)\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m) != weights length (\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mlen\u001b[39m(weights)\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m)\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 222\u001b[39m )\n\u001b[32m--> \u001b[39m\u001b[32m223\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mFillN\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mlen\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mweights\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 224\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m 225\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m._Fill(*args)\n", + "\u001b[31mTypeError\u001b[39m: none of the 2 overloaded methods succeeded. Full details:\n void TH1::FillN(Int_t, const Double_t*, const Double_t*, const Double_t*, Int_t) =>\n TypeError: takes at least 5 arguments (3 given)\n could not convert argument to buffer or nullptr" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: Ch. 14 BL subt. (Potential memory leak).\n", + "Warning in : Replacing existing TH1: Ch. 2 No correction (Potential memory leak).\n" + ] + }, { "data": { "text/html": [ "\n", "\n", - "
\n", + "
\n", "
\n", "\n", "\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "if ROOT.gROOT.FindObject('cv_4plots'):\n", " ROOT.gROOT.FindObject('cv_4plots').Close()\n", @@ -364,36 +317,14 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "46bd074b-d5aa-48a4-b027-f966314995cd", "metadata": { "jupyter": { "source_hidden": true } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "ADC channels range:\n", - "--------------------------------------------------\n", - "Channel\t\tBL-subtracted\tNo correction\n", - "--------------------------------------------------\n", - "Summed\t\t0.0 - 3965.0\t935.5 - 4880.5\n", - "Ch. 2\t\t0.0 - 3951.0\t51.0 - 4095.0\n", - "Ch. 14\t\t0.0 - 3979.0\t0.0 - 4095.0\n", - "\n", - "\n", - "\n", - "Number of values above single channel range in 'Summed' channel:\n", - "------------------------------\n", - "BL-subtracted\tNo correction\n", - "------------------------------\n", - " 0 \t 22 \n" - ] - } - ], + "outputs": [], "source": [ "print(\"ADC channels range:\\n\"\n", " f\"{\"\":->50}\\n\"\n", @@ -441,7 +372,7 @@ }, { "cell_type": "code", - "execution_count": 76, + "execution_count": null, "id": "a4f39db9-ad0d-4f5b-b3b5-ea67040593a6", "metadata": {}, "outputs": [], @@ -464,36 +395,14 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "0b44ed9c-6f45-4174-8f81-1407bcdacc05", "metadata": { "jupyter": { "source_hidden": true } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "ADC channels range after baseline subtraction:\n", - "CHANNEL 14\n", - "------------------------------\n", - "Data set\tRange\n", - "------------------------------\n", - "Background\t0.0 - 3924.0\n", - "Source\t\t0.0 - 3797.0\n", - "\n", - "SUMMED CHANNEL\n", - "------------------------------\n", - "Data set\tRange\n", - "------------------------------\n", - "Background\t0.0 - 3924.0\n", - "Source\t\t0.0 - 3797.0\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "print(\"ADC channels range after baseline subtraction:\\n\"\n", " \"CHANNEL 14\\n\"\n", @@ -520,24 +429,14 @@ }, { "cell_type": "code", - "execution_count": 43, + "execution_count": null, "id": "2093b302-1f12-4e19-86f6-74506d4e5ff6", "metadata": { - "collapsed": true, "jupyter": { - "outputs_hidden": true + "source_hidden": true } }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning in : Replacing existing TH1: Background (Potential memory leak).\n", - "Warning in : Replacing existing TH1: Signal 137Cs (Potential memory leak).\n" - ] - } - ], + "outputs": [], "source": [ "# Raw histograms\n", "hist_bg = ROOT.TH1F(\"Background\", \"\", BITS9, 0, BITS12)\n", @@ -557,84 +456,10 @@ }, { "cell_type": "code", - "execution_count": 62, + "execution_count": null, "id": "eacf6495-684b-42fa-9eab-8ca73fb2afef", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Error in : Cannot set Y axis to log scale\n" - ] - }, - { - "data": { - "text/html": [ - "\n", - "\n", - "
\n", - "
\n", - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "# Plot all the histograms\n", "if ROOT.gROOT.FindObject('cv1'):\n", @@ -673,119 +498,14 @@ }, { "cell_type": "code", - "execution_count": 72, + "execution_count": null, "id": "1b94d6b2-ea7d-41ea-acb2-a7e820e9dd34", "metadata": { "jupyter": { "source_hidden": true } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "****************************************\n", - "Minimizer is Minuit2 / Migrad\n", - "MinFCN = 3.07992\n", - "Chi2 = 6.15985\n", - "NDf = 45\n", - "Edm = 4.33544e-06\n", - "NCalls = 378\n", - "Const = 967.676 +/- 315.198 \n", - "Denom = 204.714 +/- 16.8531 \n", - "Norm = 26.0236 +/- 2.92104 \n", - "Mean = 748.998 +/- 10.0065 \n", - "Sigma = 69.1169 +/- 10.4622 \n", - "****************************************\n", - "Minimizer is Minuit2 / Migrad\n", - "MinFCN = 4.22284\n", - "Chi2 = 8.44567\n", - "NDf = 45\n", - "Edm = 1.92878e-06\n", - "NCalls = 607\n", - "Const = 65.7641 +/- 57.5374 \n", - "Denom = 441.486 +/- 153.51 \n", - "Norm = 1.79752 +/- 0.994549 \n", - "Mean = 1240.31 +/- 25.9607 \n", - "Sigma = 44.4995 +/- 28.325 \n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning in : Deleting canvas with same name: cv_fit\n", - "Error in : Cannot set Y axis to log scale\n" - ] - }, - { - "data": { - "text/html": [ - "\n", - "\n", - "
\n", - "
\n", - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "# Fit background peaks\n", "xrange_k40 = (600,1000)\n", @@ -809,25 +529,10 @@ }, { "cell_type": "code", - "execution_count": 103, + "execution_count": null, "id": "3e0906bf-c280-4c5e-b555-871735d01a54", - "metadata": { - "collapsed": true, - "jupyter": { - "outputs_hidden": true, - "source_hidden": true - } - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning in : Replacing existing TH1: Calibrated background (Potential memory leak).\n", - "Warning in : Replacing existing TH1: Calibrated signal (Potential memory leak).\n" - ] - } - ], + "metadata": {}, + "outputs": [], "source": [ "# Calibrate data and get calibrated histograms\n", "K40FIT = fit_fn_k40.GetParameter('Mean')\n", @@ -859,88 +564,10 @@ }, { "cell_type": "code", - "execution_count": 101, + "execution_count": null, "id": "5416087a-578e-4d42-8157-2516db3c6b39", - "metadata": { - "jupyter": { - "source_hidden": true - } - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Error in : Cannot set Y axis to log scale\n" - ] - }, - { - "data": { - "text/html": [ - "\n", - "\n", - "
\n", - "
\n", - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "metadata": {}, + "outputs": [], "source": [ "# Plot all the histograms\n", "if ROOT.gROOT.FindObject('cv1'):\n", @@ -972,118 +599,20 @@ }, { "cell_type": "code", - "execution_count": 115, + "execution_count": null, "id": "dcbad448-c208-4054-a4ab-e2b71b7ed2a4", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "****************************************\n", - "Minimizer is Minuit2 / Migrad\n", - "MinFCN = 5628.78\n", - "Chi2 = 11257.6\n", - "NDf = 112\n", - "Edm = 1.88824e-07\n", - "NCalls = 124\n", - "Constant = 4002.29 +/- 17.4275 \n", - "Mean = 0.479418 +/- 0.000489701 \n", - "Sigma = 0.136089 +/- 0.000348374 \t (limited)\n" - ] - } - ], + "outputs": [], "source": [ "res = hist_clean_clb.Fit(\"gaus\", \"L S\", \"\", 0, 2)" ] }, { "cell_type": "code", - "execution_count": 129, + "execution_count": null, "id": "6d7f6090-75e8-4779-8a1c-35f1cad0a43d", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Error in the location of Cs peak: 27.5%\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning in : Deleting canvas with same name: cv_fit\n" - ] - }, - { - "data": { - "text/html": [ - "\n", - "\n", - "
\n", - "
\n", - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "cv_fit = ROOT.TCanvas(\"cv_fit\", \"cv_fit\", 800,600)\n", "# hist_bg_fit.SetFillStyle(0)\n", @@ -1099,18 +628,10 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "58f9c3b4-03c0-453a-9d37-cd9cd565ad5f", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "FWHM=0.32 Energy resolution at 480MeV: 66.84%\n" - ] - } - ], + "outputs": [], "source": [ "import numpy as np\n", "sigma = 0.136089\n", @@ -1121,9 +642,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python (ROOT in Conda)", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "comcube" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1135,7 +656,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.2" + "version": "3.14.3" } }, "nbformat": 4, diff --git a/noteboks/energyResolution-Copy1.ipynb b/noteboks/energyResolution-Copy1.ipynb new file mode 100644 index 0000000..d68f13a --- /dev/null +++ b/noteboks/energyResolution-Copy1.ipynb @@ -0,0 +1,893 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4a3b1bc5-ec73-4d69-a840-07d057f4e2eb", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "project_root=Path(project_root)\n", + "\n", + "from processing import SiphraAcquisition, MetadataLoader\n", + "from analysis import fit_peak_expbg\n", + "\n", + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "\n", + "# Energy emission peaks in MeV\n", + "Cs137_MeV = 0.662\n", + "Na22_MeV = [0.511, 1.275, 1.786]\n", + "Co60_MeV = [1.173, 1.332, 2,505]\n", + "Am241_MeV = 0.0595\n", + "# Background emission\n", + "K40_MeV = 1.460\n", + "Tl208_MEV = 2.614\n", + "\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "75493f51-2089-4f0a-8057-0f2c2bf7466f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[SIPHRAAcquisition(File: 'SUBT_01_EnergyResolution_thr30_gain1over100_1over3_Background'), SIPHRAAcquisition(File: 'SUBT_02_EnergyResolution_thr30_gain1over100_1over3_Cs137'), SIPHRAAcquisition(File: 'SUBT_03_EnergyResolution_thr30_gain1over100_1over3_Background2'), SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22'), SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4'), SIPHRAAcquisition(File: 'SUBT_06_EnergyResolution_thr30_gain1over100_1over3_Co60'), SIPHRAAcquisition(File: 'SUBT_07_EnergyResolution_thr30_gain1over100_1over3_Background6'), SIPHRAAcquisition(File: 'SUBT_08_EnergyResolution_thr30_gain1over100_1over3_Am241'), SIPHRAAcquisition(File: 'SUBT_09_EnergyResolution_thr30_gain1over100_1over3_Background8')]\n" + ] + } + ], + "source": [ + "files = sorted((project_root/'data/260323').glob('SUBT_*'))\n", + "acqs = [SiphraAcquisition(file) for file in files]\n", + "print(acqs)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "88f4c856-78c0-4ca6-ad58-cc10f3b53649", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cs-137\n", + "Na22\n", + "Co-60\n", + "Am-241\n" + ] + } + ], + "source": [ + "# histograms\n", + "hists = {}\n", + "sources = []\n", + "for sgnl, bg in zip(acqs[1::2], acqs[2::2]):\n", + " filepath = sgnl.filepath.stem\n", + " src = (MetadataLoader.load(sgnl.metadataFile)).source\n", + " sources.append(src)\n", + " print(src)\n", + " # Signal and Background\n", + " hist_sgnl = ROOT.TH1F(f\"{src} Signal\", \"\", BITS12, 0, BITS12)\n", + " hist_bg = ROOT.TH1F(f\"{src} Background\", \"\", BITS12, 0 , BITS12)\n", + " hist_sgnl.Fill(sgnl['s']/len(sgnl.active_chs))\n", + " hist_bg.Fill(bg['s']/len(bg.active_chs))\n", + " hist_bg.Scale(sgnl.exposure/bg.exposure)\n", + " # Clean spectrum\n", + " hist_clean = hist_sgnl.Clone(f\"{src} Clean\")\n", + " hist_clean.Add(hist_bg, -1)\n", + " for hist in [hist_sgnl, hist_bg, hist_clean]:\n", + " # hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Normalized counts\")\n", + " hists[src] = hist_sgnl\n", + " hists[f\"{src}_BG\"] = hist_bg\n", + " hists[f\"{src}_clean\"] = hist_clean\n", + " del hist_sgnl, hist_bg\n", + "#print(hists)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "36f654a6-bd64-471f-a0e3-01927db52e3e", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1600, 1200)\n", + "cv.Divide(2,2)\n", + "\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lgds = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "for idx, (src, lg) in enumerate(zip(sources, lgds)): \n", + " cv.cd(idx+1)\n", + " \n", + " sgnl = hists[src]\n", + " bg = hists[src+'_BG']\n", + " clean = hists[src+'_clean']\n", + " \n", + " lg.AddEntry(sgnl, \"Signal\")\n", + " lg.AddEntry(bg, \"Background\")\n", + " lg.AddEntry(clean, \"Subtracted\")\n", + " sgnl.SetLineColor(colors[0])\n", + " bg.SetLineColor(colors[1])\n", + " clean.SetLineColor(colors[2])\n", + " sgnl.SetTitle(src)\n", + " sgnl.Draw(\"hist\")\n", + " bg.Draw(\"sames hist\")\n", + " clean.Draw(\"sames hist\")\n", + " ROOT.gPad.Update()\n", + " for i, sp in enumerate([sgnl, bg, clean]):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + " lg.Draw()\n", + " cv.cd(idx+1).SetLogy()\n", + " cv.Draw()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6db7ba78-ddd9-4d76-8bba-bf310f8f1221", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.004280158214429103\n", + "-0.2655932438428593\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 196.518\n", + "NDf = 127\n", + "Edm = 2.29967e-07\n", + "NCalls = 65\n", + "Constant = 2459.66 +/- 7.74946 \n", + "Mean = 183.993 +/- 0.078652 \n", + "Sigma = 29.2438 +/- 0.0713573 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 356.977\n", + "NDf = 147\n", + "Edm = 1.20533e-06\n", + "NCalls = 64\n", + "Constant = 1361.97 +/- 6.38965 \n", + "Mean = 275.54 +/- 0.144272 \n", + "Sigma = 32.7382 +/- 0.141925 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 52.0578\n", + "NDf = 47\n", + "Edm = 1.53695e-06\n", + "NCalls = 107\n", + "Constant = 255.901 +/- 3.65957 \n", + "Mean = 478.856 +/- 0.861944 \n", + "Sigma = 32.2107 +/- 1.8725 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 50.431\n", + "NDf = 47\n", + "Edm = 6.26293e-06\n", + "NCalls = 110\n", + "Constant = 187.956 +/- 4.04943 \n", + "Mean = 177.384 +/- 1.67095 \n", + "Sigma = 37.5292 +/- 4.35406 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 152.171\n", + "NDf = 97\n", + "Edm = 2.80563e-08\n", + "NCalls = 96\n", + "Constant = 606.675 +/- 4.80363 \n", + "Mean = 371.93 +/- 0.452134 \n", + "Sigma = 45.3973 +/- 0.782852 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 63.6754\n", + "NDf = 67\n", + "Edm = 1.95162e-06\n", + "NCalls = 89\n", + "Constant = 880.588 +/- 5.38234 \n", + "Mean = 471.391 +/- 0.435139 \n", + "Sigma = 32.0034 +/- 0.470993 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 0.00408946\n", + "NDf = 1\n", + "p0 = -0.265593 +/- 0.109303 \n", + "p1 = 0.00428016 +/- 0.000302376 \n" + ] + } + ], + "source": [ + "#Fitting all the visible peaks. Values around peak input manually\n", + "\n", + "Cs137_fit = ROOT.TF1(\"Cs137_fit\", \"gaus\", 120,250)\n", + "\n", + "Co60a_fit = ROOT.TF1(\"Co60a_fit\", \"gaus\", 200,350)\n", + "Co60b_fit = ROOT.TF1(\"Co60b_fit\", \"gaus\", 450,500)\n", + "\n", + "Na22a_fit = ROOT.TF1(\"Na22a_fit\", \"gaus\", 150,200)\n", + "Na22b_fit = ROOT.TF1(\"Na22b_fit\", \"gaus\", 320,420)\n", + "Na22c_fit = ROOT.TF1(\"Na22c_fit\", \"gaus\", 450,520)\n", + "\n", + "fits = [Cs137_fit, Co60a_fit, Co60b_fit, Na22a_fit, Na22b_fit, Na22c_fit]\n", + "\n", + "#Fit to histogram and find mean value of peak\n", + "hists['Cs-137_clean'].Fit(Cs137_fit, \"R+S\")\n", + "hists['Co-60_clean'].Fit(Co60a_fit, \"R+S\")\n", + "hists['Co-60_clean'].Fit(Co60b_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22a_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22b_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22c_fit, \"R+S\")\n", + "\n", + "Cs137_mean = Cs137_fit.GetParameter('Mean')\n", + "Co60a_mean = Co60a_fit.GetParameter('Mean')\n", + "Co60b_mean = Co60b_fit.GetParameter('Mean')\n", + "Na22a_mean = Na22a_fit.GetParameter('Mean')\n", + "Na22b_mean = Na22b_fit.GetParameter('Mean')\n", + "Na22c_mean = Na22c_fit.GetParameter('Mean')\n", + "\n", + "channels = [Na22a_mean, Na22b_mean, Na22c_mean]\n", + "energies = Na22_MeV\n", + "channels = np.array(channels, dtype='float64')\n", + "energies = np.array(energies, dtype='float64')\n", + "\n", + "g = ROOT.TGraph(len(channels), channels, energies)\n", + "\n", + "f = ROOT.TF1(\"calib\", \"pol1\", min(channels), max(channels))\n", + "g.Fit(f)\n", + "\n", + "a = f.GetParameter(1)\n", + "b = f.GetParameter(0)\n", + "print(a)\n", + "print(b)\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()\n", + "\n", + "#lg.Draw()\n", + "#cv.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "6b5a5060-fe8f-48ba-af47-c12bca5ea672", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-0.265593239861784\n", + "17.265934723735494\n", + "17.265934723735494\n", + "-0.265593239861784\n", + "SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22')\n", + "SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4')\n", + "[0.96495224 0.7338237 0.70760773 ... 0.76057469 0.40532156 0.49787998]\n", + "16.08942624009443\n", + "-0.11364762396634862\n", + "4096\n", + "4096\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n", + "Warning in : Replacing existing TH1: Calibrated background (Potential memory leak).\n", + "Warning in : Replacing existing TH1: Calibrated signal (Potential memory leak).\n" + ] + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv_cal'):\n", + " ROOT.gROOT.FindObject('cv_cal').Close()\n", + "\n", + "\n", + " \n", + "hist = hists['Na22_clean']\n", + "\n", + "emin = a * hist.GetXaxis().GetXmin() + b\n", + "emax = a * hist.GetXaxis().GetXmax() + b\n", + "\n", + "emax_2 = a * BITS12 + b\n", + "emin_2 = a * 0 + b\n", + "\n", + "print(emin)\n", + "print(emin_2)\n", + "print(emax)\n", + "print(emax_2)\n", + "\n", + "print(acqs[3])\n", + "print(acqs[4])\n", + "\n", + "data_bg_clb = a * (acqs[4]['s']/len(acqs[4].active_chs)) + b\n", + "data_src_clb = a * (acqs[3]['s']/len(acqs[4].active_chs)) + b\n", + "\n", + "print(data_bg_clb)\n", + "\n", + "print(data_bg_clb.max())\n", + "print(data_bg_clb.min())\n", + "\n", + "h_cal = ROOT.TH1F(\"h_cal\", \"Calibrated Spectrum\",\n", + " hist.GetNbinsX(), emin, emax)\n", + "\n", + "print(hist.GetNbinsX())\n", + "print(BITS12)\n", + "\n", + "#for i in range(1, hist.GetNbinsX() + 1):\n", + "# counts = hist.GetBinContent(i)\n", + "# h_cal.SetBinContent(i, counts)\n", + "\n", + "hist_bg_clb = ROOT.TH1F(\"Calibrated background\", \"Calibrated Histograms\", BITS12, emin, emax)\n", + "hist_src_clb = ROOT.TH1F(\"Calibrated signal\", \"Calibrated Histograms\", BITS12, emin, emax)\n", + "\n", + "\n", + "hist_bg_clb.Fill(data_bg_clb)\n", + "#hist_bg_clb.Scale(NORMFACTOR)\n", + "\n", + "hist_src_clb.Fill(data_src_clb)\n", + "\n", + "hist_clean_clb = hist_src_clb.Clone(\"Calibrated signal no BG\")\n", + "hist_clean_clb.Add(hist_bg_clb, -1)\n", + "\n", + "#for hist, clr in zip([hist_bg_clb, hist_src_clb, hist_clean_clb], colors):\n", + "# hist.GetXaxis().SetTitle(\"Energy (MeV)\")\n", + " # hist.GetYaxis().SetTitle(\"Normalized counts\")\n", + " # hist.SetLineColor(clr) \n", + "#\n", + "#cv_cal.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "57b1d3c2-9d5b-43c3-8b49-fd081e3fc77f", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Error in : Cannot set Y axis to log scale\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot all the histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "hist_bg_clb.SetLineColor(colors[0])\n", + "hist_src_clb.SetLineColor(colors[1])\n", + "lg1.AddEntry(hist_bg_clb, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_src_clb, r\"Signal ^{137}Cs\", \"l\")\n", + "hist_bg_clb.Draw(\"hist\")\n", + "hist_src_clb.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "hist_bg_clb.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "# lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_clean_clb.SetLineColor(colors[2])\n", + "hist_clean_clb.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "79a219b6-cf1a-419e-825f-272accfa28eb", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 137.425\n", + "NDf = 67\n", + "Edm = 2.14133e-06\n", + "NCalls = 84\n", + "Constant = 190.012 +/- 2.40871 \n", + "Mean = 0.492756 +/- 0.00341755 \n", + "Sigma = 0.142388 +/- 0.00438193 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 465.001\n", + "NDf = 91\n", + "Edm = 1.22703e-07\n", + "NCalls = 75\n", + "Constant = 611.78 +/- 3.97344 \n", + "Mean = 1.31966 +/- 0.00146494 \n", + "Sigma = 0.178741 +/- 0.0023569 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 742.699\n", + "NDf = 90\n", + "Edm = 7.47939e-06\n", + "NCalls = 95\n", + "Constant = 827.73 +/- 4.69558 \n", + "Mean = 1.74335 +/- 0.00131288 \n", + "Sigma = 0.172477 +/- 0.00195917 \t (limited)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Error in : Cannot set Y axis to log scale\n" + ] + } + ], + "source": [ + "# Plot all the histograms\n", + "\n", + "\n", + "\n", + "\n", + "Na22a_fit_cal = ROOT.TF1(\"Na22a_fit_cal\", \"gaus\", 0.3, 0.6)\n", + "Na22b_fit_cal = ROOT.TF1(\"Na22b_fit_cal\", \"gaus\", 1.1 , 1.5)\n", + "Na22c_fit_cal = ROOT.TF1(\"Na22c_fit_cal\", \"gaus\", 1.5, 1.9)\n", + "hist_clean_clb.Fit(Na22a_fit_cal, \"R+S\")\n", + "hist_clean_clb.Fit(Na22b_fit_cal, \"R+S\")\n", + "hist_clean_clb.Fit(Na22c_fit_cal, \"R+S\")\n", + "\n", + "hist_clean_clb.Draw(\"hist\")\n", + "cv1.cd(idx+1).SetLogy()\n", + "\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "8947e352-94db-4b13-8666-c1f2d7eca282", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.0, 0.0]\n", + "0.0 0.0\n", + "0.0 0.0\n", + "SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4')\n", + "(array([0., 0., 0., ..., 0., 0., 0.], shape=(46706,)), 0.0, 0.0)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 50.431\n", + "NDf = 47\n", + "Edm = 6.26293e-06\n", + "NCalls = 110\n", + "Constant = 187.956 +/- 4.04943 \n", + "Mean = 177.384 +/- 1.67095 \n", + "Sigma = 37.5292 +/- 4.35406 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 152.171\n", + "NDf = 97\n", + "Edm = 2.80563e-08\n", + "NCalls = 96\n", + "Constant = 606.675 +/- 4.80363 \n", + "Mean = 371.93 +/- 0.452134 \n", + "Sigma = 45.3973 +/- 0.782852 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 63.6754\n", + "NDf = 67\n", + "Edm = 1.95162e-06\n", + "NCalls = 89\n", + "Constant = 880.588 +/- 5.38234 \n", + "Mean = 471.391 +/- 0.435139 \n", + "Sigma = 32.0034 +/- 0.470993 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 0.00408946\n", + "NDf = 1\n", + "p0 = -0.265593 +/- 0.109303 \n", + "p1 = 0.00428016 +/- 0.000302376 \n" + ] + } + ], + "source": [ + "def calibration_fit(histogram, energy_ranges, energies):\n", + " \"\"\"Function to create a linear calibration fit based on a histogram. Returns slope and constant of linear fit. \n", + " Inputs: Histogram to base calibration on.\n", + " Ranges within which the peaks of the histogram are located.\n", + " Known energies of these peaks, in MeV.\"\"\"\n", + " \n", + " channels = []\n", + " for i in range(len(energy_ranges)):\n", + " cal_fit=ROOT.TF1(\"cal_fit_\" + str(i), \"gaus\", energy_ranges[i][0], energy_ranges[i][1])\n", + " hist.Fit(cal_fit, \"R+S\")\n", + " channels.append(cal_fit.GetParameter('Mean'))\n", + "\n", + " channels = np.array(channels, dtype='float64')\n", + " energies = np.array(energies, dtype='float64')\n", + "\n", + " graph = ROOT.TGraph(len(channels), channels, energies)\n", + " fit = ROOT.TF1(\"calib\", \"pol1\", min(channels), max(channels))\n", + " graph.Fit(f)\n", + " \n", + " a = fit.GetParameter(1)\n", + " b = fit.GetParameter(0)\n", + "\n", + " return [a, b]\n", + "\n", + "def calibrated_histogram(linear_fit, acquisition, n_of_bins):\n", + " a = linear_fit[0]\n", + " b = linear_fit[1]\n", + " data_cal = a * (acquisition['s']/len(acquisition.active_chs)) + b\n", + "\n", + " emax = a * n_of_bins + b\n", + " emin = b\n", + "\n", + " hist_cal = ROOT.TH1F(\"h_cal\", \"Calibrated Spectrum\", n_of_bins, emin, emax)\n", + " hist_cal.Fill(data_cal)\n", + "\n", + " return(hist_cal)\n", + "\n", + "def calibrated_acquisition(linear_fit, acquisition):\n", + " a = linear_fit[0]\n", + " b = linear_fit[1]\n", + " print(a, b)\n", + " data_cal = a * (acquisition['s']/len(acquisition.active_chs)) + b\n", + "\n", + " emax = a * BITS12 + b\n", + " emin = b\n", + "\n", + " return(data_cal, emin, emax)\n", + "\n", + " \n", + " \n", + "hist = hists['Na22_clean']\n", + "energy_ranges = [(150, 200), (320, 420), (450, 520)]\n", + "energies = Na22_MeV\n", + "\n", + "c_fit = calibration_fit(hist, energy_ranges,energies)\n", + "hist_cal_bg = calibrated_histogram(c_fit, acqs[4])\n", + "hist_cal_Na22 = calibrated_histogram(c_fit, acqs[3])\n", + "hist_cal_clean = hist_cal_Na22.Clone(\"Calibrated signal no BG\")\n", + "\n", + "#hist_cal = ROOT.TH1F(\"h_cal\", \"Calibrated Spectrum\", BITS12, emin, emax)\n", + "#hist_cal.Fill(data_cal)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "a9b090a7-1c19-4ea1-8676-176a371a3f62", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot all the histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "calibration_hist.SetLineColor(colors[0])\n", + "hist_src_clb.SetLineColor(colors[1])\n", + "lg1.AddEntry(calibration_hist, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_src_clb, r\"Signal ^{137}Cs\", \"l\")\n", + "calibration_hist.Draw(\"hist\")\n", + "hist_src_clb.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "#hist_bg_clb.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_clean_clb.SetLineColor(colors[2])\n", + "hist_clean_clb.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f60a5136-480e-48d5-b6e2-8b32ff255bbb", + "metadata": {}, + "outputs": [], + "source": [ + "def energy_Resolution(histogram, peak_range" + ] + } + ], + "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.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/noteboks/root_tests/hist010.ipynb b/noteboks/root_tests/hist010.ipynb deleted file mode 100644 index 9e8293b..0000000 --- a/noteboks/root_tests/hist010.ipynb +++ /dev/null @@ -1,142 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "id": "initial_id", - "metadata": { - "collapsed": true, - "ExecuteTime": { - "end_time": "2026-02-09T15:16:23.015658Z", - "start_time": "2026-02-09T15:16:21.512166Z" - } - }, - "source": [ - "import ROOT\n", - "import numpy as np\n", - "\n", - "c1 = ROOT.TCanvas(\"c1\",\"hists with different scales\",600,400)\n", - "\n", - "ROOT.gStyle.SetOptStat(False)\n", - "\n", - "h1 = ROOT.TH1F(\"h1\",\"my histogram\",100,-3,3)\n", - "\n", - "h1.Fill(np.array([ROOT.gRandom.Gaus(0, 1) for _ in range(10000)]))\n", - "\n", - "h1.Draw()\n", - "c1.Update()\n", - "\n", - "hint1 = ROOT.TH1F(\"hint1\",\"h1 bins integral\",100,-3,3)\n", - "\n", - "sum = 0\n", - "for i in range(1,101) :\n", - " sum += h1.GetBinContent(i)\n", - " hint1.SetBinContent(i,sum)\n", - "\n", - "rightmax = 1.1*hint1.GetMaximum()\n", - "scale = ROOT.gPad.GetUymax()/rightmax\n", - "hint1.SetLineColor(\"kRed\")\n", - "hint1.Scale(scale)\n", - "hint1.Draw(\"same\")\n", - "\n", - "axis = ROOT.TGaxis(ROOT.gPad.GetUxmax(),ROOT.gPad.GetUymin(),\n", - " ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(),0,rightmax,510,\"+L\")\n", - "axis.SetLineColor(\"kRed\")\n", - "axis.SetLabelColor(\"kRed\")\n", - "axis.Draw()" - ], - "outputs": [ - { - "data": { - "text/plain": [ - "" - ], - "text/html": [ - "\n", - "\n", - "
\n", - "
\n", - "\n", - "\n" - ] - }, - "metadata": {}, - "output_type": "display_data", - "jetTransient": { - "display_id": "7d405be2723d680b0e953166868d69b6" - } - } - ], - "execution_count": 1 - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.6" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/noteboks/root_tests/hist010_TH1_two_scales.py b/noteboks/root_tests/hist010_TH1_two_scales.py deleted file mode 100644 index bb08855..0000000 --- a/noteboks/root_tests/hist010_TH1_two_scales.py +++ /dev/null @@ -1,44 +0,0 @@ -## \file -## \ingroup tutorial_hist -## \notebook -## Example of macro illustrating how to superimpose two histograms -## with different scales in the "same" pad. -## Inspired by work of Rene Brun. -## -## \macro_image -## \macro_code -## -## \author Alberto Ferro - -import ROOT -import numpy as np - -c1 = ROOT.TCanvas("c1","hists with different scales",600,400) - -ROOT.gStyle.SetOptStat(False) - -h1 = ROOT.TH1F("h1","my histogram",100,-3,3) - -h1.Fill(np.array([ROOT.gRandom.Gaus(0, 1) for _ in range(10000)])) - -h1.Draw() -c1.Update() - -hint1 = ROOT.TH1F("hint1","h1 bins integral",100,-3,3) - -sum = 0 -for i in range(1,101) : - sum += h1.GetBinContent(i) - hint1.SetBinContent(i,sum) - -rightmax = 1.1*hint1.GetMaximum() -scale = ROOT.gPad.GetUymax()/rightmax -hint1.SetLineColor("kRed") -hint1.Scale(scale) -hint1.Draw("same") - -axis = ROOT.TGaxis(ROOT.gPad.GetUxmax(),ROOT.gPad.GetUymin(), - ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(),0,rightmax,510,"+L") -axis.SetLineColor("kRed") -axis.SetLabelColor("kRed") -axis.Draw() diff --git a/noteboks/22Na_60Co_bg-subtraction_calib.ipynb b/notebooks/22Na_60Co_bg-subtraction_calib.ipynb similarity index 100% rename from noteboks/22Na_60Co_bg-subtraction_calib.ipynb rename to notebooks/22Na_60Co_bg-subtraction_calib.ipynb diff --git a/noteboks/260204_CC_BG_comparison.ipynb b/notebooks/260204_CC_BG_comparison.ipynb similarity index 100% rename from noteboks/260204_CC_BG_comparison.ipynb rename to notebooks/260204_CC_BG_comparison.ipynb diff --git a/noteboks/260205.ipynb b/notebooks/260205.ipynb similarity index 100% rename from noteboks/260205.ipynb rename to notebooks/260205.ipynb diff --git a/noteboks/260206_CMIS_voffset.ipynb b/notebooks/260206_CMIS_voffset.ipynb similarity index 100% rename from noteboks/260206_CMIS_voffset.ipynb rename to notebooks/260206_CMIS_voffset.ipynb diff --git a/noteboks/260209_Cs137_BGsubtraction.ipynb b/notebooks/260209_Cs137_BGsubtraction.ipynb similarity index 100% rename from noteboks/260209_Cs137_BGsubtraction.ipynb rename to notebooks/260209_Cs137_BGsubtraction.ipynb diff --git a/noteboks/260212.ipynb b/notebooks/260212.ipynb similarity index 100% rename from noteboks/260212.ipynb rename to notebooks/260212.ipynb diff --git a/noteboks/CT_Cs137_calib.ipynb b/notebooks/CT_Cs137_calib.ipynb similarity index 100% rename from noteboks/CT_Cs137_calib.ipynb rename to notebooks/CT_Cs137_calib.ipynb diff --git a/notebooks/Calibration_EnergyResolution.ipynb b/notebooks/Calibration_EnergyResolution.ipynb new file mode 100644 index 0000000..e32e9af --- /dev/null +++ b/notebooks/Calibration_EnergyResolution.ipynb @@ -0,0 +1,783 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "286e580e-9b51-47e9-95b4-bbae76e7242a", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "project_root=Path(project_root)\n", + "\n", + "from processing import SiphraAcquisition, MetadataLoader\n", + "from analysis import fit_peak_expbg\n", + "\n", + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "\n", + "# Energy emission peaks in MeV\n", + "Cs137_MeV = 0.662\n", + "Na22_MeV = [0.511, 1.275, 1.786]\n", + "Co60_MeV = [1.173, 1.332, 2,505]\n", + "Am241_MeV = 0.0595\n", + "# Background emission\n", + "K40_MeV = 1.460\n", + "Tl208_MEV = 2.614\n", + "\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "bcb11ff6-e91a-4253-af2b-bc3427c5876e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[SIPHRAAcquisition(File: 'SUBT_01_EnergyResolution_thr30_gain1over100_1over3_Background'), SIPHRAAcquisition(File: 'SUBT_02_EnergyResolution_thr30_gain1over100_1over3_Cs137'), SIPHRAAcquisition(File: 'SUBT_03_EnergyResolution_thr30_gain1over100_1over3_Background2'), SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22'), SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4'), SIPHRAAcquisition(File: 'SUBT_06_EnergyResolution_thr30_gain1over100_1over3_Co60'), SIPHRAAcquisition(File: 'SUBT_07_EnergyResolution_thr30_gain1over100_1over3_Background6'), SIPHRAAcquisition(File: 'SUBT_08_EnergyResolution_thr30_gain1over100_1over3_Am241'), SIPHRAAcquisition(File: 'SUBT_09_EnergyResolution_thr30_gain1over100_1over3_Background8')]\n" + ] + } + ], + "source": [ + "files = sorted((project_root/'data/260323').glob('SUBT_*'))\n", + "acqs = [SiphraAcquisition(file) for file in files]\n", + "print(acqs)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "3b56d168-d4f0-4a4f-b917-97835015c884", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cs-137\n", + "Na22\n", + "Co-60\n", + "Am-241\n" + ] + } + ], + "source": [ + "# histograms\n", + "hists = {}\n", + "sources = []\n", + "for sgnl, bg in zip(acqs[1::2], acqs[2::2]):\n", + " filepath = sgnl.filepath.stem\n", + " src = (MetadataLoader.load(sgnl.metadataFile)).source\n", + " print(src)\n", + " sources.append(src)\n", + " # Signal and Background\n", + " hist_sgnl = ROOT.TH1F(f\"{src} Signal\", \"\", BITS12, 0, BITS12)\n", + " hist_bg = ROOT.TH1F(f\"{src} Background\", \"\", BITS12, 0 , BITS12)\n", + " hist_sgnl.Fill(sgnl['s']/len(sgnl.active_chs))\n", + " hist_bg.Fill(bg['s']/len(bg.active_chs))\n", + " hist_bg.Scale(sgnl.exposure/bg.exposure)\n", + " # Clean spectrum\n", + " hist_clean = hist_sgnl.Clone(f\"{src} Clean\")\n", + " hist_clean.Add(hist_bg, -1)\n", + " for hist in [hist_sgnl, hist_bg, hist_clean]:\n", + " # hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Normalized counts\")\n", + " hists[src] = hist_sgnl\n", + " hists[f\"{src}_BG\"] = hist_bg\n", + " hists[f\"{src}_clean\"] = hist_clean\n", + " del hist_sgnl, hist_bg" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "00444e1b-2030-452a-aeca-9ef6f89b66dd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1600, 1200)\n", + "cv.Divide(2,2)\n", + "\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lgds = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "for idx, (src, lg) in enumerate(zip(sources, lgds)): \n", + " cv.cd(idx+1)\n", + " \n", + " sgnl = hists[src]\n", + " bg = hists[src+'_BG']\n", + " clean = hists[src+'_clean']\n", + " \n", + " lg.AddEntry(sgnl, \"Signal\")\n", + " lg.AddEntry(bg, \"Background\")\n", + " lg.AddEntry(clean, \"Subtracted\")\n", + " sgnl.SetLineColor(colors[0])\n", + " bg.SetLineColor(colors[1])\n", + " clean.SetLineColor(colors[2])\n", + " sgnl.SetTitle(src)\n", + " sgnl.Draw(\"hist\")\n", + " bg.Draw(\"sames hist\")\n", + " clean.Draw(\"sames hist\")\n", + " ROOT.gPad.Update()\n", + " for i, sp in enumerate([sgnl, bg, clean]):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + " lg.Draw()\n", + " cv.cd(idx+1).SetLogy()\n", + " cv.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "9ecbda3c-0946-404a-afff-9e8014ef5a2c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.0035990386193816273, 0.0]\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 219.546\n", + "NDf = 147\n", + "Edm = 2.30957e-11\n", + "NCalls = 73\n", + "Constant = 2454.4 +/- 7.52715 \n", + "Mean = 183.938 +/- 0.0762262 \n", + "Sigma = 29.3406 +/- 0.0632264 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 0\n", + "NDf = 0\n", + "p0 = 0 +/- 1 \n", + "p1 = 0.00359904 +/- 0.00768853 \n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Fit data is empty \n" + ] + } + ], + "source": [ + "from analysis import *\n", + "\n", + "#Calibration based on the 3 peaks of the Na-22 spectrum\n", + "\n", + "\n", + "hist = hists['Cs-137_clean']\n", + "energy_ranges = [(0,1), (100, 250)]\n", + "energies = [0, Cs137_MeV]\n", + "\n", + "c_fit = calibration_fit(hist, energy_ranges, energies)\n", + "print(c_fit)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "22724987-e04a-4617-8c06-9fa511199544", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.3761890893184299]\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 518.581\n", + "NDf = 164\n", + "Edm = 1.52714e-08\n", + "NCalls = 73\n", + "Constant = 25.7054 +/- 0.0782673 \n", + "Mean = 0.66185 +/- 0.000273606 \n", + "Sigma = 0.105748 +/- 0.000220684 \t (limited)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n", + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#Calibrating Cs137 histograms based on Na22 calibration fit\n", + "hist_cal_bg_Cs137 = calibrated_histogram(c_fit, acqs[2], BITS12)\n", + "hist_cal_bg_Cs137.Scale(1/acqs[2].exposure) \n", + "hist_cal_src_Cs137 = calibrated_histogram(c_fit, acqs[1], BITS12)\n", + "hist_cal_src_Cs137.Scale(1/acqs[1].exposure) \n", + "hist_cal_clean_Cs137 = hist_cal_src_Cs137.Clone(\"Calibrated signal no BG\")\n", + "hist_cal_clean_Cs137.Add(hist_cal_bg_Cs137, -1)\n", + "\n", + "\n", + "# Plot the Cs-137 histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "hist_cal_bg_Cs137.SetLineColor(colors[0])\n", + "hist_cal_src_Cs137.SetLineColor(colors[1])\n", + "lg1.AddEntry(hist_cal_bg_Cs137, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_cal_src_Cs137, r\"Signal ^{137}Cs\", \"l\")\n", + "hist_cal_src_Cs137.Draw(\"hist\")\n", + "hist_cal_bg_Cs137.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "#hist_cal_bg.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_cal_clean_Cs137.SetLineColor(colors[2])\n", + "hist_cal_clean_Cs137.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()\n", + "\n", + "#Calculate energy resolution\n", + "res_Cs137 = energy_resolution(hist_cal_clean_Cs137, [(0.3, 0.9)], [Cs137_MeV])\n", + "print(res_Cs137)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "83d50dbc-584b-4327-8dca-120cc181bb0e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.6098052497160661, 0.3310662458402524, 0.19985352745576823]\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 114.793\n", + "NDf = 107\n", + "Edm = 1.05832e-07\n", + "NCalls = 78\n", + "Constant = 0.657598 +/- 0.0098944 \n", + "Mean = 0.487964 +/- 0.00391232 \n", + "Sigma = 0.132319 +/- 0.00277695 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 288.638\n", + "NDf = 91\n", + "Edm = 1.06692e-07\n", + "NCalls = 90\n", + "Constant = 2.06376 +/- 0.0165861 \n", + "Mean = 1.31918 +/- 0.00187749 \n", + "Sigma = 0.17924 +/- 0.0030605 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 748.173\n", + "NDf = 137\n", + "Edm = 5.20097e-06\n", + "NCalls = 67\n", + "Constant = 2.9031 +/- 0.0151085 \n", + "Mean = 1.73526 +/- 0.00089493 \n", + "Sigma = 0.151566 +/- 0.000800909 \t (limited)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n", + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#Calibration based on the 3 peaks of the Na-22 spectrum\n", + "hist = hists['Na22_clean']\n", + "energy_ranges = [(150, 200), (320, 420), (450, 520)]\n", + "energies = Na22_MeV\n", + "\n", + "c_fit = calibration_fit(hist, energy_ranges, energies)\n", + "print(c_fit)\n", + "\n", + "#Calibrating Na22 histograms based on Na22 calibration fit\n", + "hist_cal_bg_Na22 = calibrated_histogram(c_fit, acqs[4], BITS12)\n", + "hist_cal_bg_Na22.Scale(1/acqs[4].exposure) \n", + "hist_cal_src_Na22 = calibrated_histogram(c_fit, acqs[3], BITS12)\n", + "hist_cal_src_Na22.Scale(1/acqs[3].exposure) \n", + "hist_cal_clean_Na22 = hist_cal_src_Na22.Clone(\"Calibrated signal no BG\")\n", + "hist_cal_clean_Na22.Add(hist_cal_bg_Na22, -1)\n", + "\n", + "\n", + "# Plot the Cs-137 histograms\n", + "if ROOT.gROOT.FindObject('cv2'):\n", + " ROOT.gROOT.FindObject('cv2').Close()\n", + "cv2 = ROOT.TCanvas(\"cv2\", \"cv2\", 1600, 800)\n", + "cv2.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv2.cd(1)\n", + "hist_cal_bg_Na22.SetLineColor(colors[0])\n", + "hist_cal_src_Na22.SetLineColor(colors[1])\n", + "lg1.AddEntry(hist_cal_bg_Na22, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_cal_src_Na22, r\"Signal ^{22}Na\", \"l\")\n", + "hist_cal_src_Na22.Draw(\"hist\")\n", + "hist_cal_bg_Na22.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv2.cd(1).SetLogy()\n", + "#hist_cal_bg.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv2.cd(2)\n", + "hist_cal_clean_Na22.SetLineColor(colors[2])\n", + "hist_cal_clean_Na22.Draw(\"hist\")\n", + "cv2.cd(2).SetLogy()\n", + "\n", + "cv2.Draw()\n", + "\n", + "#Calculate energy resolution\n", + "res_Na22 = energy_resolution(hist_cal_clean_Na22, [(0.13, 0.6), (1.1, 1.5), (1.5 , 2.1)], Na22_MeV)\n", + "print(res_Na22)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "d2a4c5d9-9542-4795-92a3-16bc1d663518", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.0029093405113661537, -0.0009214807963730011]\n", + "[0.5258062473970452]\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 8.89298e-09\n", + "NDf = 0\n", + "Edm = 8.89072e-09\n", + "NCalls = 57\n", + "Constant = 3.63942 +/- 76.2026 \n", + "Mean = 0.316732 +/- 5.06048 \n", + "Sigma = 0.251842 +/- 1.78535 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 26.1098\n", + "NDf = 22\n", + "Edm = 1.46194e-05\n", + "NCalls = 71\n", + "Constant = 491.938 +/- 33.4685 \n", + "Mean = 20.7681 +/- 0.42253 \n", + "Sigma = 4.70087 +/- 0.220973 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 6.9514e-35\n", + "NDf = 0\n", + "p0 = -0.000921481 +/- 1.01561 \n", + "p1 = 0.00290934 +/- 0.0691501 \n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 32.6911\n", + "NDf = 24\n", + "Edm = 7.8732e-06\n", + "NCalls = 74\n", + "Constant = 5.93532 +/- 0.399578 \n", + "Mean = 0.0588274 +/- 0.00116134 \n", + "Sigma = 0.0132847 +/- 0.000589967 \t (limited)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n", + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#Calibration based on the 3 peaks of the Na-22 spectrum\n", + "hist = hists['Am-241_clean']\n", + "energy_ranges = [(0, 1), (5, 30)]\n", + "energies = [0, Am241_MeV]\n", + "\n", + "c_fit = calibration_fit(hist, energy_ranges, energies)\n", + "print(c_fit)\n", + "\n", + "#Calibrating Na22 histograms based on Na22 calibration fit\n", + "hist_cal_bg_Am241 = calibrated_histogram(c_fit, acqs[8], BITS12)\n", + "hist_cal_bg_Am241.Scale(1/acqs[8].exposure) \n", + "hist_cal_src_Am241 = calibrated_histogram(c_fit, acqs[7], BITS12)\n", + "hist_cal_src_Am241.Scale(1/acqs[7].exposure) \n", + "hist_cal_clean_Am241 = hist_cal_src_Am241.Clone(\"Calibrated signal no BG\")\n", + "hist_cal_clean_Am241.Add(hist_cal_bg_Am241, -1)\n", + "\n", + "\n", + "# Plot the Cs-137 histograms\n", + "if ROOT.gROOT.FindObject('cv4'):\n", + " ROOT.gROOT.FindObject('cv4').Close()\n", + "cv4 = ROOT.TCanvas(\"cv4\", \"cv4\", 1600, 800)\n", + "cv4.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv4.cd(1)\n", + "hist_cal_bg_Am241.SetLineColor(colors[0])\n", + "hist_cal_src_Am241.SetLineColor(colors[1])\n", + "lg1.AddEntry(hist_cal_bg_Am241, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_cal_src_Am241, r\"Signal ^{241}Am\", \"l\")\n", + "hist_cal_src_Am241.Draw(\"hist\")\n", + "hist_cal_bg_Am241.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv4.cd(1).SetLogy()\n", + "#hist_cal_bg.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv4.cd(2)\n", + "hist_cal_clean_Am241.SetLineColor(colors[2])\n", + "hist_cal_clean_Am241.Draw(\"hist\")\n", + "cv4.cd(2).SetLogy()\n", + "\n", + "cv4.Draw()\n", + "\n", + "#Calculate energy resolution\n", + "res_Am241 = energy_resolution(hist_cal_clean_Am241, [(0.01, 0.09)], [Am241_MeV])\n", + "print(res_Am241)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8dc817d8-9e6c-45c6-8045-f89abe669bb2", + "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.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/noteboks/QT_compare_channels.ipynb b/notebooks/QT_compare_channels.ipynb similarity index 100% rename from noteboks/QT_compare_channels.ipynb rename to notebooks/QT_compare_channels.ipynb diff --git a/noteboks/Untitled.ipynb b/notebooks/Untitled.ipynb similarity index 100% rename from noteboks/Untitled.ipynb rename to notebooks/Untitled.ipynb diff --git a/notebooks/Voltage_dependency.ipynb b/notebooks/Voltage_dependency.ipynb new file mode 100644 index 0000000..a5290db --- /dev/null +++ b/notebooks/Voltage_dependency.ipynb @@ -0,0 +1,695 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "248d5111-9526-4226-b6ef-a67dc0c71445", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "\n", + "from processing import SiphraAcquisition\n", + "from analysis import fit_peak_expbg\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c9421d1e-23fe-4585-905c-2f56603eafca", + "metadata": {}, + "outputs": [], + "source": [ + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "# Energy emission peaks in MeV\n", + "Cs237_MEV = 0.662" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a1ebd8ed-1b60-42dd-b3d9-2372f7596aad", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['Ch2', 'Ch4', 'Ch6', 'Ch8', 'Ch10', 'Ch12', 'Ch14', 'Ch16']\n" + ] + } + ], + "source": [ + "#Make sure json files have matching names to the dat files. \n", + "\n", + "channels = [f'Ch{n}' for n in range(2,17,2)]\n", + "print(channels)\n", + "files_lsts = {}\n", + "files_lsts[channels[0]] = sorted(Path(project_root+'/data/260311').glob('SUBTRACTED*Ch2*.csv'))\n", + "files_lsts[channels[0]].pop(3)\n", + "files_lsts[channels[0]].pop(3)\n", + "files_lsts[channels[0]].pop(3)\n", + "for ch in channels[1:]:\n", + " files_lsts[ch] = sorted(Path(project_root+'/data/260311').glob(f'SUBTRACTED*{ch}*.csv'))\n", + "\n", + "vbiases = [28.0, 28.5, 29.0]\n", + "# Datasets\n", + "datasets = {}\n", + "for ch, lst in files_lsts.items():\n", + " datasets[ch] = []\n", + " datasets[ch].append(SiphraAcquisition(lst[2]))\n", + " datasets[ch].append(SiphraAcquisition(lst[0]))\n", + " datasets[ch].append(SiphraAcquisition(lst[1]))\n", + "\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kCyan, ROOT.kYellow, ROOT.kSpring, ROOT.kViolet]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "4f82a8be-e799-4405-9464-bb085d20cdf3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "canvas_raw = []\n", + "cv_names = [\"canvas_raw_\"+str(_) for _ in range(1, 9)]\n", + "for name in cv_names:\n", + " if ROOT.gROOT.FindObject(name):\n", + " ROOT.gROOT.FindObject(name).Close()\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.15)\n", + "\n", + "legends_raw = [ROOT.TLegend(0.65, 0.7, 0.8, 0.88) for _ in range(1,9)]\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "hists_raw_dict = {}\n", + "hists_raw_names = [\"Ch2\", \"Ch4\", \"Ch6\", \"Ch8\", \"Ch10\", \"Ch12\", \"Ch14\", \"Ch16\"]\n", + "voltages = [\"28.0V\", \"28.5V\", \"29.0V\"]\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "for name, cv, lg in zip(hists_raw_names, cv_names, legends_raw):\n", + " dataset = datasets[name]\n", + " hists_raw_dict[name] = []\n", + " current_hist_list = hists_raw_dict[name]\n", + " canvas_raw.append(ROOT.TCanvas(cv, cv, 800, 600))\n", + " lg.SetHeader('Voltage')\n", + " # legends_raw.append(ROOT.TLegend(0.65, 0.7, 0.8, 0.88))\n", + " for idx, (acq, color) in enumerate(zip(dataset, colors)):\n", + " # print(f\"Current file: {acq.filepath.name}\")\n", + " ch = acq.active_chs[0]\n", + " hist = ROOT.TH1F(f\"{name}, {voltages[idx]}\", \"\", BITS9, 0, BITS12)\n", + " hist.Fill(acq[ch])\n", + " hist.Scale(1/acq.exposure)\n", + " #Preeting up..\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Count rate (s^{-1})\")\n", + " hist.SetLineColor(color+2)\n", + " hist.SetTitle(name)\n", + " lg.AddEntry(hist, f\"{voltages[idx]}\", 'l')\n", + " current_hist_list.append(hist)\n", + " current_hist_list[-1].Draw('sames hist')\n", + " lg.Draw()\n", + " canvas_raw[-1].SetLogy()\n", + " ROOT.gPad.Update()\n", + " for i, sp in enumerate(current_hist_list):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + " #legends_raw[-1].Draw()\n", + " canvas_raw[-1].Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1931d8ae-6b41-4b2d-85cf-33da23d47698", + "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.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/noteboks/banana.ipynb b/notebooks/banana.ipynb similarity index 100% rename from noteboks/banana.ipynb rename to notebooks/banana.ipynb diff --git a/notebooks/baseline_subtraction.ipynb b/notebooks/baseline_subtraction.ipynb new file mode 100644 index 0000000..ced9272 --- /dev/null +++ b/notebooks/baseline_subtraction.ipynb @@ -0,0 +1,1143 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "cda6a9c1-c4c0-4e69-b50d-7dcb1b8434f2", + "metadata": {}, + "source": [ + "# Performing baselie subtraction\n", + "\n", + "The baseline-subtracted histograms of an acquisition with 2 active channels were obtained. The baseline value for each channel is the mean of the readings from all the events triggered ONLY by external HOLD assertion. This value is subtracted from the values of each \"real\" event (triggered from internal HOLD).\n", + "\n", + "The dataset for this analysis is a background acquisition of 10 000 events. SIPHRA channels 2 and 14 were connected to SiPM channels 1-4 and 5-8, respectively, using charge comparator, with a threshold value of 20." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "initial_id", + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-10T09:37:32.276396Z", + "start_time": "2026-02-10T09:37:32.169808Z" + } + }, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "\n", + "from processing import SiphraAcquisition\n", + "from analysis import fit_peak_expbg\n", + "\n", + "# ROOT.gROOT.SetBatch(False)\n", + "# ROOT.gEnv.SetValue(\"Canvas.ShowEventStatus\", 1)\n", + "# ROOT.gEnv.SetValue(\"Canvas.ShowToolBar\", 1)\n", + "# ROOT.gEnv.SetValue(\"Canvas.UseGL\", 0)\n", + "# ROOT.gROOT.ProcessLine(\"gVirtualX = new TGX11();\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7aa83758-b8f5-4a20-96c3-c8381ea71536", + "metadata": {}, + "outputs": [], + "source": [ + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "# Energy emission peaks in MeV\n", + "K40_MEV = 1.460\n", + "TL208_MEV = 2.614\n", + "CS137_MEV = 0.661\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "68a9cd08dc972f65", + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-10T09:37:35.457562Z", + "start_time": "2026-02-10T09:37:33.669215Z" + } + }, + "outputs": [], + "source": [ + "# Datasets\n", + "acq_notSubtracted = SiphraAcquisition(project_root+'/data/260203/1_SiPM_ChannelsTest_Ch1-4_Ch2_Ch5-8_Ch14_QT_Thr20_Hys0_Background.csv',\n", + " active_chs=[2,14],\n", + " exposure_sec=1)\n", + "acq_subtracted = SiphraAcquisition(project_root+'/data/260203/1_datdecodertest_subtractBaseline.csv',\n", + " active_chs=[2,14],\n", + " exposure_sec=1)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "cae21a3ee684d69f", + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-10T09:43:25.401070Z", + "start_time": "2026-02-10T09:43:24.995723Z" + }, + "jupyter": { + "source_hidden": true + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "SUMMED_XR=5000 # Right limit of the summed channel's x-axis\n", + "\n", + "c = ROOT.TCanvas('cv', 'cv', 1200, 800)\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "names=['BL subt.', 'No correction']\n", + "hists = []\n", + "\n", + "lg = ROOT.TLegend(0.48, 0.61, 0.75, 0.83)\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "for idx, (acq, name) in enumerate(zip([acq_subtracted, acq_notSubtracted], names)):\n", + " # print(f\"Current file: {acq.filepath.name}\")\n", + " ch = 'Summed' #acq.active_chs[0]\n", + " hist = ROOT.TH1F(f\"{ch} {name}\", \"\", int(BITS9*SUMMED_XR/BITS12), 0, SUMMED_XR)\n", + " hist_singlech = ROOT.TH1F(f\"Ch. {acq.active_chs[0]} {name}\", \"\", BITS9, 0, BITS12)\n", + " hist.Fill(acq[ch]/len(acq.active_chs))\n", + " hist_singlech.Fill(acq[acq.active_chs[0]])\n", + " hist.Scale(1/acq.exposure)\n", + " hist_singlech.Scale(1/acq.exposure)\n", + "\n", + " hist_singlech2 = ROOT.TH1F(f\"Ch. {acq.active_chs[1]} {name}\", \"\", BITS9, 0, BITS12)\n", + " hist_singlech2.Fill(acq[acq.active_chs[1]])\n", + " hist_singlech2.Scale(1/acq.exposure)\n", + " #Preeting up..\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Counts\")\n", + " hist.SetLineColor(colors[idx])\n", + " hist.SetTitle(\"Baseline subtraction comparison\")\n", + " # hist.GetYaxis().SetRangeUser(0, 1e4)\n", + " lg.SetHeader(\"SIPHRA Channel\")\n", + " lg.AddEntry(hist, f\"{ch} {name}\", 'l')\n", + " hists.append(hist)\n", + " hists[-1].Draw('sames hist')\n", + " # hist_singlech.GetXaxis().SetTitle(\"ADC channel number\")\n", + " # hist_singlech.GetYaxis().SetTitle(r\"Count rate (s^{-1})\")\n", + " hist_singlech.SetLineColor(colors[idx + 2]+2)\n", + " hist_singlech2.SetLineColor(colors[idx + 4]+1)\n", + " # hist_singlech.SetTitle(\"Baseline subtraction comparison\")\n", + " lg.AddEntry(hist_singlech, f\"Ch. {acq.active_chs[0]} {name}\", 'l')\n", + " lg.AddEntry(hist_singlech2, f\"Ch. {acq.active_chs[1]} {name}\", 'l')\n", + " hists.append(hist_singlech)\n", + " hists[-1].Draw('sames hist')\n", + " hists.append(hist_singlech2)\n", + " hists[-1].Draw('sames hist')\n", + "c.SetLogy()\n", + "ROOT.gPad.Update()\n", + "for i, sp in enumerate(hists):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + "lg.Draw()\n", + "c.Draw()\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "4fe40546-908d-46cd-868c-f4c3e9bd0902", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv_4plots'):\n", + " ROOT.gROOT.FindObject('cv_4plots').Close()\n", + "\n", + "canvas4 = ROOT.TCanvas('cv_4plots', 'cv_4plots', 1600, 1200)\n", + "canvas4.Divide(2,2)\n", + "\n", + "legends = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "canvas4.cd(1)\n", + "hists[0].Draw('hist')\n", + "hists[3].Draw('sames hist')\n", + "legends[0].AddEntry(hists[0], \"Baseline subtracted\", 'l')\n", + "legends[0].AddEntry(hists[3], \"Not corrected\", 'l')\n", + "legends[0].SetHeader(\"\\'Summed\\' channel spectra\")\n", + "legends[0].Draw()\n", + "canvas4.cd(1).SetLogy()\n", + "\n", + "canvas4.cd(2)\n", + "hists[1].SetXTitle(\"ADC channel number\")\n", + "hists[1].SetYTitle(\"Counts\")\n", + "hists[1].Draw('hist')\n", + "hists[4].Draw('sames hist')\n", + "legends[1].AddEntry(hists[1], \"Baseline subtracted\", 'l')\n", + "legends[1].AddEntry(hists[4], \"Not corrected\", 'l')\n", + "legends[1].SetHeader(\"Ch. 2 spectra\")\n", + "legends[1].Draw()\n", + "canvas4.cd(2).SetLogy()\n", + "\n", + "canvas4.cd(3)\n", + "hists[2].SetXTitle(\"ADC channel number\")\n", + "hists[2].SetYTitle(\"Counts\")\n", + "hists[2].Draw('hist')\n", + "hists[5].Draw('sames hist')\n", + "legends[2].AddEntry(hists[2], \"Baseline subtracted\", 'l')\n", + "legends[2].AddEntry(hists[5], \"Not corrected\", 'l')\n", + "legends[2].SetHeader(\"Ch. 14 spectra\")\n", + "legends[2].Draw()\n", + "canvas4.cd(3).SetLogy()\n", + "\n", + "canvas4.cd(4)\n", + "h0_zoom = hists[0].Clone(\"hists_0_zoomed\")\n", + "h0_zoom.GetXaxis().SetRangeUser(0,BITS12)\n", + "h0_zoom.Draw('hist')\n", + "hists[1].Draw('sames hist')\n", + "hists[2].Draw('sames hist')\n", + "legends[3].AddEntry(h0_zoom, \"\\'Summed\\'\", 'l')\n", + "legends[3].AddEntry(hists[1], \"Ch. 2\", 'l')\n", + "legends[3].AddEntry(hists[2], \"Ch. 14\", 'l')\n", + "legends[3].SetHeader(\"Corrected channels\")\n", + "legends[3].Draw()\n", + "canvas4.cd(4).SetLogy()\n", + "\n", + "canvas4.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "46bd074b-d5aa-48a4-b027-f966314995cd", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ADC channels range:\n", + "--------------------------------------------------\n", + "Channel\t\tBL-subtracted\tNo correction\n", + "--------------------------------------------------\n", + "Summed\t\t0.0 - 3965.0\t935.5 - 4880.5\n", + "Ch. 2\t\t0.0 - 3951.0\t51.0 - 4095.0\n", + "Ch. 14\t\t0.0 - 3979.0\t0.0 - 4095.0\n", + "\n", + "\n", + "\n", + "Number of values above single channel range in 'Summed' channel:\n", + "------------------------------\n", + "BL-subtracted\tNo correction\n", + "------------------------------\n", + " 0 \t 22 \n" + ] + } + ], + "source": [ + "print(\"ADC channels range:\\n\"\n", + " f\"{\"\":->50}\\n\"\n", + " \"Channel\\t\\tBL-subtracted\\tNo correction\\n\"\n", + " f\"{\"\":->50}\\n\"\n", + " f\"Summed\\t\\t{acq_subtracted['s'].min()/2} - {acq_subtracted['s'].max()/2}\\t{acq_notSubtracted['s'].min()/2} - {acq_notSubtracted['s'].max()/2}\\n\"\n", + " f\"Ch. 2\\t\\t{acq_subtracted[2].min()} - {acq_subtracted[2].max()}\\t{acq_notSubtracted[2].min()} - {acq_notSubtracted[2].max()}\\n\"\n", + " f\"Ch. 14\\t\\t{acq_subtracted[14].min()} - {acq_subtracted[14].max()}\\t{acq_notSubtracted[14].min()} - {acq_notSubtracted[14].max()}\\n\")\n", + "\n", + "print(\"\\n\\n\"\n", + " \"Number of values above single channel range in \\'Summed\\' channel:\\n\"\n", + " f\"{\"\":->30}\\n\"\n", + " \"BL-subtracted\\tNo correction\\n\"\n", + " f\"{\"\":->30}\\n\"\n", + " f\"{len(acq_subtracted['s'][acq_subtracted['s'][:]/2>BITS12]):^13}\\t{len(acq_notSubtracted['s'][acq_notSubtracted['s'][:]/2>BITS12]):^13}\")" + ] + }, + { + "cell_type": "markdown", + "id": "ae86c5f0-762f-4cca-a8bb-b0c00694959e", + "metadata": {}, + "source": [ + "The baseline subtraction is performed directly at the conversion script. No averaging over the number of channels is carried out by the conversion script, but the above spectra and data correspond to values averaged over the number of active channels (2).\n", + "\n", + "The values of the baselines of all the channels are between 98 and 150.\n", + "\n", + "The range of the summed channel does not exceed the maximum single-channel value (4096) after baseline correction.\n", + "\n", + "**Questions:**\n", + "- Is this baseline correction enough for pedestal subtraction?\n", + "- What about the large peak at the lower range of the spectrum? Is it still originating from noise at this stage?" + ] + }, + { + "cell_type": "markdown", + "id": "19798694-0ef8-43fb-a09c-6b87cc8b03f5", + "metadata": {}, + "source": [ + "# Calibration of corrected spectrum with the backgrond $^{40}$K and $^{208}$Tl peaks\n", + "\n", + "The current process of calibration to actual energies is shown here. The two peaks of $^{40}$K and $^{208}$Tl present in the background are used for this purpose. The emission lines of these nuclides are located at 1.460 MeV and 2.614 MeV, respectively. First, both peaks are fitted to a Gaussian with exponentially-decaying background in the baseline-subtracted background histogram. The Gaussian means are taken as the actual locations of the emission lines. The raw data is linearly transformed to aligh with these two points.\n", + "\n", + "To verify the calibration, we use the same calibration parameters to fit a spectrum taken in the pressence of a sample of $^{137}$Cs. After background subtraction, it is expected that the mean of the Gaussian fit of the resulting peak, lands at 0.661MeV." + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "id": "a4f39db9-ad0d-4f5b-b3b5-ea67040593a6", + "metadata": {}, + "outputs": [], + "source": [ + "# Datasets\n", + "folder = Path(project_root)/'data/260209'\n", + "times_lst = [0, 119.634, 161.062, 148.920, 164.828, 34.836, 80.792, 54.386, 90.837, 24.645, 24.144, 24.606, 23.381, 24.720, 24.916, 24.848, 24.934] # Index coincide with starting number of file name\n", + "\n", + "# Select SiPM channels 5-8, SIPHRA channel 14 -> gives cleaner signal\n", + "# Background is file with #4, source is file with #16, with preffix 'BLsubtr_'\n", + "BGIDX = 4; SRCIDX=16; ACTIVECH=14\n", + "bgfile = sorted(folder.glob(f'BLsubtr_{BGIDX}*.csv'))[0]\n", + "srcfile = sorted(folder.glob(f'BLsubtr_{SRCIDX}*.csv'))[0]\n", + "\n", + "acq_bg = SiphraAcquisition(bgfile, active_chs=ACTIVECH, exposure_sec=times_lst[BGIDX])\n", + "acq_src = SiphraAcquisition(srcfile, active_chs=ACTIVECH, exposure_sec=times_lst[SRCIDX])\n", + "colors_clb = [colors[1]-1, colors[3]-1, colors[0]-2] # colors for BG, signal and clean, respectively\n", + "NORMFACTOR = acq_src.exposure/(acq_bg.exposure)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "0b44ed9c-6f45-4174-8f81-1407bcdacc05", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ADC channels range after baseline subtraction:\n", + "CHANNEL 14\n", + "------------------------------\n", + "Data set\tRange\n", + "------------------------------\n", + "Background\t0.0 - 3924.0\n", + "Source\t\t0.0 - 3797.0\n", + "\n", + "SUMMED CHANNEL\n", + "------------------------------\n", + "Data set\tRange\n", + "------------------------------\n", + "Background\t0.0 - 3924.0\n", + "Source\t\t0.0 - 3797.0\n", + "\n" + ] + } + ], + "source": [ + "print(\"ADC channels range after baseline subtraction:\\n\"\n", + " \"CHANNEL 14\\n\"\n", + " f\"{\"\":->30}\\n\"\n", + " \"Data set\\tRange\\n\"\n", + " f\"{\"\":->30}\\n\"\n", + " f\"Background\\t{acq_bg[ACTIVECH].min()} - {acq_bg[ACTIVECH].max()}\\n\"\n", + " f\"Source\\t\\t{acq_src[ACTIVECH].min()} - {acq_src[ACTIVECH].max()}\\n\\n\"\n", + " \"SUMMED CHANNEL\\n\"\n", + " f\"{\"\":->30}\\n\"\n", + " \"Data set\\tRange\\n\"\n", + " f\"{\"\":->30}\\n\"\n", + " f\"Background\\t{acq_bg['s'].min()} - {acq_bg['s'].max()}\\n\"\n", + " f\"Source\\t\\t{acq_src['s'].min()} - {acq_src['s'].max()}\\n\")" + ] + }, + { + "cell_type": "markdown", + "id": "cddfadfd-d1a7-4257-9712-15c6daf6d565", + "metadata": {}, + "source": [ + "Baseline subtraction is effectively filtering out noise from inactive channels. It is noteworthy that, when subtracting single channels, the upper ADC channels will be empty, hence there is a loss of dynamic range due to baseline subtraction." + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "2093b302-1f12-4e19-86f6-74506d4e5ff6", + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: Background (Potential memory leak).\n", + "Warning in : Replacing existing TH1: Signal 137Cs (Potential memory leak).\n" + ] + } + ], + "source": [ + "# Raw histograms\n", + "hist_bg = ROOT.TH1F(\"Background\", \"\", BITS9, 0, BITS12)\n", + "hist_bg.Fill(acq_bg[ACTIVECH])\n", + "hist_bg.Scale(NORMFACTOR) # Normalize counts to the source exposure time\n", + "\n", + "hist_src = ROOT.TH1F(\"Signal 137Cs\", \"\", BITS9, 0, BITS12)\n", + "hist_src.Fill(acq_src[ACTIVECH])\n", + "\n", + "# Bacground subtraction\n", + "hist_clean = hist_src.Clone(\"137Cs no BG\")\n", + "hist_clean.Add(hist_bg, -1)\n", + "# Removing negative bins\n", + "for i in range(1, hist_clean.GetNbinsX() + 1): # bin 0 is the underflow\n", + " if hist_clean.GetBinContent(i) < 0: hist_clean.SetBinContent(i, 0)" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "id": "eacf6495-684b-42fa-9eab-8ca73fb2afef", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Error in : Cannot set Y axis to log scale\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot all the histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "hist_bg.GetXaxis().SetTitle(\"ADC Channel\")\n", + "hist_bg.GetYaxis().SetTitle(\"Normalized counts\")\n", + "hist_bg.SetLineColor(colors_clb[0])\n", + "hist_bg.SetFillColor(colors_clb[0])\n", + "hist_src.SetLineColor(colors_clb[1])\n", + "hist_src.SetFillColorAlpha(colors_clb[1], 0.8)\n", + "lg1.AddEntry(hist_bg, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_src, r\"Signal ^{137}Cs\", \"l\")\n", + "hist_bg.Draw(\"hist\")\n", + "hist_src.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "hist_bg.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "# lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_clean.GetXaxis().SetTitle(\"ADC Channel\")\n", + "hist_clean.GetYaxis().SetTitle(\"Normalized counts\")\n", + "hist_clean.SetLineColor(colors_clb[2])\n", + "hist_clean.SetFillColor(colors_clb[2])\n", + "hist_clean.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "id": "1b94d6b2-ea7d-41ea-acb2-a7e820e9dd34", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "MinFCN = 3.07992\n", + "Chi2 = 6.15985\n", + "NDf = 45\n", + "Edm = 4.33544e-06\n", + "NCalls = 378\n", + "Const = 967.676 +/- 315.198 \n", + "Denom = 204.714 +/- 16.8531 \n", + "Norm = 26.0236 +/- 2.92104 \n", + "Mean = 748.998 +/- 10.0065 \n", + "Sigma = 69.1169 +/- 10.4622 \n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "MinFCN = 4.22284\n", + "Chi2 = 8.44567\n", + "NDf = 45\n", + "Edm = 1.92878e-06\n", + "NCalls = 607\n", + "Const = 65.7641 +/- 57.5374 \n", + "Denom = 441.486 +/- 153.51 \n", + "Norm = 1.79752 +/- 0.994549 \n", + "Mean = 1240.31 +/- 25.9607 \n", + "Sigma = 44.4995 +/- 28.325 \n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Deleting canvas with same name: cv_fit\n", + "Error in : Cannot set Y axis to log scale\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Fit background peaks\n", + "xrange_k40 = (600,1000)\n", + "nrm_k40_i = 55\n", + "mean_k40_i = 740\n", + "xrange_tl208 = (1000,1400)\n", + "nrm_tl208_i = 7\n", + "mean_tl208_i = 1210\n", + "hist_bg_fit = hist_bg.Clone()\n", + "hist_bg_fit.GetYaxis().SetRangeUser(0,1e3)\n", + "hist_bg.SetTitle(r\"^{40}K and ^{208}Tl peaks fitting\")\n", + "fit_fn_k40, fit_res_k40 = fit_peak_expbg(hist_bg_fit, name=\"K40_peak\", xl=xrange_k40[0], xr=xrange_k40[1], norm=nrm_k40_i, mean=mean_k40_i, showFit=True)\n", + "fit_fn_tl208, fit_res_tl208 = fit_peak_expbg(hist_bg_fit, name=\"Tl208_peak\", xl=xrange_tl208[0], xr=xrange_tl208[1], norm=nrm_tl208_i, mean=mean_tl208_i, showFit=True)\n", + "\n", + "cv_fit = ROOT.TCanvas(\"cv_fit\", \"cv_fit\", 800,600)\n", + "hist_bg_fit.SetFillStyle(0)\n", + "hist_bg_fit.Draw(\"hist\")\n", + "cv_fit.SetLogy()\n", + "cv_fit.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 103, + "id": "3e0906bf-c280-4c5e-b555-871735d01a54", + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true, + "source_hidden": true + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: Calibrated background (Potential memory leak).\n", + "Warning in : Replacing existing TH1: Calibrated signal (Potential memory leak).\n" + ] + } + ], + "source": [ + "# Calibrate data and get calibrated histograms\n", + "K40FIT = fit_fn_k40.GetParameter('Mean')\n", + "TL208FIT = fit_fn_tl208.GetParameter('Mean')\n", + "clb_slope = (TL208_MEV - K40_MEV)/(TL208FIT - K40FIT)\n", + "clb_const = K40_MEV - clb_slope * K40FIT\n", + "\n", + "data_bg_clb = acq_bg[ACTIVECH]*clb_slope + clb_const\n", + "data_src_clb = acq_src[ACTIVECH]*clb_slope + clb_const\n", + "\n", + "CLBHISTS_XR = np.max([data_bg_clb.max(), data_src_clb.max()])\n", + "\n", + "hist_bg_clb = ROOT.TH1F(\"Calibrated background\", \"Calibrated Histograms\", BITS9, 0, CLBHISTS_XR)\n", + "hist_src_clb = ROOT.TH1F(\"Calibrated signal\", \"Calibrated Histograms\", BITS9, 0, CLBHISTS_XR)\n", + "\n", + "hist_bg_clb.Fill(data_bg_clb)\n", + "hist_bg_clb.Scale(NORMFACTOR)\n", + "\n", + "hist_src_clb.Fill(data_src_clb)\n", + "\n", + "hist_clean_clb = hist_src_clb.Clone(\"Calibrated signal no BG\")\n", + "hist_clean_clb.Add(hist_bg_clb, -1)\n", + "\n", + "for hist, clr in zip([hist_bg_clb, hist_src_clb, hist_clean_clb], colors_clb):\n", + " hist.GetXaxis().SetTitle(\"Energy (MeV)\")\n", + " hist.GetYaxis().SetTitle(\"Normalized counts\")\n", + " hist.SetLineColor(clr)" + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "id": "5416087a-578e-4d42-8157-2516db3c6b39", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Error in : Cannot set Y axis to log scale\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot all the histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "hist_bg_clb.SetFillColor(colors_clb[0])\n", + "hist_src_clb.SetFillColorAlpha(colors_clb[1], 0.8)\n", + "lg1.AddEntry(hist_bg_clb, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_src_clb, r\"Signal ^{137}Cs\", \"l\")\n", + "hist_bg_clb.Draw(\"hist\")\n", + "hist_src_clb.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "hist_bg_clb.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "# lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_clean_clb.SetFillColor(colors_clb[2])\n", + "hist_clean_clb.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 115, + "id": "dcbad448-c208-4054-a4ab-e2b71b7ed2a4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "MinFCN = 5628.78\n", + "Chi2 = 11257.6\n", + "NDf = 112\n", + "Edm = 1.88824e-07\n", + "NCalls = 124\n", + "Constant = 4002.29 +/- 17.4275 \n", + "Mean = 0.479418 +/- 0.000489701 \n", + "Sigma = 0.136089 +/- 0.000348374 \t (limited)\n" + ] + } + ], + "source": [ + "res = hist_clean_clb.Fit(\"gaus\", \"L S\", \"\", 0, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 129, + "id": "6d7f6090-75e8-4779-8a1c-35f1cad0a43d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Error in the location of Cs peak: 27.5%\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Deleting canvas with same name: cv_fit\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "cv_fit = ROOT.TCanvas(\"cv_fit\", \"cv_fit\", 800,600)\n", + "# hist_bg_fit.SetFillStyle(0)\n", + "hist_clean_clb.Draw(\"hist\")\n", + "hist_clean_clb.SetFillColor(colors_clb[2]-6)\n", + "cv_fit.SetLogy()\n", + "cv_fit.Draw()\n", + "\n", + "\n", + "err = (CS137_MEV - hist_clean_clb.GetFunction(\"gaus\").GetParameter(\"Mean\"))/CS137_MEV\n", + "print(f\"Error in the location of Cs peak: {err*100:.1f}%\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "58f9c3b4-03c0-453a-9d37-cd9cd565ad5f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "FWHM=0.32 Energy resolution at 480MeV: 66.84%\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "sigma = 0.136089\n", + "FWHM = np.sqrt(8*np.log(2))*sigma\n", + "print(f\"{FWHM=:.2f} Energy resolution at 480MeV: {FWHM*100/0.479418:.2f}%\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (ROOT in Conda)", + "language": "python", + "name": "comcube" + }, + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/energyResolution-Copy1.ipynb b/notebooks/energyResolution-Copy1.ipynb new file mode 100644 index 0000000..9eddcb7 --- /dev/null +++ b/notebooks/energyResolution-Copy1.ipynb @@ -0,0 +1,539 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4a3b1bc5-ec73-4d69-a840-07d057f4e2eb", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "project_root=Path(project_root)\n", + "\n", + "from processing import SiphraAcquisition, MetadataLoader\n", + "from analysis import fit_peak_expbg\n", + "\n", + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "\n", + "# Energy emission peaks in MeV\n", + "Cs137_MeV = 0.662\n", + "Na22_MeV = [0.511, 1.275, 1.786]\n", + "Co60_MeV = [1.173, 1.332, 2,505]\n", + "Am241_MeV = 0.0595\n", + "# Background emission\n", + "K40_MeV = 1.460\n", + "Tl208_MEV = 2.614\n", + "\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "75493f51-2089-4f0a-8057-0f2c2bf7466f", + "metadata": {}, + "outputs": [], + "source": [ + "files = sorted((project_root/'data/260323').glob('SUBT_*'))\n", + "acqs = [SiphraAcquisition(file) for file in files]" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "88f4c856-78c0-4ca6-ad58-cc10f3b53649", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cs-137\n", + "Na22\n", + "Co-60\n", + "Am-241\n" + ] + } + ], + "source": [ + "# histograms\n", + "hists = {}\n", + "sources = []\n", + "for sgnl, bg in zip(acqs[1::2], acqs[2::2]):\n", + " filepath = sgnl.filepath.stem\n", + " src = (MetadataLoader.load(sgnl.metadataFile)).source\n", + " sources.append(src)\n", + " print(src)\n", + " # Signal and Background\n", + " hist_sgnl = ROOT.TH1F(f\"{src} Signal\", \"\", BITS12, 0, BITS12)\n", + " hist_bg = ROOT.TH1F(f\"{src} Background\", \"\", BITS12, 0 , BITS12)\n", + " hist_sgnl.Fill(sgnl['s']/len(sgnl.active_chs))\n", + " hist_bg.Fill(bg['s']/len(bg.active_chs))\n", + " hist_bg.Scale(sgnl.exposure/bg.exposure)\n", + " # Clean spectrum\n", + " hist_clean = hist_sgnl.Clone(f\"{src} Clean\")\n", + " hist_clean.Add(hist_bg, -1)\n", + " for hist in [hist_sgnl, hist_bg, hist_clean]:\n", + " # hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Normalized counts\")\n", + " hists[src] = hist_sgnl\n", + " hists[f\"{src}_BG\"] = hist_bg\n", + " hists[f\"{src}_clean\"] = hist_clean\n", + " del hist_sgnl, hist_bg\n", + "#print(hists)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "36f654a6-bd64-471f-a0e3-01927db52e3e", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1600, 1200)\n", + "cv.Divide(2,2)\n", + "\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lgds = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "for idx, (src, lg) in enumerate(zip(sources, lgds)): \n", + " cv.cd(idx+1)\n", + " \n", + " sgnl = hists[src]\n", + " bg = hists[src+'_BG']\n", + " clean = hists[src+'_clean']\n", + " \n", + " lg.AddEntry(sgnl, \"Signal\")\n", + " lg.AddEntry(bg, \"Background\")\n", + " lg.AddEntry(clean, \"Subtracted\")\n", + " sgnl.SetLineColor(colors[0])\n", + " bg.SetLineColor(colors[1])\n", + " clean.SetLineColor(colors[2])\n", + " sgnl.SetTitle(src)\n", + " sgnl.Draw(\"hist\")\n", + " bg.Draw(\"sames hist\")\n", + " clean.Draw(\"sames hist\")\n", + " ROOT.gPad.Update()\n", + " for i, sp in enumerate([sgnl, bg, clean]):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + " lg.Draw()\n", + " cv.cd(idx+1).SetLogy()\n", + " cv.Draw()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "6db7ba78-ddd9-4d76-8bba-bf310f8f1221", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.004280158214429103\n", + "-0.2655932438428593\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 196.518\n", + "NDf = 127\n", + "Edm = 2.29967e-07\n", + "NCalls = 65\n", + "Constant = 2459.66 +/- 7.74946 \n", + "Mean = 183.993 +/- 0.078652 \n", + "Sigma = 29.2438 +/- 0.0713573 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 356.977\n", + "NDf = 147\n", + "Edm = 1.20533e-06\n", + "NCalls = 64\n", + "Constant = 1361.97 +/- 6.38965 \n", + "Mean = 275.54 +/- 0.144272 \n", + "Sigma = 32.7382 +/- 0.141925 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 52.0578\n", + "NDf = 47\n", + "Edm = 1.53695e-06\n", + "NCalls = 107\n", + "Constant = 255.901 +/- 3.65957 \n", + "Mean = 478.856 +/- 0.861944 \n", + "Sigma = 32.2107 +/- 1.8725 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 50.431\n", + "NDf = 47\n", + "Edm = 6.26293e-06\n", + "NCalls = 110\n", + "Constant = 187.956 +/- 4.04943 \n", + "Mean = 177.384 +/- 1.67095 \n", + "Sigma = 37.5292 +/- 4.35406 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 152.171\n", + "NDf = 97\n", + "Edm = 2.80563e-08\n", + "NCalls = 96\n", + "Constant = 606.675 +/- 4.80363 \n", + "Mean = 371.93 +/- 0.452134 \n", + "Sigma = 45.3973 +/- 0.782852 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 63.6754\n", + "NDf = 67\n", + "Edm = 1.95162e-06\n", + "NCalls = 89\n", + "Constant = 880.588 +/- 5.38234 \n", + "Mean = 471.391 +/- 0.435139 \n", + "Sigma = 32.0034 +/- 0.470993 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 0.00408946\n", + "NDf = 1\n", + "p0 = -0.265593 +/- 0.109303 \n", + "p1 = 0.00428016 +/- 0.000302376 \n" + ] + } + ], + "source": [ + "#Fitting all the visible peaks. Values around peak input manually\n", + "\n", + "Cs137_fit = ROOT.TF1(\"Cs137_fit\", \"gaus\", 120,250)\n", + "\n", + "Co60a_fit = ROOT.TF1(\"Co60a_fit\", \"gaus\", 200,350)\n", + "Co60b_fit = ROOT.TF1(\"Co60b_fit\", \"gaus\", 450,500)\n", + "\n", + "Na22a_fit = ROOT.TF1(\"Na22a_fit\", \"gaus\", 150,200)\n", + "Na22b_fit = ROOT.TF1(\"Na22b_fit\", \"gaus\", 320,420)\n", + "Na22c_fit = ROOT.TF1(\"Na22c_fit\", \"gaus\", 450,520)\n", + "\n", + "fits = [Cs137_fit, Co60a_fit, Co60b_fit, Na22a_fit, Na22b_fit, Na22c_fit]\n", + "\n", + "#Fit to histogram and find mean value of peak\n", + "hists['Cs-137_clean'].Fit(Cs137_fit, \"R+S\")\n", + "hists['Co-60_clean'].Fit(Co60a_fit, \"R+S\")\n", + "hists['Co-60_clean'].Fit(Co60b_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22a_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22b_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22c_fit, \"R+S\")\n", + "\n", + "Cs137_mean = Cs137_fit.GetParameter('Mean')\n", + "Co60a_mean = Co60a_fit.GetParameter('Mean')\n", + "Co60b_mean = Co60b_fit.GetParameter('Mean')\n", + "Na22a_mean = Na22a_fit.GetParameter('Mean')\n", + "Na22b_mean = Na22b_fit.GetParameter('Mean')\n", + "Na22c_mean = Na22c_fit.GetParameter('Mean')\n", + "\n", + "channels = [Na22a_mean, Na22b_mean, Na22c_mean]\n", + "energies = Na22_MeV\n", + "channels = np.array(channels, dtype='float64')\n", + "energies = np.array(energies, dtype='float64')\n", + "\n", + "g = ROOT.TGraph(len(channels), channels, energies)\n", + "\n", + "f = ROOT.TF1(\"calib\", \"pol1\", min(channels), max(channels))\n", + "g.Fit(f)\n", + "\n", + "a = f.GetParameter(1)\n", + "b = f.GetParameter(0)\n", + "print(a)\n", + "print(b)\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()\n", + "\n", + "#lg.Draw()\n", + "#cv.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "6b5a5060-fe8f-48ba-af47-c12bca5ea672", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 137.647\n", + "NDf = 67\n", + "Edm = 2.07898e-06\n", + "NCalls = 84\n", + "Constant = 189.854 +/- 2.40817 \n", + "Mean = 0.492681 +/- 0.00341376 \n", + "Sigma = 0.142338 +/- 0.00437932 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 155.228\n", + "NDf = 91\n", + "Edm = 3.6702e-08\n", + "NCalls = 96\n", + "Constant = 613.942 +/- 3.98488 \n", + "Mean = 1.32007 +/- 0.00148807 \n", + "Sigma = 0.180551 +/- 0.00242551 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 372.407\n", + "NDf = 90\n", + "Edm = 6.00943e-09\n", + "NCalls = 92\n", + "Constant = 833.135 +/- 4.70932 \n", + "Mean = 1.74383 +/- 0.00131477 \n", + "Sigma = 0.172367 +/- 0.0019528 \t (limited)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: h_cal (Potential memory leak).\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv_cal'):\n", + " ROOT.gROOT.FindObject('cv_cal').Close()\n", + "\n", + " \n", + "hist = hists['Na22_clean']\n", + "\n", + "emin = a * hist.GetXaxis().GetXmin() + b\n", + "emax = a * hist.GetXaxis().GetXmax() + b\n", + "\n", + "h_cal = ROOT.TH1F(\"h_cal\", \"Calibrated Spectrum\",\n", + " hist.GetNbinsX(), emin, emax)\n", + "\n", + "for i in range(1, hist.GetNbinsX() + 1):\n", + " counts = hist.GetBinContent(i)\n", + " h_cal.SetBinContent(i, counts)\n", + "\n", + "h_cal.GetXaxis().SetTitle(\"Energy (MeV)\")\n", + "\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "cv_cal = ROOT.TCanvas('cv_cal', 'cv_cal', 1600, 1200)\n", + "\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lgds = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "\n", + "\n", + "Na22a_fit_cal = ROOT.TF1(\"Na22a_fit_cal\", \"gaus\", 0.3, 0.6)\n", + "Na22b_fit_cal = ROOT.TF1(\"Na22b_fit_cal\", \"gaus\", 1.1 , 1.5)\n", + "Na22c_fit_cal = ROOT.TF1(\"Na22c_fit_cal\", \"gaus\", 1.5, 1.9)\n", + "h_cal.Fit(Na22a_fit_cal, \"R+S\")\n", + "h_cal.Fit(Na22b_fit_cal, \"R+S\")\n", + "h_cal.Fit(Na22c_fit_cal, \"R+S\")\n", + "\n", + "h_cal.Draw(\"hist\")\n", + "cv_cal.cd(idx+1).SetLogy()\n", + "\n", + "cv_cal.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8947e352-94db-4b13-8666-c1f2d7eca282", + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/energyResolution-Copy2.ipynb b/notebooks/energyResolution-Copy2.ipynb new file mode 100644 index 0000000..9e22d77 --- /dev/null +++ b/notebooks/energyResolution-Copy2.ipynb @@ -0,0 +1,929 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4a3b1bc5-ec73-4d69-a840-07d057f4e2eb", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "project_root=Path(project_root)\n", + "\n", + "from processing import SiphraAcquisition, MetadataLoader\n", + "from analysis import fit_peak_expbg\n", + "\n", + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "\n", + "# Energy emission peaks in MeV\n", + "Cs137_MeV = 0.662\n", + "Na22_MeV = [0.511, 1.275, 1.786]\n", + "Co60_MeV = [1.173, 1.332, 2,505]\n", + "Am241_MeV = 0.0595\n", + "# Background emission\n", + "K40_MeV = 1.460\n", + "Tl208_MEV = 2.614\n", + "\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "75493f51-2089-4f0a-8057-0f2c2bf7466f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[SIPHRAAcquisition(File: 'SUBT_01_EnergyResolution_thr30_gain1over100_1over3_Background'), SIPHRAAcquisition(File: 'SUBT_02_EnergyResolution_thr30_gain1over100_1over3_Cs137'), SIPHRAAcquisition(File: 'SUBT_03_EnergyResolution_thr30_gain1over100_1over3_Background2'), SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22'), SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4'), SIPHRAAcquisition(File: 'SUBT_06_EnergyResolution_thr30_gain1over100_1over3_Co60'), SIPHRAAcquisition(File: 'SUBT_07_EnergyResolution_thr30_gain1over100_1over3_Background6'), SIPHRAAcquisition(File: 'SUBT_08_EnergyResolution_thr30_gain1over100_1over3_Am241'), SIPHRAAcquisition(File: 'SUBT_09_EnergyResolution_thr30_gain1over100_1over3_Background8')]\n" + ] + } + ], + "source": [ + "files = sorted((project_root/'data/260323').glob('SUBT_*'))\n", + "acqs = [SiphraAcquisition(file) for file in files]\n", + "print(acqs)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "88f4c856-78c0-4ca6-ad58-cc10f3b53649", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cs-137\n", + "Na22\n", + "Co-60\n", + "Am-241\n" + ] + } + ], + "source": [ + "# histograms\n", + "hists = {}\n", + "sources = []\n", + "for sgnl, bg in zip(acqs[1::2], acqs[2::2]):\n", + " filepath = sgnl.filepath.stem\n", + " src = (MetadataLoader.load(sgnl.metadataFile)).source\n", + " sources.append(src)\n", + " print(src)\n", + " # Signal and Background\n", + " hist_sgnl = ROOT.TH1F(f\"{src} Signal\", \"\", BITS12, 0, BITS12)\n", + " hist_bg = ROOT.TH1F(f\"{src} Background\", \"\", BITS12, 0 , BITS12)\n", + " hist_sgnl.Fill(sgnl['s']/len(sgnl.active_chs))\n", + " hist_bg.Fill(bg['s']/len(bg.active_chs))\n", + " hist_bg.Scale(sgnl.exposure/bg.exposure)\n", + " # Clean spectrum\n", + " hist_clean = hist_sgnl.Clone(f\"{src} Clean\")\n", + " hist_clean.Add(hist_bg, -1)\n", + " for hist in [hist_sgnl, hist_bg, hist_clean]:\n", + " # hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Normalized counts\")\n", + " hists[src] = hist_sgnl\n", + " hists[f\"{src}_BG\"] = hist_bg\n", + " hists[f\"{src}_clean\"] = hist_clean\n", + " del hist_sgnl, hist_bg\n", + "#print(hists)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "36f654a6-bd64-471f-a0e3-01927db52e3e", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1600, 1200)\n", + "cv.Divide(2,2)\n", + "\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lgds = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "for idx, (src, lg) in enumerate(zip(sources, lgds)): \n", + " cv.cd(idx+1)\n", + " \n", + " sgnl = hists[src]\n", + " bg = hists[src+'_BG']\n", + " clean = hists[src+'_clean']\n", + " \n", + " lg.AddEntry(sgnl, \"Signal\")\n", + " lg.AddEntry(bg, \"Background\")\n", + " lg.AddEntry(clean, \"Subtracted\")\n", + " sgnl.SetLineColor(colors[0])\n", + " bg.SetLineColor(colors[1])\n", + " clean.SetLineColor(colors[2])\n", + " sgnl.SetTitle(src)\n", + " sgnl.Draw(\"hist\")\n", + " bg.Draw(\"sames hist\")\n", + " clean.Draw(\"sames hist\")\n", + " ROOT.gPad.Update()\n", + " for i, sp in enumerate([sgnl, bg, clean]):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + " lg.Draw()\n", + " cv.cd(idx+1).SetLogy()\n", + " cv.Draw()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6db7ba78-ddd9-4d76-8bba-bf310f8f1221", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.004280158214429103\n", + "-0.2655932438428593\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 196.518\n", + "NDf = 127\n", + "Edm = 2.29967e-07\n", + "NCalls = 65\n", + "Constant = 2459.66 +/- 7.74946 \n", + "Mean = 183.993 +/- 0.078652 \n", + "Sigma = 29.2438 +/- 0.0713573 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 356.977\n", + "NDf = 147\n", + "Edm = 1.20533e-06\n", + "NCalls = 64\n", + "Constant = 1361.97 +/- 6.38965 \n", + "Mean = 275.54 +/- 0.144272 \n", + "Sigma = 32.7382 +/- 0.141925 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 52.0578\n", + "NDf = 47\n", + "Edm = 1.53695e-06\n", + "NCalls = 107\n", + "Constant = 255.901 +/- 3.65957 \n", + "Mean = 478.856 +/- 0.861944 \n", + "Sigma = 32.2107 +/- 1.8725 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 50.431\n", + "NDf = 47\n", + "Edm = 6.26293e-06\n", + "NCalls = 110\n", + "Constant = 187.956 +/- 4.04943 \n", + "Mean = 177.384 +/- 1.67095 \n", + "Sigma = 37.5292 +/- 4.35406 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 152.171\n", + "NDf = 97\n", + "Edm = 2.80563e-08\n", + "NCalls = 96\n", + "Constant = 606.675 +/- 4.80363 \n", + "Mean = 371.93 +/- 0.452134 \n", + "Sigma = 45.3973 +/- 0.782852 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 63.6754\n", + "NDf = 67\n", + "Edm = 1.95162e-06\n", + "NCalls = 89\n", + "Constant = 880.588 +/- 5.38234 \n", + "Mean = 471.391 +/- 0.435139 \n", + "Sigma = 32.0034 +/- 0.470993 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 0.00408946\n", + "NDf = 1\n", + "p0 = -0.265593 +/- 0.109303 \n", + "p1 = 0.00428016 +/- 0.000302376 \n" + ] + } + ], + "source": [ + "#Fitting all the visible peaks. Values around peak input manually\n", + "\n", + "Cs137_fit = ROOT.TF1(\"Cs137_fit\", \"gaus\", 120,250)\n", + "\n", + "Co60a_fit = ROOT.TF1(\"Co60a_fit\", \"gaus\", 200,350)\n", + "Co60b_fit = ROOT.TF1(\"Co60b_fit\", \"gaus\", 450,500)\n", + "\n", + "Na22a_fit = ROOT.TF1(\"Na22a_fit\", \"gaus\", 150,200)\n", + "Na22b_fit = ROOT.TF1(\"Na22b_fit\", \"gaus\", 320,420)\n", + "Na22c_fit = ROOT.TF1(\"Na22c_fit\", \"gaus\", 450,520)\n", + "\n", + "fits = [Cs137_fit, Co60a_fit, Co60b_fit, Na22a_fit, Na22b_fit, Na22c_fit]\n", + "\n", + "#Fit to histogram and find mean value of peak\n", + "hists['Cs-137_clean'].Fit(Cs137_fit, \"R+S\")\n", + "hists['Co-60_clean'].Fit(Co60a_fit, \"R+S\")\n", + "hists['Co-60_clean'].Fit(Co60b_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22a_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22b_fit, \"R+S\")\n", + "hists['Na22_clean'].Fit(Na22c_fit, \"R+S\")\n", + "\n", + "Cs137_mean = Cs137_fit.GetParameter('Mean')\n", + "Co60a_mean = Co60a_fit.GetParameter('Mean')\n", + "Co60b_mean = Co60b_fit.GetParameter('Mean')\n", + "Na22a_mean = Na22a_fit.GetParameter('Mean')\n", + "Na22b_mean = Na22b_fit.GetParameter('Mean')\n", + "Na22c_mean = Na22c_fit.GetParameter('Mean')\n", + "\n", + "channels = [Na22a_mean, Na22b_mean, Na22c_mean]\n", + "energies = Na22_MeV\n", + "channels = np.array(channels, dtype='float64')\n", + "energies = np.array(energies, dtype='float64')\n", + "\n", + "g = ROOT.TGraph(len(channels), channels, energies)\n", + "\n", + "f = ROOT.TF1(\"calib\", \"pol1\", min(channels), max(channels))\n", + "g.Fit(f)\n", + "\n", + "a = f.GetParameter(1)\n", + "b = f.GetParameter(0)\n", + "print(a)\n", + "print(b)\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()\n", + "\n", + "#lg.Draw()\n", + "#cv.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "6b5a5060-fe8f-48ba-af47-c12bca5ea672", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-0.2655932438428593\n", + "-0.2655932438428593\n", + "17.265934802458748\n", + "17.265934802458748\n", + "SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22')\n", + "SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4')\n", + "[0.96495224 0.7338237 0.70760773 ... 0.76057469 0.40532156 0.49787998]\n", + "16.08942631326755\n", + "-0.11364762723062616\n", + "4096\n", + "4096\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv_cal'):\n", + " ROOT.gROOT.FindObject('cv_cal').Close()\n", + "\n", + "\n", + " \n", + "hist = hists['Na22_clean']\n", + "\n", + "emin = a * hist.GetXaxis().GetXmin() + b\n", + "emax = a * hist.GetXaxis().GetXmax() + b\n", + "\n", + "emax_2 = a * BITS12 + b\n", + "emin_2 = a * 0 + b\n", + "\n", + "print(emin)\n", + "print(emin_2)\n", + "print(emax)\n", + "print(emax_2)\n", + "\n", + "print(acqs[3])\n", + "print(acqs[4])\n", + "\n", + "data_bg_clb = a * (acqs[4]['s']/len(acqs[4].active_chs)) + b\n", + "data_src_clb = a * (acqs[3]['s']/len(acqs[4].active_chs)) + b\n", + "\n", + "print(data_bg_clb)\n", + "\n", + "print(data_bg_clb.max())\n", + "print(data_bg_clb.min())\n", + "\n", + "h_cal = ROOT.TH1F(\"h_cal\", \"Calibrated Spectrum\",\n", + " hist.GetNbinsX(), emin, emax)\n", + "\n", + "print(hist.GetNbinsX())\n", + "print(BITS12)\n", + "\n", + "#for i in range(1, hist.GetNbinsX() + 1):\n", + "# counts = hist.GetBinContent(i)\n", + "# h_cal.SetBinContent(i, counts)\n", + "\n", + "hist_bg_clb = ROOT.TH1F(\"Calibrated background\", \"Calibrated Histograms\", BITS12, emin, emax)\n", + "hist_src_clb = ROOT.TH1F(\"Calibrated signal\", \"Calibrated Histograms\", BITS12, emin, emax)\n", + "\n", + "\n", + "hist_bg_clb.Fill(data_bg_clb)\n", + "#hist_bg_clb.Scale(NORMFACTOR)\n", + "\n", + "hist_src_clb.Fill(data_src_clb)\n", + "\n", + "hist_clean_clb = hist_src_clb.Clone(\"Calibrated signal no BG\")\n", + "hist_clean_clb.Add(hist_bg_clb, -1)\n", + "\n", + "#for hist, clr in zip([hist_bg_clb, hist_src_clb, hist_clean_clb], colors):\n", + "# hist.GetXaxis().SetTitle(\"Energy (MeV)\")\n", + " # hist.GetYaxis().SetTitle(\"Normalized counts\")\n", + " # hist.SetLineColor(clr) \n", + "#\n", + "#cv_cal.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "57b1d3c2-9d5b-43c3-8b49-fd081e3fc77f", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Error in : Cannot set Y axis to log scale\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot all the histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "hist_bg_clb.SetLineColor(colors[0])\n", + "hist_src_clb.SetLineColor(colors[1])\n", + "lg1.AddEntry(hist_bg_clb, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_src_clb, r\"Signal ^{137}Cs\", \"l\")\n", + "hist_bg_clb.Draw(\"hist\")\n", + "hist_src_clb.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "hist_bg_clb.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "# lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_clean_clb.SetLineColor(colors[2])\n", + "hist_clean_clb.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "79a219b6-cf1a-419e-825f-272accfa28eb", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Error in : Cannot set Y axis to log scale\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 137.425\n", + "NDf = 67\n", + "Edm = 2.14133e-06\n", + "NCalls = 84\n", + "Constant = 190.012 +/- 2.40871 \n", + "Mean = 0.492756 +/- 0.00341755 \n", + "Sigma = 0.142388 +/- 0.00438193 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 465.001\n", + "NDf = 91\n", + "Edm = 1.22703e-07\n", + "NCalls = 75\n", + "Constant = 611.78 +/- 3.97344 \n", + "Mean = 1.31966 +/- 0.00146494 \n", + "Sigma = 0.178741 +/- 0.0023569 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 742.699\n", + "NDf = 90\n", + "Edm = 7.47939e-06\n", + "NCalls = 95\n", + "Constant = 827.73 +/- 4.69558 \n", + "Mean = 1.74335 +/- 0.00131288 \n", + "Sigma = 0.172477 +/- 0.00195917 \t (limited)\n" + ] + } + ], + "source": [ + "# Plot all the histograms\n", + "\n", + "\n", + "\n", + "\n", + "Na22a_fit_cal = ROOT.TF1(\"Na22a_fit_cal\", \"gaus\", 0.3, 0.6)\n", + "Na22b_fit_cal = ROOT.TF1(\"Na22b_fit_cal\", \"gaus\", 1.1 , 1.5)\n", + "Na22c_fit_cal = ROOT.TF1(\"Na22c_fit_cal\", \"gaus\", 1.5, 1.9)\n", + "hist_clean_clb.Fit(Na22a_fit_cal, \"R+S\")\n", + "hist_clean_clb.Fit(Na22b_fit_cal, \"R+S\")\n", + "hist_clean_clb.Fit(Na22c_fit_cal, \"R+S\")\n", + "\n", + "hist_clean_clb.Draw(\"hist\")\n", + "cv1.cd(idx+1).SetLogy()\n", + "\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "8947e352-94db-4b13-8666-c1f2d7eca282", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[177.38354136 371.93022564 471.39142955]\n", + "[0.511 1.275 1.786]\n", + "[0.004280158214429103, -0.2655932438428593]\n", + "-0.2655932438428593 17.265934802458748\n", + "-0.2655932438428593 17.265934802458748\n", + "SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22')\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 50.431\n", + "NDf = 47\n", + "Edm = 6.26293e-06\n", + "NCalls = 110\n", + "Constant = 187.956 +/- 4.04943 \n", + "Mean = 177.384 +/- 1.67095 \n", + "Sigma = 37.5292 +/- 4.35406 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 152.171\n", + "NDf = 97\n", + "Edm = 2.80563e-08\n", + "NCalls = 96\n", + "Constant = 606.675 +/- 4.80363 \n", + "Mean = 371.93 +/- 0.452134 \n", + "Sigma = 45.3973 +/- 0.782852 \t (limited)\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 63.6754\n", + "NDf = 67\n", + "Edm = 1.95162e-06\n", + "NCalls = 89\n", + "Constant = 880.588 +/- 5.38234 \n", + "Mean = 471.391 +/- 0.435139 \n", + "Sigma = 32.0034 +/- 0.470993 \t (limited)\n", + "****************************************\n", + "Minimizer is Linear / Migrad\n", + "Chi2 = 0.00408946\n", + "NDf = 1\n", + "p0 = -0.265593 +/- 0.109303 \n", + "p1 = 0.00428016 +/- 0.000302376 \n" + ] + } + ], + "source": [ + "def calibration_fit(histogram, energy_ranges, energies):\n", + " \"\"\"Function to create a linear calibration fit based on a histogram. Returns slope and constant of linear fit. \n", + " Inputs: Histogram to base calibration on.\n", + " Ranges within which the peaks of the histogram are located.\n", + " Known energies of these peaks, in MeV.\"\"\"\n", + " \n", + " channels = []\n", + " for i in range(len(energy_ranges)):\n", + " cal_fit=ROOT.TF1(\"cal_fit_\" + str(i), \"gaus\", energy_ranges[i][0], energy_ranges[i][1])\n", + " hist.Fit(cal_fit, \"R+S\")\n", + " channels.append(cal_fit.GetParameter('Mean'))\n", + "\n", + " channels = np.array(channels, dtype='float64')\n", + " energies = np.array(energies, dtype='float64')\n", + "\n", + " print(channels)\n", + " print(energies)\n", + "\n", + " graph = ROOT.TGraph(len(channels), channels, energies)\n", + " fit = ROOT.TF1(\"calib\", \"pol1\", min(channels), max(channels))\n", + " graph.Fit(fit)\n", + "\n", + " \n", + " a = fit.GetParameter(1)\n", + " b = fit.GetParameter(0)\n", + "\n", + " return [a, b]\n", + "\n", + "def calibrated_histogram(linear_fit, acquisition, n_of_bins, name):\n", + " a = linear_fit[0]\n", + " b = linear_fit[1]\n", + " data_cal = a * (acquisition['s']/len(acquisition.active_chs)) + b\n", + "\n", + " emax = a * n_of_bins + b\n", + " emin = b\n", + "\n", + " print(emin, emax)\n", + " hist_cal = ROOT.TH1F(name, \"Calibrated Spectrum\", n_of_bins, emin, emax)\n", + " hist_cal.Fill(data_cal)\n", + "\n", + " return(hist_cal)\n", + "\n", + "def calibrated_acquisition(linear_fit, acquisition):\n", + " a = linear_fit[0]\n", + " b = linear_fit[1]\n", + " print(a, b)\n", + " data_cal = a * (acquisition['s']/len(acquisition.active_chs)) + b\n", + "\n", + " emax = a * BITS12 + b\n", + " emin = b\n", + "\n", + " return(data_cal, emin, emax)\n", + "\n", + " \n", + " \n", + "hist = hists['Na22_clean']\n", + "energy_ranges = [(150, 200), (320, 420), (450, 520)]\n", + "energies = Na22_MeV\n", + "\n", + "c_fit = calibration_fit(hist, energy_ranges, energies)\n", + "print(c_fit)\n", + "hist_cal_bg = calibrated_histogram(c_fit, acqs[4], BITS12, \"hist_cal_bg\")\n", + "hist_cal_bg.Scale(1/acqs[4].exposure) \n", + "hist_cal_Na22 = calibrated_histogram(c_fit, acqs[3], BITS12, \"hist_cal_Na22\")\n", + "hist_cal_Na22.Scale(1/acqs[3].exposure) \n", + "print(acqs[3])\n", + "hist_cal_clean = hist_cal_Na22.Clone(\"Calibrated signal no BG\")\n", + "\n", + "#hist_cal = ROOT.TH1F(\"h_cal\", \"Calibrated Spectrum\", BITS12, emin, emax)\n", + "#hist_cal.Fill(data_cal)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "a9b090a7-1c19-4ea1-8676-176a371a3f62", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot all the histograms\n", + "if ROOT.gROOT.FindObject('cv1'):\n", + " ROOT.gROOT.FindObject('cv1').Close()\n", + "cv1 = ROOT.TCanvas(\"cv1\", \"cv1\", 1600, 800)\n", + "cv1.Divide(2,1)\n", + "\n", + "lg1 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(1)\n", + "hist_cal_bg.SetLineColor(colors[0])\n", + "hist_cal_Na22.SetLineColor(colors[1])\n", + "lg1.AddEntry(hist_cal_bg, \"Background\", \"l\")\n", + "lg1.AddEntry(hist_cal_Na22, r\"Signal ^{22}Na\", \"l\")\n", + "hist_cal_bg.Draw(\"hist\")\n", + "hist_cal_Na22.Draw(\"sames hist\")\n", + "lg1.Draw()\n", + "cv1.cd(1).SetLogy()\n", + "#hist_cal_bg.GetYaxis().SetRangeUser(0,1e4)\n", + "\n", + "\n", + "lg2 = ROOT.TLegend(0.48, 0.71, 0.75, 0.83)\n", + "cv1.cd(2)\n", + "hist_cal_clean.SetLineColor(colors[2])\n", + "hist_cal_clean.Draw(\"hist\")\n", + "cv1.cd(2).SetLogy()\n", + "\n", + "cv1.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "f60a5136-480e-48d5-b6e2-8b32ff255bbb", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.45246703226380525\n", + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "Chi2 = 811.832\n", + "NDf = 90\n", + "Edm = 1.49727e-08\n", + "NCalls = 102\n", + "Constant = 2.92324 +/- 0.0161931 \n", + "Mean = 1.72903 +/- 0.00140192 \n", + "Sigma = 0.19213 +/- 0.00251407 \t (limited)\n" + ] + } + ], + "source": [ + "def energy_Resolution(hist, peak_range):\n", + " resolution_fit = ROOT.TF1(\"res_fit_\" + str(i), \"gaus\", peak_range[0], peak_range[1])\n", + " hist.Fit(resolution_fit, \"R+S\")\n", + " std_dev = resolution_fit.GetParameter('Sigma')\n", + " resolution = 2.355 * std_dev\n", + "\n", + " return resolution \n", + "\n", + "r = energy_Resolution(hist_cal_clean ,(1.5 , 1.9))\n", + "\n", + "print(r)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f972dcb-27f4-4b3b-866d-1c311a9f29fd", + "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.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/energyResolution.ipynb b/notebooks/energyResolution.ipynb new file mode 100644 index 0000000..3656e54 --- /dev/null +++ b/notebooks/energyResolution.ipynb @@ -0,0 +1,307 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4a3b1bc5-ec73-4d69-a840-07d057f4e2eb", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "project_root=Path(project_root)\n", + "\n", + "from processing import SiphraAcquisition, MetadataLoader\n", + "from analysis import fit_peak_expbg\n", + "\n", + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "# Energy emission peaks in MeV\n", + "K40_MEV = [1.460]\n", + "TL208_MEV = [2.614]\n", + "CS137_MEV = [0.661]\n", + "ANNIHIL_MEV = 0.511\n", + "NA22_MEV = [ANNIHIL_MEV, 1.2745, 1.2745+ANNIHIL_MEV]\n", + "CO60_MEV = [1.1732, 1.3325, 1.1732+1.3325]\n", + "\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "75493f51-2089-4f0a-8057-0f2c2bf7466f", + "metadata": {}, + "outputs": [], + "source": [ + "files = sorted((project_root/'data/260323').glob('SUBT_*'))\n", + "acqs = [SiphraAcquisition(file) for file in files]" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "5ee51f11-cd24-44c1-bfe5-f0e3b7d8a85c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[SIPHRAAcquisition(File: 'SUBT_01_EnergyResolution_thr30_gain1over100_1over3_Background'),\n", + " SIPHRAAcquisition(File: 'SUBT_02_EnergyResolution_thr30_gain1over100_1over3_Cs137'),\n", + " SIPHRAAcquisition(File: 'SUBT_03_EnergyResolution_thr30_gain1over100_1over3_Background2'),\n", + " SIPHRAAcquisition(File: 'SUBT_04_EnergyResolution_thr30_gain1over100_1over3_Na22'),\n", + " SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4'),\n", + " SIPHRAAcquisition(File: 'SUBT_06_EnergyResolution_thr30_gain1over100_1over3_Co60'),\n", + " SIPHRAAcquisition(File: 'SUBT_07_EnergyResolution_thr30_gain1over100_1over3_Background6'),\n", + " SIPHRAAcquisition(File: 'SUBT_08_EnergyResolution_thr30_gain1over100_1over3_Am241'),\n", + " SIPHRAAcquisition(File: 'SUBT_09_EnergyResolution_thr30_gain1over100_1over3_Background8')]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "acqs" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "88f4c856-78c0-4ca6-ad58-cc10f3b53649", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cs-137\n", + "Na-22\n", + "Co-60\n", + "Am-241\n", + "{'Cs-137': , 'Cs-137_BG': , 'Cs-137_clean': , 'Na-22': , 'Na-22_BG': , 'Na-22_clean': , 'Co-60': , 'Co-60_BG': , 'Co-60_clean': , 'Am-241': , 'Am-241_BG': , 'Am-241_clean': }\n" + ] + } + ], + "source": [ + "# histograms\n", + "hists = {}\n", + "sources = []\n", + "for sgnl, bg in zip(acqs[1::2], acqs[2::2]):\n", + " filepath = sgnl.filepath.stem\n", + " src = (MetadataLoader.load(sgnl.metadataFile)).source\n", + " sources.append(src)\n", + " print(src)\n", + " # Signal and Background\n", + " hist_sgnl = ROOT.TH1F(f\"{src} Signal\", \"\", BITS12, 0, BITS12)\n", + " hist_bg = ROOT.TH1F(f\"{src} Background\", \"\", BITS12, 0 , BITS12)\n", + " hist_sgnl.Fill(sgnl['s']/len(sgnl.active_chs))\n", + " hist_bg.Fill(bg['s']/len(bg.active_chs))\n", + " hist_bg.Scale(sgnl.exposure/bg.exposure)\n", + " # Clean spectrum\n", + " hist_clean = hist_sgnl.Clone(f\"{src} Clean\")\n", + " hist_clean.Add(hist_bg, -1)\n", + " for hist in [hist_sgnl, hist_bg, hist_clean]:\n", + " # hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Normalized counts\")\n", + " hists[src] = hist_sgnl\n", + " hists[f\"{src}_BG\"] = hist_bg\n", + " hists[f\"{src}_clean\"] = hist_clean\n", + " del hist_sgnl, hist_bg\n", + "print(hists)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "36f654a6-bd64-471f-a0e3-01927db52e3e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "Yinit = 0.82 # For stat boxes\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1600, 1200)\n", + "cv.Divide(2,2)\n", + "\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lgds = [ROOT.TLegend(0.48, 0.61, 0.75, 0.83) for _ in range(4)]\n", + "\n", + "for idx, (src, lg) in enumerate(zip(sources, lgds)): \n", + " cv.cd(idx+1)\n", + " \n", + " sgnl = hists[src]\n", + " bg = hists[src+'_BG']\n", + " clean = hists[src+'_clean']\n", + " \n", + " lg.AddEntry(sgnl, \"Signal\")\n", + " lg.AddEntry(bg, \"Background\")\n", + " lg.AddEntry(clean, \"Subtracted\")\n", + " sgnl.SetLineColor(colors[0])\n", + " bg.SetLineColor(colors[1])\n", + " clean.SetLineColor(colors[2])\n", + " sgnl.SetTitle(src)\n", + " sgnl.Draw(\"hist\")\n", + " bg.Draw(\"sames hist\")\n", + " clean.Draw(\"sames hist\")\n", + " ROOT.gPad.Update()\n", + " for i, sp in enumerate([sgnl, bg, clean]):\n", + " st = sp.FindObject(\"stats\")\n", + " # print(type(st))\n", + " st.SetY1NDC(Yinit - i*0.08)\n", + " st.SetY2NDC(0.1)\n", + " lg.Draw()\n", + " cv.cd(idx+1).SetLogy()\n", + " cv.Draw()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "6db7ba78-ddd9-4d76-8bba-bf310f8f1221", + "metadata": {}, + "outputs": [], + "source": [ + "fit_xrange = {\n", + " 'Cs-137': [(100, 270)],\n", + " 'Na-22': [(325, 408), (436,508)],\n", + " 'Co-60': [(200, 361), (419, 520)],\n", + " 'Am-241': [(8, 30)]\n", + "}\n", + "\n", + "# Fit Cs-137\n" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "6b5a5060-fe8f-48ba-af47-c12bca5ea672", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[SIPHRAAcquisition(File: 'SUBT_03_EnergyResolution_thr30_gain1over100_1over3_Background2'),\n", + " SIPHRAAcquisition(File: 'SUBT_05_EnergyResolution_thr30_gain1over100_1over3_Background4'),\n", + " SIPHRAAcquisition(File: 'SUBT_07_EnergyResolution_thr30_gain1over100_1over3_Background6'),\n", + " SIPHRAAcquisition(File: 'SUBT_09_EnergyResolution_thr30_gain1over100_1over3_Background8')]" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "acqs[2::2]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (ROOT in Conda)", + "language": "python", + "name": "comcube" + }, + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/noteboks/etcTests/converters_explore.ipynb b/notebooks/etcTests/converters_explore.ipynb similarity index 100% rename from noteboks/etcTests/converters_explore.ipynb rename to notebooks/etcTests/converters_explore.ipynb diff --git a/noteboks/etcTests/siphraAcq_test_json.ipynb b/notebooks/etcTests/siphraAcq_test_json.ipynb similarity index 100% rename from noteboks/etcTests/siphraAcq_test_json.ipynb rename to notebooks/etcTests/siphraAcq_test_json.ipynb diff --git a/notebooks/gainTesting.ipynb b/notebooks/gainTesting.ipynb new file mode 100644 index 0000000..7aa8f5b --- /dev/null +++ b/notebooks/gainTesting.ipynb @@ -0,0 +1,687 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 66, + "id": "cc81fa7c-d43e-4c17-854d-8bc51e7bdf91", + "metadata": {}, + "outputs": [], + "source": [ + "import ROOT\n", + "import pandas as pd\n", + "import numpy as np\n", + "import os\n", + "import sys\n", + "import ipynbname\n", + "from pathlib import Path\n", + "\n", + "project_root = str(ipynbname.path().parent.parent)\n", + "sys.path.append(project_root)\n", + "project_root=Path(project_root)\n", + "\n", + "from processing import SiphraAcquisition\n", + "from analysis import fit_peak_expbg" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "id": "f14c60fa-8e5a-4c1b-bd05-6d8f0a23851d", + "metadata": {}, + "outputs": [], + "source": [ + "# Constants\n", + "BITS12 = 2**12\n", + "BITS9 = 2**9 # 512 typical number of bins used\n", + "# Energy emission peaks in MeV\n", + "K40_MEV = 1.460\n", + "TL208_MEV = 2.614\n", + "CS137_MEV = 0.661\n", + "colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kViolet, ROOT.kYellow, ROOT.kSpring, ROOT.kCyan,]" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "id": "9d972c1c-3934-4a10-9e7c-4e7b9caa4b3f", + "metadata": {}, + "outputs": [], + "source": [ + "files1 = sorted((project_root/'data/260317').glob('1*.csv'))\n", + "acqs1 = [SiphraAcquisition(file) for file in files1]" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "id": "2696fa6e-8416-4df7-9bf9-62d00adcee3b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[SIPHRAAcquisition(File: '10_GainTesting_SiPM_Vbias29,0V_All8Chs_Na22_SIPHRA1'),\n", + " SIPHRAAcquisition(File: '11_GainTesting_SiPM_Vbias29,0V_All8Chs_Na22_SIPHRA1'),\n", + " SIPHRAAcquisition(File: '12_GainTesting_SiPM_Vbias29,0V_All8Chs_Na22_SIPHRA1'),\n", + " SIPHRAAcquisition(File: '15_GainTesting_SiPM_Vbias29,0V_All8Chs_Na22_SIPHRA1'),\n", + " SIPHRAAcquisition(File: '16_GainTesting_SiPM_Vbias29,0V_All8Chs_Background_SIPHRA1'),\n", + " SIPHRAAcquisition(File: '17_GainTesting_SiPM_Vbias29,0V_All8Chs_Background_SIPHRA1'),\n", + " SIPHRAAcquisition(File: '18_GainTesting_SiPM_Vbias29,0V_All8Chs_Background_SIPHRA1')]" + ] + }, + "execution_count": 69, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "acqs1.pop(3)\n", + "acqs1.pop(3)\n", + "acqs1" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "id": "af88b5d2-3d88-4310-ad81-90986efb6ffd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'10'" + ] + }, + "execution_count": 70, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "acqs1[0].filepath.stem[:2]" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "id": "8d9cd696-eb5e-46c8-9673-73e2afc2ebc5", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: cmis_gain 1/10 (Potential memory leak).\n", + "Warning in : Replacing existing TH1: cmis_gain 1/100 (Potential memory leak).\n", + "Warning in : Replacing existing TH1: cmis_gain 1/200 (Potential memory leak).\n", + "Warning in : Replacing existing TH1: cmis_gain 1/400 (Potential memory leak).\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1200, 800)\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lg = ROOT.TLegend(0.58, 0.61, 0.75, 0.83)\n", + "lg.SetHeader('CMIS gain')\n", + "\n", + "ci_gain='1V/0.75nC'\n", + "cmis_gains=['1/10', '1/100', '1/200', '1/400']\n", + "hists={}\n", + "\n", + "for idx, (acq, cmsg) in enumerate(zip(acqs1[:4], cmis_gains)):\n", + " hist = ROOT.TH1F(f\"cmis_gain {cmsg}\", f\"{cmsg}\", BITS12, 0, BITS12)\n", + " hist.Fill(acq['s']/len(acq.active_chs))\n", + " hist.Scale(1/acq.exposure)\n", + " # Preeting up ...\n", + " hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Count rate (s^{-1})\")\n", + " hist.SetLineColor(colors[idx]+2)\n", + " # legend and storing\n", + " lg.AddEntry(hist, f\"{cmsg}\")\n", + " hists[f\"Acq{acq.filepath.stem[:2]}\"] = hist\n", + " del hist\n", + " hists[f\"Acq{acq.filepath.stem[:2]}\"].Draw('sames hist')\n", + "cv.SetLogy()\n", + "lg.Draw()\n", + "cv.Draw()\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "id": "65e8172d-a2bc-4306-87b7-66027afc15f4", + "metadata": {}, + "outputs": [], + "source": [ + "# Calibration with acquisition 15(Signal) and 16 (BG)\n", + "files2 = sorted([p for p in (project_root/'data/260317').iterdir() if p.stem[11:13] in ['15', '16']])\n", + "acqs2 = [SiphraAcquisition(file) for file in files2]" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "b831a930-27c1-442f-8dd0-e30f2d8a84a1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[PosixPath('/home/oscar/Projects/SIPHRA_Testing/data/260317/SUBTRACTED_15_GainTesting_SiPM_Vbias29,0V_All8Chs_Na22_SIPHRA1.csv'),\n", + " PosixPath('/home/oscar/Projects/SIPHRA_Testing/data/260317/SUBTRACTED_16_GainTesting_SiPM_Vbias29,0V_All8Chs_Background_SIPHRA1.csv')]" + ] + }, + "execution_count": 73, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "files2" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "id": "0fcebaeb-0ace-4f57-b151-3a80a51535e9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "8" + ] + }, + "execution_count": 74, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(acqs2[0].active_chs)" + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "id": "f6a83f88-b203-4af0-9f0b-18a7b5a0338f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exposures:\n", + " Signal:\t220.216 s\n", + " Background:\t625.392 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning in : Replacing existing TH1: Signal (Potential memory leak).\n", + "Warning in : Replacing existing TH1: Background (Potential memory leak).\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1200, 800)\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lg = ROOT.TLegend(0.58, 0.61, 0.75, 0.83)\n", + "# lg.SetHeader('')\n", + "\n", + "ci_gain='1V/0.75nC'\n", + "labels=['Signal', 'Background']\n", + "\n", + "ref_exposure = acqs2[0].exposure\n", + "\n", + "for idx, (acq, label) in enumerate(zip(acqs2, labels)):\n", + " hist = ROOT.TH1F(label, label, BITS12, 0, BITS12)\n", + " hist.Fill(acq['s']/len(acq.active_chs))\n", + " # hist.Scale(ref_exposure/acq.exposure)\n", + " # Preeting up ...\n", + " hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC, CMIS gain = 1/400\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Counts\")\n", + " hist.SetLineColor(colors[idx]+2)\n", + " # legend and storing\n", + " lg.AddEntry(hist, label)\n", + " hists[f\"Acq{acq.filepath.stem[:13]}\"] = hist\n", + " del hist\n", + " hists[f\"Acq{acq.filepath.stem[:13]}\"].Draw('sames hist')\n", + "\n", + "# hist_subt = hists['AcqSUBTRACTED_15'].Clone('Subtracted')\n", + "# hist_subt.Add(hists['AcqSUBTRACTED_16'], -1)\n", + "# hist_subt.Draw('sames hist')\n", + "# hist_subt.SetLineColor(colors[2])\n", + "# lg.AddEntry(hist_subt, 'Subtracted')\n", + "\n", + "cv.SetLogy()\n", + "lg.Draw()\n", + "cv.Draw()\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()\n", + "print(f\"Exposures:\\n Signal:\\t{acqs2[0].exposure:>.3f} s\\n Background:\\t{acqs2[1].exposure:>.3f} s\")" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "05919cfa-b075-4898-8f7b-56039d4c8da1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'Acq10': ,\n", + " 'Acq11': ,\n", + " 'Acq12': ,\n", + " 'Acq15': ,\n", + " 'AcqSUBTRACTED_10': ,\n", + " 'AcqSUBTRACTED_11': ,\n", + " 'AcqSUBTRACTED_16': ,\n", + " 'AcqSUBTRACTED_15': }" + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hists" + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "id": "8e7a42e8-c3a4-451c-b635-5d37c1167faf", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "****************************************\n", + "Minimizer is Minuit2 / Migrad\n", + "MinFCN = 7605.91\n", + "Chi2 = 15211.8\n", + "NDf = 197\n", + "Edm = 1.10191e-06\n", + "NCalls = 138\n", + "Constant = 1267.88 +/- 5.7826 \n", + "Mean = 172.076 +/- 0.0854468 \n", + "Sigma = 22.9615 +/- 0.0616869 \t (limited)\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "
\n", + "
\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "files2 = sorted((project_root/'data/260317').glob('SUBTRACTED_*'))\n", + "acqs2 = [SiphraAcquisition(file) for file in files2]\n", + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1200, 800)\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lg = ROOT.TLegend(0.58, 0.61, 0.75, 0.83)\n", + "\n", + "hist = hist_subt.Clone('hist_subt_calib')\n", + "hist.SetTitle(r'^{22}Na subtracted spectrum, CI gain = 1/0.75 pC, CMIS gain = 1/400')\n", + "fit1 = hist.Fit('gaus', 'L S', \"\", 100, 300)\n", + "hist.Draw()\n", + "cv.SetLogy()\n", + "cv.Draw()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "daf8a1de-c392-4879-a4a9-619db60f669e", + "metadata": {}, + "outputs": [], + "source": [ + "# Calibration with acquisition 15(Signal) and 16 (BG)\n", + "files2 = sorted((project_root/'data/260317').glob('SUBTRACTED_*'))\n", + "acqs2 = [SiphraAcquisition(file) for file in files2]\n", + "if ROOT.gROOT.FindObject('cv'):\n", + " ROOT.gROOT.FindObject('cv').Close()\n", + "\n", + "cv = ROOT.TCanvas('cv', 'cv', 1200, 800)\n", + "ROOT.gStyle.SetOptStat(11)\n", + "ROOT.gStyle.SetStatFontSize(0.03)\n", + "ROOT.gStyle.SetStatW(0.16)\n", + "\n", + "lg = ROOT.TLegend(0.58, 0.61, 0.75, 0.83)\n", + "# lg.SetHeader('CMIS gain')\n", + "\n", + "ci_gain='1V/0.75nC'\n", + "labels=['Signal', 'Background']\n", + "hists={}\n", + "\n", + "for idx, (acq, label) in enumerate(zip(acqs2, labels)):\n", + " hist = ROOT.TH1F(label, label, BITS12, 0, BITS12)\n", + " hist.Fill(acq['s']/len(acq.active_chs))\n", + " hist.Scale(1/acq.exposure)\n", + " # Preeting up ...\n", + " hist.SetTitle(r\"^{22}Na spectra - CI gain = 1/0.75 pC, CMIS gain = 1/400\")\n", + " hist.GetXaxis().SetTitle(\"ADC channel number\")\n", + " hist.GetYaxis().SetTitle(r\"Count rate (s^{-1})\")\n", + " hist.SetLineColor(colors[idx]+2)\n", + " # legend and storing\n", + " lg.AddEntry(hist, label)\n", + " hists[f\"Acq{acq.filepath.stem[:13]}\"] = hist\n", + " del hist\n", + " hists[f\"Acq{acq.filepath.stem[:13]}\"].Draw('sames hist')\n", + "\n", + "\n", + "\n", + "cv.SetLogy()\n", + "lg.Draw()\n", + "cv.Draw()\n", + "ROOT.gPad.Modified()\n", + "ROOT.gPad.Update()" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "2028d39f-0e60-4bbe-b3c9-55cbbc63b14b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2004" + ] + }, + "execution_count": 56, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "2278 - 274" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "f3576c1a-50d1-4778-920f-1dd754dd473d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1204" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "1841 - 637" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "a74af979-2cdc-4da3-8164-59ab0a3be28d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "560" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "984-424" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (ROOT in Conda)", + "language": "python", + "name": "comcube" + }, + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/noteboks/histogram_summed.ipynb b/notebooks/histogram_summed.ipynb similarity index 100% rename from noteboks/histogram_summed.ipynb rename to notebooks/histogram_summed.ipynb diff --git a/noteboks/img/etc/CMIS_internal.png b/notebooks/img/etc/CMIS_internal.png similarity index 100% rename from noteboks/img/etc/CMIS_internal.png rename to notebooks/img/etc/CMIS_internal.png diff --git a/noteboks/img/etc/SIPHRA_arch.png b/notebooks/img/etc/SIPHRA_arch.png similarity index 100% rename from noteboks/img/etc/SIPHRA_arch.png rename to notebooks/img/etc/SIPHRA_arch.png diff --git a/noteboks/include_project_root.py b/notebooks/include_project_root.py similarity index 100% rename from noteboks/include_project_root.py rename to notebooks/include_project_root.py diff --git a/noteboks/spectra_current_thrsh.ipynb b/notebooks/spectra_current_thrsh.ipynb similarity index 100% rename from noteboks/spectra_current_thrsh.ipynb rename to notebooks/spectra_current_thrsh.ipynb diff --git a/noteboks/spectrum_plotter.ipynb b/notebooks/spectrum_plotter.ipynb similarity index 100% rename from noteboks/spectrum_plotter.ipynb rename to notebooks/spectrum_plotter.ipynb diff --git a/processing/__pycache__/__init__.cpython-314.pyc b/processing/__pycache__/__init__.cpython-314.pyc deleted file mode 100644 index 2e2ac3f..0000000 Binary files a/processing/__pycache__/__init__.cpython-314.pyc and /dev/null differ diff --git a/processing/__pycache__/metadata.cpython-314.pyc b/processing/__pycache__/metadata.cpython-314.pyc deleted file mode 100644 index 7e1887b..0000000 Binary files a/processing/__pycache__/metadata.cpython-314.pyc and /dev/null differ diff --git a/processing/__pycache__/siphraacquisition.cpython-314.pyc b/processing/__pycache__/siphraacquisition.cpython-314.pyc deleted file mode 100644 index 54d32b4..0000000 Binary files a/processing/__pycache__/siphraacquisition.cpython-314.pyc and /dev/null differ diff --git a/processing/metadata.py b/processing/metadata.py index 84c3501..8321d42 100644 --- a/processing/metadata.py +++ b/processing/metadata.py @@ -11,6 +11,7 @@ class Metadata: active_chs: list[int] sipm_chs: str | list[str] n_events: int + source: str | None class MetadataLoader: @@ -29,10 +30,12 @@ def load(path: PathLike) -> Metadata: @staticmethod def _parse_v1(raw: dict) -> Metadata: acq = raw["acquisition"] + source = raw["source"]["type"] return Metadata( exposure_sec=acq["exposure_sec"], active_chs=[int(ch) for ch in acq["active_chs"]], sipm_chs=acq["sipm_chs"], n_events=acq["counts"], + source=source, ) \ No newline at end of file diff --git a/siphractl/__pycache__/d2a_lib.cpython-314.pyc b/siphractl/__pycache__/d2a_lib.cpython-314.pyc deleted file mode 100644 index ffc7aa0..0000000 Binary files a/siphractl/__pycache__/d2a_lib.cpython-314.pyc and /dev/null differ diff --git a/siphractl/siphra_controller.py b/siphractl/siphra_controller.py deleted file mode 100644 index 7c2cbbb..0000000 --- a/siphractl/siphra_controller.py +++ /dev/null @@ -1,112 +0,0 @@ -from .d2a_lib import * -from .regs_bit_structure import * - -CH_ADDRS = [] -SIPHRA_RETURN_SIZE = 15 # bytes. Returned by SIPHRA when D2a.readSIPHRA is called -# TODO: Check the size of the packet returned by D2a.readSIPHRA -REG_SIZE = 4 # bytes. Size of the bytearray passed to D2a.writeSIPHRA - -# Register map classes - -class SIPHRARegister: - def __init__(self, addr, structure): - self.addr = addr - self.structure = structure - - def __getitem__(self, param_name): - pass - - def __contains__(self, param_name): - '''True if the register has a parameter named parameter.''' - params = self.structure.subcon.subcons - param_names = [param.name for param in params if param.name] - return param_name in param_names - - def parse(self, value): - return self.structure.parse(value) - - def set_param(self, param_name, value, reg_current_value): - old_register = parse(reg_current_value) - - -class SIPHRARegMap: - def __init__(self): - self._registers = {} - - for i in range(1,17): - self._registers[f"ch{i}"] = SIPHRARegister(i, CHANNEL) - - self._registers["summ"] = SIPHRARegister(0x10, SUMM_CHANNEL) - self._registers["ch_config"] = SIPHRARegister(0x11, CHANNEL_CONFIG) - - def __getitem__(self, name): - return self._registers[name] - - def __getattr__(self, name): - return self._registers[name] - - def __iter__(self): - return iter(self._registers.values()) - - def __len__(self): - return len(self._registers) - - def __contains__(self, name): - return name in self._registers - - def items(self): - return self._registers.items() - - def keys(self): - return self._registers.keys() - - def values(self): - return self._registers.values() - - -class SIPHRA: - def __init__(self, d2a: D2a): - self._d2a = d2a - self._regs = SIPHRARegMap() - - def get_reg_value(self, reg_name, chip='A'): - addr = self._regs[reg_name].addr - reg_value = self._d2a.readSIPHRA(addr, chip)[SIPHRA_RETURN_SIZE - REG_SIZE:] - return reg_value - - def read_register(self, reg_name, chip='A'): - reg_value = self.get_reg_value(reg_name, chip) - strct = self._regs[reg_name].structure - parsed_struct = strct.parse(reg_value) - return parsed_struct - - def _find_reg_containing_param(self, parameter, ch: int=0): - ''' - Returns the address of the register containing the given parameter. - :param parameter: Name of the desired parameter. - :param ch: If the parameter is in one of the 16+1 channel addresses, this argument is used for disambiguation. ``ch`` is a number between 1 and 16 if the parameter belongs to one of the input channels, it is 17 if the parameter belongs to the summing channel and is defaulted to 0 if the parameter belongs to any other register. - :return: The address of the register containing the given parameter. - ''' - reg_addr = None - if not 0 < ch < 17: - raise ValueError(f"Channel {ch} is out of range. Use channels 1-16 and 17 for summing channel") - if ch == 0: # Not a channel register - for reg in self._regs: - if parameter in reg: - reg_addr = reg.addr - else: # Channel register - # Verify that the given register exists - register = self._regs[f"ch{ch}"] - if parameter in register: reg_addr = register.addr - if not reg_addr: - raise NameError(f"Parameter {parameter} does not exist or the register containing it is not implemented.") - return reg_addr - - - - def read_param(self, parameter, ch=0, reg_name=None, chip='A'): - if reg_name: - reg_addr = self._regs[reg_name].addr - else: - reg_addr = self._find_reg_containing_param(parameter, ch) - diff --git a/siphractl/ComptonCam_DMA_To_RAW_File/Makefile b/siphractrl/ComptonCam_DMA_To_RAW_File/Makefile similarity index 100% rename from siphractl/ComptonCam_DMA_To_RAW_File/Makefile rename to siphractrl/ComptonCam_DMA_To_RAW_File/Makefile diff --git a/siphractl/ComptonCam_DMA_To_RAW_File/src/dma_to_raw_file.c b/siphractrl/ComptonCam_DMA_To_RAW_File/src/dma_to_raw_file.c similarity index 100% rename from siphractl/ComptonCam_DMA_To_RAW_File/src/dma_to_raw_file.c rename to siphractrl/ComptonCam_DMA_To_RAW_File/src/dma_to_raw_file.c diff --git a/siphractl/MK_WriteConfiguration_251210.py b/siphractrl/MK_WriteConfiguration_251210.py similarity index 100% rename from siphractl/MK_WriteConfiguration_251210.py rename to siphractrl/MK_WriteConfiguration_251210.py diff --git a/siphractl/acquire.py b/siphractrl/acquire.py similarity index 91% rename from siphractl/acquire.py rename to siphractrl/acquire.py index fc40d14..8eeca60 100644 --- a/siphractl/acquire.py +++ b/siphractrl/acquire.py @@ -30,13 +30,14 @@ def parse_args(): help="A single string containing \'Background\' or any radioactive source present") parser.add_argument("--source-description", type=str, default='[NOT SPECIFIED]', help="A single string.") - parser.add_argument("--notes", type=str, default='[NONE]') + parser.add_argument("--notes", type=str, default='[NONE]', + help="A single string.") parser.add_argument("--siphra-config-file", type=str, default='D2a/Ongoing.txt', help="By default, it takes the \'Ongoing.txt\' file in the \'D2a\' directory and writes its contents to the metadata file. If any other configuration is used, the path to the text file containing it must be specified here.") parser.add_argument("-s", "--size", type=int, default=4095, - help="([DO NOT CHANGE UNLESS CHANGING dma_to_raw_file SETTINGS]. Block size. Default is 4095.") + help="[DO NOT CHANGE UNLESS CHANGING dma_to_raw_file SETTINGS]. Block size. Default is 4095.") parser.add_argument("--device", default="/dev/D2A_DMA", - help = "([DO NOT CHANGE UNLESS CHANGING dma_to_raw_file SETTINGS].") + help = "[DO NOT CHANGE UNLESS CHANGING dma_to_raw_file SETTINGS].") return parser.parse_args() diff --git a/siphractl/d2a_lib.py b/siphractrl/d2a_lib.py similarity index 99% rename from siphractl/d2a_lib.py rename to siphractrl/d2a_lib.py index 1dee38e..2dce62d 100644 --- a/siphractl/d2a_lib.py +++ b/siphractrl/d2a_lib.py @@ -65,7 +65,7 @@ def compareUpTo(a, b, len): class D2a: ''' - Contains all fucntions required to interface with hardware devices through OS. + Contains all functions required to interface with hardware devices through OS. ''' def __init__(self, uio='/dev/uio0', spi='/dev/spi'): try: diff --git a/siphractl/regs_bit_structure.py b/siphractrl/regs_bit_structure.py similarity index 68% rename from siphractl/regs_bit_structure.py rename to siphractrl/regs_bit_structure.py index bdfe3a8..8fa76e1 100644 --- a/siphractl/regs_bit_structure.py +++ b/siphractrl/regs_bit_structure.py @@ -16,7 +16,7 @@ ) # Register 0x10: Summing channel configuration registers -SUMM_CHANNEL = BitStruct( +SUM_CHANNEL = BitStruct( Padding(18), 'cal_select_channel' / Flag, 'qc_threshold' / BitsInteger(8), @@ -88,4 +88,46 @@ 'pu_mb' / Flag, ) -# \ No newline at end of file +# Register 0x16 +CAL_CTRL = BitStruct( + Padding(26), + 'cal_cv_range_low' / Flag, + 'cal_cv_range_mid' / Flag, + 'cal_cv_range_high' / Flag, + 'cal_trigger_select' / Flag, + 'cal_mode' / BitsInteger(2) +) + +# Register 0x17 +ch_names = [f'CH_{n}_rofl' for n in range(16, 0, -1)] +READOUT_FIXED_LIST = BitStruct( + Padding(13), + 'ADC_IN_rofl' / Flag, + 'SUM_CH_rofl' / Flag, + *[name / Flag for name in ch_names], + 'ADC_PT100_SENSE_rofl' / Flag, +) + +# Register 0x18 +READOUT_MODE = BitStruct( + Padding(17), + 'int_hold_tune' / BitsInteger(5), + 'int_hold_delay' / BitsInteger(4), + 'int_hold_source' / BitsInteger(2), + 'readout_list_mode' / BitsInteger(1), + 'readout_en_spi_forced_start' / Flag, + 'readout_en_ext_hold_start' / Flag, + 'readout_en_int_hold_start' / Flag, +) + +reg_str_lst = [*[CHANNEL]*16, + SUM_CHANNEL, + CHANNEL_CONFIG, + CHANNEL_CONTROL, + ADC_CONFIG, + CAL_DAC, + PD_MODULES, + CAL_CTRL, + READOUT_FIXED_LIST, + READOUT_MODE, + ] \ No newline at end of file diff --git a/siphractrl/siphra_controller.py b/siphractrl/siphra_controller.py new file mode 100644 index 0000000..1ec194d --- /dev/null +++ b/siphractrl/siphra_controller.py @@ -0,0 +1,149 @@ +from unittest import case + +from .d2a_lib import * +from .regs_bit_structure import * +from collections import namedtuple + +CH_ADDRS = [] +SIPHRA_RETURN_SIZE = 15 # bytes. Returned by SIPHRA when D2a.readSIPHRA is called +# TODO: Check the size of the packet returned by D2a.readSIPHRA +REG_SIZE = 4 # bytes. Size of the bytearray passed to D2a.writeSIPHRA + +# Register map classes +RegField = namedtuple('RegField', ['name', 'size']) + +class SIPHRARegister: + def __init__(self, addr, structure): + self.addr = addr + self.structure = structure + self.fields = self.structure.subcon.subcons # Includes padding + self.field_names = [field.name for field in self.fields if field.name] + self.size = sum([field.sizeof() for field in self.fields if field.name]) # Number of significant bits in the register + + def __getitem__(self, idx): + '''Returns the name and bit-size of a given field. Index 0 is the padding.''' + if idx >= (n_fields:=len(self.fields)) or idx < 0: + raise IndexError(f"Index {idx} is out of range for register with {n_fields - 1} fields + padding") + if idx == 0: + return RegField(name='Padding', size=self.fields[0].sizeof()) + return RegField(name=self.fields[idx].name, size=self.fields[idx].sizeof()) + + def __len__(self): + '''Number of fields in this register (excluding padding).''' + return len(self.fields) - 1 + + def __contains__(self, field_name): + '''True if the register has a field named ``field_name``.''' + return field_name in self.field_names + + def parse(self, value): + return self.structure.parse(value) + + def set_param(self, param_name, value, current_content): + reg_content = parse(current_content) + reg_content[param_name] = value + return self.structure.build(reg_content) + +class SIPHRA: + def __init__(self): + self._d2a = D2a() + self._regs = {} + + for i in range(1, 17): + self._regs[f"ch{i}"] = SIPHRARegister(i, CHANNEL) + + self._regs["ch17"] = SIPHRARegister(0x10, SUM_CHANNEL) + self._regs["ch_config"] = SIPHRARegister(0x11, CHANNEL_CONFIG) + self._regs["ch_control"] = SIPHRARegister(0x12, CHANNEL_CONTROL) + self._regs["adc_config"] = SIPHRARegister(0x13, ADC_CONFIG) + self._regs["cal_dac"] = SIPHRARegister(0x14, CAL_DAC) + self._regs["pd_modules"] = SIPHRARegister(0x15, PD_MODULES) + self._regs["cal_ctrl"] = SIPHRARegister(0x16, CAL_CTRL) + self._regs["readout_list"] = SIPHRARegister(0x17, READOUT_FIXED_LIST) + self._regs["readout_mode"] = SIPHRARegister(0x18, READOUT_MODE) + + def _resolve_reg_id(self, reg_id: str | int) -> tuple[str, int]: + """Resolves a register name or address to a (name, addr) tuple.""" + match reg_id: + case str(): + if reg_id in self._regs: + return reg_id, self._regs[reg_id].addr + case int(): + for name, reg in self._regs.items(): + if reg.addr == reg_id: + return name, reg_id + case _: + raise TypeError(f"reg_id must be str or int, got {type(reg_id).__name__}") + raise ValueError(f"{reg_id!r} is not a valid register id") + + def _find_reg_containing_param(self, parameter, ch: int=0): + ''' + Returns the address of the register containing the given parameter. + :param parameter: Name of the desired parameter. + :param ch: If the parameter is in one of the 16+1 channel addresses, this argument is used for disambiguation. ``ch`` is a number between 1 and 16 if the parameter belongs to one of the ``ctrl_ch`` registers; it is 17 if the parameter belongs to the summing channel control register, and is defaulted to 0 if the parameter belongs to any other register. + :return: The address of the register containing the given parameter. + ''' + addr = None + name = None + if not 0 <= ch <= 17: + raise ValueError(f"Channel {ch} is out of range. Use channels 1-16 and 17 for summing channel") + if ch == 0: # Not a channel register + for n, reg in self._regs.items(): + if parameter in reg: + name = n + addr = reg.addr + else: # Channel register + # Verify that the given field exists + reg = self._regs[f"ch{ch}"] + if parameter in reg: + addr = reg.addr + name = f"ch{ch}" + if not addr: + raise NameError(f"Parameter {parameter} does not exist or the register containing it is not implemented.") + return name, addr + + def get_reg_value(self, reg_id, chip='A'): + _, addr = self._resolve_reg_id(reg_id) + return self._d2a.readSIPHRA(addr, chip) + + def read_register(self, reg_id, chip='A'): + name, addr = self._resolve_reg_id(reg_id) + reg_value = self._d2a.readSIPHRA(addr, chip) + return self._regs[name].parse(reg_value) + + def write_register(self, reg_id, value, chip='A'): + _, addr = self._resolve_reg_id(reg_id) + try: + self._d2a.writeSIPHRAwithCheck(addr, chip, value) + except Exception as e: + print(e) + return -1 + return 0 + + def read_param(self, parameter, ch=0, reg_id=None, chip='A'): + '''``reg_id`` can be the register name or its address''' + if reg_id: + name, addr = self._resolve_reg_id(reg_id) + if not parameter in _regs[name]: + raise NameError(f"Parameter {parameter} does not exist in register {reg_id}.") + else: + name, addr = self._find_reg_containing_param(parameter, ch) + + return self.read_register(name, chip=chip)[parameter] + + def write_param(self, parameter, value, ch=0, reg_id=None, chip='A'): + if reg_id: + name, addr = self._resolve_reg_id(reg_id) + if not parameter in _regs[name]: + raise NameError(f"Parameter {parameter} does not exist in register {reg_id}.") + else: + name, addr = self._find_reg_containing_param(parameter, ch) + + reg = self._regs[name] + old_value = self.get_reg_value(reg_id, chip=chip) + new_value = reg.set_param(parameter, value, old_value) + self.write_register(name, new_value, chip=chip) + + + +