← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/fks_singular.f

  • Committer: olivier Mattelaer
  • Date: 2016-05-12 11:00:18 UTC
  • mfrom: (262.1.150 2.3.4)
  • Revision ID: olivier.mattelaer@uclouvain.be-20160512110018-sevb79f0wm4g8mpp
pass to 2.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/fks_singular.f
38
38
      include 'run.inc'
39
39
      include 'timing_variables.inc'
40
40
      double precision wgt1,wgt2,wgt3,bsv_wgt,virt_wgt,born_wgt,pi,g2
41
 
     &     ,g22
 
41
     &     ,g22,wgt4
42
42
      parameter (pi=3.1415926535897932385d0)
43
43
      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
44
44
     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
61
61
      g2=g**(nint(2*wgtbpower))
62
62
      g22=g**(nint(2*wgtbpower+2))
63
63
      wgt1=wgtnstmp*f_nb/g22
 
64
      wgt4=wgtnstmp_avgvirt*f_nb/g22
64
65
      if (ickkw.eq.3 .and. fxfx_exp_rewgt.ne.0d0) then
65
66
         wgt1=wgt1 - fxfx_exp_rewgt*born_wgt*f_nb/g2/(4d0*pi)
66
67
      elseif (ickkw.eq.-1) then
78
79
      wgt2=wgtwnstmpmur*f_nb/g22
79
80
      wgt3=wgtwnstmpmuf*f_nb/g22
80
81
      call add_wgt(3,wgt1,wgt2,wgt3)
 
82
      call add_wgt(15,wgt4,0d0,0d0)
81
83
c Special for the soft-virtual needed for the virt-tricks. The
82
84
c *_wgt_mint variable should be directly passed to the mint-integrator
83
85
c and not be part of the plots nor computation of the cross section.
955
957
c     type=12: MC subtraction with n-body kin.
956
958
c     type=13: MC subtraction with n+1-body kin.
957
959
c     type=14: virtual corrections
 
960
c     type=15: virt-trick: average born contribution
958
961
c     wgt1 : weight of the contribution not multiplying a scale log
959
962
c     wgt2 : coefficient of the weight multiplying the log[mu_R^2/Q^2]
960
963
c     wgt3 : coefficient of the weight multiplying the log[mu_F^2/Q^2]
990
993
c        PDG=21 (gluon) code.
991
994
c     If the contribution belongs to an H-event or S-event:
992
995
c        H_event(icontr)
 
996
c     The weight of the born or real-emission matrix element
 
997
c        corresponding to this contribution: wgt_ME_tree. This weight does
 
998
c        include the 'ngluon' correction factor for the Born.
993
999
c
994
1000
c Not set in this subroutine, but included in the c_weights common block
995
1001
c are the
1033
1039
      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
1034
1040
     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
1035
1041
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
 
1042
      double precision         fkssymmetryfactor,fkssymmetryfactorBorn,
 
1043
     &                         fkssymmetryfactorDeg
 
1044
      integer                                      ngluons,nquarks(-6:6)
 
1045
      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
 
1046
     &                         fkssymmetryfactorDeg,ngluons,nquarks
 
1047
      double precision       wgt_ME_born,wgt_ME_real
 
1048
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
 
1049
      double precision     iden_comp
 
1050
      common /c_iden_comp/ iden_comp
1036
1051
      if (wgt1.eq.0d0 .and. wgt2.eq.0d0 .and. wgt3.eq.0d0) return
1037
1052
c Check for NaN's and INF's. Simply skip the contribution
1038
1053
      if (wgt1.ne.wgt1) return
1064
1079
c     Anything else
1065
1080
         QCDpower(icontr)=nint(2*wgtbpower+2)
1066
1081
      endif
 
1082
c Compensate for the fact that in the Born matrix elements, we use the
 
1083
c identical particle symmetry factor of the corresponding real emission
 
1084
c matrix elements
 
1085
c IDEN_COMP STUFF NEEDS TO BE UPDATED WHEN MERGING WITH 'FKS_EW' STUFF
 
1086
      wgt_ME_tree(1,icontr)=wgt_me_born
 
1087
      wgt_ME_tree(2,icontr)=wgt_me_real
 
1088
      do i=1,nexternal
 
1089
         do j=0,3
 
1090
            if (p1_cnt(0,1,0).gt.0d0) then
 
1091
               momenta_m(j,i,1,icontr)=p1_cnt(j,i,0)
 
1092
            elseif (p1_cnt(0,1,1).gt.0d0) then
 
1093
               momenta_m(j,i,1,icontr)=p1_cnt(j,i,1)
 
1094
            elseif (p1_cnt(0,1,2).gt.0d0) then
 
1095
               momenta_m(j,i,1,icontr)=p1_cnt(j,i,2)
 
1096
            else
 
1097
               if (i.lt.fks_i_d(nFKSprocess)) then
 
1098
                  momenta_m(j,i,1,icontr)=p_born(j,i)
 
1099
               elseif(i.eq.fks_i_d(nFKSprocess)) then
 
1100
                  momenta_m(j,i,1,icontr)=0d0
 
1101
               else
 
1102
                  momenta_m(j,i,1,icontr)=p_born(j,i-1)
 
1103
               endif
 
1104
            endif
 
1105
            momenta_m(j,i,2,icontr)=p_ev(j,i)
 
1106
         enddo
 
1107
      enddo
1067
1108
      if(type.eq.1 .or. type.eq. 8 .or. type.eq.9 .or. type.eq.10 .or.
1068
1109
     &     type.eq.13) then
1069
1110
c real emission and n+1-body kin. contributions to counter terms and MC
1070
1111
c subtr term
1071
1112
         do i=1,nexternal
1072
1113
            do j=0,3
1073
 
c     'momenta' is the momentum configuration for this contribution;
1074
 
c     'momenta_m' is the momentum configuration that was used in the
1075
 
c     matrix elements of this contribution.
1076
 
               momenta(j,i,icontr)=p_ev(j,i)
1077
 
               if (type.eq.1) then
1078
 
                  momenta_m(j,i,icontr)=momenta(j,i,icontr)
1079
 
               else
1080
 
                  if (p1_cnt(0,1,0).gt.0d0) then
1081
 
                     momenta_m(j,i,icontr)=p1_cnt(j,i,0)
1082
 
                  elseif (p1_cnt(0,1,1).gt.0d0) then
1083
 
                     momenta_m(j,i,icontr)=p1_cnt(j,i,1)
1084
 
                  elseif (p1_cnt(0,1,2).gt.0d0) then
1085
 
                     momenta_m(j,i,icontr)=p1_cnt(j,i,2)
1086
 
                  else
1087
 
                     write (*,*) 'ERROR in add_wgt: no valid momenta'
1088
 
                     stop 1
1089
 
                  endif
1090
 
               endif
 
1114
               momenta(j,i,icontr)=momenta_m(j,i,2,icontr)
1091
1115
            enddo
1092
1116
         enddo
1093
1117
         H_event(icontr)=.true.
1094
1118
      elseif(type.ge.2 .and. type.le.7 .or. type.eq.11 .or. type.eq.12
1095
 
     $        .or. type.eq.14)then
 
1119
     $        .or. type.eq.14 .or. type.eq.15)then
1096
1120
c Born, counter term, soft-virtual, or n-body kin. contributions to real
1097
1121
c and MC subtraction terms.
1098
1122
         do i=1,nexternal
1099
1123
            do j=0,3
1100
 
c     'momenta' is the momentum configuration for this contribution;
1101
 
c     'momenta_m' is the momentum configuration that was used in the
1102
 
c     matrix elements of this contribution.
1103
 
               if (p1_cnt(0,1,0).gt.0d0) then
1104
 
                  momenta(j,i,icontr)=p1_cnt(j,i,0)
1105
 
               elseif (p1_cnt(0,1,1).gt.0d0) then
1106
 
                  momenta(j,i,icontr)=p1_cnt(j,i,1)
1107
 
               elseif (p1_cnt(0,1,2).gt.0d0) then
1108
 
                  momenta(j,i,icontr)=p1_cnt(j,i,2)
1109
 
               else
1110
 
                  write (*,*) 'ERROR in add_wgt: no valid momenta'
1111
 
                  stop 1
1112
 
               endif
1113
 
               if (type.eq.11) then
1114
 
                  momenta_m(j,i,icontr)=p_ev(j,i)
1115
 
               else
1116
 
                  momenta_m(j,i,icontr)=momenta(j,i,icontr)
1117
 
               endif
 
1124
               momenta(j,i,icontr)=momenta_m(j,i,1,icontr)
1118
1125
            enddo
1119
1126
         enddo
1120
1127
         H_event(icontr)=.false.
1276
1283
      include 'reweight.inc'
1277
1284
      include 'reweightNLO.inc'
1278
1285
      include 'timing_variables.inc'
1279
 
      integer i,kr,kf,iwgt_save
1280
 
      double precision xlum(maxscales),dlum,pi,mu2_r(maxscales)
1281
 
     &     ,mu2_f(maxscales),mu2_q,alphas,g(maxscales),rwgt_muR_dep_fac
 
1286
      integer i,kr,kf,iwgt_save,dd
 
1287
      double precision xlum(maxscales),dlum,pi,mu2_r(maxscales),c_mu2_r
 
1288
     $     ,c_mu2_f,mu2_f(maxscales),mu2_q,alphas,g(maxscales)
 
1289
     $     ,rwgt_muR_dep_fac
1282
1290
      external rwgt_muR_dep_fac
1283
1291
      parameter (pi=3.1415926535897932385d0)
1284
1292
      external dlum,alphas
1297
1305
         xbk(1) = bjx(1,i)
1298
1306
         xbk(2) = bjx(2,i)
1299
1307
         mu2_q=scales2(1,i)
 
1308
c Loop over the dynamical_scale_choices
 
1309
         do dd=1,dyn_scale(0)
1300
1310
c renormalisation scale variation (requires recomputation of the strong
1301
1311
c coupling)
1302
 
         do kr=1,numscales
1303
 
            mu2_r(kr)=scales2(2,i)*yfactR(kr)**2
1304
 
            g(kr)=sqrt(4d0*pi*alphas(sqrt(mu2_r(kr))))
1305
 
         enddo
 
1312
            call set_mu_central(i,dd,c_mu2_r,c_mu2_f)
 
1313
            do kr=1,nint(scalevarR(0))
 
1314
               if ((.not. lscalevar(dd)) .and. kr.ne.1) exit
 
1315
               mu2_r(kr)=c_mu2_r*scalevarR(kr)**2
 
1316
               g(kr)=sqrt(4d0*pi*alphas(sqrt(mu2_r(kr))))
 
1317
            enddo
1306
1318
c factorisation scale variation (require recomputation of the PDFs)
1307
 
         do kf=1,numscales
1308
 
            mu2_f(kf)=scales2(3,i)*yfactF(kf)**2
1309
 
            q2fact(1)=mu2_f(kf)
1310
 
            q2fact(2)=mu2_f(kf)
1311
 
            xlum(kf) = dlum()
1312
 
         enddo
1313
 
         do kr=1,numscales
1314
 
            do kf=1,numscales
1315
 
               iwgt=iwgt+1 ! increment the iwgt for the wgts() array
1316
 
               if (iwgt.gt.max_wgt) then
1317
 
                  write (*,*) 'ERROR too many weights in reweight_scale'
1318
 
     &                 ,iwgt,max_wgt
1319
 
                  stop 1
1320
 
               endif
 
1319
            do kf=1,nint(scalevarF(0))
 
1320
               if ((.not. lscalevar(dd)) .and. kf.ne.1) exit
 
1321
               mu2_f(kf)=c_mu2_f*scalevarF(kf)**2
 
1322
               q2fact(1)=mu2_f(kf)
 
1323
               q2fact(2)=mu2_f(kf)
 
1324
               xlum(kf) = dlum()
 
1325
            enddo
 
1326
            do kf=1,nint(scalevarF(0))
 
1327
               if ((.not. lscalevar(dd)) .and. kf.ne.1) exit
 
1328
               do kr=1,nint(scalevarR(0))
 
1329
                  if ((.not. lscalevar(dd)) .and. kr.ne.1) exit
 
1330
                  iwgt=iwgt+1   ! increment the iwgt for the wgts() array
 
1331
                  if (iwgt.gt.max_wgt) then
 
1332
                     write (*,*) 'ERROR too many weights in '/
 
1333
     $                    /'reweight_scale',iwgt,max_wgt
 
1334
                     stop 1
 
1335
                  endif
1321
1336
c add the weights to the array
1322
 
               wgts(iwgt,i)=xlum(kf) * (wgt(1,i)+wgt(2,i)*log(mu2_r(kr)
1323
 
     &              /mu2_q)+wgt(3,i)*log(mu2_f(kf)/mu2_q))*g(kr)
1324
 
     &              **QCDpower(i)
1325
 
               wgts(iwgt,i)=wgts(iwgt,i)
1326
 
     &              *rwgt_muR_dep_fac(sqrt(mu2_r(kr)),sqrt(mu2_r(1)))
 
1337
                  wgts(iwgt,i)=xlum(kf) * (wgt(1,i)+wgt(2,i)
 
1338
     $                 *log(mu2_r(kr)/mu2_q)+wgt(3,i)*log(mu2_f(kf)
 
1339
     $                 /mu2_q))*g(kr)**QCDpower(i)
 
1340
                  wgts(iwgt,i)=wgts(iwgt,i)*rwgt_muR_dep_fac(
 
1341
     &                 sqrt(mu2_r(kr)),sqrt(scales2(2,i)))
 
1342
               enddo
1327
1343
            enddo
1328
1344
         enddo
1329
1345
      enddo
1355
1371
      common/c_nFKSprocess/nFKSprocess
1356
1372
      call cpu_time(tBefore)
1357
1373
      if (icontr.eq.0) return
 
1374
      if (dyn_scale(0).gt.1) then
 
1375
         write (*,*) "When doing NNLL+NLO veto, "/
 
1376
     $        /"can only do one dynamical_scale_choice",dyn_scale(0)
 
1377
         stop
 
1378
      endif
 
1379
 
1358
1380
c currently we have 'iwgt' weights in the wgts() array.
1359
1381
      iwgt_save=iwgt
1360
1382
c compute the new veto multiplier factor      
1361
 
      do ks=1,numscales
1362
 
         do kh=1,numscales
 
1383
      do ks=1,nint(scalevarR(0))
 
1384
         if ((.not. lscalevar(1)) .and. ks.ne.1) exit
 
1385
         do kh=1,nint(scalevarF(0))
 
1386
            if ((.not. lscalevar(1)) .and. kh.ne.1) exit
1363
1387
            if (H1_factor_virt.ne.0d0) then
1364
 
               call compute_veto_multiplier(H1_factor_virt,yfactR(ks)
1365
 
     &              ,yfactF(kh),veto_multiplier_new(ks,kh))
 
1388
               call compute_veto_multiplier(H1_factor_virt,scalevarR(ks)
 
1389
     $              ,scalevarF(kh),veto_multiplier_new(ks,kh))
1366
1390
               veto_multiplier_new(ks,kh)=veto_multiplier_new(ks,kh)
1367
1391
     &              /veto_multiplier
1368
1392
            else
1377
1401
         xbk(1) = bjx(1,i)
1378
1402
         xbk(2) = bjx(2,i)
1379
1403
         mu2_q=scales2(1,i)
 
1404
c Hard scale variation
 
1405
         do kh=1,nint(scalevarF(0))
 
1406
            if ((.not. lscalevar(1)) .and. kh.ne.1) exit
1380
1407
c soft scale variation
1381
 
         do ks=1,numscales
1382
 
            mu2_r(ks)=scales2(2,i)*yfactR(ks)**2
1383
 
            g(ks)=sqrt(4d0*pi*alphas(sqrt(mu2_r(ks))))
1384
 
            mu2_f(ks)=scales2(2,i)*yfactR(ks)**2
1385
 
            q2fact(1)=mu2_f(ks)
1386
 
            q2fact(2)=mu2_f(ks)
1387
 
            xlum(ks) = dlum()
1388
 
c Hard scale variation
1389
 
            do kh=1,numscales
1390
 
               iwgt=iwgt+1 ! increment the iwgt for the wgts() array
 
1408
            do ks=1,nint(scalevarR(0))
 
1409
               if ((.not. lscalevar(1)) .and. ks.ne.1) exit
 
1410
               mu2_r(ks)=scales2(2,i)*scalevarR(ks)**2
 
1411
               g(ks)=sqrt(4d0*pi*alphas(sqrt(mu2_r(ks))))
 
1412
               mu2_f(ks)=scales2(2,i)*scalevarR(ks)**2
 
1413
               q2fact(1)=mu2_f(ks)
 
1414
               q2fact(2)=mu2_f(ks)
 
1415
               xlum(ks) = dlum()
 
1416
               iwgt=iwgt+1      ! increment the iwgt for the wgts() array
1391
1417
               if (iwgt.gt.max_wgt) then
1392
1418
                  write (*,*) 'ERROR too many weights in reweight_scale'
1393
1419
     &                 ,iwgt,max_wgt
1401
1427
               else
1402
1428
c special for the itype=7 (i.e, the veto-compensating factor)                  
1403
1429
                  call compute_veto_compensating_factor(H1_factor_virt
1404
 
     &                 ,born_wgt_veto,yfactR(ks),yfactF(kh)
 
1430
     &                 ,born_wgt_veto,scalevarR(ks),scalevarF(kh)
1405
1431
     &                 ,veto_compensating_factor_new)
1406
1432
                  wgts(iwgt,i)=xlum(ks) * wgt(1,i)*g(ks)**QCDpower(i)
1407
1433
     &                 /veto_compensating_factor
1408
1434
     &                 *veto_compensating_factor_new
1409
1435
               endif
1410
 
               wgts(iwgt,i)=wgts(iwgt,i)
1411
 
     &              *rwgt_muR_dep_fac(sqrt(mu2_r(ks)),sqrt(mu2_r(1)))
 
1436
               wgts(iwgt,i)=wgts(iwgt,i)*rwgt_muR_dep_fac(
 
1437
     &              sqrt(mu2_r(ks)),sqrt(scales2(2,i)))
1412
1438
               wgts(iwgt,i)=wgts(iwgt,i)*veto_multiplier_new(ks,kh)
1413
1439
            enddo
1414
1440
         enddo
1428
1454
      include 'reweight.inc'
1429
1455
      include 'reweightNLO.inc'
1430
1456
      include 'timing_variables.inc'
1431
 
      integer n,izero,i
1432
 
      parameter (izero=0)
1433
 
      double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac
 
1457
      integer n,i,nn
 
1458
      double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac,g
 
1459
     &     ,alphas
1434
1460
      external rwgt_muR_dep_fac
1435
1461
      parameter (pi=3.1415926535897932385d0)
1436
1462
      external dlum,alphas
1438
1464
      common/c_nFKSprocess/nFKSprocess
1439
1465
      call cpu_time(tBefore)
1440
1466
      if (icontr.eq.0) return
 
1467
      do nn=1,lhaPDFid(0)
1441
1468
c Use as external loop the one over the PDF sets and as internal the one
1442
1469
c over the icontr. This reduces the number of calls to InitPDF and
1443
1470
c allows for better caching of the PDFs
1444
 
      do n=1,numPDFs-1
1445
 
         iwgt=iwgt+1
1446
 
         if (iwgt.gt.max_wgt) then
1447
 
            write (*,*) 'ERROR too many weights in reweight_pdf',iwgt
1448
 
     &           ,max_wgt
1449
 
            stop 1
1450
 
         endif
1451
 
         call InitPDF(n)
1452
 
         do i=1,icontr
1453
 
            nFKSprocess=nFKS(i)
1454
 
            xbk(1) = bjx(1,i)
1455
 
            xbk(2) = bjx(2,i)
1456
 
            mu2_q=scales2(1,i)
1457
 
            mu2_r=scales2(2,i)
1458
 
            mu2_f=scales2(3,i)
1459
 
            q2fact(1)=mu2_f
1460
 
            q2fact(2)=mu2_f
1461
 
            xlum = dlum()
 
1471
         do n=0,nmemPDF(nn)
 
1472
            iwgt=iwgt+1
 
1473
            if (iwgt.gt.max_wgt) then
 
1474
               write (*,*) 'ERROR too many weights in reweight_pdf',iwgt
 
1475
     &              ,max_wgt
 
1476
               stop 1
 
1477
            endif
 
1478
            call InitPDFm(nn,n)
 
1479
            do i=1,icontr
 
1480
               nFKSprocess=nFKS(i)
 
1481
               xbk(1) = bjx(1,i)
 
1482
               xbk(2) = bjx(2,i)
 
1483
               mu2_q=scales2(1,i)
 
1484
               mu2_r=scales2(2,i)
 
1485
               mu2_f=scales2(3,i)
 
1486
               q2fact(1)=mu2_f
 
1487
               q2fact(2)=mu2_f
 
1488
c Compute the luminosity
 
1489
               xlum = dlum()
 
1490
c Recompute the strong coupling: alpha_s in the PDF might change
 
1491
               g=sqrt(4d0*pi*alphas(sqrt(mu2_r)))
1462
1492
c add the weights to the array
1463
 
            wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q) +
1464
 
     &           wgt(3,i)*log(mu2_f/mu2_q))*g_strong(i)**QCDpower(i)
1465
 
            wgts(iwgt,i)=wgts(iwgt,i)*
1466
 
     &                       rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r))
 
1493
               wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q)
 
1494
     $              +wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
 
1495
               wgts(iwgt,i)=wgts(iwgt,i)*
 
1496
     &              rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r))
 
1497
            enddo
1467
1498
         enddo
1468
1499
      enddo
1469
 
      call InitPDF(izero)
 
1500
      call InitPDFm(1,0)
1470
1501
      call cpu_time(tAfter)
1471
1502
      tr_pdf=tr_pdf+(tAfter-tBefore)
1472
1503
      return
1490
1521
      integer iproc_save(fks_configs),eto(maxproc,fks_configs),
1491
1522
     &     etoi(maxproc,fks_configs),maxproc_found
1492
1523
      common/cproc_combination/iproc_save,eto,etoi,maxproc_found
1493
 
      if (icontr.gt.7) then
 
1524
      if (icontr.gt.8) then
1494
1525
         write (*,*) 'ERROR: too many applgrid weights. '/
1495
1526
     &        /'Should have at most one of each itype.',icontr
1496
1527
         stop 1
1532
1563
            appl_QES2(2)=scales2(1,i)
1533
1564
            appl_muR2(2)=scales2(2,i)
1534
1565
            appl_muF2(2)=scales2(3,i)
1535
 
         elseif (itype(i).eq.3 .or. itype(i).eq.4 .or. itype(i).eq.14)
1536
 
     $           then
 
1566
         elseif (itype(i).eq.3 .or. itype(i).eq.4 .or. itype(i).eq.14
 
1567
     &           .or. itype(i).eq.15)then
1537
1568
c     virtual, soft-virtual or soft-counter
1538
1569
            appl_w0(2)=appl_w0(2)+wgt(1,i)*final_state_rescaling
1539
1570
            appl_wR(2)=appl_wR(2)+wgt(2,i)*final_state_rescaling
1584
1615
      integer    maxflow
1585
1616
      parameter (maxflow=999)
1586
1617
      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
1587
 
     $     icolup(2,nexternal,maxflow)
1588
 
      common /c_leshouche_inc/idup,mothup,icolup
 
1618
     $     icolup(2,nexternal,maxflow),niprocs
 
1619
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
1589
1620
      do k=1,nexternal
1590
1621
         pdg(k,ict)=idup(k,1)
1591
1622
      enddo
1631
1662
      if (icontr.eq.0) return
1632
1663
      do i=1,icontr
1633
1664
         if (itype(i).eq.2 .or. itype(i).eq.3 .or. itype(i).eq.14 .or.
1634
 
     &        itype(i).eq.7) then
 
1665
     &        itype(i).eq.7 .or. itype(i).eq.15) then
1635
1666
            sig=sig+wgts(1,i)
1636
1667
         endif
1637
1668
      enddo
1650
1681
      if (icontr.eq.0) return
1651
1682
      do i=1,icontr
1652
1683
         if (itype(i).ne.2 .and. itype(i).ne.3 .and. itype(i).ne.14
1653
 
     &        .and. itype(i).ne.7) then
 
1684
     &        .and. itype(i).ne.7 .and. itype(i).ne.15) then
1654
1685
            sig=sig+wgts(1,i)
1655
1686
         endif
1656
1687
      enddo
2050
2081
            tmp_wgt=0d0
2051
2082
            do j=1,icontr_sum(0,i)
2052
2083
               ict=icontr_sum(j,i)
2053
 
               if (itype(ict).ne.2 .and. itype(ict).ne.3 .and.
2054
 
     $             itype(ict).ne.14) tmp_wgt=tmp_wgt+wgts(1,ict)
 
2084
               if ( itype(ict).ne.2  .and. itype(ict).ne.3 .and.
 
2085
     $              itype(ict).ne.14 .and. itype(ict).ne.15)
 
2086
     $                              tmp_wgt=tmp_wgt+wgts(1,ict)
2055
2087
            enddo
2056
2088
            n1body_wgt=n1body_wgt+abs(tmp_wgt)
2057
2089
         enddo
2158
2190
      include 'reweight0.inc'
2159
2191
      include 'genps.inc'
2160
2192
      include 'nFKSconfigs.inc'
2161
 
      integer i,ii,j,jj,ict,ipr,momenta_conf
 
2193
      include 'fks_info.inc'
 
2194
      integer k,i,ii,j,jj,ict,ipr,momenta_conf(2)
2162
2195
      logical momenta_equal,found
2163
2196
      double precision conv,momenta_str_l(0:3,nexternal,max_n_ctr)
2164
2197
      external momenta_equal
2179
2212
c Check if the current set of momenta are already available in the
2180
2213
c momenta_str_l array. If not, add it.
2181
2214
         found=.false.
2182
 
         do j=1,n_mom_conf
2183
 
            if (momenta_equal(momenta_str_l(0,1,j),momenta_m(0,1,ict)))
2184
 
     &           then
2185
 
               momenta_conf=j
2186
 
               found=.true.
2187
 
               exit
 
2215
         do k=1,2
 
2216
            do j=1,n_mom_conf
 
2217
               if (momenta_m(0,1,k,ict).le.0d0) then
 
2218
                  momenta_conf(k)=0
 
2219
                  cycle
 
2220
               endif
 
2221
               if (momenta_equal(momenta_str_l(0,1,j),
 
2222
     &                           momenta_m(0,1,k,ict))) then
 
2223
                  momenta_conf(k)=j
 
2224
                  found=.true.
 
2225
                  exit
 
2226
               endif
 
2227
            enddo
 
2228
            if (.not. found) then
 
2229
               n_mom_conf=n_mom_conf+1
 
2230
               do ii=1,nexternal
 
2231
                  do jj=0,3
 
2232
                     momenta_str(jj,ii,n_mom_conf)=
 
2233
     &                                      momenta_m(jj,ii,k,ict)
 
2234
                     momenta_str_l(jj,ii,n_mom_conf)=
 
2235
     &                                      momenta_m(jj,ii,k,ict)
 
2236
                  enddo
 
2237
               enddo
 
2238
               momenta_conf(k)=n_mom_conf
2188
2239
            endif
2189
2240
         enddo
2190
 
         if (.not. found) then
2191
 
            n_mom_conf=n_mom_conf+1
2192
 
            do ii=1,nexternal
2193
 
               do jj=0,3
2194
 
                  momenta_str(jj,ii,n_mom_conf)=momenta_m(jj,ii,ict)
2195
 
                  momenta_str_l(jj,ii,n_mom_conf)=momenta_m(jj,ii,ict)
2196
 
               enddo
2197
 
            enddo
2198
 
            momenta_conf=n_mom_conf
2199
 
         endif
2200
2241
         if (.not. Hevents) then
2201
2242
c For S-events, be careful to take all the IPROC that contribute to the
2202
2243
c iproc_picked:
2204
2245
            do ii=1,iproc_save(nFKS(ict))
2205
2246
               if (eto(ii,nFKS(ict)).ne.ipr) cycle
2206
2247
               n_ctr_found=n_ctr_found+1
 
2248
 
2207
2249
               if (nincoming.eq.2) then
2208
 
                  write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
2209
 
     &                 (wgt(j,ict)*conv,j=1,3),
 
2250
                  write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
 
2251
     &                 (wgt(j,ict)*conv,j=1,3),(wgt_me_tree(j,ict),j=1,2),
2210
2252
     &                 nexternal
2211
2253
               else
2212
 
                  write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
2213
 
     &                 (wgt(j,ict),j=1,3),
 
2254
                  write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
 
2255
     &                 (wgt(j,ict),j=1,3),(wgt_me_tree(j,ict),j=1,2), 
2214
2256
     &                 nexternal
2215
2257
               endif
 
2258
 
2216
2259
               procid=''
2217
2260
               do j=1,nexternal
2218
2261
                  write (str_temp,*) parton_pdg(j,ii,ict)
2222
2265
               n_ctr_str(n_ctr_found) =
2223
2266
     &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
2224
2267
     &              //trim(adjustl(procid))
2225
 
               write (str_temp,'(i2,5(1x,d14.8),3(1x,i2),1x,d18.12)')
 
2268
 
 
2269
               write (str_temp,'(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
2226
2270
     &              QCDpower(ict),
2227
2271
     &              (bjx(j,ict),j=1,2),
2228
2272
     &              (scales2(j,ict),j=1,3),
2229
 
     &              momenta_conf,
 
2273
     &              g_strong(ict),
 
2274
     &              (momenta_conf(j),j=1,2),
2230
2275
     &              itype(ict),
2231
2276
     &              nFKS(ict),
 
2277
     &              fks_i_d(nFKS(ict)),
 
2278
     &              fks_j_d(nFKS(ict)),
 
2279
     &              parton_pdg_uborn(fks_j_d(nFKS(ict)),ii,ict),
2232
2280
     &              parton_iproc(ii,ict)
2233
2281
               n_ctr_str(n_ctr_found) =
2234
2282
     &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
2238
2286
c H-event
2239
2287
            ipr=iproc_picked
2240
2288
            n_ctr_found=n_ctr_found+1
 
2289
 
2241
2290
            if (nincoming.eq.2) then
2242
 
               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
2243
 
     &              (wgt(j,ict)*conv,j=1,3),
 
2291
               write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
 
2292
     &              (wgt(j,ict)*conv,j=1,3),(wgt_me_tree(j,ict),j=1,2),
2244
2293
     &              nexternal
2245
2294
            else
2246
 
               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
2247
 
     &              (wgt(j,ict),j=1,3),
 
2295
               write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
 
2296
     &              (wgt(j,ict),j=1,3),(wgt_me_tree(j,ict),j=1,2),
2248
2297
     &              nexternal
2249
2298
            endif
 
2299
 
2250
2300
            procid=''
2251
2301
            do j=1,nexternal
2252
2302
               write (str_temp,*) parton_pdg(j,ipr,ict)
2256
2306
            n_ctr_str(n_ctr_found) =
2257
2307
     &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
2258
2308
     &           //trim(adjustl(procid))
2259
 
            write (str_temp,'(i2,5(1x,d14.8),3(1x,i2),1x,d18.12)')
 
2309
 
 
2310
            write (str_temp,'(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
2260
2311
     &           QCDpower(ict),
2261
2312
     &           (bjx(j,ict),j=1,2),
2262
2313
     &           (scales2(j,ict),j=1,3),
2263
 
     &           momenta_conf,
2264
 
     &              itype(ict),
2265
 
     &              nFKS(ict),
 
2314
     &           g_strong(ict),
 
2315
     &           (momenta_conf(j),j=1,2),
 
2316
     &           itype(ict),
 
2317
     &           nFKS(ict),
 
2318
     &           fks_i_d(nFKS(ict)),
 
2319
     &           fks_j_d(nFKS(ict)),
 
2320
     &           parton_pdg_uborn(fks_j_d(nFKS(ict)),ipr,ict),
2266
2321
     &           parton_iproc(ipr,ict)
2267
2322
            n_ctr_str(n_ctr_found) =
2268
2323
     &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
2269
2324
     &           //trim(adjustl(str_temp))
 
2325
 
 
2326
 
2270
2327
         endif
2271
2328
         if (n_ctr_found.ge.max_n_ctr) then
2272
2329
            write (*,*) 'ERROR: too many contributions in <rwgt>'
2848
2905
      parameter (vtiny=1d-8)
2849
2906
      double complex ximag
2850
2907
      parameter (ximag=(0.d0,1.d0))
 
2908
      double precision iden_comp
 
2909
      common /c_iden_comp/iden_comp
2851
2910
C  
2852
2911
      if(p_born(0,1).le.0.d0)then
2853
2912
c Unphysical kinematics: set matrix elements equal to zero
2919
2978
         write(*,*) 'FATAL ERROR in sborncol_fsr',i_type,j_type,i_fks,j_fks
2920
2979
         stop
2921
2980
      endif
2922
 
      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)
 
2981
      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)*iden_comp
2923
2982
      return
2924
2983
      end
2925
2984
 
2967
3026
      parameter (vtiny=1d-8)
2968
3027
      double complex ximag
2969
3028
      parameter (ximag=(0.d0,1.d0))
 
3029
      double precision iden_comp
 
3030
      common /c_iden_comp/iden_comp
2970
3031
C  
2971
3032
      if(p_born(0,1).le.0.d0)then
2972
3033
c Unphysical kinematics: set matrix elements equal to zero
3043
3104
     #             wgt1(2) * dconjg(azifact)
3044
3105
         call Qterms_reduced_spacelike(m_type, i_type, t, z, Q)
3045
3106
      endif
3046
 
      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)
 
3107
      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)*iden_comp
3047
3108
      return
3048
3109
      end
3049
3110
 
3325
3386
 
3326
3387
      double precision zero,pmass(nexternal)
3327
3388
      parameter(zero=0d0)
 
3389
      double precision iden_comp
 
3390
      common /c_iden_comp/iden_comp
3328
3391
      include "pmass.inc"
3329
3392
c
3330
3393
c Call the Born to be sure that 'CalculatedBorn' is done correctly. This
3349
3412
            endif
3350
3413
         enddo
3351
3414
      enddo
3352
 
 
3353
 
      wgt=softcontr
 
3415
      wgt=softcontr*iden_comp
3354
3416
c Add minus sign to compensate the minus in the color factor
3355
3417
c of the color-linked Borns (b_sf_0??.f)
3356
3418
c Factor two to fix the limits.
3483
3545
      double precision one,pi
3484
3546
      parameter (one=1.d0)
3485
3547
      parameter (pi=3.1415926535897932385d0)
 
3548
      double precision iden_comp
 
3549
      common /c_iden_comp/iden_comp
3486
3550
 
3487
3551
      if(j_fks.gt.nincoming)then
3488
3552
c Do not include this contribution for final-state branchings
3540
3604
c The partonic flux 1/(2*s) is inserted in genps. Thus, an extra 
3541
3605
c factor z (implicit in the flux of the reduced Born in FKS) 
3542
3606
c has to be inserted here
3543
 
      xnorm=1.d0/z
 
3607
      xnorm=1.d0/z *iden_comp
3544
3608
 
3545
3609
      collrem_xi=oo2pi * born_wgt * collrem_xi * xnorm
3546
3610
      collrem_lxi=oo2pi * born_wgt * collrem_lxi * xnorm
4339
4403
      include "run.inc"
4340
4404
      include "fks_powers.inc"
4341
4405
      include 'reweight.inc'
4342
 
      double precision p(0:3,nexternal),bsv_wgt,born_wgt
 
4406
      double precision p(0:3,nexternal),bsv_wgt,born_wgt,avv_wgt
4343
4407
      double precision pp(0:3,nexternal)
4344
4408
      
4345
4409
      double complex wgt1(2)
4437
4501
         bsv_wgt=dble(wgt1(1))
4438
4502
         born_wgt=dble(wgt1(1))
4439
4503
         virt_wgt=0d0
 
4504
         avv_wgt=0d0 
4440
4505
 
4441
4506
         if (abrv.eq.'born' .or. abrv.eq.'grid') goto 549
4442
4507
         if (abrv.eq.'virt' .or. abrv.eq.'viSC' .or.
4590
4655
c$$$            bsv_wgt=bsv_wgt+virt_wgt_save
4591
4656
         endif
4592
4657
         if (abrv(1:4).ne.'virt' .and. ickkw.ne.-1)
4593
 
     &        bsv_wgt=bsv_wgt+average_virtual*born_wgt*ao2pi
 
4658
     &        avv_wgt=average_virtual*born_wgt*ao2pi
4594
4659
 
4595
4660
c eq.(MadFKS.C.13)
4596
4661
         if(abrv.eq.'viSA'.or.abrv.eq.'viSB')then
4638
4703
            wgtnstmp=bsv_wgt-born_wgt-
4639
4704
     #                wgtwnstmpmuf*log(q2fact(1)/QES2)-
4640
4705
     #                wgtwnstmpmur*log(scale**2/QES2)
 
4706
            wgtnstmp_avgvirt = avv_wgt
4641
4707
         else
4642
4708
            wgtnstmp=0d0
4643
4709
            wgtwnstmpmur=0.d0
 
4710
            wgtnstmp_avgvirt = 0d0
4644
4711
         endif
4645
4712
 
4646
4713
         if (abrv(1:2).eq.'vi') then
5372
5439
      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
5373
5440
     &                         fkssymmetryfactorDeg,ngluons,nquarks
5374
5441
 
 
5442
      double precision iden_comp
 
5443
      common /c_iden_comp/iden_comp
 
5444
 
5375
5445
      include 'coupl.inc'
5376
5446
      include 'genps.inc'
5377
5447
      include 'nexternal.inc'
5378
5448
      include 'fks_powers.inc'
5379
5449
      include 'nFKSconfigs.inc'
 
5450
      include 'c_weight.inc'
5380
5451
      integer fks_j_from_i(nexternal,0:nexternal)
5381
5452
     &     ,particle_type(nexternal),pdg_type(nexternal)
5382
5453
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
5412
5483
      character*1 integrate
5413
5484
      integer i_fks,j_fks
5414
5485
      common/fks_indices/i_fks,j_fks
5415
 
      integer fac_i,fac_j,i_fks_pdg,j_fks_pdg
 
5486
      integer fac_i,fac_j,i_fks_pdg,j_fks_pdg,iden(nexternal)
5416
5487
 
5417
5488
      integer fac_i_FKS(fks_configs),fac_j_FKS(fks_configs)
5418
 
     $     ,i_type_FKS(fks_configs),j_type_FKS(fks_configs)
5419
 
     $     ,m_type_FKS(fks_configs),ngluons_FKS(fks_configs)
 
5489
     &     ,i_type_FKS(fks_configs),j_type_FKS(fks_configs)
 
5490
     &     ,m_type_FKS(fks_configs),ngluons_FKS(fks_configs)
 
5491
     &     ,iden_real_FKS(fks_configs),iden_born_FKS(fks_configs)
5420
5492
      save fac_i_FKS,fac_j_FKS,i_type_FKS,j_type_FKS,m_type_FKS
5421
 
     $     ,ngluons_FKS
 
5493
     $     ,ngluons_FKS,iden_real_FKS,iden_born_FKS
5422
5494
 
5423
5495
      character*13 filename
5424
5496
 
5594
5666
         i_type_FKS(nFKSprocess)=i_type
5595
5667
         j_type_FKS(nFKSprocess)=j_type
5596
5668
         m_type_FKS(nFKSprocess)=m_type
 
5669
 
 
5670
 
 
5671
c Compute the identical particle symmetry factor that is in the
 
5672
c real-emission matrix elements.
 
5673
         iden_real_FKS(nFKSprocess)=1
 
5674
         do i=1,nexternal
 
5675
            iden(i)=1
 
5676
         enddo
 
5677
         do i=nincoming+2,nexternal
 
5678
            do j=nincoming+1,i-1
 
5679
               if (pdg_type(j).eq.pdg_type(i)) then
 
5680
                  iden(j)=iden(j)+1
 
5681
                  iden_real_FKS(nFKSprocess)=
 
5682
     &                 iden_real_FKS(nFKSprocess)*iden(j)
 
5683
                  exit
 
5684
               endif
 
5685
            enddo
 
5686
         enddo
 
5687
c Compute the identical particle symmetry factor that is in the
 
5688
c Born matrix elements.
 
5689
         iden_born_FKS(nFKSprocess)=1
 
5690
         call set_pdg(0,nFKSprocess)
 
5691
         do i=1,nexternal
 
5692
            iden(i)=1
 
5693
         enddo
 
5694
         do i=nincoming+2,nexternal-1
 
5695
            do j=nincoming+1,i-1
 
5696
               if (pdg_uborn(j,0).eq.pdg_uborn(i,0)) then
 
5697
                  iden(j)=iden(j)+1
 
5698
                  iden_born_FKS(nFKSprocess)=
 
5699
     &                 iden_born_FKS(nFKSprocess)*iden(j)
 
5700
                  exit
 
5701
               endif
 
5702
            enddo
 
5703
         enddo
5597
5704
      endif
5598
5705
 
5599
5706
      i_type=i_type_FKS(nFKSprocess)
5600
5707
      j_type=j_type_FKS(nFKSprocess)
5601
5708
      m_type=m_type_FKS(nFKSprocess)
5602
5709
 
 
5710
c Difference in identical particle factor in the Born and real emission
 
5711
c matrix elements. To define wgt_ME_tree for the Born, we need to
 
5712
c include this factor, because in the current Born the symmetry factor
 
5713
c for the real is used. THIS NEEDS TO BE CHANGED WHEN MERGING WITH THE
 
5714
c 'FKS_EW' STUFF
 
5715
      iden_comp=dble(iden_born_FKS(nFKSprocess))/
 
5716
     &          dble(iden_real_FKS(nFKSprocess))
 
5717
      
 
5718
      
5603
5719
c Set matrices used by MC counterterms
5604
5720
      if (match_to_shower) call set_mc_matrices
5605
5721
 
5610
5726
      if (nbody.and.pdg_type(i_fks).eq.21) then
5611
5727
         fkssymmetryfactor=dble(ngluons)
5612
5728
         fkssymmetryfactorDeg=dble(ngluons)
5613
 
         fkssymmetryfactorBorn=dble(ngluons)
 
5729
         fkssymmetryfactorBorn=1d0
5614
5730
      elseif(pdg_type(i_fks).eq.-21) then
5615
5731
         fkssymmetryfactor=1d0
5616
5732
         fkssymmetryfactorDeg=1d0
6006
6122
      FK88RANDOM = SEED*MINV
6007
6123
      END
6008
6124
 
 
6125
 
 
6126
      subroutine set_mu_central(ic,dd,c_mu2_r,c_mu2_f)
 
6127
      implicit none
 
6128
      include 'nexternal.inc'
 
6129
      include 'c_weight.inc'
 
6130
      include 'reweight0.inc'
 
6131
      include 'run.inc'
 
6132
      integer ic,dd,i,j
 
6133
      double precision c_mu2_r,c_mu2_f,muR,muF,pp(0:3,nexternal)
 
6134
      if (dd.eq.1) then
 
6135
         c_mu2_r=scales2(2,ic)
 
6136
         c_mu2_f=scales2(3,ic)
 
6137
      else
 
6138
c need to recompute the scales using the momenta
 
6139
         dynamical_scale_choice=dyn_scale(dd)
 
6140
         do i=1,nexternal
 
6141
            do j=0,3
 
6142
               pp(j,i)=momenta(j,i,ic)
 
6143
            enddo
 
6144
         enddo
 
6145
         call set_ren_scale(pp,muR)
 
6146
         c_mu2_r=muR**2
 
6147
         call set_fac_scale(pp,muF)
 
6148
         c_mu2_f=muF**2
 
6149
c     reset the default dynamical_scale_choice
 
6150
         dynamical_scale_choice=dyn_scale(1)
 
6151
      endif
 
6152
      return
 
6153
      end