diff --git a/notebooks/example.ipynb b/notebooks/example.ipynb new file mode 100644 index 0000000..b78f36d --- /dev/null +++ b/notebooks/example.ipynb @@ -0,0 +1,383 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8efad73c-6dd9-46bc-bf55-d58dc54dc63c", + "metadata": {}, + "source": [ + "This notebook is a basic example for applying salting method to study ambience interference." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "232ecce7-d687-4eab-ab7e-561f285899cc", + "metadata": {}, + "outputs": [], + "source": [ + "import saltax\n", + "import straxen\n", + "import numpy as np\n", + "from tqdm import tqdm\n", + "import matplotlib.pyplot as plt\n", + "\n", + "straxen.print_versions()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e86c4c17-a60d-46c0-8aea-560558ba18b3", + "metadata": {}, + "outputs": [], + "source": [ + "rn_runs = ['022510','022735','022723','022717']#,'022685'] bad run for match?\n", + "bk_runs = ['026577','028031','031147','027705','025579','030079','029139','028047']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3e8dae8-ec60-4a2b-bb70-4298ee1080d4", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "st_salt = saltax.contexts.sxenonnt(runid=21952,\n", + " saltax_mode='salt',\n", + " output_folder='/project/lgrandi/yuanlq/salt',\n", + " faxconf_version=\"sr0_v4\",\n", + " generator_name='flat',\n", + " recoil=8,\n", + " mode='all')\n", + "st_simu = saltax.contexts.sxenonnt(runid=21952,\n", + " saltax_mode='simu',\n", + " output_folder='/project/lgrandi/yuanlq/salt',\n", + " faxconf_version=\"sr0_v4\",\n", + " generator_name='flat',\n", + " recoil=8,\n", + " mode='all')\n", + "st_data = saltax.contexts.sxenonnt(runid=21952,\n", + " saltax_mode='data',\n", + " output_folder='/project/lgrandi/yuanlq/salt',\n", + " faxconf_version=\"sr0_v4\",\n", + " generator_name='flat',\n", + " recoil=8,\n", + " mode='all')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c64b4064-041d-4e90-a0f3-df39bfbaa701", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "for i, run in enumerate(rn_runs):\n", + " events_simu_i = st_simu.get_array(run, 'event_info')\n", + " events_salt_i = st_salt.get_array(run, 'event_info')\n", + " events_data_i = st_data.get_array(run, 'event_info')\n", + " truth_i = st_simu.get_array(run, 'truth')\n", + " match_i = st_simu.get_array(run, 'match_acceptance_extended')\n", + " events_salt_matched_to_simu_i, \\\n", + " events_simu_matched_to_salt_i = saltax.match(truth_i[(truth_i['z']<-13)&(truth_i['z']>-145)&\n", + " (np.sqrt(truth_i['x']**2 + truth_i['y']**2)<63)], \n", + " match_i[(truth_i['z']<-13)&(truth_i['z']>-145)&\n", + " (np.sqrt(truth_i['x']**2 + truth_i['y']**2)<63)],\n", + " events_simu_i[(events_simu_i['z']<-13)&(events_simu_i['z']>-145)], \n", + " events_salt_i[(events_salt_i['z']<-13)&(events_salt_i['z']>-145)])\n", + " if i==0:\n", + " events_simu = events_simu_i\n", + " events_salt = events_salt_i\n", + " events_data = events_data_i\n", + " truth = truth_i\n", + " match = match_i\n", + " events_salt_matched_to_simu = events_salt_matched_to_simu_i\n", + " events_simu_matched_to_salt = events_simu_matched_to_salt_i\n", + " else:\n", + " events_simu = np.concatenate((events_simu,events_simu_i))\n", + " events_salt = np.concatenate((events_salt,events_salt_i))\n", + " events_data = np.concatenate((events_data,events_data_i))\n", + " truth = np.concatenate((truth,truth_i))\n", + " match = np.concatenate((match,match_i))\n", + " events_salt_matched_to_simu = np.concatenate((events_salt_matched_to_simu,\n", + " events_salt_matched_to_simu_i))\n", + " events_simu_matched_to_salt = np.concatenate((events_simu_matched_to_salt,\n", + " events_simu_matched_to_salt_i))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c344665-7059-4404-af19-4436e2ad64ee", + "metadata": {}, + "outputs": [], + "source": [ + "truth_s1 = []\n", + "truth_s2 = []\n", + "for i in range(np.max(truth_i['event_number'])):\n", + " selected_truth = truth_i[truth_i['event_number']==i]\n", + " if len(selected_truth) == 2 and len(np.unique(selected_truth['type']))==2:\n", + " if np.sqrt(selected_truth['x']**2 + selected_truth['y']**2)[0] <= 63:\n", + " truth_s1.append(selected_truth[selected_truth['type']==1][0]['raw_area_trigger'])\n", + " truth_s2.append(selected_truth[selected_truth['type']==2][0]['raw_area_trigger'])\n", + "\n", + "plt.scatter(truth_s1, truth_s2, s=1)\n", + "plt.xlim(0,150)\n", + "plt.ylim(0,6000)\n", + "plt.xlabel('True S1 [PE]')\n", + "plt.ylabel('True S2 [PE]')\n", + "plt.title('WFSim Truth ER Band')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f203e22-2502-4ea4-bd83-e9153853a043", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "for i, run in enumerate(bk_runs):\n", + " events_simu_i = st_simu.get_array(run, 'event_info')\n", + " events_salt_i = st_salt.get_array(run, 'event_info')\n", + " events_data_i = st_data.get_array(run, 'event_info')\n", + " truth_i = st_simu.get_array(run, 'truth')\n", + " match_i = st_simu.get_array(run, 'match_acceptance_extended')\n", + " events_salt_matched_to_simu_i, \\\n", + " events_simu_matched_to_salt_i = saltax.match(truth_i[(truth_i['z']<-13)&(truth_i['z']>-145)], \n", + " match_i[(truth_i['z']<-13)&(truth_i['z']>-145)],\n", + " events_simu_i[(events_simu_i['z']<-13)&(events_simu_i['z']>-145)], \n", + " events_salt_i[(events_salt_i['z']<-13)&(events_salt_i['z']>-145)])\n", + " if i == 0:\n", + " events_simub = events_simu_i\n", + " events_saltb = events_salt_i\n", + " events_datab = events_data_i\n", + " truthb = truth_i\n", + " matchb = match_i\n", + " events_salt_matched_to_simub = events_salt_matched_to_simu_i\n", + " events_simu_matched_to_saltb = events_simu_matched_to_salt_i\n", + " else:\n", + " events_simub = np.concatenate((events_simub,events_simu_i))\n", + " events_saltb = np.concatenate((events_saltb,events_salt_i))\n", + " events_datab = np.concatenate((events_datab,events_data_i))\n", + " truthb = np.concatenate((truthb,truth_i))\n", + " matchb = np.concatenate((matchb,match_i))\n", + " events_salt_matched_to_simub = np.concatenate((events_salt_matched_to_simub,\n", + " events_salt_matched_to_simu_i))\n", + " events_simu_matched_to_saltb = np.concatenate((events_simu_matched_to_saltb,\n", + " events_simu_matched_to_salt_i))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5903dc56-6cb1-4b25-9dfc-8e884ea79c81", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(dpi=150)\n", + "plt.scatter(events_simu['cs1'],events_simu['cs2'],s=1, label='Salted')\n", + "plt.scatter(events_data['cs1'],events_data['cs2'],s=1, label='Data')\n", + "plt.legend()\n", + "plt.xlim(0,100)\n", + "plt.ylim(0,6000)\n", + "plt.xlabel('CS1 [PE]')\n", + "plt.ylabel('CS2 [PE]')\n", + "plt.title('Salted Rn220 into SR0 Rn220 Calibration')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7891b6cc-2fd7-4e21-b9dd-a1aca03308f1", + "metadata": {}, + "outputs": [], + "source": [ + "n_bins = 31\n", + "cs1_bins = np.linspace(0,100,n_bins)\n", + "salt_med = []\n", + "simu_med = []\n", + "salt_1sig_u = []\n", + "simu_1sig_u = []\n", + "salt_1sig_l = []\n", + "simu_1sig_l = []\n", + "salt_2sig_u = []\n", + "simu_2sig_u = []\n", + "salt_2sig_l = []\n", + "simu_2sig_l = []\n", + "\n", + "for i in range(n_bins-1):\n", + " selected_salt = events_salt_matched_to_simu[(events_salt_matched_to_simu['cs1']>=cs1_bins[i])*\n", + " (events_salt_matched_to_simu['cs1']0]))\n", + " salt_1sig_l.append(np.percentile(selected_salt['cs2'][selected_salt['cs2']>0], 16.5))\n", + " salt_1sig_u.append(np.percentile(selected_salt['cs2'][selected_salt['cs2']>0], 83.5))\n", + " salt_2sig_l.append(np.percentile(selected_salt['cs2'][selected_salt['cs2']>0], 2.5))\n", + " salt_2sig_u.append(np.percentile(selected_salt['cs2'][selected_salt['cs2']>0], 97.5))\n", + " \n", + " selected_simu = events_simu_matched_to_salt[(events_simu_matched_to_salt['cs1']>=cs1_bins[i])*\n", + " (events_simu_matched_to_salt['cs1']0]))\n", + " simu_1sig_l.append(np.percentile(selected_simu['cs2'][selected_simu['cs2']>0], 16.5))\n", + " simu_1sig_u.append(np.percentile(selected_simu['cs2'][selected_simu['cs2']>0], 83.5))\n", + " simu_2sig_l.append(np.percentile(selected_simu['cs2'][selected_simu['cs2']>0], 2.5))\n", + " simu_2sig_u.append(np.percentile(selected_simu['cs2'][selected_simu['cs2']>0], 97.5))\n", + "\n", + "cs1_coord = np.linspace(0,100,n_bins-1) + 100/2/(n_bins-1)\n", + "\n", + "plt.figure(dpi=150)\n", + "plt.scatter(events_salt_matched_to_simu['cs1'],events_salt_matched_to_simu['cs2'],s=0.5, label='Salted', alpha=0.1)\n", + "plt.scatter(events_simu_matched_to_salt['cs1'],events_simu_matched_to_salt['cs2'],s=0.5, label='Simulated', alpha=0.1)\n", + "plt.plot(cs1_coord, salt_med, color='b', label='Salt Median')\n", + "plt.plot(cs1_coord, simu_med, color='r', label='Simu Median')\n", + "plt.plot(cs1_coord, salt_1sig_l, color='b', label='Salt 1sig',linestyle='dashed')\n", + "plt.plot(cs1_coord, salt_1sig_u, color='b', linestyle='dashed')\n", + "plt.plot(cs1_coord, salt_2sig_l, color='b', label='Salt 2sig',linestyle='dashed', alpha=0.5)\n", + "plt.plot(cs1_coord, salt_2sig_u, color='b', linestyle='dashed', alpha=0.5)\n", + "plt.plot(cs1_coord, simu_1sig_l, color='r', label='Simu 1sig',linestyle='dashed')\n", + "plt.plot(cs1_coord, simu_1sig_u, color='r', linestyle='dashed')\n", + "plt.plot(cs1_coord, simu_2sig_l, color='r', label='Simu 2sig',linestyle='dashed', alpha=0.5)\n", + "plt.plot(cs1_coord, simu_2sig_u, color='r', linestyle='dashed', alpha=0.5)\n", + "\n", + "\n", + "plt.legend()\n", + "plt.xlim(0,100)\n", + "plt.ylim(0,6000)\n", + "plt.xlabel('CS1 [PE]')\n", + "plt.ylabel('CS2 [PE]')\n", + "plt.title('Ambience Interference in SR0 Rn220')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9637ca13-f562-4c0a-9aec-dd7cdd252115", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(dpi=150)\n", + "plt.scatter(events_simu_matched_to_salt['cs1'],events_salt_matched_to_simu['cs1'], s=1, label='Rn220')\n", + "plt.scatter(events_simu_matched_to_saltb['cs1'],events_salt_matched_to_simub['cs1'], s=1, label='BKG')\n", + "plt.legend()\n", + "plt.xlim(0,150)\n", + "plt.ylim(0,150)\n", + "plt.title('Salted Rn220 into SR0 Rn220/BKG')\n", + "plt.xlabel('Simu CS1 [PE]')\n", + "plt.ylabel('Salt CS1 [PE]')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "313dd015-e521-4064-8b4a-54bd5e4f07b5", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(dpi=150)\n", + "plt.scatter(events_simu_matched_to_salt['cs2'],events_salt_matched_to_simu['cs2'], s=1, label='Rn220')\n", + "plt.scatter(events_simu_matched_to_saltb['cs2'],events_salt_matched_to_simub['cs2'], s=1, label='BKG')\n", + "plt.legend()\n", + "plt.xlim(0,6000)\n", + "plt.ylim(0,6000)\n", + "plt.title('Salted Rn220 into SR0 Rn220/BKG')\n", + "plt.xlabel('Simu CS2 [PE]')\n", + "plt.ylabel('Salt CS2 [PE]')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13e9f630-14f7-4d0b-a056-e804f99130ab", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(dpi=150)\n", + "plt.hist(events_salt_matched_to_simu['cs1'] - events_simu_matched_to_salt['cs1'], \n", + " bins=np.linspace(-200,200,50), alpha=0.5, label='Rn220', density=True)\n", + "plt.hist(events_salt_matched_to_simub['cs1'] - events_simu_matched_to_saltb['cs1'], \n", + " bins=np.linspace(-200,200,50), alpha=0.5, label='BKG', density=True)\n", + "plt.legend()\n", + "plt.title('Salted Rn220 into SR0 Rn220/BKG')\n", + "plt.xlabel('Ambience Contribution in CS1 [PE]')\n", + "plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc505f30-a053-40f2-af51-bc6434da8afb", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(dpi=150)\n", + "plt.hist(events_salt_matched_to_simu['cs2'] - events_simu_matched_to_salt['cs2'], \n", + " bins=np.linspace(-200,200,50), alpha=0.5, label='Rn220', density=True)\n", + "plt.hist(events_salt_matched_to_simub['cs2'] - events_simu_matched_to_saltb['cs2'], \n", + " bins=np.linspace(-200,200,50), alpha=0.5, label='BKG', density=True)\n", + "plt.title('Salted Rn220 into SR0 Rn220/BKG')\n", + "plt.legend()\n", + "plt.xlabel('Ambience Contribution in CS2 [PE]')\n", + "plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c91e880-28e0-4776-85d8-6ab6f5f88e3b", + "metadata": {}, + "outputs": [], + "source": [ + "inds_big_change_in_s2 = np.where((events_salt_matched_to_simu['cs2'] - events_simu_matched_to_salt['cs2']>100)&\n", + " (events_salt_matched_to_simu['cs2'] - events_simu_matched_to_salt['cs2']<200))[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f26b7651-b235-4d7e-9c10-b0a7fb81766b", + "metadata": {}, + "outputs": [], + "source": [ + "saltax.plot_wf(ind=inds_big_change_in_s2[7], st_salt=st_salt, st_simu=st_simu, st_data=st_data, \n", + " runid='022685', matched_simu=events_simu_matched_to_salt, \n", + " event_ext_window_ns=2.4e6, \n", + " s1_ext_window_samples=25, \n", + " s2_ext_window_samples=100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d8b48e1-116f-48d6-bb14-7cfba27ad2a8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:anaconda-XENONnT_development]", + "language": "python", + "name": "conda-env-anaconda-XENONnT_development-py" + }, + "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.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/saltax/match/match.py b/saltax/match/match.py index 9a97a9e..72450be 100644 --- a/saltax/match/match.py +++ b/saltax/match/match.py @@ -2,10 +2,11 @@ from tqdm import tqdm -def filter_out_missing_s1_s2(truth): +def filter_out_missing_s1_s2(truth, match): """ Filter out simulated events that have no S1 or S2 or both :param truth: truth from wfsim + :param match: match_acceptance_extended from pema :return: filtered truth """ bad_mask = np.zeros(len(truth), dtype=bool) @@ -29,13 +30,14 @@ def filter_out_missing_s1_s2(truth): print("Filter out %s percent of events due to missing S1 or S2 or both"\ % (np.sum(bad_mask)/len(truth)*100)) - return truth[~bad_mask] + return truth[~bad_mask], match[~bad_mask] -def filter_out_multiple_s1_s2(truth): +def filter_out_multiple_s1_s2(truth, match): """ Filter out simulated events that have multiple S1 or S2 reconstructed from wfsim. :param truth: truth from wfsim + :param match: match_acceptance_extended from pema :return: filtered truth """ bad_mask = np.zeros(len(truth), dtype=bool) @@ -52,7 +54,7 @@ def filter_out_multiple_s1_s2(truth): print("Filter out %s percent of events due to multiple S1 or S2"\ % (np.sum(bad_mask)/len(truth)*100)) - return truth[~bad_mask] + return truth[~bad_mask], match[~bad_mask] def filter_out_not_found(truth, match): @@ -66,9 +68,11 @@ def filter_out_not_found(truth, match): max_event_number = np.max(np.unique(truth['event_number'])) for event_id in np.arange(max_event_number)+1: - selected_truth = truth[truth['event_number']==event_id] + # Temporarily only match S1 because of the timing bug in wfsim + selected_truth = truth[(truth['event_number']==event_id)&(truth['type']==1)] indices = np.where(truth['event_number']==event_id) - selected_match = match[indices] + indices_s1 = np.where((truth['event_number']==event_id)&(truth['type']==1)) + selected_match = match[indices_s1] # The only outcome should be "found", otherwise remove the event if len(selected_match['outcome']) == 1: @@ -154,18 +158,18 @@ def pair_events_to_matched_simu(matched_simu, events): return matched_to -def pair_salt_to_simu(truth, events_simu, events_salt): +def pair_salt_to_simu(truth, match, events_simu, events_salt): """ Filter out bad simulation truth and then pair salted events to matched simulation events. :param truth: filtered truth + :param match: match_acceptance_extended from pema :param events_simu: events from wfsim :param events_salt: events from saltax after reconstruction - :return: ind_salt_matched_to_simu, ind_simu_matched_to_truth, truth_filtered + :return: ind_salt_matched_to_simu, ind_simu_matched_to_truth, truth_filtered, match_filtered """ - truth_filtered = filter_out_missing_s1_s2(truth) - truth_filtered = filter_out_multiple_s1_s2(truth_filtered) - # Temporarily turn off this filter because of wfsim bug in s2 timing - #truth_filtered, match_filtered = filter_out_not_found(truth_filtered, match_filtered) + truth_filtered, match_filtered = filter_out_missing_s1_s2(truth, match) + truth_filtered, match_filtered = filter_out_multiple_s1_s2(truth_filtered, match_filtered) + truth_filtered, match_filtered = filter_out_not_found(truth_filtered, match_filtered) ind_simu_matched_to_truth = pair_events_to_filtered_truth(truth_filtered, events_simu) events_simu_matched_to_truth = events_simu[ind_simu_matched_to_truth[ind_simu_matched_to_truth>=0]] @@ -173,10 +177,10 @@ def pair_salt_to_simu(truth, events_simu, events_salt): ind_salt_matched_to_simu = pair_events_to_matched_simu(events_simu_matched_to_truth, events_salt) - return ind_salt_matched_to_simu, ind_simu_matched_to_truth, truth_filtered + return ind_salt_matched_to_simu, ind_simu_matched_to_truth, truth_filtered, match_filtered -def match(truth, events_simu, events_salt): +def match(truth, match, events_simu, events_salt): """ Match salted events to simulation events. :param truth: truth from wfsim @@ -187,7 +191,7 @@ def match(truth, events_simu, events_salt): """ ind_salt_matched_to_simu, \ ind_simu_matched_to_truth, \ - _ = pair_salt_to_simu(truth, events_simu, events_salt) + _, _ = pair_salt_to_simu(truth, match, events_simu, events_salt) events_simu_matched_to_truth = events_simu[ind_simu_matched_to_truth[ind_simu_matched_to_truth>=0]] events_salt_matched_to_simu = events_salt[ind_salt_matched_to_simu[ind_salt_matched_to_simu>=0]]