source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MC_SimulationScripting.ipf @ 950

Last change on this file since 950 was 950, checked in by srkline, 8 years ago

adding new procedures for testing -- that will enable:

(1) automation of SANS data reduction, at least in some of the more standard cases. The more consistently the data files are named, the better this works.
(2) Much simpler scripting of experiment simulation. now a simulated experiment can be set up is the same way that a real experiemtn can -- by setting up a list of "runs"

Help files are to follow for all of these. Prelimiary help for the automation has been added to the SANS Reduction Help file. Videos to follow. Loaders for these two items have been added to the Macros menu so that all of the dependencies are satisfied.

More testing is still necessary to make sure that nothing has been broken, and that sufficient error catching has been done so that meaningful testing can be done.

File size: 34.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4
5
6//************
7// un-comment the Macros menu declaration to provide easier access to the functions
8// -- remove the "x" from "xMenu" in the section following these comments.  The macros menu
9//    will then have easier access to the scripting functions.
10//
11// **To write your own scripts, follow the examples in:
12//              Example_1DSim()
13//              Example_2DSim()
14//              Example_Loop_1DSim()
15//              Example_Loop_2DSim()
16//
17// and also the instructions below in the "basic cycle"
18//
19//************
20
21
22
23// this is a file with all of the functions that can be used to
24// script the MC simulation of data
25//
26// there are a lot of utility functions to "move" the parts of an instrument that
27// are "moveable" into different configurations. The aim is to have this appear as
28// much like running an experiment as possible.
29//
30//
31//  The basic cycle:
32//
33// 1) plot the model function, and set the coefficients to what you want
34//
35// 2) set (and save) as many configurations as you like -picking ones that are most appropriate for the sample scattering
36//
37//              Sim_SaveConfiguration(waveStr)
38//
39// 3) simulate the sample scattering at your configurations for specified times. The 1D count rate
40//               should match fairly well to the 2D count rate. (Check this for low/high incoherent levels)
41//-or-
42// 4) run everything for a short time (say 1s to 10s) to get the count rates
43//  - this will actually save the files - just like a real experiment. This is like step 6, but with short count times.
44//
45// 5) with the sample count rates, find the optimal count times
46//
47//              OptimalCount(samCR,empCR)
48//
49// 6) Write a function to "run" the samples at configurations and count times. Don't forget to simulate
50//               the transmissions. (See the examples below). Copy the function to a new procedure window, rename it,
51//              and edit it there as you wish.
52//
53//              ListSASCALCConfigs() will print out the list of configurations
54//
55// 7) Then with the optimal count times, do a dry run to see how long the whole simulation would take.
56//
57//              Sim_DryRun(funcStr)
58//
59// 8) Then if that's OK, run the full simulation
60//
61//              Your function here()
62//
63// 9) If it's 1D data, use NSORT to combine. If it's 2D data, reduce as usual. Be sure to skip the detector
64//              efficiency correction, and set the dead time to a tiny value so that it won't have any
65//    effect (since the detector CR is fictionally high, and the detector is perfect):
66//
67//              Sim_SetDeadTimeTiny()
68//
69// When you have a lot of 1D waves to combine, and they are not numbered like in a real reduction
70// experiment, see:
71//
72//  MakeCombineTable_byName()
73//  DoCombineFiles_byName(lowQfile,medQfile,hiQfile,saveName) (in NSORT.ipf)
74//
75// this works with 1, 2, or 3 data files
76//
77
78
79// In general, when setting up a simulation, it's easier to set up sample conditions for a particular
80// sample, and loop through the configurations. If you want your 2D data to "look" like a typical
81// experiment, then you'll need to simulate each different sample at one configuration, then "move"
82// to a different configuration and loop through the samples again. Somewhat more cumbersome, all to
83// get the file catalog to be "grouped" like a real SANS experiment.
84
85
86//
87// Important - in the reduction process, set the  dead time correction to something
88// really small, say 10^-15
89//
90// see Sim_SetDeadTimeTiny()
91
92//
93// TODO:
94// x- need to update the 1D writer to allow a name to be passed in, rather than
95// always using the dialog
96// x- fill in the "RunSample" functions
97//
98//  x- I need a little panel to control all of this, and get the information, just
99//     like setting up a real experiment. Or maybe not. Maybe better to keep some of this
100//     hidden.
101//  -- It would be nice to be able to automatically do the sum_model to add the 2D empty cell contribution
102//     For some experiments it's quite useful to see the level of this contribution, rather than
103//     the completely clean simulation.
104//  x- in this file, clearly separate the examples from the utility functions
105//  -- get a function list of everything, and document what each does (if not obvious)
106//  -- step by step comments of one of the examples (1D and 2D) are probably the most useful
107//  x- "dry run" function to estimate the total simulation time for 2D. (? set all monitor counts to 1000
108//      or whatever I use to get other estimates, and then run through everything to get estimated times)
109//
110// set the global:
111//      NVAR g_estimateOnly = root:Packages:NIST:SAS:g_estimateOnly             // == 1 for just a time estimate, == 0 (default) to do the simulation
112// and then do the dry run
113//
114//
115// ----- transmission really in the simulation...
116//  x- add function to simulate transmission in 2D by saving a "Trans" configuration
117//     and automatically using that. Currently the trans is already written to the file
118//  x- would also need a function to simulate an empty beam measurement to rescale the
119//     transmissions. (use the model given, set scale (w[0] to something really small)
120//  x- or better, use the empty cell function that is guaranteed to have coefficients that are
121//     well behaved - and I can set/reset them to get a proper "empty beam"
122//
123//
124// -- scattering from an empty cell is NOT provided in 2D
125//
126
127
128
129/// un-comment this by removing the "x" from the word xMenu. Then compile.
130// - note that no menu items will appear unless BOTH the Reduction and Analysis are loaded
131Menu "Macros"
132                Submenu "Simulation Scripting - Beta"
133                        ScriptItem(0),ShowRunListPanel()
134//                      ScriptItem(0),Sim_StoreConfProc()
135//                      ScriptItem(1),Sim_moveConfProc()
136                        ScriptItem(1)
137//                      ScriptItem(2),ListSASCALCConfigs()
138                        ScriptItem(3),DryRunProc_1D()
139//                      ScriptItem(4),DryRunProc_2D()
140                        ScriptItem(5),OptimalCountProc()
141                        ScriptItem(6),MakeCombineTable_byName()
142                        ScriptItem(7),DoCombineFiles_byName(lowQfile,medQfile,hiQfile,saveName)
143                        ScriptItem(8),Sim_SetDeadTimeTiny()
144                        ScriptItem(9)
145                        ScriptItem(10),Setup_Sim_Example()
146//                      ScriptItem(11),Example_1DSim()
147//                      ScriptItem(12),Example_2DSim()
148                        ScriptItem(13)
149                        ScriptItem(14),DisplayProcedure "Example_1DSim"
150                End     
151End
152
153Function/S ScriptItem(num)
154        Variable num
155       
156        String str=""
157       
158        if(exists("root:SANS_RED_VERSION") && exists("root:Packages:NIST:SANS_ANA_VERSION"))
159                switch(num)     
160                        case 0:
161                                str = "Run List Panel"
162                                //str = "Save Configuration"
163                                break
164                        case 1:
165                                //str = "Move to Configuration"
166                                str = "-"
167                                break
168                        case 2:
169                                str = "List Configurations"
170                                break
171                        case 3:
172                                str = "1D Count Rates"
173                                break
174                        case 4:
175                                str = "2D Dry Run"
176                                break
177                        case 5:
178                                str = "Optimal Count Times"
179                                break
180                        case 6:
181                                str = "Make Table to Combine By Name"
182                                break
183                        case 7:
184                                str = "Combine by Name"
185                                break
186                        case 8:
187                                str = "Turn Off Dead Time Correction"
188                                break
189                        case 9:
190                                str = "-"
191                                break
192                        case 10:
193                                str = "Setup Sim Example"
194                                break
195                        case 11:
196                                str = "Run 1D Sim Example"
197                                break
198                        case 12:
199                                str = "Run 2D Sim Example"
200                                break
201                        case 13:
202                                str = "-"
203                                break
204                        case 14:
205                                str = "Display Example Code"
206                                break   
207                endswitch
208        endif
209       
210        return(str)
211end
212
213////////// --- START OF EXAMPLE SCRIPTS ---  ////////////////
214
215//
216// run this before either the 1D or 2D example to make sure that the proper named configurations and function exist.
217// this function will overwrite any same-named configurations
218//
219Function Setup_Sim_Example()
220
221        SetDataFolder root:Packages:NIST:SAS:
222       
223        Make/O/T/N=20 Config_1m,Config_4m,Config_13m
224        Config_1m = {"checkNG7","8","5.08 cm","6","0.115","checkChamber","1/2\"","0","100","25","","","","","","","","","",""}
225        Config_4m = {"checkNG7","5","5.08 cm","6","0.115","checkChamber","1/2\"","0","400","0","","","","","","","","","",""}
226        Config_13m = {"checkNG7","1","5.08 cm","6","0.115","checkChamber","1/2\"","0","1300","0","","","","","","","","","",""}
227 
228// include the model and plot it, so that it will exist. Post to queue so they execute in order
229        Execute/P "INSERTINCLUDE \"SchulzSpheres_Sq_v40\""
230        Execute/P "INSERTINCLUDE \"DAB_Model_v40\""
231        Execute/P "COMPILEPROCEDURES "
232        Execute/P "PlotSchulzSpheres_SC(256,0.001,0.7)"
233        Execute/P "PlotDAB_Model(256,0.001,0.7)"
234       
235        Execute/P "SASCALC()"
236                               
237        SetDataFolder root:
238
239End
240
241Function Example_1DSim()
242
243        String confList,ctTimeList,saveNameList,funcStr,titleStrList
244       
245
246        Sim_SetSimulationType(0)                //kill the simulation panel
247        Sim_SetSimulationType(1)                //open the 1D simulation panel
248       
249//(1)   enter the (unsmeared) function name
250        funcStr = "SchulzSpheres_SC"
251        Wave cw = $("root:"+getFunctionCoef(funcStr))
252        Sim_SetModelFunction(funcStr)                                           // model function name
253               
254//(2) model coefficients here, if needed. Wave name is "cw"
255//   then set the sample thickness. Be sure to use an appropriate incoherent background in the model
256        Sim_Set1D_Transmission(0.75)                                            // For 1D, I need to supply the transmission
257        Sim_SetThickness(0.2)                                                           // thickness
258        Sim_Set1D_ABS(1)                                                                                // absolute scaling (1== yes)
259        Sim_Set1D_Noise(1)                                                                      // noise (1== yes, add statistical noise)
260
261        cw[0] = 0.1
262        // as needed - look at the parameter list for the model
263       
264       
265//(3) set the configuration list, times, and saved names
266// -- the mumber of listed configurations must match the number of discrete count times and save names
267// titleStr is the label and is the same for each run of the same sample
268        confList = "Config_1m;Config_4m;Config_13m;"
269        ctTimeList = "100;300;900;"
270        saveNameList = "sim_1m.abs;sim_4m.abs;sim_13m.abs;"
271        titleStrList = "MySample 1 at 1m;MySample 1 at 4m;MySample 1 at 13m;"
272
273        // then this runs the samples as listed
274        Sim_RunSample_1D(confList,ctTimeList,titleStrList,saveNameList)
275
276        // no transmissions or empty beam measurements to make for 1D simulation
277
278        return(0)       
279End
280
281//
282Function Example_2DSim()
283
284        String confList,ctTimeList,titleStrList,transConfList,transCtTimeList
285        Variable runIndex,val,totalTime
286        String funcStr
287
288tic()
289
290        Sim_SetSimulationType(0)                //kill the simulation panel
291        Sim_SetSimulationType(2)                //open the 2D simulation panel
292        Sim_SetSimTimeWarning(36000)                    //sets the threshold for the warning dialog to 10 hours
293        totalTime = 0
294
295
296//(1)   enter the (unsmeared) function name
297        funcStr = "SchulzSpheres_SC"
298        Wave cw = $("root:"+getFunctionCoef(funcStr))
299        Sim_SetModelFunction(funcStr)                                           // model function name
300
301//(2) set the standard sample cell size (1" diam banjo cell)
302// and set the conditions for beam stop in, and raw counts
303        Sim_SetSampleRadius(1.27)                                                       // sam radius (cm)
304        Sim_SetRawCountsCheck(1)                                                        // raw cts? 1== yes
305        Sim_SetBeamStopInOut(1)                                                         // BS in? 1==yes
306
307//(3) model coefficients here, if needed. Wave name is "cw"
308//   then set the sample thickness and incoherent cross section
309
310        cw[0] = 0.1
311        // as needed - look at the parameter list for the model
312
313        Sim_SetThickness(0.2)                                                           // thickness (cm)
314        Sim_SetIncohXS(1.3)                                                                     // incoh XS
315       
316//(4) starting run index for the saved raw data files. this will automatically increment
317//    as the sample is "Run"
318        runIndex = 400
319       
320//(5) set the configuration list, times, a single sample label, and the starting run index
321// -- the mumber of listed configurations must match the number of discrete count times
322        confList = "Config_1m;Config_4m;Config_13m;"
323        ctTimeList = "100;300;900;"
324        transConfList = "Config_13m;"           // trans only @ 13m
325        transCtTimeList = "1;"                          // trans count time = 1s
326        titleStrList = "MySample 1 at 1m;MySample 1 at 4m;MySample 1 at 13m;"
327
328        // runIndex is PBR and updates as the number of files are written
329       
330        // this runs the sample
331        totalTime += Sim_RunSample_2D(confList,ctTimeList,titleStrList,runIndex)
332        // this runs the transmissions
333        totalTime += Sim_RunTrans_2D(transConfList,transCtTimeList,titleStrList,runIndex)
334       
335        // run the empty beam at all configurations
336        // This will automatically change the function to "EC_Empirical" and "empty beam" conditions
337        transConfList = "Config_1m;Config_4m;Config_13m;"
338        transCtTimeList = "1;1;1;"     
339        titleStrList = "Empty Beam at 1m;Empty Beam at 4m;Empty Beam at 13m;"
340        totalTime += Sim_RunEmptyBeamTrans_2D(transConfList,transCtTimeList,titleStrList,runIndex)
341       
342        Print "runIndex = ",runIndex
343
344        Sim_SetSimTimeWarning(10)
345
346toc()
347       
348        return(totalTime)
349End
350
351// this example was really used, run repeatedly with different count times to get
352// replicate data sets to test the overlap
353Function Example_Loop_1DSim()
354
355        String confList,ctTimeList,saveNameList,funcStr,titleStrList
356       
357
358        Sim_SetSimulationType(0)                //kill the simulation panel
359        Sim_SetSimulationType(1)                //open the 1D simulation panel
360       
361//(1)   enter the (unsmeared) function name
362        funcStr = "DAB_model"
363        Wave cw = $("root:"+getFunctionCoef(funcStr))
364        Sim_SetModelFunction(funcStr)                                           // model function name
365               
366//(2) model coefficients here, if needed. Wave name is "cw"
367//   then set the sample thickness. Be sure to use an appropriate incoherent background in the model
368        Sim_Set1D_Transmission(0.8)                                             // For 1D, I need to supply the transmission
369        Sim_SetThickness(0.1)                                                           // thickness
370        Sim_Set1D_ABS(1)                                                                                // absolute scaling (1== yes)
371        Sim_Set1D_Noise(1)                                                                      // noise (1== yes, add statistical noise)
372
373
374        cw = {1e-05,200,0}
375       
376       
377//(3) set the configuration list, times, and saved names
378// -- the mumber of listed configurations must match the number of discrete count times and save names
379// titleStr is the label and is the same for each run of the same sample
380        confList = "Config_4m;"
381        ctTimeList = ""         //filled in the loop
382        saveNameList = ""
383        titleStrList = "DAB versus count time t = "             //overwritten in loop
384
385        Make/O/D ctTime = {5,11,16,21,27,32,37,43,48,53,107,160,214,267,321,374,428,481,535,1604,5348,21390,53476}
386
387        Variable ii=0,nTrials=10,jj
388       
389        for(jj=0;jj<numpnts(ctTime);jj+=1)     
390                for(ii=0;ii<nTrials;ii+=1)
391                        titleStrList = "DAB versus count time t = "+num2str(ctTime[jj])
392                        saveNameList = "DAB_4m_t"+num2str(ctTime[jj])+"_"+num2str(ii)+".abs;"
393                        ctTimeList = num2str(ctTime[jj])+";"
394                        // then this runs the samples as listed
395                        Sim_RunSample_1D(confList,ctTimeList,titleStrList,saveNameList)
396                endfor
397        endfor
398
399        // no transmissions or empty beam measurements to make for 1D simulation
400
401        return(0)       
402End
403
404
405//
406//
407// This example will run the same sample with three different thicknesses, at
408// 1m, 4m, 13m
409//
410// empty beam measurements at all three distances, sample transmission at 13m
411//
412// total simulation time is < 600 seconds on my old machine...
413// do a dry run first to see how long it'll take.
414//
415//
416Function Example_Loop_2DSim()
417
418        String confList,ctTimeList,titleStrList,transConfList,transCtTimeList
419        Variable runIndex,val,totalTime
420        String funcStr
421
422tic()
423
424        Sim_SetSimulationType(0)                //kill the simulation panel
425        Sim_SetSimulationType(2)                //open the 2D simulation panel
426        Sim_SetSimTimeWarning(36000)                    //sets the threshold for the warning dialog to 10 hours
427        totalTime = 0
428
429
430//(1)   determine the (unsmeared) function name (we'll set this right before the simulation)
431        funcStr = "DAB_model"
432        Wave cw = $("root:"+getFunctionCoef(funcStr))
433
434//(2) set the standard sample cell size (1" diam banjo cell)
435// and set the conditions for beam stop in, and raw counts
436        Sim_SetSampleRadius(1.27)                                                       // sam radius (cm)
437        Sim_SetRawCountsCheck(1)                                                        // raw cts? 1== yes
438        Sim_SetBeamStopInOut(1)                                                         // BS in? 1==yes
439
440//(3) model coefficients here, if needed. Wave name is "cw"
441//   then set the sample thickness and incoherent cross section
442
443        cw = {1e-05,200,0.1}
444       
445        // as needed - look at the parameter list for the model
446
447        Sim_SetThickness(0.2)                                                           // thickness (cm)
448        Sim_SetIncohXS(1.3)                                                                     // incoh XS
449       
450//(4) starting run index for the saved raw data files. this will automatically increment
451//    as the sample is "Run"
452        runIndex = 500
453
454
455//(5)  run the transmissions and empty beam first, before you forget them
456
457        // run the empty beam at all configurations
458        // This will automatically change the function to "EC_Empirical" and "empty beam" conditions
459        transConfList = "Config_1m;Config_4m;Config_13m;"
460        transCtTimeList = "1;1;1;"     
461        titleStrList = "Empty Beam at 1m;Empty Beam at 4m;Empty Beam at 13m;"
462        totalTime += Sim_RunEmptyBeamTrans_2D(transConfList,transCtTimeList,titleStrList,runIndex)
463       
464       
465//(6) set the configuration list, times, a single sample label, and the starting run index
466// -- the mumber of listed configurations must match the number of discrete count times
467        confList = ""                   // these will be filled in the loop
468        ctTimeList = ""
469        transConfList = "Config_13m"            // trans only @ 13m
470        transCtTimeList = "1;"                          // trans count time = 1s
471        titleStrList = "MySample 1"
472
473        // runIndex is PBR and updates as the number of files are written
474
475        Variable ii,jj
476
477//      any, all, or more settings can be set up to change in the loop
478//  -- be sure these waves are the same length and the values correspond.
479        Make/O/D thick = {0.1,0.2,0.5,0.1,0.2,0.5,0.1,0.2,0.5}
480        Make/O/D ctTime = {100,100,100,300,300,300,900,900,900}
481        Make/O/D/T conf = {"Config_1m","Config_1m","Config_1m","Config_4m","Config_4m","Config_4m","Config_13m","Config_13m","Config_13m"}
482
483
484        Sim_SetModelFunction(funcStr)                                           // model function name
485       
486        for(ii=0;ii<numpnts(thick);ii+=1)
487                Sim_SetThickness(thick[ii])                                                             // thickness (cm)
488
489                confList = conf[ii] +";"
490                titleStrList = "DAB simulation, thick = "+num2str(thick[ii])
491                ctTimeList = num2str(ctTime[ii])+";"
492               
493                // this runs the transmissions (only at 13m)
494                if(cmpstr(conf[ii],"Config_13m")==0)
495                        totalTime += Sim_RunTrans_2D(transConfList,transCtTimeList,titleStrList,runIndex)
496                endif
497               
498                // this runs the sample
499                totalTime += Sim_RunSample_2D(confList,ctTimeList,titleStrList,runIndex)
500
501        endfor
502       
503        Print "runIndex = ",runIndex
504
505        Sim_SetSimTimeWarning(10)
506
507toc()
508       
509        return(totalTime)
510End
511
512
513
514
515
516//////////////// ---- END OF EXAMPLE SCRIPTS --- ////////////////////
517
518
519
520
521
522// pass in a semicolon delimited list of configurations + corresponding count times + saved names
523Function Sim_RunSample_1D(confList,ctTimeList,titleStrList,saveNameList)
524        String confList,ctTimeList,titleStrList,saveNameList
525       
526        Variable ii,num,ct,cr,numPt
527        String twStr,fname,type,titleStr
528        NVAR g_estimateOnly = root:Packages:NIST:SAS:g_estimateOnly             // == 1 for just count rate, == 0 (default) to do the simulation and save
529        WAVE/Z crWave = root:CR_1D
530        WAVE/Z/T fileWave = root:Files_1D
531
532        type = "SAS"            // since this is a simulation
533        num=ItemsInList(confList)
534       
535        for(ii=0;ii<num;ii+=1)
536                twStr = StringFromList(ii, confList,";")
537                Wave/T tw=$("root:Packages:NIST:SAS:"+twStr)
538                ct = str2num(StringFromList(ii, ctTimeList,";"))
539                fname = StringFromList(ii, saveNameList,";")
540                titleStr = StringFromList(ii, titleStrList,";")
541               
542                Sim_MoveToConfiguration(tw)
543                Sim_SetCountTime(ct)
544                cr = Sim_Do1DSimulation()
545               
546                // either save it out, or return a table of the count rates
547                if(g_estimateOnly)
548                        numPt = numpnts(crWave)
549                        InsertPoints numPt, 1, crWave,fileWave
550                        crWave[numPt] = cr
551                        fileWave[numPt] = fname
552                else
553                        Sim_Save1D_wName(type,titleStr,fname)           //this function will increment the runIndex
554                endif
555       
556        endfor
557
558        return(0)
559End
560
561// fakes an empty beam transmission by setting the conditions of the
562// empty cell temporarily to something that doesn't scatter
563Function Sim_RunEmptyBeamTrans_2D(confList,ctTimeList,titleStr,runIndex)
564        String confList,ctTimeList,titleStr
565        Variable &runIndex
566       
567        Wave cw = root:Packages:NIST:SAS:coef_ECEmp             //EC coefs are in the SAS folder, not root:
568        Variable totalTime
569       
570        // change the function and coefficients
571        Sim_SetModelFunction("EC_Empirical")
572        cw = {1e-18,3.346,0,9,0}
573        Sim_SetIncohXS(0)                                                                       // incoh XS = 0 for empty beam
574        Sim_SetThickness(0.1)                                                   // thickness (cm) to give proper level of scattering
575
576               
577        totalTime = Sim_RunTrans_2D(confList,ctTimeList,titleStr,runIndex)
578               
579        // change them back
580        cw = {2.2e-08,3.346,0.0065,9,0.016}
581       
582        return(totalTime)
583End
584
585 
586// should be called with each list having only one item
587// and a count time of 1 s
588//
589// -- for an empty beam, before calling, set the incoh XS = 0 and set the scale
590// of the model to something tiny so that there is no coherent scattering
591//
592Function Sim_RunTrans_2D(confList,ctTimeList,titleStrList,runIndex)
593        String confList,ctTimeList,titleStrList
594        Variable &runIndex
595       
596        Variable ii,num,ct,index,totalTime
597        String twStr,type,titleStr
598
599        NVAR g_estimateOnly = root:Packages:NIST:SAS:g_estimateOnly             // == 1 for just a time estimate, == 0 (default) to do the simulation
600        NVAR g_estimatedMCTime = root:Packages:NIST:SAS:g_estimatedMCTime               // estimated MC sim time
601        totalTime = 0
602               
603        type = "SAS"            // since this is a simulation
604        num=ItemsInList(confList)
605       
606        for(ii=0;ii<num;ii+=1)
607                twStr = StringFromList(ii, confList,";")
608                titleStr = StringFromList(ii, titleStrList,";")
609                Wave/T tw=$("root:Packages:NIST:SAS:"+twStr)
610                ct = str2num(StringFromList(ii, ctTimeList,";"))
611               
612                Sim_MoveToConfiguration(tw)
613                Sim_SetBeamStopInOut(0)                                                         // BS out for Trans
614                Sim_SetCountTime(ct)
615                Sim_DoMCSimulation()
616               
617                if(g_estimateOnly)
618                        totalTime += g_estimatedMCTime          //      don't save, don't increment. time passed back as global after each MC simulation
619                else
620                        SaveAsVAXButtonProc("",runIndex=runIndex,simLabel=(titleStr+" at "+twStr))
621                        runIndex += 1
622                endif
623
624        endfor
625       
626        Sim_SetBeamStopInOut(1)                                                         // put the BS back in
627
628        return(totalTime)
629End
630
631
632
633// if just an estimate, I can also get the count rate too
634// WAVE results = root:Packages:NIST:SAS:results
635//      Print "Sample Simulation (2D) CR = ",results[9]/ctTime
636// -- for estimates, iMon is set to 1000, so time=1000/(root:Packages:NIST:SAS:gImon)
637Function Sim_RunSample_2D(confList,ctTimeList,titleStrList,runIndex)
638        String confList,ctTimeList,titleStrList
639        Variable &runIndex
640
641        Variable ii,num,ct,index,totalTime
642        String twStr,type,titleStr
643
644        NVAR g_estimateOnly = root:Packages:NIST:SAS:g_estimateOnly             // == 1 for just a time estimate, == 0 (default) to do the simulation
645        NVAR g_estimatedMCTime = root:Packages:NIST:SAS:g_estimatedMCTime               // estimated MC sim time
646        totalTime = 0
647       
648        type = "SAS"            // since this is a simulation
649        num=ItemsInList(confList)
650       
651       
652        for(ii=0;ii<num;ii+=1)
653                twStr = StringFromList(ii, confList,";")
654                titleStr = StringFromList(ii, titleStrList,";")
655                Wave/T tw=$("root:Packages:NIST:SAS:"+twStr)
656                ct = str2num(StringFromList(ii, ctTimeList,";"))
657               
658                Sim_MoveToConfiguration(tw)
659                Sim_SetBeamStopInOut(1)                                                         // BS in?
660                Sim_SetCountTime(ct)
661                Sim_DoMCSimulation()
662               
663                if(g_estimateOnly)
664                        totalTime += g_estimatedMCTime          //      don't save, don't increment. time passed back as global after each MC simulation
665                else
666                        SaveAsVAXButtonProc("",runIndex=runIndex,simLabel=(titleStr+" at "+twStr))
667                        runIndex += 1
668                endif
669               
670        endfor
671
672        return(totalTime)
673End
674
675// a prototype function with no parameters
676// for a real experiment - it must return the total time, if you want a proper estimate
677Function Sim_Expt_Proto()
678        Print "In the Sim_Expt_Proto -- error"
679        return(0)
680End
681
682Proc DryRunProc_2D(funcStr)
683        String funcStr
684        Sim_2DDryRun(funcStr)
685end
686
687// pass the function string, no parameters
688Function Sim_2DDryRun(funcStr)
689        String funcStr
690       
691        FUNCREF Sim_Expt_Proto func=$funcStr
692       
693        NVAR g_estimateOnly = root:Packages:NIST:SAS:g_estimateOnly
694        g_estimateOnly = 1
695       
696        Variable totalTime
697//      totalTime = Sim_CountTime_2D()  // your "exeriment" here, returning the totalTime
698       
699        totalTime = func()
700        g_estimateOnly = 0
701       
702        Printf "Total Estimated Time = %g s or %g h\r",totalTime,totalTime/3600
703end
704
705Proc DryRunProc_1D(funcStr)
706        String funcStr
707        Sim_1DDryRun(funcStr)
708end
709
710// pass the function string, no parameters
711// makes a (new) table of the files and CR
712Function Sim_1DDryRun(funcStr)
713        String funcStr
714       
715        FUNCREF Sim_Expt_Proto func=$funcStr
716       
717        NVAR g_estimateOnly = root:Packages:NIST:SAS:g_estimateOnly
718        g_estimateOnly = 1
719       
720        Variable totalTime
721       
722        Make/O/D/N=0 root:CR_1D
723        Make/O/T/N=0 root:Files_1D
724       
725        totalTime = func()
726        g_estimateOnly = 0
727       
728        Edit Files_1D,CR_1D
729       
730//      Printf "Total Estimated Time = %g s or %g h\r",totalTime,totalTime/3600
731        return(0)
732end
733
734
735
736
737
738Proc OptimalCountProc(samCountRate,emptyCountRate)
739        Variable samCountRate,emptyCountRate
740        OptimalCount(samCountRate,emptyCountRate)
741End
742
743Function OptimalCount(samCR,empCR)
744        Variable samCR,empCR
745       
746        String str
747        Variable ratio,target
748        ratio = sqrt(samCR/empCR)
749       
750        str = "For %g Sample counts, t sam = %g s and t emp = %g s\r"
751       
752        target = 1e6
753        printf str,target,round(target/samCR),round((target/samCR)/ratio)
754       
755        target = 1e5
756        printf str,target,round(target/samCR),round((target/samCR)/ratio)       
757
758        target = 1e4
759        printf str,target,round(target/samCR),round((target/samCR)/ratio)       
760        return(0)
761end
762
763
764//
765// Important - in the reduction process, set the dead time correction to something
766// really small, say 10^-15 so that  it won't have any effect (since the detector CR is fictionally high)
767//
768Function Sim_SetDeadTimeTiny()
769
770//      NVAR DeadtimeNG3_ILL = root:myGlobals:DeadtimeNG3_ILL           //pixel resolution in cm
771//      NVAR DeadtimeNG5_ILL = root:myGlobals:DeadtimeNG5_ILL
772//      NVAR DeadtimeNG7_ILL = root:myGlobals:DeadtimeNG7_ILL
773//      NVAR DeadtimeNGB_ILL = root:myGlobals:DeadtimeNGB_ILL
774        NVAR DeadtimeNG3_ORNL_VAX = root:myGlobals:DeadtimeNG3_ORNL_VAX
775//      NVAR DeadtimeNG3_ORNL_ICE = root:myGlobals:DeadtimeNG3_ORNL_ICE
776//      NVAR DeadtimeNG5_ORNL = root:myGlobals:DeadtimeNG5_ORNL
777        NVAR DeadtimeNG7_ORNL_VAX = root:myGlobals:DeadtimeNG7_ORNL_VAX
778//      NVAR DeadtimeNG7_ORNL_ICE = root:myGlobals:DeadtimeNG7_ORNL_ICE
779        NVAR DeadtimeNGB_ORNL_ICE = root:myGlobals:DeadtimeNGB_ORNL_ICE
780        NVAR DeadtimeDefault = root:myGlobals:DeadtimeDefault
781
782        DeadtimeNG3_ORNL_VAX = 1e-15
783        DeadtimeNG7_ORNL_VAX = 1e-15
784        DeadtimeNGB_ORNL_ICE = 1e-15
785        DeadtimeDefault = 1e-15
786       
787        return(0)
788End
789
790
791
792/////////////////////////////////////////////////////////////////////
793///////////// All of the utility functions to make things "move"
794
795
796
797// set the wavelength
798Function Sim_SetLambda(val)
799        Variable val
800
801        NVAR lam = root:Packages:NIST:SAS:gLambda
802        lam=val
803        LambdaSetVarProc("",val,num2str(val),"")
804        return(0)
805end
806
807// set the wavelength spread
808// ?? what are the allowable choices
809// there are 3 choices for each instrument
810//
811// TO DO ?? do this in a better way
812// currently there is no error checking...
813//
814Function Sim_SetDeltaLambda(strVal)
815        String strVal
816       
817       
818        DeltaLambdaPopMenuProc("",1,strVal)             // recalculates intensity if 2nd param==1, skip if == 0
819        PopupMenu popup0_2 win=SASCALC,popmatch=strVal
820//      ControlUpdate/W=SASCALC popup0_2
821       
822        return(0)
823End
824
825// if state == 0, lenses out (box un-checked, no parameters change)
826// if state == 1, lenses in, changes guides, aperture, SDD and wavelength appropriately
827// if prisms desired, follow this call with a set wavelength=17.2
828Function Sim_SetLenses(state)
829        Variable state
830       
831        LensCheckProc("",state)
832        CheckBox checkLens,win=SASCALC,value=state
833       
834        return(0)
835End
836
837// instrName = "checkCGB" or "checkNG7" or "checkNGB" (checkNG3 has been removed)
838// these are the only allowable choices
839Function Sim_SetInstrument(instrName)
840        String instrName
841       
842        SelectInstrumentCheckProc(instrName,0)
843        return(0)
844End
845
846
847//two steps to change the number of guides
848// - first the number of guides, then the source aperture (if 0g, 7g where there are choices)
849//
850// strVal is "5 cm" or "1.43 cm" (MUST be identical to the popup string, MUST have the units)
851///
852Function Sim_SetNumGuides_SrcAp(ng,strVal)
853        Variable ng
854        String strVal
855
856        NVAR gNg=root:Packages:NIST:SAS:gNg
857        gNg = ng
858       
859        Slider SC_Slider,win=SASCALC,variable=root:Packages:NIST:SAS:gNg                                //guides
860
861        // with the new number of guides, update the allowable source aperture popups
862        DoWindow/F SASCALC
863        UpdateControls()
864       
865        // TO DO -- this needs some logic added, or an additional parameter to properly set the
866        // source aperture
867       
868        popupmenu popup0 win=SASCALC,popmatch=strVal
869        Variable dum=sourceApertureDiam()               //sets the proper value in the header wave
870        return(0)
871End
872
873// strVal is "1/2\"" or similar. note the extra backslash so that " is part of the string
874//
875Function Sim_SetSampleApertureDiam(strVal)
876        String strVal
877       
878        popupmenu popup0_1 win=SASCALC,popmatch=strVal
879        Variable dum=sampleApertureDiam()               //sets the proper value in the header wave
880        return(0)
881End     
882
883Function Sim_SetOffset(offset)
884        Variable offset
885       
886        NVAR detOffset = root:Packages:NIST:SAS:gOffset
887        detOffset = offset      //in cm
888        detectorOffset()                //changes header and estimates (x,y) beamcenter
889        return(0)
890end
891
892Function Sim_SetSDD(sdd)
893        Variable sdd
894       
895        NVAR detDist = root:Packages:NIST:SAS:gDetDist
896
897        detDist = sdd
898        sampleToDetectorDist()          //changes the SDD and wave (DetDist is the global)
899        return(0)
900end
901
902
903// change the huber/sample chamber position
904// if str="checkHuber", the pos is set to huber
905// otherwise, it's set to the sample chamber position
906//
907Function Sim_SetTablePos(strVal)
908        String strVal
909       
910        TableCheckProc(strVal,0)
911        return(0)
912end
913
914
915//
916// ------------ specifc for the simulation
917//
918//
919// type = 0 = none, kills panel
920// type = 1 = 1D
921// type = 2 = 2D
922//
923// -- DON'T leave the do2D global set to 1 - be sure to set back to 0 before exiting
924//
925Function Sim_SetSimulationType(type)
926        Variable type
927
928        STRUCT WMCheckboxAction CB_Struct
929        Variable checkState
930        if(type == 1 || type == 2)
931                CB_Struct.checked = 1
932                checkState = 1
933        else
934                CB_Struct.checked = 0
935                checkState = 0
936        endif
937
938        NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo
939
940        if(type == 1)
941                do2D = 0
942        else
943                do2D = 1
944        endif
945       
946        SimCheckProc(CB_Struct)
947
948        // Very important
949        do2D = 0
950       
951        CheckBox checkSim win=SASCALC,value=checkState
952        return(0)
953End
954
955Function Sim_DoMCSimulation()
956
957        STRUCT WMButtonAction ba
958        ba.eventCode = 2                        //fake mouse click on button
959
960        MC_DoItButtonProc(ba)
961        return(0)
962End
963
964Function Sim_Do1DSimulation()
965
966        ReCalculateInten(1)
967        NVAR estCR = root:Packages:NIST:SAS:g_1DEstDetCR
968        return(estCR)
969End
970
971// counting time (set this last - after all of the instrument moves are done AND the
972// intensity has been updated)
973//
974// ctTime is in seconds
975// sets the global variable, and fakes an "entry" in that field
976Function Sim_SetCountTime(ctTime)
977        Variable ctTime
978       
979        STRUCT WMSetVariableAction sva
980        NVAR ct = root:Packages:NIST:SAS:gCntTime
981        ct=ctTime
982        sva.dval = ctTime                       //seconds
983        sva.eventCode = 3               //live update code
984        CountTimeSetVarProc(sva)
985        return(0)
986End
987
988
989Function Sim_SetIncohXS(val)
990        Variable val
991       
992        NVAR xs = root:Packages:NIST:SAS:gSig_incoh
993        xs = val
994        return(0)
995End
996
997Function Sim_SetThickness(val)
998        Variable val
999       
1000        NVAR th = root:Packages:NIST:SAS:gThick
1001        th = val
1002        return(0)
1003End
1004
1005Function Sim_SetSampleRadius(val)
1006        Variable val
1007       
1008        NVAR r2 = root:Packages:NIST:SAS:gR2
1009        r2 = val
1010       
1011        return(0)
1012End
1013
1014Function Sim_SetNumberOfNeutrons(val)
1015        Variable val
1016       
1017        NVAR num = root:Packages:NIST:SAS:gImon
1018        num = val
1019       
1020        return(0)
1021End
1022
1023// 1 == beam stop in
1024// 0 == beam stop out == transmission
1025//
1026Function Sim_SetBeamStopInOut(val)
1027        Variable val
1028       
1029        NVAR BS = root:Packages:NIST:SAS:gBeamStopIn                    //set to zero for beamstop out (transmission)
1030        WAVE rw=root:Packages:NIST:SAS:realsRead
1031       
1032        BS = val                        //set the global
1033       
1034        // fake the header
1035        if(BS == 1)
1036                rw[37] = 0                              // make sure BS X = 0 if BS is in
1037        else
1038                rw[37] = -10                    // fake BS out as X = -10 cm
1039        endif
1040       
1041        return(0)
1042End
1043
1044Function Sim_SetRawCountsCheck(val)
1045        Variable val
1046       
1047        NVAR rawCt = root:Packages:NIST:SAS:gRawCounts
1048        rawCt = val
1049        return(0)
1050End
1051
1052Function Sim_Set1D_ABS(val)
1053        Variable val
1054       
1055        NVAR gAbs = root:Packages:NIST:SAS:g_1D_DoABS
1056        gAbs = val
1057       
1058        return(0)
1059End
1060
1061Function Sim_Set1D_Noise(val)
1062        Variable val
1063       
1064        NVAR doNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1065        doNoise = val
1066        return(0)
1067End
1068
1069Function Sim_Set1D_Transmission(val)
1070        Variable val
1071       
1072        NVAR trans = root:Packages:NIST:SAS:gSamTrans
1073        trans = val
1074        return(0)
1075End
1076
1077Function Sim_SetModelFunction(strVal)
1078        String strVal
1079       
1080        DoWindow MC_SASCALC
1081        if(V_flag==1)
1082                //then 2D
1083                PopupMenu MC_popup0 win=MC_SASCALC,popmatch=strVal
1084        endif
1085       
1086        DoWindow        Sim_1D_Panel
1087        if(V_flag==1)
1088                //then 1D
1089                PopupMenu MC_popup0 win=Sim_1D_Panel,popmatch=strVal
1090        endif
1091       
1092        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1093        gStr = strVal
1094
1095        return(0)
1096End
1097
1098
1099// ------------
1100
1101//
1102// set the threshold for the time warning dialog to some really large value (say 36000 s)
1103// while doing the simulation to keep the dialog from popping up and stopping the simulation.
1104// then be sure to set it back to 10 s for "normal" simulations that are atttended
1105//
1106// input val is in seconds
1107Function Sim_SetSimTimeWarning(val)
1108        Variable val
1109
1110        NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
1111        SimTimeWarn = val                       //sets the threshold for the warning dialog     
1112
1113        return(0)
1114end
1115
1116// if you need to eliminate the high count data to save the rest of the set
1117Function Sim_TrimOverflowCounts()
1118        Wave w = root:Packages:NIST:SAS:linear_data
1119        w = (w > 2767000) ? 0 : w
1120        return(0)
1121End
1122
1123// if home path exists, save there, otherwise present a dialog
1124// (no sense to use catPathName, since this is simulation, not real data
1125//
1126Function Sim_Save1D_wName(type,titleStr,fname)
1127        String type,titleStr,fname
1128       
1129        String fullPath
1130        NVAR autoSaveIndex = root:Packages:NIST:SAS:gAutoSaveIndex
1131
1132        //now save the data     
1133        PathInfo home
1134        fullPath = S_path + fname
1135       
1136        Save_1DSimData("",runIndex=autoSaveIndex,simLabel=titleStr,saveName=fullPath)
1137       
1138        autoSaveIndex += 1
1139        return(0)
1140End
1141
1142Proc Sim_StoreConfProc(waveStr)
1143        String waveStr
1144       
1145        Sim_StoreConfiguration(waveStr)
1146End
1147//
1148// just make a wave and fill it
1149// I'll keep track of what's in each element
1150// -- not very user friendly, but nobody needs to see this
1151//
1152// poll the panel to fill it in
1153//
1154// 0 = instrument
1155// 1 = number of guides
1156// 2 = source aperture
1157// 3 = lambda
1158// 4 = delta lambda
1159// 5 = huber/chamber
1160// 6 = sample aperture
1161// 7 = lenses/prism/none
1162// 8 = SDD
1163// 9 = offset
1164//
1165Function Sim_StoreConfiguration(waveStr)
1166        String waveStr
1167       
1168        String str=""
1169       
1170        SetDataFolder root:Packages:NIST:SAS
1171       
1172        waveStr = CleanupName("Config_"+waveStr,0)
1173
1174        Make/O/T/N=20 $(waveStr)
1175        Wave/T tw=$(waveStr)
1176        tw = ""
1177       
1178        // 0 = instrument
1179        SVAR instStr = root:Packages:NIST:SAS:gInstStr
1180        tw[0] = "check"+instStr
1181       
1182        // 1 = number of guides
1183        NVAR ng = root:Packages:NIST:SAS:gNg
1184        tw[1] = num2str(ng)
1185       
1186        // 2 = source aperture
1187        ControlInfo/W=SASCALC popup0
1188        tw[2] = S_Value
1189       
1190        // 3 = lambda
1191        NVAR lam = root:Packages:NIST:SAS:gLambda
1192        tw[3] = num2str(lam)
1193       
1194        // 4 = delta lambda
1195//      NVAR dlam = root:Packages:NIST:SAS:gDeltaLambda
1196//      tw[4] = num2str(dlam)
1197       
1198        //or
1199        ControlInfo/W=SASCALC popup0_2
1200        tw[4] = S_Value
1201       
1202       
1203        // 5 = huber/chamber
1204        NVAR sampleTable = root:Packages:NIST:SAS:gTable                //2=chamber, 1=table
1205        if(sampleTable==1)
1206                str = "checkHuber"
1207        else
1208                str = "checkChamber"
1209        endif
1210        tw[5] = str
1211       
1212        // 6 = sample aperture
1213        ControlInfo/W=SASCALC popup0_1
1214        tw[6] = S_Value
1215       
1216        // 7 = lenses/prism/none
1217        NVAR lens = root:Packages:NIST:SAS:gUsingLenses         //==0 for no lenses, 1 for lenses(or prisms)
1218        tw[7] = num2str(lens)
1219       
1220        // 8 = SDD
1221        NVAR sdd = root:Packages:NIST:SAS:gDetDist
1222        tw[8] = num2str(sdd)
1223       
1224        // 9 = offset
1225        NVAR offset = root:Packages:NIST:SAS:gOffset
1226        tw[9] = num2str(offset)
1227       
1228        //
1229       
1230        SetDataFolder root:
1231        return(0)
1232End
1233
1234Proc Sim_moveConfProc(waveStr)
1235        String waveStr
1236        Prompt waveStr,"Select Configuration",popup,ListSASCALCConfigs()
1237       
1238        Sim_MoveToConfiguration($("root:Packages:NIST:SAS:"+waveStr))
1239End
1240// restore the configuration given a wave of information
1241//
1242// be sure to recalculate the intensity after all is set
1243// -- some changes recalculate, some do not...
1244//
1245Function Sim_MoveToConfiguration(tw)
1246        WAVE/T tw
1247
1248
1249// 0 = instrument
1250        Sim_SetInstrument(tw[0])
1251         
1252// 1 = number of guides
1253// 2 = source aperture
1254        Sim_SetNumGuides_SrcAp(str2num(tw[1]),tw[2])
1255       
1256// 3 = lambda
1257        Sim_SetLambda(str2num(tw[3]))
1258       
1259// 4 = delta lambda
1260        Sim_SetDeltaLambda(tw[4])
1261
1262// 5 = huber/chamber
1263        Sim_SetTablePos(tw[5])
1264
1265// 6 = sample aperture
1266        Sim_SetSampleApertureDiam(tw[6])
1267
1268// 7 = lenses/prism/none
1269        Sim_SetLenses(str2num(tw[7]))
1270       
1271// 8 = SDD
1272        Sim_SetSDD(str2num(tw[8]))
1273
1274// 9 = offset
1275        Sim_SetOffset(str2num(tw[9]))
1276
1277
1278        ReCalculateInten(1)
1279        return(0)
1280end
1281
1282
Note: See TracBrowser for help on using the repository browser.