From: Ketil Rĝed (ketil.roed@student.uib.no)
Date: Thu Jun 21 2007 - 21:28:41 CEST
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
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
*$ 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 - 09:32:53 CEST