Changeset 552 for sans/XOP_Dev


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

Correct reporting of results wave accounting for multiple coherent scattering

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.c

    r458 r552  
    6868        long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 
    6969        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 
     70        long NMultipleCoherent,NCoherentEvents; 
    7071 
    7172        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
     
    217218        NMultipleScatter = 0; 
    218219        NScatterEvents = 0; 
     220        NMultipleCoherent = 0; 
     221        NCoherentEvents = 0; 
     222 
    219223        isOn = 0; 
    220224         
     
    281285                                        //sprintf(buf,"neutron scatters coherently\r"); 
    282286                                        //XOPNotice(buf); 
    283                                         coherentEvent = 1; 
     287                                        coherentEvent += 1; 
    284288                                        find_theta = 0;                 //false 
    285289                                        do { 
     
    290294                                                 
    291295                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3(&seed))*delta; 
    292                                                 theta = q0/2/pi*wavelength;             //SAS approximation 
     296                                                theta = q0/2/pi*wavelength;             //SAS approximation. 1% error at theta=30 degrees 
    293297                                                 
    294298                                                find_theta = 1;         //always accept 
     
    305309                                        //sprintf(buf,"neutron scatters incoherent\r"); 
    306310                                        //XOPNotice(buf); 
    307                                         incoherentEvent = 1; 
     311                                        incoherentEvent += 1; 
    308312                                  // phi and theta are random over the entire sphere of scattering 
    309313                                        // !can't just choose random theta and phi, won't be random over sphere solid angle 
     
    390394                                                NMultipleScatter += 1; 
    391395                                        } 
     396                                        if(coherentEvent >= 1 && incoherentEvent == 0) { 
     397                                                NCoherentEvents += 1; 
     398                                        } 
     399                                        if(coherentEvent > 1 && incoherentEvent == 0) { 
     400                                                NMultipleCoherent += 1; 
     401                                        } 
     402 
    392403                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data 
    393404                                                isOn += 1; 
    394405                                                nt[0] += 1; 
    395406                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    396                                                 indices[0] = xCtr_long; 
    397                                                 indices[1] = yCtr_long; 
     407                                                //indices[0] = xCtr_long;               //don't put everything in one pixel 
     408                                                //indices[1] = yCtr_long; 
     409                                                indices[0] = (long)round(xCtr+xx/pixSize); 
     410                                                indices[1] = (long)round(yCtr+yy/pixSize); 
     411                                                // check for valid indices - got an XOP error, probably from here 
     412                                                if(indices[0] > 127) indices[0] = 127; 
     413                                                if(indices[0] < 0) indices[0] = 0; 
     414                                                if(indices[1] > 127) indices[1] = 127; 
     415                                                if(indices[1] < 0) indices[1] = 0; 
     416                                                 
    398417                                                if (result = MDGetNumericWavePointValue(wavH, indices, value)) 
    399418                                                        return result; 
     
    429448        if (result = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    430449                return result; 
    431         value[0] = (double)NDoubleCoherent; 
     450        value[0] = (double)NMultipleCoherent; 
    432451        indices[0] = 5; 
    433452        if (result = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
     
    435454        value[0] = (double)NMultipleScatter; 
    436455        indices[0] = 6; 
     456        if (result = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
     457                return result;   
     458        value[0] = (double)NCoherentEvents; 
     459        indices[0] = 7; 
    437460        if (result = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    438461                return result;   
Note: See TracChangeset for help on using the changeset viewer.