2 *=======================================================================
5 *-----------------------------------------------------------------------
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)
14 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
17 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
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
27 COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
30 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
31 DOUBLE PRECISION PARP,PARI
33 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
36 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
37 DOUBLE PRECISION PARU,PARJ
39 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
40 DOUBLE PRECISION PMAS,PARF,VCKM
42 COMMON/PYDAT4/CHAF(500,2)
45 C MXSS = maximum number of modes
46 C NSSMOD = number of modes
47 C ISSMOD = initial particle
48 C JSSMOD = final particles
50 C BSSMOD = branching ratio
53 COMMON/SSMODE/NSSMOD,ISSMOD(MXSS),JSSMOD(5,MXSS),GSSMOD(MXSS)
55 INTEGER NSSMOD,ISSMOD,JSSMOD
59 C AMGLSS = gluino mass
60 C AMQKSS = squark mass
61 C AMLLSS = left-slepton mass
62 C AMLRSS = right-slepton mass
63 C AMNLSS = sneutrino mass
64 C TWOM1 = Higgsino mass = - mu
65 C RV2V1 = ratio v2/v1 of vev's
66 C AMTLSS,AMTRSS = left,right stop masses
67 C AMT1SS,AMT2SS = light,heavy stop masses
68 C AMBLSS,AMBRSS = left,right sbottom masses
69 C AMB1SS,AMB2SS = light,heavy sbottom masses
70 C AMZiSS = signed mass of Zi
71 C ZMIXSS = Zi mixing matrix
72 C AMWiSS = signed Wi mass
73 C GAMMAL,GAMMAR = Wi left, right mixing angles
74 C AMHL,AMHH,AMHA = neutral Higgs h0, H0, A0 masses
75 C AMHC = charged Higgs H+ mass
76 C ALFAH = Higgs mixing angle
77 C AAT = stop trilinear term
78 C THETAT = stop mixing angle
79 C AAB = sbottom trilinear term
80 C 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)
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
95 $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
96 $,AAB,THETAB,AAL,THETAL,AMGVSS
98 EQUIVALENCE (AMZISS(1),AMZ1SS)
100 REAL xpar(10),ypar(6)
107 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 * Check bounds on 2 photon mass
111 IF (iproc .EQ. 1) THEN
112 IF (amin .LE. 2.31) THEN
114 & ' %GGINIT: minimal 2-photon mass is too low, ',
119 ELSE IF (iproc .EQ. 2) THEN
120 kc = pycomp (kf_onium)
122 IF (kf_onium .EQ. 441) THEN
127 c. xggres = ggwid(kf_onium)
128 IF (xgtres .EQ. 0.) xgtres = 0.1
129 IF (xmres.LE.amin .OR. xmres.GE.amax) THEN
131 & ' %GGINIT: Xmres is out of the Luminosity mass range'
135 ELSE IF (iproc .EQ. 3) THEN
137 amprt(1) = pmas(kc,1)
138 IF (amin .LT. 2*amprt(1)) THEN
140 & ' %GGINIT: minimal 2-photon mass is too low, ',
141 & 'set to 2m(fermion)'
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
149 & ' %GGINIT: minimal 2-photon mass is too low, ',
154 ELSE IF (iproc .EQ. 5) THEN
155 WRITE (6,*) 'Process 5 is not implemented yet'
157 C CALL ssmssm (xm1, xm2, xmg, xms, xmtl, xmtr, xmll, xmlr, xmnl,
158 C & xtanb, xmha, xmu, xmt, xat, xmbr, xab, iallow)
159 C IF (iallow .NE. 0) THEN
160 C WRITE (6,*) 'Lightest neutralino is not LSP, ',
161 C & 'input other MSSM parameters'
164 C amprt(1) = abs (amw1ss)
165 C amprt(2) = abs (amz1ss)
166 C WRITE (6,'('' %GGINIT'','//
167 C & ''', M(chi+) = '',f5.1,'//
168 C & ''', M(chi0) = '',f5.1)') amprt(1), amprt(2)
170 C IF (amin .LT. 2*amprt(1)) THEN
172 C & ' %GGINIT: minimal 2-photon mass is too low, ',
173 C & 'set to 2M(chi+)'
177 ELSE IF (iproc .EQ. 6) THEN
179 amprt(1) = pmas(kc1,1) + pmas(kc1,3)
181 amprt(2) = pmas(kc2,1) + pmas(kc2,3)
182 IF (amin .LT. amprt(1)+amprt(2)) THEN
184 & ' %GGINIT: minimal 2-photon mass is too low, ',
185 & 'set to m(kv1)+m(kv2)'
186 amin = amprt(1)+amprt(2)
191 IF (amin .GT. amax) THEN
192 WRITE (6,*) ' %GGINIT: AMIN > AMAX, execution stopped'
196 * ----- TWOGAM Initialization -----
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
204 ipar(1)= 1 ! Unweighted events
205 ipar(2)= 0 ! Continuum generation
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
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
221 CALL twoini (6,ipar,xpar,ypar)
223 * The Luminosity function initialization
225 * ilumf = +1 ! Reading the Luminosity from file
226 * ilumf = -1 ! New calculation of the Luminosity function
228 CALL lumfun (ilumf,lumfil,ymin,ymax,ny,amin,amax,nmas)
230 * Integral of lum. function over the (M,y) region
232 hmas = (amax-amin)/nmas
236 f2 = dldmx(xm-hmas/2.,ny)
238 xlumint = xlumint +(f1 + 4.*f2 + f3)/6.
241 xlumint = xlumint * hmas
242 * ----- end of TWOGAM initialization -----
244 IF (iproc .EQ. 1) THEN
245 * ----- PYTHIA initialization -----
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)
261 CALL pyinit ('5MOM','gamma','gamma',0.0D0)
262 * ----- end of PYTHIA initialization -----
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
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
283 pmas(kc,1) = amprt(1)
286 * Neutralino storing in JETSET
292 pmas(kc,1) = amprt(2)
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)
300 IF (issmod(idc) .EQ. 39) THEN
301 IF (moddcy .EQ. 1) THEN
303 IF (jssmod(jmod,idc) .EQ. -12) THEN
305 xsbra = xsbra*xsbra*4
311 IF (jssmod(jmod,idc) .EQ. -12) THEN ! e nu_e
312 ssbr(1) = bssmod(idc)
314 ELSE IF (jssmod(jmod,idc) .EQ. -14) THEN ! mu nu_mu
315 ssbr(2) = bssmod(idc)
317 ELSE IF (jssmod(jmod,idc) .EQ. -16) THEN ! tau nu_tau
318 ssbr(3) = bssmod(idc)
320 ELSE IF (jssmod(jmod,idc) .EQ. -2) THEN ! u dbar
321 ssbr(4) = bssmod(idc)
323 ELSE IF (jssmod(jmod,idc) .EQ. 4) THEN ! c sbar
324 ssbr(5) = bssmod(idc)
334 ELSE IF (iproc .EQ. 6) THEN
335 * Fix parameters of process gg -> V1V2
336 IF (kv1 .GT. kv2) THEN
341 IF (kv1.EQ.113 .AND. kv2.EQ.113) THEN
343 ELSE IF (kv1.EQ.113 .AND. kv2.EQ.223) THEN
345 ELSE IF (kv1.EQ.113 .AND. kv2.EQ.333) THEN
347 ELSE IF (kv1.EQ.113 .AND. kv2.EQ.443) THEN
349 ELSE IF (kv1.EQ.223 .AND. kv2.EQ.223) THEN
351 ELSE IF (kv1.EQ.223 .AND. kv2.EQ.333) THEN
353 ELSE IF (kv1.EQ.223 .AND. kv2.EQ.443) THEN
355 ELSE IF (kv1.EQ.333 .AND. kv2.EQ.333) THEN
357 ELSE IF (kv1.EQ.333 .AND. kv2.EQ.443) THEN
359 ELSE IF (kv1.EQ.443 .AND. kv2.EQ.443) THEN
362 WRITE (6,*) ' %GGINIT: unknown (V1 V2) combination'
363 WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
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
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
380 xm = amin + ggrnd(0) * (amax-amin)
391 *=======================================================================
393 FUNCTION ggwid (kf_res)
394 *-----------------------------------------------------------------------
395 * Routine to define 2-gamma width of a resonance
396 * with JETSET code KF_RES
398 *-----------------------------------------------------------------------
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
403 ELSE IF (kf_res .EQ. 221) THEN
405 ELSE IF (kf_res .EQ. 331) THEN
407 ELSE IF (kf_res .EQ. 441
408 & .OR. kf_res .EQ. 10441
409 & .OR. kf_res .EQ. 445) THEN
411 ELSE IF (kf_res .EQ. 551
412 & .OR. kf_res .EQ. 10551
413 & .OR. kf_res .EQ. 555) THEN
419 *=======================================================================
422 *-----------------------------------------------------------------------
423 * Single event generation
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)
431 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
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,
438 COMMON /ggkin/ xkin(10)
439 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
442 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
443 DOUBLE PRECISION PARU,PARJ
445 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
446 DOUBLE PRECISION PMAS,PARF,VCKM
448 REAL pgam1(4),pgam2(4)
453 1 CALL gg2gam (amas,wsq,ygg,q1sq,q2sq,
454 & pgam1,pgam2,ptag1,ptag2,weight)
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
460 xssum = xssum + xscur
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
466 xstote = xstot / sqrt(float(ntry))
472 dbetaz = dtanh(1.d0*ygg)
474 * calculation of the 4-vector of two-photon system
475 * and filling arrays pgg(1:2,1:5), kgg(1:2)
477 p2g(ip) = pgam1(ip) + pgam2(ip)
478 pgg(1,ip) = pgam1(ip)
479 pgg(2,ip) = pgam2(ip)
488 * =======> PROCESS # 1
489 IF (iproc .EQ. 1) THEN
491 p(1,ip) = dble (pgam1(ip))
492 p(2,ip) = dble (pgam2(ip))
494 p(1,5) = dble (-xmg1)
495 p(2,5) = dble (-xmg2)
498 * =======> PROCESS # 2
499 ELSE IF (iproc .EQ. 2) THEN
502 pgg(ngg,ip) = p2g(ip)
505 * Fill the JETSET common block
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
513 IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR. iproc.EQ.5) THEN
521 ELSE IF (iproc.EQ.6) THEN
529 sinthe = sqrt ((1.-costhe)*(1.+costhe))
530 p3(1) = p3abs * cos(phi) * sinthe
531 p3(2) = p3abs * sin(phi) * sinthe
532 p3(3) = p3abs * costhe
541 pgg(igg1,ip) = p3(ip)
542 pgg(igg2,ip) = p4(ip)
548 IF (iproc .EQ. 3) THEN
552 ELSE IF (iproc .EQ. 4) THEN
556 ELSE IF (iproc .EQ. 5) THEN
559 * Decay of charginos into W+neutralino
566 ELSE IF (iproc.EQ.6) THEN
572 * Boost event into initial frame
574 CALL pyrobo(3,ngg, 0.d0,0.d0, 0.d0,0.d0, dbetaz)
582 *=======================================================================
585 *-----------------------------------------------------------------------
586 * Cross section calculation and end of run and
587 * print out cross section and related variables
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)
595 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
597 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
601 WRITE (6,'(/1x,79(''=''))')
602 WRITE (6,'(1x,''I'',36x,''TPHIC'',36x, ''I'')')
603 WRITE (6,'(1x,''I'',34x,''Process '',i1,34x,''I'')') iproc
604 WRITE (6,'(1x,''I'',34x,''---------'', 34x,''I'')')
605 WRITE (6,'(1x,''I'',77x, ''I'')')
606 WRITE (6,'(1x,''I'',20x,''Events thrown : '',i8,31x,''I'')')
608 WRITE (6,'(1x,''I'',20x,''Events accepted : '',i8,31x,''I'')')
610 WRITE (6,'(1x,''I'',20x,''Cross section : '',e10.3,'' +- '','//
611 & 'e10.3,'' nb'',12x ''I'')') xstot, xstote
612 WRITE (6,'(1x,''I'',77x, ''I'')')
613 WRITE (6,'(1x,79(''=''))')
615 IF (iproc .EQ. 1) THEN
616 * --- Alternative cross section for minimum bias events ---
618 hmas = (amax-amin)/nmas
619 f1 = dldmx(amin ,ny) * dchad(amin)
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.
627 sigtt = totxsc * hmas
628 sigmx = sigtt/sqrt(float(nevent))
629 WRITE (6,*) ' CrosSec =',sigtt,' +/-',sigmx,' mb'
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
647 *=======================================================================
650 *-----------------------------------------------------------------------
651 * Cross section calculation in current event for iproc=1,3,4,5,6
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)
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,
664 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
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
670 PARAMETER (epsilon = 0.0808, eta = -0.4525)
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)
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) -
688 xden = (1 - beta2*costh2)**2
690 ffmax = (1. + beta2) / (1. - beta2)
692 IF (ffcur .LE. rff*ffmax) GOTO 1
694 small = 1.0 - cos(thetamin)
695 2 zz = exp (ggrnd(0) * alog(2./small))
696 C sgn = sign (1., 2.*ggrnd(0)-1.)
699 costhe = sgn * (zz-1.) / (zz+1.)
700 ffcur = (1.+costhe*costhe) / (1.-beta2*costhe*costhe)
701 ffmax = 2. / (1.- costhe*costhe)
703 IF (ffcur .LE. rff*ffmax) GOTO 2
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))-
716 xscur0 = 2. * alog (sqrt(wsq)/amprt(1)) - 1.0
718 xscur0 = totfac * xscur0 * xsbra
719 xscur = xscur0 * xlumint
720 ELSE IF (iproc .EQ. 4) THEN
722 beta2 = 1. - 4.*xmw**2/wsq
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
731 ffmax = (19. + 10.*beta2 + beta2**2) / (1. - beta2)**2
733 IF (ffcur .LE. rff*ffmax) GOTO 3
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
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)
756 IF (fmax .GT. 1.e-10) THEN
757 tvar = 1./bb * dlog (fmax - ggrnd(0)*(fmax-fmin))
759 tvar = tmax + 1./bb * dlog (1. - ggrnd(0))
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
772 IF (kv1 .EQ. 443 .AND. kv2 .EQ. 443) THEN
773 factor = 1.0 / alog(xmgg/3.097)
777 xscur0 = factor * gvpar(1) * xmgg**gvpar(2) / bb * exp(bb*tmax)
778 xscur = xscur0 * xlumint
782 * Find max cross section of gg-process
784 IF (xscur0 .GT. xsmax0) xsmax0 = xscur0
789 *=======================================================================
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)
802 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
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,
809 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
810 DOUBLE PRECISION PMAS,PARF,VCKM
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)
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-amprt(2))*(eb+amprt(2)))
831 costhe = 2*ggrnd(0) - 1.
832 sinthe = sqrt((1.-costhe)*(1.+costhe))
834 pcm(1) = pb * sinthe * cos(phi)
835 pcm(2) = pb * sinthe * sin(phi)
838 CALL lorenb (amprt(1),p1,pcm,p2)
843 CALL lorenb (amprt(1),p1,pcm,p3)
845 pgg(ngg+1,ip) = p2(ip)
846 pgg(ngg+2,ip) = p3(ip)
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))
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
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))
868 * Random mass of (2+3) according to Breit-Wigner
870 xm23 = xmw**2 + xmw*gwt *
871 * tan (atmin + ggrnd(0)*(atmax - atmin))
872 xm23 = max (0., xm23)
874 IF (xm23 .GT. xm23max) GOTO 1 ! to avoid boundary effect
875 e23 = (am1**2 + xm23**2 - am4**2) / 2. / am1
877 * Rest frame of (2+3)
879 phi = 2.*pi * ggrnd(0)
880 costhe = 2.*ggrnd(0) - 1.
881 sinthe = sqrt ((1.-costhe)*(1.+costhe))
884 p2(1) = pp * sinthe * cos(phi)
885 p2(2) = pp * sinthe * sin(phi)
894 * Transformation from rest frame of (2+3) to rest frame of (1),
895 * 3-momenta p2, p3 are along z-axis
898 gambe = sqrt ((gamma-1.)*(gamma+1.))
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)
910 * Arbitrary rotation of rest frame of (1)
912 phi = 2.*pi * ggrnd(0)
913 costhe = 2.*ggrnd(0) - 1.
914 sinthe = sqrt ((1.-costhe)*(1.+costhe))
941 ptmp(ix) = ptmp(ix) + o1(ix,jx) * p2(jx)
950 ptmp(ix) = ptmp(ix) + o2(ix,jx) * p2(jx)
960 ptmp(ix) = ptmp(ix) + o1(ix,jx) * p3(jx)
969 ptmp(ix) = ptmp(ix) + o2(ix,jx) * p3(jx)
979 pp =sqrt ((e4-am4)*(e4+am4))
980 p4(1) = pp * sinthe * cos(phi)
981 p4(2) = pp * sinthe * sin(phi)
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
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)
999 * Select current phase space configuration
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
1006 * Pick up final lepton flavour
1008 IF (moddcy .EQ. 1) THEN
1009 IF (rflav .LE. 0.5) THEN
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
1025 ELSE IF (rflav .LE. rr2) THEN
1028 ELSE IF (rflav .LE. rr3) THEN
1031 ELSE IF (rflav .LE. rr4) THEN
1034 ELSE IF (rflav .LE. rr5) THEN
1042 * Transformation to CMS of 2-gamma
1043 CALL lorenb (amprt(1),p1,p2,pcm)
1045 pgg(ngg+1,ip) = pcm(ip)
1048 kgg(ngg+1) = lflav1 * sign(1.,1.*kgg(igg))
1051 CALL lorenb (amprt(1),p1,p3,pcm)
1053 pgg(ngg+1,ip) = pcm(ip)
1056 kgg(ngg+1) = lflav2 * sign(1.,1.*kgg(igg))
1059 CALL lorenb (amprt(1),p1,p4,pcm)
1061 pgg(ngg+1,ip) = pcm(ip)
1063 pgg(ngg+1,5) = amprt(2)
1064 kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
1071 *=======================================================================
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
1085 IF (igg .LE. 2) THEN
1092 p(igg,ip) = pgg(igg,ip)
1099 *=======================================================================
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
1116 *=======================================================================
1118 SUBROUTINE gg2gam (amas,wsq,y3,qs1,qs2,pgam1,pgam2,ptag1,ptag2,wt)
1119 C***********************************************************************
1121 C Generator of the Direct Two Photon Processes in S.A.Sadovsky
1122 C coherent heavy ion collisions: A+B --> 1+2+3 17.01.1995
1123 C with double differencial gg-luminosity function Yu.Kharlov
1126 C LUMFUN must be called first to initialise the generator
1128 C***********************************************************************
1130 C DETAILED DESCRIPTION
1131 C ======== ===========
1133 C INPUTS: only through common /TWOGENC/
1136 C OUTPUTS: AMAS The mass of ion.
1137 C ======== WSQ Square mass of the gamma-gamma system
1138 C Y3 The rapidity of the gamma-gamma system
1139 C QS1 Virtual mass squared of photon 1.
1140 C QS2 Virtual mass squared of photon 2.
1141 C PGAM1 Four-vector of photon 1
1142 C PGAM2 Four-vector of photon 2
1143 C PTAG1 Four-vector of tag 1
1144 C PTAG2 Four-vector of tag 2
1145 C WT The weight of the event (dummy so far)
1147 C The remaining variables are internal to the program.
1149 C***********************************************************************
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
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
1164 PARAMETER (alpha=1./137.036,pi=3.14159265)
1166 REAL ptag1(10),ptag2(10),pgam1(10),pgam2(10),wt
1168 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1169 REAL*4 y3,xm,q0,qs,qs1,qs2
1175 DATA q0,pi2,lstr / 0.060, 6.28318530718 d0, .false. /
1179 icnter(1) = icnter(1) + 1
1180 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1187 pa(3)= dsqrt(1D0*((eb-amas)*(eb+amas)))
1198 * Generate 2gamma rapidity y3 and mass xm
1199 CALL ranlum(y3,xm,wt)
1203 * Energy and z-momentum of 2gamma system
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?
1213 * Energies and z-momenta of the photons
1214 xfun12 = wsq + xm1 - xm2
1216 eg1 = (egg + pgg*sqrt(1 - 4*wsq*xm1/(xfun12*xfun12))) *
1219 pg1 = sqrt (eg1**2 - xm1)
1220 pg2 = -sqrt (eg2**2 - xm2)
1222 * 4-momenta of the photons
1232 * 4-momenta of the recoil ions
1234 ptag1(ip) = pa(ip) - pgam1(ip)
1235 ptag2(ip) = pb(ip) - pgam2(ip)
1241 *=======================================================================
1243 SUBROUTINE ranlum (ry,rm,wlum)
1246 C Corrected: Yu.Kharlov
1250 C generate two random numbers (ry and rm) with probability proportional
1251 C to luminosity function in ymin-ymax rapidity and xmin-xmax mass range
1253 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
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
1266 DATA ifl/0/, stpy/0./, stpm/0./
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'
1279 stpy = (ymx-ymn)/nny
1280 stpm = (amx-amn)/nnms
1281 rmin = exp(cexp*(xmax))
1282 rmax = exp(cexp*(xmin))
1289 rm = ggrnd(0) ! random mass
1291 rm = alog(rm)/cexp ! exponetial law
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
1299 ry = ggrnd(0) ! random Y
1301 * Find y-limits for generated gg-mass: y_max = log(E_max/M_gg)
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.
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
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
1322 *=======================================================================
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,
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)
1337 COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
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
1345 PARAMETER (al = 1./128.)
1346 DATA pi /3.14159276/, alpha /al/, gev2nb /389379./,
1348 DATA iproc, nevent, ilumf /1, 10000, -1/
1349 DATA lumfil /'gglum.dat'/
1351 DATA amprt /110.,10.,150.,-1.,0.106/
1352 DATA thetamin /1.e-03/
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
1374 *=======================================================================
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
1383 REAL ptag1(4),ptag2(4),p3(5)
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)
1391 ichg=pychge(k(i,2))/3
1392 WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5)
1396 100 FORMAT(i6,i3,3e16.8,2x,2i8)
1397 101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8)
1400 *=======================================================================
1403 SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR) TWOG0965
1404 C***********************************************************************TWOG0966
1406 C Initialises the gamma-gamma generator of Ion collisions TWOGAM. TWOG0968
1408 C INPUTS: LOUTX Logical unit for print-out TWOG0970
1409 C ======= XPAR(1) EBM Beam energy GeV TWOG0971
1410 C XPAR(2) WMN Minimum two-gamma mass GeV TWOG0972
1411 C XPAR(3) WMX Maximum two-gamma mass GeV TWOG0973
1413 C XPAR(8) WGHTMX Maximum weight TWOG0978
1414 C XPAR(9) XMRES Mass of resonance GeV TWOG0979
1415 C XPAR(10) XGRES Total width of resonance GeV TWOG0980
1417 C IPAR(1) 0 Unweighted events TWOG0982
1418 C 1 Weighted events TWOG0983
1419 C IPAR(2) 0 Continuum production TWOG0984
1420 C 1 Resonance formation TWOG0985
1422 C EXTERNALS: none TWOG0987
1423 C ========== TWOG0988
1425 C AUTHOR/DATE: S.A. Sadovsky, 17 January 1995 TWOG0990
1426 C =========== TWOG0991
1428 C***********************************************************************TWOG0993
1429 C 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
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
1447 REAL XPAR(11),YPAR(6) TWOG1012
1448 INTEGER LOUTX,IPAR(10),I TWOG1013
1450 EB = XPAR(1) TWOG1015
1451 WMIN = XPAR(2) TWOG1016
1452 WMAX = XPAR(3) TWOG1017
1453 WGHTMX= XPAR(8) TWOG1022
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
1462 C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2)) TWOG1031
1463 C rather than 1/2 COT(th/2) to be able to start theta from zero. TWOG1032
1465 S = (2.*EB)**2 TWOG1054
1466 OMMIN = WMIN**2/(4.*EB) TWOG1055
1467 OMRAT = EB/OMMIN TWOG1056
1469 C Print out initial conditions. TWOG1058
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
1476 WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX,
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
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
1508 IF (LWATE) THEN TWOG1090
1509 WRITE(LOUT,1001) TWOG1091
1510 WRITE(LOUT,1004) TWOG1092
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
1519 WRITE(LOUT,1001) TWOG1103
1520 WRITE(LOUT,1000) TWOG1104
1522 DO 10 I=1,NCNTRS TWOG1106
1523 10 ICNTER(I)=0 TWOG1107
1524 WGHTGT = 0. TWOG1108
1525 WGHTSM = 0. TWOG1109
1528 C=========================================================================
1530 FUNCTION DLDMX(Rm,Ny)
1531 PARAMETER (isize=100000)
1533 C Calculate the Luminosity integral over dY 5.02.1995
1534 C in the rapidity range ymin-ymax
1535 C using luminosity interpolation table
1536 C Ny/2 - number of subintervals for Simpson integration
1537 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1538 COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1539 + xmin,xmax,amn,amx,nnms,alfun(isize)
1541 stpm=(amx-amn)/nnms ! Interpolation steps
1544 q =(rm-im*stpm-amn)/stpm
1547 Hy =(ymax-ymin)/My ! Integration step over Y
1549 iy =(ry-ymn)/stpy+1.
1550 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 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)
1559 iy =(ry-ymn)/stpy+1.
1560 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 3
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)
1566 iy =(ry-ymn)/stpy+1.
1567 p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 2
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)
1572 SIM = SIM + (F1 + 4.*F2 + F3)/6.
1575 DLDMX = 1000.*SIM*HY ! dL/dMx in GeV^-1
1578 *=======================================================================
1580 SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms)
1581 C-----------------------------------------------------------------------
1582 C preparation of luminosity function for generator
1584 C ifl=-1 - calculate function + save in file 'file'
1585 C ifl= 0 - calculate function
1586 C ifl= 1 - read function from file
1587 C ymin,ymax,ny - rapidity range & number of points
1588 C xmn,xmx,nms - mass range & number of points
1589 C Attention: dimensions of alfun are alfun(ny+1,nms+1)
1593 C Author: G.V.Khaoustov 12.12.1994
1594 C Modified: Yu.Kharlov 18-JUN-1997
1595 C-----------------------------------------------------------------------
1596 COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1599 COMMON/D2LParam/lz,la,Rion,Gamma
1602 PARAMETER (isize=100000)
1605 COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny,
1606 + wmn,wmx,amn,amx,nnms,alfun(isize)
1608 umn = ymin ! parameters for luminosity generation
1613 C 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)
1618 * Fing max energy of gg-system (with safity factor 2)
1620 egg_max = 2.*hbarc * gamma / rion
1622 IF(ifl.EQ.1) THEN ! read lum. function from file
1623 OPEN(unit=10,err=1,file=file,status='old',
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))
1629 PRINT *,'The Luminosity has been read from File '
1633 IF (ny*nms .GT. isize) GOTO 2
1634 PRINT *,'Start the Luminosity calculation'
1643 stepm=(amx-amn)/nnms
1645 * Lum. function calculation
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)
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
1661 alfun(iy+(nny+1)*im) = DL
1666 DO i=1,nnms+1 ! find max values for exp. parameters calcul.
1669 IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1))
1674 s=(s-work(nnms+1))*stepm
1676 smax=0. ! calculation optimized parameters for exp. generator
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
1686 amp=am*exp(c*(k-1)*stepm)
1687 IF((work(k)-amp)/amp.GT.0.00001) THEN ! Accuracy problem
1698 PRINT *,'Fault to find best exponential approximation'
1702 IF(ifl.EQ.0) RETURN ! save lum. function
1703 OPEN(unit=10,err=1,file=file,status='unknown',
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))
1709 PRINT *,'The Luminosity calculation is finished'
1711 1 PRINT *,' file open error, file:',file
1713 2 PRINT *,' array "ALFUN" to small to contain luminosity function'
1715 3 PRINT *,' file read/writr error, file:',file
1718 99 FORMAT(2e13.5,i5,2e13.5,i5)
1720 C======================================================================
1722 REAL FUNCTION D2LDMDY(M,Y)
1723 C-----------------------------------------------------------------------
1724 C Calculate the Luminosity-function d^2L/dMdY
1725 C where M is the invariant mass and Y the rapidity
1726 C of the photon-photon-system
1728 C Calculations are based on the formulae of
1729 C N. Baron, G. Baur, Phys. Rev. C
1730 C C. Cahn, J. D. Jackson, Nucl. Phys.
1732 C Kai Hencken, 29. 9. 1994
1734 C M : invariant mass (MeV)
1737 C D2LDMDY : diff. Lum. (MeV**-1)
1740 C-----------------------------------------------------------------------
1742 COMMON/D2LParam/lz,la,Rion,Gamma
1747 COMMON/D2LIParam/RM,W1,W2
1750 EXTERNAL NClass,DeltaL
1753 C physical constants hbar*c and 1/alpha
1756 PARAMETER (hbarc=197.327053,falphai=137.0359895)
1760 C calculate photon frequencies from M,Y
1762 W1 = M/2.0 * EXP ( Y)
1763 W2 = M/2.0 * EXP (-Y)
1765 C calculate Rho in MeV^-1 instead of fm
1769 C calculate the "classical" value first.
1772 LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma)
1775 C substract from this the Delta-value
1779 D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL)
1783 *=======================================================================
1785 REAL FUNCTION NClass (W, R, Gamma)
1787 C the classical photon-number without impact parameter cutoff:
1793 PARAMETER (Pi = 3.141592)
1795 REAL BesselK0,BesselK1
1796 EXTERNAL BesselK0,BesselK1
1802 NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) -
1803 & Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2))
1807 *=======================================================================
1809 REAL FUNCTION DeltaL ()
1811 C the difference to the classical value
1815 COMMON/D2LParam/lz,la,Rion,Gamma
1819 COMMON/D2LIParam/RM,W1,W2
1821 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1822 COMMON /GAUSSD/XGAUSS,WGAUSS
1834 PARAMETER (Pi = 3.141592,Accur = 1D -2)
1837 C integrate first over b1
1841 C Loop incrementing the boundary
1849 C Loop for the Gauss integration
1855 DO I=2**N-1,2**(N+1)-2
1856 t = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0
1858 Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2
1860 Int = (tmax-tmin)/2.0*Int
1861 IF (Int .NE. 0.) THEN
1862 IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200
1865 IF (int .LT. 1.D-20) THEN
1868 WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE'
1873 IF (sum .NE. 0.0) THEN
1874 IF (ABS (Int2/Sum) .LT. Accur) GOTO 300
1876 IF (b1 .GT. 1.E+06) GOTO 300
1884 DeltaL = 4.0*Pi * Sum
1889 REAL FUNCTION Itgb1(b)
1891 C integration then over b2
1894 C Corrected by S.A.Sadovsky
1897 c +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997
1901 COMMON/D2LIParam/RM,W1,W2
1902 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1903 COMMON /GAUSSD/XGAUSS,WGAUSS
1909 PARAMETER (Accur = 1.D-2)
1913 IF (RM .GT. bmin) THEN
1916 bmax = b1 + 2.0 * RM
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
1925 Int = (bmax-bmin)/2.0*Int
1927 IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100
1930 IF (int .LT. 1.D-20) THEN
1933 WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int
1941 REAL*8 FUNCTION Itgb2 (b1,b2)
1943 C The function to be integrated over
1946 C Corrected by S.A.Sadovsky
1949 COMMON/D2LParam/lz,la,Rion,Gamma
1954 COMMON/D2LIParam/RM,W1,W2
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) *
1969 REAL FUNCTION Nphoton (W,Rho,Gamma)
1971 C The differential photonnumber for a nucleus
1972 C (without form factor)
1983 PARAMETER (PI=3.141592)
1987 C without form factor
1989 Wphib = WGamma*BESSELK1 (WGamma*Rho)
1991 Nphoton = 1.0/PI**2 * Wphib**2
1996 CDECK ID>, BESSELK0.
1997 REAL FUNCTION BESSELK0(X)
1999 C Special functions I0,K0,I1,K1
2000 C see Abramowitz&Stegun
2008 IF (X .LT. 2.0) THEN
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))))))
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))))))
2032 CDECK ID>, BESSELI0.
2033 REAL FUNCTION BESSELI0(X)
2040 IF (AX .LT. 3.75) THEN
2045 & 1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492
2046 & +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2)))))
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))))))))
2064 CDECK ID>, BESSELK1.
2065 REAL FUNCTION BESSELK1(X)
2073 IF (X .LT. 2.0) THEN
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)))))))
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)))))))
2094 CDECK ID>, BESSELI1.
2095 REAL FUNCTION BESSELI1(X)
2102 IF (AX .LT. 3.75) THEN
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))))))
2116 & 0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1
2120 & 0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2
2121 & +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1))))
2123 BESSELI1 = BESSELI1 *
2129 BESSELI1 = -BESSELI1
2135 CDECK ID>, GAUSSDAT.
2136 C DATA for the Gauss integration, x-values and weight for a
2137 C N=2,4,8,16,32,64 point Gauss integration.
2139 C based on program GAUSS64 of D. Trautmann
2140 C data from Abramowitz & Stegun
2142 C KAI HENCKEN, FEBRUAR 1993
2146 DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
2147 COMMON /GAUSSD/XGAUSS,WGAUSS
2149 DATA XGAUSS(1)/ .57735026918962576D0/
2150 DATA XGAUSS(2)/-.57735026918962576D0/
2151 DATA WGAUSS(1)/ 1.00000000000000000D0/
2152 DATA WGAUSS(2)/ 1.00000000000000000D0/
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/
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/
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/
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/
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/
2408 *=======================================================================
2409 SUBROUTINE LORENB (U,PS,PI,PF)
2411 C CERN PROGLIB# U102 LORENB .VERSION KERNFOR 4.04 821124
2412 C ORIG. 20/08/75 L.PAPE
2414 DOUBLE PRECISION PF4, FN
2415 DIMENSION PS(4),PI(4),PF(4)
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)