]> git.uio.no Git - u/mrichter/AliRoot.git/blame - DIME/dimemcv1.05.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / DIME / dimemcv1.05.f
CommitLineData
de4b745b 1 subroutine dimeinit
2
3 implicit double precision(a-y)
4 implicit complex*16(z)
5 double precision pt1(2),pt2(2),ptx(2),q(4,20)
6 double precision pboo(4),pcm(4),plb(4)
7 double precision tvec(4),uvec(4),svec(4),vv1(4),vv2(4)
8 integer d,h,i,j,k,l,m,n,o,p,r,iset,af,nhist,ntotal,eflag,ichi,
9 &icut,ncut,id(20),id0,id1,id2,pdf,num,ptgen,mode,iord,ii,ip,ll,nev
10 &,hh,lmax
11 character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
12 &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
13 integer iin,nch
14 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
15 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
16 common/ipars/nch
17 common/mom/q
18 common/int/ip
19 common/mompt/pt1,pt2,ptx
20 common/vars/s,rts,mmes,yx
21ccc new
22 common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp,
23 & mwidth, pi, rt2, ebeam, sum, sum1, weightm,
24 & bjac, bp
25 common /ivars/ ncut, ll, icut, nev, num
26ccc new
27 common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
28 common/flags/iin, pflag, fsi, ppbar, output, cuts, unw
29 common/gluons/g1,g2,kt1,kt2
30 common/levicivita/epsilon
31 common/ang/cost,costa,costb
32 common/alph/alfas
33 common/hvars/sh,th,uh
34 common/wvars/sig0,bb,t1,t2
35 common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
36 common/sec/cpom,cf,crho,aff,ar
37 integer iii
38 common/ivar/iii
39 common/ff/formf
40 common/ffpars/bexp,ao,bo,aoo,boo
41 common/regge/mregge
42ccccc hepevt output
43 integer nmxhep,kk
44 parameter (nmxhep=4000)
45 integer nevhep,nhep,isthep,idhep,jmohep,jdahep
46 double precision phep,vhep
47 common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
48 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
49ccccc Les Houches Event Common Block
50 INTEGER MAXNUP
51 PARAMETER (MAXNUP=500)
52 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
53 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
54 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
55 & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
56 & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
57 & SPINUP(MAXNUP)
58ccccc rambo variables
59 integer npart,nrun
60 double precision pin(4),am(100),pout(4,100),wt,ein
61 common/ramboc/ pin, am, pout, wt, ein,
62 &npart, nrun
63ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
64c c
65c Dime MC for the central exclusive production of c
66c meson pairs by double Pomeron exchange: c
67c c
68c p(1) p(2) --> p(3) + M_1(5) M_2(6) + p(4) c
69c c
70c Momenta for each event in array q(i,j), where j is c
71c the particle label and i is the 4-momentum c
72c component, with: c
73c c
74c 1,2 = transverse components c
75c 3 = beam component c
76c 4 = energy c
77c c
78c PDG number for ith particle in arrary ID(i) c
79c c
80c Also gives event record in HEPEVT or LHE format c
81c (others are available upon request) c
82c c
83ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
84c c
85c Particles generated: c
86c c
87c pflag='pipm' : M_1=pi^+ M_2=pi^- c
88c pflag='pi0' : M_1=M_2=pi^0 c
89c pflag='kpkm' : M_1=K^+ M_2=K^- c
90c pflag='ks' : M_1=M_2=K_0 c
91c pflag='rho' : M_1=M_2=rho_0 c
92c c
93c with decay: rho(5) --> pi^+(7)pi^-(8) c
94c rho(6) --> pi^+(9)pi^-(10) c
95c according to phase space, with finite rho c
96c width included c
97c c
98ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
99c c
100c User defined cuts can readily be implemented c
101c in subroutine 'icut' c
102c c
103ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
104c c
105c This is version 1.05 : 17 Sep 2014 c
106c c
107c Comments etc to: lucian.harland-lang@durham.ac.uk c
108c c
109c If you use the code please reference (and see for c
110c details of model used): c
111c c
112c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c
113c and W.J. Stirling arXiv:1105.1626 c
114c c
115c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c
116c and W.J. Stirling arXiv:1204.4803 c
117c c
118c L.A. Harland-Lang, V.A. Khoze and M.G. Ryskin c
119c arXiv:1312.4553 c
120c c
121ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
122
123ccAM open(35,file='exerec1.dat') ! event record
124
125ccAM rts=1.96d3 ! centre of mass energy
126 write(6,*) "energy", rts
127cccc Some basic cuts, imposed in subtroutine 'icut'. Other user defined cuts can readily be implemented in subroutine
128cccc note: in rhorho case these cuts are imposed on rho's, and not their decay productions. Cuts on the decay products
129cccc can be imposed by editing 'icut'
130
131ccAM rmax=1.8d1 ! max meson rapidity
132ccAM rmin=-1.8d1 ! min meson rapidity
133ccAM ecut=2d0 ! min meson p_t
134
135c -- physics parameters
136
137ccAM pflag='rho' ! Process generated - see preamble for options
138ccAM fsi='true' ! phenomenological model for exclusive suppression in Pom Pom --> meson pair. To turn on/off --> 'true/false'
139ccAM formf='orexp' ! meson - pomeron form factor.
140ccAM ppbar='false' ! set true if ppbar collisions
141 output='lhe' ! Event record style, HEPEVT or LHE
142ccAM cuts='true' ! Impose cuts or not
143ccAM unw='true' ! Set = 'true' for unweighted events
144ccAM iin=1 ! Model for soft survival factor, as described in arXiv:1306.2149. Default = 1
145
146cccccccc
147
148ccAM ntotal=1000000 ! no. of runs for weighted events
149ccAM nev=100 ! no. of unweighted events generated to event record
150
151ccccccccc
152
153 if(formf.eq.'exp')then ! Set parameters for Pomeron-Meson form factor in function 'fpi(x)'
154 bexp=1d0/2.2d0
155 elseif(formf.eq.'orexp')then
156 bo=1d0/1.1d0
157 ao=dsqrt(0.5d0)
158 elseif(formf.eq.'power')then
159 aoo=1.7d0
160 endif
161
162ccccccccccccccccccccccccccccccccccccccc
163
164ccccccccc Start of main body of code
165
166cccccccccccccccccccccccccccccccccccccc
167
168 call initpars(iin) ! Initialise soft survival parameters
169 call calcop ! proton opacity
170 call calcscreen ! screening amplitude
171
172 if(pflag.eq.'pipm')then
173 mmes=0.13957018d0 ! pi+/- mass, PDG 2011 value
174 sig0=13.63d0
175 alphapm=0.7d0
176 alpha0m=-0.7d0*mmes**2
177 cf=31.79d0/0.389d0
178 crho=4.23d0/0.389d0
179
180 nup=6
181 istup(5)=1
182 idup(5)=211
183 mothup(1,5)=1
184 mothup(2,5)=2
185 icolup(1,5)=0
186 icolup(2,5)=0
187 vtimup(5)=0
188 spinup(5)=9.
189
190 istup(6)=1
191 idup(6)=-211
192 mothup(1,6)=1
193 mothup(2,6)=2
194 icolup(1,6)=0
195 icolup(2,6)=0
196 vtimup(6)=0
197 spinup(6)=9.
198
199 elseif(pflag.eq.'pi0')then
200 mmes=0.1349766d0 ! pi0 mass, PDG 2011 value
201 sig0=13.63d0
202 alphapm=0.7d0
203 alpha0m=-0.7d0*mmes**2
204 cf=31.79d0/0.389d0
205 crho=0d0
206
207 nup=6
208 istup(5)=1
209 idup(5)=111
210 mothup(1,5)=1
211 mothup(2,5)=2
212 icolup(1,5)=0
213 icolup(2,5)=0
214 vtimup(5)=0
215 spinup(5)=9.
216
217 istup(6)=1
218 idup(6)=-111
219 mothup(1,6)=1
220 mothup(2,6)=2
221 icolup(1,6)=0
222 icolup(2,6)=0
223 vtimup(6)=0
224 spinup(6)=9.
225
226 elseif(pflag.eq.'kpkm')then
227 mmes=0.493677d0 ! K+/- mass, PDG 2011 value
228 sig0=11.82d0
229 cf=17.255d0/0.389d0
230 crho=9.105d0/0.389d0
231
232 nup=6
233 istup(5)=1
234 idup(5)=321
235 mothup(1,5)=1
236 mothup(2,5)=2
237 icolup(1,5)=0
238 icolup(2,5)=0
239 vtimup(5)=0
240 spinup(5)=9.
241
242 istup(6)=1
243 idup(6)=-321
244 mothup(1,6)=1
245 mothup(2,6)=2
246 icolup(1,6)=0
247 icolup(2,6)=0
248 vtimup(6)=0
249 spinup(6)=9.
250
251 elseif(pflag.eq.'ks')then
252 mmes=0.497614d0 ! K_0 mass, PDG 2011 value
253 sig0=11.82d0
254 cf=17.255d0/0.389d0
255 crho=0d0
256
257 nup=6
258 istup(5)=1
259 idup(5)=310
260 mothup(1,5)=1
261 mothup(2,5)=2
262 icolup(1,5)=0
263 icolup(2,5)=0
264 vtimup(5)=0
265 spinup(5)=9.
266
267 istup(6)=1
268 idup(6)=310
269 mothup(1,6)=1
270 mothup(2,6)=2
271 icolup(1,6)=0
272 icolup(2,6)=0
273 vtimup(6)=0
274 spinup(6)=9.
275
276 elseif(pflag.eq.'rho')then
277 mmes0=0.77549d0 ! rho mass, PDG 2013 value
278 mwidth=0.1491d0 ! rho width PDG 2013 value
279 sig0=10d0
280 cf=0d0
281 crho=0d0
282
283 nup=10
284 istup(5)=2
285 idup(5)=113
286 mothup(1,5)=1
287 mothup(2,5)=2
288 icolup(1,5)=0
289 icolup(2,5)=0
290 vtimup(5)=0
291 spinup(5)=9.
292
293 istup(6)=2
294 idup(6)=113
295 mothup(1,6)=1
296 mothup(2,6)=2
297 icolup(1,6)=0
298 icolup(2,6)=0
299 vtimup(6)=0
300 spinup(6)=9
301
302 do k=7,10
303 istup(k)=1
304 mothup(2,k)=0
305 icolup(1,k)=0
306 icolup(2,k)=0
307 vtimup(k)=0
308 spinup(k)=9.
309 enddo
310
311 idup(7)=211
312 idup(8)=-211
313 idup(9)=211
314 idup(10)=-211
315 mothup(1,7)=5
316 mothup(1,8)=5
317 mothup(1,9)=6
318 mothup(1,10)=6
319
320 endif
321
322c -- other parameters
323
324 ebeam=rts/2d0
325 s=rts*rts
326 zi=(0d0,1d0)
327 rt2=dsqrt(2d0)
328 pi=dacos(-1d0)
329 bp=rts/dsqrt(2.0d0)
330 mp=0.93827d0
331 beta=dsqrt(1d0-4d0*mp**2/s)
332 s0=1d0
333
334c Pomeron + t-slope
335
336 bb=4d0
337 bjac=6d0
338 bjac1=2d0
339
340 alphap=0.25d0 ! D-L '92 fit
341 alpha0=1.0808d0
342 alphapr=0.93d0
343 alpha0r=0.55d0
344
345 mf127=1.275d0
346 mf1525=1.525d0
347
348 cpom=sig0/0.389d0
349 aff=-0.860895d0
350 ar=-1.16158d0
351
352 mmes1=mmes
353 mmes2=mmes
354
355ccccc Initialise rambo (rho0 decay)
356
357 if(pflag.eq.'rho')then
358
359 mmes=mmes0
360 npart=2
361 do j=1,4
362 am(j)=0.13957018d0
363 enddo
364 endif
365
366ccccccc
367
368c initialise counters etc
369
370c do hh=1,20
371
372 nhist=1
373 sum=0d0
374 sum1=0d0
375 ncut=0
376
377 weightm=0d0
378
379c initialise histograms
380
381 do i=1,nhist
382 call histo3(i)
383 enddo
384
385 num=0
386
387cccccccccccccccccccccccccccccccccccccccccccccccccccccc
388
389c Incoming protons to event record array
390
391cccccccccccccccccccccccccccccccccccccccccccccccccccccc
392
393 ID(1)=2212
394 q(1,1)=0d0
395 q(2,1)=0d0
396 q(3,1)=ebeam*beta
397 q(4,1)=ebeam
398
399 istup(1)=-1
400 idup(1)=2212
401 mothup(1,1)=0
402 mothup(2,1)=0
403 icolup(1,1)=0
404 icolup(2,1)=0
405 do j=1,4
406 pup(j,1)=q(j,1)
407 enddo
408 pup(5,1)=dsqrt(q(4,1)**2-q(3,1)**2-q(2,1)**2-q(1,1)**2)
409 vtimup(1)=0
410 spinup(1)=9.
411
412 q(1,2)=0d0
413 q(2,2)=0d0
414 q(3,2)=-ebeam*beta
415 q(4,2)=ebeam
416
417 istup(2)=-1
418 if(ppbar.eq.'true')then
419 idup(2)=-2212
420 else
421 idup(2)=2212
422 endif
423 mothup(1,2)=0
424 mothup(2,2)=0
425 icolup(1,2)=0
426 icolup(2,2)=0
427 do j=1,4
428 pup(j,2)=q(j,2)
429 enddo
430 pup(5,2)=dsqrt(q(4,2)**2-q(3,2)**2-q(2,2)**2-q(1,2)**2)
431 vtimup(2)=0
432 spinup(2)=9.
433
434cccc outgoing initial info
435
436 istup(3)=1
437 if(ppbar.eq.'true')then
438 idup(2)=-2212
439 else
440 idup(2)=2212
441 endif
442 idup(3)=2212
443 mothup(1,3)=1
444 mothup(2,3)=2
445 icolup(1,3)=0
446 icolup(2,3)=0
447 vtimup(3)=0
448 spinup(3)=9
449
450 istup(4)=1
451 idup(4)=2212
452 mothup(1,4)=1
453 mothup(2,4)=2
454 icolup(1,4)=0
455 icolup(2,4)=0
456 vtimup(4)=0
457 spinup(4)=9
458
459ccc HEPEVT
460 if(output.eq.'hepevt')then
461
462 nhep=nup
463 write(6,*) "hepevt", nhep, nup
464 do k=1,5
465 phep(k,1)=pup(k,1)
466 phep(k,2)=pup(k,2)
467 enddo
468
469 do k=1,nup
470 isthep(k)=istup(k)
471 idhep(k)=idup(k)
472 jmohep(1,k)=mothup(1,k)
473 jmohep(2,k)=jmohep(1,k)
474 jdahep(1,k)=0
475 jdahep(2,k)=0
476 vhep(1,k)=0d0
477 vhep(2,k)=0d0
478 vhep(3,k)=0d0
479 vhep(4,k)=0d0
480 enddo
481
482 if(pflag.eq.'rho')then
483 jdahep(1,5)=7
484 jdahep(2,5)=8
485 jdahep(1,6)=9
486 jdahep(2,6)=10
487 endif
488
489 jmohep(2,5)=2
490 jmohep(2,6)=2
491
492 endif
493
494 if(unw.eq.'true')then
495 lmax=2
496 else
497 lmax=1
498 endif
499
500c do ll=1,lmax
501
502 if(ll.eq.2)then
503c ntotal=nev*10
504 ntotal=1000000000
505 endif
506
507 ip=ntotal+1
508
509ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
510c c
511c Start of event loop c
512c c
513ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
514
515 if(ll.eq.1)then
516 write(6,*)'Generating weighted events...'
517 else
518 write(6,*)'Generating unweighted events...'
519 endif
520 return
521 end
522
523 subroutine dimegenerate
524 implicit none
525
526ccccc hepevt output
527 integer nmxhep,kk
528 parameter (nmxhep=4000)
529 integer nevhep,nhep,isthep,idhep,jmohep,jdahep
530 double precision phep,vhep
531 common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
532 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
533ccccc Les Houches Event Common Block
534 INTEGER MAXNUP
535 PARAMETER (MAXNUP=500)
536 INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
537 DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
538 COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
539 & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
540 & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
541 & SPINUP(MAXNUP)
542ccccc rambo variables
543 integer npart,nrun
544 double precision pin(4), am(100), pout(4,100), wt, ein
545 common/ramboc/ pin, am, pout, wt, ein,
546 &npart, nrun
547
548ccccc
549 double precision etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,
550 &rmin,mcut
551 common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
552ccccc
553 double precision q(4,20)
554 common /mom/ q
555ccccc
556 double precision pt1(2), pt2(2), ptx(2)
557 common/mompt/pt1,pt2,ptx
558ccccc
559 double precision s, rts, mmes, yx
560 double precision bjac, bp
561 double precision mf127, mf1525
562
563 double precision m0, mmes0, mmes1, mmes2, mp, mwidth, pi, rt2,
564 &ebeam, sh, th, uh, sum, sum1, weightm
565
566 integer id(20), ncut, ll, icut, nev, num, iin
567
568 common/vars/s,rts,mmes,yx
569 common/hvars/sh,th,uh
570 common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp,
571 & mwidth, pi, rt2, ebeam, sum, sum1, weightm,
572 & bjac, bp
573 common /ivars/ ncut, ll, icut, nev, num
574ccccc
575 character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
576 &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
577 common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw
578 common/ff/formf
579
580ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
581cccccc local variables ...
582 double precision aa1, aa2, al1, al2, almax, almin, c, cc1, cc2,
583 & msmax, msmin, cont, exn, mmesp, mwidth1, mwidth2, norm,
584 & p1m, p1p, p2m, p2p, phi1, phi2, phix1, pt1sq, pt1x, pt1y,
585 & pt2sq, pt2x, pt2y, ptxsq1, ptxsq2, ptxsqmax, ptxsqmin,
586 & ran0, ran1, ran2, ran3, ran4, ran5, ran6, ranhist, ranx1,
587 & ranx2, ranx3, ranx4, rm1, rm2, rmx1, rmx2,
588 & root1sq, root2sq, snp, t1, t2, weight, weight0, weight1,
589 & weight2, wthist, x1, x2, ymax1, ymax2, ymin1, ymin2, yx1,
590 & yx2
591 double precision pboo(4), pcm(4), plb(4)
592 double precision svec(4)
593 complex*16 zmat
594 integer h, i, j, k, ntotal
595 weight=0d0
596 call r2455(ran0)
597 call r2455(ran1)
598 call r2455(ran2)
599 call r2455(ran3)
600 call r2455(ran4)
601 call r2455(ran5)
602 call r2455(ran6)
603
604ccccccccccccccccccccccccccccccccccccccccccccccccccc
605
606c incoming protons
607
608 ID(1)=2212
609 q(1,1)=0d0
610 q(2,1)=0d0
611 q(3,1)=ebeam
612 q(4,1)=ebeam
613
614 ID(2)=2212
615 q(1,2)=0d0
616 q(2,2)=0d0
617 q(3,2)=-ebeam
618 q(4,2)=ebeam
619
620 phi1=2d0*pi*ran0
621 phi2=2d0*pi*ran1
622
623 pt1sq=-dlog(ran2)/(bjac)
624 pt2sq=-dlog(ran3)/(bjac)
625
626 weight=(dexp(bjac*pt1sq)*dexp(bjac*pt2sq))/bjac**2
627
628 pt1(1)=dsqrt(pt1sq)*dsin(phi1)
629 pt1(2)=dsqrt(pt1sq)*dcos(phi1)
630 pt2(1)=dsqrt(pt2sq)*dsin(phi2)
631 pt2(2)=dsqrt(pt2sq)*dcos(phi2)
632
633 pt1x=pt1(1)
634 pt1y=pt1(2)
635 pt2x=pt2(1)
636 pt2y=pt2(2)
637
638cccccccccccccccccccccccccccccccc
639
640ccc non-zero rho width
641
642cccccccccccccccccccccccccccccccc
643
644 if(pflag.eq.'rho')then ! new part here
645
646 677 call r2455(rm1)
647 call r2455(rm2)
648
649 msmax=mmes0+4d0*mwidth
650 msmin=2d0*0.13957018d0
651
652cccccccccccccccc
653
654 almin=datan(-(mmes0**2-msmin**2)/mwidth/mmes0)
655 almax=datan(-(mmes0**2-msmax**2)/mwidth/mmes0)
656
657 al1=almin+(almax-almin)*rm1
658 al2=almin+(almax-almin)*rm2
659
660 mmes1=dsqrt(dtan(al1)*mmes0*mwidth+mmes0**2)
661 mmes2=dsqrt(dtan(al2)*mmes0*mwidth+mmes0**2)
662
663 weight=weight*(almax-almin)**2
664 weight=weight*mwidth**2*mmes0**2
665 weight=weight*(1d0+dtan(al1)**2)*(1d0+dtan(al2)**2)
666 weight=weight/4d0/mmes1/mmes2
667
668ccccccccccccccccc
669
670 mwidth1=mwidth*((1d0-4d0*0.13957018d0**2/mmes1**2)/
671 & (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0
672 mwidth2=mwidth*((1d0-4d0*0.13957018d0**2/mmes2**2)/
673 & (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0
674
675 if(mmes1.lt.2d0*0.13957018d0) goto 677
676 if(mmes2.lt.2d0*0.13957018d0) goto 677
677
678 weight=weight*2d0*mmes0*mmes1*mwidth1/pi
679 weight=weight*2d0*mmes0*mmes2*mwidth2/pi
680 weight=weight*mmes1**2*mmes2**2/mmes0**4
681 weight=weight/((mmes0**2-mmes1**2)**2+mmes1**4*
682 & mwidth1**2/mmes0**2)
683 weight=weight/((mmes0**2-mmes2**2)**2+mmes2**4*
684 & mwidth2**2/mmes0**2)
685
686 endif
687
688
689cccccccccccccccccccccccccccccccccccccccccccccccccccc
690
691c pi rapidities
692
693 call r2455(ranx1)
694 call r2455(ranx2)
695 call r2455(ranx3)
696 call r2455(ranx4)
697
698 phix1=2d0*pi*ranx1
699
700c r2max=1d0
701c r2=r2max*ranx1
702c ptxsq1=-dlog(r2)/bjac1 ! generate exp p_t^2 (can be more efficient)
703
704 ptxsqmin=0d0
705 ptxsqmax=10d0 ! generate linear p_t^2
706
707 ptxsq1=ptxsqmin+(ptxsqmax-ptxsqmin)*ranx2
708
709 q(1,5)=dsqrt(ptxsq1)*dcos(phix1)
710 q(2,5)=dsqrt(ptxsq1)*dsin(phix1)
711 q(1,6)=-pt1(1)-pt2(1)-q(1,5)
712 q(2,6)=-pt1(2)-pt2(2)-q(2,5)
713 ptxsq2=q(1,6)**2+q(2,6)**2
714
715 rmx1=dsqrt(mmes1**2+ptxsq1)
716 rmx2=dsqrt(mmes2**2+ptxsq2)
717
718 ymax1=dlog(rts/rmx1)
719
720 if(cuts.eq.'true')then
721 if(rmax.lt.ymax1)then
722 ymax1=rmax
723 endif
724 endif
725
726 ymin1=-ymax1
727
728 yx1=ymin1+(ymax1-ymin1)*ranx3
729
730 ymax2=(dlog((rts-rmx1*dexp(yx1))/rmx2))
731 ymin2=-(dlog((rts-rmx1*dexp(-yx1))/rmx2))
732
733 if(cuts.eq.'true')then
734 if(ymax2.gt.rmax)then
735 ymax2=rmax
736 endif
737 if(ymin2.lt.-rmax)then
738 ymin2=-rmax
739 endif
740 endif
741
742 if(ymax2.lt.ymin2)then
743 weight=0d0
744 goto 700
745 endif
746
747 yx2=ymin2+(ymax2-ymin2)*ranx4
748 x1=(rmx2*dexp(yx2)+rmx1*dexp(yx1))/rts
749 x2=(rmx2*dexp(-yx2)+rmx1*dexp(-yx1))/rts
750
751 weight=weight*(ptxsqmax-ptxsqmin)
752c weight=weight*dexp(bjac1*ptxsq1)*r2max/bjac1
753 weight=weight*(ymax1-ymin1)*(ymax2-ymin2)
754
755
756 q(3,5)=rmx1*(dexp(yx1)-dexp(-yx1))/2d0
757 q(4,5)=rmx1*(dexp(yx1)+dexp(-yx1))/2d0
758
759 q(3,6)=rmx2*(dexp(yx2)-dexp(-yx2))/2d0
760 q(4,6)=rmx2*(dexp(yx2)+dexp(-yx2))/2d0
761
762ccccccccccccccccccccccccccccccccccccccccccccccccccccc
763
764c impose massive on-shell condition by solving
765c p1+ + cc1/p2- = aa1
766c p2- + cc2/p1+ = aa2
767
768 aa1=bp*(1d0-x1)
769 aa2=bp*(1d0-x2)
770 cc1=0.5d0*(pt2sq+mp**2)
771 cc2=0.5d0*(pt1sq+mp**2)
772
773 root1sq=(cc1-cc2-aa1*aa2)**2-4d0*cc2*aa1*aa2
774 root2sq=(cc2-cc1-aa1*aa2)**2-4d0*cc1*aa1*aa2
775 if(root1sq.le.0d0.or.root2sq.le.0d0)then
776 weight=0d0
777 goto 700
778 endif
779 p1p=(cc2-cc1+aa1*aa2+dsqrt(root1sq))/(2d0*aa2)
780 p2m=(cc1-cc2+aa1*aa2+dsqrt(root2sq))/(2d0*aa1)
781 p1m=(pt1sq+mp**2)/(2d0*p1p)
782 p2p=(pt2sq+mp**2)/(2d0*p2m)
783
784 if(p1p.lt.0d0.or.p1m.lt.0d0.or.p2p.lt.0d0.or.p2m.lt.0d0)then
785 weight=0d0
786 goto 700
787 endif
788
789 t1=-rts*p1m*rt2
790 t2=-rts*p2p*rt2
791
792 q(1,3)=pt1(1)
793 q(2,3)=pt1(2)
794 q(3,3)=(p1p-p1m)/rt2
795 q(4,3)=(p1p+p1m)/rt2
796
797 q(1,4)=pt2(1)
798 q(2,4)=pt2(2)
799 q(3,4)=(p2p-p2m)/rt2
800 q(4,4)=(p2p+p2m)/rt2
801
802ccccccccccccccccccccccccccccccccccccccccccccccccc
803
804 do k=1,4
805 svec(k)=q(k,5)+q(k,6)
806 enddo
807
808 call dot(svec,svec,sh)
809
810ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
811
812c rho0 decays
813
814ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
815
816 if(pflag.eq.'rho')then
817
818 kk=6
819
820 do k=5,6
821
822 if(k.eq.5)then
823 mmesp=mmes1
824 else
825 mmesp=mmes2
826 endif
827 call rambo(npart,mmesp,am,pout,wt) ! call RAMBO
828
829 do j=1,4
830 pboo(j)=q(j,k)
831 enddo
832
833 do h=1,2
834 do j=1,4
835 pcm(j)=pout(j,h)
836 enddo
837 call boost(mmesp,pboo,pcm,plb)
838 kk=kk+1
839 do j=1,4
840 q(j,kk)=plb(j)
841 enddo
842 enddo
843
844 enddo
845
846 endif
847
848
849ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
850
851c place cuts
852
853 if(cuts.eq.'true')then
854
855 call cut(icut)
856
857 if(icut.eq.0)then
858 weight=0d0
859 weight0=0d0
860 weight1=0d0
861 weight2=0d0
862 ncut=ncut+1
863 goto 700
864 endif
865
866 endif
867
868c print*,i
869
870cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
871
872c Write 4-momenta to event record array
873
874ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
875
876 if(output.eq.'lhe')then
877
878 do k=3,4
879 do j=1,4
880 pup(j,k)=q(j,k)
881 enddo
882 pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
883 enddo
884
885 if(pflag.eq.'rho')then
886
887 do k=5,10
888 do j=1,4
889 pup(j,k)=q(j,k)
890 enddo
891 pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
892 enddo
893
894 else
895
896 do k=5,6
897 do j=1,4
898 pup(j,k)=q(j,k)
899 enddo
900 pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
901 enddo
902
903 endif
904
905 else
906
907 do k=3,4
908 do j=1,4
909 phep(j,k)=q(j,k)
910 enddo
911 phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
912 enddo
913
914 if(pflag.eq.'rho')then
915
916 do k=5,10
917 do j=1,4
918 phep(j,k)=q(j,k)
919 enddo
920 phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
921 enddo
922
923 else
924
925 do k=5,6
926 do j=1,4
927 phep(j,k)=q(j,k)
928 enddo
929 phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
930 enddo
931
932 endif
933
934 endif
935
936ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
937
938c Weight evaluation
939
940ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
941
942c pp-->pp(M_1 M_2) matrix element
943
944ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
945
946
947 call schimc(pt1x,pt1y,pt2x,pt2y,zmat)
948
949 weight=weight*cdabs(zmat)**2
950 weight=weight/norm**2
951
952ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
953
954 weight=weight/(s**2*(16d0)**3*pi**5)*389d3
955
956ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
957
958ccc Probability of no additional particles produced in production subprocess
959
960cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
961
962 if(fsi.eq.'true')then
963
964 if(pflag.eq.'pipm'.or.pflag.eq.'pi0')then
965
966 if(dsqrt(sh).gt.mf127)then
967 m0=mf127
968 c=0.7d0
969 cont=1d0-(m0-2d0*mmes)**2/m0**2
970 exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
971 snp=dexp(-exn)
972 else
973 snp=1d0
974 endif
975
976 elseif(pflag.eq.'kpkm'.or.pflag.eq.'ks')then
977
978 if(dsqrt(sh).gt.mf1525)then
979 m0=mf1525
980 c=0.7d0
981 cont=1d0-(m0-2d0*mmes)**2/m0**2
982 exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
983 snp=dexp(-exn)
984 else
985 snp=1d0
986 endif
987
988 else
989
990 snp=1d0
991
992 endif
993
994 else
995
996 snp=1d0
997
998 endif
999
1000 if(snp.gt.1d0)then
1001 snp=1d0
1002 endif
1003
1004 weight=weight*snp
1005
1006ccccccccccccccc
1007
1008 if(pflag.eq.'pi0'.or.pflag.eq.'ks'.or.pflag.eq.'rho')then
1009 weight=weight/2d0 ! symmetry factor
1010 endif
1011
1012cccccccccccccccccccccccccccccccccccccccccccccccccccccc
1013
1014 666 weight=weight
1015
1016 if(weightm.lt.weight)then
1017 weightm=weight
1018 endif
1019
1020
1021 if(ll.eq.1)then
1022
1023 xwgtup=weight
1024
1025 wthist=weight/dfloat(ntotal)
1026
1027 call binit(wthist)
1028
1029 else
1030
1031 call r2455(ranhist)
1032
1033 if(ranhist.lt.weight/weightm)then
1034
1035 xwgtup=1d0
1036
1037 num=num+1
1038
1039
1040c call binit(sumt/nev)
1041
1042 write(35,302)num,nup
1043
1044 if(output.eq.'lhe')then
1045
1046 do j=1,nup
1047 write(35,300)j,istup(j),idup(j),mothup(1,j),mothup(2,j),
1048 & icolup(1,j),icolup(2,j),pup(1,j),pup(2,j),
1049 & pup(3,j),pup(4,j),pup(5,j),vtimup(j),spinup(j)
1050 enddo
1051
1052 else
1053
1054 do j=1,nup
1055 write(35,301)j,idhep(j),isthep(j),jmohep(1,j),
1056 & jmohep(2,j),jdahep(1,j),jdahep(2,j),
1057 & phep(1,j),phep(2,j),phep(3,j),phep(4,j),phep(5,j)
1058 & ,vhep(1,j),vhep(2,j),vhep(3,j),vhep(4,j)
1059 enddo
1060
1061 endif
1062
1063
1064 endif
1065
1066 endif
1067
1068 if(ll.eq.2)then
1069
1070 if(num.gt.nev-1)then ! exit loop once required no. of unweighted events generated
1071
1072 ntotal=i
1073! goto 888
1074 call terminate()
1075
1076 endif
1077 endif
1078
1079 700 sum=sum+weight
1080 sum1=sum1+weight**2
1081
1082cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1083c c
1084c End of event loop c
1085c c
1086cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1087 call terminate()
1088 300 format(i4,1x,i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,
1089 &E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6)
1090 301 format(i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,E13.6
1091 &,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x
1092 &,E13.6)
1093 302 format(1x,i8,1x,i5,1x,E13.6)
1094
1095 return
1096 end
1097
1098 subroutine terminate
1099 character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
1100 &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
1101
1102 888 continue
1103 sum=sum/dfloat(ntotal)
1104 sum1=sum1/dfloat(ntotal)
1105 var=dsqrt((sum1-sum**2)/dfloat(ntotal))
1106 EFF=1d0*(ntotal-ncut)/ntotal
1107
1108 if(ll.eq.1)then
1109 print*,'...done!'
1110 write(6,*)
1111 sumt=sum
1112 else
1113 print*,'...done!'
1114 write(6,*)
1115 endif
1116
1117
1118
1119ccccc create histograms
1120
1121 if(ll.eq.1)then
1122
1123 do i=1,nhist
1124 call histo2(i,0)
1125 enddo
1126
1127 endif
1128
1129 if(ll.eq.1)then
1130 if(pflag.eq.'pipm')then
1131 write(6,*)
1132 write(6,*)'pi+pi- production'
1133 elseif(pflag.eq.'pi0')then
1134 write(6,*)
1135 write(6,*)'pi0pi0 production'
1136 elseif(pflag.eq.'kpkm')then
1137 write(6,*)
1138 write(6,*)'K+K- production'
1139 elseif(pflag.eq.'ks')then
1140 write(6,*)
1141 write(6,*)'K0K0 production'
1142 elseif(pflag.eq.'rho')then
1143 write(6,*)
1144 write(6,*)'rho0rho0 production'
1145 endif
1146 endif
1147
1148 if(ll.eq.1)then
1149 if(pflag.eq.'rho')then
1150 write(6,221)rts,ntotal,sum,var
1151 else
1152 write(6,222)rts,mmes,ntotal,sum,var
1153 endif
1154 else
1155 if(output.eq.'lhe')then
1156 print*,'LHE output'
1157 else
1158 print*,'HEPEVT output'
1159 endif
1160 write(6,223)nev,sum,var
1161 endif
1162
1163! enddo
1164
1165 close(35)
1166
1167 221 format(/,
1168 . 3x,'collider energy (GeV) : ',f10.3,/,
1169 . 3x,'number of events : ',i9,/,
1170 . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5,
1171 . /)
1172
1173 222 format(/,
1174 . 3x,'collider energy (GeV) : ',f10.3,/,
1175 . 3x,'meson mass (GeV) : ',f10.5,/,
1176 . 3x,'number of events : ',i9,/,
1177 . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5,
1178 . /)
1179
1180 223 format(3x,'number of events : ',i9,/,
1181 . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5,
1182 . /)
1183
1184 call cpu_time(t2)
1185 print*,'time elapsed = ', t2, ' s'
1186
1187 return
1188 end
1189
1190ccccc Pomeron -- (off-shell) meson form factor
1191
1192 function fpi(x)
1193 implicit double precision(a-y)
1194 character formf*10
1195 common/vars/s,rts,mmes,yx
1196 common/ff/formf
1197 common/ffpars/bexp,ao,bo,aoo,boo
1198
1199 if(formf.eq.'exp')then
1200 fpi=exp((x-mmes**2)*bexp)
1201 elseif(formf.eq.'orexp')then
1202 fpi=exp(-(dsqrt(-x+mmes**2+ao**2)-ao)*bo)
1203 elseif(formf.eq.'power')then
1204 fpi=1d0/(1d0-(x-mmes**2)/aoo)
1205 endif
1206
1207 return
1208 end
1209
1210ccccccc Pom Pom --> meson pair amplitude
1211
1212 subroutine wev(matt,v1,v2)
1213 implicit double precision(a-y)
1214 implicit complex*16(z)
1215 double precision q(4,20),svec(4),tvec(4),uvec(4),v1(4),v2(4)
1216 integer k
1217 character*10 mregge,pflag
1218 complex*16 matt
1219 common/mom/q
1220 common/vars/s,rts,mmes,yx
1221 common/wvars/sig0,bb,t1,t2
1222 common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1223 common/sec/cpom,cf,crho,aff,ar
1224 common/hvars/sh,th,uh
1225 common/regge/mregge
1226 common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw
1227
1228
1229 zi=(0d0,1d0)
1230 pi=dacos(-1d0)
1231
1232 do k=1,4
1233 q(k,10)=v1(k)
1234 q(k,11)=v2(k)
1235 enddo
1236
1237 do k=1,4
1238 tvec(k)=q(k,10)-q(k,5)
1239 uvec(k)=q(k,10)-q(k,6)
1240 enddo
1241
1242 call dot(tvec,tvec,th)
1243 call dot(uvec,uvec,uh)
1244
1245 sht1=2d0*(q(4,5)*q(4,3)-q(3,5)*q(3,3)-q(2,5)*q(2,3)-
1246 &q(1,5)*q(1,3))+mmes**2
1247
1248 sht2=2d0*(q(4,6)*q(4,4)-q(3,6)*q(3,4)-q(2,6)*q(2,4)-
1249 &q(1,6)*q(1,4))+mmes**2
1250
1251 shu1=2d0*(q(4,6)*q(4,3)-q(3,6)*q(3,3)-q(2,6)*q(2,3)-
1252 &q(1,6)*q(1,3))+mmes**2
1253
1254 shu2=2d0*(q(4,5)*q(4,4)-q(3,5)*q(3,4)-q(2,5)*q(2,4)-
1255 &q(1,5)*q(1,4))+mmes**2
1256
1257 bt1=(2d0*alphap*dlog(sht1))
1258 bt2=(2d0*alphap*dlog(sht2))
1259 bu1=(2d0*alphap*dlog(shu1))
1260 bu2=(2d0*alphap*dlog(shu2))
1261
1262 bt1r=(2d0*alphapr*dlog(sht1))
1263 bt2r=(2d0*alphapr*dlog(sht2))
1264 bu1r=(2d0*alphapr*dlog(shu1))
1265 bu2r=(2d0*alphapr*dlog(shu2))
1266
1267 sig0t1=sig0*sht1**0.0808d0/0.389d0
1268 sig0t2=sig0*sht2**0.0808d0/0.389d0
1269 sig0u1=sig0*shu1**0.0808d0/0.389d0
1270 sig0u2=sig0*shu2**0.0808d0/0.389d0
1271
1272 t11=v1(1)**2+v1(2)**2
1273 t22=v2(1)**2+v2(2)**2
1274
1275 zmt=(zi*dexp(-bt1*t11/2d0)*cpom*sht1**alpha0
1276 & +((aff+zi)*cf*sht1**alpha0r+(ar-zi)*crho*sht1**alpha0r)
1277 & *dexp(-bt1r*t11/2d0))
1278 zmt=zmt*(zi*dexp(-bt2*t22/2d0)*cpom*sht2**alpha0
1279 & +((aff+zi)*cf*sht2**alpha0r-(ar-zi)*crho*sht2**alpha0r)
1280 & *dexp(-bt2r*t11/2d0))
1281 zmt=zmt*fpi(th)**2/(mmes**2-th)
1282
1283 zmu=(zi*dexp(-bu1*t11/2d0)*cpom*shu1**alpha0
1284 & +((aff+zi)*cf*shu1**alpha0r-(ar-zi)*crho*shu1**alpha0r)
1285 & *dexp(-bu1r*t11/2d0))
1286 zmu=zmu*(zi*dexp(-bu2*t22/2d0)*cpom*shu2**alpha0
1287 & +((aff+zi)*cf*shu2**alpha0r+(ar-zi)*crho*shu2**alpha0r)
1288 & *dexp(-bu2r*t22/2d0))
1289 zmu=zmu*fpi(uh)**2/(mmes**2-uh)
1290
1291 matt=zmu+zmt
1292
1293 return
1294 end
1295
1296c binning subroutine
1297
1298 subroutine binit(wt)
1299 implicit double precision(a-y)
1300 double precision q(4,20),pt1(2),pt2(2),ptx(2)
1301 common/mom/q
1302 common/mompt/pt1,pt2,ptx
1303 common/vars/s,rts,mmes,yx
1304 common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2
1305 common/vars2/ptpi1,etapi1,ptpi2,etapi2
1306 common/ang/cost,costa,costb
1307 common/hvars/sh,th,uh
1308
1309 ypisum=0.5d0*dlog((q(4,5)+q(4,6)+q(3,5)+q(3,6))
1310 & /(q(4,5)+q(4,6)-q(3,5)-q(3,6)))
1311
1312 ypi1=0.5d0*dlog((q(4,5)+q(3,5))
1313 & /(q(4,5)-q(3,5)))
1314 ypi2=0.5d0*dlog((q(4,6)+q(3,6))
1315 & /(q(4,6)-q(3,6)))
1316
1317 call histo1(1,30,0d0,4.5d0,dsqrt(sh),wt)
1318
1319 return
1320 end
1321
1322ccccc Cutting subroutine
1323ccccc
1324ccccc User can define their own cuts on any particle momenta, stored in array q(i,j)
1325ccccc
1326
1327 subroutine cut(icut)
1328 implicit double precision(a-y)
1329 double precision q(4,20)
1330 integer icut
1331 common/mom/q
1332 common/vars/s,rts,mmes,yx
1333 common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2
1334 common/vars2/ptpi1,etapi1,ptpi2,etapi2
1335 common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
1336 common/hvars/sh,th,uh
1337 integer iii
1338 common/ivar/iii
1339
1340 icut=0
1341
1342c -- meson 1 variables
1343 ptpi1=dsqrt(q(1,5)**2+q(2,5)**2)
1344 pmodpi1=dsqrt(q(1,5)**2+q(2,5)**2+q(3,5)**2)
1345 etapi1=.5d0*dlog((pmodpi1+q(3,5))
1346 . /(pmodpi1-q(3,5)))
1347
1348 ypi1=.5d0*dlog((q(4,5)+q(3,5))
1349 . /(q(4,5)-q(3,5)))
1350
1351 etpi1=q(4,5)*dsqrt(pmodpi1**2-q(3,5)**2)/pmodpi1
1352
1353c -- meson 2 variables
1354 ptpi2=dsqrt(q(1,6)**2+q(2,6)**2)
1355 pmodpi2=dsqrt(q(1,6)**2+q(2,6)**2+q(3,6)**2)
1356 etapi2=.5d0*dlog((pmodpi2+q(3,6))
1357 . /(pmodpi2-q(3,6)))
1358
1359 ypi2=.5d0*dlog((q(4,6)+q(3,6))
1360 . /(q(4,6)-q(3,6)))
1361
1362 etpi2=q(4,6)*dsqrt(pmodpi2**2-q(3,6)**2)/pmodpi2
1363
1364c print*,ptpi1,ptpi2
1365
1366ccccc example cuts....
1367
1368c if(dsqrt(sh).lt.0.8d0)return
1369
1370 xf1=2d0*q(3,3)/rts
1371 xf2=2d0*q(3,4)/rts
1372
1373c if(dabs(xf1).lt.0.9d0)return
1374c if(dabs(xf2).lt.0.9d0)return
1375
1376ccccccccc Basic cuts described at start of code
1377
1378 if(ptpi1.lt.ecut)return
1379 if(ptpi2.lt.ecut)return
1380 if(ypi1.gt.rmax)return
1381 if(ypi2.gt.rmax)return
1382 if(ypi1.lt.rmin)return
1383 if(ypi2.lt.rmin)return
1384
1385 icut=1
1386
1387 return
1388 end
1389
1390c prints histograms
1391
1392 subroutine histo1(ih,ib,x0,x1,x,w)
1393 implicit real*8(a-h,o-z)
1394 character*1 regel(30),blank,star
1395 dimension h(20,100),hx(20),io(20),iu(20),ii(20)
1396 dimension y0(20),y1(20),ic(20)
1397 data regel / 30*' ' /,blank /' ' /,star /'*'/
1398 save
1399 open(10,file="output.dat")
1400 y0(ih)=x0
1401 y1(ih)=x1
1402 ic(ih)=ib
1403 if(x.lt.x0) goto 11
1404 if(x.gt.x1) goto 12
1405 ix=idint((x-x0)/(x1-x0)*dble(ib))+1
1406 h(ih,ix)=h(ih,ix)+w
1407 if(h(ih,ix).gt.hx(ih)) hx(ih)=h(ih,ix)
1408 ii(ih)=ii(ih)+1
1409 return
1410 11 iu(ih)=iu(ih)+1
1411 return
1412 12 io(ih)=io(ih)+1
1413 return
1414 entry histo2(ih,il)
1415 ib1=ic(ih)
1416 x01=y0(ih)
1417 x11=y1(ih)
1418 bsize=(x11-x01)/dble(ib1)
1419 hx(ih)=hx(ih)*(1.d0+1.d-06)
1420 if(il.eq.0) write(6,21) ih,ii(ih),iu(ih),io(ih)
1421 if(il.eq.1) write(6,22) ih,ii(ih),iu(ih),io(ih)
1422 21 format(' no.',i3,' lin : inside,under,over ',3i6)
1423 22 format(' no.',i3,' log : inside,under,over ',3i6)
1424 if(ii(ih).eq.0) goto 28
1425 write(6,23)
1426 23 format(35(1h ),3(10h----+----i))
1427 do 27 iv=1,ib1
1428 z=(dble(iv)-0.5d0)/dble(ib1)*(x11-x01)+x01
1429 if(il.eq.1) goto 24
1430 iz=idint(h(ih,iv)/hx(ih)*30.)+1
1431 goto 25
1432 24 iz=-1
1433 if(h(ih,iv).gt.0.d0)
1434 .iz=idint(dlog(h(ih,iv))/dlog(hx(ih))*30.)+1
1435 25 if(iz.gt.0.and.iz.le.30) regel(iz)=star
1436 write(6,26) z,h(ih,iv)/bsize,(regel(i),i=1,30)
1437c write(10,*)z-0.125d0/2d0,h(ih,iv)/bsize ! Print histogram to file
1438 write(10,*)z,h(ih,iv)/bsize ! Print histogram to file
1439c write(10,*)z,h(ih,iv)/bsize ! Print histogram to file
1440 26 format(1h ,2g15.6,4h i,30a1,1hi)
1441 36 format(1h ,2g15.6)
1442 if(iz.gt.0.and.iz.le.30) regel(iz)=blank
1443 27 continue
1444 write(6,23)
1445 return
1446 28 write(6,29)
1447 29 format(' no entries inside histogram')
1448 return
1449 entry histo3(ih)
1450 do 31 i=1,100
1451 31 h(ih,i)=0.
1452 hx(ih)=0.
1453 io(ih)=0
1454 iu(ih)=0
1455 ii(ih)=0
1456 close(10)
1457 return
1458 end
1459
1460
1461cccc Initializes soft model parameters
1462
1463 subroutine initpars(in)
1464 implicit double precision(a-y)
1465 integer in,i1,i2
1466 integer nch
1467 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1468 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1469 common/ipars/nch
1470 common/vars/s,rts,mmes,yx
1471
1472 pi=dacos(-1d0)
1473
1474 if(in.eq.1)then
1475
1476 ep=0.13d0
1477 asp=0.08d0
1478 ep1=0d0
1479 sigo=23d0/0.39d0
1480 gaa(3)=0.4d0
1481 gd2=0.3d0
1482 nch=2
1483 ntf=0
1484
1485 cc0(1)=0.45d0
1486 bm(1)=3d0
1487 bb0(1)=0.1d0
1488 pp0(1)=0.92d0
1489 bex(1)=8.5d0
1490
1491 cc0(2)=0.45d0
1492 bm(2)=1.5d0
1493 bb0(2)=0.5d0
1494 pp0(2)=0.1d0
1495 bex(2)=4.5d0
1496
1497 cc0(3)=1d0
1498 bm(3)=0d0
1499 bb0(3)=0.8d0
1500 pp0(3)=0.5d0
1501 bex(3)=0.5d0
1502
1503 elseif(in.eq.2)then
1504
1505 ep=0.115d0
1506 asp=0.11d0
1507 ep1=0d0
1508 sigo=33d0/0.39d0
1509 gaa(3)=0.6d0
1510 gd2=0.16d0
1511 nch=2d0
1512 ntf=0d0
1513
1514 cc0(1)=0.63d0
1515 bm(1)=3d0
1516 bb0(1)=0.1d0
1517 pp0(1)=0.5d0
1518 bex(1)=8d0
1519
1520 cc0(2)=0.47d0
1521 bm(2)=1.5d0
1522 bb0(2)=0.5d0
1523 pp0(2)=0.1d0
1524 bex(2)=6d0
1525
1526 cc0(3)=1d0
1527 bm(3)=0d0
1528 bb0(3)=0.8d0
1529 pp0(3)=0.5d0
1530 bex(3)=0.5d0
1531
1532 elseif(in.eq.3)then
1533
1534 ep=0.093d0
1535 asp=0.075d0
1536 ep1=0d0
1537 sigo=60d0/0.39d0
1538 gaa3=4.8d0
1539 gd2=1.03d0
1540 nch=2
1541 nga=1
1542 cc0(1)=0.55d0
1543 bm(1)=3d0
1544 bb0(1)=0.27d0
1545 pp0(1)=0.48d0
1546 bex(1)=5.3d0
1547
1548 cc0(2)=0.48d0
1549 bm(2)=1.5d0
1550 bb0(2)=0.1d0
1551 pp0(2)=1d0
1552 bex(2)=3.8d0
1553
1554 cc0(3)=0.24d0
1555
1556 elseif(in.eq.4)then
1557
1558 ep=0.11d0 !!! Capital delta
1559 asp=0.06d0 !!! alpha'
1560 ep1=0d0 !!! Zero in all models, matters (?)
1561 sigo=50d0/0.39d0 !!! sigma_0 (GeV^-2)
1562 gaa3=6d0 !!! k2/k(1.8 TeV)
1563 gd2=1.3d0 !!! k1/k(1.8 TeV)
1564 nch=2
1565 nga=1 !!! A flag (?)
1566 cc0(1)=0.6d0 !!! d1
1567 bm(1)=3d0 !!! Doesn't matter
1568 bb0(1)=0.45d0 !!! c1-0.08 (added back later)
1569 pp0(1)=0.5d0 !!! 2*|a_1|^2
1570 bex(1)=7.2d0 !!! b1
1571
1572 cc0(2)=0.48d0 !!! d2
1573 bm(2)=1.5d0 !!! Doesn't matter
1574 bb0(2)=0.16d0 !!! c2-0.08 (added back later)
1575 pp0(2)=1d0 !!! |a_2|^2 is set later
1576 bex(2)=4.2d0 !!! b2
1577
1578 cc0(3)=0.12d0 !!! Beta : k^2_min ~ s^Beta
1579
1580 endif
1581
1582 if(nch.eq.3) pp0(3)=3d0-pp0(2)-pp0(1)
1583 if(nch.eq.2) pp0(2)=2d0-pp0(1) !!! Set |a_2|^2
1584 if(nch.eq.1) pp0(1)=1d0
1585
1586 if(in.eq.3.or.in.eq.4)then
1587
1588 gamm=(1800d0/rts)**cc0(3)
1589 ga1=1d0/(1d0+gamm*gd2)
1590 ga2=1d0/(1d0+gamm*gaa3)
1591 gaa(1)=2d0*ga1/(ga1+ga2) !!! gamma_1 (+) Eq (15)
1592 gaa(2)=2d0*ga2/(ga1+ga2) !!! gamma_2 (-) Eq (15)
1593
1594 elseif(in.eq.1.or.in.eq.2)then
1595
1596 if(nch.eq.2)then
1597 gaa(1)=1d0+dsqrt(gd2)
1598 gaa(2)=1d0-dsqrt(gd2)
1599 gaa(3)=0.
1600 elseif(nch.eq.1)then
1601 gaa(1)=1d0
1602 gaa(2)=0d0
1603 gaa(3)=0d0
1604 endif
1605
1606 endif
1607
1608 sum=0d0
1609
1610 do i1=1,nch
1611 do i2=1,nch
1612 sum=sum+gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1613 enddo
1614 enddo
1615
1616 norm=sum
1617
1618 return
1619 end
1620
1621ccccc Calculates proton opacity
1622
1623 subroutine calcop
1624 implicit double precision(a-y)
1625 integer nb,i1,i2,nch,ib
1626 common/ipars/nch
1627 double precision op(5,5,10000,2),oph(5,5,10000,2)
1628 common/opac/op,oph
1629
1630 nb=900
1631 hb=100d0/dble(nb)
1632
1633 print*,'Calculating opacity'
1634
1635
1636 do ib=1,nb+1
1637 bt=dble(ib-1)*hb
1638
1639 do i1=1,nch
1640 do i2=1,nch
1641
1642 call opacity(i1,i2,bt,fr,fr1)
1643
1644
1645 op(i1,i2,ib,1)=bt
1646 op(i1,i2,ib,2)=fr
1647 oph(i1,i2,ib,1)=bt
1648 oph(i1,i2,ib,2)=fr1
1649
1650 write(40,*)bt,fr,fr1
1651
1652 enddo
1653 enddo
1654 enddo
1655
1656 return
1657 end
1658
1659cccc Calculates screening amplitude (in k_t space)
1660
1661 subroutine calcscreen
1662 implicit double precision(a-y)
1663 integer ns,i1,i2,nch,ib
1664 common/ipars/nch
1665 double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1666 common/spac/sca,sca1
1667
1668 ns=900
1669 ksqma=8.2d0
1670 inck=ksqma/dble(ns)
1671 ksqmin=0.001d0
1672 lginck=(dlog(ksqma/ksqmin))/dble(ns)
1673
1674 print*,'Calculating screening amplitude'
1675
1676 do ib=0,ns+1
1677
1678 ksq=dble(ib-1)*inck
1679 lgksq=ksq
1680
1681 if(ib.eq.0)then
1682 ksq=0d0
1683 lgksq=0d0
1684 else
1685 lgksq=dble(ib-1)*lginck+dlog(ksqmin)
1686 ksq=dexp(lgksq)
1687 endif
1688
1689 do i1=1,nch
1690 do i2=1,nch
1691
1692 call screening(i1,i2,ksq,sc,sc1)
1693
1694 sca(i1,i2,ib,1)=lgksq
1695 sca(i1,i2,ib,2)=sc
1696 sca1(i1,i2,ib,1)=lgksq
1697 sca1(i1,i2,ib,2)=sc1
1698
1699 enddo
1700 enddo
1701
1702 enddo
1703
1704 return
1705 end
1706
1707ccccc Interpolation....
1708
1709 subroutine opacityint(i,j,bt,out,out1)
1710 implicit double precision(a-y)
1711 double precision op(5,5,10000,2),oph(5,5,10000,2)
1712 integer i,j,it
1713 common/opac/op,oph
1714
1715 incbt=op(1,1,2,1)-op(1,1,1,1)
1716 it=nint(bt/incbt)
1717 if(dble(it).gt.(bt/incbt))then
1718 it=it-1
1719 endif
1720
1721 m=(op(i,j,it+2,2)-op(i,j,it+1,2))
1722 &/(op(i,j,it+2,1)-op(i,j,it+1,1))
1723 del=bt-op(1,1,it+1,1)
1724 mh=(oph(i,j,it+2,2)-oph(i,j,it+1,2))
1725 &/(oph(i,j,it+2,1)-oph(i,j,it+1,1))
1726 delh=bt-oph(1,1,it+1,1)
1727
1728 out=m*del+op(i,j,it+1,2)
1729 out1=mh*delh+oph(i,j,it+1,2)
1730
1731 return
1732 end
1733
1734 subroutine screeningint(i,j,ktsq,out,out1)
1735 implicit double precision(a-y)
1736 double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1737 integer i,j,it,ns
1738 common/spac/sca,sca1
1739
1740 if(ktsq.lt.0.001d0)then
1741
1742 it=0
1743
1744 m=(sca(i,j,it+1,2)-sca(i,j,it,2))
1745 & /(dexp(sca(i,j,it+1,1))-sca(i,j,it,1))
1746 del=ktsq-sca(1,1,it,1)
1747 m1=sca1(i,j,it+1,2)-sca1(i,j,it,2)
1748 & /(dexp(sca1(i,j,it+1,1))-sca1(i,j,it,1))
1749 del1=ktsq-sca1(1,1,it,1)
1750
1751 out=m*del+sca(i,j,it,2)
1752 out1=m1*del1+sca1(i,j,it,2)
1753
1754 elseif(ktsq.lt.8d0)then
1755
1756 ktmin=sca(1,1,1,1)
1757 inckt=sca(1,1,2,1)-sca(1,1,1,1)
1758 it=nint((dlog(ktsq)-ktmin)/inckt)
1759 if(dble(it).gt.((dlog(ktsq)-ktmin)/inckt))then
1760 it=it-1
1761 endif
1762
1763 m=(sca(i,j,it+2,2)-sca(i,j,it+1,2))
1764 & /(sca(i,j,it+2,1)-sca(i,j,it+1,1))
1765 del=dlog(ktsq)-sca(1,1,it+1,1)
1766 m1=sca1(i,j,it+2,2)-sca1(i,j,it+1,2)
1767 & /(sca1(i,j,it+2,1)-sca1(i,j,it+1,1))
1768 del1=dlog(ktsq)-sca1(1,1,it+1,1)
1769
1770 out=m*del+sca(i,j,it+1,2)
1771 out1=m1*del1+sca1(i,j,it+1,2)
1772
1773 else
1774
1775 out=0d0
1776 out1=0d0
1777
1778 endif
1779
1780 return
1781 end
1782
1783cccc Integrates round Pomeron loop (to calculate screened amplitude)
1784
1785 subroutine schimc(p1x,p1y,p2x,p2y,out)
1786 implicit double precision(a-y)
1787 integer nx,jx,jy,nch,i1,i2,it,k
1788 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1789 double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1790 double precision a1(4),a2(4),q(4,20)
1791 complex*16 out,x0,x00
1792 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1793 common/ipars/nch
1794 common/spac/sca,sca1
1795 common/mom/q
1796 common/wvars/sig0,bb,t1,t2
1797 common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1798
1799 do k=1,4
1800 a1(k)=q(k,1)-q(k,3)
1801 a2(k)=q(k,2)-q(k,4)
1802 enddo
1803
1804 nx=10
1805 hx=2d0/dble(nx)
1806
1807 out=0d0
1808
1809 do jx=-nx,nx
1810 do jy=-nx,nx
1811
1812 tpx=dble(jx)*hx
1813 tpy=dble(jy)*hx
1814 tp2=tpx**2+tpy**2
1815 wt=hx*hx
1816
1817 p1xp=p1x-tpx
1818 p1yp=p1y-tpy
1819 t12=p1xp**2+p1yp**2
1820 p2xp=tpx+p2x
1821 p2yp=tpy+p2y
1822 t22=p2xp**2+p2yp**2
1823
1824 a1(1)=-p1xp
1825 a1(2)=-p1yp
1826 a2(1)=-p2xp
1827 a2(2)=-p2yp
1828
1829 do i1=1,nch
1830 do i2=1,nch
1831
1832 call screeningint(i1,i2,tp2,sc,sc1)
1833 call wev(x0,a1,a2)
1834
1835 x0=x0*pp0(i1)*pp0(i2)/dble(nch)**2
1836 x0=x0*dexp(-((t12+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+
1837 & (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1))
1838 x0=x0*dexp(-((t22+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+
1839 & (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2))
1840
1841 out=out+x0*wt*sc
1842
1843 enddo
1844 enddo
1845
1846
1847 enddo
1848 enddo
1849
1850 a1(1)=-p1x
1851 a1(2)=-p1y
1852 a2(1)=-p2x
1853 a2(2)=-p2y
1854
1855 t11=p1x**2+p1y**2
1856 t22=p2x**2+p2y**2
1857
1858 call formfac(t11,t22,x00p)
1859
1860 call wev(x00,a1,a2)
1861
1862
1863
1864 out=out+x00*x00p
1865
1866 return
1867 end
1868
1869ccccc Pomeron -- diffrative eignstate i form factor
1870
1871 subroutine formfac(t1,t2,out)
1872 implicit double precision(a-y)
1873 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1874 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1875 integer nch,i1,i2
1876 common/ipars/nch
1877
1878 out=0d0
1879
1880 delta1=dexp(-t1/2d0)
1881 delta2=dexp(-t2/2d0)
1882
1883 do i1=1,nch
1884 do i2=1,nch
1885
1886 wt=gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1887 wt=wt*dexp(-((t1+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+
1888 & (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1))
1889 wt=wt*dexp(-((t2+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+
1890 & (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2))
1891
1892 out=out+wt
1893
1894 enddo
1895 enddo
1896
1897 return
1898 end
1899
1900 subroutine screening(i,j,ktsq,out,out1)
1901 implicit double precision(a-y)
1902 integer i,j,ib,nb,it,nt,nch
1903 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1904 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1905 common/ipars/nch
1906 common/vars/s,rts,mmes,yx
1907
1908 pi=dacos(-1d0)
1909
1910 nb=5000
1911 hb=99d0/dble(nb)
1912
1913 out=0d0
1914 out1=0d0
1915
1916 do ib=1,nb
1917 bt=dble(ib)*hb
1918 wt=-bt/2d0/pi*hb
1919
1920 call opacityint(i,j,bt,fr,fr1)
1921
1922 sige=sigo*dexp(dlog(rts)*2d0*ep)
1923 fr=fr*gaa(i)*gaa(j)*sige
1924 fr1=fr1*gaa(i)*gaa(j)*sige
1925
1926 out=out+wt*(1d0-dexp(-fr/2d0))*besj0(bt*dsqrt(ktsq))
1927 & *gaa(i)*gaa(j)
1928 out1=out1+wt*(1d0-dexp(-fr1/2d0))*besj0(bt*dsqrt(ktsq))
1929 & *gaa(i)*gaa(j)
1930
1931 enddo
1932
1933 return
1934 end
1935
1936
1937 subroutine opacity(i,j,bt,out,out1)
1938 implicit double precision(a-y)
1939 integer i,j,ib,nb,it,nt,nch
1940 double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1941 common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1942 common/ipars/nch
1943 common/vars/s,rts,mmes,yx
1944
1945 pi=dacos(-1d0)
1946
1947 ampi=0.02d0
1948 amro=1d0
1949 a4=4d0*ampi
1950 alo=dlog(amro/ampi)
1951 bpol=2.4d0
1952
1953 nt=6000
1954 htt=6d0/dble(nt)**2
1955
1956 out=0d0
1957 out1=0d0
1958
1959 do it=0,nt
1960 t=dble(it)**2*htt
1961 if(it.eq.0) t=1d-8
1962 wt=htt*2d0*dble(it)/4d0/pi
1963 if(it.eq.0) wt=wt/2d0
1964 bes0=besj0(bt*dsqrt(t))
1965
1966 ffi=dexp(-((t+0.08d0+bb0(i))*bex(i))**cc0(i)+
1967 & (bex(i)*(bb0(i)+0.08d0))**cc0(i))
1968 ffj=dexp(-((t+0.08d0+bb0(j))*bex(j))**cc0(j)+
1969 & (bex(j)*(bb0(j)+0.08d0))**cc0(j))
1970
1971 asp1=asp
1972 form1=dlog(ffi*ffj)-2d0*t*asp1*dlog(rts)
1973cccccccccccccccc
1974 ah=dsqrt(1d0+a4/t)
1975 h1pi=2d0*a4+t*(alo-ah*ah*ah*dlog((1d0+ah)/(ah-1d0))) !!! Pion loop insertion
1976 h1pi=h1pi*sigo/(72d0*pi**3*(1d0+t/bpol)**2)
1977ccccccccccccccc
1978 ww=bes0*dexp(form1-2d0*h1pi*dlog(rts))
1979 aspt=t*asp+h1pi
1980
1981 out=out+ww*wt
1982 out1=out1+bes0*wt*ffi*ffj
1983
1984 enddo
1985
1986
1987 return
1988 end
1989
1990 function rgauss(r1,r2)
1991 implicit double precision (a-y)
1992
1993 pi=dacos(-1d0)
1994 rgauss=dsqrt(-2d0*dlog(r1))*dsin(2d0*pi*r2)
1995
1996 return
1997 end
1998
1999c boosting subroutine
2000
2001 subroutine boost(p,pboo,pcm,plb)
2002 real*8 pboo(4),pcm(4),plb(4),p,fact
2003 plb(4)=(pboo(4)*pcm(4)+pboo(3)*pcm(3)
2004 & +pboo(2)*pcm(2)+pboo(1)*pcm(1))/p
2005 fact=(plb(4)+pcm(4))/(p+pboo(4))
2006 do 10 j=1,3
2007 10 plb(j)=pcm(j)+fact*pboo(j)
2008 return
2009 end
2010
2011c calculates lorentzian dot product for real 4-vectors
2012
2013 subroutine dot(v1,v2,dt)
2014 double precision v1(4),v2(4),dt
2015
2016 dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
2017 &-v1(4)*v2(4))
2018
2019 return
2020 end
2021c
2022
2023*
2024* subtractive mitchell-moore generator
2025* ronald kleiss - october 2, 1987
2026*
2027* the algorithm is N(i)=[ N(i-24) - N(i-55) ]mod M,
2028* implemented in a cirucular array with identifcation
2029* of NR(i+55) and nr(i), such that effectively:
2030* N(1) <--- N(32) - N(1)
2031* N(2) <--- N(33) - N(2) ....
2032* .... N(24) <--- N(55) - N(24)
2033* N(25) <--- N(1) - N(25) ....
2034* .... N(54) <--- N(30) - N(54)
2035* N(55) <--- N(31) - N(55)
2036*
2037* in this version M =2**30 and RM=1/M=2.D0**(-30.D0)
2038*
2039* the array NR has been initialized by putting NR(i)=i
2040* and subsequently running the algorithm 100,000 times.
2041*
2042
2043 subroutine R2455(Ran)
2044 IMPLICIT REAL*8(A-H,O-Z)
2045 DIMENSION N(55)
2046 DATA N/
2047 . 980629335, 889272121, 422278310,1042669295, 531256381,
2048 . 335028099, 47160432, 788808135, 660624592, 793263632,
2049 . 998900570, 470796980, 327436767, 287473989, 119515078,
2050 . 575143087, 922274831, 21914605, 923291707, 753782759,
2051 . 254480986, 816423843, 931542684, 993691006, 343157264,
2052 . 272972469, 733687879, 468941742, 444207473, 896089285,
2053 . 629371118, 892845902, 163581912, 861580190, 85601059,
2054 . 899226806, 438711780, 921057966, 794646776, 417139730,
2055 . 343610085, 737162282,1024718389, 65196680, 954338580,
2056 . 642649958, 240238978, 722544540, 281483031,1024570269,
2057 . 602730138, 915220349, 651571385, 405259519, 145115737/
2058 DATA M/1073741824/
2059 DATA RM/0.9313225746154785D-09/
2060 DATA K/55/,L/31/
2061 IF(K.EQ.55) THEN
2062 K=1
2063 ELSE
2064 K=K+1
2065 ENDIF
2066 IF(L.EQ.55) THEN
2067 L=1
2068 ELSE
2069 L=L+1
2070 ENDIF
2071 J=N(L)-N(K)
2072 IF(J.LT.0) J=J+M
2073 N(K)=J
2074 RAN=J*RM
2075 END
2076
2077
2078
2079 SUBROUTINE RAMBO(N,ET,XM,P,WT)
2080C------------------------------------------------------
2081C
2082C RAMBO
2083C
2084C RA(NDOM) M(OMENTA) B(EAUTIFULLY) O(RGANIZED)
2085C
2086C A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
2087C AUTHORS: S.D. ELLIS, R. KLEISS, W.J. STIRLING
2088C THIS IS VERSION 1.0 - WRITTEN BY R. KLEISS
2089C
2090C N = NUMBER OF PARTICLES (>1, IN THIS VERSION <101)
2091C ET = TOTAL CENTRE-OF-MASS ENERGY
2092C XM = PARTICLE MASSES ( DIM=100 )
2093C P = PARTICLE MOMENTA ( DIM=(4,100) )
2094C WT = WEIGHT OF THE EVENT
2095C
2096C------------------------------------------------------
2097 implicit none
2098! boost arguments
2099 integer n
2100 double precision et,xm(100),p(4,100),wt
2101! function
2102 double precision rn
2103! internal variables
2104 double precision q(4,100),z(100),r(4),
2105 & b(3),p2(100),xm2(100),e(100),v(100)
2106 integer iwarn(5)
2107 double precision acc
2108 integer itmax,ibegin
2109 data acc/1.0d-14/,itmax/6/,ibegin/0/,iwarn/5*0/
2110
2111 integer i,k,iter
2112 double precision twopi,po2log,xmt,c,s,f,rmas,g,a,x,bq,
2113 & accu,f0,g0,wt2,wt3,wtm,x2,xmax
2114 integer nm
2115 save
2116
2117
2118C
2119C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
2120 IF(IBEGIN.NE.0) GOTO 103
2121 IBEGIN=1
2122 TWOPI=8.*DATAN(1.D0)
2123 PO2LOG=DLOG(TWOPI/4.)
2124 Z(2)=PO2LOG
2125 DO 101 K=3,100
2126 101 Z(K)=Z(K-1)+PO2LOG-2.*DLOG(DFLOAT(K-2))
2127 DO 102 K=3,100
2128 102 Z(K)=(Z(K)-DLOG(DFLOAT(K-1)))
2129C
2130C CHECK ON THE NUMBER OF PARTICLES
2131 103 IF(N.GT.1.AND.N.LT.101) GOTO 104
2132 PRINT 1001,N
2133 STOP
2134C
2135C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
2136 104 XMT=0.
2137 NM=0
2138 DO 105 I=1,N
2139 IF(XM(I).NE.0.D0) NM=NM+1
2140 105 XMT=XMT+DABS(XM(I))
2141 IF(XMT.LE.ET) GOTO 201
2142 PRINT 1002,XMT,ET
2143 STOP
2144C
2145C THE PARAMETER VALUES ARE NOW ACCEPTED
2146C
2147C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
2148 201 DO 202 I=1,N
2149 C=2.*RN(1)-1.
2150 S=DSQRT(1.-C*C)
2151 F=TWOPI*RN(2)
2152 Q(4,I)=-DLOG(RN(3)*RN(4))
2153 Q(3,I)=Q(4,I)*C
2154 Q(2,I)=Q(4,I)*S*DCOS(F)
2155 202 Q(1,I)=Q(4,I)*S*DSIN(F)
2156
2157C
2158C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
2159 DO 203 I=1,4
2160 203 R(I)=0.
2161 DO 204 I=1,N
2162 DO 204 K=1,4
2163 204 R(K)=R(K)+Q(K,I)
2164 RMAS=DSQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
2165 DO 205 K=1,3
2166 205 B(K)=-R(K)/RMAS
2167 G=R(4)/RMAS
2168 A=1./(1.+G)
2169 X=ET/RMAS
2170C
2171C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
2172 DO 207 I=1,N
2173 BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
2174 DO 206 K=1,3
2175 206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
2176 207 P(4,I)=X*(G*Q(4,I)+BQ)
2177C
2178C CALCULATE WEIGHT AND POSSIBLE WARNINGS
2179 WT=PO2LOG
2180 IF(N.NE.2) WT=(2.*N-4.)*DLOG(ET)+Z(N)
2181 IF(WT.GE.-180.D0) GOTO 208
2182 IF(IWARN(1).LE.5) PRINT 1004,WT
2183 IWARN(1)=IWARN(1)+1
2184 208 IF(WT.LE. 174.D0) GOTO 209
2185 IF(IWARN(2).LE.5) PRINT 1005,WT
2186 IWARN(2)=IWARN(2)+1
2187C
2188C RETURN FOR WEIGHTED MASSLESS MOMENTA
2189 209 IF(NM.NE.0) GOTO 210
2190 WT=DEXP(WT)
2191 RETURN
2192C
2193C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
2194 210 XMAX=DSQRT(1.-(XMT/ET)**2)
2195 DO 301 I=1,N
2196 XM2(I)=XM(I)**2
2197 301 P2(I)=P(4,I)**2
2198 ITER=0
2199 X=XMAX
2200 ACCU=ET*ACC
2201 302 F0=-ET
2202 G0=0.
2203 X2=X*X
2204 DO 303 I=1,N
2205 E(I)=DSQRT(XM2(I)+X2*P2(I))
2206 F0=F0+E(I)
2207 303 G0=G0+P2(I)/E(I)
2208 IF(DABS(F0).LE.ACCU) GOTO 305
2209 ITER=ITER+1
2210 IF(ITER.LE.ITMAX) GOTO 304
2211C PRINT 1006,ITMAX,ACCU,DABS(F0)
2212 WRITE(99,1006)ITMAX,ACCU,DABS(F0)
2213 IF (DABS(F0).GT.1.0E-6) THEN
2214 WRITE(*,1007)DABS(F0)
2215 ENDIF
2216 GOTO 305
2217 304 X=X-F0/(X*G0)
2218 GOTO 302
2219 305 DO 307 I=1,N
2220 V(I)=X*P(4,I)
2221 DO 306 K=1,3
2222 306 P(K,I)=X*P(K,I)
2223 307 P(4,I)=E(I)
2224C
2225C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
2226 WT2=1.
2227 WT3=0.
2228 DO 308 I=1,N
2229 WT2=WT2*V(I)/E(I)
2230 308 WT3=WT3+V(I)**2/E(I)
2231 WTM=(2.*N-3.)*DLOG(X)+DLOG(WT2/WT3*ET)
2232C
2233C RETURN FOR WEIGHTED MASSIVE MOMENTA
2234 WT=WT+WTM
2235 IF(WT.GE.-180.D0) GOTO 309
2236 IF(IWARN(3).LE.5) PRINT 1004,WT
2237 IWARN(3)=IWARN(3)+1
2238 309 IF(WT.LE. 174.D0) GOTO 310
2239 IF(IWARN(4).LE.5) PRINT 1005,WT
2240 IWARN(4)=IWARN(4)+1
2241 310 WT=DEXP(WT)
2242 RETURN
2243C
2244 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
2245 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
2246 . ' SMALLER THAN TOTAL ENERGY =',D15.6)
2247 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
2248 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY OVERFLOW')
2249 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
2250 . ' DESIRED ACCURACY =',D10.3,' (STOPS AT',D10.3,')')
2251 1007 FORMAT(' RAMBO WARNS: END OF ITERATION ACCURACY TOO LOW =',D10.3)
2252 END
2253C
2254
2255 function rn(idum)
2256*
2257* SUBTRACTIVE MITCHELL-MOORE GENERATOR
2258* RONALD KLEISS - OCTOBER 2, 1987
2259*
2260* THE ALGORITHM IS N(I)=[ N(I-24) - N(I-55) ]MOD M,
2261* IMPLEMENTED IN A CIRUCULAR ARRAY WITH IDENTIFCATION
2262* OF NR(I+55) AND NR(I), SUCH THAT EFFECTIVELY:
2263* N(1) <--- N(32) - N(1)
2264* N(2) <--- N(33) - N(2) ....
2265* .... N(24) <--- N(55) - N(24)
2266* N(25) <--- N(1) - N(25) ....
2267* .... N(54) <--- N(30) - N(54)
2268* N(55) <--- N(31) - N(55)
2269*
2270* IN THIS VERSION M =2**30 AND RM=1/M=2.D0**(-30D0)
2271*
2272* THE ARRAY NR HAS BEEN INITIALIZED BY PUTTING NR(I)=I
2273* AND SUBSEQUENTLY RUNNING THE ALGORITHM 100,000 TIMES.
2274*
2275
2276 implicit none
2277 double precision rn
2278 integer idum
2279 integer n(55)
2280 data n/
2281 . 980629335, 889272121, 422278310,1042669295, 531256381,
2282 . 335028099, 47160432, 788808135, 660624592, 793263632,
2283 . 998900570, 470796980, 327436767, 287473989, 119515078,
2284 . 575143087, 922274831, 21914605, 923291707, 753782759,
2285 . 254480986, 816423843, 931542684, 993691006, 343157264,
2286 . 272972469, 733687879, 468941742, 444207473, 896089285,
2287 . 629371118, 892845902, 163581912, 861580190, 85601059,
2288 . 899226806, 438711780, 921057966, 794646776, 417139730,
2289 . 343610085, 737162282,1024718389, 65196680, 954338580,
2290 . 642649958, 240238978, 722544540, 281483031,1024570269,
2291 . 602730138, 915220349, 651571385, 405259519, 145115737/
2292 double precision eps
2293 double precision rm
2294 integer j,k,l,m
2295
2296 data eps/1D-9/
2297 data m/1073741824/
2298 data rm/0.9313225746154785D-09/
2299 data k/55/,l/31/
2300
2301
2302 1 CONTINUE
2303 IF(K.EQ.55) THEN
2304 K=1
2305 ELSE
2306 K=K+1
2307 ENDIF
2308 IF(L.EQ.55) THEN
2309 L=1
2310 ELSE
2311 L=L+1
2312 ENDIF
2313 J=N(L)-N(K)
2314 IF(J.LT.0) J=J+M
2315 N(K)=J
2316 RN=J*RM
2317 IF(RN.LT.EPS) GOTO 1
2318 IF(RN.GT.1D0-EPS) GOTO 1
2319 RETURN
2320 END