#include "isajet/pilot.h"
      SUBROUTINE PTFUN
C
C          Calculate an envelope
C          D(SIGMA)/D(PT**2)D(Y1)D(Y2) < PTFUN1*PT**PTFUN2
C          used to generate initial PT values.
C
#if defined(CERNLIB_IMPNONE)
      IMPLICIT NONE
#endif
#include "isajet/itapes.inc"
#include "isajet/keys.inc"
#include "isajet/const.inc"
#include "isajet/jetlim.inc"
#include "isajet/ptpar.inc"
#include "isajet/jetpar.inc"
#include "isajet/jetsig.inc"
C
      REAL PCPY(24)
      EQUIVALENCE(P(1),PCPY(1))
      REAL PTS(51),SIGSAV(51),STOR(24),DPT,DPTMIN,A,B,DEVMAX,DEV
      REAL DY1,DY2,B1
      INTEGER I,NPT,NDIV1,NDIV2,I1,I2
C
      DATA DPTMIN/0.2/
C
C          Initialize
      DO 89 I=1,24
   89 STOR(I)=PCPY(I)
      YJ(1)=0
      YJ(2)=0
      TH(1)=PI/2.
      TH(2)=PI/2.
      STH(1)=1.
      STH(2)=1.
      CTH(1)=0.
      CTH(2)=0.
      PHI(1)=0.
      PHI(2)=PI
      IF(FIXPT(1).OR.FIXPT(2)) GOTO 300
      DPT=(PTMAX(1)-PTMIN(1))/25.
      IF(DPT.LT.DPTMIN) DPT=DPTMIN
      NPT=(PTMAX(1)-PTMIN(1))/DPT+1
      IF(NPT.GT.51) NPT=50
      IF(NPT.LE.1) NPT=2
C
C          Calculate sigma vs PT at Y1=Y2=0
      DO 100 I=1,NPT
        PT(1)=PTMIN(1)+DPT*(I-1)
        PT(2)=PT(1)
        P(1)=PT(1)
        P(2)=PT(2)
        IF(KEYS(1)) THEN
          CALL SIGQCD
        ELSEIF(KEYS(5)) THEN
          CALL SIGSSY
        ELSEIF(KEYS(6)) THEN
          CALL SIGWW
        ELSEIF(KEYS(8)) THEN
          CALL SIGGAM
        ELSEIF(KEYS(10)) THEN
          CALL SIGWH
        ENDIF
        IF(SIGMA.EQ.0.) GO TO 9999
        SIGSAV(I)=ALOG(SIGMA)
        PTS(I)=ALOG(PT(1))
  100 CONTINUE
C
C          Fit to power and shift to get envelope
C
      CALL LSTSQ(PTS,SIGSAV,NPT,A,B)
      DEVMAX=0.
      DO 101 I=1,NPT
        DEV=SIGSAV(I)-A-B*PTS(I)
        IF(DEV.GT.DEVMAX) DEVMAX=DEV
  101 CONTINUE
C
C        Scan in Y1, Y2 for 3 PT values
C
      DO 104 I=1,3
        IF(I.EQ.1) PT(1)=PTMIN(1)
        IF(I.EQ.2) PT(1)=(PTMIN(1)+PTMAX(1))/2.
        IF(I.EQ.3) PT(1)=PTMAX(1)
        PT(2)=PT(1)
        NDIV1=YJMAX(1)-YJMIN(1) 
        IF(NDIV1.GT.20) NDIV1=20
        NDIV2=YJMAX(2)-YJMIN(2) 
        IF(NDIV2.GT.20) NDIV2=20
        IF(NDIV1.LE.1) NDIV1=2
        IF(NDIV2.LE.1) NDIV2=2
        DY1=(YJMAX(1)-YJMIN(1))/(NDIV1-1)
        DY2=(YJMAX(2)-YJMIN(2))/(NDIV2-1)
        IF(FIXYJ(1)) NDIV1=1
        IF(FIXYJ(2)) NDIV2=1
C
        DO 103 I1=1,NDIV1
          YJ(1)=YJMIN(1)+(I1-1)*DY1
          CTH(1)=TANH(YJ(1))
          STH(1)=SQRT(1.-CTH(1)**2)
          IF(STH(1).EQ.0) GOTO 103
          TH(1)=ACOS(CTH(1))
          P(1)=PT(1)/STH(1)
C
          DO 102 I2=1,NDIV2
            YJ(2)=YJMIN(2)+(I2-1)*DY2
            CTH(2)=TANH(YJ(2))
            STH(2)=SQRT(1.-CTH(2)**2)
            IF(STH(2).EQ.0) GOTO 103
            TH(2)=ACOS(CTH(2))
            P(2)=PT(2)/STH(2)
            IF(KEYS(1)) THEN
              CALL SIGQCD
            ELSEIF(KEYS(5)) THEN
              CALL SIGSSY
            ELSEIF(KEYS(6)) THEN
              CALL SIGWW
            ELSEIF(KEYS(8)) THEN
              CALL SIGGAM
            ELSEIF(KEYS(10)) THEN
              CALL SIGWH
            ENDIF
            IF(SIGMA.EQ.0.) GO TO 102
            DEV=ALOG(SIGMA)-A-B*ALOG(PT(1))
            IF(DEV.GT.DEVMAX) DEVMAX=DEV
  102     CONTINUE
  103   CONTINUE
  104 CONTINUE
C
      A=A+DEVMAX
      B1=B+2.
      PTFUN1=EXP(A)
      PTFUN2=B
C
C          Use envelope to generate initial PT values
C
      PTGEN1=PTMIN(1)**B1
      PTGEN2=PTMAX(1)**B1-PTGEN1
      PTGEN3=1./B1
      DO 109 I=1,24
  109 PCPY(I)=STOR(I)
C
C          Write envelope parameters on listing
C
      WRITE(ITLIS,200) PTFUN1,PTFUN2,PTGEN1,PTGEN2,PTGEN3
200   FORMAT(//10X,'FIT AT Y1=Y2=0 IS D(SIGMA)/D(PT**2)D(Y1)D(Y2)='
     C  ,E11.5,'*PT**(',E11.5,')'//
     C  10X,'PT FIRST GENERATED BY PT=(',E11.5,'+',E11.5,'*RANF)**(',
     C  E11.5,')')
C
      RETURN
C
C          Fixed PT
C
  300 CONTINUE
      IF(FIXPT(1)) PT(2)=PT(1)
      IF(FIXPT(2)) PT(1)=PT(2)
      P(1)=PT(1)
      P(2)=PT(2)
      IF(KEYS(1)) THEN
        CALL SIGQCD
      ELSEIF(KEYS(5)) THEN
        CALL SIGSSY
      ELSEIF(KEYS(6)) THEN
        CALL SIGWW
      ELSEIF(KEYS(8)) THEN
        CALL SIGGAM
      ELSEIF(KEYS(10)) THEN
        CALL SIGWH
      ENDIF
      SIGMAX=SIGMA
      DO 301 I=1,24
301   PCPY(I)=STOR(I)
C
      RETURN
C
C          Fit fails if SIGMA=0 in specified range
9999  WRITE(ITLIS,1010) PT(1)
1010  FORMAT(//' ERROR IN PTFUN...SIGMA=0 FOR PT = ',E12.4/
     1' CHECK YOUR LIMITS.')
      STOP 99
      END