← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to vendor/CutTools/src/cts/cts_kinematics.f90

merge with 2.0.0beta4 revision 216

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/cts/cts_kinematics.f90
255
255
  end interface!build_l
256
256
  contains
257
257
!
258
 
  subroutine dp_build_l(k1,k2,l1,l2,l3,l4,al1,al2,bet,ga)  
 
258
  subroutine dp_build_l(k1,k2in,l1,l2,l3,l4,al1,al2,bet,ga)  
 
259
   use scale
259
260
   include 'cts_dpr.h' 
260
261
    :: p
261
262
   include 'cts_dpr.h'
262
 
    , intent(in), dimension(0:3) :: k1,k2
 
263
    , intent(in), dimension(0:3) :: k1,k2in
 
264
   include 'cts_dpr.h'
 
265
    , dimension(0:3) :: k2
263
266
   include 'cts_dpc.h'
264
267
    , intent(out), dimension(0:3) :: l1,l2,l3,l4
265
268
   include 'cts_dpc.h'
266
269
    , intent(out) :: al1,al2,bet,ga 
267
270
   include 'cts_dpr.h'
268
271
    :: k1k1,k1k2,k2k2
 
272
   include 'cts_dpr.h'
 
273
    :: d12a,d12b,d12c,d12
 
274
   include 'cts_dpr.h'
 
275
    , dimension(1:3) :: s10,s20,d10,d20
269
276
   include 'cts_dpc.h' 
270
 
    :: del12,b1p,b1m,b2p,b2m,c1p,c1m,c2p,c2m,ausp,ausm
 
277
    :: del12,rdel12,b1p,b1m,b2p,b2m,c1p,c1m,c2p,c2m,ausp,ausm
271
278
   integer :: k
 
279
   k2= k2in
 
280
 1 do k= 1,3
 
281
    s10(k)= k1(0)+k1(k); d10(k)= k1(0)-k1(k) 
 
282
    s20(k)= k2(0)+k2(k); d20(k)= k2(0)-k2(k) 
 
283
   enddo
 
284
   d12a= (k2(1)*d10(3)-k1(1)*d20(3))*(k2(1)*s10(3)-k1(1)*s20(3))
 
285
   d12b= (k2(2)*s10(1)-k1(2)*s20(1))*(k2(2)*d10(1)-k1(2)*d20(1))
 
286
   d12c= (k2(3)*s10(2)-k1(3)*s20(2))*(k2(3)*d10(2)-k1(3)*d20(2))
 
287
   d12 =  d12a+d12b+d12c
 
288
   if (abs(d12/roots**4).lt.1.d-60) then
 
289
     k2(0)= k2(0)*1.0000000000001d0 
 
290
     goto 1 
 
291
   endif
 
292
   del12= d12*c1(p)
 
293
   rdel12= sqrt(del12)
272
294
   call contr(k1,k1,k1k1)
273
295
   call contr(k1,k2,k1k2)
274
296
   call contr(k2,k2,k2k2)
275
 
   del12= (k1k2*k1k2-k1k1*k2k2)*c1(p)
276
297
   if (k1k2.gt.0.d0) then
277
 
    ga= k1k2+sqrt(del12)
 
298
    k = 1 
278
299
   else
279
 
    ga= k1k2-sqrt(del12)
 
300
    k =-1 
280
301
   endif
 
302
   ga= k1k2+k*rdel12
281
303
   al1= k1k1/ga
282
304
   al2= k2k2/ga
283
 
   bet= c1(p)/(c1(p)-al1*al2)
 
305
   if (dreal(al1*al2).lt.0.d0) then
 
306
     bet= c1(p)/(c1(p)-al1*al2)
 
307
   else
 
308
     bet= k*ga**2*(c1(p)+al1*al2)/4.d0/k1k2/rdel12
 
309
   endif
284
310
   do k= 0,3
285
311
    l1(k)= bet*(k1(k)-al1*k2(k))
286
312
    l2(k)= bet*(k2(k)-al2*k1(k))
321
347
   l4(3)=     b2m*b1p  - c2m*c1p
322
348
  end subroutine dp_build_l
323
349
!
324
 
  subroutine mp_build_l(k1,k2,l1,l2,l3,l4,al1,al2,bet,ga)  
 
350
  subroutine mp_build_l(k1,k2in,l1,l2,l3,l4,al1,al2,bet,ga)  
 
351
   use scale
325
352
   include 'cts_mpr.h' 
326
353
    :: p
327
354
   include 'cts_mpr.h'
328
 
    , intent(in), dimension(0:3) :: k1,k2
 
355
    , intent(in), dimension(0:3) :: k1,k2in
 
356
   include 'cts_mpr.h'
 
357
    , dimension(0:3) :: k2
329
358
   include 'cts_mpc.h'
330
359
    , intent(out), dimension(0:3) :: l1,l2,l3,l4
331
360
   include 'cts_mpc.h'
332
361
    , intent(out) :: al1,al2,bet,ga 
333
362
   include 'cts_mpr.h'
334
 
    :: k1k1,k1k2,k2k2
 
363
    :: k1k1,k1k2,k2k2,r12
 
364
   include 'cts_mpr.h'
 
365
    :: d12a,d12b,d12c,d12
 
366
   include 'cts_mpr.h'
 
367
    , dimension(1:3) :: s10,s20,d10,d20
335
368
   include 'cts_mpc.h' 
336
 
    :: del12,b1p,b1m,b2p,b2m,c1p,c1m,c2p,c2m,ausp,ausm
 
369
    :: del12,rdel12,b1p,b1m,b2p,b2m,c1p,c1m,c2p,c2m,ausp,ausm
337
370
   integer :: k
 
371
   k2= k2in
 
372
 1 do k= 1,3
 
373
    s10(k)= k1(0)+k1(k); d10(k)= k1(0)-k1(k) 
 
374
    s20(k)= k2(0)+k2(k); d20(k)= k2(0)-k2(k) 
 
375
   enddo
 
376
   d12a= (k2(1)*d10(3)-k1(1)*d20(3))*(k2(1)*s10(3)-k1(1)*s20(3))
 
377
   d12b= (k2(2)*s10(1)-k1(2)*s20(1))*(k2(2)*d10(1)-k1(2)*d20(1))
 
378
   d12c= (k2(3)*s10(2)-k1(3)*s20(2))*(k2(3)*d10(2)-k1(3)*d20(2))
 
379
   d12 =  d12a+d12b+d12c
 
380
   if (abs(d12/roots**4).lt.1.d-120) then
 
381
     k2(0)= k2(0)*1.0000000000001d0 
 
382
     goto 1 
 
383
   endif
 
384
   del12= d12*c1(p)
 
385
   rdel12= sqrt(del12)
338
386
   call contr(k1,k1,k1k1)
339
387
   call contr(k1,k2,k1k2)
340
388
   call contr(k2,k2,k2k2)
341
 
   del12= (k1k2*k1k2-k1k1*k2k2)*c1(p)
342
389
   if (k1k2.gt.0.d0) then
343
 
    ga= k1k2+sqrt(del12)
 
390
    k = 1
344
391
   else
345
 
    ga= k1k2-sqrt(del12)
 
392
    k =-1
346
393
   endif
 
394
   ga= k1k2+k*rdel12
347
395
   al1= k1k1/ga
348
396
   al2= k2k2/ga
349
 
   bet= c1(p)/(c1(p)-al1*al2)
 
397
   r12= (al1*al2+conjg(al1*al2))
 
398
   if (r12.lt.0.d0) then
 
399
     bet= c1(p)/(c1(p)-al1*al2)
 
400
   else
 
401
     bet= k*ga**2*(c1(p)+al1*al2)/4.d0/k1k2/rdel12
 
402
   endif
350
403
   do k= 0,3
351
404
    l1(k)= bet*(k1(k)-al1*k2(k))
352
405
    l2(k)= bet*(k2(k)-al2*k1(k))
549
602
   x10= z*(dd2-al2*dd1-dd0*(1.d0-al2))
550
603
   x20= z*(dd1-al1*dd2-dd0*(1.d0-al1))
551
604
   cc = 0.25d0*(x10*x20-dd0/gm)
 
605
!  comment: the next line to avoid underflows 
 
606
   if (abs(cc).lt.1.d-60) cc= 0.d0
552
607
   cut3%gm= gm
553
608
   cut3%cc= cc
554
609
   do i= 0,3
3736
3791
       l4vec(k,ib) =  cut3%l4(k)
3737
3792
       p0vecc(k,ib) =  den(bbn3(1,ib))%p(k) 
3738
3793
      enddo
3739
 
      cc = cut3%cc; if (abs(cc).lt.1.d-60) cc= 0.d0;! fix 20.12.2012
 
3794
      cc = cut3%cc
3740
3795
      tau= cut3%tau
3741
3796
      c34= -2.d0*cut3%gm
3742
3797
      cden1= cc*tau