Ticket #7548: Sr122_individual_components_mg.py

File Sr122_individual_components_mg.py, 6.9 KB (added by Martyn Gigg, 7 years ago)
Line 
1try:
2    from mantidplottests import screenshot
3except ImportError:
4    def screenshot(*args,**kwargs):
5        pass
6
7from mantid.simpleapi import *
8import math
9import sys
10
11#=================== helper functions ===============================
12
13def create_cuboid_xml(xlength,ylength,zlength):
14    xml = """<cuboid id="sample0">
15<left-front-bottom-point x="%(xpt)f" y="-%(ypt)f" z="-%(zpt)f"  />
16<left-front-top-point  x="%(xpt)f" y="-%(ypt)f" z="%(zpt)f"  />
17<left-back-bottom-point  x="-%(xpt)f" y="-%(ypt)f" z="-%(zpt)f"  />
18<right-front-bottom-point  x="%(xpt)f" y="%(ypt)f" z="-%(zpt)f"  />
19</cuboid>
20<algebra val="sample0" />
21"""
22    return xml % {"xpt": xlength/2.0,"ypt":ylength/2.0,"zpt":zlength/2.0}
23 
24def plot_slice(in_ws, out_ws):
25        BinMD(InputWorkspace=in_ws, OutputWorkspace=out_ws,
26                AlignedDim0='[H,0,0], -12.000000, 9.000000, 100',
27                AlignedDim1='[0,K,0], -6.000000, 7.000000, 100',
28                AlignedDim2='[0,0,L], 0.000000, 6.000000, 1',
29                AlignedDim3='DeltaE, 100.000000, 150.000000, 1')
30        plotSlice(out_ws,colormax=180,normalization=0)
31
32#============================================================
33
34ei = 300.
35bins = [-30,3,279]
36temperature = 6.
37chopper_speed = 600.
38 
39# Oriented lattice & goniometer.
40alatt = 5.57
41blatt = 5.51
42clatt = 12.298
43uvec = [9.700000e-03,9.800000e-03,9.996000e-01]
44vvec = [9.992000e-01,-3.460000e-02,-4.580000e-02]
45
46omega = 0.0
47alpha = 0.0
48beta = 0.0
49gamma = 0.0
50 
51# sample dimensions
52sx = 0.05 # Perp
53sy = 0.025 # Up direction
54sz = 0.04 # Beam direction
55 
56# Crystal mosaic
57eta_sig = 4
58 
59data_dir = '/data/Software/VATES_testing/Sr122_tests/'
60
61 
62#Use Iliad processed workspace for simulation
63# Workspaces from Iliad have had their source component moved to redefine the location of t0 to the first monitor
64# It needs to be put back in the proper place here
65base_moderator = mtd['w1'].getInstrument().getBaseInstrument().getSource()
66mod_pos = base_moderator.getPos()
67MoveInstrumentComponent('w1',ComponentName=base_moderator.getName(),X=mod_pos.getX(), 
68                                             Y=mod_pos.getY(),Z=mod_pos.getZ(),RelativePosition=False)
69
70
71##
72## Required log entries, can be taken from real ones by placing an instrument parameter of the same
73## name pointing to the log name
74##
75AddSampleLog(Workspace='w1', LogName='Ei',LogText=str(ei), LogType="Number")
76AddSampleLog(Workspace='w1', LogName='temperature_log',LogText=str(temperature), LogType="Number")
77AddSampleLog(Workspace='w1', LogName='chopper_speed_log',LogText=str(chopper_speed), LogType="Number")
78AddSampleLog(Workspace='w1', LogName='eta_sigma',LogText=str(eta_sig), LogType="Number")
79 
80##
81## Sample shape
82##
83CreateSampleShape(InputWorkspace='w1', ShapeXML=create_cuboid_xml(sx,sy,sz))
84 
85##
86## Chopper & Moderator models.
87##
88CreateModeratorModel(Workspace='w1',ModelType='IkedaCarpenterModerator',
89                     Parameters="TiltAngle=32,TauF=2.7,TauS=0,R=0")
90CreateChopperModel(Workspace='w1',ModelType='FermiChopperModel',
91                   Parameters="AngularVelocity=chopper_speed_log,ChopperRadius=0.049,SlitThickness=0.0023,SlitRadius=1.3,Ei=Ei,JitterSigma=0.0")
92 
93##
94## UB matrix
95##
96SetUB(Workspace='w1',a=alatt,b=blatt,c=clatt,u=uvec,v=vvec)
97 
98##
99## Sample rotation. Simulate 1 run at zero degrees psi
100##
101
102psi = 0.0
103AddSampleLog(Workspace='w1',LogName='psi',LogText=str(psi),LogType='Number')
104SetGoniometer(Workspace='w1',Axis0="psi,0,1,0,1")
105
106# Create the MD workspace
107qscale = 'Q in A^-1'
108sr122_md = ConvertToMD(InputWorkspace='w1', QDimensions="Q3D", QConversionScales=qscale,SplitInto=[3], SplitThreshold=100,
109                      MinValues="-15,-15,-15,-30", MaxValues="25,25,25,279",OverwriteExisting=True)
110 
111 #Modified version of the simulation, with greater control of which contributions to the resolution model we have:
112resol_model = "TobyFitResolutionModel"
113fg_model = "Strontium122"
114                                                                 
115                                                 
116fg_parameters = "Seff=0.7,J1a=38.7,J1b=-5,J2=27.3,SJc=10,GammaSlope=0.08,MultEps=0,TwinType=1,FormFactorIon=0"
117xsec_model = "Strontium122"
118# Just forground model, no resolution convolution. Should reproduce exactly what Tobyfit does...
119# Sequentially work through the invidiual contributions to the resolution:
120res_parameters = "MCLoopMin=10,MCLoopMax=10,MCType=1,ForegroundOnly=1" 
121parameters = res_parameters +  "," + fg_parameters
122
123# Sequentially work through the invidiual contributions to the resolution:
124
125## Moderator ##
126res_parameters = "MCLoopMin=10,MCLoopMax=10,MCType=1,Moderator=1,Chopper=0,ChopperJitter=0,Aperture=0,SampleVolume=0,DetectorDepth=0,DetectorArea=0,DetectionTime=0,CrystalMosaic=0" 
127parameters = res_parameters +  "," + fg_parameters
128simulated_mod= SimulateResolutionConvolvedModel(InputWorkspace="sr122_md",
129                                            ResolutionFunction=resol_model,
130                                            ForegroundModel=xsec_model,Parameters=parameters)
131BinMD(InputWorkspace="simulated_mod", OutputWorkspace="simulated_mod_cut",
132                AlignedDim0='[H,0,0], 1, 1.3, 1',
133                AlignedDim1='[0,K,0], -1, 1, 41',
134                AlignedDim2='[0,0,L], -100, 100, 1',
135                AlignedDim3='DeltaE, 65, 75.000000, 1')
136mt_mod_table=QueryMDWorkspace('simulated_mod_cut')
137mt_mod_sim=ConvertTableToMatrixWorkspace('mt_mod_table',ColumnX='[0,K,0]/A^-1',ColumnY='Signal/none')
138
139
140## Sample Volume ##
141res_parameters = "MCLoopMin=10,MCLoopMax=10,MCType=1,Moderator=0,Chopper=0,ChopperJitter=0,Aperture=0,SampleVolume=1,DetectorDepth=0,DetectorArea=0,DetectionTime=0,CrystalMosaic=0" 
142parameters = res_parameters +  "," + fg_parameters
143simulated_vol= SimulateResolutionConvolvedModel(InputWorkspace="sr122_md",
144                                             ResolutionFunction=resol_model,
145                                             ForegroundModel=xsec_model,
146                                             Parameters=parameters)
147
148BinMD(InputWorkspace="simulated_vol", OutputWorkspace="simulated_vol_cut",
149                AlignedDim0='[H,0,0], 1, 1.3, 1',
150                AlignedDim1='[0,K,0], -1, 1, 41',
151                AlignedDim2='[0,0,L], -100, 100, 1',
152                AlignedDim3='DeltaE, 65, 75.000000, 1')
153mt_volume_table=QueryMDWorkspace('simulated_vol_cut')
154mt_volume_sim=ConvertTableToMatrixWorkspace('mt_volume_table',ColumnX='[0,K,0]/A^-1',ColumnY='Signal/none')
155
156## Mosaic ##
157res_parameters = "MCLoopMin=10,MCLoopMax=10,MCType=1,Moderator=0,Chopper=0,ChopperJitter=0,Aperture=0,SampleVolume=0,DetectorDepth=0,DetectorArea=0,DetectionTime=0,CrystalMosaic=1" 
158parameters = res_parameters +  "," + fg_parameters
159simulated_mosaic= SimulateResolutionConvolvedModel(InputWorkspace="sr122_md",
160                                             ResolutionFunction=resol_model,
161                                             ForegroundModel=xsec_model,
162                                             Parameters=parameters)
163
164BinMD(InputWorkspace="simulated_mosaic", OutputWorkspace="simulated_mosaic_cut",
165                AlignedDim0='[H,0,0], 1, 1.3, 1',
166                AlignedDim1='[0,K,0], -1, 1, 41',
167                AlignedDim2='[0,0,L], -100, 100, 1',
168                AlignedDim3='DeltaE, 65, 75.000000, 1')
169mt_mosaic_table=QueryMDWorkspace('simulated_mosaic_cut')
170mt_mosaic_sim=ConvertTableToMatrixWorkspace('mt_mosaic_table',ColumnX='[0,K,0]/A^-1',ColumnY='Signal/none')