]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/pythia/pymaxi.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Precision parameter for pT sampling plus corresponding getter introduced.
[u/mrichter/AliRoot.git] / PYTHIA / pythia / pymaxi.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE PYMAXI
5
6C...Finds optimal set of coefficients for kinematical variable selection
7C...and the maximum of the part of the differential cross-section used
8C...in the event weighting.
9 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
10 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
11 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
12 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
13 COMMON/PYINT1/MINT(400),VINT(400)
14 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
15 COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
16 COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
17 COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
18 COMMON/PYINT6/PROC(0:200)
19 CHARACTER PROC*28
20 COMMON/PYINT7/SIGT(0:6,0:6,0:5)
21 SAVE /LUDAT1/,/LUDAT2/
22 SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,
23 &/PYINT5/,/PYINT6/,/PYINT7/
24 CHARACTER CVAR(4)*4
25 DIMENSION NPTS(4),MVARPT(500,4),VINTPT(500,30),SIGSPT(500),
26 &NAREL(7),WTREL(7),WTMAT(7,7),WTRELN(7),COEFU(7),COEFO(7),
27 &IACCMX(4),SIGSMX(4),SIGSSM(3)
28 DATA CVAR/'tau ','tau''','y* ','cth '/
29 DATA SIGSSM/3*0./
30
31C...Select subprocess to study: skip cases not applicable.
32 NPOSI=0
33 VINT(143)=1.
34 VINT(144)=1.
35 XSEC(0,1)=0.
36 DO 440 ISUB=1,200
37 MINT(51)=0
38 IF(ISET(ISUB).EQ.11) THEN
39 XSEC(ISUB,1)=1.00001*COEF(ISUB,1)
40 NPOSI=NPOSI+1
41 GOTO 430
42 ELSEIF(ISUB.GE.91.AND.ISUB.LE.95) THEN
43 XSEC(ISUB,1)=SIGT(0,0,ISUB-90)
44 IF(MSUB(ISUB).NE.1) GOTO 440
45 NPOSI=NPOSI+1
46 GOTO 430
47 ELSEIF(ISUB.EQ.96) THEN
48 IF(MINT(50).EQ.0) GOTO 440
49 IF(MSUB(95).NE.1.AND.MSTP(81).LE.0.AND.MSTP(131).LE.0) GOTO 440
50 IF(MINT(49).EQ.0.AND.MSTP(131).EQ.0) GOTO 440
51 ELSEIF(ISUB.EQ.11.OR.ISUB.EQ.12.OR.ISUB.EQ.13.OR.ISUB.EQ.28.OR.
52 &ISUB.EQ.53.OR.ISUB.EQ.68) THEN
53 IF(MSUB(ISUB).NE.1.OR.MSUB(95).EQ.1) GOTO 440
54 ELSE
55 IF(MSUB(ISUB).NE.1) GOTO 440
56 ENDIF
57 MINT(1)=ISUB
58 ISTSB=ISET(ISUB)
59 IF(ISUB.EQ.96) ISTSB=2
60 IF(MSTP(122).GE.2) WRITE(MSTU(11),5000) ISUB
61 MWTXS=0
62 IF(MSTP(142).GE.1.AND.ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+
63 &MSUB(94)+MSUB(95).EQ.0) MWTXS=1
64
65C...Find resonances (explicit or implicit in cross-section).
66 MINT(72)=0
67 KFR1=0
68 IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
69 KFR1=KFPR(ISUB,1)
70 ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
71 &ISUB.EQ.171.OR.ISUB.EQ.176) THEN
72 KFR1=23
73 ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
74 &ISUB.EQ.177) THEN
75 KFR1=24
76 ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
77 KFR1=25
78 IF(MSTP(46).EQ.5) THEN
79 KFR1=30
80 PMAS(30,1)=PARP(45)
81 PMAS(30,2)=PARP(45)**3/(96.*PARU(1)*PARP(47)**2)
82 ENDIF
83 ENDIF
84 CKMX=CKIN(2)
85 IF(CKMX.LE.0.) CKMX=VINT(1)
86 IF(KFR1.NE.0) THEN
87 IF(CKIN(1).GT.PMAS(KFR1,1)+20.*PMAS(KFR1,2).OR.
88 & CKMX.LT.PMAS(KFR1,1)-20.*PMAS(KFR1,2)) KFR1=0
89 ENDIF
90 IF(KFR1.NE.0) THEN
91 TAUR1=PMAS(KFR1,1)**2/VINT(2)
92 GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
93 MINT(72)=1
94 MINT(73)=KFR1
95 VINT(73)=TAUR1
96 VINT(74)=GAMR1
97 ENDIF
98 IF(ISUB.EQ.141) THEN
99 KFR2=23
100 TAUR2=PMAS(KFR2,1)**2/VINT(2)
101 GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
102 IF(CKIN(1).GT.PMAS(KFR2,1)+20.*PMAS(KFR2,2).OR.
103 & CKMX.LT.PMAS(KFR2,1)-20.*PMAS(KFR2,2)) KFR2=0
104 IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
105 MINT(72)=2
106 MINT(74)=KFR2
107 VINT(75)=TAUR2
108 VINT(76)=GAMR2
109 ELSEIF(KFR2.NE.0) THEN
110 KFR1=KFR2
111 TAUR1=TAUR2
112 GAMR1=GAMR2
113 MINT(72)=1
114 MINT(73)=KFR1
115 VINT(73)=TAUR1
116 VINT(74)=GAMR1
117 ENDIF
118 ENDIF
119
120C...Find product masses and minimum pT of process.
121 SQM3=0.
122 SQM4=0.
123 MINT(71)=0
124 VINT(71)=CKIN(3)
125 VINT(80)=1.
126 IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
127 NBW=0
128 DO 100 I=1,2
129 IF(KFPR(ISUB,I).EQ.0) THEN
130 ELSEIF(MSTP(42).LE.0.OR.PMAS(LUCOMP(KFPR(ISUB,I)),2).LT.
131 & PARP(41)) THEN
132 IF(I.EQ.1) SQM3=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
133 IF(I.EQ.2) SQM4=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
134 ELSE
135 NBW=NBW+1
136 ENDIF
137 100 CONTINUE
138 IF(NBW.GE.1) THEN
139 CALL PYOFSH(3,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
140 IF(MINT(51).EQ.1) THEN
141 WRITE(MSTU(11),5100) ISUB
142 MSUB(ISUB)=0
143 GOTO 440
144 ENDIF
145 SQM3=PQM3**2
146 SQM4=PQM4**2
147 ENDIF
148 IF(MIN(SQM3,SQM4).LT.CKIN(6)**2) MINT(71)=1
149 IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
150 IF(ISUB.EQ.96.AND.MSTP(82).LE.1) VINT(71)=PARP(81)
151 IF(ISUB.EQ.96.AND.MSTP(82).GE.2) VINT(71)=0.08*PARP(82)
152 ELSEIF(ISTSB.EQ.6) THEN
153 CALL PYOFSH(5,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
154 IF(MINT(51).EQ.1) THEN
155 WRITE(MSTU(11),5100) ISUB
156 MSUB(ISUB)=0
157 GOTO 440
158 ENDIF
159 SQM3=PQM3**2
160 SQM4=PQM4**2
161 ENDIF
162 VINT(63)=SQM3
163 VINT(64)=SQM4
164
165C...Prepare for additional variable choices in 2 -> 3.
166 IF(ISTSB.EQ.5) THEN
167 VINT(201)=0.
168 IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(KFPR(ISUB,2),1)
169 VINT(206)=VINT(201)
170 VINT(204)=PMAS(23,1)
171 IF(ISUB.EQ.124) VINT(204)=PMAS(24,1)
172 IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
173 & ISUB.EQ.186.OR.ISUB.EQ.187) VINT(204)=VINT(201)
174 VINT(209)=VINT(204)
175 ENDIF
176
177C...Number of points for each variable: tau, tau', y*, cos(theta-hat).
178 NPTS(1)=2+2*MINT(72)
179 IF(MINT(47).EQ.1) THEN
180 IF(ISTSB.EQ.1.OR.ISTSB.EQ.2.OR.ISTSB.EQ.6) NPTS(1)=1
181 ELSEIF(MINT(47).EQ.5) THEN
182 IF(ISTSB.LE.2.OR.ISTSB.GE.6) NPTS(1)=NPTS(1)+1
183 ENDIF
184 NPTS(2)=1
185 IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
186 IF(MINT(47).GE.2) NPTS(2)=2
187 IF(MINT(47).EQ.5) NPTS(2)=3
188 ENDIF
189 NPTS(3)=1
190 IF(MINT(47).GE.4) NPTS(3)=3
191 IF(MINT(45).EQ.3) NPTS(3)=NPTS(3)+1
192 IF(MINT(46).EQ.3) NPTS(3)=NPTS(3)+1
193 NPTS(4)=1
194 IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) NPTS(4)=5
195 NTRY=NPTS(1)*NPTS(2)*NPTS(3)*NPTS(4)
196
197C...Reset coefficients of cross-section weighting.
198 DO 110 J=1,20
199 COEF(ISUB,J)=0.
200 110 CONTINUE
201 COEF(ISUB,1)=1.
202 COEF(ISUB,8)=0.5
203 COEF(ISUB,9)=0.5
204 COEF(ISUB,13)=1.
205 COEF(ISUB,18)=1.
206 MCTH=0
207 MTAUP=0
208 METAUP=0
209 VINT(23)=0.
210 VINT(26)=0.
211 SIGSAM=0.
212
213C...Find limits and select tau, y*, cos(theta-hat) and tau' values,
214C...in grid of phase space points.
215 CALL PYKLIM(1)
216 METAU=MINT(51)
217 NACC=0
218 DO 140 ITRY=1,NTRY
219 MINT(51)=0
220 IF(METAU.EQ.1) GOTO 140
221 IF(MOD(ITRY-1,NPTS(2)*NPTS(3)*NPTS(4)).EQ.0) THEN
222 MTAU=1+(ITRY-1)/(NPTS(2)*NPTS(3)*NPTS(4))
223 IF(MTAU.GT.2+2*MINT(72)) MTAU=7
224 CALL PYKMAP(1,MTAU,0.5)
225 IF(ISTSB.GE.3.AND.ISTSB.LE.5) CALL PYKLIM(4)
226 METAUP=MINT(51)
227 ENDIF
228 IF(METAUP.EQ.1) GOTO 140
229 IF(ISTSB.GE.3.AND.ISTSB.LE.5.AND.MOD(ITRY-1,NPTS(3)*NPTS(4))
230 &.EQ.0) THEN
231 MTAUP=1+MOD((ITRY-1)/(NPTS(3)*NPTS(4)),NPTS(2))
232 CALL PYKMAP(4,MTAUP,0.5)
233 ENDIF
234 IF(MOD(ITRY-1,NPTS(3)*NPTS(4)).EQ.0) THEN
235 CALL PYKLIM(2)
236 MEYST=MINT(51)
237 ENDIF
238 IF(MEYST.EQ.1) GOTO 140
239 IF(MOD(ITRY-1,NPTS(4)).EQ.0) THEN
240 MYST=1+MOD((ITRY-1)/NPTS(4),NPTS(3))
241 IF(MYST.EQ.4.AND.MINT(45).NE.3) MYST=5
242 CALL PYKMAP(2,MYST,0.5)
243 CALL PYKLIM(3)
244 MECTH=MINT(51)
245 ENDIF
246 IF(MECTH.EQ.1) GOTO 140
247 IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
248 MCTH=1+MOD(ITRY-1,NPTS(4))
249 CALL PYKMAP(3,MCTH,0.5)
250 ENDIF
251 IF(ISUB.EQ.96) VINT(25)=VINT(21)*(1.-VINT(23)**2)
252
253C...Store position and limits.
254 MINT(51)=0
255 CALL PYKLIM(0)
256 IF(MINT(51).EQ.1) GOTO 140
257 NACC=NACC+1
258 MVARPT(NACC,1)=MTAU
259 MVARPT(NACC,2)=MTAUP
260 MVARPT(NACC,3)=MYST
261 MVARPT(NACC,4)=MCTH
262 DO 120 J=1,30
263 VINTPT(NACC,J)=VINT(10+J)
264 120 CONTINUE
265
266C...Normal case: calculate cross-section.
267 IF(ISTSB.NE.5) THEN
268 CALL PYSIGH(NCHN,SIGS)
269 IF(MWTXS.EQ.1) THEN
270 CALL PYEVWT(WTXS)
271 SIGS=WTXS*SIGS
272 ENDIF
273
274C..2 -> 3: find highest value out of a number of tries.
275 ELSE
276 SIGS=0.
277 DO 130 IKIN3=1,MSTP(129)
278 CALL PYKMAP(5,0,0.)
279 IF(MINT(51).EQ.1) GOTO 130
280 CALL PYSIGH(NCHN,SIGTMP)
281 IF(MWTXS.EQ.1) THEN
282 CALL PYEVWT(WTXS)
283 SIGTMP=WTXS*SIGTMP
284 ENDIF
285 IF(SIGTMP.GT.SIGS) SIGS=SIGTMP
286 130 CONTINUE
287 ENDIF
288
289C...Store cross-section.
290 SIGSPT(NACC)=SIGS
291 IF(SIGS.GT.SIGSAM) SIGSAM=SIGS
292 IF(MSTP(122).GE.2) WRITE(MSTU(11),5200) MTAU,MYST,MCTH,MTAUP,
293 &VINT(21),VINT(22),VINT(23),VINT(26),SIGS
294 140 CONTINUE
295 IF(NACC.EQ.0) THEN
296 WRITE(MSTU(11),5100) ISUB
297 MSUB(ISUB)=0
298 GOTO 440
299 ELSEIF(SIGSAM.EQ.0.) THEN
300 WRITE(MSTU(11),5300) ISUB
301 MSUB(ISUB)=0
302 GOTO 440
303 ENDIF
304 IF(ISUB.NE.96) NPOSI=NPOSI+1
305
306C...Calculate integrals in tau over maximal phase space limits.
307 TAUMIN=VINT(11)
308 TAUMAX=VINT(31)
309 ATAU1=LOG(TAUMAX/TAUMIN)
310 IF(NPTS(1).GE.2) THEN
311 ATAU2=(TAUMAX-TAUMIN)/(TAUMAX*TAUMIN)
312 ENDIF
313 IF(NPTS(1).GE.4) THEN
314 ATAU3=LOG(TAUMAX/TAUMIN*(TAUMIN+TAUR1)/(TAUMAX+TAUR1))/TAUR1
315 ATAU4=(ATAN((TAUMAX-TAUR1)/GAMR1)-ATAN((TAUMIN-TAUR1)/GAMR1))/
316 & GAMR1
317 ENDIF
318 IF(NPTS(1).GE.6) THEN
319 ATAU5=LOG(TAUMAX/TAUMIN*(TAUMIN+TAUR2)/(TAUMAX+TAUR2))/TAUR2
320 ATAU6=(ATAN((TAUMAX-TAUR2)/GAMR2)-ATAN((TAUMIN-TAUR2)/GAMR2))/
321 & GAMR2
322 ENDIF
323 IF(NPTS(1).GT.2+2*MINT(72)) THEN
324 ATAU7=LOG(MAX(2E-6,1.-TAUMIN)/MAX(2E-6,1.-TAUMAX))
325 ENDIF
326
327C...Reset. Sum up cross-sections in points calculated.
328 DO 300 IVAR=1,4
329 IF(NPTS(IVAR).EQ.1) GOTO 300
330 IF(ISUB.EQ.96.AND.IVAR.EQ.4) GOTO 300
331 NBIN=NPTS(IVAR)
332 DO 160 J1=1,NBIN
333 NAREL(J1)=0
334 WTREL(J1)=0.
335 COEFU(J1)=0.
336 DO 150 J2=1,NBIN
337 WTMAT(J1,J2)=0.
338 150 CONTINUE
339 160 CONTINUE
340 DO 170 IACC=1,NACC
341 IBIN=MVARPT(IACC,IVAR)
342 IF(IVAR.EQ.1.AND.IBIN.EQ.7) IBIN=3+2*MINT(72)
343 IF(IVAR.EQ.3.AND.IBIN.EQ.5.AND.MINT(45).NE.3) IBIN=4
344 NAREL(IBIN)=NAREL(IBIN)+1
345 WTREL(IBIN)=WTREL(IBIN)+SIGSPT(IACC)
346
347C...Sum up tau cross-section pieces in points used.
348 IF(IVAR.EQ.1) THEN
349 TAU=VINTPT(IACC,11)
350 WTMAT(IBIN,1)=WTMAT(IBIN,1)+1.
351 WTMAT(IBIN,2)=WTMAT(IBIN,2)+(ATAU1/ATAU2)/TAU
352 IF(NBIN.GE.4) THEN
353 WTMAT(IBIN,3)=WTMAT(IBIN,3)+(ATAU1/ATAU3)/(TAU+TAUR1)
354 WTMAT(IBIN,4)=WTMAT(IBIN,4)+(ATAU1/ATAU4)*TAU/
355 & ((TAU-TAUR1)**2+GAMR1**2)
356 ENDIF
357 IF(NBIN.GE.6) THEN
358 WTMAT(IBIN,5)=WTMAT(IBIN,5)+(ATAU1/ATAU5)/(TAU+TAUR2)
359 WTMAT(IBIN,6)=WTMAT(IBIN,6)+(ATAU1/ATAU6)*TAU/
360 & ((TAU-TAUR2)**2+GAMR2**2)
361 ENDIF
362 IF(NBIN.GT.2+2*MINT(72)) THEN
363 WTMAT(IBIN,NBIN)=WTMAT(IBIN,NBIN)+(ATAU1/ATAU7)*
364 & TAU/MAX(2E-6,1.-TAU)
365 ENDIF
366
367C...Sum up tau' cross-section pieces in points used.
368 ELSEIF(IVAR.EQ.2) THEN
369 TAU=VINTPT(IACC,11)
370 TAUP=VINTPT(IACC,16)
371 TAUPMN=VINTPT(IACC,6)
372 TAUPMX=VINTPT(IACC,26)
373 ATAUP1=LOG(TAUPMX/TAUPMN)
374 ATAUP2=((1.-TAU/TAUPMX)**4-(1.-TAU/TAUPMN)**4)/(4.*TAU)
375 WTMAT(IBIN,1)=WTMAT(IBIN,1)+1.
376 WTMAT(IBIN,2)=WTMAT(IBIN,2)+(ATAUP1/ATAUP2)*(1.-TAU/TAUP)**3/
377 & TAUP
378 IF(NBIN.GE.3) THEN
379 ATAUP3=LOG(MAX(2E-6,1.-TAUPMN)/MAX(2E-6,1.-TAUPMX))
380 WTMAT(IBIN,3)=WTMAT(IBIN,3)+(ATAUP1/ATAUP3)*
381 & TAUP/MAX(2E-6,1.-TAUP)
382 ENDIF
383
384C...Sum up y* cross-section pieces in points used.
385 ELSEIF(IVAR.EQ.3) THEN
386 YST=VINTPT(IACC,12)
387 YSTMIN=VINTPT(IACC,2)
388 YSTMAX=VINTPT(IACC,22)
389 AYST0=YSTMAX-YSTMIN
390 AYST1=0.5*(YSTMAX-YSTMIN)**2
391 AYST2=AYST1
392 AYST3=2.*(ATAN(EXP(YSTMAX))-ATAN(EXP(YSTMIN)))
393 WTMAT(IBIN,1)=WTMAT(IBIN,1)+(AYST0/AYST1)*(YST-YSTMIN)
394 WTMAT(IBIN,2)=WTMAT(IBIN,2)+(AYST0/AYST2)*(YSTMAX-YST)
395 WTMAT(IBIN,3)=WTMAT(IBIN,3)+(AYST0/AYST3)/COSH(YST)
396 IF(MINT(45).EQ.3) THEN
397 TAUE=VINTPT(IACC,11)
398 IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINTPT(IACC,16)
399 YST0=-0.5*LOG(TAUE)
400 AYST4=LOG(MAX(1E-6,EXP(YST0-YSTMIN)-1.)/
401 & MAX(1E-6,EXP(YST0-YSTMAX)-1.))
402 WTMAT(IBIN,4)=WTMAT(IBIN,4)+(AYST0/AYST4)/
403 & MAX(1E-6,1.-EXP(YST-YST0))
404 ENDIF
405 IF(MINT(46).EQ.3) THEN
406 TAUE=VINTPT(IACC,11)
407 IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINTPT(IACC,16)
408 YST0=-0.5*LOG(TAUE)
409 AYST5=LOG(MAX(1E-6,EXP(YST0+YSTMAX)-1.)/
410 & MAX(1E-6,EXP(YST0+YSTMIN)-1.))
411 WTMAT(IBIN,NBIN)=WTMAT(IBIN,NBIN)+(AYST0/AYST5)/
412 & MAX(1E-6,1.-EXP(-YST-YST0))
413 ENDIF
414
415C...Sum up cos(theta-hat) cross-section pieces in points used.
416 ELSE
417 RM34=MAX(1E-20,2.*SQM3*SQM4/(VINTPT(IACC,11)*VINT(2))**2)
418 RSQM=1.+RM34
419 CTHMAX=SQRT(1.-4.*VINT(71)**2/(TAUMAX*VINT(2)))
420 CTHMIN=-CTHMAX
421 IF(CTHMAX.GT.0.9999) RM34=MAX(RM34,2.*VINT(71)**2/
422 & (TAUMAX*VINT(2)))
423 ACTH1=CTHMAX-CTHMIN
424 ACTH2=LOG(MAX(RM34,RSQM-CTHMIN)/MAX(RM34,RSQM-CTHMAX))
425 ACTH3=LOG(MAX(RM34,RSQM+CTHMAX)/MAX(RM34,RSQM+CTHMIN))
426 ACTH4=1./MAX(RM34,RSQM-CTHMAX)-1./MAX(RM34,RSQM-CTHMIN)
427 ACTH5=1./MAX(RM34,RSQM+CTHMIN)-1./MAX(RM34,RSQM+CTHMAX)
428 CTH=VINTPT(IACC,13)
429 WTMAT(IBIN,1)=WTMAT(IBIN,1)+1.
430 WTMAT(IBIN,2)=WTMAT(IBIN,2)+(ACTH1/ACTH2)/MAX(RM34,RSQM-CTH)
431 WTMAT(IBIN,3)=WTMAT(IBIN,3)+(ACTH1/ACTH3)/MAX(RM34,RSQM+CTH)
432 WTMAT(IBIN,4)=WTMAT(IBIN,4)+(ACTH1/ACTH4)/MAX(RM34,RSQM-CTH)**2
433 WTMAT(IBIN,5)=WTMAT(IBIN,5)+(ACTH1/ACTH5)/MAX(RM34,RSQM+CTH)**2
434 ENDIF
435 170 CONTINUE
436
437C...Check that equation system solvable; else trivial way out.
438 IF(MSTP(122).GE.2) WRITE(MSTU(11),5400) CVAR(IVAR)
439 MSOLV=1
440 WTRELS=0.
441 DO 180 IBIN=1,NBIN
442 IF(MSTP(122).GE.2) WRITE(MSTU(11),5500) (WTMAT(IBIN,IRED),
443 &IRED=1,NBIN),WTREL(IBIN)
444 IF(NAREL(IBIN).EQ.0) MSOLV=0
445 WTRELS=WTRELS+WTREL(IBIN)
446 180 CONTINUE
447 IF(MSOLV.EQ.0) THEN
448 DO 190 IBIN=1,NBIN
449 COEFU(IBIN)=1.
450 WTRELN(IBIN)=0.1
451 IF(WTRELS.GT.0.) WTRELN(IBIN)=MAX(0.1,WTREL(IBIN)/WTRELS)
452 190 CONTINUE
453
454C...Solve to find relative importance of cross-section pieces.
455 ELSE
456 DO 200 IBIN=1,NBIN
457 WTRELN(IBIN)=MAX(0.1,WTREL(IBIN)/WTRELS)
458 200 CONTINUE
459 DO 230 IRED=1,NBIN-1
460 DO 220 IBIN=IRED+1,NBIN
461 RQT=WTMAT(IBIN,IRED)/WTMAT(IRED,IRED)
462 WTREL(IBIN)=WTREL(IBIN)-RQT*WTREL(IRED)
463 DO 210 ICOE=IRED,NBIN
464 WTMAT(IBIN,ICOE)=WTMAT(IBIN,ICOE)-RQT*WTMAT(IRED,ICOE)
465 210 CONTINUE
466 220 CONTINUE
467 230 CONTINUE
468 DO 250 IRED=NBIN,1,-1
469 DO 240 ICOE=IRED+1,NBIN
470 WTREL(IRED)=WTREL(IRED)-WTMAT(IRED,ICOE)*COEFU(ICOE)
471 240 CONTINUE
472 COEFU(IRED)=WTREL(IRED)/WTMAT(IRED,IRED)
473 250 CONTINUE
474 ENDIF
475
476C...Normalize coefficients, with piece shared democratically.
477 COEFSU=0.
478 WTRELS=0.
479 DO 260 IBIN=1,NBIN
480 COEFU(IBIN)=MAX(0.,COEFU(IBIN))
481 COEFSU=COEFSU+COEFU(IBIN)
482 WTRELS=WTRELS+WTRELN(IBIN)
483 260 CONTINUE
484 IF(COEFSU.GT.0.) THEN
485 DO 270 IBIN=1,NBIN
486 COEFO(IBIN)=PARP(122)/NBIN+(1.-PARP(122))*0.5*
487 & (COEFU(IBIN)/COEFSU+WTRELN(IBIN)/WTRELS)
488 270 CONTINUE
489 ELSE
490 DO 280 IBIN=1,NBIN
491 COEFO(IBIN)=1./NBIN
492 280 CONTINUE
493 ENDIF
494 IF(IVAR.EQ.1) IOFF=0
495 IF(IVAR.EQ.2) IOFF=17
496 IF(IVAR.EQ.3) IOFF=7
497 IF(IVAR.EQ.4) IOFF=12
498 DO 290 IBIN=1,NBIN
499 ICOF=IOFF+IBIN
500 IF(IVAR.EQ.1.AND.IBIN.GT.2+2*MINT(72)) ICOF=7
501 IF(IVAR.EQ.3.AND.IBIN.EQ.4.AND.MINT(45).NE.3) ICOF=ICOF+1
502 COEF(ISUB,ICOF)=COEFO(IBIN)
503 290 CONTINUE
504 IF(MSTP(122).GE.2) WRITE(MSTU(11),5600) CVAR(IVAR),
505 &(COEFO(IBIN),IBIN=1,NBIN)
506 300 CONTINUE
507
508C...Find two most promising maxima among points previously determined.
509 DO 310 J=1,4
510 IACCMX(J)=0
511 SIGSMX(J)=0.
512 310 CONTINUE
513 NMAX=0
514 DO 370 IACC=1,NACC
515 DO 320 J=1,30
516 VINT(10+J)=VINTPT(IACC,J)
517 320 CONTINUE
518 IF(ISTSB.NE.5) THEN
519 CALL PYSIGH(NCHN,SIGS)
520 IF(MWTXS.EQ.1) THEN
521 CALL PYEVWT(WTXS)
522 SIGS=WTXS*SIGS
523 ENDIF
524 ELSE
525 SIGS=0.
526 DO 330 IKIN3=1,MSTP(129)
527 CALL PYKMAP(5,0,0.)
528 IF(MINT(51).EQ.1) GOTO 330
529 CALL PYSIGH(NCHN,SIGTMP)
530 IF(MWTXS.EQ.1) THEN
531 CALL PYEVWT(WTXS)
532 SIGTMP=WTXS*SIGTMP
533 ENDIF
534 IF(SIGTMP.GT.SIGS) SIGS=SIGTMP
535 330 CONTINUE
536 ENDIF
537 IEQ=0
538 DO 340 IMV=1,NMAX
539 IF(ABS(SIGS-SIGSMX(IMV)).LT.1E-4*(SIGS+SIGSMX(IMV))) IEQ=IMV
540 340 CONTINUE
541 IF(IEQ.EQ.0) THEN
542 DO 350 IMV=NMAX,1,-1
543 IIN=IMV+1
544 IF(SIGS.LE.SIGSMX(IMV)) GOTO 360
545 IACCMX(IMV+1)=IACCMX(IMV)
546 SIGSMX(IMV+1)=SIGSMX(IMV)
547 350 CONTINUE
548 IIN=1
549 360 IACCMX(IIN)=IACC
550 SIGSMX(IIN)=SIGS
551 IF(NMAX.LE.1) NMAX=NMAX+1
552 ENDIF
553 370 CONTINUE
554
555C...Read out starting position for search.
556 IF(MSTP(122).GE.2) WRITE(MSTU(11),5700)
557 SIGSAM=SIGSMX(1)
558 DO 420 IMAX=1,NMAX
559 IACC=IACCMX(IMAX)
560 MTAU=MVARPT(IACC,1)
561 MTAUP=MVARPT(IACC,2)
562 MYST=MVARPT(IACC,3)
563 MCTH=MVARPT(IACC,4)
564 VTAU=0.5
565 VYST=0.5
566 VCTH=0.5
567 VTAUP=0.5
568
569C...Starting point and step size in parameter space.
570 DO 410 IRPT=1,2
571 DO 400 IVAR=1,4
572 IF(NPTS(IVAR).EQ.1) GOTO 400
573 IF(IVAR.EQ.1) VVAR=VTAU
574 IF(IVAR.EQ.2) VVAR=VTAUP
575 IF(IVAR.EQ.3) VVAR=VYST
576 IF(IVAR.EQ.4) VVAR=VCTH
577 IF(IVAR.EQ.1) MVAR=MTAU
578 IF(IVAR.EQ.2) MVAR=MTAUP
579 IF(IVAR.EQ.3) MVAR=MYST
580 IF(IVAR.EQ.4) MVAR=MCTH
581 IF(IRPT.EQ.1) VDEL=0.1
582 IF(IRPT.EQ.2) VDEL=MAX(0.01,MIN(0.05,VVAR-0.02,0.98-VVAR))
583 IF(IRPT.EQ.1) VMAR=0.02
584 IF(IRPT.EQ.2) VMAR=0.002
585 IMOV0=1
586 IF(IRPT.EQ.1.AND.IVAR.EQ.1) IMOV0=0
587 DO 390 IMOV=IMOV0,8
588
589C...Define new point in parameter space.
590 IF(IMOV.EQ.0) THEN
591 INEW=2
592 VNEW=VVAR
593 ELSEIF(IMOV.EQ.1) THEN
594 INEW=3
595 VNEW=VVAR+VDEL
596 ELSEIF(IMOV.EQ.2) THEN
597 INEW=1
598 VNEW=VVAR-VDEL
599 ELSEIF(SIGSSM(3).GE.MAX(SIGSSM(1),SIGSSM(2)).AND.
600 &VVAR+2.*VDEL.LT.1.-VMAR) THEN
601 VVAR=VVAR+VDEL
602 SIGSSM(1)=SIGSSM(2)
603 SIGSSM(2)=SIGSSM(3)
604 INEW=3
605 VNEW=VVAR+VDEL
606 ELSEIF(SIGSSM(1).GE.MAX(SIGSSM(2),SIGSSM(3)).AND.
607 &VVAR-2.*VDEL.GT.VMAR) THEN
608 VVAR=VVAR-VDEL
609 SIGSSM(3)=SIGSSM(2)
610 SIGSSM(2)=SIGSSM(1)
611 INEW=1
612 VNEW=VVAR-VDEL
613 ELSEIF(SIGSSM(3).GE.SIGSSM(1)) THEN
614 VDEL=0.5*VDEL
615 VVAR=VVAR+VDEL
616 SIGSSM(1)=SIGSSM(2)
617 INEW=2
618 VNEW=VVAR
619 ELSE
620 VDEL=0.5*VDEL
621 VVAR=VVAR-VDEL
622 SIGSSM(3)=SIGSSM(2)
623 INEW=2
624 VNEW=VVAR
625 ENDIF
626
627C...Convert to relevant variables and find derived new limits.
628 IF(IVAR.EQ.1) THEN
629 VTAU=VNEW
630 CALL PYKMAP(1,MTAU,VTAU)
631 IF(ISTSB.GE.3.AND.ISTSB.LE.5) CALL PYKLIM(4)
632 ENDIF
633 IF(IVAR.LE.2.AND.ISTSB.GE.3.AND.ISTSB.LE.5) THEN
634 IF(IVAR.EQ.2) VTAUP=VNEW
635 CALL PYKMAP(4,MTAUP,VTAUP)
636 ENDIF
637 IF(IVAR.LE.2) CALL PYKLIM(2)
638 IF(IVAR.LE.3) THEN
639 IF(IVAR.EQ.3) VYST=VNEW
640 CALL PYKMAP(2,MYST,VYST)
641 CALL PYKLIM(3)
642 ENDIF
643 IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
644 IF(IVAR.EQ.4) VCTH=VNEW
645 CALL PYKMAP(3,MCTH,VCTH)
646 ENDIF
647 IF(ISUB.EQ.96) VINT(25)=VINT(21)*(1.-VINT(23)**2)
648
649C...Evaluate cross-section. Save new maximum. Final maximum.
650 IF(ISTSB.NE.5) THEN
651 CALL PYSIGH(NCHN,SIGS)
652 IF(MWTXS.EQ.1) THEN
653 CALL PYEVWT(WTXS)
654 SIGS=WTXS*SIGS
655 ENDIF
656 ELSE
657 SIGS=0.
658 DO 380 IKIN3=1,MSTP(129)
659 CALL PYKMAP(5,0,0.)
660 IF(MINT(51).EQ.1) GOTO 380
661 CALL PYSIGH(NCHN,SIGTMP)
662 IF(MWTXS.EQ.1) THEN
663 CALL PYEVWT(WTXS)
664 SIGTMP=WTXS*SIGTMP
665 ENDIF
666 IF(SIGTMP.GT.SIGS) SIGS=SIGTMP
667 380 CONTINUE
668 ENDIF
669 SIGSSM(INEW)=SIGS
670 IF(SIGS.GT.SIGSAM) SIGSAM=SIGS
671 IF(MSTP(122).GE.2) WRITE(MSTU(11),5800) IMAX,IVAR,MVAR,IMOV,
672 &VNEW,VINT(21),VINT(22),VINT(23),VINT(26),SIGS
673 390 CONTINUE
674 400 CONTINUE
675 410 CONTINUE
676 420 CONTINUE
677 IF(MSTP(121).EQ.1) SIGSAM=PARP(121)*SIGSAM
678 XSEC(ISUB,1)=1.05*SIGSAM
679 430 CONTINUE
680 IF(MSTP(173).EQ.1.AND.ISUB.NE.96) XSEC(ISUB,1)=
681 &PARP(174)*XSEC(ISUB,1)
682 IF(ISUB.NE.96) XSEC(0,1)=XSEC(0,1)+XSEC(ISUB,1)
683 440 CONTINUE
684 MINT(51)=0
685
686C...Print summary table.
687 IF(NPOSI.EQ.0) THEN
688 WRITE(MSTU(11),5900)
689 STOP
690 ENDIF
691 IF(MSTP(122).GE.1) THEN
692 WRITE(MSTU(11),6000)
693 WRITE(MSTU(11),6100)
694 DO 450 ISUB=1,200
695 IF(MSUB(ISUB).NE.1.AND.ISUB.NE.96) GOTO 450
696 IF(ISUB.EQ.96.AND.MINT(50).EQ.0) GOTO 450
697 IF(ISUB.EQ.96.AND.MSUB(95).NE.1.AND.MSTP(81).LE.0) GOTO 450
698 IF(ISUB.EQ.96.AND.MINT(49).EQ.0.AND.MSTP(131).EQ.0) GOTO 450
699 IF(MSUB(95).EQ.1.AND.(ISUB.EQ.11.OR.ISUB.EQ.12.OR.ISUB.EQ.13.OR.
700 & ISUB.EQ.28.OR.ISUB.EQ.53.OR.ISUB.EQ.68)) GOTO 450
701 WRITE(MSTU(11),6200) ISUB,PROC(ISUB),XSEC(ISUB,1)
702 450 CONTINUE
703 WRITE(MSTU(11),6300)
704 ENDIF
705
706C...Format statements for maximization results.
707 5000 FORMAT(/1X,'Coefficient optimization and maximum search for ',
708 &'subprocess no',I4/1X,'Coefficient modes tau',10X,'y*',9X,
709 &'cth',9X,'tau''',7X,'sigma')
710 5100 FORMAT(1X,'Warning: requested subprocess ',I3,' has no allowed ',
711 &'phase space.'/1X,'Process switched off!')
712 5200 FORMAT(1X,4I4,F12.8,F12.6,F12.7,F12.8,1P,E12.4)
713 5300 FORMAT(1X,'Warning: requested subprocess ',I3,' has vanishing ',
714 &'cross-section.'/1X,'Process switched off!')
715 5400 FORMAT(1X,'Coefficients of equation system to be solved for ',A4)
716 5500 FORMAT(1X,1P,8E11.3)
717 5600 FORMAT(1X,'Result for ',A4,':',7F9.4)
718 5700 FORMAT(1X,'Maximum search for given coefficients'/2X,'MAX VAR ',
719 &'MOD MOV VNEW',7X,'tau',7X,'y*',8X,'cth',7X,'tau''',7X,'sigma')
720 5800 FORMAT(1X,4I4,F8.4,F11.7,F9.3,F11.6,F11.7,1P,E12.4)
721 5900 FORMAT(1X,'Error: no requested process has non-vanishing ',
722 &'cross-section.'/1X,'Execution stopped!')
723 6000 FORMAT(/1X,8('*'),1X,'PYMAXI: summary of differential ',
724 &'cross-section maximum search',1X,8('*'))
725 6100 FORMAT(/11X,58('=')/11X,'I',38X,'I',17X,'I'/11X,'I ISUB ',
726 &'Subprocess name',15X,'I Maximum value I'/11X,'I',38X,'I',
727 &17X,'I'/11X,58('=')/11X,'I',38X,'I',17X,'I')
728 6200 FORMAT(11X,'I',2X,I3,3X,A28,2X,'I',2X,1P,E12.4,3X,'I')
729 6300 FORMAT(11X,'I',38X,'I',17X,'I'/11X,58('='))
730
731 RETURN
732 END