1 | #!/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | """ |
---|
4 | Fit model using multiple scattering. |
---|
5 | |
---|
6 | As of this writing, multiscattering isn't integrated into sasmodels, and a |
---|
7 | number of hacks are required to get it running as a fit. |
---|
8 | |
---|
9 | The appropriate items need to be on the python path. These include |
---|
10 | sasview (for reading the data), sasmodels and bumps. The multiscat module |
---|
11 | (currently in the sasmodels/explore directory) is also needed, either beside |
---|
12 | this example fit file, or by putting sasmdoels/explore on the python path. |
---|
13 | |
---|
14 | On Unix/Mac running as developer I do:: |
---|
15 | |
---|
16 | # Show the model without fitting |
---|
17 | PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py |
---|
18 | |
---|
19 | # Run the fit |
---|
20 | PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \ |
---|
21 | multiscatfit.py --store=/tmp/t1 |
---|
22 | |
---|
23 | You may be able to run multiscatfit.py against the distributed sasview |
---|
24 | application (if it is new enough, and if you have multiscat.py in the |
---|
25 | same directory). You probably need a command such as:: |
---|
26 | |
---|
27 | sasview.exe bumps.cli multiscatfit.py --store=t1 |
---|
28 | """ |
---|
29 | |
---|
30 | import sys |
---|
31 | from bumps.names import * |
---|
32 | from sasmodels.core import load_model |
---|
33 | from sasmodels.bumps_model import Model, Experiment |
---|
34 | from sasmodels.data import load_data, set_beam_stop, set_top |
---|
35 | |
---|
36 | from multiscat import MultipleScattering |
---|
37 | |
---|
38 | ## Load the data |
---|
39 | #data = load_data('DEC07267.DAT') |
---|
40 | #set_beam_stop(data, 0.003, outer=0.025) |
---|
41 | data = load_data('latex_smeared.xml', index=1) |
---|
42 | |
---|
43 | ## Define the model |
---|
44 | kernel = load_model("ellipsoid") |
---|
45 | |
---|
46 | model = Model( |
---|
47 | kernel, |
---|
48 | scale=0.005, background=0.05, |
---|
49 | radius_polar=2200, radius_equatorial=2200, |
---|
50 | sld=.291, sld_solvent=7.105, |
---|
51 | #theta=90, theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, |
---|
52 | #phi=90, phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, |
---|
53 | radius_polar_pd=0.222296, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, |
---|
54 | radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0, |
---|
55 | ) |
---|
56 | |
---|
57 | # SET THE FITTING PARAMETERS |
---|
58 | model.radius_polar.range(15, 3000) |
---|
59 | model.radius_equatorial.range(15, 3000) |
---|
60 | #model.theta.range(0, 90) |
---|
61 | #model.theta_pd.range(0,10) |
---|
62 | #model.phi_pd.range(0,20) |
---|
63 | #model.phi.range(0, 180) |
---|
64 | model.background.range(0,1000) |
---|
65 | model.scale.range(0, 0.1) |
---|
66 | |
---|
67 | # Mulitple scattering probability parameter |
---|
68 | # HACK: the parameter is assigned to model.theta, which is otherwise unused |
---|
69 | # since the dataset is 1D; this won't work for 2D data |
---|
70 | probability = Parameter(name="probability", value=0.0) |
---|
71 | probability.range(0.0, 0.9) |
---|
72 | model.phi = probability |
---|
73 | |
---|
74 | M = Experiment(data=data, model=model) |
---|
75 | |
---|
76 | # Stack mulitple scattering on top of the existing resolution function. |
---|
77 | # Because resolution functions in sasview don't have fitting parameters, |
---|
78 | # we instead allow the multiple scattering calculator to take a function |
---|
79 | # instead of a probability. This function returns the current value of |
---|
80 | # the parameter. ** THIS IS TEMPORARY ** when multiple scattering is |
---|
81 | # properly integrated into sasmodels and sasview, its fittable parameter |
---|
82 | # will be treated like the model parameters. |
---|
83 | M.resolution = MultipleScattering(resolution=M.resolution, |
---|
84 | probability=lambda: probability.value, |
---|
85 | ) |
---|
86 | M._kernel_inputs = M.resolution.q_calc |
---|
87 | problem = FitProblem(M) |
---|
88 | |
---|
89 | if __name__ == "__main__": |
---|
90 | #M.theory() |
---|
91 | M.plot() |
---|
92 | import pylab; pylab.show() |
---|