PEARL Procedures rev-distro-1.6.0-0-gcf1399e-dirty Igor procedures for the analysis of PEARL data

- Main Page
- Related Pages
- Packages
- Namespaces
- Files

- File List
- Globals

pearl-scienta-preprocess.ipf
Go to the documentation of this file.
1 #pragma rtGlobals=3// Use modern global access method and strict wave access.
2 #pragma IgorVersion = 6.1
3 #pragma ModuleName = PearlScientaPreprocess
4 #pragma version = 1.03
5 #include "mm-fitfuncs"
6
7 // $Id$
8 // author: matthias.muntwiler@psi.ch
9 // Copyright (c) 2013-17 Paul Scherrer Institut
10
11 // Licensed under the Apache License, Version 2.0 (the "License");
12 // you may not use this file except in compliance with the License.
13 // You may obtain a copy of the License at
14 // http://www.apache.org/licenses/LICENSE-2.0
15
30
35
36 variable prompt_int_linbg_reduction(string* param){
37 string &param
38
39 variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
40 variable Lsize = NumberByKey("Lsize", param, "=", ";")
41 variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
42 variable Hsize = NumberByKey("Hsize", param, "=", ";")
43 variable Cpos = NumberByKey("Cpos", param, "=", ";")
44 variable Csize = NumberByKey("Csize", param, "=", ";")
45
46 prompt Lcrop, "Lower cropping region"
47 prompt Hcrop, "Upper cropping region"
48 prompt Lsize, "Lower background region"
49 prompt Hsize, "Upper background region"
50 prompt Cpos, "Center position"
51 prompt Csize, "Center integration region"
52
53 doprompt "int_linbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
54 if (v_flag == 0)
55 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
56 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
57 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
58 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
59 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
60 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
61 endif
62
63 return v_flag
64 };
65
66 string capture_int_linbg_cursors(){
67 // this function is for testing only, until we implement a proper interface
68 string param = csr_int_linbg_reduction("")
69 svar /z global_params = root:packages:pearl_explorer:s_reduction_params
70 if (svar_exists(global_params))
71 global_params = param
72 endif
73 return param
74 };
75
76 string csr_int_linbg_reduction(string win){
77 // PRELIMINARY - function arguments may change
78
79 // sets reduction parameters from cursors in a graph.
80 // an even number of cursors (2 or more) must be set on the image.
81 // cursor names and order do not matter,
82 // except that the alphabetically first cursor which is attached to an image selects the image.
83 // the cursors mark the following positions, from innermost to outermost pair:
84 // 1) low and high limits of peak region.
85 // 2) peak-side boundary of lower and upper background region.
86 // 3) lower and upper cropping region.
87
88 string win
89
90 // read all cursor positions
91 variable ic
92 string sc
93 variable nc = 10
94 make /n=(nc) /free positions
95 variable np = 0
96 wave /z image = $""
97 string imagename = ""
98 string tracename = ""
99 string info
100
101 for (ic = 0; ic < nc; ic += 1)
102 sc = num2char(char2num("A") + ic)
103 wave /z wc = CsrWaveRef($sc, win)
104 info = CsrInfo($sc, win)
105 tracename = StringByKey("TNAME", info, ":", ";")
106 if (waveexists(wc) && (wavedims(wc) == 2))
107 if (!waveexists(image))
108 wave image = wc
109 imagename = tracename
110 endif
111 if (cmpstr(tracename, imagename) == 0)
112 positions[np] = pcsr($sc, win)
113 np += 1
114 endif
115 endif
116 endfor
117
118 np = floor(np / 2) * 2// ignore odd cursor
119 redimension /n=(np) positions
120 sort positions, positions
121 // shift upper positions by one so that the rightmost pixel becomes 1.0
122 positions = p >= np / 2 ? positions + 1 : positions
123 positions = positions / dimsize(image, 0)
124
125 // map innermost cursor pair to peak center and size
126 variable ip2 = np / 2
127 variable ip1 = ip2 - 1
128 variable Cpos = (positions[ip1] + positions[ip2]) / 2
129 variable Csize = positions[ip2] - positions[ip1]
130 if (ip1 >= 0)
131 Cpos = (positions[ip1] + positions[ip2]) / 2
132 Csize = positions[ip2] - positions[ip1]
133 else
134 // default: a small region in the center
135 Cpos = 0.5
136 Csize = 0.2
137 endif
138
139 // background region
140 ip1 -= 1
141 ip2 += 1
142 variable Lsize
143 variable Hsize
144 if (ip1 >= 0)
145 Lsize = positions[ip1]
146 Hsize = 1 - positions[ip2]
147 else
148 // default: everything outside the peak region
149 Lsize = Cpos - Csize / 2
150 Hsize = 1 - (Cpos + Csize / 2)
151 endif
152
153 // crop region
154 ip1 -= 1
155 ip2 += 1
156 variable Lcrop
157 variable Hcrop
158 if (ip1 >= 0)
159 Lcrop = positions[ip1]
160 Hcrop = 1 - positions[ip2]
161 else
162 // default: dead corners of the EW4000 at PEARL
163 Lcrop = 0.11
164 Hcrop = 0.11
165 endif
166 Lsize = max(Lsize - Lcrop, 0)
167 Hsize = max(Hsize - Hcrop, 0)
168
169 string param = ""
170 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
171 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
172 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
173 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
174 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
175 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
176
177 return param
178 };
179
180 variable test_int_linbg(wave image){
181 // useful for testing or manual processing
182 // since the int_linbg_reduction function cannot be called from the command line directly.
183 wave image
184
185 string param = ""
186 prompt_int_linbg_reduction(param)
187
188 string data1_name = "test_data1"
189 string data2_name = "test_data2"
190 duplicate /o image, $data1_name, $data2_name
191 wave w_data1 = $data1_name
192 wave w_data2 = $data2_name
193
194 int_linbg_reduction(image, w_data1, w_data2, param)
195 };
196
197 threadsafe variable int_linbg_reduction(wave source, wave dest1, wave dest2, string* param){
198 // data reduction function for adh5_load_reduced_detector
199 // calculates the average pixel value of each angular slice
200 // in one center and two background intervals.
201 // a background value is calculated at the center position
202 // by linear interpolation from the two background values.
203 // returns the center minus linear background in dest1.
204 // returns the Poisson one-sigma error in dest2.
205 wave source// source wave
206 // Scienta detector image, energy axis along X, angle axis along Y
207 wave dest1, dest2// destination waves
208 // each wave is a one-dimensional intensity distribution
209 // the function may redimension these waves to one of the image dimensions
210 // (it must be clear to the user which dimension this is).
211 // the meaning of dest1 and dest2 is up to the particular function,
212 // e.g. dest1 could hold the mean value and dest2 the one-sigma error,
213 // or dest1 could hold the X-profile, and dest2 the Y-profile.
214 string &param// parameters in a key1=value1;key2=value2;... list
215 // all region parameters are relative to the image size (0...1)
216 // Lcrop = size of the lower cropping region
217 // Hcrop = size of the upper cropping region
218 // Lsize = size of the lower background integration region
219 // Hsize = size of the upper background integration region
220 // Cpos = center position of the of the peak integration region
221 // Csize = size of the peak integration region
222
223 // typical values (peak centered on detector, FWHM ~ 20 % of image)
224 // Lcrop=0.11;Hcrop=0.11;Lsize=0.2;Hsize=0.2;Cpos=0.5;Csize=0.2
225
226 variable nx = dimsize(source, 0)
227 variable ny = dimsize(source, 1)
228
229 // read parameters
230 variable lcrop = NumberByKey("Lcrop", param, "=", ";")
231 variable lsize = NumberByKey("Lsize", param, "=", ";")
232 variable hcrop = NumberByKey("Hcrop", param, "=", ";")
233 variable hsize = NumberByKey("Hsize", param, "=", ";")
234 variable cpos = NumberByKey("Cpos", param, "=", ";")
235 variable csize = NumberByKey("Csize", param, "=", ";")
236
237 // validate parameters
238 // background parameters are optional, center parameter is required.
239 if (numtype(lcrop) != 0)
240 lcrop = 0
241 endif
242 if (numtype(lsize) != 0)
243 lsize = 0
244 endif
245 if (numtype(hcrop) != 0)
246 hcrop = 0
247 endif
248 if (numtype(hsize) != 0)
249 hsize = 0
250 endif
251 if (numtype(Cpos) != 0)
252 return 1// Cpos parameter missing
253 endif
254 if (numtype(Csize) != 0)
255 return 2// Csize parameter missing
256 endif
257
258 variable lpos = lcrop + lsize / 2
259 variable hpos = 1 - (hcrop + hsize / 2)
260
261 variable p0
262 variable p1
263
264 adh5_setup_profile(source, dest1, 1)
265 adh5_setup_profile(source, dest2, 1)
266
267 duplicate /free dest1, lbg, hbg
268 if (lsize > 0)
269 p0 = round(lcrop * nx)
270 p1 = round((lcrop + lsize) * nx)
271 ad_profile_y_w(source, p0, p1, lbg)
272 else
273 lbg = 0
274 endif
275 if (hsize > 0)
276 p0 = round((1 - hcrop - hsize) * nx)
277 p1 = round((1 - hcrop) * nx)
278 ad_profile_y_w(source, p0, p1, hbg)
279 else
280 hbg = 0
281 endif
282 if (csize > 0)
283 p0 = round((cpos - csize/2) * nx)
284 p1 = round((cpos + csize/2) * nx)
285 ad_profile_y_w(source, p0, p1, dest1)
286 else
287 dest1 = 0
288 endif
289
290 variable scale = (cpos - lpos) / (hpos - lpos)
291 dest2 = dest1
292 dest1 -= scale * (hbg - lbg) + lbg
293 dest2 = sqrt(dest2 + scale^2 * (hbg + lbg))
294 return 0// return zero if successful, non-zero if an error occurs
295 };
296
297 variable test_shockley_anglefit(wave image, variable branch){
298 // apply the Shockley_anglefit function to a single image
299 // useful for testing or manual processing
300 // since the Shockley_anglefit function cannot be called from the command line directly.
301 wave image
302 variable branch// +1 or -1
303
304 string param = ""
305 param = ReplaceStringByKey("branch", param, num2str(branch), "=", ";")
306
307 string s_branch
308 if (branch >= 0)
309 s_branch = "p"
310 else
311 s_branch = "n"
312 endif
313 string pkpos_name = "saf_pkpos_" + s_branch
314 string pkwid_name = "saf_pkwid_" + s_branch
315 duplicate /o image, $pkpos_name, $pkwid_name
316 wave w_pkpos = $pkpos_name
317 wave w_pkwid = $pkwid_name
318
319 shockley_anglefit(image, w_pkpos, w_pkwid, param)
320 };
321
322 variable prompt_Shockley_anglefit(string* param){
323 string &param
324
325 variable branch = NumberByKey("branch", param, "=", ";")
326
327 prompt branch, "Branch (-1 or +1)"
328
329 doprompt "Shockley_anglefit_reduction Parameters", branch
330 if (v_flag == 0)
331 param = ReplaceNumberByKey("branch", param, branch, "=", ";")
332 endif
333
334 return v_flag
335 };
336
337 threadsafe variable Shockley_anglefit(wave source, wave dest1, wave dest2, string* param){
338 // data reduction function for adh5_load_reduced_detector
339 // specialized for analysing the Cu(111) Shockley surface state
340 // do curve fitting of one branch of the surface state
341 // the result is peak position versus energy
342 // TODO: this function contains hard-coded parameters. please generalize as necessary.
343 wave source// source wave
344 // Scienta detector image, energy axis along X, angle axis along Y
345 // the apex of the surface state must be at angle 0
346 wave dest1, dest2// destination waves
347 // dest1: peak position
348 // dest2: peak width (sigma)
349 string &param// parameters in a key1=value1;key2=value2;... list
350 // branch=-1 or +1: select negative (positive) angles for the fit interval, respectively
351
352 variable nx = dimsize(source, 0)
353 variable ny = dimsize(source, 1)
354
355 // read parameters
356 variable branch = NumberByKey("branch", param, "=", ";")
357
358 // validate parameters
359 if (numtype(branch) != 0)
360 branch = +1
361 endif
362
363 // prepare output
364 adh5_setup_profile(source, dest1, 0)
365 adh5_setup_profile(source, dest2, 0)
366 dest1 = nan
367 dest2 = nan
368
369 // select angle range
370 // hard-coded for a particular measurement series
371 variable y0
372 variable y1
373 if (branch < 0)
374 y0 = -5
375 y1 = 0
376 else
377 y0 = 0
378 y1 = 5
379 endif
380
381 // select energy range
382 // start at the point of highest intensity and go up 0.45 eV
383 variable p0
384 variable p1
385 variable q0
386 variable q1
387 duplicate /free dest1, center
388 q0 = round((y0 - dimoffset(source, 1)) / dimdelta(source, 1))
389 q1 = round((y1 - dimoffset(source, 1)) / dimdelta(source, 1))
390 ad_profile_x_w(source, q0, q1, center)
391 wavestats /q/m=1 center
392 p0 = round((v_maxloc - dimoffset(source, 0)) / dimdelta(source, 0))
393 p1 = round((v_maxloc + 0.4 - dimoffset(source, 0)) / dimdelta(source, 0))
394
395 // prepare intermediate data buffer
396 make /n=(ny)/d/free profile
397 setscale /p x dimoffset(source,1), dimdelta(source,1), waveunits(source,1), profile
398
399 variable pp
400 for (pp = p0; pp <= p1; pp += 1)
401 profile = source[pp][p]
402 curvefit /Q /NTHR=1 /W=2 gauss profile(y0,y1)
403 wave w_coef
404 dest1[pp] = w_coef[2]
405 dest2[pp] = w_coef[3]
406 endfor
407 return 0// return zero if successful, non-zero if an error occurs
408 };
409
410 variable prompt_int_quadbg_reduction(string* param){
411 string &param
412
413 variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
414 variable Lsize = NumberByKey("Lsize", param, "=", ";")
415 variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
416 variable Hsize = NumberByKey("Hsize", param, "=", ";")
417 variable Cpos = NumberByKey("Cpos", param, "=", ";")
418 variable Csize = NumberByKey("Csize", param, "=", ";")
419
420 prompt Lcrop, "Lower cropping region"
421 prompt Hcrop, "Upper cropping region"
422 prompt Lsize, "Lower background region"
423 prompt Hsize, "Upper background region"
424 prompt Cpos, "Center position"
425 prompt Csize, "Center integration region"
426
427 doprompt "int_quadbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
428 if (v_flag == 0)
429 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
430 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
431 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
432 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
433 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
434 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
435 endif
436
437 return v_flag
438 };
439
440 variable test_int_quadbg(wave image){
441 // useful for testing or manual processing
442 // since the int_quadbg_reduction function cannot be called from the command line directly.
443 wave image
444
445 string param = ""
446 prompt_int_quadbg_reduction(param)
447
448 string data1_name = "test_data1"
449 string data2_name = "test_data2"
450 duplicate /o image, $data1_name, $data2_name
451 wave w_data1 = $data1_name
452 wave w_data2 = $data2_name
453
454 int_quadbg_reduction(image, w_data1, w_data2, param)
455 };
456
457 threadsafe variable int_quadbg_reduction(wave source, wave dest1, wave dest2, string* param){
458 // data reduction function for adh5_load_reduced_detector
459 // integrates peak area minus a quadratic backgrouind
460 wave source// source wave
461 // Scienta detector image, energy axis along X, angle axis along Y
462 wave dest1, dest2// destination waves
463 string &param// parameters in a key1=value1;key2=value2;... list
464 // all region parameters are relative to the image size (0...1)
465 // Lcrop = size of the lower cropping region
466 // Hcrop = size of the upper cropping region
467 // Lsize = size of the lower background integration region
468 // Hsize = size of the upper background integration region
469 // Cpos = center position of the of the peak integration region
470 // Csize = size of the peak integration region
471
472 // typical values (peak centered on detector, FWHM ~ 20 % of image)
473 // Lcrop=0.11;Hcrop=0.11;Lsize=0.2;Hsize=0.2;Cpos=0.5;Csize=0.2
474
475 variable nx = dimsize(source, 0)
476 variable ny = dimsize(source, 1)
477
478 // read parameters
479 variable lcrop = NumberByKey("Lcrop", param, "=", ";")
480 variable lsize = NumberByKey("Lsize", param, "=", ";")
481 variable hcrop = NumberByKey("Hcrop", param, "=", ";")
482 variable hsize = NumberByKey("Hsize", param, "=", ";")
483 variable cpos = NumberByKey("Cpos", param, "=", ";")
484 variable csize = NumberByKey("Csize", param, "=", ";")
485
486 // validate parameters
487 // background parameters are optional, center parameter is required.
488 if (numtype(lcrop) != 0)
489 lcrop = 0
490 endif
491 if (numtype(lsize) != 0)
492 lsize = 0
493 endif
494 if (numtype(hcrop) != 0)
495 hcrop = 0
496 endif
497 if (numtype(hsize) != 0)
498 hsize = 0
499 endif
500 if (numtype(Cpos) != 0)
501 return 1// Cpos parameter missing
502 endif
503 if (numtype(Csize) != 0)
504 return 2// Csize parameter missing
505 endif
506
507 // crop boundaries
508 variable pcl = round(lcrop * nx)
509 variable pch = round((1 - hcrop) * nx)
510 // fit boundaries
511 variable pfl = round((lcrop + lsize) * nx)
512 variable pfh = round((1 - hcrop - hsize) * nx)
513 // integration boundaries
514 variable pil = round((cpos - csize/2) * nx)
515 variable pih = round((cpos + csize/2) * nx)
516
517 adh5_setup_profile(source, dest1, 0)
518 adh5_setup_profile(source, dest2, 0)
519
520 // prepare intermediate data buffer
521 make /n=(nx) /d /free profile, mask, fit
522 setscale /p x dimoffset(source,0), dimdelta(source,0), waveunits(source,0), profile, mask, fit
523 mask = ((p >= pcl) && (p < pfl)) || ((p >= pfh) && (p < pch))
524
525 variable qq
526 variable sp, sf
527 variable xil = x2pnt(profile, pil)
528 variable xih = x2pnt(profile, pih)
529
530 make /n=3 /free /d w_coef
531 for (qq = 0; qq < ny; qq += 1)
532 profile = source[p][qq]
533 curvefit /Q /NTHR=1 /W=2 poly 3, kwCWave=w_coef, profile /M=mask
534 fit = poly(w_coef, x)
535 sp = sum(profile, xil, xih)
536 sf = sum(fit, xil, xih)
537 dest1[qq] = sp - sf
538 dest2[qq] = sqrt(sp)
539 endfor
540
541 return 0// return zero if successful, non-zero if an error occurs
542 };
543
544 variable scienta_norm(wave w, variable x){
545 wave w
546 variable x
547
548 return w[0] * (x^2 - w[1]^2)
549 };
550
551 wave fit_scienta_ang_transm(wave data, wave params){
552 wave data// measured angular distribution (1D)
553 wave /z params
554
555 if (!waveexists(params))
556 make /n=12 /o params
557 endif
558 redimension /n=12/d params
559
560 variable h = wavemax(data) - wavemin(data)
561 params[0] = h / 2
562 params[1] = 0
563 params[2] = 7
564 params[3] = h / 4
565 params[4] = -23
566 params[5] = 4
567 params[6] = h / 4
568 params[7] = +23
569 params[8] = 4
570 params[9] = h / 2
571 params[10] = 0
572 params[11] = -0.001
573 FuncFit /NTHR=0 /q scienta_ang_transm params data
574
575 return params
576 };
577
578 threadsafe variable scienta_ang_transm(wave w, variable x){
579 // parameterized angular transmission function of the analyser
580 wave w// coefficients
581 // w[0] = amplitude gauss 1
582 // w[1] = position gauss 1
583 // w[2] = width gauss 1
584 // w[3] = amplitude gauss 2
585 // w[4] = position gauss 2
586 // w[5] = width gauss 2
587 // w[6] = amplitude gauss 3
588 // w[7] = position gauss 3
589 // w[8] = width gauss 3
590 // w[9] = constant background
591 // w[10] = linear background
592 // w[11] = quadratic background
593 variable x
594
595 make /free /n=4 /d w_int
596 w_int[0] = 0
597 w_int[1,] = w[p - 1]
598 variable pk1 = gauss1d(w_int, x)
599 w_int[1,] = w[p + 2]
600 variable pk2 = gauss1d(w_int, x)
601 w_int[1,] = w[p + 5]
602 variable pk3 = gauss1d(w_int, x)
603 w_int[0,2] = w[9 + p]
604 w_int[3] = 0
605 variable bg = poly(w_int, x)
606
607 return bg + pk1 + pk2 + pk3
608 };
609
610 wave fit_scienta_poly_bg(wave data, wave params, variable bgterms){
611 wave data// measured distribution (2D)
612 wave /z params// wave, will be redimensioned for the correct size
613 variable bgterms// number of terms in the polynomial background: 2, 3, or 4
614
615 if (!waveexists(params))
616 make /n=15 /o params
617 endif
618 redimension /n=15 /d params
619
620 variable wmax = wavemax(data)
621 variable wmin = wavemin(data)
622 params[0] = 0
623 params[1] = 7
624 params[2] = 1 / 2
625 params[3] = -23
626 params[4] = 4
627 params[5] = 1 / 2
628 params[6] = +23
629 params[7] = 4
630 params[8] = 1
631 params[9] = 0
632 params[10] = -0.001
633 params[11] = wmin
634 params[12] = (wmax - wmin) / dimdelta(data,1) / dimsize(data,1)
635 params[13] = 0
636 params[14] = 0
637
638 string h = "0000000000000"
639 if (bgterms < 3)
640 h = h + "1"
641 else
642 h = h + "0"
643 endif
644 if (bgterms < 4)
645 h = h + "1"
646 else
647 h = h + "0"
648 endif
649 FuncFitMD /NTHR=1 /q /h=h scienta_poly_bg params data
650
651 return params
652 };
653
654 variable scienta_poly_bg(wave w, variable e, variable a){
655 // polynomial background with
656 // parameterized angular transmission function of the analyser
657 wave w// coefficients
658 // angular transmission, varies with a
659 // amplitude of gauss 1 = 1 constant
660 // other peak amplitudes and linear terms are relative to gauss 1
661 // w[0] = position gauss 1
662 // w[1] = width gauss 1
663 // w[2] = amplitude gauss 2, relative to gauss 1
664 // w[3] = position gauss 2
665 // w[4] = width gauss 2
666 // w[5] = amplitude gauss 3, relative to gauss 1
667 // w[6] = position gauss 3
668 // w[7] = width gauss 3
669 // w[8] = constant term
670 // w[9] = linear term
671 // w[10] = quadratic term
672 // spectral background, varies with e
673 // w[11] = constant term
674 // w[12] = linear term
675 // w[13] = quadratic term
676 // w[14] = cubic term
677 variable a// detection angle
678 variable e// electron energy
679
680 make /free /n=4 /d w_int
681 variable p0 = 0
682
683 w_int[0] = 0
684 w_int[1] = 1
685 w_int[2,] = w[p0 + p - 2]
686 variable pk1 = gauss1d(w_int, a)
687 p0 += 2
688
689 w_int[1,] = w[p0 + p - 1]
690 variable pk2 = gauss1d(w_int, a)
691 p0 += 3
692
693 w_int[1,] = w[p0 + p - 1]
694 variable pk3 = gauss1d(w_int, a)
695 p0 += 3
696
697 w_int[0,2] = w[p0 + p]
698 w_int[3] = 0
699 variable base = poly(w_int, a)
700 p0 += 3
701
702 w_int[0,3] = w[p0 + p]
703 variable bg = poly(w_int, e)
704
705 return bg * (base + pk1 + pk2 + pk3)
706 };
707
716 variable prompt_redim_linbg_reduction(string* param){
717 string &param
718
719 variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
720 variable Lsize = NumberByKey("Lsize", param, "=", ";")
721 variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
722 variable Hsize = NumberByKey("Hsize", param, "=", ";")
723 variable Cpos = NumberByKey("Cpos", param, "=", ";")
724 variable Csize = NumberByKey("Csize", param, "=", ";")
725
726 prompt Lcrop, "Lower cropping region"
727 prompt Hcrop, "Upper cropping region"
728 prompt Lsize, "Lower background region"
729 prompt Hsize, "Upper background region"
730 prompt Cpos, "Center position"
731 prompt Csize, "Center integration region"
732
733 doprompt "redim_linbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
734 if (v_flag == 0)
735 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
736 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
737 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
738 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
739 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
740 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
741 endif
742
743 return v_flag
744 };
745
787 threadsafe variable redim_linbg_reduction(wave source, wave dest1, wave dest2, string* param){
788 wave source
789 wave dest1, dest2
790 string &param
791
792 variable nx = dimsize(source, 0)
793 variable ny = dimsize(source, 1)
794
795 duplicate /free source, source_redim
796 redimension /n=(nx * ny) source_redim
797 nx += 1
798 redimension /n=(nx, ny) source_redim
799
800 return int_linbg_reduction(source_redim, dest1, dest2, param)
801 };
802
803 // DblGaussLinBG function fit reduction
804
810 variable test_gauss2_reduction(wave image){
811 wave image
812
813 string param = ""
814 param = ReplaceNumberByKey("rngl", param, -inf, "=", ";")
815 param = ReplaceNumberByKey("rngh", param, inf, "=", ";")
816 param = ReplaceNumberByKey("pos1", param, 113.70, "=", ";")
817 param = ReplaceNumberByKey("wid1", param, 1.46, "=", ";")
818 param = ReplaceNumberByKey("pos2", param, 126.53, "=", ";")
819 param = ReplaceNumberByKey("wid2", param, 1.63, "=", ";")
820 param = ReplaceNumberByKey("ybox", param, 1, "=", ";")
821 param = ReplaceStringByKey("return", param, "amp1", "=", ";")
822
823 make /o test1
824 make /o test2
825
826 gauss2_reduction(image, test1, test2, param)
827 };
828
834 variable prompt_gauss2_reduction(string* param){
835 string &param
836
837 variable rngl = NumberByKey("rngl", param, "=", ";")
838 variable rngh = NumberByKey("rngh", param, "=", ";")
839 variable pos1 = NumberByKey("pos1", param, "=", ";")
840 variable wid1 = NumberByKey("wid1", param, "=", ";")
841 variable pos2 = NumberByKey("pos2", param, "=", ";")
842 variable wid2 = NumberByKey("wid2", param, "=", ";")
843 variable ybox = NumberByKey("ybox", param, "=", ";")
844 string retpar = StringByKey("return", param, "=", ";")
845
846 prompt rngl, "range low"
847 prompt rngh, "range high"
848 prompt pos1, "position 1"
849 prompt wid1, "width 1"
850 prompt pos2, "position 2"
851 prompt wid2, "width 2"
852 prompt ybox, "ybox (1 or 3)"
853 prompt retpar, "return", popup "amp1;amp2;"
854
855 doprompt "gauss2_reduction reduction parameters", rngl, rngh, pos1, wid1, pos2, wid2, ybox, retpar
856 if (v_flag == 0)
857 param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
858 param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
859 param = ReplaceNumberByKey("pos1", param, pos1, "=", ";")
860 param = ReplaceNumberByKey("wid1", param, wid1, "=", ";")
861 param = ReplaceNumberByKey("pos2", param, pos2, "=", ";")
862 param = ReplaceNumberByKey("wid2", param, wid2, "=", ";")
863 param = ReplaceStringByKey("return", param, retpar, "=", ";")
864 param = ReplaceNumberByKey("ybox", param, ybox, "=", ";")
865 endif
866
867 return v_flag
868 };
869
903 threadsafe variable gauss2_reduction(wave source, wave dest1, wave dest2, string* param){
904 wave source
905 wave dest1, dest2
906 string &param
907
908 variable nx = dimsize(source, 0)
909 variable ny = dimsize(source, 1)
910
911 // prepare output
912 adh5_setup_profile(source, dest1, 1)
913 adh5_setup_profile(source, dest2, 1)
914 dest1 = 0
915 dest2 = 0
916
917 // read parameters
918 variable rngl = NumberByKey("rngl", param, "=", ";")
919 variable rngh = NumberByKey("rngh", param, "=", ";")
920 variable pos1 = NumberByKey("pos1", param, "=", ";")
921 variable wid1 = NumberByKey("wid1", param, "=", ";")
922 variable pos2 = NumberByKey("pos2", param, "=", ";")
923 variable wid2 = NumberByKey("wid2", param, "=", ";")
924 variable ybox = NumberByKey("ybox", param, "=", ";")
925 string retpar = StringByKey("return", param, "=", ";")
926
927 variable idestcoef
928 if (cmpstr(retpar, "amp2") == 0)
929 idestcoef = 3
930 else
931 idestcoef = 0
932 endif
933
934 // prepare curve fit
935 make /free xprof
936 adh5_setup_profile(source, xprof, 0)
937 duplicate /free xprof, xprof_sig
938 make /free /d /n=8 w_coef, W_sigma
939 w_coef[0] = {1, pos1, wid1, 1, pos2, wid2, 0, 0}
940
941 variable pl = max(x2pnt(xprof, rngl), 0)
942 variable ph = min(x2pnt(xprof, rngh), numpnts(xprof) - 1)
943
944 // text constraints cannot be used in threadsafe functions.
945 // the following matrix-vector forumlation is equivalent to:
946 // make /free /T /N=4 constraints
947 // constraints[0] = {"K0 >= 0", "K3 >= 0", "K6 <= 0", "K7 => 0"}
948 make /free /n=(4,8) cmat
949 make /free /n=4 cvec
950 cmat = 0
951 cmat[0][0] = -1
952 cmat[1][3] = -1
953 cmat[2][6] = 1
954 cmat[3][7] = -1
955 cvec = 0
956
957 // loop over angle scale
958 variable p0 = 0
959 variable p1 = numpnts(dest1) - 1
960 variable pp
961 variable wmin
962 variable wmax
963 if (ybox > 1)
964 p0 += ceil((ybox - 1) / 2)
965 p1 -= ceil((ybox - 1) / 2)
966 endif
967 variable V_FitNumIters
968
969 for (pp = p0; pp <= p1; pp += 1)
970 xprof = source[p][pp]
971 if (ybox > 1)
972 xprof += source[p][pp-1] + source[p][pp+1]
973 endif
974 xprof_sig = max(sqrt(xprof), 1)
975 xprof /= ybox
976 xprof_sig /= ybox
977
978 wmin = wavemin(xprof)
979 wmax = wavemax(xprof)
980 w_coef[0] = wmax - wmin
981 w_coef[3] = w_coef[0]
982 w_coef[6] = 0
983 w_coef[7] = wmin
984 FuncFit /H="01101100" /Q /NTHR=1 /N /W=2 DblGaussLinBG w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
985 wave w_sigma
986
987 if (V_FitNumIters < 40)
988 dest1[pp] = max(w_coef[idestcoef], 0)
989 dest2[pp] = max(w_sigma[idestcoef], 0)
990 endif
991 endfor
992
993 dest1 *= w_coef[idestcoef+2] * sqrt(pi)
994 dest2 *= w_coef[idestcoef+2] * sqrt(pi)
995
996 return 0
997 };
998
fit_scienta_ang_transm
wave fit_scienta_ang_transm(wave data, wave params)
Definition: pearl-scienta-preprocess.ipf:551
test_int_linbg
variable test_int_linbg(wave image)
Definition: pearl-scienta-preprocess.ipf:180
Shockley_anglefit
threadsafe variable Shockley_anglefit(wave source, wave dest1, wave dest2, string *param)
Definition: pearl-scienta-preprocess.ipf:337
prompt_gauss2_reduction
variable prompt_gauss2_reduction(string *param)
prompt for the gauss2_reduction parameters
Definition: pearl-scienta-preprocess.ipf:834
redim_linbg_reduction
threadsafe variable redim_linbg_reduction(wave source, wave dest1, wave dest2, string *param)
linear background reduction function for incorrectly dimensioned scienta image
Definition: pearl-scienta-preprocess.ipf:787
scienta_poly_bg
variable scienta_poly_bg(wave w, variable e, variable a)
Definition: pearl-scienta-preprocess.ipf:654
ad_profile_x_w
threadsafe wave ad_profile_x_w(wave dataset, variable q1, variable q2, wave destwave, variable noavg=defaultValue)
1D cut through 2D dataset along X dimension, existing destination wave.
Definition: pearl-area-profiles.ipf:504
fit_scienta_poly_bg
wave fit_scienta_poly_bg(wave data, wave params, variable bgterms)
Definition: pearl-scienta-preprocess.ipf:610
prompt_Shockley_anglefit
variable prompt_Shockley_anglefit(string *param)
Definition: pearl-scienta-preprocess.ipf:322
adh5_setup_profile
threadsafe variable adh5_setup_profile(wave image, wave profile, variable dim)
set up a one-dimensional wave for a line profile based on a 2D original wave.
Definition: pearl-area-import.ipf:1113
prompt_int_quadbg_reduction
variable prompt_int_quadbg_reduction(string *param)
Definition: pearl-scienta-preprocess.ipf:410
scienta_ang_transm
threadsafe variable scienta_ang_transm(wave w, variable x)
Definition: pearl-scienta-preprocess.ipf:578
int_quadbg_reduction
threadsafe variable int_quadbg_reduction(wave source, wave dest1, wave dest2, string *param)
Definition: pearl-scienta-preprocess.ipf:457
gauss2_reduction
threadsafe variable gauss2_reduction(wave source, wave dest1, wave dest2, string *param)
fit horizontal cuts of an image with two gaussian peaks on a linear background
Definition: pearl-scienta-preprocess.ipf:903
test_shockley_anglefit
variable test_shockley_anglefit(wave image, variable branch)
Definition: pearl-scienta-preprocess.ipf:297
int_linbg_reduction
threadsafe variable int_linbg_reduction(wave source, wave dest1, wave dest2, string *param)
Definition: pearl-scienta-preprocess.ipf:197
scienta_norm
variable scienta_norm(wave w, variable x)
Definition: pearl-scienta-preprocess.ipf:544
prompt_redim_linbg_reduction
variable prompt_redim_linbg_reduction(string *param)
parameter dialog for the redim_linbg_reduction() function
Definition: pearl-scienta-preprocess.ipf:716
ad_profile_y_w
threadsafe wave ad_profile_y_w(wave dataset, variable p1, variable p2, wave destwave, variable noavg=defaultValue)
1D cut through 2D dataset along X dimension, existing destination wave.
Definition: pearl-area-profiles.ipf:567
prompt_int_linbg_reduction
variable prompt_int_linbg_reduction(string *param)
Definition: pearl-scienta-preprocess.ipf:36
test_int_quadbg
variable test_int_quadbg(wave image)
Definition: pearl-scienta-preprocess.ipf:440
csr_int_linbg_reduction
string csr_int_linbg_reduction(string win)
Definition: pearl-scienta-preprocess.ipf:76
capture_int_linbg_cursors
string capture_int_linbg_cursors()
Definition: pearl-scienta-preprocess.ipf:66
test_gauss2_reduction
variable test_gauss2_reduction(wave image)
apply the gauss2_reduction function to a single image
Definition: pearl-scienta-preprocess.ipf:810