C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C  File        : RGSCIN.FOR
C  Author      : Jean Roy
C  Date        : 1-JUL-1990  (Oh Canada !)
C  Description : Function to return the range of a pion or muon in
C                scintillator as a function of momentum or kinetic energy.
C                The value is determined by interpolating the NASA table
C                values.
C                   Now possible to give a range in cm and have the
C                momentum returned by RGSCIN.
C
C  Input arguments :   EPR  Kinetic energy, momentum or range.
C                     IEPR  1 = kinetic energy  2 = momentum  3 = range
C                    ITYPE  1 = pion  2 = muon
C
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      REAL FUNCTION RGSCIN (EPRIN,IEPR,ITYPE)
C
      IMPLICIT NONE
C
      REAL     EPRIN,EPR,E(35),RGPI(35),RGMU(35),RHO,DIVDIF
      REAL     RMPI,RMMU,RMASS
      INTEGER  IEPR,ITYPE
C
      DATA     RHO   /1.032/
      DATA     RMPI  /139.56755/ , RMMU /105.65839/
C
      DATA     E /   2.0,  4.0,  6.0,  8.0, 10.0, 14.0, 18.0, 22.0,
     &              26.0, 30.0, 34.0, 38.0, 42.0, 46.0, 50.0, 60.0,
     &              70.0, 80.0, 90.0,100.0,110.0,120.0,130.0,140.0,
     &             150.0,160.0,170.0,180.0,190.0,200.0,220.0,240.0,
     &             260.0,280.0,300.0  /
C
      DATA     RGPI /   0.0313,  0.110,  0.229,  0.385,  0.576,
     &                  1.053 ,  1.646,  2.343,  3.134,  4.011,
     &                  4.967 ,  5.995,  7.089,  8.245,  9.458,
     &                 12.711 , 16.236, 19.989, 23.933, 28.042,
     &                 32.292 , 36.665, 41.143, 45.712, 50.361,
     &                 55.078 , 59.856, 64.687, 69.563, 74.480,
     &                 84.417 , 94.464,104.60 ,114.79 ,125.04   /
C
      DATA     RGMU /   0.0393,  0.138,  0.287,  0.482,  0.719,
     &                  1.307 ,  2.031,  2.875,  3.825,  4.870,
     &                  6.000 ,  7.207,  8.483,  9.821, 11.215,
     &                 14.920 , 18.877, 23.045, 27.386, 31.869,
     &                 36.469 , 41.166, 45.944, 50.791, 55.694,
     &                 60.646 , 65.639, 70.665, 75.721, 80.801,
     &                 91.018 ,101.29 ,111.60 ,121.94 ,132.29   /
C
      SAVE
C
C***********************************************************************
C  Get the appropriate rest mass.
C  If the momentum is given as argument, compute the kinetic energy.
C  Interpolate to find value of range and divide by scintillator density
C  to get the value in cm.
C  If the range (in cm) was given as an argument, return the
C  corresponding momentum.
C***********************************************************************
C
      EPR = EPRIN
      RGSCIN = 0.0
      IF (ITYPE .EQ. 1) THEN
        RMASS = RMPI
        IF (IEPR .EQ. 3) THEN
          EPR = EPR * RHO
          IF (EPR .LT. RGPI(1)) RETURN
          RGSCIN = DIVDIF(E,RGPI,35,EPR,2)
          RGSCIN = SQRT(RGSCIN**2 + 2.0*RMASS*RGSCIN)
        ELSE
          IF (EPR .LT. E(1)) RETURN
          IF (IEPR .EQ. 2) EPR = SQRT(EPR**2 + RMASS**2) - RMASS
          RGSCIN = DIVDIF(RGPI,E,35,EPR,2) / RHO
        ENDIF
      ELSE IF (ITYPE .EQ. 2) THEN
        RMASS = RMMU
        IF (IEPR .EQ. 3) THEN
          EPR = EPR * RHO
          IF (EPR .LT. RGMU(1)) RETURN
          RGSCIN = DIVDIF(E,RGMU,35,EPR,2)
          RGSCIN = SQRT(RGSCIN**2 + 2.0*RMASS*RGSCIN)
        ELSE
          IF (EPR .LT. E(1)) RETURN
          IF (IEPR .EQ. 2) EPR = SQRT(EPR**2 + RMASS**2) - RMASS
          RGSCIN = DIVDIF(RGMU,E,35,EPR,2) / RHO
        ENDIF
      ENDIF
C
C
      RETURN
C
      END