Changeset 729bcf6 in sasview


Ignore:
Timestamp:
Mar 10, 2011 11:37:13 AM (14 years ago)
Author:
Jae Cho <jhjcho@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
510e7ad
Parents:
a234820
Message:

for circular and sectot average, calculate qvalues and dq values more accurately.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • DataLoader/manipulations.py

    r342a506 r729bcf6  
    210210             
    211211            #TODO: find better definition of x[i_q] based on q_data 
    212             x[i_q] = min + (i_q + 1) * self.bin_width / 2.0 
     212            x[i_q] += frac * q_value#min + (i_q + 1) * self.bin_width / 2.0 
    213213            y[i_q] += frac * data[npts] 
    214214             
    … …  
    226226        err_y = err_y / y_counts 
    227227        y    = y / y_counts 
    228          
     228        x    = x / y_counts 
    229229        idx = (numpy.isfinite(y) & numpy.isfinite(x))  
    230230         
    … …  
    426426        :return: Data1D object 
    427427        """ 
    428         # Get data 
     428        # Get data W/ finite values 
    429429        data = data2D.data[numpy.isfinite(data2D.data)] 
    430430        q_data = data2D.q_data[numpy.isfinite(data2D.data)] 
     431        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
    431432        err_data = data2D.err_data[numpy.isfinite(data2D.data)] 
     433         
    432434        dq_data = None 
     435         
     436        # Get the dq for resolution averaging 
    433437        if data2D.dqx_data != None and data2D.dqy_data != None: 
    434             dq_data = data2D.dqx_data[numpy.isfinite(data2D.data)] * \ 
    435                         data2D.dqx_data[numpy.isfinite(data2D.data)] 
    436             dq_data += data2D.dqy_data[numpy.isfinite(data2D.data)] * \ 
    437                         data2D.dqy_data[numpy.isfinite(data2D.data)] 
     438            # The pinholes and det. pix contribution present  
     439            # in both direction of the 2D which must be subtracted when 
     440            # converting to 1D: dq_overlap should calculated ideally at 
     441            # q = 0. Note This method works on only pinhole geometry.  
     442            # Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. 
     443            z_max = max(data2D.q_data) 
     444            z_min = min(data2D.q_data) 
     445            x_max = data2D.dqx_data[data2D.q_data[z_max]] 
     446            x_min = data2D.dqx_data[data2D.q_data[z_min]]     
     447            y_max = data2D.dqy_data[data2D.q_data[z_max]] 
     448            y_min = data2D.dqy_data[data2D.q_data[z_min]]           
     449            # Find qdx at q = 0 
     450            dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
     451            # when extrapolation goes wrong 
     452            if dq_overlap_x > min(data2D.dqx_data): 
     453                dq_overlap_x = min(data2D.dqx_data) 
     454            dq_overlap_x *= dq_overlap_x  
     455            # Find qdx at q = 0 
     456            dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
     457            # when extrapolation goes wrong 
     458            if dq_overlap_y > min(data2D.dqy_data): 
     459                dq_overlap_y = min(data2D.dqy_data) 
     460            # get dq at q=0. 
     461            dq_overlap_y *= dq_overlap_y 
     462 
     463            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y)/2.0) 
     464            # Final protection of dq 
     465            if dq_overlap < 0: 
     466                dq_overlap = y_min 
     467            dqx_data = data2D.dqx_data[numpy.isfinite(data2D.data)] 
     468            dqy_data = data2D.dqy_data[numpy.isfinite(data2D.data)] - dq_overlap 
     469            # def; dqx_data = dq_r dqy_data = dq_phi 
     470            # Convert dq 2D to 1D here 
     471            dqx = dqx_data * dqx_data             
     472            dqy = dqy_data * dqy_data 
     473            dq_data = numpy.add(dqx, dqy) 
    438474            dq_data = numpy.sqrt(dq_data) 
    439475             
    440         q_data_max = numpy.max(q_data) 
    441  
     476        #q_data_max = numpy.max(q_data) 
    442477        if len(data2D.q_data) == None: 
    443478            msg = "Circular averaging: invalid q_data: %g" % data2D.q_data 
    … …  
    475510                i_q = nbins -1               
    476511            y[i_q] += frac * data_n 
    477  
     512            # Take dqs from data to get the q_average 
     513            x[i_q] += frac * q_value 
    478514            if err_data == None or err_data[npt] == 0.0: 
    479515                if data_n < 0: 
    … …  
    483519                err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] 
    484520            if dq_data != None: 
    485                 err_x[i_q] += frac * frac * dq_data[npt] * dq_data[npt] 
     521                # To be consistent with dq calculation in 1d reduction,  
     522                # we need just the averages (not quadratures) because  
     523                # it should not depend on the number of the q points  
     524                # in the qr bins. 
     525                err_x[i_q] += frac * dq_data[npt] 
    486526            else: 
    487527                err_x = None 
    488528            y_counts[i_q]  += frac 
    489529         
    490         # We take the center of ring area, not radius.   
    491         # This is more accurate than taking the radial center of ring. 
    492         x = (qbins + self.bin_width) * (qbins + self.bin_width) 
    493         x += qbins * qbins 
    494         x = x / 2.0 
    495         x = numpy.sqrt(x) 
    496530        # Average the sums        
    497531        for n in range(nbins): 
    498532            if err_y[n] < 0: err_y[n] = -err_y[n] 
    499533            err_y[n] = math.sqrt(err_y[n]) 
    500             if err_x != None: 
    501                 err_x[n] = math.sqrt(err_x[n]) 
     534            #if err_x != None: 
     535            #    err_x[n] = math.sqrt(err_x[n]) 
    502536             
    503537        err_y = err_y / y_counts 
    504538        err_y[err_y==0] = numpy.average(err_y) 
    505539        y    = y / y_counts 
     540        x    = x / y_counts 
    506541        idx = (numpy.isfinite(y)) & (numpy.isfinite(x))  
    507542        if err_x != None: 
    … …  
    514549            raise ValueError, msg 
    515550         
    516         #elif len(y[idx])!= nbins: 
    517         #    print "resulted",nbins- len(y[idx]) 
    518         #,"empty bin(s) due to tight binning..." 
    519551        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) 
    520552     
    … …  
    746778        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 
    747779        dq_data = None 
    748         # dx (smear) data 
     780             
     781        # Get the dq for resolution averaging 
    749782        if data2D.dqx_data != None and data2D.dqy_data != None: 
    750             dq_data = data2D.dqx_data[numpy.isfinite(data2D.data)] * \ 
    751                         data2D.dqx_data[numpy.isfinite(data2D.data)] 
    752             dq_data += data2D.dqy_data[numpy.isfinite(data2D.data)] * \ 
    753                         data2D.dqy_data[numpy.isfinite(data2D.data)] 
     783            # The pinholes and det. pix contribution present  
     784            # in both direction of the 2D which must be subtracted when 
     785            # converting to 1D: dq_overlap should calculated ideally at 
     786            # q = 0.  
     787            # Extrapolate dqy(perp) at q = 0 
     788            z_max = max(data2D.q_data) 
     789            z_min = min(data2D.q_data) 
     790            x_max = data2D.dqx_data[data2D.q_data[z_max]] 
     791            x_min = data2D.dqx_data[data2D.q_data[z_min]]     
     792            y_max = data2D.dqy_data[data2D.q_data[z_max]] 
     793            y_min = data2D.dqy_data[data2D.q_data[z_min]]           
     794            # Find qdx at q = 0 
     795            dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
     796            # when extrapolation goes wrong 
     797            if dq_overlap_x > min(data2D.dqx_data): 
     798                dq_overlap_x = min(data2D.dqx_data) 
     799            dq_overlap_x *= dq_overlap_x  
     800            # Find qdx at q = 0 
     801            dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
     802            # when extrapolation goes wrong 
     803            if dq_overlap_y > min(data2D.dqy_data): 
     804                dq_overlap_y = min(data2D.dqy_data) 
     805            # get dq at q=0. 
     806            dq_overlap_y *= dq_overlap_y 
     807 
     808            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0)   
     809            if dq_overlap < 0: 
     810                dq_overlap = y_min 
     811            dqx_data = data2D.dqx_data[numpy.isfinite(data2D.data)] 
     812            dqy_data = data2D.dqy_data[numpy.isfinite(data2D.data)] - dq_overlap 
     813            # def; dqx_data = dq_r dqy_data = dq_phi 
     814            # Convert dq 2D to 1D here 
     815            dqx = dqx_data * dqx_data             
     816            dqy = dqy_data * dqy_data 
     817            dq_data = numpy.add(dqx, dqy) 
    754818            dq_data = numpy.sqrt(dq_data) 
    755819             
    … …  
    828892            ## Get the total y           
    829893            y[i_bin] += frac * data_n 
    830  
     894            x[i_bin] += frac * q_value 
    831895            if err_data[n] == None or err_data[n] == 0.0: 
    832896                if data_n < 0: 
    … …  
    837901                 
    838902            if dq_data != None: 
    839                 x_err[i_bin] += frac * frac * dq_data[n] * dq_data[n] 
     903                # To be consistent with dq calculation in 1d reduction,  
     904                # we need just the averages (not quadratures) because  
     905                # it should not depend on the number of the q points  
     906                # in the qr bins. 
     907                x_err[i_bin] += frac * dq_data[n] 
    840908            else: 
    841909                x_err = None 
    … …  
    846914            y[i] = y[i] / y_counts[i] 
    847915            y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 
    848             if x_err != None: 
    849                 x_err[i] = math.sqrt(x_err[i]) / y_counts[i] 
     916 
    850917            # The type of averaging: phi,q2, or q 
    851918            # Calculate x[i]should be at the center of the bin 
    … …  
    856923                # We take the center of ring area, not radius.   
    857924                # This is more accurate than taking the radial center of ring. 
    858                 delta_r = (self.r_max - self.r_min) / self.nbins 
    859                 r_inner = self.r_min + delta_r * i 
    860                 r_outer = r_inner + delta_r 
    861                 x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) 
    862                  
     925                #delta_r = (self.r_max - self.r_min) / self.nbins 
     926                #r_inner = self.r_min + delta_r * i 
     927                #r_outer = r_inner + delta_r 
     928                #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) 
     929                x[i] = x[i] / y_counts[i] 
    863930        y_err[y_err==0] = numpy.average(y_err)         
    864931        idx = (numpy.isfinite(y) & numpy.isfinite(y_err)) 
    865932        if x_err != None: 
    866             d_x = x_err[idx] 
     933            d_x = x_err[idx] / y_counts[idx] 
    867934        else: 
    868935            d_x = None 
Note: See TracChangeset for help on using the changeset viewer.