*$ 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
 


