***
*     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
*KEND.
cn..          Define of HBOOK word
      integer   istat
      character chtags(7)*5
      data      chtags/'I','Npad1','Npad2','Nstr1','Nstr2','Npmd','Nx'/
      character ckovtags(4)*5
      data      ckovtags/'I','N','x','y'/
      character emtags(6)*5
      data      emtags/'I','E','x','y','r','sita'/
      character pitags(6)*5
      data      pitags/'I','E1','E2','Pt','phi','M'/
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
      IPFIT=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=1 ! 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
**
*
*
*             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,'PHENIX','phenix.nt','n',1024,istat)
**
*
*             Booking Histogram & Ntuple file
*
      CALL HBOOKN(10,'Number of gamma',7,'PHENIX',5000,chtags)
c      CALL HBOOKN(20,'Cherenkov Number',4,'PHENIX',5000,ckovtags)
c      CALL HBOOKN(30,'Energy of gamma',6,'PHENIX',5000,emtags)
c      CALL HBOOKN(40,'Energy of pion0',6,'PHENIX',5000,pitags)
C
      CALL HBOOK1(100,'Pt-distribution(MeV/c)',100,0.,2000.,0.)
      CALL HBOOK1(110,'dN/(Pt*dPt)(MeV/c)',100,0.,2000.,0.)
      CALL HBOOK1(120,'rapidity-distribution',20,-2.,8.,0.)
C
      CALL HBOOK1(210,'PAD1 flame (MeV)',100,0.,2000.,0.)
      CALL HBOOK1(220,'PAD2 flame (MeV)',100,0.,2000.,0.)
      CALL HBOOK1(230,'STR1 flame (MeV)',100,0.,2000.,0.)
      CALL HBOOK1(240,'STR2 flame (MeV)',100,0.,2000.,0.)
C
      CALL HBOOK1(300,'ratio of hitting 2nd arm (%)',200,0.,40.,0.)
*
      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('PHENIX')
**
*
      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/uNUMBER/Npad1,Npad2,Nstr1,Nstr2,Narm,Npmd,Nx,
     +     IT_pad1(2000),IS_pad1(2000),IT_pad2(2000),IS_pad2(2000),
     +     IT_str1(2000),IS_str1(2000),IT_str2(2000),IS_str2(2000),
     +     IT_arm(2000),IS_arm(2000),IT_pmd(2000),IS_pmd(2000)
C
*KEND.
C
      REAL      UBUF,PLAB(3),VERT(3)
      INTEGER   NVTX,NT,IPART
C
      real      pt, p4(4), Init_pt,Init_num, y_mean, y_sigma, y
      integer   Iseed, Init_PART,In_part
      real      Nckov,Nck(96,48)
      integer   inum,i
      character eventname*20
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_pt,Iptswit
*
      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,*)' '
c         WRITE(6,*)' Init_PART = 7 : pion0 '
         READ(5,*)Init_PART
*
         WRITE(6,*)' '
         WRITE(6,*)'Number of paeticles  ..integer'
c         WRITE(6,*)'  dN/dy = 375'
         READ(5,*)Init_num
      ENDIF
*
      VERT(1) = XTARG
      VERT(2) = YTARG
      VERT(3) = ZTARG
      CALL GSVERT(VERT,0,0,UBUF,1,NVTX)
*
      y_mean  = 2.90
      y_sigma = 1.31
c      Init_PART = 7
c      Init_num = 375*6          ! dN/dy * (y_max - y_min)
C
      Iptswit = 0
      Iseed = NRNDM(2)
      do inum = 1,Init_num
         call pion_gen(y_mean,y_sigma,y,p4,Iseed)
         pt = sqrt(p4(1)**2 + p4(2)**2)
*
         PLAB(1) = p4(1)
         PLAB(2) = p4(2)
         PLAB(3) = p4(3)
         CALL GSKINE(PLAB,Init_PART,NVTX,UBUF,1,NT)
*
c         WRITE(6,*) IEVENT,Init_PART,pt
         CALL HFILL(100,pt*1.0E+3,0.,1.) ! Fill of Pt-distribution
         CALL HFILL(110,pt*1.0E+3,0.,1./pt) ! Fill of Pt-distribution
         CALL HFILL(120,y,0.,1.) ! Fill of rapidity-distribution
      enddo
*
      Npad1 = 0
      Npad2 = 0
      Nstr1 = 0
      Nstr2 = 0
      Narm  = 0
      Npmd  = 0
      Nx    = 0
      do i = 1,2000
         IT_pad1(i) = 0
         IS_pad1(i) = 0
         IT_pad2(i) = 0
         IS_pad2(i) = 0
         IT_str1(i) = 0
         IS_str1(i) = 0
         IT_str2(i) = 0
         IS_str2(i) = 0
         IT_arm(i)  = 0
         IS_arm(i)  = 0
         IT_pmd(i)  = 0
         IS_pmd(i)  = 0
      enddo
*
*
*
*             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
C
      DIMENSION X(3),F(3)
C.
C.    ------------------------------------------------------------------
*
      F(1)=0.
      F(2)=0.
      F(3)=0.
C
      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,user common block
      common/uNUMBER/Npad1,Npad2,Nstr1,Nstr2,Narm,Npmd,Nx,
     +     IT_pad1(2000),IS_pad1(2000),IT_pad2(2000),IS_pad2(2000),
     +     IT_str1(2000),IS_str1(2000),IT_str2(2000),IS_str2(2000),
     +     IT_arm(2000),IS_arm(2000),IT_pmd(2000),IS_pmd(2000)
C
*KEND.
C
      integer   Npad1,Npad2,Nstr1,Nstr2,Npmd,Nx
     +     IT_pad1(2000),IS_pad1(2000),IT_pad2(2000),IS_pad2(2000),
     +     IT_str1(2000),IS_str1(2000),IT_str2(2000),IS_str2(2000),
     +     IT_arm(2000),IS_arm(2000),IT_pmd(2000),IS_pmd(2000)
*     -----------------------------------------------------------------
*
      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((INWVOL.eq.1).and.(IPART.eq.1).and.(GETOT .ge. 0.005))then
         if((NAMES(3) .eq. 'JPAD').and.(NAMES(4) .ne. 'JPDW')
     +        .and.(NUMBER(3).eq.1))then
            CALL HFILL(210,GETOT*1.0E3,0.,1.)
            Npad1 = Npad1 + 1
            IT_pad1(Npad1) = ITRA
            IS_pad1(Npad1) = ISTAK
c            write(6,*)'Npad1=',Npad1,ITRA,ISTAK
         elseif((NAMES(3) .eq. 'JPAD').and.(NAMES(4) .ne. 'JPDW')
     +           .and.(NUMBER(3).eq.2))then
            CALL HFILL(220,GETOT*1.0E3,0.,1.)
            Npad2 = Npad2 + 1
            IT_pad2(Npad2) = ITRA
            IS_pad2(Npad2) = ISTAK
c            write(6,*)'Npad2=',Npad2,ITRA,ISTAK
         elseif((NAMES(3) .eq. 'JSTR').and.(NAMES(4) .ne. 'JSTW')
     +           .and.(NUMBER(3).eq.1))then
            CALL HFILL(230,GETOT*1.0E3,0.,1.)
            Nstr1 = Nstr1 + 1
            IT_str1(Nstr1) = ITRA
            IS_str1(Nstr1) = ISTAK
c            write(6,*)'Nstr1=',Nstr1,ITRA,ISTAK
         elseif((NAMES(3) .eq. 'JSTR').and.(NAMES(4) .ne. 'JSTW')
     +           .and.(NUMBER(3).eq.2))then
            CALL HFILL(240,GETOT*1.0E3,0.,1.)
            Nstr2 = Nstr2 + 1
            IT_str2(Nstr2) = ITRA
            IS_str2(Nstr2) = ISTAK
c            write(6,*)'Nstr2=',Nstr2,ITRA,ISTAK
         endif
         if((NAMES(2) .eq. 'PMD ').and.(NAMES(3) .ne. 'PHOL'))then
            Npmd = Npmd + 1
            IT_pmd(Npmd) = ITRA
            IS_pmd(Npmd) = ISTAK
c            write(6,*)'Npmd=',Npmd,ITRA,ISTAK
         endif
      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,user common block
      common/uNUMBER/Npad1,Npad2,Nstr1,Nstr2,Narm,Npmd,Nx,
     +     IT_pad1(2000),IS_pad1(2000),IT_pad2(2000),IS_pad2(2000),
     +     IT_str1(2000),IS_str1(2000),IT_str2(2000),IS_str2(2000),
     +     IT_arm(2000),IS_arm(2000),IT_pmd(2000),IS_pmd(2000)
C
*KEND.
C
      real      event(7), ratio
      integer   i, j, Mpmd(2000)
c      parameter (maxN = 1000)
c      real      Nem(maxN),E_em(maxN),Rem(maxN),sita_em(maxN),Rx,Ry
c      integer   xem(maxN),yem(maxN)
c      real      E_pi,Pt_pi,phi_pi,M_pi
*     ------------------------------------------------------------------
*
      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)
**
*
*             Filling Histogram file
*
      do i = 1,2000
         Mpmd(i) = 0
      enddo
      do i = 1,Npmd
         do j = 1,Npad1
            if((IT_pmd(i).eq.IT_pad1(j)).and.
     +           (IS_pmd(i).eq.IS_pad1(j)))then
               Mpmd(i) = Mpmd(i) + 1
c               write(6,*)Mpmd(i),IT_pmd(i),IS_pmd(i)
            endif
         enddo
         do j = 1,Npad2
            if((IT_pmd(i).eq.IT_pad2(j)).and.
     +           (IS_pmd(i).eq.IS_pad2(j)))then
               Mpmd(i) = Mpmd(i) + 1
c               write(6,*)Mpmd(i),IT_pmd(i),IS_pmd(i)
            endif
         enddo
         do j = 1,Nstr1
            if((IT_pmd(i).eq.IT_str1(j)).and.
     +           (IS_pmd(i).eq.IS_str1(j)))then
               Mpmd(i) = Mpmd(i) + 1
c               write(6,*)Mpmd(i),IT_pmd(i),IS_pmd(i)
            endif
         enddo
         do j = 1,Nstr2
            if((IT_pmd(i).eq.IT_str2(j)).and.
     +           (IS_pmd(i).eq.IS_str2(j)))then
               Mpmd(i) = Mpmd(i) + 1
c               write(6,*)Mpmd(i),IT_pmd(i),IS_pmd(i)
            endif
         enddo
         if(Mpmd(i) .ne. 0)then
            Nx = Nx + 1
         endif
      enddo
C
      ratio = float(Nx)/float(Npmd)*100.
      write(6,*)IEVENT,Nx,Npmd,ratio
      CALL HFILL(300,ratio,0.,1.)
      event(1) = IEVENT
      event(2) = Npad1
      event(3) = Npad2
      event(4) = Nstr1
      event(5) = Nstr2
      event(6) = Npmd
      event(7) = Nx
      CALL HFN(10,event)
**
*
*             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
*
*
