Error in Calculation of absolute Flux (Source_flux) ???

Dr. Georg Artus Georg_Artus at Physik.TU-Muenchen.DE
Fri Jan 22 10:28:03 CET 1999


Hello Kristian,

a few days ago my colleague Ralph Gilles (he was also at the Riso
meeting Nov. 98) started with McStas. Due to 'learning' we detected some
problems in calculating the absolute flux with McStas using the
component Source_flux!

Look at the following simple example: It is just a straight guide with a
monitor at the beginning and at the end:

DEFINE INSTRUMENT TEST()

TRACE

COMPONENT source = Source_flux(
        radius = 0.050,
        dist = 2.0,
        xw = 0.02, yh = 0.1,
        E0 = 81.804,    /* = 1A*/
        dE = 0.818,
        flux = 2.542E+13)
  AT (0,0,0) ABSOLUTE

COMPONENT monitor1 = Monitor(
        xmin = -0.010,
        xmax = 0.010,
        ymin = -0.050,
        ymax = 0.050)
  AT (0,0,2) ABSOLUTE

COMPONENT part1 = Guide2(
        w1 = 0.02,
        h1 = 0.1,
        w2 = 0.02,
        h2 = 0.01,
        l = 10.0,
        R0 = 1.0,
        Qc = 0.0214,
        alpha = 5.617,
        mh = 3,
        mv = 2,
        W = 0.0033)
  AT (0,0,2.0) ABSOLUTE

COMPONENT monitor2 = Monitor(
        xmin = -0.010,
        xmax = 0.010,
        ymin = -0.050,
        ymax = 0.050)
  AT (0,0,12) ABSOLUTE

END

Let's compare the results with the same instrument, but now the guide is
focusing (exit height is now 8cm instead of 10cm):

DEFINE INSTRUMENT TEST()

TRACE

COMPONENT source = Source_flux(
        radius = 0.050,
        dist = 2.0,
        xw = 0.02, yh = 0.1,
        E0 = 81.804,    /* = 1A*/
        dE = 0.818,
        flux = 2.542E+13)
  AT (0,0,0) ABSOLUTE

COMPONENT monitor1 = Monitor(
        xmin = -0.010,
        xmax = 0.010,
        ymin = -0.050,
        ymax = 0.050)
  AT (0,0,2) ABSOLUTE

COMPONENT part1 = Guide2(
        w1 = 0.02,
        h1 = 0.1,
        w2 = 0.02,
        h2 = 0.008,
        l = 10.0,
        R0 = 1.0,
        Qc = 0.0214,
        alpha = 5.617,
        mh = 3,
        mv = 2,
        W = 0.0033)
  AT (0,0,2.0) ABSOLUTE

COMPONENT monitor2 = Monitor(
        xmin = -0.010,
        xmax = 0.010,
        ymin = -0.040,
        ymax = 0.040)
  AT (0,0,12) ABSOLUTE

END


McStas:     Transmission (%)   Flux at Monitor2 (units ?)

straight     3.43                3.42E+08
focusing     2.87                2.87E+08

Transmission = Flux at Monitor2 / Flux at Monitor1  
Flux at Monitor1 = 9.977E+09

Going from a straight to a focusing guide the transmission should
decrease, because the average angle between a neutron trajectory an the
top and bottom face of the guide increases and reflection gets less
probable. But the flux at the exit of the guide should increase due to
focusing! It seems as if McStas calculates the flux proportional to the
transmission.
For comparison look at the results for the same guide system calculated
with Beamline (J. Cook):

Beamline:   Transmission (%)   Flux at Monitor2 (n/cm**2/AA/s)

straight     3.5                 1.75E+09
focusing     2.9                 1.82E+09

The transmissions fit well (as we have verified with several other guide
systems). Small differences can be explained by the different handling
of the mirror reflectivity. Here the absolute flux increases going from
straight to focusing guide. Also the fluxes computed by both programs
don't match. This is why I am not sure, what the units of the fluxes
calculated with McStas are and how they are calculated. Due to an erlier
discussion about this matter on this newsgroup (11/23/98) I expected it
to be n/cm**2/s/AA/ster.

"
> What are the units of I and M2 for any monitor when I use source_flux?
> I assume it is n/cm**2/s/AA?

It is n/cm**2/time/st/AA, where the 'st' means steradian (units of solid
angle), and 'time' is the same unit as the one supplied for the 'flux'
parameter of the component (McStas does not really have any comcept of
the duration of an experiment).
"

To get more insight, I tried a simple test without guide: A small source
with radius 1cm and monitor1 of 10cm*10cm in a distance of 1m. For this
example the flux at the monitor can be calculated by hand:

Fi = Q * o * Aq / Ai

Fi = flux at monitor1 (n/cm**2/s/AA)
Q = brightness of source (n/cm**2/s/AA/ster) (I used a wavelength
dependent Q which was calculated for our beamline at FRMII)
o = solid angle (here 0.01 ster)
Aq = Area of source (here pi*cm**2)
Ai = Area of detector (100 cm**2)

Here the numbers:

Wavelength  (AA)       0.5       1.0       2.0       3.0

Flux calc. by hand     7.5E+08   8.0E+09   2.4E+09   4.2E+08

Flux Beamline          7.5E+08   8.0E+09   2.4E+09   4.2E+08

Flux McStas            3.7E+08   8.0E+09   4.8E+09   1.3E+09

While the fluxes calculated by hand and calulated by Beamline fit well,
from the numbers it seems as if McStas calculates fluxes with units
n/cm**2/s!

If we now add a guide to the system above the flux Ff (n/cm**2/s/AA) and
the end of the guide can be calculated as

Ff = nf / Af  
nf = t * ni 
ni = Q * o * Aq   =>

Ff = t * Q * o * Aq / Af

Af = Area of guide exit (=area of monitor2)
nf = neutrons reaching monitor2
ni = neutrons reaching monitor1
t = transmission of guide

The above formula reproduces the fluxes calculated by Beamline
perfectly. 

It seems as if McStas calculates the absolute flux proportional to the
overall transmission, which is not always correct. Also I'm slightly
confused about the units.

Best wishes,

Georg

-- 
*********************************************

Dr. Georg Artus
Technische Universitaet Muenchen
FRM-II Reaktorstation
D-85747 Garching

Tel: +49 (0)89/289-14675
Fax: +49 (0)89/289-14666 or 12112
E-mail: gartus at ph.tum.de



More information about the mcstas-users mailing list