Ticket #773: LOQ_analysis_new.py

File LOQ_analysis_new.py, 10.9 KB (added by Russell Taylor, 11 years ago)
Line 
1import datetime
2#LOQ data analysis script
3#path = "C:/MantidInstall/"
4#path = "C:/Mantid/Test/"
5path="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
8def 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#######################
18print datetime.datetime.now()
19for 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
244print datetime.datetime.now()