PEARL Procedures rev-distro-3.1.0-0-gea838b3-dirty
Igor procedures for the analysis of PEARL data
Loading...
Searching...
No Matches
pearl-scienta-preprocess.ipf
Go to the documentation of this file.
1#pragma TextEncoding = "UTF-8"
2#pragma rtGlobals=3 // Use modern global access method and strict wave access.
3#pragma IgorVersion = 6.1
4#pragma ModuleName = PearlScientaPreprocess
5#include "pearl-fitfuncs"
6#include "pearl-area-import"
7
8// Copyright (c) 2013-18 Paul Scherrer Institut
9//
10// Licensed under the Apache License, Version 2.0 (the "License");
11// you may not use this file except in compliance with the License.
12// You may obtain a copy of the License at
13// http://www.apache.org/licenses/LICENSE-2.0
14
29
34
37function prompt_int_linbg_reduction(param)
38 string &param
39
40 variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
41 variable Lsize = NumberByKey("Lsize", param, "=", ";")
42 variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
43 variable Hsize = NumberByKey("Hsize", param, "=", ";")
44 variable Cpos = NumberByKey("Cpos", param, "=", ";")
45 variable Csize = NumberByKey("Csize", param, "=", ";")
46
47 prompt Lcrop, "Lower cropping region"
48 prompt Hcrop, "Upper cropping region"
49 prompt Lsize, "Lower background region"
50 prompt Hsize, "Upper background region"
51 prompt Cpos, "Center position"
52 prompt Csize, "Center integration region"
53
54 doprompt "int_linbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
55 if (v_flag == 0)
56 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
57 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
58 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
59 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
60 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
61 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
62 endif
63
64 return v_flag
65end
66
95function /s capture_int_linbg_cursors()
96 string param = csr_int_linbg_reduction("")
97 svar /z global_params = root:packages:pearl_explorer:s_reduction_params
98 if (svar_exists(global_params))
99 global_params = param
100 endif
101 return param
102end
103
135function /s csr_int_linbg_reduction(win)
136 string win
137
138 // read all cursor positions
139 variable ic
140 string sc
141 variable nc = 10
142 make /n=(nc) /free positions
143 variable np = 0
144 wave /z image = $""
145 string imagename = ""
146 string tracename = ""
147 string info
148
149 for (ic = 0; ic < nc; ic += 1)
150 sc = num2char(char2num("A") + ic)
151 wave /z wc = CsrWaveRef($sc, win)
152 info = CsrInfo($sc, win)
153 tracename = StringByKey("TNAME", info, ":", ";")
154 if (waveexists(wc) && (wavedims(wc) == 2))
155 if (!waveexists(image))
156 wave image = wc
157 imagename = tracename
158 endif
159 if (cmpstr(tracename, imagename) == 0)
160 positions[np] = pcsr($sc, win)
161 np += 1
162 endif
163 endif
164 endfor
165
166 np = floor(np / 2) * 2 // ignore odd cursor
167 redimension /n=(np) positions
168 sort positions, positions
169 // shift upper positions by one so that the rightmost pixel becomes 1.0
170 positions = p >= np / 2 ? positions + 1 : positions
171 positions = positions / dimsize(image, 0)
172
173 // map innermost cursor pair to peak center and size
174 variable ip2 = np / 2
175 variable ip1 = ip2 - 1
176 variable Cpos
177 variable Csize
178 if (ip1 >= 0)
179 Cpos = (positions[ip1] + positions[ip2]) / 2
180 Csize = positions[ip2] - positions[ip1]
181 else
182 // default: a small region in the center
183 Cpos = 0.5
184 Csize = 0.2
185 endif
186
187 // background region
188 ip1 -= 1
189 ip2 += 1
190 variable Lsize
191 variable Hsize
192 if (ip1 >= 0)
193 Lsize = positions[ip1]
194 Hsize = 1 - positions[ip2]
195 else
196 // default: everything outside the peak region
197 Lsize = Cpos - Csize / 2
198 Hsize = 1 - (Cpos + Csize / 2)
199 endif
200
201 // crop region
202 ip1 -= 1
203 ip2 += 1
204 variable Lcrop = 0
205 variable Hcrop = 0
206 if (ip1 >= 0)
207 Lcrop = positions[ip1]
208 Hcrop = 1 - positions[ip2]
209 else
210 // default: in fixed mode: dark corners of the EW4000 at PEARL, 0 otherwise
211 if (dimsize(image, 0) >= 992)
212 Lcrop = 0.11
213 Hcrop = 0.11
214 endif
215 endif
216 Lsize = max(Lsize - Lcrop, 0)
217 Hsize = max(Hsize - Hcrop, 0)
218
219 string param = ""
220 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
221 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
222 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
223 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
224 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
225 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
226
227 return param
228end
229
262threadsafe function /wave int_linbg_reduction(source, param)
263 wave source
264 string &param
265
266 variable nx = dimsize(source, 0)
267 variable ny = dimsize(source, 1)
268
269 // read parameters
270 variable lcrop = NumberByKey("Lcrop", param, "=", ";")
271 variable lsize = NumberByKey("Lsize", param, "=", ";")
272 variable hcrop = NumberByKey("Hcrop", param, "=", ";")
273 variable hsize = NumberByKey("Hsize", param, "=", ";")
274 variable cpos = NumberByKey("Cpos", param, "=", ";")
275 variable csize = NumberByKey("Csize", param, "=", ";")
276
277 make /wave /free /n=2 result_waves
278 make /free /n=0 dest1, dest2
279 result_waves[0] = dest1
280 result_waves[1] = dest2
281 adh5_setup_profile(source, dest1, 1)
282 adh5_setup_profile(source, dest2, 1)
283
284 // validate parameters
285 // background parameters are optional, center parameter is required.
286 if (numtype(lcrop) != 0)
287 lcrop = 0
288 endif
289 if (numtype(lsize) != 0)
290 lsize = 0
291 endif
292 if (numtype(hcrop) != 0)
293 hcrop = 0
294 endif
295 if (numtype(hsize) != 0)
296 hsize = 0
297 endif
298 if (numtype(Cpos) != 0)
299 redimension /n=0 result_waves
300 return result_waves // Cpos parameter missing
301 endif
302 if (numtype(Csize) != 0)
303 redimension /n=0 result_waves
304 return result_waves // Csize parameter missing
305 endif
306
307 variable lpos = lcrop + lsize / 2
308 variable hpos = 1 - (hcrop + hsize / 2)
309
310 variable p0
311 variable p1
312
313 duplicate /free dest1, lbg, hbg
314 if (lsize > 0)
315 p0 = round(lcrop * nx)
316 p1 = round((lcrop + lsize) * nx)
317 ad_profile_y_w(source, p0, p1, lbg)
318 else
319 lbg = 0
320 endif
321 if (hsize > 0)
322 p0 = round((1 - hcrop - hsize) * nx)
323 p1 = round((1 - hcrop) * nx)
324 ad_profile_y_w(source, p0, p1, hbg)
325 else
326 hbg = 0
327 endif
328 if (csize > 0)
329 p0 = round((cpos - csize/2) * nx)
330 p1 = round((cpos + csize/2) * nx)
331 ad_profile_y_w(source, p0, p1, dest1)
332 else
333 dest1 = 0
334 endif
335
336 variable scale = (cpos - lpos) / (hpos - lpos)
337 dest2 = dest1
338 dest1 -= scale * (hbg - lbg) + lbg
339 dest2 = sqrt(dest2 + scale^2 * (hbg + lbg))
340
341 string s_note1
342 string s_note2
343 sprintf s_note1, "AxisLabelD=peak integral"
344 sprintf s_note2, "KineticEnergy=%.3f", cpos * nx * dimdelta(source, 0) + dimoffset(source, 0)
345 Note dest1, s_note1
346 Note dest1, s_note2
347 Note dest2, s_note1
348 Note dest2, s_note2
349
350 return result_waves
351end
352
353function prompt_int_quadbg_reduction(param)
354 string &param
355
356 variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
357 variable Lsize = NumberByKey("Lsize", param, "=", ";")
358 variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
359 variable Hsize = NumberByKey("Hsize", param, "=", ";")
360 variable Cpos = NumberByKey("Cpos", param, "=", ";")
361 variable Csize = NumberByKey("Csize", param, "=", ";")
362
363 prompt Lcrop, "Lower cropping region"
364 prompt Hcrop, "Upper cropping region"
365 prompt Lsize, "Lower background region"
366 prompt Hsize, "Upper background region"
367 prompt Cpos, "Center position"
368 prompt Csize, "Center integration region"
369
370 doprompt "int_quadbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
371 if (v_flag == 0)
372 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
373 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
374 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
375 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
376 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
377 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
378 endif
379
380 return v_flag
381end
382
415threadsafe function /wave int_quadbg_reduction(source, param)
416 wave source
417 string &param
418
419 variable nx = dimsize(source, 0)
420 variable ny = dimsize(source, 1)
421
422 // read parameters
423 variable lcrop = NumberByKey("Lcrop", param, "=", ";")
424 variable lsize = NumberByKey("Lsize", param, "=", ";")
425 variable hcrop = NumberByKey("Hcrop", param, "=", ";")
426 variable hsize = NumberByKey("Hsize", param, "=", ";")
427 variable cpos = NumberByKey("Cpos", param, "=", ";")
428 variable csize = NumberByKey("Csize", param, "=", ";")
429
430 make /wave /free /n=2 result_waves
431 make /free /n=0 dest1, dest2
432 result_waves[0] = dest1
433 result_waves[1] = dest2
434 adh5_setup_profile(source, dest1, 1)
435 adh5_setup_profile(source, dest2, 1)
436
437 // validate parameters
438 // background parameters are optional, center parameter is required.
439 if (numtype(lcrop) != 0)
440 lcrop = 0
441 endif
442 if (numtype(lsize) != 0)
443 lsize = 0
444 endif
445 if (numtype(hcrop) != 0)
446 hcrop = 0
447 endif
448 if (numtype(hsize) != 0)
449 hsize = 0
450 endif
451 if (numtype(Cpos) != 0)
452 redimension /n=0 result_waves
453 return result_waves // Cpos parameter missing
454 endif
455 if (numtype(Csize) != 0)
456 redimension /n=0 result_waves
457 return result_waves // Csize parameter missing
458 endif
459
460 // crop boundaries
461 variable pcl = round(lcrop * nx)
462 variable pch = round((1 - hcrop) * nx)
463 // fit boundaries
464 variable pfl = round((lcrop + lsize) * nx)
465 variable pfh = round((1 - hcrop - hsize) * nx)
466 // integration boundaries
467 variable pil = round((cpos - csize/2) * nx)
468 variable pih = round((cpos + csize/2) * nx)
469
470 // prepare intermediate data buffer
471 make /n=(nx) /d /free profile, mask, fit
472 setscale /p x dimoffset(source,0), dimdelta(source,0), waveunits(source,0), profile, mask, fit
473 mask = ((p >= pcl) && (p < pfl)) || ((p >= pfh) && (p < pch))
474
475 variable qq
476 variable sp, sf
477 variable xil = x2pnt(profile, pil)
478 variable xih = x2pnt(profile, pih)
479
480 make /n=3 /free /d w_coef
481 for (qq = 0; qq < ny; qq += 1)
482 profile = source[p][qq]
483 curvefit /Q /NTHR=1 /W=2 poly 3, kwCWave=w_coef, profile /M=mask
484 fit = poly(w_coef, x)
485 sp = sum(profile, xil, xih)
486 sf = sum(fit, xil, xih)
487 dest1[qq] = sp - sf
488 dest2[qq] = sqrt(sp)
489 endfor
490
491 string s_note1
492 string s_note2
493 sprintf s_note1, "AxisLabelD=peak integral"
494 sprintf s_note2, "KineticEnergy=%.3f", cpos * nx * dimdelta(source, 0) + dimoffset(source, 0)
495 Note dest1, s_note1
496 Note dest1, s_note2
497 Note dest2, s_note1
498 Note dest2, s_note2
499
500 return result_waves
501end
502
511function prompt_redim_linbg_reduction(param)
512 string &param
513
514 variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
515 variable Lsize = NumberByKey("Lsize", param, "=", ";")
516 variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
517 variable Hsize = NumberByKey("Hsize", param, "=", ";")
518 variable Cpos = NumberByKey("Cpos", param, "=", ";")
519 variable Csize = NumberByKey("Csize", param, "=", ";")
520
521 prompt Lcrop, "Lower cropping region"
522 prompt Hcrop, "Upper cropping region"
523 prompt Lsize, "Lower background region"
524 prompt Hsize, "Upper background region"
525 prompt Cpos, "Center position"
526 prompt Csize, "Center integration region"
527
528 doprompt "redim_linbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
529 if (v_flag == 0)
530 param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
531 param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
532 param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
533 param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
534 param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
535 param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
536 endif
537
538 return v_flag
539end
540
574threadsafe function /wave redim_linbg_reduction(source, param)
575 wave source
576 string &param
577
578 variable nx = dimsize(source, 0)
579 variable ny = dimsize(source, 1)
580
581 duplicate /free source, source_redim
582 redimension /n=(nx * ny) source_redim
583 nx += 1
584 redimension /n=(nx, ny) source_redim
585
586 return int_linbg_reduction(source_redim, param)
587end
588
594function test_gauss4_reduction(image)
595 wave image
596
597 string param = ""
598
599 param = ReplaceNumberByKey("rngl", param, -inf, "=", ";")
600 param = ReplaceNumberByKey("rngh", param, inf, "=", ";")
601 param = ReplaceNumberByKey("npeaks", param, 4, "=", ";")
602 param = ReplaceNumberByKey("ybox", param, 1, "=", ";")
603 param = ReplaceNumberByKey("pos1", param, 11, "=", ";")
604 param = ReplaceNumberByKey("wid1", param, 0.1, "=", ";")
605 param = ReplaceNumberByKey("pos2", param, 12, "=", ";")
606 param = ReplaceNumberByKey("wid2", param, 0.2, "=", ";")
607 param = ReplaceNumberByKey("pos3", param, 13, "=", ";")
608 param = ReplaceNumberByKey("wid3", param, 0.3, "=", ";")
609 param = ReplaceNumberByKey("pos4", param, 14, "=", ";")
610 param = ReplaceNumberByKey("wid4", param, 0.4, "=", ";")
611
612 wave /wave results = gauss4_reduction(image, param)
613
614 variable npk = numpnts(results) / 2
615 variable ipk
616 string sw
617 for (ipk = 0; ipk < npk; ipk += 1)
618 sw = "test_int_" + num2str(ipk + 1)
619 duplicate /o results[ipk], $sw
620 sw = "test_sig_" + num2str(ipk + 1)
621 duplicate /o results[ipk + npk], $sw
622 endfor
623end
624
628function prompt_gauss4_reduction(param)
629 string &param
630
631 variable rngl = NumberByKey("rngl", param, "=", ";")
632 variable rngh = NumberByKey("rngh", param, "=", ";")
633 variable pos1 = NumberByKey("pos1", param, "=", ";")
634 variable wid1 = NumberByKey("wid1", param, "=", ";")
635 variable pos2 = NumberByKey("pos2", param, "=", ";")
636 variable wid2 = NumberByKey("wid2", param, "=", ";")
637 variable pos3 = NumberByKey("pos3", param, "=", ";")
638 variable wid3 = NumberByKey("wid3", param, "=", ";")
639 variable pos4 = NumberByKey("pos4", param, "=", ";")
640 variable wid4 = NumberByKey("wid4", param, "=", ";")
641 variable npeaks = NumberByKey("npeaks", param, "=", ";")
642 variable ybox = NumberByKey("ybox", param, "=", ";")
643
644 prompt rngl, "range low"
645 prompt rngh, "range high"
646 prompt pos1, "position 1"
647 prompt wid1, "width 1"
648 prompt pos2, "position 2"
649 prompt wid2, "width 2"
650 prompt pos3, "position 3"
651 prompt wid3, "width 3"
652 prompt pos4, "position 4"
653 prompt wid4, "width 4"
654 prompt npeaks, "number of peaks"
655 prompt ybox, "ybox (1 or 3)"
656
657 doprompt "gauss4_reduction reduction parameters (1/2)", rngl, rngh, npeaks, ybox
658 if (v_flag == 0)
659 param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
660 param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
661 param = ReplaceNumberByKey("npeaks", param, npeaks, "=", ";")
662 param = ReplaceNumberByKey("ybox", param, ybox, "=", ";")
663
664 doprompt "gauss4_reduction reduction parameters (2/2)", pos1, wid1, pos2, wid2, pos3, wid3, pos4, wid4
665 if (v_flag == 0)
666 param = ReplaceNumberByKey("pos1", param, pos1, "=", ";")
667 param = ReplaceNumberByKey("wid1", param, wid1, "=", ";")
668 param = ReplaceNumberByKey("pos2", param, pos2, "=", ";")
669 param = ReplaceNumberByKey("wid2", param, wid2, "=", ";")
670 param = ReplaceNumberByKey("pos3", param, pos3, "=", ";")
671 param = ReplaceNumberByKey("wid3", param, wid3, "=", ";")
672 param = ReplaceNumberByKey("pos4", param, pos4, "=", ";")
673 param = ReplaceNumberByKey("wid4", param, wid4, "=", ";")
674 endif
675 endif
676
677 return v_flag
678end
679
720threadsafe function /wave gauss4_reduction(source, param)
721 wave source
722 string &param
723
724 variable nx = dimsize(source, 0)
725 variable ny = dimsize(source, 1)
726
727 // read parameters
728 variable rngl = NumberByKey("rngl", param, "=", ";")
729 variable rngh = NumberByKey("rngh", param, "=", ";")
730 variable pos1 = NumberByKey("pos1", param, "=", ";")
731 variable wid1 = NumberByKey("wid1", param, "=", ";")
732 variable pos2 = NumberByKey("pos2", param, "=", ";")
733 variable wid2 = NumberByKey("wid2", param, "=", ";")
734 variable pos3 = NumberByKey("pos3", param, "=", ";")
735 variable wid3 = NumberByKey("wid3", param, "=", ";")
736 variable pos4 = NumberByKey("pos4", param, "=", ";")
737 variable wid4 = NumberByKey("wid4", param, "=", ";")
738 variable npeaks = NumberByKey("npeaks", param, "=", ";")
739 variable ybox = NumberByKey("ybox", param, "=", ";")
740
741 // prepare curve fit
742 variable ipk
743 make /free xprof
744 adh5_setup_profile(source, xprof, 0)
745 duplicate /free xprof, xprof_sig
746 variable pl = max(x2pnt(xprof, rngl), 0)
747 variable ph = min(x2pnt(xprof, rngh), numpnts(xprof) - 1)
748
749 make /free /n=(npeaks) peak_coef
750 peak_coef = p * 3 + 2
751 variable n_coef = npeaks * 3 + 2
752 make /free /d /n=14 w_coef, W_sigma
753 w_coef[0] = {0, 0, 1, pos1, wid1, 1, pos2, wid2, 1, pos3, wid3, 1, pos4, wid4}
754 redimension /n=(n_coef) w_coef, w_sigma
755
756 // text constraints cannot be used in threadsafe functions.
757 // the following matrix-vector forumlation is equivalent to:
758 // make /free /T /N=6 constraints
759 // constraints[0] = {"K2 >= 0", "K5 >= 0", "K8 >= 0", "K11 >= 0", "K1 <= 0", "K0 => 0"}
760 make /free /n=(npeaks + 2, numpnts(w_coef)) cmat
761 make /free /n=(npeaks + 2) cvec
762 cmat = 0
763 cmat[0][0] = -1
764 cmat[1][1] = 1
765 cvec = 0
766
767 string hold = "00"
768 for (ipk=0; ipk < npeaks; ipk += 1)
769 hold += "011"
770 cmat[2 + ipk][2 + ipk*3] = -1
771 endfor
772
773 // prepare output
774 make /free /n=(npeaks * 2) /wave result_waves
775 string s_note
776 for (ipk = 0; ipk < npeaks; ipk += 1)
777 make /free /n=0 pk_int
778 adh5_setup_profile(source, pk_int, 1)
779 pk_int = nan
780 sprintf s_note, "AxisLabelD=peak %u integral", ipk+1
781 Note pk_int, s_note
782 sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * 3]
783 Note pk_int, s_note
784 result_waves[ipk] = pk_int
785
786 make /free /n=0 pk_sig
787 adh5_setup_profile(source, pk_sig, 1)
788 pk_sig = nan
789 sprintf s_note, "AxisLabelD=peak %u sigma", ipk+1
790 Note pk_sig, s_note
791 sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * 3]
792 Note pk_sig, s_note
793 result_waves[ipk + npeaks] = pk_sig
794
795 waveclear pk_int, pk_sig
796 endfor
797
798 // loop over angle scale
799 variable p0 = 0
800 variable p1 = dimsize(source, 1) - 1
801 variable pp
802 variable wmin
803 variable wmax
804 if (ybox > 1)
805 p0 += ceil((ybox - 1) / 2)
806 p1 -= ceil((ybox - 1) / 2)
807 endif
808 variable V_FitNumIters
809 variable V_FitError
810
811 for (pp = p0; pp <= p1; pp += 1)
812 // box average
813 xprof = source[p][pp]
814 if (ybox > 1)
815 xprof += source[p][pp-1] + source[p][pp+1]
816 endif
817 xprof_sig = max(sqrt(xprof), 1)
818 xprof /= ybox
819 xprof_sig /= ybox
820
821 // generate guess
822 wmin = wavemin(xprof)
823 wmax = wavemax(xprof)
824 w_coef[0] = wmin
825 w_coef[1] = 0
826 for (ipk=0; ipk < npeaks; ipk += 1)
827 w_coef[2 + ipk*3] = wmax - wmin
828 endfor
829
830 V_FitError = 0
831 FuncFit /H=hold /Q /NTHR=1 /N /W=2 MultiGaussLinBG_AO w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
832 wave w_sigma
833
834 // retrieve results, leave them at nan if the fit did not converge
835 if (V_FitNumIters < 40)
836 for (ipk = 0; ipk < npeaks; ipk += 1)
837 wave val = result_waves[ipk]
838 wave sig = result_waves[ipk + npeaks]
839 val[pp] = max(w_coef[peak_coef[ipk]], 0)
840 sig[pp] = max(w_sigma[peak_coef[ipk]], 0)
841 endfor
842 endif
843 endfor
844
845 // calculate integral
846 for (ipk = 0; ipk < npeaks; ipk += 1)
847 wave val = result_waves[ipk]
848 wave sig = result_waves[ipk + npeaks]
849 val *= w_coef[peak_coef[ipk] + 2] * sqrt(pi)
850 sig *= w_coef[peak_coef[ipk] + 2] * sqrt(pi)
851 endfor
852
853 return result_waves
854end
855
856threadsafe function /wave gauss6_reduction(source, param)
857 wave source
858 string &param
859
860 variable nx = dimsize(source, 0)
861 variable ny = dimsize(source, 1)
862
863 // read parameters
864 variable rngl = NumberByKey("rngl", param, "=", ";")
865 variable rngh = NumberByKey("rngh", param, "=", ";")
866 variable pos1 = NumberByKey("pos1", param, "=", ";")
867 variable wid1 = NumberByKey("wid1", param, "=", ";")
868 variable pos2 = NumberByKey("pos2", param, "=", ";")
869 variable wid2 = NumberByKey("wid2", param, "=", ";")
870 variable pos3 = NumberByKey("pos3", param, "=", ";")
871 variable wid3 = NumberByKey("wid3", param, "=", ";")
872 variable pos4 = NumberByKey("pos4", param, "=", ";")
873 variable wid4 = NumberByKey("wid4", param, "=", ";")
874 variable pos5 = NumberByKey("pos5", param, "=", ";")
875 variable wid5 = NumberByKey("wid5", param, "=", ";")
876 variable pos6 = NumberByKey("pos6", param, "=", ";")
877 variable wid6 = NumberByKey("wid6", param, "=", ";")
878 variable npeaks = NumberByKey("npeaks", param, "=", ";")
879 variable ybox = NumberByKey("ybox", param, "=", ";")
880
881 // prepare curve fit
882 variable ipk
883 make /free xprof
884 adh5_setup_profile(source, xprof, 0)
885 duplicate /free xprof, xprof_sig
886 variable pl = max(x2pnt(xprof, rngl), 0)
887 variable ph = min(x2pnt(xprof, rngh), numpnts(xprof) - 1)
888
889 make /free /n=(npeaks) peak_coef
890 peak_coef = p * 3 + 2
891 variable n_coef = npeaks * 3 + 2
892 make /free /d /n=(n_coef) w_coef, W_sigma
893 w_coef[0] = {0, 0, 1, pos1, wid1, 1, pos2, wid2, 1, pos3, wid3, 1, pos4, wid4, 1, pos5, wid5, 1, pos6, wid6}
894 redimension /n=(n_coef) w_coef, w_sigma
895
896 // text constraints cannot be used in threadsafe functions.
897 // the following matrix-vector forumlation is equivalent to:
898 // make /free /T /N=6 constraints
899 // constraints[0] = {"K2 >= 0", "K5 >= 0", "K8 >= 0", "K11 >= 0", "K1 <= 0", "K0 => 0"}
900 make /free /n=(npeaks + 2, numpnts(w_coef)) cmat
901 make /free /n=(npeaks + 2) cvec
902 cmat = 0
903 cmat[0][0] = -1
904 cmat[1][1] = 1
905 cvec = 0
906
907 string hold = "00"
908 for (ipk=0; ipk < npeaks; ipk += 1)
909 hold += "011"
910 cmat[2 + ipk][2 + ipk*3] = -1
911 endfor
912
913 // prepare output
914 make /free /n=(npeaks * 2) /wave result_waves
915 string s_note
916 for (ipk = 0; ipk < npeaks; ipk += 1)
917 make /free /n=0 pk_int
918 adh5_setup_profile(source, pk_int, 1)
919 pk_int = nan
920 sprintf s_note, "AxisLabelD=peak %u integral", ipk+1
921 Note pk_int, s_note
922 sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * 3]
923 Note pk_int, s_note
924 result_waves[ipk] = pk_int
925
926 make /free /n=0 pk_sig
927 adh5_setup_profile(source, pk_sig, 1)
928 pk_sig = nan
929 sprintf s_note, "AxisLabelD=peak %u sigma", ipk+1
930 Note pk_sig, s_note
931 sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * 3]
932 Note pk_sig, s_note
933 result_waves[ipk + npeaks] = pk_sig
934
935 waveclear pk_int, pk_sig
936 endfor
937
938 // loop over angle scale
939 variable p0 = 0
940 variable p1 = dimsize(source, 1) - 1
941 variable pp
942 variable wmin
943 variable wmax
944 if (ybox > 1)
945 p0 += ceil((ybox - 1) / 2)
946 p1 -= ceil((ybox - 1) / 2)
947 endif
948 variable V_FitNumIters
949 variable V_FitError
950
951 for (pp = p0; pp <= p1; pp += 1)
952 // box average
953 xprof = source[p][pp]
954 if (ybox > 1)
955 xprof += source[p][pp-1] + source[p][pp+1]
956 endif
957 xprof_sig = max(sqrt(xprof), 1)
958 xprof /= ybox
959 xprof_sig /= ybox
960
961 // generate guess
962 wmin = wavemin(xprof)
963 wmax = wavemax(xprof)
964 w_coef[0] = wmin
965 w_coef[1] = 0
966 for (ipk=0; ipk < npeaks; ipk += 1)
967 w_coef[2 + ipk*3] = wmax - wmin
968 endfor
969
970 V_FitError = 0
971 FuncFit /H=hold /Q /NTHR=1 /N /W=2 MultiGaussLinBG_AO w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
972 wave w_sigma
973
974 // retrieve results, leave them at nan if the fit did not converge
975 if (V_FitNumIters < 40)
976 for (ipk = 0; ipk < npeaks; ipk += 1)
977 wave val = result_waves[ipk]
978 wave sig = result_waves[ipk + npeaks]
979 val[pp] = max(w_coef[peak_coef[ipk]], 0)
980 sig[pp] = max(w_sigma[peak_coef[ipk]], 0)
981 endfor
982 endif
983 endfor
984
985 // calculate integral
986 for (ipk = 0; ipk < npeaks; ipk += 1)
987 wave val = result_waves[ipk]
988 wave sig = result_waves[ipk + npeaks]
989 val *= w_coef[peak_coef[ipk] + 2] * sqrt(pi)
990 sig *= w_coef[peak_coef[ipk] + 2] * sqrt(pi)
991 endfor
992
993 return result_waves
994end
995
996
997
998
1030function /s find_gauss4_reduction_params(spectrum, peakpos)
1031 wave spectrum
1032 wave peakpos
1033 string param = ""
1034
1035 variable wmin = wavemin(spectrum)
1036 variable wmax = wavemax(spectrum)
1037
1038 // read parameters
1039 variable rngl = dimoffset(spectrum, 0)
1040 variable rngh = dimoffset(spectrum, 0) + dimdelta(spectrum, 0) * (dimsize(spectrum, 0) - 1)
1041 make /n=4 /free positions, widths
1042 variable npeaks = numpnts(peakpos)
1043 variable ybox = 1
1044 positions = 0
1045 positions[0, npeaks-1] = peakpos[p]
1046 widths = 0.2
1047
1048 variable n_coef = npeaks * 3 + 2
1049 make /free /d /n=(n_coef) w_coef
1050 w_coef = 0
1051 w_coef[0] = wmin
1052 w_coef[1] = 0
1053
1054 make /free /n=(2+npeaks, numpnts(w_coef)) cmat
1055 make /free /n=(2+npeaks) cvec
1056 cmat = 0
1057 cmat[0][0] = -1
1058 cmat[1][1] = 1
1059 cvec = 0
1060
1061 variable ip
1062 for (ip=0; ip < npeaks; ip += 1)
1063 cmat[2 + ip][2 + ip*3] = -1
1064 w_coef[2 + ip*3] = wmax - wmin
1065 w_coef[3 + ip*3] = peakpos[ip]
1066 w_coef[4 + ip*3] = widths[ip]
1067 endfor
1068
1069 variable V_FitNumIters
1070 FuncFit /Q /NTHR=1 /N MultiGaussLinBG w_coef spectrum /C={cmat, cvec}
1071
1072 for (ip=0; ip < npeaks; ip += 1)
1073 positions[ip] = w_coef[3 + ip * 3]
1074 widths[ip ] = abs(w_coef[4 + ip * 3])
1075 endfor
1076 for (ip=npeaks; ip < 4; ip += 1)
1077 positions[ip] = 0
1078 widths[ip] = 0.2
1079 endfor
1080
1081 param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
1082 param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
1083 param = ReplaceNumberByKey("npeaks", param, npeaks, "=", ";")
1084 param = ReplaceNumberByKey("ybox", param, ybox, "=", ";")
1085 param = ReplaceNumberByKey("pos1", param, positions[0], "=", ";")
1086 param = ReplaceNumberByKey("wid1", param, widths[0], "=", ";")
1087 param = ReplaceNumberByKey("pos2", param, positions[1], "=", ";")
1088 param = ReplaceNumberByKey("wid2", param, widths[1], "=", ";")
1089 param = ReplaceNumberByKey("pos3", param, positions[2], "=", ";")
1090 param = ReplaceNumberByKey("wid3", param, widths[2], "=", ";")
1091 param = ReplaceNumberByKey("pos4", param, positions[3], "=", ";")
1092 param = ReplaceNumberByKey("wid4", param, widths[3], "=", ";")
1093
1094 return param
1095end
1096
1138threadsafe function /wave voigt4_reduction(source, param)
1139 wave source
1140 string &param
1141
1142 dfref orig_dfr = GetDataFolderDFR()
1143
1144 variable nx = dimsize(source, 0)
1145 variable ny = dimsize(source, 1)
1146
1147 // read parameters
1148 variable rngl = NumberByKey("rngl", param, "=", ";")
1149 variable rngh = NumberByKey("rngh", param, "=", ";")
1150 variable pos1 = NumberByKey("pos1", param, "=", ";")
1151 variable wid1 = NumberByKey("wid1", param, "=", ";")
1152 variable shp1 = NumberByKey("shp1", param, "=", ";")
1153 variable pos2 = NumberByKey("pos2", param, "=", ";")
1154 variable wid2 = NumberByKey("wid2", param, "=", ";")
1155 variable shp2 = NumberByKey("shp2", param, "=", ";")
1156 variable pos3 = NumberByKey("pos3", param, "=", ";")
1157 variable wid3 = NumberByKey("wid3", param, "=", ";")
1158 variable shp3 = NumberByKey("shp3", param, "=", ";")
1159 variable pos4 = NumberByKey("pos4", param, "=", ";")
1160 variable wid4 = NumberByKey("wid4", param, "=", ";")
1161 variable shp4 = NumberByKey("shp4", param, "=", ";")
1162 variable npeaks = NumberByKey("npeaks", param, "=", ";")
1163
1164 // prepare curve fit
1165 variable ipk
1166 make /free xprof
1167 adh5_setup_profile(source, xprof, 0)
1168 duplicate /free xprof, xprof_sig
1169 variable pl = max(x2pnt(xprof, rngl), 0)
1170 variable ph = min(x2pnt(xprof, rngh), numpnts(xprof) - 1)
1171
1172 make /free /n=(npeaks) peak_coef
1173 variable coef_per_peak = 4
1174 peak_coef = p * coef_per_peak + 2
1175 variable n_coef = npeaks * coef_per_peak + 2
1176 make /free /d /n=18 w_coef, W_sigma
1177 w_coef[0] = {0, 0, 1, pos1, wid1, shp1, 1, pos2, wid2, shp2, 1, pos3, wid3, shp3, 1, pos4, wid4, shp4}
1178 redimension /n=(n_coef) w_coef, w_sigma
1179
1180 // text constraints cannot be used in threadsafe functions.
1181 // the following matrix-vector formulation enforces all peak amplitudes to be positive.
1182 make /free /n=(npeaks, numpnts(w_coef)) cmat
1183 make /free /n=(npeaks) cvec
1184 cmat = 0
1185 cvec = 0
1186
1187 string hold = "00"
1188 for (ipk=0; ipk < npeaks; ipk += 1)
1189 hold += "0111"
1190 cmat[ipk][2 + ipk * coef_per_peak] = -1
1191 endfor
1192
1193 // prepare output
1194 make /free /n=(npeaks * 2) /wave result_waves
1195 string s_note
1196 for (ipk = 0; ipk < npeaks; ipk += 1)
1197 make /free /n=0 pk_int
1198 adh5_setup_profile(source, pk_int, 1)
1199 pk_int = nan
1200 sprintf s_note, "AxisLabelD=peak %u integral", ipk+1
1201 Note pk_int, s_note
1202 sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * coef_per_peak]
1203 Note pk_int, s_note
1204 result_waves[ipk] = pk_int
1205
1206 make /free /n=0 pk_sig
1207 adh5_setup_profile(source, pk_sig, 1)
1208 pk_sig = nan
1209 sprintf s_note, "AxisLabelD=peak %u sigma", ipk+1
1210 Note pk_sig, s_note
1211 sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * coef_per_peak]
1212 Note pk_sig, s_note
1213 result_waves[ipk + npeaks] = pk_sig
1214
1215 waveclear pk_int, pk_sig
1216 endfor
1217
1218 // loop over angle scale
1219 variable p0 = 0
1220 variable p1 = dimsize(source, 1) - 1
1221 variable pp
1222 variable wmin
1223 variable wmax
1224 variable V_FitNumIters
1225 variable V_FitError
1226
1227 for (pp = p0; pp <= p1; pp += 1)
1228 xprof = source[p][pp]
1229 xprof_sig = max(sqrt(xprof), 1)
1230
1231 // generate guess
1232 wmin = wavemin(xprof)
1233 wmax = wavemax(xprof)
1234 w_coef[0] = wmin
1235 w_coef[1] = 0
1236
1237 for (ipk=0; ipk < npeaks; ipk += 1)
1238 w_coef[2 + ipk * coef_per_peak] = wmax - wmin
1239 endfor
1240
1241 V_FitError = 0
1242 FuncFit /H=hold /Q /N /W=2 MultiVoigtLinBG_AO w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
1243 wave w_sigma
1244
1245 // retrieve results, leave them at nan if the fit did not converge
1246 if (V_FitNumIters < 40)
1247 for (ipk = 0; ipk < npeaks; ipk += 1)
1248 wave val = result_waves[ipk]
1249 wave sig = result_waves[ipk + npeaks]
1250 val[pp] = max(w_coef[peak_coef[ipk]], 0)
1251 sig[pp] = max(w_sigma[peak_coef[ipk]], 0)
1252 endfor
1253 endif
1254 endfor
1255
1256 SetDataFolder orig_dfr
1257 return result_waves
1258end
1259
1260
1264function prompt_voigt4_reduction(param)
1265 string &param
1266
1267 variable rngl = NumberByKey("rngl", param, "=", ";")
1268 variable rngh = NumberByKey("rngh", param, "=", ";")
1269 variable pos1 = NumberByKey("pos1", param, "=", ";")
1270 variable wid1 = NumberByKey("wid1", param, "=", ";")
1271 variable shp1 = NumberByKey("shp1", param, "=", ";")
1272 variable pos2 = NumberByKey("pos2", param, "=", ";")
1273 variable wid2 = NumberByKey("wid2", param, "=", ";")
1274 variable shp2 = NumberByKey("shp2", param, "=", ";")
1275 variable pos3 = NumberByKey("pos3", param, "=", ";")
1276 variable wid3 = NumberByKey("wid3", param, "=", ";")
1277 variable shp3 = NumberByKey("shp3", param, "=", ";")
1278 variable pos4 = NumberByKey("pos4", param, "=", ";")
1279 variable wid4 = NumberByKey("wid4", param, "=", ";")
1280 variable shp4 = NumberByKey("shp4", param, "=", ";")
1281 variable npeaks = NumberByKey("npeaks", param, "=", ";")
1282 variable dummy = nan
1283
1284 prompt rngl, "range low"
1285 prompt rngh, "range high"
1286 prompt pos1, "position 1"
1287 prompt wid1, "width 1"
1288 prompt shp1, "shape 1"
1289 prompt pos2, "position 2"
1290 prompt wid2, "width 2"
1291 prompt shp2, "shape 2"
1292 prompt pos3, "position 3"
1293 prompt wid3, "width 3"
1294 prompt shp3, "shape 3"
1295 prompt pos4, "position 4"
1296 prompt wid4, "width 4"
1297 prompt shp4, "shape 4"
1298 prompt npeaks, "number of peaks"
1299 prompt dummy, "(not used)"
1300
1301 doprompt "voigt4_reduction reduction parameters (1/2)", rngl, rngh, npeaks, dummy, pos1, pos2, pos3, pos4
1302 if (v_flag == 0)
1303 param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
1304 param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
1305 param = ReplaceNumberByKey("npeaks", param, npeaks, "=", ";")
1306 param = ReplaceNumberByKey("pos1", param, pos1, "=", ";")
1307 param = ReplaceNumberByKey("pos2", param, pos2, "=", ";")
1308 param = ReplaceNumberByKey("pos3", param, pos3, "=", ";")
1309 param = ReplaceNumberByKey("pos4", param, pos4, "=", ";")
1310
1311 doprompt "voigt4_reduction reduction parameters (2/2)", wid1, shp1, wid2, shp2, wid3, shp3, wid4, shp4
1312 if (v_flag == 0)
1313 param = ReplaceNumberByKey("wid1", param, wid1, "=", ";")
1314 param = ReplaceNumberByKey("shp1", param, shp1, "=", ";")
1315 param = ReplaceNumberByKey("wid2", param, wid2, "=", ";")
1316 param = ReplaceNumberByKey("shp2", param, shp2, "=", ";")
1317 param = ReplaceNumberByKey("wid3", param, wid3, "=", ";")
1318 param = ReplaceNumberByKey("shp3", param, shp3, "=", ";")
1319 param = ReplaceNumberByKey("wid4", param, wid4, "=", ";")
1320 param = ReplaceNumberByKey("shp4", param, shp4, "=", ";")
1321 endif
1322 endif
1323
1324 return v_flag
1325end
1326
1327
1335function test_shockley_anglefit(image, branch)
1336 wave image
1337 variable branch
1338
1339 string param = ""
1340 param = ReplaceStringByKey("branch", param, num2str(branch), "=", ";")
1341
1342 string s_branch
1343 if (branch >= 0)
1344 s_branch = "p"
1345 else
1346 s_branch = "n"
1347 endif
1348 string pkpos_name = "saf_pkpos_" + s_branch
1349 string pkwid_name = "saf_pkwid_" + s_branch
1350
1351 wave /wave results = shockley_anglefit(image, param)
1352 duplicate results[0], $pkpos_name
1353 duplicate results[1], $pkwid_name
1354end
1355
1356function prompt_Shockley_anglefit(param)
1357 string &param
1358
1359 variable branch = NumberByKey("branch", param, "=", ";")
1360
1361 prompt branch, "Branch (-1 or +1)"
1362
1363 doprompt "Shockley_anglefit_reduction Parameters", branch
1364 if (v_flag == 0)
1365 param = ReplaceNumberByKey("branch", param, branch, "=", ";")
1366 endif
1367
1368 return v_flag
1369end
1370
1393threadsafe function /wave Shockley_anglefit(source, param)
1394 wave source
1395 string &param
1396
1397 variable nx = dimsize(source, 0)
1398 variable ny = dimsize(source, 1)
1399
1400 // read parameters
1401 variable branch = NumberByKey("branch", param, "=", ";")
1402
1403 // validate parameters
1404 if (numtype(branch) != 0)
1405 branch = +1
1406 endif
1407
1408 // prepare output
1409 make /wave /free /n=2 result_waves
1410 make /free /n=0 dest1, dest2
1411 result_waves[0] = dest1
1412 result_waves[1] = dest2
1413 adh5_setup_profile(source, dest1, 0)
1414 adh5_setup_profile(source, dest2, 0)
1415 dest1 = nan
1416 dest2 = nan
1417
1418 // select angle range
1419 // hard-coded for a particular measurement series
1420 variable y0
1421 variable y1
1422 if (branch < 0)
1423 y0 = -5
1424 y1 = 0
1425 else
1426 y0 = 0
1427 y1 = 5
1428 endif
1429
1430 // select energy range
1431 // start at the point of highest intensity and go up 0.45 eV
1432 variable p0
1433 variable p1
1434 variable q0
1435 variable q1
1436 duplicate /free dest1, center
1437 q0 = round((y0 - dimoffset(source, 1)) / dimdelta(source, 1))
1438 q1 = round((y1 - dimoffset(source, 1)) / dimdelta(source, 1))
1439 ad_profile_x_w(source, q0, q1, center)
1440 wavestats /q/m=1 center
1441 p0 = round((v_maxloc - dimoffset(source, 0)) / dimdelta(source, 0))
1442 p1 = round((v_maxloc + 0.4 - dimoffset(source, 0)) / dimdelta(source, 0))
1443
1444 // prepare intermediate data buffer
1445 make /n=(ny)/d/free profile
1446 setscale /p x dimoffset(source,1), dimdelta(source,1), waveunits(source,1), profile
1447
1448 variable pp
1449 for (pp = p0; pp <= p1; pp += 1)
1450 profile = source[pp][p]
1451 curvefit /Q /NTHR=1 /W=2 gauss profile(y0,y1)
1452 wave w_coef
1453 dest1[pp] = w_coef[2]
1454 dest2[pp] = w_coef[3]
1455 endfor
1456
1457 return result_waves
1458end
1459
1460function scienta_norm(w, x): fitfunc
1461 wave w
1462 variable x
1463
1464 return w[0] * (x^2 - w[1]^2)
1465end
1466
1467function /wave fit_scienta_ang_transm(data, params)
1468 wave data // measured angular distribution (1D)
1469 wave /z params
1470
1471 if (!waveexists(params))
1472 make /n=12 /o params
1473 endif
1474 redimension /n=12/d params
1475
1476 variable h = wavemax(data) - wavemin(data)
1477 params[0] = h / 2
1478 params[1] = 0
1479 params[2] = 7
1480 params[3] = h / 4
1481 params[4] = -23
1482 params[5] = 4
1483 params[6] = h / 4
1484 params[7] = +23
1485 params[8] = 4
1486 params[9] = h / 2
1487 params[10] = 0
1488 params[11] = -0.001
1489 FuncFit /NTHR=0 /q scienta_ang_transm params data
1490
1491 return params
1492end
1493
1494threadsafe function scienta_ang_transm(w, x): fitfunc
1495 // parameterized angular transmission function of the analyser
1496 wave w // coefficients
1497 // w[0] = amplitude gauss 1
1498 // w[1] = position gauss 1
1499 // w[2] = width gauss 1
1500 // w[3] = amplitude gauss 2
1501 // w[4] = position gauss 2
1502 // w[5] = width gauss 2
1503 // w[6] = amplitude gauss 3
1504 // w[7] = position gauss 3
1505 // w[8] = width gauss 3
1506 // w[9] = constant background
1507 // w[10] = linear background
1508 // w[11] = quadratic background
1509 variable x
1510
1511 make /free /n=4 /d w_int
1512 w_int[0] = 0
1513 w_int[1,] = w[p - 1]
1514 variable pk1 = gauss1d(w_int, x)
1515 w_int[1,] = w[p + 2]
1516 variable pk2 = gauss1d(w_int, x)
1517 w_int[1,] = w[p + 5]
1518 variable pk3 = gauss1d(w_int, x)
1519 w_int[0,2] = w[9 + p]
1520 w_int[3] = 0
1521 variable bg = poly(w_int, x)
1522
1523 return bg + pk1 + pk2 + pk3
1524end
1525
1526function /wave fit_scienta_poly_bg(data, params, bgterms)
1527 wave data // measured distribution (2D)
1528 wave /z params // wave, will be redimensioned for the correct size
1529 variable bgterms // number of terms in the polynomial background: 2, 3, or 4
1530
1531 if (!waveexists(params))
1532 make /n=15 /o params
1533 endif
1534 redimension /n=15 /d params
1535
1536 variable wmax = wavemax(data)
1537 variable wmin = wavemin(data)
1538 params[0] = 0
1539 params[1] = 7
1540 params[2] = 1 / 2
1541 params[3] = -23
1542 params[4] = 4
1543 params[5] = 1 / 2
1544 params[6] = +23
1545 params[7] = 4
1546 params[8] = 1
1547 params[9] = 0
1548 params[10] = -0.001
1549 params[11] = wmin
1550 params[12] = (wmax - wmin) / dimdelta(data,1) / dimsize(data,1)
1551 params[13] = 0
1552 params[14] = 0
1553
1554 string h = "0000000000000"
1555 if (bgterms < 3)
1556 h = h + "1"
1557 else
1558 h = h + "0"
1559 endif
1560 if (bgterms < 4)
1561 h = h + "1"
1562 else
1563 h = h + "0"
1564 endif
1565 FuncFitMD /NTHR=1 /q /h=h scienta_poly_bg params data
1566
1567 return params
1568end
1569
1570function scienta_poly_bg(w, e, a): fitfunc
1571 // polynomial background with
1572 // parameterized angular transmission function of the analyser
1573 wave w // coefficients
1574 // angular transmission, varies with a
1575 // amplitude of gauss 1 = 1 constant
1576 // other peak amplitudes and linear terms are relative to gauss 1
1577 // w[0] = position gauss 1
1578 // w[1] = width gauss 1
1579 // w[2] = amplitude gauss 2, relative to gauss 1
1580 // w[3] = position gauss 2
1581 // w[4] = width gauss 2
1582 // w[5] = amplitude gauss 3, relative to gauss 1
1583 // w[6] = position gauss 3
1584 // w[7] = width gauss 3
1585 // w[8] = constant term
1586 // w[9] = linear term
1587 // w[10] = quadratic term
1588 // spectral background, varies with e
1589 // w[11] = constant term
1590 // w[12] = linear term
1591 // w[13] = quadratic term
1592 // w[14] = cubic term
1593 variable a // detection angle
1594 variable e // electron energy
1595
1596 make /free /n=4 /d w_int
1597 variable p0 = 0
1598
1599 w_int[0] = 0
1600 w_int[1] = 1
1601 w_int[2,] = w[p0 + p - 2]
1602 variable pk1 = gauss1d(w_int, a)
1603 p0 += 2
1604
1605 w_int[1,] = w[p0 + p - 1]
1606 variable pk2 = gauss1d(w_int, a)
1607 p0 += 3
1608
1609 w_int[1,] = w[p0 + p - 1]
1610 variable pk3 = gauss1d(w_int, a)
1611 p0 += 3
1612
1613 w_int[0,2] = w[p0 + p]
1614 w_int[3] = 0
1615 variable base = poly(w_int, a)
1616 p0 += 3
1617
1618 w_int[0,3] = w[p0 + p]
1619 variable bg = poly(w_int, e)
1620
1621 return bg * (base + pk1 + pk2 + pk3)
1622end
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:1138
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:568
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:505
MultiGaussLinBG_AO
threadsafe variable MultiGaussLinBG_AO(wave pw, wave yw, wave xw)
multiple gaussian peaks on a linear background fit function (all at once).
Definition pearl-fitfuncs.ipf:80
MultiVoigtLinBG_AO
threadsafe variable MultiVoigtLinBG_AO(wave pw, wave yw, wave xw)
multiple voigt peaks on a linear background fit function.
Definition pearl-fitfuncs.ipf:224
MultiGaussLinBG
threadsafe variable MultiGaussLinBG(wave w, variable x)
multiple gaussian peaks on a linear background fit function.
Definition pearl-fitfuncs.ipf:45
gauss6_reduction
threadsafe wave gauss6_reduction(wave source, string *param)
Definition pearl-scienta-preprocess.ipf:856
prompt_int_linbg_reduction
variable prompt_int_linbg_reduction(string *param)
prompt the user for integrate on linear background reduction parameters.
Definition pearl-scienta-preprocess.ipf:37
prompt_gauss4_reduction
variable prompt_gauss4_reduction(string *param)
prompt for the gauss4_reduction parameters
Definition pearl-scienta-preprocess.ipf:628
int_linbg_reduction
threadsafe wave int_linbg_reduction(wave source, string *param)
linear-background subtracted integration reduction function.
Definition pearl-scienta-preprocess.ipf:262
prompt_int_quadbg_reduction
variable prompt_int_quadbg_reduction(string *param)
Definition pearl-scienta-preprocess.ipf:353
prompt_redim_linbg_reduction
variable prompt_redim_linbg_reduction(string *param)
parameter dialog for the redim_linbg_reduction() function
Definition pearl-scienta-preprocess.ipf:511
gauss4_reduction
threadsafe wave gauss4_reduction(wave source, string *param)
fit horizontal cuts of an image with up to four gaussian peaks on a linear background
Definition pearl-scienta-preprocess.ipf:720
redim_linbg_reduction
threadsafe wave redim_linbg_reduction(wave source, string *param)
linear background reduction function for incorrectly dimensioned scienta image
Definition pearl-scienta-preprocess.ipf:574
csr_int_linbg_reduction
string csr_int_linbg_reduction(string win)
calculate linear background reduction parameters from cursors in a graph.
Definition pearl-scienta-preprocess.ipf:135
int_quadbg_reduction
threadsafe wave int_quadbg_reduction(wave source, string *param)
integrate peak area minus a quadratic background
Definition pearl-scienta-preprocess.ipf:415
test_gauss4_reduction
variable test_gauss4_reduction(wave image)
apply the gauss4_reduction function to a single image
Definition pearl-scienta-preprocess.ipf:594
capture_int_linbg_cursors
string capture_int_linbg_cursors()
capture linear background reduction parameters from cursors in a graph.
Definition pearl-scienta-preprocess.ipf:95