From: Alberto Fasso' (fasso@SLAC.Stanford.EDU)
Date: Mon Jan 02 2006 - 19:51:26 CET
As already pointed out by Giuseppe Battistoni, the error
is not due to your source routine, but to the lack of a
command BEAM in the input file defining the maximum energy
used in the present run.
However, I would like to call your attention on the fact that there
is no need for a source routine in order to sample uniformly
a starting position between two values in x and y, and an energy
between two values.
With a BEAM card (I copy from the manual):
WHAT(2) > 0.0 : beam momentum spread in GeV/c. The momentum
distribution is assumed to be rectangular
[OK, this is a uniform MOMENTUM distribution, and not an
ENERGY distribution, but at the energies you have mentioned there is
no real difference]
WHAT(4) > 0.0 : beam width in x-direction in cm, if WHAT(6) > 0.0.
The beam profile is assumed to be rectangular.
WHAT(5) > 0.0 : beam width in y-direction in cm, if WHAT(6) > 0.0.
The beam profile is assumed to be rectangular.
A few more, small remarks:
- to "start the random number generator" you must use the command
RANDOMIZe in the input file. The argument of FLRNDM is a dummy which
is not used anywhere, not an initialization value.
- the values passed to subroutine source by command SOURCE are
called WHASOU(1), WHASOU(2), etc., not wha1, wha2 etc.
- Subroutine source is called with only one argument (NOMORE)
and calling it with 3 arguments can only create trouble. (I am
surprised that your compiler did not complain). The variables
WHASOU(1), WHASOU(2) etc. - see above - are passed via COMMON
SOURCM, that you have correctly INCLUDEd.
Alberto
On Sun, 1 Jan 2006, gin chen wrote:
> Hi,
> I need to generate particles ( say protons) randomly
> in x direction from -500.0 c m to 1000.00 cm and in y
> direction the same. Also, the energy of the particle
> needs to be chosen randomly from 500.0 GeV to 20000.0
> GeV.
> Since *.inp file can choose one value for BEAM and
> BEAMPOS, I thought I should use source.f file for this
> purpose. Reading the fluka manual and getting the
> program from src/ area I create the following source.f
> file.
>
> Then I compiled it using the command
> $FLUPRO/flutil/lfluka -m fluka source.f -o myfluka
>
> To run, I use
> $FLUPRO/flutil/rfluka -e myfluka -N0 -M1 mine &
>
> mine is for mine.inp file.
>
> Also I include the line in mine.inp file.
>
> SOURCE 22.0 33.0 45.0
>
> where 22.0 33.0 45.0 are some numbers to start the
> random number generator FLRNDM(xxx).
>
> Then I get the following error message and program
> stops.
> Could someone tell me how to use the source.f file in
> fluka for the above purpose. I can not find anything
> helpful in the manaual or the discussions.
>
> I have attached the error message and source.f file I
> have used.
>
> Thank you,
>
> Gin
>
> =======================ERROR=========================
>
> 1NUMBER OF BEAM NUMBER OF BEAM
> APPROXIMATE NUMBER AVERAGE TIME USED TIME LEFT
> (RESERVED NUMBER OF STARS
> PARTICLES HANDLED PARTICLES LEFT OF BEAM
> PARTICLES BY A BEAM PARTICLE 250.0 SECONDS
> CREATED
> THAT CAN
> STILL BE FOR PRINTOUT)
> HANDLED
>
> NEXT SEEDS: 0 0 0 0 0
> 0 181CD 3039 0 0
> Stepop, dp/dx=<0!! ij, z, a, mmat, ekin, pthrmx,
> dpdx0 1 1 1 26 1343.05945
> 545.843114 -4.66927858E-05
> Abort called from STEPOP reason dp/dx=<0 Run stopped!
> STOP dp/dx=<0
>
>
> ======================= source.f ( modified for x,y,z
> and energy)
> *$ CREATE SOURCE.FOR
> *COPY SOURCE
> *
> *=== source
> ===========================================================*
> *
> SUBROUTINE SOURCE (real wha1,real wha2,real
> wha3)
>
> INCLUDE '(DBLPRC)'
> INCLUDE '(DIMPAR)'
> INCLUDE '(IOUNIT)'
> *
> *----------------------------------------------------------------------*
> *
> *
> * Copyright (C) 1990-2005 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 17-jun-05 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)'
> *
> LOGICAL LFIRST
> *
> SAVE LFIRST
> DATA LFIRST / .TRUE. /
> REAL ENSAMP
> *======================================================================*
> *
> *
> * 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
> * Lstack 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
> *
> +-------------------------------------------------------------------*
> * | Heavy ion:
> IF ( IJBEAM .EQ. -2 ) THEN
> IJHION = IPROZ * 1000 + IPROA
> IJHION = IJHION * 100 + KXHEAV
> IONID = IJHION
> CALL DCDION ( IONID )
> CALL SETION ( IONID )
> ILOFLK (NPFLKA) = IJHION
> * |
> *
> +-------------------------------------------------------------------*
> * | Normal hadron:
> ELSE
> IONID = IJBEAM
> ILOFLK (NPFLKA) = IJBEAM
> 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)
> c TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM
> (IONID)**2 ) - AM (IONID)
> * Particle momentum
> c PMOFLK (NPFLKA) = PBEAM
> * PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * (
> TKEFLK (NPFLKA)
> * & + TWOTWO * AM
> (ILOFLK(NPFLKA)) ) )
> ENSAMP = 500.0 + 19500.0*FLRNDM(wha1)
> TKEFLK(NPFLKA) = ENSAMP
> PMOFLK (NPFLKA) = SQRT(ENSAMP*(ENSAMP +
> TWOTWO*AM(IJBEAM)))
> * Cosines (tx,ty,tz)
> c TXFLK (NPFLKA) = UBEAM
> c TYFLK (NPFLKA) = VBEAM
> c TZFLK (NPFLKA) = WBEAM
> TXFLK (NPFLKA) = 0.0
> TYFLK (NPFLKA) = 0.0
> TZFLK (NPFLKA) = -1.0
>
> * TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK
> (NPFLKA)**2
> * & - TYFLK (NPFLKA)**2 )
> * Polarization cosines:
> TXPOL (NPFLKA) = -TWOTWO
> TYPOL (NPFLKA) = +ZERZER
> TZPOL (NPFLKA) = +ZERZER
> * Particle coordinates
> c XFLK (NPFLKA) = XBEAM
> c YFLK (NPFLKA) = YBEAM
> c ZFLK (NPFLKA) = ZBEAM
> XFLK (NPFLKA) = -500.0 + 1000.0*FLRNDM(wha2)
> YFLK (NPFLKA) = -500.0 + 1000.0*FLRNDM(wha3)
> ZFLK (NPFLKA) = +650.0
> * 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
> * Flag this is prompt radiation
> LRADDC (NPFLKA) = .FALSE.
> 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
>
>
>
>
> __________________________________________
> Yahoo! DSL – Something to write home about.
> Just $16.99/mo. or less.
> dsl.yahoo.com
>
>
>
-- Alberto Fassò SLAC-RP, MS 48, 2575 Sand Hill Road, Menlo Park CA 94025 Phone: (1 650) 926 4762 Fax: (1 650) 926 3569 fasso@slac.stanford.edu
This archive was generated by hypermail 2.1.6 : Mon Jan 02 2006 - 20:22:18 CET