Ignore:
Timestamp:
Nov 3, 2008 12:48:40 PM (14 years ago)
Author:
srkline
Message:

More tweaking of the MonteCarlo?, trying in vain to thread the XOP. As is, works fine with the XOP and single threading.

File:
1 edited

Legend:

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

    r431 r434  
    2929// - settings for "count for X seconds" or "how long to 1E6 cts on detector" (run short sim, then multiply) 
    3030// - add quartz window scattering to the simulation somehow 
    31 // - do smeared models make any sense?? 
     31// - do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation 
    3232// - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition 
    3333//   or the volume fraction of solvent. 
     
    4242//   effects on the absolute scale can be seen? 
    4343// 
    44 // - why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back? 
     44// X why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back? 
    4545// - can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering 
    46 // - ask John how to verify what is going on 
     46// X ask John how to verify what is going on 
    4747// - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up. 
    4848//   a warning has been added - but the models are better fixed with the limiting value. 
    4949// 
    5050// 
     51 
     52// threaded call to the main function, adds up the individual runs, and returns what is to be displayed 
     53// results is calculated and sent back for display 
     54Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     55        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 
     56 
     57        //initialize ran1 in the XOP by passing a negative integer 
     58        // does nothing in the Igor code 
     59        Duplicate/O results retWave 
     60        //results[0] = -1*(datetime) 
     61 
     62        Variable NNeutron=inputWave[0] 
     63        Variable i,nthreads= ThreadProcessorCount 
     64        if(nthreads>2)          //only support 2 processors until I can figure out how to properly thread the XOP and to loop it 
     65                nthreads=2 
     66        endif 
     67         
     68//      nthreads = 1 
     69         
     70        variable mt= ThreadGroupCreate(nthreads) 
     71         
     72        inputWave[0] = NNeutron/nthreads                //split up the number of neutrons 
     73         
     74        for(i=0;i<nthreads;i+=1) 
     75                Duplicate/O nt $("nt"+num2istr(i))              //new instance for each thread 
     76                Duplicate/O j1 $("j1"+num2istr(i)) 
     77                Duplicate/O j2 $("j2"+num2istr(i)) 
     78                Duplicate/O nn $("nn"+num2istr(i)) 
     79                Duplicate/O linear_data $("linear_data"+num2istr(i)) 
     80                Duplicate/O retWave $("retWave"+num2istr(i)) 
     81                Duplicate/O inputWave $("inputWave"+num2istr(i)) 
     82                Duplicate/O ran_dev $("ran_dev"+num2istr(i)) 
     83                 
     84                // ?? I need explicit wave references? 
     85                if(i==0) 
     86                        WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0 
     87                        retWave0[0] = -1*datetime               //to initialize ran3 
     88                        ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0) 
     89                        Print "started thread 0" 
     90                endif 
     91                if(i==1) 
     92                        WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1 
     93                        //retWave1[0] = -1*datetime             //to initialize ran3 
     94                        ThreadStart mt,i,Monte_SANS_W1(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1) 
     95                        Print "started thread 1" 
     96                endif 
     97//              if(i==2) 
     98//                      WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2 
     99//                      retWave2[0] = -1*datetime               //to initialize ran3 
     100//                      ThreadStart mt,i,Monte_SANS_W(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2) 
     101//              endif 
     102//              if(i==3) 
     103//                      WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3 
     104//                      retWave3[0] = -1*datetime               //to initialize ran3 
     105//                      ThreadStart mt,i,Monte_SANS_W(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3) 
     106//              endif 
     107        endfor 
     108 
     109// wait until done 
     110        do 
     111                variable tgs= ThreadGroupWait(mt,100) 
     112        while( tgs != 0 ) 
     113        variable dummy= ThreadGroupRelease(mt) 
     114        Print "done with all threads" 
     115 
     116        // calculate all of the bits for the results 
     117        if(nthreads == 1) 
     118                nt = nt0                // add up each instance 
     119                j1 = j10 
     120                j2 = j20 
     121                nn = nn0 
     122                linear_data = linear_data0 
     123                retWave = retWave0 
     124        endif 
     125        if(nthreads == 2) 
     126                nt = nt0+nt1            // add up each instance 
     127                j1 = j10+j11 
     128                j2 = j20+j21 
     129                nn = nn0+nn1 
     130                linear_data = linear_data0+linear_data1 
     131                retWave = retWave0+retWave1 
     132        endif 
     133//      if(nthreads == 3) 
     134//              nt = nt0+nt1+nt2                // add up each instance 
     135//              j1 = j10+j11+j12 
     136//              j2 = j20+j21+j22 
     137//              nn = nn0+nn1+nn2 
     138//              linear_data = linear_data0+linear_data1+linear_data2 
     139//              retWave = retWave0+retWave1+retWave2 
     140//      endif 
     141//      if(nthreads == 4) 
     142//              nt = nt0+nt1+nt2+nt3            // add up each instance 
     143//              j1 = j10+j11+j12+j13 
     144//              j2 = j20+j21+j22+j23 
     145//              nn = nn0+nn1+nn2+nn3 
     146//              linear_data = linear_data0+linear_data1+linear_data2+linear_data3 
     147//              retWave = retWave0+retWave1+retWave2+retWave3 
     148//      endif 
     149         
     150        // fill up the results wave 
     151        Variable xc,yc 
     152        xc=inputWave[3] 
     153        yc=inputWave[4] 
     154        results[0] = inputWave[9]+inputWave[10]         //total XS 
     155        results[1] = inputWave[10]                                              //SAS XS 
     156        results[2] = retWave[1]                                                 //number that interact n2 
     157        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0) 
     158        results[4] = retWave[3]/retWave[1]                              //avg# times scattered 
     159        results[5] = retWave[4]/retWave[1]                                              //single coherent fraction 
     160        results[6] = retWave[5]/retWave[1]                              //double coherent fraction 
     161        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction 
     162        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //trasnmitted fraction 
     163         
     164        return(0) 
     165End 
     166 
     167// worker function for threads, does nothing except switch between XOP and Igor versions 
     168ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     169        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 
     170         
     171#if exists("Monte_SANSX") 
     172        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     173#else 
     174        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     175#endif 
     176 
     177        return (0) 
     178End 
     179// worker function for threads, does nothing except switch between XOP and Igor versions 
     180ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     181        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 
     182         
     183#if exists("xxxxMonte_SANSX") 
     184        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     185#else 
     186        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     187#endif 
     188 
     189        return (0) 
     190End 
     191 
     192// NON-threaded call to the main function returns what is to be displayed 
     193// results is calculated and sent back for display 
     194Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     195        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 
     196 
     197        //initialize ran1 in the XOP by passing a negative integer 
     198        // does nothing in the Igor code, enoise is already initialized 
     199        Duplicate/O results retWave 
     200        WAVE retWave 
     201        retWave[0] = -1*abs(trunc(100000*enoise(1))) 
     202         
     203#if exists("Monte_SANSX") 
     204        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave) 
     205#else 
     206        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave) 
     207#endif 
     208 
     209        // fill up the results wave 
     210        Variable xc,yc 
     211        xc=inputWave[3] 
     212        yc=inputWave[4] 
     213        results[0] = inputWave[9]+inputWave[10]         //total XS 
     214        results[1] = inputWave[10]                                              //SAS XS 
     215        results[2] = retWave[1]                                                 //number that interact n2 
     216        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0) 
     217        results[4] = retWave[3]/retWave[1]                              //avg# times scattered 
     218        results[5] = retWave[4]/retWave[1]                                              //single coherent fraction 
     219        results[6] = retWave[5]/retWave[1]                              //double coherent fraction 
     220        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction 
     221        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //trasnmitted fraction 
     222         
     223        return(0) 
     224End 
    51225 
    52226 
     
    72246// 
    73247 
    74 Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results) 
     248ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results) 
    75249        WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results 
    76250 
     
    82256        Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG 
    83257        Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI 
    84         //in John's implementation, he dimensioned the indices of the arrays to begin 
    85         // at 0, making things much easier for me... 
    86         //DIMENSION  NT(0:5000),J1(0:5000),J2(0:5000),NN(0:100) 
    87  
    88258        Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat 
    89259        Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single 
     
    93263        Variable Sig_scat,Sig_abs,Ratio,Sig_total 
    94264        Variable isOn=0,testQ,testPhi,xPixel,yPixel 
     265        Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent 
     266        Variable NDoubleCoherent,NMultipleScatter 
    95267         
    96268        imon = inputWave[0] 
     
    122294         
    123295//c       total SAS cross-section 
    124 // 
    125 // input a test value for the incoherent scattering from water 
    126 // 
    127 //      sig_incoh = 5.6                 //[1/cm] as calculated for H2O, now a parameter 
    128 // 
    129296//      SIG_SAS = zpow/thick 
    130297        zpow = sig_sas*thick                    //since I now calculate the sig_sas from the model 
     
    133300//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total 
    134301//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas 
    135         results[0] = sig_total 
    136         results[1] = sig_sas 
     302//      results[0] = sig_total          //assign these after everything's done 
     303//      results[1] = sig_sas 
    137304//      RATIO = SIG_ABS / SIG_TOTAL 
    138305        RATIO = sig_incoh / SIG_TOTAL 
    139 //!!!! the ratio is not yet properly weighted for the volume fractions of each component!!! 
    140306         
    141307//c       assuming theta = sin(theta)...OK 
     
    149315        N2 = 0 
    150316        N3 = 0 
     317        NSingleIncoherent = 0 
     318        NSingleCoherent = 0 
     319        NDoubleCoherent = 0 
     320        NMultipleScatter = 0 
     321        NScatterEvents = 0 
    151322 
    152323//C     INITIALIZE ARRAYS. 
     
    171342                INDEX = 0                       //      Set counter for number of scattering events. 
    172343                zz = 0.0                        //      Set entering dimension of sample. 
     344                incoherentEvent = 0 
     345                coherentEvent = 0 
    173346                 
    174347                do                                      //      Makes sure position is within circle. 
     
    203376                                //Split neutron interactions into scattering and absorption events 
    204377                                IF(ran > ratio )                //C             NEUTRON SCATTERED coherently 
     378                                        coherentEvent = 1 
    205379                                        FIND_THETA = 0                  //false 
    206380                                        DO 
     
    228402                  // phi and theta are random over the entire sphere of scattering 
    229403                  // !can't just choose random theta and phi, won't be random over sphere solid angle 
     404                        incoherentEvent = 1 
    230405                         
    231406                        ran = abs(enoise(1))            //[0,1] 
     
    255430                                FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
    256431                                 
    257 //                              if(xPixel != xCtr && yPixel != yCtr) 
    258 //                                      Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel 
    259 //                              endif 
    260                                  
    261432                                if(xPixel != -1 && yPixel != -1) 
    262                                         isOn += 1 
    263433                                        //if(index==1)  // only the single scattering events 
    264434                                                MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering 
    265435                                        //endif 
     436                                                isOn += 1               // scattered neutron that lands on detector 
    266437                                endif 
    267438 
     
    281452                                EndIf 
    282453                                 
     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 
    283468                        ENDIF 
    284469                while (!done) 
     
    286471 
    287472//      Print "Monte Carlo Done" 
    288         trans_th = exp(-sig_total*thick) 
     473        results[0] = n1 
     474        results[1] = n2 
     475        results[2] = isOn 
     476        results[3] = NScatterEvents             //sum of # of times that neutrons scattered 
     477        results[4] = NSingleCoherent            //# of events that are single, coherent 
     478        results[5] = NDoubleCoherent 
     479        results[6] = NMultipleScatter           //# of multiple scattering events 
     480         
     481         
     482//      trans_th = exp(-sig_total*thick) 
    289483//      TRANS_exp = (N1-N2) / N1                        // Transmission 
    290484        // dsigma/domega assuming isotropic scattering, with no absorption. 
     
    292486//      Print "total # of neutrons reaching 2D detector",isOn 
    293487//      Print "fraction of incident neutrons reaching detector ",isOn/iMon 
    294         results[2] = isOn 
    295         results[3] = isOn/iMon 
    296          
    297                                  
    298 // OUTPUT of the 1D data, not necessary now since I want the 2D 
    299 //      Make/O/N=(num_bins) qvals,int_ms,sig_ms,int_sing,int_doub 
    300 //      ii=1 
    301 //      Print "binning" 
    302 //      do 
    303 //              //CALCULATE MEAN THETA IN BIN. 
    304 //              THETA_z = (ii-0.5)*DTH                  // Mean scattering angle of bin. 
    305 //              //Solid angle of Ith bin. 
    306 //              D_OMEGA = 2*PI*ABS( COS(ii*DTH) - COS((ii-1)*DTH) ) 
    307 //              //SOLID ANGLE CORRECTION: YIELDING CROSS-SECTION. 
    308 //              dS_dW = NT[ii]/(N1*THICK*Trans_th*D_OMEGA) 
    309 //              SIG = NT[ii]^0.5/(N1*THICK*Trans_th*D_OMEGA) 
    310 //              ds_dw_single = J1[ii]/(N1*THICK*Trans_th*D_OMEGA) 
    311 //              ds_dw_double = J2[ii]/(N1*THICK*Trans_th*D_OMEGA) 
    312 //              //Deviation from isotropic model. 
    313 //              qq = 4*pi*sin(0.5*theta_z)/wavelength 
    314 //              qvals[ii-1] = qq 
    315 //              int_ms[ii-1] = dS_dW 
    316 //              sig_ms[ii-1] = sig 
    317 //              int_sing[ii-1] = ds_dw_single 
    318 //              int_doub[ii-1] = ds_dw_double 
    319 //              //140     WRITE (7,145) qq,dS_dW,SIG,ds_dw_single,ds_dw_double 
    320 //              ii+=1 
    321 //      while(ii<=num_bins) 
    322 //      Print "done binning" 
    323  
    324  
    325 //        Write(7,*) 
    326 //        Write(7,*) '#Times Sc.   #Events    ' 
    327 //      DO 150 I = 1,N_INDEX 
    328 //150    WRITE (7,146) I,NN(I) 
    329 //146   Format (I5,T10,F8.0) 
    330  
    331 ///       Write(7,171) N1 
    332 //        Write(7,172) N2 
    333 //        Write(7,173) N3 
    334 //        Write(7,174) N2-N3 
    335  
    336 //171     Format('Total number of neutrons:         N1= ',E10.5) 
    337 ///172     Format('Number of neutrons that interact: N2= ',E10.5) 
    338 //173     Format('Number of absorption events:      N3= ',E10.5) 
    339 //174     Format('# of neutrons that scatter out:(N2-N3)= ',E10.5) 
    340  
     488         
    341489//      Print "Total number of neutrons = ",N1 
    342490//      Print "Total number of neutrons that interact = ",N2 
    343491//      Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2 
    344         results[4] = N2 
    345         results[5] = sum(j1,-inf,inf)/N2 
    346          
    347         Tabs = (N1-N3)/N1 
    348         Ttot = (N1-N2)/N1 
    349         I1_sumI = NN[0]/(N2-N3) 
     492//      results[2] = N2                                         //number that scatter 
     493//      results[3] = isOn - MC_linear_data[xCtr][yCtr]                  //# scattered reaching detector minus zero angle 
     494 
     495         
     496//      Tabs = (N1-N3)/N1 
     497//      Ttot = (N1-N2)/N1 
     498//      I1_sumI = NN[0]/(N2-N3) 
    350499//      Print "Tabs = ",Tabs 
    351500//      Print "Transmitted neutrons = ",Ttot 
    352         results[6] = Ttot 
     501//      results[8] = Ttot 
    353502//      Print "I1 / all I1 = ", I1_sumI 
    354 //      Print "DONE!" 
    355503 
    356504End 
     
    391539End 
    392540 
    393 Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
     541ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
    394542        Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel 
    395543 
     
    484632 
    485633        // SANS Reduction bits 
    486         tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;" 
     634        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;" 
    487635        list = RemoveFromList(tmp, list  ,";") 
    488636        list = RemoveFromList("Monte_SANS", list) 
     
    507655         
    508656        list = SortList(list) 
     657         
     658        list = "default;"+list 
    509659        return(list) 
    510660End               
     
    527677//calculates new direction (xyz) from an old direction 
    528678//theta and phi don't change 
    529 Function NewDirection(vx,vy,vz,theta,phi) 
     679ThreadSafe Function NewDirection(vx,vy,vz,theta,phi) 
    530680        Variable &vx,&vy,&vz 
    531681        Variable theta,phi 
     
    563713End 
    564714 
    565 Function path_len(aval,sig_tot) 
     715ThreadSafe Function path_len(aval,sig_tot) 
    566716        Variable aval,sig_tot 
    567717         
     
    642792End 
    643793 
    644 /////UNUSED, testing routines that have note been updated to work with SASCALC 
     794/////UNUSED, testing routines that have not been updated to work with SASCALC 
    645795// 
    646796//Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr) 
Note: See TracChangeset for help on using the changeset viewer.