*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
      SUBROUTINE SOURCE ( NOMORE )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2021      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-Apr-21    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)'

      INCLUDE '(CASLIM)'                  ! NCASES ....
*
      LOGICAL LFIRST, LISNUT
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*
*  ....... ==== START of my content =========
*
*  from ../Prod_Real/  : target : W + Iron refl/  
*                        5*5*30 W-core     4*5*30 - moderator
*  summary : rSUM/  ...      work10-stend/run4/Prod_Real/rSUMdat/
*
*                    C3 (0.5 m) :  BIG-C3-81-all.dat - 52167506  lines
*                    C4 (5.5 m) :  BIG-C4-91-all.dat -   641597  lines
*   i.e.
*          ../Prod_Real/rSUMdat/BIG-C3-81-all.dat
*          ../Prod_Real/rSUMdat/BIG-C4-91-all.dat


      PARAMETER (NMAX=10000000)           ! - for BIG-input, total
      COMMON /ICOM/ NNN, MYIJ(NMAX)
      COMMON /COM/ ERG(NMAX), TOTL(NMAX), ORIGL(NMAX), 
     &             TT(NMAX), XXX(NMAX), YYY(NMAX), ZZZ(NMAX),
     &             COSX(NMAX), COSY(NMAX), COSZ(NMAX)

c      INTEGER NMAX
c      PARAMETER (NMAX=10000000)           ! - for BIG-input, total
c      INTEGER NNN, MYIJ(NMAX)
c      REAL   ERG(NMAX), TOTL(NMAX), ORIGL(NMAX)
c      REAL   TT(NMAX), XXX(NMAX), YYY(NMAX), ZZZ(NMAX)
c      REAL COSX(NMAX), COSY(NMAX), COSZ(NMAX)
ccc      REAL*8 COSX(NMAX), COSY(NMAX), COSZ(NMAX)
c      COMMON /ICOM/ NNN, MYIJ
c      COMMON /COM/ ERG, TOTL,ORIGL, TT, XXX,YYY,ZZZ,
c     &              COSX,COSY,COSZ

      SAVE /ICOM/
      SAVE /COM/
*
*    NNN  - # of lines
*    MYIJ -  - part. type                 - replace  IJBEAM here
*
*     ERG - energy, GeV 
*    TOTL -   Ltot, cm                  ...for C4 Ltot points to C3 !
*    ORIGL-  Lorig, cm
*      TT -  time, sec
*   XXX, YYY, ZZZ - coordinates,  cm  
*   COSX, COSY, COSZ - normalized angles (cos)
*

      INTEGER NENTRY                   ! counter
      SAVE    NENTRY

      CHARACTER*250 LINE
      CHARACTER*80  FNAME



*   ....... ==== END OF my content =========


*
*  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 ***   ======...my entry : =========

         CALL READ_FILE          ! --- read file and fill common block

*
* check input parameters from SOURCE card :
*

        IPAR1 = NINT (WHASOU(1))
        IPAR2 = NINT (WHASOU(2))
        IPAR3 = NINT (WHASOU(3))

        WRITE(LUNOUT,*)         
        WRITE(LUNOUT,*) " ------             my  output ------------"         
        WRITE(LUNOUT,*) "----- PARAMETERS : "         
        WRITE(LUNOUT,*) "------------------- "         
        WRITE(LUNOUT,*) "- 1 - : ", IPAR1         
        WRITE(LUNOUT,*) "- 2 - : ", IPAR2         
        WRITE(LUNOUT,*) "- 3 - : ", IPAR3         
        WRITE(LUNOUT,*)         
        WRITE(LUNOUT,*) " - NCASES - : ", NCASES        
        WRITE(LUNOUT,*) " - NCASE  - : ", NCASE        
        WRITE(LUNOUT,*) " - MCASE  - : ", MCASE        
        WRITE(LUNOUT,*) " - NSOUIT - : ", NSOUIT        
        WRITE(LUNOUT,*) " - NCOINC - : ", NCOINC        
        WRITE(LUNOUT,*) " - MCOINC - : ", MCOINC        
        WRITE(LUNOUT,*)         
        WRITE(LUNOUT,*)         

        WRITE(LUNOUT,*)         !  put empty line in *out  file

        NENTRY = 0              ! my event counter

*        ==========... end of my ENTRY .. ==========

      END IF
*  |


*  +-------------------------------------------------------------------*
*  ----------- prepare my next event -----------
*   conver my variables into    FLUKA  ones .....
*
******   using RNDM  ...................
*      ISAMPL = FLRNDM (DUMMY) * NLINES 
*      ISAMPL = ISAMPL + 1 
*

*  Problem ! each cycle start from the 1-st event !
*        ...how to find cycle number ?!!! 
*  when I use consequtive method : evnet by event ....

*  ==> only RNDM method !?
*           ===========

******   using RNDM  ...................
c      NENTRY = FLRNDM (DUMMY) * NNN
c      NENTRY = NENTRY + 1 
c      IF (NENTRY .GT. NNN) NENTRY = NNN

******   use another variable ! using RNDM  ...................
c      IF (IGEN .LT. NGEN) GO TO 11
cccccc      IF (IGEN.NE.0 .AND. IGEN.LT.NGEN) GO TO 11
c      IF (IGEN.GT.0 .AND. IGEN.LE.NGEN) GO TO 11

c 10   CONTINUE                         ! start new generation...

*  **************************************
*  ...select input particle and
* ... generate new particles for the NEW selected  particle...

cccc      IENTRY = FLRNDM (DUMMY) * NNN
cccc      IENTRY = IENTRY + 1 
cccc      IF (IENTRY .GT. NNN) IENTRY = NNN

c      CALL GENER (NENTRY)          ! fill new array
c      CALL GENER (IENTRY)          ! fill new generated array
cccc      IFL= IGENER (IENTRY)          ! fill new generated array
                                   ! = 0 - simple generation
                                   ! = 1 - general generation

cccc      IGEN = 0                     ! calculator of used generated part.

cccc      IF (NENTRY .LT. 21) THEN      ! ...check variables...
cccc         WRITE(LUNOUT,*)
c         WRITE(LUNOUT,"(A,I7,1x,A,I7,1x,A,I7)") 
cccc         WRITE(LUNOUT,"(4(A,I7.1x))") 
cccc     &      " ==> check: Ient= ",IENTRY,"  Nent= ",NENTRY,
cccc     &      " NNN= ",NNN," IFL= ",IFL
c         WRITE(LUNOUT,*)
ccc      END IF
*
*  ****************************************
* ... and now use generated particles and NENTRY variable !
*
cccc 11   CONTINUE                  ! select next generation...
cccc      
cccc      NENTRY = NENTRY + 1        ! my calculator of events
cccc      IGEN = IGEN + 1            ! gener. calculator
cccc
cccc      IJBEAM = MYIJ(IENTRY)                   ! particle index
cccc      
cccc      XBEAM = XGEN(IGEN)
cccc      YBEAM = YGEN(IGEN)
cccc      ZBEAM = ZGEN(IGEN)
cccc
cccc      UBEAM = CXGEN(IGEN)
cccc      VBEAM = CYGEN(IGEN)
cccc      WBEAM = CZGEN(IGEN)
*      UBEAM = SQRT ( ONEONE - YC**2 - ZC**2 )  ! normalise vector: XC
cccc
cccc      TIMEOLD= TOLD(IGEN)
cccc      TIMEH = TGEN(IGEN)
cccc
cccc      EKIN  = EGEN(IGEN)
cccc      OL = ALGEN(IGEN)
cccc
cccc      IF (NENTRY .LT. 21) THEN      ! ...check variables...
cccc         WRITE(LUNOUT,*)
c         WRITE(LUNOUT,"(A,I7,1x,A,I7,1x,A,I7,1x,A,I7)") 
cccc         WRITE(LUNOUT,"(6(A,I7,1x))") 
cccc     &      " ...> check2: Ient= ",IENTRY,"  Nent= ",NENTRY,
cccc     &            " NPFLKA= ",NPFLKA," NCASE= ",NCASE,
cccc     &      " IGEN= ",IGEN," IJ= ",IJBEAM
*         WRITE(LUNOUT,"(A,E12.4,1x,A,E12.4,1x,A,E12.4,1x,A,E12.4)") 
cccc         WRITE(LUNOUT,"(5(A,E12.4,1x))") 
cccc     &      "        Ekin= ",EKIN,"  T= ",TIMEH," Told= ",TIMEOLD,
cccc     &               " L= ",OL," X= ",XBEAM
cccc         WRITE(LUNOUT,*)
cccc      END IF


*  ****************************************
* ... old varaible definitions...
*
cccccc      NENTRY =  NENTRY + 1
c
      NENTRY = FLRNDM (DUMMY) * NNN
      NENTRY = NENTRY + 1 
      IF (NENTRY .GT. NNN) NENTRY = NNN

      IJBEAM = MYIJ(NENTRY)                   ! particle index
      
      XBEAM = XXX(NENTRY)
      YBEAM = YYY(NENTRY)
      ZBEAM = ZZZ(NENTRY)

      UBEAM = COSX(NENTRY)
      VBEAM = COSY(NENTRY)
      WBEAM = COSZ(NENTRY)
*      UBEAM = SQRT ( ONEONE - YC**2 - ZC**2 )  ! normalise vector: XC

      TIMEH = TT(NENTRY)
ccc      TIMEOLD = TT(NENTRY)

      EKIN  = ERG(NENTRY)
      OL = ORIGL(NENTRY)

* ----------
* NCASE     = number of primaries treated so far (including the current one)

**      IF ( (NENTRY.EQ.1) .OR. (MOD (NENTRY, 1000) .EQ. 0) ) THEN
      IF ( (NCASE.EQ.1) .OR. (MOD (NCASE, 1000) .EQ. 0) ) THEN
        WRITE(LUNOUT,*)
        WRITE(LUNOUT,*) " ==== SOURCE: Nentry=  ",NENTRY," NPFLKA= ",
     &                 "  NPFLKA= ",NPFLKA, " NCASE= ",NCASE
c        WRITE(LUNOUT,*) "      IJ, Ekin, L =   ",IJBEAM,"..",EKIN,
c     &                 "..",OL
        WRITE(LUNOUT,"(A,I7,1x,A,E12.4,1x,A,E12.4)") 
     &   "      IJ, Ekin, L =   ",IJBEAM,"..",EKIN,"..",OL
      END IF  
      

*
* ===>  now use FLUKA counter  NPFLKA, stack counter 

*  ---------   next event ready ----------------
*  +-------------------------------------------------------------------*
*
*  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)

*  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 or lepton:    -  ! it is my case : N  and G ....
*       try to use FLUKA variables ...
*

      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
         AGESTK (NPFLKA) = TIMEH                 !  my time def.


*  |  Kinetic energy of the particle (GeV)
***         TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 )
***     &                   - AM (IONID)
         TKEFLK (NPFLKA) = EKIN                  !  my energy def.

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

*  Cosines (tx,ty,tz)
      TXFLK  (NPFLKA) = UBEAM
      TYFLK  (NPFLKA) = VBEAM
      TZFLK  (NPFLKA) = WBEAM
*     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) = 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.
     &          ABS (ILOFLK (NPFLKA)) .GE. KPIOMN ) 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



***************  MY entry ***********************************
*=======================================================
*=======================================================


      SUBROUTINE READ_FILE 
*               ---------
*       use my common and fill it with input data
* --------------------------------

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'

      INCLUDE '(CASLIM)'                  ! NCASES ....

* ------------------------------------------------
*  from ../Prod_Real/  : target : W + Iron refl/  
*                        5*5*30 W-core     4*5*30 - moderator
*  summary : rSUM/  ...      work10-stend/run4/Prod_Real/rSUMdat/
*
*                    C3 (0.5 m) :  BIG-C3-81-all.dat - 52167506  lines
*                    C4 (5.5 m) :  BIG-C4-91-all.dat -   641597  lines
*   i.e.
*          ../Prod_Real/rSUMdat/BIG-C3-81-all.dat
*          ../Prod_Real/rSUMdat/BIG-C4-91-all.dat


      PARAMETER (NMAX=10000000)           ! - for BIG-input, total
      COMMON /ICOM/ NNN, MYIJ(NMAX)
      COMMON /COM/ ERG(NMAX), TOTL(NMAX), ORIGL(NMAX), 
     &             TT(NMAX), XXX(NMAX), YYY(NMAX), ZZZ(NMAX),
     &             COSX(NMAX), COSY(NMAX), COSZ(NMAX)

c      INTEGER NMAX
c      PARAMETER (NMAX=10000000)           ! - for BIG-input, total
c      INTEGER NNN, MYIJ(NMAX)
c      REAL   ERG(NMAX), TOTL(NMAX), ORIGL(NMAX)
c      REAL   TT(NMAX), XXX(NMAX), YYY(NMAX), ZZZ(NMAX)
c      REAL COSX(NMAX), COSY(NMAX), COSZ(NMAX)
ccc      REAL*8 COSX(NMAX), COSY(NMAX), COSZ(NMAX)
c      COMMON /ICOM/ NNN, MYIJ
c      COMMON /COM/ ERG, TOTL,ORIGL, TT, XXX,YYY,ZZZ,
c     &              COSX,COSY,COSZ

      SAVE /ICOM/
      SAVE /COM/

*
*    NNN  - # of lines
*    MYIJ -  - part. type                 - replace  IJBEAM here
*
*     ERG - energy, GeV 
*    TOTL -   Ltot, cm                  ...for C4 Ltot points to C3 !
*  TOTORIGL-  Lorig, cm
*      TT -  time, sec
*   XXX, YYY, ZZZ - coordinates,  cm  
*   COSX, COSY, COSZ - normalized angles (cos)
*

* ----------------------------------------------


c      CHARACTER*250 LINE
      CHARACTER*160 LINE
      CHARACTER*80  FNAME

*======================================================================

*======================================================================
*======================================================================
*
*                     get input file methods :
*
* - here we use standart method : define open channel from OPEN card
*
* - open file right here :
*        WRITE(LUNOUT,*)
*        CALL OAUXFI('phsp3.dat',LUNRDB,'OLD',IERR)
*        OPEN(UNIT=89,FILE="testsource.dat",STATUS="UNKNOWN")
*        WRITE(LUNOUT,*)
*         WRITE(LUNOUT,*)
*         WRITE(LUNOUT,'(A,132A)') ("*",I=1,132)
*         WRITE(LUNOUT,*)
*
* OR
*
*         IF(NINT(WHASOU(1)) .EQ. 1) THEN
*            OPEN(UNIT = 81, FILE = "../sr.spectrum", 
*     &         STATUS = "OLD")
*            WRITE(LUNOUT,'(A,A,A)') "    sr white spectrum ",
*     &      "(file sr.spectrum)" 
*         END IF
*         WRITE(LUNOUT,*)
*         WRITE(LUNOUT,*)
*         WRITE(LUNOUT,'(A,132A)') ("*",I=1,132)
*         WRITE(LUNOUT,*)
* 


*** NOTE:  from 3 methods to read input file with data 
***           ONLY the 3-d one works !

*======================================================================

*
*                    .... my start - read file ...
*                      -------------------------
****  |  WHASOU(1) = UNIT number (File name is linked with OPEN card)
***         LUNRD = NINT(WHASOU(1))
****         POSMIN = WHASOU(2)
****         POSMAX = WHASOU(3)

*-----------------------------
**          
**            WRITE(LUNOUT,*) "-my- ",LUNRD
**
***            LUNOUT= 91    ! ...see parameter
**
***            OPEN(UNIT = 81, FILE = "../sr.spectrum", 
**            OPEN(LUNRD, FILE = "test.dat",   IOSTAT= IOS,
**     &         STATUS = "OLD")
**
**            WRITE(LUNOUT,*) "-my- open: iostat = ",IOS
**
***            WRITE(LUNOUT,'(A,A,A)') "    sr white spectrum ",
***     &      "(file sr.spectrum)" 
**
**
**
*----------------------------
*               another method to open file 
*

        WRITE(LUNOUT,*)         !  put empty line in *out  file

c        CALL OAUXFI(' ../Prod_Real/rSUMdat/BIG-C4-91-all.dat', 
c     &                LUNRD, 'OLD', IERR)
****        CALL OAUXFI(' ../BIG-C4-91-all.dat', LUNRD, 'OLD', IERR)

c        FNAME = "../BIG-C4-91-all.dat"   
        FNAME = "../BIG-C4-91-all.dat"    ! Direct part. for 50*50 window
        CALL OAUXFI(FNAME, LUNRD, 'OLD', IERR)


        WRITE(LUNOUT,*)
        WRITE(LUNOUT,*) " ---- open input file: events : ",FNAME
        WRITE(LUNOUT,*)
        WRITE(LUNOUT,*) " -my- : ===== OPEN file :  ierr= ", IERR
        WRITE(LUNOUT,'(A,132A)') ("*",I=1,132)   ! line with **/...
        WRITE(LUNOUT,*)

*----------------------------

*         WTOT = ZERZER               ! internal  - total weight
         NNN = 0                     ! internal ... counter

 10      CONTINUE                   ! start loop
************************** 

*            READ( LUNRD, '(A)', ERR=9999, END=20 ) LINE
            READ( LUNRD, '(A)', ERR=9999, END=20, IOSTAT= IOS) LINE

            IF (NNN .LT. 10) THEN
              WRITE(LUNOUT,*) "-......my- read : ios =",IOS
              WRITE(LUNOUT,'(A)') LINE
***            WRITE(LUNOUT,*) "***", LINE
            ENDIF

***   **************************************
***    input format :
***    Det - Type -  Nev - Ekin - Ltot- L.Orig - Time - w - X - Y - Z  - 
***                                                    cosX - cosY - cosZ 

***            now read line ... :
***
***      ...ignore                     Det,     Nev,     w 
***  ****************************
*            READ (LINE,*,ERR=10) 
            READ (LINE,*, IOSTAT= IOS,ERR=10) 
     &         IDD, IP, NN, EK, TL, TO, T, W, X, Y, Z, XC, YC, ZC
*     &         I, IDD, DUM, DUM, IP, EK, AL, T, W, X, Y, Z, XC, YC, ZC

*            WRITE(LUNOUT,*) "***", IOS
*   **********************************8

*           ... now apply some cuts :  X = [598.9 - 599.1}
*                                      XC > 0.9  (forward direction only)
*
*            IF ( X.LT.598.9 .OR. X.GT.599.1 ) GOTO 10
*            IF ( XC .LT. 0.9 ) GOTO 10

*
*           ... now save good event ....
*
            NNN = NNN + 1

            IF (NNN.GT.NMAX) CALL FLABRT('SOURCE','Increase NMAX')

            MYIJ(NNN)  = IP                ! part. type
            ERG(NNN)   = EK                !  kinetic energy, GeV  
            TOTL(NNN)  = TL                !  Ltot, cm 
           ORIGL(NNN)  = TO                !  Lorig, cm 
              TT(NNN)  = T                 ! time, sec 

            XXX(NNN) = X                   ! coordinates
            YYY(NNN) = Y
            ZZZ(NNN) = Z

*** C4 has X= 549.6 - 540.8, ... so to fill C4 set X= 549.4

c            XNEW = 549.4
            XNEW = 549.0
            XXX(NNN) = XNEW


*           ... make normalization of forward direction vector:
*               input angles are :   YC  and ZC,,,, XC is calculated !
*
            COSY(NNN) = YC
            COSZ(NNN) = ZC
            COSX(NNN) = SQRT ( ONEONE - YC**2 - ZC**2)

*  -----------------------------

            IF (NNN .LT. 10) THEN
          WRITE(LUNOUT,* ) "-N- :", NNN, MYIJ(NNN), ERG(NNN), TT(NNN), 
     &                    XXX(NNN), YYY(NNN), ZZZ(NNN), COSX(NNN)
            ENDIF


         GOTO 10                    ! read next ev.
*********************************

 20      CONTINUE                   ! end of ev. loop

         IF (NNN.EQ.0) CALL FLABRT('SOURCE','Error reading file')
         WRITE (LUNOUT,*)
         WRITE (LUNOUT,*) '*** input file:  ',NNN,' particles loaded'
         WRITE (LUNOUT,*)

*                        ---------------------------
*                    .... my start - end of read file ...
*
      CLOSE( LUNRD )

      RETURN

*=== End of subroutine READ_FILE =========================================*
*                  ... my : end of job .....
 9999 CONTINUE
      CALL FLABRT('SOURCE-READ_FILE','Error reading source file')
*                  ..........................

      END
*===================================================
*=======================================================

