*$ CREATE SOURCE.FOR *COPY SOURCE * *=== source ===========================================================* * SUBROUTINE SOURCE (NOMORE) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Created on 07 january 1990 by Alfredo Ferrari & Paola Sala * * Infn - Milan * * * * Last change on 21-jun-98 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 '(AACOLL)' INCLUDE '(BEAM)' INCLUDE '(EPISOR)' INCLUDE '(FHEAVY)' INCLUDE '(PAPROP)' INCLUDE '(LTCLCM)' INCLUDE '(STACK)' INCLUDE '(STARS)' * LOGICAL LFIRST * SAVE LFIRST DATA LFIRST / .TRUE. / COMMON /XXXEVENTSTACKXXX/IGLPRTCNTXX,JTRACKSTRXX(5000) COMMON /XXEVENTSTACKXX/ENGYXX(5000),AGETXX(5000), & XTRCKXX(5000),YTRCKXX(5000),ZTRCKXX(5000), & DXTRCKXX(5000),DYTRCKXX(5000),DZTRCKXX(5000) * *======================================================================* * * * BASIC VERSION * * * *======================================================================* NOMORE = 0 * +-------------------------------------------------------------------* * | First call initializations: IF ( LFIRST ) THEN * | *** The following 3 cards are mandatory *** c OPEN(99,FILE='TEST.DAT',STATUS='NEW') OPEN(98,FILE='/home/clem/myFlukaArea/TEST.DAT',STATUS='OLD') TKESUM = ZERZER LFIRST = .FALSE. LUSSRC = .TRUE. IGLPRTCNTXX=0 NUMEVENTS2=0 c iseed=6161431 c iseed2=1431 read(98,*) iseed,iseed2 close(98) write(99,*) iseed,iseed2 iseed3=iseed2/2 -1 izeroxx=0 * | *** User initialization *** END IF IF (NUMEVENTS2.GT.0 ) THEN IF(IGLPRTCNTXX.GT.0) THEN WRITE(99,500) IJHION,xxpbeam,tinx,tiny,tinz,IGLPRTCNTXX, & izeroxx,NUMEVENTS2 DO LLL=1,IGLPRTCNTXX WRITE(99,501) JTRACKSTRXX(LLL), & ENGYXX(LLL),AGETXX(LLL), & XTRCKXX(LLL),YTRCKXX(LLL),ZTRCKXX(LLL), & DXTRCKXX(LLL),DYTRCKXX(LLL),DZTRCKXX(LLL) ENDDO IGLPRTCNTXX=0 ENDIF 500 FORMAT (1X,I7,1X,4(1X,E13.6),3(1X,I9)) 501 FORMAT (1X,I7,1X,2(1X,E13.7),3(1X,F11.0),3(1X,E11.5)) endif NUMEVENTS2=NUMEVENTS2+1 * | * +-------------------------------------------------------------------* * 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 * Lstack is the stack counter: of course any time source is called it * must be =0 LSTACK = LSTACK + 1 * Wt is the weight of the particle WT (LSTACK) = ONEONE WEIPRI = WEIPRI + WT (LSTACK) * Particle type (1=proton.....). Ijbeam is the type set by the BEAM * card IF ( IJBEAM .EQ. -2 ) THEN c XXSAMPXX = xxRAN3xx(ISEED3) c IF(XXSAMPXX.LT.0.91) THEN IPROZ=2 IPROA=4 c ELSEIF(XXSAMPXX.GE.0.91.AND.XXSAMPXX.LT.0.973) THEN c IPROZ=7 c IPROA=14 c ELSEIF(XXSAMPXX.LT.0.99.AND.XXSAMPXX.GE.0.973) THEN c IPROZ=14 c IPROA=28 c ELSEIF(XXSAMPXX.GE.0.99) THEN c IPROZ=26 c IPROA=56 c ENDIF IJHION = IPROZ * 1000 + IPROA IJHION = IJHION * 100 + KXHEAV c write(99,*) iproa,iproz,kxheav,ijhion IONID = IJHION CALL DCDION ( IONID ) CALL SETION ( IONID ) ILO (LSTACK) = IJHION * | * +-------------------------------------------------------------------* * | Normal hadron: ELSE ILO (LSTACK) = IJBEAM END IF * From this point ..... * Particle generation (1 for primaries) LO (LSTACK) = 1 * User dependent flag: LOUSE (LSTACK) = 0 * User dependent spare variables: DO 100 ISPR = 1, MKBMX1 SPAREK (ISPR,LSTACK) = ZERZER 100 CONTINUE * User dependent spare flags: DO 200 ISPR = 1, MKBMX2 ISPARK (ISPR,LSTACK) = 0 200 CONTINUE * Save the track number of the stack particle: ISPARK (MKBMX2,LSTACK) = LSTACK NPARMA = NPARMA + 1 NUMPAR (LSTACK) = NPARMA NEVENT (LSTACK) = 0 DFNEAR (LSTACK) = +ZERZER * ... to this point: don't change anything * Particle age (s) AGESTK (LSTACK) = +ZERZER AKNSHR (LSTACK) = -TWOTWO * Group number for "low" energy neutrons, set to 0 anyway IGROUP (LSTACK) = 0 * Kinetic energy of the particle (GeV) call PSP_LSI(ISEED,xxpbeam) xxpbeamxx = xxpbeam c write(99,*) pbeam c xxpbeam = float(ichrge(IJBEAM)) * xxpbeam xxpbeam = float(IPROZ) * xxpbeam c WRITE(99,*) IPROZ,IPROA,AM(IJbEAM) TKE (LSTACK) = SQRT ( xxPBEAM**2 + AM (IJBEAM)**2 ) - AM (IJBEAM) * Particle momentum PMOM (LSTACK) = xxPBEAM c PMOM (LSTACK) = PBEAM * PMOM (LSTACK) = SQRT ( TKE (LSTACK) * ( TKE (LSTACK) + TWOTWO * & * AM (ILO(LSTACK)) ) ) * Cosines (tx,ty,tz) c tinx=0.0 c tiny=0. c tinz=-1.0 COSTH = SQRT(ONEONE-0.9*xxRAN2xx(ISEED2)) IF(COSTH.GT.1.0) COSTH=1.0 ANG = TWOPIP*xxRAN2xx(ISEED2) TINX=SQRT(ONEONE-COSTH*COSTH)*COS(ANG) TINY=SQRT(ONEONE-COSTH*COSTH)*SIN(ANG) TINZ=-COSTH TX (LSTACK) = TINX TY (LSTACK) = TINY TZ (LSTACK) = TINZ * TZ (LSTACK) = SQRT ( ONEONE - TX(LSTACK)**2 - TY(LSTACK)**2 ) * Polarization cosines: TXPOL (LSTACK) = -TWOTWO TYPOL (LSTACK) = +ZERZER TZPOL (LSTACK) = +ZERZER * Particle coordinates xina=0.0 yina=0.0 ZINA=644814000. XA (LSTACK) = XINA YA (LSTACK) = YINA ZA (LSTACK) = ZINA c WRITE(99,500) IJHION,ichrge(IJBEAM),AM (IJBEAM),xxpbeamxx, c & TKE(LSTACK),tinx,tiny,tinz c 500 FORMAT ('P',1X,I7,1X,I2,6(1X,E13.6)) c WRITE(98,502) IJHION,ichrge(IJBEAM),pbeam,TKE(LSTACK),AM (IJBEAM), c & tinx,tiny,tinz c 502 FORMAT ('P',1X,I7,1X,I2,6(1X,E13.6)) c WRITE(*,*) IJBEAM,XINA,YINA,ZINA * Calculate the total kinetic energy of the primaries: don't change IF ( ILO(LSTACK) .EQ. -2 .OR. ILO(LSTACK) .GT. 100000 ) THEN TKESUM = TKESUM + TKE (LSTACK) * WT (LSTACK) ELSE IF ( ILO(LSTACK) .NE. 0 ) THEN TKESUM = TKESUM + ( TKE (LSTACK) + AMDISC (ILO(LSTACK)) ) & * WT (LSTACK) ELSE TKESUM = TKESUM + TKE (LSTACK) * WT (LSTACK) END IF * Here we ask for the region number of the hitting point. * NREG (LSTACK) = ... * The following line makes the starting region search much more * robust if particles are starting very close to a boundary: CALL GEOCRS ( TX (LSTACK), TY (LSTACK), TZ (LSTACK) ) CALL GEOREG ( XA (LSTACK), YA (LSTACK), ZA (LSTACK), & NREG (LSTACK), IDISC ) * Do not change these cards: CALL GEOHSM ( NHSPNT (LSTACK), 1, -11, MLATTC ) NLATTC(LSTACK) = MLATTC CALL SOEVSV CMPATH (LSTACK) = ZERZER RETURN *=== End of subroutine Source =========================================* END FUNCTION xxran1xx(idum) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL*8 xxran1xx real AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum xxran1xx=min(AM*iy,RNMX) return END C (C) Copr. 1986-92 Numerical Recipes Software 0(Z'-. FUNCTION xxran2xx(idum) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL*8 xxran2xx real AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum xxran2xx=min(AM*iy,RNMX) return END C (C) Copr. 1986-92 Numerical Recipes Software 0(Z'-. FUNCTION xxran3xx(idum) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL*8 xxran3xx real AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum xxran3xx=min(AM*iy,RNMX) return END C (C) Copr. 1986-92 Numerical Recipes Software 0(Z'-. SUBROUTINE PSP_LSI(ISEED,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL FIRST DATA FIRST/.TRUE./ IF(FIRST) THEN FIRST=.FALSE. G=2.65 C PMX=2000. PMX=19999. PC=1.0 AMX = 1./PMX**(G-1.) C=(1.-G)*(AMX - 1./PC**(G-1.) ) C=1./C PINV = 1./(G-1.) ENDIF PR=xxRAN1xx(ISEED) Y=1./(AMX+PR/(C*(G-1.)))**PINV RETURN END