Ticket #6863: t6863_refine.py

File t6863_refine.py, 9.8 KB (added by Wenduo Zhou, 7 years ago)
Line 
1######################################################################
2#
3# This is a partial copy from LeBailFitScript.py
4#
5# Python Script as Step 1 of Le Bail Fitting to
6# 1. Load file
7# 2. Create LeBailFitInput
8# 3. Fit Peaks
9#
10# Step 1:   Load data, model (HKL list) and starting instrument parameters values,
11#           and do an initial LeBailFit/calculation to see how much the starting
12#           values are off;
13#
14######################################################################
15from Calibration_ImportInformation import *
16
17global numsteps
18
19numsteps = 2002  # Step 4
20numprofsteps = 20000 # Step 4.1 (step 5)
21
22#--------------  Definition of Global Variables ---------------------
23bankid = 0
24
25datafilename = "" 
26hklfilename = ""
27irffilename = ""
28
29datawsname = ""
30instrparamwsname = ""
31braggpeakparamwsname = ""
32
33# Range for Le Bail Fit of all peaks
34startx = -1
35endx =  -1
36
37# Range for fitting single peaks for step 1~3
38tofmin_singlepeaks = -1
39tofmax_singlepeaks = -1
40
41# Background related
42backgroundtype = "Polynomial"
43backgroundorder = 6
44bkgdtablewsname = ""
45bkgdwsname = ""
46bkgdfilename = ""
47usrbkgdpoints = ''
48
49latticesize = 4.1568899999999998
50#--------------------------------------------------------------------
51
52def setupGlobals(infofilename):
53    """ Set up globals values
54    """
55
56    global datafilename, hklfilename, irffilename 
57    global datawsname, instrparamwsname, braggpeakparamwsname
58    global bkgdtablewsname, bkgdwsname, bkgdfilename, backgroundorder
59    global tofmin_singlepeaks, tofmax_singlepeaks, startx, endx
60    global bankid, latticesize
61    global usrbkgdpoints
62
63    bankid, calibDict = importCalibrationInformation(infofilename)
64    bankid = int(bankid)
65
66    datafilename = calibDict["DataFileDir"] + calibDict[bankid]["DataFileName"] 
67    hklfilename  = calibDict["HKLFileDir"]  + calibDict[bankid]["HKLFileName"]
68    irffilename  = calibDict["IrfFileDir"]  + calibDict[bankid]["IrfFileName"]
69   
70    startx = float(calibDict[bankid]["LeBailFitMinTOF"])
71    endx   = float(calibDict[bankid]["LeBailFitMaxTOF"])
72
73    # Name of workspaces
74    datawsname = calibDict[bankid]["DataWorkspace"]
75    instrparamwsname     = "Bank%sInstrumentParameterTable" % (bankid)
76    braggpeakparamwsname = 'BraggPeakParameterTable'
77
78    # Background related
79    usrbkgdpoints   = calibDict[bankid]["UserSpecifiedBkgdPts"]
80    bkgdwsname      = datawsname+"_Background"
81    backgroundtype  = calibDict["BackgroundType"]
82    backgroundorder = int(calibDict["BackgroundOrder"])
83    bkgdfilename    = calibDict["WorkingDir"] + datawsname + "_Parameters.bak'"
84    bkgdwsname      = datawsname + "_Background"
85    bkgdtablewsname = datawsname + "_Background_Parameters"
86
87    # Other constants
88    latticesize   = calibDict[bankid]["LatticeSize"]
89
90    return
91
92def generateMCSetupTable(wsname):
93    """ Generate a Le Bail fit Monte Carlo random walk setup table
94    """
95    import mantid.simpleapi as api
96   
97    tablews = api.CreateEmptyTableWorkspace(OutputWorkspace=str(wsname))
98   
99    tablews.addColumn("str", "Name")
100    tablews.addColumn("double", "A0")
101    tablews.addColumn("double", "A1")
102    tablews.addColumn("int", "NonNegative")
103    tablews.addColumn("int", "Group")
104   
105    group = 0
106    tablews.addRow(["Dtt1"  , 5.0, 0.0, 0, group]) 
107    tablews.addRow(["Dtt1t" , 5.0, 0.0, 0, group])
108    tablews.addRow(["Dtt2t" , 0.1, 1.0, 0, group])
109    tablews.addRow(["Zero"  , 5.0, 0.0, 0, group])
110    tablews.addRow(["Zerot" , 5.0, 0.0, 0, group])
111    tablews.addRow(["Width" , 0.0, 0.1, 1, group])
112    tablews.addRow(["Tcross", 0.0, 1.0, 1, group])
113   
114    group = 1
115    tablews.addRow(["Beta0" , 0.50, 1.0, 0, group]) 
116    tablews.addRow(["Beta1" , 0.05, 1.0, 0, group])
117    tablews.addRow(["Beta0t", 0.50, 1.0, 0, group])
118    tablews.addRow(["Beta1t", 0.05, 1.0, 0, group])
119   
120    group = 2
121    tablews.addRow(["Alph0" , 0.05, 1.0, 0, group])
122    tablews.addRow(["Alph1" , 0.02, 1.0, 0, group])
123    tablews.addRow(["Alph0t", 0.10, 1.0, 0, group])
124    tablews.addRow(["Alph1t", 0.05, 1.0, 0, group])
125   
126    group = 3
127    tablews.addRow(["Sig0", 2.0, 1.0, 1, group])
128    tablews.addRow(["Sig1", 2.0, 1.0, 1, group])
129    tablews.addRow(["Sig2", 2.0, 1.0, 1, group])
130   
131    group = 4
132    tablews.addRow(["Gam0", 2.0, 1.0, 0, group])
133    tablews.addRow(["Gam1", 2.0, 1.0, 0, group])
134    tablews.addRow(["Gam2", 2.0, 1.0, 0, group])
135
136    return tablews
137
138def breakParametersGroups(tablews):
139    """ Break the parameter groups.  Such that each parameter/row has an individual group
140    """
141    numrows = tablews.rowCount()
142    for ir in xrange(numrows):
143        tablews.setCell(ir, 4, ir)
144
145    return
146
147def resetParametersGroups(tablews):
148    """ Set the group number to original setup
149    """
150    numrows = tablews.rowCount()
151    for ir in xrange(numrows):
152        parname = tablews.cell(ir, 0)
153        if parname in ["Dtt1", "Dtt1t", "Dtt2t", "Zero", "Zerot", "Width", "Tcross"]:
154            group = 0
155        elif parname in ["Beta0", "Beta1", "Beta0t", "Beta1t"]:
156            group = 1
157        elif parname in ["Alph0", "Alph1", "Alph0t", "Alph1t"]:
158            group = 2
159        elif parname in ["Sig0", "Sig1", "Sig2"]:
160            group = 3
161        else:
162            group = 4
163        tablews.setCell(ir, 4, group)
164        return
165
166def doStep4(numsteps):
167    """ Use Monte Carlo random/drunken walk for solution
168            FitRegion                   = '8500,49000',
169            BackgroundParametersWorkspace   ='PG3_10808_Background_Parameters',
170    """ 
171    global datawsname
172    global instrparamwsname
173    global braggpeakparamwsname
174    global startx, endx
175    global backgroundorder, bkgdtablewsname
176
177    outputwsname = '%s_MC_%s' % (datawsname, numsteps)
178    outputwsname2 = '%s_MC_%s' % (datawsname, numprofsteps)
179   
180    print instrparamwsname   
181    print braggpeakparamwsname
182    print startx, endx
183
184    instwsname = 'Bank1InstrumentParametersTable'
185
186    # Set up the parameters to refine
187    UpdatePeakParameterTableValue(
188            InputWorkspace=instwsname,
189            Column="FitOrTie",
190            NewStringValue="tie")
191    UpdatePeakParameterTableValue( 
192            InputWorkspace=instwsname, 
193            Column="FitOrTie",
194            ParameterNames=["Dtt1", "Dtt1t", "Zerot", "Width"],
195            #ParameterNames=["Dtt1", "Dtt1t", "Zerot"],
196            NewStringValue="fit")
197
198    # Limit the range of MC
199    cwl = 0.533
200    UpdatePeakParameterTableValue( 
201            InputWorkspace=instwsname, 
202            Column="Min",
203            ParameterNames=["Width"],
204            NewFloatValue=0.50) #cwl*0.25)
205           
206    UpdatePeakParameterTableValue( 
207            InputWorkspace=instwsname, 
208            Column="Max",
209            ParameterNames=["Width"],
210            NewFloatValue=1.25) #cwl*4.0)
211           
212    wsname = "MCSetupParameterTable"
213    tablews = generateMCSetupTable(wsname)
214    #breakParametersGroups(tablews)
215
216    print "Fit range: %f , %f" % (startx, endx)
217    LeBailFit(
218            InputWorkspace                  = datawsname,
219            OutputWorkspace                 = outputwsname,
220            InputParameterWorkspace         = instwsname,
221            OutputParameterWorkspace        = 'Bank1RefinedProfileTable',
222            InputHKLWorkspace               = 'BraggPeakParameterTable1',
223            OutputPeaksWorkspace            = 'BraggPeakRefinedParameterTable',
224            FitRegion                       = '%f, %f' % (startx, endx),
225            Function                        = 'MonteCarlo', 
226            NumberMinimizeSteps             = numsteps, 
227            BackgroundParametersWorkspace   = bkgdtablewsname,
228            UseInputPeakHeights             = False, 
229            PeakRadius                      ='8',
230            Minimizer                       = 'Levenberg-Marquardt',
231            MCSetupWorkspace        = str(wsname),
232            Damping                         = '5.0',
233            RandomSeed                      = 0,
234            AnnealingTemperature            = 100.0,
235            DrunkenWalk                     = True)         
236   
237    # Set up the parameters to refine
238   
239    instwsname = str('Bank1RefinedProfileTable')
240    UpdatePeakParameterTableValue(
241            InputWorkspace=instwsname,
242            Column="FitOrTie",
243            NewStringValue="tie")
244    UpdatePeakParameterTableValue( 
245            InputWorkspace=instwsname, 
246            Column="FitOrTie",
247            ParameterNames=["Alph0", "Alph0t", "Alph1", "Beta0", "Beta1", "Beta0t", "Beta1t", "Sig1", "Sig2"],
248            NewStringValue="fit")
249           
250    #breakParametersGroups(tablews)         
251   
252    print "Fit range: %f , %f" % (startx, endx)
253    LeBailFit(
254            InputWorkspace                  = datawsname,
255            OutputWorkspace                 = outputwsname2,
256            InputParameterWorkspace         = instwsname,
257            OutputParameterWorkspace        = 'Bank1RefinedProfileTable2',
258            InputHKLWorkspace               = 'BraggPeakParameterTable1',
259            OutputPeaksWorkspace            = 'BraggPeakParameterTable4',
260            FitRegion                       = '%f, %f' % (startx, endx),
261            Function                        = 'MonteCarlo', 
262            NumberMinimizeSteps             = numprofsteps/10, 
263            BackgroundParametersWorkspace   = bkgdtablewsname,
264            UseInputPeakHeights             = False, 
265            PeakRadius                      ='8',
266            Minimizer                       = 'Levenberg-Marquardt',
267            MCSetupWorkspace        = str(wsname),
268            Damping                         = '1.0',
269            RandomSeed                      = 0,
270            AnnealingTemperature            = 100.0,
271            DrunkenWalk                     = True)
272    return
273
274
275def main(argv):
276    """ Main
277    """
278    global numsteps
279    if numsteps <= 0:
280            numsteps = 200
281            print "Using default number of steps 200"
282
283    setupGlobals("/home/wzz/Projects/MantidTests/LeBailFit/Test2013B/Bank1/Calibration_Information.config")
284
285    doStep4(numsteps)
286
287    return
288
289if __name__=="__main__":
290    main([])