Removing warnings (Sun)
[u/mrichter/AliRoot.git] / TPHIC / tphic17.f
CommitLineData
0b5dd071 1CDECK ID>, GGCDES.
2*=======================================================================
3CDECK ID>, GGINIT.
4 SUBROUTINE gginit
5*-----------------------------------------------------------------------
6* GGHIC initialization
7* Author: Yu.Kharlov
8*-----------------------------------------------------------------------
9 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
10 & amin,amax,ymin,ymax,nmas,ny, kferm,
11 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
12 & thetamin, costhv1, kv1,kv2,gvpar(4)
13 CHARACTER lumfil*80
14 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
15 & gvconst(4,10)
16 REAL nc
17 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
18 & xstote, ssbr(10)
19 COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
20 & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
21 COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
22 & xmll, xmlr, xmnl, xtanb, xmha, xmu,
23 & xmt, xat, xmbr, xab, u11, v11
24 COMMON/D2LParam/lz,la,Rion,Gamma
25 INTEGER lz,la
26 REAL Rion,Gamma
27 COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
28 DOUBLE PRECISION CKIN
29 SAVE /PYSUBS/
30 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
31 DOUBLE PRECISION PARP,PARI
32 SAVE /PYPARS/
33 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
34 DOUBLE PRECISION P,V
35 SAVE /PYJETS/
36 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
37 DOUBLE PRECISION PARU,PARJ
38 SAVE /PYDAT1/
39 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
40 DOUBLE PRECISION PMAS,PARF,VCKM
41 SAVE /PYDAT2/
42 COMMON/PYDAT4/CHAF(500,2)
43 CHARACTER CHAF*16
44 SAVE /PYDAT4/
45C MXSS = maximum number of modes
46C NSSMOD = number of modes
47C ISSMOD = initial particle
48C JSSMOD = final particles
49C GSSMOD = width
50C BSSMOD = branching ratio
51 INTEGER MXSS
52 PARAMETER (MXSS=1000)
53 COMMON/SSMODE/NSSMOD,ISSMOD(MXSS),JSSMOD(5,MXSS),GSSMOD(MXSS)
54 $,BSSMOD(MXSS)
55 INTEGER NSSMOD,ISSMOD,JSSMOD
56 REAL GSSMOD,BSSMOD
57 SAVE /SSMODE/
58C SUSY parameters
59C AMGLSS = gluino mass
60C AMQKSS = squark mass
61C AMLLSS = left-slepton mass
62C AMLRSS = right-slepton mass
63C AMNLSS = sneutrino mass
64C TWOM1 = Higgsino mass = - mu
65C RV2V1 = ratio v2/v1 of vev's
66C AMTLSS,AMTRSS = left,right stop masses
67C AMT1SS,AMT2SS = light,heavy stop masses
68C AMBLSS,AMBRSS = left,right sbottom masses
69C AMB1SS,AMB2SS = light,heavy sbottom masses
70C AMZiSS = signed mass of Zi
71C ZMIXSS = Zi mixing matrix
72C AMWiSS = signed Wi mass
73C GAMMAL,GAMMAR = Wi left, right mixing angles
74C AMHL,AMHH,AMHA = neutral Higgs h0, H0, A0 masses
75C AMHC = charged Higgs H+ mass
76C ALFAH = Higgs mixing angle
77C AAT = stop trilinear term
78C THETAT = stop mixing angle
79C AAB = sbottom trilinear term
80C THETAB = sbottom mixing angle
81 COMMON/SSPAR/AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS
82 $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS
83 $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS
84 $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS
85 $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS(4,4)
86 $,AMW1SS,AMW2SS
87 $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
88 $,AAB,THETAB,AAL,THETAL,AMGVSS
89 REAL AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS
90 $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS
91 $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS
92 $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS
93 $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS
94 $,AMW1SS,AMW2SS
95 $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
96 $,AAB,THETAB,AAL,THETAL,AMGVSS
97 REAL AMZISS(4)
98 EQUIVALENCE (AMZISS(1),AMZ1SS)
99 SAVE /SSPAR/
82764a9c 100 REAL xpar(10),ypar(6)
0b5dd071 101 INTEGER ipar(2)
102 INTEGER pycomp
103 REAL*8 ggrnd
104 EXTERNAL pydata
105 la = ia
106 lz = iz
107* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108*
109* Check bounds on 2 photon mass
110*
111 IF (iproc .EQ. 1) THEN
112 IF (amin .LE. 2.31) THEN
113 WRITE (6,*)
114 & ' %GGINIT: minimal 2-photon mass is too low, ',
115 & 'set to 3m(pho)'
116 amin = 2.31
117 END IF
118*
119 ELSE IF (iproc .EQ. 2) THEN
120 kc = pycomp (kf_onium)
121 xmres = pmas(kc,1)
122 IF (kf_onium .EQ. 441) THEN
123 xgtres = 0.013
124 ELSE
125 xgtres = pmas(kc,2)
126 END IF
127c. xggres = ggwid(kf_onium)
128 IF (xgtres .EQ. 0.) xgtres = 0.1
129 IF (xmres.LE.amin .OR. xmres.GE.amax) THEN
130 WRITE(6,*)
131 & ' %GGINIT: Xmres is out of the Luminosity mass range'
132 STOP
133 ENDIF
134*
135 ELSE IF (iproc .EQ. 3) THEN
136 kc = pycomp (kferm)
137 amprt(1) = pmas(kc,1)
138 IF (amin .LT. 2*amprt(1)) THEN
139 WRITE (6,*)
140 & ' %GGINIT: minimal 2-photon mass is too low, ',
141 & 'set to 2m(fermion)'
142 amin = 2*amprt(1)
143 END IF
144*
145 ELSE IF (iproc .EQ. 4) THEN
146 IF (amprt(4) .LE. 0.) amprt(4) = pmas(pycomp(24),1)
147 IF (amin .LT. 2*amprt(4)) then
148 WRITE (6,*)
149 & ' %GGINIT: minimal 2-photon mass is too low, ',
150 & 'set to 2m(W+-)'
151 amin = 2*amprt(4)
152 END IF
153*
154 ELSE IF (iproc .EQ. 5) THEN
155 WRITE (6,*) 'Process 5 is not implemented yet'
156 STOP
157C CALL ssmssm (xm1, xm2, xmg, xms, xmtl, xmtr, xmll, xmlr, xmnl,
158C & xtanb, xmha, xmu, xmt, xat, xmbr, xab, iallow)
159C IF (iallow .NE. 0) THEN
160C WRITE (6,*) 'Lightest neutralino is not LSP, ',
161C & 'input other MSSM parameters'
162C STOP
163C END IF
164C amprt(1) = abs (amw1ss)
165C amprt(2) = abs (amz1ss)
166C WRITE (6,'('' %GGINIT'','//
167C & ''', M(chi+) = '',f5.1,'//
168C & ''', M(chi0) = '',f5.1)') amprt(1), amprt(2)
169C *
170C IF (amin .LT. 2*amprt(1)) THEN
171C WRITE (6,*)
172C & ' %GGINIT: minimal 2-photon mass is too low, ',
173C & 'set to 2M(chi+)'
174C amin = 2*amprt(1)
175C END IF
176*
177 ELSE IF (iproc .EQ. 6) THEN
178 kc1 = pycomp (kv1)
179 amprt(1) = pmas(kc1,1) + pmas(kc1,3)
180 kc2 = pycomp (kv2)
181 amprt(2) = pmas(kc2,1) + pmas(kc2,3)
182 IF (amin .LT. amprt(1)+amprt(2)) THEN
183 WRITE (6,*)
184 & ' %GGINIT: minimal 2-photon mass is too low, ',
185 & 'set to m(kv1)+m(kv2)'
186 amin = amprt(1)+amprt(2)
187 END IF
188*
189 END IF
190*
191 IF (amin .GT. amax) THEN
192 WRITE (6,*) ' %GGINIT: AMIN > AMAX, execution stopped'
193 STOP
194 END IF
195*
196* ----- TWOGAM Initialization -----
197*
198 gamma = ebmn*ia/amas
199 xpar(1)= ebmn*ia ! Beam energy
200 xpar(2)= amin ! Minimum gg mass
201 xpar(3)= amax ! Maximum gg mass
202 xpar(8)= 100. ! Weight, but not too long
203*
204 ipar(1)= 1 ! Unweighted events
205 ipar(2)= 0 ! Continuum generation
206*
207 IF (iproc .EQ. 2) THEN
208 xpar(9) = xmres ! Resonance mass
209 xpar(10)= xgtres ! Resonance total width
210 ipar(2) = 1 ! Resonance generation
211 END IF
212*
213 ypar(1)= iz ! Ion charge
214 ypar(2)= ia ! Ion atomic number
215 ypar(3)= amas ! Ion mass
216 ypar(4)= gamma ! Gamma factor of ion
217 ypar(5)= ymin ! Minimum gg-rapidity
218 ypar(6)= ymax ! Maximum gg-rapidity
219 eb = xpar(1)
220*
221 CALL twoini (6,ipar,xpar,ypar)
222*
223* The Luminosity function initialization
224*
225* ilumf = +1 ! Reading the Luminosity from file
226* ilumf = -1 ! New calculation of the Luminosity function
227*
228 CALL lumfun (ilumf,lumfil,ymin,ymax,ny,amin,amax,nmas)
229*
230* Integral of lum. function over the (M,y) region
231 xlumint = 0.0
232 hmas = (amax-amin)/nmas
233 f1 = dldmx(amin ,ny)
234 DO i=1,nmas
235 xm = amin + i*hmas
236 f2 = dldmx(xm-hmas/2.,ny)
237 f3 = dldmx(xm ,ny)
238 xlumint = xlumint +(f1 + 4.*f2 + f3)/6.
239 f1 = f3
240 END DO
241 xlumint = xlumint * hmas
242* ----- end of TWOGAM initialization -----
243*
244 IF (iproc .EQ. 1) THEN
245* ----- PYTHIA initialization -----
246 parp(2) = dble(amin)
247 mstp(14) = 20
248 mstp(82) = 1
249 mstp(171) = 1
250 mstp(172) = 1
251 DO ip = 1,2
252 p(1,ip) = 0.D0
253 p(2,ip) = 0.D0
254 END DO
255 p(1,3) = dble( amax/2)
256 p(2,3) = dble(-amax/2)
257 p(1,4) = dble( amax/2)
258 p(2,4) = dble( amax/2)
259 p(1,5) = 0.D0
260 p(2,5) = 0.D0
261 CALL pyinit ('5MOM','gamma','gamma',0.0D0)
262* ----- end of PYTHIA initialization -----
263*
264 ELSE IF (iproc .EQ. 2) THEN
265* --- Cross section for narrow Resonance production ---
266 totxsc = dldmx(xmres,ny) ! dL/dMx in GeV^-1
267 njres = kf_onium / 10
268 njres = kf_onium - 10 * njres
269 xnorm = 4. * pi**2 * njres * xggres * gev2nb / xmres**2
270 xstot = xnorm * totxsc
271 xstote = 0.
272*
273 ELSE IF (iproc .EQ. 3) THEN
274 qf = abs(kchg(kc,1))/3.
275 nc = abs(kchg(kc,2))*2. + 1.
276 ELSE IF (iproc .EQ. 5) THEN
277* Chargino storing in JETSET
278 kc = pycomp (41)
279 chaf(kc,1) = 'chi_1'
280 kchg(kc,1) = -3
281 kchg(kc,2) = 0
282 kchg(kc,3) = 1
283 pmas(kc,1) = amprt(1)
284 qf = 1.
285 nc = 1.
286* Neutralino storing in JETSET
287 kc = pycomp (42)
288 chaf(kc,1) = 'chi_1'
289 kchg(kc,1) = 0
290 kchg(kc,2) = 0
291 kchg(kc,3) = 1
292 pmas(kc,1) = amprt(2)
293*
294* Gaugino mixing parameters
295 u11 = -zmixss(4,1)*sin(gammar)/sqrt(2.) +
296 & zmixss(2,1)*cos(gammar)
297 v11 = zmixss(3,1)*sin(gammal)/sqrt(2.) +
298 & zmixss(2,1)*cos(gammal)
299 DO idc = 1,nssmod
300 IF (issmod(idc) .EQ. 39) THEN
301 IF (moddcy .EQ. 1) THEN
302 DO jmod = 1,5
303 IF (jssmod(jmod,idc) .EQ. -12) THEN
304 xsbra = bssmod(idc)
305 xsbra = xsbra*xsbra*4
306 GOTO 10
307 END IF
308 END DO
309 ELSE
310 DO jmod = 1,5
311 IF (jssmod(jmod,idc) .EQ. -12) THEN ! e nu_e
312 ssbr(1) = bssmod(idc)
313 GOTO 20
314 ELSE IF (jssmod(jmod,idc) .EQ. -14) THEN ! mu nu_mu
315 ssbr(2) = bssmod(idc)
316 GOTO 20
317 ELSE IF (jssmod(jmod,idc) .EQ. -16) THEN ! tau nu_tau
318 ssbr(3) = bssmod(idc)
319 GOTO 20
320 ELSE IF (jssmod(jmod,idc) .EQ. -2) THEN ! u dbar
321 ssbr(4) = bssmod(idc)
322 GOTO 20
323 ELSE IF (jssmod(jmod,idc) .EQ. 4) THEN ! c sbar
324 ssbr(5) = bssmod(idc)
325 GOTO 20
326 END IF
327 END DO
328 20 CONTINUE
329 END IF
330 END IF
331 END DO
332*
333 10 CONTINUE
334 ELSE IF (iproc .EQ. 6) THEN
335* Fix parameters of process gg -> V1V2
336 IF (kv1 .GT. kv2) THEN
337 kvtmp = kv1
338 kv1 = kv2
339 kv2 = kvtmp
340 END IF
341 IF (kv1.EQ.113 .AND. kv2.EQ.113) THEN
342 idx_vv = 1
343 ELSE IF (kv1.EQ.113 .AND. kv2.EQ.223) THEN
344 idx_vv = 2
345 ELSE IF (kv1.EQ.113 .AND. kv2.EQ.333) THEN
346 idx_vv = 3
347 ELSE IF (kv1.EQ.113 .AND. kv2.EQ.443) THEN
348 idx_vv = 4
349 ELSE IF (kv1.EQ.223 .AND. kv2.EQ.223) THEN
350 idx_vv = 5
351 ELSE IF (kv1.EQ.223 .AND. kv2.EQ.333) THEN
352 idx_vv = 6
353 ELSE IF (kv1.EQ.223 .AND. kv2.EQ.443) THEN
354 idx_vv = 7
355 ELSE IF (kv1.EQ.333 .AND. kv2.EQ.333) THEN
356 idx_vv = 8
357 ELSE IF (kv1.EQ.333 .AND. kv2.EQ.443) THEN
358 idx_vv = 9
359 ELSE IF (kv1.EQ.443 .AND. kv2.EQ.443) THEN
360 idx_vv = 10
361 ELSE
362 WRITE (6,*) ' %GGINIT: unknown (V1 V2) combination'
363 WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
364 STOP
365 END IF
366 DO 101 igvpar=1,4
367 101 gvpar(igvpar) = gvconst(igvpar,idx_vv)
368 IF (gvpar(1) .EQ. 0.0) THEN
369 WRITE (6,*) ' %GGINIT: (V1 V2) combination not implemented'
370 WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
371 STOP
372 END IF
373 END IF
374*
375 IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR.
376 & iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN
377* Max cross section estimation
378 xsmax0 = 0.
379 DO nest = 1,100
380 xm = amin + ggrnd(0) * (amax-amin)
381 wsq = xm*xm
382 CALL ggxsec
383 END DO
384 xssum = 0.
385 ntry = 0
386 END IF
387*
388 RETURN
389 END
390*
391*=======================================================================
392CDECK ID>, GGWID.
393 FUNCTION ggwid (kf_res)
394*-----------------------------------------------------------------------
395* Routine to define 2-gamma width of a resonance
396* with JETSET code KF_RES
397* Author: Yu.Kharlov
398*-----------------------------------------------------------------------
399 REAL wid(5)
400 DATA wid /9.E-09, 1.E-06, 5.E-06, 6.3E-06, 0.41E-06/
401 IF (kf_res .EQ. 111) THEN
402 ggwid = wid(1)
403 ELSE IF (kf_res .EQ. 221) THEN
404 ggwid = wid(2)
405 ELSE IF (kf_res .EQ. 331) THEN
406 ggwid = wid(3)
407 ELSE IF (kf_res .EQ. 441
408 & .OR. kf_res .EQ. 10441
409 & .OR. kf_res .EQ. 445) THEN
410 ggwid = wid(4)
411 ELSE IF (kf_res .EQ. 551
412 & .OR. kf_res .EQ. 10551
413 & .OR. kf_res .EQ. 555) THEN
414 ggwid = wid(5)
415 END IF
416 RETURN
417 END
418*
419*=======================================================================
420CDECK ID>, GGRUN.
421 SUBROUTINE ggrun
422*-----------------------------------------------------------------------
423* Single event generation
424* Author: Yu.Kharlov
425*-----------------------------------------------------------------------
426 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
427 & amin,amax,ymin,ymax,nmas,ny, kferm,
428 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
429 & thetamin, costhv1, kv1,kv2,gvpar(4)
430 CHARACTER lumfil*80
431 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
432 & gvconst(4,10)
433 REAL nc
434 COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
435 & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
436 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
437 & xstote, ssbr(10)
438 COMMON /ggkin/ xkin(10)
439 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
440 DOUBLE PRECISION P,V
441 SAVE /PYJETS/
442 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
443 DOUBLE PRECISION PARU,PARJ
444 SAVE /PYDAT1/
445 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
446 DOUBLE PRECISION PMAS,PARF,VCKM
447 SAVE /PYDAT2/
82764a9c 448 REAL pgam1(4),pgam2(4)
0b5dd071 449 REAL p3(4), p4(4)
450 REAL*8 dbetaz
451 REAL*8 ggrnd
452*
453 1 CALL gg2gam (amas,wsq,ygg,q1sq,q2sq,
454 & pgam1,pgam2,ptag1,ptag2,weight)
455*
456* Cross section calculation
457 IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR.
458 & iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN
459 CALL ggxsec
460 xssum = xssum + xscur
461 ntry = ntry + 1
462* Select event on cross section of process gg -> XX
463 IF (xscur0/xsmax0 .LT. sngl (ggrnd(0))) GOTO 1
464* Total cross section calculation due to current statistics
465 xstot = xssum/ntry
466 xstote = xstot / sqrt(float(ntry))
467 END IF
468*
469 xmgg = sqrt (wsq)
470 xmg1 = sqrt (q1sq)
471 xmg2 = sqrt (q2sq)
472 dbetaz = dtanh(1.d0*ygg)
473*
474* calculation of the 4-vector of two-photon system
475* and filling arrays pgg(1:2,1:5), kgg(1:2)
476 DO ip=1,4
477 p2g(ip) = pgam1(ip) + pgam2(ip)
478 pgg(1,ip) = pgam1(ip)
479 pgg(2,ip) = pgam2(ip)
480 END DO
481 p2g(5) = xmgg
482 pgg(1,5) = -xmg1
483 pgg(2,5) = -xmg2
484 kgg(1) = 22
485 kgg(2) = 22
486 ngg = 2
487*
488* =======> PROCESS # 1
489 IF (iproc .EQ. 1) THEN
490 DO ip = 1,4
491 p(1,ip) = dble (pgam1(ip))
492 p(2,ip) = dble (pgam2(ip))
493 END DO
494 p(1,5) = dble (-xmg1)
495 p(2,5) = dble (-xmg2)
496 CALL pyevnt
497*
498* =======> PROCESS # 2
499 ELSE IF (iproc .EQ. 2) THEN
500 ngg = ngg + 1
501 DO ip = 1,5
502 pgg(ngg,ip) = p2g(ip)
503 END DO
504 kgg(ngg) = kf_onium
505* Fill the JETSET common block
506 CALL gg2py (ngg)
507 CALL pyexec
508*
509* =======> PROCESS # 3, # 4, # 5, # 6
510 ELSE IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR.
511 & iproc.EQ.5 .OR. iproc.EQ.6) THEN
512 phi = 2*pi*ggrnd(0)
513 IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR. iproc.EQ.5) THEN
514 e3 = 0.5*sqrt(wsq)
515 e4 = e3
516 beta = xkin(1)
517 costhe = xkin(2)
518 am1 = xkin(3)
519 am2 = xkin(3)
520 p3abs = e3 * beta
521 ELSE IF (iproc.EQ.6) THEN
522 am1 = xkin(1)
523 am2 = xkin(2)
524 costhe = xkin(3)
525 p3abs = xkin(4)
526 e3 = xkin(5)
527 e4 = xmgg - e3
528 END IF
529 sinthe = sqrt (1 - costhe*costhe)
530 p3(1) = p3abs * cos(phi) * sinthe
531 p3(2) = p3abs * sin(phi) * sinthe
532 p3(3) = p3abs * costhe
533 p3(4) = e3
534 p4(4) = e4
535 DO ip = 1,3
536 p4(ip) = -p3(ip)
537 END DO
538 igg1 = ngg + 1
539 igg2 = ngg + 2
540 DO ip = 1,4
541 pgg(igg1,ip) = p3(ip)
542 pgg(igg2,ip) = p4(ip)
543 END DO
544 pgg(igg1,5) = am1
545 pgg(igg2,5) = am2
546 ngg = ngg + 2
547*
548 IF (iproc .EQ. 3) THEN
549 kgg(igg1) = kferm
550 kgg(igg2) = -kferm
551 CALL gg2py (ngg)
552 ELSE IF (iproc .EQ. 4) THEN
553 kgg(igg1) = 24
554 kgg(igg2) = -24
555 CALL gg2py (ngg)
556 ELSE IF (iproc .EQ. 5) THEN
557 kgg(igg1) = 41
558 kgg(igg2) = -41
559* Decay of charginos into W+neutralino
560 CALL ggdecy (igg1)
561 CALL ggdecy (igg2)
562*
563 CALL gg2py (ngg)
564 k(igg1,1) = 11
565 k(igg2,1) = 11
566 ELSE IF (iproc.EQ.6) THEN
567 kgg(igg1) = kv1
568 kgg(igg2) = kv2
569 CALL gg2py (ngg)
570 END IF
571*
572* Boost event into initial frame
573 mstu(33) = 1
574 CALL pyrobo(3,ngg, 0.d0,0.d0, 0.d0,0.d0, dbetaz)
575 CALL pyexec
576*
577 END IF
578*
579 RETURN
580 END
581*
582*=======================================================================
583CDECK ID>, GGEXIT.
584 SUBROUTINE ggexit
585*-----------------------------------------------------------------------
586* Cross section calculation and end of run and
587* print out cross section and related variables
588* Author: Yu.Kharlov
589*-----------------------------------------------------------------------
590 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
591 & amin,amax,ymin,ymax,nmas,ny, kferm,
592 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
593 & thetamin, costhv1, kv1,kv2,gvpar(4)
594 CHARACTER lumfil*80
595 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
596 & xstote, ssbr(10)
597 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
598 & gvconst(4,10)
599 REAL nc
600*
601 WRITE (6,'(/x,79(''=''))')
602 WRITE (6,'(x,''I'',36x,''TPHIC'',36x, ''I'')')
603 WRITE (6,'(x,''I'',34x,''Process '',i1,34x,''I'')') iproc
604 WRITE (6,'(x,''I'',34x,''---------'', 34x,''I'')')
605 WRITE (6,'(x,''I'',77x, ''I'')')
606 WRITE (6,'(x,''I'',20x,''Events thrown : '',i8,31x,''I'')')
607 & ntry
608 WRITE (6,'(x,''I'',20x,''Events accepted : '',i8,31x,''I'')')
609 & nevent
610 WRITE (6,'(x,''I'',20x,''Cross section : '',e10.3,'' +- '','//
611 & 'e10.3,'' nb'',12x ''I'')') xstot, xstote
612 WRITE (6,'(x,''I'',77x, ''I'')')
613 WRITE (6,'(x,79(''=''))')
614*
615 IF (iproc .EQ. 1) THEN
616* --- Alternative cross section for minimum bias events ---
617 totxsc = 0.0
618 hmas = (amax-amin)/nmas
619 f1 = dldmx(amin ,ny) * dchad(amin)
620 DO i=1,nmas
621 xm = amin + i*hmas
622 f2 = dldmx(xm-hmas/2.,ny) * dchad(xm-hmas/2.)
623 f3 = dldmx(xm ,ny) * dchad(xm)
624 totxsc = totxsc +(f1 + 4.*f2 + f3)/6.
625 f1 = f3
626 END DO
627 sigtt = totxsc * hmas
628 sigmx = sigtt/sqrt(float(nevent))
629 WRITE (6,*) ' CrosSec =',sigtt,' +/-',sigmx,' mb'
630 END IF
631*
632 RETURN
633 END
634*
635 REAL FUNCTION dchad(xm)
636*-----------------------------------------------------------------------
637* Total cross section for gamma-gamma collisions vs. Xm in mb
638* Parametrization is given by G.A.Schuler, T.Sjostrand,
639* CERN-TH.7193/94, formula (9).
640* Coded by Yu.Kharlov 20.03.1995.
641*-----------------------------------------------------------------------
642 DATA epsilon, eta /0.0808, 0.4525/
643 dchad = 0.211e-3 * (xm**2)**epsilon + 0.297e-3/(xm**2)**eta
644 RETURN
645 END
646*
647*=======================================================================
648CDECK ID>, GGXSEC.
649 SUBROUTINE ggxsec
650*-----------------------------------------------------------------------
651* Cross section calculation in current event for iproc=1,3,4,5,6
652* Author: Yu.Kharlov
653*-----------------------------------------------------------------------
654 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
655 & amin,amax,ymin,ymax,nmas,ny, kferm,
656 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
657 & thetamin, costhv1, kv1,kv2,gvpar(4)
658 CHARACTER lumfil*80
659 COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
660 & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
661 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
662 & gvconst(4,10)
663 REAL nc
664 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
665 & xstote, ssbr(10)
666 COMMON /ggkin/ xkin(10)
667 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
668 DOUBLE PRECISION PMAS,PARF,VCKM
669 SAVE /PYDAT2/
670 PARAMETER (epsilon = 0.0808, eta = -0.4525)
671 REAL*8 ggrnd
672 REAL*8 pymass
673*
674 IF (iproc .EQ. 1) THEN
675* Parametrization by G.A.Schuler, T.Sjostrand,
676* CERN-TH.7193/94, formula (9).
677 xscur0 = 211. * wsq**epsilon + 297. * wsq**eta
678 xscur = xscur0 * xlumint
679 ELSE IF (iproc .EQ. 3 .OR. iproc .EQ. 5) THEN
680 beta2 = 1.0 - 4.0*amprt(1)**2/wsq
681 beta2 = min (beta2, 1.0)
682 beta = sqrt (beta2)
683 IF (beta2 .LT. 0.9) THEN
684 1 costhe = 2.*ggrnd(0) - 1.
685 costh2 = costhe*costhe
686 xnum = 1. + 2.*beta2 * (1.-beta2) * (1.-costh2) -
687 & (beta2*costh2)**2
688 xden = (1 - beta2*costh2)**2
689 ffcur = xnum / xden
690 ffmax = (1. + beta2) / (1. - beta2)
691 rff = ggrnd(0)
692 IF (ffcur .LE. rff*ffmax) GOTO 1
693 ELSE
694 small = 1.0 - cos(thetamin)
695 2 zz = exp (ggrnd(0) * alog(2./small))
696C sgn = sign (1., 2.*ggrnd(0)-1.)
697 sgn = 2.*ggrnd(0)-1.
698 sgn = sign (1., sgn)
699 costhe = sgn * (zz-1.) / (zz+1.)
700 ffcur = (1.+costhe*costhe) / (1.-beta2*costhe*costhe)
701 ffmax = 2. / (1.- costhe*costhe)
702 rff = ggrnd(0)
703 IF (ffcur .LE. rff*ffmax) GOTO 2
704 END IF
705*
706 xkin(1) = beta
707 xkin(2) = costhe
708 xkin(3) = amprt(1)
709*
710* Yu.Kharlov, 3.06.1995
711 totfac = 4. * pi * alpha**2 * gev2nb * qf**4 * nc / wsq
712 IF (beta2 .LT. 0.99) THEN
713 xscur0 =((3.-beta2**2)/(2.*beta) * alog ((1.+beta)/(1.-beta))-
714 & 2. + beta2) * beta
715 ELSE
716 xscur0 = 2. * alog (sqrt(wsq)/amprt(1)) - 1.0
717 END IF
718 xscur0 = totfac * xscur0 * xsbra
719 xscur = xscur0 * xlumint
720 ELSE IF (iproc .EQ. 4) THEN
721 xmw = amprt(4)
722 beta2 = 1. - 4.*xmw**2/wsq
723 beta = sqrt (beta2)
724 3 costhe = 2.*ggrnd(0) - 1.
725 costh2 = costhe*costhe
726 xnum = 19. - 6.*beta2 * (1.-beta2) +
727 & 2.*beta2 * (8. - 3.*beta2) * costh2 +
728 & 3.*(beta2*costh2)**2
729 xden = (1. - beta2*costh2)**2
730 ffcur = xnum / xden
731 ffmax = (19. + 10.*beta2 + beta2**2) / (1. - beta2)**2
732 rff = ggrnd(0)
733 IF (ffcur .LE. rff*ffmax) GOTO 3
734*
735 xkin(1) = beta
736 xkin(2) = costhe
737 xkin(3) = xmw
738*
739* Yu.Kharlov, 3.06.1995
740 totfac = pi * alpha**2 * beta * gev2nb/ wsq
741 xscur0 = 2. * (22. - 9.*beta2 + 3.*beta2**2) / (1. - beta2) -
742 & 3. * (1.-beta2**2)/beta * alog ((1.+beta)/(1.-beta))
743 xscur0 = totfac * xscur0 * xsbra
744 xscur = xscur0 * xlumint
745 ELSE IF (iproc .EQ. 6) THEN
746 xmgg = sqrt (wsq)
747 bb = gvpar(3) + gvpar(4)*alog(xmgg/50.0)
748 am1 = sngl (pymass(kv1))
749 am2 = sngl (pymass(kv2))
750 p1 = pfork (xmgg,am1,am2)
751 e1 = sqrt (p1*p1 + am1*am1)
752 tmin = am1*am1 - xmgg*(e1+p1*costhv1)
753 tmax = am1*am1 - xmgg*(e1-p1*costhv1)
754 fmin = exp (bb*tmin)
755 fmax = exp (bb*tmax)
756 IF (fmax .GT. 1.e-10) THEN
757 tvar = 1./bb * dlog (fmax - ggrnd(0)*(fmax-fmin))
758 ELSE
759 tvar = tmax + 1./bb * dlog (1. - ggrnd(0))
760 END IF
761 costhe = (e1*xmgg - am1*am1 + tvar) / xmgg / p1
762 IF (costhe .GT. 1.0) costhe = 1.0
763 IF (costhe .LT.-1.0) costhe = -1.0
764 IF (ggrnd(0) .LT. 0.5) costhe = -costhe
765*
766 xkin(1) = am1
767 xkin(2) = am2
768 xkin(3) = costhe
769 xkin(4) = p1
770 xkin(5) = e1
771*
772 IF (kv1 .EQ. 443 .AND. kv2 .EQ. 443) THEN
773 factor = 1.0 / alog(xmgg/3.097)
774 ELSE
775 factor = 1.0
776 END IF
777 xscur0 = factor * gvpar(1) * xmgg**gvpar(2) / bb * exp(bb*tmax)
778 xscur = xscur0 * xlumint
779*
780 END IF
781*
782* Find max cross section of gg-process
783*
784 IF (xscur0 .GT. xsmax0) xsmax0 = xscur0
785*
786 RETURN
787 END
788*
789*=======================================================================
790CDECK ID>, GGDECY.
791 SUBROUTINE ggdecy (igg)
792*-----------------------------------------------------------------------
793* Routine to decay a particle number I in GGEVNT common block
794* Input : IGG - particle number in common GGEVNT arrays
795* Author : Yu.Kharlov
796*-----------------------------------------------------------------------
797 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
798 & amin,amax,ymin,ymax,nmas,ny, kferm,
799 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
800 & thetamin, costhv1, kv1,kv2,gvpar(4)
801 CHARACTER lumfil*80
802 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
803 & gvconst(4,10)
804 REAL nc
805 COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
806 & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
807 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
808 & xstote, ssbr(10)
809 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
810 DOUBLE PRECISION PMAS,PARF,VCKM
811 SAVE /PYDAT2/
812 COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
813 & xmll, xmlr, xmnl, xtanb, xmha, xmu,
814 & xmt, xat, xmbr, xab, u11, v11
815 REAL p1(4), pcm(4), p2(4), p3(4), p4(4), ptmp(3)
816 REAL o1(3,3), o2(3,3)
817 INTEGER pycomp
818 DATA ifl /0/
819 REAL*8 ggrnd
820*
821 DO ip = 1,4
822 p1(ip) = pgg(igg,ip)
823 END DO
824*
825 IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 3) THEN
826* Decay chi_1^+- --> chi_1^0 W^+- in CMS, uniformly
827 eb = (amprt(1)**2 + amprt(2)**2 - amprt(4)**2) / (2 * amprt(1))
828 ec = (amprt(1)**2 - amprt(2)**2 + amprt(4)**2) / (2 * amprt(1))
829 pb = sqrt (eb*eb - amprt(2)*amprt(2))
830 pc = pb
831 costhe = 2*ggrnd(0) - 1.
832 sinthe = sqrt (1 - costhe*costhe)
833 phi = 2*pi*ggrnd(0)
834 pcm(1) = pb * sinthe * cos(phi)
835 pcm(2) = pb * sinthe * sin(phi)
836 pcm(3) = pb * costhe
837 pcm(4) = eb
838 CALL lorenb (amprt(1),p1,pcm,p2)
839 DO ip = 1,3
840 pcm(ip) = -pcm(ip)
841 END DO
842 pcm(ip) = ec
843 CALL lorenb (amprt(1),p1,pcm,p3)
844 DO ip = 1,4
845 pgg(ngg+1,ip) = p2(ip)
846 pgg(ngg+2,ip) = p3(ip)
847 END DO
848 pgg(ngg+1,5) = amprt(2)
849 pgg(ngg+2,5) = amprt(4)
850 kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
851 kgg(ngg+2) = 24 * sign(1.,1.*kgg(igg))
852 ngg = ngg + 2
853*
854 ELSE IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 5) THEN
855* Decay chi_1^+- --> e^+- nu_e chi_1^0 in CMS,
856* according to the exact matrix element
857*
858 xm23max = amprt(1) - amprt(2)
859 xmw = sngl (pmas(pycomp(24),1))
860 gwt = sngl (pmas(pycomp(24),2))
861 atmin = - atan (xmw/gwt)
862 atmax = atan ((xm23max**2 - xmw**2) / (xmw*gwt))
863 am1 = amprt(1)
864 am4 = amprt(2)
865*
866 1 CONTINUE
867*
868* Random mass of (2+3) according to Breit-Wigner
869*
870 xm23 = xmw**2 + xmw*gwt *
871 * tan (atmin + ggrnd(0)*(atmax - atmin))
872 xm23 = max (0., xm23)
873 xm23 = sqrt (xm23)
874 IF (xm23 .GT. xm23max) GOTO 1 ! to avoid boundary effect
875 e23 = (am1**2 + xm23**2 - am4**2) / 2. / am1
876*
877* Rest frame of (2+3)
878*
879 phi = 2.*pi * ggrnd(0)
880 costhe = 2.*ggrnd(0) - 1.
881 sinthe = sqrt (1. - costhe*costhe)
882*
883 pp = xm23 / 2.
884 p2(1) = pp * sinthe * cos(phi)
885 p2(2) = pp * sinthe * sin(phi)
886 p2(3) = pp * costhe
887 p2(4) = pp
888*
889 DO ip = 1,3
890 p3(ip) = -p2(ip)
891 END DO
892 p3(4) = pp
893*
894* Transformation from rest frame of (2+3) to rest frame of (1),
895* 3-momenta p2, p3 are along z-axis
896*
897 gamma = e23 / xm23
898 gambe = sqrt (gamma*gamma - 1.)
899*
900 p2z = gamma*p2(3) + gambe*p2(4)
901 e2 = gamma*p2(4) + gambe*p2(3)
902 p3z = gamma*p3(3) + gambe*p3(4)
903 e3 = gamma*p3(4) + gambe*p3(3)
904*
905 p2(3) = p2z
906 p2(4) = e2
907 p3(3) = p3z
908 p3(4) = e3
909*
910* Arbitrary rotation of rest frame of (1)
911*
912 phi = 2.*pi * ggrnd(0)
913 costhe = 2.*ggrnd(0) - 1.
914 sinthe = sqrt (1. - costhe*costhe)
915 cosphi = cos(phi)
916 sinphi = sin(phi)
917*
918 o1(1,1) = costhe
919 o1(1,2) = 0.
920 o1(1,3) = -sinthe
921 o1(2,1) = 0.
922 o1(2,2) = 1.
923 o1(2,3) = 0.
924 o1(3,1) = sinthe
925 o1(3,2) = 0.
926 o1(3,3) = costhe
927*
928 o2(1,1) = cosphi
929 o2(1,2) = -sinphi
930 o2(1,3) = 0.
931 o2(2,1) = sinphi
932 o2(2,2) = cosphi
933 o2(2,3) = 0.
934 o2(3,1) = 0.
935 o2(3,2) = 0.
936 o2(3,3) = 1.
937*
938 DO ix = 1,3
939 ptmp(ix) = 0.
940 DO jx = 1,3
941 ptmp(ix) = ptmp(ix) + o1(ix,jx) * p2(jx)
942 END DO
943 END DO
944 DO ix = 1,3
945 p2(ix) = ptmp(ix)
946 END DO
947 DO ix = 1,3
948 ptmp(ix) = 0.
949 DO jx = 1,3
950 ptmp(ix) = ptmp(ix) + o2(ix,jx) * p2(jx)
951 END DO
952 END DO
953 DO ix = 1,3
954 p2(ix) = ptmp(ix)
955 END DO
956*
957 DO ix = 1,3
958 ptmp(ix) = 0.
959 DO jx = 1,3
960 ptmp(ix) = ptmp(ix) + o1(ix,jx) * p3(jx)
961 END DO
962 END DO
963 DO ix = 1,3
964 p3(ix) = ptmp(ix)
965 END DO
966 DO ix = 1,3
967 ptmp(ix) = 0.
968 DO jx = 1,3
969 ptmp(ix) = ptmp(ix) + o2(ix,jx) * p3(jx)
970 END DO
971 END DO
972 DO ix = 1,3
973 p3(ix) = ptmp(ix)
974 END DO
975*
976 p12 = am1 * p2(4)
977*
978 e4 = am1 - e23
979 pp =sqrt (e4*e4 - am4*am4)
980 p4(1) = pp * sinthe * cos(phi)
981 p4(2) = pp * sinthe * sin(phi)
982 p4(3) = -pp * costhe
983 p4(4) = e4
984*
985 width = (- 4.*p12**2*u11**2 - 4.*p12**2*v11**2 +
986 & 4.*p12*xm23**2*u11**2 - 2.*p12*am4**2*u11**2 -
987 & 2.*p12*am4**2*v11**2 + 2.*p12*am1**2*u11**2 +
988 & 2.*p12*am1**2*v11**2 - xm23**4*u11**2 +
989 & xm23**2*am4**2*u11**2 - 2.*xm23**2*am4*am1*u11*v11 -
990 & xm23**2*am1**2*u11**2) /
991 & ((xm23**2-xmw**2)**2 + (xmw*gwt)**2)
992 width = width * pp * xm23
993*
994 widmax = (2.*am1*xm23max *
995 & (2.*(xm23max*u11)**2 + am1**2*(u11**2 + v11**2)) +
996 & (xm23max*am4*u11)**2) * (am1**2 - am4**2) / 2./am1 *
997 & xm23 / ((xm23**2-xmw**2)**2 + (xmw*gwt)**2)
998*
999* Select current phase space configuration
1000*
1001 IF (width .GT. widmax) WRITE (6,*)
1002 & '%GGDECY: Maximum chi+ width violated'
1003 rwid = widmax * ggrnd(0)
1004 IF (width .LE. rwid) GOTO 1
1005*
1006* Pick up final lepton flavour
1007 rflav = ggrnd(0)
1008 IF (moddcy .EQ. 1) THEN
1009 IF (rflav .LE. 0.5) THEN
1010 lflav1 = 11
1011 lflav2 = -12
1012 ELSE
1013 lflav1 = 13
1014 lflav2 = -14
1015 END IF
1016 ELSE
1017 rr1 = ssbr(1)
1018 rr2 = ssbr(1)+ssbr(2)
1019 rr3 = ssbr(1)+ssbr(2)+ssbr(3)
1020 rr4 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)
1021 rr5 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)+ssbr(5)
1022 IF (rflav .LE. rr1) THEN
1023 lflav1 = 11
1024 lflav2 = -12
1025 ELSE IF (rflav .LE. rr2) THEN
1026 lflav1 = 13
1027 lflav2 = -14
1028 ELSE IF (rflav .LE. rr3) THEN
1029 lflav1 = 15
1030 lflav2 = -16
1031 ELSE IF (rflav .LE. rr4) THEN
1032 lflav1 = -2
1033 lflav2 = 1
1034 ELSE IF (rflav .LE. rr5) THEN
1035 lflav1 = -4
1036 lflav2 = 3
1037 ELSE
1038 GOTO 1
1039 END IF
1040 END IF
1041*
1042* Transformation to CMS of 2-gamma
1043 CALL lorenb (amprt(1),p1,p2,pcm)
1044 DO ip = 1,4
1045 pgg(ngg+1,ip) = pcm(ip)
1046 END DO
1047 pgg(ngg+1,5) = 0.
1048 kgg(ngg+1) = lflav1 * sign(1.,1.*kgg(igg))
1049 ngg = ngg + 1
1050*
1051 CALL lorenb (amprt(1),p1,p3,pcm)
1052 DO ip = 1,4
1053 pgg(ngg+1,ip) = pcm(ip)
1054 END DO
1055 pgg(ngg+1,5) = 0.
1056 kgg(ngg+1) = lflav2 * sign(1.,1.*kgg(igg))
1057 ngg = ngg + 1
1058*
1059 CALL lorenb (amprt(1),p1,p4,pcm)
1060 DO ip = 1,4
1061 pgg(ngg+1,ip) = pcm(ip)
1062 END DO
1063 pgg(ngg+1,5) = amprt(2)
1064 kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
1065 ngg = ngg + 1
1066 END IF
1067*
1068 RETURN
1069 END
1070*
1071*=======================================================================
1072CDECK ID>, GG2PY.
1073 SUBROUTINE gg2py (nngg)
1074*-----------------------------------------------------------------------
1075* Routine to fill PYTHIA event record by particles from GG
1076* Input: NNGG : number of particles
1077* Author: Yu.Kharlov
1078*-----------------------------------------------------------------------
1079 COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
1080 & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
1081 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1082 DOUBLE PRECISION P,V
1083 SAVE /PYJETS/
1084 DO igg = 1,nngg
1085 IF (igg .LE. 2) THEN
1086 k(igg,1) = 21
1087 ELSE
1088 k(igg,1)=1
1089 END IF
1090 k(igg,2) = kgg(igg)
1091 DO ip = 1,5
1092 p(igg,ip) = pgg(igg,ip)
1093 END DO
1094 END DO
1095 n = nngg
1096 RETURN
1097 END
1098*
1099*=======================================================================
1100CDECK ID>, PFORK.
1101 FUNCTION pfork (am0,am1,am2)
1102*-----------------------------------------------------------------------
1103* Routine to find a momentum of a secondary particle in the decay
1104* of a particle with mass AM0 in the rest to 2 particles with
1105* masses AM1 and AM2.
1106*-----------------------------------------------------------------------
1107 denom = am0**4 + am1**4 + am2**4 -
1108 & 2.*(am0*am1)**2 - 2.*(am0*am2)**2 - 2.*(am1*am2)**2
1109 IF (denom .GT. 0.) THEN
1110 pfork = sqrt (denom) / 2./am0
1111 ELSE
1112 pfork = -1.
1113 END IF
1114 RETURN
1115 END
1116*=======================================================================
1117CDECK ID>, GG2GAM.
1118 SUBROUTINE gg2gam (amas,wsq,y3,qs1,qs2,pgam1,pgam2,ptag1,ptag2,wt)
1119C***********************************************************************
1120C
1121C Generator of the Direct Two Photon Processes in S.A.Sadovsky
1122C coherent heavy ion collisions: A+B --> 1+2+3 17.01.1995
1123C with double differencial gg-luminosity function Yu.Kharlov
1124C 20.03.1995
1125C
1126C LUMFUN must be called first to initialise the generator
1127C
1128C***********************************************************************
1129C
1130C DETAILED DESCRIPTION
1131C ======== ===========
1132C
1133C INPUTS: only through common /TWOGENC/
1134C =======
1135C
1136C OUTPUTS: AMAS The mass of ion.
1137C ======== WSQ Square mass of the gamma-gamma system
1138C Y3 The rapidity of the gamma-gamma system
1139C QS1 Virtual mass squared of photon 1.
1140C QS2 Virtual mass squared of photon 2.
1141C PGAM1 Four-vector of photon 1
1142C PGAM2 Four-vector of photon 2
1143C PTAG1 Four-vector of tag 1
1144C PTAG2 Four-vector of tag 2
1145C WT The weight of the event (dummy so far)
1146C
1147C The remaining variables are internal to the program.
1148C
1149C***********************************************************************
1150 INTEGER icnter,ncntrs,nevts,lout,debug
1151 PARAMETER (ncntrs =7)
1152 COMMON /twgenc/s,eb,wmin,wmax,th1pmn,th1pmx,th2pmn,th2pmx,
1153 + st12mn,st12mx,st12rt,st22mn,st22mx,st22rt,ommin,omrat,eps,
1154 + wghtsm,wghtgt,wghtmx,wtmxfd,xmres,xgres,
1155 + icnter(ncntrs),nevts,lout,lres,lwate,debug
1156 SAVE /twgenc/
1157 REAL st12mn,st12mx,st12rt,st22mn,st22mx,st22rt
1158 REAL wmin,wmax,ommin,omrat,eps,th1pmn,th1pmx
1159 REAL th2pmn,th2pmx,s,eb,wghtsm,wghtmx,wghtgt,wtmxfd
1160 REAL xmres,xgres
1161 LOGICAL lres,lwate
1162C
82764a9c 1163 REAL alpha,pi
0b5dd071 1164 PARAMETER (alpha=1./137.036,pi=3.14159265)
1165C
1166 REAL ptag1(10),ptag2(10),pgam1(10),pgam2(10),wt
82764a9c 1167 REAL amas,wsq
0b5dd071 1168C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82764a9c 1169 REAL*4 y3,xm,q0,qs,qs1,qs2
1170 REAL*8 pa(5),pb(5)
0b5dd071 1171C
82764a9c 1172 REAL*8 pi2
0b5dd071 1173C
1174 LOGICAL lstr
1175 DATA q0,pi2,lstr / 0.060, 6.28318530718 d0, .false. /
1176 REAL*8 ggrnd
1177C
1178 999 CONTINUE
1179 icnter(1) = icnter(1) + 1
1180C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1181 IF (lstr) go to 10
1182 lstr =.true.
1183 qs = q0/sqrt(2.)
1184C
1185 pa(1)= 0D0
1186 pa(2)= 0D0
1187 pa(3)= dsqrt(1D0*(eb*eb-amas*amas))
1188 pa(4)= eb*1D0
1189 pa(5)= amas*1D0
1190C
1191 pb(1)= 0D0
1192 pb(2)= 0D0
1193 pb(3)=-pa(3)
1194 pb(4)= pa(4)
1195 pb(5)= amas*1D0
1196C
1197 10 xm = xmres
1198* Generate 2gamma rapidity y3 and mass xm
1199 CALL ranlum(y3,xm,wt)
1200C
1201 wsq = xm*xm
1202C
1203* Energy and z-momentum of 2gamma system
1204 egg = xm*cosh(y3)
1205 pgg = xm*sinh(y3)
1206*
1207* Squared masses of the photons
1208 qs1 = - qs**2 * dlog(ggrnd(0)) ! is it true that squared mass
1209 qs2 = - qs**2 * dlog(ggrnd(0)) ! is distributed in exponent?
1210 xm1 = -qs1
1211 xm2 = -qs2
1212*
1213* Energies and z-momenta of the photons
1214 xfun12 = wsq + xm1 - xm2
1215*
1216 eg1 = (egg + pgg*sqrt(1 - 4*wsq*xm1/(xfun12*xfun12))) *
1217 & xfun12/(2*wsq)
1218 eg2 = egg - eg1
1219 pg1 = sqrt (eg1**2 - xm1)
1220 pg2 = -sqrt (eg2**2 - xm2)
1221*
1222* 4-momenta of the photons
1223 DO ip=1,2
1224 pgam1(ip) = 0.
1225 pgam2(ip) = 0.
1226 END DO
1227 pgam1(3) = pg1
1228 pgam1(4) = eg1
1229 pgam2(3) = pg2
1230 pgam2(4) = eg2
1231*
1232* 4-momenta of the recoil ions
1233 DO ip=1,4
1234 ptag1(ip) = pa(ip) - pgam1(ip)
1235 ptag2(ip) = pb(ip) - pgam2(ip)
1236 END DO
1237*
1238 RETURN
1239 END
1240*
1241*=======================================================================
1242CDECK ID>, RANLUM.
1243 SUBROUTINE ranlum (ry,rm,wlum)
1244C G.V.Khaoustov
1245C 12.12.1994
1246C Corrected: Yu.Kharlov
1247C 23.05.1995
1248C 24-JUN-1997
1249C
1250C generate two random numbers (ry and rm) with probability proportional
1251C to luminosity function in ymin-ymax rapidity and xmin-xmax mass range
1252C
1253 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1254 & gvconst(4,10)
1255 REAL nc
1256 PARAMETER (isize=100000)
1257 COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1258 + xmin,xmax,amn,amx,nnms,alfun(isize)
1259 PARAMETER (NCNTRS =7)
1260 COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX,
1261 + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS,
1262 + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES,
1263 + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG
1264 LOGICAL lres,lwate
1265 REAL*4 wlum
1266 DATA ifl/0/, stpy/0./, stpm/0./
1267 REAL*8 ggrnd
1268C
1269 IF (ifl.EQ.0) THEN
1270 IF (ymin.LT.ymn .OR. ymin.GT.ymx .OR.
1271 * ymax.LT.ymn .OR. ymax.GT.ymx .OR.
1272 * xmin.LT.amn .OR. xmin.GT.amx .OR.
1273 * xmax.LT.amn .OR. xmax.GT.amx) THEN
1274 WRITE (6,*) 'Luminosity function is not defined in this range'
1275 STOP
1276 ENDIF
1277 aky = ymax-ymin
1278 akm = xmax-xmin
1279 stpy = (ymx-ymn)/nny
1280 stpm = (amx-amn)/nnms
1281 rmin = exp(cexp*(xmax))
1282 rmax = exp(cexp*(xmin))
1283 dr = rmax-rmin
1284 ifl = 1
1285 ENDIF
1286C
1287 1 CONTINUE
1288 IF (.NOT.lres) THEN
1289 rm = ggrnd(0) ! random mass
1290 rm = rmin + rm*dr
1291 rm = alog(rm)/cexp ! exponetial law
1292 ELSE
1293 wsq = xmres**2 + xmres*xgres * tan(pi*(ggrnd(0)-0.5))
1294 IF (wsq .LT. xmin*xmin) GOTO 1
1295 IF (wsq .GT. xmax*xmax) GOTO 1
1296 rm = sqrt (wsq)
1297 END IF
1298C
1299 ry = ggrnd(0) ! random Y
1300*
1301* Find y-limits for generated gg-mass: y_max = log(E_max/M_gg)
1302*
1303 ym_max = alog (egg_max/rm)
1304 ygg_max = min (ymax, ym_max)
1305 ygg_min = max (ymin,-ym_max)
1306 ry = ygg_min + ry*(ygg_max-ygg_min)
1307 iy = (ry-ymn)/stpy+1.
1308C
1309 im = (rm-amn)/stpm
1310 r = anorm*exp(cexp*rm)*ggrnd(0) ! normalization
1311 p = (ry-(iy-1)*stpy-ymn)/stpy ! four point interpolation
1312 q = (rm-im*stpm-amn)/stpm
1313 k = iy+im*(nny+1)
1314 f = (1.-p)*(1.-q)*alfun(k) + p*(1.-q)*alfun(k+1) +
1315 * q*(1.-p)*alfun(k+nny+1) + q*p*alfun(k+2+nny)
1316 IF (r .GT. f) GOTO 1
1317 wlum = 1000.*f ! differential luminosity in GeV^-1
1318*
1319 RETURN
1320 END
1321*
1322*=======================================================================
1323CDECK ID>, GGDATA.
1324 BLOCK DATA ggdata
1325*-----------------------------------------------------------------------
1326* Some useful constants and masses of MSSM particles
1327* Author: Yu.Kharlov
1328*-----------------------------------------------------------------------
1329 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1330 & gvconst(4,10)
1331 REAL nc
1332 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
1333 & amin,amax,ymin,ymax,nmas,ny, kferm,
1334 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
1335 & thetamin, costhv1, kv1,kv2,gvpar(4)
1336 CHARACTER lumfil*80
1337 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
1338 & xstote, ssbr(10)
1339 COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
1340 & xmll, xmlr, xmnl, xtanb, xmha, xmu,
1341 & xmt, xat, xmbr, xab, u11, v11
1342 COMMON/D2LParam/lz,la,Rion,Gamma
1343 INTEGER lz,la
1344 REAL Rion,Gamma
1345 PARAMETER (al = 1./128.)
1346 DATA pi /3.14159276/, alpha /al/, gev2nb /389379./,
1347 & hbarc /0.197327/
1348 DATA iproc, nevent, ilumf /1, 10000, -1/
1349 DATA lumfil /'gglum.dat'/
1350 DATA Rion /0.0/
1351 DATA amprt /110.,10.,150.,-1.,0.106/
1352 DATA thetamin /1.e-03/
1353 DATA costhv1 /1.0/
1354 DATA xsbra /1./
1355 DATA ny, nmas /50, 50/
1356 DATA xm1, xm2, xmg, xms, xmtl, xmtr,
1357 & xmll, xmlr, xmnl, xtanb, xmha, xmu,
1358 & xmt, xat, xmbr, xab
1359 & /30., 50., 250., 700., 600., 600.,
1360 & 400., 400., 400., 4., 300., -300.,
1361 & 174., 300., 300., 300./
1362 DATA gvconst / 63.0, 0.22, 6.0, 1.0, ! rho rho
1363 & 14.0, 0.22, 6.0, 1.0, ! rho omega
1364 & 7.0, 0.22, 4.0, 1.0, ! rho phi
1365 & 0.05, 0.80, 2.5, 0.0, ! rho J/psi
1366 & 0.0, 0.00, 0.0, 0.0, ! omega omega
1367 & 0.8, 0.22, 4.0, 1.0, ! omega phi
1368 & 6.E-3, 0.80, 2.5, 0.0, ! omega J/psi
1369 & 0.2, 0.22, 2.0, 1.0, ! phi phi
1370 & 3.E-3, 0.80, 1.5, 0.0, ! phi J/psi
1371 & 6.E-5, 2.00, 1.5, 0.0/ ! J/psi J/psi
1372 END
1373*
1374*=======================================================================
1375CDECK ID>, EVPRINT.
1376 SUBROUTINE evprint(Lun,Iev,P3,Ygg,ifl,Eb,ia,iz,amas,ptag1,ptag2)
1377*-----------------------------------------------------------------------
1378* Event output for internal use by GGHIC authors
1379*-----------------------------------------------------------------------
1380 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1381 DOUBLE PRECISION P,V
1382 SAVE /PYJETS/
1383 REAL ptag1(4),ptag2(4),p3(5)
1384C
1385 WRITE(lun,100) iev,n,Eb,p3(5),ygg,ifl,-ifl
1386 WRITE(lun,101) 90000+ia,iz,ptag1,amas
1387 WRITE(lun,101) 90000+ia,iz,ptag2,amas
1388 WRITE(lun,101) 90000+ 3, 0,(p3(j),j=1,5)
1389 IF (n.NE.0) THEN
1390 DO i=1,n
1391 ichg=pychge(k(i,2))/3
1392 WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5)
1393 ENDDO
1394 ENDIF
1395 RETURN
1396 100 FORMAT(i6,i3,3e16.8,2x,2i8)
1397 101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8)
1398 END
1399*
1400*=======================================================================
1401C
1402CDECK ID>, TWOINI.
1403 SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR) TWOG0965
1404C***********************************************************************TWOG0966
1405C TWOG0967
1406C Initialises the gamma-gamma generator of Ion collisions TWOGAM. TWOG0968
1407C TWOG0969
1408C INPUTS: LOUTX Logical unit for print-out TWOG0970
1409C ======= XPAR(1) EBM Beam energy GeV TWOG0971
1410C XPAR(2) WMN Minimum two-gamma mass GeV TWOG0972
1411C XPAR(3) WMX Maximum two-gamma mass GeV TWOG0973
1412C TWOG0974
1413C XPAR(8) WGHTMX Maximum weight TWOG0978
1414C XPAR(9) XMRES Mass of resonance GeV TWOG0979
1415C XPAR(10) XGRES Total width of resonance GeV TWOG0980
1416C TWOG0981
1417C IPAR(1) 0 Unweighted events TWOG0982
1418C 1 Weighted events TWOG0983
1419C IPAR(2) 0 Continuum production TWOG0984
1420C 1 Resonance formation TWOG0985
1421C TWOG0986
1422C EXTERNALS: none TWOG0987
1423C ========== TWOG0988
1424C TWOG0989
1425C AUTHOR/DATE: S.A. Sadovsky, 17 January 1995 TWOG0990
1426C =========== TWOG0991
1427C TWOG0992
1428C***********************************************************************TWOG0993
1429C IMPLICIT NONE TWOG0994
1430 INTEGER ICNTER,NCNTRS,NEVTS,LOUT,DEBUG TWOG0995
1431 PARAMETER (NCNTRS =7) TWOG0996
1432 COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX, TWOG0997
1433 + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS, TWOG0998
1434 + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES, TWOG0999
1435 + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG TWOG1000
1436 SAVE /TWGENC/ TWOG1001
1437 REAL ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT TWOG1002
1438 REAL WMIN,WMAX,OMMIN,OMRAT,EPS,TH1PMN,TH1PMX TWOG1003
1439 REAL TH2PMN,TH2PMX,S,EB,WGHTSM,WGHTMX,WGHTGT,WTMXFD TWOG1004
1440 REAL XMRES,XGRES TWOG1005
1441 LOGICAL LRES,LWATE TWOG1006
1442C TWOG1007
1443 REAL ALPHA,PI,TWOPI,XME,XMU TWOG1008
1444 PARAMETER (ALPHA=1./137.036,PI=3.14159265) TWOG1009
1445 PARAMETER (TWOPI=2.*PI,XME=0.000511,XMU=0.105658) TWOG1010
1446C TWOG1011
1447 REAL XPAR(11),YPAR(6) TWOG1012
1448 INTEGER LOUTX,IPAR(10),I TWOG1013
1449C TWOG1014
1450 EB = XPAR(1) TWOG1015
1451 WMIN = XPAR(2) TWOG1016
1452 WMAX = XPAR(3) TWOG1017
1453 WGHTMX= XPAR(8) TWOG1022
1454 WTMXFD= 0 TWOG1023
1455 XMRES = XPAR(9) TWOG1024
1456 XGRES = XPAR(10) TWOG1025
1457 LOUT = LOUTX TWOG1026
1458 LWATE = IPAR(1) .EQ. 1 TWOG1027
1459 LRES = IPAR(2) .EQ. 1 TWOG1028
1460 EPS = 1./64.*(XME*WMIN**2/EB**3)**2/100. TWOG1029
1461C+ TWOG1030
1462C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2)) TWOG1031
1463C rather than 1/2 COT(th/2) to be able to start theta from zero. TWOG1032
1464C- TWOG1033
1465 S = (2.*EB)**2 TWOG1054
1466 OMMIN = WMIN**2/(4.*EB) TWOG1055
1467 OMRAT = EB/OMMIN TWOG1056
1468C+ TWOG1057
1469C Print out initial conditions. TWOG1058
1470C- TWOG1059
1471 WRITE(LOUT,1000) TWOG1060
1472 1000 FORMAT(1H ,10X,27('-'),' TWINIT ',27('-')) TWOG1061
1473 WRITE(LOUT,1001) TWOG1062
1474 1001 FORMAT(1H ,10X,'I',60X,'I') TWOG1063
1475C
1476 WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX,
1477 + YPAR(5),YPAR(6)
1478C TWOG1065
1479 1002 FORMAT TWOG1066
1480 + (1H ,10X,'I Beam energy ',E15.5, TWOG1067
1481 + ' GeV',11X,'I',/, TWOG1068
1482 + 1H ,10X,'I Mass of ion ',F12.5, TWOG1067
1483 + ' GeV',11X,'I',/, TWOG1068
1484 + 1H ,10X,'I Gamma of ion ',F12.5, TWOG1067
1485 + ' ',11X,'I',/, TWOG1068
1486 + 1H ,10X,'I Charge of ion ',F12.5, TWOG1067
1487 + ' ',11X,'I',/, TWOG1068
1488 + 1H ,10X,'I Atomic number of ion ',F12.5, TWOG1067
1489 + ' ',11X,'I',/, TWOG1068
1490 + 1H ,10X,'I Minimum gamma-gamma mass ',F12.5, TWOG1069
1491 + ' GeV',11X,'I',/, TWOG1070
1492 + 1H ,10X,'I Maximum gamma-gamma mass ',F12.5, TWOG1071
1493 + ' GeV',11X,'I',/, TWOG1080
1494 + 1H ,10X,'I Minimum rapidity ',F12.5, TWOG1069
1495 + ' ',11X,'I',/, TWOG1070
1496 + 1H ,10X,'I Maximum rapidity ',F12.5, TWOG1071
1497 + ' ',11X,'I' ) TWOG1080
1498 IF (LRES) THEN TWOG1081
1499 WRITE(LOUT,1001) TWOG1082
1500 WRITE(LOUT,1003)XMRES,XGRES TWOG1083
1501 ENDIF TWOG1084
1502 1003 FORMAT TWOG1085
1503 + (1H ,10X,'I Resonance mass ',F12.5, TWOG1086
1504 + ' GeV',11X,'I',/, TWOG1087
1505 + 1H ,10X,'I Resonance total width ',F12.5, TWOG1088
1506 + ' GeV',11X,'I') TWOG1089
1507C
1508 IF (LWATE) THEN TWOG1090
1509 WRITE(LOUT,1001) TWOG1091
1510 WRITE(LOUT,1004) TWOG1092
1511 ELSE TWOG1095
1512 WRITE(LOUT,1001) TWOG1096
1513 WRITE(LOUT,1005)WGHTMX TWOG1097
1514 1004 FORMAT(1H ,10X,'I Produce weighted events ',27X,'I') TWOG1094
1515 1005 FORMAT(1H ,10X,'I Maximum weight ',F12.5, TWOG1099
1516 + ' ',11X,'I') TWOG1100
1517 ENDIF TWOG1101
1518C TWOG1102
1519 WRITE(LOUT,1001) TWOG1103
1520 WRITE(LOUT,1000) TWOG1104
1521C TWOG1105
1522 DO 10 I=1,NCNTRS TWOG1106
1523 10 ICNTER(I)=0 TWOG1107
1524 WGHTGT = 0. TWOG1108
1525 WGHTSM = 0. TWOG1109
1526 RETURN TWOG1110
1527 END TWOG1111
1528C=========================================================================
1529CDECK ID>, DLDMX.
1530 FUNCTION DLDMX(Rm,Ny)
1531 PARAMETER (isize=100000)
1532C S.A.Sadovsky
1533C Calculate the Luminosity integral over dY 5.02.1995
1534C in the rapidity range ymin-ymax
1535C using luminosity interpolation table
1536C Ny/2 - number of subintervals for Simpson integration
1537C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1538 COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1539 + xmin,xmax,amn,amx,nnms,alfun(isize)
1540 stpy=(ymx-ymn)/nny
1541 stpm=(amx-amn)/nnms ! Interpolation steps
1542C
1543 im=(rm-amn)/stpm
1544 q =(rm-im*stpm-amn)/stpm
1545C
1546 My = Ny/2
1547 Hy =(ymax-ymin)/My ! Integration step over Y
1548 ry = ymin
1549 iy =(ry-ymn)/stpy+1.
1550 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 1
1551 k = iy+im*(nny+1)
1552 F1 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1553 * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1554C
1555 SIM = 0.0
1556 DO 10 i=1,My
1557C
1558 ry = ymin + Hy*i
1559 iy =(ry-ymn)/stpy+1.
1560 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 3
1561 k = iy+im*(nny+1)
1562 F3 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1563 * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1564C
1565 ry = ry - Hy*0.5
1566 iy =(ry-ymn)/stpy+1.
1567 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 2
1568 k = iy+im*(nny+1)
1569 F2 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1570 * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1571C
1572 SIM = SIM + (F1 + 4.*F2 + F3)/6.
1573 10 F1 = F3
1574C
1575 DLDMX = 1000.*SIM*HY ! dL/dMx in GeV^-1
1576 RETURN
1577 END
1578*=======================================================================
1579CDECK ID>, LUMFUN.
1580 SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms)
1581C-----------------------------------------------------------------------
1582C preparation of luminosity function for generator
1583C
1584C ifl=-1 - calculate function + save in file 'file'
1585C ifl= 0 - calculate function
1586C ifl= 1 - read function from file
1587C ymin,ymax,ny - rapidity range & number of points
1588C xmn,xmx,nms - mass range & number of points
1589C Attention: dimensions of alfun are alfun(ny+1,nms+1)
1590C
1591C Called by GGINIT
1592C
1593C Author: G.V.Khaoustov 12.12.1994
1594C Modified: Yu.Kharlov 18-JUN-1997
1595C-----------------------------------------------------------------------
1596 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1597 & gvconst(4,10)
1598 REAL nc
1599 COMMON/D2LParam/lz,la,Rion,Gamma
1600 INTEGER lz,la
1601 REAL Rion,Gamma
1602 PARAMETER (isize=100000)
1603 CHARACTER*(*) file
1604 REAL work(1000)
1605 COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny,
1606 + wmn,wmx,amn,amx,nnms,alfun(isize)
1607C
1608 umn = ymin ! parameters for luminosity generation
1609 umx = ymax
1610 wmn = xmn
1611 wmx = xmx
1612C
1613C Define R : radii of the ions (in fm) (if set non-zero, use it)
1614 IF (Rion .LE. 0) THEN
1615 rion = 1.2*float(la)**(1.0/3.0)
1616 ENDIF
1617*
1618* Fing max energy of gg-system (with safity factor 2)
1619*
1620 egg_max = 2.*hbarc * gamma / rion
1621*
1622 IF(ifl.EQ.1) THEN ! read lum. function from file
1623 OPEN(unit=10,err=1,file=file,status='old',
1624 * form='formatted')
1625 READ(10, *,err=3) anorm,cexp
1626 READ(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1627 READ(10,100,err=3) (alfun(I),i=1,(nny+1)*(nnms+1))
1628 CLOSE(unit=10)
1629 PRINT *,'The Luminosity has been read from File '
1630 RETURN
1631 ENDIF
1632c
1633 IF (ny*nms .GT. isize) GOTO 2
1634 PRINT *,'Start the Luminosity calculation'
1635 nny = ny
1636 ymn = ymin
1637 ymx = ymax
1638 nnms= nms
1639 amn = xmn
1640 amx = xmx
1641c
1642 stepy=(ymx-ymn)/nny
1643 stepm=(amx-amn)/nnms
1644*
1645* Lum. function calculation
1646*
1647 DO IM = 0,nnms
1648 am=amn+im*stepm
1649* Find Y-limits for mass AM (Yu.Kharlov)
1650 ym_max = alog (egg_max/am)
1651 ym_max = max(ym_max,5.*stepy)
1652 ygg_max = min (ymax, ym_max)
1653 ygg_min = max (ymin,-ym_max)
1654 DO IY = 1,nny+1
1655 Y = ymn+(iy-1)*stepy
1656 IF (Y.LT.ygg_max .AND. Y.GT.ygg_min) THEN
1657 dl = D2LDMDY(1000.*am,Y) ! Differential Luminosity in MeV^-1
1658 ELSE
1659 dl = 0.0
1660 END IF
1661 alfun(iy+(nny+1)*im) = DL
1662 ENDDO
1663 ENDDO
1664c
1665 s=0.
1666 DO i=1,nnms+1 ! find max values for exp. parameters calcul.
1667 am=0.
1668 DO j=1,nny+1
1669 IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1))
1670 ENDDO
1671 work(i)=am
1672 s=s+am
1673 ENDDO
1674 s=(s-work(nnms+1))*stepm
1675c
1676 smax=0. ! calculation optimized parameters for exp. generator
1677 iflag=-1
1678 DO 101 i=1,nnms
1679 DO 102 j=i+1,nnms+1
1680 IF (work(j) .LE. 0.0) GOTO 102
1681 c=alog(work(i)/work(j))/((i-j)*stepm)
1682 am=work(i)/exp(c*(i-1)*stepm)
1683 se=am/c*(exp(c*(amx-amn))-1.)
1684 IF(s/se.lt.smax) GOTO 102
1685 DO k=1,nnms+1
1686 amp=am*exp(c*(k-1)*stepm)
1687 IF((work(k)-amp)/amp.GT.0.00001) THEN ! Accuracy problem
1688 GOTO 102
1689 ENDIF
1690 ENDDO
1691 smax =s/se
1692 anorm=am
1693 cexp =c
1694 iflag=1
1695 102 CONTINUE
1696 101 CONTINUE
1697 IF(iflag.LT.0) THEN
1698 PRINT *,'Fault to find best exponential approximation'
1699 STOP
1700 ENDIF
1701c
1702 IF(ifl.EQ.0) RETURN ! save lum. function
1703 OPEN(unit=10,err=1,file=file,status='unknown',
1704 * form='formatted')
1705 WRITE(10, *,err=3) anorm,cexp
1706 WRITE(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1707 WRITE(10,100,err=3) (alfun(i),i=1,(nny+1)*(nnms+1))
1708 CLOSE(unit=10)
1709 PRINT *,'The Luminosity calculation is finished'
1710 RETURN
1711 1 PRINT *,' file open error, file:',file
1712 STOP
1713 2 PRINT *,' array "ALFUN" to small to contain luminosity function'
1714 STOP
1715 3 PRINT *,' file read/writr error, file:',file
1716 STOP
1717 100 FORMAT(5e13.5)
1718 99 FORMAT(2e13.5,i5,2e13.5,i5)
1719 END
1720C======================================================================
1721CDECK ID>, D2LDMDY.
1722 REAL FUNCTION D2LDMDY(M,Y)
1723C-----------------------------------------------------------------------
1724C Calculate the Luminosity-function d^2L/dMdY
1725C where M is the invariant mass and Y the rapidity
1726C of the photon-photon-system
1727C
1728C Calculations are based on the formulae of
1729C N. Baron, G. Baur, Phys. Rev. C
1730C C. Cahn, J. D. Jackson, Nucl. Phys.
1731C
1732C Kai Hencken, 29. 9. 1994
1733C
1734C M : invariant mass (MeV)
1735C Y : rapidity
1736C
1737C D2LDMDY : diff. Lum. (MeV**-1)
1738C
1739C Called by LUMFUN
1740C-----------------------------------------------------------------------
1741 IMPLICIT NONE
1742 COMMON/D2LParam/lz,la,Rion,Gamma
1743 INTEGER lz,la
1744 REAL Rion,Gamma
1745 REAL M,Y
1746 REAL RM,W1,W2
1747 COMMON/D2LIParam/RM,W1,W2
1748
1749 REAL NClass,DeltaL
1750 EXTERNAL NClass,DeltaL
1751
1752C
1753C physical constants hbar*c and 1/alpha
1754C
1755 REAL hbarc,falphai
1756 PARAMETER (hbarc=197.327053,falphai=137.0359895)
1757
1758 REAL LClass,DL
1759C
1760C calculate photon frequencies from M,Y
1761C
1762 W1 = M/2.0 * EXP ( Y)
1763 W2 = M/2.0 * EXP (-Y)
1764C
1765C calculate Rho in MeV^-1 instead of fm
1766C
1767 RM = Rion / hbarc
1768C
1769C calculate the "classical" value first.
1770C
1771
1772 LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma)
1773
1774C
1775C substract from this the Delta-value
1776C
1777 DL = DeltaL ()
1778
1779 D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL)
1780
1781 RETURN
1782 END
1783*=======================================================================
1784CDECK ID>, NCLASS.
1785 REAL FUNCTION NClass (W, R, Gamma)
1786C
1787C the classical photon-number without impact parameter cutoff:
1788C
1789 IMPLICIT NONE
1790 REAL W,R, Gamma
1791
1792 REAL Pi
1793 PARAMETER (Pi = 3.141592)
1794
1795 REAL BesselK0,BesselK1
1796 EXTERNAL BesselK0,BesselK1
1797
1798 REAL Xi
1799
1800 Xi = W*R/Gamma
1801
1802 NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) -
1803 & Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2))
1804
1805 RETURN
1806 END
1807*=======================================================================
1808CDECK ID>, DELTAL.
1809 REAL FUNCTION DeltaL ()
1810C
1811C the difference to the classical value
1812C Called by d2LdMdY
1813C
1814 IMPLICIT NONE
1815 COMMON/D2LParam/lz,la,Rion,Gamma
1816 INTEGER lz,la
1817 REAL Rion,Gamma
1818 REAL RM,W1,W2
1819 COMMON/D2LIParam/RM,W1,W2
1820
1821 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1822 COMMON /GAUSSD/XGAUSS,WGAUSS
1823
1824 INTEGER N,I
1825 REAL b1
1826 REAL*8 t,tmin,tmax
1827 REAL*8 Sum,Int,Int2
1828
1829 REAL Itgb1
1830 EXTERNAL Itgb1
1831
1832 REAL Pi
1833 REAL*8 Accur
1834 PARAMETER (Pi = 3.141592,Accur = 1D -2)
1835
1836C
1837C integrate first over b1
1838C
1839
1840C
1841C Loop incrementing the boundary
1842C
1843 tmin = 0.00
1844 tmax = 0.25
1845 Sum = 0.00
1846
1847 100 CONTINUE
1848C
1849C Loop for the Gauss integration
1850C
1851 Int=0.0
1852 DO N=1,6
1853 Int2 = Int
1854 Int=0.0
1855 DO I=2**N-1,2**(N+1)-2
1856 t = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0
1857 b1 = RM * DEXP (t)
1858 Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2
1859 ENDDO
1860 Int = (tmax-tmin)/2.0*Int
1861 IF (Int .NE. 0.) THEN
1862 IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200
1863 END IF
1864 ENDDO
1865 IF (int .LT. 1.D-20) THEN
1866 int = 0.0D0
1867 ELSE
1868 WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE'
1869 END IF
1870
1871 200 CONTINUE
1872 Sum = Sum + Int
1873 IF (sum .NE. 0.0) THEN
1874 IF (ABS (Int2/Sum) .LT. Accur) GOTO 300
1875 ELSE
1876 IF (b1 .GT. 1.E+06) GOTO 300
1877 END IF
1878 tmin = tmax
1879 tmax = tmax + 0.5
1880 GOTO 100
1881
1882 300 CONTINUE
1883
1884 DeltaL = 4.0*Pi * Sum
1885 RETURN
1886 END
1887
1888CDECK ID>, ITGB1.
1889 REAL FUNCTION Itgb1(b)
1890C
1891C integration then over b2
1892C Called by DeltaL
1893C
1894C Corrected by S.A.Sadovsky
1895C 25.10.1994
1896 IMPLICIT NONE
1897c +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997
1898 REAL b
1899 REAL*8 b1,bmin,bmax
1900 REAL RM,W1,W2
1901 COMMON/D2LIParam/RM,W1,W2
1902 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1903 COMMON /GAUSSD/XGAUSS,WGAUSS
1904 INTEGER N,I
1905 REAL*8 b2,Int,Int2
1906 REAL*8 Itgb2
1907 EXTERNAL Itgb2
1908 REAL*8 Accur
1909 PARAMETER (Accur = 1.D-2)
1910
1911 b1 = b
1912 bmin = b1 - 2.0*RM
1913 IF (RM .GT. bmin) THEN
1914 bmin = RM
1915 ENDIF
1916 bmax = b1 + 2.0 * RM
1917 Int = 0.0
1918 DO N=1,6
1919 Int2 = Int
1920 Int = 0.0
1921 DO I=2**N-1,2**(N+1)-2
1922 b2 = (bmax-bmin)/2.0*XGAUSS(I) + (bmax+bmin)/2.0
1923 Int = Int + WGAUSS(I) * Itgb2(b1,b2) * b2
1924 END DO
1925 Int = (bmax-bmin)/2.0*Int
1926 IF (Int.NE.0.) THEN
1927 IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100
1928 ENDIF
1929 ENDDO
1930 IF (int .LT. 1.D-20) THEN
1931 int = 0.0D0
1932 ELSE
1933 WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int
1934 END IF
1935 100 CONTINUE
1936 Itgb1 = Int
1937 RETURN
1938 END
1939
1940CDECK ID>, ITGB2.
1941 REAL*8 FUNCTION Itgb2 (b1,b2)
1942C
1943C The function to be integrated over
1944C Called by Itgb1
1945C
1946C Corrected by S.A.Sadovsky
1947C 25.10.1994
1948 IMPLICIT NONE
1949 COMMON/D2LParam/lz,la,Rion,Gamma
1950 INTEGER lz,la
1951 REAL Rion,Gamma
1952 REAL*8 b1,b2,arg
1953 REAL RM,W1,W2
1954 COMMON/D2LIParam/RM,W1,W2
1955
1956 REAL Nphoton
1957 EXTERNAL Nphoton
1958
1959 arg = (b1**2 + b2**2 - 4.0 * RM**2) / (2.0*b1*b2)
1960 IF (arg .GT. 1.D0) arg = 1.D0
1961 Itgb2 = Nphoton(W1,sngl(b1),Gamma) *
1962 & Nphoton(W2,sngl(b2),Gamma) *
1963 & DACOS (arg)
1964
1965 RETURN
1966 END
1967
1968CDECK ID>, NPHOTON.
1969 REAL FUNCTION Nphoton (W,Rho,Gamma)
1970C
1971C The differential photonnumber for a nucleus
1972C (without form factor)
1973C
1974 IMPLICIT NONE
1975 REAL W,Rho,Gamma
1976
1977 REAL BESSELK1
1978 EXTERNAL BESSELK1
1979
1980 REAL Wphib,WGamma
1981
1982 REAL PI
1983 PARAMETER (PI=3.141592)
1984
1985 WGamma = W/Gamma
1986C
1987C without form factor
1988C
1989 Wphib = WGamma*BESSELK1 (WGamma*Rho)
1990
1991 Nphoton = 1.0/PI**2 * Wphib**2
1992
1993 RETURN
1994 END
1995
1996CDECK ID>, BESSELK0.
1997 REAL FUNCTION BESSELK0(X)
1998C
1999C Special functions I0,K0,I1,K1
2000C see Abramowitz&Stegun
2001C
2002 IMPLICIT NONE
2003 REAL X,Y
2004
2005 REAL BESSELI0
2006 EXTERNAL BESSELI0
2007
2008 IF (X .LT. 2.0) THEN
2009
2010 Y = X**2/4.0
2011
2012 BESSELK0 =
2013 & (-LOG(X/2.0)*BESSELI0(X))+(-.57721566+Y*(0.42278420
2014 & +Y*(0.23069756+Y*(0.3488590E-1+Y*(0.262698E-2
2015 & +Y*(0.10750E-3+Y*0.740E-5))))))
2016
2017 ELSE
2018
2019 Y = 2.0/X
2020
2021 BESSELK0 =
2022 & (EXP(-X)/SQRT(X))*(1.25331414+Y*(-0.7832358E-1
2023 & +Y*(0.2189568E-1+Y*(-0.1062446E-1+Y*(0.587872E-2
2024 & +Y*(-0.251540E-2+Y*0.53208E-3))))))
2025
2026 ENDIF
2027
2028 RETURN
2029 END
2030
2031
2032CDECK ID>, BESSELI0.
2033 REAL FUNCTION BESSELI0(X)
2034 IMPLICIT NONE
2035
2036 REAL X,Y,AX
2037
2038 AX = ABS(X)
2039
2040 IF (AX .LT. 3.75) THEN
2041
2042 Y = (X/3.75)**2
2043
2044 BESSELI0 =
2045 & 1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492
2046 & +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2)))))
2047
2048 ELSE
2049
2050 Y = 3.75/AX
2051
2052 BESSELI0 =
2053 & (EXP(AX)/SQRT(AX))*(0.39894228+Y*(0.1328592E-1
2054 & +Y*(0.225319E-2+Y*(-0.157565E-2+Y*(0.916281E-2
2055 & +Y*(-0.2057706E-1+Y*(0.2635537E-1+Y*(-0.1647633E-1
2056 & +Y*0.392377E-2))))))))
2057
2058 ENDIF
2059
2060 RETURN
2061 END
2062
2063
2064CDECK ID>, BESSELK1.
2065 REAL FUNCTION BESSELK1(X)
2066 IMPLICIT NONE
2067
2068 REAL X,Y
2069
2070 REAL BESSELI1
2071 EXTERNAL BESSELI1
2072
2073 IF (X .LT. 2.0) THEN
2074
2075 Y = X**2/4.0
2076 BESSELK1 =
2077 & (LOG(X/2.0)*BESSELI1(X))+(1.0/X)*(1.0+Y*(0.15443144
2078 & +Y*(-0.67278579+Y*(-0.18156897+Y*(-0.1919402E-1
2079 & +Y*(-0.110404E-2+Y*(-0.4686E-4)))))))
2080
2081 ELSE
2082
2083 Y=2.0/X
2084 BESSELK1 =
2085 & (EXP(-X)/SQRT(X))*(1.25331414+Y*(0.23498619
2086 & +Y*(-0.3655620E-1+Y*(0.1504268E-1+Y*(-0.780353E-2
2087 & +Y*(0.325614E-2+Y*(-0.68245E-3)))))))
2088
2089 ENDIF
2090
2091 RETURN
2092 END
2093
2094CDECK ID>, BESSELI1.
2095 REAL FUNCTION BESSELI1(X)
2096 IMPLICIT NONE
2097
2098 REAL X,Y,AX
2099
2100 AX = ABS(X)
2101
2102 IF (AX .LT. 3.75) THEN
2103
2104 Y = (X/3.75)**2
2105
2106 BESSELI1 =
2107 & AX*(0.5+Y*(0.87890594+Y*(0.51498869+Y*(0.15084934
2108 & +Y*(0.2658733E-1+Y*(0.301532E-2+Y*0.32411E-3))))))
2109
2110
2111 ELSE
2112
2113 Y = 3.75/AX
2114
2115 BESSELI1 =
2116 & 0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1
2117 & -Y*0.420059E-2))
2118
2119 BESSELI1 =
2120 & 0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2
2121 & +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1))))
2122
2123 BESSELI1 = BESSELI1 *
2124 & EXP(AX)/SQRT(AX)
2125
2126 ENDIF
2127
2128 IF (X .LT. 0) THEN
2129 BESSELI1 = -BESSELI1
2130 ENDIF
2131
2132 RETURN
2133 END
2134
2135CDECK ID>, GAUSSDAT.
2136C DATA for the Gauss integration, x-values and weight for a
2137C N=2,4,8,16,32,64 point Gauss integration.
2138C
2139C based on program GAUSS64 of D. Trautmann
2140C data from Abramowitz & Stegun
2141C
2142C KAI HENCKEN, FEBRUAR 1993
2143C
2144 BLOCKDATA GAUSSDATA
2145 IMPLICIT NONE
2146 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
2147 COMMON /GAUSSD/XGAUSS,WGAUSS
2148
2149 DATA XGAUSS(1)/ .57735026918962576D0/
2150 DATA XGAUSS(2)/-.57735026918962576D0/
2151 DATA WGAUSS(1)/ 1.00000000000000000D0/
2152 DATA WGAUSS(2)/ 1.00000000000000000D0/
2153
2154 DATA XGAUSS(3)/ .33998104358485627D0/
2155 DATA XGAUSS(4)/ .86113631159405258D0/
2156 DATA XGAUSS(5)/-.33998104358485627D0/
2157 DATA XGAUSS(6)/-.86113631159405258D0/
2158 DATA WGAUSS(3)/ .65214515486254613D0/
2159 DATA WGAUSS(4)/ .34785484513745385D0/
2160 DATA WGAUSS(5)/ .65214515486254613D0/
2161 DATA WGAUSS(6)/ .34785484513745385D0/
2162
2163 DATA XGAUSS(7)/ .18343464249564981D0/
2164 DATA XGAUSS(8)/ .52553240991632899D0/
2165 DATA XGAUSS(9)/ .79666647741362674D0/
2166 DATA XGAUSS(10)/ .96028985649753623D0/
2167 DATA XGAUSS(11)/-.18343464249564981D0/
2168 DATA XGAUSS(12)/-.52553240991632899D0/
2169 DATA XGAUSS(13)/-.79666647741362674D0/
2170 DATA XGAUSS(14)/-.96028985649753623D0/
2171 DATA WGAUSS(7)/ .36268378337836198D0/
2172 DATA WGAUSS(8)/ .31370664587788727D0/
2173 DATA WGAUSS(9)/ .22238103445337448D0/
2174 DATA WGAUSS(10)/ .10122853629037627D0/
2175 DATA WGAUSS(11)/ .36268378337836198D0/
2176 DATA WGAUSS(12)/ .31370664587788727D0/
2177 DATA WGAUSS(13)/ .22238103445337448D0/
2178 DATA WGAUSS(14)/ .10122853629037627D0/
2179
2180 DATA XGAUSS(15)/ .0950125098376374402D0/
2181 DATA XGAUSS(16)/ .281603550779258913D0/
2182 DATA XGAUSS(17)/ .458016777657227386D0/
2183 DATA XGAUSS(18)/ .617876244402643748D0/
2184 DATA XGAUSS(19)/ .755404408355003034D0/
2185 DATA XGAUSS(20)/ .865631202387831744D0/
2186 DATA XGAUSS(21)/ .944575023073232576D0/
2187 DATA XGAUSS(22)/ .989400934991649933D0/
2188 DATA XGAUSS(23)/-.0950125098376374402D0/
2189 DATA XGAUSS(24)/-.281603550779258913D0/
2190 DATA XGAUSS(25)/-.458016777657227386D0/
2191 DATA XGAUSS(26)/-.617876244402643748D0/
2192 DATA XGAUSS(27)/-.755404408355003034D0/
2193 DATA XGAUSS(28)/-.865631202387831744D0/
2194 DATA XGAUSS(29)/-.944575023073232576D0/
2195 DATA XGAUSS(30)/-.989400934991649933D0/
2196 DATA WGAUSS(15)/ .189450610455068496D0/
2197 DATA WGAUSS(16)/ .182603415044923589D0/
2198 DATA WGAUSS(17)/ .169156519395002538D0/
2199 DATA WGAUSS(18)/ .149595988816576732D0/
2200 DATA WGAUSS(19)/ .124628971255533872D0/
2201 DATA WGAUSS(20)/ .0951585116824927848D0/
2202 DATA WGAUSS(21)/ .0622535239386478929D0/
2203 DATA WGAUSS(22)/ .0271524594117540949D0/
2204 DATA WGAUSS(23)/ .189450610455068496D0/
2205 DATA WGAUSS(24)/ .182603415044923589D0/
2206 DATA WGAUSS(25)/ .169156519395002538D0/
2207 DATA WGAUSS(26)/ .149595988816576732D0/
2208 DATA WGAUSS(27)/ .124628971255533872D0/
2209 DATA WGAUSS(28)/ .0951585116824927848D0/
2210 DATA WGAUSS(29)/ .0622535239386478929D0/
2211 DATA WGAUSS(30)/ .0271524594117540949D0/
2212
2213 DATA XGAUSS(31)/ .0483076656877383162D0/
2214 DATA XGAUSS(32)/ .144471961582796493D0/
2215 DATA XGAUSS(33)/ .239287362252137075D0/
2216 DATA XGAUSS(34)/ .331868602282127650D0/
2217 DATA XGAUSS(35)/ .421351276130635345D0/
2218 DATA XGAUSS(36)/ .506899908932229390D0/
2219 DATA XGAUSS(37)/ .587715757240762329D0/
2220 DATA XGAUSS(38)/ .663044266930215201D0/
2221 DATA XGAUSS(39)/ .732182118740289680D0/
2222 DATA XGAUSS(40)/ .794483795967942407D0/
2223 DATA XGAUSS(41)/ .849367613732569970D0/
2224 DATA XGAUSS(42)/ .896321155766052124D0/
2225 DATA XGAUSS(43)/ .934906075937739689D0/
2226 DATA XGAUSS(44)/ .964762255587506430D0/
2227 DATA XGAUSS(45)/ .985611511545268335D0/
2228 DATA XGAUSS(46)/ .997263861849481564D0/
2229 DATA XGAUSS(47)/-.0483076656877383162D0/
2230 DATA XGAUSS(48)/-.144471961582796493D0/
2231 DATA XGAUSS(49)/-.239287362252137075D0/
2232 DATA XGAUSS(50)/-.331868602282127650D0/
2233 DATA XGAUSS(51)/-.421351276130635345D0/
2234 DATA XGAUSS(52)/-.506899908932229390D0/
2235 DATA XGAUSS(53)/-.587715757240762329D0/
2236 DATA XGAUSS(54)/-.663044266930215201D0/
2237 DATA XGAUSS(55)/-.732182118740289680D0/
2238 DATA XGAUSS(56)/-.794483795967942407D0/
2239 DATA XGAUSS(57)/-.849367613732569970D0/
2240 DATA XGAUSS(58)/-.896321155766052124D0/
2241 DATA XGAUSS(59)/-.934906075937739689D0/
2242 DATA XGAUSS(60)/-.964762255587506430D0/
2243 DATA XGAUSS(61)/-.985611511545268335D0/
2244 DATA XGAUSS(62)/-.997263861849481564D0/
2245 DATA WGAUSS(31)/ .0965400885147278006D0/
2246 DATA WGAUSS(32)/ .0956387200792748594D0/
2247 DATA WGAUSS(33)/ .0938443990808045654D0/
2248 DATA WGAUSS(34)/ .0911738786957638847D0/
2249 DATA WGAUSS(35)/ .0876520930044038111D0/
2250 DATA WGAUSS(36)/ .0833119242269467552D0/
2251 DATA WGAUSS(37)/ .0781938957870703065D0/
2252 DATA WGAUSS(38)/ .0723457941088485062D0/
2253 DATA WGAUSS(39)/ .0658222227763618468D0/
2254 DATA WGAUSS(40)/ .0586840934785355471D0/
2255 DATA WGAUSS(41)/ .0509980592623761762D0/
2256 DATA WGAUSS(42)/ .0428358980222266807D0/
2257 DATA WGAUSS(43)/ .0342738629130214331D0/
2258 DATA WGAUSS(44)/ .0253920653092620595D0/
2259 DATA WGAUSS(45)/ .0162743947309056706D0/
2260 DATA WGAUSS(46)/ .00701861000947009660D0/
2261 DATA WGAUSS(47)/ .0965400885147278006D0/
2262 DATA WGAUSS(48)/ .0956387200792748594D0/
2263 DATA WGAUSS(49)/ .0938443990808045654D0/
2264 DATA WGAUSS(50)/ .0911738786957638847D0/
2265 DATA WGAUSS(51)/ .0876520930044038111D0/
2266 DATA WGAUSS(52)/ .0833119242269467552D0/
2267 DATA WGAUSS(53)/ .0781938957870703065D0/
2268 DATA WGAUSS(54)/ .0723457941088485062D0/
2269 DATA WGAUSS(55)/ .0658222227763618468D0/
2270 DATA WGAUSS(56)/ .0586840934785355471D0/
2271 DATA WGAUSS(57)/ .0509980592623761762D0/
2272 DATA WGAUSS(58)/ .0428358980222266807D0/
2273 DATA WGAUSS(59)/ .0342738629130214331D0/
2274 DATA WGAUSS(60)/ .0253920653092620595D0/
2275 DATA WGAUSS(61)/ .0162743947309056706D0/
2276 DATA WGAUSS(62)/ .00701861000947009660D0/
2277
2278 DATA XGAUSS(63)/ .02435029266342443250D0/
2279 DATA XGAUSS(64)/ .0729931217877990394D0/
2280 DATA XGAUSS(65)/ .121462819296120554D0/
2281 DATA XGAUSS(66)/ .169644420423992818D0/
2282 DATA XGAUSS(67)/ .217423643740007084D0/
2283 DATA XGAUSS(68)/ .264687162208767416D0/
2284 DATA XGAUSS(69)/ .311322871990210956D0/
2285 DATA XGAUSS(70)/ .357220158337668116D0/
2286 DATA XGAUSS(71)/ .402270157963991604D0/
2287 DATA XGAUSS(72)/ .446366017253464088D0/
2288 DATA XGAUSS(73)/ .489403145707052957D0/
2289 DATA XGAUSS(74)/ .531279464019894546D0/
2290 DATA XGAUSS(75)/ .571895646202634034D0/
2291 DATA XGAUSS(76)/ .611155355172393250D0/
2292 DATA XGAUSS(77)/ .648965471254657340D0/
2293 DATA XGAUSS(78)/ .685236313054233243D0/
2294 DATA XGAUSS(79)/ .719881850171610827D0/
2295 DATA XGAUSS(80)/ .752819907260531897D0/
2296 DATA XGAUSS(81)/ .783972358943341408D0/
2297 DATA XGAUSS(82)/ .813265315122797560D0/
2298 DATA XGAUSS(83)/ .840629296252580363D0/
2299 DATA XGAUSS(84)/ .865999398154092820D0/
2300 DATA XGAUSS(85)/ .889315445995114106D0/
2301 DATA XGAUSS(86)/ .910522137078502806D0/
2302 DATA XGAUSS(87)/ .929569172131939576D0/
2303 DATA XGAUSS(88)/ .946411374858402816D0/
2304 DATA XGAUSS(89)/ .961008799652053719D0/
2305 DATA XGAUSS(90)/ .973326827789910964D0/
2306 DATA XGAUSS(91)/ .983336253884625957D0/
2307 DATA XGAUSS(92)/ .991013371476744321D0/
2308 DATA XGAUSS(93)/ .996340116771955279D0/
2309 DATA XGAUSS(94)/ .999305041735772139D0/
2310 DATA XGAUSS(95)/-.02435029266342443250D0/
2311 DATA XGAUSS(96)/-.0729931217877990394D0/
2312 DATA XGAUSS(97)/-.121462819296120554D0/
2313 DATA XGAUSS(98)/-.169644420423992818D0/
2314 DATA XGAUSS(99)/-.217423643740007084D0/
2315 DATA XGAUSS(100)/-.264687162208767416D0/
2316 DATA XGAUSS(101)/-.311322871990210956D0/
2317 DATA XGAUSS(102)/-.357220158337668116D0/
2318 DATA XGAUSS(103)/-.402270157963991604D0/
2319 DATA XGAUSS(104)/-.446366017253464088D0/
2320 DATA XGAUSS(105)/-.489403145707052957D0/
2321 DATA XGAUSS(106)/-.531279464019894546D0/
2322 DATA XGAUSS(107)/-.571895646202634034D0/
2323 DATA XGAUSS(108)/-.611155355172393250D0/
2324 DATA XGAUSS(109)/-.648965471254657340D0/
2325 DATA XGAUSS(110)/-.685236313054233243D0/
2326 DATA XGAUSS(111)/-.719881850171610827D0/
2327 DATA XGAUSS(112)/-.752819907260531897D0/
2328 DATA XGAUSS(113)/-.783972358943341408D0/
2329 DATA XGAUSS(114)/-.813265315122797560D0/
2330 DATA XGAUSS(115)/-.840629296252580363D0/
2331 DATA XGAUSS(116)/-.865999398154092820D0/
2332 DATA XGAUSS(117)/-.889315445995114106D0/
2333 DATA XGAUSS(118)/-.910522137078502806D0/
2334 DATA XGAUSS(119)/-.929569172131939576D0/
2335 DATA XGAUSS(120)/-.946411374858402816D0/
2336 DATA XGAUSS(121)/-.961008799652053719D0/
2337 DATA XGAUSS(122)/-.973326827789910964D0/
2338 DATA XGAUSS(123)/-.983336253884625957D0/
2339 DATA XGAUSS(124)/-.991013371476744321D0/
2340 DATA XGAUSS(125)/-.996340116771955279D0/
2341 DATA XGAUSS(126)/-.999305041735772139D0/
2342 DATA WGAUSS(63)/ .0486909570091397204D0/
2343 DATA WGAUSS(64)/ .0485754674415034269D0/
2344 DATA WGAUSS(65)/ .0483447622348029572D0/
2345 DATA WGAUSS(66)/ .0479993885964583077D0/
2346 DATA WGAUSS(67)/ .0475401657148303087D0/
2347 DATA WGAUSS(68)/ .0469681828162100173D0/
2348 DATA WGAUSS(69)/ .0462847965813144172D0/
2349 DATA WGAUSS(70)/ .0454916279274181445D0/
2350 DATA WGAUSS(71)/ .0445905581637565631D0/
2351 DATA WGAUSS(72)/ .0435837245293234534D0/
2352 DATA WGAUSS(73)/ .0424735151236535890D0/
2353 DATA WGAUSS(74)/ .0412625632426235286D0/
2354 DATA WGAUSS(75)/ .0399537411327203414D0/
2355 DATA WGAUSS(76)/ .0385501531786156291D0/
2356 DATA WGAUSS(77)/ .0370551285402400460D0/
2357 DATA WGAUSS(78)/ .0354722132568823838D0/
2358 DATA WGAUSS(79)/ .0338051618371416094D0/
2359 DATA WGAUSS(80)/ .0320579283548515535D0/
2360 DATA WGAUSS(81)/ .0302346570724024789D0/
2361 DATA WGAUSS(82)/ .0283396726142594832D0/
2362 DATA WGAUSS(83)/ .0263774697150546587D0/
2363 DATA WGAUSS(84)/ .0243527025687108733D0/
2364 DATA WGAUSS(85)/ .0222701738083832542D0/
2365 DATA WGAUSS(86)/ .0201348231535302094D0/
2366 DATA WGAUSS(87)/ .0179517157756973431D0/
2367 DATA WGAUSS(88)/ .0157260304760247193D0/
2368 DATA WGAUSS(89)/ .0134630478967186426D0/
2369 DATA WGAUSS(90)/ .0111681394601311288D0/
2370 DATA WGAUSS(91)/ .00884675982636394772D0/
2371 DATA WGAUSS(92)/ .00650445796897836286D0/
2372 DATA WGAUSS(93)/ .00414703326056246764D0/
2373 DATA WGAUSS(94)/ .00178328072169643295D0/
2374 DATA WGAUSS(95)/ .0486909570091397204D0/
2375 DATA WGAUSS(96)/ .0485754674415034269D0/
2376 DATA WGAUSS(97)/ .0483447622348029572D0/
2377 DATA WGAUSS(98)/ .0479993885964583077D0/
2378 DATA WGAUSS(99)/ .0475401657148303087D0/
2379 DATA WGAUSS(100)/ .0469681828162100173D0/
2380 DATA WGAUSS(101)/ .0462847965813144172D0/
2381 DATA WGAUSS(102)/ .0454916279274181445D0/
2382 DATA WGAUSS(103)/ .0445905581637565631D0/
2383 DATA WGAUSS(104)/ .0435837245293234534D0/
2384 DATA WGAUSS(105)/ .0424735151236535890D0/
2385 DATA WGAUSS(106)/ .0412625632426235286D0/
2386 DATA WGAUSS(107)/ .0399537411327203414D0/
2387 DATA WGAUSS(108)/ .0385501531786156291D0/
2388 DATA WGAUSS(109)/ .0370551285402400460D0/
2389 DATA WGAUSS(110)/ .0354722132568823838D0/
2390 DATA WGAUSS(111)/ .0338051618371416094D0/
2391 DATA WGAUSS(112)/ .0320579283548515535D0/
2392 DATA WGAUSS(113)/ .0302346570724024789D0/
2393 DATA WGAUSS(114)/ .0283396726142594832D0/
2394 DATA WGAUSS(115)/ .0263774697150546587D0/
2395 DATA WGAUSS(116)/ .0243527025687108733D0/
2396 DATA WGAUSS(117)/ .0222701738083832542D0/
2397 DATA WGAUSS(118)/ .0201348231535302094D0/
2398 DATA WGAUSS(119)/ .0179517157756973431D0/
2399 DATA WGAUSS(120)/ .0157260304760247193D0/
2400 DATA WGAUSS(121)/ .0134630478967186426D0/
2401 DATA WGAUSS(122)/ .0111681394601311288D0/
2402 DATA WGAUSS(123)/ .00884675982636394772D0/
2403 DATA WGAUSS(124)/ .00650445796897836286D0/
2404 DATA WGAUSS(125)/ .00414703326056246764D0/
2405 DATA WGAUSS(126)/ .00178328072169643295D0/
2406
2407 END
2408*=======================================================================
2409 SUBROUTINE LORENB (U,PS,PI,PF)
2410C
2411C CERN PROGLIB# U102 LORENB .VERSION KERNFOR 4.04 821124
2412C ORIG. 20/08/75 L.PAPE
2413C
2414 DOUBLE PRECISION PF4, FN
2415 DIMENSION PS(4),PI(4),PF(4)
2416
2417 IF (PS(4).EQ.U) GO TO 17
2418 PF4 = (PI(4)*PS(4)+PI(3)*PS(3)+PI(2)*PS(2)+PI(1)*PS(1)) / U
2419 FN = (PF4+PI(4)) / (PS(4)+U)
2420 PF(1)= PI(1) + FN*PS(1)
2421 PF(2)= PI(2) + FN*PS(2)
2422 PF(3)= PI(3) + FN*PS(3)
2423 PF(4)= PF4
2424 GO TO 18
2425C
2426 17 PF(1)= PI(1)
2427 PF(2)= PI(2)
2428 PF(3)= PI(3)
2429 PF(4)= PI(4)
2430C
2431 18 CONTINUE
2432C
2433 RETURN
2434C
2435 END