LUXcode Dump

Added code for LUX

CodeForLUX/Matlab/MyFunctions/EstimateField.m

@@ -0,0 +1,97 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Purpose: 
+%   This function takes s2radius and s2drift inputs in the form of 1D vectors.  It
+%   uses Scott Hertel's field simulation to determine the field at the location of each data point.
+%   The output of this function is a 1D vector of the same length of the original data which 
+%   has corresponding field values for each data point in the units of V/cm.
+% 
+%   We use a cubic interpolation of Scott's map within the bounds of his simulated data
+%   and use a nearest neighbor interpolation outside of the bounds.  The bounds are defined by 
+%   fitting a spline to a rvZ scatterplot in the simulated data, with the additional constraint
+%   that any simulated data with drift_time<20 uSec is thrown out since the simualtion seems 
+%   inaccurate in that region.  Additionally, any real data that has radius < 1 cm is 
+%   considered to be outside of the simulation bounds, since a cubic interpolation fails in that region
+%   despite the data being within the simulation bounds.
+%
+%  Inputs:
+%       - s2radius             - 1xN Event radii after golden selection
+%       - s2drift              - 1xN Event drift times after golden selection
+%       - year                 - A string defining which of Scott's simulations to use.  ('2014' or '2015)
+%
+%  Outputs:
+%       - FieldValues          - 1xN Estimated field values in the vicinity of each event
+% 
+%  Author:
+%       - Richard Knoche
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ FieldValues ] = EstimateField( s2radius, s2drift, year )
+
+%% Load simualted field data from Scott
+
+switch(year)
+    case '2014'
+        load('C:\Program Files\MATLAB\R2012a\bin\Run04_ERBand\Kr2p22\sep2014_field_map.mat');
+    case '2015'
+        load('C:\Program Files\MATLAB\R2012a\bin\Run04_ERBand\Kr2p22\sep2015_field_map.mat')
+end
+
+
+%% Clean up simulated field data
+
+%reshape the data
+sim_radius  = reshape(sep2014_field_map.r_end_cm,1,size(sep2014_field_map.r_end_cm,1)*size(sep2014_field_map.r_end_cm,2));
+sim_drift   = reshape(sep2014_field_map.t_end_us,1,size(sep2014_field_map.t_end_us,1)*size(sep2014_field_map.t_end_us,2));
+sim_field   = reshape(sep2014_field_map.E_start_Vcm,1,size(sep2014_field_map.E_start_Vcm,1)*size(sep2014_field_map.E_start_Vcm,2));
+
+% Remove NaN from simulated field data
+nan_cut     = (~isnan(sim_radius) & ~isnan(sim_drift) & ~isnan(sim_field));
+sim_radius  = sim_radius(nan_cut);
+sim_drift   = sim_drift(nan_cut);
+sim_field   = sim_field(nan_cut);
+
+%% Make a 2D map of the field, with radius_bin and drift_bin vectors
+
+% Can use griddata to interpolate the Field Map, but it doesn't work well for points outside the boundary
+% Instead, we first find the boundary, use nearest neighbor interp outside of it, and cubic interp inside of it
+
+%Finds the boundary of the radius, drift_time scatter plot
+clear r_center max_z
+i       = 1;
+r_step  = 0.5;
+
+for r_max=[0.45:r_step:25];
+    r_center(i) = mean(sim_radius(inrange(sim_radius,[r_max-r_step,r_max]))); %#ok<AGROW>
+    max_z(i)    = max(sim_drift(inrange(sim_radius,[r_max-r_step,r_max]))); %#ok<AGROW>
+    i=i+1;
+end
+
+%Construct the boundary from a spline fit to the maximum values found above
+boundary_spline                  = csapi(r_center,max_z);
+
+%Evaluates if data is inside or outside the boundary
+data_boundary                    = fnval(boundary_spline, s2radius); 
+
+%A cut that is one in inside, zero if outside.  Has extra specifications to avoid bad simulation data at low drift and radius
+data_inside_bound                = (s2drift<=data_boundary & s2drift>=20 & s2radius>=1); 
+
+%To get field value for the data, first split into inside/outside bound data vectors
+data_inbound_radius              = s2radius(data_inside_bound);
+data_inbound_drift               = s2drift(data_inside_bound);
+data_outbound_radius             = s2radius(~data_inside_bound);
+data_outbound_drift              = s2drift(~data_inside_bound);
+
+%Interpolate: outside bound use nearest neighbor, inside bound use cubic
+field_inbound_value              = griddata(sim_drift,sim_radius,sim_field,data_inbound_drift,data_inbound_radius,'cubic');
+
+%only uses nearest for sim_drift>=20 uSec to avoid problem at low drift in sim
+field_outbound_value             = griddata(sim_drift(sim_drift>=20),sim_radius(sim_drift>=20),sim_field(sim_drift>=20),data_outbound_drift,data_outbound_radius,'nearest'); 
+
+%Patch inbound and outbound field values together
+FieldValues                      = zeros(length(field_inbound_value)+length(field_outbound_value),1);
+FieldValues(data_inside_bound)   = field_inbound_value;
+FieldValues(~data_inside_bound)  = field_outbound_value;
+
+end

CodeForLUX/Matlab/MyFunctions/GetGoldenCuts.m

@@ -0,0 +1,53 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Purpose: 
+%       This function produces golden event selection cuts between S1 and S2 area bounds
+%
+%  Inputs:
+%       - d                  - the data structure to be analysed.
+%       - s1area_bound_min   - Minimum S1 size
+%       - s1area_bound_max   - Maximum S1 size
+%       - s2area_bound_min   - Minimum S2 size
+%       - s2area_bound_max   - Maximum S2 size
+%
+%  Outputs:
+%       - s1_single_cut      -Golden event cut for S1 pulses
+%       - s2_single_cut      - Golden event cut for S2 pulses
+%
+%  Author:
+%       - Richard Knoche
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ s1_single_cut, s2_single_cut ] = GetGoldenCuts( d, s1area_bound_min, s1area_bound_max, s2area_bound_min, s2area_bound_max )
+
+%Pulse area cuts
+cut_s2_area = inrange(d.pulse_area_phe, [s2area_bound_min, s2area_bound_max]);
+cut_s1_area = inrange(d.pulse_area_phe, [s1area_bound_min, s1area_bound_max]);
+
+%Pulse classification cuts (S2 cut is defined above 100 phe, and paired with an s1)
+cut_pulse_s1 = d.pulse_classification == 1 | d.pulse_classification == 9;
+cut_pulse_s2 = d.pulse_classification == 2;
+cut_s2_with_threshold = d.pulse_area_phe.*cut_pulse_s2 > 100; % subset of cut_pulse_s2
+cut_legit_s2_in_legit_event = d.s1s2_pairing.*cut_s2_with_threshold; % this should be used as s2 classification cuts
+
+%Golden cut defines golden events to be events which have one and only one paired S2 above the threshold of 100 phe
+%Note there can be multiple S1s in the golden event
+cut_golden_event = sum(cut_legit_s2_in_legit_event) == 1; 
+cut_s2_in_golden_events = logical(repmat(cut_golden_event,[10,1]).*cut_legit_s2_in_legit_event); %Selects S2 that is in a golden event
+cut_s1_in_golden_events = logical(repmat(cut_golden_event,[10,1]).*cut_pulse_s1.*d.s1s2_pairing); %Selects first S1 that is in a golden event
+
+%Combine area cuts and golden cuts
+cut_s2_for = cut_s2_in_golden_events.*cut_s2_area; %Selects S2 that is in a golden event and in Kr area bounds
+cut_s1_for = cut_s1_in_golden_events.*cut_s1_area; %Selects first S1 that is in a golden event and in Kr area bounds
+
+%Require that "good" golden events have only one S1, that the S1 be within area bounds, and the S2 be within area bounds
+%Note sum(cut_s1_for) == 1 above only requires that the first of the S1 in an event be within area bounds, since the S1S2pairing part of cut_s1_in_golden_events is 0 for all subsequent S1s in the events
+cut_selected_events = sum(cut_s2_for) == 1 & sum(cut_s1_for) == 1 & sum(d.pulse_classification==1)==1; 
+
+%Final event selection cuts
+s1_single_cut = logical(repmat(cut_selected_events,[10,1]).*cut_s1_in_golden_events);
+s2_single_cut = logical(repmat(cut_selected_events,[10,1]).*cut_s2_in_golden_events);
+
+end
+

CodeForLUX/Matlab/MyFunctions/GetS1aS1bCuts.m

@@ -0,0 +1,33 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Purpose: 
+%       This function produces golden event selection cuts for Kr S1a and S1b pulses 
+%
+%  Inputs:
+%       - d                  - the data structure to be analysed.
+%       - s1_single_cut      - Golden event cut for merged Kr S1 pulses
+%       - s2_single_cut      - Golden event cut for merged Kr S2 pulses
+%
+%  Outputs:
+%       - s1ab_cut           - Golden event cut for S1a and S1b pulses
+%
+%  Author:
+%       - Richard Knoche
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ s1ab_cut ] = GetS1aS1bCuts(d, s1_single_cut, s2_single_cut)
+
+%Make timing variables
+s1ab_timing                        = d.Kr83fit_dt_samples;
+
+%Set up cuts
+s1ab_golden_cut                    = logical(sum(s1_single_cut(:,:),1));
+s1ab_pulsearea_cut                 = d.Kr83fit_s1b_area_phe > 30;
+s1ab_timing_cut                    = s1ab_timing > 13;  
+s1ab_radius_cut_temp               = (sqrt(d.x_cm.^2 + d.y_cm.^2) < 25) & s2_single_cut;
+s1ab_radius_cut                    = logical(sum(s1ab_radius_cut_temp(:,:),1));
+s1ab_cut                           = s1ab_golden_cut & s1ab_timing_cut & s1ab_pulsearea_cut & s1ab_radius_cut;
+
+end
+

CodeForLUX/Matlab/MyFunctions/GetS1aS1bFiducialCuts.m

@@ -0,0 +1,46 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Purpose: 
+%       This function produces fiducial volume cuts for Kr S1a and S1b pulses 
+%
+%  Inputs:
+%       - s1ab_z             - S1a/S1b z coordinates after golden cut is applied
+%       - s2radius           - S1a/S1b radius coordinates after golden cut is applied
+%
+%  Outputs:
+%       - s1ab_z_cut           - Fidiucal cut for making XY dependence maps with S1a/S1b
+%       - s1ab_r_cut           - Fidiucal cut for making Z dependence maps with S1a/S1b
+%       - s1ab_z_cut_upper     - Upper limit on S1a/S1b drift time histogram 
+%
+%  Author:
+%       - Richard Knoche
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ s1ab_z_cut, s1ab_r_cut, s1ab_z_cut_upper ] = GetS1aS1bFiducialCuts(s1ab_z, s1ab_radius)
+
+%Additional Z cut for s1ab_xy measurement
+s1ab_z_hist_binning                = (0:1:400); %Drift time range
+s1ab_z_hist_cut                    = inrange(s1ab_z,[0,400]);
+[s1ab_z_hist s1ab_z_hist_bins]     = hist(s1ab_z(s1ab_z_hist_cut),s1ab_z_hist_binning);
+frac_hist                          = cumsum(s1ab_z_hist)./sum(s1ab_z_hist);
+s1ab_z_cut_lower                   = min(s1ab_z_hist_bins(frac_hist>.10));
+s1ab_z_cut_upper                   = max(s1ab_z_hist_bins(frac_hist<.90));
+s1ab_z_cut                         = inrange(s1ab_z,[s1ab_z_cut_lower,s1ab_z_cut_upper]);
+
+%Additional radius cut for s1ab_Z measurement
+s1ab_r_hist_binning                = (0:1:26);
+s1ab_r_upperhist_cut               = (s1ab_z>s1ab_z_cut_upper & s1ab_radius<=26);
+s1ab_r_lowerhist_cut               = (s1ab_z<s1ab_z_cut_lower & s1ab_radius<=26);
+[upperz_r_hist upperz_r_hist_bins] = hist(s1ab_radius(s1ab_r_upperhist_cut),s1ab_r_hist_binning);
+[lowerz_r_hist lowerz_r_hist_bins] = hist(s1ab_radius(s1ab_r_lowerhist_cut),s1ab_r_hist_binning);
+frac_lowerz_r_hist                 = cumsum(lowerz_r_hist)./sum(lowerz_r_hist); %Find 80% radial edge of bottom 10% drift time
+lowerz_rbound                      = max(lowerz_r_hist_bins(frac_lowerz_r_hist<0.80));
+frac_upperz_r_hist                 = cumsum(upperz_r_hist)./sum(upperz_r_hist); %Find 90% radial edge of top 90% drift time
+upperz_rbound                      = max(upperz_r_hist_bins(frac_upperz_r_hist<0.90));   
+rslope                             = (s1ab_z_cut_lower-s1ab_z_cut_upper)/(lowerz_rbound-upperz_rbound); %find slope and intercept of radial cut line
+rintercept                         = s1ab_z_cut_lower-rslope.*lowerz_rbound;
+s1ab_r_cut                         = (s1ab_z-rslope.*s1ab_radius)<rintercept;    
+
+end
+

CodeForLUX/Matlab/MyFunctions/GetSECuts.m

@@ -0,0 +1,32 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Purpose: 
+%       This function produces single electron selection cuts from golden S1 and S2 events
+%
+%  Inputs:
+%       - d                  - the data structure to be analysed.
+%       - s1_single_cut      - 10xN Golden S1 selection cut
+%       - s2_single_cut      - 10xN Golden S2 selection cut
+%
+%  Outputs:
+%       - SE_good_cut      - Golden event cut for S1 pulses
+%
+%  Author:
+%       - Richard Knoche
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ SE_good_cut ] = GetSECuts( d, s1_single_cut, s2_single_cut)
+
+    %Pulse classification cut
+    s4_class=(d.pulse_classification==4); 
+
+    %Cut to remove S1 tails that improperly classified as SE
+    d = se_wrt_first_s1_func(d);
+    cut_se_s1_z = d.t_btw_se_first_s1 > 100*10; %Cut is based on timing after S1
+
+    %Final cut must be between the Kr S1 and S2, and meet the cuts above
+    SE_good_cut =logical( cut_se_s1_z & s4_class & ( cumsum(s1_single_cut)-cumsum(s2_single_cut) ) );
+
+end
+