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