EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lqcdpr.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lqcdpr.F
1 C **********************************************************************
2 
3  SUBROUTINE lqcdpr(QG,QQB)
4 
5  IMPLICIT NONE
6 
7 C...Probabilities for hard QCD events, qg or qqb, from integration of
8 C...QCD matrix elements event-by event or interpolation on x-W grid.
9 
10 *
11 * to avoid variable conflictions, a second keep element is necessary
12 * with the same common block name (see LPTOU2)
13 *
14  COMMON /leptou/ cut(14),lst(40),parl(30),
15  & x,y,w2,q2,u
16  REAL cut,parl,x,y,w2,q2,u
17  INTEGER lst
18  SAVE /leptou/
19 
20  COMMON /linter/ pari(50),ewqc(2,2,8),qc(8),zl(2,4),zq(2,8),pq(17)
21  REAL pari,ewqc,qc,zl,zq,pq
22  SAVE /linter/
23 
24  COMMON /lgrid/ nxx,nww,xx(31),ww(21),pqg(31,21,3),pqqb(31,21,2),
25  &qgmax(31,21,3),qqbmax(31,21,2),ycut(31,21),xtot(31,21),np
26  REAL xx,ww,pqg,pqqb,qgmax,qqbmax,ycut,xtot
27  INTEGER nxx,nww,np
28  SAVE /lgrid/
29 
30 
31 
32  common/debug1/iw,ix,ih
33 
34 
35  dimension pqsave(17)
36  EXTERNAL dsigma,dsigm2
37  INTEGER nout,nabove,nwarn,ih,ix,iw,ip
38  INTEGER i,iycut
39  REAL eps,xpmax,rq,pqsave,rtot,rqqb,rqg,ulalps
40  REAL xpmin
41  REAL qg,qqb,w,xp,wd,xd,x1p,x2p,qgip,pq17,qtot,qqbip,qgqqb
42  REAL p27max,ycmin,ycmax
43  DATA nout,nabove/2*0/,nwarn/10/
44 
45  LOGICAL zoom
46 
47 C...Get ycut from grid
48  IF(lst(19).GE.0.OR.lst(19).EQ.-10) THEN
49 C...
50 C...qg and qqb event probabilities from interpolation on grid
51  qg=0.
52  qqb=0.
53 C...QCD weight zero for x->1 above grid and W small below grid
54  IF(x.GT.xx(nxx).AND.x.GT.0.999) RETURN
55  IF(lst(19).LT.10.AND.sqrt(w2).LT.ww(1).AND.ww(1).LT.6.) RETURN
56 
57  xp=x
58 C...Local variable W is W or y
59  w=sqrt(w2)
60  IF(lst(19).GE.10.OR.lst(19).EQ.-10) w=y
61  IF(x.LT.xx(1).OR.x.GT.xx(nxx).OR.
62  &w.LT.ww(1).OR.w.GT.ww(nww)) THEN
63 C...x and/or W/y outside limits of grid, write warning NWARN first times
64  nout=nout+1
65  IF(lst(3).GE.1.AND.nout.LE.nwarn)
66  & WRITE(6,1000) x,w,int(pari(29)),nwarn
67  IF(x.LT.xx(1)) xp=xx(1)
68  IF(x.GT.xx(nxx)) xp=xx(nxx)
69  IF(w.LT.ww(1)) w=ww(1)
70  IF(w.GT.ww(nww)) w=ww(nww)
71  ENDIF
72 
73  ih=1
74  IF(lst(30).EQ.1) ih=2
75 *Pepsi>>
76  IF(lst(40).NE.0) ih=1
77 *Pepsi<<
78  ix=0
79  100 ix=ix+1
80  IF(xp.GT.xx(ix+1)) goto 100
81  iw=0
82  200 iw=iw+1
83  IF(w.GT.ww(iw+1)) goto 200
84  wd=(w-ww(iw))/(ww(iw+1)-ww(iw))
85  xd=(xp-xx(ix))/(xx(ix+1)-xx(ix))
86 
87  DO 500 ip=1,np
88  x1p=(pqg(ix+1,iw,ip)-pqg(ix,iw,ip))*xd+pqg(ix,iw,ip)
89  x2p=(pqg(ix+1,iw+1,ip)-pqg(ix,iw+1,ip))*xd+pqg(ix,iw+1,ip)
90  qgip=(x2p-x1p)*wd+x1p
91  IF(np.EQ.1) THEN
92  qg=qgip
93 
94  pari(15)=max(qgmax(ix,iw,ih),qgmax(ix+1,iw+1,ih),
95  & qgmax(ix+1,iw,ih),qgmax(ix,iw+1,ih))
96  ELSE
97  qg=qg+pari(23+ip)*qgip
98  pari(14+ip)=max(qgmax(ix,iw,ip),qgmax(ix+1,iw+1,ip),
99  & qgmax(ix+1,iw,ip),qgmax(ix,iw+1,ip))
100  ENDIF
101  IF(ip.EQ.3) goto 500
102  x1p=(pqqb(ix+1,iw,ip)-pqqb(ix,iw,ip))*xd+pqqb(ix,iw,ip)
103  x2p=(pqqb(ix+1,iw+1,ip)-pqqb(ix,iw+1,ip))*xd+pqqb(ix,iw+1,ip)
104  qqbip=(x2p-x1p)*wd+x1p
105  IF(np.EQ.1) THEN
106  qqb=qqbip
107  pari(18)=max(qqbmax(ix,iw,ih),qqbmax(ix+1,iw+1,ih),
108  & qqbmax(ix+1,iw,ih),qqbmax(ix,iw+1,ih))
109  ELSE
110  qqb=qqb+pari(23+ip)*qqbip
111  pari(17+ip)=max(qqbmax(ix,iw,ip),qqbmax(ix+1,iw+1,ip),
112  & qqbmax(ix+1,iw,ip),qqbmax(ix,iw+1,ip))
113  ENDIF
114  500 CONTINUE
115 
116  IF(np.NE.1) THEN
117 C...Get total x-section from interpolation to be used for normalization.
118  x1p=(xtot(ix+1,iw)-xtot(ix,iw))*xd+xtot(ix,iw)
119  x2p=(xtot(ix+1,iw+1)-xtot(ix,iw+1))*xd+xtot(ix,iw+1)
120  pq17=(x2p-x1p)*wd+x1p
121  qg=qg/pq17
122  qqb=qqb/pq17
123  ENDIF
124 
125 C..Interpolate in the grid
126  x1p=(ycut(ix+1,iw)-ycut(ix,iw))*xd+ycut(ix,iw)
127  x2p=(ycut(ix+1,iw+1)-ycut(ix,iw+1))*xd+ycut(ix,iw+1)
128  parl(27)=(x2p-x1p)*wd+x1p
129 C...Include alpha-strong in weight.
130  qg=qg*parl(25)
131  qqb=qqb*parl(25)
132 C...Get value of y-cut,
133  IF(lst(19).GE.0) THEN
134  IF(lst(33).EQ.-91) THEN
135 C...Include 3-jet cross section in denominator
136  qtot=1.+qg+qqb
137  qg =qg/qtot
138  qqb=qqb/qtot
139  ENDIF
140  IF(qg+qqb.GT.1) THEN
141 C...Sum of QCD event probabilities larger than unity, rescale to unity
142 C...and print warning for first NWARN cases.
143  nabove=nabove+1
144  IF(lst(3).GE.1.AND.nabove.LE.nwarn)
145  & WRITE(6,1100) qg,qqb,x,w,int(pari(29)),nwarn
146  qgqqb=qg+qqb
147  qg=qg/qgqqb
148  qqb=qqb/qgqqb
149  ENDIF
150  ELSE
151  IF(max(ycut(ix,iw),ycut(ix+1,iw+1),
152  & ycut(ix+1,iw),ycut(ix,iw+1))-
153  & min(ycut(ix,iw),ycut(ix+1,iw+1),
154  & ycut(ix+1,iw),ycut(ix,iw+1)).EQ.0.0) THEN
155  RETURN
156  ELSE
157 C...Get the minimum from the grid
158  parl(27)=min(ycut(ix,iw),ycut(ix+1,iw+1),
159  & ycut(ix+1,iw),ycut(ix,iw+1))
160  ENDIF
161  ENDIF
162 
163 C...Grid
164  ENDIF
165 
166 C...Calculate probabilities directly or refine value from grid
167  IF(lst(19).LE.0) THEN
168 
169 C...qg and qqbar event probabilities (and max values for simulation)
170 C...obtained by integrating QCD matrix elements for each event.
171 C LST2=LST(2)
172 C LST(2)=-3
173 C NP=1
174  lst(32)=1
175 
176  DO 1 i=1,17
177  1 pqsave(i)=pq(i)
178 
179  parl(25)=ulalps(q2)
180  pari(20)=pq(17)
181  IF(lst(19).GT.-10) THEN
182  IF(lst(20).LE.1) THEN
183  parl(27)=max(parl(9)**2/w2,parl(8))
184  p27max=1.0
185  ELSEIF(lst(20).EQ.2) THEN
186  parl(27)=max(parl(9)**2/q2,parl(8))
187  p27max=w2/q2
188  ELSEIF(lst(20).GE.3.AND.lst(20).LE.5) THEN
189  parl(27)=parl(8)
190  p27max=0.5
191  ELSEIF(lst(20).EQ.6) THEN
192  parl(27)=parl(9)
193  p27max=w2
194  ENDIF
195  ELSE
196  IF(lst(20).LE.1) THEN
197  p27max=1.0
198  ELSEIF(lst(20).EQ.2) THEN
199  p27max=w2/q2
200  ELSEIF(lst(20).GE.3.AND.lst(20).LE.5) THEN
201  p27max=0.5
202  ELSEIF(lst(20).EQ.6) THEN
203  p27max=w2
204  ENDIF
205  ENDIF
206 
207  zoom=.false.
208  iycut=0
209  ycmin=parl(27)
210  ycmax=parl(27)
211  10 iycut=iycut+1
212  rqg=0.
213  rqqb=0.
214 CAE.Scheme for ME cutoff: W2, Q2, mixed
215  IF(lst(20).LE.1) THEN
216  xpmin=dble(x)/(1.d0-2.d0*(1.d0-dble(x))*dble(parl(27)))
217  xpmax=dble(x)/(dble(x)+(1.d0-dble(x))*dble(parl(27)))
218  ELSEIF(lst(20).EQ.2) THEN
219  xpmin=dble(x)/(1.d0-2.d0*dble(x)*dble(parl(27)))
220  xpmax=1.d0/(1.d0+dble(parl(27)))
221  ELSEIF(lst(20).EQ.3.OR.lst(20).EQ.4) THEN
222  xpmin=x
223  xpmax=1./(1.+parl(9))
224  ELSEIF(lst(20).EQ.5) THEN
225  xpmin=x
226  xpmax=q2/(q2+parl(9))
227  ELSEIF(lst(20).EQ.6) THEN
228  xpmin=x
229  xpmax=q2/(q2+parl(27))
230  ELSE
231  WRITE(6,*) 'LQCDPR: No such jet scheme!'
232  ENDIF
233 CAE
234  IF(xpmin.LT.x.OR.xpmin.GT.1.) goto 40
235  IF(xpmin.GE.xpmax) goto 40
236 
237  pari(15)=0.
238  pari(16)=0.
239  pari(18)=0.
240  pari(19)=0.
241 C...QCD-Compton -> qg-event
242  lst(24)=2
243  eps=parl(11)
244 CAE CALL GADAP(XPMIN,XPMAX,DSIGMA,EPS,RQG)
245  CALL gadap(log(1.0-xpmax),log(1.0-xpmin),dsigm2,eps,rqg)
246 C...QCD-fusion -> qq-event
247  lst(24)=3
248  eps=parl(11)
249 CAE CALL GADAP(XPMIN,XPMAX,DSIGMA,EPS,RQQB)
250  CALL gadap(log(1.0-xpmax),log(1.0-xpmin),dsigm2,eps,rqqb)
251 C...q-event
252  rq=1.-rqg-rqqb
253 CAE WRITE(6,*) IYCUT,RQ,PARL(27),YCMIN,YCMAX
254  IF(.NOT.zoom) THEN
255 CAE.First find interval so that RQ>0
256  IF(rq.LT.0.AND.iycut.LT.10) THEN
257  parl(27)=min(1.1*exp(-2.0*rq)*parl(27),p27max)
258  ycmin=ycmax
259  ycmax=parl(27)
260  ELSEIF(rq.LT.0.AND.iycut.GE.10) THEN
261 C...Terminate procedure after some iterations
262  rtot=(rqg+rqqb)*1.05
263  rqg=rqg/rtot
264  rqqb=rqqb/rtot
265  rq=1.-rqg-rqqb
266 C IF(LST(3).GE.1) THEN
267 C WRITE(6,*) 'Warning! sigma>tot for x,q2,cut=',X,Q2,PARL(27)
268 C WRITE(6,*) 'Weights set to=',RQ,RQG,RQQB
269 C ENDIF
270 C...Break loop
271  goto 40
272  ELSEIF(iycut.GE.2.AND.rq.GT.parl(13)) THEN
273 C...If RQ>PARL(13), then ycut was increased to much
274  zoom=.true.
275  parl(27)=(ycmin+ycmax)/2.
276  ELSE
277 C...correct ycut found
278  goto 40
279  ENDIF
280  ELSE
281 C...Zoom in on ycut so that 0<RQ<PARL(13)
282  IF(rq.LT.0.AND.iycut.LT.40) THEN
283  ycmin=parl(27)
284  parl(27)=(ycmin+ycmax)/2.
285  ELSEIF(rq.GT.parl(13).AND.iycut.LT.40) THEN
286  ycmax=parl(27)
287  parl(27)=(ycmin+ycmax)/2.
288 C...Catch infinite loop
289  ELSEIF(iycut.GE.40) THEN
290  IF(lst(3).GE.1) THEN
291  WRITE(6,*) 'LQCDPR: Warning, PARL(27) not found.'
292  ENDIF
293  rtot=(rqg+rqqb)*1.05
294  rqg=rqg/max(1.0,rtot)
295  rqqb=rqqb/max(1.0,rtot)
296  rq=1.-rqg-rqqb
297 C...Break loop
298  goto 40
299  ELSE
300 C...ycut found, break loop
301  goto 40
302  ENDIF
303  ENDIF
304 C...Loop until correct weights found
305  goto 10
306 
307  40 CONTINUE
308 CAE
309  IF(lst(33).EQ.-91) THEN
310 C...Include 3-jet cross section in denominator
311  qtot=1.+rqg+rqqb
312  rqg =rqg/qtot
313  rqqb=rqqb/qtot
314  rq=1.-rqg-rqqb
315  ENDIF
316 
317 C LST(2)=LST2
318  lst(32)=0
319  DO 90 i=1,17
320  90 pq(i)=pqsave(i)
321  qg=rqg
322  qqb=rqqb
323 
324 C...Refine
325  ENDIF
326 
327  1000 FORMAT(' Warning: x=',f7.4,' or W/y=',f10.4,' outside QCD grid',
328  &' in event no.',i8,/,10x,
329  &'weight on limit of grid used. Only first',i5,' warnings printed')
330  1100 FORMAT(' Warning: Sum of QCD probabilities larger than unity ',
331  &' QG, QQB =',2f8.4,/10x,'occurs at x, W/y =',2f10.4,
332  &' in event no.',i8,/,10x,
333  &'Weights rescaled to unit sum. Only first',i5,' warnings printed')
334  RETURN
335  END