EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lweits.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lweits.F
1 
2 C **********************************************************************
3 
4  SUBROUTINE lweits(LFILE)
5 
6  IMPLICIT NONE
7 
8 C...Integrates the QCD matrix elements to obtain probabilities for
9 C...qg- and qq-events as a function of (x,W). Also finds various
10 C...maximum values to be used for the QCD simulation. Results stored
11 C...in common LGRID and optionally written to logical file LFILE.
12 
13 *
14 * to avoid variable conflictions, a second keep element is necessary
15 * with the same common block name (see LPTOU2)
16 *
17 
18  include "leptou.inc"
19  include "linter.inc"
20  include "ludat1.inc"
21 
22  COMMON /loptim/ optx(4),opty(4),optq2(4),optw2(4),comfac
23  REAL optx,opty,optq2,optw2,comfac
24  SAVE /loptim/
25 
26  COMMON /lgrid/ nxx,nww,xx(31),ww(21),pqg(31,21,3),pqqb(31,21,2),
27  &qgmax(31,21,3),qqbmax(31,21,2),ycut(31,21),xtot(31,21),np
28  REAL xx,ww,pqg,pqqb,qgmax,qqbmax,ycut,xtot
29  INTEGER nxx,nww,np
30  SAVE /lgrid/
31 
32 
33  COMMON /lintrl/ psave(3,4,5),ksave(4),xmin,xmax,ymin,ymax,
34  &q2min,q2max,w2min,w2max,ilep,inu,ig,iz
35  REAL psave,xmin,xmax,ymin,ymax,q2min,q2max,w2min,w2max
36  INTEGER ksave,ilep,inu,ig,iz
37  SAVE /lintrl/
38 
39  INTEGER lfile,ncall,lst2,ipmax,iw,ix,lw,lx,iycut,ip
40  REAL wwi,xxi,wmax,pqcom,rqg,xpmin,xpmax,eps,rq,
41  &qtot,rqqb,w,result
42  REAL ysplit,xsplit,rtot,ulalps
43  INTEGER ny,nx,nl,i
44  REAL p27max,ycmin,ycmax
45  dimension wwi(21,4),xxi(31,4)
46  EXTERNAL dsigma,dsigm2
47  LOGICAL zoom
48 
49  DATA wwi/
50  1 5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,
51  1 20.,21.,22.,23.,24.,25.,
52  2 5.,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,
53  2 6.,6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.,7.3,
54 * 2 16.,17.,18.,19.,20.,21.,22.,23.,24.,25.,
55 * 2 5.,7.,9.,11.,13.,15.,17.,19.,21.,23.,25.,27.,29.,31.,33.,
56 * 2 35.,37.,39.,41.,43.,45.,
57  3 5.,7.,10.,15.,20.,25.,30.,40.,50.,75.,100.,125.,150.,175.,
58  3 200.,225.,250.,275.,300.,325.,350.,
59  4 5.,10.,15.,20.,30.,50.,75.,100.,150.,200.,250.,300.,400.,
60  4 500.,700.,900.,1200.,1500.,1800.,2100.,2500./
61  DATA xxi/
62  1 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,8.e-3,
63  1 1.e-2,2.e-2,3.e-2,4.e-2,5.e-2,6.e-2,8.e-2,
64  1 .1,.15,.2,.25,.3,.35,.4,.45,.5,.55,.6,.65,.72,.8,.87,.94,.999,
65 C.HI 1 1.E-2,2.E-2,2.5E-2,3.E-2,3.5E-2,4.E-2,4.5E-2,5.E-2,
66 C.HI 1 5.5E-2, 6.E-2,6.5E-2,7.E-2, 7.5E-2,8.E-2, 8.5E-2,
67 C.HI 1 .1,.15,.2,.25,.3,.35,.4,.45,.5,.55,.6,.65,.72,.8,.9,.999,
68  2 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,8.e-3,
69  2 1.e-2,2.e-2,3.e-2,4.e-2,5.e-2,6.e-2,8.e-2,
70  2 .1,.15,.2,.25,.3,.35,.4,.45,.5,.55,.6,.65,.72,.8,.87,.94,.999,
71  3 1.e-5,2.e-5,4.e-5,6.e-5,8.e-5,1.e-4,2.e-4,4.e-4,6.e-4,8.e-4,
72  3 1.e-3,2.e-3,4.e-3,6.e-3,8.e-3,1.e-2,2.e-2,4.e-2,6.e-2,8.e-2,
73  3 .1,.2,.3,.4,.5,.6,.7,.8,.87,.94,.999,
74  4 1.e-5,2.e-5,4.e-5,6.e-5,8.e-5,1.e-4,2.e-4,4.e-4,6.e-4,8.e-4,
75  4 1.e-3,2.e-3,4.e-3,6.e-3,8.e-3,1.e-2,2.e-2,4.e-2,6.e-2,8.e-2,
76  4 .1,.2,.3,.4,.5,.6,.7,.8,.87,.94,.999/
77  DATA ncall/0/,xsplit/0.1/,ysplit/2./
78 
79  ncall=ncall+1
80  lst2=lst(2)
81  lst(2)=-3
82  wmax=sqrt(parl(21))
83 
84  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
85  &WRITE(6,1000) parl(11),lst(13),mstu(112),paru(112),
86  &parl(8),parl(9),parl(12),parl(13)
87  IF(lst(17).EQ.0) THEN
88  np=1
89  ipmax=2
90  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,1010)
91  ELSE
92  np=3
93  IF(lst(23).EQ.1) np=2
94  ipmax=3
95  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) THEN
96  IF(lst(19).GE.0.AND.lst(19).LT.10) WRITE(6,1020)
97  IF(lst(19).GE.10.OR.lst(19).EQ.-10) WRITE(6,2020)
98  ENDIF
99  ENDIF
100 
101  IF(lst(19).GE.0.AND.lst(19).LT.10) THEN
102 C...Fixed grid in x and W
103  IF(lst(19).EQ.0)THEN
104 C...Grid specified by user.
105  READ(5,*) nww,nxx
106  READ(5,*) (ww(iw),iw=1,nww)
107  READ(5,*) (xx(ix),ix=1,nxx)
108  IF(xx(nxx).GT..99) xx(nxx)=.99
109  ELSEIF(lst(19).GE.1.AND.lst(19).LE.4) THEN
110 C...Grid taken from data in arrays WWI, XXI.
111  DO 10 iw=1,nww
112  10 ww(iw)=wwi(iw,lst(19))
113  DO 20 ix=1,nxx
114  20 xx(ix)=xxi(ix,lst(19))
115  ENDIF
116  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
117  & WRITE(6,1030) lst(19),nxx,nww,xx,ww
118  IF(wmax.GT.ww(nww)) THEN
119  IF(lst(3).GE.1) WRITE(6,1040) wmax,ww(nww)
120  IF(lst(3).GE.2) THEN
121  WRITE(6,1900)
122  stop
123  ENDIF
124  ENDIF
125  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,1100)
126  ELSEIF(lst(19).GE.10.OR.lst(19).EQ.-10) THEN
127 C...Grid in x,y automatically defined from available x,y-ranges
128  nx=nxx
129  ny=nww
130  IF(xmin.GE.xsplit) THEN
131  DO 30 i=1,nx
132  30 xx(i)=min(0.999,xmin+(xmax-xmin)*(i-1)/float(nx-1))
133  ELSEIF(xmax.GT.xsplit) THEN
134  nl=min(2.*nx/3.,max(nx/3.,nx*log(xsplit/xmin)/log(xmax/xmin)))
135  DO 40 i=1,nl
136  40 xx(i)=min(0.999,
137  & 10.**(log10(xmin)+(log10(xsplit)-log10(xmin))*(i-1)/float(nl)))
138  DO 41 i=nl+1,nx
139  41 xx(i)=min(0.999,xsplit+(xmax-xsplit)*(i-nl-1)/float(nx-nl-1))
140  ELSE
141  DO 50 i=1,nx
142  50 xx(i)=min(0.999,
143  & 10.**(log10(xmin)+(log10(xmax)-log10(xmin))*(i-1)/float(nx-1)))
144  ENDIF
145 C...y-variable stored in same array as W
146  IF(ymin.GE.ysplit) THEN
147  DO 60 i=1,ny
148  60 ww(i)=min(0.999,ymin+(ymax-ymin)*(i-1)/float(ny-1))
149  ELSEIF(ymax.GT.ysplit) THEN
150  nl=min(2.*ny/3.,max(ny/3.,ny*log(ysplit/ymin)/log(ymax/ymin)))
151  DO 70 i=1,nl
152  70 ww(i)=min(0.999,
153  & 10.**(log10(ymin)+(log10(ysplit)-log10(ymin))*(i-1)/float(nl)))
154  DO 71 i=nl+1,ny
155  71 ww(i)=min(0.999,ysplit+(ymax-ysplit)*(i-nl-1)/float(ny-nl-1))
156  ELSE
157  DO 80 i=1,ny
158  80 ww(i)=min(0.999,
159  & 10.**(log10(ymin)+(log10(ymax)-log10(ymin))*(i-1)/float(ny-1)))
160  ENDIF
161  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
162  & WRITE(6,2030) lst(19),nxx,nww,xx,ww
163  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,2100)
164  ENDIF
165 
166  lw=0
167  DO 500 iw=1,nww
168  w=ww(iw)
169  y=ww(iw)
170  IF(lw.GT.0) goto 600
171  IF(lst(19).LT.10.AND.w.GT.wmax) lw=lw+1
172  lx=0
173  DO 400 ix=1,nxx
174  x=xx(ix)
175  IF(lst(19).GE.0.AND.lst(19).LT.10) THEN
176 C...x,W given.
177  w2=w**2
178  u=(w2-psave(3,2,5)**2)/(2.*psave(3,2,5)*(1.-x))
179  q2=2.*psave(3,2,5)*u*x
180  y=q2/(parl(21)*x)
181  ELSEIF(lst(19).GE.10.OR.lst(19).EQ.-10) THEN
182 C...x,y given.
183  parl(22)=y*parl(21)
184  q2=x*parl(22)
185  u=parl(22)/(2.*psave(3,2,5))
186  w2=parl(22)*(1.-x)+psave(3,2,5)**2
187  w=sqrt(w2)
188  ENDIF
189 C...Protection against too small Q2 in structure functions
190 Ctest IF(Q2.LT.2.098) GOTO 400
191  IF(lx.GT.0) goto 500
192  IF(lst(19).GE.0.AND.lst(19).LT.10.AND.x.GT.1.-w2/wmax**2) lx=lx+1
193  CALL lepto
194  pqcom=pari(31)*pq(17)*comfac
195  parl(25)=ulalps(q2)
196  pari(20)=pq(17)
197  xtot(ix,iw)=pq(17)
198  IF(lst(20).LE.1) THEN
199  parl(27)=max(parl(9)**2/w2,parl(8))
200  p27max=1.0
201  ELSEIF(lst(20).EQ.2) THEN
202  parl(27)=max(parl(9)**2/q2,parl(8))
203  p27max=w2/q2
204  ELSEIF(lst(20).EQ.3.OR.lst(20).EQ.4) THEN
205  parl(27)=parl(8)
206  p27max=0.5
207  ELSEIF(lst(20).EQ.5) THEN
208  parl(27)=parl(8)
209  p27max=0.5
210  ELSEIF(lst(20).EQ.6) THEN
211  parl(27)=parl(9)
212  p27max=w2
213  ENDIF
214 
215  zoom=.false.
216  ycmin=parl(27)
217  ycmax=parl(27)
218  iycut=0
219  100 iycut=iycut+1
220  rqg=0.
221  rqqb=0.
222 CAE.Scheme for ME cutoff: W2, Q2, mixed
223  IF(lst(20).LE.1) THEN
224  xpmin=dble(x)/(1.d0-2.d0*(1.d0-dble(x))*dble(parl(27)))
225  xpmax=dble(x)/(dble(x)+(1.d0-dble(x))*dble(parl(27)))
226  ELSEIF(lst(20).EQ.2) THEN
227  xpmin=dble(x)/(1.d0-2.d0*dble(x)*dble(parl(27)))
228  xpmax=1.d0/(1.d0+dble(parl(27)))
229  ELSEIF(lst(20).EQ.3.OR.lst(20).EQ.4) THEN
230  xpmin=x
231  xpmax=1.d0/(1.d0+dble(parl(9)))
232  ELSEIF(lst(20).EQ.5) THEN
233  xpmin=x
234  xpmax=q2/(q2+parl(9))
235  ELSEIF(lst(20).EQ.6) THEN
236  xpmin=x
237  xpmax=q2/(q2+parl(27))
238  ELSE
239  WRITE(6,*) 'LWEITS: No such jet scheme!'
240  ENDIF
241 CAE
242  IF(xpmin.GE.xpmax.OR.xpmin.LE.0.) goto 210
243  DO 200 ip=1,np
244  IF(lst(17).EQ.0) THEN
245  pari(15)=0.
246  pari(16)=0.
247  pari(18)=0.
248  pari(19)=0.
249  ELSE
250  pari(14+ip)=0.
251  IF(ip.LE.2) pari(17+ip)=0.
252  ENDIF
253  lst(32)=ip
254  lst(24)=2
255  eps=parl(11)
256 CAE CALL GADAP(XPMIN,XPMAX,DSIGMA,EPS,RESULT)
257  CALL gadap(log(1.0-xpmax),log(1.0-xpmin),dsigm2,eps,result)
258  rqg=rqg+result
259  pqg(ix,iw,ip)=result/parl(25)
260  IF(lst(17).EQ.0) THEN
261  qgmax(ix,iw,1)=pari(15)
262  qgmax(ix,iw,2)=pari(16)
263  ELSE
264  pqg(ix,iw,ip)=result*pari(20)/pari(23+ip)/parl(25)
265  qgmax(ix,iw,ip)=pari(14+ip)
266  ENDIF
267  IF(ip.EQ.3) goto 200
268  lst(24)=3
269  eps=parl(11)
270 CAE CALL GADAP(XPMIN,XPMAX,DSIGMA,EPS,RESULT)
271  CALL gadap(log(1.0-xpmax),log(1.0-xpmin),dsigm2,eps,result)
272  rqqb=rqqb+result
273  pqqb(ix,iw,ip)=result/parl(25)
274  IF(lst(17).EQ.0) THEN
275  qqbmax(ix,iw,1)=pari(18)
276  qqbmax(ix,iw,2)=pari(19)
277  ELSE
278  pqqb(ix,iw,ip)=result*pari(20)/pari(23+ip)/parl(25)
279  qqbmax(ix,iw,ip)=pari(17+ip)
280  ENDIF
281  200 CONTINUE
282  210 CONTINUE
283  rq=1.-rqg-rqqb
284 *HI>>
285  IF(rq.lt.-20) rq=-20
286 *HI<<
287  IF(.NOT.zoom) THEN
288 CAE.First find interval so that RQ>0
289  IF(rq.LT.0.AND.iycut.LT.10) THEN
290  parl(27)=min(1.1*exp(-2.0*rq)*parl(27),p27max)
291  ycmin=ycmax
292  ycmax=parl(27)
293  ELSEIF(rq.LT.0.AND.iycut.GE.10) THEN
294 C...Terminate procedure after some iterations
295  WRITE(6,*) 'Warning! sigma>tot for x,q2,cut=',x,q2,parl(27)
296  WRITE(6,*) 'Weights=',rq,rqg,rqqb
297  rtot=(rqg+rqqb)*1.05
298  rqg=rqg/rtot
299  rqqb=rqqb/rtot
300  rq=1.-rqg-rqqb
301  DO 220 ip=1,3
302  pqg(ix,iw,ip)=pqg(ix,iw,ip)/rtot
303  qgmax(ix,iw,ip)=qgmax(ix,iw,ip)/rtot
304 220 CONTINUE
305  DO 230 ip=1,2
306  pqqb(ix,iw,ip)=pqqb(ix,iw,ip)/rtot
307  qqbmax(ix,iw,ip)=qqbmax(ix,iw,ip)/rtot
308 230 CONTINUE
309 C...Break loop
310  goto 250
311  ELSEIF(iycut.GE.2.AND.rq.GT.parl(13)) THEN
312 C...If RQ>PARL(13), then ycut was increased to much
313  zoom=.true.
314  parl(27)=(ycmin+ycmax)/2.
315  ELSE
316 C...correct ycut found
317  goto 250
318  ENDIF
319  ELSE
320 C...Zoom in on ycut so that 0<RQ<PARL(13)
321  IF(rq.LT.0.AND.iycut.LT.100) THEN
322  ycmin=parl(27)
323  parl(27)=(ycmin+ycmax)/2.
324  ELSEIF(rq.GT.parl(13).AND.iycut.LT.100) THEN
325  ycmax=parl(27)
326  parl(27)=(ycmin+ycmax)/2.
327  ELSEIF(iycut.GE.100) THEN
328  IF(lst(3).GE.1) THEN
329  WRITE(6,*) 'LWEITS: Warning, PARL(27) not found.'
330  ENDIF
331  rtot=(rqg+rqqb)*1.05
332  rqg=rqg/max(1.0,rtot)
333  rqqb=rqqb/max(1.0,rtot)
334  rq=1.-rqg-rqqb
335 C...Break loop
336  goto 250
337  ELSE
338 C...ycut found, break loop
339  goto 250
340  ENDIF
341  ENDIF
342 C...Loop until correct weights found
343  goto 100
344 
345  250 CONTINUE
346 CAE
347  ycut(ix,iw)=parl(27)
348  IF(lst(33).EQ.-91) THEN
349 Ctest...Include 3-jet cross section in denominator
350  qtot=1.+rqg+rqqb
351  rqg =rqg/qtot
352  rqqb=rqqb/qtot
353  rq=1.-rqg-rqqb
354  ENDIF
355  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) THEN
356  IF(lst(19).LT.10) WRITE(6,1200) w,x,y,q2,parl(25),pqcom,
357  & parl(27),iycut,rq,rqg,rqqb,(qgmax(ix,iw,ip),ip=1,ipmax),
358  & (qqbmax(ix,iw,ip),ip=1,min(2,ipmax))
359  IF(lst(19).GE.10) WRITE(6,2200) x,y,q2,w,parl(25),pqcom,
360  & parl(27),iycut,rq,rqg,rqqb,(qgmax(ix,iw,ip),ip=1,ipmax),
361  & (qqbmax(ix,iw,ip),ip=1,min(2,ipmax))
362  ENDIF
363  400 CONTINUE
364  500 CONTINUE
365  600 CONTINUE
366 
367  lst(2)=lst2
368  IF(lfile.LT.0) THEN
369 C...Write results on logical file number IABS(LFILE)
370  WRITE(iabs(lfile)) lst,parl,nxx,nww,np,xx,ww
371  WRITE(iabs(lfile))(((pqg(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,np),
372  & (((pqqb(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,np),
373  & (((qgmax(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,ipmax),
374  & (((qqbmax(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,min(2,ipmax)),
375  & ycut
376  IF(np.NE.1) WRITE(iabs(lfile)) xtot
377  CLOSE(iabs(lfile))
378  ENDIF
379  RETURN
380 
381  1000 FORMAT('1',/,5x,'Integration of 1st order QCD matrix elements',
382  & /,5x,'============================================',
383  &/,' for gluon radiation (qg-event) and boson-gluon fusion ',
384  &'(qq-event) probability.',
385  &//,' Required precision in integration, PARL(11) =',f8.4,
386  &//,' Heaviest flavour produced in boson-gluon fusion, LST(13) =',
387  &i5,//,' Alpha-strong parameters: # flavours, MSTU(112) =',i3,
388  &/,25x,' QCD lambda, PARU(112) =',f6.3,' GeV',
389  &//,' Cuts on matrix elements:',
390  &/,' PARL(8), PARL(9), PARL(12), PARL(13) =',4f8.4,/)
391  1010 FORMAT(' Lepton energy not allowed to vary in simulation.',/)
392  1020 FORMAT(' Lepton energy allowed to vary in simulation, ',/,
393  &' y in table below calculated assuming max energy.',/)
394  1030 FORMAT(' Grid choice, LST(19) =',i3,5x,'# grid points in x, W =',
395  &2i5,/,' x-values in array XX:',/,10f8.5,/,10f8.5,/,11f8.5,
396  & /,' W-values in array WW:',/,10f7.1,/,11f7.1,/)
397  1040 FORMAT(' Warning: max W outside grid, Wmax, grid-max =',2f12.1)
398  1100 FORMAT(
399  &9x,'cut',' it',2x,'q-event',1x,'qg-event',
400  &1x,'qq-event',' max of matrix elements qg & qq; L,R or T,S,I',
401  &/,1x,132(1h-),/)
402  1200 FORMAT(f7.1,2f8.4,1pg10.3,0pf6.2,1pg11.3,0pf8.4,i3,3f9.4,1p,5e9.2)
403  1900 FORMAT(' Execution stopped ',/)
404  2020 FORMAT(' Lepton energy allowed to vary in simulation, ',/,
405  &' W in table below calculated assuming max energy.',/)
406  2030 FORMAT(' Grid choice, LST(19) =',i3,5x,'# grid points in x, y =',
407  &2i5,/,' x-values in array XX:',/,10f8.5,/,10f8.5,/,11f8.5,
408  & /,' y-values in array WW:',/,10f7.4,/,11f7.4,/)
409  2100 FORMAT(
410  &9x,'cut',' it',2x,'q-event',1x,'qg-event',
411  &1x,'qq-event',' max of matrix elements qg & qq; L,R or T,S,I',
412  &/,1x,132(1h-),/)
413  2200 FORMAT(2f8.5,1pg10.3,0pf7.1,f6.2,1pg11.3,0pf8.4,i3,3f9.4,1p,5e9.2)
414  END