Ignore:
Timestamp:
Nov 25, 2008 4:35:13 PM (14 years ago)
Author:
srkline
Message:

minor changes to MonteCarlo? procedures to add in a realistic beamstop, plus some tests of including absorption and detector efficiency to the simulation. These have been commented out as is seems that they do little to improve the agreement betwen Mc countrate and real data countrate. MC simulation is 2.7 or 3.7x higher CR on detector (only two data points). Need to verify flux estimates first.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS
Files:
2 edited

Legend:

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

    r455 r456  
    1010// 
    1111// 
     12// - Why am I off by a factor of 2.7 - 3.7 (MC too high) relative to real data? 
     13//   I need to include efficiency (70%?) - do I knock these off be fore the simulation or do I  
     14//    really simulate that some fraction of neutrons on the detector don't actually get counted? 
     15//   Is the flux estimate up-to-date? 
    1216// - Most importantly, this needs to be checked for correctness of the MC simulation 
    1317// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation 
    1418// X why does my integrated tau not match up with John's analytical calculations? where are the assumptions? 
    1519// - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles 
    16  
     20// - my simulated transmission is larger than what is measured, even after correcting for the quartz cell. 
     21//   Why? Do I need to include absorption? Just inherent problems with incoherent cross sections? 
     22// 
    1723// X at the larger angles, is the "flat" detector being properly accounted for - in terms of 
    1824//   the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector 
     
    268274        Variable isOn=0,testQ,testPhi,xPixel,yPixel 
    269275        Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent 
    270         Variable NDoubleCoherent,NMultipleScatter 
     276        Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency 
     277         
     278        detEfficiency = 1.0             //70% counting efficiency = 0.7 
    271279         
    272280        imon = inputWave[0] 
     
    291299        num_bins = 200          //number of 1-D bins (not really used) 
    292300         
    293 // my additions - calculate the randome deviate function as needed 
    294 // and calculate the scattering power from the model function 
     301// my additions - calculate the random deviate function as needed 
     302// and calculate the scattering power from the model function (passed in as a wave) 
    295303// 
    296304        Variable left = leftx(ran_dev) 
     
    301309        zpow = sig_sas*thick                    //since I now calculate the sig_sas from the model 
    302310        SIG_ABS = SIGABS_0 * WAVElength 
     311        sig_abs = 0.0           //cm-1 
    303312        SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh 
    304313//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total 
     
    306315//      results[0] = sig_total          //assign these after everything's done 
    307316//      results[1] = sig_sas 
    308 //      RATIO = SIG_ABS / SIG_TOTAL 
     317//      variable ratio1,ratio2 
     318//      ratio1 = sig_abs/sig_total 
     319//      ratio2 = sig_incoh/sig_total 
     320//      // 0->ratio1 = abs 
     321//      // ratio1 -> ratio2 = incoh 
     322//      // > ratio2 = coherent 
    309323        RATIO = sig_incoh / SIG_TOTAL 
    310324         
     
    373387                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2)) 
    374388                                //NEUTRON INTERACTED. 
     389                                ran = abs(enoise(1))            //[0,1] 
     390                                 
     391//                              if(ran<ratio1) 
     392//                                      //absorption event 
     393//                                      n3 +=1 
     394//                                      done=1 
     395//                              else 
     396 
    375397                                INDEX += 1                      //Increment counter of scattering events. 
    376398                                IF(INDEX == 1) 
    377399                                        N2 += 1                 //Increment # of scat. neutrons 
    378400                                Endif 
    379                                 ran = abs(enoise(1))            //[0,1] 
    380401                                //Split neutron interactions into scattering and absorption events 
    381                                 IF(ran > ratio )                //C             NEUTRON SCATTERED coherently 
     402//                              IF(ran > (ratio1+ratio2) )              //C             NEUTRON SCATTERED coherently 
     403                                IF(ran > ratio)         //C             NEUTRON SCATTERED coherently 
    382404                                        coherentEvent = 1 
    383405                                        FIND_THETA = 0                  //false 
     
    414436                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle. 
    415437                                ENDIF           //(ran > ratio) 
     438//                              endif           // event was absorption 
    416439                        ELSE 
    417440                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere 
    418441                                DONE = 1                //done = true, will exit from loop 
     442                                 
     443//                              countIt = 1 
     444//                              if(abs(enoise(1)) > detEfficiency)              //efficiency of 70% wired @top 
     445//                                      countIt = 0                                     //detector does not register 
     446//                              endif 
     447                                 
    419448                                //Increment #scattering events array 
    420449                                If (index <= N_Index) 
     
    430459                                        // since the scattering is isotropic, I can safely pick a new, random value 
    431460                                        // this would not be true if simulating anisotropic scattering. 
    432                                         testPhi = abs(enoise(1))*2*Pi 
     461                                        //testPhi = abs(enoise(1))*2*Pi 
     462                                        testPhi = FindPhi(Vx,Vy)                //use the exiting phi value as defined by Vx and Vy 
     463                                         
    433464                                        // is it on the detector?        
    434465                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
     
    471502                                        endif 
    472503                                        //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel 
     504                                         
    473505                                else    // if neutron escaped without interacting 
     506                                 
    474507                                        // then it must be a transmitted neutron 
    475508                                        // don't need to calculate, just increment the proper counters 
     
    478511                                        nt[0] += 1 
    479512                                         
    480                                         // (Don't do this - it's a waste of time)  
    481                                         // re-initialize and go back and count this neutron again 
    482                                         // go back to xy position 
    483 //                                      xx -= ll*vx 
    484 //                                      yy -= ll*vy 
    485 //                                      zz -= ll*vz 
    486 //                                      RR = sqrt(xx*xx+yy*yy)          //radial position of scattering event. 
    487 //                       
    488 //                                      Vx = 0.0                        // Initialize incident direction vector. 
    489 //                                      Vy = 0.0 
    490 //                                      Vz = 1.0 
    491 //                                      zz = 0 
    492 //                                      theta = 0 
    493 //                                      phi = 0 
    494                                 endif 
    495                                  
     513                                endif           //if interacted 
    496514                        ENDIF 
    497515                while (!done) 
     
    507525        results[6] = NMultipleScatter           //# of multiple scattering events 
    508526         
    509          
     527//      Print "# absorbed = ",n3 
     528 
    510529//      trans_th = exp(-sig_total*thick) 
    511530//      TRANS_exp = (N1-N2) / N1                        // Transmission 
     
    581600        return(0) 
    582601End 
     602 
     603//phi is defined from +x axis, proceeding CCW around [0,2Pi] 
     604Threadsafe Function FindPhi(vx,vy) 
     605        variable vx,vy 
     606         
     607        variable phi 
     608         
     609        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2 
     610         
     611        // special cases 
     612        if(vx==0 && vy > 0) 
     613                return(pi/2) 
     614        endif 
     615        if(vx==0 && vy < 0) 
     616                return(3*pi/2) 
     617        endif 
     618        if(vx >= 0 && vy == 0) 
     619                return(0) 
     620        endif 
     621        if(vx < 0 && vy == 0) 
     622                return(pi) 
     623        endif 
     624         
     625         
     626         
     627        if(vx > 0 && vy > 0) 
     628                return(phi) 
     629        endif 
     630        if(vx < 0 && vy > 0) 
     631                return(abs(phi) + pi/2) 
     632        endif 
     633        if(vx < 0 && vy < 0) 
     634                return(phi + pi) 
     635        endif 
     636        if( vx > 0 && vy < 0) 
     637                return(abs(phi) + 3*pi/2) 
     638        endif 
     639         
     640        return(phi) 
     641end 
     642 
    583643 
    584644ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf

    r455 r456  
    116116        Variable/G root:Packages:NIST:SAS:gDoMonteCarlo = 0 
    117117        Make/O/D/N=10 root:Packages:NIST:SAS:results = 0 
    118         Make/O/T/N=10 root:Packages:NIST:SAS:results_desc = {"total X-section (1/cm)","SAS X-section (1/cm)","number that scatter","number that reach detector","avg # times scattered","fraction single coherent","fraction double coherent","fraction multiple scattered","fraction transmitted","-"} 
     118        Make/O/T/N=10 root:Packages:NIST:SAS:results_desc = {"total X-section (1/cm)","SAS X-section (1/cm)","number that scatter","number that reach detector","avg # times scattered","fraction single coherent","fraction double coherent","fraction multiple scattered","fraction transmitted","detector counts w/beamstop"} 
    119119         
    120120        //tick labels for SDD slider 
     
    773773                endif 
    774774                 
     775                linear_data = 0         //initialize 
    775776// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...           
    776777                t0 = stopMStimer(-2) 
     
    787788                endif 
    788789 
     790                Print "counts on detector = ",sum(linear_data,-inf,inf) 
     791                 
     792                linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike 
     793                Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf) 
     794 
     795                // or simulate a beamstop 
     796                Variable rad=beamstopDiam()/2           //beamstop radius in cm 
     797                rad /= 0.5                              //convert cm to pixels 
     798                rad += 1                                        // add an extra pixel to each side to account for edge 
     799                Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask 
     800                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask 
     801                tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1 
     802                 
     803                linear_data *= tmp_mask 
     804                Print "counts on detector not behind beamstop = ",sum(linear_data,-inf,inf) 
     805                results[9] = sum(linear_data,-inf,inf) 
     806                 
    789807                // convert to absolute scale 
    790                 Variable kappa,beaminten = beamIntensity() 
     808                Variable kappa          //,beaminten = beamIntensity() 
    791809//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten) 
    792810                kappa = thick*(pixSize/sdd)^2*trans*iMon 
    793811 
    794812                linear_data = linear_data / kappa 
    795                 linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike 
     813                 
    796814                data = linear_data 
    797815                 
     
    799817                S_CircularAverageTo1D("SAS") 
    800818                // multiply either estimate by beamstop shadowing 
    801                 aveint *= fSubS 
     819                 
     820//              aveint *= fSubS 
     821 
    802822                // put the new result into the simulation folder 
    803823                Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")      
Note: See TracChangeset for help on using the changeset viewer.