Re: about BXDRAW or region problem

From: Vasilis Vlachoudis (Vasilis.Vlachoudis@cern.ch)
Date: Tue May 22 2007 - 11:27:14 CEST

  • Next message: Yung-Shun Yeh: "Re: about BXDRAW or region problem"

    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
    >
    >




  • Next message: Yung-Shun Yeh: "Re: about BXDRAW or region problem"

    This archive was generated by hypermail 2.1.6 : Tue May 22 2007 - 12:56:26 CEST