forked from npshub/mantid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
testCenterREF_M.py
48 lines (40 loc) · 2.26 KB
/
testCenterREF_M.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
# Mantid Repository : https://github.com/mantidproject/mantid
#
# Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
# NScD Oak Ridge National Laboratory, European Spallation Source,
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
#pylint: disable=invalid-name
import numpy as np
import mantid
def calculateCenter(ws):
intensities=ws.extractY()[:,0]
detector_table=mantid.simpleapi.PreprocessDetectorsToMD(ws,GetMaskState="1")
l2=np.array(detector_table.column(1)) #in meters
tt=np.array(detector_table.column(2)) #in radians
phi=np.array(detector_table.column(3)) #in radians
masked=1-np.array(detector_table.column(7))
intensities=intensities*masked
x=l2*np.sin(tt)*np.cos(phi)
y=l2*np.sin(tt)*np.sin(phi)
z=l2*np.cos(tt)
avex=(x*intensities).sum()/intensities.sum()
avey=(y*intensities).sum()/intensities.sum()
avez=(z*intensities).sum()/intensities.sum()
rotation010=np.degrees(mantid.kernel.V3D(avex,0,avez).angle(mantid.kernel.V3D(0,0,1)))
return (-avey,-rotation010)
central=mantid.simpleapi.LoadEventNexus('REF_M_22715',NXentryName='entry-Off_Off')
original=mantid.simpleapi.CloneWorkspace(central)
translation,rotation=calculateCenter(original)
mantid.simpleapi.MoveInstrumentComponent(Workspace=central,ComponentName="DetectorArm",X=0,Y=translation,Z=1,RelativePosition=1)
mantid.simpleapi.RotateInstrumentComponent(Workspace=central,ComponentName="DetectorArm",X=0,Y=1,Z=0,Angle=rotation,RelativeRotation=1)
atangle=mantid.simpleapi.LoadEventNexus('REF_M_22710',NXentryName='entry-Off_Off')
mantid.simpleapi.MoveInstrumentComponent(Workspace=atangle,ComponentName="DetectorArm",X=0,Y=translation,Z=1,RelativePosition=1)
mantid.simpleapi.RotateInstrumentComponent(Workspace=atangle,ComponentName="DetectorArm",X=0,Y=1,Z=0,Angle=rotation,RelativeRotation=1)
central=mantid.simpleapi.ConvertUnits(central,Target="Wavelength",EMode="Elastic")
central=mantid.simpleapi.Rebin(central,"3.5,0.1,7")
atangle=mantid.simpleapi.ConvertUnits(atangle,Target="Wavelength",EMode="Elastic")
atangle=mantid.simpleapi.Rebin(atangle,"3.5,0.1,7")
csum=mantid.simpleapi.SumSpectra(central)
asum=mantid.simpleapi.SumSpectra(atangle)
normalized=asum/csum