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