[3fdb4b6] | 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 |
---|