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