Re: need help for source.f file

From: Alberto Fasso' (fasso@SLAC.Stanford.EDU)
Date: Mon Jan 02 2006 - 19:51:26 CET

  • Next message: Burge, F \(Frances\): "RE: Abort from GTEMFF"

    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
    

  • Next message: Burge, F \(Frances\): "RE: Abort from GTEMFF"

    This archive was generated by hypermail 2.1.6 : Mon Jan 02 2006 - 20:22:18 CET