1
subroutine add_write_info(p_born,pp,ybst_til_tolab,iconfig,Hevents
2
& ,putonshell,ndim,ipole,x,jpart,npart,pb,shower_scale)
3
c Computes all the info needed to write out the events including the
4
c intermediate resonances. It also boosts the events to the lab frame
7
include "nexternal.inc"
8
include "born_nhel.inc"
9
include "coloramps.inc"
10
include "reweight0.inc"
11
include "nFKSconfigs.inc"
12
include "leshouche_info.inc"
16
double precision p_born(0:3,nexternal-1),pp(0:3,nexternal)
17
double precision ybst_til_tolab,shower_scale
19
logical Hevents,putonshell
20
integer ndim,ipole,jpart(7,-nexternal+3:2*nexternal-3),npart
21
double precision pb(0:4,-nexternal+3:2*nexternal-3)
24
integer np,i,j,iBornGraph,ida(2),size,nexpart,idum,ires,nres,ns
25
integer icolalt(2,-nexternal+3:2*nexternal-3),ip,iflow
26
integer ito(-nexternal+3:nexternal)
27
double precision xtarget,sumborn,jampsum
28
double complex wgt1(2)
29
logical firsttime,firsttime2
30
data firsttime/.true./
31
data firsttime2/.true./
33
c The process chosen to write
35
common/c_addwrite/i_process
41
c Jamp amplitudes of the Born (to be filled with a call the sborn())
42
double Precision amp2(maxamps), jamp2(0:maxamps)
43
common/to_amps/ amp2, jamp2
45
C iforest and other configuration info. Read once and saved.
46
integer itree_S_t(2,-max_branch:-1),sprop_tree_S_t(-max_branch:-1)
47
$ ,itree_H_t(2,-max_branch:-1),sprop_tree_H_t(-max_branch:-1)
48
integer itree_S(2,-max_branch:-1,fks_configs),sprop_tree_S(
49
$ -max_branch:-1,fks_configs),itree_H(2,-max_branch:-1
50
$ ,fks_configs),sprop_tree_H(-max_branch:-1,fks_configs)
51
save itree_S,sprop_tree_S,itree_H,sprop_tree_H
53
c Masses and widths of the (internal) propagators. Read once and saved.
54
double precision pmass_tree_S_t(-nexternal:0) ,pwidth_tree_S_t(
55
$ -nexternal:0),pmass_tree_H_t(-nexternal:0),pwidth_tree_H_t(
57
double precision pmass_tree_S(-nexternal:0,fks_configs)
58
$ ,pwidth_tree_S(-nexternal:0,fks_configs),pmass_tree_H(
59
$ -nexternal:0,fks_configs),pwidth_tree_H(-nexternal:0
61
save pmass_tree_S,pwidth_tree_S,pmass_tree_H,pwidth_tree_H
63
c tree, s-channel props, masses and widths, copied from the save values
65
integer itree(2,-max_branch:-1),sprop_tree(-max_branch:-1)
66
double precision pmass_tree(-nexternal:0)
67
double precision pwidth_tree(-nexternal:0)
70
logical OnBW(-nexternal:0)
74
parameter (maxflow=999)
75
integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
76
& icolup(2,nexternal,maxflow)
77
c include "leshouche.inc"
78
common /c_leshouche_inc/idup,mothup,icolup
80
c Common block to check if we are doing MC over helicities.
83
common/to_matrix/isum_hel, multi_channel
85
c Masses of external (n+1)-body particles
86
real*8 emass(nexternal)
91
common/fks_indices/i_fks,j_fks
93
c For the boost to the lab frame
94
double precision chy,shy,chymo,p1(0:3,nexternal),xdir(3)
95
data (xdir(i),i=1,3) /0,0,1/
97
c For (n+1)-body this is the configuration mapping
98
integer mapconfig(0:lmaxconfigs), this_config
99
common/to_mconfigs/mapconfig, this_config
101
c For shifting QCD partons from zero to their mass-shell
102
double precision x(99),p(0:3,99)
104
double precision xmi,xmj,xm1,xm2,emsum,tmpecm,dot,wgt
105
double precision p1_cnt(0:3,nexternal,-2:2)
106
double precision wgt_cnt(-2:2)
107
double precision pswgt_cnt(-2:2)
108
double precision jac_cnt(-2:2)
109
common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
110
double precision xmcmass(nexternal)
111
common/cxmcmass/xmcmass
113
parameter (mohdr=-100)
117
COMMON/C_NFKSPROCESS/NFKSPROCESS
118
integer save_nFKSprocess
119
double precision SCALUP(fks_configs*2)
120
common /cshowerscale/SCALUP
121
integer iSorH_lhe,ifks_lhe(fks_configs) ,jfks_lhe(fks_configs)
122
& ,fksfather_lhe(fks_configs) ,ipartner_lhe(fks_configs)
123
common/cto_LHE1/iSorH_lhe,ifks_lhe,jfks_lhe,
124
# fksfather_lhe,ipartner_lhe
126
c Set the leshouche info and fks info
128
call fks_inc_chooser()
129
call leshouche_inc_chooser()
132
c Set the number of external particles and overwrite the iSorH_lhe value
142
c Determine which Born graph was used for multi-channeling
147
c Fill the itree, sprop, pmass and pwidth of this configuration. Needed
148
c to determine possible intermediate s-channel resonances. Note that the
149
c set_itree subroutine does not properly set the t-channel info.
152
save_nFKSprocess=nFKSprocess
153
do nFKSprocess=1,FKS_configs
154
call fks_inc_chooser()
156
call set_itree(iconfig,.false.,itree_S_t,sprop_tree_S_t
157
$ ,pmass_tree_S_t,pwidth_tree_S_t)
159
call set_itree(iconfig,.true.,itree_H_t,sprop_tree_H_t
160
$ ,pmass_tree_H_t,pwidth_tree_H_t)
162
itree_H(1,j,nFKSprocess)=itree_H_t(1,j)
163
itree_H(2,j,nFKSprocess)=itree_H_t(2,j)
164
sprop_tree_H(j,nFKSprocess)=sprop_tree_H_t(j)
165
pmass_tree_H(j,nFKSprocess)=pmass_tree_H_t(j)
166
pwidth_tree_H(j,nFKSprocess)=pwidth_tree_H_t(j)
167
itree_S(1,j,nFKSprocess)=itree_S_T(1,j)
168
itree_S(2,j,nFKSprocess)=itree_S_t(2,j)
169
sprop_tree_S(j,nFKSprocess)=sprop_tree_S_t(j)
170
pmass_tree_S(j,nFKSprocess)=pmass_tree_S_t(j)
171
pwidth_tree_S(j,nFKSprocess)=pwidth_tree_S_t(j)
175
nFKSprocess=save_nFKSprocess
176
call fks_inc_chooser()
178
c Copy the saved information to the arrays actually used
180
do j=-(nexternal-3),-1
182
itree(i,j)=itree_H(i,j,nFKSprocess)
184
sprop_tree(j)=sprop_tree_H(j,nFKSprocess)
185
pmass_tree(j)=pmass_tree_H(j,nFKSprocess)
186
pwidth_tree(j)=pwidth_tree_H(j,nFKSprocess)
189
do j=-(nexternal-4),-1
191
itree(i,j)=itree_S(i,j,nFKSprocess)
193
sprop_tree(j)=sprop_tree_S(j,nFKSprocess)
194
pmass_tree(j)=pmass_tree_S(j,nFKSprocess)
195
pwidth_tree(j)=pwidth_tree_S(j,nFKSprocess)
198
c Set the shower scale
200
shower_scale=SCALUP(nFKSprocess*2)
202
shower_scale=SCALUP(nFKSprocess*2-1)
205
c This is an (n+1)-body process (see update_unwgt_table in
206
c driver_mintMC.f). For S events it corresponds to the underlying Born
209
if (ip.lt.1 .or. ip.gt.maxproc_used) then
210
write (*,*)'ERROR #12 in add_write_info,'/
211
& /' not a well-defined process',ip,Hevents
216
c Fill jpart particle info for the final state particles of
217
c the (n+1)-body events. Color is done below.
220
jpart(1,i) = idup(i,ip)
221
jpart(2,i) = mothup(1,i,ip)
222
jpart(3,i) = mothup(2,i,ip)
223
if (i.le.nincoming) then
229
c Assume helicity summed
233
if (firsttime2 .and. isum_hel.ne.0) then
234
write (*,*) 'WARNING: for writing the events, no helicity '//
235
& 'info is used even though some info could be available.'
238
c Can be filled when doing MC over helicities...
239
c$$$ read(hel_buf,'(15i5)') (jpart(7,i),i=1,nexternal)
242
c Get color flow that is consistent with iconfig from Born
244
call sborn(p_born,wgt1)
247
if (icolamp(i,iBornGraph,1)) then
248
sumborn=sumborn+jamp2(i)
251
if (sumborn.eq.0d0) then
252
write (*,*) 'Error #1 in add_write_info:'
253
write (*,*) 'in MadFKS, sumborn should always be larger'//
254
& ' than zero, because always QCD partons around',sumborn,max_bcol
256
write (*,*) i,iBornGraph,icolamp(i,iBornGraph,1),jamp2(i)
260
xtarget=ran2()*sumborn
263
if (icolamp(1,iBornGraph,1)) then
268
do while (jampsum .lt. xtarget)
270
if (icolamp(iflow,iBornGraph,1)) then
271
jampsum=jampsum+jamp2(iflow)
274
if (iflow.gt.max_bcol) then
275
write (*,*) 'ERROR #2 in add_write_info',iflow,max_bcol
280
c Shift particle momenta to put them on the mass shell as given in the
281
c subroutine fill_MC_mshell().
285
c fills the common block with the MC masses (special treatment of i_fks
286
c and j_fks, because phase-space generation won't work in all cases).
287
call put_on_MC_mshell(Hevents,jpart,xmi,xmj,xm1,xm2,emsum)
288
c Prevents the code from crashing in the extremely rare case in which
289
c the cm energy is smaller than sum of masses - keep massless partons
290
tmpecm=min(dot(pp(0,1),pp(0,2)),
291
# dot(p_born(0,1),p_born(0,2)))
292
tmpecm=sqrt(2d0*tmpecm)
293
if(tmpecm.lt.0.99*emsum)then
294
write (*,*) 'Momenta generation for put_on_MC_mshell failed'
299
c generate a phase-space point with the MC masses
300
call generate_momenta(ndim,iconfig,wgt,x,p)
302
call set_cms_stuff(mohdr)
303
c special treament here for i_fks and j_fks masses
304
call put_on_MC_mshell_Hev(p,xmi,xmj,xm1,xm2,mfail)
305
c include initial state masses
306
if(j_fks.gt.nincoming.and.mfail.eq.0)
307
& call put_on_MC_mshell_in(p,xm1,xm2,mfail)
309
c include initial state masses
310
call set_cms_stuff(izero)
311
call put_on_MC_mshell_in(p1_cnt(0,1,0),xm1,xm2,mfail)
314
c restore the common block for the masses to the original MG masses
315
call put_on_MG_mshell()
317
c all went fine and we can copy the new momenta onto the old ones.
322
elseif(.not.Hevents .and. i.lt.max(i_fks,j_fks)) then
323
p_born(j,i)=p1_cnt(j,i,0)
324
elseif(.not.Hevents .and. i.gt.max(i_fks,j_fks)) then
325
p_born(j,i-1)=p1_cnt(j,i,0)
329
elseif(mfail.eq.1)then
330
c Probably not needed, but just to make sure: fill the momenta common
331
c blocks again by call generate momenta again.
333
call generate_momenta(ndim,iconfig,wgt,x,p)
335
call set_cms_stuff(mohdr)
337
call set_cms_stuff(izero)
339
c Also, set the masses that need to written in the event file to the MG
341
call write_masses_lhe_MG()
342
elseif(mfail.eq.-1)then
343
write(*,*)'Error in driver: mfail not set',mfail
347
c Use the MadGraph masses in the event file
348
call write_masses_lhe_MG()
352
c Derive the n-body from the (n+1)-body if we are doing S-events
354
if (.not.Hevents) then
356
if(i.lt.min(i_fks,j_fks)) then
357
jpart(1,i) = jpart(1,i)
358
jpart(2,i) = jpart(2,i)
359
jpart(3,i) = jpart(3,i)
360
elseif(i.eq.min(i_fks,j_fks)) then
361
if(jpart(1,i_fks).eq.-jpart(1,j_fks)
362
& .and.j_fks.gt.nincoming) then
364
elseif(jpart(1,i_fks).eq.jpart(1,j_fks)
365
& .and.j_fks.le.nincoming) then
367
elseif(jpart(1,i_fks).eq.21) then
368
jpart(1,i)=jpart(1,j_fks)
369
elseif(jpart(1,j_fks).eq.21.and.j_fks.le.nincoming) then
370
jpart(1,i)=-jpart(1,i_fks)
372
write (*,*) 'ERROR #5 in add_write_info()',
373
& i_fks,j_fks,jpart(1,i_fks),jpart(1,j_fks)
376
jpart(2,i) = jpart(2,i)
377
jpart(3,i) = jpart(3,i)
378
elseif(i.lt.max(i_fks,j_fks)) then
379
jpart(1,i) = jpart(1,i)
380
jpart(2,i) = jpart(2,i)
381
jpart(3,i) = jpart(3,i)
383
jpart(1,i) = jpart(1,i+1)
384
jpart(2,i) = jpart(2,i+1)
385
jpart(3,i) = jpart(3,i+1)
391
c Fill color of external particles
394
call fill_icolor_H(iflow,jpart)
396
call fill_icolor_S(iflow,jpart,idum)
399
icolalt(1,i)=jpart(4,i)
400
icolalt(2,i)=jpart(5,i)
404
c Set-up the external momenta that should be written in event file
405
c Also boost them to the lab frame.
407
if (abs(ybst_til_tolab).le.1d-7) then
416
if (Hevents .or. i.lt.i_fks) then
423
chy=cosh(ybst_til_tolab)
424
shy=sinh(ybst_til_tolab)
428
call boostwdir2(chy,shy,chymo,xdir,
431
call boostwdir2(chy,shy,chymo,xdir,
432
& p_born(0,i),p1(0,i))
437
if (Hevents .or. i.lt.i_fks) then
443
c In some rare cases (due to numerical inaccuracies in the boost) the
444
c energy of the incoming partons can be larger than the beam energy This
445
c is, of course, non-physical, so we use the non-boosted events
446
if (pb(0,1).gt.ebeam(1).or.pb(0,2).gt.ebeam(2)) then
447
write (*,*) 'WARNING: boost from center-of-momentum to '/
448
& /'laboratory frame too extreme. Use the center-of-mo'/
449
& /'mentum momenta instead.',pb(0,1),pb(0,2),ebeam(1)
463
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
464
C ALL EXTERNAL PARTICLE INFO SET-UP. CHECK WHICH RESONANCES SHOULD BE WRITTEN C
465
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
467
c Fill the OnBW array to determine which resonances should be written
469
call OnBreitWigner(pp,p_born,Hevents,itree,sprop_tree,pmass_tree,
472
c First check number of resonant s-channel propagators
476
c Loop over propagators to find mother-daughter information
477
do i=-1,-nexpart+3,-1
481
c Skip the t-channels
482
if ( itree(1,i) .eq. 1 .or. itree(2,i) .eq. 1 .or.
483
& itree(1,i) .eq. 2 .or. itree(2,i) .eq. 2 ) exit
484
jpart(1,i)=sprop_tree(i)
486
c Set status codes for propagator
488
c Resonance whose mass should be preserved
492
c Tag the s-channel internal propagators that we'll remove later
495
c Calculate momentum (p1+p2 for s-channel)
497
pb(j,i)=pb(j,ida(1))+pb(j,ida(2))
500
& sqrt(max(0d0,pb(0,i)**2-pb(1,i)**2-pb(2,i)**2-pb(3,i)**2))
502
c Set color info for all s-channels
503
c Fist set "safe" color info
504
if(icolalt(1,ida(1))+icolalt(1,ida(2))-
505
$ icolalt(2,ida(1))-icolalt(2,ida(2)).eq.0) then ! color singlet
508
elseif(icolalt(1,ida(1))-icolalt(2,ida(2)).eq.0) then
509
icolalt(1,i) = icolalt(1,ida(2))
510
icolalt(2,i) = icolalt(2,ida(1))
511
else if(icolalt(1,ida(2))-icolalt(2,ida(1)).eq.0) then
512
icolalt(1,i) = icolalt(1,ida(1))
513
icolalt(2,i) = icolalt(2,ida(2))
514
else if(jpart(6,i).ge.3) then ! Don't need to match
515
icolalt(1,i) = icolalt(1,ida(1))+icolalt(1,ida(2))
516
icolalt(2,i) = icolalt(2,ida(1))+icolalt(2,ida(2))
518
c Erraneous color assignment for propagator
519
write(*,*) 'ERROR: Safe color assignment wrong!'/
520
& /' This should never happen !',i,icolalt(1,ida(1))
521
& ,icolalt(1,ida(2)),icolalt(2,ida(1)),icolalt(2,ida(2))
525
c Set initial state as tentative mothers
528
c Set mother info for daughters
533
c Just zero helicity info for intermediate states
535
enddo ! do i (loop over internal propagators)
537
c Remove non-resonant mothers, set position of particles
540
jpart(4,i)=icolalt(1,i)
541
jpart(5,i)=icolalt(2,i)
542
if(i.eq.1.or.i.eq.2) then
543
ito(i)=i ! initial state particle
545
ito(i)=i+nres ! final state particle
546
elseif(i.le.-1.and.jpart(6,i).eq.2) then
548
ito(i)=2+ires ! s-channel resonances
553
if(jpart(2,i).lt.0.and.jpart(6,jpart(2,i)).ne.2) then
554
jpart(2,i)=jpart(2,jpart(2,i))
555
jpart(3,i)=jpart(3,jpart(3,i))
560
c Shift particles to right place
563
if(ito(i).le.0) cycle
565
jpart(j,ito(i))=jpart(j,i)
567
if(jpart(2,ito(i)).lt.0) then
568
jpart(2,ito(i))=ito(jpart(2,ito(i)))
569
jpart(3,ito(i))=ito(jpart(3,ito(i)))
576
c Set the number of particles that needs to be written in event file
584
subroutine set_itree(iconfig,Hevents,itree,sprop_tree,pmass_tree
590
include 'nexternal.inc'
592
double precision ZERO
594
integer itree(2,-max_branch:-1)
595
integer sprop_tree(-max_branch:-1)
596
integer iforest(2,-max_branch:-1,lmaxconfigs)
597
integer sprop(-max_branch:-1,lmaxconfigs),mapconfig(0:lmaxconfigs)
598
integer tprid(-max_branch:-1,lmaxconfigs)
599
include "born_conf.inc"
601
double precision pmass(-nexternal:0,lmaxconfigs)
602
double precision pwidth(-nexternal:0,lmaxconfigs)
603
integer pow(-nexternal:0,lmaxconfigs)
604
double precision pmass_tree(-nexternal:0)
605
double precision pwidth_tree(-nexternal:0)
607
common/fks_indices/i_fks,j_fks
608
include "born_props.inc"
610
c Do not really care about t-channels: loop should just go to
611
c nexternal-4, even though there is one more when there are t-channels
613
do j=-(nexternal-4),-1
615
itree(i,j)=iforest(i,j,iconfig)
617
sprop_tree(j)=sprop(j,iconfig)
618
pmass_tree(j)=pmass(j,iconfig)
619
pwidth_tree(j)=pwidth(j,iconfig)
622
c When we are doing H-events, we need to add --when j_fks is final
623
c state-- the s-channel branching fks_mother -> j_fks + i_fks. When
624
c j_fks is initial state, we need to add a 'bogus' t-channel splitting
625
c to make sure that all the loops stop properly
628
c must re-label the external particles to get the correct daughters
629
if (i_fks.le.j_fks) then
630
write (*,*) 'ERROR: i_fks should be greater than j_fks'
633
do j=-(nexternal-4),-1
635
if ( itree(i,j).ge.i_fks ) then
636
itree(i,j)=itree(i,j)+1
641
if (j_fks.gt.nincoming) then
642
c we must add an extra s-channel. Easiest is to add it all the way at
643
c the beginning. Therefore, relabel everything else first. Use the
644
c original ones to make sure we are doing it correctly (and not
645
c accidentally overwriting something).
646
do j=-(nexternal-4),-1
647
sprop_tree(j-1)=sprop(j,iconfig)
648
pmass_tree(j-1)=pmass(j,iconfig)
649
pwidth_tree(j-1)=pwidth(j,iconfig)
651
itree(i,j-1)=iforest(i,j,iconfig)
652
c Also update the internal references
653
if ( itree(i,j-1).lt. 0 ) then
654
itree(i,j-1)=itree(i,j-1)-1
659
c Add the new s-channel
663
c This will never be an on-shell s-channel. Hence, give it some bogus values:
668
c We have to make sure that the fks_mother (which is equal to the j_fks
669
c label of the Born) is replaced by the new s-channel
671
do j=-(nexternal-3),-2 ! do not include the new s-channel
673
if ( itree(i,j).eq. j_fks ) then
674
itree(i,j)=-1 ! reference to the new s-channel
680
c j_fks is initial state
681
jj=-(nexternal-3) ! Just add it at the end
682
c setting itree to 1 (or 2) makes sure that the loops over s-channel
683
c propagators will exit
698
subroutine OnBreitWigner(p,p_born,Hevents,itree,sprop_tree,
700
c*****************************************************************************
701
c Decides if internal s-channel propagator is on-shell
702
c*****************************************************************************
708
include 'nexternal.inc'
713
double precision p(0:3,nexternal),p_born(0:3,nexternal-1)
714
logical OnBW(-nexternal:0),Hevents
715
integer itree(2,-max_branch:-1),sprop_tree(-max_branch:-1)
716
double precision pmass(-nexternal:0)
717
double precision pwidth(-nexternal:0)
721
double precision xp(0:3,-nexternal:nexternal)
724
double precision xmass
725
integer ida(2),idenpart
726
integer IDUP(nexternal)
746
if (i.ne.nexternal) then
748
elseif (i.eq.nexternal) then
755
do i=-1,-iloop,-1 !Loop over propagators
757
c Skip the t-channels
758
if ( itree(1,i) .eq. 1 .or. itree(2,i) .eq. 1 .or.
759
& itree(1,i) .eq. 2 .or. itree(2,i) .eq. 2 ) exit
761
xp(j,i) = xp(j,itree(1,i))+xp(j,itree(2,i))
763
if (pwidth(i) .gt. 0d0) then !This is B.W.
765
c If the invariant mass is close to pole mass, set OnBW to true
767
xmass = sqrt(dot(xp(0,i),xp(0,i)))
768
onshell = ( abs(xmass-pmass(i)) .lt. bwcutoff*pwidth(i) )
771
c If mother and daughter have the same ID, remove one of them
776
if (sprop_tree(i).eq.sprop_tree(ida(j))) then
777
idenpart=ida(j) ! mother and daugher have same ID
779
elseif (ida(j).gt.0) then
780
if (sprop_tree(i).eq.IDUP(ida(j))) then
781
idenpart=ida(j) ! mother and daugher have same ID
785
c Always remove if daughter final-state (and identical)
786
if(idenpart.gt.0) then
788
c Else remove either this resonance or daughter,
789
c whichever is closer to mass shell
790
elseif(idenpart.lt.0.and.abs(xmass-pmass(i)).gt.
791
$ abs(sqrt(dot(xp(0,idenpart),xp(0,idenpart)))-
793
OnBW(i)=.false. ! mother off-shell
794
elseif(idenpart.lt.0) then
795
OnBW(idenpart)=.false. ! daughter off-shell
802
subroutine get_ID_H(IDUP_tmp)
805
include 'nexternal.inc'
807
parameter (maxflow=999)
808
integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
809
& icolup(2,nexternal,maxflow)
810
c include 'leshouche.inc'
811
common /c_leshouche_inc/idup,mothup,icolup
812
integer IDUP_tmp(nexternal),i
815
IDUP_tmp(i)=IDUP(i,1)
821
subroutine get_ID_S(IDUP_tmp)
824
include 'nexternal.inc'
826
parameter (maxflow=999)
827
integer idup(nexternal,maxproc)
828
integer mothup(2,nexternal,maxproc)
829
integer icolup(2,nexternal,maxflow)
830
include 'born_leshouche.inc'
831
integer IDUP_tmp(nexternal),i
834
IDUP_tmp(i)=IDUP(i,1)
836
IDUP_tmp(nexternal)=0
842
subroutine fill_icolor_H(iflow,jpart)
844
include "nexternal.inc"
846
integer fks_j_from_i(nexternal,0:nexternal)
847
& ,particle_type(nexternal),pdg_type(nexternal)
848
common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
851
common/fks_indices/i_fks,j_fks
852
double precision xtarget,ran2
854
integer jpart(7,-nexternal+3:2*nexternal-3),iflow
855
integer i_part,j_part,imother,lc
857
call fill_icolor_S(iflow,jpart,lc)
859
j_part = particle_type(j_fks)
860
i_part = particle_type(i_fks)
863
if (i.gt.max(i_fks,j_fks)) then
864
jpart(4,i)=jpart(4,i-1)
865
jpart(5,i)=jpart(5,i-1)
869
imother=min(i_fks,j_fks)
870
c The following works only if i_fks is always greater than j_fks.
871
if (j_fks.gt.nincoming) then
872
if (j_part.eq.3.and.i_part.eq.-3) then
873
if(jpart(4,imother).eq.0 .or. jpart(5,imother).eq.0)then
874
write (*,*) 'Error #3 in fill_icolor_H',
875
& jpart(4,imother),jpart(5,imother)
879
jpart(5,i_fks)=jpart(5,imother)
880
jpart(4,j_fks)=jpart(4,imother)
882
elseif (j_part.eq.-3.and.i_part.eq.3) then
883
if(jpart(4,imother).eq.0 .or. jpart(5,imother).eq.0)then
884
write (*,*) 'Error #4 in fill_icolor_H',
885
& jpart(4,imother),jpart(5,imother)
888
jpart(4,i_fks)=jpart(4,imother)
891
jpart(5,j_fks)=jpart(5,imother)
892
elseif (j_part.eq.3.and.i_part.eq.8) then
893
if(jpart(4,imother).eq.0 .or. jpart(5,imother).ne.0)then
894
write (*,*) 'Error #5 in fill_icolor_H',
895
& jpart(4,imother),jpart(5,imother)
898
jpart(4,i_fks)=jpart(4,imother)
902
elseif (j_part.eq.-3.and.i_part.eq.8) then
903
if(jpart(4,imother).ne.0 .or. jpart(5,imother).eq.0)then
904
write (*,*) 'Error #6 in fill_icolor_H',
905
& jpart(4,imother),jpart(5,imother)
909
jpart(5,i_fks)=jpart(5,imother)
912
elseif (j_part.eq.8.and.i_part.eq.8) then
913
if(jpart(4,imother).eq.0 .or. jpart(5,imother).eq.0)then
914
write (*,*) 'Error #7 in fill_icolor_H',
915
& jpart(4,imother),jpart(5,imother)
918
if (ran2().gt.0.5d0) then
920
jpart(5,i_fks)=jpart(5,imother)
921
jpart(4,j_fks)=jpart(4,imother)
924
jpart(4,i_fks)=jpart(4,imother)
927
jpart(5,j_fks)=jpart(5,imother)
930
write (*,*) 'Error #1 in fill_icolor_H',i_part,j_part
934
if (j_part.eq.-3.and.i_part.eq.-3) then
935
if(jpart(4,imother).eq.0 .or. jpart(5,imother).eq.0)then
936
write (*,*) 'Error #8 in fill_icolor_H',
937
& jpart(4,imother),jpart(5,imother)
941
jpart(5,i_fks)=jpart(4,imother)
943
jpart(5,j_fks)=jpart(5,imother)
944
elseif (j_part.eq.3.and.i_part.eq.3) then
945
if(jpart(4,imother).eq.0 .or. jpart(5,imother).eq.0)then
946
write (*,*) 'Error #9 in fill_icolor_H',
947
& jpart(4,imother),jpart(5,imother),'HERE'
950
jpart(4,i_fks)=jpart(5,imother)
952
jpart(4,j_fks)=jpart(4,imother)
954
elseif (j_part.eq.3.and.i_part.eq.8) then
955
if(jpart(4,imother).eq.0 .or. jpart(5,imother).ne.0)then
956
write (*,*) 'Error #10 in fill_icolor_H',
957
& jpart(4,imother),jpart(5,imother)
961
jpart(5,i_fks)=jpart(4,imother)
964
elseif (j_part.eq.-3.and.i_part.eq.8) then
965
if(jpart(4,imother).ne.0 .or. jpart(5,imother).eq.0)then
966
write (*,*) 'Error #11 in fill_icolor_H',
967
& jpart(4,imother),jpart(5,imother)
970
jpart(4,i_fks)=jpart(5,imother)
974
elseif (j_part.eq.8.and.i_part.eq.3) then
975
if(jpart(4,imother).ne.0 .or. jpart(5,imother).eq.0)then
976
write (*,*) 'Error #12 in fill_icolor_H',
977
& jpart(4,imother),jpart(5,imother)
983
jpart(5,j_fks)=jpart(5,imother)
984
elseif (j_part.eq.8.and.i_part.eq.-3) then
985
if(jpart(4,imother).eq.0 .or. jpart(5,imother).ne.0)then
986
write (*,*) 'Error #13 in fill_icolor_H',
987
& jpart(4,imother),jpart(5,imother)
992
jpart(4,j_fks)=jpart(4,imother)
994
elseif (j_part.eq.8.and.i_part.eq.8) then
995
if(jpart(4,imother).eq.0 .or. jpart(5,imother).eq.0)then
996
write (*,*) 'Error #14 in fill_icolor_H',
997
& jpart(4,imother),jpart(5,imother)
1000
if (ran2().gt.0.5d0) then
1002
jpart(5,i_fks)=jpart(4,imother)
1004
jpart(5,j_fks)=jpart(5,imother)
1006
jpart(4,i_fks)=jpart(5,imother)
1008
jpart(4,j_fks)=jpart(4,imother)
1012
write (*,*) 'Error #2 in fill_icolor_H',i_part,j_part
1021
subroutine fill_icolor_S(iflow,jpart,lc)
1024
include 'nexternal.inc'
1026
parameter (maxflow=999)
1028
integer idup(nexternal,maxproc)
1029
integer mothup(2,nexternal,maxproc)
1030
integer icolup(2,nexternal,maxflow)
1031
include "born_leshouche.inc"
1032
integer jpart(7,-nexternal+3:2*nexternal-3),lc,iflow
1036
jpart(4,i)=ICOLUP(1,i,iflow)
1037
lc=max(lc,jpart(4,i))
1038
jpart(5,i)=ICOLUP(2,i,iflow)
1039
lc=max(lc,jpart(5,i))
1047
subroutine put_on_MC_mshell(Hevents,jpart,xmi,xmj,xm1,xm2,emsum)
1048
c Sets the common block /to_mass/emass, to be passed to one_tree
1049
c to generate a massive kinematics with efficiency 1.
1051
c ip is the number of the partonic process chosen at random in the case
1052
c of multiple possibilities.
1053
c This routines assumes that the mother of (i_fks,j_fks) has
1054
c label min(i_fks,j_fks)
1057
double precision xmi,xmj,xm1,xm2,emsum
1060
integer i,j,idpart,idparti,idpartj
1061
double precision zero,tmpmass
1062
parameter (zero=0.d0)
1063
c jpart contains ID of particles
1064
include 'nexternal.inc'
1065
integer jpart(7,-nexternal+3:2*nexternal-3)
1067
c Masses to be given to genps_fks.f
1068
double precision emass(nexternal)
1069
common/to_mass/ emass
1071
c Monte Carlo masses: use PDG conventions
1072
double precision mcmass(-5:21)
1073
common/cmcmass/mcmass
1075
c Masses used by write_events_lhe
1076
double precision xmcmass(nexternal)
1077
common/cxmcmass/xmcmass
1079
c cuts.inc contains maxjetflavor
1083
common/fks_indices/i_fks,j_fks
1085
c Masses of the real process, as set by MadGraph
1087
double precision pmass(nexternal)
1090
if(j_fks.ne.min(i_fks,j_fks))then
1091
write(*,*)'j_fks#min(i_fks,j_fks) in put_on_MC_mshell'
1101
if(pmass(i).ne.0.d0)then
1102
write(*,*)'Fatal error in put_on_MC_mshell',i_fks,i,pmass(i)
1106
elseif(i.eq.j_fks)then
1107
idpartj=jpart(1,j_fks)
1108
idparti=jpart(1,i_fks)
1109
if(.not.Hevents)then
1112
if(idparti.eq.-idpartj.and.j_fks.gt.nincoming) then
1114
elseif(idparti.eq.idpartj.and.j_fks.le.nincoming) then
1116
elseif(idparti.eq.21) then
1118
elseif(idpartj.eq.21.and.j_fks.le.nincoming) then
1121
write(*,*)'Error #2 in put_on_MC_mshell',
1122
& i_fks,j_fks,idparti,idpartj
1125
if( (abs(idpart).ge.1.and.abs(idpart).le.3) .or.
1127
# ( (abs(idpart).eq.4.or.abs(idpart).eq.5) .and.
1128
# pmass(j_fks).eq.0.d0 ) )then
1131
c j_fks is an heavy particle
1132
if(idparti.ne.21)then
1133
write(*,*)'Error #3 in put_on_MC_mshell',
1134
& i_fks,j_fks,i,pmass(j_fks)
1137
c If MadFKS has a non-zero mass for a c or b quark, one probably wants
1138
c to use that in the shower phase as well (ie bottom or charm production)
1141
if(j_fks.gt.nincoming)then
1144
c one_tree assumes massless incoming QCD particles
1156
if( (abs(idpartj).ge.1.and.abs(idpartj).le.3) .or.
1157
# idpartj.eq.21 .or.
1158
# ( (abs(idpartj).eq.4.or.abs(idpartj).eq.5) .and.
1159
# pmass(j_fks).eq.0.d0 ) )then
1162
c j_fks is an heavy particle
1163
if(idparti.ne.21)then
1164
write(*,*)'Error #3 in put_on_MC_mshell',
1165
& i_fks,j_fks,i,pmass(j_fks),(jpart(1,j),j=1,nexternal)
1168
c If MadFKS has a non-zero mass for a c or b quark, one probably wants
1169
c to use that in the shower phase as well (ie bottom or charm production)
1172
if(j_fks.gt.nincoming)then
1173
emass(j_fks)=xmi+xmj
1175
c one_tree assumes massless incoming QCD particles
1188
if( idpart.eq.21 .or.
1189
# (abs(idpart).ge.1.and.abs(idpart).le.5) )then
1190
if(pmass(i).eq.0.d0)then
1191
tmpmass=mcmass(idpart)
1193
c If MadFKS has a non-zero mass for a "light" quark, one probably wants
1194
c to use that in the shower phase as well (ie bottom or charm production).
1195
c One may use equivalently a condition on maxjetflavor
1201
if(i.gt.nincoming)then
1215
do i=nincoming+1,nexternal
1216
emsum=emsum+emass(i)
1219
if( xmi.eq.-1.d0.or.xmj.eq.-1.d0 .or.
1220
# xm1.eq.-1.d0.or.xm2.eq.-1.d0 )then
1221
write(*,*)'Error #4 in put_on_MC_mshell',i_fks,j_fks
1222
write(*,*)xmi,xmj,xm1,xm2
1230
subroutine put_on_MC_mshell_Hev(p,xmi,xmj,xm1,xm2,mfail)
1232
include 'nexternal.inc'
1233
double precision p(0:3,99),xmi,xmj,xm1,xm2
1236
common/fks_indices/i_fks,j_fks
1238
if (p(0,1).lt.0d0) then
1240
write (*,*) 'Momenta generation for put_on_MC_mshell failed'
1244
if(j_fks.le.nincoming)then
1245
call put_on_MC_mshell_Hevin(p,xmi,xm1,xm2,mfail)
1247
call put_on_MC_mshell_Hevout(p,xmi,xmj,mfail)
1254
subroutine put_on_MC_mshell_in(p,xm1,xm2,mfail)
1256
include 'nexternal.inc'
1257
double precision p(0:3,nexternal),xm1,xm2
1259
double precision xm1_r,xm2_r
1261
double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
1262
common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
1264
double precision xmcmass(nexternal)
1265
common/cxmcmass/xmcmass
1267
if (p(0,1).lt.0d0) then
1269
write (*,*) 'Momenta generation for put_on_MC_mshell failed'
1273
if(abs(p(3,1)+p(3,2)).gt.1.d-10)then
1274
write(*,*)'Error #1 in put_on_MC_mshell_in',p(3,1),p(3,2)
1277
if(shat.le.(xm1+xm2)**2)then
1283
c$$$CHECK AGAIN USE OF ybst_til_tolab IN getxmss.
1284
c$$$MUST BE THE SAME BOOST AS WHEN WRITING EVENTS
1285
call getxmss_madfks(shat,ybst_til_tolab,
1286
# p(3,1),xm1_r,p(3,2),xm2_r,
1287
# p(3,1),p(3,2),mfail)
1289
p(0,1)=sqrt(xm1_r**2+p(3,1)**2)
1290
p(0,2)=sqrt(xm2_r**2+p(3,2)**2)
1299
subroutine getxmss_madfks(shat,ycm,p13cm,xm1,p23cm,xm2,
1301
c This routine is taken from the MC@NLO package. It is identical
1302
c to the routine getxmss() there, except for the fact that here,
1303
c by setting donotforce=.true., the partons are left massless if
1304
c putting them on the MC mass shell implies that they travel in
1305
c the same direction.
1306
c NOTE: the convention on the sign of the boost is opposite to that
1307
c used in MC@NLO, hence the different definition of ytmp here
1309
c Original comment in getxmss():
1310
c After putting the momenta on shell, the two incoming partons may
1311
c travel in the same direction. This routine prevents this to happen,
1312
c redefining Monte Carlo masses if necessary
1314
real*8 shat,ycm,p13cm,xm1,p23cm,xm2,p13,p23
1316
real*8 tiny,fact,sqs,xm1s,xm2s,xkp2prime_norm2,xkp2prime_norm,
1317
# ytmp,e1,e2,p13p,p23p,s1p,s2p,xif,sol
1318
integer iflag,idone,ileg
1319
parameter (fact=0.98d0)
1320
parameter (tiny=1.d-6)
1322
parameter (donotforce=.false.)
1327
c NOTE: was ytmp=-ycm in MC@NLO owing to different conventions
1331
xkp2prime_norm2=( shat-2*(xm1**2+xm2**2)+
1332
# (xm1**2-xm2**2)**2/shat )/4.d0
1333
xkp2prime_norm=sqrt(xkp2prime_norm2)
1334
if(sign(1.d0,p13cm).ne.1.d0.or.sign(1.d0,p23cm).ne.-1.d0)then
1335
write(*,*)'Error # 0 in getxmss_madfks'
1340
e1=sqrt(p13**2+xm1**2)
1341
e2=sqrt(p23**2+xm2**2)
1342
p13p=p13*cosh(ytmp)-e1*sinh(ytmp)
1343
p23p=p23*cosh(ytmp)-e2*sinh(ytmp)
1347
if(s1p.eq.1.d0 .and. s2p.eq.-1.d0)then
1349
elseif(s1p.eq.-1.d0 .and. s2p.eq.-1.d0)then
1350
if(ytmp.lt.0.d0)then
1351
write(*,*)'getxmss_madfks: wrong y sign, # 1'
1356
elseif(s1p.eq.1.d0 .and. s2p.eq.1.d0)then
1357
if(ytmp.gt.0.d0)then
1358
write(*,*)'getxmss_madfks: wrong y sign, # 2'
1364
write(*,*)'Error # 1 in getxmss_madfks',s1p,s2p
1365
write(*,*)shat,xm1,xm2
1366
write(*,*)p13,e1,p23,e2,ytmp
1369
if(iflag.eq.0.and.(.not.donotforce))then
1370
sol=xif+cosh(2*ytmp)-
1371
# sqrt(2.d0)*cosh(ytmp)*sqrt(cosh(2*ytmp)-1+2*xif)
1372
if(sol.le.0.d0.or.idone.eq.1)then
1373
c The procedure failed; pass the massless event to the Monte Carlo,
1374
c and let the Monte Carlo deal with it
1383
xm1=fact*sqrt(sol*shat)
1385
write(*,*)'Mass # 1 too large in getxmss_madfks'
1388
elseif(ileg.eq.2)then
1389
xm2=fact*sqrt(sol*shat)
1391
write(*,*)'Mass # 2 too large in getxmss_madfks'
1395
write(*,*)'Error # 2 in getxmss_madfks'
1406
subroutine put_on_MC_mshell_Hevout(p,xmi,xmj,mfail)
1408
double precision p(0:3,99),xmi,xmj
1411
double precision dirbst(1:3),dirpj(1:3)
1412
double precision pipj(0:3),pibst(0:3),pjbst(0:3)
1414
double precision p1(0:3),p1R(0:3),pipjR(0:3)
1416
double precision E1o(-1:1),p1o(-1:1),E2o(-1:1),p2o(-1:1)
1419
integer i,i1vec,i2vec,isol
1420
double precision dot,rho,threedot,Q,Q2,Q0,Qv,xm1,xm2,
1421
# cosphi_pipj,cosphi_p1R,costh_pipj,costh_p1R,projj,proji,
1422
# phi_p1R,phi_pipj,sinth_p1R,sinphi_pipj,sinphi_p1R,sinth_pipj,
1424
external dot,rho,threedot
1425
include 'nexternal.inc'
1427
common/fks_indices/i_fks,j_fks
1429
if(j_fks.le.nincoming)then
1430
write(*,*)'put_on_MC_mshell_Hevout must not be called for ISR'
1435
pipj(i)=p(i,i_fks)+p(i,j_fks)
1439
if(Q2.lt.(xmi+xmj)**2)then
1440
write(*,*)'Error in put_on_MC_mshell_Hevout',pipj(0),xmi,xmj
1445
proji=threedot(pipj,p(0,i_fks))
1446
projj=threedot(pipj,p(0,j_fks))
1447
if(proji.gt.projj)then
1461
call getangles(pipj,th_pipj,costh_pipj,sinth_pipj,
1462
# phi_pipj,cosphi_pipj,sinphi_pipj)
1463
call trp_rotate_invar(pipj,pipjR,
1464
# costh_pipj,sinth_pipj,
1465
# cosphi_pipj,sinphi_pipj)
1466
if( abs(pipjR(0)-Q0).gt.max(Q0,1.d0)*1.d-8 .or.
1467
# abs(pipjR(3)-Qv).gt.max(Qv,1.d0)*1.d-8 )then
1468
write(*,*)'Error #1 in put_on_MC_mshell_Hevout',pipj,Q0,Qv
1471
call trp_rotate_invar(p1,p1R,
1472
# costh_pipj,sinth_pipj,
1473
# cosphi_pipj,sinphi_pipj)
1474
call getangles(p1R,th_p1R,costh_p1R,sinth_p1R,
1475
# phi_p1R,cosphi_p1R,sinphi_p1R)
1476
call xkin_2body(Q0,Qv,xm1,xm2,costh_p1R,E1o,p1o,E2o,p2o,ifail)
1477
if(ifail(1).eq.1.and.ifail(-1).eq.1)then
1480
elseif(ifail(1).eq.1.and.ifail(-1).eq.0)then
1482
elseif(ifail(1).eq.0.and.ifail(-1).eq.1)then
1485
if(abs(p1R(0)-E1o(1)).lt.abs(p1R(0)-E1o(-1)))then
1496
call rotate_invar(p1R,p1,
1497
# costh_p1R,sinth_p1R,
1498
# cosphi_p1R,sinphi_p1R)
1499
call rotate_invar(p1,p1,
1500
# costh_pipj,sinth_pipj,
1501
# cosphi_pipj,sinphi_pipj)
1504
p(i,i2vec)=pipj(i)-p1(i)
1512
subroutine xkin_2body(Q0,Qv,xm1,xm2,cth1,E1,p1,E2,p2,ifail)
1513
c Returns energies and moduli of 3-momenta of the two final-state particles
1514
c in a two-body phase space, in a frame where the incoming momentum has
1515
c energy equal to Q0 and modulus of 3-momentum equal to Qv.
1516
c cth1 is the cosine of the angle between \vec{Q} and \vec{p1}
1518
c There are two possible solutions, returned in arrays E*(#) and p*(#)
1519
c with #=-1,1. If a given solution is acceptable, ifail(#)=0,
1520
c and ifail(#)=1 otherwise
1522
double precision Q0,Qv,xm1,xm2,cth1
1523
double precision E1(-1:1),p1(-1:1),E2(-1:1),p2(-1:1)
1525
double precision Q02,Qv2,delta,den,arg,xA,xB
1532
delta=(Q02-Qv2)+xm1**2-xm2**2
1534
arg=delta**2-4*xm1**2*den
1541
xB=Qv*abs(cth1)*sqrt(arg)/(2*den)
1544
if(cth1.gt.0.d0)then
1545
if( (delta-2*Q0*E1( 1)).gt.0.d0 )ifail( 1)=1
1546
if( (delta-2*Q0*E1(-1)).gt.0.d0 )ifail(-1)=1
1547
elseif(cth1.gt.0.d0)then
1548
if( (delta-2*Q0*E1( 1)).lt.0.d0 )ifail( 1)=1
1549
if( (delta-2*Q0*E1(-1)).lt.0.d0 )ifail(-1)=1
1551
if(ifail(1).eq.1.and.ifail(-1).eq.1)return
1555
if(E1(i).lt.xm1)then
1559
p1(i)=E1(i)**2-xm1**2
1560
if(p1(i).lt.0.d0)then
1567
if(E2(i).lt.xm2)then
1571
p2(i)=E2(i)**2-xm2**2
1572
if(p2(i).lt.0.d0)then
1584
subroutine put_on_MC_mshell_Hevin(p,xmi,xm1,xm2,mfail)
1586
double precision p(0:3,99),xmi,xm1,xm2
1589
include 'nexternal.inc'
1590
double precision chy,shy,chymo,shatp,q0,stot,xm1_r,xm2_r,
1592
double precision p1(0:3),p2(0:3),xk(0:3)
1596
common/fks_indices/i_fks,j_fks
1598
double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
1599
common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
1602
double precision xmcmass(nexternal)
1603
common/cxmcmass/xmcmass
1605
double precision zaxis(1:3)
1608
double precision rho
1611
if(j_fks.gt.nincoming)then
1612
write(*,*)'put_on_MC_mshell_Hevin must not be called for FSR'
1616
stot=4d0*ebeam(1)*ebeam(2)
1617
chy=cosh(ybst_til_tocm)
1618
shy=sinh(ybst_til_tocm)
1625
c write(*,*) 'p1:',(p1(i),i=0,3)
1626
c write(*,*) 'p2:',(p2(i),i=0,3)
1627
c write(*,*) 'xk:',(xk(i),i=0,3)
1628
call boostwdir2(chy,shy,chymo,zaxis,p1,p1)
1629
call boostwdir2(chy,shy,chymo,zaxis,p2,p2)
1630
call boostwdir2(chy,shy,chymo,zaxis,xk,xk)
1631
if(abs(p1(3)+p2(3)).gt.1.d-6)then
1632
write(*,*)'Error #1 in put_on_MC_mshell_Hevin',p1(3),p2(3)
1633
write(*,*) 'p1:',(p1(i),i=0,3)
1634
write(*,*) 'p2:',(p2(i),i=0,3)
1635
write(*,*) 'xk:',(xk(i),i=0,3)
1638
q0=p1(0)+p2(0)-xk(0)
1639
xk(0)=sqrt(xmi**2+rho(xk)**2)
1641
if(shatp.lt.shat*0.9999d0)then
1642
write(*,*)'Error #2 in put_on_MC_mshell_Hevin',shat,shatp
1644
elseif(shatp.le.shat)then
1647
if(shatp.ge.stot)then
1651
ybst_cm_tolab=ybst_til_tolab-ybst_til_tocm
1654
call getxmss_madfks(shatp,ybst_cm_tolab,
1655
# p1(3),xm1_r,p2(3),xm2_r,
1656
# p1(3),p2(3),mfail)
1657
if(mfail.eq.1)return
1658
p1(0)=sqrt(xm1_r**2+p1(3)**2)
1659
p2(0)=sqrt(xm2_r**2+p2(3)**2)
1662
call boostwdir2(chy,-shy,chymo,zaxis,p1,p1)
1663
call boostwdir2(chy,-shy,chymo,zaxis,p2,p2)
1664
call boostwdir2(chy,-shy,chymo,zaxis,xk,xk)
1675
subroutine put_on_MG_mshell()
1676
c Puts particles back on the MadGraph mass shell
1678
include 'nexternal.inc'
1680
double precision zero
1681
parameter (zero=0.d0)
1682
c Masses to be given to genps_fks.f
1683
double precision emass(nexternal)
1684
common/to_mass/ emass
1685
c Masses of the real process, as set by MadGraph
1687
double precision pmass(nexternal)
1698
subroutine write_masses_lhe_MG()
1699
c Set masses used by MC equal to MG ones
1701
include 'nexternal.inc'
1702
double precision xmcmass(nexternal)
1703
common/cxmcmass/xmcmass
1706
double precision zero
1707
parameter (ZERO = 0d0)
1709
double precision pmass(nexternal)
1720
subroutine get3space(pin,xlength,xdir)
1722
real*8 pin(0:3),xlength,xdir(1:3)
1725
xlength=pin(1)**2+pin(2)**2+pin(3)**2
1726
if(xlength.eq.0)then
1731
xlength=sqrt(xlength)
1733
xdir(i)=pin(i)/xlength
1740
subroutine put_on_MC_mshell_Hevout_old(p,xmi,xmj)
1742
double precision p(0:3,99),xmi,xmj
1743
double precision dirbst(1:3),dirpj(1:3)
1744
double precision pipj(0:3),pibst(0:3),pjbst(0:3)
1747
double precision dot,Q2,Q,expy,chybst,shybst,chybstmo,tmp,
1748
# Emo,xpmo,pij,xmi2,xmj2
1751
include 'nexternal.inc'
1754
common/fks_indices/i_fks,j_fks
1757
pipj(i)=p(i,i_fks)+p(i,j_fks)
1761
if(Q2.lt.(xmi+xmj)**2)then
1762
write(*,*)'Error in put_on_MC_mshell_Hevout',Emo,xmi,xmj
1766
call get3space(pipj,xpmo,dirbst)
1767
expy=sqrt((Emo+xpmo)/(Emo-xpmo))
1768
chybst=0.5d0*(expy+1.d0/expy)
1769
shybst=0.5d0*(expy-1.d0/expy)
1770
chybstmo=chybst-1.d0
1771
call boostwdir2(chybst,shybst,chybstmo,dirbst,p(0,j_fks),pjbst)
1772
call boostwdir2(chybst,shybst,chybstmo,dirbst,p(0,i_fks),pibst)
1773
call get3space(pjbst,tmp,dirpj)
1775
write(*,*)'Parton j soft in put_on_MC_mshell_Hevout'
1781
pjbst(0)=Q/2.d0*(1+(xmj2-xmi2)/Q2)
1782
pibst(0)=Q/2.d0*(1-(xmj2-xmi2)/Q2)
1783
pij=Q2**2-2*Q2*(xmi2+xmj2)+(xmi2-xmj2)**2
1784
pij=1/2.d0*sqrt( pij/Q2 )
1786
pjbst(i)= dirpj(i) * pij
1787
pibst(i)=-dirpj(i) * pij
1790
call boostwdir2(chybst,-shybst,chybstmo,dirbst,pibst,p(0,i_fks))
1791
call boostwdir2(chybst,-shybst,chybstmo,dirbst,pjbst,p(0,j_fks))
1793
if(j_fks.le.nincoming)then
1794
write(*,*)'ISR not yet implemented in put_on_MC_mshell_Hevout'
added
removed
Template/NLO/SubProcesses/add_write_info.f
