Ignore:
Timestamp:
Nov 3, 2008 5:32:02 PM (14 years ago)
Author:
srkline
Message:

Change to MC calculation to avoid unnecessary calculations for transmitted neutrons. Marginal speedup.

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

    r434 r436  
    1212// - Most importantly, this needs to be checked for correctness of the MC simulation 
    1313// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation 
    14 // - why does my integrated tau not match up with John's analytical calculations? where are the assumptions? 
     14// X why does my integrated tau not match up with John's analytical calculations? where are the assumptions? 
    1515// - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles 
    16 // - what is magical about Qu? Is this an assumpution? 
    17  
    18 // - at the larger angles, is the "flat" detector being properly accounted for - in terms of 
     16 
     17// X at the larger angles, is the "flat" detector being properly accounted for - in terms of 
    1918//   the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector 
    2019//   so that what I see is already "corrected"? 
     
    2726// X if no MC desired, still use the selected model 
    2827// X better display of MC results on panel 
    29 // - settings for "count for X seconds" or "how long to 1E6 cts on detector" (run short sim, then multiply) 
     28// - settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...) 
    3029// - add quartz window scattering to the simulation somehow 
    3130// - do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation 
     
    3332//   or the volume fraction of solvent. 
    3433// 
    35 // - add to the results the fraction of coherently scattered neutrons that are singly scattered, different than 
     34// X add to the results the fraction of coherently scattered neutrons that are singly scattered, different than 
    3635//   the overall fraction of singly scattered, and maybe more important to know. 
    3736// 
    38 // - change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons 
     37// X change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons 
    3938//   aren't counted. Is the # that interact a better number? 
    4039// 
     
    160159        results[6] = retWave[5]/retWave[1]                              //double coherent fraction 
    161160        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction 
    162         results[8] = (retWave[0]-retWave[1])/retWave[0]                 //trasnmitted fraction 
     161        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction 
    163162         
    164163        return(0) 
     
    219218        results[6] = retWave[5]/retWave[1]                              //double coherent fraction 
    220219        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction 
    221         results[8] = (retWave[0]-retWave[1])/retWave[0]                 //trasnmitted fraction 
     220        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction 
    222221         
    223222        return(0) 
     
    339338                Phi = 0.0                       //      Intialize azimuthal angle. 
    340339                N1 += 1                 //      Increment total number neutrons counter. 
    341                 DONE = 0                        //      True when neutron is absorbed or when  scattered out of the sample. 
     340                DONE = 0                        //      True when neutron is scattered out of the sample. 
    342341                INDEX = 0                       //      Set counter for number of scattering events. 
    343342                zz = 0.0                        //      Set entering dimension of sample. 
     
    417416                                        NN[INDEX] += 1 
    418417                                Endif 
    419                                 //IF (VZ > 1.0)         // FIX INVALID ARGUMENT 
    420                                         //VZ = 1.0 - 1.2e-7 
    421                                 //ENDIF 
    422                                 Theta_z = acos(Vz)              // Angle WITH respect to z axis. 
    423                                 testQ = 2*pi*sin(theta_z)/wavelength 
    424418                                 
    425                                 // pick a random phi angle, and see if it lands on the detector 
    426                                 // since the scattering is isotropic, I can safely pick a new, random value 
    427                                 // this would not be true if simulating anisotropic scattering. 
    428                                 testPhi = abs(enoise(1))*2*Pi 
    429                                 // is it on the detector?        
    430                                 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
     419                                if(index != 0)          //the neutron interacted at least once, figure out where it ends up 
     420 
     421                                        Theta_z = acos(Vz)              // Angle WITH respect to z axis. 
     422                                        testQ = 2*pi*sin(theta_z)/wavelength 
     423                                         
     424                                        // pick a random phi angle, and see if it lands on the detector 
     425                                        // since the scattering is isotropic, I can safely pick a new, random value 
     426                                        // this would not be true if simulating anisotropic scattering. 
     427                                        testPhi = abs(enoise(1))*2*Pi 
     428                                        // is it on the detector?        
     429                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
     430                                         
     431                                        if(xPixel != -1 && yPixel != -1) 
     432                                                //if(index==1)  // only the single scattering events 
     433                                                        MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering 
     434                                                //endif 
     435                                                        isOn += 1               // neutron that lands on detector 
     436                                        endif 
     437         
     438                                        If(theta_z < theta_max) 
     439                                                //Choose index for scattering angle array. 
     440                                                //IND = NINT(THETA_z/DTH + 0.4999999) 
     441                                                ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint() 
     442                                                NT[ind] += 1                    //Increment bin for angle. 
     443                                                //Increment angle array for single scattering events. 
     444                                                IF(INDEX == 1) 
     445                                                        j1[ind] += 1 
     446                                                Endif 
     447                                                //Increment angle array for double scattering events. 
     448                                                IF (INDEX == 2) 
     449                                                        j2[ind] += 1 
     450                                                Endif 
     451                                        EndIf 
     452                                         
     453                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 
     454                                        NScatterEvents += index         //total number of scattering events 
     455                                        if(index == 1 && incoherentEvent == 1) 
     456                                                NSingleIncoherent += 1 
     457                                        endif 
     458                                        if(index == 1 && coherentEvent == 1) 
     459                                                NSingleCoherent += 1 
     460                                        endif 
     461                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) 
     462                                                NDoubleCoherent += 1 
     463                                        endif 
     464                                        if(index > 1) 
     465                                                NMultipleScatter += 1 
     466                                        endif 
     467                                        //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel 
     468                                else    // if neutron escaped without interacting 
     469                                        // then it must be a transmitted neutron 
     470                                        // don't need to calculate, just increment the proper counters 
     471                                        MC_linear_data[xCtr][yCtr] += 1 
     472                                        isOn += 1 
     473                                        nt[0] += 1 
     474                                         
     475                                        // (Don't do this - it's a waste of time)  
     476                                        // re-initialize and go back and count this neutron again 
     477                                        // go back to xy position 
     478//                                      xx -= ll*vx 
     479//                                      yy -= ll*vy 
     480//                                      zz -= ll*vz 
     481//                                      RR = sqrt(xx*xx+yy*yy)          //radial position of scattering event. 
     482//                       
     483//                                      Vx = 0.0                        // Initialize incident direction vector. 
     484//                                      Vy = 0.0 
     485//                                      Vz = 1.0 
     486//                                      zz = 0 
     487//                                      theta = 0 
     488//                                      phi = 0 
     489                                endif 
    431490                                 
    432                                 if(xPixel != -1 && yPixel != -1) 
    433                                         //if(index==1)  // only the single scattering events 
    434                                                 MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering 
    435                                         //endif 
    436                                                 isOn += 1               // scattered neutron that lands on detector 
    437                                 endif 
    438  
    439                                 If(theta_z < theta_max) 
    440                                         //Choose index for scattering angle array. 
    441                                         //IND = NINT(THETA_z/DTH + 0.4999999) 
    442                                         ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint() 
    443                                         NT[ind] += 1                    //Increment bin for angle. 
    444                                         //Increment angle array for single scattering events. 
    445                                         IF(INDEX == 1) 
    446                                                 j1[ind] += 1 
    447                                         Endif 
    448                                         //Increment angle array for double scattering events. 
    449                                         IF (INDEX == 2) 
    450                                                 j2[ind] += 1 
    451                                         Endif 
    452                                 EndIf 
    453                                  
    454                                 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 
    455                                 NScatterEvents += index         //total number of scattering events 
    456                                 if(index == 1 && incoherentEvent == 1) 
    457                                         NSingleIncoherent += 1 
    458                                 endif 
    459                                 if(index == 1 && coherentEvent == 1) 
    460                                         NSingleCoherent += 1 
    461                                 endif 
    462                                 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) 
    463                                         NDoubleCoherent += 1 
    464                                 endif 
    465                                 if(index > 1) 
    466                                         NMultipleScatter += 1 
    467                                 endif 
    468491                        ENDIF 
    469492                while (!done) 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf

    r434 r436  
    711711                inputWave[10] = sig_sas 
    712712 
     713                linear_data = 0         //initialize 
     714                 
    713715                Variable t0 = stopMStimer(-2) 
    714716         
     
    722724                Variable trans 
    723725                trans = results[8]                      //(n1-n2)/n1 
     726                if(trans == 0) 
     727                        trans = 1 
     728                endif 
    724729 
    725730                // convert to absolute scale 
Note: See TracChangeset for help on using the changeset viewer.