Added geometry correction. Work still in progress. This refs #4303

Code/Mantid/scripts/reduction/instruments/reflectometer/ref_reduction.py

@@ -26,7 +26,7 @@
 Xrange = [115, 210]
 
 #nexus_path = '/mnt/hgfs/j35/results/'
-nexus_path = '/home/m2d/data/ref_l/'
+nexus_path = '/home/j35/results/'
 list_data = ['66421',
              '66420',
              '66422',
@@ -153,9 +153,30 @@
 _tof_axis = mt2.readX(0)[:]
 ########## This was used to test the R(Q) 
 ##Convert the data without background subtraction to R(Q)
-q_array = wks_utility.convertToRvsQ(dMD=dMD,
-                                    theta=theta,
-                                    tof=_tof_axis)
+
+
+##without geometry correction
+#q_array = wks_utility.convertToRvsQ(dMD=dMD,
+#                                    theta=theta,
+#                                    tof=_tof_axis)
+
+##without geometry correction
+yrange = arange(int(list_data_peak[0][1]) - int(list_data_peak[0][0])) + int(list_data_peak[0][0])
+q_array = wks_utility.convertToRvsQWithCorrection(mt, 
+                                                  dMD=dMD,
+                                                  theta=theta,
+                                                  tof=_tof_axis,
+                                                  yrange = yrange,
+                                                  cpix=data_cpix)
+
+
+
+
+imshow(q_array)
+legend()
+show()
+
+sys.exit("stop")
 
 q_array_reversed = q_array[::-1]
 #if bPlot:

Code/Mantid/scripts/reduction/instruments/reflectometer/test.py

@@ -4,18 +4,85 @@
 from pylab import *
 import wks_utility
 
-#I'm testing the divide algorithm
-x_axis = [0,1,2,3,4]
-y_axis = [10,20,30,40,50,110,120,130,140,150]
-CreateWorkspace('wp1', DataX=x_axis, DataY=y_axis, Nspec=2)
-
-x_axis = [0,1,2,3,4]
-y_axis = [2,2,2,2,2,4,4,4,4,4]
-CreateWorkspace('wp2', DataX=x_axis, DataY=y_axis, Nspec=2)
-
-#divide
-Divide('wp1','wp2','wp3')
-mt3 = mtd['wp3']
-print mt3.readX(0)[:]
-print mt3.readY(0)[:]
-print mt3.readY(1)[:]
+data_file = '/mnt/hgfs/j35/results/REF_L_66420_event.nxs'
+LoadEventNexus(Filename=data_file, OutputWorkspace='EventData')
+rebin_parameters = '0,200,200000'  #microS
+rebin(InputWorkspace='EventData', OutputWorkspace='HistoData', Params=rebin_parameters)
+mt = mtd['HistoData']
+
+cpix = 256 / 2
+
+#dimension of the detector (256 by 304 pixels)
+maxX = 304
+maxY = 256
+
+sample = mt.getInstrument().getSample()
+source = mt.getInstrument().getSource()
+dSM = sample.getDistance(source)
+#create array of distances pixel->sample
+dPS_array = zeros((maxY, maxX))
+for x in range(maxX):
+    for y in range(maxY):
+        _index = maxY * x + y
+        detector = mt.getDetector(_index)
+        dPS_array[y, x] = sample.getDistance(detector)
+#array of distances pixel->source
+dMP_array = dPS_array + dSM
+#distance sample->center of detector
+dSD = dPS_array[maxY / 2, maxX / 2]
+
+#print wks_utility.ref_beamdiv_correct(cpix, mt, dSD)
+
+x = maxX/2.
+dPS_array = zeros(maxY)
+for y in range(maxY):
+    _index = maxY*x + y
+    detector = mt.getDetector(int(_index))
+    dPS_array[y] = sample.getDistance(detector)
+
+dSD = dPS_array[maxY / 2]
+dMD = dSD + dSM
+value = wks_utility.ref_beamdiv_correct(cpix, mt, dSD, 128)
+
+#Get metadata
+mt_run = mt.getRun()
+##get angles value
+ths_value = mt_run.getProperty('ths').value[0]
+ths_units = mt_run.getProperty('ths').units
+tthd_value = mt_run.getProperty('tthd').value[0]
+tthd_units = mt_run.getProperty('tthd').units
+ths_rad = wks_utility.angleUnitConversion(value=ths_value,
+                              from_units=ths_units,
+                              to_units='rad')
+tthd_rad = wks_utility.angleUnitConversion(value=tthd_value,
+                               from_units=tthd_units,
+                               to_units='rad')
+
+theta = tthd_rad - ths_rad
+_tof_axis = mt.readX(0)[:]
+
+q_array_before = wks_utility.convertToRvsQ(dMD=dMD,
+                                           theta=theta,
+                                           tof=_tof_axis)
+
+print 'q_array before'
+print q_array_before
+
+list_data_peak = [126, 134]
+yrange = arange(134-126+1) + 126
+
+q_array_after = wks_utility.convertToRvsQWithCorrection(mt,
+                                                         dMD=dMD,
+                                                         theta=theta,
+                                                         tof=_tof_axis,
+                                                         yrange=yrange,
+                                                         cpix=cpix)
+
+
+print 'q_array after'
+print q_array_after
+
+
+
+
+

Code/Mantid/scripts/reduction/instruments/reflectometer/wks_utility.py

@@ -1,4 +1,4 @@
-from numpy import zeros, arctan2, arange
+from numpy import zeros, arctan2, arange, shape
 from mantidsimple import *
 import math
 
@@ -235,6 +235,8 @@ def _convertToRvsQ(_tof_axis,
     return q_array
 
 
+
+
 def convertToRvsQ(dMD=-1,theta=-1,tof=None):
     """
     This function converts the pixel/TOF array to the R(Q) array
@@ -243,6 +245,7 @@ def convertToRvsQ(dMD=-1,theta=-1,tof=None):
             TOF: TOF of pixel
             theta: angle of detector
     """
+
     _const = float(4) * math.pi * m * dMD / h
     sz_tof = numpy.shape(tof)[0]
     q_array = zeros(sz_tof-1)
@@ -252,7 +255,55 @@ def convertToRvsQ(dMD=-1,theta=-1,tof=None):
         tofm = (tof1+tof2)/2.
         _Q = _const * math.sin(theta) / (tofm*1e-6)
         q_array[t] = _Q*1e-10
+    return q_array
+
+
+
+
+def convertToRvsQWithCorrection(mt, dMD=-1,theta=-1,tof=None, yrange=None, cpix=None):
+    """
+    This function converts the pixel/TOF array to the R(Q) array
+    using Q = (4.Pi.Mn)/h  *  L.sin(theta/2)/TOF
+    with    L: distance central_pixel->source
+            TOF: TOF of pixel
+            theta: angle of detector
+    """
+
+    sample = mt.getInstrument().getSample()
+    source = mt.getInstrument().getSource()
+    dSM = sample.getDistance(source)
+
+    maxX = 304
+    maxY = 256
 
+    dPS_array = zeros((maxY, maxX))
+    for x in range(maxX):
+        for y in range(maxY):
+            _index = maxY * x + y
+            detector = mt.getDetector(_index)
+            dPS_array[y, x] = sample.getDistance(detector)
+    #array of distances pixel->source
+    dMP_array = dPS_array + dSM
+    #distance sample->center of detector
+    dSD = dPS_array[maxY / 2, maxX / 2]
+
+    _const = float(4) * math.pi * m * dMD / h
+    sz_tof = len(tof)
+    q_array = zeros((len(yrange), sz_tof-1))
+
+    xrange = range(len(yrange))
+    for x in xrange:
+        _px = yrange[x]
+        dangle = ref_beamdiv_correct(cpix, mt, dSD, _px)
+        #theta += (2.0 * dangle)
+        theta += dangle
+        for t in range(sz_tof-1):
+            tof1 = tof[t]
+            tof2 = tof[t+1]
+            tofm = (tof1+tof2)/2.
+            _Q = _const * math.sin(theta) / (tofm*1e-6)
+            q_array[x,t] = _Q*1e-10
+
     return q_array
 
 def getQHisto(source_to_detector, theta, tof_array):
@@ -264,4 +315,169 @@ def getQHisto(source_to_detector, theta, tof_array):
         q_array[t] = _Q*1e-10
 
     return q_array
-
+
+
+def ref_beamdiv_correct(cpix, mt, det_secondary, 
+                        pixel_index, 
+                        pixel_width=0.0007):
+    """
+    This function calculates the acceptance diagram, determines pixel overlap
+    and computes the offset to the scattering angle.
+    """
+
+    # This is currently set to the same number for both REF_L and REF_M
+    epsilon = 0.5 * 1.3 * 1.0e-3
+
+    # Set the center pixel
+    if cpix is None:
+        cpix = 133.5
+
+    first_slit_size = getS1h(mt)
+    last_slit_size = getS2h(mt)
+    last_slit_dist = 0.654 #m
+    slit_dist = 2.193 #m
+
+    _y = 0.5 * (first_slit_size[0] + last_slit_size[0])
+    _x = slit_dist
+    gamma_plus = math.atan2( _y, _x)
+
+    _y = 0.5 * (first_slit_size[0] - last_slit_size[0])
+    _x = slit_dist 
+    gamma_minus = math.atan2( _y, _x)
+
+    half_last_aperture = 0.5 * last_slit_size[0]
+    neg_half_last_aperture = -1.0 * half_last_aperture
+
+    last_slit_to_det = last_slit_dist + det_secondary
+    dist_last_aper_det_sin_gamma_plus = last_slit_to_det * math.sin(gamma_plus)
+    dist_last_aper_det_sin_gamma_minus = last_slit_to_det * math.sin(gamma_minus)
+
+    #set the delta theta coordinates of the acceptance polygon
+    accept_poly_x = []
+    accept_poly_x.append(-1.0 * gamma_minus)
+    accept_poly_x.append(gamma_plus)
+    accept_poly_x.append(gamma_plus)
+    accept_poly_x.append(gamma_minus)
+    accept_poly_x.append(-1.0 * gamma_plus)
+    accept_poly_x.append(-1.0 * gamma_plus)
+    accept_poly_x.append(accept_poly_x[0])
+
+    #set the z coordinates of the acceptance polygon
+    accept_poly_y = []
+    accept_poly_y.append(half_last_aperture - dist_last_aper_det_sin_gamma_minus + epsilon)
+    accept_poly_y.append(half_last_aperture + dist_last_aper_det_sin_gamma_plus + epsilon)
+    accept_poly_y.append(half_last_aperture + dist_last_aper_det_sin_gamma_plus - epsilon)
+    accept_poly_y.append(neg_half_last_aperture + dist_last_aper_det_sin_gamma_minus - epsilon)
+    accept_poly_y.append(neg_half_last_aperture - dist_last_aper_det_sin_gamma_plus - epsilon)
+    accept_poly_y.append(neg_half_last_aperture - dist_last_aper_det_sin_gamma_plus + epsilon)
+    accept_poly_y.append(accept_poly_y[0])
+
+    cur_index = pixel_index
+
+    #set the z band for the pixel
+    xMinus = (cur_index - cpix - 0.5) * pixel_width
+    xPlus = (cur_index - cpix + 0.5) * pixel_width
+
+    #calculate the intersection
+    yLeftCross = -1
+    yRightCross = -1
+
+    xI = accept_poly_x[0]
+    yI = accept_poly_y[0]
+
+    int_poly_x = []
+    int_poly_y = []
+
+    for i in range(len(accept_poly_x)):
+
+        xF = accept_poly_y[i]
+        yF = accept_poly_x[i]
+
+        if xI < xF:
+
+            if xI < xMinus and xF > xMinus:
+                yLeftCross = yI + (yF - yI) * (xMinus - xI) / (xF - xI)
+                int_poly_x.append(yLeftCross)
+                int_poly_y.append(xMinus)
+
+            if xI < xPlus and xF >= xPlus:
+                yRightCross = yI + (yF - yI) * (xPlus - xI) / (xF - xI)
+                int_poly_x.append(yRightCross)
+                int_poly_y.append(xPlus)
+
+        else:
+
+            if xF < xPlus and xI >= xPlus:
+                yRightCross = yI + (yF - yI) * (xPlus - xI) / (xF - xI)
+                int_poly_x.append(yRightCross)
+                int_poly_y.append(xPlus)
+
+            if xF < xMinus and xI >= xMinus:
+                yLeftCross = yI + (yF - yI) * (xMinus - xI) / (xF - xI)
+                int_poly_x.append(yLeftCross)
+                int_poly_y.append(xMinus)
+
+        #This catches points on the polygon inside the range of interest
+        if xF >= xMinus and xF < xPlus:
+            int_poly_x.append(yF)
+            int_poly_y.append(xF)
+
+        xI = xF
+        yI = yF
+
+    if (len(int_poly_x) > 2):
+        int_poly_x.append(int_poly_x[0])
+        int_poly_y.append(int_poly_y[0])
+        int_poly_x.append(int_poly_x[1])
+        int_poly_y.append(int_poly_y[1])
+    else:
+        #Intersection polygon is null, point or line, so has no area
+        #therefore there is no angle corrction
+        return None
+
+    #Calculate intersection polygon area
+    area = calc_area_2D_polygon(int_poly_x, 
+                                int_poly_y,
+                                len(int_poly_x) - 2)
+
+
+    center_of_mass = calc_center_of_mass(int_poly_x, 
+                                         int_poly_y, 
+                                         area)
+
+    return center_of_mass
+
+def calc_area_2D_polygon(x_coord, y_coord, size_poly):
+    """
+    Calculation of the area defined by the 2D polygon
+    """
+    _range = arange(size_poly) + 1
+    area = 0
+    for i in _range:
+        area += (x_coord[i] * (y_coord[i+1] - y_coord[i-1]))
+    return area/2.
+
+def calc_center_of_mass(arr_x, arr_y, A):
+    """
+    Function that calculates the center-of-mass for the given polygon
+    
+    @param arr_x: The array of polygon x coordinates
+    @param arr_y: The array of polygon y coordinates
+    @param A: The signed area of the polygon
+    
+    @return: The polygon center-of-mass
+    """
+
+    center_of_mass = 0.0
+    SIXTH = 1./6.
+    for j in arange(len(arr_x)-2):
+        center_of_mass += (arr_x[j] + arr_x[j+1]) \
+                * ((arr_x[j] * arr_y[j+1]) - \
+                   (arr_x[j+1] * arr_y[j]))
+
+    if A != 0.0:
+        return (SIXTH * center_of_mass) / A
+    else:
+        return 0.0
+
+