Ticket #7: red_noCell_samples_pcNormalize.py

File red_noCell_samples_pcNormalize.py, 13.8 KB (added by wangzy, 6 years ago)
Line 
1# reduction without sensitivity and efficiency
2from mantid.simpleapi import *
3import math
4from xml.etree import ElementTree as ET
5import glob
6import re
7import os
8import sys
9
10def mergeWS(prepath, wsList,output):
11    for ix in wsList:
12        name='RUN'+str(ix).zfill(7)
13        LoadNexus(Filename=prepath+'\\'+name+'\\'+output+'.nxs',OutputWorkspace=str(ix))
14    for i in range(len(wsList)):
15        if len(wsList)== 1:
16            CloneWorkspace(InputWorkspace=str(wsList[i]), OutputWorkspace=output)
17            DeleteWorkspace(Workspace=str(wsList[i]))
18        elif wsList[i+1]!=wsList[-1]:
19            Plus(LHSWorkspace=str(wsList[i]),RHSWorkspace=str(wsList[i+1]), OutputWorkspace=str(wsList[i+1]))
20            DeleteWorkspace(Workspace=str(wsList[i]))
21        else:
22            Plus(LHSWorkspace=str(wsList[i]),RHSWorkspace=str(wsList[i+1]), OutputWorkspace=output)
23            DeleteWorkspace(Workspace=str(wsList[i]))
24            DeleteWorkspace(Workspace=str(wsList[i+1]))
25            break
26
27
28def getMask(xPids, yPids):
29    index=[]
30    xnum=xPids[1]-xPids[0]+1
31    ynum=yPids[1]-yPids[0]+1
32    sta=yPids[0]*125+xPids[0]
33    for i in range(xnum):
34        for j in range(ynum):
35            index.append(sta+j*125+i)
36    print index
37    return index
38
39
40def findCmaskWS(wsname):
41    x=0.0009
42    y=-0.00229
43    MoveInstrumentComponent(Workspace=wsname, ComponentName='Bank01', X=x, Y=y, RelativePosition=False)
44    index=getMask([0,18],[0,119])
45    MaskDetectors(Workspace=wsname, WorkspaceIndexList=index)
46    index=getMask([105,124],[0,119])
47    MaskDetectors(Workspace=wsname, WorkspaceIndexList=index)
48    inner=0.048
49    xml_inner='<infinite-cylinder id="beam_stop"><centre x="0" y="0" z="0.0" /><axis x="0" y="0" z="1" /><radius val= "'+str(inner)+'" /></infinite-cylinder><algebra val="beam_stop"/>'
50    MaskDetectorsInShape(Workspace=wsname, ShapeXML=xml_inner)
51
52    index=getMask([60,63],[0,57])
53    MaskDetectors(Workspace=wsname, WorkspaceIndexList=index)
54
55def calResolution(M2sample,wavelength, Q):
56    L1=11.90-7.45
57    #L1, source to sample distance, unit: m
58    L2=4.0
59    #L2, sample to detector distance, unit: m
60    d1=1.9/100.0
61    #d1,source aperture diameter, unit: m
62    d2=0.6/100.0
63    #d2, sample aperture diameter, unit: m
64    d3=1.6/100.0
65    #d3, detector res, unit: m   
66    detector_min=0.048
67    #detector minimum position
68    detector_max=0.7
69    #detector maximum position
70    delta_lamda=0.0057
71    #variance of lamda
72   
73    #Q=mtd['Q'].readX(0)
74   
75    #wavelength=mtd['wavelength'].readX(0)
76    efficiency=[]
77    for efficiency_index in range(len(wavelength)):
78        efficiency.append(4.7613/1000000.0+0.00054333*wavelength[efficiency_index])
79    #M2 efficiency
80
81   
82    Intensity=M2sample
83    Intensity_eff=[]
84    for Intensity_index in range(len(Intensity)):
85        Intensity_eff.append(Intensity[Intensity_index]/efficiency[Intensity_index])
86    #Intensity=Intensity/efficiency
87   
88    sigmaQ=[]
89    for Q_value in Q:
90        total_intensity=0
91        lamda_value_index=0
92        sigmaQ_sum=0
93        for lamda_value in wavelength:
94            Ldet=L2*Q_value*lamda_value/6.28
95           
96            if detector_min < Ldet < detector_max:
97                total_intensity=total_intensity+Intensity_eff[lamda_value_index]
98                sigmaQ_sum=sigmaQ_sum+Intensity_eff[lamda_value_index]*math.sqrt((delta_lamda)**2.0+0.25*((L2/L1)**2.0)*((d1/2.0/Ldet)**2.0)+0.25*((1+L2/L1)**2.0)*((d2/2.0/Ldet)**2.0)+1.0/12.0*((d3/Ldet)**2.0))*Q_value
99            lamda_value_index=lamda_value_index+1
100                     
101        sigmaQ.append(sigmaQ_sum/total_intensity)
102           
103       
104       
105    return sigmaQ
106       
107       
108def outputIq(prepath, runNo, transRunNo, qName, wName, m2Sample):
109
110    Q=mtd[qName].readX(0)
111    Q_shift=[]
112    for Q_shift_index in range(len(Q)-1):
113        Q_shift.append((Q[Q_shift_index]+Q[Q_shift_index+1])/2.0)
114       
115    iQ=mtd[qName].readY(0)
116   
117    sIQ=mtd[qName].readE(0)
118     
119    wavelength=mtd[wName].readX(0)
120    wavelength_shift=[]
121    for wavelength_shift_index in range(len(wavelength)-1):
122        wavelength_shift.append((wavelength[wavelength_shift_index]+wavelength[wavelength_shift_index+1])/2.0)
123    #change wavelength from bin to point
124   
125    M2sample=mtd[m2Sample].readY(0)
126     
127    sQ=calResolution(M2sample,wavelength_shift,Q)
128     
129    f_fp=prepath+'/SANS_'+str(runNo)+'.dat'     
130    ff=open(f_fp, "w")
131   
132    ff.write('#Run '+str(runNo)+' and Trans '+str(transRunNo)+' from SANS@CSNS\n')
133    #ff.write('%d',% (runNo))
134    ff.write('Q (1/A) | IQ (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ (1/A) \n')
135    #ff.write('%s', "RUN000001 from SANS@CSNS 2018----")
136     
137    for i in range(len(Q_shift)):
138        #strtmp='%.6f %.6f %.6f %.6f' % (Q[i],iQ[i],sIQ[i],sQ[i])
139        ff.write('%.6f %.6f %.6f %.6f \n' % (Q_shift[i],iQ[i],sIQ[i],sQ[i]))
140       
141       
142def readProtonCharge(prepath,wslist):           
143    pc=[]       
144    for ix in wslist:
145        runnum='RUN'+str(ix).zfill(7)
146        rawpath=prepath+'\\'+runnum+'\\'
147        xmllist=glob.glob(os.path.join(rawpath, '*.xml'))
148        xmlfile=xmllist[0]
149        root = ET.parse(xmlfile)
150        item=root.find('NXentry')   
151        pc.append(float(item.findtext('proton_charge')))
152    pc_final=sum(pc)/len(pc)
153
154    return pc_final   
155
156
157
158
159
160#Al_Mg_ht_1
161'''
162sample_run=[1095]
163transam_run=[1155]
164volume=3.14*0.34*0.34*0.09
165trans_run=[1154]
166Np=80
167
168
169trans_spline=True
170'''
171
172#JZFM_1
173'''
174sample_run=[1233]
175transam_run=[1230]
176volume=3.14*0.34*0.34*0.128
177trans_run=[1226]
178Np=80
179
180
181'''
182#J80
183'''
184sample_run=[1239]
185transam_run=[1173]
186volume=3.14*0.34*0.34*0.096
187trans_run=[1171]
188Np=80
189
190
191'''
192
193#J82
194'''
195sample_run=[1240]
196transam_run=[1176]
197volume=3.14*0.34*0.34*0.11
198trans_run=[1171]
199Np=80
200'''
201#J83
202
203sample_run=[1257]
204transam_run=[1177]
205volume=3.14*0.34*0.34*0.14
206trans_run=[1171]
207Np=80
208
209
210
211
212#BP_62.5K_4#
213'''
214sample_run=[1248]
215transam_run=[1247]
216volume=3.14*0.34*0.34*0.114
217trans_run=[1246]
218Np=80
219'''
220
221
222
223
224empty_run=[960,961]
225
226
227
228
229
230monSUM=True
231prepath=r'C:\Users\wzy\Documents\python\CSNSSANS'
232
233#readDx
234#empty
235'''
236mergeWS(prepath, empty_run, 'sample')
237CloneWorkspace(InputWorkspace='sample',OutputWorkspace='empty_tmp')
238LoadInstrument(Workspace='empty_tmp', Filename=prepath+'/paramData/detector_mod.xml', RewriteSpectraMap=True)
239#CloneWorkspace(InputWorkspace='bg',OutputWorkspace='empty_tmp')
240findCmaskWS('empty_tmp')
241ConvertUnits(InputWorkspace='empty_tmp', OutputWorkspace='empty_tmp0', Target='Wavelength', AlignBins=True)
242
243
244#sample
245
246mergeWS(prepath, sample_run, 'sample')
247CloneWorkspace(InputWorkspace='sample',OutputWorkspace='scattering')
248CloneWorkspace(InputWorkspace='scattering',OutputWorkspace='sam_tmp')
249
250LoadInstrument(Workspace='sam_tmp', Filename=prepath+'/paramData/detector_mod.xml', RewriteSpectraMap=True)
251findCmaskWS('sam_tmp')
252ConvertUnits(InputWorkspace='sam_tmp', OutputWorkspace='sam_tmp0', Target='Wavelength', AlignBins=True)
253'''
254
255
256
257Np=80
258wave_rebin='4.6,0.02,7.9'
259wave_rebin_all='2,0.02,7.9'
260Np=Np/4
261
262
263
264
265CreateSingleValuedWorkspace(OutputWorkspace='v', DataValue=volume)
266#calculate sample transmission
267mergeWS(prepath, transam_run, 'monitors')
268RenameWorkspace(InputWorkspace='monitors',OutputWorkspace='monTranSam')
269mergeWS(prepath, trans_run,'monitors')
270RenameWorkspace(InputWorkspace='monitors',OutputWorkspace='monTran')
271
272
273
274LoadInstrument(Workspace='monTran', Filename=prepath+'/paramData/monitors.xml', RewriteSpectraMap=True)
275LoadInstrument(Workspace='monTranSam', Filename=prepath+'/paramData/monitors.xml', RewriteSpectraMap=True)
276
277
278
279ConvertUnits(InputWorkspace='monTranSam', OutputWorkspace='monTranSam', Target='Wavelength', AlignBins=True)
280ConvertUnits(InputWorkspace='monTran', OutputWorkspace='monTran', Target='Wavelength', AlignBins=True)
281
282Rebin(InputWorkspace='monTranSam',OutputWorkspace='monTranSam',Params=wave_rebin)
283name=mtd['monTranSam']
284sample_trans_m3=(name.readY(1))*1.0
285
286Rebin(InputWorkspace='monTran',OutputWorkspace='monTran',Params=wave_rebin)
287name=mtd['monTran']
288empty_m3=(name.readY(1))*1.0
289
290wavelength_trans_M3_sample_to_empty=(name.readX(0))*1.0
291
292pc_sample=readProtonCharge(prepath, sample_run)
293pc_sample_trans=readProtonCharge(prepath, transam_run)
294pc_empty=readProtonCharge(prepath, trans_run)
295
296trans_M3_sample_to_empty=sample_trans_m3/empty_m3
297
298CreateWorkspace(OutputWorkspace='trans_M3_ratio_unfitted', DataX=wavelength_trans_M3_sample_to_empty, DataY=trans_M3_sample_to_empty)
299
300Fit("name=Polynomial,n=1",InputWorkspace='trans_M3_ratio_unfitted',Output='trans_M3_ratio')
301
302name=mtd['trans_M3_ratio_Workspace']
303trans_M3_ratio=(name.readY(1))*1.0
304
305CreateWorkspace(OutputWorkspace='trans_M3_ratio', DataX=wavelength_trans_M3_sample_to_empty, DataY=trans_M3_ratio)
306
307
308
309
310scaling_pc_trans=pc_sample_trans/pc_empty;
311#scaling_from_m2Counts_to_pc_trans=1;
312
313CreateSingleValuedWorkspace(OutputWorkspace='scaling_trans_pc', DataValue=scaling_pc_trans)
314
315Divide(LHSWorkspace='trans_M3_ratio', RHSWorkspace='scaling_trans_pc', OutputWorkspace='trans_sam', AllowDifferentNumberSpectra=True)
316Rebin(InputWorkspace='trans_sam', OutputWorkspace="trans_sam", Params=wave_rebin)
317
318
319dlambda=0.02
320CreateSingleValuedWorkspace(OutputWorkspace='dlambda', DataValue=dlambda)
321Multiply(LHSWorkspace='trans_sam', RHSWorkspace='dlambda', OutputWorkspace='trans_sam', AllowDifferentNumberSpectra=True)
322ConvertToDistribution(Workspace='trans_sam')
323
324Multiply(LHSWorkspace='trans_M3_ratio_unfitted', RHSWorkspace='dlambda', OutputWorkspace='trans_sam_unfitted', AllowDifferentNumberSpectra=True)
325Divide(LHSWorkspace='trans_sam_unfitted', RHSWorkspace='scaling_trans_pc', OutputWorkspace='trans_sam_unfitted', AllowDifferentNumberSpectra=True)
326Rebin(InputWorkspace='trans_sam_unfitted', OutputWorkspace='trans_sam_unfitted', Params=wave_rebin)
327ConvertToDistribution(Workspace='trans_sam_unfitted')
328
329
330
331
332DeleteWorkspace(Workspace='monTran')
333
334'''
335#load calibration
336#==========================================
337LoadAscii(Filename=prepath+'/paramData/cal_bp_1M.dat', OutputWorkspace='cali', Separator='Automatic', Unit='Wavelength')
338
339ConvertToHistogram(InputWorkspace='cali', OutputWorkspace='cali')
340
341Rebin(InputWorkspace='cali', OutputWorkspace='cali', Params=wave_rebin)
342
343
344#empty
345
346
347Rebin(InputWorkspace='empty_tmp0',OutputWorkspace='empty_tmp',Params=wave_rebin)
348
349mergeWS(prepath, empty_run,'monitor2')
350RenameWorkspace(InputWorkspace='monitor2',OutputWorkspace='mon_empty')
351LoadInstrument(Workspace='mon_empty', Filename=prepath+'/paramData/monitor2_based_on_monitors.xml', RewriteSpectraMap=True)
352
353ConvertUnits(InputWorkspace='mon_empty', OutputWorkspace='empty_m2', Target='Wavelength', AlignBins=True)
354Rebin(InputWorkspace='empty_m2',OutputWorkspace='empty_m2',Params=wave_rebin_all)
355
356
357if monSUM:
358    name=mtd['empty_m2']
359    counts=sum(name.readY(0))*1.0
360    CreateSingleValuedWorkspace(OutputWorkspace='m2',DataValue=counts)
361    Divide(LHSWorkspace='empty_tmp', RHSWorkspace='m2', OutputWorkspace='empty_tmp', AllowDifferentNumberSpectra=True)
362else:
363    Divide(LHSWorkspace='empty_tmp', RHSWorkspace='empty_m2', OutputWorkspace='empty_tmp', AllowDifferentNumberSpectra=True)
364
365Divide(LHSWorkspace='empty_tmp', RHSWorkspace='cali', OutputWorkspace='empty_tmp', AllowDifferentNumberSpectra=True)
366
367
368
369# sample and emptycell
370#========================================
371#sample
372
373
374Rebin(InputWorkspace='sam_tmp0',OutputWorkspace='sam_tmp',Params=wave_rebin)
375
376
377mergeWS(prepath, sample_run,'monitor2')
378RenameWorkspace(InputWorkspace='monitor2',OutputWorkspace='sam_m2')
379LoadInstrument(Workspace='sam_m2', Filename=prepath+'/paramData/monitor2_based_on_monitors.xml', RewriteSpectraMap=True)
380
381ConvertUnits(InputWorkspace='sam_m2', OutputWorkspace='sam_m2', Target='Wavelength', AlignBins=True)
382Rebin(InputWorkspace='sam_m2',OutputWorkspace='sam_m2',Params=wave_rebin_all)
383
384#
385if monSUM:
386    name=mtd['sam_m2']
387    counts=sum(name.readY(0))*1.0
388    CreateSingleValuedWorkspace(OutputWorkspace='m2',DataValue=counts)
389    Divide(LHSWorkspace='sam_tmp', RHSWorkspace='m2', OutputWorkspace='sam_tmp', AllowDifferentNumberSpectra=True)
390    CloneWorkspace(InputWorkspace='sam_tmp',OutputWorkspace='sam_raw')
391else:
392    Divide(LHSWorkspace='sam_tmp', RHSWorkspace='sam_m2', OutputWorkspace='sam_tmp', AllowDifferentNumberSpectra=True)
393
394
395
396scaling_pc_sample_scat=pc_sample/pc_empty;
397#scaling_from_m2Counts_to_pc_sample_scat=1;
398
399CreateSingleValuedWorkspace(OutputWorkspace='scaling_scat_sample', DataValue=scaling_pc_sample_scat)
400
401Multiply(LHSWorkspace='sam_tmp', RHSWorkspace='scaling_scat_sample', OutputWorkspace='sam_tmp', AllowDifferentNumberSpectra=True)
402
403Divide(LHSWorkspace='sam_tmp', RHSWorkspace='cali', OutputWorkspace='sam_tmp', AllowDifferentNumberSpectra=True)
404
405ApplyTransmissionCorrection(InputWorkspace='sam_tmp',TransmissionWorkspace='trans_sam',OutputWorkspace='sam_tmp',ThetaDependent=True)
406
407q_rebin='0.01,-0.08,0.4'
408
409Q1D(DetBankWorkspace='empty_tmp', OutputWorkspace='IQ_air', OutputBinning=q_rebin, AccountForGravity=False, SolidAngleWeighting=True)
410ReplaceSpecialValues(InputWorkspace='IQ_air', OutputWorkspace='IQ_air', NaNValue=0, InfinityValue=0)
411
412Q1D(DetBankWorkspace='sam_tmp', OutputWorkspace='IQ_sam_raw', OutputBinning=q_rebin, AccountForGravity=False, SolidAngleWeighting=True)
413ReplaceSpecialValues(InputWorkspace='IQ_sam_raw', OutputWorkspace='IQ_sam_raw', NaNValue=0, InfinityValue=0)
414
415Minus(LHSWorkspace='IQ_sam_raw', RHSWorkspace='IQ_air',OutputWorkspace='IQ_v')
416
417Divide(LHSWorkspace='IQ_v', RHSWorkspace='v', OutputWorkspace='IQ_v', AllowDifferentNumberSpectra=True)
418
419dlamda=0.5
420CreateSingleValuedWorkspace(OutputWorkspace='dlamda_sample', DataValue=dlamda)
421Multiply(LHSWorkspace='IQ_v', RHSWorkspace='dlamda_sample', OutputWorkspace='IQ_v', AllowDifferentNumberSpectra=True)
422'''