← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/Source/invarients.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/Source/invarients.f
 
1
      subroutine set_invarients(nfinal,ninvar)
 
2
c***************************************************************************
 
3
c     Calculates all of the invarients for a 2->n process
 
4
c***************************************************************************
 
5
      implicit none
 
6
c     
 
7
c     Constants
 
8
c
 
9
      include 'genps.inc'
 
10
c
 
11
c     Arguments
 
12
c
 
13
      integer nfinal,ninvar
 
14
c
 
15
c     Local
 
16
c     
 
17
      integer ip1,ip2,ipstart,ipstop,np,i
 
18
      integer ncycle
 
19
      character*10 buff
 
20
c
 
21
c     Global
 
22
c
 
23
      integer              imom(maxinvar),ninvarients
 
24
      common/to_invarients/imom          ,ninvarients
 
25
c-----
 
26
c  Begin Code
 
27
c-----
 
28
 
 
29
      do i=1,nfinal
 
30
         imom(i)=i
 
31
      enddo
 
32
      ipstart=1
 
33
      ipstop =nfinal
 
34
      np     =nfinal
 
35
c
 
36
c     First do all the s-channel
 
37
c
 
38
      do ncycle=2,nfinal-1
 
39
         do ip1 = ipstart,ipstop-1
 
40
            do ip2=int((real(imom(ip1))/10.-imom(ip1)/10)*10+.1)+1,
 
41
     $           nfinal
 
42
               np=np+1
 
43
               if (np .gt. maxinvar) then
 
44
                  print*,'Sorry too many invarients',np,ip1,ip2,ncycle
 
45
                  stop
 
46
               endif
 
47
               imom(np)=imom(ip1)*10+imom(ip2)
 
48
               if (imom(np) .lt. 10) then
 
49
                  write(buff,'(a2,i1)') 'S?',imom(np)
 
50
               elseif (imom(np) .lt. 100) then
 
51
                  write(buff,'(a2,i2)') 'S?',imom(np)
 
52
               elseif (imom(np) .lt. 1000) then
 
53
                  write(buff,'(a2,i3)') 'S?',imom(np)
 
54
               elseif (imom(np) .lt. 10000) then
 
55
                  write(buff,'(a2,i4)') 'S?',imom(np)
 
56
               elseif (imom(np) .lt. 100000) then
 
57
                  write(buff,'(a2,i5)') 'S?',imom(np)
 
58
               else
 
59
                  write(buff,'(a2,i6)') 'S?',imom(ip1)
 
60
               endif
 
61
c               call hbook1(100+np-nfinal,buff,100,0.,1.,0.)
 
62
c               write(*,'(i4,i6)') np-nfinal,imom(np)
 
63
               write(*,'(i4,a1,a6)') np-nfinal,'=',buff
 
64
               if ((np-nfinal)/7 .eq. real(np-nfinal)/7.) write(*,*)' '
 
65
            enddo
 
66
         enddo
 
67
         ipstart=ipstop+1
 
68
         ipstop = np
 
69
      enddo
 
70
c
 
71
c     Now do the t-channel
 
72
c
 
73
      ipstop = np
 
74
      do ip1 = 1,ipstop
 
75
c         write(*,'(i4,a2,i6)') np-nfinal+ip1,'a-',imom(ip1)
 
76
         if (imom(ip1) .lt. 10) then
 
77
            write(buff,'(a2,i1)') 'T?',imom(ip1)
 
78
         elseif (imom(ip1) .lt. 100) then
 
79
            write(buff,'(a2,i2)') 'T?',imom(ip1)
 
80
         elseif (imom(ip1) .lt. 1000) then
 
81
            write(buff,'(a2,i3)') 'T?',imom(ip1)
 
82
         elseif (imom(ip1) .lt. 10000) then
 
83
            write(buff,'(a2,i4)') 'T?',imom(ip1)
 
84
         elseif (imom(ip1) .lt. 100000) then
 
85
            write(buff,'(a2,i5)') 'T?',imom(ip1)
 
86
         else
 
87
            write(buff,'(a2,i6)') 'T?',imom(ip1)
 
88
         endif
 
89
c         call hbook1(100+np-nfinal+ip1,buff,100,0.,1.,0.)
 
90
c         write(*,*) np-nfinal+ip1,buff
 
91
         write(*,'(i4,a1,a6)') np-nfinal+ip1,'=',buff
 
92
         if ((np-nfinal+ip1)/7 .eq. real(np-nfinal+ip1)/7.) write(*,*)
 
93
      enddo
 
94
      write(*,*)
 
95
      print*,'Particles, Invarients',nfinal,np-nfinal+np
 
96
      ninvarients=np-nfinal+np
 
97
      ninvar=ninvarients
 
98
      if (ninvarients .gt. maxinvar) then
 
99
         print*,'Error too many invarients to map'
 
100
c         stop
 
101
      endif
 
102
      end
 
103
 
 
104
 
 
105
      subroutine fill_invarients(nfinal,p1,s,xx)
 
106
c***************************************************************************
 
107
c     Calculates all of the invarients for a 2->n process
 
108
c***************************************************************************
 
109
      implicit none
 
110
c     
 
111
c     Constants
 
112
c     
 
113
      include 'genps.inc'
 
114
c
 
115
c     Arguments
 
116
c
 
117
      integer nfinal
 
118
      double precision p1(0:3,nfinal+2),s,xx(55)
 
119
c
 
120
c     Local
 
121
c     
 
122
      integer ip1,ip2,ipstart,ipstop,np,i,j
 
123
      integer imom(maxinvar)
 
124
      integer ncycle
 
125
      character*10 buff
 
126
      double precision p(0:3,maxinvar),ptemp(0:3)
 
127
c
 
128
c     External
 
129
c
 
130
      double precision dot
 
131
      external         dot
 
132
c-----
 
133
c  Begin Code
 
134
c-----
 
135
 
 
136
      do i=1,nfinal
 
137
         imom(i) = i
 
138
         do j=0,3
 
139
            p(j,i)=p1(j,i+2)
 
140
         enddo
 
141
c         write(*,'(i3,4f17.8)') i,(p(j,i),j=0,3)
 
142
      enddo
 
143
      ipstart=1
 
144
      ipstop =nfinal
 
145
      np     =nfinal
 
146
c
 
147
c     First do all the s-channel
 
148
c
 
149
      do ncycle=2,nfinal-1
 
150
         do ip1 = ipstart,ipstop-1
 
151
            do ip2=int((real(imom(ip1))/10.-imom(ip1)/10)*10+.1)+1
 
152
     $           ,nfinal
 
153
               np=np+1
 
154
               if (np .gt. maxinvar) then
 
155
                  print*,'Sorry too many invarients',np,ip1,ip2,ncycle
 
156
                  stop
 
157
               endif
 
158
               imom(np)=imom(ip1)*10+imom(ip2)
 
159
               do j=0,3
 
160
                  p(j,np) = p(j,ip1)+p(j,ip2)
 
161
               enddo
 
162
               xx(np-nfinal) = dot(p(0,np),p(0,np))/s
 
163
c               call hfill(100+np-nfinal,
 
164
c     &              real(dot(p(0,np),p(0,np))/s),0.,wgt)
 
165
c               write(*,'(i4,3f20.8)') np-nfinal,
 
166
c     &              real(dot(p(0,np),p(0,np))/s)
 
167
            enddo
 
168
         enddo
 
169
         ipstart=ipstop+1
 
170
         ipstop = np
 
171
      enddo
 
172
c
 
173
c     Now do the t-channel
 
174
c
 
175
      ipstop = np
 
176
      do ip1 = 1,ipstop
 
177
         do j = 0,3
 
178
            ptemp(j)=p1(j,1)-p(j,ip1)
 
179
         enddo
 
180
         xx(np-nfinal+ip1)= .5d0*(dot(ptemp,ptemp)/s+1d0)
 
181
c         call hfill(100+np-nfinal+ip1,real(-dot(ptemp,ptemp)/s),0.,wgt)
 
182
c         write(*,'(i4,3f20.8)') np-nfinal+ip1,
 
183
c     &         real(-dot(ptemp,ptemp)/s)
 
184
      enddo
 
185
      end
 
186
 
 
187
 
 
188
      subroutine map_invarients(Minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
 
189
c****************************************************************************
 
190
c     Determines mappings for each structure of invarients onto integration
 
191
c     variables.  Input: Ninvar, iforest.  Output: Minvar, ninvar
 
192
c****************************************************************************
 
193
      implicit none
 
194
c
 
195
c     Constants
 
196
c
 
197
      include 'genps.inc'
 
198
c
 
199
c     Arguments
 
200
c
 
201
      integer Minvar(maxdim,lmaxconfigs),nconfigs,ninvar,nexternal,nincoming
 
202
      integer mincfig,maxcfig
 
203
c
 
204
c     Local
 
205
c
 
206
      integer iconfig, jgrid,j, nbranch
 
207
      logical found,tchannel
 
208
c
 
209
c     Global
 
210
c
 
211
      integer              imom(maxinvar),ninvarients
 
212
      common/to_invarients/imom          ,ninvarients
 
213
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
214
      common/to_forest/ iforest
 
215
 
 
216
c-----
 
217
c  Begin Code
 
218
c----
 
219
 
 
220
      nbranch = nexternal-2
 
221
      jgrid=0
 
222
c
 
223
c     
 
224
c     Try simple mapping if nconfigs = 1
 
225
c
 
226
 
 
227
      if (nconfigs .eq. 1) then
 
228
         do j=1,3*nbranch-4+2
 
229
            minvar(j,mincfig)=j
 
230
         enddo
 
231
         jgrid=j-1
 
232
      else
 
233
c      if (ep) jgrid=1
 
234
c      if (pp) jgrid=2
 
235
      do iconfig=mincfig,maxcfig
 
236
         tchannel = .false.         
 
237
         do j=1,nbranch-1
 
238
            if (iforest(1,-j,iconfig) .eq. 1) then
 
239
               tchannel=.true.
 
240
            endif
 
241
            jgrid=jgrid+1
 
242
            minvar(j,iconfig) = jgrid
 
243
            if (tchannel .and. j .lt. nbranch-1) then
 
244
               jgrid=jgrid+1            
 
245
               minvar(nbranch-1+2*j,iconfig)=jgrid
 
246
            endif
 
247
         enddo             !Each Branch
 
248
         if (.not. tchannel .and. nincoming.eq.2) then          !Don't need last s-channel
 
249
            jgrid=jgrid-1
 
250
            minvar(nbranch-1,iconfig)=0
 
251
         endif
 
252
c         if (pp) then
 
253
c            jgrid=jgrid+1
 
254
c            minvar(3*nbranch-3,iconfig)=jgrid
 
255
c            jgrid=jgrid+1
 
256
c            minvar(3*nbranch-2,iconfig)=jgrid
 
257
c         elseif (ep) then
 
258
c            jgrid=jgrid+1
 
259
c            minvar(3*nbranch-3,iconfig)=jgrid
 
260
c         endif
 
261
      enddo  !Each configurations
 
262
      endif
 
263
      ninvar = jgrid
 
264
      end
 
265
 
 
266
      subroutine sortint(n,ra)
 
267
      integer ra(n)
 
268
      l=n/2+1
 
269
      ir=n
 
270
10    continue
 
271
        if(l.gt.1)then
 
272
          l=l-1
 
273
          rra=ra(l)
 
274
        else
 
275
          rra=ra(ir)
 
276
          ra(ir)=ra(1)
 
277
          ir=ir-1
 
278
          if(ir.eq.1)then
 
279
            ra(1)=rra
 
280
            return
 
281
          endif
 
282
        endif
 
283
        i=l
 
284
        j=l+l
 
285
20      if(j.le.ir)then
 
286
          if(j.lt.ir)then
 
287
            if(ra(j).lt.ra(j+1))j=j+1
 
288
          endif
 
289
          if(rra.lt.ra(j))then
 
290
            ra(i)=ra(j)
 
291
            i=j
 
292
            j=j+j
 
293
          else
 
294
            j=ir+1
 
295
          endif
 
296
        go to 20
 
297
        endif
 
298
        ra(i)=rra
 
299
      go to 10
 
300
      end
 
301
 
 
302