]>
Commit | Line | Data |
---|---|---|
0795afa3 | 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 |