Re: 245Cm

From: Alberto Fasso' <fasso_at_SLAC.Stanford.EDU>
Date: Mon, 31 Aug 2009 12:03:52 -0700 (PDT)

Dear Francesca,

first, a minor point: you sent a log file truncated because it was too long,
and that was right. Only, the information part of the log file is at
its beginning: you should have kept the first 5 lines and delete
everything else. But it doesn't matter, in this case there was no
information, as it can be seen from the other log file.

Now coming to the serious matter: your problem is certainly due to your
source routine.
The two sources you sent are practically identical, except for
replacing Cm instead of Am. Both contain more than one error,
so that it is pure chance that one of the two works (or better: seems
to work - its results are probably wrong).

It is always difficult to check a program written by somebody else,
but I will try to indicate both the parts which are certainly wrong and
parts which are not written in the best way we recommend for FLUKA.
Some of the latter maybe cannot do any harm, but it is difficult to
say with certainty: I recommend that you change them to make sure.

Certainly wrong:
          AHS1=(A1+A2)/2
          AHS=INT(AH1S)
neither AHS1 nor AH1S appear anywhere else, so one of the two is mistyped.

          W = ((WA - ANSYM(AMEZZI))/(1-WA*(ASYM(AHS))))
AMEZZI, AHS are floating point variables, although their decimal part
has been set = 0: they cannot be used as array indeces (I am surprised that
you didn't get a compilation error)

          ZPEAK = MASS * 95/241 + DELTAZ
Dividing 95 (integer) by 241 (integer) gives zero. The above statement
is equivalent to ZPEAK = DELTAZ

In the whole routine no attention is paid to the variable types:
there can be other errors as the one above.
In particular:
most floating point numbers are written as reals, while in FLUKA they
must all be written in double precision:
not -0.45 but -0.45D0
not 0.6 but 0.6D0
not 241./2 but 241.d0/TWOTWO
etc.
Z = INT(Z) (many statements of this type)
Fortran is very tolerant and accepts statements like the above one.
Actually it means mixing a mathematical operation (taking the integer
part of a floating point number) with a type conversion (from integer
to floating point). It would be more clear (and less likely to lead
to an error) to write
IZ = INT(Z)
Z = DBLE(IZ)
But there is an intrinsic function to do this:
Z = DINT(Z) (notice the D! Z is in double precision)

             ASYM(I) = EXP ( -((I-AMEZZI)**2)/(2*(SIGMA**2)))
             ANSYM(I) = EXP(-((I-A2)**2)/(2*SIGMA2**2)) +
      & EXP(-((I-(AINI-A2))**2)/(2*SIGMA2**2)) + HLFHLF *
      & (EXP(-((I-A1)**2)/(2*SIGMA1**2)) +EXP(-((I-(AINI-A1))
      & **2)/(2*SIGMA1**2)) )
Very confusing: I is integer, AMEZZI, A2, AINI, A1 are double precision.
Probably there is no error, but an error can easily creep in in this
kind of programming, and it is difficult for me to check.

SAMPLE = I+((XI-CUMUL(I))/(CUMUL(I+1)-CUMUL(I)))
this is an integer I plus a double precision number < 1, that is
followed by:
MASS = INT(SAMPLE)
Taking the integer part, is equivalent to have MASS = I!

CALL FLNRRN(RGAUSS)
CALL FLNRRN(RGAUSS)
What is the purpose of the double call? Nothing wrong, but only the second
is relevant (waste of time)

* You want a 4PI isotropi angular distribution
Not an error, but the following 7 statements to sample isotropic cosines
could be replaced by a single one:
CALL RACO ( UBEAM, VBEAM, ZBEAM )
Unfortunately this routine is not reported in the Manual, but I recommend
you to use it because it is optimized in various ways.

  99 X1 = 0.D0
       X2 = 4.D0
       SEGNOX = FLRNDM (SSSS)
       IF (SEGNOX.LT.0.5) XFLK (NPFLKA) = -( X2 - X1 ) * FLRNDM(XXXX)
       IF (SEGNOX.GE.0.5) XFLK (NPFLKA) = +( X2 - X1 ) * FLRNDM(XXXX)
Why not simply
       XFLK (NPFLKA) = -X2 + FLRNDM(XXXX) * TWOTWO * X2 ?
(same for Y in the following)

Kind regards,

Alberto

-- 
Alberto Fasso`
SLAC-RP, MS 48, 2575 Sand Hill Road, Menlo Park CA 94025
Phone: (1 650) 926 4762   Fax: (1 650) 926 3569
fasso_at_slac.stanford.edu
Received on Tue Sep 01 2009 - 09:38:44 CEST

This archive was generated by hypermail 2.2.0 : Tue Sep 01 2009 - 09:38:48 CEST