source: sasview/sansmodels/src/c_smearer/smearer.cpp @ d560a37

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since d560a37 was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Reorganizing models in preparation of cpp cleanup

  • Property mode set to 100644
File size: 11.4 KB
Line 
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
6        If you use DANSE applications to do scientific research that leads to
7        publication, we ask that you acknowledge the use of the software with the
8        following sentence:
9
10        "This work benefited from DANSE software developed under NSF award DMR-0520547."
11
12        copyright 2009, University of Tennessee
13 */
14#include "smearer.hh"
15#include <stdio.h>
16#include <math.h>
17using namespace std;
18
19#if defined(_MSC_VER)
20extern "C" {
21#include "winFuncs.h"
22}
23#endif
24/**
25 * Constructor for BaseSmearer
26 *
27 * @param qmin: minimum Q value
28 * @param qmax: maximum Q value
29 * @param nbins: number of Q bins
30 */
31BaseSmearer :: BaseSmearer(double qmin, double qmax, int nbins) {
32        // Number of bins
33        this->nbins = nbins;
34        this->qmin = qmin;
35        this->qmax = qmax;
36        // Flag to keep track of whether we have a smearing matrix or
37        // whether we need to compute one
38        has_matrix = false;
39        even_binning = true;
40};
41
42/**
43 * Constructor for BaseSmearer
44 *
45 * Used for uneven binning
46 * @param q: array of Q values
47 * @param nbins: number of Q bins
48 */
49BaseSmearer :: BaseSmearer(double* q, int nbins) {
50        // Number of bins
51        this->nbins = nbins;
52        this->q_values = q;
53        // Flag to keep track of whether we have a smearing matrix or
54        // whether we need to compute one
55        has_matrix = false;
56        even_binning = false;
57};
58
59/**
60 * Constructor for SlitSmearer
61 *
62 * @param width: slit width in Q units
63 * @param height: slit height in Q units
64 * @param qmin: minimum Q value
65 * @param qmax: maximum Q value
66 * @param nbins: number of Q bins
67 */
68SlitSmearer :: SlitSmearer(double width, double height, double qmin, double qmax, int nbins) :
69        BaseSmearer(qmin, qmax, nbins){
70        this->height = height;
71        this->width = width;
72};
73
74/**
75 * Constructor for SlitSmearer
76 *
77 * @param width: slit width in Q units
78 * @param height: slit height in Q units
79 * @param q: array of Q values
80 * @param nbins: number of Q bins
81 */
82SlitSmearer :: SlitSmearer(double width, double height, double* q, int nbins) :
83        BaseSmearer(q, nbins){
84        this->height = height;
85        this->width = width;
86};
87
88/**
89 * Constructor for QSmearer
90 *
91 * @param width: array slit widths for each Q point, in Q units
92 * @param qmin: minimum Q value
93 * @param qmax: maximum Q value
94 * @param nbins: number of Q bins
95 */
96QSmearer :: QSmearer(double* width, double qmin, double qmax, int nbins) :
97        BaseSmearer(qmin, qmax, nbins){
98        this->width = width;
99};
100
101/**
102 * Constructor for QSmearer
103 *
104 * @param width: array slit widths for each Q point, in Q units
105 * @param q: array of Q values
106 * @param nbins: number of Q bins
107 */
108QSmearer :: QSmearer(double* width, double* q, int nbins) :
109        BaseSmearer(q, nbins){
110        this->width = width;
111};
112
113/**
114 * Compute the slit smearing matrix
115 *
116 * For even binning (q_min to q_max with nbins):
117 *
118 *   step = (q_max-q_min)/(nbins-1)
119 *   first bin goes from q_min to q_min+step
120 *   last bin goes from q_max to q_max+step
121 *
122 * For binning according to q array:
123 *
124 * Each q point represents a bin going from half the distance between it
125 * and the previous point to half the distance between it and the next point.
126 *
127 *    Example: bin i goes from (q_values[i-1]+q_values[i])/2 to (q_values[i]+q_values[i+1])/2
128 *
129 * The exceptions are the first and last bins, which are centered at the first and
130 * last q-values, respectively. The width of the first and last bins is the distance between
131 * their respective neighboring q-value.
132 */
133void SlitSmearer :: compute_matrix(){
134
135        weights = new vector<double>(nbins*nbins,0);
136
137        // Check the length of the data
138        if (nbins<2) return;
139        int npts_h = height>0.0 ? npts : 1;
140        int npts_w = width>0.0 ? npts : 1;
141
142        // If both height and width are great than zero,
143        // modify the number of points in each direction so
144        // that the total number of points is still what
145        // the user would expect (downgrade resolution)
146        //if(npts_h>1 && npts_w>1){
147        //      npts_h = (int)ceil(sqrt((double)npts));
148        //      npts_w = npts_h;
149        //}
150        double shift_h, shift_w, hbin_size, wbin_size;
151        // Make sure height and width are all positive (FWMH/2)
152        // Assumption; height and width are all same for all q points
153        if(npts_h == 1){
154                shift_h = 0.0;
155        } else {
156                shift_h = fabs(height);
157        }
158        if(npts_w == 1){
159                shift_w = 0.0;
160        } else {
161                shift_w = fabs(width);
162        }
163        // size of the h bin and w bin
164        hbin_size = shift_h / nbins;
165        wbin_size = shift_w / nbins;
166
167        // Loop over all q-values
168        double q, q_min, q_max, q_0=0.0;
169        for(int i=0; i<nbins; i++) {
170                // Find Weights
171                // Find q where the resolution smearing calculation of I(q) occurs
172                get_bin_range(i, &q, &q_min, &q_max);
173                // Block q becomes <=0
174                if (q <= 0){
175                        continue;
176                }
177                // Find q[0] value to normalize the weight later,
178                //  otherwise, we will have a precision problem.
179                if (i == 0){
180                        q_0 = q;
181                }
182                // Loop over all qj-values
183                for(int j=0; j<nbins; j++) {
184                        double q_j, q_high, q_low;
185                        // Calculate bin size of q_j
186                        get_bin_range(j, &q_j, &q_low, &q_high);
187                        // Block q_j becomes <=0
188                        if (q_j <= 0){
189                                continue;
190                        }
191                        // Check q_low that can not be negative.
192                        if (q_low < 0.0){
193                                q_low = 0.0;
194                        }
195                        // default parameter values
196                        (*weights)[i*nbins+j] = 0.0;
197                        // protect for negative q
198                        if (q <= 0.0 || q_j <= 0.0){
199                                        continue;
200                        }
201                        double shift_w = 0.0;
202                        // Condition: zero slit smear.
203                        if (npts_w == 1 && npts_h == 1){
204                                if(q_j == q) {
205                                        (*weights)[i*nbins+j] = 1.0;
206                                }
207                        }
208                        //Condition:Smear weight integration for width >0 when the height (=0) does not present.
209                        //Or height << width.
210                        else if((npts_w!=1 && npts_h == 1)|| (npts_w!=1 && npts_h != 1 && width/height > 100.0)){
211                                shift_w = width;
212                                //del_w = width/((double)npts_w-1.0);
213                                double q_shifted_low = q - shift_w;
214                                // High limit of the resolution range
215                                double q_shifted_high = q + shift_w;
216                                // Go through all the q_js for weighting those points
217                                if(q_j >= q_shifted_low && q_j <= q_shifted_high) {
218                                        // The weighting factor comes,
219                                        // Give some weight (delq_bin) for the q_j within the resolution range
220                                        // Weight should be same for all qs except
221                                        // for the q bin size at j.
222                                        // Note that the division by q_0 is only due to the precision problem
223                                        // where q_high - q_low gets to very small.
224                                        // Later, it will be normalized again.
225                                        (*weights)[i*nbins+j] += (q_high - q_low)/q_0 ;
226                                }
227                        }
228                        else{
229                                // Loop for width (;Height is analytical.)
230                                // Condition: height >>> width, otherwise, below is not accurate enough.
231                                // Smear weight numerical iteration for width >0 when the height (>0) presents.
232                                // When width = 0, the numerical iteration will be skipped.
233                                // The resolution calculation for the height is done by direct integration,
234                                // assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within a q' bin, [q_high, q_low].
235                                // In general, this weight numerical iteration for width >0 might be a rough approximation,
236                                // but it must be good enough when height >>> width.
237                                for(int k=(-npts_w + 1); k<npts_w; k++){
238                                        if(npts_w!=1){
239                                                shift_w = width/((double)npts_w-1.0)*(double)k;
240                                        }
241                                        // For each q-value, compute the weight of each other q-bin
242                                        // in the I(q) array
243                                        // Low limit of the resolution range
244                                        double q_shift = q + shift_w;
245                                        if (q_shift < 0.0){
246                                                q_shift = 0.0;
247                                        }
248                                        double q_shifted_low = q_shift;
249                                        // High limit of the resolution range
250                                        double q_shifted_high = sqrt(q_shift * q_shift + shift_h * shift_h);
251
252
253                                        // Go through all the q_js for weighting those points
254                                        if(q_j >= q_shifted_low && q_j <= q_shifted_high) {
255                                                // The weighting factor comes,
256                                                // Give some weight (delq_bin) for the q_j within the resolution range
257                                                // Weight should be same for all qs except
258                                                // for the q bin size at j.
259                                                // Note that the division by q_0 is only due to the precision problem
260                                                // where q_high - q_low gets to very small.
261                                                // Later, it will be normalized again.
262
263                                                // The fabs below are not necessary but in case: the weight should never be imaginary.
264                                                // At the edge of each sub_width. weight += u(at q_high bin) - u(0), where u(0) = 0,
265                                                // and weighted by (2.0* npts_w -1.0)once for each q.
266                                                //if (q == q_j) {
267                                                if (q_low <= q_shift && q_high > q_shift) {
268                                                        //if (k==0)
269                                                                (*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)-q_shift * q_shift)))/q_0;// * (2.0*double(npts_w)-1.0);
270                                                }
271                                                // For the rest of sub_width. weight += u(at q_high bin) - u(at q_low bin)
272                                                else{// if (u > 0.0){
273                                                        (*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)- q_shift * q_shift))-sqrt(fabs((q_low)*(q_low)- q_shift * q_shift)))/q_0 ;
274                                                }
275                                        }
276                                }
277                        }
278                }
279        }
280};
281
282/**
283 * Compute the point smearing matrix
284 */
285void QSmearer :: compute_matrix(){
286        weights = new vector<double>(nbins*nbins,0);
287
288        // Loop over all q-values
289        double q, q_min, q_max;
290        double q_j, q_jmax, q_jmin;
291        for(int i=0; i<nbins; i++) {
292                get_bin_range(i, &q, &q_min, &q_max);
293
294                for(int j=0; j<nbins; j++) {
295                        get_bin_range(j, &q_j, &q_jmin, &q_jmax);
296
297                        // Compute the fraction of the Gaussian contributing
298                        // to the q_j bin between q_jmin and q_jmax
299                        long double value =  erf( (q_jmax-q)/(sqrt(2.0)*width[i]) );
300                value -= erf( (q_jmin-q)/(sqrt(2.0)*width[i]) );
301                (*weights)[i*nbins+j] += value;
302                }
303        }
304}
305
306/**
307 * Computes the Q range of a given bin of the Q distribution.
308 * The range is computed according the the data distribution that
309 * was given to the object at initialization.
310 *
311 * @param i: number of the bin in the distribution
312 * @param q: q-value of bin i
313 * @param q_min: lower bound of the bin
314 * @param q_max: higher bound of the bin
315 *
316 */
317int BaseSmearer :: get_bin_range(int i, double* q, double* q_min, double* q_max) {
318        if (even_binning) {
319                double step = (qmax-qmin)/((double)nbins-1.0);
320                *q = qmin + (double)i*step;
321                *q_min = *q - 0.5*step;
322                *q_max = *q + 0.5*step;
323                return 1;
324        } else if (i>=0 && i<nbins) {
325                *q = q_values[i];
326                if (i==0) {
327                        double step = (q_values[1]-q_values[0])/2.0;
328                        *q_min = *q - step;
329                        *q_max = *q + step;
330                } else if (i==nbins-1) {
331                        double step = (q_values[i]-q_values[i-1])/2.0;
332                        *q_min = *q - step;
333                        *q_max = *q + step;
334                } else {
335                        *q_min = *q - (q_values[i]-q_values[i-1])/2.0;
336                        *q_max = *q + (q_values[i+1]-q_values[i])/2.0;
337                }
338                return 1;
339        }
340        return -1;
341}
342
343/**
344 * Perform smearing by applying the smearing matrix to the input Q array
345 */
346void BaseSmearer :: smear(double *iq_in, double *iq_out, int first_bin, int last_bin){
347
348        // If we haven't computed the smearing matrix, do it now
349        if(!has_matrix) {
350                compute_matrix();
351                has_matrix = true;
352        }
353
354        // Loop over q-values and multiply apply matrix
355        for(int q_i=first_bin; q_i<=last_bin; q_i++){
356                double sum = 0.0;
357                double counts = 0.0;
358
359                for(int i=first_bin; i<=last_bin; i++){
360                        // Skip if weight is less than 1e-03(this value is much smaller than
361                        // the weight at the 3*sigma distance
362                        // Will speed up a little bit...
363                        if ((*weights)[q_i*nbins+i] < 1.0e-003){
364                                continue;
365                        }
366                        sum += iq_in[i] * (*weights)[q_i*nbins+i];
367                        counts += (*weights)[q_i*nbins+i];
368                }
369
370                // Normalize counts
371                iq_out[q_i] = (counts>0.0) ? sum/counts : 0.0;
372        }
373}
Note: See TracBrowser for help on using the repository browser.