1 | """ |
---|
2 | #This software was developed by the University of Tennessee as part of the |
---|
3 | #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
4 | #project funded by the US National Science Foundation. |
---|
5 | #See the license text in license.txt |
---|
6 | """ |
---|
7 | from __future__ import division |
---|
8 | |
---|
9 | import numpy as np |
---|
10 | from numpy import pi, cos, sin, sqrt |
---|
11 | |
---|
12 | from .resolution import Resolution |
---|
13 | |
---|
14 | ## Singular point |
---|
15 | SIGMA_ZERO = 1.0e-010 |
---|
16 | ## Limit of how many sigmas to be covered for the Gaussian smearing |
---|
17 | # default: 2.5 to cover 98.7% of Gaussian |
---|
18 | NSIGMA = 3.0 |
---|
19 | ## Defaults |
---|
20 | NR = {'xhigh':10, 'high':5, 'med':5, 'low':3} |
---|
21 | NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4} |
---|
22 | |
---|
23 | class Pinhole2D(Resolution): |
---|
24 | """ |
---|
25 | Gaussian Q smearing class for SAS 2d data |
---|
26 | """ |
---|
27 | |
---|
28 | def __init__(self, data=None, index=None, |
---|
29 | nsigma=NSIGMA, accuracy='Low', coords='polar'): |
---|
30 | """ |
---|
31 | Assumption: equally spaced bins in dq_r, dq_phi space. |
---|
32 | |
---|
33 | :param data: 2d data used to set the smearing parameters |
---|
34 | :param index: 1d array with len(data) to define the range |
---|
35 | of the calculation: elements are given as True or False |
---|
36 | :param nr: number of bins in dq_r-axis |
---|
37 | :param nphi: number of bins in dq_phi-axis |
---|
38 | :param coord: coordinates [string], 'polar' or 'cartesian' |
---|
39 | """ |
---|
40 | ## Accuracy: Higher stands for more sampling points in both directions |
---|
41 | ## of r and phi. |
---|
42 | ## number of bins in r axis for over-sampling |
---|
43 | self.nr = NR[accuracy.lower()] |
---|
44 | ## number of bins in phi axis for over-sampling |
---|
45 | self.nphi = NPHI[accuracy.lower()] |
---|
46 | ## maximum nsigmas |
---|
47 | self.nsigma = nsigma |
---|
48 | self.coords = coords |
---|
49 | self._init_data(data, index) |
---|
50 | |
---|
51 | def _init_data(self, data, index): |
---|
52 | """ |
---|
53 | Get qx_data, qy_data, dqx_data,dqy_data, |
---|
54 | and calculate phi_data=arctan(qx_data/qy_data) |
---|
55 | """ |
---|
56 | # TODO: maybe don't need to hold copy of qx,qy,dqx,dqy,data,index |
---|
57 | # just need q_calc and weights |
---|
58 | self.data = data |
---|
59 | self.index = index |
---|
60 | |
---|
61 | self.qx_data = data.qx_data[index] |
---|
62 | self.qy_data = data.qy_data[index] |
---|
63 | self.q_data = data.q_data[index] |
---|
64 | |
---|
65 | dqx = getattr(data, 'dqx_data', None) |
---|
66 | dqy = getattr(data, 'dqy_data', None) |
---|
67 | if dqx is not None and dqy is not None: |
---|
68 | # Here dqx and dqy mean dq_parr and dq_perp |
---|
69 | self.dqx_data = dqx[index] |
---|
70 | self.dqy_data = dqy[index] |
---|
71 | ## Remove singular points if exists |
---|
72 | self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO |
---|
73 | self.dqy_data[self.dqy_data < SIGMA_ZERO] = SIGMA_ZERO |
---|
74 | qx_calc, qy_calc, weights = self._calc_res() |
---|
75 | self.q_calc = [qx_calc, qy_calc] |
---|
76 | self.q_calc_weights = weights |
---|
77 | else: |
---|
78 | # No resolution information |
---|
79 | self.dqx_data = self.dqy_data = None |
---|
80 | self.q_calc = [self.qx_data, self.qy_data] |
---|
81 | self.q_calc_weights = None |
---|
82 | |
---|
83 | #self.phi_data = np.arctan(self.qx_data / self.qy_data) |
---|
84 | |
---|
85 | def _calc_res(self): |
---|
86 | """ |
---|
87 | Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, |
---|
88 | then find smeared intensity |
---|
89 | """ |
---|
90 | nr, nphi = self.nr, self.nphi |
---|
91 | # Total number of bins = # of bins |
---|
92 | nbins = nr * nphi |
---|
93 | # Number of bins in the dqr direction (polar coordinate of dqx and dqy) |
---|
94 | bin_size = self.nsigma / nr |
---|
95 | # in dq_r-direction times # of bins in dq_phi-direction |
---|
96 | # data length in the range of self.index |
---|
97 | nq = len(self.qx_data) |
---|
98 | |
---|
99 | # Mean values of dqr at each bins |
---|
100 | # starting from the half of bin size |
---|
101 | r = bin_size / 2.0 + np.arange(nr) * bin_size |
---|
102 | # mean values of qphi at each bines |
---|
103 | phi = np.arange(nphi) |
---|
104 | dphi = phi * 2.0 * pi / nphi |
---|
105 | dphi = dphi.repeat(nr) |
---|
106 | |
---|
107 | ## Transform to polar coordinate, |
---|
108 | # and set dphi at each data points ; 1d array |
---|
109 | dphi = dphi.repeat(nq) |
---|
110 | q_phi = self.qy_data / self.qx_data |
---|
111 | |
---|
112 | # Starting angle is different between polar |
---|
113 | # and cartesian coordinates. |
---|
114 | #if self.coords != 'polar': |
---|
115 | # dphi += np.arctan( q_phi * self.dqx_data/ \ |
---|
116 | # self.dqy_data).repeat(nbins).reshape(nq,\ |
---|
117 | # nbins).transpose().flatten() |
---|
118 | |
---|
119 | # The angle (phi) of the original q point |
---|
120 | q_phi = np.arctan(q_phi).repeat(nbins)\ |
---|
121 | .reshape(nq, nbins).transpose().flatten() |
---|
122 | ## Find Gaussian weight for each dq bins: The weight depends only |
---|
123 | # on r-direction (The integration may not need) |
---|
124 | weight_res = (np.exp(-0.5 * (r - bin_size / 2.0)**2) - |
---|
125 | np.exp(-0.5 * (r + bin_size / 2.0)**2)) |
---|
126 | # No needs of normalization here. |
---|
127 | #weight_res /= np.sum(weight_res) |
---|
128 | weight_res = weight_res.repeat(nphi).reshape(nr, nphi) |
---|
129 | weight_res = weight_res.transpose().flatten() |
---|
130 | |
---|
131 | ## Set dr for all dq bins for averaging |
---|
132 | dr = r.repeat(nphi).reshape(nr, nphi).transpose().flatten() |
---|
133 | ## Set dqr for all data points |
---|
134 | dqx = np.outer(dr, self.dqx_data).flatten() |
---|
135 | dqy = np.outer(dr, self.dqy_data).flatten() |
---|
136 | |
---|
137 | qx = self.qx_data.repeat(nbins)\ |
---|
138 | .reshape(nq, nbins).transpose().flatten() |
---|
139 | qy = self.qy_data.repeat(nbins)\ |
---|
140 | .reshape(nq, nbins).transpose().flatten() |
---|
141 | |
---|
142 | # The polar needs rotation by -q_phi |
---|
143 | if self.coords == 'polar': |
---|
144 | q_r = sqrt(qx**2 + qy**2) |
---|
145 | qx_res = ((dqx*cos(dphi) + q_r) * cos(-q_phi) |
---|
146 | + dqy*sin(dphi) * sin(-q_phi)) |
---|
147 | qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) |
---|
148 | + dqy*sin(dphi) * cos(-q_phi)) |
---|
149 | else: |
---|
150 | qx_res = qx + dqx*cos(dphi) |
---|
151 | qy_res = qy + dqy*sin(dphi) |
---|
152 | |
---|
153 | |
---|
154 | return qx_res, qy_res, weight_res |
---|
155 | |
---|
156 | def apply(self, theory): |
---|
157 | if self.q_calc_weights is not None: |
---|
158 | # TODO: interpolate rather than recomputing all the different qx,qy |
---|
159 | # Resolution needs to be applied |
---|
160 | nq, nbins = len(self.qx_data), self.nr * self.nphi |
---|
161 | ## Reshape into 2d array to use np weighted averaging |
---|
162 | theory = np.reshape(theory, (nbins, nq)) |
---|
163 | ## Averaging with Gaussian weighting: normalization included. |
---|
164 | value = np.average(theory, axis=0, weights=self.q_calc_weights) |
---|
165 | ## Return the smeared values in the range of self.index |
---|
166 | return value |
---|
167 | else: |
---|
168 | return theory |
---|