*$ CREATE SOURCE.FOR *COPY SOURCE * *=== source ===========================================================* * SUBROUTINE SOURCE ( NOMORE ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2020 by Alfredo Ferrari & Paola Sala * * All Rights Reserved. * * * * * * New source for FLUKA9x-FLUKA20xy: * * * * Created on 07 January 1990 by Alfredo Ferrari & Paola Sala * * Infn - Milan * * * * Last change on 03-Jan-20 by Paola Sala * * Infn - Milan * * * * 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. * * * * Output variables: * * * * Nomore = if > 0 the run will be terminated * * * *----------------------------------------------------------------------* * INCLUDE '(BEAMCM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(IOIOCM)' INCLUDE '(LTCLCM)' INCLUDE '(PAPROP)' INCLUDE '(SOURCM)' INCLUDE '(SUMCOU)' * real*8 ParRndm !My Code********* integer,parameter::neutron_dimension=1012,gamma_dimension=1001 real*8 energy_low, & energy_high,spec_count, & neutron_spec_count(neutron_dimension), & gamma_spec_count(gamma_dimension), & increment_count character(len=1000)::temp,temp1 integer i real*8 neu_spec_pro(1012),gamma_spec_pro(1001) real*8 energy real*8 DirRndm real*8 EneRndm real*8 thi * LOGICAL LFIRST, LISNUT SAVE LFIRST DATA LFIRST / .TRUE. / * Statement function: LISNUT (IJ) = INDEX ( PRNAME (IJ), 'NEUTRI' ) .GT. 0 *======================================================================* * * * 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 * | * +-------------------------------------------------------------------* * 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) neutron_spec_count(1)=0 gamma_spec_count(1)=0 ParRndm=FLRNDM(XDUMMY) !My Code********* IF(ParRndm.LT.1.076E-3/(1.076E-3+6.91E-6))THEN !My Code********* IJBEAM=8 !My Code********* ELSE !My Code********* IJBEAM=7 !My Code********* END IF !My Code********* * 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 = IPROM * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA IJHION = IJHION * 100 + KXHEAV IONID = IJHION CALL DCDION ( IONID ) CALL SETION ( IONID ) EEXSTK (NPFLKA) = EXENRG (IONID) TMNSTK (NPFLKA) = TMNLF (IONID) ILVSTK (NPFLKA) = IEXLVL (IONID) LFRPHN (NPFLKA) = .FALSE. * | * +-------------------------------------------------------------------* * | Heavy ion: ELSE IF ( IJBEAM .EQ. -2 ) THEN IJHION = IPROM * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA IJHION = IJHION * 100 + KXHEAV IONID = IJHION CALL DCDION ( IONID ) CALL SETION ( IONID ) EEXSTK (NPFLKA) = EXENRG (IONID) TMNSTK (NPFLKA) = TMNLF (IONID) ILVSTK (NPFLKA) = IEXLVL (IONID) ILOFLK (NPFLKA) = IJHION * | Flag this is prompt radiation LRADDC (NPFLKA) = .FALSE. * | Group number for "low" energy neutrons, set to 0 anyway IGROUP (NPFLKA) = 0 * | Parent radioactive isotope: IRDAZM (NPFLKA) = 0 * | Particle age (s) AGESTK (NPFLKA) = +ZERZER * | 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 (IONID) ) ) LFRPHN (NPFLKA) = .FALSE. * | * +-------------------------------------------------------------------* * | Normal hadron: ELSE IONID = IJBEAM EEXSTK (NPFLKA) = EXENRG (IJBEAM) TMNSTK (NPFLKA) = TMNLF (IJBEAM) ILVSTK (NPFLKA) = IEXLVL (IJBEAM) ILOFLK (NPFLKA) = IJBEAM * | Flag this is prompt radiation LRADDC (NPFLKA) = .FALSE. * | Group number for "low" energy neutrons, set to 0 anyway IGROUP (NPFLKA) = 0 * | Parent radioactive isotope: IRDAZM (NPFLKA) = 0 * | Particle age (s) AGESTK (NPFLKA) = +ZERZER * | Kinetic energy of the particle (GeV) OPEN(UNIT=11,FILE="/home/llh/flukaworks/np_spectrum/80M/n31") OPEN(UNIT=12,FILE="/home/llh/flukaworks/np_spectrum/80M/n32") OPEN(UNIT=13,FILE="/home/llh/flukaworks/np_spectrum/80M/n33") OPEN(UNIT=14,FILE="/home/llh/flukaworks/np_spectrum/80M/n34") OPEN(UNIT=15,FILE="/home/llh/flukaworks/np_spectrum/80M/n35") OPEN(UNIT=16,FILE="/home/llh/flukaworks/np_spectrum/80M/n36") OPEN(UNIT=17,FILE="/home/llh/flukaworks/np_spectrum/80M/n37") OPEN(UNIT=21,FILE="/home/llh/flukaworks/np_spectrum/80M/p51") OPEN(UNIT=22,FILE="/home/llh/flukaworks/np_spectrum/80M/p52") * Cosines (tx,ty,tz) ****************************************************************** * when the particle is neutron ****************************************************************** IF(IJBEAM.eq.8)THEN DirRndm=FLRNDM(XDUMMY) ***************************************************************** *neutron moment angle is 0-2.5 ***************************************************************** IF(DirRndm.le.3.093E-5/1.076E-3)THEN WBEAM=cos(FLRNDM(XDUMMY)*2.5/180*3.1415926) read(11,"(A1000)")temp1 read(11,"(A1000)")temp1 do i=1,1011 read(11,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(11) read(11,"(A1000)")temp1 read(11,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(11,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************ *the neutron moment angle is 0-2.5 ************************************************************************ ************************************************************************* *the neutron moment angle is 2.5-5 ************************************************************************ else if(DirRndm.gt.3.093E-5/1.076E-3.and. & DirRndm.le.1.146e-4/1.076E-3)then increment_count=0 WBEAM=cos((FLRNDM(XDUMMY)*2.5+2.5)/180*3.1415926) read(12,"(A1000)")temp1 read(12,"(A1000)")temp1 do i=1,1011 read(12,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(12) read(12,"(A1000)")temp1 read(12,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(12,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************* *the neutron moment angle is 2.5-5 ************************************************************************ ************************************************************************* *the neutron moment angle is 5-7.5 ************************************************************************ else IF(DirRndm.gt.1.146e-4/1.076E-3.and. & DirRndm.le.2.645e-4/1.076E-3)THEN WBEAM=cos((FLRNDM(XDUMMY)*2.5+5)/180*3.1415926) increment_count=0 read(13,"(A1000)")temp1 read(13,"(A1000)")temp1 do i=1,1011 read(13,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(13) read(13,"(A1000)")temp1 read(13,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(13,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************* *the neutron moment angle is 5-7.5 ************************************************************************ ************************************************************************* *the neutron moment angle is 7.5-10 ************************************************************************ else IF(DirRndm.gt.2.645e-4/1.076E-3.and. & DirRndm.le.4.17e-4/1.076E-3)THEN WBEAM=cos((FLRNDM(XDUMMY)*2.5+7.5)/180*3.1415926) increment_count=0 read(14,"(A1000)")temp1 read(14,"(A1000)")temp1 do i=1,1011 read(14,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(14) read(14,"(A1000)")temp1 read(14,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(14,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************* *the neutron moment angle is 7.5-10 ************************************************************************ ************************************************************************* *the neutron moment angle is 10-12.5 ************************************************************************ else IF(DirRndm.gt.4.17e-4/1.076E-3.and. & DirRndm.le.6.173e-4/1.076E-3)THEN WBEAM=cos((FLRNDM(XDUMMY)*2.5+10)/180*3.1415926) increment_count=0 read(15,"(A1000)")temp1 read(15,"(A1000)")temp1 do i=1,1011 read(15,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(15) read(15,"(A1000)")temp1 read(15,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(15,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************* *the neutron moment angle is 10-12.5 ************************************************************************ ************************************************************************* *the neutron moment angle is 12.5-15 ************************************************************************ else IF(DirRndm.gt.6.173e-4/1.076E-3.and. & DirRndm.le.8.385e-4/1.076E-3)THEN WBEAM=cos((FLRNDM(XDUMMY)*2.5+12.5)/180*3.1415926) increment_count=0 read(16,"(A1000)")temp1 read(16,"(A1000)")temp1 do i=1,1011 read(16,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(16) read(16,"(A1000)")temp1 read(16,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(16,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************* *the neutron moment angle is 12.5-15 ************************************************************************ ************************************************************************* *the neutron moment angle is 15-17.5 ************************************************************************ else IF(DirRndm.gt.8.385e-4/1.076E-3)THEN WBEAM=cos((FLRNDM(XDUMMY)*2.5+15)/180*3.1415926) increment_count=0 read(17,"(A1000)")temp1 read(17,"(A1000)")temp1 do i=1,1011 read(17,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) neutron_spec_count(i+1)=increment_count end do do i=1,1012 neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012) end do rewind(17) read(17,"(A1000)")temp1 read(17,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1011 read(17,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1)) & then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do end if ************************************************************************* *the neutron moment angle is 15-17.5 ************************************************************************ ************************************************************************ * when the particle is neutron ************************************************************************ ************************************************************************ * when the particle is photon ************************************************************************ ************************************************************************ *the photon moment angle is 0-15 ************************************************************************ else if(ijbeam.eq.7)then DirRndm=FLRNDM(XDUMMY) if(DirRndm.le.1.762/6.906)then wbeam=cos(Flrndm(xdummy)*15/180*3.1415926) increment_count=0 read(21,"(A1000)")temp1 read(21,"(A1000)")temp1 do i=1,1000 read(21,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) gamma_spec_count(i+1)=increment_count end do do i=1,1001 gamma_spec_pro(i)=gamma_spec_count(i)/gamma_spec_count(1001) end do rewind(21) read(21,"(A1000)")temp1 read(21,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1000 read(21,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.gamma_spec_pro(i).and.EneRndm.le. & gamma_spec_pro(i+1))then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do ************************************************************************ *the photon moment angle is 0-15 ************************************************************************ ************************************************************************ *the photon moment angle is 15-30 ************************************************************************ else if(DirRndm.gt.1.762/6.906)then wbeam=cos((Flrndm(xdummy)*15+15)/180*3.1415926) increment_count=0 read(22,"(A1000)")temp1 read(22,"(A1000)")temp1 do i=1,1000 read(22,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high read(temp(25:33),*)spec_count increment_count=increment_count+spec_count*(energy_high & -energy_low) gamma_spec_count(i+1)=increment_count end do do i=1,1001 gamma_spec_pro(i)=gamma_spec_count(i)/gamma_spec_count(1001) end do rewind(22) read(22,"(A1000)")temp1 read(22,"(A1000)")temp1 EneRndm=FLRNDM(XDUMMY) do i=1,1000 read(22,"(A1000)")temp read(temp(3:11),*)energy_low read(temp(14:22),*)energy_high if(EneRndm.gt.gamma_spec_pro(i).and.EneRndm.le. & gamma_spec_pro(i+1))then energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY) end if end do end if end if ************************************************************************ *the photon moment angle is 15-30 ************************************************************************ ************************************************************************ * when the particle is photon ************************************************************************ PBEAM=SQRT(energy*(energy+2*am(ionid))) * PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA) * & + TWOTWO * AM (IONID) ) ) thi=flrndm(xdummy)*360 ubeam=sin(thi/180.*3.1415926)*sqrt(1-wbeam*wbeam) vbeam=cos(thi/180.*3.1415926)*sqrt(1-wbeam*wbeam) TXFLK (NPFLKA) = UBEAM TYFLK (NPFLKA) = VBEAM TZFLK (NPFLKA) = WBEAM * TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2 * & - TYFLK (NPFLKA)**2 ) TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) & - AM (IONID) * | Particle momentum PMOFLK (NPFLKA) = PBEAM * PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA) * & + TWOTWO * AM (IONID) ) ) * | +----------------------------------------------------------------* * | | Check if it is a neutrino, if so force the interaction * | | (unless the relevant flag has been disabled) IF ( LISNUT (IJBEAM) .AND. LNUFIN ) THEN LFRPHN (NPFLKA) = .TRUE. * | | * | +----------------------------------------------------------------* * | | Not a neutrino ELSE LFRPHN (NPFLKA) = .FALSE. END IF * | | * | +----------------------------------------------------------------* END IF * | * +-------------------------------------------------------------------* * From this point ..... * Particle generation (1 for primaries) LOFLK (NPFLKA) = 1 * User dependent flag: LOUSE (NPFLKA) = 0 * No channeling: KCHFLK (NPFLKA) = 0 ECRFLK (NPFLKA) = ZERZER * Extra infos: INFSTK (NPFLKA) = 0 LNFSTK (NPFLKA) = 0 ANFSTK (NPFLKA) = ZERZER * Parent variables: IPRSTK (NPFLKA) = 0 EKPSTK (NPFLKA) = ZERZER * 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 AKNSHR (NPFLKA) = -TWOTWO * ahead particle Moment here * Polarization cosines: TXPOL (NPFLKA) = -TWOTWO TYPOL (NPFLKA) = +ZERZER TZPOL (NPFLKA) = +ZERZER * Particle coordinates xbeam=0 ybeam=0 zbeam=0 XFLK (NPFLKA) = XBEAM YFLK (NPFLKA) = YBEAM ZFLK (NPFLKA) = ZBEAM * Calculate the total kinetic energy of the primaries: don't change * +-------------------------------------------------------------------* * | (Radioactive) isotope: IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN * | * +-------------------------------------------------------------------* * | Heavy ion: ELSE IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. & ILOFLK (NPFLKA) .GT. 100000 ) THEN TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA) * | * +-------------------------------------------------------------------* * | Standard particle: 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