Context Navigation

VC_DetectorBinning_Utils.ipf @ 1242

Last change on this file since 1242 was 1242, checked in by srkline, 2 years ago
updating the IgorVersion? pragma to v7.0 for all files to be consistent.
File size: 47.6 KB

Line
1	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
2	#pragma IgorVersion = 7.00
3
4	/////////////////////////
5	//
6	// Utility functions to:
7	// calculate Q, Qx, Qy, Qz
8	// fill the detector panels with simulated data (the model functions are here)
9	// bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
10	//
11	/////////////////////////
12
13
14
15
16	// x- hard wired for a sphere - change this to allow minimal selections and altering of coefficients
17	// x- add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
18	//
19	// NOTE - this is a VCALC only routine, so it's not been made completely generic
20	//
21	Function FillPanel_wModelData(det,qTot,type)
22	Wave det,qTot
23	String type
24
25	// SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
26
27	// q-values and detector arrays already allocated and calculated
28	Duplicate/O det tmpInten,tmpSig,prob_i
29
30	Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
31
32	//imon = VC_BeamIntensity()*CountTime
33	imon = VCALC_getImon() //TODO: currently from the panel, not calculated
34	trans = 0.8
35	thick = 0.1
36
37	// need SDD
38	// need pixel dimensions
39	// nominal sdd in cm, offset in cm, want result in cm !
40
41	sdd = VC_getSDD(type) // setback is already included VCALC_getTopBottomSDDSetback(type) // result is sdd in [cm]
42
43	pixSizeX = VCALC_getPixSizeX(type) // cm
44	pixSizeY = VCALC_getPixSizeY(type)
45
46
47	//?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
48	// keep the panel to keep people used to using it.)
49	// peak @ 0.1 ~ AgBeh
50	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}
51	//
52	// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
53	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}
54	//
55	//
56	Variable addEmpBgd=0
57
58	// TODOHIGHRES
59	// this is a slow step - try to figure out how to multithread this efficiently. simply adding the
60	// keyword has little effect. maybe only do this for "B", maybe rewrite the calculation to not use pq indexing
61	//
62	String funcStr = VCALC_getModelFunctionStr()
63	strswitch(funcStr)
64	case "Big Debye":
65	tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
66	break
67	case "Big Sphere":
68	tmpInten = VC_SphereForm(1,2000,1e-6,0.01,qTot[p][q])
69	break
70	case "Debye":
71	MultiThread tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
72	break
73	case "Sphere":
74	tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])
75	break
76	case "AgBeh":
77	tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
78	break
79	case "Vycor":
80	tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
81	break
82	case "Empty Cell":
83	tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
84	break
85	case "Blocked Beam":
86	tmpInten = VC_BlockedBeam(0.01,qTot[p][q])
87	break
88	case "Debye +":
89	tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
90	addEmpBgd = 1
91	break
92	case "AgBeh +":
93	tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
94	addEmpBgd = 1
95	break
96	case "Empty Cell +":
97	tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
98	tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
99	break
100	default:
101	tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
102	endswitch
103
104
105	if(addEmpBgd == 1)
106	tmpInten += VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
107	tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
108	endif
109
110
111	// x- this is faked to get around the singularity at the center of the back detector
112	//
113	//
114	if(cmpstr(type,"B") == 0)
115	Variable nx,ny,px,py
116	nx = VCALC_get_nPix_X(type)
117	ny = VCALC_get_nPix_Y(type)
118	px = trunc(nx/2)
119	py = trunc(ny/2)
120
121	tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
122	endif
123
124
125
126	///////////////
127	// // calculate the scattering cross section simply to be able to estimate the transmission
128	// Variable sig_sas=0
129	//
130	// // remember that the random deviate is the coherent portion ONLY - the incoherent background is
131	// // subtracted before the calculation.
132	// CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
133	//
134	// if(sig_sas > 100)
135	// DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
136	// endif
137	//
138	// // calculate the multiple scattering fraction for display (10/2009)
139	// Variable ii,nMax=10,tau
140	// mScat=0
141	// tau = thick*sig_sas
142	// // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
143	// // neutrons out of those that were scattered
144	// for(ii=2;ii<nMax;ii+=1)
145	// mScat += tau^(ii)/factorial(ii)
146	//// print tau^(ii)/factorial(ii)
147	// endfor
148	// estTrans = exp(-1thicksig_sas) //thickness and sigma both in units of cm
149	// mscat *= (estTrans)/(1-estTrans)
150	//
151	//// if(mScat > 0.1) // Display warning
152	//
153	// Print "Sig_sas = ",sig_sas
154	////////////////////
155
156	MultiThread prob_i = transthickpixSizeXpixSizeY/(sdd)^2tmpInten //probability of a neutron in q-bin(i)
157
158	tmpInten = (imon)*prob_i //tmpInten is not the model calculation anymore!!
159
160
161	/// **** can I safely assume a Gaussian error in the count rate??
162	MultiThread tmpSig = sqrt(tmpInten) // corrected based on John's memo, from 8/9/99
163
164	MultiThread tmpInten += gnoise(tmpSig)
165	MultiThread tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q] // MAR 2013 -- is this the right thing to do
166	tmpInten = trunc(tmpInten)
167
168
169	det = tmpInten
170
171	// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
172	// det /= transthickpixSizeXpixSizeY/(sdd)^2imon
173
174
175	KillWaves/Z tmpInten,tmpSig,prob_i
176	SetDataFolder root:
177
178	return(0)
179	End
180
181
182	// For a given detector panel, calculate the q-values
183	// -work with everything as arrays
184	// Input needed:
185	// detector data
186	// detector type (LRTB?)
187	// beam center (may be off the detector)
188	// SDD
189	// lambda
190	//
191	// pixel dimensions for detector type (global constants)
192	// - data dimensions read directly from array
193	//
194	// --What is calculated:
195	// array of Q
196	// array of qx,qy,qz
197	// array of error already exists
198	//
199	//
200	// -- sdd in cm
201	// -- lambda in Angstroms
202	Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
203	Wave data,qTot,qx,qy,qz
204	Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
205
206	// loop over the array and calculate the values - this is done as a wave assignment
207	// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf
208	MultiThread qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
209	MultiThread qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
210	MultiThread qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
211	MultiThread qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
212
213	return(0)
214	End
215
216
217	// for testing, a version that will calculate the q-arrays for VCALC based on whatever nonlinear coefficients
218	// exist in the RAW data folder
219	//
220	// reverts to the "regular" linear detector if waves not found or a flag is set
221	//
222	// need to call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
223	//
224	// ***** everything passed in is [cm], except for wavelength [A]
225	//
226	// **** TODO :: calibration constants are still in [mm]
227	//
228	//
229	// TODO:
230	// -- tube width is hard-wired in
231	//
232	//
233	Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
234	Wave data,qTot,qx,qy,qz
235	Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
236	String detStr
237
238	String destPath = "root:Packages:NIST:VSANS:VCALC"
239
240	// be sure that the real distance waves exist
241	// TODO -- this may not be the best location?
242
243	Variable tube_width = 8.4 // TODO: UNITS!!! Hard-wired value in [mm]
244
245	String orientation
246	Variable dimX,dimY
247	dimX = DimSize(data,0)
248	dimY = DimSize(data,1)
249	if(dimX > dimY)
250	orientation = "horizontal"
251	else
252	orientation = "vertical"
253	endif
254
255
256	Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
257	Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
258	NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
259
260	if(kBCTR_CM && cmpstr(detStr,"B") != 0)
261	if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
262	// beam ctr is in cm already
263
264	// calculate all of the q-values
265	qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
266	qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
267	qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
268	qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
269
270	// Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
271	// Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
272	else
273	// do the q-calculation using linear detector
274	//VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
275	qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
276	qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
277	qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
278	qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
279	endif
280
281
282	else
283	// using the old calculation with beam center in pixels
284	if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
285	// convert the beam centers to mm
286	// String orientation
287	Variable newX,newY
288	dimX = DimSize(data_realDistX,0)
289	dimY = DimSize(data_realDistX,1)
290	if(dimX > dimY)
291	orientation = "horizontal"
292	else
293	orientation = "vertical"
294	endif
295
296
297	//
298	if(cmpstr(orientation,"vertical")==0)
299	// this is data dimensioned as (Ntubes,Npix)
300	newX = tube_width*xCtr
301	newY = data_realDistY[0][yCtr]
302	else
303	// this is data (horizontal) dimensioned as (Npix,Ntubes)
304	newX = data_realDistX[xCtr][0]
305	newY = tube_width*yCtr
306	endif
307
308	//if detector "B", different calculation for the centers (not tubes)
309	if(cmpstr(detStr,"B")==0)
310	newX = data_realDistX[xCtr][0]
311	newY = data_realDistY[0][yCtr]
312	//newX = xCtr
313	//newY = yCtr
314	endif
315
316	// calculate all of the q-values
317	qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
318	qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
319	qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
320	qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
321
322	// Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
323	// Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
324	else
325	// do the q-calculation using linear detector
326	VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
327	endif
328
329	endif
330
331	KillWaves/Z tmpCalib,tmpCalibX,tmpCalibY
332
333	return(0)
334	End
335
336
337
338	// make the data_realDistX,Y Waves that are needed for the calculation of q
339	Function VC_MakeRealDistXYWaves(data,detStr)
340	Wave data
341	String detStr
342
343	String destPath = "root:Packages:NIST:VSANS:VCALC"
344
345	// calibration waves do not exist yet, so make some fake ones '
346	// do I count on the orientation as an input, or do I just figure it out on my own?
347	String orientation
348	Variable dimX,dimY
349	dimX = DimSize(data,0)
350	dimY = DimSize(data,1)
351	if(dimX > dimY)
352	orientation = "horizontal"
353	else
354	orientation = "vertical"
355	endif
356
357	if(cmpstr(orientation,"vertical")==0)
358	Make/O/D/N=(3,48) tmpCalib
359	// for the "tall" L/R banks
360	tmpCalib[0][] = -512
361	tmpCalib[1][] = 8
362	tmpCalib[2][] = 0
363	else
364	Make/O/D/N=(3,48) tmpCalib
365	// for the "short" T/B banks
366	tmpCalib[0][] = -256
367	tmpCalib[1][] = 4
368	tmpCalib[2][] = 0
369	endif
370	// override if back panel
371	if(cmpstr(detStr,"B") == 0)
372	// and for the back detector "B"
373	Make/O/D/N=3 tmpCalibX,tmpCalibY
374	tmpCalibX[0] = VCALC_getPixSizeX(detStr) // pixel size in [cm] VCALC_getPixSizeX(detStr) is [cm]
375	tmpCalibX[1] = 1
376	tmpcalibX[2] = 10000
377	tmpCalibY[0] = VCALC_getPixSizeY(detStr) // pixel size in [cm] VCALC_getPixSizeX(detStr) is [cm]
378	tmpCalibY[1] = 1
379	tmpcalibY[2] = 10000
380	endif
381
382
383	// Wave w_calib = V_getDetTube_spatialCalib("VCALC",detStr)
384	Variable tube_width = 8.4 // TODO: UNITS!!! Hard-wired value in [mm]
385	if(cmpstr(detStr,"B") == 0)
386	V_NonLinearCorrection_B("VCALC",data,tmpCalibX,tmpCalibY,detStr,destPath)
387	// beam center is in pixels, so use the old routine
388	V_ConvertBeamCtrPix_to_mmB("VCALC","B",destPath)
389	else
390	V_NonLinearCorrection("VCALC",data,tmpCalib,tube_width,detStr,destPath)
391	endif
392
393
394	return(0)
395	End
396
397
398
399	//////////////////////
400	// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
401	// accout the different x/y pixel sizes and the beam center not being on the detector -
402	// off a different edge for each LRTB type
403	/////////////////////
404
405	//function to calculate the overall q-value, given all of the necesary trig inputs
406	//and are in detector coordinates (1,128) rather than axis values
407	//the pixel locations need not be integers, reals are ok inputs
408	//sdd is in [cm]
409	//wavelength is in Angstroms
410	//
411	//returned magnitude of Q is in 1/Angstroms
412	//
413	//
414	Threadsafe Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
415	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
416
417	Variable dx,dy,qval,two_theta,dist
418
419
420	dx = (xaxval - xctr)*pixSizeX //delta x in cm
421	dy = (yaxval - yctr)*pixSizeY //delta y in cm
422	dist = sqrt(dx^2 + dy^2)
423
424	two_theta = atan(dist/sdd)
425
426	qval = 4Pi/lamsin(two_theta/2)
427
428	return qval
429	End
430
431	//calculates just the q-value in the x-direction on the detector
432	//input/output is the same as CalcQval()
433	//ALL inputs are in detector coordinates
434	//
435	//sdd is in [cm]
436	//wavelength is in Angstroms
437	//
438	// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
439	// now properly accounts for qz
440	//
441	Threadsafe Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
442	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
443
444	Variable qx,qval,phi,dx,dy,dist,two_theta
445
446	qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
447
448	// sdd *=100 //convert to cm
449	dx = (xaxval - xctr)*pixSizeX //delta x in cm
450	dy = (yaxval - yctr)*pixSizeY //delta y in cm
451	phi = V_FindPhi(dx,dy)
452
453	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
454	dist = sqrt(dx^2 + dy^2)
455	two_theta = atan(dist/sdd)
456
457	qx = qvalcos(two_theta/2)cos(phi)
458
459	return qx
460	End
461
462	//calculates just the q-value in the y-direction on the detector
463	//input/output is the same as CalcQval()
464	//ALL inputs are in detector coordinates
465	//sdd is in [cm]
466	//wavelength is in Angstroms
467	//
468	// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
469	// now properly accounts for qz
470	//
471	Threadsafe Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
472	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
473
474	Variable dy,qval,dx,phi,qy,dist,two_theta
475
476	qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
477
478	// sdd *=100 //convert to cm
479	dx = (xaxval - xctr)*pixSizeX //delta x in cm
480	dy = (yaxval - yctr)*pixSizeY //delta y in cm
481	phi = V_FindPhi(dx,dy)
482
483	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
484	dist = sqrt(dx^2 + dy^2)
485	two_theta = atan(dist/sdd)
486
487	qy = qvalcos(two_theta/2)sin(phi)
488
489	return qy
490	End
491
492	//calculates just the z-component of the q-vector, not measured on the detector
493	//input/output is the same as CalcQval()
494	//ALL inputs are in detector coordinates
495	//sdd is in [cm]
496	//wavelength is in Angstroms
497	//
498	// not actually used, but here for completeness if anyone asks
499	//
500	Threadsafe Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
501	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
502
503	Variable dy,qval,dx,phi,qz,dist,two_theta
504
505	qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
506
507	// sdd *=100 //convert to cm
508
509	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
510	dx = (xaxval - xctr)*pixSizeX //delta x in cm
511	dy = (yaxval - yctr)*pixSizeY //delta y in cm
512	dist = sqrt(dx^2 + dy^2)
513	two_theta = atan(dist/sdd)
514
515	qz = qval*sin(two_theta/2)
516
517	return qz
518	End
519
520	//phi is defined from +x axis, proceeding CCW around [0,2Pi]
521	Threadsafe Function V_FindPhi(vx,vy)
522	variable vx,vy
523
524	variable phi
525
526	phi = atan(vy/vx) //returns a value from -pi/2 to pi/2
527
528	// special cases
529	if(vx==0 && vy > 0)
530	return(pi/2)
531	endif
532	if(vx==0 && vy < 0)
533	return(3*pi/2)
534	endif
535	if(vx >= 0 && vy == 0)
536	return(0)
537	endif
538	if(vx < 0 && vy == 0)
539	return(pi)
540	endif
541
542
543	if(vx > 0 && vy > 0)
544	return(phi)
545	endif
546	if(vx < 0 && vy > 0)
547	return(phi + pi)
548	endif
549	if(vx < 0 && vy < 0)
550	return(phi + pi)
551	endif
552	if( vx > 0 && vy < 0)
553	return(phi + 2*pi)
554	endif
555
556	return(phi)
557	end
558
559	Function VC_SphereForm(scale,radius,delrho,bkg,x)
560	Variable scale,radius,delrho,bkg
561	Variable x
562
563	// variables are:
564	//[0] scale
565	//[1] radius (A)
566	//[2] delrho (A-2)
567	//[3] background (cm-1)
568
569	// Variable scale,radius,delrho,bkg
570	// scale = w[0]
571	// radius = w[1]
572	// delrho = w[2]
573	// bkg = w[3]
574
575
576	// calculates scale * f^2/Vol where f=Vol3delrho*((sin(qr)-qrcos(qr))/qr^3
577	// and is rescaled to give [=] cm^-1
578
579	Variable bes,f,vol,f2
580	//
581	//handle q==0 separately
582	If(x==0)
583	f = 4/3piradius^3delrhodelrhoscale1e8 + bkg
584	return(f)
585	Endif
586
587	// bes = 3(sin(xradius)-xradiuscos(x*radius))/x^3/radius^3
588
589	bes = 3sqrt(pi/(2xradius))BesselJ(1.5,xradius)/(xradius)
590
591	vol = 4pi/3radius^3
592	f = volbesdelrho // [=] A
593	// normalize to single particle volume, convert to 1/cm
594	f2 = f * f / vol * 1.0e8 // [=] 1/cm
595
596	return (scale*f2+bkg) // Scale, then add in the background
597
598	End
599
600	ThreadSafe Function VC_Debye(scale,rg,bkg,x)
601	Variable scale,rg,bkg
602	Variable x
603
604	// variables are:
605	//[0] scale factor
606	//[1] radius of gyration [A]
607	//[2] background [cm-1]
608
609	// calculates (scale*debye)+bkg
610	Variable Pq,qr2
611
612	qr2=(x*rg)^2
613	Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
614
615	//scale
616	Pq *= scale
617	// then add in the background
618	return (Pq+bkg)
619	End
620
621	// a sum of a power law and debye to approximate the scattering from a real empty cell
622	//
623	// make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
624	//
625	Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
626	Variable aa,mm,scale,rg,bkg
627	Variable x
628
629	// variables are:
630	//[0] = A
631	//[1] = power m
632	//[2] scale factor
633	//[3] radius of gyration [A]
634	//[4] background [cm-1]
635
636	Variable Iq
637
638	// calculates (scale*debye)+bkg
639	Variable Pq,qr2
640
641	// if(x*Rg < 1e-3) //added Oct 2008 to avoid numerical errors at low arg values
642	// return(scale+bkg)
643	// endif
644
645	Iq = aa*x^-mm
646
647	qr2=(x*rg)^2
648	Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
649
650	//scale
651	Pq *= scale
652	// then add the terms up
653	return (Iq + Pq + bkg)
654	End
655
656	// blocked beam
657	//
658	Function VC_BlockedBeam(bkg,x)
659	Variable bkg
660	Variable x
661
662	return (bkg)
663	End
664
665
666	//
667	// a broad peak to simulate silver behenate or vycor
668	//
669	// peak @ 0.1 ~ AgBeh
670	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}
671	//
672	//
673	// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
674	// Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}
675	//
676	//
677	Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
678	Variable aa,nn,cc,LL,Qzero,mm,bgd
679	Variable x
680
681	// variables are:
682	//[0] Porod term scaling
683	//[1] Porod exponent
684	//[2] Lorentzian term scaling
685	//[3] Lorentzian screening length [A]
686	//[4] peak location [1/A]
687	//[5] Lorentzian exponent
688	//[6] background
689
690	// local variables
691	Variable inten, qval
692	// x is the q-value for the calculation
693	qval = x
694	// do the calculation and return the function value
695
696	inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
697
698	Return (inten)
699
700	End
701
702	//
703	// updated to new folder structure Feb 2016
704	// folderStr = RAW,SAM, VCALC or other
705	// detStr is the panel identifer "ML", etc.
706	//
707	Function V_SetDeltaQ(folderStr,detStr)
708	String folderStr,detStr
709
710	NVAR binWidth = root:Packages:NIST:VSANS:Globals:gBinWidth
711
712	Variable isVCALC
713	if(cmpstr(folderStr,"VCALC") == 0)
714	isVCALC = 1
715	endif
716
717	String folderPath = "root:Packages:NIST:VSANS:"+folderStr
718	String instPath = ":entry:instrument:detector_"
719
720	if(isVCALC)
721	WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) // 2D detector data
722	else
723	Wave inten = V_getDetectorDataW(folderStr,detStr)
724	endif
725
726	Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
727	Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
728
729	Variable xDim,yDim,delQ
730
731	xDim=DimSize(inten,0)
732	yDim=DimSize(inten,1)
733
734	if(xDim<yDim)
735	delQ = abs(qx[0][0] - qx[1][0])/2
736	else
737	delQ = abs(qy[0][1] - qy[0][0])/2
738	endif
739
740	// TODOHIGHRES
741	// -- is this how I want to handle the too-fine resolution of 1x1 binning?
742	NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
743
744	if(cmpstr(detStr,"B") == 0 && gHighResBinning == 1)
745	delQ = 4*delQ
746	Print "Reset delta Q for binning the back detector to 4x pix = ",delQ
747	endif
748
749	// multiply the deltaQ by the binWidth (=multiple of pixels)
750	// this defaults to 1, and is set in VSANS preferences
751	delQ *= binWidth
752
753	// set the global
754	Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
755	// Print "SET delQ = ",delQ," for ",type
756
757	return(delQ)
758	end
759
760
761	//TODO -- need a switch here to dispatch to the averaging type
762	Proc VC_BinQxQy_to_1D(folderStr,type)
763	String folderStr
764	String type
765	// Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
766	// Prompt type,"detector identifier"
767
768
769	VC_fDoBinning_QxQy2D(folderStr, type)
770
771
772	/// this is for a tall, narrow slit mode
773	// VC_fBinDetector_byRows(folderStr,type)
774
775	End
776
777
778	// folderStr is RAW, VCALC, SAM, etc.
779	// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
780	//
781	Proc VC_Graph_1D_detType(folderStr,type)
782	String folderStr,type
783
784	SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
785
786	Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
787	ModifyGraph mirror=2,grid=1,log=1
788	ModifyGraph mode=4,marker=19,msize=2
789	// ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy) // for simulations, I don't have 2D uncertainty
790	ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
791	legend
792
793	SetDataFolder root:
794
795	End
796
797
798
799	//////////
800	//
801	// Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
802	// - each pixel QxQyQz has been calculated beforehand
803	// - if multiple panels are selected to be combined, it is done here during the binning
804	// - the setting of deltaQ step is still a little suspect (TODO)
805	//
806	//
807	// see the equivalent function in PlotUtils2D_v40.ipf
808	//
809	//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
810	//
811	// this has been modified to accept different detector panels and to take arrays
812	// -- type = FL or FR or...other panel identifiers
813	//
814	// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
815	//
816	//
817	// updated Feb2016 to take new folder structure
818	// TODO
819	// -- VERIFY
820	// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
821	// but this is a lousy choice.
822	// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
823	// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
824	// a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
825	// step where data is added to a WORK file (see Raw_to_Work())
826	//
827	//
828	// TODO:
829	// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
830	// hard-wired. these need to be corrected before even the pinhole resolution is correct
831	// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
832	// and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
833	// is separate functions that are resolution-aware.
834	//
835	//
836	// folderStr = WORK folder, type = the binning type (may include multiple detectors)
837	Function VC_fDoBinning_QxQy2D(folderStr,type,collimationStr)
838	String folderStr,type,collimationStr
839
840	Variable nSets = 0
841	Variable xDim,yDim
842	Variable ii,jj
843	Variable qVal,nq,var,avesq,aveisq
844	Variable binIndex,val,isVCALC=0,maskMissing
845
846	String folderPath = "root:Packages:NIST:VSANS:"+folderStr
847	String instPath = ":entry:instrument:detector_"
848	String detStr
849
850	if(cmpstr(folderStr,"VCALC") == 0)
851	isVCALC = 1
852	endif
853
854	detStr = type
855
856	// now switch on the type to determine which waves to declare and create
857	// since there may be more than one panel to step through. There may be two, there may be four
858	//
859
860	// TODO:
861	// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
862	//
863	// assume that the mask files are missing unless we can find them. If VCALC data,
864	// then the Mask is missing by definition
865	maskMissing = 1
866
867	strswitch(type) // string switch
868	// case "FL": // execute if case matches expression
869	// case "FR":
870	// detStr = type
871	// if(isVCALC)
872	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
873	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
874	// else
875	// Wave inten = V_getDetectorDataW(folderStr,detStr)
876	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
877	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
878	// if(WaveExists(mask) == 1)
879	// maskMissing = 0
880	// endif
881	//
882	// endif
883	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
884	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
885	// nSets = 1
886	// break
887
888	// case "FT":
889	// case "FB":
890	// detStr = type
891	// if(isVCALC)
892	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
893	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
894	// else
895	// Wave inten = V_getDetectorDataW(folderStr,detStr)
896	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
897	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
898	// if(WaveExists(mask) == 1)
899	// maskMissing = 0
900	// endif
901	// endif
902	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
903	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
904	// nSets = 1
905	// break
906
907	// case "ML":
908	// case "MR":
909	// detStr = type
910	// if(isVCALC)
911	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
912	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
913	// else
914	// Wave inten = V_getDetectorDataW(folderStr,detStr)
915	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
916	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
917	// if(WaveExists(mask) == 1)
918	// maskMissing = 0
919	// endif
920	// endif
921	// //TODO:
922	// // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
923	// // determined specifically for that panel (currently using one tube width as deltaQ)
924	// // -- this is repeated multiple times in this switch
925	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
926	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
927	// nSets = 1
928	// break
929
930	// case "MT":
931	// case "MB":
932	// detStr = type
933	// if(isVCALC)
934	// WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
935	// WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
936	// else
937	// Wave inten = V_getDetectorDataW(folderStr,detStr)
938	// Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
939	// Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
940	// if(WaveExists(mask) == 1)
941	// maskMissing = 0
942	// endif
943	// endif
944	// NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
945	// Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
946	// nSets = 1
947	// break
948
949	// only one panel, simply pick that panel and move on out of the switch
950	case "FL":
951	case "FR":
952	case "FT":
953	case "FB":
954	case "ML":
955	case "MR":
956	case "MT":
957	case "MB":
958	case "B":
959	if(isVCALC)
960	WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
961	WAVE/Z iErr = $("iErr_"+detStr) // 2D errors -- may not exist, especially for simulation
962	else
963	Wave inten = V_getDetectorDataW(folderStr,detStr)
964	Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
965	endif
966	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
967	if(WaveExists(mask) == 1)
968	maskMissing = 0
969	endif
970	NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
971	Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr) // 2D q-values
972	nSets = 1
973	break
974
975	case "FLR":
976	// detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
977	// TODO
978	// -- see if I can un-hard-wire some of this below when more than one panel is combined
979	if(isVCALC)
980	WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
981	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
982	WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
983	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
984	else
985	Wave inten = V_getDetectorDataW(folderStr,"FL")
986	Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
987	Wave inten2 = V_getDetectorDataW(folderStr,"FR")
988	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
989	endif
990	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
991	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
992	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
993	maskMissing = 0
994	endif
995
996	NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
997
998	Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values
999	Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values
1000
1001	nSets = 2
1002	break
1003
1004	case "FTB":
1005	if(isVCALC)
1006	WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
1007	WAVE/Z iErr = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
1008	WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
1009	WAVE/Z iErr2 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
1010	else
1011	Wave inten = V_getDetectorDataW(folderStr,"FT")
1012	Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
1013	Wave inten2 = V_getDetectorDataW(folderStr,"FB")
1014	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
1015	endif
1016	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1017	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1018	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1019	maskMissing = 0
1020	endif
1021
1022	NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
1023
1024	Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values
1025	Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values
1026
1027	nSets = 2
1028	break
1029
1030	case "FLRTB":
1031	if(isVCALC)
1032	WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
1033	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
1034	WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
1035	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
1036	WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
1037	WAVE/Z iErr3 = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
1038	WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
1039	WAVE/Z iErr4 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
1040	else
1041	Wave inten = V_getDetectorDataW(folderStr,"FL")
1042	Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
1043	Wave inten2 = V_getDetectorDataW(folderStr,"FR")
1044	Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
1045	Wave inten3 = V_getDetectorDataW(folderStr,"FT")
1046	Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
1047	Wave inten4 = V_getDetectorDataW(folderStr,"FB")
1048	Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
1049
1050	endif
1051	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
1052	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
1053	Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1054	Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1055	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1056	maskMissing = 0
1057	endif
1058
1059	NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1060
1061	Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL") // 2D q-values
1062	Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR") // 2D q-values
1063	Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT") // 2D q-values
1064	Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB") // 2D q-values
1065
1066	nSets = 4
1067	break
1068
1069	case "MLR":
1070	if(isVCALC)
1071	WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1072	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
1073	WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1074	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
1075	else
1076	Wave inten = V_getDetectorDataW(folderStr,"ML")
1077	Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1078	Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1079	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1080
1081	endif
1082	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1083	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1084	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1085	maskMissing = 0
1086	endif
1087
1088	NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1089
1090	Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values
1091	Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values
1092
1093	nSets = 2
1094	break
1095
1096	case "MTB":
1097	if(isVCALC)
1098	WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1099	WAVE/Z iErr = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
1100	WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1101	WAVE/Z iErr2 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
1102	else
1103	Wave inten = V_getDetectorDataW(folderStr,"MT")
1104	Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1105	Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1106	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1107
1108	endif
1109	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1110	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1111	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1112	maskMissing = 0
1113	endif
1114
1115	NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1116
1117	Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values
1118	Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values
1119
1120	nSets = 2
1121	break
1122
1123	case "MLRTB":
1124	if(isVCALC)
1125	WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1126	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
1127	WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1128	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
1129	WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1130	WAVE/Z iErr3 = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
1131	WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1132	WAVE/Z iErr4 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
1133	else
1134	Wave inten = V_getDetectorDataW(folderStr,"ML")
1135	Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1136	Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1137	Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1138	Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1139	Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1140	Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1141	Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1142
1143	endif
1144	Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1145	Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1146	Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1147	Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1148	if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1149	maskMissing = 0
1150	endif
1151
1152	NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1153
1154	Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML") // 2D q-values
1155	Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR") // 2D q-values
1156	Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT") // 2D q-values
1157	Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB") // 2D q-values
1158
1159	nSets = 4
1160	break
1161
1162	default:
1163	nSets = 0
1164	Print "ERROR ---- type is not recognized "
1165	endswitch
1166
1167	// Print "delQ = ",delQ," for ",type
1168
1169	if(nSets == 0)
1170	SetDataFolder root:
1171	return(0)
1172	endif
1173
1174
1175	// RAW data is currently read in and the 2D error wave is correctly generated
1176	// 2D error is propagated through all reduction steps, but I have not
1177	// verified that it is an exact duplication of the 1D error
1178	//
1179	//
1180	//
1181	// IF ther is no 2D error wave present for some reason, make a fake one
1182	if(WaveExists(iErr)==0 && WaveExists(inten) != 0)
1183	Duplicate/O inten,iErr
1184	Wave iErr=iErr
1185	// iErr = 1+sqrt(inten+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1186	iErr = sqrt(inten+0.75) // TODO -- here I'm just using some fictional value
1187	endif
1188	if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1189	Duplicate/O inten2,iErr2
1190	Wave iErr2=iErr2
1191	// iErr2 = 1+sqrt(inten2+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1192	iErr2 = sqrt(inten2+0.75) // TODO -- here I'm just using some fictional value
1193	endif
1194	if(WaveExists(iErr3)==0 && WaveExists(inten3) != 0)
1195	Duplicate/O inten3,iErr3
1196	Wave iErr3=iErr3
1197	// iErr3 = 1+sqrt(inten3+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1198	iErr3 = sqrt(inten3+0.75) // TODO -- here I'm just using some fictional value
1199	endif
1200	if(WaveExists(iErr4)==0 && WaveExists(inten4) != 0)
1201	Duplicate/O inten4,iErr4
1202	Wave iErr4=iErr4
1203	// iErr4 = 1+sqrt(inten4+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
1204	iErr4 = sqrt(inten4+0.75) // TODO -- here I'm just using some fictional value
1205	endif
1206
1207	// TODO -- nq will need to be larger, once the back detector is installed
1208	//
1209	// note that the back panel of 320x320 (1mm res) results in 447 data points!
1210	// - so I upped nq to 600
1211
1212	if(cmpstr(type,"B") == 0)
1213	nq = 8000
1214	else
1215	nq=600
1216	endif
1217
1218	//****TODO**** -- where to put the averaged data -- right now, folderStr is forced to ""
1219	// SetDataFolder $("root:"+folderStr) //should already be here, but make sure...
1220	Make/O/D/N=(nq) $(folderPath+":"+"iBin_qxqy"+"_"+type)
1221	Make/O/D/N=(nq) $(folderPath+":"+"qBin_qxqy"+"_"+type)
1222	Make/O/D/N=(nq) $(folderPath+":"+"nBin_qxqy"+"_"+type)
1223	Make/O/D/N=(nq) $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1224	Make/O/D/N=(nq) $(folderPath+":"+"eBin_qxqy"+"_"+type)
1225	Make/O/D/N=(nq) $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1226
1227	Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1228	Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1229	Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1230	Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1231	Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1232	Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1233
1234
1235	// delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2)) //use bins of 1 pixel width
1236	// TODO: not sure if I want to set dQ in x or y direction...
1237	// the short dimension is the 8mm tubes, use this direction as dQ?
1238	// but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0]
1239	// WRT the beam center. use qx or qy directly. Still not happy with this way...
1240
1241
1242	qBin_qxqy[] = p*delQ
1243	SetScale/P x,0,delQ,"",qBin_qxqy //allows easy binning
1244
1245	iBin_qxqy = 0
1246	iBin2_qxqy = 0
1247	eBin_qxqy = 0
1248	eBin2D_qxqy = 0
1249	nBin_qxqy = 0 //number of intensities added to each bin
1250
1251	// now there are situations of:
1252	// 1 panel
1253	// 2 panels
1254	// 4 panels
1255	//
1256	// this needs to be a double loop now...
1257	// TODO:
1258	// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1259	// -- verify the 2D error propagation by reducing it to 1D error
1260	//
1261	//
1262	// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1263	//
1264	// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
1265
1266	//
1267	// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1268	if(maskMissing == 1)
1269	Print "Mask file not found for at least one detector - so all data is used"
1270	endif
1271
1272	NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1273	if(gIgnoreDetB && cmpstr(type,"B") == 0)
1274	maskMissing = 1
1275	Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)"
1276	endif
1277
1278	Variable mask_val
1279	// use set 1 (no number) only
1280	if(nSets >= 1)
1281	xDim=DimSize(inten,0)
1282	yDim=DimSize(inten,1)
1283
1284	for(ii=0;ii<xDim;ii+=1)
1285	for(jj=0;jj<yDim;jj+=1)
1286	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1287	qVal = qTotal[ii][jj]
1288	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1289	val = inten[ii][jj]
1290
1291	// if(isVCALC \|\| maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1292	if(maskMissing) // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1293	mask_val = 0
1294	else
1295	mask_val = mask[ii][jj]
1296	endif
1297	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1298	iBin_qxqy[binIndex] += val
1299	iBin2_qxqy[binIndex] += val*val
1300	eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1301	nBin_qxqy[binIndex] += 1
1302	endif
1303	endfor
1304	endfor
1305
1306	endif
1307
1308	// add in set 2 (set 1 already done)
1309	if(nSets >= 2)
1310	xDim=DimSize(inten2,0)
1311	yDim=DimSize(inten2,1)
1312
1313	for(ii=0;ii<xDim;ii+=1)
1314	for(jj=0;jj<yDim;jj+=1)
1315	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1316	qVal = qTotal2[ii][jj]
1317	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1318	val = inten2[ii][jj]
1319
1320	// if(isVCALC \|\| maskMissing)
1321	if(maskMissing)
1322	mask_val = 0
1323	else
1324	mask_val = mask2[ii][jj]
1325	endif
1326	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1327	iBin_qxqy[binIndex] += val
1328	iBin2_qxqy[binIndex] += val*val
1329	eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1330	nBin_qxqy[binIndex] += 1
1331	endif
1332	endfor
1333	endfor
1334
1335	endif
1336
1337	// add in set 3 and 4 (set 1 and 2 already done)
1338	if(nSets == 4)
1339	xDim=DimSize(inten3,0)
1340	yDim=DimSize(inten3,1)
1341
1342	for(ii=0;ii<xDim;ii+=1)
1343	for(jj=0;jj<yDim;jj+=1)
1344	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1345	qVal = qTotal3[ii][jj]
1346	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1347	val = inten3[ii][jj]
1348
1349	// if(isVCALC \|\| maskMissing)
1350	if(maskMissing)
1351	mask_val = 0
1352	else
1353	mask_val = mask3[ii][jj]
1354	endif
1355	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1356	iBin_qxqy[binIndex] += val
1357	iBin2_qxqy[binIndex] += val*val
1358	eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1359	nBin_qxqy[binIndex] += 1
1360	endif
1361	endfor
1362	endfor
1363
1364
1365	xDim=DimSize(inten4,0)
1366	yDim=DimSize(inten4,1)
1367
1368	for(ii=0;ii<xDim;ii+=1)
1369	for(jj=0;jj<yDim;jj+=1)
1370	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1371	qVal = qTotal4[ii][jj]
1372	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1373	val = inten4[ii][jj]
1374
1375	// if(isVCALC \|\| maskMissing)
1376	if(maskMissing)
1377	mask_val = 0
1378	else
1379	mask_val = mask4[ii][jj]
1380	endif
1381	if (numType(val)==0 && mask_val == 0) //count only the good points, ignore Nan or Inf
1382	iBin_qxqy[binIndex] += val
1383	iBin2_qxqy[binIndex] += val*val
1384	eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1385	nBin_qxqy[binIndex] += 1
1386	endif
1387	endfor
1388	endfor
1389
1390	endif
1391
1392
1393	// after looping through all of the data on the panels, calculate errors on I(q),
1394	// just like in CircSectAve.ipf
1395	// TODO:
1396	// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
1397	// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1398	for(ii=0;ii<nq;ii+=1)
1399	if(nBin_qxqy[ii] == 0)
1400	//no pixels in annuli, data unknown
1401	iBin_qxqy[ii] = 0
1402	eBin_qxqy[ii] = 1
1403	eBin2D_qxqy[ii] = NaN
1404	else
1405	if(nBin_qxqy[ii] <= 1)
1406	//need more than one pixel to determine error
1407	iBin_qxqy[ii] /= nBin_qxqy[ii]
1408	eBin_qxqy[ii] = 1
1409	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1410	else
1411	//assume that the intensity in each pixel in annuli is normally distributed about mean...
1412	// -- this is correctly calculating the error as the standard error of the mean, as
1413	// was always done for SANS as well.
1414	iBin_qxqy[ii] /= nBin_qxqy[ii]
1415	avesq = iBin_qxqy[ii]^2
1416	aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1417	var = aveisq-avesq
1418	if(var<=0)
1419	eBin_qxqy[ii] = 1e-6
1420	else
1421	eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1422	endif
1423	// and calculate as it is propagated pixel-by-pixel
1424	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1425	endif
1426	endif
1427	endfor
1428
1429	eBin2D_qxqy = sqrt(eBin2D_qxqy) // as equation (3) of John's memo
1430
1431	// find the last non-zero point, working backwards
1432	val=nq
1433	do
1434	val -= 1
1435	while((nBin_qxqy[val] == 0) && val > 0)
1436
1437	// print val, nBin_qxqy[val]
1438	DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1439
1440	if(val == 0)
1441	// all the points were deleted
1442	return(0)
1443	endif
1444
1445
1446	// since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1447	// find the first non-zero point, working forwards
1448	val = -1
1449	do
1450	val += 1
1451	while(nBin_qxqy[val] == 0)
1452	DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1453
1454	// ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1455	val = numpnts(nBin_qxqy)-1
1456	do
1457	if(nBin_qxqy[val] == 0)
1458	DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1459	endif
1460	val -= 1
1461	while(val>0)
1462
1463	// utility function to remove NaN values from the waves
1464	V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
1465
1466
1467	// TODO:
1468	// -- This is where I calculate the resolution in SANS (see CircSectAve)
1469	// -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1470	// -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1471	//
1472	nq = numpnts(qBin_qxqy)
1473	Make/O/D/N=(nq) $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1474	Make/O/D/N=(nq) $(folderPath+":"+"qBar_qxqy"+"_"+type)
1475	Make/O/D/N=(nq) $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1476	Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1477	Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1478	Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1479	Variable ret1,ret2,ret3
1480
1481	// all of the different collimation conditions are handled within the V_getResolution function
1482	// which is responsible for switching based on the different collimation types (white beam, slit, Xtal, etc)
1483	// to calculate the correct resolution, or fill the waves with the correct "flags"
1484	//
1485
1486	// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
1487	// geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
1488	// as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
1489	// with the inner(lambda) calculated first, then the outer(geometry).
1490	//
1491
1492	// possible values are:
1493	//
1494	// pinhole
1495	// pinhole_whiteBeam
1496	// convergingPinholes
1497	//
1498	// *slit data should be reduced using the slit routine, not here, proceed but warn
1499	// narrowSlit
1500	// narrowSlit_whiteBeam
1501	//
1502	ii=0
1503	do
1504	V_getResolution(qBin_qxqy[ii],folderStr,type,collimationStr,ret1,ret2,ret3)
1505	sigmaq[ii] = ret1
1506	qbar[ii] = ret2
1507	fsubs[ii] = ret3
1508	ii+=1
1509	while(ii<nq)
1510
1511
1512	SetDataFolder root:
1513
1514	return(0)
1515	End
1516
1517

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

Download in other formats:

Original Format