Ignore:
Timestamp:
Sep 10, 2009 3:45:30 PM (13 years ago)
Author:
srkline
Message:

in SASCALC - changed some of the logic to allow 2D simulation

in MultiScatter_MC:

  • changed results wave to report coherent multiple scattering counters correctly (also changed XOP)
  • added save of 1D sim data

in U_CALC:

  • adjusted display for optimal view on windows
  • 7 angle ranges now (added defaults, propagated where all controls are modified)
  • added the angle axis correctly now as a free axis, with a hook linked to the bottom axis
  • removed XS and trans calculations, as integration of ran_dev isn't well

in USANS_EmptyWaves:

  • double precision waves
File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r544 r551  
    6262// 
    6363// 
     64 
     65// setting the flag to 1 == 2D simulation 
     66// any other value for flag == 1D simulation 
     67// 
     68// must remember to close/reopen the simulation control panel to get the correct panel 
     69// 
     70Function Set_2DMonteCarlo_Flag(value) 
     71        Variable value 
     72         
     73        NVAR flag=root:Packages:NIST:SAS:gDoMonteCarlo 
     74        flag=value 
     75        return(0) 
     76end 
    6477 
    6578// threaded call to the main function, adds up the individual runs, and returns what is to be displayed 
     
    175188        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0) 
    176189        results[4] = retWave[3]/retWave[1]                              //avg# times scattered 
    177         results[5] = retWave[4]/retWave[1]                                              //single coherent fraction 
    178         results[6] = retWave[5]/retWave[1]                              //double coherent fraction 
     190        results[5] = retWave[4]/retWave[7]                                              //single coherent fraction 
     191        results[6] = retWave[5]/retWave[7]                              //multiple coherent fraction 
    179192        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction 
    180193        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction 
     
    239252        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0) 
    240253        results[4] = retWave[3]/retWave[1]                              //avg# times scattered 
    241         results[5] = retWave[4]/retWave[1]                                              //single coherent fraction 
    242         results[6] = retWave[5]/retWave[1]                              //double coherent fraction 
     254        results[5] = retWave[4]/retWave[7]                                              //single coherent fraction 
     255        results[6] = retWave[5]/retWave[7]                              //double coherent fraction 
    243256        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction 
    244257        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction 
     
    288301        Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent 
    289302        Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency 
     303        Variable NMultipleCoherent,NCoherentEvents 
    290304         
    291305        detEfficiency = 1.0             //70% counting efficiency = 0.7 
     
    351365        NMultipleScatter = 0 
    352366        NScatterEvents = 0 
     367        NMultipleCoherent = 0 
     368        NCoherentEvents = 0 
    353369 
    354370//C     INITIALIZE ARRAYS. 
     
    415431//                              IF(ran > (ratio1+ratio2) )              //C             NEUTRON SCATTERED coherently 
    416432                                IF(ran > ratio)         //C             NEUTRON SCATTERED coherently 
    417                                         coherentEvent = 1 
     433                                        coherentEvent += 1 
    418434                                        FIND_THETA = 0                  //false 
    419435                                        DO 
     
    441457                  // phi and theta are random over the entire sphere of scattering 
    442458                  // !can't just choose random theta and phi, won't be random over sphere solid angle 
    443                         incoherentEvent = 1 
     459                        incoherentEvent += 1 
    444460                         
    445461                        ran = abs(enoise(1))            //[0,1] 
     
    514530                                                NMultipleScatter += 1 
    515531                                        endif 
     532                                        if(coherentEvent >= 1 && incoherentEvent == 0) 
     533                                                NCoherentEvents += 1 
     534                                        endif 
     535                                        if(coherentEvent > 1 && incoherentEvent == 0) 
     536                                                NMultipleCoherent += 1 
     537                                        endif 
     538                                         
     539                                         
    516540                                        //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel 
    517541                                         
     
    534558        results[1] = n2 
    535559        results[2] = isOn 
    536         results[3] = NScatterEvents             //sum of # of times that neutrons scattered 
     560        results[3] = NScatterEvents             //sum of # of times that neutrons scattered (coh+incoh) 
    537561        results[4] = NSingleCoherent            //# of events that are single, coherent 
    538         results[5] = NDoubleCoherent 
    539         results[6] = NMultipleScatter           //# of multiple scattering events 
     562        results[5] = NMultipleCoherent  //# of scattered neutrons that are coherently scattered more than once 
     563        results[6] = NMultipleScatter           //# of scattered neutrons that are scattered more than once (coh and/or incoh) 
     564        results[7] = NCoherentEvents            //# of scattered neutrons that are scattered coherently one or more times 
    540565         
    541566//      Print "# absorbed = ",n3 
     
    843868//      RenameWindow #,T_results 
    844869//      SetActiveSubwindow ## 
    845 EndMacro 
     870//EndMacro 
    846871 
    847872// as a stand-alone panel, extra control bar  (right) and subwindow implementations don't work right  
    848873// for various reasons... 
    849874Window MC_SASCALC() : Panel 
     875 
     876        // when opening the panel, set the raw counts check to 1 
     877        root:Packages:NIST:SAS:gRawCounts = 1 
     878         
    850879        PauseUpdate; Silent 1           // building window... 
    851880        NewPanel /W=(92,556,713,818)/K=1 as "SANS Simulator" 
     
    866895        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1 
    867896        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo" 
    868         SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)" 
     897        SetVariable cntVar,pos={185,73},size={90,15},proc=CountTimeSetVarProc,title="time(s)" 
    869898        SetVariable cntVar,format="%d" 
    870899        SetVariable cntVar,limits={1,600,1},value= root:Packages:NIST:SAS:gCntTime 
     
    12641293        Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation" 
    12651294        Button MC_button0,fColor=(3,52428,1) 
     1295        Button MC_button1,pos={17,211},size={150,20},proc=Save_1DSimData,title="Save Simulated Data" 
    12661296        GroupBox group0,pos={15,42},size={280,130},title="Sample Setup" 
    12671297        CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset 
     
    13691399        return 0 
    13701400End 
     1401 
     1402 
     1403// 
     1404// 
     1405// 
     1406Function Save_1DSimData(ctrlName) : ButtonControl 
     1407        String ctrlName 
     1408 
     1409        String type="SAS",fullpath="" 
     1410        Variable dialog=1               //=1 will present dialog for name 
     1411         
     1412        String destStr="" 
     1413        destStr = "root:Packages:NIST:"+type 
     1414         
     1415        Variable refNum 
     1416        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n" 
     1417        String fname,ave="C",hdrStr1="",hdrStr2="" 
     1418        Variable step=1 
     1419         
     1420        If(1) 
     1421                //setup a "fake protocol" wave, sice I have no idea of the current state of the data 
     1422                Make/O/T/N=8 root:myGlobals:Protocols:SIMProtocol 
     1423                Wave/T SIMProtocol = $"root:myGlobals:Protocols:SIMProtocol" 
     1424                String junk="****SIMULATED DATA****" 
     1425                //stick in the fake protocol... 
     1426                NVAR ctTime = root:Packages:NIST:SAS:gCntTime 
     1427                NVAR totalCts = root:Packages:NIST:SAS:g_1DTotCts                       //summed counts (simulated) 
     1428                NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR                // estimated detector count rate 
     1429                NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt 
     1430         
     1431                SIMProtocol[0] = junk 
     1432                SIMProtocol[1] = "\tCounting time (s) = "+num2str(ctTime) 
     1433                SIMProtocol[2] = "\tTotal detector counts = "+num2str(totalCts) 
     1434                SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR) 
     1435                SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat) 
     1436                SIMProtocol[5] = junk 
     1437                SIMProtocol[6] = "" 
     1438                SIMProtocol[7] = "" 
     1439                //set the global 
     1440                String/G root:myGlobals:Protocols:gProtoStr = "SIMProtocol" 
     1441        Endif 
     1442         
     1443         
     1444        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error**** 
     1445        WAVE intw=$(destStr + ":integersRead") 
     1446        WAVE rw=$(destStr + ":realsRead") 
     1447        WAVE/T textw=$(destStr + ":textRead") 
     1448        WAVE qvals =$(destStr + ":qval") 
     1449        WAVE inten=$(destStr + ":aveint") 
     1450        WAVE sig=$(destStr + ":sigave") 
     1451        WAVE qbar = $(destStr + ":QBar") 
     1452        WAVE sigmaq = $(destStr + ":SigmaQ") 
     1453        WAVE fsubs = $(destStr + ":fSubS") 
     1454 
     1455        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr 
     1456        Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr) 
     1457         
     1458        //check each wave 
     1459        If(!(WaveExists(intw))) 
     1460                Abort "intw DNExist Save_1DSimData()" 
     1461        Endif 
     1462        If(!(WaveExists(rw))) 
     1463                Abort "rw DNExist Save_1DSimData()" 
     1464        Endif 
     1465        If(!(WaveExists(textw))) 
     1466                Abort "textw DNExist Save_1DSimData()" 
     1467        Endif 
     1468        If(!(WaveExists(qvals))) 
     1469                Abort "qvals DNExist Save_1DSimData()" 
     1470        Endif 
     1471        If(!(WaveExists(inten))) 
     1472                Abort "inten DNExist Save_1DSimData()" 
     1473        Endif 
     1474        If(!(WaveExists(sig))) 
     1475                Abort "sig DNExist Save_1DSimData()" 
     1476        Endif 
     1477        If(!(WaveExists(qbar))) 
     1478                Abort "qbar DNExist Save_1DSimData()" 
     1479        Endif 
     1480        If(!(WaveExists(sigmaq))) 
     1481                Abort "sigmaq DNExist Save_1DSimData()" 
     1482        Endif 
     1483        If(!(WaveExists(fsubs))) 
     1484                Abort "fsubs DNExist Save_1DSimData()" 
     1485        Endif 
     1486        If(!(WaveExists(proto))) 
     1487                Abort "current protocol wave DNExist Save_1DSimData()" 
     1488        Endif 
     1489 
     1490        //strings can be too long to print-- must trim to 255 chars 
     1491        Variable ii,num=8 
     1492        Make/O/T/N=(num) tempShortProto 
     1493        for(ii=0;ii<num;ii+=1) 
     1494                tempShortProto[ii] = (proto[ii])[0,240] 
     1495        endfor 
     1496         
     1497        if(dialog) 
     1498                PathInfo/S catPathName 
     1499                fullPath = DoSaveFileDialog("Save data as") 
     1500                If(cmpstr(fullPath,"")==0) 
     1501                        //user cancel, don't write out a file 
     1502                        Close/A 
     1503                        Abort "no data file was written" 
     1504                Endif 
     1505                //Print "dialog fullpath = ",fullpath 
     1506        Endif 
     1507         
     1508        NVAR monCt = root:Packages:NIST:SAS:gImon 
     1509        NVAR thick = root:Packages:NIST:SAS:gThick 
     1510        NVAR trans = root:Packages:NIST:SAS:gSamTrans                   //for 1D, default value 
     1511         
     1512 
     1513         
     1514        hdrStr1 = num2str(monCt)+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18]) 
     1515        hdrStr1 += "     "+num2str(trans)+"     "+num2str(thick) + ave +"   "+num2str(step) + "\r\n" 
     1516 
     1517        hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    " 
     1518        hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+"ORNL  " + "\r\n" 
     1519         
     1520        //actually open the file here 
     1521        Open refNum as fullpath 
     1522         
     1523        //write out the standard header information 
     1524        fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n","SIMULATED DATA",(date() +"  "+ time()) 
     1525        fprintf refnum,"LABEL: %s\r\n","SIMULATED DATA" 
     1526        fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n" 
     1527        fprintf refnum,hdrStr1 
     1528        fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n" 
     1529        fprintf refnum,hdrStr2 
     1530//      fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step 
     1531 
     1532        //insert protocol information here 
     1533        //-1 list of sample files 
     1534        //0 - bkg 
     1535        //1 - emp 
     1536        //2 - div 
     1537        //3 - mask 
     1538        //4 - abs params c2-c5 
     1539        //5 - average params 
     1540        fprintf refnum, "SAM: %s\r\n",tempShortProto[0] 
     1541        fprintf refnum, "BGD: %s\r\n",tempShortProto[1] 
     1542        fprintf refnum, "EMP: %s\r\n",tempShortProto[2] 
     1543        fprintf refnum, "DIV: %s\r\n",tempShortProto[3] 
     1544        fprintf refnum, "MASK: %s\r\n",tempShortProto[4] 
     1545        fprintf refnum, "ABS: %s\r\n",tempShortProto[5] 
     1546        fprintf refnum, "Average Choices: %s\r\n",tempShortProto[6] 
     1547         
     1548        //write out the data columns 
     1549        fprintf refnum,"The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n" 
     1550        wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs 
     1551         
     1552        Close refnum 
     1553         
     1554        SetDataFolder root:             //(redundant) 
     1555         
     1556        //write confirmation of write operation to history area 
     1557        Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath) 
     1558        KillWaves/Z tempShortProto 
     1559 
     1560        //clear the stuff that was created for case of saving files 
     1561        If(1) 
     1562                Killwaves/Z root:myGlobals:Protocols:SIMProtocol 
     1563                String/G root:myGlobals:Protocols:gProtoStr = "" 
     1564        Endif 
     1565         
     1566         
     1567        return(0) 
     1568         
     1569End 
Note: See TracChangeset for help on using the changeset viewer.