EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pollinit.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pollinit.F
1 
2 
3 *72*********************************************************************
4  SUBROUTINE pollinit(LFILE,LEPIN,PLZ,PPZ,INTER)
5 
6  IMPLICIT NONE
7 
8 C...Initialize for an incoming lepton (type LEPIN, momentum pz=PLZ)
9 C...and target nucleon (momentum pz=PPZ) to interact via INTER.
10 C...Find maximum of differential cross section, calculate QCD event
11 C...probabilities or read them from logical file LFILE (if >0).
12 C... Numerical integration to obtain total cross-section.
13 
14  COMMON /lintrl/ psave(3,4,5),ksave(4),xmin,xmax,ymin,ymax,
15  &q2min,q2max,w2min,w2max,ilep,inu,ig,iz
16  REAL psave,xmin,xmax,ymin,ymax,q2min,q2max,w2min,w2max
17  INTEGER ksave,ilep,inu,ig,iz
18  SAVE /lintrl/
19 
20 *
21 * to avoid variable conflictions, a second keep element is necessary
22 * with the same common block name (see LPTOU2)
23 *
24  COMMON /leptou/ cut(14),lst(40),parl(30),
25  & x,y,w2,q2,u
26  REAL cut,parl,x,y,w2,q2,u
27  INTEGER lst
28  SAVE /leptou/
29 
30  COMMON /linter/ pari(50),ewqc(2,2,8),qc(8),zl(2,4),zq(2,8),pq(17)
31  REAL pari,ewqc,qc,zl,zq,pq
32  SAVE /linter/
33 
34  COMMON /lgrid/ nxx,nww,xx(31),ww(21),pqg(31,21,3),pqqb(31,21,2),
35  &qgmax(31,21,3),qqbmax(31,21,2),ycut(31,21),xtot(31,21),np
36  REAL xx,ww,pqg,pqqb,qgmax,qqbmax,ycut,xtot
37  INTEGER nxx,nww,np
38  SAVE /lgrid/
39 
40 
41  COMMON /loptim/ optx(4),opty(4),optq2(4),optw2(4),comfac
42  REAL optx,opty,optq2,optw2,comfac
43  SAVE /loptim/
44 
45  INTEGER nlupdm,nplbuf
46  parameter(nlupdm=4000,nplbuf=5)
47  common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
48  INTEGER n,k
49  REAL p,v
50  SAVE /lujets/
51 
52  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
53  INTEGER mstu,mstj
54  REAL paru,parj
55  SAVE /ludat1/
56 
57  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
58  INTEGER kchg
59  REAL pmas,parf,vckm
60  SAVE /ludat2/
61 
62  COMMON /lboost/ dbeta(2,3),stheta(2),sphi(2),pb(5),phir
63  DOUBLE PRECISION dbeta
64  REAL stheta,sphi,pb,phir
65  SAVE /lboost/
66 
67  COMMON /lminui/ xkin(4),ukin(4),wkin(4),ain(4),bin(4),
68  &maxfin,relup,relerr,reler2,fcnmax
69  REAL xkin,ukin,wkin,ain,bin,relerr,relup,reler2,fcnmax
70  INTEGER maxfin
71  SAVE /lminui/
72 
73  COMMON /lminuc/ namkin(4),nam(30)
74  CHARACTER*10 namkin,nam
75  SAVE /lminuc/
76 
77  COMMON /lpflag/ lst3
78  INTEGER lst3
79  SAVE /lpflag/
80 
81  COMMON /pypara/ ipy(80),pypar(80),pyvar(80)
82  REAL pypar,pyvar
83  INTEGER ipy
84  SAVE /pypara/
85 
86  common/pypars/mstp(200),parp(200),msti(200),paripy(200)
87  INTEGER mstp,msti
88  REAL parp,paripy
89  SAVE /pypars/
90 
91  common/maxfromhand/ fmaxfh
92  REAL fmaxfh
93  SAVE/maxfromhand/
94 
95 
96 *
97  INTEGER imxpdf
98  parameter(imxpdf=40)
99  COMMON /pepadm/cpdfnam(2,imxpdf),ipdfnam(2,imxpdf),
100  & iplst(10),cunpol,cpol
101  CHARACTER*256 cpdfnam,cunpol,cpol
102  INTEGER iplst,ipdfnam
103  SAVE /pepadm/
104 
105 **************************************************************
106 *
107 * IPLST(1) = 0 (default) : number of PDF warnings
108 * IPLST(2) = 11 (default) : unit -1- for pdf files
109 * IPLST(3) = 12 (default) : unit -2- for pdf files
110 *
111 **************************************************************
112 
113 
114  INTEGER inter,lfile,lepin,ncall,i,j,ifl
115  REAL plz,ppz,pi,pm2,s,roots
116  REAL ulmass,ulangl,umin,umax
117  INTEGER lqcd,ltm,ipmax,ip,iw,ix
118  REAL ti1,ti2
119 
120  INTEGER lstw(40)
121  REAL parlw(30)
122 
123  DATA pi/3.1415927/,ncall/0/
124 
125  ncall=ncall+1
126  lst3=lst(3)
127 
128 *MM,HI >>
129  DO i=1,imxpdf
130  IF(lst(15).eq.ipdfnam(1,i)) THEN
131  cunpol=cpdfnam(1,ipdfnam(2,i))
132  cpol =cpdfnam(2,ipdfnam(2,i))
133  ENDIF
134  ENDDO
135  IF (lst(15).GE.113.AND.lst(15).LE.115) CALL polini
136  IF (lst(15).GE.107.AND.lst(15).LE.109) CALL nloini
137  IF ((lst(15).GE.124).AND.(lst(15).LE.129)) CALL inideflo
138  IF (lst(40).EQ.0) RETURN
139 *MM,HI<<
140 
141  IF(lst(18).GE.1) THEN
142 C...W, Z masses from theta-Weinberg, Fermi constant GF and rad. corr.
143  pmas(24,1)=sqrt(pi*parl(16)/(sqrt(2.)*parl(17)*parl(5)*
144  & (1.-parl(18))))
145  pmas(23,1)=pmas(24,1)/sqrt(1.-parl(5))
146  ENDIF
147 C...Couplings between Z0 and left/right-handed leptons and quarks.
148  zl(1,1)=-.5+parl(5)
149  zl(1,2)=parl(5)
150  zl(2,1)=zl(1,2)
151  zl(2,2)=zl(1,1)
152  zl(1,3)=0.5
153  zl(2,3)=0.
154  zl(1,4)=0.
155  zl(2,4)=0.5
156  DO 10 ifl=1,8
157  zq(1,ifl)=sign(0.5,qc(ifl))-qc(ifl)*parl(5)
158  10 zq(2,ifl)=-qc(ifl)*parl(5)
159 
160 C...Set initial state.
161  lst(23)=inter
162  ksave(1)=lepin
163  ksave(2)=2212
164  k(1,1)=21
165  k(1,2)=ksave(1)
166  k(1,3)=0
167  k(1,4)=0
168  k(1,5)=0
169  k(2,1)=21
170  k(2,2)=ksave(2)
171  k(2,3)=0
172  k(2,4)=0
173  k(2,5)=0
174  p(1,1)=0.
175  p(1,2)=0.
176  p(1,3)=plz
177  p(1,5)=ulmass(ksave(1))
178  p(1,4)=sqrt(p(1,3)**2+p(1,5)**2)
179  p(2,1)=0.
180  p(2,2)=0.
181  p(2,3)=ppz
182  p(2,5)=ulmass(ksave(2))
183  p(2,4)=sqrt(p(2,3)**2+p(2,5)**2)
184  n=2
185  lst(28)=3
186 C...Save momentum vectors of incoming particles
187  DO 20 i=1,2
188  DO 20 j=1,5
189  20 psave(3,i,j)=p(i,j)
190 C...Dot-product of initial particles, cms energy
191  parl(21)=2.*(dble(p(1,4))*dble(p(2,4))-dble(p(1,3))*dble(p(2,3)))
192  roots=sqrt((dble(p(1,4))+dble(p(2,4)))**2
193  & -(dble(p(1,3))+dble(p(2,3)))**2)
194  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,1000)
195  &lepin,(p(1,j),j=1,3),parl(1),parl(2),(p(2,j),j=1,3),inter,roots
196 *...Write out also polarization state and structure function
197  IF(lst(40).NE.0)WRITE(6,1005) lst(40),lst(15)
198 *---
199  IF(plz*ppz.GT.0.1) THEN
200  WRITE(6,1010)
201  stop
202  ENDIF
203 
204 C...Reduced header for Jetset/Pythia
205  mstu(12)=0
206  mstp(122)=0
207  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
208  &WRITE(6,1020) mstu(181),mstu(182),mstp(181),mstp(182)
209 C...If JETSET version before 7.402, problem with azimuthal dependence
210 C...in LUSHOW solved by chosing flat azimuthal dependence.
211  IF(mstu(181).LE.7.AND.mstu(182).LT.402) THEN
212  mstj(46)=0
213  WRITE(6,1030) mstj(46)
214  ENDIF
215 C...Initialize PYTHIA for parton densities.
216  IF(lst(15).GT.0) THEN
217 C...Set switches and parameters for parton densities in LYSTFU.
218  mstp(51)=lst(15)
219  mstp(52)=lst(16)
220  mstp(58)=lst(12)
221  ENDIF
222 C CALL PYINIT('NONE','e-','p',ROOTS)
223  parl(26)=parp(1)
224 CAE-- use Lambda from parton densities in initial cascade
225  pypar(21)=parp(1)
226 C...Reset PYTHIA 4.8 parameters from LEPTO parameters.
227  IF(mod(lst(8),10).EQ.3.OR.mod(lst(8),10).EQ.5) ipy(13)=0
228  IF(lst(35).NE.1.AND.
229  &(mod(lst(8),10).EQ.4.OR.mod(lst(8),10).EQ.5)) ipy(14)=0
230  ipy(8)=lst(12)
231 
232 
233 *...Special restrictions apply for polarized version: Only first
234 *...order QCD. Only EM interaction. Beam energy not allowed to vary.
235 *...Parton cascades, longitudinal SF and higher twist not allowed.
236 C IF((LST(40).NE.0.AND.LST(23).NE.1).OR.
237  IF((lst(40).NE.0.AND.lst(23).GE.3).OR.
238  &(lst(40).NE.0.AND.lst(8).GE.2).OR.
239  &(lst(40).NE.0.AND.lst(17).NE.0).OR.
240 * &(LST(40).NE.0.AND.LST(11).NE.0).OR.
241  &(lst(40).NE.0.AND.iabs(lepin).NE.11)) THEN
242  WRITE(6,1025)
243  stop
244  ENDIF
245 *---
246  IF(psave(3,1,3).LT.0.) THEN
247 C...Flip event to have initial lepton along +z axis
248  p(1,3)=-p(1,3)
249  p(2,3)=-p(2,3)
250  ENDIF
251 C...Boost parameters to cms of incoming particles
252  dbeta(1,1)=0.d0
253  dbeta(1,2)=0.d0
254  dbeta(1,3)=(dble(p(1,3))+dble(p(2,3)))/(dble(p(1,4))+dble(p(2,4)))
255  sphi(1)=0.d0
256  stheta(1)=0.d0
257  IF(lst(17).NE.0) THEN
258 C...For varying beam energies, transform to cms, lepton along +z axis.
259  CALL ludbrb(0,0,0.,0.,0.d0,0.d0,-dbeta(1,3))
260  sphi(1)=ulangl(p(1,1),p(1,2))
261  CALL ludbrb(0,0,0.,-sphi(1),0.d0,0.d0,0.d0)
262  stheta(1)=ulangl(p(1,3),p(1,1))
263  CALL ludbrb(0,0,-stheta(1),0.,0.d0,0.d0,0.d0)
264  lst(28)=2
265  ENDIF
266 
267 C...Effective limits on kinematic variables x, y, Q**2, W**2
268  pm2=p(2,5)**2
269  s=parl(21)
270  xmin=max(cut(1),0.)
271  xmax=min(cut(2),1.)
272  ymin=max(cut(3),0.)
273  ymax=min(cut(4),1.)
274  q2min=max(cut(5),0.)
275  q2max=min(cut(6),s)
276  w2min=max(cut(7),0.)
277  w2max=min(cut(8),s+pm2)
278  umin=max(cut(9),0.)
279  umax=min(cut(10),s/(2.*p(2,5)))
280  DO 40 i=1,2
281  xmin=max(xmin,q2min/(s*ymax),q2min/(2.*p(2,5)*umax),
282  &1.-(w2max-pm2)/max(s*ymin,1.e-22),
283  &1.-(w2max-pm2)/max(2.*p(2,5)*umin,1.e-22))
284  xmax=min(xmax,q2max/max(s*ymin,1.e-22),
285  &q2max/max(2.*p(2,5)*umin,1.e-22),
286  &1.-(w2min-pm2)/(s*ymax),1.-(w2min-pm2)/(2.*p(2,5)*umax))
287  ymin=max(ymin,q2min/(s*xmax),(w2min-pm2)/(s*(1.-xmin)),
288  &(w2min-pm2+q2min)/s,2.*p(2,5)*umin/s)
289  ymax=min(ymax,q2max/max(s*xmin,1.e-22),
290  &(w2max-pm2)/max(s*(1.-xmax),1.e-22),
291  &(w2max-pm2+q2max)/s,2.*p(2,5)*umax/s)
292  q2min=max(q2min,s*xmin*ymin,s*ymin-w2max+pm2,
293  &2.*p(2,5)*umin*xmin,(w2min-pm2)*xmin/(1.-xmin))
294  q2max=min(q2max,s*xmax*ymax,s*ymax-w2min+pm2,
295  &2.*p(2,5)*umax*xmax,(w2max-pm2)*xmax/max(1.-xmax,1.e-22))
296  w2min=max(w2min,s*(1.-xmax)*ymin+pm2,q2min*(1.-xmax)/xmax+pm2,
297  &s*ymin-q2max+pm2,2.*p(2,5)*umin*(1.-xmax)+pm2)
298  w2max=min(w2max,s*(1.-xmin)*ymax+pm2,
299  &q2max*(1.-xmin)/max(xmin,1.e-22)+pm2,
300  &s*ymax-q2min+pm2,2.*p(2,5)*umax*(1.-xmin)+pm2)
301 C UMIN=MAX(UMIN,....)
302 C UMAX=MIN(UMAX,....)
303  40 CONTINUE
304  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,1050)
305  &cut,xmin,xmax,ymin,ymax,q2min,q2max,w2min,w2max,umin,umax
306  IF(xmax.LT.xmin.OR.ymax.LT.ymin.OR.q2max.LT.q2min.OR.
307  &w2max.LT.w2min) THEN
308  IF(lst(3).GE.1) WRITE(6,1100)
309  IF(lst(3).GE.2) THEN
310  WRITE(6,1900)
311  stop
312  ENDIF
313  ENDIF
314  IF(xmin.LT.1.e-10.OR.q2min.LT.1.e-01) THEN
315  IF(lst(3).GE.1) WRITE(6,1110)
316  IF(lst(3).GE.2) THEN
317  WRITE(6,1900)
318  stop
319  ENDIF
320  ENDIF
321 
322  pari(11)=(parl(1)-parl(2))/parl(1)
323  ksave(4)=lepin
324  ilep=1
325  IF(lepin.LT.0) ilep=2
326  inu=0
327  IF(iabs(lepin).EQ.12.OR.iabs(lepin).EQ.14
328  &.OR.iabs(lepin).EQ.16) inu=1
329  IF(inu.EQ.1) THEN
330 C...Set full polarisation for incoming neutrino.
331  parl(6)=-1.
332  IF(lepin.LT.0) parl(6)=1.
333  ENDIF
334  IF(lst(23).EQ.1.AND.inu.EQ.0) THEN
335 C...Electromagnetic interaction.
336  ksave(3)=22
337  ig=1
338  iz=0
339  ELSEIF(lst(23).EQ.2) THEN
340 C...Weak charged current, only one helicity state contributes.
341  IF(ksave(1).LT.0.AND.parl(6).LT.-0.99
342  & .OR.ksave(1).GT.0.AND.parl(6).GT.0.99) THEN
343  IF(lst(3).GE.1) WRITE(6,1150) lepin,parl(6)
344  IF(lst(3).GE.2) THEN
345  WRITE(6,1900)
346  stop
347  ENDIF
348  ENDIF
349  IF(mod(iabs(lepin),2).EQ.0) THEN
350  ksave(3)=isign(24,lepin)
351  ksave(4)=isign(iabs(lepin)-1,lepin)
352  ELSE
353  ksave(3)=isign(24,-lepin)
354  ksave(4)=isign(iabs(lepin)+1,lepin)
355  ENDIF
356  ELSEIF(lst(23).EQ.3.OR.(lst(23).EQ.4.AND.inu.EQ.1)) THEN
357 C...Weak neutral current.
358  ksave(3)=23
359  ig=0
360  iz=1
361  ELSEIF(lst(23).EQ.4.AND.inu.EQ.0) THEN
362 C...Neutral current, electromagnetic and weak with interference.
363  ksave(3)=23
364  ig=1
365  iz=1
366  ELSE
367  IF(lst(3).GE.1) WRITE(6,1200) inter,lepin
368  IF(lst(3).GE.2) THEN
369  WRITE(6,1900)
370  stop
371  ENDIF
372  ENDIF
373 
374 C...Choice of independent variables.
375  IF(lst(1).EQ.0) THEN
376  lst(31)=1
377  IF(inter.EQ.2.OR.inter.EQ.3) lst(31)=2
378  ELSE
379  lst(31)=iabs(lst(1))
380  ENDIF
381  IF(lst(31).LT.1.OR.lst(31).GT.3) THEN
382  IF(lst(3).GE.1) WRITE(6,1210) lst(1),lst(31)
383  IF(lst(3).GE.2) THEN
384  WRITE(6,1900)
385  stop
386  ENDIF
387  ENDIF
388  IF(lst(1).LT.0) THEN
389 C...User-defined optimization parameters.
390  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
391  & WRITE(6,1220) optx,opty,optq2,optw2
392  ELSE
393 C...Set optimization parameters.
394  DO 50 i=1,4
395  optx(i)=0.
396  opty(i)=0.
397  optq2(i)=0.
398  50 optw2(i)=0.
399  IF(inter.EQ.1) THEN
400  optx(2)=1.
401  opty(1)=1.
402  optq2(3)=1.
403  optw2(3)=1.
404  ELSEIF(inter.EQ.4) THEN
405  optx(1)=0.1
406  optx(2)=1.
407  opty(1)=1.
408  optq2(1)=0.5
409  optq2(2)=0.5
410  optq2(3)=1.
411  optw2(1)=0.5
412  optw2(2)=0.5
413  optw2(3)=1.
414  ELSE
415  optx(1)=1.
416  opty(1)=1.
417  optq2(1)=1.
418  optw2(1)=1.
419  ENDIF
420  ENDIF
421 
422 C...Initialize Monte Carlo estimate of cross section.
423  parl(24)=0.
424  pari(27)=0.
425  pari(28)=0.
426  pari(29)=0.
427  pari(30)=0.
428  pari(32)=0.
429  IF(lst(23).EQ.2) THEN
430 C...Constant factor GF**2/pi for CC, transformation to picobarn.
431  pari(31)=parl(17)**2/pi*0.39e+09
432  ELSE
433 C...Constant factor 2*pi*alpha**2 for NC, transformation to picobarn.
434  pari(31)=2.*pi*parl(16)**2*0.39e+09
435  ENDIF
436  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
437  &WRITE(6,1250) (i,lst(i),lst(i+10),parl(i),parl(i+10),i=1,10)
438 
439 C...Set up grid with longitudinal structure function, QCD & target mass;
440 C...only when photon exchange is included
441  lqcd=mod(lst(11),10)
442  ltm=mod(lst(11)/10,10)
443  IF(lst(11).NE.0.AND.(inter.EQ.1.OR.inter.EQ.4)) CALL fltabl
444 
445 C...Get integrated cross-section.
446  parl(23)=0.
447  IF(lst(10).GT.0) CALL lxsect
448  IF(lqcd.EQ.2.OR.ltm.EQ.2) THEN
449  WRITE(6,1300)
450  IF(lqcd.EQ.2) WRITE(6,1310)
451  IF(ltm .EQ.2) WRITE(6,1320)
452  WRITE(6,1330)
453  ENDIF
454 
455  IF(lst(2).EQ.1) THEN
456 C...Find max value of differential cross section for rejection.
457  ukin(1)=(xmax+xmin)/2.
458  wkin(1)=0.8*(xmax-xmin)/2.
459  ain(1)=xmin
460  bin(1)=xmax
461  IF(lst(31).EQ.1) THEN
462  ukin(2)=(q2max+q2min)/2.
463  wkin(2)=0.8*(q2max-q2min)/2.
464  ain(2)=q2min
465  bin(2)=q2max
466  namkin(2)=' Q**2'
467  ELSEIF(lst(31).EQ.2) THEN
468  ukin(2)=(ymax+ymin)/2.
469  wkin(2)=0.8*(ymax-ymin)/2.
470  ain(2)=ymin
471  bin(2)=ymax
472  namkin(2)=' y'
473  ELSEIF(lst(31).EQ.3) THEN
474  ukin(2)=(w2max+w2min)/2.
475  wkin(2)=0.8*(w2max-w2min)/2.
476  ain(2)=w2min
477  bin(2)=w2max
478  namkin(2)=' W**2'
479  ENDIF
480 C...Maximum obtained by minimizing -(diff. x-section).
481  CALL ltimex(ti1)
482  CALL lminew
483  CALL ltimex(ti2)
484  pari(lst(23))=fcnmax*1.1
485 CMM.. Maximum inserted by hand if MINUIT fails
486  IF(pari(lst(23)).EQ.0.0) pari(lst(23))=fmaxfh
487  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
488  & WRITE(6,1400) pari(lst(23)),ti2-ti1
489  ENDIF
490 
491  IF(lfile.GT.0.AND.lst(19).GE.0) THEN
492 C...Read QCD weights from file.
493  READ(lfile) lstw,parlw,nxx,nww,np,xx,ww
494  ipmax=2
495  IF(lstw(17).NE.0) ipmax=3
496  READ(lfile) (((pqg(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,np),
497  & (((pqqb(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,np),
498  & (((qgmax(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,ipmax),
499  & (((qqbmax(ix,iw,ip),ix=1,nxx),iw=1,nww),ip=1,min(2,ipmax)),
500  & ycut
501  IF(np.NE.1) READ(lfile) xtot
502  CLOSE(lfile)
503 C...Reset parameters for matrix element integration.
504  parl(8)=parlw(8)
505  parl(9)=parlw(9)
506  parl(11)=parlw(11)
507  parl(12)=parlw(12)
508  parl(13)=parlw(13)
509 C...Check current parameter values against those used when
510 C...calculating weights.
511 *...Added check also for polarization flag
512  IF(lst(12).NE.lstw(12).OR.lst(13).NE.lstw(13)
513  & .OR.lst(15).NE.lstw(15).OR.lst(16).NE.lstw(16)
514  & .OR.lst(17).NE.lstw(17).OR.lst(23).NE.lstw(23)
515  & .OR.lst(40).NE.lstw(40)
516  & .OR.abs(parl(1)-parlw(1)).GT.0.1.OR.abs(parl(2)-parlw(2)).GT.0.1
517  & .OR.abs(parl(5)-parlw(5)).GT.0.01
518  & .OR.abs(parl(6)-parlw(6)).GT.0.1) THEN
519  IF(lst(3).GE.1)
520  & WRITE(6,1500) lst(12),lstw(12),lst(13),lstw(13),lst(15),
521  & lstw(15),lst(16),lstw(16),lst(17),lstw(17),lst(23),lstw(23),
522  & lst(40),lstw(40),
523  & parl(1),parlw(1),parl(2),parlw(2),parl(5),parlw(5),parl(6),
524  & parlw(6)
525 *---
526  IF(lst(3).GE.2) THEN
527  WRITE(6,1900)
528  stop
529  ENDIF
530  ENDIF
531  ELSEIF((lst(19).GE.0.OR.lst(19).EQ.-10).AND.
532  &(lst(8).EQ.1.OR.lst(8)/10.EQ.1.OR.mod(lst(8),10).EQ.9)) THEN
533 C...Calculate weights if 1st order QCD from grid is requested.
534  CALL ltimex(ti1)
535  CALL lweits(lfile)
536  CALL ltimex(ti2)
537  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
538  & WRITE(6,1510) ti2-ti1
539  ENDIF
540 
541 C...Reset counters to zero for Monte Carlo estimate of cross section.
542  pari(27)=0.
543  pari(28)=0.
544  pari(29)=0.
545  pari(30)=0.
546  lst(32)=0
547  RETURN
548 
549  1000 FORMAT(' ',68('='),//,
550  &'A MONTE CARLO GENERATOR FOR DEEP INELASTIC LEPTON-'
551  &,'NUCLEON SCATTERING',//,
552  &'LEPTO version 6.5, April 20, 1996',/,68('='),//,
553  &'Lepton: type =',i3,5x,'momentum (px,py,pz) =',3f8.1,
554  &' GeV',//,'Target: A, Z =',2f3.0,2x,
555  &'momentum (px,py,pz) =',3f8.1,' GeV',//,
556  &'Interaction :',i3,14x,' CMS energy =',1pg12.4,' GeV',/68('='))
557 *...Write also relative polarization state
558  1005 FORMAT(' Realtive e-p polarization state:',i3,/,
559  &' Polarized structure function set:',i3,/)
560 *---
561  1010 FORMAT(' Warning: lepton and nucleon momenta in same direction',
562  &' not allowed.',/,10x,'Execution stopped.')
563 *...Stop execution if inconsistency is detected in polarized case
564  1020 FORMAT(/,' JETSET version ',i3,'.',i3,' is used.',/,
565  &' Parton densities in PYTHIA version ',i3,'.',i3,' are used.',/)
566  1025 FORMAT(' Warning: unallowed type of interaction.',/,10x,
567  &' Execution stopped.')
568  1030 FORMAT(' Warning (LINIT): JETSET version before 7.402, MSTJ(46)'
569  & ,' set to',i4,/,18x,'to avoid mismatch LEPTO<-->LUSHOW.',/)
570 *---
571  1050 FORMAT(/,' User applied cuts (+ phase space) : ',1p,
572  & g12.4,' < x < ',g12.4,
573  &/,37x,g12.4,' < y < ',g12.4,
574  &/,37x,g12.4,' < Q**2 < ',g12.4,
575  &/,37x,g12.4,' < W**2 < ',g12.4,
576  &/,37x,g12.4,' < nu < ',g12.4,
577  &/,37x,g12.4,' < E'' < ',g12.4,
578  &/,37x,g12.4,' < theta < ',g12.4,/,
579  &/, ' Effective ranges (from above cuts): ',
580  & g12.4,' < x < ',g12.4,
581  &/,37x,g12.4,' < y < ',g12.4,
582  &/,37x,g12.4,' < Q**2 < ',g12.4,
583  &/,37x,g12.4,' < W**2 < ',g12.4,
584  &/,37x,g12.4,' < nu < ',g12.4)
585  1100 FORMAT(' Warning: effective upper limit of kinematical ',
586  &'variable(s) smaller than corresponding lower limit.')
587  1110 FORMAT(' Warning: lower limit in x and/or Q2 too small for ',
588  &'DIS formalism.')
589  1150 FORMAT(' Warning: weak charged current cross section zero for ',
590  &'specified lepton helicity; LEPIN, PARL(6) =',i3,f5.2)
591  1200 FORMAT(' Warning: unrecognized interaction in LINIT call: ',
592  &'INTER = ',i5,' for lepton LEPIN =',i5)
593  1210 FORMAT(' Warning: unallowed value of LST(1) =',i3,
594  &' and/or LST(31) =',i3)
595  1220 FORMAT(/,' User-defined optimization parameters:',
596  &/,5x,'OPTX(1...4) =',4g11.3,/,5x,'OPTY(1...4) =',4g11.3,
597  &/,5x,'OPYQ2(1...4) =',4g11.3,/,5x,'OPTW2(1...4) =',4g11.3,/)
598  1250 FORMAT(/,' Parameter values:',
599  &'LST(I+10)',8x,'PARL(I)',5x,'PARL(I+10)',1p,
600  &/,5x,55('-'),10(/,3i10,2g15.4),/)
601  1300 FORMAT(' Warning: cross section, PARL(23), excludes FL (see ',
602  &'LST(11)) from:')
603  1310 FORMAT(10x,'QCD, since evaluated event by event for LQCD=2')
604  1320 FORMAT(10x,'TM , since evaluated event by event for LTM =2')
605  1330 FORMAT(' Cross section in PARL(24) includes these contributions.')
606  1400 FORMAT(' Max of differential cross section (for weighting) =',
607  &e12.4,/,' obtained in ',f7.2,' seconds.',/)
608 *...Added also LST(40)
609  1500 FORMAT(
610  &'with those used when calculating QCD weights.',
611  &'current value value for weights',/,
612  &/,' LST(12) ',i12,10x,i12,
613  &/,' LST(13) ',i12,10x,i12,
614  &/,' LST(15) ',i12,10x,i12,
615  &/,' LST(16) ',i12,10x,i12,
616  &/,' LST(17) ',i12,10x,i12,
617  &/,' LST(23) ',i12,10x,i12,
618  &/,' LST(40) ',i12,10x,i12,
619  &/,' PARL(1) ',e12.4,10x,e12.4,
620  &/,' PARL(2) ',e12.4,10x,e12.4,
621  &/,' PARL(5) ',e12.4,10x,e12.4,
622  &/,' PARL(6) ',e12.4,10x,e12.4)
623 *---
624  1510 FORMAT(/,' Time for calculating QCD weights =',f5.1,' seconds',/)
625  1900 FORMAT(' Execution stopped ',/)
626  END