Context Navigation

VC_DetectorBinning_Slit.ipf @ 1077

Last change on this file since 1077 was 1077, checked in by srkline, 5 years ago
more attempts to correct the slit binning - to get the absolute scaling correct - still in progress.
File size: 12.0 KB

Line
1	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
2
3
4	//////////////////
5	//
6	// Procedures to average data taken in "slit" aperture geometry. The fake data on the detector panels
7	// is as would be collected in PINHOLE geometry - I do not currently have a simulation for slit
8	// apertures (this would need to be MonteCarlo) - the I(q) averaging here gives the I(q) that you
9	// would measure in 1D using slit geometry. It is done by simply summing the columns of data on each detector.
10	//
11	//////////////////
12	//
13	//
14	// TODO - big question about averaging in this way...
15	// can the T/B panels really be used at all for slit mode? - since there's a big "hole" in the scattering data
16	// collected -- you're not getting the full column of data covering a wide range of Qy. L/R panels should be fine,
17	// although the different carriages cover very different "chunks" of Qy
18	//
19	// - best answer so far is to skip the T/B panels, and simply not use them
20	//
21	// TODO -- be sure that the absolute scaling of this is correct. Right now, something is off.
22	// -- so write up the documentaiton for this operation, and try to find the error in the process...
23	//
24	//
25	/////////////////
26
27
28	// TODO:
29	// -- verify the error calculation
30	// -- ? add in functionality to handle FLR and MLR cases (2 panels of data)
31	//
32	// seems backwards to call this "byRows", but this is the way that Igor indexes
33	// LR banks are defined as (48,256) (n,m), sumRows gives sum w/ dimension (n x 1)
34	//
35	// updated to new folder structure Feb 2016
36	// folderStr = RAW,SAM, VCALC or other
37	// detStr is the panel identifer "ML", etc.
38	Function VC_fBinDetector_byRows(folderStr,detStr)
39	String folderStr,detStr
40
41	// SetDataFolder root:Packages:NIST:VSANS:VCALC
42
43	Variable pixSizeX,pixSizeY,delQy
44	Variable isVCALC=0,maskMissing,nsets=0
45
46	maskMissing = 1 // set to zero if a mask is actually present
47
48	if(cmpstr(folderStr,"VCALC") == 0)
49	isVCALC = 1
50	endif
51
52	String folderPath = "root:Packages:NIST:VSANS:"+folderStr
53	String instPath = ":entry:instrument:detector_"
54
55	strswitch(detStr) // string switch
56	// only one panel, simply pick that panel and move on out of the switch
57	case "FL":
58	case "FR":
59	case "ML":
60	case "MR":
61	case "B":
62	if(isVCALC)
63	WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) // 2D detector data
64	WAVE/Z iErr = $("asdf_iErr_"+detStr) // TODO: 2D errors -- may not exist, especially for simulation
65	else
66	Wave inten = V_getDetectorDataW(folderStr,detStr)
67	Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
68	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
69	if(WaveExists(mask) == 1)
70	maskMissing = 0
71	endif
72	endif
73
74	nsets = 1
75	break
76
77	// case "FLR":
78	// // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
79	// if(isVCALC)
80	// WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
81	// WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
82	// WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
83	// WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
84	// else
85	// Wave inten = V_getDetectorDataW(folderStr,"FL")
86	// Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
87	// Wave inten2 = V_getDetectorDataW(folderStr,"FR")
88	// Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
89	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
90	// Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
91	// if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
92	// maskMissing = 0
93	// endif
94	// endif
95	//
96	// nsets = 2
97	// break
98	//
99	// case "MLR":
100	// if(isVCALC)
101	// WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
102	// WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
103	// WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
104	// WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
105	// else
106	// Wave inten = V_getDetectorDataW(folderStr,"ML")
107	// Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
108	// Wave inten2 = V_getDetectorDataW(folderStr,"MR")
109	// Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
110	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
111	// Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
112	// if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
113	// maskMissing = 0
114	// endif
115	// endif
116	//
117	// nSets = 2
118	// break
119
120	default:
121	nSets = 0
122	Print "ERROR ---- type is not recognized "
123	endswitch
124
125	if(nSets == 0)
126	SetDataFolder root:
127	return(0)
128	endif
129
130
131	Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
132	Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
133	Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
134
135
136	// delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution
137	// delQx = abs(qx[0][0] - qx[1][0])
138
139	Variable nq,val
140	nq = DimSize(inten,0) //nq == the number of columns (x dimension)
141
142	// SetDataFolder $(folderPath+instPath+detStr)
143
144	Make/O/D/N=(nq) $(folderPath+":"+"iBin_qxqy_"+detStr)
145	Make/O/D/N=(nq) $(folderPath+":"+"qBin_qxqy_"+detStr)
146	Make/O/D/N=(nq) $(folderPath+":"+"nBin_qxqy_"+detStr)
147	Make/O/D/N=(nq) $(folderPath+":"+"iBin2_qxqy_"+detStr)
148	Make/O/D/N=(nq) $(folderPath+":"+"eBin_qxqy_"+detStr)
149	Make/O/D/N=(nq) $(folderPath+":"+"eBin2D_qxqy_"+detStr)
150
151	Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+detStr)
152	Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy_"+detStr)
153	Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy_"+detStr)
154	Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy_"+detStr)
155	Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy_"+detStr)
156	Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy_"+detStr)
157
158
159	iBin_qxqy = 0
160	iBin2_qxqy = 0
161	eBin_qxqy = 0
162	eBin2D_qxqy = 0
163	nBin_qxqy = 0
164
165
166	// sum the rows
167
168	// MatrixOp would be fast, but I can't figure out how to apply the mask with sumRows???
169	// MatrixOp/O iBin_qxqy = sumRows(inten) //automatically generates the destination
170	//
171	//// how to properly calculate the error?
172	//// This MatrixOp gives values way too large -- larger than the intensity (be sure to correct *delQy below...
173	//
174	// MatrixOp/O tmp = sqrt(varCols(inten^t)) // variance: no varRows operation, so use the transpose of the matrix
175	// eBin_qxqy = tmp[0][p]
176
177	Variable ii,jj,ntube,npix,sum_inten, sum_n, sum_inten2,avesq,aveisq,var,mask_val
178	ntube = DimSize(inten,0)
179	npix = DimSize(inten,1)
180
181	for(ii=0;ii<ntube;ii+=1)
182	sum_inten = 0 // initialize the sum
183	sum_n = 0
184	sum_inten2 = 0
185
186	for(jj=0;jj<npix;jj+=1)
187	val = inten[ii][jj]
188	if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
189	mask_val = 0
190	else
191	mask_val = mask[ii][jj]
192	endif
193	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
194	sum_inten += val
195	sum_n += 1
196	sum_inten2 += val*val
197	endif
198	endfor
199	iBin_qxqy[ii] = sum_inten/sum_n //the average value
200
201	avesq = sum_inten
202	aveisq = sum_inten2/sum_n
203	var = aveisq-avesq
204	if(var<=0)
205	eBin_qxqy[ii] = 1e-6
206	else
207	eBin_qxqy[ii] = sqrt(var/(sum_n - 1))
208	endif
209
210	endfor
211
212	//TODO: use only the Qx component in the y-center of the detector, not Qtotal
213	qBin_qxqy = abs(qx[p][npix/2])
214
215	// if the detectors are "L", then the values are all negative...
216	qBin_qxqy = abs(qBin_qxqy)
217
218	// for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed
219	// won't hurt to sort the R data
220	Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy
221
222	// not needed to corrrect data, but put in the file???
223	delQy = abs(qy[0][npix-1] - qy[0][0]) // TODO: do I use dQ for the height of the panel?
224	// iBin_qxqy *= delQy
225	// eBin_qxqy *= delQy
226
227
228	/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty"
229	// clean out the near-zero Q point in the T/B and Back detectors
230	// qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]
231
232
233	// clear out zero data values before exiting...
234	// find the last non-zero point, working backwards
235	val = numpnts(iBin_qxqy)
236	do
237	val -= 1
238	while((iBin_qxqy[val] == 0) && val > 0)
239
240	DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy
241
242	// work forwards? this doesn't work...
243	// val = -1
244	// do
245	// val += 1
246	// while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)
247	// DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy
248
249
250	// TODO:
251	// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS
252	// since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian
253
254	// TODO:
255	// -- This is where I calculate the resolution in SANS (see CircSectAve)
256	// -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
257	// -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging
258	//
259
260	nq = numpnts(qBin_qxqy)
261	Make/O/D/N=(nq) $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr)
262	Make/O/D/N=(nq) $(folderPath+":"+"qBar_qxqy"+"_"+detStr)
263	Make/O/D/N=(nq) $(folderPath+":"+"fSubS_qxqy"+"_"+detStr)
264	Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr)
265	Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr)
266	Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr)
267
268	// TODO
269	// -- these are DUMMY VALUES!!!
270	sigmaq = -delQy
271	qbar = -delQy
272	fsubs = -delQy
273
274
275	SetDataFolder root:
276
277	return(0)
278	End
279
280
281	// TODO -- update to new folder structure
282	// unused -- update if necessary
283	Proc CopyIQWaves()
284
285	SetDataFolder root:Packages:NIST:VSANS:VCALC
286
287	duplicate/O iBin_qxqy_B iBin_qxqy_B_pin
288	duplicate/O qBin_qxqy_B qBin_qxqy_B_pin
289
290	duplicate/O iBin_qxqy_MR iBin_qxqy_MR_pin
291	duplicate/O qBin_qxqy_MR qBin_qxqy_MR_pin
292	duplicate/O iBin_qxqy_MT iBin_qxqy_MT_pin
293	duplicate/O qBin_qxqy_MT qBin_qxqy_MT_pin
294	duplicate/O iBin_qxqy_ML iBin_qxqy_ML_pin
295	duplicate/O qBin_qxqy_ML qBin_qxqy_ML_pin
296	duplicate/O iBin_qxqy_MB iBin_qxqy_MB_pin
297	duplicate/O qBin_qxqy_MB qBin_qxqy_MB_pin
298
299	duplicate/O iBin_qxqy_FR iBin_qxqy_FR_pin
300	duplicate/O qBin_qxqy_FR qBin_qxqy_FR_pin
301	duplicate/O iBin_qxqy_FT iBin_qxqy_FT_pin
302	duplicate/O qBin_qxqy_FT qBin_qxqy_FT_pin
303	duplicate/O iBin_qxqy_FL iBin_qxqy_FL_pin
304	duplicate/O qBin_qxqy_FL qBin_qxqy_FL_pin
305	duplicate/O iBin_qxqy_FB iBin_qxqy_FB_pin
306	duplicate/O qBin_qxqy_FB qBin_qxqy_FB_pin
307
308	SetDataFolder root:
309	End
310
311	// TODO -- update to new folder structure
312	// unused -- update if necessary
313	Window slit_vs_pin_graph() : Graph
314	PauseUpdate; Silent 1 // building window...
315	String fldrSav0= GetDataFolder(1)
316	SetDataFolder root:Packages:NIST:VSANS:VCALC:
317	Display /W=(1296,44,1976,696) iBin_qxqy_B vs qBin_qxqy_B
318	AppendToGraph iBin_qxqy_B_pin vs qBin_qxqy_B_pin
319	AppendToGraph iBin_qxqy_MR vs qBin_qxqy_MR
320	AppendToGraph iBin_qxqy_MR_pin vs qBin_qxqy_MR_pin
321	AppendToGraph iBin_qxqy_MT vs qBin_qxqy_MT
322	AppendToGraph iBin_qxqy_MT_pin vs qBin_qxqy_MT_pin
323	AppendToGraph iBin_qxqy_FR vs qBin_qxqy_FR
324	AppendToGraph iBin_qxqy_FR_pin vs qBin_qxqy_FR_pin
325	AppendToGraph iBin_qxqy_FT vs qBin_qxqy_FT
326	AppendToGraph iBin_qxqy_FT_pin vs qBin_qxqy_FT_pin
327	SetDataFolder fldrSav0
328	ModifyGraph mode=4
329	ModifyGraph marker=19
330	ModifyGraph lSize=2
331	ModifyGraph rgb(iBin_qxqy_B)=(0,0,0),rgb(iBin_qxqy_B_pin)=(65535,16385,16385),rgb(iBin_qxqy_MR)=(2,39321,1)
332	ModifyGraph rgb(iBin_qxqy_MR_pin)=(0,0,65535),rgb(iBin_qxqy_MT)=(39321,1,31457)
333	ModifyGraph rgb(iBin_qxqy_MT_pin)=(48059,48059,48059),rgb(iBin_qxqy_FR)=(65535,32768,32768)
334	ModifyGraph rgb(iBin_qxqy_FR_pin)=(0,65535,0),rgb(iBin_qxqy_FT)=(16385,65535,65535)
335	ModifyGraph rgb(iBin_qxqy_FT_pin)=(65535,32768,58981)
336	ModifyGraph msize=2
337	ModifyGraph grid=1
338	ModifyGraph log=1
339	ModifyGraph mirror=1
340	SetAxis bottom 1e-05,0.05287132
341	Legend/C/N=text0/J/X=64.73/Y=7.25 "\\Z12\\s(iBin_qxqy_B) iBin_qxqy_B\r\\s(iBin_qxqy_B_pin) iBin_qxqy_B_pin\r\\s(iBin_qxqy_MR) iBin_qxqy_MR"
342	AppendText "\\s(iBin_qxqy_MR_pin) iBin_qxqy_MR_pin\r\\s(iBin_qxqy_MT) iBin_qxqy_MT\r\\s(iBin_qxqy_MT_pin) iBin_qxqy_MT_pin\r\\s(iBin_qxqy_FR) iBin_qxqy_FR"
343	AppendText "\\s(iBin_qxqy_FR_pin) iBin_qxqy_FR_pin\r\\s(iBin_qxqy_FT) iBin_qxqy_FT\r\\s(iBin_qxqy_FT_pin) iBin_qxqy_FT_pin"
344	EndMacro

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format