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