***
*     gxuser.f
***
*
      SUBROUTINE UGINIT
*
************************************************************************
*                                                                      *
*             To initialise GEANT3.21 program and read data cards      *
*                                                                      *
************************************************************************
*
      PARAMETER (NWGEAN=9000000,NWPAW=2000000,NWKUIP=30000)
      COMMON/GCBANK/GEANT(NWGEAN)
      COMMON/PAWC/PAW(NWPAW)
C
*KEEP,GCKINE.
      COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP
     +      ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD
C
      COMMON/GCUNIT/LIN,LOUT,NUNITS,LUNITS(5)
      INTEGER LIN,LOUT,NUNITS,LUNITS
      COMMON/GCMAIL/CHMAIL
      CHARACTER*132 CHMAIL
C
*KEEP,GCCUTS.
      COMMON/GCCUTS/CUTGAM,CUTELE,CUTNEU,CUTHAD,CUTMUO,BCUTE,BCUTM
     +             ,DCUTE ,DCUTM ,PPCUTM,TOFMAX,GCUTS(5)
C
*KEEP,GCPHYS.
      COMMON/GCPHYS/IPAIR,SPAIR,SLPAIR,ZINTPA,STEPPA
     +             ,ICOMP,SCOMP,SLCOMP,ZINTCO,STEPCO
     +             ,IPHOT,SPHOT,SLPHOT,ZINTPH,STEPPH
     +             ,IPFIS,SPFIS,SLPFIS,ZINTPF,STEPPF
     +             ,IDRAY,SDRAY,SLDRAY,ZINTDR,STEPDR
     +             ,IANNI,SANNI,SLANNI,ZINTAN,STEPAN
     +             ,IBREM,SBREM,SLBREM,ZINTBR,STEPBR
     +             ,IHADR,SHADR,SLHADR,ZINTHA,STEPHA
     +             ,IMUNU,SMUNU,SLMUNU,ZINTMU,STEPMU
     +             ,IDCAY,SDCAY,SLIFE ,SUMLIF,DPHYS1
     +             ,ILOSS,SLOSS,SOLOSS,STLOSS,DPHYS2
     +             ,IMULS,SMULS,SOMULS,STMULS,DPHYS3
     +             ,IRAYL,SRAYL,SLRAYL,ZINTRA,STEPRA
      COMMON/GCPHLT/ILABS,SLABS,SLLABS,ZINTLA,STEPLA
     +             ,ISYNC
     +             ,ISTRA
C
*KEEP,GCFLAG.
      COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN
     +        ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2)
      COMMON/GCFLAX/BATCH, NOLOG
      LOGICAL BATCH, NOLOG
C
*KEEP,user common block
      COMMON/MKMAG/RMKF
      integer RMKF
*KEND.
cn..          Define of HBOOK word
      integer   istat
      character chtags(5)*5
      data      chtags/'I','NHIT','P','theta','phi'/
      character tagkine(12)*6
      data      tagkine/'Pc','Pxc','Pyc','Pzc','Thec','Phic',
     +                   'P','Px','Py','Pz','theta','phi'/
C
***   -----------------------------------------------------------------
*
*
*             Initialize GEANT
*
*
      CALL GINIT           ! Initialize GEANT system
*
* --- Set the number of significant characters for FFREAD keys to 5 ---

      CALL FFSET('SIZE',5)

* --- Define GWA98 specific FFREAD cards ---

      CALL SXKEY

* --- Read all FFREAD data cards ---

CHK We comment this out for now.
CHK      OPEN(25,FILE='RUNCARD.DEF',STATUS='OLD')
c      CALL GFFGO
CHK      CLOSE(25)

*
*
*             Input tracking thresholds
*
*
      CUTGAM=0.001 ! threshold for gammas transport
      CUTELE=0.001 ! threshold for electrons transport
      CUTHAD=0.01  ! threshold for hadrons   transport
      CUTNEU=0.01  ! threshold for neutral hadrons transport
      CUTMUO=0.01  ! threshold for muons transport
      BCUTE=CUTGAM ! threshold for photons produced by electron Brems.
      BCUTM=CUTGAM ! threshold for photons produced by muon Brems.
      DCUTE=CUTELE ! threshold for electron produced by electron&delta-rays
      DCUTM=CUTELE ! threshold for electron produced by muon delta-rays
      PPCUTM=0.002 ! threshold for e+e- pair production by muons
C      TOFMAX=1010  ! threshold on time of flight counted
C                                             from primary interaction time
*
*
*             Input physics processes
*
*
      IPAIR=0 ! pair production with generation of electron/positron
      ICOMP=0 ! compton scattering with generation of electron
      IPHOT=0 ! photo-electric effect with generation of electron
      IPFIS=0 ! photo-fission effect with generation of secondaries
      IDRAY=0 ! delta ray production with generation of electron
      IANNI=0 ! positron annihilation with generation of photons
      IBREM=0 ! bremsstrahlung with generation of gamma
      IHADR=0 ! hadronic interactions with generation of secondaries
      IMUNU=0 ! muon nuclear interaction with generation of secondaries
      IDCAY=0 ! decay in flight with generation of secondaries
      ILOSS=0 ! continuous energy loss without generation of delta ray 
C                    and Landau-Vaviliov-Gauss fluctuations
      IMULS=0 ! multiple scattering according to Moli'ere theory
      IRAYL=0 ! Rayleigh scattering
      ILABS=0 ! absorbtion of Cerenkov photons with posssible detection
      ISYNC=0 ! To control synchrotron radiation
      ISTRA=0 ! To control energy loss fluctuation model
**
c      write(6,*)' '
c      write(6,*)'      Type in !'
c      write(6,*)' '
*
**
*
*
*             Initialize GEANT/ZBOOK data structures
*
*
      CALL GZINIT          ! Initialize mother banks
*          
*
*             Initialize graphics package
*
*
      CALL GDINIT          ! Initialize drawing packag
*
      CALL GPART           ! Store the standard particle constants
c      CALL GMATE           ! Initialize material constants
*
      CALL UGEOM           ! Define user geometry
      CALL GPHYSI          ! Initialize GEANT physics processes
*
c      CALL GPTMED          ! Prints the tracking medium parameters
*
**
*
*             Open HBOOK
*
      CALL HROPEN(1,'WA98','wa98.nt','n',1024,istat)
**
*
*             Booking Histogram & Ntuple file
*
      CALL HBOOKN(10,'Momentum of pion+',5,'WA98',5000,chtags)
      CALL HBOOKN(50,'Check of kinematics',12,'WA98',5000,tagkine)
C
*
      END
*
***
*
      SUBROUTINE UGLAST
*
************************************************************************
*                                                                      *
*             Termination routine to print histograms and statistics   *
*                                                                      *
************************************************************************
*
      COMMON/GCOMIS/ICOMIS,JUINIT,JUGEOM,JUKINE,JUSTEP,JUOUT,JULAST
C
*KEND.
cn..          Define of HBOOK word
      integer icycle
C
*     -----------------------------------------------------------------
*
      IF(JULAST.NE.0)THEN
         CALL CSJCAL(JULAST,0,X,X,X,X,X,X,X,X,X,X)
         GO TO 99
      ENDIF
*
      CALL GLAST
**
*
*
*             Close HIGZ
*
*
      CALL IGEND
**
*
*             Close HBOOK
*
      CALL HROUT(0,icycle,' ')
      CALL HREND('WA98')
**
*
      RETURN
 99   END
*
***
*
      SUBROUTINE GUKINE
*
************************************************************************
*                                                                      *
*             Generates Kinematics for primary tracks                  *
*                                                                      *
************************************************************************
*
*KEEP,GCFLAG.
      COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN
     +     ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2)
      COMMON/GCFLAX/BATCH, NOLOG
      LOGICAL BATCH, NOLOG
C
*KEEP,GCKINE.
      COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP
     +     ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD
C
*KEEP,GCONST.
      COMMON/GCONST/PI,TWOPI ,PIBY2,DEGRAD,RADDEG,CLIGHT ,BIG,EMASS
      COMMON/GCONSX/EMMU,PMASS,AVO
C
*KEEP,GCSCAN.
      PARAMETER (MSLIST=32,MAXMDT=3)
      COMMON/GCSCAN/SCANFL,NPHI,PHIMIN,PHIMAX,NTETA,TETMIN,TETMAX,
     +              MODTET,IPHIMI,IPHIMA,IPHI1,IPHIL,NSLMAX,
     +              NSLIST,ISLIST(MSLIST),VSCAN(3),FACTX0,FACTL,
     +              FACTR,IPHI,ITETA,ISCUR,SX0,SABS,TETMID(MAXMDT),
     +              TETMAD(MAXMDT)
      LOGICAL SCANFL
      COMMON/GCSCAC/SFIN,SFOUT
      CHARACTER*80 SFIN,SFOUT
C
*KEEP,GCOMIS.
      COMMON/GCOMIS/ICOMIS,JUINIT,JUGEOM,JUKINE,JUSTEP,JUOUT,JULAST
C
*KEEP,SCXGEO.
C --- Common which contains some general geometry parameters ---
      REAL DWA98,DWMAG,DTARG
     $,    DWA98X,DWA98Y,DWA98Z
     $,    DMAGNX,DMAGNY,DMAGNZ
     $,    DTARGX,DTARGY,DTARGZ
     $,    XTARG,YTARG,ZTARG,XMAG,YMAG,ZMAG
      COMMON /SCXGEO/ DWA98(3),DWMAG(3),DTARG(3)
     $,               DWA98X,DWA98Y,DWA98Z
     $,               DMAGNX,DMAGNY,DMAGNZ
     $,               DTARGX,DTARGY,DTARGZ
     $,               XTARG,YTARG,ZTARG,XMAG,YMAG,ZMAG
C
*KEEP,user common block
      COMMON/MKMAG/RMKF
      real       RMKF
      INTEGER    NHIT,NDETNUM,Num_EVENT/1/
      PARAMETER  (NDETNUM=10)
      REAL       DETPOS(NDETNUM,8),INITP(6)
      common/uTRACK/NHIT,Num_EVENT,DETPOS(NDETNUM,8),INITP(6),pl
C
*KEND.
C
      REAL      UBUF,PLAB(3),VERT(3)
      INTEGER   NVTX,NT,IPART
C
      real      pt, p4(4),Init_num, y_mean, y_sigma, y
      real      Init_pz, theta, phi, pl, p_theta, p_phi
      integer   Iseed, Init_PART,In_part
      integer   inum,i,j
C
*     -----------------------------------------------------------------
*
      IF(JUKINE.NE.0)THEN
         CALL CSJCAL(JUKINE,0,X,X,X,X,X,X,X,X,X,X)
         GO TO 99
      ENDIF
*
*
*             Define type, momentum and vertex of particle
*
*
      save Init_PART,Init_num,Init_pz,p_theta,p_phi
*
      IF(IEVENT.EQ.1)THEN
         WRITE(6,*)' '
         WRITE(6,*)'      Type in !'
*     
         WRITE(6,*)' '
         WRITE(6,*)'Initial particle number ..integer'
         WRITE(6,*)' <example>  1,Gamma   3,Electron   6,Muon-'
         WRITE(6,*)'            8,Pion+   9,Pion-     11,Kaon+'
         WRITE(6,*)'           12,Kaon-  14,Proton     etc.   '
         WRITE(6,*)' '
         READ(5,*)Init_PART
*
         WRITE(6,*)'Magnetic field is 1 tesla (-10 kGauss)'
         WRITE(6,*)' '
*
         WRITE(6,*)' '
         WRITE(6,*)'Initial momentum (GeV) ..real'
         READ(5,*)Init_pz
*
c         WRITE(6,*)' '
c         WRITE(6,*)'Initial angle of theta and phy (degree) ..real'
c         READ(5,*)theta,phi
c         p_theta = theta*3.1415926/180.
c         p_phi = phi*3.1415926/180.
*
c         WRITE(6,*)' '
c         WRITE(6,*)'Magnetic field (T:tesla) ..real'
c         WRITE(6,*)' 1 tesla (-10 kGauss)'
c         WRITE(6,*)' '
c         READ(5,*) RMKF
      ENDIF
      RMKF = 1.
*
      Iseed = NRNDM(2)
      VERT(1) = XTARG
      VERT(2) = YTARG
      VERT(3) = ZTARG
c      VERT(3) = ZTARG - 0.02 + 0.04*ran(Iseed)
      CALL GSVERT(VERT,0,0,UBUF,1,NVTX)
*
      if((Num_EVENT.ge.0).and.(Num_EVENT.le.100)) then
         pl = Init_pz 
      elseif((Num_EVENT.gt.10).and.(Num_EVENT.le.200)) then
         pl = Init_pz + 0.5
      elseif((Num_EVENT.gt.200).and.(Num_EVENT.le.300)) then
         pl = Init_pz + 1.0
      elseif((Num_EVENT.gt.300).and.(Num_EVENT.le.400)) then
         pl = Init_pz + 1.5
      elseif((Num_EVENT.gt.400).and.(Num_EVENT.le.500)) then
         pl = Init_pz + 2.0
      else
         pl = Init_pz + 2.5
      endif
c      pl = Init_pz
*
      theta = -5. + 10.*ran(Iseed)
      phi = -90. + 180.*ran(Iseed)
      p_theta = theta*3.1415926/180.
      p_phi = phi*3.1415926/180.
*
      PLAB(1) = pl*sin(p_theta)*cos(p_phi)
      PLAB(2) = pl*sin(p_theta)*sin(p_phi)
      PLAB(3) = pl*cos(p_theta)
      CALL GSKINE(PLAB,Init_PART,NVTX,UBUF,1,NT)
*
*     fill for initial 
*
      INITP(1)=pl
      INITP(2)=PLAB(1)
      INITP(3)=PLAB(2)
      INITP(4)=PLAB(3)
      INITP(5)=p_theta
      INITP(6)=p_phi
*
      NHIT = 0
      do i = 1,10
         do j = 1,8
            DETPOS(i,j) = 0.
         enddo
      enddo
C
*
*
*
*             Kinematic debug (controled by ISWIT(1))
*
*
      IF((IDEBUG.EQ.1).AND.(ISWIT(1).EQ.1))THEN
         CALL GPRINT('VERT',0)
         CALL GPRINT('KINE',0)
      ENDIF
*     
*
 99   END
*
***
*
      SUBROUTINE GUTREV 
C.  
C.    ******************************************************************    
C.    *                                                                *    
C.    *       User routine to control tracking of one event            *    
C.    *       Called by GRUN                                           *    
C.    *                                                                *    
C.    ******************************************************************    
C.  
C.
      CALL GTREVE
      RETURN    
      END
*
***
*
      SUBROUTINE GUTRAK
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       User routine to control tracking of one track            *
C.    *                                                                *
C.    *    ==>Called by : GTREVE                                       *
C.    *                                                                *
C.    ******************************************************************
C.
*KEEP,GCFLAG
      COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN
     +     ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2)
      COMMON/GCFLAX/BATCH, NOLOG
      LOGICAL BATCH, NOLOG
C
*KEEP,GCKINE.
      COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP
     +     ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD
      CHARACTER*4 NAPART                       
C
*KEEP,GCTRAK.
      PARAMETER (MAXMEC=30)
      COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC)
     +     ,NAMEC(MAXMEC),NSTEP ,MAXNST,DESTEP,DESTEL,SAFETY,SLENG
     +     ,STEP  ,SNEXT ,SFIELD,TOFG  ,GEKRAT,UPWGHT,IGNEXT,INWVOL
     +     ,ISTOP ,IGAUTO,IEKBIN, ILOSL, IMULL,INGOTO,NLDOWN,NLEVIN
     +     ,NLVSAV,ISTORY
      CHARACTER*4 NAMEC
      PARAMETER (MAXME1=30)
      COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1)
*KEEP,user common block

C
      real Init_ENERGY
C.
C.    ------------------------------------------------------------------
C.
      if((VECT(1) .le. 0.1).and.(VECT(2) .le. 0.1).and.
     +     (VECT(3) .le. 0.1).and.(IPART .eq. In_part))then
         Init_ENERGY = GETOT*1.0E+3   ! GETOT:(GeV)--> Init_ENERGY:(MeV)
      endif
*
      CALL GTRACK
      RETURN
      END
*
***
*
      SUBROUTINE GUFLD(X,F) 
C.  
C.    ******************************************************************    
C.    *                                                                *    
C.    *       User routine to compute the magnetic field F             *    
C.    *       at space point X                                         *    
C.    *       Called by GHELIX,GRKUTA                                  *    
C.    *                                                                *    
C.    ******************************************************************    
C.  
*KEEP,GCVOLU.
      COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15),    
     +     LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15),  
     +     GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3) 
      CHARACTER*4 NAMES
*KEEP,user common block
      COMMON/MKMAG/RMKF
      SAVE
C
      DIMENSION X(3),F(3)
C.
C.    ------------------------------------------------------------------
*
      F(1)=0.
      F(2)=-RMKF*10             ! 1 tesla = 10kGauss
      F(3)=0.
*
      RETURN
      END
*
***
*
      SUBROUTINE GUSTEP
*
************************************************************************
*                                                                      *
*             User routine called at the end of each tracking step     *
*             MEC   is the mechanism origin of the step                *
*             INWVOL is different from 0 when the track has reached    *
*                    a volume boundary                                 *
*             ISTOP is different from 0 if the track has stopped       *
*                                                                      *
*             Called by GTRACK                                         *
*                                                                      *
************************************************************************
*
*KEEP,GCFLAG
      COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN
     +     ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2)
      COMMON/GCFLAX/BATCH, NOLOG
      LOGICAL BATCH, NOLOG
C
*KEEP,GCKINE.
      COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP
     +     ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD
      CHARACTER*4 NAPART
C
*KEEP,GCKING.
      INTEGER MXGKIN
      PARAMETER (MXGKIN=100)
      COMMON/GCKING/KCASE,NGKINE,GKIN(5,MXGKIN),
     +                           TOFD(MXGKIN),IFLGK(MXGKIN)
      INTEGER       KCASE,NGKINE ,IFLGK,MXPHOT,NGPHOT
      REAL          GKIN,TOFD,XPHOT
C
      PARAMETER (MXPHOT=800)
      COMMON/GCKIN2/NGPHOT,XPHOT(11,MXPHOT)
C
      COMMON/GCKIN3/GPOS(3,MXGKIN)
      REAL          GPOS
C
*KEEP,GCMATE.
      COMMON/GCMATE/NMAT,NAMATE(5),A,Z,DENS,RADL,ABSL
C
*KEEP,GCTMED.
      COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFIELD,FIELDM,TMAXFD,STEMAX
     +     ,DEEMAX,EPSIL,STMIN,CFIELD,CMULS,IUPD,ISTPAR,NUMOLD
      COMMON/GCTLIT/THRIND,PMIN,DP,DNDL,JMIN,ITCKOV,IMCKOV,NPCKOV
C
*KEEP,GCTRAK.
      PARAMETER (MAXMEC=30)
      COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC)
     +     ,NAMEC(MAXMEC),NSTEP ,MAXNST,DESTEP,DESTEL,SAFETY,SLENG
     +     ,STEP  ,SNEXT ,SFIELD,TOFG  ,GEKRAT,UPWGHT,IGNEXT,INWVOL
     +     ,ISTOP ,IGAUTO,IEKBIN, ILOSL, IMULL,INGOTO,NLDOWN,NLEVIN
     +     ,NLVSAV,ISTORY
      CHARACTER*4 NAMEC
      PARAMETER (MAXME1=30)
      COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1)
C   
*KEEP,GCVOLU.
      COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15),    
     +     LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15),  
     +     GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3) 
      CHARACTER*4 NAMES
C
      COMMON/GCOMIS/ICOMIS,JUINIT,JUGEOM,JUKINE,JUSTEP,JUOUT,JULAST
C
*KEEP,GCONST.
      COMMON/GCONST/PI,TWOPI ,PIBY2,DEGRAD,RADDEG,CLIGHT ,BIG,EMASS
      COMMON/GCONSX/EMMU,PMASS,AVO
C
*KEEP,SCXGEO.
C --- Common which contains some general geometry parameters ---
      REAL DWA98,DWMAG,DTARG
     $,    DWA98X,DWA98Y,DWA98Z
     $,    DMAGNX,DMAGNY,DMAGNZ
     $,    DTARGX,DTARGY,DTARGZ
     $,    XTARG,YTARG,ZTARG,XMAG,YMAG,ZMAG
      COMMON /SCXGEO/ DWA98(3),DWMAG(3),DTARG(3)
     $,               DWA98X,DWA98Y,DWA98Z
     $,               DMAGNX,DMAGNY,DMAGNZ
     $,               DTARGX,DTARGY,DTARGZ
     $,               XTARG,YTARG,ZTARG,XMAG,YMAG,ZMAG
C
*KEEP,user common block
      INTEGER    NHIT,NDETNUM
      PARAMETER  (NDETNUM=10)
      REAL       DETPOS(NDETNUM,8),INITP(6)
      common/uTRACK/NHIT,Num_EVENT,DETPOS(NDETNUM,8),INITP(6),pl
*KEND.
C
      LOGICAL LH(10)/10*.FALSE./
      INTEGER MI,MMACS
*     -----------------------------------------------------------------
*
      IF(JUSTEP.NE.0)THEN
         CALL CSJCAL(JUSTEP,0,X,X,X,X,X,X,X,X,X,X)
         GO TO 99
      ENDIF
*
*
*             Print the tracking and physics
*
*
      IF(ISWIT(1).EQ.1)THEN
         CALL GPCXYZ
      ENDIF
*
*
*             Save in the stack particle generated during the current step
*
*
      IF(NGKINE.NE.0)THEN
         CALL GSKING(0)
      ENDIF
*
*
*             Stores current space point into JXYZ
*
*
      CALL GSXYZ
**
*
*             Filling Physics parameter
*
      if(Num_EVENT.ge.12001) goto 99
      if((INWVOL.eq.1).and.(CHARGE.ne.0).and.(VECT(3).ge.ZMAG))then
         if((NAMES(5) .eq. 'JPDA').and.(NUMBER(3).eq.1))then
            LH(6) = .TRUE.
            NHIT = NHIT + 1
c            write(6,*)'PAD1 = ',VECT(1) ,VECT(2) ,VECT(3) - ZMAG 
         elseif((NAMES(5) .eq. 'JPDA').and.(NUMBER(3).eq.2))then
            LH(7) = .TRUE.
            NHIT = NHIT + 1
c            write(6,*)'PAD2 = ',VECT(1) ,VECT(2) ,VECT(3) - ZMAG 
         elseif((NAMES(5) .eq. 'JSTA').and.(NUMBER(3).eq.1))then
            LH(8) = .TRUE.
            NHIT = NHIT + 1
c            write(6,*)'STR1 = ',VECT(1) ,VECT(2) ,VECT(3) - ZMAG 
         elseif((NAMES(5) .eq. 'JSTA').and.(NUMBER(3).eq.2))then
            LH(9) = .TRUE.
            NHIT = NHIT + 1
c            write(6,*)'STR2 = ',VECT(1) ,VECT(2) ,VECT(3) - ZMAG 
         endif
c..   save position for new hit
         do I=1,NDETNUM
            if(LH(I)) THEN
               DETPOS(I,1)=VECT(1) ! POS X
               DETPOS(I,2)=VECT(2) ! POS Y
               DETPOS(I,3)=VECT(3) - ZMAG ! POS Z
               DETPOS(I,4)=TOFG*1.0e+9 ! TOF (ns)
               DETPOS(I,5)=VECT(1) ! POS X with resolution
               DETPOS(I,6)=VECT(2) ! POS Y with resolution
               DETPOS(I,7)=VECT(3) - ZMAG ! POS Z with resolution
               DETPOS(I,8)=TOFG*1.0e+9 ! TOF (ns) with resolution
               LH(I) = .FALSE.
            endif
         enddo
      endif
C
      if(INWVOL .eq. 2)then
         NAMES(2)='    '
         NAMES(3)='    '
         NAMES(4)='    '
         NAMES(5)='    '
      endif
*
  99  CONTINUE
      END
*
*
***
*
      SUBROUTINE GUOUT
*
************************************************************************
*                                                                      *
*             User routine called at the end of each event.            *
*                                                                      *
************************************************************************
*
*KEEP,GCFLAG.
      COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN
     +        ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2)
      COMMON/GCFLAX/BATCH, NOLOG
      LOGICAL BATCH, NOLOG
C
      COMMON/GCOMIS/ICOMIS,JUINIT,JUGEOM,JUKINE,JUSTEP,JUOUT,JULAST
C
*KEEP,GCVOLU.
      COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15),
     +     LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15),
     +     GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3)
      CHARACTER*4 NAMES
C
*KEEP,SCXGEO.
C --- Common which contains some general geometry parameters ---
      REAL DWA98,DWMAG,DTARG
     $,    DWA98X,DWA98Y,DWA98Z
     $,    DMAGNX,DMAGNY,DMAGNZ
     $,    DTARGX,DTARGY,DTARGZ
     $,    XTARG,YTARG,ZTARG,XMAG,YMAG,ZMAG
      COMMON /SCXGEO/ DWA98(3),DWMAG(3),DTARG(3)
     $,               DWA98X,DWA98Y,DWA98Z
     $,               DMAGNX,DMAGNY,DMAGNZ
     $,               DTARGX,DTARGY,DTARGZ
     $,               XTARG,YTARG,ZTARG,XMAG,YMAG,ZMAG
C
*KEEP,SCXFF.
C --- Common which contains FFREAD stuff for the GWA98 package ---
C --- as well as the total processed event counter ---
      INTEGER NSXDET,NLUNS,NPLANS
      REAL SXGATE,SXMAGN
      INTEGER ISXDET,ISXDB1,ISXDB2,JDETF,ISXTRA,JTRAF,ISXDRG,JDRGF
     $,       ISXDRT,JDRTF,IFDRAT,ISXDCH,JDCHM,JDCHN,JDCHP,ISXLUN
     $,       ISXOUT,JOUTF,ISXPW1,ISXPW2,IMXTAR,ISXTAR,ISXTIM
     $,       ISXPLN,JPLNF,ISXWKS,ISXMAP,ISXVAC,ISXTUB,NSXEVT,ISXBG1
     $,       ISXBG2,IFPART,IFVOLU,IFMATE,IFTMED,IFVERT,IFKINE,IFSETS
     $,       IFHITS,IFDIGI,IFSECS,IFXSEC,IFLOSS,ISXEVT,ISXCHM
     $,       MAXNSTEP
     $,       ISXINP,ISXDIN,ISXCOOR

      PARAMETER (NSXDET=15,NLUNS=8,NPLANS=15)

      COMMON /SCXFF/ ISXDET(NSXDET+1),ISXDB1(NSXDET+1)
     $,              ISXDB2(NSXDET+1),JDETF(NSXDET+1)
     $,              ISXTRA(NSXDET+1),JTRAF(NSXDET+1)
     $,              ISXDRG(NSXDET+1),JDRGF(NSXDET+1)
     $,              ISXDRT(NSXDET+1),JDRTF(NSXDET+1),IFDRAT
     $,              ISXDCH,JDCHM,JDCHN,JDCHP
     $,              ISXLUN(NLUNS),ISXOUT(NSXDET+1),JOUTF(NSXDET+1)
     $,              ISXPW1(NSXDET+1),ISXPW2(NSXDET+1),IMXTAR
     $,              ISXTAR,ISXTIM,SXGATE(NSXDET)
     $,              ISXPLN(NPLANS),JPLNF(NPLANS)
     $,              ISXWKS,SXMAGN,ISXMAP,ISXVAC,ISXTUB,NSXEVT,ISXBG1
     $,              ISXBG2,IFPART,IFVOLU,IFMATE,IFTMED,IFVERT,IFKINE
     $,              IFSETS,IFHITS,IFDIGI,IFSECS
     $,              IFXSEC,IFLOSS
     $,              ISXEVT,ISXCHM
     $,              MAXNSTEP
     $,              ISXINP,ISXDIN,ISXCOOR
C
*KEEP,user common block
      INTEGER    NHIT,NDETNUM
      PARAMETER  (NDETNUM=10)
      REAL       DETPOS(NDETNUM,8),INITP(6)
      common/uTRACK/NHIT,Num_EVENT,DETPOS(NDETNUM,8),INITP(6),pl
C
*KEND.
C normal one
        REAL    XRES(10), YRES(10), TRES(10)
        DATA XRES /0.15, 0.15, 0.15, 0.15, 0.800,
     &             0.02, 0.02, 0.40, 0.40, 0.433/  ! cm
        DATA YRES /0.15, 0.15, 0.15, 0.15, 1.0, 
     &             0.40, 0.40, 0.80, 0.80, 1.0/    ! cm
        DATA TRES /0.09, 0.09, 0.09, 0.09, 0.09, 
     &             0.09, 0.09, 0.09, 0.09, 0.09/   ! nsec
C
      real      XSIGMA4(4),YSIGMA4(4),POSX4(4),POSY4(4),POSZ4(4)
      real      XSIGMA3(3),YSIGMA3(3),POSX3(3),POSY3(3),POSZ3(3)
      real      FISLPX(2),FISLPY(2),MOM(8),theta,phi
      real      FITX4(4),FITY4(4),FITX3(3),FITY3(3),CHI2X,CHI2Y
      integer   IDET,ISTAT,Iseed
      integer   i, j, k
      real      event(5),Kine(12)
*     ------------------------------------------------------------------
*
      IF(JUOUT.NE.0)THEN
         CALL CSJCAL(JUOUT,0,X,X,X,X,X,X,X,X,X,X)
         GO TO 99
      ENDIF
*
      CALL ICLRWK(0,0)
**
*
*           calcurate momentum
*
      if(Num_EVENT.ge.12001) goto 99
      Iseed = NRNDM(2)
      do IDET = 1,10
c         DETPOS(IDET,5)=DETPOS(IDET,1)+rangaus(0.,XRES(IDET),Iseed)
c         DETPOS(IDET,6)=DETPOS(IDET,2)+rangaus(0.,YRES(IDET),Iseed)
c         DETPOS(IDET,7)=DETPOS(IDET,3)
c         DETPOS(IDET,8)=DETPOS(IDET,4)+rangaus(0.,TRES(IDET),Iseed)
      enddo
C
      if(NHIT .eq. 4)then
         POSX4(1) = DETPOS(6,5)
         POSX4(2) = DETPOS(7,5)
         POSX4(3) = DETPOS(8,5)
         POSX4(4) = DETPOS(9,5)
C
         POSY4(1) = DETPOS(6,6)
         POSY4(2) = DETPOS(7,6)
         POSY4(3) = DETPOS(8,6)
         POSY4(4) = DETPOS(9,6)
C
         POSZ4(1) = DETPOS(6,7)
         POSZ4(2) = DETPOS(7,7)
         POSZ4(3) = DETPOS(8,7)
         POSZ4(4) = DETPOS(9,7)
C
         XSIGMA4(1)=XRES(6)
         XSIGMA4(2)=XRES(7)
         XSIGMA4(3)=XRES(8)
         XSIGMA4(4)=XRES(9)
C
         YSIGMA4(1)=YRES(6)
         YSIGMA4(2)=YRES(7)
         YSIGMA4(3)=YRES(8)
         YSIGMA4(4)=YRES(9)
C
C.. GET CHISQ (X-Z)
         CALL TRGTCH4(XSIGMA4,POSX4,POSZ4,FITX4,CHI2X,ISTAT,FISLPX)
c         write(6,*)'X = ',FISLPX(1),'*Z +',FISLPX(2),CHI2X,ISTAT
C.. GET CHISQ (Y-Z)
         CALL TRGTCH4(YSIGMA4,POSY4,POSZ4,FITY4,CHI2Y,ISTAT,FISLPY)
c         write(6,*)'Y = ',FISLPY(1),'*Z +',FISLPY(2),CHI2Y,ISTAT
cmk.. calcurate momentum and angle
         CALL MCMOM(FISLPX,FISLPY,MOM)
         theta = MOM(5)*180/3.14159265
         phi   = MOM(6)*180/3.14159265
         write(6,*)'P =',MOM(1),' theta =',theta,' phi =',phi
C
      elseif(NHIT .eq. 3)then
         if(DETPOS(6,5) .le. 80.)then
            i = 7
            j = 8
            k = 9
         elseif(DETPOS(7,5) .le. 80.)then
            i = 6
            j = 8
            k = 9
         elseif(DETPOS(8,5) .le. 80.)then
            i = 6
            j = 7
            k = 9
         elseif(DETPOS(9,5) .le. 80.)then
            i = 6
            j = 7
            k = 8
         else
            goto 200
         endif
         POSX3(1) = DETPOS(i,5)
         POSX3(2) = DETPOS(j,5)
         POSX3(3) = DETPOS(k,5)
C     
         POSY3(1) = DETPOS(i,6)
         POSY3(2) = DETPOS(j,6)
         POSY3(3) = DETPOS(k,6)
C
         POSZ3(1) = DETPOS(i,7)
         POSZ3(2) = DETPOS(j,7)
         POSZ3(3) = DETPOS(k,7)
C
         XSIGMA3(1)=XRES(i)
         XSIGMA3(2)=XRES(j)
         XSIGMA3(3)=XRES(k)
C
         YSIGMA3(1)=YRES(i)
         YSIGMA3(2)=YRES(j)
         YSIGMA3(3)=YRES(k)
C
C.. GET CHISQ (X-Z)
         CALL TRGTCH3(XSIGMA3,POSX3,POSZ3,FITX3,CHI2X,ISTAT,FISLPX)
c         write(6,*)'X = ',FISLPX(1),'*Z +',FISLPX(2),CHI2X,ISTAT
C.. GET CHISQ (Y-Z)
         CALL TRGTCH3(YSIGMA3,POSY3,POSZ3,FITY3,CHI2Y,ISTAT,FISLPY)
c         write(6,*)'Y = ',FISLPY(1),'*Z +',FISLPY(2),CHI2Y,ISTAT
cmk.. calcurate momentum and angle
         CALL MCMOM(FISLPX,FISLPY,MOM)
         theta = MOM(5)*180/3.14159265
         phi   = MOM(6)*180/3.14159265
         write(6,*)'P =',MOM(1),' theta =',theta,' phi =',phi
      else
 200     continue
         write(6,*)'Do not calcurate momentum'
      endif
**
*
*             Filling Histogram file
*
      if((NHIT .eq. 4).or.(NHIT .eq.3))then
         Num_EVENT = Num_EVENT + 1
         event(1) = IEVENT
         event(2) = NHIT
         event(3) = MOM(1)
         event(4) = theta
         event(5) = phi
         CALL HFN(10,event)
         kine(1) = MOM(1)
         Kine(2) = MOM(2)
         Kine(3) = MOM(3)
         Kine(4) = MOM(4)
         Kine(5) = MOM(5)
         Kine(6) = MOM(6)
         Kine(7) = INITP(1)
         Kine(8) = INITP(2)
         Kine(9) = INITP(3)
         Kine(10)= INITP(4)
         Kine(11)= INITP(5)
         Kine(12)= INITP(6)
         CALL HFN(50,kine)
      endif
**
*
*             Draw track of particle and detector
*
*
      if( ISWIT(5) .eq. 0 )then
         if( ISWIT(2) .eq. 2 )then
            CALL GDRAW('WA98',90.,180.,0.,17.,10.,0.008,0.008)
         elseif( ISWIT(2) .eq. 3)then
            CALL GDRAW('WA98',90.,90.,0.,17.,10.,0.008,0.008)
         elseif( ISWIT(2) .eq. 4)then
            CALL GDRAW('WA98',40.,60.,0.,6.,5.,0.012,0.012)
         elseif( ISWIT(2) .eq. 5)then
            CALL GDRAW('WA98',150.,90.,0.,10.,15.,0.012,0.012)
         elseif( ISWIT(2) .eq. 5)then
            CALL GDRAW('WA98',30.,90.,0.,10.,5.,0.012,0.012)
         else
            CALL GDRAW('WA98',140.,120.,0.,14.,15.,0.012,0.012)
         endif
*
         CALL GDHEAD(101,'jTOF second arm$',.5)
         CALL GDXYZ(0)
         CALL GDAXIS(XTARG,YTARG,ZTARG,100.)
         CALL GDPART(0,11,0.25)
         CALL GDMAN(18.,4.0)
         CALL GDSCAL(10.,2.)
      endif
*
      RETURN
 99      END
*
*
**
*
      subroutine pion_gen(y_mean,y_sigma,y,p4,Iseed)
*
************************************************************************
*  this routine generates pions with constant rapidity density
*  and with the slope of 150 MeV/c.  pion0 = 139 MeV
*     y_mean  ; INPUT  : Mean value of pion rapidity
*     y_sigma ; INPUT  : Sigma of pion rapidity
*     y       ; OUTPUT : pion rapidity
*     p4      ; OUTPUT : Four Momentum of the particles
*     Iseed   ; IN/OUT : INTEGER*4  Seed of random number
************************************************************************
*
      real y, y_mean, y_sigma, mass, pt, p4(4), slope, sita, arg
C
c      y = min(y_min,y_max) + abs(y_max-y_min)*ran(Iseed)  ! rapidity
      y =  rangaus(y_mean,y_sigma,Iseed) ! rapidity
C
      slope = 1./6.             ! 1/6 [GeV/c] ==> 150 [MeV/c]
      pt = ran_pt(slope, Iseed)
C
      arg = 2.*3.1415926*ran(Iseed)
C
      p4(1) = pt*cos(arg)
      p4(2) = pt*sin(arg)
      mass = 0.139
      p4(3) = sqrt(pt**2 + mass**2)*(exp(y) - exp(-y))/2. ! p_para
      p4(4) = sqrt(mass**2 + p4(1)**2 + p4(2)**2 + p4(3)**2)
C
      return
      end
*
***
*
      real function rangaus(Mean,Sigma,Iseed)
*
************************************************************************
*  Random number generation of Gauss distr.
*       Amean ; INPUT  : Mean value of Gaussian
*       Sigma ; INPUT  : Sigma      of Gaussian
*       Iseed ; INPUT  : INTEGER*4  Seed of random number
************************************************************************
*
      real Mean, Sigma, arg
C
      arg  = Sigma * sqrt(abs(2.*log(ran(Iseed))))
      rangaus = arg*cos(6.283185307*ran(Iseed)) + Mean
*
      RETURN
      END
*
***
*
      real function ranexp(Slope,Iseed)
*
************************************************************************
*  Random number generation of Exponential distribution.
*       Slope ; INPUT  : Slope of exponentialdistribution
*       Iseed ; INPUT  : INTEGER*4  Seed of random number
************************************************************************
*
      ranexp = -Slope * Alog(1.-RAN(Iseed))
*
      RETURN
      END
*
**
*
      real function ran_pt(Slope,Iseed)
*
************************************************************************
*  Generate pt when the invariant cross section follows
*  Exponential distribution.
*       Slope ; INPUT  : Slope of exponential distribution
*       Iseed ; INPUT  : INTEGER*4  Seed of random number
************************************************************************
*       dN
*      ---- == Pt*exp(-Pt/S)   ;  S = slope
*       dPt
************************************************************************
*
      common/com_ran_pt/prob, pt_slope
      external pt_func
      real prob, pt_slope
C
      pt_slope = Slope
      prob = ran(Iseed)
C
      CALL RZERO(0., 5., ran_pt, R, 0.01, 100, pt_func)
      return
      end
**
      real function pt_func(x)
      common/com_ran_pt/prob, pt_slope
C
      pt_func = 1.- (1.+ x/pt_slope)*exp(-x/pt_slope) - prob
      return
      end
*
*
************************************************************************
      SUBROUTINE   TRGTCH4(SIGPS,POSC1,POSC2,POSF1,CHI2,ISTAT6,FISLP)
cmk                       (sigma,px(py),pz  ,fitX(Y),            )
C       GET Chi2 & positions.
************************************************************************

      IMPLICIT        NONE
c.. define DETPARAMS.
      INTEGER         MAXHODS
      PARAMETER      (MAXHODS=4)
c.. define subroutine parameters
      INTEGER         ISTAT6
      INTEGER         IHOD, IP1
      REAL            SIGPS(MAXHODS)
      REAL            POSC1(MAXHODS), POSC2(MAXHODS)
      REAL            POSF1(MAXHODS)
      REAL            CHI2
c.. define local variables
      REAL            SUM(5), PAR(2), NORM, D
cmk+
      REAL            FISLP(2)
cmk+
C--------1---------2---------3---------4---------5---------6---------7--
c.. initialize
      ISTAT6=0
      DO IP1=1,5
         SUM(IP1)=0.00
      ENDDO

      DO IHOD=1, MAXHODS
         NORM=SIGPS(IHOD)*SIGPS(IHOD)
         SUM(1)=SUM(1)+1.00/NORM
         SUM(2)=SUM(2)+POSC2(IHOD)/NORM
         SUM(3)=SUM(3)+POSC1(IHOD)/NORM
         SUM(4)=SUM(4)+POSC2(IHOD)*POSC2(IHOD)/NORM
         SUM(5)=SUM(5)+POSC1(IHOD)*POSC2(IHOD)/NORM
       ENDDO
       D=SUM(1)*SUM(4)-SUM(2)*SUM(2)
       IF(D.NE.0.) THEN
          PAR(1)=(SUM(1)*SUM(5)-SUM(2)*SUM(3))/D
          PAR(2)=(SUM(3)*SUM(4)-SUM(2)*SUM(5))/D
cmk linear equation x=par(1)*z+par(2)
       ELSE
          ISTAT6=2
          WRITE(6,1400)
 1400     FORMAT(1H ,' % TRGTCH2: ERROR. DEVIDED BY ZERO(D=0) IN
     &         sub.TRGTCH2.')
          RETURN
       ENDIF
       
        CHI2=0.00
        DO IHOD=1,MAXHODS
           POSF1(IHOD)=PAR(1)*POSC2(IHOD)+PAR(2)
           CHI2=CHI2+(POSC1(IHOD)-POSF1(IHOD))/SIGPS(IHOD)*
     &              (POSC1(IHOD)-POSF1(IHOD))/SIGPS(IHOD)
        ENDDO
cmk+
        FISLP(1)=PAR(1)
        FISLP(2)=PAR(2)
cmk+
        RETURN
        END
*
*
************************************************************************
      SUBROUTINE   TRGTCH3(SIGPS,POSC1,POSC2,POSF1,CHI2,ISTAT6,FISLP)
cmk                       (sigma,px(py),pz  ,fitX(Y),            )
C       GET Chi2 & positions.
************************************************************************

      IMPLICIT        NONE
c.. define DETPARAMS.
      INTEGER         MAXHODS
      PARAMETER      (MAXHODS=3)
c.. define subroutine parameters
      INTEGER         ISTAT6
      INTEGER         IHOD, IP1
      REAL            SIGPS(MAXHODS)
      REAL            POSC1(MAXHODS), POSC2(MAXHODS)
      REAL            POSF1(MAXHODS)
      REAL            CHI2
c.. define local variables
      REAL            SUM(5), PAR(2), NORM, D
cmk+
      REAL            FISLP(2)
cmk+
C--------1---------2---------3---------4---------5---------6---------7--
c.. initialize
      ISTAT6=0
      DO IP1=1,5
         SUM(IP1)=0.00
      ENDDO

      DO IHOD=1, MAXHODS
         NORM=SIGPS(IHOD)*SIGPS(IHOD)
         SUM(1)=SUM(1)+1.00/NORM
         SUM(2)=SUM(2)+POSC2(IHOD)/NORM
         SUM(3)=SUM(3)+POSC1(IHOD)/NORM
         SUM(4)=SUM(4)+POSC2(IHOD)*POSC2(IHOD)/NORM
         SUM(5)=SUM(5)+POSC1(IHOD)*POSC2(IHOD)/NORM
       ENDDO
       D=SUM(1)*SUM(4)-SUM(2)*SUM(2)
       IF(D.NE.0.) THEN
          PAR(1)=(SUM(1)*SUM(5)-SUM(2)*SUM(3))/D
          PAR(2)=(SUM(3)*SUM(4)-SUM(2)*SUM(5))/D
cmk linear equation x=par(1)*z+par(2)
       ELSE
          ISTAT6=2
          WRITE(6,1400)
 1400     FORMAT(1H ,' % TRGTCH2: ERROR. DEVIDED BY ZERO(D=0) IN
     &         sub.TRGTCH2.')
          RETURN
       ENDIF
       
        CHI2=0.00
        DO IHOD=1,MAXHODS
           POSF1(IHOD)=PAR(1)*POSC2(IHOD)+PAR(2)
           CHI2=CHI2+(POSC1(IHOD)-POSF1(IHOD))/SIGPS(IHOD)*
     &              (POSC1(IHOD)-POSF1(IHOD))/SIGPS(IHOD)
        ENDDO
cmk+
        FISLP(1)=PAR(1)
        FISLP(2)=PAR(2)
cmk+
        RETURN
        END
*
*
c#########################################################################
c     Calcurate a momentun from a tracking information      '95/9/22
c     written by     Mizuki Kurata
c#########################################################################
c
c*********-*********-*********-*********-*********-*********-*********-****
      SUBROUTINE MCMOM(SLPX,SLPY,MOM)
c*********-*********-*********-*********-*********-*********-*********-***
c
c     input :X=SLPX(1)*Z+SLPX(2), Y=SLPY(1)*Z+SLPY(2)
c     output:MOM 1;P,2;Px,3;Py,4;Pz,5;theta,6;phi
c*********-*********-*********-*********-*********-*********-*********-****
c
*KEEP,user common block
      COMMON/MKMAG/RMKF
*KEND.
      REAL SLPX(2),SLPY(2),MOM(6),TOTP
      REAL SLPX1S,SLPX1S1,SLPX2S,SLPX12,R0,R0S,POSP
      REAL OGZ,OGX,OGY,CRZ,CRX,HCRZ,HCRX,HCR_SLP
      REAL OGPRB,HB
      REAL PCZ,PCX,CDZ,CDX
      REAL PROU
      REAL MAG,MOMXZ
      REAL THXZ1,THXZ2,THXZ3,THOUT,DTHOUT
      REAL DST_CD,DST_CB,FAN,DST_TM,PLGXZ
      REAL THYZ1
      PARAMETER (R0=328.)
c      PARAMETER (R0=200.)
c
c---------*---------*---------*---------*---------*---------*---------*----
c
C     FOR ON ZX PLANE
c     definition of line and cross point
c     line #1;X=SLPX(1)*Z+SLPX(2)
c     line #2;X=-1./SLPX(1)*Z+OGPRB
c     line #3;X=HCRX/HCRZ*Z (HCRX/HCRZ=HCR_SLP)
c     line #4;X=-1./HCR_SLP*Z+HB
c
c     point A;cross of line #1 and circle with radius of R0 and line #4 
c     point B;cross of line #3 and line #4
c     point C;cross of line #2 and line #3
c     point D;cross of line #2 and line #4
c     point E;cross of line #1 and magnet edge
c
c     useful parameter
      SLPX1S =SLPX(1)**2.
      SLPX1S1=1+SLPX1S
      SLPX2S =SLPX(2)**2.
      SLPX12 =SLPX(1)*SLPX(2)
      R0S =R0**2.
c
c starting point;
c     line #1 X=SLPX(1)*Z+SLPX(2) is given by tracking information
c
c step 1;cross point E between line #1 and the magnet edge.
c (OGZ,OGX)
      POSP=-SLPX2S+10000*SLPX1S1
      IF(POSP.GT.0) THEN
         OGZ=(-SLPX12+SQRT(POSP))/SLPX1S1
         OGX=SLPX(1)*OGZ+SLPX(2)
cdb
c         WRITE(6,*) 'X =',OGX,' and ','Z =',OGZ
c         WRITE(6,*) 'Function is   X=',SLPX(1),'Z+',SLPX(2)
cdb
      ELSE
*         WRITE(6,*) 'WARNING SQR MK',POSP
      ENDIF
c
c step 2;getting line #2
      OGPRB=OGX+1./SLPX(1)*OGZ
c     line #2 X=-1./SLPX(1)*Z+OGPRB is available
cdbg
c      TYPE*,'line2;X=',-1./SLPX(1),'*Z+',OGPRB
cdbg
c
c step 3;cross point A between x**2+y**2=R0**2 and line #1
c (CRZ,CRX)
      CRZ=(-SLPX12+SQRT(-SLPX2S+R0S*SLPX1S1))/SLPX1S1
      CRX=SLPX(1)*CRZ+SLPX(2)
c 
c step 4;getting half point B between point A and the traget.
c        and line #3
c (HCRZ,HCRX)
      HCRZ=(CRZ-R0)/2.
      HCRX=CRX/2.
      HCR_SLP=HCRX/HCRZ
c     line #3 X=HCRX/HCRZ*Z (HCRX/HCRZ=HCR_SLP) is available
cdbg
c      TYPE*,'line3;X=',HCR_SLP,'*Z+'
cdbg
c
c step 5; getting line #4
      HB=CRX+CRZ/HCR_SLP
c     line #4 X=-1./HCR_SLP*Z+HB is available
cdbg
c      TYPE*,'line4;X=',-1./HCR_SLP,'*Z+',HB
cdbg
cdb
c      write(6,*) 'cross point z=',crz
c      write(6,*) 'cross point x=',crx
c      write(6,*) 'half point  z=',hcrz      
c      write(6,*) 'half point  x=',hcrx
c      write(6,*) 'half slope   =',hcr_slp
cdb
c
c step 6;getting cross point C between line #2 and line #3
c        (PCZ,PCX)
      PCZ=OGPRB/(HCR_SLP+1./SLPX(1))
      PCX=HCR_SLP*PCZ
c
c step 7;getting cross point D between line #2 and line #4
c        (CDZ,CDX)
      CDZ=(-OGPRB+HB)/(1./HCR_SLP-1./SLPX(1))
      CDX=-CDZ/SLPX(1)+OGPRB
c
c** for Momentum calcuration ******
c   
c     radius of curvature on XZ plane
      PROU=(PCZ-OGZ)**2.+(PCX-OGX)**2.
      PROU=SQRT(PROU)
c
c     momentum
      MAG = RMKF
      MOMXZ=3.*PROU*MAG/1000.       !Pxz
cdbg{
c      TYPE*,'cross point C',PCZ,PCX
c      TYPE*,'cross point E',OGZ,OGX
c      TYPE*,'PROU',PROU
c      TYPE*,'PXZ',MOMXZ
cdbg}
cdbg{
c      MOM(1)=MOMXZ
cdbg}
c
c*** for out-going angle on xz plane **
      THXZ1=ATAN(SLPX(1))          !angle of line #1
      THXZ2=ATAN(-1./HCR_SLP)      !angle of line #4
      THOUT=THXZ1-2.*THXZ2         !out-going angle in rad.
      DTHOUT=THOUT*180./3.141592   !out-going angle in deg.
cdbg{
c      TYPE*,'THXZ1 ON XZ',THXZ1*180./3.141592
c      TYPE*,'THXZ2 ON XZ',THXZ2*180./3.141592
c      TYPE*,'THOUT ON XZ',DTHOUT
cdbg}
c
c**   for Px and Pz
      MOM(2)=-MOMXZ*SIN(THOUT) !Px
      MOM(4)= MOMXZ*COS(THOUT) !Pz
c**   for Py
c     path integral in magnetic field
      DST_CD=SQRT((PCZ-CDZ)**2.+(PCX-CDX)**2.)
      DST_CB=SQRT((PCZ-HCRZ)**2.+(PCX-HCRX)**2.)
      THXZ3=ACOS(DST_CB/DST_CD)
c
      FAN=PROU*2.*THXZ3  ! fan shaped
c     path integral a straght line from target to the magnet
c     which is jast same as the distance A and E      
      DST_TM=SQRT((CRZ-OGZ)**2.+(CRX-OGX)**2.)
c     totaly
      PLGXZ=FAN+DST_TM
c    
c     y-position
      OGY=SLPY(1)*OGZ+SLPY(2)
c     out-going angle on YZ plane
      THYZ1=ATAN(OGY/PLGXZ)
c
c     momentum P
      TOTP=MOMXZ/COS(THYZ1)    !p
      IF(HB.GT.SLPX(2)) THEN
         MOM(1)=TOTP           !positive charged
      ELSEIF(HB.LT.SLPX(2)) THEN
         MOM(1)=-TOTP          !negative charged
      ELSE
         MOM(1)=0.
         WRITE(6,*) 'Neutral ??'
      ENDIF
      MOM(3)=MOMXZ*OGY/PLGXZ   !Py
c
c*** for angle of theta and phi
      MOM(5)=ACOS(MOM(4)/TOTP) ! theta
      IF(MOM(2).GT.0) THEN
         MOM(6)=ATAN(MOM(3)/MOM(2)) ! phi
      ELSEIF(MOM(2).LT.0) THEN
         MOM(6)=(3.141592/2.+ATAN(ABS(MOM(2)/MOM(3))))
     +        *(MOM(3)/ABS(MOM(3)))
      ELSE
         MOM(6)=3.141592/2.*(MOM(3)/ABS(MOM(3)))
      ENDIF
c
c
      RETURN
      END
c
c   f_gen_GAUSS.for
c
c     generate gaussian distribution
c         
c
      REAL FUNCTION F_GEN_GAUSS()
c
      implicit none
c
      real*4 hrndm1, dummy
      integer*4 idms
      logical first
      external f_GAUSS
c
      data idms/21003/
      data first/.true./
c
      if(first) then            ! prepare function
         first = .false.
         call HBFUN1(idms,'  Gaussian Distribution $',
     $        200, -5., +5., f_GAUSS)
      endif
c
      f_gen_GAUSS = HRNDM1(idms,dummy)
c
      return
      end
c
      function f_GAUSS(x)
      f_GAUSS = exp(-0.5*x**2)
      return
      end
