From: Paola Sala (paola.sala@cern.ch)
Date: Fri Jun 22 2007 - 12:15:21 CEST
Dear Ketil
your input options look fine.
As for scoring:
There is no built-in scoring that can give you
double differential distributions in position and energy.
However, I would like to suggest you to consider two possibilities:
- set-up an "hard" position binning, by segmenting your target or your
dummy region in several concentric ones , and score with bxdraw on each
sub-boundary
- score the particle fluence with a cartesian userbin, and the spectra
with usrbdx. This of course does dot give you the correlation between
energy and position.
Moreover, I have a question: does the position really matter? I mean, if
you have downstream detectors, the only observable will be the angular
distibution, not the exit position on your small target.
For mgdraw : what you coded seems correct. The warning about low energy
neutrons is also correct! You must be careful to "bin" your results in a
way that is consistent with the low energy neutron group structure,
otherwise you can get artifacts. You must also take the particle weight
into account, since with the NEW-DEFA card you have biasing in the
thermal neutron group.
In general, built-in scoring is recommended, it is well tested,
reliable, avoids large intermediate files, post-processing is already
available.
For the use of region names in MGDRAW:
Internally, fluka works with region numbers only. However, there are
routines that can be called to convert from name to number and vice
versa. These are:
GEOR2N ( NREG, REGNAM, IERR )
* Input variable: *
* Nreg = region number *
* *
* Output variables: *
* Regnam = region name *
* Ierr = error code (0 on success, 1 on failure) *
* *
*
GEON2R ( REGNAM, NREG, IERR )
* Input variable: *
* Regnam = region name *
* *
* Output variables: *
* Nreg = region number *
* Ierr = error code (0 on success, 1 on failure) *
* *
*----------------------------------------------------------------------*
For both, regnam is a
character*8 variable.
WARNING : geon2r "GEOmetry: Name 2 Region " can be very slow, expecially
on big geometries, since it loops on all region names and compares with
the input name: it should be used to initialize the needed variables and
avoided runtime.
Greetings
Paola
On Thu, 2007-06-21 at 15:28 -0400, Ketil Røed wrote:
>
> Dear All,
>
> I want to score the type of particle, its energy, its position and its
> direction when crossing a boundary between two regions.
> With USRYIELD and USRBDX cards it is possible to score some of this
> information but I can not find the option to score its position in the
> boundary plane.
>
> My geometry setup is simply a 1 cm thick and rectangular silicon target
> with a dummy region at the end for boundary scoring (see geometry in
> .inp file for more details)
>
> For this example I am using a 0.1x0.1 cm2 proton beam of energy 0.1 GeV.
> Later I plan to also use neutron and Heavy Ion beams of different energies.
>
> The type of particles I am intersted in is of course the primare beam
> particle itself (being in this example protons) and some of the
> secondaries produced in the spallation reactions with the Silicon
> (mainly protons, neutrons, Heavy Ions (Z > 1)). At the beam energy in
> this example I dont expect the Heavy Ion secondaries to be energetic
> enough to induce a second spallation reactions, thus I have not enable
> the DPMJET options. However I have enable transportation of all the
> secondaries. At least this is what I believe I have done.
> Can anyone confirm that my settings in the EVENTYPE and PHYSICS cards is
> correct given the above statements. Are there any other settings I have
> forgotten if I want to transport all secondaries produced?
>
> As mentioned in the beginning I have not found a standard fluka card
> that gives me the position of the particle crossing a boundary. However
> I have looked into the possibility of using the mgdraw.f routine. By
> searching through the discussion list and the manual I have found only a
> few examples of how to edit this file.
>
> For boundary crossing I have edited the BXDRAW entry as given below(see
> also enclosed mgdraw.f file). This produces a file that seems to contain
> the information I want. Can anyone comment on the validity of this method?
>
>
> ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
>
> IF(MREG.EQ.3.AND.NEWREG.EQ.4) THEN ! select boundary of interest
> IF(JTRACK.EQ.1) THEN ! select proton
> IF(ETRACK.GT.AM(JTRACK)) THEN ! proton has survived
> WRITE(41.0,*) ETRACK-AM(JTRACK),
> XSCO,YSCO,ZSCO,CXTRCK,CYRCK,CZTRCK,WTRACK
> ENDIF
> ENDIF
> ENDIF
>
> IF(MREG.EQ.3.AND.NEWREG.EQ.4) THEN
> IF(JTRACK.EQ.8) THEN ! select neutron
> IF(ETRACK.GT.AM(JTRACK)) THEN
> WRITE(42.0,*) ETRACK-AM(JTRACK),
> XSCO,YSCO,ZSCO,CXTRCK,CYRCK,CZTRCK,WTRACK,KTRACK
> ENDIF
> ENDIF
> ENDIF
>
> RETURN
>
> Where:
> ETRACK-AM(JTRACK) is the kinetic energy of the particle
> Cx,y,ztrck are the direction cosines of the current particle
> xco,yco,zco are the coordinates of the particle at boundary crossing
> WTRACK is the particle weight
> Ktrack = if > 0 neutron group of the particle (neutron)
>
>
> I have also noticed the following posting in the discussion list:
> http://www.fluka.org/web_archive/earchive/new-fluka-discuss/0807.html
> Thus it seems like using the mgdraw.f routine with low energy neutrons
> is not a good idea. Can anyone please tell me if it is possible to score
> the position and direction of neutrons below 20 MeV when crossing a
> boundary without using the mgdraw.f routine?
>
> I will also like to point out that I had to change from using names to
> numbers on the regions of interest. Using
> "IF(MREG.EQ.TARG1.AND.NEWREG.EQ.TARG2)" did not work. How can I enable
> the use of region names in the mgdraw.f routine?
>
>
> All comments and help are welcome. Thanks!
>
> Best regards,
>
> Ketil Røed
> PhD student
> University of Bergen
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> plain text document attachment (example.inp)
> TITLE
> Proton beam
> * FLAIR: Use names
> GLOBAL 1.0 1.0
> * ...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
> DEFAULTS NEW-DEFA
> BEAM -0.1 -0.0 0.0 0.1 0.1 1.0PROTON
> BEAMPOS 0.0 0.0 -0.5 0.0 0.0
> EVENTYPE 2.0 EVAP
> PHYSICS 3.0 EVAPORAT
> PHYSICS 1.0 COALESCE
> GEOBEGIN COMBNAME
> 0 0 Cylindrical Target
> * *AAA*IIII_________+_________+_________+_________+_________+_________+
> SPH BLK 0.0 0.0 0.0 10000.0
> * vacuum box
> RPP VOI -1000.0 1000.0 -1000.0 1000.0 -1000.0 1000.0
> * target
> ZCC TARG 0.0 0.0 10.0
> * planes limiting the target
> XYP START 0.0
> XYP END 1.0
> YZP TOP 5.0
> YZP BOTTOM -5.0
> XZP LEFT -5.0
> XZP RIGHT 5.0
> * Dummy region before target (for boundary scoring)
> XYP P1 -0.1
> * Dummy region after target (for boundary scoring)
> XYP P2 1.1
> END
> * Regions
> * _AAA.....AA.....AA.....AA.....AA.....AA.....AA.....AA.....AA.....AA.....
> * Black Hole
> BLKHOLE 5 +BLK -VOI
> * Dummy
> DUMMY1 5 +TARG -P1 +START -LEFT +RIGHT -BOTTOM +TOP
> * Target1
> 3.0 5 +TARG -START +END -LEFT +RIGHT -BOTTOM +TOP
> * Target2
> 4.0 5 +TARG +P2 -END -LEFT +RIGHT -BOTTOM +TOP
> * Vac
> VAC 5 +VOI -( +TARG -P1 +P2 -LEFT +RIGHT -BOTTOM +TOP )
> END
> * AAAAAAAxx_________+_________+_________+_________+_________+_________+
> GEOEND
> * ...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
> ASSIGNMA BLCKHOLE BLKHOLE
> ASSIGNMA VACUUM DUMMY1
> ASSIGNMA SILICON 3.0
> ASSIGNMA SILICON 4.0
> ASSIGNMA VACUUM VAC
> PART-THR -0.0001 PROTON PROTON 0.0
> USRBIN 10.0 ENERGY -21.0 5.0 5.0 1.1Energy
> USRBIN -5.0 -5.0 0.0 200.0 200.0 110.0&
> USERDUMP 100.0 37.0 0.0 1.0
> START 5000.0 0.0
> OPEN 42.0
> bdx_neutron2
> OPEN 41.0
> bdx_proton2
> RANDOMIZ 1.0
> STOP
> plain text document attachment (mgdraw.f)
> *$ 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
> * IF ( MTRACK .GT. 0 ) THEN
> * RULLL = ZERZER
> * CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
> * WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)), I = 1, MTRACK ),
> * & JBK = 1, NQEMGD )
> * END IF
> * 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.4) THEN ! select boundary of interest
> IF(JTRACK.EQ.1) THEN ! select proton
> IF(ETRACK.GT.AM(JTRACK)) THEN ! proton has survived
> WRITE(41.0,*) ETRACK-AM(JTRACK), XSCO,YSCO,ZSCO,CXTRCK,CYRCK,CZTRCK,WTRACK
> ENDIF
> ENDIF
> ENDIF
>
> IF(MREG.EQ.3.AND.NEWREG.EQ.4) THEN
> IF(JTRACK.EQ.8) THEN ! select neutron
> IF(ETRACK.GT.AM(JTRACK)) THEN
> WRITE(42.0,*) ETRACK-AM(JTRACK), XSCO,YSCO,ZSCO,CXTRCK,CYRCK,CZTRCK,WTRACK,KTRACK
> 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
> * |
> * +-------------------------------------------------------------------*
> 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:
>
> RETURN
> *=== End of subrutine Mgdraw ==========================================*
> END
>
This archive was generated by hypermail 2.1.6 : Fri Jun 22 2007 - 12:31:55 CEST