From: Vasilis Vlachoudis (Vasilis.Vlachoudis@cern.ch)
Date: Tue May 22 2007 - 11:27:14 CEST
Dear Yung-Shun
you have a problem with your input. You use in the GEOBEGIN card
SDUM=COMBNAME which means to activate the free format with names for the
geometry, and then you are using a Fixed format with numbers which
should be defined by SDUM=COMBINAT. The two formats should not be mixed!
Once I have changed the SDUM (and the region END card which is
misaligned) it seems that works.
Regards
Vasilis
Yung-Shun Yeh wrote:
> Dear all,
>
> I am simulating neutron propagation in a water ball and getting the
> properties of neutrons across some boundaries. But I encountered some
> problem when I used the subroutine bxdraw in mgdraw.f to get the
> informations in boundaries. Let me describe the detail, there are two
> homocentric balls made of water. One's radius is 5 cm, the other's is
> 10 cm. And the neutron source is put on the center and is isotropic. I
> can take data in 10cm boundary, but can not in 5cm boundary. I don't
> know where I am wrong ? Could anyone give me instrucions? Thanks very
> much!
> Attached is my input file and mgdraw.f and source.f.
>
> Best regards,
> Yung-Shun
> ------------------------------------------------------------------------
>
> *$ 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)'
> *
> LOGICAL LFIRST
> *
> SAVE LFIRST
> DATA LFIRST / .TRUE. /
> PARAMETER (PI=3.141592653589793238462643383279D+00 )
> *======================================================================*
> * *
> * 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
> ***********************************************************************
> *//* Initial beam position
> *// XFLK(NPFLKA)=0.0
> *// YFLK(NPFLKA)=0.0
> *// ZFLK(NPFLKA)=0.0
> *//
> *//* Random direction
> *// THETA=PI*FLRNDM(XXX)
> *// PHI=2*PI*FLRNDM(XXX)
> *// TXFLK(NPFLKA)=SIN(THETA)*COS(PHI)
> *// TYFLK(NPFLKA)=SIN(THETA)*SIN(PHI)
> *// TZFLK(NPFLKA)=COS(THETA)
> * MOMENTUM
> * PBEAM=0.44458
> PBEAM=100
> * PARTICLE TYPE
> IJBEAM=8.0
> ***********************************************************************
> * |
> * +-------------------------------------------------------------------*
> * 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)
> TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
> * Particle momentum
> PMOFLK (NPFLKA) = PBEAM
> * PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
> * & + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
> * Cosines (tx,ty,tz)
> THETA=ACOS(-1+2*FLRNDM(XXX))
> PHI=2*PI*FLRNDM(XXX)
> TXFLK(NPFLKA)=SIN(THETA)*COS(PHI)
> TYFLK(NPFLKA)=SIN(THETA)*SIN(PHI)
> TZFLK(NPFLKA)=COS(THETA)
> *// TXFLK (NPFLKA) = UBEAM
> *// TYFLK (NPFLKA) = VBEAM
> *// TZFLK (NPFLKA) = WBEAM
> * TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
> * & - TYFLK (NPFLKA)**2 )
> * Polarization cosines:
> TXPOL (NPFLKA) = -TWOTWO
> TYPOL (NPFLKA) = +ZERZER
> TZPOL (NPFLKA) = +ZERZER
> * Particle coordinates
> XFLK (NPFLKA) = XBEAM
> YFLK (NPFLKA) = YBEAM
> ZFLK (NPFLKA) = ZBEAM
> * 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
>
>
> ------------------------------------------------------------------------
>
> *$ CREATE MGDRAW.FOR
> *COPY MGDRAW
> * *
> *=== mgdraw ===========================================================*
> * *
> SUBROUTINE MGDRAW ( ICODE, MREG )
>
> INCLUDE '(DBLPRC)'
> INCLUDE '(DIMPAR)'
> INCLUDE '(IOUNIT)'
> *
> *----------------------------------------------------------------------*
> * *
> * Copyright (C) 1990-2006 by Alfredo Ferrari *
> * All Rights Reserved. *
> * *
> * *
> * MaGnetic field trajectory DRAWing: actually this entry manages *
> * all trajectory dumping for *
> * drawing *
> * *
> * Created on 01 march 1990 by Alfredo Ferrari *
> * INFN - Milan *
> * Last change 05-may-06 by Alfredo Ferrari *
> * INFN - Milan *
> * *
> *----------------------------------------------------------------------*
> *
> INCLUDE '(CASLIM)'
> INCLUDE '(COMPUT)'
> INCLUDE '(SOURCM)'
> INCLUDE '(FHEAVY)'
> INCLUDE '(FLKSTK)'
> INCLUDE '(GENSTK)'
> INCLUDE '(MGDDCM)'
> INCLUDE '(PAPROP)'
> INCLUDE '(QUEMGD)'
> INCLUDE '(SUMCOU)'
> INCLUDE '(TRACKR)'
> *
> DIMENSION DTQUEN ( MXTRCK, MAXQMG )
> *
> CHARACTER*20 FILNAM
> LOGICAL LFCOPE
> SAVE LFCOPE
> DATA LFCOPE / .FALSE. /
> *
> *----------------------------------------------------------------------*
> * *
> * Icode = 1: call from Kaskad *
> * Icode = 2: call from Emfsco *
> * Icode = 3: call from Kasneu *
> * Icode = 4: call from Kashea *
> * Icode = 5: call from Kasoph *
> * *
> *----------------------------------------------------------------------*
> * *
> IF ( .NOT. LFCOPE ) THEN
> LFCOPE = .TRUE.
> IF ( KOMPUT .EQ. 2 ) THEN
> FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
> ELSE
> FILNAM = CFDRAW
> END IF
> OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
> & 'UNFORMATTED' )
> END IF
> WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
> & SNGL (WTRACK)
> WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
> & SNGL (ZTRACK (I)), I = 0, NTRACK ),
> & ( SNGL (DTRACK (I)), I = 1, MTRACK ),
> & SNGL (CTRACK)
> * +-------------------------------------------------------------------*
> * | Quenching is activated
> IF ( LQEMGD ) THEN
> RULLL = ZERZER
> CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
> WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)), I = 1, MTRACK ),
> & JBK = 1, NQEMGD )
> END IF
> * | End of quenching
> * +-------------------------------------------------------------------*
> RETURN
> *
> *======================================================================*
> * *
> * Boundary-(X)crossing DRAWing: *
> * *
> * Icode = 1x: call from Kaskad *
> * 19: boundary crossing *
> * Icode = 2x: call from Emfsco *
> * 29: boundary crossing *
> * Icode = 3x: call from Kasneu *
> * 39: boundary crossing *
> * Icode = 4x: call from Kashea *
> * 49: boundary crossing *
> * Icode = 5x: call from Kasoph *
> * 59: boundary crossing *
> * *
> *======================================================================*
> * *
> ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
>
> IF(MREG.EQ.3.AND.NEWREG.EQ.2) THEN
> IF(JTRACK.EQ.8) THEN
> IF(ETRACK.GT.AM(JTRACK)) THEN
> WRITE(42.0,*) ETRACK-AM(JTRACK),
> + XSCO,YSCO,ZSCO
> c WRITE(41.0,*) ETRACK-AM(JTRACK)
> ENDIF
> ENDIF
> ENDIF
>
> IF(MREG.EQ.4.AND.NEWREG.EQ.3) THEN
> IF(JTRACK.EQ.8) THEN
> IF(ETRACK.GT.AM(JTRACK)) THEN
> WRITE(41.0,*) ETRACK-AM(JTRACK),
> + XSCO,YSCO,ZSCO
> c WRITE(41.0,*) ETRACK-AM(JTRACK)
> ENDIF
> ENDIF
> ENDIF
>
> RETURN
> *
> *======================================================================*
> * *
> * Event End DRAWing: *
> * *
> *======================================================================*
> * *
> ENTRY EEDRAW ( ICODE )
> RETURN
> *
> *======================================================================*
> * *
> * ENergy deposition DRAWing: *
> * *
> * Icode = 1x: call from Kaskad *
> * 10: elastic interaction recoil *
> * 11: inelastic interaction recoil *
> * 12: stopping particle *
> * 13: pseudo-neutron deposition *
> * 14: escape *
> * 15: time kill *
> * Icode = 2x: call from Emfsco *
> * 20: local energy deposition (i.e. photoelectric) *
> * 21: below threshold, iarg=1 *
> * 22: below threshold, iarg=2 *
> * 23: escape *
> * 24: time kill *
> * Icode = 3x: call from Kasneu *
> * 30: target recoil *
> * 31: below threshold *
> * 32: escape *
> * 33: time kill *
> * Icode = 4x: call from Kashea *
> * 40: escape *
> * 41: time kill *
> * 42: delta ray stack overflow *
> * Icode = 5x: call from Kasoph *
> * 50: optical photon absorption *
> * 51: escape *
> * 52: time kill *
> * *
> *======================================================================*
> * *
> ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO )
> IF ( .NOT. LFCOPE ) THEN
> LFCOPE = .TRUE.
> IF ( KOMPUT .EQ. 2 ) THEN
> FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
> ELSE
> FILNAM = CFDRAW
> END IF
> OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
> & 'UNFORMATTED' )
> END IF
> WRITE (IODRAW) 0, ICODE, JTRACK, SNGL (ETRACK), SNGL (WTRACK)
> WRITE (IODRAW) SNGL (XSCO), SNGL (YSCO), SNGL (ZSCO), SNGL (RULL)
> * +-------------------------------------------------------------------*
> * | Quenching is activated : calculate quenching factor
> * | and store quenched energy in DTQUEN(1, jbk)
> IF ( LQEMGD ) THEN
> RULLL = RULL
> CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
> WRITE (IODRAW) ( SNGL (DTQUEN(1, JBK)), JBK = 1, NQEMGD )
> END IF
> * | end quenching
> * +-------------------------------------------------------------------*
> RETURN
> *
> *======================================================================*
> * *
> * SOurce particle DRAWing: *
> * *
> *======================================================================*
> *
> ENTRY SODRAW
> IF ( .NOT. LFCOPE ) THEN
> LFCOPE = .TRUE.
> IF ( KOMPUT .EQ. 2 ) THEN
> FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
> ELSE
> FILNAM = CFDRAW
> END IF
> OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
> & 'UNFORMATTED' )
> END IF
> WRITE (IODRAW) -NCASE, NPFLKA, NSTMAX, SNGL (TKESUM),
> & SNGL (WEIPRI)
> * +-------------------------------------------------------------------*
> * | (Radioactive) isotope: it works only for 1 source particle on
> * | the stack for the time being
> IF ( ILOFLK (NPFLKA) .GE. 100000 .AND. LRADDC (NPFLKA) ) THEN
> IARES = MOD ( ILOFLK (NPFLKA), 100000 ) / 100
> IZRES = MOD ( ILOFLK (NPFLKA), 10000000 ) / 100000
> IISRES = ILOFLK (NPFLKA) / 10000000
> IONID = ILOFLK (NPFLKA)
> WRITE (IODRAW) ( IONID,SNGL(-TKEFLK(I)),
> & SNGL (WTFLK(I)), SNGL (XFLK (I)),
> & SNGL (YFLK (I)), SNGL (ZFLK (I)),
> & SNGL (TXFLK(I)), SNGL (TYFLK(I)),
> & SNGL (TZFLK(I)), I = 1, NPFLKA )
> * |
> * +-------------------------------------------------------------------*
> * | Patch for heavy ions: it works only for 1 source particle on
> * | the stack for the time being
> ELSE IF ( ABS (ILOFLK (NPFLKA)) .GE. 10000 ) THEN
> IONID = ILOFLK (NPFLKA)
> CALL DCDION ( IONID )
> WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-IONID)),
> & SNGL (WTFLK(I)), SNGL (XFLK (I)),
> & SNGL (YFLK (I)), SNGL (ZFLK (I)),
> & SNGL (TXFLK(I)), SNGL (TYFLK(I)),
> & SNGL (TZFLK(I)), I = 1, NPFLKA )
> * |
> * +-------------------------------------------------------------------*
> * | Patch for heavy ions: ???
> ELSE IF ( ILOFLK (NPFLKA) .LT. -6 ) THEN
> WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-ILOFLK(NPFLKA))),
> & SNGL (WTFLK(I)), SNGL (XFLK (I)),
> & SNGL (YFLK (I)), SNGL (ZFLK (I)),
> & SNGL (TXFLK(I)), SNGL (TYFLK(I)),
> & SNGL (TZFLK(I)), I = 1, NPFLKA )
> * |
> * +-------------------------------------------------------------------*
> * |
> ELSE
> WRITE (IODRAW) ( ILOFLK(I), SNGL (TKEFLK(I)+AM(ILOFLK(I))),
> & SNGL (WTFLK(I)), SNGL (XFLK (I)),
> & SNGL (YFLK (I)), SNGL (ZFLK (I)),
> & SNGL (TXFLK(I)), SNGL (TYFLK(I)),
> & SNGL (TZFLK(I)), I = 1, NPFLKA )
> END IF
> * |
> * +-------------------------------------------------------------------*
>
> write(42.0) txflk(i), tyflk(i), tzflk(i)
> RETURN
> *
> *======================================================================*
> * *
> * USer dependent DRAWing: *
> * *
> * Icode = 10x: call from Kaskad *
> * 100: elastic interaction secondaries *
> * 101: inelastic interaction secondaries *
> * 102: particle decay secondaries *
> * 103: delta ray generation secondaries *
> * 104: pair production secondaries *
> * 105: bremsstrahlung secondaries *
> * 110: decay products *
> * Icode = 20x: call from Emfsco *
> * 208: bremsstrahlung secondaries *
> * 210: Moller secondaries *
> * 212: Bhabha secondaries *
> * 214: in-flight annihilation secondaries *
> * 215: annihilation at rest secondaries *
> * 217: pair production secondaries *
> * 219: Compton scattering secondaries *
> * 221: photoelectric secondaries *
> * 225: Rayleigh scattering secondaries *
> * Icode = 30x: call from Kasneu *
> * 300: interaction secondaries *
> * Icode = 40x: call from Kashea *
> * 400: delta ray generation secondaries *
> * For all interactions secondaries are put on GENSTK common (kp=1,np) *
> * but for KASHEA delta ray generation where only the secondary elec- *
> * tron is present and stacked on FLKSTK common for kp=npflka *
> * *
> *======================================================================*
> *
> ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
> IF ( .NOT. LFCOPE ) THEN
> LFCOPE = .TRUE.
> IF ( KOMPUT .EQ. 2 ) THEN
> FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
> ELSE
> FILNAM = CFDRAW
> END IF
> OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
> & 'UNFORMATTED' )
> END IF
> * No output by default:
> * write(41.0,*) mreg, sqrt(xsco*xsco+ysco*ysco+zsco*zsco)
>
> RETURN
> *=== End of subrutine Mgdraw ==========================================*
> END
>
>
This archive was generated by hypermail 2.1.6 : Tue May 22 2007 - 12:56:26 CEST