1212import numpy as np
1313from scipy .optimize import brentq
1414
15+ from lr_reduction .gravity_correction import GravityDirection , gravity_correction
1516from lr_reduction .instrument_settings import InstrumentSettings
1617from lr_reduction .utils import mantid_algorithm_exec
1718
@@ -326,6 +327,8 @@ class EventReflectivity:
326327 If ``UseDeadTimeThreshold`` is True, this is the upper limit for dead-time correction ratios
327328 use_emmission_time : bool
328329 If True, the emission time delay will be computed
330+ gravity_direction : GravityDirection
331+ Direction of gravity correction to apply
329332 """
330333
331334 QX_VS_QZ = 0
@@ -363,6 +366,7 @@ def __init__(
363366 dead_time_threshold : Optional [float ] = 1.5 ,
364367 instrument_settings : InstrumentSettings = None ,
365368 use_emission_time = True ,
369+ gravity_direction = None # undetermined
366370 ):
367371 if instrument in [self .INSTRUMENT_4A , self .INSTRUMENT_4B ]:
368372 self .instrument = instrument
@@ -393,6 +397,7 @@ def __init__(
393397 self .dead_time_threshold = dead_time_threshold
394398 self .instrument_settings = instrument_settings
395399 self .use_emission_time = use_emission_time
400+ self .grav_direction = gravity_direction
396401
397402 # Turn on functional background estimation
398403 self .use_functional_bck = functional_background
@@ -941,9 +946,17 @@ def _reflectivity(
941946 # convert tof values to wavelength.
942947 wl_list = tofs / self .constant
943948
944- # Gravity correction
945- d_theta = self .gravity_correction (ws , wl_list )
946- # collect weighted events
949+ # Gravity correction (check on value in template for gravity correction):
950+ if self .grav_direction is not None :
951+ d_theta = gravity_correction (ws , wl_list , gravity_direction = self .grav_direction )
952+ else :
953+ if peak_position == 0 : # is direct beam
954+ # for direct beam, find the direction using the `ths` value of the associated reflectivity run
955+ d_theta = gravity_correction (
956+ ws , wl_list , gravity_direction = GravityDirection .find_direction (self ._ws_sc ))
957+ else : # is reflectivity run
958+ d_theta = gravity_correction (ws , wl_list , gravity_direction = None )
959+
947960 event_weights = evt_list .getWeights ()
948961
949962 # Calculate per-spectum offset in theta for q-summing.
@@ -955,6 +968,7 @@ def _reflectivity(
955968 ths_value = ws .getRun ()["ths" ].value [- 1 ]
956969 delta_theta_f *= np .sign (ths_value )
957970
971+ # TODO: Check in code for any other calls of gravity_correction.
958972 # convert wavelengths into qz. This could be separated to enable output in lam and q.
959973 qz = 4.0 * np .pi / wl_list * np .sin (theta + delta_theta_f - d_theta )
960974 # Remove unfeasible negative values:
@@ -1224,75 +1238,6 @@ def emission_time_correction(self, ws, tofs):
12241238 tofs -= t_off + t_mult * tofs / self .constant
12251239 return tofs
12261240
1227- def gravity_correction (self , ws , wl_list ):
1228- """
1229- Gravity correction for each event
1230- Think this works on an array of wavelengths so could work for non-event list too.
1231-
1232- Parameters
1233- ----------
1234- ws : mantid.api.Workspace
1235- Mantid workspace to extract correction meta-data from.
1236- wl_list : numpy.ndarray
1237- Array of wavelengths for each event.
1238-
1239- Returns
1240- -------
1241- numpy.ndarray
1242- Array of gravity-corrected theta values for each event, in radians.
1243- """
1244- # Xi reference would be the position of xi if the si slit were to be positioned
1245- # at the sample. The distance from the sample to si is then xi_reference - xi.
1246- xi_reference = 445
1247- if ws .getInstrument ().hasParameter ("xi-reference" ):
1248- xi_reference = ws .getInstrument ().getNumberParameter ("xi-reference" )[0 ]
1249-
1250- # Distance between the s1 and the sample
1251- s1_sample_distance = 1485
1252- if ws .getInstrument ().hasParameter ("s1-sample-distance" ):
1253- s1_sample_distance = ws .getInstrument ().getNumberParameter ("s1-sample-distance" )[0 ] * 1000
1254-
1255- xi = 310
1256- if ws .getInstrument ().hasParameter ("BL4B:Mot:xi.RBV" ):
1257- xi = abs (ws .getRun ().getProperty ("BL4B:Mot:xi.RBV" ).value [0 ])
1258-
1259- sample_si_distance = xi_reference - xi
1260- slit_distance = s1_sample_distance - sample_si_distance
1261-
1262- # Angle of the incident beam on a horizontal sample
1263- # TODO: this will need to be updated with logged value.
1264- theta_in = - 4.0
1265-
1266- # Calculation from the ILL paper. This works for inclined beams.
1267- # Calculated theta is the angle on the sample
1268-
1269- g = 9.8067 # m/s^2
1270- h = 6.6260715e-34 # Js=kg m^2/s
1271- mn = 1.67492749804e-27 # kg
1272-
1273- v = h / (mn * wl_list * 1e-10 )
1274- k = g / (2 * v ** 2 )
1275-
1276- # Define the sample position as x=0, y=0. increasing x is towards moderator
1277- xs = 0
1278-
1279- # positions of slits
1280- x1 = sample_si_distance / 1000
1281- x2 = (sample_si_distance + slit_distance ) / 1000
1282-
1283- # height of slits determined by incident theta, y=0 is the sample height
1284- y1 = x1 * np .tan (theta_in * np .pi / 180 )
1285- y2 = x2 * np .tan (theta_in * np .pi / 180 )
1286-
1287- # This is the location of the top of the parabola
1288- x0 = (y1 - y2 + k * (x1 ** 2 - x2 ** 2 )) / (2 * k * (x1 - x2 ))
1289-
1290- # Angle is arctan(dy/dx) at sample
1291- theta_sample = np .arctan (2 * k * (x0 - xs )) * 180 / np .pi
1292-
1293- return (theta_sample - theta_in ) * np .pi / 180.0
1294-
1295-
12961241def compute_resolution (ws , default_dq = 0.027 , theta = None , q_summing = False ):
12971242 """
12981243 Compute the Q resolution from the meta data.
0 commit comments