*$ 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 '(CASLIM)' INCLUDE '(BEAMCM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(IOIOCM)' INCLUDE '(LTCLCM)' INCLUDE '(PAPROP)' INCLUDE '(SOURCM)' INCLUDE '(SUMCOU)' * LOGICAL LFIRST * SAVE LFIRST DATA LFIRST / .TRUE. / * *======================================================================* * * * BASIC VERSION * * * *======================================================================* NOMORE = 0 * +-------------------------------------------------------------------* * | First call initializations: IF ( LFIRST ) THEN * | *** The following 3 cards are mandatory *** TKESUM = ZERZER LFIRST = .FALSE. LUSSRC = .TRUE. * | *** User initialization *** rt = 5.0D0 rc = 2.0D0 pathL = 700.0D0 tanmax = (rt+rc)/pathL p3_low = sqrt(1.0D0/(tanmax*tanmax+1.0D0)) p3_diff = 1.0D0 - p3_low aln10 = log(10.0D0) x1 = 3.0D0/(2.0D0*aln10) x2 = 3.0D0/aln10 A1 = (3.0D7)*exp(2.0D0*aln10/3.0D0) A2 = (3.0D5)*exp(4.0D0*aln10/3.0D0) area1 = 2.7D7 area2 = (3.0D5)*99.0D0*x1 area3 = (3.0D5)*x2 S = area1 + area2 + area3 y0 = area1/S y1 = area2/S aa1 = 1.0D0 aa2 = 4.0D0 Eth = 0.1D0 xstart = rand(97531) 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 * 40 continue phibm = 0.0D0 phi_start = 0.0D0 thetabm = 0.0D0 rbm = 0.0D0 rbefore = 0.0D0 xstart = 1.82D0*rand(0) - 0.91D0 ystart = 0.2D0*rand(0) - 0.1D0 50 continue phi_start = TWOPIP*rand(0) p1 = rand(0) p3 = p3_low + p3_diff*rand(0) rbm = rt*sqrt(p1) thetbm = acos(p3) phibm = TWOPIP*rand(0) t = pathL/p3 stheta = sqrt(ONEONE - p3*p3) dx = stheta*cos(phibm)*t dy = stheta*sin(phibm)*t xbef = xstart + dx ybef = ystart + dy rbefore = sqrt(xbef*xbef + ybef*ybef) if (rbefore.gt.2.0D0) go to 50 cost = cos(thetbm) sint = sin(thetbm) XPART = xstart YPART = ystart ZPART = 0.0D0 UPART = sint*cos(phibm) VPART = sint*sin(phibm) WPART = sqrt(ONEONE-VPART*VPART-UPART*UPART) * R = rand(0) C = (y0-R)*S/(A1*x1) + exp(-aa1/x1) D = (y0+y1-R)*S/(A2*x2) + exp(-aa2/x2) Ekin = 0.0D0 E1 = Eth + S*exp(aa1/x1)*R/A1 E2 = -x1*log(C) E3 = -x2*log(D) if (R.lt.y0) then Ekin = E1 elseif (R.lt.(y0+y1)) then Ekin = E2 else Ekin = E3 endif Energ = Ekin + AM(IONID) PPART = sqrt(Ekin*Ekin + TWOTWO*Ekin*AM(IONID)) * if ((ncase.ge.515).and.(ncase.le.520)) go to 40 if (ncase.gt.500) 1 write(*,'(i5,6e12.4)') ncase,xpart,ypart,upart,vpart,ppart,energ * PPART = sqrt(Energ*Energ - AM(IONID)**2) * 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 ( PPART**2 + AM (IONID)**2 ) - AM (IONID) * Particle momentum PMOFLK (NPFLKA) = PPART * PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA) * & + TWOTWO * AM (IM) ) ) * Cosines (tx,ty,tz) TXFLK (NPFLKA) = UPART TYFLK (NPFLKA) = VPART TZFLK (NPFLKA) = WPART * 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) = XPART YFLK (NPFLKA) = YPART ZFLK (NPFLKA) = ZPART * 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