From: Alfredo Ferrari (alfredo.ferrari@cern.ch)
Date: Thu Mar 15 2007 - 16:43:12 CET
Hi all
you are not forced to bin the spectrum into the 72 groups, however the
code will do it for you as soon as a neutron is generated by the source
routine (that is it the corresponding group will be computed and used).
At least you can space some efforts (results will be identical as you
bin it into the 72 group structure yourself).
Ciao
Alfredo
On Thu, 15 Mar 2007, Sebastien WURTH wrote:
> Hello Sunil,
>
> I did almost the same work with Nisy's basis work on Cf252 (thanks again!).
> I tried to sample an Am-Be source with respect to the ISO 8529 norm.
> I respected the 72 neutron groups and to be honnest did some "internal
> kitchen work" to obtain the spectrum I wished.
> Indeed the neutron groups intervalls are a bit to large to define properly
> what really happens between 2 and 11 MeV, but I found my sampling satisfying.
> I leave your questions to the experts, but I would like to show the
> differences between the spectrum obtained with the proper definition of Am-Be
> in terms of ISO 8259 norm and what I managed to sample for FLUKA. See the pdf
> page attached.
> In terms of doses obtained, it matched with real measurements with less than
> 5% error.
>
> Regards.
> Sebastien.
>
> Sunil C wrote :
>
>> Hi All
>>
>> I have the Am-Be neutron spectrum modeled in the source.f attached here,
>> based on what Nisy had done for Cf-252. I have few questions.
>> 1. Is it necessary to bin the source spectrum strictly in accordance to the
>> group structure? What happens if I use the spectrum as given in IAEA tecdoc
>> 413 (2001), the lethargy converted to fluence and normalized to one, but
>> energy remaining exactly as in the book?
>> Such as
>>
>> En (MeV) Fluence
>> 7.94E-02 -
>> 1.00E-01 -
>> 1.25E-01 --
>> 1.58E-01 -
>> 1.99E-01 -
>> 2.51E-01 -
>> 3.16E-01 -
>> ......
>> ....
>> ....
>> ...
>>
>> Should all the 72 groups be listed in the cumulative index? In this case
>> the numbers are zero below 7e-2 MeV and above 14.9 MeV. The cumulative
>> index will show lots of zero below 0.07 MeV and lots of 1.0 above 14.9MeV.
>> In the source attached I have re-binned the spectrum from the IAEA
>> reference to match the group structure and have retained only the first
>> groups with 0.0 and 1.0 as cumulative index.
>>
>> 2. Does the neutron group 'NEUGRP' needs to be faithfully passed on from
>> the source subroutine to the main program? IS the cross section invoked by
>> this information or is there an internal check?
>> In the source.f here I am starting from group 3 by setting the 'NEUGRP'
>> index to 'NEUGRP+2'
>>
>> And Thank you very much Nisy.
>>
>>
>> Best Regards
>>
>>
>> Sunil C
>>
>>
>> ------------------------------------------------------------------------
>>
>> *$ CREATE SOURCE.FOR
>> *COPY SOURCE
>> *
>> *=== source ===========================================================*
>> *
>> SUBROUTINE SOURCE ( NOMORE )
>>
>> INCLUDE '(DBLPRC)'
>> INCLUDE '(DIMPAR)'
>> INCLUDE '(IOUNIT)'
>> *----------------------------------------------------------------------*
>> * Version for isotropic Cf source.
>> *
>> c DIMENSION CFCUM(0:72), ENEDGE(73), ANG(0:40), CUMANG(0:40)
>> DIMENSION CFCUM(0:40), ENEDGE(41), ANG(0:40), CUMANG(0:40)
>> *----------------------------------------------------------------------*
>> *
>> INCLUDE '(BEAMCM)'
>> INCLUDE '(CASLIM)'
>> INCLUDE '(FHEAVY)'
>> INCLUDE '(FLKSTK)'
>> INCLUDE '(IOIOCM)'
>> INCLUDE '(LTCLCM)'
>> INCLUDE '(PAPROP)'
>> INCLUDE '(SOURCM)'
>> INCLUDE '(SUMCOU)'
>> *
>> LOGICAL LFIRST
>> SAVE LFIRST
>> DATA LFIRST / .TRUE. /
>> *----------------------------------------------------------------------*
>> * Neutron energy group boundaries
>> DATA ENEDGE /
>> c Am-Be neutron spectrum
>> & 1.4918E-02, 1.3499E-02, 1.2214E-02, & 1.1052E-02, 1.0000E-02,
>> 9.0484E-03, 8.1873E-03, 7.4082E-03, & 6.7032E-03, 6.0653E-03,
>> 5.4881E-03, 4.9659E-03, 4.4933E-03, & 4.0657E-03, 3.6788E-03,
>> 3.3287E-03, 3.0119E-03, 2.7253E-03, & 2.4660E-03, 2.2313E-03,
>> 2.0190E-03, 1.8268E-03, 1.6530E-03, & 1.4957E-03, 1.3534E-03,
>> 1.2246E-03, 1.1080E-03, 1.0026E-03, & 9.0718E-04, 8.2085E-04,
>> 7.4274E-04, 6.0810E-04, 4.9787E-04, & 4.0762E-04, 3.3373E-04,
>> 2.7324E-04, 2.2371E-04, 1.8316E-04, & 1.4996E-04, 1.2277E-04,
>> 8.6517E-05/ *
>> * Cumulative Am-Be spectrum in group strcuture
>> DATA CFCUM / 0.D0,
>> &1.871941772E-03,4.460515263E-03,
>> &7.681224425E-03,1.166466537E-02,1.654782161E-02,2.244066112E-02,
>> &2.947992860E-02,3.777423905E-02,4.725383086E-02,5.774037795E-02,
>> &6.890897217E-02,8.058700982E-02,9.281815275E-02,1.055467203E-01,
>> &1.186954005E-01,1.327763145E-01,1.492876327E-01,1.684029556E-01,
>> &1.901719153E-01,2.149087836E-01,2.455909710E-01,2.844440519E-01,
>> &3.313109456E-01,3.847880155E-01,4.455136331E-01,5.112235084E-01,
>> &5.824383965E-01,6.541031298E-01,7.240992794E-01,7.917544886E-01,
>> &8.561530386E-01,9.152071332E-01,9.564046266E-01,9.777466352E-01,
>> &9.901464874E-01,9.956112046E-01,9.972629018E-01,9.985685291E-01,
>> &9.994966257E-01,1.000000000E+00/
>>
>> *...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
>> * Angle cosines
>> DATA ANG/ -1.0D0, -0.95D0, -0.9D0, -0.85D0, -0.8D0, -0.75D0,
>> & -0.7D0, -0.65D0, -0.6D0, -0.55D0, -0.5D0, -0.45D0,
>> & -0.4D0, -0.35D0, -0.3D0, -0.25D0, -0.2D0, -0.15D0,
>> & -0.1D0, -5.D-02, 0.0D0, 5.D-02, 0.1D0, 0.15D0,
>> & 0.2D0, 0.25D0, 0.3D0, 0.35D0, 0.4D0, 0.45D0,
>> & 0.5D0, 0.55D0, 0.6D0, 0.65D0, 0.7D0, 0.75D0,
>> & 0.8D0, 0.85D0, 0.9D0, 0.95D0, 1.0D0 /
>> * Cumulative angular distribution
>> DATA CUMANG / 0.0D0, 4.3787904871D-4, 9.309108266D-4,
>> & 1.5080318037D-3, 2.1835838416D-3, 2.9743548826D-3,
>> & 3.8909576476D-3, 4.8923191479D-3, 5.9862770681D-3,
>> & 7.1813938468D-3, 8.4870236975D-3, 9.9133858209D-3,
>> & 1.1471644393D-2, 1.3097707545D-2, 1.4776914040D-2,
>> & 1.6511000723D-2, 1.8301761200D-2, 2.0151047694D-2,
>> & 2.2060772967D-2, 2.4032912292D-2, 2.6069505497D-2,
>> & 2.9084203728D-2, 3.3546756876D-2, 4.0152519402D-2,
>> & 4.9930800378D-2, 6.4405250301D-2, 8.5831276094D-2,
>> & 0.1175474784D0, 0.1569148680D0, 0.1981130258D0,
>> & 0.2412270912D0, 0.2863461629D0, 0.3335634830D0,
>> & 0.3829766299D0, 0.4346877199D0, 0.4919529577D0,
>> & 0.5559660349D0, 0.6275220820D0, 0.7075099237D0,
>> & 0.8189694928D0, 1.0D0 /
>> *----------------------------------------------------------------------*
>> *...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
>> *----------------------------------------------------------------------*
>> NOMORE = 0
>> * +-------------------------------------------------------------------*
>> * | First call initializations:
>> IF ( LFIRST ) THEN
>> * | *** The following 3 cards are mandatory ***
>> TKESUM = ZERZER
>> LFIRST = .FALSE.
>> LUSSRC = .TRUE.
>> * | *** User initialization ***
>> write(lunout,*)
>> write(lunout,'(a,132a)') ("*",i=1,132)
>> write(lunout,*)
>> write(lunout,*)
>> write(lunout,'(a)') "Californium source. Version corrected"
>> write(lunout,*)
>> write(lunout,*)
>> write(lunout,'(a,132a)') ("*",i=1,132)
>> write(lunout,*)
>> END IF
>> * |
>> * +-------------------------------------------------------------------*
>> * Sample the energy group
>> XI = FLRNDM(DUMMY)
>> DO 500 K = 1, 40
>> IF(XI .LE. CFCUM(K)) THEN
>> NEUGRP = K+2
>> ENERGY = ENEDGE(K) - &
>> (XI-CFCUM(K-1))*(ENEDGE(K)-ENEDGE(K+1))/(CFCUM(K)-CFCUM(K-1))
>> GO TO 600 END IF
>> 500 CONTINUE
>> STOP ' Failed to sample the energy group'
>> 600 CONTINUE
>> kount=kount+1
>> *
>> * Sample the cosine of the polar angle
>> XI = FLRNDM(XI)
>> DO 300 K = 1, 40
>> IF(XI .LE. CUMANG(K)) THEN
>> COSTHE = ANG(K-1) + & (XI -
>> CUMANG(K-1))*(ANG(K)-ANG(K-1))/(CUMANG(K)-CUMANG(K-1))
>> GO TO 400 END IF
>> 300 CONTINUE
>> STOP ' Failed to sample the polar angle cosine'
>> 400 CONTINUE
>> SINTHE = SQRT(ONEONE - COSTHE**2)
>> IF(FLRNDM(COSTHE) .LE. HLFHLF) SINTHE = -SINTHE
>> * Samples the azimuthal angle
>> PHI = TWOPIP * FLRNDM(COSTHE)
>> COSPHI = COS(PHI)
>> * 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
>> * +-------------------------------------------------------------------*
>> * | 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
>> IGROUP (NPFLKA) = NEUGRP
>> * Kinetic energy of the particle (GeV)
>> TKEFLK (NPFLKA) = ENERGY
>> * Particle momentum
>> PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
>> * + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
>> * Cosines (tx,ty,tz) (make sure they are properly normalized)
>> TZFLK (NPFLKA) = COSTHE
>> TYFLK (NPFLKA) = SINTHE * COSPHI
>> TXFLK (NPFLKA) = SQRT(ONEONE-TZFLK(NPFLKA)**2-TYFLK(NPFLKA)**2)
>> IF(FLRNDM(COSPHI) .LE. HLFHLF) TXFLK (NPFLKA) = -TXFLK (NPFLKA)
>> * 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
>> * 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
>>
>>
>
>
-- +----------------------------------------------------------------------------+ | Alfredo Ferrari || Tel.: +41.22.767.6119 | | CERN-AB || Fax.: +41.22.767.7555 | | 1211 Geneva 23 || e-mail: Alfredo.Ferrari@cern.ch | | Switzerland || Alfredo.Ferrari@mi.infn.it | +----------------------------------------------------------------------------+
This archive was generated by hypermail 2.1.6 : Thu Mar 15 2007 - 18:15:03 CET