Transition to NewIO
[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/
100 REAL pgam1(4),pgam2(4),xpar(10),ypar(6)
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/
448 REAL pgam1(4),pgam2(4),pev(4,10)
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
1163 REAL alpha,pi,wlum,weight
1164 PARAMETER (alpha=1./137.036,pi=3.14159265)
1165C
1166 REAL ptag1(10),ptag2(10),pgam1(10),pgam2(10),wt
1167 REAL amas,wsq,q1sq,q2sq
1168C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1169 REAL*4 xmin,xmax,ymin,ymax,gmma,y3,xm,q0,qs,qs1,qs2
1170 REAL*8 pa(5),pb(5),p1(5),p2(5),p3(5)
1171C
1172 REAL*8 ph1,ph2,ak1,ak2,qk1,qk2,q1,q2,q3,e1,e2,e3
1173 REAL*8 y,y1,y2,phi,pi2,pae,aeq,beq,wps,tm3,dts,wpt
1174C
1175 LOGICAL lstr
1176 DATA q0,pi2,lstr / 0.060, 6.28318530718 d0, .false. /
1177 REAL*8 ggrnd
1178C
1179 999 CONTINUE
1180 icnter(1) = icnter(1) + 1
1181C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1182 IF (lstr) go to 10
1183 lstr =.true.
1184 qs = q0/sqrt(2.)
1185C
1186 pa(1)= 0D0
1187 pa(2)= 0D0
1188 pa(3)= dsqrt(1D0*(eb*eb-amas*amas))
1189 pa(4)= eb*1D0
1190 pa(5)= amas*1D0
1191C
1192 pb(1)= 0D0
1193 pb(2)= 0D0
1194 pb(3)=-pa(3)
1195 pb(4)= pa(4)
1196 pb(5)= amas*1D0
1197C
1198 10 xm = xmres
1199* Generate 2gamma rapidity y3 and mass xm
1200 CALL ranlum(y3,xm,wt)
1201C
1202 wsq = xm*xm
1203C
1204* Energy and z-momentum of 2gamma system
1205 egg = xm*cosh(y3)
1206 pgg = xm*sinh(y3)
1207*
1208* Squared masses of the photons
1209 qs1 = - qs**2 * dlog(ggrnd(0)) ! is it true that squared mass
1210 qs2 = - qs**2 * dlog(ggrnd(0)) ! is distributed in exponent?
1211 xm1 = -qs1
1212 xm2 = -qs2
1213*
1214* Energies and z-momenta of the photons
1215 xfun12 = wsq + xm1 - xm2
1216*
1217 eg1 = (egg + pgg*sqrt(1 - 4*wsq*xm1/(xfun12*xfun12))) *
1218 & xfun12/(2*wsq)
1219 eg2 = egg - eg1
1220 pg1 = sqrt (eg1**2 - xm1)
1221 pg2 = -sqrt (eg2**2 - xm2)
1222*
1223* 4-momenta of the photons
1224 DO ip=1,2
1225 pgam1(ip) = 0.
1226 pgam2(ip) = 0.
1227 END DO
1228 pgam1(3) = pg1
1229 pgam1(4) = eg1
1230 pgam2(3) = pg2
1231 pgam2(4) = eg2
1232*
1233* 4-momenta of the recoil ions
1234 DO ip=1,4
1235 ptag1(ip) = pa(ip) - pgam1(ip)
1236 ptag2(ip) = pb(ip) - pgam2(ip)
1237 END DO
1238*
1239 RETURN
1240 END
1241*
1242*=======================================================================
1243CDECK ID>, RANLUM.
1244 SUBROUTINE ranlum (ry,rm,wlum)
1245C G.V.Khaoustov
1246C 12.12.1994
1247C Corrected: Yu.Kharlov
1248C 23.05.1995
1249C 24-JUN-1997
1250C
1251C generate two random numbers (ry and rm) with probability proportional
1252C to luminosity function in ymin-ymax rapidity and xmin-xmax mass range
1253C
1254 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1255 & gvconst(4,10)
1256 REAL nc
1257 PARAMETER (isize=100000)
1258 COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1259 + xmin,xmax,amn,amx,nnms,alfun(isize)
1260 PARAMETER (NCNTRS =7)
1261 COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX,
1262 + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS,
1263 + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES,
1264 + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG
1265 LOGICAL lres,lwate
1266 REAL*4 wlum
1267 DATA ifl/0/, stpy/0./, stpm/0./
1268 REAL*8 ggrnd
1269C
1270 IF (ifl.EQ.0) THEN
1271 IF (ymin.LT.ymn .OR. ymin.GT.ymx .OR.
1272 * ymax.LT.ymn .OR. ymax.GT.ymx .OR.
1273 * xmin.LT.amn .OR. xmin.GT.amx .OR.
1274 * xmax.LT.amn .OR. xmax.GT.amx) THEN
1275 WRITE (6,*) 'Luminosity function is not defined in this range'
1276 STOP
1277 ENDIF
1278 aky = ymax-ymin
1279 akm = xmax-xmin
1280 stpy = (ymx-ymn)/nny
1281 stpm = (amx-amn)/nnms
1282 rmin = exp(cexp*(xmax))
1283 rmax = exp(cexp*(xmin))
1284 dr = rmax-rmin
1285 ifl = 1
1286 ENDIF
1287C
1288 1 CONTINUE
1289 IF (.NOT.lres) THEN
1290 rm = ggrnd(0) ! random mass
1291 rm = rmin + rm*dr
1292 rm = alog(rm)/cexp ! exponetial law
1293 ELSE
1294 wsq = xmres**2 + xmres*xgres * tan(pi*(ggrnd(0)-0.5))
1295 IF (wsq .LT. xmin*xmin) GOTO 1
1296 IF (wsq .GT. xmax*xmax) GOTO 1
1297 rm = sqrt (wsq)
1298 END IF
1299C
1300 ry = ggrnd(0) ! random Y
1301*
1302* Find y-limits for generated gg-mass: y_max = log(E_max/M_gg)
1303*
1304 ym_max = alog (egg_max/rm)
1305 ygg_max = min (ymax, ym_max)
1306 ygg_min = max (ymin,-ym_max)
1307 ry = ygg_min + ry*(ygg_max-ygg_min)
1308 iy = (ry-ymn)/stpy+1.
1309C
1310 im = (rm-amn)/stpm
1311 r = anorm*exp(cexp*rm)*ggrnd(0) ! normalization
1312 p = (ry-(iy-1)*stpy-ymn)/stpy ! four point interpolation
1313 q = (rm-im*stpm-amn)/stpm
1314 k = iy+im*(nny+1)
1315 f = (1.-p)*(1.-q)*alfun(k) + p*(1.-q)*alfun(k+1) +
1316 * q*(1.-p)*alfun(k+nny+1) + q*p*alfun(k+2+nny)
1317 IF (r .GT. f) GOTO 1
1318 wlum = 1000.*f ! differential luminosity in GeV^-1
1319*
1320 RETURN
1321 END
1322*
1323*=======================================================================
1324CDECK ID>, GGDATA.
1325 BLOCK DATA ggdata
1326*-----------------------------------------------------------------------
1327* Some useful constants and masses of MSSM particles
1328* Author: Yu.Kharlov
1329*-----------------------------------------------------------------------
1330 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1331 & gvconst(4,10)
1332 REAL nc
1333 COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
1334 & amin,amax,ymin,ymax,nmas,ny, kferm,
1335 & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
1336 & thetamin, costhv1, kv1,kv2,gvpar(4)
1337 CHARACTER lumfil*80
1338 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
1339 & xstote, ssbr(10)
1340 COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
1341 & xmll, xmlr, xmnl, xtanb, xmha, xmu,
1342 & xmt, xat, xmbr, xab, u11, v11
1343 COMMON/D2LParam/lz,la,Rion,Gamma
1344 INTEGER lz,la
1345 REAL Rion,Gamma
1346 PARAMETER (al = 1./128.)
1347 DATA pi /3.14159276/, alpha /al/, gev2nb /389379./,
1348 & hbarc /0.197327/
1349 DATA iproc, nevent, ilumf /1, 10000, -1/
1350 DATA lumfil /'gglum.dat'/
1351 DATA Rion /0.0/
1352 DATA amprt /110.,10.,150.,-1.,0.106/
1353 DATA thetamin /1.e-03/
1354 DATA costhv1 /1.0/
1355 DATA xsbra /1./
1356 DATA ny, nmas /50, 50/
1357 DATA xm1, xm2, xmg, xms, xmtl, xmtr,
1358 & xmll, xmlr, xmnl, xtanb, xmha, xmu,
1359 & xmt, xat, xmbr, xab
1360 & /30., 50., 250., 700., 600., 600.,
1361 & 400., 400., 400., 4., 300., -300.,
1362 & 174., 300., 300., 300./
1363 DATA gvconst / 63.0, 0.22, 6.0, 1.0, ! rho rho
1364 & 14.0, 0.22, 6.0, 1.0, ! rho omega
1365 & 7.0, 0.22, 4.0, 1.0, ! rho phi
1366 & 0.05, 0.80, 2.5, 0.0, ! rho J/psi
1367 & 0.0, 0.00, 0.0, 0.0, ! omega omega
1368 & 0.8, 0.22, 4.0, 1.0, ! omega phi
1369 & 6.E-3, 0.80, 2.5, 0.0, ! omega J/psi
1370 & 0.2, 0.22, 2.0, 1.0, ! phi phi
1371 & 3.E-3, 0.80, 1.5, 0.0, ! phi J/psi
1372 & 6.E-5, 2.00, 1.5, 0.0/ ! J/psi J/psi
1373 END
1374*
1375*=======================================================================
1376CDECK ID>, EVPRINT.
1377 SUBROUTINE evprint(Lun,Iev,P3,Ygg,ifl,Eb,ia,iz,amas,ptag1,ptag2)
1378*-----------------------------------------------------------------------
1379* Event output for internal use by GGHIC authors
1380*-----------------------------------------------------------------------
1381 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1382 DOUBLE PRECISION P,V
1383 SAVE /PYJETS/
1384 REAL ptag1(4),ptag2(4),p3(5)
1385C
1386 WRITE(lun,100) iev,n,Eb,p3(5),ygg,ifl,-ifl
1387 WRITE(lun,101) 90000+ia,iz,ptag1,amas
1388 WRITE(lun,101) 90000+ia,iz,ptag2,amas
1389 WRITE(lun,101) 90000+ 3, 0,(p3(j),j=1,5)
1390 IF (n.NE.0) THEN
1391 DO i=1,n
1392 ichg=pychge(k(i,2))/3
1393 WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5)
1394 ENDDO
1395 ENDIF
1396 RETURN
1397 100 FORMAT(i6,i3,3e16.8,2x,2i8)
1398 101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8)
1399 END
1400*
1401*=======================================================================
1402C
1403CDECK ID>, TWOINI.
1404 SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR) TWOG0965
1405C***********************************************************************TWOG0966
1406C TWOG0967
1407C Initialises the gamma-gamma generator of Ion collisions TWOGAM. TWOG0968
1408C TWOG0969
1409C INPUTS: LOUTX Logical unit for print-out TWOG0970
1410C ======= XPAR(1) EBM Beam energy GeV TWOG0971
1411C XPAR(2) WMN Minimum two-gamma mass GeV TWOG0972
1412C XPAR(3) WMX Maximum two-gamma mass GeV TWOG0973
1413C TWOG0974
1414C XPAR(8) WGHTMX Maximum weight TWOG0978
1415C XPAR(9) XMRES Mass of resonance GeV TWOG0979
1416C XPAR(10) XGRES Total width of resonance GeV TWOG0980
1417C TWOG0981
1418C IPAR(1) 0 Unweighted events TWOG0982
1419C 1 Weighted events TWOG0983
1420C IPAR(2) 0 Continuum production TWOG0984
1421C 1 Resonance formation TWOG0985
1422C TWOG0986
1423C EXTERNALS: none TWOG0987
1424C ========== TWOG0988
1425C TWOG0989
1426C AUTHOR/DATE: S.A. Sadovsky, 17 January 1995 TWOG0990
1427C =========== TWOG0991
1428C TWOG0992
1429C***********************************************************************TWOG0993
1430C IMPLICIT NONE TWOG0994
1431 INTEGER ICNTER,NCNTRS,NEVTS,LOUT,DEBUG TWOG0995
1432 PARAMETER (NCNTRS =7) TWOG0996
1433 COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX, TWOG0997
1434 + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS, TWOG0998
1435 + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES, TWOG0999
1436 + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG TWOG1000
1437 SAVE /TWGENC/ TWOG1001
1438 REAL ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT TWOG1002
1439 REAL WMIN,WMAX,OMMIN,OMRAT,EPS,TH1PMN,TH1PMX TWOG1003
1440 REAL TH2PMN,TH2PMX,S,EB,WGHTSM,WGHTMX,WGHTGT,WTMXFD TWOG1004
1441 REAL XMRES,XGRES TWOG1005
1442 LOGICAL LRES,LWATE TWOG1006
1443C TWOG1007
1444 REAL ALPHA,PI,TWOPI,XME,XMU TWOG1008
1445 PARAMETER (ALPHA=1./137.036,PI=3.14159265) TWOG1009
1446 PARAMETER (TWOPI=2.*PI,XME=0.000511,XMU=0.105658) TWOG1010
1447C TWOG1011
1448 REAL XPAR(11),YPAR(6) TWOG1012
1449 INTEGER LOUTX,IPAR(10),I TWOG1013
1450C TWOG1014
1451 EB = XPAR(1) TWOG1015
1452 WMIN = XPAR(2) TWOG1016
1453 WMAX = XPAR(3) TWOG1017
1454 WGHTMX= XPAR(8) TWOG1022
1455 WTMXFD= 0 TWOG1023
1456 XMRES = XPAR(9) TWOG1024
1457 XGRES = XPAR(10) TWOG1025
1458 LOUT = LOUTX TWOG1026
1459 LWATE = IPAR(1) .EQ. 1 TWOG1027
1460 LRES = IPAR(2) .EQ. 1 TWOG1028
1461 EPS = 1./64.*(XME*WMIN**2/EB**3)**2/100. TWOG1029
1462C+ TWOG1030
1463C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2)) TWOG1031
1464C rather than 1/2 COT(th/2) to be able to start theta from zero. TWOG1032
1465C- TWOG1033
1466 S = (2.*EB)**2 TWOG1054
1467 OMMIN = WMIN**2/(4.*EB) TWOG1055
1468 OMRAT = EB/OMMIN TWOG1056
1469C+ TWOG1057
1470C Print out initial conditions. TWOG1058
1471C- TWOG1059
1472 WRITE(LOUT,1000) TWOG1060
1473 1000 FORMAT(1H ,10X,27('-'),' TWINIT ',27('-')) TWOG1061
1474 WRITE(LOUT,1001) TWOG1062
1475 1001 FORMAT(1H ,10X,'I',60X,'I') TWOG1063
1476C
1477 WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX,
1478 + YPAR(5),YPAR(6)
1479C TWOG1065
1480 1002 FORMAT TWOG1066
1481 + (1H ,10X,'I Beam energy ',E15.5, TWOG1067
1482 + ' GeV',11X,'I',/, TWOG1068
1483 + 1H ,10X,'I Mass of ion ',F12.5, TWOG1067
1484 + ' GeV',11X,'I',/, TWOG1068
1485 + 1H ,10X,'I Gamma of ion ',F12.5, TWOG1067
1486 + ' ',11X,'I',/, TWOG1068
1487 + 1H ,10X,'I Charge of ion ',F12.5, TWOG1067
1488 + ' ',11X,'I',/, TWOG1068
1489 + 1H ,10X,'I Atomic number of ion ',F12.5, TWOG1067
1490 + ' ',11X,'I',/, TWOG1068
1491 + 1H ,10X,'I Minimum gamma-gamma mass ',F12.5, TWOG1069
1492 + ' GeV',11X,'I',/, TWOG1070
1493 + 1H ,10X,'I Maximum gamma-gamma mass ',F12.5, TWOG1071
1494 + ' GeV',11X,'I',/, TWOG1080
1495 + 1H ,10X,'I Minimum rapidity ',F12.5, TWOG1069
1496 + ' ',11X,'I',/, TWOG1070
1497 + 1H ,10X,'I Maximum rapidity ',F12.5, TWOG1071
1498 + ' ',11X,'I' ) TWOG1080
1499 IF (LRES) THEN TWOG1081
1500 WRITE(LOUT,1001) TWOG1082
1501 WRITE(LOUT,1003)XMRES,XGRES TWOG1083
1502 ENDIF TWOG1084
1503 1003 FORMAT TWOG1085
1504 + (1H ,10X,'I Resonance mass ',F12.5, TWOG1086
1505 + ' GeV',11X,'I',/, TWOG1087
1506 + 1H ,10X,'I Resonance total width ',F12.5, TWOG1088
1507 + ' GeV',11X,'I') TWOG1089
1508C
1509 IF (LWATE) THEN TWOG1090
1510 WRITE(LOUT,1001) TWOG1091
1511 WRITE(LOUT,1004) TWOG1092
1512 ELSE TWOG1095
1513 WRITE(LOUT,1001) TWOG1096
1514 WRITE(LOUT,1005)WGHTMX TWOG1097
1515 1004 FORMAT(1H ,10X,'I Produce weighted events ',27X,'I') TWOG1094
1516 1005 FORMAT(1H ,10X,'I Maximum weight ',F12.5, TWOG1099
1517 + ' ',11X,'I') TWOG1100
1518 ENDIF TWOG1101
1519C TWOG1102
1520 WRITE(LOUT,1001) TWOG1103
1521 WRITE(LOUT,1000) TWOG1104
1522C TWOG1105
1523 DO 10 I=1,NCNTRS TWOG1106
1524 10 ICNTER(I)=0 TWOG1107
1525 WGHTGT = 0. TWOG1108
1526 WGHTSM = 0. TWOG1109
1527 RETURN TWOG1110
1528 END TWOG1111
1529C=========================================================================
1530CDECK ID>, DLDMX.
1531 FUNCTION DLDMX(Rm,Ny)
1532 PARAMETER (isize=100000)
1533C S.A.Sadovsky
1534C Calculate the Luminosity integral over dY 5.02.1995
1535C in the rapidity range ymin-ymax
1536C using luminosity interpolation table
1537C Ny/2 - number of subintervals for Simpson integration
1538C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1539 COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1540 + xmin,xmax,amn,amx,nnms,alfun(isize)
1541 stpy=(ymx-ymn)/nny
1542 stpm=(amx-amn)/nnms ! Interpolation steps
1543C
1544 im=(rm-amn)/stpm
1545 q =(rm-im*stpm-amn)/stpm
1546C
1547 My = Ny/2
1548 Hy =(ymax-ymin)/My ! Integration step over Y
1549 ry = ymin
1550 iy =(ry-ymn)/stpy+1.
1551 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 1
1552 k = iy+im*(nny+1)
1553 F1 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1554 * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1555C
1556 SIM = 0.0
1557 DO 10 i=1,My
1558C
1559 ry = ymin + Hy*i
1560 iy =(ry-ymn)/stpy+1.
1561 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 3
1562 k = iy+im*(nny+1)
1563 F3 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1564 * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1565C
1566 ry = ry - Hy*0.5
1567 iy =(ry-ymn)/stpy+1.
1568 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 2
1569 k = iy+im*(nny+1)
1570 F2 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1571 * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1572C
1573 SIM = SIM + (F1 + 4.*F2 + F3)/6.
1574 10 F1 = F3
1575C
1576 DLDMX = 1000.*SIM*HY ! dL/dMx in GeV^-1
1577 RETURN
1578 END
1579*=======================================================================
1580CDECK ID>, LUMFUN.
1581 SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms)
1582C-----------------------------------------------------------------------
1583C preparation of luminosity function for generator
1584C
1585C ifl=-1 - calculate function + save in file 'file'
1586C ifl= 0 - calculate function
1587C ifl= 1 - read function from file
1588C ymin,ymax,ny - rapidity range & number of points
1589C xmn,xmx,nms - mass range & number of points
1590C Attention: dimensions of alfun are alfun(ny+1,nms+1)
1591C
1592C Called by GGINIT
1593C
1594C Author: G.V.Khaoustov 12.12.1994
1595C Modified: Yu.Kharlov 18-JUN-1997
1596C-----------------------------------------------------------------------
1597 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1598 & gvconst(4,10)
1599 REAL nc
1600 COMMON/D2LParam/lz,la,Rion,Gamma
1601 INTEGER lz,la
1602 REAL Rion,Gamma
1603 PARAMETER (isize=100000)
1604 CHARACTER*(*) file
1605 REAL work(1000)
1606 COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny,
1607 + wmn,wmx,amn,amx,nnms,alfun(isize)
1608C
1609 umn = ymin ! parameters for luminosity generation
1610 umx = ymax
1611 wmn = xmn
1612 wmx = xmx
1613C
1614C Define R : radii of the ions (in fm) (if set non-zero, use it)
1615 IF (Rion .LE. 0) THEN
1616 rion = 1.2*float(la)**(1.0/3.0)
1617 ENDIF
1618*
1619* Fing max energy of gg-system (with safity factor 2)
1620*
1621 egg_max = 2.*hbarc * gamma / rion
1622*
1623 IF(ifl.EQ.1) THEN ! read lum. function from file
1624 OPEN(unit=10,err=1,file=file,status='old',
1625 * form='formatted')
1626 READ(10, *,err=3) anorm,cexp
1627 READ(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1628 READ(10,100,err=3) (alfun(I),i=1,(nny+1)*(nnms+1))
1629 CLOSE(unit=10)
1630 PRINT *,'The Luminosity has been read from File '
1631 RETURN
1632 ENDIF
1633c
1634 IF (ny*nms .GT. isize) GOTO 2
1635 PRINT *,'Start the Luminosity calculation'
1636 nny = ny
1637 ymn = ymin
1638 ymx = ymax
1639 nnms= nms
1640 amn = xmn
1641 amx = xmx
1642c
1643 stepy=(ymx-ymn)/nny
1644 stepm=(amx-amn)/nnms
1645*
1646* Lum. function calculation
1647*
1648 DO IM = 0,nnms
1649 am=amn+im*stepm
1650* Find Y-limits for mass AM (Yu.Kharlov)
1651 ym_max = alog (egg_max/am)
1652 ym_max = max(ym_max,5.*stepy)
1653 ygg_max = min (ymax, ym_max)
1654 ygg_min = max (ymin,-ym_max)
1655 DO IY = 1,nny+1
1656 Y = ymn+(iy-1)*stepy
1657 IF (Y.LT.ygg_max .AND. Y.GT.ygg_min) THEN
1658 dl = D2LDMDY(1000.*am,Y) ! Differential Luminosity in MeV^-1
1659 ELSE
1660 dl = 0.0
1661 END IF
1662 alfun(iy+(nny+1)*im) = DL
1663 ENDDO
1664 ENDDO
1665c
1666 s=0.
1667 DO i=1,nnms+1 ! find max values for exp. parameters calcul.
1668 am=0.
1669 DO j=1,nny+1
1670 IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1))
1671 ENDDO
1672 work(i)=am
1673 s=s+am
1674 ENDDO
1675 s=(s-work(nnms+1))*stepm
1676c
1677 smax=0. ! calculation optimized parameters for exp. generator
1678 iflag=-1
1679 DO 101 i=1,nnms
1680 DO 102 j=i+1,nnms+1
1681 IF (work(j) .LE. 0.0) GOTO 102
1682 c=alog(work(i)/work(j))/((i-j)*stepm)
1683 am=work(i)/exp(c*(i-1)*stepm)
1684 se=am/c*(exp(c*(amx-amn))-1.)
1685 IF(s/se.lt.smax) GOTO 102
1686 DO k=1,nnms+1
1687 amp=am*exp(c*(k-1)*stepm)
1688 IF((work(k)-amp)/amp.GT.0.00001) THEN ! Accuracy problem
1689 GOTO 102
1690 ENDIF
1691 ENDDO
1692 smax =s/se
1693 anorm=am
1694 cexp =c
1695 iflag=1
1696 102 CONTINUE
1697 101 CONTINUE
1698 IF(iflag.LT.0) THEN
1699 PRINT *,'Fault to find best exponential approximation'
1700 STOP
1701 ENDIF
1702c
1703 IF(ifl.EQ.0) RETURN ! save lum. function
1704 OPEN(unit=10,err=1,file=file,status='unknown',
1705 * form='formatted')
1706 WRITE(10, *,err=3) anorm,cexp
1707 WRITE(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1708 WRITE(10,100,err=3) (alfun(i),i=1,(nny+1)*(nnms+1))
1709 CLOSE(unit=10)
1710 PRINT *,'The Luminosity calculation is finished'
1711 RETURN
1712 1 PRINT *,' file open error, file:',file
1713 STOP
1714 2 PRINT *,' array "ALFUN" to small to contain luminosity function'
1715 STOP
1716 3 PRINT *,' file read/writr error, file:',file
1717 STOP
1718 100 FORMAT(5e13.5)
1719 99 FORMAT(2e13.5,i5,2e13.5,i5)
1720 END
1721C======================================================================
1722CDECK ID>, D2LDMDY.
1723 REAL FUNCTION D2LDMDY(M,Y)
1724C-----------------------------------------------------------------------
1725C Calculate the Luminosity-function d^2L/dMdY
1726C where M is the invariant mass and Y the rapidity
1727C of the photon-photon-system
1728C
1729C Calculations are based on the formulae of
1730C N. Baron, G. Baur, Phys. Rev. C
1731C C. Cahn, J. D. Jackson, Nucl. Phys.
1732C
1733C Kai Hencken, 29. 9. 1994
1734C
1735C M : invariant mass (MeV)
1736C Y : rapidity
1737C
1738C D2LDMDY : diff. Lum. (MeV**-1)
1739C
1740C Called by LUMFUN
1741C-----------------------------------------------------------------------
1742 IMPLICIT NONE
1743 COMMON/D2LParam/lz,la,Rion,Gamma
1744 INTEGER lz,la
1745 REAL Rion,Gamma
1746 REAL M,Y
1747 REAL RM,W1,W2
1748 COMMON/D2LIParam/RM,W1,W2
1749
1750 REAL NClass,DeltaL
1751 EXTERNAL NClass,DeltaL
1752
1753C
1754C physical constants hbar*c and 1/alpha
1755C
1756 REAL hbarc,falphai
1757 PARAMETER (hbarc=197.327053,falphai=137.0359895)
1758
1759 REAL LClass,DL
1760C
1761C calculate photon frequencies from M,Y
1762C
1763 W1 = M/2.0 * EXP ( Y)
1764 W2 = M/2.0 * EXP (-Y)
1765C
1766C calculate Rho in MeV^-1 instead of fm
1767C
1768 RM = Rion / hbarc
1769C
1770C calculate the "classical" value first.
1771C
1772
1773 LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma)
1774
1775C
1776C substract from this the Delta-value
1777C
1778 DL = DeltaL ()
1779
1780 D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL)
1781
1782 RETURN
1783 END
1784*=======================================================================
1785CDECK ID>, NCLASS.
1786 REAL FUNCTION NClass (W, R, Gamma)
1787C
1788C the classical photon-number without impact parameter cutoff:
1789C
1790 IMPLICIT NONE
1791 REAL W,R, Gamma
1792
1793 REAL Pi
1794 PARAMETER (Pi = 3.141592)
1795
1796 REAL BesselK0,BesselK1
1797 EXTERNAL BesselK0,BesselK1
1798
1799 REAL Xi
1800
1801 Xi = W*R/Gamma
1802
1803 NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) -
1804 & Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2))
1805
1806 RETURN
1807 END
1808*=======================================================================
1809CDECK ID>, DELTAL.
1810 REAL FUNCTION DeltaL ()
1811C
1812C the difference to the classical value
1813C Called by d2LdMdY
1814C
1815 IMPLICIT NONE
1816 COMMON/D2LParam/lz,la,Rion,Gamma
1817 INTEGER lz,la
1818 REAL Rion,Gamma
1819 REAL RM,W1,W2
1820 COMMON/D2LIParam/RM,W1,W2
1821
1822 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1823 COMMON /GAUSSD/XGAUSS,WGAUSS
1824
1825 INTEGER N,I
1826 REAL b1
1827 REAL*8 t,tmin,tmax
1828 REAL*8 Sum,Int,Int2
1829
1830 REAL Itgb1
1831 EXTERNAL Itgb1
1832
1833 REAL Pi
1834 REAL*8 Accur
1835 PARAMETER (Pi = 3.141592,Accur = 1D -2)
1836
1837C
1838C integrate first over b1
1839C
1840
1841C
1842C Loop incrementing the boundary
1843C
1844 tmin = 0.00
1845 tmax = 0.25
1846 Sum = 0.00
1847
1848 100 CONTINUE
1849C
1850C Loop for the Gauss integration
1851C
1852 Int=0.0
1853 DO N=1,6
1854 Int2 = Int
1855 Int=0.0
1856 DO I=2**N-1,2**(N+1)-2
1857 t = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0
1858 b1 = RM * DEXP (t)
1859 Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2
1860 ENDDO
1861 Int = (tmax-tmin)/2.0*Int
1862 IF (Int .NE. 0.) THEN
1863 IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200
1864 END IF
1865 ENDDO
1866 IF (int .LT. 1.D-20) THEN
1867 int = 0.0D0
1868 ELSE
1869 WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE'
1870 END IF
1871
1872 200 CONTINUE
1873 Sum = Sum + Int
1874 IF (sum .NE. 0.0) THEN
1875 IF (ABS (Int2/Sum) .LT. Accur) GOTO 300
1876 ELSE
1877 IF (b1 .GT. 1.E+06) GOTO 300
1878 END IF
1879 tmin = tmax
1880 tmax = tmax + 0.5
1881 GOTO 100
1882
1883 300 CONTINUE
1884
1885 DeltaL = 4.0*Pi * Sum
1886 RETURN
1887 END
1888
1889CDECK ID>, ITGB1.
1890 REAL FUNCTION Itgb1(b)
1891C
1892C integration then over b2
1893C Called by DeltaL
1894C
1895C Corrected by S.A.Sadovsky
1896C 25.10.1994
1897 IMPLICIT NONE
1898c +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997
1899 REAL b
1900 REAL*8 b1,bmin,bmax
1901 REAL RM,W1,W2
1902 COMMON/D2LIParam/RM,W1,W2
1903 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1904 COMMON /GAUSSD/XGAUSS,WGAUSS
1905 INTEGER N,I
1906 REAL*8 b2,Int,Int2
1907 REAL*8 Itgb2
1908 EXTERNAL Itgb2
1909 REAL*8 Accur
1910 PARAMETER (Accur = 1.D-2)
1911
1912 b1 = b
1913 bmin = b1 - 2.0*RM
1914 IF (RM .GT. bmin) THEN
1915 bmin = RM
1916 ENDIF
1917 bmax = b1 + 2.0 * RM
1918 Int = 0.0
1919 DO N=1,6
1920 Int2 = Int
1921 Int = 0.0
1922 DO I=2**N-1,2**(N+1)-2
1923 b2 = (bmax-bmin)/2.0*XGAUSS(I) + (bmax+bmin)/2.0
1924 Int = Int + WGAUSS(I) * Itgb2(b1,b2) * b2
1925 END DO
1926 Int = (bmax-bmin)/2.0*Int
1927 IF (Int.NE.0.) THEN
1928 IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100
1929 ENDIF
1930 ENDDO
1931 IF (int .LT. 1.D-20) THEN
1932 int = 0.0D0
1933 ELSE
1934 WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int
1935 END IF
1936 100 CONTINUE
1937 Itgb1 = Int
1938 RETURN
1939 END
1940
1941CDECK ID>, ITGB2.
1942 REAL*8 FUNCTION Itgb2 (b1,b2)
1943C
1944C The function to be integrated over
1945C Called by Itgb1
1946C
1947C Corrected by S.A.Sadovsky
1948C 25.10.1994
1949 IMPLICIT NONE
1950 COMMON/D2LParam/lz,la,Rion,Gamma
1951 INTEGER lz,la
1952 REAL Rion,Gamma
1953 REAL*8 b1,b2,arg
1954 REAL RM,W1,W2
1955 COMMON/D2LIParam/RM,W1,W2
1956
1957 REAL Nphoton
1958 EXTERNAL Nphoton
1959
1960 arg = (b1**2 + b2**2 - 4.0 * RM**2) / (2.0*b1*b2)
1961 IF (arg .GT. 1.D0) arg = 1.D0
1962 Itgb2 = Nphoton(W1,sngl(b1),Gamma) *
1963 & Nphoton(W2,sngl(b2),Gamma) *
1964 & DACOS (arg)
1965
1966 RETURN
1967 END
1968
1969CDECK ID>, NPHOTON.
1970 REAL FUNCTION Nphoton (W,Rho,Gamma)
1971C
1972C The differential photonnumber for a nucleus
1973C (without form factor)
1974C
1975 IMPLICIT NONE
1976 REAL W,Rho,Gamma
1977
1978 REAL BESSELK1
1979 EXTERNAL BESSELK1
1980
1981 REAL Wphib,WGamma
1982
1983 REAL PI
1984 PARAMETER (PI=3.141592)
1985
1986 WGamma = W/Gamma
1987C
1988C without form factor
1989C
1990 Wphib = WGamma*BESSELK1 (WGamma*Rho)
1991
1992 Nphoton = 1.0/PI**2 * Wphib**2
1993
1994 RETURN
1995 END
1996
1997CDECK ID>, BESSELK0.
1998 REAL FUNCTION BESSELK0(X)
1999C
2000C Special functions I0,K0,I1,K1
2001C see Abramowitz&Stegun
2002C
2003 IMPLICIT NONE
2004 REAL X,Y
2005
2006 REAL BESSELI0
2007 EXTERNAL BESSELI0
2008
2009 IF (X .LT. 2.0) THEN
2010
2011 Y = X**2/4.0
2012
2013 BESSELK0 =
2014 & (-LOG(X/2.0)*BESSELI0(X))+(-.57721566+Y*(0.42278420
2015 & +Y*(0.23069756+Y*(0.3488590E-1+Y*(0.262698E-2
2016 & +Y*(0.10750E-3+Y*0.740E-5))))))
2017
2018 ELSE
2019
2020 Y = 2.0/X
2021
2022 BESSELK0 =
2023 & (EXP(-X)/SQRT(X))*(1.25331414+Y*(-0.7832358E-1
2024 & +Y*(0.2189568E-1+Y*(-0.1062446E-1+Y*(0.587872E-2
2025 & +Y*(-0.251540E-2+Y*0.53208E-3))))))
2026
2027 ENDIF
2028
2029 RETURN
2030 END
2031
2032
2033CDECK ID>, BESSELI0.
2034 REAL FUNCTION BESSELI0(X)
2035 IMPLICIT NONE
2036
2037 REAL X,Y,AX
2038
2039 AX = ABS(X)
2040
2041 IF (AX .LT. 3.75) THEN
2042
2043 Y = (X/3.75)**2
2044
2045 BESSELI0 =
2046 & 1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492
2047 & +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2)))))
2048
2049 ELSE
2050
2051 Y = 3.75/AX
2052
2053 BESSELI0 =
2054 & (EXP(AX)/SQRT(AX))*(0.39894228+Y*(0.1328592E-1
2055 & +Y*(0.225319E-2+Y*(-0.157565E-2+Y*(0.916281E-2
2056 & +Y*(-0.2057706E-1+Y*(0.2635537E-1+Y*(-0.1647633E-1
2057 & +Y*0.392377E-2))))))))
2058
2059 ENDIF
2060
2061 RETURN
2062 END
2063
2064
2065CDECK ID>, BESSELK1.
2066 REAL FUNCTION BESSELK1(X)
2067 IMPLICIT NONE
2068
2069 REAL X,Y
2070
2071 REAL BESSELI1
2072 EXTERNAL BESSELI1
2073
2074 IF (X .LT. 2.0) THEN
2075
2076 Y = X**2/4.0
2077 BESSELK1 =
2078 & (LOG(X/2.0)*BESSELI1(X))+(1.0/X)*(1.0+Y*(0.15443144
2079 & +Y*(-0.67278579+Y*(-0.18156897+Y*(-0.1919402E-1
2080 & +Y*(-0.110404E-2+Y*(-0.4686E-4)))))))
2081
2082 ELSE
2083
2084 Y=2.0/X
2085 BESSELK1 =
2086 & (EXP(-X)/SQRT(X))*(1.25331414+Y*(0.23498619
2087 & +Y*(-0.3655620E-1+Y*(0.1504268E-1+Y*(-0.780353E-2
2088 & +Y*(0.325614E-2+Y*(-0.68245E-3)))))))
2089
2090 ENDIF
2091
2092 RETURN
2093 END
2094
2095CDECK ID>, BESSELI1.
2096 REAL FUNCTION BESSELI1(X)
2097 IMPLICIT NONE
2098
2099 REAL X,Y,AX
2100
2101 AX = ABS(X)
2102
2103 IF (AX .LT. 3.75) THEN
2104
2105 Y = (X/3.75)**2
2106
2107 BESSELI1 =
2108 & AX*(0.5+Y*(0.87890594+Y*(0.51498869+Y*(0.15084934
2109 & +Y*(0.2658733E-1+Y*(0.301532E-2+Y*0.32411E-3))))))
2110
2111
2112 ELSE
2113
2114 Y = 3.75/AX
2115
2116 BESSELI1 =
2117 & 0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1
2118 & -Y*0.420059E-2))
2119
2120 BESSELI1 =
2121 & 0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2
2122 & +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1))))
2123
2124 BESSELI1 = BESSELI1 *
2125 & EXP(AX)/SQRT(AX)
2126
2127 ENDIF
2128
2129 IF (X .LT. 0) THEN
2130 BESSELI1 = -BESSELI1
2131 ENDIF
2132
2133 RETURN
2134 END
2135
2136CDECK ID>, GAUSSDAT.
2137C DATA for the Gauss integration, x-values and weight for a
2138C N=2,4,8,16,32,64 point Gauss integration.
2139C
2140C based on program GAUSS64 of D. Trautmann
2141C data from Abramowitz & Stegun
2142C
2143C KAI HENCKEN, FEBRUAR 1993
2144C
2145 BLOCKDATA GAUSSDATA
2146 IMPLICIT NONE
2147 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
2148 COMMON /GAUSSD/XGAUSS,WGAUSS
2149
2150 DATA XGAUSS(1)/ .57735026918962576D0/
2151 DATA XGAUSS(2)/-.57735026918962576D0/
2152 DATA WGAUSS(1)/ 1.00000000000000000D0/
2153 DATA WGAUSS(2)/ 1.00000000000000000D0/
2154
2155 DATA XGAUSS(3)/ .33998104358485627D0/
2156 DATA XGAUSS(4)/ .86113631159405258D0/
2157 DATA XGAUSS(5)/-.33998104358485627D0/
2158 DATA XGAUSS(6)/-.86113631159405258D0/
2159 DATA WGAUSS(3)/ .65214515486254613D0/
2160 DATA WGAUSS(4)/ .34785484513745385D0/
2161 DATA WGAUSS(5)/ .65214515486254613D0/
2162 DATA WGAUSS(6)/ .34785484513745385D0/
2163
2164 DATA XGAUSS(7)/ .18343464249564981D0/
2165 DATA XGAUSS(8)/ .52553240991632899D0/
2166 DATA XGAUSS(9)/ .79666647741362674D0/
2167 DATA XGAUSS(10)/ .96028985649753623D0/
2168 DATA XGAUSS(11)/-.18343464249564981D0/
2169 DATA XGAUSS(12)/-.52553240991632899D0/
2170 DATA XGAUSS(13)/-.79666647741362674D0/
2171 DATA XGAUSS(14)/-.96028985649753623D0/
2172 DATA WGAUSS(7)/ .36268378337836198D0/
2173 DATA WGAUSS(8)/ .31370664587788727D0/
2174 DATA WGAUSS(9)/ .22238103445337448D0/
2175 DATA WGAUSS(10)/ .10122853629037627D0/
2176 DATA WGAUSS(11)/ .36268378337836198D0/
2177 DATA WGAUSS(12)/ .31370664587788727D0/
2178 DATA WGAUSS(13)/ .22238103445337448D0/
2179 DATA WGAUSS(14)/ .10122853629037627D0/
2180
2181 DATA XGAUSS(15)/ .0950125098376374402D0/
2182 DATA XGAUSS(16)/ .281603550779258913D0/
2183 DATA XGAUSS(17)/ .458016777657227386D0/
2184 DATA XGAUSS(18)/ .617876244402643748D0/
2185 DATA XGAUSS(19)/ .755404408355003034D0/
2186 DATA XGAUSS(20)/ .865631202387831744D0/
2187 DATA XGAUSS(21)/ .944575023073232576D0/
2188 DATA XGAUSS(22)/ .989400934991649933D0/
2189 DATA XGAUSS(23)/-.0950125098376374402D0/
2190 DATA XGAUSS(24)/-.281603550779258913D0/
2191 DATA XGAUSS(25)/-.458016777657227386D0/
2192 DATA XGAUSS(26)/-.617876244402643748D0/
2193 DATA XGAUSS(27)/-.755404408355003034D0/
2194 DATA XGAUSS(28)/-.865631202387831744D0/
2195 DATA XGAUSS(29)/-.944575023073232576D0/
2196 DATA XGAUSS(30)/-.989400934991649933D0/
2197 DATA WGAUSS(15)/ .189450610455068496D0/
2198 DATA WGAUSS(16)/ .182603415044923589D0/
2199 DATA WGAUSS(17)/ .169156519395002538D0/
2200 DATA WGAUSS(18)/ .149595988816576732D0/
2201 DATA WGAUSS(19)/ .124628971255533872D0/
2202 DATA WGAUSS(20)/ .0951585116824927848D0/
2203 DATA WGAUSS(21)/ .0622535239386478929D0/
2204 DATA WGAUSS(22)/ .0271524594117540949D0/
2205 DATA WGAUSS(23)/ .189450610455068496D0/
2206 DATA WGAUSS(24)/ .182603415044923589D0/
2207 DATA WGAUSS(25)/ .169156519395002538D0/
2208 DATA WGAUSS(26)/ .149595988816576732D0/
2209 DATA WGAUSS(27)/ .124628971255533872D0/
2210 DATA WGAUSS(28)/ .0951585116824927848D0/
2211 DATA WGAUSS(29)/ .0622535239386478929D0/
2212 DATA WGAUSS(30)/ .0271524594117540949D0/
2213
2214 DATA XGAUSS(31)/ .0483076656877383162D0/
2215 DATA XGAUSS(32)/ .144471961582796493D0/
2216 DATA XGAUSS(33)/ .239287362252137075D0/
2217 DATA XGAUSS(34)/ .331868602282127650D0/
2218 DATA XGAUSS(35)/ .421351276130635345D0/
2219 DATA XGAUSS(36)/ .506899908932229390D0/
2220 DATA XGAUSS(37)/ .587715757240762329D0/
2221 DATA XGAUSS(38)/ .663044266930215201D0/
2222 DATA XGAUSS(39)/ .732182118740289680D0/
2223 DATA XGAUSS(40)/ .794483795967942407D0/
2224 DATA XGAUSS(41)/ .849367613732569970D0/
2225 DATA XGAUSS(42)/ .896321155766052124D0/
2226 DATA XGAUSS(43)/ .934906075937739689D0/
2227 DATA XGAUSS(44)/ .964762255587506430D0/
2228 DATA XGAUSS(45)/ .985611511545268335D0/
2229 DATA XGAUSS(46)/ .997263861849481564D0/
2230 DATA XGAUSS(47)/-.0483076656877383162D0/
2231 DATA XGAUSS(48)/-.144471961582796493D0/
2232 DATA XGAUSS(49)/-.239287362252137075D0/
2233 DATA XGAUSS(50)/-.331868602282127650D0/
2234 DATA XGAUSS(51)/-.421351276130635345D0/
2235 DATA XGAUSS(52)/-.506899908932229390D0/
2236 DATA XGAUSS(53)/-.587715757240762329D0/
2237 DATA XGAUSS(54)/-.663044266930215201D0/
2238 DATA XGAUSS(55)/-.732182118740289680D0/
2239 DATA XGAUSS(56)/-.794483795967942407D0/
2240 DATA XGAUSS(57)/-.849367613732569970D0/
2241 DATA XGAUSS(58)/-.896321155766052124D0/
2242 DATA XGAUSS(59)/-.934906075937739689D0/
2243 DATA XGAUSS(60)/-.964762255587506430D0/
2244 DATA XGAUSS(61)/-.985611511545268335D0/
2245 DATA XGAUSS(62)/-.997263861849481564D0/
2246 DATA WGAUSS(31)/ .0965400885147278006D0/
2247 DATA WGAUSS(32)/ .0956387200792748594D0/
2248 DATA WGAUSS(33)/ .0938443990808045654D0/
2249 DATA WGAUSS(34)/ .0911738786957638847D0/
2250 DATA WGAUSS(35)/ .0876520930044038111D0/
2251 DATA WGAUSS(36)/ .0833119242269467552D0/
2252 DATA WGAUSS(37)/ .0781938957870703065D0/
2253 DATA WGAUSS(38)/ .0723457941088485062D0/
2254 DATA WGAUSS(39)/ .0658222227763618468D0/
2255 DATA WGAUSS(40)/ .0586840934785355471D0/
2256 DATA WGAUSS(41)/ .0509980592623761762D0/
2257 DATA WGAUSS(42)/ .0428358980222266807D0/
2258 DATA WGAUSS(43)/ .0342738629130214331D0/
2259 DATA WGAUSS(44)/ .0253920653092620595D0/
2260 DATA WGAUSS(45)/ .0162743947309056706D0/
2261 DATA WGAUSS(46)/ .00701861000947009660D0/
2262 DATA WGAUSS(47)/ .0965400885147278006D0/
2263 DATA WGAUSS(48)/ .0956387200792748594D0/
2264 DATA WGAUSS(49)/ .0938443990808045654D0/
2265 DATA WGAUSS(50)/ .0911738786957638847D0/
2266 DATA WGAUSS(51)/ .0876520930044038111D0/
2267 DATA WGAUSS(52)/ .0833119242269467552D0/
2268 DATA WGAUSS(53)/ .0781938957870703065D0/
2269 DATA WGAUSS(54)/ .0723457941088485062D0/
2270 DATA WGAUSS(55)/ .0658222227763618468D0/
2271 DATA WGAUSS(56)/ .0586840934785355471D0/
2272 DATA WGAUSS(57)/ .0509980592623761762D0/
2273 DATA WGAUSS(58)/ .0428358980222266807D0/
2274 DATA WGAUSS(59)/ .0342738629130214331D0/
2275 DATA WGAUSS(60)/ .0253920653092620595D0/
2276 DATA WGAUSS(61)/ .0162743947309056706D0/
2277 DATA WGAUSS(62)/ .00701861000947009660D0/
2278
2279 DATA XGAUSS(63)/ .02435029266342443250D0/
2280 DATA XGAUSS(64)/ .0729931217877990394D0/
2281 DATA XGAUSS(65)/ .121462819296120554D0/
2282 DATA XGAUSS(66)/ .169644420423992818D0/
2283 DATA XGAUSS(67)/ .217423643740007084D0/
2284 DATA XGAUSS(68)/ .264687162208767416D0/
2285 DATA XGAUSS(69)/ .311322871990210956D0/
2286 DATA XGAUSS(70)/ .357220158337668116D0/
2287 DATA XGAUSS(71)/ .402270157963991604D0/
2288 DATA XGAUSS(72)/ .446366017253464088D0/
2289 DATA XGAUSS(73)/ .489403145707052957D0/
2290 DATA XGAUSS(74)/ .531279464019894546D0/
2291 DATA XGAUSS(75)/ .571895646202634034D0/
2292 DATA XGAUSS(76)/ .611155355172393250D0/
2293 DATA XGAUSS(77)/ .648965471254657340D0/
2294 DATA XGAUSS(78)/ .685236313054233243D0/
2295 DATA XGAUSS(79)/ .719881850171610827D0/
2296 DATA XGAUSS(80)/ .752819907260531897D0/
2297 DATA XGAUSS(81)/ .783972358943341408D0/
2298 DATA XGAUSS(82)/ .813265315122797560D0/
2299 DATA XGAUSS(83)/ .840629296252580363D0/
2300 DATA XGAUSS(84)/ .865999398154092820D0/
2301 DATA XGAUSS(85)/ .889315445995114106D0/
2302 DATA XGAUSS(86)/ .910522137078502806D0/
2303 DATA XGAUSS(87)/ .929569172131939576D0/
2304 DATA XGAUSS(88)/ .946411374858402816D0/
2305 DATA XGAUSS(89)/ .961008799652053719D0/
2306 DATA XGAUSS(90)/ .973326827789910964D0/
2307 DATA XGAUSS(91)/ .983336253884625957D0/
2308 DATA XGAUSS(92)/ .991013371476744321D0/
2309 DATA XGAUSS(93)/ .996340116771955279D0/
2310 DATA XGAUSS(94)/ .999305041735772139D0/
2311 DATA XGAUSS(95)/-.02435029266342443250D0/
2312 DATA XGAUSS(96)/-.0729931217877990394D0/
2313 DATA XGAUSS(97)/-.121462819296120554D0/
2314 DATA XGAUSS(98)/-.169644420423992818D0/
2315 DATA XGAUSS(99)/-.217423643740007084D0/
2316 DATA XGAUSS(100)/-.264687162208767416D0/
2317 DATA XGAUSS(101)/-.311322871990210956D0/
2318 DATA XGAUSS(102)/-.357220158337668116D0/
2319 DATA XGAUSS(103)/-.402270157963991604D0/
2320 DATA XGAUSS(104)/-.446366017253464088D0/
2321 DATA XGAUSS(105)/-.489403145707052957D0/
2322 DATA XGAUSS(106)/-.531279464019894546D0/
2323 DATA XGAUSS(107)/-.571895646202634034D0/
2324 DATA XGAUSS(108)/-.611155355172393250D0/
2325 DATA XGAUSS(109)/-.648965471254657340D0/
2326 DATA XGAUSS(110)/-.685236313054233243D0/
2327 DATA XGAUSS(111)/-.719881850171610827D0/
2328 DATA XGAUSS(112)/-.752819907260531897D0/
2329 DATA XGAUSS(113)/-.783972358943341408D0/
2330 DATA XGAUSS(114)/-.813265315122797560D0/
2331 DATA XGAUSS(115)/-.840629296252580363D0/
2332 DATA XGAUSS(116)/-.865999398154092820D0/
2333 DATA XGAUSS(117)/-.889315445995114106D0/
2334 DATA XGAUSS(118)/-.910522137078502806D0/
2335 DATA XGAUSS(119)/-.929569172131939576D0/
2336 DATA XGAUSS(120)/-.946411374858402816D0/
2337 DATA XGAUSS(121)/-.961008799652053719D0/
2338 DATA XGAUSS(122)/-.973326827789910964D0/
2339 DATA XGAUSS(123)/-.983336253884625957D0/
2340 DATA XGAUSS(124)/-.991013371476744321D0/
2341 DATA XGAUSS(125)/-.996340116771955279D0/
2342 DATA XGAUSS(126)/-.999305041735772139D0/
2343 DATA WGAUSS(63)/ .0486909570091397204D0/
2344 DATA WGAUSS(64)/ .0485754674415034269D0/
2345 DATA WGAUSS(65)/ .0483447622348029572D0/
2346 DATA WGAUSS(66)/ .0479993885964583077D0/
2347 DATA WGAUSS(67)/ .0475401657148303087D0/
2348 DATA WGAUSS(68)/ .0469681828162100173D0/
2349 DATA WGAUSS(69)/ .0462847965813144172D0/
2350 DATA WGAUSS(70)/ .0454916279274181445D0/
2351 DATA WGAUSS(71)/ .0445905581637565631D0/
2352 DATA WGAUSS(72)/ .0435837245293234534D0/
2353 DATA WGAUSS(73)/ .0424735151236535890D0/
2354 DATA WGAUSS(74)/ .0412625632426235286D0/
2355 DATA WGAUSS(75)/ .0399537411327203414D0/
2356 DATA WGAUSS(76)/ .0385501531786156291D0/
2357 DATA WGAUSS(77)/ .0370551285402400460D0/
2358 DATA WGAUSS(78)/ .0354722132568823838D0/
2359 DATA WGAUSS(79)/ .0338051618371416094D0/
2360 DATA WGAUSS(80)/ .0320579283548515535D0/
2361 DATA WGAUSS(81)/ .0302346570724024789D0/
2362 DATA WGAUSS(82)/ .0283396726142594832D0/
2363 DATA WGAUSS(83)/ .0263774697150546587D0/
2364 DATA WGAUSS(84)/ .0243527025687108733D0/
2365 DATA WGAUSS(85)/ .0222701738083832542D0/
2366 DATA WGAUSS(86)/ .0201348231535302094D0/
2367 DATA WGAUSS(87)/ .0179517157756973431D0/
2368 DATA WGAUSS(88)/ .0157260304760247193D0/
2369 DATA WGAUSS(89)/ .0134630478967186426D0/
2370 DATA WGAUSS(90)/ .0111681394601311288D0/
2371 DATA WGAUSS(91)/ .00884675982636394772D0/
2372 DATA WGAUSS(92)/ .00650445796897836286D0/
2373 DATA WGAUSS(93)/ .00414703326056246764D0/
2374 DATA WGAUSS(94)/ .00178328072169643295D0/
2375 DATA WGAUSS(95)/ .0486909570091397204D0/
2376 DATA WGAUSS(96)/ .0485754674415034269D0/
2377 DATA WGAUSS(97)/ .0483447622348029572D0/
2378 DATA WGAUSS(98)/ .0479993885964583077D0/
2379 DATA WGAUSS(99)/ .0475401657148303087D0/
2380 DATA WGAUSS(100)/ .0469681828162100173D0/
2381 DATA WGAUSS(101)/ .0462847965813144172D0/
2382 DATA WGAUSS(102)/ .0454916279274181445D0/
2383 DATA WGAUSS(103)/ .0445905581637565631D0/
2384 DATA WGAUSS(104)/ .0435837245293234534D0/
2385 DATA WGAUSS(105)/ .0424735151236535890D0/
2386 DATA WGAUSS(106)/ .0412625632426235286D0/
2387 DATA WGAUSS(107)/ .0399537411327203414D0/
2388 DATA WGAUSS(108)/ .0385501531786156291D0/
2389 DATA WGAUSS(109)/ .0370551285402400460D0/
2390 DATA WGAUSS(110)/ .0354722132568823838D0/
2391 DATA WGAUSS(111)/ .0338051618371416094D0/
2392 DATA WGAUSS(112)/ .0320579283548515535D0/
2393 DATA WGAUSS(113)/ .0302346570724024789D0/
2394 DATA WGAUSS(114)/ .0283396726142594832D0/
2395 DATA WGAUSS(115)/ .0263774697150546587D0/
2396 DATA WGAUSS(116)/ .0243527025687108733D0/
2397 DATA WGAUSS(117)/ .0222701738083832542D0/
2398 DATA WGAUSS(118)/ .0201348231535302094D0/
2399 DATA WGAUSS(119)/ .0179517157756973431D0/
2400 DATA WGAUSS(120)/ .0157260304760247193D0/
2401 DATA WGAUSS(121)/ .0134630478967186426D0/
2402 DATA WGAUSS(122)/ .0111681394601311288D0/
2403 DATA WGAUSS(123)/ .00884675982636394772D0/
2404 DATA WGAUSS(124)/ .00650445796897836286D0/
2405 DATA WGAUSS(125)/ .00414703326056246764D0/
2406 DATA WGAUSS(126)/ .00178328072169643295D0/
2407
2408 END
2409*=======================================================================
2410 SUBROUTINE LORENB (U,PS,PI,PF)
2411C
2412C CERN PROGLIB# U102 LORENB .VERSION KERNFOR 4.04 821124
2413C ORIG. 20/08/75 L.PAPE
2414C
2415 DOUBLE PRECISION PF4, FN
2416 DIMENSION PS(4),PI(4),PF(4)
2417
2418 IF (PS(4).EQ.U) GO TO 17
2419 PF4 = (PI(4)*PS(4)+PI(3)*PS(3)+PI(2)*PS(2)+PI(1)*PS(1)) / U
2420 FN = (PF4+PI(4)) / (PS(4)+U)
2421 PF(1)= PI(1) + FN*PS(1)
2422 PF(2)= PI(2) + FN*PS(2)
2423 PF(3)= PI(3) + FN*PS(3)
2424 PF(4)= PF4
2425 GO TO 18
2426C
2427 17 PF(1)= PI(1)
2428 PF(2)= PI(2)
2429 PF(3)= PI(3)
2430 PF(4)= PI(4)
2431C
2432 18 CONTINUE
2433C
2434 RETURN
2435C
2436 END