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

Last change on this file since 942 was 942, checked in by srkline, 9 years ago

Some cleanup of the FFT routines to be more exact in declaring mat

Rearranged the Event mode panel so that it's a little more obvious what to do, and in what order

Cleaned up the examples for the Simulation Scripting.

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