***
*     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
c      character chtags(3)*5
c      data      chtags/'E','ToF','N'/
c      character ckovtags(4)*5
c      data      ckovtags/'I','N','x','y'/
      character emtags(6)*5
      data      emtags/'I','E','x','y','r','sita'/
c      character pitags(6)*5
c      data      pitags/'I','E1','E2','Pt','phi','M'/
      character caltags(40)*6
      data caltags/'I','N','E0','Eless','Etot',
     +     'E1','IP1','x1','y1','i1','j1','E2','IP2','x2','y2','i2',
     +     'j2','E3','IP3','x3','y3','i3','j3','Nem','Eem1','iem1',
     +     'jem1','Ntof','Etof1','IPtof1','itof1','jtof1',
     +     'Etof2','IPtof2','itof2','jtof2',
     +     'Etof3','IPtof3','itof3','jtof3'/
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,*)' '
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)
      CALL HBOOKN(30,'Energy of gamma',6,'PHENIX',5000,emtags)
C
      CALL HBOOKN(60,'Energy(Pt= 250MeV)',40,'PHENIX',5000,caltags)
      CALL HBOOKN(70,'Energy(Pt= 500MeV)',40,'PHENIX',5000,caltags)
      CALL HBOOKN(80,'Energy(Pt=1000MeV)',40,'PHENIX',5000,caltags)
      CALL HBOOKN(90,'Energy(Pt=1500MeV)',40,'PHENIX',5000,caltags)
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.(P=250MeV/c)',
     +        100,-200.,200.,50,-100.,100.,0.)
         CALL HBOOK2(260,'ToF Pair pro. and Compton scat.(P=500MeV/c)',
     +        100,-200.,200.,50,-100.,100.,0.)
c         CALL HBOOK2(270,'ToF Pair pro. and Compton scat.(P=1000MeV/c)',
c     +        100,-200.,200.,50,-100.,100.,0.)
c         CALL HBOOK2(280,'ToF Pair pro. and Compton scat.(P=1500MeV/c)',
c     +        100,-200.,200.,50,-100.,100.,0.)
      endif
C
      CALL HBOOK1(300,'EMcal (P=250MeV/c)',100,0.,400,0.)
      CALL HBOOK1(310,'EMcal (P=250MeV/c)',200,0.,400,0.)
      CALL HBOOK1(320,'EMcal (P=500MeV/c)',100,0.,800,0.)
      CALL HBOOK1(330,'EMcal (P=500MeV/c)',200,0.,800,0.)
      CALL HBOOK1(340,'EMcal (P=1000MeV/c)',100,0.,1600,0.)
      CALL HBOOK1(350,'EMcal (P=1000MeV/c)',200,0.,1600,0.)
      CALL HBOOK1(360,'EMcal (P=1500MeV/c)',100,0.,2400,0.)
      CALL HBOOK1(370,'EMcal (P=1500MeV/c)',200,0.,2400,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/Etot,ENEcal(10),xcal(10),ycal(10),ical(10),
     +     jcal(10),Ncal,Ntof,ENEtof(10),itof(10),jtof(10),IP(10)
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,*)'Number of paeticles  ..integer'
         READ(5,*)Init_num
      ENDIF
*
      VERT(1) = 0.
      VERT(2) = 0.
      VERT(3) = 0.
      CALL GSVERT(VERT,0,0,UBUF,1,NVTX)
*
      if(IEVENT.EQ.1)then
         Init_pt = 0.25
      endif
      if(IEVENT.EQ.10001)then
         Init_pt = 0.50
      endif
      if(IEVENT.EQ.20001)then
         Init_pt = 1.0
      endif
      if(IEVENT.EQ.30001)then
         Init_pt = 1.5
      endif
*
      Iseed = NRNDM(2)+1812
      do inum = 1,Init_num
         thita = 0.400          !  = arctan(200/500)
         p_theta = 3.1415926/2 - thita + 2.*thita*ran(Iseed)
         phi = 0.200            !  = arctan(100/500)
         p_phi = -phi + 2.*phi*ran(Iseed)
*
         pt = Init_pt
         PLAB(1) = pt*sin(p_theta)*cos(p_phi)
         PLAB(2) = pt*sin(p_theta)*sin(p_phi)
         PLAB(3) = pt*cos(p_theta)
         CALL GSKINE(PLAB,Init_PART,NVTX,UBUF,1,NT)
      enddo
*
      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
      Etot = 0.
      Ncal = 0
      Ntof = 0
      do i = 1,10
         ENEcal(i) = 0.
         xcal(i) = 0.
         ycal(i) = 0.
         ical(i) = 0
         jcal(i) = 0
         ENEtof(i) = 0.
         itof(i) = 0
         jtof(i) = 0
         IP(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/Etot,ENEcal(10),xcal(10),ycal(10),ical(10),
     +     jcal(10),Ncal,Ntof,ENEtof(10),itof(10),jtof(10),IP(10)
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      Etot/0./,ENEcal(10),xcal(10),ycal(10)
      integer   ical(10),jcal(10),IP(10)
*     -----------------------------------------------------------------
*
      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
         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
c               write(6,*)NAMES(3),NAMEC(LMEC(NMEC)),ITRA,ISTAK,DESTEP
               if((pt .ge. 0.240) .and. (pt .lt. 0.260))then
                  CALL HFILL(250,VECT(3),VECT(2),1.)
               endif
               if((pt .ge. 0.490) .and. (pt .lt. 0.510))then
                  CALL HFILL(260,VECT(3),VECT(2),1.)
               endif
               if((pt .ge. 0.990) .and. (pt .lt. 1.010))then
                  CALL HFILL(270,VECT(3),VECT(2),1.)
               endif
               if((pt .ge. 1.490) .and. (pt .lt. 1.510))then
                  CALL HFILL(280,VECT(3),VECT(2),1.)
               endif
            endif
         endif
      endif
      if(((NAMES(3).eq.'SCTS').or.(NAMES(3).eq.'SCTL')).and.
     +     (Ncal.lt.10).and.(INWVOL.eq.1).and.
     +     (itf.lt.10).and.(jtf.lt.100))then
         Etot = Etot + GETOT*1.0E3
         Ncal = Ncal + 1
         ENEcal(Ncal) = GETOT*1.0E3
         IP(Ncal)   = IPART
         xcal(Ncal) = VECT(3)
         ycal(Ncal) = VECT(2)
         ical(Ncal) = itf
         jcal(Ncal) = jtf
      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
         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/Etot,ENEcal(10),xcal(10),ycal(10),ical(10),
     +     jcal(10),Ncal,Ntof,ENEtof(10),itof(10),jtof(10),IP(10)
C
*KEND.
C
      real      range_EM,event_em(6),event_cal(40)
      integer   level_cal(10)/1,2,3,4,5,6,7,8,9,10/
      integer   level_em(10)/1,2,3,4,5,6,7,8,9,10/
      integer   level_tof(10)/1,2,3,4,5,6,7,8,9,10/
      real      ENEtof(10)
      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
      parameter (maxN = 1000)
      real      Nem(maxN),E_em(maxN),Rem(maxN),sita_em(maxN)
      integer   xem(maxN),yem(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
*
c      write(6,*) IEVENT,Init_ENERGY,'(MeV) ToF ',ELOSStof,' EMcal',Nckov
C
      do i = 1,10
         IPtof(i) = 0
      enddo
      do itf = 1,8
         do jtf = 1,96
            if(ELOSS(itf,jtf) .ge. 0.1)then
               Ntof = Ntof + 1
               ENEtof(Ntof) = ELOSS(itf,jtf)
               itof(Ntof) = itf
               jtof(Ntof) = jtf
               do i = 1,10
                  if((itf .eq. ical(i)).and.(jtf .eq. jcal(i)))then
                     IPtof(Ntof) = IP(i)
                  endif
               enddo
               CALL HFILL(200,ELOSS(itf,jtf),0.,1.)
            endif
         enddo
      enddo
C
      if(ISWIT(7) .eq. 3)then
         range_EM = 60.
      else
         range_EM = 530.
      endif
C
      NN = 0
      do ick = 2,95
         do jck = 2,47
            if((Nck(ick,jck) .ne. 0.).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)
               sita_em(NN) = atan(Rem(NN)/range_EM)*180./3.1415926
            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
            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)
            CALL HFN(30,event_em)
C
            if(E_em(i) .gt. 20.)then
               if((pt .ge. 0.240) .and. (pt .lt. 0.260))then
                  CALL HFILL(300,E_em(i),0.,1.)
                  CALL HFILL(310,E_em(i),0.,1.)
               endif
               if((pt .ge. 0.490) .and. (pt .lt. 0.510))then
                  CALL HFILL(320,E_em(i),0.,1.)
                  CALL HFILL(330,E_em(i),0.,1.)
               endif
               if((pt .ge. 0.990) .and. (pt .lt. 1.010))then
                  CALL HFILL(340,E_em(i),0.,1.)
                  CALL HFILL(350,E_em(i),0.,1.)
               endif
               if((pt .ge. 1.490) .and. (pt .lt. 1.510))then
                  CALL HFILL(360,E_em(i),0.,1.)
                  CALL HFILL(370,E_em(i),0.,1.)
               endif
               CALL HFILL(300,E_em(i),0.,1.)
               CALL HFILL(310,sita_em(i),0.,1.)
            endif
         enddo
      endif
      do i = NN+1,10
         E_em(i) = 0.
         xem(i) = 0.
         yem(i) = 0.
      enddo
C
      do i = 1,9
         do j = i+1,10
            if(ENEcal(level_cal(i)) .lt. ENEcal(level_cal(j)))then
               temp = level_cal(i)
               level_cal(i) = level_cal(j)
               level_cal(j) = temp
            endif
            if(E_em(level_em(i)) .lt. E_em(level_em(j)))then
               temp = level_em(i)
               level_em(i) = level_em(j)
               level_em(j) = temp
            endif
            if(ENEtof(level_tof(i)) .lt. ENEtof(level_tof(j)))then
               temp = level_tof(i)
               level_tof(i) = level_tof(j)
               level_tof(j) = temp
            endif
         enddo
      enddo
C
      write(6,*) IEVENT,Init_ENERGY,'(MeV) N_ToF ',Ntof,' N_EMcal',Ncal
C
      event_cal(1) = IEVENT
      event_cal(2) = Ncal
      event_cal(3) = Init_ENERGY
      if((Ncal .eq. 0.).or.(Init_ENERGY - Etot .lt. 0.))then
         event_cal(4) = 0.
      else
         event_cal(4) = Init_ENERGY - Etot
      endif
      event_cal(5) = Etot
      event_cal(6) = ENEcal(level_cal(1))
      event_cal(7) = IP(level_cal(1))
      event_cal(8) = xcal(level_cal(1))
      event_cal(9) = ycal(level_cal(1))
      event_cal(10) = ical(level_cal(1))
      event_cal(11) = jcal(level_cal(1))
      event_cal(12) = ENEcal(level_cal(2))
      event_cal(13) = IP(level_cal(2))
      event_cal(14) = xcal(level_cal(2))
      event_cal(15) = ycal(level_cal(2))
      event_cal(16) = ical(level_cal(2))
      event_cal(17) = jcal(level_cal(2))
      event_cal(18) = ENEcal(level_cal(3))
      event_cal(19) = IP(level_cal(3))
      event_cal(20) = xcal(level_cal(3))
      event_cal(21) = ycal(level_cal(3))
      event_cal(22) = ical(level_cal(3))
      event_cal(23) = jcal(level_cal(3))
      event_cal(24) = NN
      event_cal(25) = E_em(level_em(1))
      event_cal(26) = xem(level_em(1))
      event_cal(27) = yem(level_em(1))
      event_cal(28) = Ntof
      event_cal(29) = ENEtof(level_tof(1))
      event_cal(30) = IPtof(level_tof(1))
      event_cal(31) = itof(level_tof(1))
      event_cal(32) = jtof(level_tof(1))
      event_cal(33) = ENEtof(level_tof(2))
      event_cal(34) = IPtof(level_tof(2))
      event_cal(35) = itof(level_tof(2))
      event_cal(36) = jtof(level_tof(2))
      event_cal(37) = ENEtof(level_tof(3))
      event_cal(38) = IPtof(level_tof(3))
      event_cal(39) = itof(level_tof(3))
      event_cal(40) = jtof(level_tof(3))
      
C
      if((pt .ge. 0.240) .and. (pt .lt. 0.260))then
         CALL HFN(60,event_cal)
      elseif((pt .ge. 0.490) .and. (pt .lt. 0.510))then
         CALL HFN(70,event_cal)
      elseif((pt .ge. 0.990) .and. (pt .lt. 1.010))then
         CALL HFN(80,event_cal)
      elseif((pt .ge. 1.490) .and. (pt .lt. 1.510))then
         CALL HFN(90,event_cal)
      endif
**
*
*             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)
            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
*
*
