1
subroutine addmothers(ip,jpart,pb,isym,jsym,rscale,aqcd,aqed,buff,npart)
5
include 'nexternal.inc'
12
integer jpart(7,-nexternal+3:2*nexternal-3),npart,ip
13
double precision pb(0:4,-nexternal+3:2*nexternal-3)
14
double precision rscale,aqcd,aqed,targetamp(maxamps)
17
integer isym(nexternal,99), jsym
18
integer i,j,k,ida(2),ns,nres,ires,icl,ito2,idenpart,nc,ic
19
integer ito(-nexternal+3:nexternal),iseed
20
integer icolalt(2,-nexternal+3:2*nexternal-3)
21
double precision qicl(-nexternal+3:2*nexternal-3), factpm
22
double precision xtarget
27
double precision pmass(-nexternal:0,lmaxconfigs)
28
double precision pwidth(-nexternal:0,lmaxconfigs)
29
integer pow(-nexternal:0,lmaxconfigs)
32
data first_time /.true./
34
Double Precision amp2(maxamps), jamp2(0:maxflow)
35
common/to_amps/ amp2, jamp2
37
integer mincfig, maxcfig
38
common/to_configs/mincfig, maxcfig
39
integer idmap(-nexternal:nexternal),icmp
41
integer iforest(2,-max_branch:-1,lmaxconfigs)
42
common/to_forest/ iforest
43
integer sprop(-max_branch:-1,lmaxconfigs)
44
integer tprid(-max_branch:-1,lmaxconfigs)
45
common/to_sprop/sprop,tprid
46
integer mapconfig(0:lmaxconfigs), this_config
47
common/to_mconfigs/mapconfig, this_config
48
integer iconfig,igraph
50
integer idup(nexternal,maxproc)
51
integer mothup(2,nexternal,maxproc)
52
integer icolup(2,nexternal,maxflow)
53
include 'leshouche.inc'
55
include 'coloramps.inc'
57
logical OnBW(-nexternal:0) !Set if event is on B.W.
58
common/to_BWEvents/ OnBW
60
c integer ncols,ncolflow(maxamps),ncolalt(maxamps),icorg
61
c common/to_colstats/ncols,ncolflow,ncolalt,icorg
63
double precision pt,ran1
74
c Choose the config (diagram) which was actually used to produce the event
77
if(maxcfig-mincfig.ne.0) then
78
write(*,*) 'Error! Cannot cmaxcfig != mincfig: ',maxcfig,mincfig
82
c ...unless the diagram is passed in igscl(1); then use that diagram
83
c if (igscl(0).ne.0) then
84
c if (btest(mlevel,3)) then
85
c write(*,*)'unwgt.f: write out diagram ',igscl(1)
91
c Choose a color flow which is certain to work with the propagator
92
c structure of the chosen diagram and use that as an alternative
94
igraph=mapconfig(iconfig)
98
if(icolamp(igraph,1)) then
100
c print *,'Color flow 1 allowed for graph ',igraph
105
if(icolamp(igraph,ic))then
106
targetamp(ic) = jamp2(ic)+targetamp(ic-1)
107
c print *,'Color flow ',ic,' allowed for graph ',igraph
109
targetamp(ic)=targetamp(ic-1)
112
xtarget=ran1(iseed)*targetamp(nc)
115
do while (targetamp(ic) .lt. xtarget .and. ic .lt. nc)
118
c print *,'Chose color flow ',ic
120
if(targetamp(nc).eq.0) ic=0
123
icolalt(1,isym(i,jsym))=icolup(1,i,ic)
124
icolalt(2,isym(i,jsym))=icolup(2,i,ic)
126
icolalt(1,i)=jpart(4,i)
127
icolalt(2,i)=jpart(5,i)
141
c Get mother information from chosen graph
144
c First check number of resonant s-channel propagators
148
c Loop over propagators to find mother-daughter information
149
do i=-1,-nexternal+3,-1
151
ida(1)=iforest(1,i,iconfig)
152
ida(2)=iforest(2,i,iconfig)
154
if(ida(j).gt.0) ida(j)=isym(ida(j),jsym)
156
c Decide s- or t-channel
157
if(iabs(sprop(i,iconfig)).gt.0) then ! s-channel propagator
158
jpart(1,i)=sprop(i,iconfig)
161
c Don't care about t-channel propagators
164
c Set status codes for propagator
165
c if((igscl(0).ne.0.and.
166
c $ (iabs(jpart(1,i)).gt.5.and.iabs(jpart(1,i)).lt.11).or.
167
c $ (iabs(jpart(1,i)).gt.16.and.iabs(jpart(1,i)).ne.21)).or.
168
c $ (igscl(0).eq.0.and.OnBW(i))) then
170
c Resonance whose mass should be preserved
174
c Propagator for documentation only - not included
177
c Calculate momentum (p1+p2 for s-channel, p2-p1 for t-channel)
179
pb(j,i)=pb(j,ida(1))+pb(j,ida(2))
181
pb(4,i)=sqrt(max(0d0,pb(0,i)**2-pb(1,i)**2-pb(2,i)**2-pb(3,i)**2))
182
c if(jpart(6,i).eq.2.and.
183
c $ abs(pb(4,i)-pmass(i,iconfig)).gt.5d0*pwidth(i,iconfig)) then
187
c Remove propagator if "decay" to same particle (i.e. radiation of g/gamma/higgs)
188
if(jpart(6,i).eq.2 .and. (jpart(1,i).eq.jpart(1,ida(1)).or.
189
$ jpart(1,i).eq.jpart(1,ida(2)))) then
190
if(jpart(1,i).eq.jpart(1,ida(1))) idenpart=ida(1)
191
if(jpart(1,i).eq.jpart(1,ida(2))) idenpart=ida(2)
192
c Always remove if daughter final-state or already removed for this reason
193
if(idenpart.gt.0.or.jpart(6,idenpart).eq.4) then
196
c Else remove either this resonance or daughter, which is closer to mass shell
197
elseif(abs(pb(4,i)-pmass(i,iconfig)).gt.
198
$ abs(pb(4,idenpart)-pmass(i,iconfig))) then
201
else if(jpart(6,idenpart).eq.2) then
206
c Set color info for all s-channels
207
c Fist set "safe" color info
208
if(icolalt(1,ida(1))+icolalt(1,ida(2))-
209
$ icolalt(2,ida(1))-icolalt(2,ida(2)).eq.0) then ! color singlet
212
elseif(icolalt(1,ida(1))-icolalt(2,ida(2)).eq.0) then
213
icolalt(1,i) = icolalt(1,ida(2))
214
icolalt(2,i) = icolalt(2,ida(1))
215
else if(icolalt(1,ida(2))-icolalt(2,ida(1)).eq.0) then
216
icolalt(1,i) = icolalt(1,ida(1))
217
icolalt(2,i) = icolalt(2,ida(2))
218
else if(jpart(6,i).ge.3) then ! Don't need to match
219
icolalt(1,i) = icolalt(1,ida(1))+icolalt(1,ida(2))
220
icolalt(2,i) = icolalt(2,ida(1))+icolalt(2,ida(2))
222
c Erraneous color assignment for propagator
225
c Set tentative mothers
228
c Set mother info for daughters
233
c Just zero helicity info for intermediate states
236
goto 100 ! only s-channel propagators found
238
c Propagator not compatible with color flow -> remove the propagator
239
print *,'ERROR: Safe color assignment wrong! This should never happen!'
241
if(jpart(6,i).eq.2) then
246
c print *,'Diagram: ',igraph,', Color flow: ',ic
247
c print *,'colors:',icolalt(1,ida(1)),icolalt(2,ida(1)),
248
c $ icolalt(1,ida(2)),icolalt(2,ida(2))
249
c write(*,'(15i5)'),(jpart(1,j),j=1,nexternal),(jpart(1,j),j=-1,i,-1)
250
c write(*,'(15i5)'),(jpart(2,j),j=1,nexternal),(jpart(2,j),j=-1,i,-1)
251
c write(*,'(15i5)'),(jpart(3,j),j=1,nexternal),(jpart(3,j),j=-1,i,-1)
252
c write(*,'(15i5)'),(icolalt(1,j),j=1,nexternal),(icolalt(1,j),j=-1,i,-1)
253
c write(*,'(15i5)'),(icolalt(2,j),j=1,nexternal),(icolalt(2,j),j=-1,i,-1)
254
c write(*,'(15i5)'),(jpart(6,j),j=1,nexternal),(jpart(6,j),j=-1,i,-1)
255
c write(*,'(15f5.0)'),(pb(4,j),j=1,nexternal),(pb(4,j),j=-1,i,-1)
259
c Remove non-resonant mothers, set position of particles
262
jpart(4,i)=icolalt(1,i)
263
jpart(5,i)=icolalt(2,i)
264
if(i.eq.1.or.i.eq.2) then
265
ito(i)=i ! initial state particle
267
ito(i)=i+nres ! final state particle
268
else if(i.le.-1.and.jpart(6,i).eq.2) then
270
ito(i)=2+ires ! s-channel resonances
275
if(jpart(2,i).lt.0.and.jpart(6,jpart(2,i)).ne.2) then
276
jpart(2,i)=jpart(2,jpart(2,i))
277
jpart(3,i)=jpart(3,jpart(3,i))
282
c Shift particles to right place and set mothers of particles
284
do i=nexternal,-ns,-1
285
if(ito(i).le.0) cycle
287
jpart(j,ito(i))=jpart(j,i)
289
if(jpart(2,ito(i)).lt.0) then
290
jpart(2,ito(i))=ito(jpart(2,ito(i)))
291
jpart(3,ito(i))=ito(jpart(3,ito(i)))
299
$ write(buff,'(a1,9e15.7)') '#',(ptclus(i),i=3,min(nexternal,11))
300
npart = nexternal+nres
added
removed
Template/SubProcesses/addmothers.f
