***
*     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 toftags(9)*5
      data      toftags/'I','IDtof','IP','TOF','E','Eloss',
     &     'X','Y','EMcal'/
      character emtags(6)*5
      data      emtags/'I','E','x','y','r','sita'/
C
***   -----------------------------------------------------------------
*
*
*             Initialize GEANT
*
*
      CALL GINIT           ! Initialize GEANT system
*
*
*             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=1 ! pair production with generation of electron/positron
      ICOMP=1 ! compton scattering with generation of electron
      IPHOT=1 ! photo-electric effect with generation of electron
      IPFIS=0 ! photo-fission effect with generation of secondaries
      IDRAY=1 ! delta ray production with generation of electron
      IANNI=1 ! positron annihilation with generation of photons
      IBREM=1 ! bremsstrahlung with generation of gamma
      IHADR=1 ! hadronic interactions with generation of secondaries
      IMUNU=1 ! muon nuclear interaction with generation of secondaries
      IDCAY=1 ! decay in flight with generation of secondaries
      ILOSS=2 ! continuous energy loss without generation of delta ray 
C                    and Landau-Vaviliov-Gauss fluctuations
      IMULS=1 ! multiple scattering according to Moli'ere theory
      IRAYL=0 ! Rayleigh scattering
      ILABS=1 ! absorbtion of Cerenkov photons with posssible detection
      ISYNC=0 ! To control synchrotron radiation
      ISTRA=0 ! To control energy loss fluctuation model
**
      write(6,*)' '
      write(6,*)'      Type in !'
      write(6,*)' '
      write(6,*)' Select of CUT parameter'
      write(6,*)'    cut for gammas and electron (CUTGAM,CUTELE)[GeV]'
      read(5,*) CUTGAM
      CUTELE = CUTGAM 
      write(6,*)'    cut for hadrons  (CUTHAD,CUTNEU,CUTMUO)[GeV]'
      read(5,*) CUTHAD
      CUTNEU = CUTHAD
      CUTMUO = CUTHAD
*
      write(6,*)' '
      write(6,*)'Select of GAS '
      write(6,*)'      Vaccume ===>0'
      write(6,*)'      Air     ===>1'
      write(6,*)'      He(gas) ===>2'
      read(5,*) ISWIT(6)
*
      write(6,*)' '
      write(6,*)'Select of ToF status '
      write(6,*)'    non-panel ===>0'
      write(6,*)'      1-panel ===>1'
      write(6,*)'      8-panel ===>2'
      write(6,*)'    non-panel &'
      write(6,*)'    EMcal 50cm===>3'
c      read(5,*) ISWIT(7)
      ISWIT(7) = 1
*
c      BCUTE=CUTGAM
c      BCUTM=CUTGAM
c      DCUTE=CUTELE
c      DCUTM=CUTELE
*
*
*             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
      CALL GMATE           ! Initialize material constants
*
      CALL UGEOM           ! Define user geometry
      CALL GPHYSI          ! Initialize GEANT physics processes
*
c      CALL GPTMED          ! Prints the tracking medium parameters
*
      write(6,*)' ---------- USER Selection ----------'
      if(ISWIT(6) .eq. 1)then
         write(6,*)' Air'
      elseif(ISWIT(6) .eq. 2)then
         write(6,*)' Helium gas'
      else
         write(6,*) ' Vaccume'
      endif
      if(ISWIT(7) .eq. 0)then
         write(6,*)' This geom. is EMcal only'
      elseif(ISWIT(7) .eq. 1)then
         write(6,*)' This geom. is ToF 1-panel and EMcal small'
      elseif(ISWIT(7) .eq. 3)then
         write(6,*)' This geom. is EMcal only (near)'
      else
         write(6,*)' This geom. is ToF 8-panel and EMcal '
         write(6,*)'      Full Scale !!'
      endif
      write(6,*)' '
      write(6,*)' ------------------------------------'
      write(6,*)' '
*
**
*
*             Open HBOOK
*
      CALL HROPEN(1,'PHENIX','phenix.nt','n',1024,istat)
**
*
*             Booking Histogram & Ntuple file
*
      CALL HBOOKN(10,'ToF Alvedo full',9,'PHENIX',5000,toftags)
      CALL HBOOKN(20,'ToF Alvedo full (slat)',9,'PHENIX',5000,toftags)
      CALL HBOOKN(30,'ToF Alvedo eloss',9,'PHENIX',5000,toftags)
      CALL HBOOKN(40,'ToF Alvedo eloss (slat)',9,'PHENIX',5000,toftags)
C
      CALL HBOOK1(300,'EMcal',100,0.,1600,0.)
      CALL HBOOK1(310,'EMcal',200,0.,1600,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,user common block
      common/uELOSS/In_part,IDcal,pt,Init_ENERGY,ELOSStof,ELOSS(8,96),
     +              Nckov,Nck(96,48)
      common/uENERGY/IP(96),TOFcal(96),ENEcal(96),xcal(96),ycal(96)
C
*KEND.
C
      REAL      UBUF,PLAB(3),VERT(3)
      INTEGER   NVTX,NT,IPART
C
      real      pt, p4(4), Init_pt, y_min, y_max, y
      integer   Iseed, Init_PART,In_part,Iptswit, IDtof, IDcal
      real      Nckov,Nck(96,48)
      integer   inum,ick, jck
      character eventname*20
      real      theta, phi, p_theta, p_phi
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.   '
         READ(5,*)Init_PART
         In_part = Init_PART
*
         WRITE(6,*)' '
         WRITE(6,*)'Momentum of paeticles  ..real'
         READ(5,*)Init_pt
*
      ENDIF
*
      Iseed = NRNDM(2)
*
      VERT(1) = 0.
      VERT(2) = 0.
      VERT(3) = 0.
c      VERT(2) = -1. + 2.*run(Iseed)
c      VERT(3) = -1. + 2.*run(Iseed)
      CALL GSVERT(VERT,0,0,UBUF,1,NVTX)
      pt = Init_pt
      PLAB(1) = pt
      PLAB(2) = 0.
      PLAB(3) = 0.
      CALL GSKINE(PLAB,Init_PART,NVTX,UBUF,1,NT)
*
      ELOSStof = 0.
      Nckov = 0.
      do ick = 1,96
         do jck = 1,48
            Nck(ick,jck) = 0.
            if(jck .le. 8)then
               ELOSS(jck,ick) = 0.
            endif
         enddo
      enddo
C
      Ncal = 0
      Ntof = 0
      do i = 1,96
         IP(i) = 0
         TOFcal(i) = 0.
         ENEcal(i) = 0.
         xcal(i) = 0.
         ycal(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,NAMEC1
      PARAMETER (MAXME1=30)
      COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1)
*KEEP,user common block
      common/uELOSS/In_part,IDcal,pt,Init_ENERGY,ELOSStof,ELOSS(8,96),
     +              Nckov,Nck(96,48)
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/uELOSS/In_part,IDcal,pt,Init_ENERGY,ELOSStof,ELOSS(8,96),
     +              Nckov,Nck(96,48)
      common/uENERGY/IP(96),TOFcal(96),ENEcal(96),xcal(96),ycal(96)
C
*KEND.
C
      real      TOFobs, ELOSSobs, Npe, beta
      real      ELOSStof, ELOSS(8,96), Nckov, Nck(96,48)
      integer   ick, jck, itf, jtf, Ncal/0/
      real      TOFcal(96),ENEcal(96),xcal(96),ycal(96)
      integer   IP(96)
*     -----------------------------------------------------------------
*
      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((ISWIT(7) .ne. 0) .and. (ISWIT(7) .ne. 3))then
         if((NAMES(3) .eq. 'SCTS').or.(NAMES(3) .eq. 'SCTL'))then
            TOFobs = TOFG * 1.0E+9 ! TOFG:(sec)  --> TOFobs:(nsec)
            ELOSSobs = DESTEP * 1.0E+3 ! DESTEP:(GeV)--> ELOSSobs:(MeV)
            ELOSStof = ELOSStof + ELOSSobs
            if(NAMES(3) .eq. 'SCTS')then
               itf = NUMBER(2)
               if((NUMBER(3) .ge. 1).and.(NUMBER(3) .le. 16))then
                  jtf = NUMBER(3)
               else
                  jtf = 113 - NUMBER(3)
               endif
            elseif(NAMES(3) .eq. 'SCTL')then
               itf = NUMBER(2)
               if((NUMBER(3).ge. 1).and.(NUMBER(3).le.16))then
                  jtf = 32 + NUMBER(3)
               elseif((NUMBER(3).ge.17).and.(NUMBER(3).le.32))then
                  jtf = 81 - NUMBER(3)
               elseif((NUMBER(3).ge.33).and.(NUMBER(3).le.48))then
                  jtf = 32 + NUMBER(3)
               elseif((NUMBER(3).ge.49).and.(NUMBER(3).le.64))then
                  jtf = 81 - NUMBER(3)
               endif
            endif
            ELOSS(itf,jtf) = ELOSS(itf,jtf) + ELOSSobs
            if(INWVOL.eq.1)then
               IP(jtf)   = IPART
               TOFcal(jtf) = TOFobs
               ENEcal(jtf) = GETOT*1.0E3
               xcal(jtf) = VECT(3)
               ycal(jtf) = VECT(2)
            endif
         endif
      endif
      if((NAMES(3).eq.'EMCL').and.(CHARGE.ne.0).and.(VECT(7).ne.0.))then
         beta = VECT(7)/GETOT   ! VECT(7):(GeV/c) GETOT:(GeV)
         Npe = STEP*370.*1.0*(1. - 1./(1.5*beta)**2) !1.5=kusseturitu
         if(1./(beta*1.5) .lt. 1.)then
            Nckov = Nckov + Npe
         endif
      endif
C
      if(INWVOL .eq. 2)then
         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,user common block
      common/uELOSS/In_part,IDcal,pt,Init_ENERGY,ELOSStof,ELOSS(8,96),
     +              Nckov,Nck(96,48)
      common/uENERGY/IP(96),TOFcal(96),ENEcal(96),xcal(96),ycal(96)
C
*KEND.
C
      real      event_tof(9),E_em
      integer   Ntof,itof(10),jtof(10),IPtof(10),temp
      real      pt,Init_ENERGY,ELOSStof,ELOSS(8,96),Nckov,Nck(96,48),Npi
      integer   icycle, itf, jtf, ick, jck, i, j, maxN, NN
*     ------------------------------------------------------------------
*
      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
*
c      write(6,*) IEVENT,Init_ENERGY,'(MeV) ToF ',ELOSStof,' EMcal',Nckov
C
      do i = 1,10
         IPtof(i) = 0
      enddo
      CALL HFILL(200,ELOSS(1,47),0.,1.)
C
      E_em = 5.995E-2*Nckov + 0.606   ! MeV : cut= 1MeV
c      E_em = 7.067E-2*Nckov + 0.634   ! MeV : cut= 5MeV
c      E_em = 9.328E-2*Nckov + 7.427   ! MeV : cut=10MeV
C
      if(E_em .gt. 20.)then
         CALL HFILL(300,E_em,0.,1.)
         CALL HFILL(310,E_em,0.,1.)
      endif
      itf = 1
      do jtf = 1,96
         event_tof(1) = IEVENT
         event_tof(2) = jtf
         event_tof(3) = IP(jtf)
         event_tof(4) = TOFcal(jtf)
         event_tof(5) = ENEcal(jtf)
         event_tof(6) = ELOSS(itf,jtf)
         event_tof(7) = xcal(jtf)
         event_tof(8) = ycal(jtf)
         event_tof(9) = E_em
         if((IP(jtf) .ge. 0.1).and.(abs(xcal(jtf)).le.50)
     &        .and.(jtf.ne.47))then
            CALL HFN(10,event_tof)
            if((jtf.eq.62).or.(jtf.eq.63).or.(jtf.eq.94)
     &           .or.(jtf.eq.79))then
               CALL HFN(20,event_tof)
               write(6,*)IEVENT,jtf,IP(jtf),ELOSS(itf,jtf),TOFcal(jtf)
            endif
         endif
         if((ELOSS(itf,jtf) .ge. 0.1).and.(abs(xcal(jtf)).le.50)
     &        .and.(jtf.ne.47))then
            CALL HFN(30,event_tof)
            if((jtf.eq.62).or.(jtf.eq.63).or.(jtf.eq.94)
     &           .or.(jtf.eq.79))then
               CALL HFN(40,event_tof)
               write(6,*)IEVENT,jtf,IP(jtf),ELOSS(itf,jtf),TOFcal(jtf)
            endif
         endif
      enddo
C
**
*
*             Draw track of particle and detector
*
*
      if( ISWIT(5) .eq. 0 )then
         if( ISWIT(2) .eq. 2 )then
            CALL GDRAW('MOTH',0.,0.,0.,0.5,10.,0.030,0.030)
         elseif( ISWIT(2) .eq. 3)then
            CALL GDRAW('MOTH',90.,90.,270.,0.5,10.,0.030,0.030)
         elseif( ISWIT(2) .eq. 4)then
            CALL GDRAW('MOTH',90.,0.,0.,10.,10.,0.030,0.030)
         elseif( ISWIT(2) .eq. 5)then
            CALL GDRAW('MOTH',90.,0.,0.,10.,10.,0.050,0.050)
         elseif( ISWIT(2) .eq. 6)then
            CALL GDRAW('MOTH',60.,40.,0.,1.,19.,0.050,0.050)
         else
            CALL GDRAW('MOTH',60.,120.,0.,3.,5.,0.030,0.030)
            CALL GDMAN(18.,4.0)
         endif
*
         CALL GDHEAD(101,'TOF and EM-cal$',.5)
         CALL GDXYZ(0)
         CALL GDAXIS(0.,0.,0.,100.)
         CALL GDPART(0,11,0.25)
         CALL GDSCAL(10.,2.)
      endif
*
      RETURN
 99      END
*
*
**
*
      subroutine pion_gen(y_min,y_max,y,p4,Init_pt,Iptswit,Iseed)
*
************************************************************************
*  this routine generates pions with constant rapidity density
*  and with the slope of 150 MeV/c.  pion0 = 139 MeV
*     y_min  ; INPUT  : Lower limit of pion rapidity
*     y_max  ; INPUT  : Higher limit of pion rapidity
*     y      ; OUTPUT : pion rapidity
*     p4     ; OUTPUT : Four Momentum of the particles
*     Init_pt; INPUT  : Initial momemtum Pt (if Iptswit=1)
*     Iptswit; INPUT  : Pt switch 0 or 1
*     Iseed  ; IN/OUT : INTEGER*4  Seed of random number
************************************************************************
*
      real y, y_min, y_max, mass, pt, p4(4), slope, sita, arg
C
      y = min(y_min,y_max) + abs(y_max-y_min)*ran(Iseed)  ! rapidity
C
      if(Iptswit .eq. 0)then
         slope = 1./6.          ! 1/6 [GeV/c] ==> 150 [MeV/c]
         pt = ran_pt(slope, Iseed)
      elseif(Iptswit .eq. 1)then
         pt = Init_pt
      endif
C
      arg = 2.*3.1415926*ran(Iseed)
*      sita = 0.178              !  = arctan(90/500)
*      arg = -sita + 2.*sita*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 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
*
*
