source: sasview/DataLoader/extensions/smearer.cpp @ cad821b

/**
        This software was developed by the University of Tennessee as part of the
        Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
        project funded by the US National Science Foundation.

        If you use DANSE applications to do scientific research that leads to
        publication, we ask that you acknowledge the use of the software with the
        following sentence:

        "This work benefited from DANSE software developed under NSF award DMR-0520547."

        copyright 2009, University of Tennessee
 */
#include "smearer.hh"
#include <stdio.h>
#include <math.h>
using namespace std;


/**
 * Constructor for BaseSmearer
 *
 * @param qmin: minimum Q value
 * @param qmax: maximum Q value
 * @param nbins: number of Q bins
 */
BaseSmearer :: BaseSmearer(double qmin, double qmax, int nbins) {
        // Number of bins
        this->nbins = nbins;
        this->qmin = qmin;
        this->qmax = qmax;
        // Flag to keep track of whether we have a smearing matrix or
        // whether we need to compute one
        has_matrix = false;
};

/**
 * Constructor for SlitSmearer
 *
 * @param width: slit width in Q units
 * @param height: slit height in Q units
 * @param qmin: minimum Q value
 * @param qmax: maximum Q value
 * @param nbins: number of Q bins
 */
SlitSmearer :: SlitSmearer(double width, double height, double qmin, double qmax, int nbins) :
        BaseSmearer(qmin, qmax, nbins){
        this->height = height;
        this->width = width;
};

        /**
         * Constructor for SlitSmearer
         *
         * @param width: array slit widths for each Q point, in Q units
         * @param qmin: minimum Q value
         * @param qmax: maximum Q value
         * @param nbins: number of Q bins
         */
QSmearer :: QSmearer(double* width, double qmin, double qmax, int nbins) :
        BaseSmearer(qmin, qmax, nbins){
        this->width = width;
};

/**
 * Compute the slit smearing matrix
 */
void SlitSmearer :: compute_matrix(){

        weights = new vector<double>(nbins*nbins,0);

        // Loop over all q-values
        for(int i=0; i<nbins; i++) {
                double q = qmin + (double)i*(qmax-qmin)/((double)nbins-1.0);

                // For each q-value, compute the weight of each other q-bin
                // in the I(q) array
                int npts_h = height>0 ? npts : 1;
                int npts_w = width>0 ? npts : 1;

                // If both height and width are great than zero,
                // modify the number of points in each direction so
                // that the total number of points is still what
                // the user would expect (downgrade resolution)
                if(npts_h>1 && npts_w>1){
                        npts_h = (int)ceil(sqrt((double)npts));
                        npts_w = npts_h;
                }

                double shift_h, shift_w;
                for(int k=0; k<npts_h; k++){
                        if(npts_h==1){
                                shift_h = 0;
                        } else {
                                shift_h = height/((double)npts_h-1.0) * (double)k;
                        }
                        for(int j=0; j<npts_w; j++){
                                if(npts_w==1){
                                        shift_w = 0;
                                } else {
                                        shift_w = width/((double)npts_w-1.0) * (double)j;
                                }
                                double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) );
                                int q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) ));

                                // Skip the entries outside our I(q) range
                                //TODO: be careful with edge effect
                                if(q_i<nbins)
                                        (*weights)[i*nbins+q_i]++;
                        }
                }
        }
};

/**
 * Compute the point smearing matrix
 */
void QSmearer :: compute_matrix(){
        weights = new vector<double>(nbins*nbins,0);

        // Loop over all q-values
        double step = (qmax-qmin)/((double)nbins-1.0);
        for(int i=0; i<nbins; i++) {
                double q = qmin + (double)i*step;
                double q_min = q - 0.5*step;
                double q_max = q + 0.5*step;

                for(int j=0; j<nbins; j++) {
                        double q_j = qmin + (double)j*step;

                        // Compute the fraction of the Gaussian contributing
                        // to the q bin between q_min and q_max
                        double value =  erf( (q_max-q_j)/(sqrt(2.0)*width[j]) );
                value -= erf( (q_min-q_j)/(sqrt(2.0)*width[j]) );
                (*weights)[i*nbins+j] += value;
                }
        }
}

/**
 * Perform smearing by applying the smearing matrix to the input Q array
 */
void BaseSmearer :: smear(double *iq_in, double *iq_out, int first_bin, int last_bin){

        // If we haven't computed the smearing matrix, do it now
        if(!has_matrix) {
                compute_matrix();
                has_matrix = true;
        }

        // Loop over q-values and multiply apply matrix
        for(int q_i=first_bin; q_i<=last_bin; q_i++){
                double sum = 0.0;
                double counts = 0.0;

                for(int i=first_bin; i<=last_bin; i++){
                        sum += iq_in[i] * (*weights)[q_i*nbins+i];
                        counts += (*weights)[q_i*nbins+i];
                }

                // Normalize counts
                iq_out[q_i] = (counts>0.0) ? sum/counts : 0;
        }
}