Context Navigation

V_IQ_Annular.ipf @ 1095

Last change on this file since 1095 was 1095, checked in by srkline, 4 years ago

Added two model functions for white beam smearing.

Many other small changes for processing of the back detector, shuffling of VSANS menu items, and consistent naming of V_ procedures.

File size: 25.3 KB

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

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

Download in other formats:

Original Format