1 | from scipy.special import erf |
---|
2 | from numpy import sqrt |
---|
3 | import numpy as np |
---|
4 | |
---|
5 | def pinhole_resolution(q_calc, q, dq): |
---|
6 | """ |
---|
7 | Compute the convolution matrix *W* for pinhole resolution 1-D data. |
---|
8 | |
---|
9 | Each row *W[i]* determines the normalized weight that the corresponding |
---|
10 | points *q_calc* contribute to the resolution smeared point *q[i]*. Given |
---|
11 | *W*, the resolution smearing can be computed using *dot(W,q)*. |
---|
12 | |
---|
13 | *q_calc* must be increasing. |
---|
14 | """ |
---|
15 | edges = bin_edges(q_calc) |
---|
16 | edges[edges<0.] = 0. # clip edges below zero |
---|
17 | G = erf( (edges[:,None]-q[None,:]) / (sqrt(2.0)*dq)[None,:] ) |
---|
18 | weights = G[1:,:] - G[:-1,:] |
---|
19 | weights /= sum(weights, axis=1) |
---|
20 | return weights |
---|
21 | |
---|
22 | def slit_resolution(q_calc, q, qx_width, qy_width): |
---|
23 | edges = bin_edges(q_calc) # Note: requires q > 0 |
---|
24 | edges[edges<0.] = 0. # clip edges below zero |
---|
25 | |
---|
26 | weights = np.zeros((len(q),len(q_calc)),'d') |
---|
27 | # Loop for width (;Height is analytical.) |
---|
28 | # Condition: height >>> width, otherwise, below is not accurate enough. |
---|
29 | # Smear weight numerical iteration for width >0 when the height (>0) presents. |
---|
30 | # When width = 0, the numerical iteration will be skipped. |
---|
31 | # The resolution calculation for the height is done by direct integration, |
---|
32 | # assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within a q' bin, [q_high, q_low]. |
---|
33 | # In general, this weight numerical iteration for width >0 might be a rough approximation, |
---|
34 | # but it must be good enough when height >>> width. |
---|
35 | E_sq = edges**2[:,None] |
---|
36 | y_pts = 500 if np.any(qy_width>0) else 1 |
---|
37 | for k in range(-y_pts+1,y_pts): |
---|
38 | qy = q if y_pts == 1 else q + qy_width/(y_pts-1)*k |
---|
39 | qy = np.clip(qy, 0.0, edges[-1]) |
---|
40 | qx_low = qy |
---|
41 | qx_high = sqrt(qx_low**2 + qx_width**2) |
---|
42 | in_x = (q_calc[:,None]>=qx_low[None,:])*(q_calc[:,None]<=qx_high[None,:]) |
---|
43 | weights += (sqrt(E_sq[1:]-qy[None,:]**2)-sqrt(E_sq[:-1]-qy[None,:]**2))*in_x |
---|
44 | weights /= sum(weights, axis=1) |
---|
45 | # Condition: zero slit smear. |
---|
46 | if (npts_w == 1 and npts_h == 1): |
---|
47 | if(q_j == q) : |
---|
48 | weights[i,j] = 1.0 |
---|
49 | #Condition:Smear weight integration for width >0 when the height (=0) does not present. |
---|
50 | #Or height << width. |
---|
51 | elif (npts_w!=1 and npts_h==1)or(npts_w!=1 and npts_h != 1 and width/height > 100.0): |
---|
52 | shift_w = width |
---|
53 | #del_w = width/((double)npts_w-1.0); |
---|
54 | q_shifted_low = q - shift_w |
---|
55 | # High limit of the resolution range |
---|
56 | q_shifted_high = q + shift_w |
---|
57 | # Go through all the q_js for weighting those points |
---|
58 | if(q_j >= q_shifted_low and q_j <= q_shifted_high): |
---|
59 | # The weighting factor comes, |
---|
60 | # Give some weight (delq_bin) for the q_j within the resolution range |
---|
61 | # Weight should be same for all qs except |
---|
62 | # for the q bin size at j. |
---|
63 | # Note that the division by q_0 is only due to the precision problem |
---|
64 | # where q_high - q_low gets to very small. |
---|
65 | # Later, it will be normalized again. |
---|
66 | weights[i,j] += (q_high - q_low)/q_0 |
---|
67 | else: |
---|
68 | # Loop for width (;Height is analytical.) |
---|
69 | # Condition: height >>> width, otherwise, below is not accurate enough. |
---|
70 | # Smear weight numerical iteration for width >0 when the height (>0) presents. |
---|
71 | # When width = 0, the numerical iteration will be skipped. |
---|
72 | # The resolution calculation for the height is done by direct integration, |
---|
73 | # assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within a q' bin, [q_high, q_low]. |
---|
74 | # In general, this weight numerical iteration for width >0 might be a rough approximation, |
---|
75 | # but it must be good enough when height >>> width. |
---|
76 | for k in range(-npts_w + 1,npts_w+1): |
---|
77 | if(npts_w!=1): |
---|
78 | shift_w = width/(npts_w-1.0)*k |
---|
79 | # For each q-value, compute the weight of each other q-bin |
---|
80 | # in the I(q) array |
---|
81 | # Low limit of the resolution range |
---|
82 | q_shift = q + shift_w |
---|
83 | if (q_shift < 0.0): |
---|
84 | q_shift = 0.0 |
---|
85 | q_shifted_low = q_shift |
---|
86 | # High limit of the resolution range |
---|
87 | q_shifted_high = sqrt(q_shift * q_shift + shift_h * shift_h) |
---|
88 | |
---|
89 | |
---|
90 | # Go through all the q_js for weighting those points |
---|
91 | if(q_j >= q_shifted_low and q_j <= q_shifted_high) : |
---|
92 | # The weighting factor comes, |
---|
93 | # Give some weight (delq_bin) for the q_j within the resolution range |
---|
94 | # Weight should be same for all qs except |
---|
95 | # for the q bin size at j. |
---|
96 | # Note that the division by q_0 is only due to the precision problem |
---|
97 | # where q_high - q_low gets to very small. |
---|
98 | # Later, it will be normalized again. |
---|
99 | |
---|
100 | # The fabs below are not necessary but in case: the weight should never be imaginary. |
---|
101 | # At the edge of each sub_width. weight += u(at q_high bin) - u(0), where u(0) = 0, |
---|
102 | # and weighted by (2.0* npts_w -1.0)once for each q. |
---|
103 | #if (q == q_j) { |
---|
104 | if (q_low <= q_shift and q_high > q_shift) : |
---|
105 | #if (k==0) |
---|
106 | weights[i,j] += (sqrt(abs((q_high)*(q_high)-q_shift * q_shift)))/q_0# * (2.0*double(npts_w)-1.0); |
---|
107 | # For the rest of sub_width. weight += u(at q_high bin) - u(at q_low bin) |
---|
108 | else:# if (u > 0.0){ |
---|
109 | weights[i,j] += (sqrt(abs((q_high)*(q_high)- q_shift * q_shift))-sqrt(abs((q_low)*(q_low)- q_shift * q_shift)))/q_0 |
---|
110 | |
---|
111 | |
---|
112 | def bin_edges(x): |
---|
113 | if len(x) < 2 or (np.diff(x)<0).any(): |
---|
114 | raise ValueError("Expected bins to be an increasing set") |
---|
115 | edges = np.hstack([ |
---|
116 | x[0] - 0.5*(x[1] - x[0]), # first point minus half first interval |
---|
117 | 0.5*(x[1:] + x[:-1]), # mid points of all central intervals |
---|
118 | x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval |
---|
119 | ]) |
---|
120 | return edges |
---|