← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2-alpha0

« back to all changes in this revision

Viewing changes to Template/SubProcesses/check_intdip.f

Added Template and HELAS into bzr

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/SubProcesses/check_intdip.f
 
1
      PROGRAM checkintdip
 
2
C**************************************************************************
 
3
C     This is a program for testing the integrated dipoles. For a couple of
 
4
C     phase space points, chosen by RAMBO, it prints the the singular and
 
5
C     finite terms for the given process.
 
6
C     R.Frederix and N.Greiner - February 2010
 
7
C**************************************************************************
 
8
      IMPLICIT NONE
 
9
C     
 
10
C     CONSTANTS  
 
11
C     
 
12
      REAL*8 ZERO
 
13
      PARAMETER (ZERO=0D0)
 
14
C     
 
15
C     INCLUDE FILES
 
16
C     
 
17
C---  the include file with the values of the parameters and masses     
 
18
      INCLUDE "coupl.inc"
 
19
C---  integer nexternal ! number particles (incoming+outgoing) in the me 
 
20
      INCLUDE "nexternal.inc" 
 
21
C---  particle masses
 
22
      REAL*8 PMASS(NEXTERNAL)   
 
23
C---  integer    n_max_cg
 
24
      INCLUDE "ngraphs.inc"     !how many diagrams (could be useful to know...)
 
25
      INCLUDE "dipole.inc" 
 
26
 
 
27
C     
 
28
C     LOCAL
 
29
C     
 
30
      INTEGER I,J,K,l,m,l1,m1, number
 
31
      REAL*8 P(0:3,NEXTERNAL)      ! four momenta. Energy is the zeroth component.
 
32
      REAL*8 SQRTS                 ! sqrt(s)= center of mass energy 
 
33
      REAL*8 ME,SUBTRACT,SUBTRACT2 ! Matrix Element and subtraction terms
 
34
      real*8 start                 ! Starting value for invariants and the number of phase space points
 
35
      logical cut                  ! Cut to prevent going into double limits
 
36
      real*8 sum,x(2)
 
37
      real*8 Z, eps,epssq, finite
 
38
C     
 
39
C     EXTERNAL
 
40
C     
 
41
      REAL*8 DOT
 
42
      EXTERNAL DOT
 
43
      
 
44
C-----
 
45
C     BEGIN CODE
 
46
C-----
 
47
C     
 
48
C---  INITIALIZATION CALLS
 
49
C     
 
50
c---  Call to initialize the values of the couplings, masses and widths 
 
51
c     used in the evaluation of the matrix element. The primary parameters of the
 
52
c     models are read from Cards/param_card.dat. The secondary parameters are calculated
 
53
c     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
 
54
c     in coupl.inc .
 
55
 
 
56
      call setpara('param_card.dat',.true.)  !first call to setup the paramaters
 
57
      include "pmass.inc"                    !set up masses
 
58
      mu  = 91.188d0  !renormalization scale
 
59
      muf = 91.188d0  !factorization scale
 
60
 
 
61
      SQRTS=1000d0              !CMS energy in GEV
 
62
      x(1)=1d0                  !Bjorken x's
 
63
      x(2)=1d0
 
64
      number=3            !Number of phase space points
 
65
 
 
66
c For this check routine to work properly,
 
67
c the last parton should be a massless one
 
68
      if (pmass(nexternal).ne.0d0) then
 
69
         write (*,*) 'Error in process order:'
 
70
         write (*,*) 'Please, change the order of the particles in the'
 
71
         write (*,*)
 
72
     &        'process such that the last one is a massless particle'
 
73
         stop
 
74
      endif
 
75
 
 
76
c     Loop over phase space points
 
77
      do i=1,number
 
78
 
 
79
c---  Now use a simple multipurpose PS generator (RAMBO) just to get a 
 
80
c     RANDOM set of four momenta of given masses pmass(i) to be used to evaluate 
 
81
c     the integrated dipoles.
 
82
c       
 
83
         CALL GET_INT_MOMENTA(SQRTS,PMASS,P,Z)
 
84
 
 
85
         CALL INTDIPOLES(P,X,Z,1d0,EPS,EPSSQ,FINITE)
 
86
c Note that for the intdipolesfinite the phase space point from
 
87
c Rambo is, in general not correct
 
88
c         CALL INTDIPOLESFINITE(P,X,Z,1d0,EPS,EPSSQ,FINITE)
 
89
         write (*,*) ' '
 
90
         write (*,*) 'PHASE SPACE POINT',i,':'
 
91
         do j=1, nexternal-1
 
92
            write (*,*)'   ',j, P(0,j), P(1,j), P(2,j), P(3,j)
 
93
         enddo
 
94
         write (*,*) '    Z =', Z
 
95
         write (*,*) ' '
 
96
         write (*,*) '    1/eps^2: ',epssq
 
97
         write (*,*) '    1/eps:   ',eps
 
98
         write (*,*) '    Finite:  ',finite
 
99
         write (*,*) ' '
 
100
         write (*,*) ' '
 
101
      enddo
 
102
 
 
103
      end
 
104
 
 
105
 
 
106
 
 
107
      SUBROUTINE GET_INT_MOMENTA(ENERGY,PMASS,P,Z)
 
108
C     Provide four-momenta for integrated dipoles.
 
109
 
 
110
      IMPLICIT NONE
 
111
      INCLUDE "nexternal.inc"
 
112
C     ARGUMENTS
 
113
      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT,Z
 
114
C     LOCAL
 
115
      INTEGER I
 
116
      
 
117
      P(0,1)=energy/2
 
118
      P(1,1)=0d0
 
119
      P(2,1)=0d0
 
120
      P(3,1)=energy/2
 
121
      
 
122
      P(0,2)=energy/2
 
123
      P(1,2)=0d0
 
124
      P(2,2)=0d0
 
125
      P(3,2)=-energy/2
 
126
 
 
127
      call rambo_int(nexternal-3,energy,pmass(3),prambo,Z)
 
128
      DO I=3, NEXTERNAL-1
 
129
         P(0,I)=PRAMBO(4,I-2)   
 
130
         P(1,I)=PRAMBO(1,I-2)
 
131
         P(2,I)=PRAMBO(2,I-2)
 
132
         P(3,I)=PRAMBO(3,I-2)   
 
133
      ENDDO
 
134
 
 
135
      P(0,nexternal)=0d0
 
136
      P(1,nexternal)=0d0
 
137
      P(2,nexternal)=0d0
 
138
      P(3,nexternal)=0d0
 
139
 
 
140
      RETURN
 
141
      END
 
142
 
 
143
 
 
144
 
 
145
 
 
146
      double precision function dot(p1,p2)
 
147
C****************************************************************************
 
148
C     4-Vector Dot product
 
149
C****************************************************************************
 
150
      implicit none
 
151
      double precision p1(0:3),p2(0:3)
 
152
      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
 
153
      end
 
154
 
 
155
 
 
156
      SUBROUTINE RAMBO_INT(N,ET,XM,P,Zvar)
 
157
C     Modified RAMBO suited for integrated dipoles
 
158
C     R.Frederix and N.Greiner - February 2010
 
159
***********************************************************************
 
160
*                       RAMBO                                         *
 
161
*    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                  *
 
162
*                                                                     *
 
163
*    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR                *
 
164
*    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                 *
 
165
*    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                      *
 
166
*    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)    *
 
167
*                                                                     *
 
168
*    N  = NUMBER OF PARTICLES                                         *
 
169
*    ET = TOTAL CENTRE-OF-MASS ENERGY                                 *
 
170
*    XM = PARTICLE MASSES ( DIM=NEXTERNAL-2 )                         *
 
171
*    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-2) )                    *
 
172
*    WT = WEIGHT OF THE EVENT                                         *
 
173
***********************************************************************
 
174
      IMPLICIT REAL*8(A-H,O-Z)
 
175
          INCLUDE "nexternal.inc"
 
176
      DIMENSION XM(NEXTERNAL-2),P(4,NEXTERNAL-2)
 
177
      DIMENSION Q(4,NEXTERNAL-2),Z(NEXTERNAL-2),R(4),
 
178
     .   B(3),P2(NEXTERNAL-2),XM2(NEXTERNAL-2),
 
179
     .   E(NEXTERNAL-2),V(NEXTERNAL-2),IWARN(5)
 
180
      SAVE ACC,ITMAX,IBEGIN,IWARN
 
181
      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
 
182
*
 
183
* INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
 
184
      IF(IBEGIN.NE.0) GOTO 103
 
185
      IBEGIN=1
 
186
      TWOPI=8.*DATAN(1.D0)
 
187
      PO2LOG=LOG(TWOPI/4.)
 
188
      Z(2)=PO2LOG
 
189
      DO 101 K=3,10
 
190
  101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
 
191
      DO 102 K=3,10
 
192
  102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
 
193
*
 
194
* CHECK ON THE NUMBER OF PARTICLES
 
195
  103 IF(N.GT.1.AND.N.LT.101) GOTO 104
 
196
      PRINT 1001,N
 
197
      STOP
 
198
*
 
199
* CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
 
200
  104 XMT=0.
 
201
      NM=0
 
202
      DO 105 I=1,N
 
203
      IF(XM(I).NE.0.D0) NM=NM+1
 
204
  105 XMT=XMT+ABS(XM(I))
 
205
      IF(XMT.LE.ET) GOTO 201
 
206
      PRINT 1002,XMT,ET
 
207
      STOP
 
208
*
 
209
* THE PARAMETER VALUES ARE NOW ACCEPTED
 
210
*
 
211
* GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
 
212
  201 DO 202 I=1,N
 
213
         r1=rn(1)
 
214
      C=2.*r1-1.
 
215
      S=SQRT(1.-C*C)
 
216
      F=TWOPI*RN(2)
 
217
      r1=rn(3)
 
218
      r2=rn(4)
 
219
      Q(4,I)=-LOG(r1*r2)
 
220
      Q(3,I)=Q(4,I)*C
 
221
      Q(2,I)=Q(4,I)*S*COS(F)
 
222
  202 Q(1,I)=Q(4,I)*S*SIN(F)
 
223
 
 
224
      Zvar=rn(1)
 
225
*
 
226
* CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
 
227
      DO 203 I=1,4
 
228
  203 R(I)=0.
 
229
      DO 204 I=1,N
 
230
      DO 204 K=1,4
 
231
  204 R(K)=R(K)+Q(K,I)
 
232
      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
 
233
      DO 205 K=1,3
 
234
  205 B(K)=-R(K)/RMAS
 
235
      G=R(4)/RMAS
 
236
      A=1./(1.+G)
 
237
      X=ET/RMAS
 
238
*
 
239
* TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
 
240
      DO 207 I=1,N
 
241
      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
 
242
      DO 206 K=1,3
 
243
  206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
 
244
  207 P(4,I)=X*(G*Q(4,I)+BQ)
 
245
*
 
246
* CALCULATE WEIGHT AND POSSIBLE WARNINGS
 
247
      WT=PO2LOG
 
248
      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
 
249
      IF(WT.GE.-180.D0) GOTO 208
 
250
      IF(IWARN(1).LE.5) PRINT 1004,WT
 
251
      IWARN(1)=IWARN(1)+1
 
252
  208 IF(WT.LE. 174.D0) GOTO 209
 
253
      IF(IWARN(2).LE.5) PRINT 1005,WT
 
254
      IWARN(2)=IWARN(2)+1
 
255
*
 
256
* RETURN FOR WEIGHTED MASSLESS MOMENTA
 
257
  209 IF(NM.NE.0) GOTO 210
 
258
* RETURN LOG OF WEIGHT
 
259
      WT=WT
 
260
      RETURN
 
261
*
 
262
* MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
 
263
  210 XMAX=SQRT(1.-(XMT/ET)**2)
 
264
      DO 301 I=1,N
 
265
      XM2(I)=XM(I)**2
 
266
  301 P2(I)=P(4,I)**2
 
267
      ITER=0
 
268
      X=XMAX
 
269
      ACCU=ET*ACC
 
270
  302 F0=-ET
 
271
      G0=0.
 
272
      X2=X*X
 
273
      DO 303 I=1,N
 
274
      E(I)=SQRT(XM2(I)+X2*P2(I))
 
275
      F0=F0+E(I)
 
276
  303 G0=G0+P2(I)/E(I)
 
277
      IF(ABS(F0).LE.ACCU) GOTO 305
 
278
      ITER=ITER+1
 
279
      IF(ITER.LE.ITMAX) GOTO 304
 
280
      PRINT 1006,ITMAX
 
281
      GOTO 305
 
282
  304 X=X-F0/(X*G0)
 
283
      GOTO 302
 
284
  305 DO 307 I=1,N
 
285
      V(I)=X*P(4,I)
 
286
      DO 306 K=1,3
 
287
  306 P(K,I)=X*P(K,I)
 
288
  307 P(4,I)=E(I)
 
289
*
 
290
* CALCULATE THE MASS-EFFECT WEIGHT FACTOR
 
291
      WT2=1.
 
292
      WT3=0.
 
293
      DO 308 I=1,N
 
294
      WT2=WT2*V(I)/E(I)
 
295
  308 WT3=WT3+V(I)**2/E(I)
 
296
      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
 
297
*
 
298
* RETURN FOR  WEIGHTED MASSIVE MOMENTA
 
299
      WT=WT+WTM
 
300
      IF(WT.GE.-180.D0) GOTO 309
 
301
      IF(IWARN(3).LE.5) PRINT 1004,WT
 
302
      IWARN(3)=IWARN(3)+1
 
303
  309 IF(WT.LE. 174.D0) GOTO 310
 
304
      IF(IWARN(4).LE.5) PRINT 1005,WT
 
305
      IWARN(4)=IWARN(4)+1
 
306
* RETURN LOG OF WEIGHT
 
307
  310 WT=WT
 
308
      RETURN
 
309
*
 
310
 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
 
311
 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
 
312
     . ' SMALLER THAN TOTAL ENERGY =',D15.6)
 
313
 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
 
314
 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
 
315
 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
 
316
     . ' DESIRED ACCURACY =',D15.6)
 
317
      END
 
318
 
 
319
      FUNCTION RN(IDUMMY)
 
320
      REAL*8 RN,RAN
 
321
      SAVE INIT
 
322
      DATA INIT /1/
 
323
      IF (INIT.EQ.1) THEN
 
324
        INIT=0
 
325
        CALL RMARIN(1802,9373)
 
326
      END IF
 
327
*
 
328
  10  CALL RANMAR(RAN)
 
329
      IF (RAN.LT.1D-16) GOTO 10
 
330
      RN=RAN
 
331
*
 
332
      END
 
333
 
 
334
 
 
335
 
 
336
      SUBROUTINE RANMAR(RVEC)
 
337
*     -----------------
 
338
* Universal random number generator proposed by Marsaglia and Zaman
 
339
* in report FSU-SCRI-87-50
 
340
* In this version RVEC is a double precision variable.
 
341
      IMPLICIT REAL*8(A-H,O-Z)
 
342
      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
 
343
      COMMON/ RASET2 / IRANMR,JRANMR
 
344
      SAVE /RASET1/,/RASET2/
 
345
      UNI = RANU(IRANMR) - RANU(JRANMR)
 
346
      IF(UNI .LT. 0D0) UNI = UNI + 1D0
 
347
      RANU(IRANMR) = UNI
 
348
      IRANMR = IRANMR - 1
 
349
      JRANMR = JRANMR - 1
 
350
      IF(IRANMR .EQ. 0) IRANMR = 97
 
351
      IF(JRANMR .EQ. 0) JRANMR = 97
 
352
      RANC = RANC - RANCD
 
353
      IF(RANC .LT. 0D0) RANC = RANC + RANCM
 
354
      UNI = UNI - RANC
 
355
      IF(UNI .LT. 0D0) UNI = UNI + 1D0
 
356
      RVEC = UNI
 
357
      END
 
358
 
 
359
      SUBROUTINE RMARIN(IJ,KL)
 
360
*     -----------------
 
361
* Initializing routine for RANMAR, must be called before generating
 
362
* any pseudorandom numbers with RANMAR. The input values should be in
 
363
* the ranges 0<=ij<=31328 ; 0<=kl<=30081
 
364
      IMPLICIT REAL*8(A-H,O-Z)
 
365
      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
 
366
      COMMON/ RASET2 / IRANMR,JRANMR
 
367
      SAVE /RASET1/,/RASET2/
 
368
* This shows correspondence between the simplified input seeds IJ, KL
 
369
* and the original Marsaglia-Zaman seeds I,J,K,L.
 
370
* To get the standard values in the Marsaglia-Zaman paper (i=12,j=34
 
371
* k=56,l=78) put ij=1802, kl=9373
 
372
      I = MOD( IJ/177 , 177 ) + 2
 
373
      J = MOD( IJ     , 177 ) + 2
 
374
      K = MOD( KL/169 , 178 ) + 1
 
375
      L = MOD( KL     , 169 )
 
376
      DO 300 II = 1 , 97
 
377
        S =  0D0
 
378
        T = .5D0
 
379
        DO 200 JJ = 1 , 24
 
380
          M = MOD( MOD(I*J,179)*K , 179 )
 
381
          I = J
 
382
          J = K
 
383
          K = M
 
384
          L = MOD( 53*L+1 , 169 )
 
385
          IF(MOD(L*M,64) .GE. 32) S = S + T
 
386
          T = .5D0*T
 
387
  200   CONTINUE
 
388
        RANU(II) = S
 
389
  300 CONTINUE
 
390
      RANC  =   362436D0 / 16777216D0
 
391
      RANCD =  7654321D0 / 16777216D0
 
392
      RANCM = 16777213D0 / 16777216D0
 
393
      IRANMR = 97
 
394
      JRANMR = 33
 
395
      END
 
396
 
 
397
 
 
398
 
 
399
 
 
400
 
 
401