***
*     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(3)*5
      data      chtags/'E','ToF','N'/
      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
*
*
*             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
      IPFIT=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,*)' '
c      write(6,*)' Select of CUT parameter'
c      write(6,*)'    cut for gammas    (CUTGAM)[GeV]'
c      read(5,*) CUTGAM
c      write(6,*)'    cut for electron  (CUTELE)[GeV]'
c      read(5,*) CUTELE
*
      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'
      read(5,*) ISWIT(7)
*
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
*
c      CALL HBOOKN(10,'Ntuple',3,'PHENIX',5000,chtags)
c      CALL HBOOKN(20,'Cherenkov Number',4,'PHENIX',5000,ckovtags)
c      CALL HBOOKN(30,'Energy of gamma',6,'PHENIX',5000,emtags)
      CALL HBOOKN(40,'Energy of pion0',6,'PHENIX',5000,pitags)
      CALL HBOOKN(50,'Energy of pion0 (trur)',6,'PHENIX',5000,pitags)
      CALL HBOOKN(60,'Energy of pion0 (B.G.)',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',100,-5.,5.,0.)
C
      if((ISWIT(7) .ne. 0) .and. (ISWIT(7) .ne. 3))then
         CALL HBOOK1(200,'(ToF) Energy loss(MeV)',100,0.,20.,0.)
         CALL HBOOK2(250,'ToF Pair pro. and Compton scat.',
     +        100,-200.,200.,50,-100.,100.,0.)
      endif
C
      CALL HBOOK1(300,'(EMcal) GAMMA (MeV)',200,0.,2000.,0.)
      CALL HBOOK1(310,' GAMMA sita-distribution ',120,0.,40.,0.)
      CALL HBOOK1(320,' pion --> 2-GAMMA (MeV)',100,0.,2000.,0.)
      CALL HBOOK1(330,' 2-GAMMA phi-distribution ',120,0.,60.,0.)
      CALL HBOOK1(340,' mass of pion',200,0.,1000.,0.)
      CALL HBOOK1(350,' mass of pion (true)',200,0.,1000.,0.)
      CALL HBOOK1(360,' mass of pion (B.G.)',200,0.,1000.,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),ITRA_pi(96,48)
C
*KEND.
C
      REAL      UBUF,PLAB(3),VERT(3)
      INTEGER   NVTX,NT,IPART
C
      real      pt, p4(4), Init_pt,Init_num, y_min, y_max, y
      integer   Iseed, Init_PART,In_part,Iptswit, IDtof, IDcal
      real      Nckov,Nck(96,48)
      integer   inum, ick, jck, ITRA_pi(96,48)
      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,*)' '
         WRITE(6,*)' Init_PART = 7 : pion0 '
c         READ(5,*)Init_PART
*
         WRITE(6,*)' '
         WRITE(6,*)'Number of paeticles  ..integer'
         WRITE(6,*)'  dN/dy = 375'
c         READ(5,*)Init_num
      ENDIF
*
      VERT(1) = 0.
      VERT(2) = 0.
      VERT(3) = 0.
      CALL GSVERT(VERT,0,0,UBUF,1,NVTX)
*
      y_min = -3.
      y_max =  3.
      Init_PART = 7
      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_min,y_max,y,p4,Init_pt,Iptswit,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
*
      ELOSStof = 0.
      Nckov = 0.
      do ick = 1,96
         do jck = 1,48
            Nck(ick,jck) = 0.
            ITRA_pi(ick,jck) = 0
            if(jck .le. 8)then
               ELOSS(jck,ick) = 0.
            endif
         enddo
      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
      common/uELOSS/In_part,IDcal,pt,Init_ENERGY,ELOSStof,ELOSS(8,96),
     +              Nckov,Nck(96,48),ITRA_pi(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),ITRA_pi(96,48)
C
*KEND.
C
      real      TOFobs, ELOSSobs, Npe, beta
      real      ELOSStof, ELOSS(8,96), Nckov, Nck(96,48)
      integer   ick, jck
*     -----------------------------------------------------------------
*
      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(ELOSSobs .ne. 0)then
               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
            endif
         endif
         if((NAMES(3) .eq. 'SCTS').or.(NAMES(3) .eq. 'SCTL').or.
     +        (NAMES(3) .eq. 'PMT ').or.(NAMES(3) .eq. 'MMTL').or.
     +        (NAMES(3) .eq. 'CRBN').or.(NAMES(3) .eq. 'HCMB').or.
     +        (NAMES(3) .eq. 'ALM1').or.(NAMES(3) .eq. 'ALM2').or.
     +        (NAMES(3) .eq. 'ALM3').or.(NAMES(3) .eq. 'CABL'))then
            if((LMEC(NMEC) .ge. 6).and.(LMEC(NMEC) .le. 8))then
               CALL HFILL(250,VECT(3),VECT(2),1.)
            endif
         endif
      endif
      if((NAMES(5).eq.'PbGy').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
            ick = NUMBER(4)
            jck = NUMBER(5)
            Nck(ick,jck)=Nck(ick,jck) + Npe
            if(ITRA_pi(ick,jck) .eq. 0)then
               ITRA_pi(ick,jck) = ITRA
            endif
         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),ITRA_pi(96,48)
C
*KEND.
C
      real      range_EM,event(3),event_ck(4),event_em(6),event_pi(6)
      real      pt,Init_ENERGY,ELOSStof,ELOSS(8,96),Nckov,Nck(96,48),Npi
      integer   icycle, ick, jck, i, j, maxN, NN
      parameter (maxN = 1000)
      real      Nem(maxN),E_em(maxN),Rem(maxN),sita_em(maxN),Rx,Ry
      integer   xem(maxN),yem(maxN),ITRA_pion(maxN)
      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
*
      write(6,*) IEVENT,Init_ENERGY,'(MeV) ToF ',ELOSStof,' EMcal',Nckov
C
      event(1) = Init_ENERGY
      event(2) = ELOSStof
      event(3) = Nckov
c      CALL HFN(10,event)
C
      do itf = 1,8
         do jtf = 1,96
            if(ELOSS(itf,jtf) .ge. 0.1)then
               CALL HFILL(200,ELOSS(itf,jtf),0.,1.)
            endif
         enddo
      enddo
C
      event_ck(1) = IEVENT
      event_ck(2) = Nckov
      event_ck(3) = 0.
      event_ck(4) = 0.
c      CALL HFN(20,event_ck)
C
      if(ISWIT(7) .eq. 3)then
         range_EM = 60.
      else
         range_EM = 530.
      endif
C
      NN = 0
      do ick = 1,96
         do jck = 1,48
            if(Nck(ick,jck) .ne. 0.)then
               event_ck(1) = IEVENT
               event_ck(2) = Nck(ick,jck)
               event_ck(3) = ick
               event_ck(4) = jck
               CALL HFN(20,event_ck)
               if((ick .ne. 1).and.(ick .ne. 96).and.
     +              (jck .ne. 1).and.(jck .ne. 48).and.
     +              (Nck(ick,jck) .gt. Nck(ick-1,jck-1)).and.
     +              (Nck(ick,jck) .gt. Nck(ick-1,jck)).and.
     +              (Nck(ick,jck) .gt. Nck(ick-1,jck+1)).and.
     +              (Nck(ick,jck) .gt. Nck(ick,jck-1)).and.
     +              (Nck(ick,jck) .gt. Nck(ick,jck+1)).and.
     +              (Nck(ick,jck) .gt. Nck(ick+1,jck-1)).and.
     +              (Nck(ick,jck) .gt. Nck(ick+1,jck)).and.
     +              (Nck(ick,jck) .gt. Nck(ick+1,jck+1)))then
                  NN = NN + 1
                  xem(NN) = ick
                  yem(NN) = jck
                  Rem(NN) = 4.1*sqrt((ick-48.5)**2 + (jck-24.5)**2)
                  Rx = 4.1*(ick - 48.5)
                  Ry = sqrt((4.1*(jck-24.5))**2 + range_EM**2)
                  sita_em(NN) = atan(Rx/Ry)*180./3.1415926
               endif
            endif
         enddo
      enddo
C
      if(NN .ge. 1)then
         do i = 1,NN
            Nem(i) = Nck(xem(i)-1,yem(i)-1)+Nck(xem(i)-1,yem(i))
     +           +Nck(xem(i)-1,yem(i)+1)+Nck(xem(i),yem(i)-1)
     +           +Nck(xem(i),yem(i))+Nck(xem(i),yem(i)+1)
     +           +Nck(xem(i)+1,yem(i)-1)+Nck(xem(i)+1,yem(i))
     +           +Nck(xem(i)+1,yem(i)+1)
C
            E_em(i) = 5.995E-2*Nem(i) + 0.606   ! MeV : cut= 1MeV
c            E_em(i) = 7.067E-2*Nem(i) + 0.634   ! MeV : cut= 5MeV
c            E_em(i) = 9.328E-2*Nem(i) + 7.427   ! MeV : cut=10MeV
C
            ITRA_pion(i) = ITRA_pi(xem(i),yem(i))
C
            event_em(1) = IEVENT
            event_em(2) = E_em(i)
            event_em(3) = xem(i)
            event_em(4) = yem(i)
            event_em(5) = Rem(i)
            event_em(6) = sita_em(i)
            if(E_em(i) .gt. 10.)then
c               CALL HFN(30,event_em)
               CALL HFILL(300,E_em(i),0.,1.)
               CALL HFILL(310,sita_em(i),0.,1.)
            endif
         enddo
      endif
      if(NN .ge. 2)then
         do i = 1,NN-1
            do j = i+1,NN
               if((E_em(i) .gt. 250.).and.(E_em(j) .gt. 250.))then
                  Npi = Nem(i) + Nem(j)
                  E_pi = E_em(i) + E_em(j)
                  Pt_pi = E_em(i)*cos(sita_em(i)*3.1415926/180.) 
     +                 +E_em(j)*cos(sita_em(j)*3.1415926/180.)
                  phi_pi=acos((Rem(i)**2+Rem(j)**2+2*range_EM**2
     +                 -((xem(i)-xem(j))**2+(yem(i)-yem(j))**2)*4.1**2)
     +                 /(2*sqrt((Rem(i)**2+range_EM**2)
     +                 *(Rem(j)**2 + range_EM**2))))*180./3.1415926
                  M_pi = 2*sqrt(E_em(i)*E_em(j))*
     +                 sin(phi_pi/2*3.1415926/180.)
C
                  event_pi(1) = IEVENT
                  if(E_em(i) .ge. E_em(j))then
                     event_pi(2) = E_em(i)
                     event_pi(3) = E_em(j)
                  else
                     event_pi(2) = E_em(j)
                     event_pi(3) = E_em(i)
                  endif
                  event_pi(4) = Pt_pi
                  event_pi(5) = phi_pi
                  event_pi(6) = M_pi
C
                  if(phi_pi .ge. 3.)then
                     CALL HFN(40,event_pi)
                     CALL HFILL(320,E_pi,0.,1.)
                     CALL HFILL(330,phi_pi,0.,1.)
                     CALL HFILL(340,M_pi,0.,1.)
                     if(ITRA_pion(i) .eq. ITRA_pion(j))then
                        CALL HFN(50,event_pi)
                        CALL HFILL(350,M_pi,0.,1.)
                     else
                        CALL HFN(60,event_pi)
                        CALL HFILL(360,M_pi,0.,1.)
                     endif
                  endif
               endif
            enddo
         enddo
      endif
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)
         else
            CALL GDRAW('MOTH',60.,120.,0.,3.,5.,0.030,0.030)
         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 GDMAN(18.,4.0)
         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
*
*
