1 | import datetime |
---|
2 | #LOQ data analysis script |
---|
3 | #path = "C:/MantidInstall/" |
---|
4 | #path = "C:/Mantid/Test/" |
---|
5 | path="C:/Documents and Settings/wmx35332/Mantid/Test/" |
---|
6 | # this makes a list for a rectangular block (wide x high) from pixel startID [with (1,1) at bottom left in usual LOQ sense |
---|
7 | # though may need to change for SANS2d ??? ] for square array dim x dim |
---|
8 | def detBlock(startID,wide,high,dim=128): |
---|
9 | output = "" |
---|
10 | # no idea why this says zero here rather than one ! - but it works |
---|
11 | for j in range(0,high): |
---|
12 | for i in range(0,wide): |
---|
13 | output += str(startID+i+(dim*j))+"," |
---|
14 | output = output.rstrip(",") |
---|
15 | return output |
---|
16 | # |
---|
17 | ####################### |
---|
18 | print datetime.datetime.now() |
---|
19 | for x in range(0,1): |
---|
20 | # Call the LOQScriptInput algorithm to get the input data and parameters |
---|
21 | #input = LOQScriptInputDialog("48098","48094","48130","48127","48128",path+"Data/DIRECT.041","38","419","2.2","10.0","-0.035","0.008","0.28","0.002","324.95","328.02") |
---|
22 | input = LOQScriptInput("48098","48094","48130","48127","48128",path+"Data/DIRECT.041","38","419","2.2","10.0","-0.035","0.008","0.28","0.002","324.95","328.02") |
---|
23 | #,path+"Data/LOQ sans configuration/DIRECT.041") |
---|
24 | data_file = path + "Data/LOQ" + input.getPropertyValue("SampleWorkspace") + ".raw" |
---|
25 | empty_file = path + "Data/LOQ" + input.getPropertyValue("EmptyCanWorkspace") + ".raw" |
---|
26 | trans_sample_file = path + "Data/LOQ" + input.getPropertyValue("TransmissionSampleWorkspace") + ".raw" |
---|
27 | trans_direct_file = path + "Data/LOQ" + input.getPropertyValue("TransmissionDirectWorkspace") + ".raw" |
---|
28 | trans_empty_file = path + "Data/LOQ" + input.getPropertyValue("TransmissionEmptyCanWorkspace") + ".raw" |
---|
29 | min_radius = float(input.getPropertyValue("RadiusMin")) |
---|
30 | max_radius = float(input.getPropertyValue("RadiusMax")) |
---|
31 | wav1 = float(input.getPropertyValue("WavelengthMin")) |
---|
32 | wav2 = float(input.getPropertyValue("WavelengthMax")) |
---|
33 | dwav = input.getPropertyValue("WavelengthDelta") |
---|
34 | q1 = input.getPropertyValue("QMin") |
---|
35 | q2 = input.getPropertyValue("Q_max") |
---|
36 | dq = input.getPropertyValue("Q_delta") |
---|
37 | xshift=str((317.5-float(input.getPropertyValue("BeamCentreX")))/1000.0) |
---|
38 | yshift=str((317.5-float(input.getPropertyValue("BeamCentreY")))/1000.0) |
---|
39 | direct_beam_file = input.getPropertyValue("EfficiencyCorrectionFile") |
---|
40 | |
---|
41 | outputWS = "Small_Angle" |
---|
42 | |
---|
43 | ####################### |
---|
44 | #Step 1 - Load the data file |
---|
45 | LoadDataAlg = LoadRaw(data_file,OutputWorkspace="Monitor",SpectrumMin="2",SpectrumMax="2") |
---|
46 | #data_file = LoadDataAlg.getPropertyValue("Filename") |
---|
47 | #Load the small angle bank |
---|
48 | # whole det is 3 to 16386, but skip first and last two rows is 130 to 16130, note the crop on flat_cell needs these numbers |
---|
49 | firstsmall=130 |
---|
50 | lastsmall=16130 |
---|
51 | LoadRaw(Filename = data_file, OutputWorkspace=outputWS,SpectrumMin=str(firstsmall),SpectrumMax=str(lastsmall)) |
---|
52 | #Load the high angle bank |
---|
53 | #LoadRaw(Filename = data_file, OutputWorkspace="High_Angle",spectrummin="16386",spectrummax="17792") |
---|
54 | |
---|
55 | ####################### |
---|
56 | #Step 1.2 - Correct the monitor spectrum for flat background and remove prompt spike |
---|
57 | # Interpolation is just linear for now |
---|
58 | RemoveBins("Monitor","Monitor","19900","20500",Interpolation="Linear") |
---|
59 | FlatBackground("Monitor","Monitor","0","31000","39000") |
---|
60 | |
---|
61 | ####################### |
---|
62 | #Step 1.3 - Mask out unwanted detetctors |
---|
63 | # XML description of a cylinder containing detectors to remove |
---|
64 | inner = "<infinite-cylinder id='inner'> " |
---|
65 | inner += "<centre x='0.0' y='0.0' z='0.0' /> " |
---|
66 | inner += "<axis x='0.0' y='0.0' z='1' /> " |
---|
67 | inner += "<radius val='"+str(min_radius/1000.0)+"' /> " |
---|
68 | inner += "</infinite-cylinder> " |
---|
69 | inner += "<algebra val='inner' /> " |
---|
70 | |
---|
71 | # Remove the beam stop |
---|
72 | MaskDetectorsInShape(outputWS,inner) |
---|
73 | |
---|
74 | outer = "<infinite-cylinder id='outer'> " |
---|
75 | outer += "<centre x='0.0' y='0.0' z='0.0' /> " |
---|
76 | outer += "<axis x='0.0' y='0.0' z='1' /> " |
---|
77 | outer += "<radius val='"+str(max_radius/1000.0)+"' /> " |
---|
78 | outer += "</infinite-cylinder> " |
---|
79 | # The hash in front of the name says that we want to keep everything inside the cylinder |
---|
80 | outer += "<algebra val='#outer' /> " |
---|
81 | |
---|
82 | # Remove the corners |
---|
83 | MaskDetectorsInShape(outputWS,outer) |
---|
84 | |
---|
85 | # mask right hand column |
---|
86 | MaskDetectors(outputWS,DetectorList=detBlock(128+2,1,128)) |
---|
87 | |
---|
88 | ####################### |
---|
89 | #Step 1.1 - try move detector for beam centre (324.95, 328.02), |
---|
90 | #This step has to happen after the detector masking |
---|
91 | #RelativePosition="0" is a logical variable, which means an "absolute shift" |
---|
92 | # the detector centr is at X=Y=0 so does relative or absolute matter ? |
---|
93 | MoveInstrumentComponent(outputWS,"main-detector-bank",X=xshift,Y=yshift,RelativePosition="1") |
---|
94 | |
---|
95 | ####################### |
---|
96 | #Step 2 - Convert all of the files to wavelength and rebin |
---|
97 | # ConvertUnits does have a rebin option, but it's crude. In particular it rebins on linear scale. |
---|
98 | wavbin=str(wav1)+","+str(dwav)+","+str(wav2) |
---|
99 | wavname="_"+str(int(wav1))+"_"+str(int(wav2)) |
---|
100 | # |
---|
101 | ConvertUnits("Monitor","Monitor","Wavelength") |
---|
102 | Rebin("Monitor","Monitor",wavbin) |
---|
103 | ConvertUnits(outputWS,outputWS,"Wavelength") |
---|
104 | Rebin(outputWS,outputWS,wavbin) |
---|
105 | #ConvertUnits("High_Angle","High_Angle","Wavelength") |
---|
106 | #Rebin("High_Angle","High_Angle",wavbin) |
---|
107 | |
---|
108 | # FROM NOW ON JUST LOOKING AT MAIN DETECTOR DATA..... |
---|
109 | |
---|
110 | ####################### |
---|
111 | #step 2.5 remove unwanted Wavelength bins |
---|
112 | #Remove non edge bins and linear interpolate #373 |
---|
113 | #RemoveBins("High_Angle","High_Angle","3.00","3.50",Interpolation="Linear") |
---|
114 | #Future Missing - add cubic interpolation |
---|
115 | |
---|
116 | ####################### |
---|
117 | #Step 3 - Flat cell correction |
---|
118 | |
---|
119 | #OPTION 1 |
---|
120 | #calculate from raw flood source file |
---|
121 | #maybe later! |
---|
122 | |
---|
123 | #OPTION 2 |
---|
124 | #LoadRKH with scalar value for wavelength ranges |
---|
125 | #data/flat(wavelength) |
---|
126 | #LoadRKH(path+"Data/LOQ sans configuration/FLAT_CELL.061","flat","SpectraNumber") |
---|
127 | #CropWorkspace("flat","flat",StartSpectrum=str(firstsmall),EndSpectrum=str(lastsmall)) |
---|
128 | #optional correct for diffrernt detector positions |
---|
129 | #SolidAngle("flat","flatsolidangle") |
---|
130 | #SolidAngle(outputWS,"solidangle") |
---|
131 | #Divide("flatsolidangle","solidangle","SAcorrection") |
---|
132 | #Divide("flat","SAcorrection","flat") |
---|
133 | #RETRY |
---|
134 | #Divide(outputWS,"flat",outputWS) |
---|
135 | |
---|
136 | #OPTION 3 |
---|
137 | #Correction using instrument geometry |
---|
138 | #SolidAngle(outputWS,"solidangle") |
---|
139 | #Divide(outputWS,"solidangle","Small_Angle corrected by solidangle") |
---|
140 | |
---|
141 | ####################### |
---|
142 | #Step 4 - Correct by incident beam monitor |
---|
143 | # At this point need to fork off workspace name to keep a workspace containing raw counts |
---|
144 | outputWS_cor = outputWS + "_tmp" |
---|
145 | Divide(outputWS,"Monitor",outputWS_cor) |
---|
146 | mantid.deleteWorkspace("Monitor") |
---|
147 | |
---|
148 | ####################### |
---|
149 | #Step 6 - Correct by transmission |
---|
150 | # Load the run with sample |
---|
151 | #LoadRawDialog(path+"Data/LOQ48130.raw","sample") |
---|
152 | LoadRaw(Filename=trans_sample_file,OutputWorkspace="sample") |
---|
153 | # Change the instrument definition to the correct one |
---|
154 | LoadInstrument("sample",path+"Instrument/LOQ_trans_Definition.xml") |
---|
155 | # Need to remove prompt spike and, later, flat background |
---|
156 | RemoveBins("sample","sample","19900","20500",Interpolation="Linear") |
---|
157 | FlatBackground("sample","sample","1,2","31000","39000") |
---|
158 | ConvertUnits("sample","sample","Wavelength") |
---|
159 | Rebin("sample","sample",wavbin) |
---|
160 | # Now load and convert the direct run |
---|
161 | #LoadRawDialog(path+"Data/LOQ48127.raw","direct") |
---|
162 | LoadRaw(Filename=trans_direct_file,OutputWorkspace="direct") |
---|
163 | LoadInstrument("direct",path+"Instrument/LOQ_trans_Definition.xml") |
---|
164 | RemoveBins("direct","direct","19900","20500",Interpolation="Linear") |
---|
165 | FlatBackground("direct","direct","1,2","31000","39000") |
---|
166 | ConvertUnits("direct","direct","Wavelength") |
---|
167 | Rebin("direct","direct",wavbin) |
---|
168 | |
---|
169 | CalculateTransmission("sample","direct","transmission") |
---|
170 | mantid.deleteWorkspace("sample") |
---|
171 | mantid.deleteWorkspace("direct") |
---|
172 | |
---|
173 | # Now do the correction |
---|
174 | Divide(outputWS_cor,"transmission",outputWS_cor) |
---|
175 | #mantid.deleteWorkspace("transmission") |
---|
176 | |
---|
177 | ####################### |
---|
178 | #Step 7 - Correct for efficiency |
---|
179 | CorrectToFile(outputWS_cor,direct_beam_file,outputWS_cor,"Wavelength","Divide") |
---|
180 | |
---|
181 | ####################### |
---|
182 | #Step 8 - Rescale(detector) |
---|
183 | # Enter the rescale value here |
---|
184 | rescale = 1.508*100.0 |
---|
185 | |
---|
186 | ####################### |
---|
187 | #Step 10 - Correct for sample/Can volume |
---|
188 | # Could easily move this to be done after the sample-can operation |
---|
189 | # Enter the value here |
---|
190 | thickness = 1.0 |
---|
191 | area = 3.1414956*8*8/4 |
---|
192 | |
---|
193 | correction = rescale/(thickness*area) |
---|
194 | |
---|
195 | CreateSingleValuedWorkspace("scalar",str(correction),"0.0") |
---|
196 | Multiply(outputWS_cor,"scalar",outputWS_cor) |
---|
197 | mantid.deleteWorkspace("scalar") |
---|
198 | |
---|
199 | ####################### |
---|
200 | #Step 11 - Convert to Q |
---|
201 | #Convert units to Q (MomentumTransfer) |
---|
202 | ConvertUnits(outputWS_cor,outputWS_cor,"MomentumTransfer") |
---|
203 | ConvertUnits(outputWS,outputWS,"MomentumTransfer") |
---|
204 | |
---|
205 | #Need to mark the workspace as a distribution at this point to get next rebin right |
---|
206 | ws = mantid.getMatrixWorkspace(outputWS_cor) |
---|
207 | ws.isDistribution(True) |
---|
208 | |
---|
209 | # Calculate the solid angle corrections |
---|
210 | SolidAngle(outputWS_cor,"solidangle") |
---|
211 | |
---|
212 | #rebin to desired Q bins |
---|
213 | q_bins = q1+","+dq+","+q2 |
---|
214 | Rebin(outputWS,outputWS,q_bins) |
---|
215 | Rebin(outputWS_cor,outputWS_cor,q_bins) |
---|
216 | RebinPreserveValue("solidangle","solidangle",q_bins) |
---|
217 | |
---|
218 | #Sum all spectra |
---|
219 | SumSpectra(outputWS,outputWS) |
---|
220 | SumSpectra(outputWS_cor,outputWS_cor) |
---|
221 | SumSpectra("solidangle","solidangle") |
---|
222 | |
---|
223 | ####################### |
---|
224 | #correct for solidangle |
---|
225 | Divide(outputWS_cor,"solidangle",outputWS_cor) |
---|
226 | #mantid.deleteWorkspace("solidangle") |
---|
227 | |
---|
228 | ####################### |
---|
229 | # Now put back the fractional error from the raw count workspace into the result |
---|
230 | PoissonErrors(outputWS_cor,outputWS,outputWS+wavname) |
---|
231 | mantid.deleteWorkspace(outputWS_cor) |
---|
232 | mantid.deleteWorkspace(outputWS) |
---|
233 | |
---|
234 | ####################### |
---|
235 | #step 12 - Cross section (remove can scatering) |
---|
236 | #Perform steps 1-11 for the sample and can |
---|
237 | # sample-can |
---|
238 | |
---|
239 | ####################### |
---|
240 | #step 12 - Save 1D data |
---|
241 | # this writes to the /bin directory - why ? |
---|
242 | SaveRKH(outputWS+wavname,Filename=outputWS+wavname+".Q",FirstColumnValue="MomentumTransfer") |
---|
243 | # print x |
---|
244 | print datetime.datetime.now() |
---|