{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Stream Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook provides an analysis of the software triggered data of a detector module. This includes the calculation of rate and stability cuts, energy spectra and light yield plots, and efficiencies." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2021-07-04T23:31:21.371501Z", "start_time": "2021-07-04T23:31:21.367783Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Let´s start!\n" ] } ], "source": [ "print('Let´s start!')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook builds upon the hardware data analysis notebook and the trigger script. We do cuts and the energy calibration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we import packages." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import cait as ai\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import h5py\n", "import pickle\n", "from scipy.stats import norm, expon\n", "%config InlineBackend.figure_formats = ['svg'] # we need this for a suitable resolution of the plots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we define a set of constants and paths." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "RUN = ... # put an string for the number of the experiments run, e.g. '34'\n", "MODULE = ... # put a name for the detector, e.g. 'DetA'\n", "PATH_PROC_DATA = ... # path to where you want to store the HDF5 files\n", "FILE_NMBRS = [] # a list of string, the file number you want to analyse, e.g. ['001', '002', '003']\n", "RDT_CHANNELS = [] # a list of strings of the channels, e.g. [0, 1] (written in PAR file - attention, the PAR file counts from 1, Cait from 0)\n", "RECORD_LENGTH = 16384 # the number of samples within a record window (read in PAR file)\n", "SAMPLE_FREQUENCY = 25000 # the sample frequency of the measurement (read in PAR file)\n", "DOWN_SEF = 4 # the downsample rate for the standard event fit\n", "DOWN_BLF = 16 # the downsample rate for the baseline fit\n", "PROCESSES = 8 # the number of processes for parallelization\n", "PCA_COMPONENTS = 2 # the number of pca components to calculate\n", "SKIP_FNMR = [] # in case the loop crashed at some point and you want to start from a specific file number, write here the numbers to ignore, e.g. ['001', '002']\n", "\n", "# typically you do not need to change the values below\n", "\n", "FNAME_EFFICIENCY = 'efficiency_{:03d}'.format(len(FILE_NMBRS) - 1)\n", "FNAME_TRAINING = 'training_001'\n", "FNAME_STREAM = 'stream_{:03d}'.format(len(FILE_NMBRS) - 1)\n", "H5_CHANNELS = list(range(len(RDT_CHANNELS)))\n", "SEF_APP = '_down{}'.format(DOWN_SEF) if DOWN_SEF > 1 else ''" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need some values that we already calculated from the hardware data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "DETECTOR_MASS = ... # the detector mass in kg\n", "print('Detector mass in kg: ', DETECTOR_MASS)\n", "BL_RESOLUTION_OF = [] # list of the baseline resolutions, calculated with the superposition method, in mV\n", "THRESHOLDS = [(6.5 * r) * 1e-3 for r in BL_RESOLUTION_OF]\n", "print('OF resolution in V: ', BL_RESOLUTION_OF)\n", "print('OF thresholds in V: ', THRESHOLDS)\n", "FIT_BL_SIGMAS = [] # list of the baseline resolutions, calculated with the noise fit model, in V\n", "FIT_THRESHOLDS = [] # list of the trigger thresholds, calculated with the noise fit model, in V\n", "print('Fit thresholds in V: ', FIT_THRESHOLDS)\n", "TRUNCATION_LEVELS = [] # list of the truncation levels\n", "MAXIMAL_EVENT_HEIGHTS = [] # list of the maximal event heights to be included in the energy calibration\n", "PATH_PULSER_MODEL = ... # put an string of the path where you want to store the pulser model\n", "\n", "PEAK_ENERGY = 5.89 * 8/9 + 6.49 * 1/9 # in case of an iron source - change this for different source" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we assemble calculated values from further down in the notebook. Fill them up while you go." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "PEC_FACTOR = [...] # list of the values source_peak_energy/source_peak_pulse_height\n", "print('Energy resolution (eV): ', [b*p for b,p in zip(BL_RESOLUTION_OF, PEC_FACTOR)])\n", "print('Treshold (eV): ', [1000*b*p for b,p in zip(THRESHOLDS, PEC_FACTOR)])\n", "print('Treshold (eV) 1 Noise Trigger Fit: ', [1000*b*p for b,p in zip(FIT_THRESHOLDS, PEC_FACTOR)])\n", "print('Treshold from Trigger Efficiency (eV): ', 1000*np.array([...])) # The fitted trigger efficiency threshold\n", "CPE_FACTOR = [...]# list of the values source_peak_energy/source_peak_equivalent_tpa_value\n", "LIVE_TIME = ... # the total live time of the detector, calculated at the end of the file\n", "EXPOSURE = LIVE_TIME/24 * DETECTOR_MASS\n", "print('Live Time in h: ', LIVE_TIME)\n", "print('Exposure in kg * days: ', EXPOSURE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a DataHandler instance to do the processing." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dh_stream = ai.DataHandler(run=RUN,\n", " module=MODULE,\n", " channels=RDT_CHANNELS)\n", "\n", "dh_stream.set_filepath(path_h5=PATH_PROC_DATA,\n", " fname=FNAME_STREAM,\n", " appendix=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## View Events" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets have a look at the stream triggered events to check, if the triggering worked properly." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ei = ai.EventInterface(module=MODULE,\n", " run=RUN,\n", " nmbr_channels=len(H5_CHANNELS))\n", "ei.load_h5(path=PATH_PROC_DATA, \n", " fname=FNAME_STREAM, \n", " channels=RDT_CHANNELS, \n", " appendix=False, \n", " which_to_label=['events'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to include a labels file." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ei.create_labels_csv(path=PATH_PROC_DATA)\n", "# ei.load_labels_csv(path=PATH_PROC_DATA)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also load the optimal filter transfer function to view the events the way they were triggered." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ei.load_of(group_name_appendix='_tp')\n", "# ei.load_sev_par(name_appendix='_down4', group_name_appendix='_tp')\n", "# ei.load_of()\n", "# ei.load_sev_par(name_appendix='_down4')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, lets see the events." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ei.start(start_from_idx=0, print_label_list=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cuts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To clean the data from artifacts, we want to apply some quality cuts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case we have not included the SEV, NPS and OF during the trigger loop, and have not calculated the corresponding features, we can do this here." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Only run this if it wasnt run during the trigger loop!\n", "\n", "dh_hw = ai.DataHandler(run=RUN,\n", " module=MODULE,\n", " channels=RDT_CHANNELS)\n", " \n", "dh_hw.set_filepath(path_h5=PATH_PROC_DATA,\n", " fname='hw_{:03d}'.format(len(FILE_NMBRS)-1),\n", " appendix=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Only run this if it wasnt run during the trigger loop!\n", "\n", "dh_stream.include_sev(sev=dh_hw.get('stdevent','event'), \n", " fitpar=dh_hw.get('stdevent','fitpar'), \n", " mainpar=dh_hw.get('stdevent','mainpar'))\n", "\n", "dh_stream.include_nps(nps=dh_hw.get('noise','nps'))\n", "\n", "dh_stream.include_of(of_real=dh_hw.get('optimumfilter','optimumfilter_real'), \n", " of_imag=dh_hw.get('optimumfilter','optimumfilter_imag'))\n", "\n", "dh_stream.include_sev(sev=dh_hw.get('stdevent_tp','event'), \n", " fitpar=dh_hw.get('stdevent_tp','fitpar'), \n", " mainpar=dh_hw.get('stdevent_tp','mainpar'),\n", " group_name_appendix='_tp')\n", "\n", "dh_stream.include_of(of_real=dh_hw.get('optimumfilter_tp','optimumfilter_real'), \n", " of_imag=dh_hw.get('optimumfilter_tp','optimumfilter_imag'),\n", " group_name_appendix='_tp')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Only run this if it wasnt run during the trigger loop!\n", "\n", "dh_stream.apply_of(first_channel_dominant=True)\n", "dh_stream.apply_of(type='testpulses', name_appendix_group='_tp')\n", "dh_stream.apply_sev_fit(down=DOWN_SEF, name_appendix='_down{}'.format(DOWN_SEF), processes=PROCESSES,\n", " truncation_level=TRUNCATION_LEVELS, verb=True, first_channel_dominant=True)\n", "dh_stream.apply_sev_fit(type='testpulses', group_name_appendix='_tp', \n", " down=DOWN_SEF, name_appendix='_down{}'.format(DOWN_SEF), processes=PROCESSES,\n", " truncation_level=TRUNCATION_LEVELS, verb=True)\n", "\n", "# do your own cuts!\n", "clean_events = ai.cuts.LogicalCut(initial_condition=np.abs(dh_stream.get('events', 'mainpar')[0,:,8]) < 2e-6)\n", "clean_events.add_condition(np.abs(dh_stream.get('events', 'mainpar')[1,:,8]) < 2e-6)\n", "clean_events.add_condition(dh_stream.get('events', 'mainpar')[0,:,0] < 1) \n", "clean_events.add_condition(dh_stream.get('events', 'mainpar')[1,:,0] < 1.5) \n", "clean_events.add_condition(dh_stream.get('events', 'mainpar')[0,:,3] < 4500) \n", "clean_events.add_condition(dh_stream.get('events', 'mainpar')[0,:,3] > 3900)\n", "clean_events.add_condition(dh_stream.get('events', 'mainpar')[1,:,3] < 4500) \n", "clean_events.add_condition(dh_stream.get('events', 'mainpar')[1,:,3] > 3900)\n", "\n", "for c in H5_CHANNELS:\n", " dh_stream.apply_logical_cut(cut_flag=clean_events.get_flag(),\n", " naming='clean_events',\n", " channel=c,\n", " type='events',\n", " delete_old=False)\n", "\n", "dh_stream.apply_pca(nmbr_components=PCA_COMPONENTS, \n", " down=DOWN_SEF, \n", " fit_idx=clean_events.get_idx())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One cut we want to apply is a rate cut." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dh_stream.calc_rate_cut(min=0, max=30)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the histogram of the rate." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nmbr_bins = 500\n", "total_hours = 1000\n", "\n", "good_intervals = dh_stream.get('metainfo', 'rate_stable')\n", "\n", "dh_stream.show_values(group='events', key='hours', bins=nmbr_bins, \n", " xlabel='Time (h)', ylabel=f'Counts/({total_hours/nmbr_bins*60:.1f} min)', title='Event Rate',\n", " show=False)\n", "for iv in good_intervals:\n", " plt.axvspan(xmin=iv[0], xmax=iv[1], color='green', alpha=0.3, zorder=80)\n", "\n", "# plt.xlim((400,450))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do the stability cut." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cp_ph = dh_stream.get('controlpulses', 'pulse_height')\n", "\n", "bounds = []\n", "\n", "for c in H5_CHANNELS:\n", " median_cp = np.median(cp_ph[c])\n", " print('Median pulse height Channel {}: {} V'.format(c, median_cp))\n", " bounds.append((median_cp*0.8, median_cp*1.2))\n", " dh_stream.calc_controlpulse_stability(channel=c, significance=3, lb=bounds[c][0], ub=bounds[c][1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the control pulse heights." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cp_h = dh_stream.get('controlpulses', 'hours')\n", "\n", "y_lims = [None, None]\n", "\n", "for c in H5_CHANNELS:\n", " instable_iv = dh_stream.get('metainfo', f'controlpulse_instable_ch{c}')\n", " \n", " plt.close()\n", " xlims, ylims, I = ai.styles.scatter_img(x_data=cp_h, y_data=cp_ph[c], ylims=y_lims[c], height=80, width=200, alpha=0.3)\n", " ax_extent = list(xlims)+list(ylims)\n", " ai.styles.use_cait_style()\n", " plt.imshow(I,\n", " vmin=0,\n", " vmax=1, \n", " cmap=plt.get_cmap('viridis'),\n", " interpolation='lanczos',\n", " aspect='auto',\n", " extent=ax_extent,\n", " )\n", " for iv in instable_iv:\n", " plt.axvspan(iv[0], iv[1], color='red', alpha=0.2, zorder=70)\n", " plt.grid(alpha=0.2)\n", " plt.title('Control Pulse Stability Channel {}'.format(c))\n", " plt.xlabel('Time (h)')\n", " plt.ylabel('Pulse Height (V)')\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do the quality cuts further down." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Energy Calibration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can double check the CPE factors with the iron calibration peak, which is done below. Afterwards, we do the energy fine tuning with the test pulses.\n", "\n", "For as estimator for the pulse height, we always take the optimum filter." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_ranges = [(0, 1.6), (0, 0.3)]\n", "y_ranges = [(0, 200), (0, 1000)]\n", "\n", "for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):\n", " dh_stream.show_values(group='events', key='of_ph', bins=250, idx0=c, idx2=None, range=xr, yran=yr,\n", " xlabel='Pulse Height (V)', ylabel='Counts', title='Spectrum PH Channel {}'.format(c), show=False)\n", " # plt.axvline(x=calibration_peak[c], color='black', linestyle='dotted', linewidth=2.5, zorder=70) # run the cell below first!\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ph = dh_stream.get('events','of_ph')\n", "iv = [(0.3, 0.45), (0.008, 0.03)]\n", "calibration_peak = []\n", "\n", "for c in H5_CHANNELS:\n", " peak_events = ai.cuts.LogicalCut(initial_condition=ph[c]>iv[c][0])\n", " peak_events.add_condition(ph[c] 1.6)\n", "tpa_peak_cut.add_condition(tpa_equivalents[0] < 2.4)\n", "tpa_peak_cut.add_condition(tpa_equivalents[1] < 0.9)\n", "\n", "peak_pos = []\n", "cpe_factor = []\n", "\n", "for c in H5_CHANNELS:\n", " peak_pos.append(np.median(tpa_equivalents[c, tpa_peak_cut.get_flag()]))\n", " print('Peak position Channel {}: {} V_TPA'.format(c, peak_pos[-1]))\n", " cpe_factor.append(iron_peak/peak_pos[-1])\n", " print('CPE factor Channel {}: {} keV/V_TPA'.format(c, cpe_factor[-1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can overwrite the calculated recoil energies with the new ones, with a precise CPE factor." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in H5_CHANNELS:\n", " dh_stream.include_values(values=dh_stream.get('events', 'tpa_equivalent')[c]*cpe_factor[c], naming='recoil_energy', channel=c, delete_old=True)\n", " dh_stream.include_values(values=dh_stream.get('events', 'tpa_equivalent_sigma')[c]*cpe_factor[c], naming='recoil_energy_sigma', channel=c, delete_old=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Quality Cuts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We apply some logical cuts. Before that we have a look at the main parameter distribution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ph = dh_stream.get('events', 'pulse_height')\n", "ydata = dh_stream.get('events', 'rise_time')\n", "# ydata = dh_stream.get('events', 'decay_time')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xlims = [(0,1.5), (0,0.8)]\n", "ylims = [(-1,5), (-1,10)]\n", "# ylims = [(-1,40), (-1,25)]\n", "\n", "for c, x, y in zip(H5_CHANNELS, xlims, ylims):\n", " plt.close()\n", " xlims, ylims, I = ai.styles.scatter_img(x_data=ph[c], y_data=ydata[c], ylims=y, xlims=x, height=100, width=200, alpha=0.5)\n", " ax_extent = list(xlims)+list(ylims)\n", " ai.styles.use_cait_style()\n", " plt.imshow(I,\n", " vmin=0,\n", " vmax=1, \n", " cmap=plt.get_cmap('viridis'),\n", " interpolation='lanczos',\n", " aspect='auto',\n", " extent=ax_extent,\n", " )\n", " #plt.axhline(y=0.25, color='red', linewidth=2) # three event types\n", " #plt.axhline(y=0.3, color='red', linewidth=2) # three event types\n", " #plt.axvline(x=2.5, color='green', linewidth=2) # saturation\n", " #plt.axvline(x=0.5, color='green', linewidth=2) # noise\n", " plt.grid(alpha=0.2)\n", " plt.title(' Channel {}'.format(c))\n", " plt.xlabel('Pulse Height (V)')\n", " plt.ylabel('Rise Time (ms)')\n", " # plt.ylabel('Decay Time (ms)')\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in H5_CHANNELS:\n", " dh_stream.show_values(group='events', key='mainpar', bins=200, idx0=c, idx2=3, #range=(-5e-6, 5e-6),\n", " xlabel='Sample Index', ylabel='Counts', title='Maximum Position Channel {}'.format(c))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ph_cut = ai.cuts.LogicalCut(initial_condition=np.abs(dh_stream.get('events', 'slope')[0,:]) < 0.2) # slope\n", "ph_cut.add_condition(np.abs(dh_stream.get('events', 'slope')[1,:]) < 0.2) # slope\n", "ph_cut.add_condition(dh_stream.get('events', 'pulse_height')[0,:,0] < 1.6) # ph\n", "ph_cut.add_condition(dh_stream.get('events', 'pulse_height')[1,:,0] < 0.3) # ph\n", "ph_cut.add_condition(dh_stream.get('events', 'onset')[0,:,3] < 50) # max pos phonon\n", "ph_cut.add_condition(dh_stream.get('events', 'onset')[0,:,3] > 10) # max pos phonon\n", "ph_cut.add_condition(dh_stream.get('events', 'pulse_height')[0,:] < MAXIMAL_EVENT_HEIGHTS[0]) # saturation\n", "ph_cut.add_condition(dh_stream.get('events', 'of_ph')[0,:] > FIT_THRESHOLDS[0]) # threshold\n", "\n", "print('Surviving Event Ratio: ', ph_cut.counts(), np.sum(ph_cut.get_flag())/len(ph_cut.get_flag()))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in H5_CHANNELS:\n", " dh_stream.apply_logical_cut(cut_flag=ph_cut.get_flag(),\n", " naming='ph_cut',\n", " channel=c,\n", " type='events',\n", " delete_old=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cut Efficiency" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We determine the cut efficiency, that we need in the high level analysis to re-weight the event counts. This is the point at which you should run the script for the cut efficiency simulation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dh_eff = ai.DataHandler(run=RUN,\n", " module=MODULE,\n", " channels=RDT_CHANNELS)\n", "\n", "dh_eff.set_filepath(path_h5=PATH_PROC_DATA,\n", " fname=FNAME_EFFICIENCY,\n", " appendix=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dh_eff.calc_rate_cut(intervals=dh_stream.get('metainfo', 'rate_stable'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nmbr_bins = 500\n", "total_hours = 1000\n", "\n", "good_intervals = dh_eff.get('metainfo', 'rate_stable')\n", "\n", "dh_eff.show_values(group='events', key='hours', bins=nmbr_bins, \n", " xlabel='Time (h)', ylabel=f'Counts/({total_hours/nmbr_bins*60:.1f} min)', title='Event Rate',\n", " show=False)\n", "for iv in good_intervals:\n", " plt.axvspan(xmin=iv[0], xmax=iv[1], color='green', alpha=0.3, zorder=80)\n", "\n", "# plt.xlim((400,450))\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in H5_CHANNELS:\n", " dh_eff.calc_controlpulse_stability(channel=c, instable_iv=dh_stream.get('metainfo', f'controlpulse_instable_ch{c}'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "surviving = ai.cuts.LogicalCut(initial_condition=np.abs(dh_eff.get('events', 'slope')[0]) < 0.1)\n", "surviving.add_condition(np.abs(dh_eff.get('events', 'onset')[0]) < 20) # max pos phonon\n", "surviving.add_condition(np.abs(dh_eff.get('events', 'rise_time')[0]) < 2.5) # rise time\n", "surviving.add_condition(np.abs(dh_eff.get('events', 'rise_time')[0]) > 0.5) # rise time\n", "surviving.add_condition(dh_eff.get('events', 'rate_cut')) # rate cut\n", "surviving.add_condition(dh_eff.get('events', 'controlpulse_stability')[0]) # stability cut\n", "surviving.add_condition(dh_eff.get('events', 'pulse_height')[0,:] < MAXIMAL_EVENT_HEIGHTS[0]) # saturation\n", "surviving.add_condition(dh_eff.get('events', 'of_ph')[0,:] > FIT_THRESHOLDS[0]) # threshold\n", "\n", "print('Surviving Event Ratio: ', surviving.counts(), np.sum(surviving.get_flag())/len(surviving.get_flag()))\n", "\n", "for c in H5_CHANNELS:\n", " dh_eff.apply_logical_cut(cut_flag=surviving.get_flag(),\n", " naming='surviving',\n", " channel=c,\n", " type='events',\n", " delete_old=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "surviving = dh_eff.get('events', 'surviving')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "with open(PATH_PULSER_MODEL, 'rb') as f:\n", " pm = pickle.load(f)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# run this only if you havent during the efficiency loop\n", "\n", "dh_eff.calc_calibration(starts_saturation=MAXIMAL_EVENT_HEIGHTS,\n", " cpe_factor=CPE_FACTOR,\n", " poly_order=3,\n", " plot=False,\n", " method='of',\n", " pulser_models=pm,\n", " name_appendix_energy='_reconstructed',\n", " )\n", "\n", "dh_eff.calc_calibration(starts_saturation=MAXIMAL_EVENT_HEIGHTS,\n", " cpe_factor=CPE_FACTOR,\n", " poly_order=3,\n", " plot=False,\n", " method='true_ph',\n", " pulser_models=pm,\n", " name_appendix_energy='_true',\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in [H5_CHANNELS[0]]:\n", " efficiency, counts, bins = dh_eff.show_efficiency(channel=c, \n", " cut_flag=surviving[c], \n", " which_quantity='recoil_energy_true', \n", " bins=200,\n", " title='Cut Efficiency',\n", " xlabel='True Energy (keV)',\n", " ylabel='Survival Probability',\n", " range=[1e-3, 2.6],\n", " show=True,\n", " xscale='linear')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... and only for the low energy region:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in [H5_CHANNELS[0]]:\n", " efficiency, counts, bins = dh_eff.show_efficiency(channel=c, \n", " cut_flag=surviving[c], \n", " which_quantity='recoil_energy_true', \n", " bins=200,\n", " title='Cut Efficiency',\n", " xlabel='True Energy (keV)',\n", " ylabel='Survival Probability',\n", " show=True,\n", " xscale='linear',\n", " range=(1e-3, 0.5),)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit trigger efficiency." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for c in H5_CHANNELS:\n", " ai.fit.fit_trigger_efficiency(binned_energies=bins, survived_fraction=efficiency,\n", " a1_0=0.09, a0_0=1, a2_0=0.01, plot=True, \n", " title='Trigger Efficiency Channel {}'.format(c), \n", " xlim=(0,0.2),\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write the time dependent cut efficiency to a text file for limit calculations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.savetxt(PATH_PROC_DATA + 'xy_files/' + MODULE + '_cuteff.xy',\n", " np.column_stack([dh_eff.get('events', 'hours')*60*60*1e6, \n", " dh_eff.get('events', 'recoil_energy_true')[0],\n", " dh_eff.get('events', 'recoil_energy_reconstructed')[0],\n", " dh_eff.get('events', 'surviving')[0].astype(int), \n", " (dh_eff.get('events', 'of_ph')[0] > FIT_THRESHOLDS[0]).astype(int)]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Spectrum and Light Yield Plot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We plot an energy spectrum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_ranges = [(0, 8), (0, 80)]\n", "y_ranges = [(0, 250), (0, 150)]\n", "bins=200\n", "\n", "for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):\n", " dh_stream.show_values(group='events', key='recoil_energy', bins=bins, idx0=c, cut_flag=ph_cut.get_flag(), range=xr, yran=yr,\n", " xlabel='Recoil Energy (keV)', ylabel='Counts / {} eV'.format(int((xr[1]-xr[0])/bins*1000)), \n", " save_path='plots/{}/{}_Spectrum_LE_Channel{}.pdf'.format(MODULE, MODULE, c))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... and only the low-energy part." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_ranges = [(0, 0.5), (0, 80)]\n", "y_ranges = [(0, 50), (0, 150)]\n", "bins=150\n", "\n", "for c, xr, yr in zip([0], x_ranges, y_ranges):\n", " dh_stream.show_values(group='events', key='recoil_energy', bins=bins, idx0=c, cut_flag=ph_cut.get_flag(), range=xr, yran=yr,\n", " xlabel='Recoil Energy (keV)', ylabel='Counts / {} eV'.format(int((xr[1]-xr[0])/bins*1000)), \n", " save_path='plots/{}/{}_Spectrum_LE_Channel{}.pdf'.format(MODULE, MODULE, c))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sometimes we see the peaks better, when we perform density estimation with a Gaussian kernal and the individual recoil energy uncertainties." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "recoil_energy = dh_stream.get('events', 'recoil_energy')\n", "recoil_energy = recoil_energy[:, ph_cut.get_flag()]\n", "recoil_energy_sigma = dh_stream.get('events', 'recoil_energy_sigma')\n", "recoil_energy_sigma = recoil_energy_sigma[:, ph_cut.get_flag()]\n", "\n", "x_ranges = [(0, 10), (0, 100)]\n", "y_ranges = [(0, 35000), (0, 3500)]\n", "\n", "for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):\n", "\n", " # density estimation\n", " x_d = np.linspace(xr[0], xr[1], 1000)\n", " density = sum(norm(loc=m, scale=s).pdf(x_d) for m,s in zip(recoil_energy[c], recoil_energy_sigma[c]))/EXPOSURE\n", "\n", " plt.close()\n", " ai.styles.use_cait_style(dpi=150)\n", " plt.fill_between(x_d, density, alpha=0.8, zorder=25, label='Density Estimate')\n", " plt.plot(recoil_energy[c], np.full_like(recoil_energy[c], -0.1), '|k', markeredgewidth=1, label='Counts')\n", " ai.styles.make_grid()\n", " plt.xlim(xr)\n", " plt.ylim(yr)\n", " plt.xlabel('Recoil Energy (keV)')\n", " plt.ylabel('Density (1/(keV kg days))')\n", " plt.legend()\n", " plt.savefig('plots/{}/{}_Density_Channel{}.pdf'.format(MODULE, MODULE, c))\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We single out the iron lines." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "recoil_energy = dh_stream.get('events', 'recoil_energy')\n", "recoil_energy_sigma = dh_stream.get('events', 'recoil_energy_sigma')\n", "\n", "iron_lower = ai.cuts.LogicalCut(recoil_energy[0] > 5)\n", "iron_lower.add_condition(recoil_energy[0] < 6.2)\n", "iron_lower.add_condition(ph_cut.get_flag())\n", "\n", "iron_upper = ai.cuts.LogicalCut(recoil_energy[0] > 6.2)\n", "iron_upper.add_condition(recoil_energy[0] < 6.7)\n", "iron_upper.add_condition(ph_cut.get_flag())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rec_lower = recoil_energy[:, iron_lower.get_flag()]\n", "rec_sigma_lower = recoil_energy_sigma[:, iron_lower.get_flag()]\n", "rec_upper = recoil_energy[:, iron_upper.get_flag()]\n", "rec_sigma_upper = recoil_energy_sigma[:, iron_upper.get_flag()]\n", "\n", "x_ranges = [(5, 7), (0, 20)]\n", "y_ranges = [(0, 30000), (0, 2000)]\n", "\n", "for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):\n", "\n", " print('Mean lower line, channel {}: {}'.format(c, np.mean(rec_lower[c])))\n", " print('Mean upper line, channel {}: {}'.format(c, np.mean(rec_upper[c])))\n", " \n", " # density estimation lower iron\n", " x_d = np.linspace(xr[0], xr[1], 1000)\n", " density_lower = sum(norm(loc=m, scale=s).pdf(x_d) for m,s in zip(rec_lower[c], rec_sigma_lower[c]))/EXPOSURE\n", "\n", " # density estimation upper iron\n", " density_upper = sum(norm(loc=m, scale=s).pdf(x_d) for m,s in zip(rec_upper[c], rec_sigma_upper[c]))/EXPOSURE\n", "\n", " plt.close()\n", " ai.styles.use_cait_style(dpi=150)\n", " # lower iron\n", " plt.fill_between(x_d, density_lower, alpha=0.8, zorder=25, label='K alpha')\n", " plt.plot(rec_lower[c], np.full_like(rec_lower[c], -0.1), '|k', markeredgewidth=1)\n", " # upper iron\n", " plt.fill_between(x_d, density_upper, alpha=0.8, zorder=25, label='K beta')\n", " plt.plot(rec_upper[c], np.full_like(rec_upper[c], -0.1), '|k', markeredgewidth=1)\n", " ai.styles.make_grid()\n", " plt.xlim(xr)\n", " plt.ylim(yr)\n", " plt.xlabel('Recoil Energy (keV)')\n", " plt.ylabel('Density (1/(keV kg days))')\n", " plt.legend()\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we do a light yield plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "recoil_energy_corr = dh_stream.get('events', 'recoil_energy')\n", "recoil_energy_corr = recoil_energy_corr[:, ph_cut.get_flag()]\n", "\n", "plt.close()\n", "xlims, ylims, I = ai.styles.scatter_img(x_data=recoil_energy_corr[0], \n", " y_data=recoil_energy_corr[1]/recoil_energy_corr[0], \n", " ylims=(-4, 10),\n", " xlims=(0, 10),\n", " height=200, \n", " width=250,\n", " alpha=0.9)\n", "ax_extent = list(xlims)+list(ylims)\n", "ai.styles.use_cait_style()\n", "plt.imshow(I,\n", " vmin=0,\n", " vmax=1, \n", " cmap=plt.get_cmap('viridis'),\n", " interpolation='lanczos',\n", " aspect='auto',\n", " extent=ax_extent,\n", " )\n", "plt.grid(alpha=0.2)\n", "plt.title('Light Yield Plot')\n", "plt.xlabel('Recoil Energy (keV)')\n", "plt.ylabel('Light Yield (a.u.)')\n", "#plt.colorbar()\n", "plt.savefig('plots/{}/{}_LY.pdf'.format(MODULE, MODULE))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is the event rate of low energetic events?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "only_le = ai.cuts.LogicalCut(initial_condition=dh_stream.get('events', 'recoil_energy')[0] < 0.5) # cut at 0.5 keV\n", "only_le.add_condition(ph_cut.get_flag())\n", "print('Total Counts: ', ph_cut.counts())\n", "print('Counts below 0.5 keV: ', only_le.counts())\n", "\n", "dh_stream.show_values(group='events', key='hours', bins=100, \n", " xlabel='Time (h)', ylabel='Counts', title='Event Rate below 0.5 keV', cut_flag=only_le.get_flag())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write the spectrum to a xy file to process it for limit calculations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "recoil_energy_corr = dh_stream.get('events', 'recoil_energy')\n", "recoil_energy_corr = recoil_energy_corr[:, ph_cut.get_flag()]\n", "\n", "data = [recoil_energy_corr[0], \n", " recoil_energy_corr[1]/recoil_energy_corr[0]]\n", "\n", "ai.data.write_xy_file(filepath=PATH_PROC_DATA + 'xy_files/' + MODULE + '_spectrum.xy',\n", " data=data,\n", " title='Run ' + RUN + ' ' + MODULE + ' Spectrum',\n", " axis=['Recoil Energy (keV)', \n", " 'Light Yield (a.u.)'])\n", "\n", "data.append(dh_stream.get('events', 'hours')[ph_cut.get_flag()])\n", "\n", "ai.data.write_xy_file(filepath=PATH_PROC_DATA + 'xy_files/' + MODULE + '_spectrum_timedependent.xy',\n", " data=data,\n", " title='Run ' + RUN + ' ' + MODULE + 'Spectrum Time Dependent',\n", " axis=['Recoil Energy (keV)', \n", " 'Light Yield (a.u.)',\n", " 'Time (h)'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write the times at which the individual files start and stop to an xy file. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dh_stream.generate_startstop()\n", "\n", "ai.data.write_xy_file(filepath=PATH_PROC_DATA + 'xy_files/' + MODULE + '_startstop.xy',\n", " data=dh_stream.get('metainfo', 'startstop_hours').T,\n", " title='Run ' + RUN + ' ' + MODULE + ' Startstop',\n", " axis=['Start Files (hours)', \n", " 'Stop Files (hours)'])\n", "\n", "np.savetxt(PATH_PROC_DATA + 'xy_files/' + MODULE + '_startstop_tstamp.xy',\n", " np.column_stack([dh_stream.get('metainfo', 'startstop_s')[:,0].astype(int), \n", " dh_stream.get('metainfo', 'startstop_mus')[:,0].astype(int), \n", " dh_stream.get('metainfo', 'startstop_s')[:,1].astype(int), \n", " dh_stream.get('metainfo', 'startstop_mus')[:,1].astype(int)]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, what is the exposure?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dh_stream.exposure(detector_mass=DETECTOR_MASS, exclude_instable=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finished." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }