Changeset 553 for sans/XOP_Dev/MonteCarlo/MonteCarlo2.c

Timestamp:

Sep 10, 2009 3:53:23 PM (14 years ago)

Author:

srkline

Message:

Correct reporting of results wave accounting for multiple coherent scattering
(2nd version with different ran - otherwise identical to MonteCarlo?.c)

File:

• sans/XOP_Dev/MonteCarlo/MonteCarlo2.c (modified) (7 diffs)

sans/XOP_Dev/MonteCarlo/MonteCarlo2.c

@@ -70 +70 @@
         long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
         long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
+        long NMultipleCoherent,NCoherentEvents;
+
 
         // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
@@ -219 +221 @@
         NMultipleScatter = 0;
         NScatterEvents = 0;
+        NMultipleCoherent = 0;
+        NCoherentEvents = 0;
+        
         isOn = 0;
         
@@ -283 +288 @@
                                         //sprintf(buf,"neutron scatters coherently\r");
                                         //XOPNotice(buf);
-                                        coherentEvent = 1;
+                                        coherentEvent += 1;
                                         find_theta = 0;                 //false
                                         do {
@@ -307 +312 @@
                                         //sprintf(buf,"neutron scatters incoherent\r");
                                         //XOPNotice(buf);
-                                        incoherentEvent = 1;
+                                        incoherentEvent += 1;
                                   // phi and theta are random over the entire sphere of scattering
                                         // !can't just choose random theta and phi, won't be random over sphere solid angle
@@ -392 +397 @@
                                                 NMultipleScatter += 1;
                                         }
+                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
+                                                NCoherentEvents += 1;
+                                        }
+                                        if(coherentEvent > 1 && incoherentEvent == 0) {
+                                                NMultipleCoherent += 1;
+                                        }
+
                                 } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
                                                 isOn += 1;
                                                 nt[0] += 1;
                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
-                                                indices[0] = xCtr_long;
-                                                indices[1] = yCtr_long;
+                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
+                                                //indices[1] = yCtr_long;
+                                                indices[0] = (long)round(xCtr+xx/pixSize);
+                                                indices[1] = (long)round(yCtr+yy/pixSize);
+                                                // check for valid indices - got an XOP error, probably from here
+                                                if(indices[0] > 127) indices[0] = 127;
+                                                if(indices[0] < 0) indices[0] = 0;
+                                                if(indices[1] > 127) indices[1] = 127;
+                                                if(indices[1] < 0) indices[1] = 0;
+                                                
                                                 if (result = MDGetNumericWavePointValue(wavH, indices, value))
                                                         return result;
@@ -431 +451 @@
         if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
                 return result;
-        value[0] = (double)NDoubleCoherent;
+        value[0] = (double)NMultipleCoherent;
         indices[0] = 5;
         if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
@@ -439 +459 @@
         if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
                 return result;  
-
+        value[0] = (double)NCoherentEvents;
+        indices[0] = 7;
+        if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
+                return result;
 
 //      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave