ISAJET/code/ptfun.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE PTFUN
   3 C
   4 C          Calculate an envelope
   5 C          D(SIGMA)/D(PT**2)D(Y1)D(Y2) < PTFUN1*PT**PTFUN2
   6 C          used to generate initial PT values.
   7 C
   8 #if defined(CERNLIB_IMPNONE)
   9       IMPLICIT NONE
  10 #endif
  11 #include "isajet/itapes.inc"
  12 #include "isajet/keys.inc"
  13 #include "isajet/const.inc"
  14 #include "isajet/jetlim.inc"
  15 #include "isajet/ptpar.inc"
  16 #include "isajet/jetpar.inc"
  17 #include "isajet/jetsig.inc"
  18 C
  19       REAL PCPY(24)
  20       EQUIVALENCE(P(1),PCPY(1))
  21       REAL PTS(51),SIGSAV(51),STOR(24),DPT,DPTMIN,A,B,DEVMAX,DEV
  22       REAL DY1,DY2,B1
  23       INTEGER I,NPT,NDIV1,NDIV2,I1,I2
  24 C
  25       DATA DPTMIN/0.2/
  26 C
  27 C          Initialize
  28       DO 89 I=1,24
  29    89 STOR(I)=PCPY(I)
  30       YJ(1)=0
  31       YJ(2)=0
  32       TH(1)=PI/2.
  33       TH(2)=PI/2.
  34       STH(1)=1.
  35       STH(2)=1.
  36       CTH(1)=0.
  37       CTH(2)=0.
  38       PHI(1)=0.
  39       PHI(2)=PI
  40       IF(FIXPT(1).OR.FIXPT(2)) GOTO 300
  41       DPT=(PTMAX(1)-PTMIN(1))/25.
  42       IF(DPT.LT.DPTMIN) DPT=DPTMIN
  43       NPT=(PTMAX(1)-PTMIN(1))/DPT+1
  44       IF(NPT.GT.51) NPT=50
  45       IF(NPT.LE.1) NPT=2
  46 C
  47 C          Calculate sigma vs PT at Y1=Y2=0
  48       DO 100 I=1,NPT
  49         PT(1)=PTMIN(1)+DPT*(I-1)
  50         PT(2)=PT(1)
  51         P(1)=PT(1)
  52         P(2)=PT(2)
  53         IF(KEYS(1)) THEN
  54           CALL SIGQCD
  55         ELSEIF(KEYS(5)) THEN
  56           CALL SIGSSY
  57         ELSEIF(KEYS(6)) THEN
  58           CALL SIGWW
  59         ELSEIF(KEYS(8)) THEN
  60           CALL SIGGAM
  61         ELSEIF(KEYS(10)) THEN
  62           CALL SIGWH
  63         ENDIF
  64         IF(SIGMA.EQ.0.) GO TO 9999
  65         SIGSAV(I)=ALOG(SIGMA)
  66         PTS(I)=ALOG(PT(1))
  67   100 CONTINUE
  68 C
  69 C          Fit to power and shift to get envelope
  70 C
  71       CALL LSTSQ(PTS,SIGSAV,NPT,A,B)
  72       DEVMAX=0.
  73       DO 101 I=1,NPT
  74         DEV=SIGSAV(I)-A-B*PTS(I)
  75         IF(DEV.GT.DEVMAX) DEVMAX=DEV
  76   101 CONTINUE
  77 C
  78 C        Scan in Y1, Y2 for 3 PT values
  79 C
  80       DO 104 I=1,3
  81         IF(I.EQ.1) PT(1)=PTMIN(1)
  82         IF(I.EQ.2) PT(1)=(PTMIN(1)+PTMAX(1))/2.
  83         IF(I.EQ.3) PT(1)=PTMAX(1)
  84         PT(2)=PT(1)
  85         NDIV1=YJMAX(1)-YJMIN(1)
  86         IF(NDIV1.GT.20) NDIV1=20
  87         NDIV2=YJMAX(2)-YJMIN(2)
  88         IF(NDIV2.GT.20) NDIV2=20
  89         IF(NDIV1.LE.1) NDIV1=2
  90         IF(NDIV2.LE.1) NDIV2=2
  91         DY1=(YJMAX(1)-YJMIN(1))/(NDIV1-1)
  92         DY2=(YJMAX(2)-YJMIN(2))/(NDIV2-1)
  93         IF(FIXYJ(1)) NDIV1=1
  94         IF(FIXYJ(2)) NDIV2=1
  95 C
  96         DO 103 I1=1,NDIV1
  97           YJ(1)=YJMIN(1)+(I1-1)*DY1
  98           CTH(1)=TANH(YJ(1))
  99           STH(1)=SQRT(1.-CTH(1)**2)
 100           IF(STH(1).EQ.0) GOTO 103
 101           TH(1)=ACOS(CTH(1))
 102           P(1)=PT(1)/STH(1)
 103 C
 104           DO 102 I2=1,NDIV2
 105             YJ(2)=YJMIN(2)+(I2-1)*DY2
 106             CTH(2)=TANH(YJ(2))
 107             STH(2)=SQRT(1.-CTH(2)**2)
 108             IF(STH(2).EQ.0) GOTO 103
 109             TH(2)=ACOS(CTH(2))
 110             P(2)=PT(2)/STH(2)
 111             IF(KEYS(1)) THEN
 112               CALL SIGQCD
 113             ELSEIF(KEYS(5)) THEN
 114               CALL SIGSSY
 115             ELSEIF(KEYS(6)) THEN
 116               CALL SIGWW
 117             ELSEIF(KEYS(8)) THEN
 118               CALL SIGGAM
 119             ELSEIF(KEYS(10)) THEN
 120               CALL SIGWH
 121             ENDIF
 122             IF(SIGMA.EQ.0.) GO TO 102
 123             DEV=ALOG(SIGMA)-A-B*ALOG(PT(1))
 124             IF(DEV.GT.DEVMAX) DEVMAX=DEV
 125   102     CONTINUE
 126   103   CONTINUE
 127   104 CONTINUE
 128 C
 129       A=A+DEVMAX
 130       B1=B+2.
 131       PTFUN1=EXP(A)
 132       PTFUN2=B
 133 C
 134 C          Use envelope to generate initial PT values
 135 C
 136       PTGEN1=PTMIN(1)**B1
 137       PTGEN2=PTMAX(1)**B1-PTGEN1
 138       PTGEN3=1./B1
 139       DO 109 I=1,24
 140   109 PCPY(I)=STOR(I)
 141 C
 142 C          Write envelope parameters on listing
 143 C
 144       WRITE(ITLIS,200) PTFUN1,PTFUN2,PTGEN1,PTGEN2,PTGEN3
 145 200   FORMAT(//10X,'FIT AT Y1=Y2=0 IS D(SIGMA)/D(PT**2)D(Y1)D(Y2)='
 146      C  ,E11.5,'*PT**(',E11.5,')'//
 147      C  10X,'PT FIRST GENERATED BY PT=(',E11.5,'+',E11.5,'*RANF)**(',
 148      C  E11.5,')')
 149 C
 150       RETURN
 151 C
 152 C          Fixed PT
 153 C
 154   300 CONTINUE
 155       IF(FIXPT(1)) PT(2)=PT(1)
 156       IF(FIXPT(2)) PT(1)=PT(2)
 157       P(1)=PT(1)
 158       P(2)=PT(2)
 159       IF(KEYS(1)) THEN
 160         CALL SIGQCD
 161       ELSEIF(KEYS(5)) THEN
 162         CALL SIGSSY
 163       ELSEIF(KEYS(6)) THEN
 164         CALL SIGWW
 165       ELSEIF(KEYS(8)) THEN
 166         CALL SIGGAM
 167       ELSEIF(KEYS(10)) THEN
 168         CALL SIGWH
 169       ENDIF
 170       SIGMAX=SIGMA
 171       DO 301 I=1,24
 172 301   PCPY(I)=STOR(I)
 173 C
 174       RETURN
 175 C
 176 C          Fit fails if SIGMA=0 in specified range
 177 9999  WRITE(ITLIS,1010) PT(1)
 178 1010  FORMAT(//' ERROR IN PTFUN...SIGMA=0 FOR PT = ',E12.4/
 179      1' CHECK YOUR LIMITS.')
 180       STOP 99
 181       END