From: Vasilis Vlachoudis (Vasilis.Vlachoudis@cern.ch)
Date: Wed Jun 06 2007 - 09:43:14 CEST
Dear Ioannis,
your source routine is wrong!
1) The random number function yyy=FLRNDM(xxx) needs one argument (xxx)
which should be a dummy floating point variable. You are using it with
the FLRNDM(NPFLKA) where the NPFLKA is the stack pointer of FLUKA. Thus
you are modifying the stack position (which is also an integer) and then
results are unpredictable.
2) when you sample from a histogram. You have to select ONLY ONCE a
random number with the FLRNDM and then check this random number against
the distribution. Your routine on the contrary, on every IF statement is
using 2 new random numbers for checking against your distribution.
Therefore there is a big probability that none of the IF is full filled,
moreover the sampling of your spectrum will be totally nonsene.
Please find attached a corrected version of your routine.
Note that the version is not optimized (I am using the same structure
with ifs like in your program). It would be easier and preferable since
your are sampling from a linear histogram with equal width bins to find
the index of the bin with an operation like IBIN =
NINT((XRNDM-0.121)/0.003) and get the results from an array.
Best Regards
Vasilis
Ioannis Kantemiris wrote:
> Dear FLUKA experts
>
> I have a problem with my source.f . When I use :
> IF (FLRNDM(NPFLKA).LT.0.5) THEN
> XFLK (NPFLKA) = -50.D0
> YFLK (NPFLKA) = -7.D0 + 2.D0*FLRNDM(NPFLKA)
> ZFLK (NPFLKA) = 0.4D0 + 5.2D0*FLRNDM(NPFLKA)
> * Cosines (tx,ty,tz)
> TXFLK (NPFLKA) = 1.D0
> TYFLK (NPFLKA) = 0.D0
> TZFLK (NPFLKA) = 0.D0
> ELSE
> XFLK (NPFLKA) = 50.D0
> YFLK (NPFLKA) = -7.D0 + 2.D0*FLRNDM(NPFLKA)
> ZFLK (NPFLKA) = 0.4D0 + 5.2D0*FLRNDM(NPFLKA)
> * Cosines (tx,ty,tz)
> TXFLK (NPFLKA) = -1.D0
> TYFLK (NPFLKA) = 0.D0
> TZFLK (NPFLKA) = 0.D0
> ENDIF
>
> and my run lasts for ~1 hour everything is ok, but when I try to increase
> the simulation time to ~4 hours the simulation stops and I get in the .out
>
> NEXT SEEDS:2CBC78F7 0 0 0 0 0 181CD
> 3039 0 0
>
> **** Isotope tabulation data start at location 3509121 and end at
> 3510058 (I*4 addr.) ****
> 300000 999700000 3489198
> 3.9337001E-03 1.3805460E+04 524304
> NEXT SEEDS: 77B54BE 1 0 0 0 0 181CD
> 3039 0 0
> 400000 999600000 3390137
> 3.9327250E-03 1.3412480E+04 698572
> NEXT SEEDS:1DDAF90D 1 0 0 0 0 181CD
> 3039 0 0
> 1Invalid region or body in Normal. IR= 3 LAT= 0 NASC= 2
> IRRNOR= 3 IR1(IR)= 3 IR2(IR)= 3
> X,Y,Z: -9.22542135150 -8.71320759940 11.2172917844
> 3 -4 -5 -6 -7 -8 -9
> Abort called from NORML reason UNKNOWN PROBLEM Run stopped!
> STOP UNKNOWN PROBLEM
>
> I run geometry debugger and I don't get any geometry errors.
> When I change the source.f eg
>
> IF (FLRNDM(NPFLKA).LT.0.5) THEN
> XFLK (NPFLKA) = -50.D0
> YFLK (NPFLKA) = -7.D0 + 2.D0
> ZFLK (NPFLKA) = 0.4D0 + 5.2D0
> ...
> the simulation runs without any problem.
>
> I attach the input file and the source.f
>
> Thank you in advance
>
> Best regards
>
> Ioannis
> ------------------------------------------------------------------------
>
> *$ CREATE SOURCE.FOR
> *COPY SOURCE
> *
> *=== source ===========================================================*
> *
> SUBROUTINE SOURCE ( NOMORE )
>
> INCLUDE '(DBLPRC)'
> INCLUDE '(DIMPAR)'
> INCLUDE '(IOUNIT)'
> *
> *----------------------------------------------------------------------*
> * *
> * Copyright (C) 1990-2006 by Alfredo Ferrari & Paola Sala *
> * All Rights Reserved. *
> * *
> * *
> * New source for FLUKA9x-FLUKA200x: *
> * *
> * Created on 07 january 1990 by Alfredo Ferrari & Paola Sala *
> * Infn - Milan *
> * *
> * Last change on 03-mar-06 by Alfredo Ferrari *
> * *
> * This is just an example of a possible user written source routine. *
> * note that the beam card still has some meaning - in the scoring the *
> * maximum momentum used in deciding the binning is taken from the *
> * beam momentum. Other beam card parameters are obsolete. *
> * *
> *----------------------------------------------------------------------*
> *
> INCLUDE '(BEAMCM)'
> INCLUDE '(FHEAVY)'
> INCLUDE '(FLKSTK)'
> INCLUDE '(IOIOCM)'
> INCLUDE '(LTCLCM)'
> INCLUDE '(PAPROP)'
> INCLUDE '(SOURCM)'
> INCLUDE '(SUMCOU)'
> *
> INCLUDE '(CASLIM)'
> *
> LOGICAL LFIRST
> *
> SAVE LFIRST
> DATA LFIRST / .TRUE. /
> *======================================================================*
> * *
> * BASIC VERSION *
> * *
> *======================================================================*
> NOMORE = 0
> * +-------------------------------------------------------------------*
> * | First call initializations:
> IF ( LFIRST ) THEN
> * | *** The following 3 cards are mandatory ***
> TKESUM = ZERZER
> LFIRST = .FALSE.
> LUSSRC = .TRUE.
> * | *** User initialization ***
> END IF
> * |
> * +-------------------------------------------------------------------*
> * Push one source particle to the stack. Note that you could as well
> * push many but this way we reserve a maximum amount of space in the
> * stack for the secondaries to be generated
> * Npflka is the stack counter: of course any time source is called it
> * must be =0
> NPFLKA = NPFLKA + 1
> * Wt is the weight of the particle
> WTFLK (NPFLKA) = ONEONE
> WEIPRI = WEIPRI + WTFLK (NPFLKA)
> * Particle type (1=proton.....). Ijbeam is the type set by the BEAM
> * card
> * +-------------------------------------------------------------------*
> * | (Radioactive) isotope:
> IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
> IARES = IPROA
> IZRES = IPROZ
> IISRES = IPROM
> CALL STISBM ( IARES, IZRES, IISRES )
> IJHION = IPROZ * 1000 + IPROA
> IJHION = IJHION * 100 + KXHEAV
> IONID = IJHION
> CALL DCDION ( IONID )
> CALL SETION ( IONID )
> * |
> * +-------------------------------------------------------------------*
> * | Heavy ion:
> ELSE IF ( IJBEAM .EQ. -2 ) THEN
> IJHION = IPROZ * 1000 + IPROA
> IJHION = IJHION * 100 + KXHEAV
> IONID = IJHION
> CALL DCDION ( IONID )
> CALL SETION ( IONID )
> ILOFLK (NPFLKA) = IJHION
> * | Flag this is prompt radiation
> LRADDC (NPFLKA) = .FALSE.
> * |
> * +-------------------------------------------------------------------*
> * | Normal hadron:
> ELSE
> IONID = IJBEAM
> ILOFLK (NPFLKA) = IJBEAM
> * | Flag this is prompt radiation
> LRADDC (NPFLKA) = .FALSE.
> END IF
> * |
> * +-------------------------------------------------------------------*
> * From this point .....
> * Particle generation (1 for primaries)
> LOFLK (NPFLKA) = 1
> * User dependent flag:
> LOUSE (NPFLKA) = 0
> * User dependent spare variables:
> DO 100 ISPR = 1, MKBMX1
> SPAREK (ISPR,NPFLKA) = ZERZER
> 100 CONTINUE
> * User dependent spare flags:
> DO 200 ISPR = 1, MKBMX2
> ISPARK (ISPR,NPFLKA) = 0
> 200 CONTINUE
> * Save the track number of the stack particle:
> ISPARK (MKBMX2,NPFLKA) = NPFLKA
> NPARMA = NPARMA + 1
> NUMPAR (NPFLKA) = NPARMA
> NEVENT (NPFLKA) = 0
> DFNEAR (NPFLKA) = +ZERZER
> * ... to this point: don't change anything
> * Particle age (s)
> AGESTK (NPFLKA) = +ZERZER
> AKNSHR (NPFLKA) = -TWOTWO
> * Group number for "low" energy neutrons, set to 0 anyway
> IGROUP (NPFLKA) = 0
> * Kinetic energy of the particle (GeV)
> IF(FLRNDM(NPFLKA).LE.0.023) THEN
> TKEFLK (NPFLKA)=0.121
> ELSE IF(FLRNDM(NPFLKA).GT.0.023.AND.FLRNDM(NPFLKA).LE.0.049) THEN
> TKEFLK (NPFLKA)=0.124
> ELSE IF(FLRNDM(NPFLKA).GT.0.049.AND.FLRNDM(NPFLKA).LE.0.075) THEN
> TKEFLK (NPFLKA)=0.127
> ELSE IF(FLRNDM(NPFLKA).GT.0.075.AND.FLRNDM(NPFLKA).LE.0.108) THEN
> TKEFLK (NPFLKA)=0.130
> ELSE IF(FLRNDM(NPFLKA).GT.0.108.AND.FLRNDM(NPFLKA).LE.0.148) THEN
> TKEFLK (NPFLKA)=0.133
> ELSE IF(FLRNDM(NPFLKA).GT.0.148.AND.FLRNDM(NPFLKA).LE.0.197) THEN
> TKEFLK (NPFLKA)=0.136
> ELSE IF(FLRNDM(NPFLKA).GT.0.197.AND.FLRNDM(NPFLKA).LE.0.249) THEN
> TKEFLK (NPFLKA)=0.139
> ELSE IF(FLRNDM(NPFLKA).GT.0.249.AND.FLRNDM(NPFLKA).LE.0.305) THEN
> TKEFLK (NPFLKA)=0.142
> ELSE IF(FLRNDM(NPFLKA).GT.0.305.AND.FLRNDM(NPFLKA).LE.0.367) THEN
> TKEFLK (NPFLKA)=0.145
> ELSE IF(FLRNDM(NPFLKA).GT.0.367.AND.FLRNDM(NPFLKA).LE.0.436) THEN
> TKEFLK (NPFLKA)=0.148
> ELSE IF(FLRNDM(NPFLKA).GT.0.436.AND.FLRNDM(NPFLKA).LE.0.508) THEN
> TKEFLK (NPFLKA)=0.151
> ELSE IF(FLRNDM(NPFLKA).GT.0.508.AND.FLRNDM(NPFLKA).LE.0.587) THEN
> TKEFLK (NPFLKA)=0.154
> ELSE IF(FLRNDM(NPFLKA).GT.0.587.AND.FLRNDM(NPFLKA).LE.0.655) THEN
> TKEFLK (NPFLKA)=0.157
> ELSE IF(FLRNDM(NPFLKA).GT.0.655.AND.FLRNDM(NPFLKA).LE.0.744) THEN
> TKEFLK (NPFLKA)=0.160
> ELSE IF(FLRNDM(NPFLKA).GT.0.744.AND.FLRNDM(NPFLKA).LE.0.830) THEN
> TKEFLK (NPFLKA)=0.163
> ELSE IF(FLRNDM(NPFLKA).GT.0.830.AND.FLRNDM(NPFLKA).LE.0.921) THEN
> TKEFLK (NPFLKA)=0.166
> ELSE IF(FLRNDM(NPFLKA).GT.0.921.AND.FLRNDM(NPFLKA).LE.1.) THEN
> TKEFLK (NPFLKA)=0.169
> ENDIF
> * Particle momentum
> * PMOFLK (NPFLKA) = PBEAM
> PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
> & + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
> * Polarization cosines:
> TXPOL (NPFLKA) = -TWOTWO
> TYPOL (NPFLKA) = +ZERZER
> TZPOL (NPFLKA) = +ZERZER
> * Particle coordinates
> IF (FLRNDM(NPFLKA).LT.0.5) THEN
> XFLK (NPFLKA) = -50.D0
> YFLK (NPFLKA) = -7.D0 + 2.D0*FLRNDM(NPFLKA)
> ZFLK (NPFLKA) = 0.4D0 + 5.2D0*FLRNDM(NPFLKA)
> * Cosines (tx,ty,tz)
> TXFLK (NPFLKA) = 1.D0
> TYFLK (NPFLKA) = 0.D0
> TZFLK (NPFLKA) = 0.D0
> ELSE
> XFLK (NPFLKA) = 50.D0
> YFLK (NPFLKA) = -7.D0 + 2.D0*FLRNDM(NPFLKA)
> ZFLK (NPFLKA) = 0.4D0 + 5.2D0*FLRNDM(NPFLKA)
> * Cosines (tx,ty,tz)
> TXFLK (NPFLKA) = -1.D0
> TYFLK (NPFLKA) = 0.D0
> TZFLK (NPFLKA) = 0.D0
> ENDIF
> * Calculate the total kinetic energy of the primaries: don't change
> IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
> & THEN
> TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
> ELSE IF ( ILOFLK (NPFLKA) .NE. 0 ) THEN
> TKESUM = TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) )
> & * WTFLK (NPFLKA)
> ELSE
> TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
> END IF
> RADDLY (NPFLKA) = ZERZER
> * Here we ask for the region number of the hitting point.
> * NREG (NPFLKA) = ...
> * The following line makes the starting region search much more
> * robust if particles are starting very close to a boundary:
> CALL GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
> CALL GEOREG ( XFLK (NPFLKA), YFLK (NPFLKA), ZFLK (NPFLKA),
> & NRGFLK(NPFLKA), IDISC )
> * Do not change these cards:
> CALL GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
> NLATTC (NPFLKA) = MLATTC
> CMPATH (NPFLKA) = ZERZER
> CALL SOEVSV
> RETURN
> *=== End of subroutine Source =========================================*
> END
>
This archive was generated by hypermail 2.1.6 : Wed Jun 06 2007 - 10:19:17 CEST