← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_242195/PROC_242195/Source/invarients.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_242195/PROC_242195/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
      include 'maxconfigs.inc'
 
199
c
 
200
c     Arguments
 
201
c
 
202
      integer Minvar(maxdim,lmaxconfigs),nconfigs,ninvar,nexternal,nincoming
 
203
      integer mincfig,maxcfig
 
204
c
 
205
c     Local
 
206
c
 
207
      integer iconfig, jgrid,j, nbranch
 
208
      logical found,tchannel
 
209
c
 
210
c     Global
 
211
c
 
212
      integer              imom(maxinvar),ninvarients
 
213
      common/to_invarients/imom          ,ninvarients
 
214
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
215
      common/to_forest/ iforest
 
216
 
 
217
c-----
 
218
c  Begin Code
 
219
c----
 
220
 
 
221
      nbranch = nexternal-2
 
222
      jgrid=0
 
223
c
 
224
c     
 
225
c     Try simple mapping if nconfigs = 1
 
226
c
 
227
 
 
228
      if (nconfigs .eq. 1) then
 
229
c         do j=1,3*nbranch-4+2
 
230
         do j=1,maxdim
 
231
            minvar(j,mincfig)=j
 
232
         enddo
 
233
         jgrid=j-1
 
234
      else
 
235
c      if (ep) jgrid=1
 
236
c      if (pp) jgrid=2
 
237
      do iconfig=mincfig,maxcfig
 
238
         tchannel = .false.         
 
239
         do j=1,nbranch-1
 
240
            if (iforest(1,-j,iconfig) .eq. 1) then
 
241
               tchannel=.true.
 
242
            endif
 
243
            jgrid=jgrid+1
 
244
            minvar(j,iconfig) = jgrid
 
245
            if (tchannel .and. j .lt. nbranch-1) then
 
246
               jgrid=jgrid+1            
 
247
               minvar(nbranch-1+2*j,iconfig)=jgrid
 
248
            endif
 
249
         enddo             !Each Branch
 
250
         if (.not. tchannel .and. nincoming.eq.2) then          !Don't need last s-channel
 
251
            jgrid=jgrid-1
 
252
            minvar(nbranch-1,iconfig)=0
 
253
         endif
 
254
c         if (pp) then
 
255
c            jgrid=jgrid+1
 
256
c            minvar(3*nbranch-3,iconfig)=jgrid
 
257
c            jgrid=jgrid+1
 
258
c            minvar(3*nbranch-2,iconfig)=jgrid
 
259
c         elseif (ep) then
 
260
c            jgrid=jgrid+1
 
261
c            minvar(3*nbranch-3,iconfig)=jgrid
 
262
c         endif
 
263
      enddo  !Each configurations
 
264
      endif
 
265
      ninvar = jgrid
 
266
      end
 
267
 
 
268
      subroutine sortint(n,ra)
 
269
      integer ra(n)
 
270
      l=n/2+1
 
271
      ir=n
 
272
10    continue
 
273
        if(l.gt.1)then
 
274
          l=l-1
 
275
          rra=ra(l)
 
276
        else
 
277
          rra=ra(ir)
 
278
          ra(ir)=ra(1)
 
279
          ir=ir-1
 
280
          if(ir.eq.1)then
 
281
            ra(1)=rra
 
282
            return
 
283
          endif
 
284
        endif
 
285
        i=l
 
286
        j=l+l
 
287
20      if(j.le.ir)then
 
288
          if(j.lt.ir)then
 
289
            if(ra(j).lt.ra(j+1))j=j+1
 
290
          endif
 
291
          if(rra.lt.ra(j))then
 
292
            ra(i)=ra(j)
 
293
            i=j
 
294
            j=j+j
 
295
          else
 
296
            j=ir+1
 
297
          endif
 
298
        go to 20
 
299
        endif
 
300
        ra(i)=rra
 
301
      go to 10
 
302
      end
 
303
 
 
304