Context Navigation

V_IQ_Annular.ipf @ 1073

Last change on this file since 1073 was 1073, checked in by srkline, 5 years ago

lots of changes here:
many little fixes to clean up TODO items and marke them DONE

changed the handling of the panel "gap" to split the gap evenly. Q-calculations have been re-verified with this change.

re-named the list of "bin Type" values, and added a few more choices. Streamlined how the averaging and plotting works with this list so that it can be more easily modified as different combinations of binning are envisioned. This resulted in a lot of excess code being cut out and replaced with cleaner logic. This change has also been verified to work as intended.

Attenuation is now always calculated from the table. The table also by (NEW) definition has values for the white beam (one waelength) and graphite (multiple possible wavelengths) where the wavelengths are artificially scaled (*1000) or *1e6) so that the interpolations can be done internally without the need for multiple attenuator tables.

File size: 23.8 KB

Line
1	#pragma TextEncoding = "MacRoman"
2	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
3
4
5	// Procedures to do an annular binning of the data
6	//
7	// As for SANS, needs a Q-center and Q-delta to define the annular ring,
8	// and the number of bins to divide the 360 degree circle
9	//
10	// qWidth is +/- around the q-center
11	//
12	//
13	// TODO
14	// x- add error bars to the plot of phi
15	// x- data writer to export annular data (3-column)
16	// x- loader to re-read annular data (will the normal loader work?) -- yes
17	// -- integrate this with the protocol
18	// -- integrate this with a more general "average panel"
19	// -- draw the q-center and width on the image (as a contour? - if no, use a draw layer...)
20	// Can I "flag" the pixels that contribute to the annualr average, and overlay them like a
21	// Mask(translucent?), like a "thresholding" operation, but act on the Q-value, not the Z- value.
22
23	//
24	//
25	Proc Annular_Binning(folderStr,detGroup,qCtr_Ann,qWidth)
26	String folderStr="SAM",detGroup="F"
27	Variable qCtr_Ann=0.1,qWidth=0.01 // +/- in A^-1
28
29	V_QBinAllPanels_Annular(folderStr,detGroup,qCtr_Ann,qWidth)
30
31	V_Phi_Graph_Proc(folderStr,detGroup)
32
33	End
34
35
36	// TODO -- binType == 4 (slit mode) should never end up here, as it makes no sense
37	//
38	// -- really, the only binning that makes any sense is "one", treating each panel individually,
39	// so I may scrap the parameter, or ignore it. so don't count on it in the future.
40	//
41	Function V_QBinAllPanels_Annular(folderStr,detGroup,qCtr_Ann,qWidth)
42	String folderStr,detGroup
43	Variable qCtr_Ann,qWidth // +/- in A^-1
44
45	Variable ii,delQ
46	String detStr
47
48	if(cmpstr(detGroup,"F") == 0)
49	detStr = "FLRTB"
50	else
51	detStr = "MLRTB"
52	endif
53
54	// right now, use all of the detectors. There is a lot of waste in this and it could be
55	// done a lot faster, but...
56
57	// TODO
58	// detStr = "FLRTB" or "MLRTB", depending which panel the q-ring is centered on/
59	// for now, no crossing of the rings onto different carriages
60
61	V_fDoAnnularBin_QxQy2D(folderStr,detStr,qCtr_Ann,qWidth)
62
63
64	return(0)
65	End
66
67	Proc V_Phi_Graph_Proc(folderStr,detGroup)
68	String folderStr,detGroup
69
70	PauseUpdate; Silent 1 // building window...
71	String fldrSav0= GetDataFolder(1)
72	SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
73
74	DoWindow/F V_Phi_Graph
75	if(V_flag == 0)
76	Display /W=(35,45,572,419)/N=V_Phi_Graph /K=1
77	else
78	RemoveFromGraph/Z iPhiBin_qxqy_FLRTB
79	RemoveFromGraph/Z iPhiBin_qxqy_MLRTB
80	endif
81
82	if(cmpstr(detGroup,"F") == 0)
83	AppendToGraph iPhiBin_qxqy_FLRTB vs phiBin_qxqy_FLRTB
84	// Display /W=(35,45,572,419)/N=V_Phi_Graph /K=1 iPhiBin_qxqy_FLRTB vs phiBin_qxqy_FLRTB
85	ModifyGraph mode=4
86	ModifyGraph marker=19
87	ErrorBars iPhiBin_qxqy_FLRTB Y,wave=(ePhiBin_qxqy_FLRTB,ePhiBin_qxqy_FLRTB)
88	Label left "Counts"
89	Label bottom "Phi (degrees)"
90	Legend
91	else
92	AppendToGraph iPhiBin_qxqy_MLRTB vs phiBin_qxqy_MLRTB
93	ModifyGraph mode=4
94	ModifyGraph marker=19
95	ErrorBars iPhiBin_qxqy_MLRTB Y,wave=(ePhiBin_qxqy_MLRTB,ePhiBin_qxqy_MLRTB)
96	Label left "Counts"
97	Label bottom "Phi (degrees)"
98	Legend
99	endif
100
101
102	SetDataFolder fldrSav0
103	EndMacro
104
105
106	//////////
107	//
108	//
109	//
110	// TODO
111	// -- "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
112	//
113	// ** Currently, type is being passed in as "" and ignored (looping through all of the detector panels
114	// to potentially add to the annular bins)
115	//
116	// folderStr = WORK folder, type = the binning type (may include multiple detectors)
117	//
118	// Variable qCtr_Ann = 0.1
119	// Variable qWidth = 0.02 // +/- in A^-1
120	//
121	Function V_fDoAnnularBin_QxQy2D(folderStr,type,qCtr_Ann,qWidth)
122	String folderStr,type
123	Variable qCtr_Ann,qWidth
124
125	Variable nSets = 0
126	Variable xDim,yDim
127	Variable ii,jj
128	Variable qVal,nq,var,avesq,aveisq
129	Variable binIndex,val,isVCALC=0,maskMissing
130
131	String folderPath = "root:Packages:NIST:VSANS:"+folderStr
132	String instPath = ":entry:instrument:detector_"
133	String detStr
134
135	if(cmpstr(folderStr,"VCALC") == 0)
136	isVCALC = 1
137	endif
138
139	// now switch on the type to determine which waves to declare and create
140	// since there may be more than one panel to step through. There may be two, there may be four
141	//
142
143	// TODO:
144	// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
145	//
146	// assume that the mask files are missing unless we can find them. If VCALC data,
147	// then the Mask is missing by definition
148	maskMissing = 1
149
150	strswitch(type) // string switch
151
152
153	case "FLRTB":
154	if(isVCALC)
155	WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
156	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
157	WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
158	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
159	WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
160	WAVE/Z iErr3 = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
161	WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
162	WAVE/Z iErr4 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
163	else
164	Wave inten = V_getDetectorDataW(folderStr,"FL")
165	Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
166	Wave inten2 = V_getDetectorDataW(folderStr,"FR")
167	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
168	Wave inten3 = V_getDetectorDataW(folderStr,"FT")
169	Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
170	Wave inten4 = V_getDetectorDataW(folderStr,"FB")
171	Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
172	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
173	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
174	Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
175	Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
176	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
177	maskMissing = 0
178	endif
179	endif
180	// NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
181
182	Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values
183	Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values
184	Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values
185	Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values
186
187	Wave qx = $(folderPath+instPath+"FL"+":qx_"+"FL") // 2D qx-values
188	Wave qx2 = $(folderPath+instPath+"FR"+":qx_"+"FR") // 2D qx-values
189	Wave qx3 = $(folderPath+instPath+"FT"+":qx_"+"FT") // 2D qx-values
190	Wave qx4 = $(folderPath+instPath+"FB"+":qx_"+"FB") // 2D qx-values
191
192	Wave qy = $(folderPath+instPath+"FL"+":qy_"+"FL") // 2D qy-values
193	Wave qy2 = $(folderPath+instPath+"FR"+":qy_"+"FR") // 2D qy-values
194	Wave qy3 = $(folderPath+instPath+"FT"+":qy_"+"FT") // 2D qy-values
195	Wave qy4 = $(folderPath+instPath+"FB"+":qy_"+"FB") // 2D qy-values
196
197	nSets = 4
198	break
199
200
201
202	case "MLRTB":
203	if(isVCALC)
204	WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
205	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
206	WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
207	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
208	WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
209	WAVE/Z iErr3 = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
210	WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
211	WAVE/Z iErr4 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
212	else
213	Wave inten = V_getDetectorDataW(folderStr,"ML")
214	Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
215	Wave inten2 = V_getDetectorDataW(folderStr,"MR")
216	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
217	Wave inten3 = V_getDetectorDataW(folderStr,"MT")
218	Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
219	Wave inten4 = V_getDetectorDataW(folderStr,"MB")
220	Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
221	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
222	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
223	Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
224	Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
225	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
226	maskMissing = 0
227	endif
228	endif
229	// NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
230
231	Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values
232	Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values
233	Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values
234	Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values
235
236	Wave qx = $(folderPath+instPath+"ML"+":qx_"+"ML") // 2D qx-values
237	Wave qx2 = $(folderPath+instPath+"MR"+":qx_"+"MR") // 2D qx-values
238	Wave qx3 = $(folderPath+instPath+"MT"+":qx_"+"MT") // 2D qx-values
239	Wave qx4 = $(folderPath+instPath+"MB"+":qx_"+"MB") // 2D qx-values
240
241	Wave qy = $(folderPath+instPath+"ML"+":qy_"+"ML") // 2D qy-values
242	Wave qy2 = $(folderPath+instPath+"MR"+":qy_"+"MR") // 2D qy-values
243	Wave qy3 = $(folderPath+instPath+"MT"+":qy_"+"MT") // 2D qy-values
244	Wave qy4 = $(folderPath+instPath+"MB"+":qy_"+"MB") // 2D qy-values
245
246	nSets = 4
247	break
248
249	default:
250	nSets = 0 // optional default expression executed
251	Print "ERROR ---- type is not recognized "
252	endswitch
253
254	// Print "delQ = ",delQ," for ",type
255
256	if(nSets == 0)
257	SetDataFolder root:
258	return(0)
259	endif
260
261
262	//TODO: properly define the 2D errors here - I'll have this if I do the simulation
263	// -- need to propagate the 2D errors up to this point
264	//
265	if(WaveExists(iErr)==0 && WaveExists(inten) != 0)
266	Duplicate/O inten,iErr
267	Wave iErr=iErr
268	// iErr = 1+sqrt(inten+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
269	iErr = sqrt(inten+0.75) // TODO -- here I'm just using some fictional value
270	endif
271	if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
272	Duplicate/O inten2,iErr2
273	Wave iErr2=iErr2
274	// iErr2 = 1+sqrt(inten2+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
275	iErr2 = sqrt(inten2+0.75) // TODO -- here I'm just using some fictional value
276	endif
277	if(WaveExists(iErr3)==0 && WaveExists(inten3) != 0)
278	Duplicate/O inten3,iErr3
279	Wave iErr3=iErr3
280	// iErr3 = 1+sqrt(inten3+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
281	iErr3 = sqrt(inten3+0.75) // TODO -- here I'm just using some fictional value
282	endif
283	if(WaveExists(iErr4)==0 && WaveExists(inten4) != 0)
284	Duplicate/O inten4,iErr4
285	Wave iErr4=iErr4
286	// iErr4 = 1+sqrt(inten4+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
287	iErr4 = sqrt(inten4+0.75) // TODO -- here I'm just using some fictional value
288	endif
289
290	// TODO -- nq will need to be larger, once the back detector is installed
291	//
292	// note that the back panel of 320x320 (1mm res) results in 447 data points!
293	// - so I upped nq to 600
294
295	nq = 600
296
297	//****TODO**** -- where to put the averaged data -- right now, folderStr is forced to ""
298	// SetDataFolder $("root:"+folderStr) //should already be here, but make sure...
299	Make/O/D/N=(nq) $(folderPath+":"+"iPhiBin_qxqy"+"_"+type)
300	// Make/O/D/N=(nq) $(folderPath+":"+"qBin_qxqy"+"_"+type)
301	Make/O/D/N=(nq) $(folderPath+":"+"phiBin_qxqy"+"_"+type)
302	Make/O/D/N=(nq) $(folderPath+":"+"nPhiBin_qxqy"+"_"+type)
303	Make/O/D/N=(nq) $(folderPath+":"+"iPhiBin2_qxqy"+"_"+type)
304	Make/O/D/N=(nq) $(folderPath+":"+"ePhiBin_qxqy"+"_"+type)
305	Make/O/D/N=(nq) $(folderPath+":"+"ePhiBin2D_qxqy"+"_"+type)
306
307	Wave iPhiBin_qxqy = $(folderPath+":"+"iPhiBin_qxqy_"+type)
308	// Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
309	Wave phiBin_qxqy = $(folderPath+":"+"phiBin_qxqy"+"_"+type)
310	Wave nPhiBin_qxqy = $(folderPath+":"+"nPhiBin_qxqy"+"_"+type)
311	Wave iPhiBin2_qxqy = $(folderPath+":"+"iPhiBin2_qxqy"+"_"+type)
312	Wave ePhiBin_qxqy = $(folderPath+":"+"ePhiBin_qxqy"+"_"+type)
313	Wave ePhiBin2D_qxqy = $(folderPath+":"+"ePhiBin2D_qxqy"+"_"+type)
314
315
316	// TODO:
317	// -- get the q-binning values from the VSANS equivalent
318	// SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr
319	// Variable qc = NumberByKey("QCENTER",keyListStr,"=",";") // q-center
320	// Variable nw = NumberByKey("QDELTA",keyListStr,"=",";") // (in SANS - number of pixels wide)
321
322
323	Variable nphi,dphi,isIn,phiij,iphi
324
325	// TODO: define nphi
326	// dr = 1 //minimum annulus width, keep this fixed at one
327	NVAR numPhiSteps = root:Packages:NIST:VSANS:Globals:gNPhiSteps
328	nphi = numPhiSteps //number of anular sectors is set by users
329	dphi = 360/nphi
330
331
332	iPhiBin_qxqy = 0
333	iPhiBin2_qxqy = 0
334	ePhiBin_qxqy = 0
335	ePhiBin2D_qxqy = 0
336	nPhiBin_qxqy = 0 //number of intensities added to each bin
337
338
339	// 4 panels is currently the only situation
340	//
341	// this needs to be a double loop now...
342	// TODO:
343	// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
344	// -- the solid angle per pixel is not completely implemented.
345	// it will be present for WORK data other than RAW, but not for RAW
346
347	// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
348	if(maskMissing == 1)
349	Print "Mask file not found for at least one detector - so all data is used"
350	endif
351
352	Variable mask_val
353	// use set 1 (no number) only
354	if(nSets >= 1)
355	xDim=DimSize(inten,0)
356	yDim=DimSize(inten,1)
357
358	for(ii=0;ii<xDim;ii+=1)
359	for(jj=0;jj<yDim;jj+=1)
360	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
361	qVal = qTotal[ii][jj]
362
363	isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
364
365	if(isIn) // it's in the annulus somewhere, do something
366	// now I need the qx and qy to find phi
367	if (qy[ii][jj] >= 0)
368	//phiij is in degrees
369	phiij = atan2(qy[ii][jj],qx[ii][jj])*180/Pi //0 to 180 deg
370	else
371	phiij = 360 + atan2(qy[ii][jj],qx[ii][jj])*180/Pi //180 to 360 deg
372	Endif
373	if (phiij > (360-0.5*dphi))
374	phiij -= 360
375	Endif
376	iphi = trunc(phiij/dphi + 1.501) // TODO: why the value of 1.501????
377
378	val = inten[ii][jj]
379
380	if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
381	mask_val = 0
382	else
383	mask_val = mask[ii][jj]
384	endif
385	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
386	iPhiBin_qxqy[iphi-1] += val
387	iPhiBin2_qxqy[iphi-1] += val*val
388	ePhiBin2D_qxqy[iphi-1] += iErr[ii][jj]*iErr[ii][jj]
389	nPhiBin_qxqy[iphi-1] += 1
390	endif
391
392	endif // isIn
393
394	endfor
395	endfor
396
397	endif
398
399	// add in set 2 (set 1 already done)
400	if(nSets >= 2)
401	xDim=DimSize(inten2,0)
402	yDim=DimSize(inten2,1)
403
404	for(ii=0;ii<xDim;ii+=1)
405	for(jj=0;jj<yDim;jj+=1)
406	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
407	qVal = qTotal2[ii][jj]
408
409	isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
410
411	if(isIn) // it's in the annulus somewhere, do something
412	// now I need the qx and qy to find phi
413	if (qy2[ii][jj] >= 0)
414	//phiij is in degrees
415	phiij = atan2(qy2[ii][jj],qx2[ii][jj])*180/Pi //0 to 180 deg
416	else
417	phiij = 360 + atan2(qy2[ii][jj],qx2[ii][jj])*180/Pi //180 to 360 deg
418	Endif
419	if (phiij > (360-0.5*dphi))
420	phiij -= 360
421	Endif
422	iphi = trunc(phiij/dphi + 1.501) // TODO: why the value of 1.501????
423
424	val = inten2[ii][jj]
425
426	if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
427	mask_val = 0
428	else
429	mask_val = mask2[ii][jj]
430	endif
431	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
432	iPhiBin_qxqy[iphi-1] += val
433	iPhiBin2_qxqy[iphi-1] += val*val
434	ePhiBin2D_qxqy[iphi-1] += iErr2[ii][jj]*iErr2[ii][jj]
435	nPhiBin_qxqy[iphi-1] += 1
436	endif
437
438	endif // isIn
439
440	endfor //jj
441	endfor //ii
442
443	endif // set 2
444
445
446	// add in set 3 and 4 (set 1 and 2 already done)
447	if(nSets == 4)
448	xDim=DimSize(inten3,0)
449	yDim=DimSize(inten3,1)
450
451	for(ii=0;ii<xDim;ii+=1)
452	for(jj=0;jj<yDim;jj+=1)
453	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
454	qVal = qTotal3[ii][jj]
455
456	isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
457
458	if(isIn) // it's in the annulus somewhere, do something
459	// now I need the qx and qy to find phi
460	if (qy3[ii][jj] >= 0)
461	//phiij is in degrees
462	phiij = atan2(qy3[ii][jj],qx3[ii][jj])*180/Pi //0 to 180 deg
463	else
464	phiij = 360 + atan2(qy3[ii][jj],qx3[ii][jj])*180/Pi //180 to 360 deg
465	Endif
466	if (phiij > (360-0.5*dphi))
467	phiij -= 360
468	Endif
469	iphi = trunc(phiij/dphi + 1.501) // TODO: why the value of 1.501????
470
471	val = inten3[ii][jj]
472
473	if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
474	mask_val = 0
475	else
476	mask_val = mask3[ii][jj]
477	endif
478	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
479	iPhiBin_qxqy[iphi-1] += val
480	iPhiBin2_qxqy[iphi-1] += val*val
481	ePhiBin2D_qxqy[iphi-1] += iErr3[ii][jj]*iErr3[ii][jj]
482	nPhiBin_qxqy[iphi-1] += 1
483	endif
484
485	endif // isIn
486
487
488	endfor
489	endfor // end of ij loop over set 3
490
491
492	xDim=DimSize(inten4,0)
493	yDim=DimSize(inten4,1)
494
495	for(ii=0;ii<xDim;ii+=1)
496	for(jj=0;jj<yDim;jj+=1)
497	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
498	qVal = qTotal4[ii][jj]
499
500	isIn = V_CloseEnough(qVal,qCtr_Ann,qWidth)
501
502	if(isIn) // it's in the annulus somewhere, do something
503	// now I need the qx and qy to find phi
504	if (qy4[ii][jj] >= 0)
505	//phiij is in degrees
506	phiij = atan2(qy4[ii][jj],qx4[ii][jj])*180/Pi //0 to 180 deg
507	else
508	phiij = 360 + atan2(qy4[ii][jj],qx4[ii][jj])*180/Pi //180 to 360 deg
509	Endif
510	if (phiij > (360-0.5*dphi))
511	phiij -= 360
512	Endif
513	iphi = trunc(phiij/dphi + 1.501) // TODO: why the value of 1.501????
514
515	val = inten4[ii][jj]
516
517	if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
518	mask_val = 0
519	else
520	mask_val = mask4[ii][jj]
521	endif
522	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
523	iPhiBin_qxqy[iphi-1] += val
524	iPhiBin2_qxqy[iphi-1] += val*val
525	ePhiBin2D_qxqy[iphi-1] += iErr4[ii][jj]*iErr4[ii][jj]
526	nPhiBin_qxqy[iphi-1] += 1
527	endif
528
529	endif // isIn
530
531
532	endfor
533	endfor // end ij loop over set 4
534
535	endif // adding sets 3 and 4
536
537
538	// after looping through all of the data on the panels, calculate errors on I(q),
539	// just like in CircSectAve.ipf
540	// TODO:
541	// -- 2D Errors were NOT properly acculumated through reduction, so this loop of calculations is NOT MEANINGFUL (yet)
542	// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
543	for(ii=0;ii<nphi;ii+=1)
544
545	phiBin_qxqy[ii] = dphi*ii
546
547	if(nPhiBin_qxqy[ii] == 0)
548	//no pixels in annuli, data unknown
549	iPhiBin_qxqy[ii] = 0
550	ePhiBin_qxqy[ii] = 1
551	ePhiBin2D_qxqy[ii] = NaN
552	else
553	if(nPhiBin_qxqy[ii] <= 1)
554	//need more than one pixel to determine error
555	iPhiBin_qxqy[ii] /= nPhiBin_qxqy[ii]
556	ePhiBin_qxqy[ii] = 1
557	ePhiBin2D_qxqy[ii] /= (nPhiBin_qxqy[ii])^2
558	else
559	//assume that the intensity in each pixel in annuli is normally distributed about mean...
560	// -- this is correctly calculating the error as the standard error of the mean, as
561	// was always done for SANS as well.
562	iPhiBin_qxqy[ii] /= nPhiBin_qxqy[ii]
563	avesq = iPhiBin_qxqy[ii]^2
564	aveisq = iPhiBin2_qxqy[ii]/nPhiBin_qxqy[ii]
565	var = aveisq-avesq
566	if(var<=0)
567	ePhiBin_qxqy[ii] = 1e-6
568	else
569	ePhiBin_qxqy[ii] = sqrt(var/(nPhiBin_qxqy[ii] - 1))
570	endif
571	// and calculate as it is propagated pixel-by-pixel
572	ePhiBin2D_qxqy[ii] /= (nPhiBin_qxqy[ii])^2
573	endif
574	endif
575	endfor
576
577	ePhiBin2D_qxqy = sqrt(ePhiBin2D_qxqy) // as equation (3) of John's memo
578
579	// find the last non-zero point, working backwards
580	val=nq
581	do
582	val -= 1
583	while((nPhiBin_qxqy[val] == 0) && val > 0)
584
585	// print val, nPhiBin_qxqy[val]
586	DeletePoints val, nq-val, iPhiBin_qxqy,phiBin_qxqy,nPhiBin_qxqy,iPhiBin2_qxqy,ePhiBin_qxqy,ePhiBin2D_qxqy
587
588	if(val == 0)
589	// all the points were deleted
590	return(0)
591	endif
592
593
594	// since the beam center is not always on the detector, many of the low Q bins will have zero pixels
595	// find the first non-zero point, working forwards
596	val = -1
597	do
598	val += 1
599	while(nPhiBin_qxqy[val] == 0)
600	DeletePoints 0, val, iPhiBin_qxqy,phiBin_qxqy,nPhiBin_qxqy,iPhiBin2_qxqy,ePhiBin_qxqy,ePhiBin2D_qxqy
601
602	// ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
603	val = numpnts(nPhiBin_qxqy)-1
604	do
605	if(nPhiBin_qxqy[val] == 0)
606	DeletePoints val, 1, iPhiBin_qxqy,phiBin_qxqy,nPhiBin_qxqy,iPhiBin2_qxqy,ePhiBin_qxqy,ePhiBin2D_qxqy
607	endif
608	val -= 1
609	while(val>0)
610
611	// TODO:
612	// -- is this where I do the resolution calculation? This is where I calculate the resolution in SANS (see CircSectAve)
613	// -- or do I do it as a separate call?
614	// -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
615	//
616
617	SetDataFolder root:
618
619	return(0)
620	End
621
622
623	Proc V_Write1DAnnular(pathStr,folderStr,detGroup,saveName)
624	String pathStr="root:Packages:NIST:VSANS:"
625	String folderStr="SAM"
626	String detGroup = "F"
627	String saveName = "Annular_Data.dat"
628
629
630	V_fWrite1DAnnular(pathStr,folderStr,detGroup,saveName)
631
632	end
633
634	// TODO:
635	// -- this is a temporary solution before a real writer is created
636	// -- resolution is not generated here (and it shouldn't be) since resolution is not known yet.
637	// -- but a real writer will need to be aware of resolution, and there may be different forms
638	//
639	// this will bypass save dialogs
640	// -- AND WILL OVERWITE DATA WITH THE SAME NAME
641	//
642	// V_Write1DData_ITX("root:Packages:NIST:VSANS:",type,saveName,binType)
643	//
644	Function V_fWrite1DAnnular(pathStr,folderStr,detGroup,saveName)
645	String pathStr,folderStr,detGroup,saveName
646
647	String formatStr="",fullpath=""
648	Variable refnum,dialog=1
649
650	SetDataFolder $(pathStr+folderStr)
651
652	if(cmpstr(detGroup,"F") == 0)
653	Wave/Z pw = phiBin_qxqy_FLRTB
654	Wave/Z iw = iPhiBin_qxqy_FLRTB
655	Wave/Z sw = ePhiBin_qxqy_FLRTB
656	else
657	Wave/Z pw = phiBin_qxqy_MLRTB
658	Wave/Z iw = iPhiBin_qxqy_MLRTB
659	Wave/Z sw = ePhiBin_qxqy_MLRTB
660	endif
661
662
663	String dataSetFolderParent,basestr
664
665	//make sure the waves exist
666
667	if(WaveExists(pw) == 0)
668	SetDataFolder root:
669	Abort "q is missing"
670	endif
671	if(WaveExists(iw) == 0)
672	SetDataFolder root:
673	Abort "i is missing"
674	endif
675	if(WaveExists(sw) == 0)
676	SetDataFolder root:
677	Abort "s is missing"
678	endif
679	// if(WaveExists(resw) == 0)
680	// Abort "Resolution information is missing."
681	// endif
682
683	// Duplicate/O qw qbar,sigQ,fs
684	// if(dimsize(resW,1) > 4)
685	// //it's USANS put -dQv back in the last 3 columns
686	// NVAR/Z dQv = USANS_dQv
687	// if(NVAR_Exists(dQv) == 0)
688	// SetDataFolder root:
689	// Abort "It's USANS data, and I don't know what the slit height is."
690	// endif
691	// sigQ = -dQv
692	// qbar = -dQv
693	// fs = -dQv
694	// else
695	// //it's SANS
696	// sigQ = resw[p][0]
697	// qbar = resw[p][1]
698	// fs = resw[p][2]
699	// endif
700	//
701
702	PathInfo catPathName
703	fullPath = S_Path + saveName
704
705	Open refnum as fullpath
706
707	fprintf refnum,"Annular data written from folder %s on %s\r\n",folderStr,(date()+" "+time())
708
709	// TODO -- make this work for 6-columns (or??)
710	// formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
711	// fprintf refnum, "The 6 columns are \| Q (1/A) \| I(Q) (1/cm) \| std. dev. I(Q) (1/cm) \| sigmaQ \| meanQ \| ShadowFactor\|\r\n"
712	// wfprintf refnum,formatStr,qw,iw,sw,sigQ,qbar,fs
713
714	//currently, only three columns
715	formatStr = "%15.4g %15.4g %15.4g\r\n"
716	fprintf refnum, "The 3 columns are \| Phi (degrees) \| I(phi) (1/cm) \| std. dev. I(phi) (1/cm)\r\n"
717
718	wfprintf refnum,formatStr,pw,iw,sw
719	Close refnum
720
721	// KillWaves/Z sigQ,qbar,fs
722	Print "Data written to: ",fullpath
723
724	SetDataFolder root:
725	return(0)
726	End

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

Download in other formats:

Original Format