5
INCLUDE "ngraphs.inc" !how many diagrams
10
include 'nexternal.inc'
11
include 'nexternal_prod.inc'
18
REAL*8 P(0:3,NEXTERNAL_PROD)
19
REAL*8 PFULL(0:3,NEXTERNAL), Ptrial(0:3,NEXTERNAL)
20
double precision x(36), Ecollider
21
CHARACTER*120 BUFF(NEXTERNAL_PROD)
22
integer iforest(2,-nexternal:-1,N_MAX_CG)
23
integer itree(2,-nexternal:-1), iconfig
24
integer map_external2res(nexternal_prod) ! map (index in production) -> index in the full structure
25
c integer mapconfig(0:lmaxconfigs)
26
integer sprop(1,-nexternal:-1,N_MAX_CG) ! first col has one entry, since we should have group_processes=false
27
integer tprid(-nexternal:-1,N_MAX_CG)
28
integer mapconfig(0:N_MAX_CG), this_config
29
common/to_mconfigs/mapconfig, this_config
30
double precision prmass(-nexternal:0,N_MAX_CG)
31
double precision prwidth(-nexternal:0,N_MAX_CG)
32
integer pow(-nexternal:0,N_MAX_CG)
33
double precision qmass(-nexternal:0),qwidth(-nexternal:0),jac
34
double precision M_PROD, M_FULL
36
integer counter,mode,nbpoints, counter2
37
double precision mean, variance, maxweight,weight,std
39
double precision Pprod(0:3,nexternal_prod)
40
integer nb_mc_masses, indices_mc_masses(nexternal)
41
double precision values_mc_masses(nexternal)
43
! variables to keep track of the vegas numbers for the production part
44
logical keep_inv(-nexternal:-1),no_gen
46
double precision fixedinv(-nexternal:0)
47
double precision phi_tchan(-nexternal:0),m2_tchan(-nexternal:0)
48
double precision cosphi_schan(-nexternal:0), phi_schan(-nexternal:0)
49
common /to_fixed_kin/keep_inv,no_gen, ivar, fixedinv,
50
& phi_tchan,m2_tchan,cosphi_schan, phi_schan
52
double precision BWcut, maxBW
53
common /to_BWcut/BWcut, maxBW
55
c Conflicting BW stuff
56
integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
57
double precision cBW_mass(-nexternal:-1,-1:1),
58
& cBW_width(-nexternal:-1,-1:1)
59
common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
63
DOUBLE PRECISION AMP2(n_max_cg)
66
! variables associate with the PS generation
67
double precision totmassin, totmass
68
double precision shat, sqrtshat, stot, y, m(-nexternal:nexternal)
69
integer nbranch, ns_channel,nt_channel, pos_pz
71
& totmassin, totmass,shat, sqrtshat, stot,y, m,
72
& nbranch, ns_channel,nt_channel, pos_pz
74
integer*8 iseed, P_seed
78
real*8 ranu(97),ranc,rancd,rancm
80
common/ raset1 / ranu,ranc,rancd,rancm
81
common/ raset2 / iranmr,jranmr
83
c call ntuple(x,0d0,1d0,1,2) ! initialize the sequence of random
84
! numbers at the position reached
85
! at the previous termination of the
88
open(unit=56,file='seeds.dat',status='old')
91
open(unit=56,file='offset.dat',status='old')
94
iseed = iseed + P_seed
96
cccccccccccccccccccccccccccccccccccccccccccccccccccc
97
c I. read momenta for the production events
99
cccccccccccccccccccccccccccccccccccccccccccccccccccc
101
! hard-code for testing
102
c buff(1)="1 0.5000000E+03 0.0000000E+00 0.0000000E+00 0.5000000E+03 0.0000000E+00"
103
c buff(2)="2 0.5000000E+03 0.0000000E+00 0.0000000E+00 -0.5000000E+03 0.0000000E+00"
104
c buff(3)="3 0.5000000E+03 0.1040730E+03 0.4173556E+03 -0.1872274E+03 0.1730000E+03"
105
c buff(4)="4 0.5000000E+03 -0.1040730E+03 -0.4173556E+03 0.1872274E+03 0.1730000E+03"
106
c do i=1,nexternal_prod
107
c read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
113
read(*,*) mode, BWcut, Ecollider, temp
116
if (mode.eq.1) then ! calculate the maximum weight
118
elseif (mode.eq.2) then
119
maxweight=temp ! unweight decay config
120
elseif (mode.eq.3) then
121
continue ! just retrun the value of M_full
123
write(*,*) ranu,ranc,rancd,rancm,iranmr,jranmr
125
goto 2 ! and close the program
128
do i=1,nexternal_prod
129
read (*,*) P(0,i),P(1,i),P(2,i),P(3,i)
133
read(*,*) nb_mc_masses
134
if (nb_mc_masses.gt.0) then
135
read (*,*) (indices_mc_masses(k), k=1,nb_mc_masses)
136
read (*,*) (values_mc_masses(k), k=1,nb_mc_masses)
139
c write(*,*) sqrt(dot(P(0,3),P(0,3)))
140
c write(*,*) sqrt(dot(P(0,4),P(0,4)))
141
c write(*,*) dot(P(0,5),P(0,5))
142
c write(*,*) 'shat', sqrt(2d0*dot(P(0,1),P(0,2)))
143
c write(*,*) 'pt2g', sqrt(2d0*dot(P(0,4),P(0,5))+dot(P(0,4),P(0,4)))
144
cccccccccccccccccccccccccccccccccccccccccccccccccccc
145
c II. initialization of masses and widths c
146
c also load production topology information c
147
cccccccccccccccccccccccccccccccccccccccccccccccccccc
149
include 'configs_production.inc'
152
call set_parameters(p,Ecollider)
154
include 'props_production.inc'
157
! do not consider the case of conflicting BW
158
do i = -nexternal, -1
163
cccccccccccccccccccccccccccccccccccccccccccccccccccc
164
c III. compute production matrix element c
165
cccccccccccccccccccccccccccccccccccccccccccccccccccc
171
CALL SMATRIX_PROD(P,M_PROD)
172
c write(*,*) 'M_prod ', M_prod
173
cccccccccccccccccccccccccccccccccccccccccccccccccccc
174
c IV. select one topology c
175
cccccccccccccccccccccccccccccccccccccccccccccccccccc
177
call get_config(iconfig)
178
do i=-nexternal_prod+2,-1
180
itree(j,i)=iforest(j,i,iconfig)
184
do i=-nexternal_prod+3,-1
185
qmass(i)=prmass(i,iconfig)
186
qwidth(i)=prwidth(i,iconfig)
189
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
190
c V. load topology for the whole event select c
191
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
193
call merge_itree(itree,qmass,qwidth, p,map_external2res)
194
!write(*,*) keep_inv(-5)
195
!write(*,*) 'm2_tchan ',m2_tchan(-5)
196
!write(*,*) 'fixedinv', fixedinv(-5)
197
!write(*,*) 'phi_tchan', phi_tchan(-5)
199
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
200
c VIa. Calculate the max. weight c
201
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
217
do j = 1, 3*(nexternal-nexternal_prod)
218
call ntuple(x(j),0d0,1d0,j,1)
221
c do j=1,nexternal_prod
222
c write (*,*) (p(k,j), k=0,3)
225
call generate_momenta_conf(jac,x,itree,qmass,qwidth,pfull,pprod,map_external2res)
228
if (counter2.ge.3) then ! use another topology to generate PS points
229
call get_config(iconfig)
230
do k=-nexternal_prod+2,-1
232
itree(j,k)=iforest(j,k,iconfig)
236
do k=-nexternal_prod+3,-1
237
qmass(k)=prmass(k,iconfig)
238
qwidth(k)=prwidth(k,iconfig)
240
call merge_itree(itree,qmass,qwidth, p,map_external2res)
246
! write (*,*) (pfull(k,j), k=0,3)
248
call SMATRIX(pfull,M_full)
249
call SMATRIX_PROD(pprod,M_prod)
250
c write(*,*) 'M_full ', M_full
251
c write(*,*) 'jac',jac
253
weight=M_full*jac/M_prod
254
if (weight.gt.maxweight) then
260
c max_mom(j,k)=pfull(j,k)
265
c variance=variance+weight**2
267
c mean=mean/real(nbpoints)
268
c variance=variance/real(nbpoints)-mean**2
270
write (*,*) maxweight ! ,mean,std
271
c write (*,*) 'max_m',max_m
272
c write (*,*) 'max_jac', jac
273
c write (*,*) 'Extrenal masses'
275
c write(*,*) dot(max_mom(0,k), max_mom(0,k))
278
c pw1(j)=max_mom(j,4)+max_mom(j,5)
279
c pt1(j)=pw1(j)+max_mom(j,3)
280
c pw2(j)=max_mom(j,7)+max_mom(j,8)
281
c pt2(j)=pw2(j)+max_mom(j,6)
282
c pt2g(j)=pt2(j)+max_mom(j,9)
285
c write (*,*) 'm45', sqrt(2D0*dot(max_mom(0,4),max_mom(0,5)))
286
c write (*,*) 'm78', sqrt(2d0*dot(max_mom(0,7),max_mom(0,8)))
287
c write (*,*) 'mt1', sqrt(dot(pt1,pt1))
288
c write (*,*) 'mt2', sqrt(dot(pt2,pt2))
289
c write (*,*) 'mt2g', sqrt(dot(pt2g,pt2g))
290
c write (*,*) 'm9', sqrt(dot(max_mom(0,9),max_mom(0,9)))
291
c write (*,*) 'shat', sqrt(2D0*dot(max_mom(0,2),max_mom(0,1)))
292
c write(*,*) (max_mom(j,1), j=0,3)
293
c write(*,*) (max_mom(j,2), j=0,3)
294
c write(*,*) (pt1(j), j=0,3)
295
c write(*,*) (pt2(j), j=0,3)
296
c write(*,*) (max_mom(j,9), j=0,3)
302
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
303
c VIb. generate unweighted decay config c
304
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
316
do i = 1, 3*(nexternal-nexternal_prod)+1
317
call ntuple(x(i),0d0,1d0,i,1)
320
call generate_momenta_conf(jac,x,itree,qmass,qwidth,pfull,pprod,map_external2res)
323
if (counter2.ge.3) then ! use another topology to generate PS points
324
call get_config(iconfig)
325
do i=-nexternal_prod+2,-1
327
itree(j,i)=iforest(j,i,iconfig)
331
do i=-nexternal_prod+3,-1
332
qmass(i)=prmass(i,iconfig)
333
qwidth(i)=prwidth(i,iconfig)
335
call merge_itree(itree,qmass,qwidth, p,map_external2res)
340
call SMATRIX(pfull,M_full)
341
call SMATRIX_PROD(pprod,M_prod)
343
weight=M_full*jac/M_prod
345
if (weight.gt.x(3*(nexternal-nexternal_prod)+1)*maxweight) notpass=.false.
350
if (nb_mc_masses.gt.0) then
354
m(indices_mc_masses(k))=values_mc_masses(k)
356
call generate_momenta_conf(jac,x,itree,qmass,qwidth,ptrial,pprod,map_external2res)
359
write(*,*) nexternal, counter, maxBW, weight, counter2, 0
361
write (*,*) (pfull(j,i), j=0,3)
364
write(*,*) nexternal, counter, maxBW, weight, counter2, 1
366
write (*,*) (ptrial(j,i), j=0,3)
370
write(*,*) nexternal, counter, maxBW, weight, counter2, 0
372
write (*,*) (pfull(j,i), j=0,3)
380
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
381
c VIc. return M_full^2 c
382
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
393
do i = 1, 3*(nexternal-nexternal_prod)+1
394
call ntuple(x(i),0d0,1d0,i,1)
397
call generate_momenta_conf(jac,x,itree,qmass,qwidth,pfull,pprod,map_external2res)
398
if (jac.lt.0d0) cycle
400
call SMATRIX(pfull,M_full)
411
subroutine get_config(iconfig)
414
C--- integer n_max_cg
415
INCLUDE "ngraphs.inc" !how many diagrams
422
double precision cumulweight(0:n_max_cg),random
425
DOUBLE PRECISION AMP2(n_max_cg)
428
integer mapconfig(0:N_MAX_CG), this_config
429
common/to_mconfigs/mapconfig, this_config
431
if (mapconfig(0).eq.1) then
438
cumulweight(i)=amp2(mapconfig(i))+cumulweight(i-1)
442
call ntuple(random,0d0,1d0,24,1)
446
cumulweight(i)=cumulweight(i)/cumulweight(mapconfig(0))
448
!write(*,*) cumulweight(i-1)
449
!write(*,*) cumulweight(i)
450
if (random.ge.cumulweight(i-1).and.random.le.cumulweight(i)) then
456
write(*,*) 'Unable to generate iconfig ', random, mapconfig(0),
462
subroutine set_parameters(p,Ecollider)
465
double precision ZERO
467
include 'nexternal_prod.inc'
469
double precision p(0:3, nexternal_prod), Ecollider
473
double precision ptot(0:3)
475
include 'nexternal.inc'
479
! variables associate with the PS generation
480
double precision totmassin, totmass
481
double precision shat, sqrtshat, stot, y, m(-nexternal:nexternal)
482
integer nbranch, ns_channel,nt_channel, pos_pz
484
& totmassin, totmass,shat, sqrtshat, stot,y, m,
485
& nbranch, ns_channel,nt_channel, pos_pz
487
c Masses of particles. Should be filled in setcuts.f
488
double precision pmass(nexternal)
489
common /to_mass/pmass
494
include "../parameters.inc"
501
! m(i) pmass(i) already refer to the full kinematics
506
if (p(3,1).gt.p(3,2)) then
512
ptot(j)=p(j,1)+p(j,2)
516
y = 0.5*log((ptot(0)+ptot(3))/(ptot(0)-ptot(3)))
518
c write(*,*) sqrtshat
520
c write(*,*) (ptot(j), j=0,3)
523
c Make sure have enough mass for external particles
526
totmassin=totmassin+pmass(i)
533
if (stot .lt. max(totmass,totmassin)**2) then
534
write (*,*) 'Fatal error #0 in one_tree:'/
535
& /'insufficient collider energy'
542
subroutine merge_itree(itree,qmass,qwidth, p_ext,mapext2res)
546
include 'nexternal_prod.inc'
547
include 'nexternal.inc'
551
integer itree(2,-nexternal:-1) ! PS structure for the production
552
double precision qmass(-nexternal:0),qwidth(-nexternal:0)
553
double precision p_ext(0:3,nexternal_prod)
554
integer mapext2res(nexternal_prod)
558
double precision normp
559
! info for the full process
560
integer itree_full(2,-nexternal:-1),i,j
561
double precision qmass_full(-nexternal:0),qwidth_full(-nexternal:0)
562
! info for the decay part
563
double precision idecay(2,-nexternal:-1), prmass(-nexternal:-1),prwidth(-nexternal:-1)
564
integer ns_channel_decay
565
integer map_external2res(nexternal_prod) ! map (index in production) -> index in the full structure
566
double precision p(0:3,-nexternal:nexternal)
568
integer idB, id1, index_p2
569
double precision pa(0:3), pb(0:3), p1(0:3), p2(0:3),pboost(0:3)
570
double precision pb_cms(0:3), p1_cms(0:3), p1_rot(0:3)
573
! variables to keep track of the vegas numbers for the production part
574
logical keep_inv(-nexternal:-1),no_gen
576
double precision fixedinv(-nexternal:0)
577
double precision phi_tchan(-nexternal:0),m2_tchan(-nexternal:0)
578
double precision cosphi_schan(-nexternal:0), phi_schan(-nexternal:0)
579
common /to_fixed_kin/keep_inv,no_gen, ivar, fixedinv,
580
& phi_tchan,m2_tchan,cosphi_schan, phi_schan
582
! variables associate with the PS generation
583
double precision totmassin, totmass
584
double precision shat, sqrtshat, stot, y, m(-nexternal:nexternal)
585
integer nbranch, ns_channel,nt_channel, pos_pz
587
& totmassin, totmass,shat, sqrtshat, stot,y, m,
588
& nbranch, ns_channel,nt_channel, pos_pz
596
include 'configs_decay.inc'
598
mapext2res=map_external2res
599
c Determine number of s- and t-channel branches, at this point it
600
c includes the s-channel p1+p2
601
c write(*,*) (itree(i,-1), i=1,2)
602
c write(*,*) qmass(-1)
603
c write(*,*) qwidth(-1)
605
nbranch=nexternal_prod -2
607
do while(itree(1,-ns_channel).ne.1 .and.
608
& itree(1,-ns_channel).ne.2 .and. ns_channel.lt.nbranch)
610
ns_channel=ns_channel+1
612
ns_channel=ns_channel - 1
613
nt_channel=nbranch-ns_channel-1
614
c If no t-channles, ns_channels is one less, because we want to exclude
615
c the s-channel p1+p2
616
if (nt_channel .eq. 0 .and. nincoming .eq. 2) then
617
ns_channel=ns_channel-1
619
c Set one_body to true if it s a 2->1 process at the Born (i.e. 2->2 for the n+1-body)
620
if((nexternal-nincoming).eq.1)then
624
elseif((nexternal-nincoming).gt.1)then
628
write(*,*)'Error#1 in genps_madspin.f',nexternal,nincoming
632
c write(*,*) 'ns_channel ',ns_channel
633
c write(*,*) 'nt_channel ',nt_channel
635
! first fill the new itree for the the legs in the decay part
636
do i=-(ns_channel_decay),-1
637
itree_full(1,i) = idecay(1,i)
638
itree_full(2,i) = idecay(2,i)
639
qmass_full(i)=prmass(i)
640
qwidth_full(i)=prwidth(i)
644
c write(*,*) (itree_full(i,-1), i=1,2)
645
c write(*,*) (itree_full(i,-2), i=1,2)
646
c write(*,*) (itree_full(i,-3), i=1,2)
647
c write(*,*) (itree_full(i,-4), i=1,2)
648
c write(*,*) qmass_full(-1)
649
c write(*,*) qmass_full(-2)
650
c write(*,*) qmass_full(-3)
651
c write(*,*) qmass_full(-4)
652
c write(*,*) qwidth_full(-1)
653
c write(*,*) qwidth_full(-2)
654
c write(*,*) qwidth_full(-3)
655
c write(*,*) qwidth_full(-4)
657
! store the external momenta in the production event a
658
! new variable p that has the same labeling system as the new itree
659
do i=1, nexternal_prod
661
p(j,map_external2res(i)) = p_ext(j,i)
665
! fill the new itree with the kinematics associated with the production
666
do i=-1,-(ns_channel+nt_channel)-1,-1
667
c write(*,*) 'i prod',i
668
c write(*,*) 'i full',i-ns_channel_decay
669
c write(*,*) 'd1 prod',itree(1,i)
670
c write(*,*) 'd2 prod',itree(2,i)
671
if (itree(1,i).lt.0 ) then
672
c write(*,*) 'd1 full',itree(1,i)-ns_channel_decay
673
itree_full(1,i-ns_channel_decay) = itree(1,i)-ns_channel_decay
675
itree_full(1,i-ns_channel_decay) = map_external2res(itree(1,i))
677
if (itree(2,i).lt.0 ) then
678
c write(*,*) 'd2 full',itree(2,i)-ns_channel_decay
679
itree_full(2,i-ns_channel_decay) = itree(2,i)-ns_channel_decay
681
itree_full(2,i-ns_channel_decay) = map_external2res(itree(2,i))
684
! record the momentum of the intermediate leg
686
if (nt_channel.ne.0.and.i .lt.-ns_channel) then
687
p(j,i-ns_channel_decay)=p(j,itree_full(1,i-ns_channel_decay))-p(j,itree_full(2,i-ns_channel_decay))
689
p(j,i-ns_channel_decay)=p(j,itree_full(1,i-ns_channel_decay))+p(j,itree_full(2,i-ns_channel_decay))
693
keep_inv(i-ns_channel_decay)=.TRUE.
695
if (i.ne.(-ns_channel-nt_channel-1)) then
696
fixedinv(i-ns_channel_decay)=dot(p(0,i-ns_channel_decay),p(0,i-ns_channel_decay))
697
qmass_full(i-ns_channel_decay)=qmass(i)
698
qwidth_full(i-ns_channel_decay)=qwidth(i)
705
!write(*,*) -ns_channel-nt_channel-1
706
!write(*,*) map_external2res(itree(2,-ns_channel-nt_channel-1))
707
!write(*,*) itree(1,-ns_channel-nt_channel-1)
708
!write(*,*) itree(2,-ns_channel-nt_channel-1)
710
!write(*,*) nt_channel
713
!write(*,*) (itree_full(i,-1), i=1,2)
714
c write(*,*) (itree_full(i,-2), i=1,2)
715
c write(*,*) (itree_full(i,-3), i=1,2)
716
c write(*,*) (itree_full(i,-4), i=1,2)
717
c write(*,*) (itree_full(i,-5), i=1,2)
718
c c write(*,*) (itree_full(i,-6), i=1,2)
719
c c write(*,*) qmass_full(-1)
720
c write(*,*) qmass_full(-2)
721
c write(*,*) qmass_full(-3)
722
c write(*,*) qmass_full(-4)
723
c write(*,*) qmass_full(-5)
724
c write(*,*) qwidth_full(-1)
725
c write(*,*) qwidth_full(-2)
726
c write(*,*) qwidth_full(-3)
727
c write(*,*) qwidth_full(-4)
728
c write(*,*) qwidth_full(-5)
730
c write(*,*) 'p -2', (p(j,-2), j=0,3)
731
c write(*,*) 'p -4', (p(j,-4), j=0,3)
732
c write(*,*) 'p 9', (p(j,9), j=0,3)
733
c !overwrite the previous information
735
ns_channel= ns_channel+ns_channel_decay
736
c write(*,*) 'ns_channel ',ns_channel
737
c write(*,*) 'nt_channel ',nt_channel
738
do i =-(ns_channel+nt_channel)-1,-1
739
itree(1,i) =itree_full(1,i)
740
itree(2,i) =itree_full(2,i)
741
qmass(i)=qmass_full(i)
742
qwidth(i)=qwidth_full(i)
745
! extract the phi and m2 numbers for each t-channel branching
746
do i=-(ns_channel+nt_channel),-ns_channel-1
747
! the t-branching process has the structure A+B > 1 + 2
748
! with A and 1 the two daughters in itree
755
p2(j)=pa(j)+ pb(j)-p1(j)
757
m2_tchan(i)=dot(p2,p2)
758
if (m2_tchan(i).gt.0d0) then
759
m2_tchan(i)=sqrt(m2_tchan(i))
761
! might be negative because of numerical unstabilities
762
index_p2=itree(2,i-1)
763
if (index_p2.gt.0) then
764
m2_tchan(i)=m(index_p2)
766
write(*,*) 'Warning: m_2^2 is negative in t-channel branching '
771
pboost(j) = pa(j)+pb(j)
773
pboost(j) = -pboost(j)
777
c go to pa+pb cms system
779
call boostx(pb,pboost,pb_cms)
780
call boostx(p1,pboost,p1_cms)
781
c rotate such that pb is aligned with z axis
782
call invrot(p1_cms, pb_cms,p1_rot)
783
phi_tchan(i)=phi(p1_rot)
786
do i=-ns_channel-1, -1
787
if (keep_inv(i)) then
788
if (nt_channel.ne.0.and.i.eq.(ns_channel+1)) cycle
793
call boostx(p(0,itree(1,i)),pboost,p1)
794
call boostx(p(0,itree(2,i)),pboost,p2)
795
c write(*,*) 'p1 ', (p1(j), j=0,3)
796
c write(*,*) 'p2 ', (p2(j), j=0,3)
797
normp=sqrt(p1(1)**2+p1(2)**2+p1(3)**2)
798
cosphi_schan(i)=p1(3)/normp
805
subroutine generate_momenta_conf(jac,x,itree,qmass,qwidth,p,p_prod,mapext2res)
810
include 'nexternal.inc'
811
include 'nexternal_prod.inc'
814
double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal)
815
integer itree(2,-nexternal:-1)
816
double precision qmass(-nexternal:0),qwidth(-nexternal:0)
819
! variables to keep track of the vegas numbers for the production part
820
logical keep_inv(-nexternal:-1),no_gen
822
double precision fixedinv(-nexternal:0)
823
double precision phi_tchan(-nexternal:0),m2_tchan(-nexternal:0)
824
double precision cosphi_schan(-nexternal:0), phi_schan(-nexternal:0)
825
common /to_fixed_kin/keep_inv,no_gen, ivar, fixedinv,
826
& phi_tchan,m2_tchan,cosphi_schan, phi_schan
829
! variables associate with the PS generation
830
double precision totmassin, totmass
831
double precision shat, sqrtshat, stot, y, m(-nexternal:nexternal)
832
integer nbranch, ns_channel,nt_channel,pos_pz
834
& totmassin, totmass,shat, sqrtshat, stot,y, m,
835
& nbranch, ns_channel,nt_channel, pos_pz
838
c Masses of particles. Should be filled in setcuts.f
839
double precision pmass(nexternal)
840
common /to_mass/pmass
843
integer mapext2res(nexternal_prod) ! map (index in production) -> index in the full structure
846
double precision pb(0:3,-nexternal:nexternal)
847
& ,S(-nexternal:nexternal)
849
& ,pwgt,p_born_CHECK(0:3,nexternal), ptot(0:3)
852
double precision smax, smin, xm02,bwmdpl,bwmdmn,bwfmpl,bwfmmn
853
double precision bwdelf, BWshift
856
double precision lambda
858
double precision xbwmass3,bwfunc
859
external xbwmass3,bwfunc
862
parameter (pi=3.1415926535897932d0)
864
data firsttime/.true./
865
double precision zero
867
c Conflicting BW stuff
868
integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
869
double precision cBW_mass(-nexternal:-1,-1:1),
870
& cBW_width(-nexternal:-1,-1:1)
871
common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
874
double precision BWcut, maxBW
875
common /to_BWcut/BWcut, maxBW
885
c STEP 1: generate the initial momenta
887
if (nexternal_prod.gt.3) then
889
m(-nbranch) = sqrtshat
890
pb(0,-nbranch)= m(-nbranch)
894
c write(*,*) 'nbranch'
896
! P.A.: Shat needs to be generated accoding to a BW distribution
897
smax=min(stot*0.99D0,(qmass(-nbranch+1)+BWcut*qwidth(-nbranch+1))**2 )
898
smin=max(0.1d0,(qmass(-nbranch+1)-BWcut*qwidth(-nbranch+1))**2 )
899
xm02=qmass(-nbranch+1)**2
902
bwfmpl=atan(bwmdpl/(qmass(-nbranch+1)*qwidth(-nbranch+1)))
903
bwfmmn=atan(bwmdmn/(qmass(-nbranch+1)*qwidth(-nbranch+1)))
904
bwdelf=(bwfmpl+bwfmmn)/pi
906
s(-nbranch)=xbwmass3(x(ivar),xm02,qwidth(-nbranch+1),bwdelf
910
if (s(-nbranch).gt.smax.or.s(-nbranch).lt.smin) then
916
sqrtshat=dsqrt(s(-nbranch))
917
xjac=xjac*bwdelf/bwfunc(s(-nbranch),xm02,qwidth(-nbranch+1))
918
m(-nbranch) = dsqrt(s(-nbranch))
919
BWshift=abs(m(-nbranch)-qmass(-nbranch+1))/qwidth(-nbranch+1)
920
if (BWshift.gt.maxBW) maxBW=BWshift
921
pb(0,-nbranch)=m(-nbranch)
927
c write(*,*) qmass(-nbranch+1)
928
c write(*,*) qwidth(-nbranch+1)
929
c write(*,*) m(-nbranch)
935
if(nincoming.eq.2) then
937
call mom2cx(sqrtshat,m(1),m(2),1d0,0d0,pb(0,1),pb(0,2))
939
call mom2cx(sqrtshat,m(2),m(1),1d0,0d0,pb(0,2),pb(0,1))
948
ptot(0)=sqrtshat*cosh(y)
951
ptot(3)=sqrtshat*sinh(y)
952
call boostx(pb(0,1),ptot,pb(0,1))
953
call boostx(pb(0,2),ptot,pb(0,2))
956
c STEP 2: generate all the invariant masses of the s-channels
957
call generate_inv_mass_sch(ns_channel,itree,m,sqrtshat
958
& ,totmass,qwidth,qmass,cBW,cBW_mass,cBW_width,s,x,xjac,pass)
960
!write(*,*) 'jac s-chan ', xjac
967
c If only s-channels, also set the p1+p2 s-channel
968
if (nt_channel .eq. 0 .and. nincoming .eq. 2) then
969
s(-nbranch+1)=s(-nbranch)
970
m(-nbranch+1)=m(-nbranch) !Basic s-channel has s_hat
971
pb(0,-nbranch+1) = m(-nbranch+1)*cosh(y)!and 0 momentum
972
pb(1,-nbranch+1) = 0d0
973
pb(2,-nbranch+1) = 0d0
974
pb(3,-nbranch+1) = m(-nbranch+1)*sinh(y)
978
c STEP 3: do the T-channel branchings
980
if (nt_channel.ne.0) then
981
call generate_t_channel_branchings(ns_channel,nbranch,itree,m,s
982
& ,x,pb,xjac,xpswgt,pass)
988
!write(*,*) 'jac t-chan ', xjac
990
c STEP 4: generate momentum for all intermediate and final states
991
c being careful to calculate from more massive to less massive states
992
c so the last states done are the final particle states.
994
call fill_momenta(nbranch,nt_channel
995
& ,x,itree,m,s,pb,xjac,xpswgt,pass)
1000
!write(*,*) 'jac fill ', xjac
1009
do i = 1, nexternal_prod
1011
p_prod(j,i)=pb(j,mapext2res(i))
1021
subroutine fill_momenta(nbranch,nt_channel
1022
& ,x,itree,m,s,pb,xjac0,xpswgt0,pass)
1025
parameter (pi=3.1415926535897932d0)
1026
!include 'genps.inc'
1027
include 'nexternal.inc'
1028
integer nbranch,nt_channel,ionebody
1029
double precision M(-nexternal:nexternal),x(36)
1030
double precision s(-nexternal:nexternal)
1031
double precision pb(0:3,-nexternal:nexternal)
1032
integer itree(2,-nexternal:-1)
1033
double precision xjac0,xpswgt0
1034
logical pass,one_body
1036
double precision one
1038
double precision costh,phi,xa2,xb2
1040
double precision lambda,dot
1042
double precision vtiny
1043
parameter (vtiny=1d-12)
1045
! variables to keep track of the vegas numbers for the production part
1046
logical keep_inv(-nexternal:-1),no_gen
1048
double precision fixedinv(-nexternal:0)
1049
double precision phi_tchan(-nexternal:0),m2_tchan(-nexternal:0)
1050
double precision cosphi_schan(-nexternal:0), phi_schan(-nexternal:0)
1051
common /to_fixed_kin/keep_inv,no_gen, ivar, fixedinv,
1052
& phi_tchan,m2_tchan,cosphi_schan, phi_schan
1056
do i = -nbranch+nt_channel+(nincoming-1),-1
1057
if (keep_inv(i).or.no_gen) then
1058
costh=cosphi_schan(i)
1062
costh= 2d0*x(ivar)-1d0
1064
phi = 2d0*pi*x(ivar)
1066
cosphi_schan(i)=costh
1067
xjac0 = xjac0 * 4d0*pi
1069
xa2 = m(itree(1,i))*m(itree(1,i))/s(i)
1070
xb2 = m(itree(2,i))*m(itree(2,i))/s(i)
1071
if (m(itree(1,i))+m(itree(2,i)) .ge. m(i)) then
1078
c write(*,*) sqrt(s(i))
1079
c write(*,*) m(itree(1,i))
1080
c write(*,*) m(itree(2,i))
1081
c write(*,*) xa2, xb2
1083
if (.not.keep_inv(i).and..not.no_gen) then
1084
xpswgt0 = xpswgt0*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI)
1086
c write(*,*) xpswgt0
1088
call mom2cx(m(i),m(itree(1,i)),m(itree(2,i)),costh,phi,
1089
& pb(0,itree(1,i)),pb(0,itree(2,i)))
1090
c write(*,*) 'i ', i
1091
c write(*,*) 'costh, phi ', costh,phi
1092
c write(*,*) 'm init ', m(i)
1093
c write(*,*) (pb(j,itree(1,i)), j=0,3)
1094
c write(*,*) (pb(j,itree(2,i)), j=0,3)
1095
c If there is an extremely large boost needed here, skip the phase-space point
1096
c because of numerical stabilities.
1097
if (dsqrt(abs(dot(pb(0,i),pb(0,i))))/pb(0,i)
1103
call boostx(pb(0,itree(1,i)),pb(0,i),pb(0,itree(1,i)))
1104
call boostx(pb(0,itree(2,i)),pb(0,i),pb(0,itree(2,i)))
1105
c write(*,*) (pb(j,itree(1,i)), j=0,3)
1106
c write(*,*) (pb(j,itree(2,i)), j=0,3)
1111
c Special phase-space fix for the one_body
1112
c if (one_body) then
1113
c Factor due to the delta function in dphi_1
1114
c xpswgt0=pi/m(ionebody)
1115
c Kajantie s normalization of phase space (compensated below in flux)
1116
c xpswgt0=xpswgt0/(2*pi)
1118
c pb(i,3) = pb(i,1)+pb(i,2)
1125
subroutine generate_inv_mass_sch(ns_channel,itree,m,sqrtshat
1126
& ,totmass,qwidth,qmass,cBW,cBW_mass,cBW_width,s,x,xjac0,pass)
1129
parameter (pi=3.1415926535897932d0)
1130
!include 'genps.inc'
1131
include 'nexternal.inc'
1133
double precision totmass,BWshift
1134
double precision sqrtshat, m(-nexternal:nexternal)
1136
double precision qmass(-nexternal:0),qwidth(-nexternal:0)
1137
double precision x(36)
1138
double precision s(-nexternal:nexternal)
1139
double precision xjac0
1140
integer itree(2,-nexternal:-1)
1142
double precision smin,smax,xm02,bwmdpl,bwmdmn,bwfmpl,bwfmmn,bwdelf
1144
double precision xbwmass3,bwfunc
1145
external xbwmass3,bwfunc
1147
integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
1148
double precision cBW_mass(-nexternal:-1,-1:1),
1149
& cBW_width(-nexternal:-1,-1:1)
1150
double precision b(-1:1),x0
1152
! variables to keep track of the vegas numbers for the production part
1153
logical keep_inv(-nexternal:-1),no_gen
1155
double precision fixedinv(-nexternal:0)
1156
double precision phi_tchan(-nexternal:0),m2_tchan(-nexternal:0)
1157
double precision cosphi_schan(-nexternal:0), phi_schan(-nexternal:0)
1158
common /to_fixed_kin/keep_inv,no_gen, ivar, fixedinv,
1159
& phi_tchan,m2_tchan,cosphi_schan, phi_schan
1161
double precision BWcut, maxBW
1162
common /to_BWcut/BWcut, maxBW
1166
do i = -1,-ns_channel,-1
1167
c Generate invariant masses for all s-channel branchings of the Born
1168
if (keep_inv(i).or.no_gen) then
1172
smin = (m(itree(1,i))+m(itree(2,i)))**2
1173
smax = (sqrtshat-totalmass+sqrt(smin))**2
1174
if(smax.lt.smin.or.smax.lt.0.d0.or.smin.lt.0.d0)then
1175
write(*,*)'Error#13 in genps_madspin.f'
1176
write(*,*)smin,smax,i
1180
c Choose the appropriate s given our constraints smin,smax
1181
if(qwidth(i).ne.0.d0)then
1183
if (cBW(i).eq.1 .and.
1184
& cBW_width(i,1).gt.0d0 .and. cBW_width(i,-1).gt.0d0) then
1185
c conflicting BW on both sides
1187
b(j)=(cBW_mass(i,j)-qmass(i))/
1188
& (qwidth(i)+cBW_width(i,j))
1189
b(j)=qmass(i)+b(j)*qwidth(i)
1192
if (x(ivar).lt.1d0/3d0) then
1194
s(i)=(b(-1)-smin)*x0+smin
1195
xjac0=3d0*xjac0*(b(-1)-smin)
1196
elseif (x(ivar).gt.1d0/3d0 .and. x(ivar).lt.2d0/3d0) then
1201
bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
1202
bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
1203
bwdelf=(bwfmpl+bwfmmn)/pi
1204
s(i)=xbwmass3(x0,xm02,qwidth(i),bwdelf
1206
xjac0=3d0*xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
1209
s(i)=(smax-b(1))*x0+b(1)
1210
xjac0=3d0*xjac0*(smax-b(1))
1212
elseif (cBW(i).eq.1.and.cBW_width(i,1).gt.0d0) then
1213
c conflicting BW with alternative mass larger
1214
b(1)=(cBW_mass(i,1)-qmass(i))/
1215
& (qwidth(i)+cBW_width(i,1))
1216
b(1)=qmass(i)+b(1)*qwidth(i)
1218
if (x(ivar).lt.0.5d0) then
1223
bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
1224
bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
1225
bwdelf=(bwfmpl+bwfmmn)/pi
1226
s(i)=xbwmass3(x0,xm02,qwidth(i),bwdelf
1228
xjac0=2d0*xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
1231
s(i)=(smax-b(1))*x0+b(1)
1232
xjac0=2d0*xjac0*(smax-b(1))
1234
elseif (cBW(i).eq.1.and.cBW_width(i,-1).gt.0d0) then
1235
c conflicting BW with alternative mass smaller
1236
b(-1)=(cBW_mass(i,-1)-qmass(i))/
1237
& (qwidth(i)+cBW_width(i,-1)) ! b(-1) is negative here
1238
b(-1)=qmass(i)+b(-1)*qwidth(i)
1240
if (x(ivar).lt.0.5d0) then
1242
s(i)=(b(-1)-smin)*x0+smin
1243
xjac0=2d0*xjac0*(b(-1)-smin)
1249
bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
1250
bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
1251
bwdelf=(bwfmpl+bwfmmn)/pi
1252
s(i)=xbwmass3(x0,xm02,qwidth(i),bwdelf
1254
xjac0=2d0*xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
1258
! P.A.: introduce the BWcutoff here
1259
smax=min(smax,(qmass(i)+BWcut*qwidth(i))**2 )
1260
smin=max(smin,(qmass(i)-BWcut*qwidth(i))**2 )
1265
bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
1266
bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
1267
bwdelf=(bwfmpl+bwfmmn)/pi
1268
s(i)=xbwmass3(x(ivar),xm02,qwidth(i),bwdelf
1271
if (s(i).gt.smax.or.s(i).lt.smin) then
1277
xjac0=xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
1278
!write(*,*) i , sqrt(s(i))
1281
c not a Breit Wigner
1282
s(i) = (smax-smin)*x(ivar)+smin
1283
xjac0 = xjac0*(smax-smin)
1286
c If numerical inaccuracy, quit loop
1287
if (xjac0 .lt. 0d0) then
1292
if (s(i) .lt. smin) then
1299
c fill masses, update totalmass
1303
if (.not.keep_inv(i).and..not.no_gen) then
1304
BWshift=abs(m(i)-qmass(i))/qwidth(i)
1305
if (BWshift.gt.maxBW) maxBW=BWshift
1307
totalmass=totalmass+m(i)-
1308
& m(itree(1,i))-m(itree(2,i))
1309
if ( totalmass.gt.sqrtshat )then
1320
subroutine generate_t_channel_branchings(ns_channel,nbranch,itree
1321
& ,m,s,x,pb,xjac0,xpswgt0,pass)
1322
c First we need to determine the energy of the remaining particles this
1323
c is essentially in place of the cos(theta) degree of freedom we have
1324
c with the s channel decay sequence
1327
parameter (pi=3.1415926535897932d0)
1328
!include 'genps.inc'
1329
include 'nexternal.inc'
1330
double precision xjac0,xpswgt0
1331
double precision M(-nexternal:nexternal),x(36)
1332
double precision s(-nexternal:nexternal)
1333
double precision pb(0:3,-nexternal:nexternal)
1334
integer itree(2,-nexternal:-1)
1335
integer ns_channel,nbranch
1338
double precision totalmass,smin,smax,s1,ma2,mbq,m12,mnq,tmin,tmax
1341
double precision lambda,dot
1344
! variables to keep track of the vegas numbers for the production part
1345
logical keep_inv(-nexternal:-1),no_gen
1347
double precision fixedinv(-nexternal:0)
1348
double precision phi_tchan(-nexternal:0),m2_tchan(-nexternal:0)
1349
double precision cosphi_schan(-nexternal:0), phi_schan(-nexternal:0)
1350
common /to_fixed_kin/keep_inv,no_gen, ivar, fixedinv,
1351
& phi_tchan,m2_tchan,cosphi_schan, phi_schan
1356
do ibranch = -ns_channel-1,-nbranch,-1
1357
totalmass=totalmass+m(itree(2,ibranch))
1359
m(-ns_channel-1) = dsqrt(S(-nbranch))
1361
c Choose invariant masses of the pseudoparticles obtained by taking together
1362
c all final-state particles or pseudoparticles found from the current
1363
c t-channel propagator down to the initial-state particle found at the end
1364
c of the t-channel line.
1365
do ibranch = -ns_channel-1,-nbranch+2,-1
1366
totalmass=totalmass-m(itree(2,ibranch))
1368
smax = (m(ibranch) - m(itree(2,ibranch)))**2
1369
if (smin .gt. smax) then
1374
if (keep_inv(ibranch).or.no_gen) then
1375
m(ibranch-1)=m2_tchan(ibranch)
1378
m(ibranch-1)=dsqrt((smax-smin)*
1380
m2_tchan(ibranch)=m(ibranch-1)
1381
xjac0 = xjac0*(smax-smin)
1383
if (m(ibranch-1)**2.lt.smin.or.m(ibranch-1)**2.gt.smax
1384
& .or.m(ibranch-1).ne.m(ibranch-1)) then
1391
c Set m(-nbranch) equal to the mass of the particle or pseudoparticle P
1392
c attached to the vertex (P,t,p2), with t being the last t-channel propagator
1393
c in the t-channel line, and p2 the incoming particle opposite to that from
1394
c which the t-channel line starts
1395
m(-nbranch) = m(itree(2,-nbranch))
1397
c Now perform the t-channel decay sequence. Most of this comes from:
1398
c Particle Kinematics Chapter 6 section 3 page 166
1400
c From here, on we can just pretend this is a 2->2 scattering with
1401
c Pa + Pb -> P1 + P2
1402
c p(0,itree(ibranch,1)) + p(0,2) -> p(0,ibranch)+ p(0,itree(ibranch,2))
1403
c M(ibranch) is the total mass available (Pa+Pb)^2
1404
c M(ibranch-1) is the mass of P2 (all the remaining particles)
1406
do ibranch=-ns_channel-1,-nbranch+1,-1
1407
s1 = m(ibranch)**2 !Total mass available
1409
mbq = dot(pb(0,itree(1,ibranch)),pb(0,itree(1,ibranch)))
1410
m12 = m(itree(2,ibranch))**2
1411
mnq = m(ibranch-1)**2
1412
call yminmax(s1,t,m12,ma2,mbq,mnq,tmin,tmax)
1414
if (keep_inv(ibranch).or.no_gen) then
1415
t = fixedinv(ibranch)
1418
t=(tmax_temp-tmin)*x(ivar)+tmin
1420
xjac0=xjac0*(tmax_temp-tmin)
1423
if (t .lt. tmin .or. t .gt. tmax) then
1428
if (keep_inv(ibranch).or.no_gen) then
1429
phi=phi_tchan(ibranch)
1432
phi = 2d0*pi*x(ivar)
1433
phi_tchan(ibranch)=phi
1434
xjac0 = xjac0*2d0*pi
1437
c Finally generate the momentum. The call is of the form
1438
c pa+pb -> p1+ p2; t=(pa-p1)**2; pr = pa-p1
1439
c gentcms(pa,pb,t,phi,m1,m2,p1,pr)
1441
call gentcms(pb(0,itree(1,ibranch)),pb(0,2),t,phi,
1442
& m(itree(2,ibranch)),m(ibranch-1),pb(0,itree(2,ibranch)),
1443
& pb(0,ibranch),xjac0)
1445
if (xjac0 .lt. 0d0) then
1446
write(*,*) 'Failedgentcms',ibranch,xjac0
1450
if (.not.keep_inv(ibranch).and..not.no_gen) then
1451
xpswgt0 = xpswgt0/(4d0*dsqrt(lambda(s1,ma2,mbq)))
1454
c We need to get the momentum of the last external particle. This
1455
c should just be the sum of p(0,2) and the remaining momentum from our
1456
c last t channel 2->2
1458
pb(i,itree(2,-nbranch)) = pb(i,-nbranch+1)+pb(i,2)
1463
subroutine gentcms(pa,pb,t,phi,m1,m2,p1,pr,jac)
1464
c*************************************************************************
1465
c Generates 4 momentum for particle 1, and remainder pr
1466
c given the values t, and phi
1467
c Assuming incoming particles with momenta pa, pb
1468
c And outgoing particles with mass m1,m2
1469
c s = (pa+pb)^2 t=(pa-p1)^2
1470
c*************************************************************************
1475
double precision t,phi,m1,m2 !inputs
1476
double precision pa(0:3),pb(0:3),jac
1477
double precision p1(0:3),pr(0:3) !outputs
1481
double precision ptot(0:3),E_acms,p_acms,pa_cms(0:3)
1482
double precision esum,ed,pp,md2,ma2,pt,ptotm(0:3)
1487
double precision dot
1493
ptot(i) = pa(i)+pb(i)
1502
c determine magnitude of p1 in cms frame (from dhelas routine mom2cx)
1504
ESUM = sqrt(max(0d0,dot(ptot,ptot)))
1505
if (esum .eq. 0d0) then
1506
jac=-8d0 !Failed esum must be > 0
1511
IF (M1*M2.EQ.0.) THEN
1512
PP=(ESUM-ABS(ED))*0.5d0
1514
PP=(MD2/ESUM)**2-2.0d0*(M1**2+M2**2)+ESUM**2
1518
write(*,*) 'Warning#12 in genps_madspin.f',pp
1524
c Energy of pa in pa+pb cms system
1526
call boostx(pa,ptotm,pa_cms)
1528
p_acms = dsqrt(pa_cms(1)**2+pa_cms(2)**2+pa_cms(3)**2)
1530
p1(0) = MAX((ESUM+ED)*0.5d0,0.d0)
1531
p1(3) = -(m1*m1+ma2-t-2d0*p1(0)*E_acms)/(2d0*p_acms)
1532
pt = dsqrt(max(pp*pp-p1(3)*p1(3),0d0))
1536
call rotxxx(p1,pa_cms,p1) !Rotate back to pa_cms frame
1537
call boostx(p1,ptot,p1) !boost back to lab fram
1539
pr(i)=pa(i)-p1(i) !Return remainder of momentum
1546
function bwfunc(s,xm02,gah)
1547
c Returns the Breit Wigner function, normalized in such a way that
1548
c its integral in the range (-inf,inf) is one
1550
real*8 bwfunc,s,xm02,gah
1552
parameter (pi=3.1415926535897932d0)
1555
bwfunc=xm0*gah/(pi*((s-xm02)**2+xm02*gah**2))
1560
SUBROUTINE YMINMAX(X,Y,Z,U,V,W,YMIN,YMAX)
1561
C**************************************************************************
1562
C This is the G function from Particle Kinematics by
1563
C E. Byckling and K. Kajantie, Chapter 4 p. 91 eqs 5.28
1564
C It is used to determine physical limits for Y based on inputs
1565
C**************************************************************************
1570
double precision tiny
1571
parameter (tiny=1d-199)
1575
Double precision x,y,z,u,v,w !inputs y is dummy
1576
Double precision ymin,ymax !output
1580
double precision y1,y2,yr,ysqr
1584
double precision lambda
1588
ysqr = lambda(x,u,v)*lambda(x,w,z)
1589
if (ysqr .ge. 0d0) then
1592
print*,'Error in yminymax sqrt(-x)',lambda(x,u,v),lambda(x,w,z)
1595
y1 = u+w -.5d0* ((x+u-v)*(x+w-z) - yr)/(x+tiny)
1596
y2 = u+w -.5d0* ((x+u-v)*(x+w-z) + yr)/(x+tiny)
1602
function xbwmass3(t,xm02,ga,bwdelf,bwfmmn)
1603
c Returns the boson mass squared, given 0<t<1, the nominal mass (xm0),
1604
c and the mass range (implicit in bwdelf and bwfmmn). This function
1605
c is the inverse of F(M^2), where
1606
c F(M^2)=\int_{xmlow2}^{M^2} ds BW(sqrt(s),M0,Ga)
1607
c BW(M,M0,Ga)=M0 Ga/pi 1/((M^2-M0^2)^2+M0^2 Ga^2
1608
c and therefore eats up the Breit-Wigner when changing integration
1609
c variable M^2 --> t
1611
real*8 xbwmass3,t,xm02,ga,bwdelf,bwfmmn
1613
parameter (pi=3.1415926535897932d0)
1616
xbwmass3=xm02+xm0*ga*tan(pi*bwdelf*t-bwfmmn)
1621
subroutine invrot(p,q, pp)
1622
! inverse of the "rot" operation
1623
! q is the four momentum that is aligned with z in the new frame
1624
! p is the four momentum to be rotated
1625
! pp is rotated four momentum p
1628
double precision pp(0:3), p(0:3), q(0:3)
1630
double precision qt2, qt, qq
1633
qt2 = (q(1))**2 + (q(2))**2
1635
if(qt2.eq.0.0d0) then
1636
if ( q(3).gt.0d0 ) then
1646
qq = sqrt(qt2+q(3)**2)
1648
pp(2)=-q(2)/qt*p(1)+q(1)/qt*p(2)
1649
if (q(3).eq.0d0) then
1651
if (q(2).ne.0d0) then
1652
pp(3)=(p(2)-q(2)*q(3)/qq/qt-q(1)/qt*pp(2))*qq/q(2)
1654
pp(3)=(p(1)-q(1)*q(3)/qq/qt*pp(1)+q(2)/qt*pp(2))*qq/q(1)
1657
if (q(2).ne.0d0) then
1658
pp(3)=(qt**2*p(2)+q(2)*q(3)*p(3)-q(1)*qt*pp(2))/qq/q(2)
1660
pp(3)=(q(1)*p(1)+q(3)*p(3))/qq
1662
pp(1)=(-p(3)+q(3)/qq*pp(3))*qq/qt
1671
DOUBLE PRECISION FUNCTION phi(p)
1672
c************************************************************************
1673
c MODIF 16/11/06 : this subroutine defines phi angle
1674
c phi is defined from 0 to 2 pi
1675
c************************************************************************
1680
double precision p(0:3)
1685
double precision pi,zero
1686
parameter (pi=3.141592654d0,zero=0d0)
1691
if(p(1).gt.zero) then
1692
phi=datan(p(2)/p(1))
1693
else if(p(1).lt.zero) then
1694
phi=datan(p(2)/p(1))+pi
1695
else if(p(2).GE.zero) then !remind that p(1)=0
1697
else if(p(2).lt.zero) then !remind that p(1)=0
1700
if(phi.lt.zero) phi=phi+2*pi
1704
double precision function dot(p1,p2)
1705
C****************************************************************************
1706
C 4-Vector Dot product
1707
C****************************************************************************
1709
double precision p1(0:3),p2(0:3)
1710
dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
1712
if(dabs(dot).lt.1d-6)then ! solve numerical problem
1718
DOUBLE PRECISION FUNCTION LAMBDA(S,MA2,MB2)
1720
C****************************************************************************
1721
C THIS IS THE LAMBDA FUNCTION FROM VERNONS BOOK COLLIDER PHYSICS P 662
1722
C MA2 AND MB2 ARE THE MASS SQUARED OF THE FINAL STATE PARTICLES
1723
C 2-D PHASE SPACE = .5*PI*SQRT(1.,MA2/S^2,MB2/S^2)*(D(OMEGA)/4PI)
1724
C****************************************************************************
1725
DOUBLE PRECISION MA2,MB2,S,tiny,tmp,rat
1726
parameter (tiny=1.d-8)
1728
tmp=S**2+MA2**2+MB2**2-2d0*S*MA2-2d0*MA2*MB2-2d0*S*MB2
1730
if(ma2.lt.0.d0.or.mb2.lt.0.d0)then
1731
write(6,*)'Error#1 in function Lambda:',s,ma2,mb2
1734
rat=1-(sqrt(ma2)+sqrt(mb2))/s
1735
if(rat.gt.-tiny)then
1738
write(6,*)'Error#2 in function Lambda:',s,ma2,mb2
added
removed
MadSpin/src/driver.f
