← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

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

  • Committer: olivier Mattelaer
  • Date: 2015-03-05 00:14:16 UTC
  • mfrom: (258.1.9 2.3)
  • mto: (258.8.1 2.3)
  • mto: This revision was merged to the branch mainline in revision 259.
  • Revision ID: olivier.mattelaer@uclouvain.be-20150305001416-y9mzeykfzwnl9t0j
partial merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/reweight_xsec.f
2
2
c NLO and aMC@NLO results, for the computation of scale 
3
3
c and PDF uncertainties
4
4
 
5
 
      subroutine reweight_fillkin(pp,itype)
6
 
c Fills four-momenta common blocks, according to itype
7
 
      implicit none
8
 
      include "genps.inc"
9
 
      include "nexternal.inc"
10
 
      include "reweight.inc"
11
 
 
12
 
      double precision pp(0:3,nexternal)
13
 
      integer itype
14
 
 
15
 
      double precision p1_cnt(0:3,nexternal,-2:2)
16
 
      double precision wgt_cnt(-2:2)
17
 
      double precision pswgt_cnt(-2:2)
18
 
      double precision jac_cnt(-2:2)
19
 
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
20
 
 
21
 
      integer i,j,icnt,ione
22
 
      parameter (ione=1)
23
 
c
24
 
      if(itype.eq.1)then
25
 
        do i=1,nexternal
26
 
          do j=0,3
27
 
            wgtkin(j,i,ione)=pp(j,i)
28
 
          enddo
29
 
        enddo
30
 
      elseif(itype.ge.2.and.itype.le.4)then
31
 
        icnt=itype-2
32
 
        do i=1,nexternal
33
 
          do j=0,3
34
 
            wgtkin(j,i,itype)=p1_cnt(j,i,icnt)
35
 
          enddo
36
 
        enddo
37
 
      else
38
 
        write(*,*)'Error in reweight_fillkin: itype=',itype
39
 
        stop
40
 
      endif
41
 
      return
42
 
      end
43
 
 
44
 
 
45
 
      subroutine reweight_fillkin_all(pp,itype,iFKS)
46
 
c Fills four-momenta common blocks, according to itype
47
 
      implicit none
48
 
      include "genps.inc"
49
 
      include "nexternal.inc"
50
 
      include "nFKSconfigs.inc"
51
 
      include "reweight_all.inc"
52
 
 
53
 
      double precision pp(0:3,nexternal)
54
 
      integer itype,iFKS
55
 
 
56
 
      double precision p1_cnt(0:3,nexternal,-2:2)
57
 
      double precision wgt_cnt(-2:2)
58
 
      double precision pswgt_cnt(-2:2)
59
 
      double precision jac_cnt(-2:2)
60
 
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
61
 
 
62
 
      integer i,j,icnt,ione
63
 
      parameter (ione=1)
64
 
c
65
 
      if(itype.eq.1)then
66
 
        do i=1,nexternal
67
 
          do j=0,3
68
 
            wgtkin_all(j,i,ione,iFKS)=pp(j,i)
69
 
          enddo
70
 
        enddo
71
 
      elseif(itype.ge.2.and.itype.le.4)then
72
 
        icnt=2-itype
73
 
        do i=1,nexternal
74
 
          do j=0,3
75
 
            wgtkin_all(j,i,itype,iFKS)=p1_cnt(j,i,icnt)
76
 
          enddo
77
 
        enddo
78
 
      else
79
 
        write(*,*)'Error in reweight_fillkin: itype=',itype
80
 
        stop
81
 
      endif
82
 
      return
83
 
      end
84
 
 
85
 
 
86
 
      subroutine reweight_replkin_cnt()
87
 
c Fills wgtkin(*,*,2) with the content of wgtkin(*,*,3) or wgtkin(*,*,4).
88
 
c This is needed since only wgtkin(*,*,0) is used in the computations of
89
 
c cross sections through weights, and counterterms may be non zero
90
 
c only for collinear or soft-collinear configurations (eg because of
91
 
c choices of xicut and delta's). This is fine so long as the counterevent
92
 
c kinematic configurations are identical for the computations performed
93
 
c here 
94
 
      implicit none
95
 
      include "genps.inc"
96
 
      include "nexternal.inc"
97
 
      include "reweight.inc"
98
 
 
99
 
      integer i,j,icnt,itwo
100
 
      parameter (itwo=2)
101
 
c
102
 
      if (wgtkin(0,1,2).gt.0.d0 .or.
103
 
     #   (wgtkin(0,1,3).lt.0.d0.and.wgtkin(0,1,4).lt.0.d0) )then
104
 
        write(*,*)'Error in reweight_replkin_cnt'
105
 
        write(*,*)wgtkin(0,1,2),wgtkin(0,1,3),wgtkin(0,1,4)
106
 
        stop
107
 
      endif
108
 
c
109
 
      if(wgtkin(0,1,4).gt.0.d0)then
110
 
        icnt=4
111
 
      else
112
 
        icnt=3
113
 
      endif
114
 
c
115
 
      do i=1,nexternal
116
 
        do j=0,3
117
 
          wgtkin(j,i,itwo)=wgtkin(j,i,icnt)
118
 
        enddo
119
 
      enddo
120
 
c
121
 
      return
122
 
      end
123
 
 
124
 
 
125
 
      subroutine reweight_fill_extra()
126
 
c Fills arrays with dimensions maxparticles, equating them with their
127
 
c counterparts with dimensions nexternal. Needed before calling routines
128
 
c that only include reweight0.inc, such as read_lhef_event() and 
129
 
c write_lhef_event()
130
 
      implicit none
131
 
      include "genps.inc"
132
 
      include "nexternal.inc"
133
 
      include "reweight.inc"
134
 
      integer i,j,k
135
 
c
136
 
      do k=1,4
137
 
        do i=1,nexternal
138
 
          do j=0,3
139
 
            wgtkinE(j,i,k)=wgtkin(j,i,k)
140
 
          enddo
141
 
        enddo
142
 
      enddo
143
 
c
144
 
      do k=1,nexternal
145
 
        wgtwmcxsecE(k)=wgtwmcxsec(k)
146
 
        wgtmcxbjE(1,k)=wgtmcxbj(1,k)
147
 
        wgtmcxbjE(2,k)=wgtmcxbj(2,k)
148
 
      enddo
149
 
c
150
 
      return
151
 
      end
152
 
 
153
 
 
 
5
      
154
6
      subroutine reweight_fill_extra_inverse()
155
7
c Fills arrays with dimensions nexternal, equating them with their
156
8
c counterparts with dimensions maxparticles; this is thus the inverse
316
168
      subroutine fill_reweight0inc_nbody(iproc)
317
169
c Set all reweight variables equal to zero
318
170
      implicit none
 
171
      include "nexternal.inc"
319
172
      include "reweight_all.inc"
320
173
      integer iproc
321
174
      logical debug
464
317
      return
465
318
      end
466
319
 
467
 
 
468
 
      subroutine sum_reweight(iFKS_s,iFKS,iproc)
469
 
c Set all reweight variables equal to zero
470
 
      implicit none
471
 
      include "genps.inc"
472
 
      include "nexternal.inc"
473
 
      include "nFKSconfigs.inc"
474
 
      include "reweight_all.inc"
475
 
      integer i,j,k,iFKS_s,iFKS,iproc
476
 
      logical debug
477
 
      parameter (debug=.false.)
478
 
c
479
 
      if (debug) write (*,*) 'wgtref',iFKS_s,iFKS,wgtref_all(iFKS_s
480
 
     $     ,iproc),wgtref_all(iFKS,iproc)
481
 
      call reweight_sum(wgtref_all(iFKS_s,iproc),wgtref_all(iFKS,iproc))
482
 
      if (debug) write (*,*) 'wgtref_nbody',iFKS_s,iFKS
483
 
     $     ,wgtref_nbody_all(iproc),wgtref_nbody_all(iproc)
484
 
      call reweight_sum(wgtref_nbody_all(iproc),wgtref_nbody_all(iproc))
485
 
      do k=1,4
486
 
         if (debug) write (*,*) 'wgtqes2',k,iFKS_s,iFKS,wgtqes2_all(k
487
 
     &        ,iFKS_s),wgtqes2_all(k,iFKS)
488
 
         call reweight_check_equal(wgtqes2_all(k,iFKS_s),wgtqes2_all(k
489
 
     &        ,iFKS))
490
 
         if (debug) write (*,*) 'wgtxbj1',k,iFKS_s,iFKS,wgtxbj_all(1,k
491
 
     &        ,iFKS_s),wgtxbj_all(1,k,iFKS)
492
 
         if (k.eq.1.or.k.eq.2) then
493
 
            call reweight_check_equal(wgtxbj_all(1,k,iFKS_s)
494
 
     &           ,wgtxbj_all(1,k,iFKS))
495
 
         else
496
 
            call reweight_overwrite(wgtxbj_all(1,k,iFKS_s)
497
 
     &           ,wgtxbj_all(1,k,iFKS),2)
498
 
         endif
499
 
         if (debug) write (*,*) 'wgtxbj2',k,iFKS_s,iFKS,wgtxbj_all(2,k
500
 
     &        ,iFKS_s),wgtxbj_all(2,k,iFKS)
501
 
         if (k.eq.1.or.k.eq.2) then
502
 
            call reweight_check_equal(wgtxbj_all(2,k,iFKS_s)
503
 
     &           ,wgtxbj_all(2,k,iFKS))
504
 
         else
505
 
            call reweight_overwrite(wgtxbj_all(2,k,iFKS_s)
506
 
     &           ,wgtxbj_all(2,k,iFKS),2)
507
 
         endif
508
 
        do i=1,nexternal
509
 
          do j=0,3
510
 
             if (debug) write (*,*) 'wgtkin',j,i,k,iFKS_s,iFKS
511
 
     &            ,wgtkin_all(j,i,k,iFKS_s),wgtkin_all(j,i,k,iFKS)
512
 
             if (k.eq.1.or.k.eq.2) then
513
 
                call reweight_check_equal(wgtkin_all(j,i,k,iFKS_s)
514
 
     &               ,wgtkin_all(j,i,k,iFKS))
515
 
             else
516
 
                call reweight_overwrite(wgtkin_all(j,i,k,iFKS_s)
517
 
     &               ,wgtkin_all(j,i,k,iFKS),2)
518
 
             endif
519
 
          enddo
520
 
        enddo
521
 
        if (debug) write (*,*) 'wgtmuR2',k,iFKS_s,iFKS,wgtmuR2_all(k
522
 
     &       ,iFKS_s),wgtmuR2_all(k,iFKS)
523
 
        call reweight_check_equal(wgtmuR2_all(k,iFKS_s),wgtmuR2_all(k
524
 
     &       ,iFKS))
525
 
        if (debug) write (*,*) 'wgtmuF12',k,iFKS_s,iFKS,wgtmuF12_all(k
526
 
     &       ,iFKS_s),wgtmuF12_all(k,iFKS)
527
 
        call reweight_check_equal(wgtmuF12_all(k,iFKS_s),wgtmuF12_all(k
528
 
     &       ,iFKS))
529
 
        if (debug) write (*,*) 'wgtmuF22',k,iFKS_s,iFKS,wgtmuF22_all(k
530
 
     &       ,iFKS_s),wgtmuF22_all(k,iFKS)
531
 
        call reweight_check_equal(wgtmuF22_all(k,iFKS_s),wgtmuF22_all(k
532
 
     &       ,iFKS))
533
 
        if (debug) write (*,*) 'wgtwreal',k,iFKS_s,iFKS,wgtwreal_all(k
534
 
     &       ,iFKS_s),wgtwreal_all(k,iFKS)
535
 
        call reweight_sum(wgtwreal_all(k,iFKS_s),wgtwreal_all(k,iFKS))
536
 
        if (debug) write (*,*) 'wgtwdeg',k,iFKS_s,iFKS,wgtwdeg_all(k
537
 
     &       ,iFKS_s),wgtwdeg_all(k,iFKS)
538
 
        call reweight_sum(wgtwdeg_all(k,iFKS),wgtwdeg_all(k,iFKS))
539
 
        if (debug) write (*,*) 'wgtwdegmuF',k,iFKS_s,iFKS
540
 
     &       ,wgtwdegmuf_all(k,iFKS_s),wgtwdegmuf_all(k,iFKS)
541
 
        call reweight_sum(wgtwdegmuf_all(k,iFKS),wgtwdegmuf_all(k,iFKS))
542
 
c$$$        if(k.eq.2)then
543
 
c$$$           if (debug) write (*,*) 'wgtwborn',k,wgtwborn(k),wgtwborn_all
544
 
c$$$           call reweight_overwrite(wgtwborn(k),wgtwborn_all,0)
545
 
c$$$           if (debug) write (*,*) 'wgtwns',k,wgtwns(k),wgtwns_all
546
 
c$$$           call reweight_overwrite(wgtwns(k),wgtwns_all,0)
547
 
c$$$           if (debug) write (*,*) 'wgtwnsmuf',k,wgtwnsmuf(k)
548
 
c$$$     &          ,wgtwnsmuf_all
549
 
c$$$           call reweight_overwrite(wgtwnsmuf(k),wgtwnsmuf_all,0)
550
 
c$$$           if (debug) write (*,*) 'wgtwnsmur',k,wgtwnsmur(k)
551
 
c$$$     &          ,wgtwnsmur_all
552
 
c$$$          call reweight_overwrite(wgtwnsmur(k),wgtwnsmur_all,0)
553
 
c$$$        endif
554
 
      enddo
555
 
      do k=1,nexternal
556
 
         if (debug) write (*,*) 'wgtwmcxsec',k,iFKS_s,iFKS
557
 
     &        ,wgtwmcxsec_all(k,iFKS_s),wgtwmcxsec_all(k,iFKS)
558
 
         call reweight_sum(wgtwmcxsec_all(k,iFKS_s),wgtwmcxsec_all(k
559
 
     &        ,iFKS))
560
 
         if (debug) write (*,*) 'wgtwmcxbj1',k,iFKS_s,iFKS
561
 
     &        ,wgtmcxbj_all(1,k,iFKS_s),wgtmcxbj_all(1,k,iFKS)
562
 
         call reweight_check_equal(wgtmcxbj_all(1,k,iFKS_s)
563
 
     &        ,wgtmcxbj_all(1,k,iFKS))
564
 
         if (debug) write (*,*) 'wgtwmcxbj2',k,iFKS_s,iFKS
565
 
     &        ,wgtmcxbj_all(2,k,iFKS_s),wgtmcxbj_all(2,k,iFKS)
566
 
         call reweight_check_equal(wgtmcxbj_all(2,k,iFKS_s)
567
 
     &        ,wgtmcxbj_all(2,k,iFKS))
568
 
      enddo
569
 
      if (debug) write (*,*) 'iwgtnumpartn',iFKS
570
 
     &     ,iwgtnumpartn_all(iFKS_s),iwgtnumpartn_all(iFKS)
571
 
      call reweight_check_equal_int(iwgtnumpartn_all(iFKS_s)
572
 
     &     ,iwgtnumpartn_all(iFKS))
573
 
      return
574
 
      end
575
 
 
576
 
      subroutine reweight_check_equal_int(a,b)
577
 
      implicit none
578
 
      integer a,b
579
 
      if (a.ne.b) then
580
 
         write (*,*) 'Error #1 in reweight_check_equal_int',a,b
581
 
         stop
582
 
      endif
583
 
      return
584
 
      end
585
 
 
586
 
      subroutine reweight_check_equal(a,b)
587
 
      implicit none
588
 
      double precision a,b
589
 
      double precision vtiny
590
 
      parameter (vtiny=1d-10)
591
 
      if (a.ne.0) then
592
 
         if (abs(a-b)/abs(a).gt.vtiny) then
593
 
            write (*,*) 'Error #1 in reweight_check_equal',a,b
594
 
            stop
595
 
         endif
596
 
      elseif (b.ne.0) then
597
 
         if (abs(a-b)/abs(b).gt.vtiny) then
598
 
            write (*,*) 'Error #2 in reweight_check_equal',a,b
599
 
            stop
600
 
         endif
601
 
      endif
602
 
      return
603
 
      end
604
 
 
605
 
      subroutine reweight_sum(a,b)
606
 
      implicit none
607
 
      double precision a,b
608
 
      a=a+b
609
 
      return
610
 
      end
611
 
 
612
 
 
613
 
      function compute_rwgt_wgt_NLO(xmuR_over_ref,xmuF1_over_ref,
614
 
     #                              xmuF2_over_ref,xQES_over_ref,
615
 
     #                              kwgtinfo)
616
 
c Recomputes the NLO cross section using the weights saved, and compares
617
 
c with the reference weight
618
 
      implicit none
619
 
      include "genps.inc"
620
 
      include "nexternal.inc"
621
 
      include 'coupl.inc'
622
 
      include 'q_es.inc'
623
 
      include 'run.inc'
624
 
      include "reweight.inc"
625
 
      include "appl_common.inc"
626
 
      include "nFKSconfigs.inc"
627
 
 
628
 
      double precision compute_rwgt_wgt_NLO
629
 
      double precision xmuR_over_ref,xmuF1_over_ref,
630
 
     #                 xmuF2_over_ref,xQES_over_ref
631
 
      integer kwgtinfo
632
 
      double precision rwgt,xsec,xlum,dlum,xlgmuf,xlgmur,alphas
633
 
      double precision xsec11,xsec12,xsec20
634
 
      double precision QES2_local
635
 
      double precision save_murrat,save_muf1rat,save_muf2rat,save_qesrat
636
 
      double precision tiny,pi
637
 
      parameter (tiny=1.d-2)
638
 
      parameter (pi=3.14159265358979323846d0)
639
 
 
640
 
      integer k,izero,mohdr
641
 
      parameter (izero=0)
642
 
      parameter (mohdr=-100)
643
 
      INTEGER NFKSPROCESS
644
 
      COMMON/C_NFKSPROCESS/NFKSPROCESS
645
 
      integer save_nFKSprocess
646
 
c FxFx merging
647
 
      logical rewgt_mohdr_calculated,rewgt_izero_calculated
648
 
      double precision rewgt_mohdr,rewgt_izero,rewgt_exp_mohdr
649
 
     $     ,rewgt_exp_izero
650
 
      logical setclscales
651
 
      double precision rewgt
652
 
      external setclscales,rewgt
653
 
      double precision rwgt_muR_dep_fac
654
 
c APPLGRID
655
 
      integer iappl
656
 
      common /for_applgrid/ iappl
657
 
      double precision vegas_weight
658
 
      common/cvegas_weight/vegas_weight
659
 
      integer flavour_map(fks_configs)
660
 
      common/c_flavour_map/flavour_map
661
 
      double precision final_state_rescaling
662
 
      integer iproc_save(fks_configs),eto(maxproc,fks_configs),
663
 
     1        etoi(maxproc,fks_configs),maxproc_found
664
 
      common/cproc_combination/iproc_save,eto,etoi,maxproc_found
665
 
c
666
 
      save_murrat=muR_over_ref
667
 
      save_muf1rat=muF1_over_ref
668
 
      save_muf2rat=muF2_over_ref
669
 
      save_qesrat=QES_over_ref
670
 
      save_nFKSprocess=nFKSprocess
671
 
c
672
 
      muR_over_ref=xmuR_over_ref 
673
 
      muF1_over_ref=xmuF1_over_ref
674
 
      muF2_over_ref=xmuF2_over_ref
675
 
      QES_over_ref=xQES_over_ref 
676
 
c
677
 
      if(kwgtinfo.ne.4.and.kwgtinfo.ne.5) wgtbpower=rwgtbpower
678
 
c
679
 
      xsec=0.d0
680
 
      xsec11=0.d0
681
 
      xsec12=0.d0
682
 
      xsec20=0.d0
683
 
 
684
 
      if (wgtwreal(1).eq.0d0) goto 541
685
 
 
686
 
      call set_cms_stuff(mohdr)
687
 
      if( (kwgtinfo.eq.1.and.wgtmuR2(1).ne.0.d0) .or.
688
 
     $     ((kwgtinfo.ge.3.or.kwgtinfo.le.5).and.wgtkin(0,1,1).gt.0.d0)
689
 
     $     )then
690
 
        if(kwgtinfo.eq.1)then
691
 
          scale=muR_over_ref*sqrt(wgtmuR2(1))
692
 
          g=sqrt(4d0*pi*alphas(scale))
693
 
          q2fact(1)=wgtmuF12(1) * muF1_over_ref**2
694
 
          q2fact(2)=wgtmuF22(1) * muF2_over_ref**2
695
 
c Should cause the code to crash if used
696
 
          QES2=0.d0
697
 
        elseif(kwgtinfo.ge.3.or.kwgtinfo.le.5)then
698
 
          call set_cms_stuff(mohdr)
699
 
           if (ickkw.eq.3) then
700
 
              mu_r=sqrt(wgtmuR2(1))*muR_over_ref
701
 
              scale=mu_r
702
 
              g=sqrt(4d0*pi*alphas(scale))
703
 
              call update_as_param()
704
 
              q2fact(1)=muF1_over_ref**2*wgtmuF12(1)
705
 
              q2fact(2)=muF2_over_ref**2*wgtmuF22(1)
706
 
           else
707
 
              call set_alphaS(wgtkin(0,1,1))
708
 
           endif
709
 
        else
710
 
          write(*,*)'Error #0a in compute_rwgt_wgt_NLO',kwgtinfo
711
 
          stop
712
 
        endif
713
 
        xbk(1) = wgtxbj(1,1)
714
 
        xbk(2) = wgtxbj(2,1)
715
 
        if(xbk(1).le.0.d0.or.xbk(2).le.0.d0.or.
716
 
     #     xbk(1).gt.1.d0.or.xbk(2).gt.1.d0)then
717
 
          if(wgtwreal(1).ne.0.d0)then
718
 
            write(*,*)'Error #1 in compute_rwgt_wgt_NLO'
719
 
            write(*,*)xbk(1),xbk(2),wgtwreal(1)
720
 
            stop
721
 
          endif
722
 
        else
723
 
          nFKSprocess=nFKSprocess_used
724
 
          xlum = dlum()
725
 
          xsec11=xsec11+xlum*wgtwreal(1)*g**(2*wgtbpower+2.d0)
726
 
     f         * rwgt_muR_dep_fac(scale)
727
 
com-- muR-dependent fac is reweighted here
728
 
        endif
729
 
      endif
730
 
c
731
 
 541  continue
732
 
      if ( wgtwreal(2).eq.0d0 .and.
733
 
     $     wgtwreal(3).eq.0d0 .and. wgtwreal(4).eq.0d0 .and.
734
 
     $     wgtwdeg(2).eq.0d0 .and.
735
 
     $     wgtwdeg(3).eq.0d0 .and. wgtwdeg(4).eq.0d0 .and.
736
 
     $     wgtwdegmuf(2).eq.0d0 .and.
737
 
     $     wgtwdegmuf(3).eq.0d0 .and. wgtwdegmuf(4).eq.0d0 .and.
738
 
     $     wgtwborn(2).eq.0d0 .and. wgtwns(2).eq.0d0 .and.
739
 
     $     wgtwnsmuf(2).eq.0d0 .and. wgtwnsmur(2).eq.0d0) goto 542
740
 
 
741
 
      call set_cms_stuff(izero)
742
 
      if( (kwgtinfo.eq.1.and.wgtmuR2(2).ne.0.d0) .or.
743
 
     $     ((kwgtinfo.ge.3.or.kwgtinfo.le.5).and.wgtkin(0,1,2).gt.0.d0)
744
 
     $     )then
745
 
        if(kwgtinfo.eq.1)then
746
 
          scale=muR_over_ref*sqrt(wgtmuR2(2))
747
 
          g=sqrt(4d0*pi*alphas(scale))
748
 
          q2fact(1)=wgtmuF12(2) * muF1_over_ref**2
749
 
          q2fact(2)=wgtmuF22(2) * muF2_over_ref**2
750
 
c Should cause the code to crash if used
751
 
          QES2=0.d0
752
 
        elseif(kwgtinfo.ge.3.or.kwgtinfo.le.5)then
753
 
           if (ickkw.eq.3) then
754
 
              mu_r=sqrt(wgtmuR2(2))*muR_over_ref
755
 
              scale=mu_r
756
 
              g=sqrt(4d0*pi*alphas(scale))
757
 
              call update_as_param()
758
 
              q2fact(1)=muF1_over_ref**2*wgtmuF12(2)
759
 
              q2fact(2)=muF2_over_ref**2*wgtmuF22(2)
760
 
           else
761
 
              call set_alphaS(wgtkin(0,1,2))
762
 
           endif
763
 
        else
764
 
          write(*,*)'Error #0b in compute_rwgt_wgt_NLO',kwgtinfo
765
 
          stop
766
 
        endif
767
 
        QES2_local=wgtqes2(2)
768
 
        if(QES2_local.eq.0.d0)then
769
 
          if(wgtwdegmuf(3).ne.0.d0.or.
770
 
     #       wgtwdegmuf(4).ne.0.d0.or.
771
 
     #       wgtwnsmuf(2).ne.0.d0.or.
772
 
     #       wgtwnsmur(2).ne.0.d0)then
773
 
            write(*,*)'Error in compute_rwgt_wgt_NLO'
774
 
            write(*,*)' Ellis-Sexton scale was not set'
775
 
            write(*,*)wgtwdegmuf(3),wgtwdegmuf(4),
776
 
     #                wgtwnsmuf(2),wgtwnsmur(2)
777
 
            stop
778
 
          endif
779
 
          xlgmuf=0.d0
780
 
          xlgmur=0.d0
781
 
        else
782
 
           if(abs(QES2/QES2_local-1.d0).gt.tiny.and.
783
 
     &          (kwgtinfo.ge.3.or.kwgtinfo.le.5))then
784
 
              write(*,*)'Error in compute_rwgt_wgt_NLO'
785
 
              write(*,*)' Mismatch in ES scale',QES2,QES2_local
786
 
              stop
787
 
           endif
788
 
          xlgmuf=log(q2fact(1)/QES2_local)
789
 
          xlgmur=log(scale**2/QES2_local)
790
 
        endif
791
 
        do k=2,4
792
 
          xbk(1) = wgtxbj(1,k)
793
 
          xbk(2) = wgtxbj(2,k)
794
 
          if(xbk(1).le.0.d0.or.xbk(2).le.0.d0.or.
795
 
     #       xbk(1).gt.1.d0.or.xbk(2).gt.1.d0)then
796
 
            if(wgtwreal(k).ne.0.d0.or.
797
 
     #         wgtwdeg(k).ne.0.d0.or.
798
 
     #         wgtwdegmuf(k).ne.0.d0.or.
799
 
     #         (k.eq.2.and.(wgtwborn(2).ne.0.d0.or.
800
 
     #                      wgtwns(2).ne.0.d0.or.
801
 
     #                      wgtwnsmuf(2).ne.0.d0.or.
802
 
     #                      wgtwnsmur(2).ne.0.d0)))then
803
 
              write(*,*)'Error #2 in compute_rwgt_wgt_NLO'
804
 
              write(*,*)k,xbk(1),xbk(2)
805
 
              write(*,*)wgtwreal(k),wgtwdeg(k),wgtwdegmuf(k)
806
 
              if(k.eq.2)write(*,*)wgtwborn(k),wgtwns(k),
807
 
     #                            wgtwnsmuf(k),wgtwnsmur(k)
808
 
              stop
809
 
            endif
810
 
          else
811
 
            nFKSprocess=nFKSprocess_used
812
 
            xlum = dlum()
813
 
            xsec12=xsec12+xlum*( wgtwreal(k)+wgtwdeg(k)+
814
 
     #                           wgtwdegmuf(k)*xlgmuf )*
815
 
     #                g**(2*wgtbpower+2.d0) * rwgt_muR_dep_fac(scale)
816
 
com-- muR-dependent fac is reweighted here
817
 
            if(k.eq.2)then
818
 
               nFKSprocess=nFKSprocess_used_born
819
 
               xlum = dlum()
820
 
              if(wgtbpower.gt.0)then
821
 
                xsec20=xsec20+xlum*wgtwborn(k)*g**(2*wgtbpower)
822
 
     f                * rwgt_muR_dep_fac(scale)
823
 
com-- muR-dependent fac is reweighted here
824
 
              else
825
 
                xsec20=xsec20+xlum*wgtwborn(k) * rwgt_muR_dep_fac(scale)
826
 
com-- muR-dependent fac is reweighted here
827
 
              endif
828
 
              xsec12=xsec12+xsec20
829
 
              xsec12=xsec12+xlum*( wgtwns(k)+
830
 
     #                         wgtwnsmuf(k)*xlgmuf+
831
 
     #                         wgtwnsmur(k)*xlgmur )*
832
 
     #                  g**(2*wgtbpower+2.d0) * rwgt_muR_dep_fac(scale)
833
 
com-- muR-dependent fac is reweighted here
834
 
            endif
835
 
          endif
836
 
        enddo
837
 
      endif
838
 
 
839
 
 542  continue
840
 
 
841
 
c
842
 
      muR_over_ref=save_murrat
843
 
      muF1_over_ref=save_muf1rat
844
 
      muF2_over_ref=save_muf2rat
845
 
      QES_over_ref=save_qesrat
846
 
      nFKSprocess=save_nFKSprocess
847
 
c
848
 
      wgtNLO11=xsec11
849
 
      wgtNLO12=xsec12-xsec20
850
 
      wgtNLO20=xsec20
851
 
      xsec=xsec11+xsec12
852
 
      compute_rwgt_wgt_NLO=xsec
853
 
c
854
 
 
855
 
************************************************************************
856
 
*     APPLgrid (Refer to arXiv:1110.4738).
857
 
************************************************************************
858
 
      if(iappl.ne.0)then
859
 
         do k=1,4
860
 
c     Bjorken-x
861
 
            appl_x1(k) = wgtxbj(1,k)
862
 
            appl_x2(k) = wgtxbj(2,k)
863
 
c     Scales
864
 
            if(k.eq.1)then
865
 
               appl_muF2(k) = wgtmuF12(1) * muF1_over_ref**2
866
 
               appl_muR2(k) = wgtmuR2(1)  * muR_over_ref**2
867
 
               appl_QES2(k) = wgtqes2(2)
868
 
            elseif(k.ge.2.and.k.le.4)then
869
 
               appl_muF2(k) = wgtmuF12(2) * muF1_over_ref**2
870
 
               appl_muR2(k) = wgtmuR2(2)  * muR_over_ref**2
871
 
               appl_QES2(k) = wgtqes2(2)
872
 
            endif
873
 
c     Initialize weights
874
 
            appl_w0(k) = 0d0
875
 
            appl_wR(k) = 0d0
876
 
            appl_wF(k) = 0d0
877
 
            appl_wB(k) = 0d0
878
 
c     Fill weights (Same notation as in Eq. (2.17) of the paper above).
879
 
c     Flavor map: integer from 1 to nproc, with nproc the number of the
880
 
c     parton luminosities, defined at generation in
881
 
c     initial_states_map.dat.
882
 
            if(k.eq.1)then
883
 
c     Events (only W0 is non-zero)
884
 
               if(wgtwreal(k).ne.0d0)then
885
 
c     Weights
886
 
                  appl_w0(k) = wgtwreal(k)
887
 
c     Flavour map
888
 
                  appl_flavmap(k) = flavour_map(nFKSprocess_used)
889
 
               endif
890
 
            else
891
 
               if ( wgtwreal(2).ne.0d0 .or. wgtwreal(3).ne.0d0 .or.
892
 
     $              wgtwreal(4).ne.0d0 .or. wgtwdeg(2).ne.0d0 .or.
893
 
     $              wgtwdeg(3).ne.0d0 .or. wgtwdeg(4).ne.0d0 .or.
894
 
     $              wgtwdegmuf(2).ne.0d0 .or. wgtwdegmuf(3).ne.0d0 .or.
895
 
     $              wgtwdegmuf(4).ne.0d0 .or. wgtwborn(2).ne.0d0 .or.
896
 
     $              wgtwns(2).ne.0d0 .or. wgtwnsmuf(2).ne.0d0 .or.
897
 
     $              wgtwnsmur(2).ne.0d0) then
898
 
c     Soft events.
899
 
c     Note that for k=2 it is the Born flavour map that is always used.
900
 
c     So there is a unique set of weights for k=2, all of which uses
901
 
c     flavour_map(nFKSprocess_used_Born)
902
 
                  if(k.eq.2)then
903
 
c     Weights
904
 
                     appl_w0(k) = wgtwreal(k) + wgtwdeg(k) + wgtwns(k)
905
 
                     appl_wR(k) = wgtwnsmur(k)
906
 
                     appl_wF(k) = wgtwdegmuf(k) + wgtwnsmuf(k)
907
 
                     appl_wB(k) = wgtwborn(k)
908
 
c     Flavour map
909
 
                     appl_flavmap(k)=flavour_map(nFKSprocess_used_Born)
910
 
c     Collinear and sof-collinear events
911
 
                  elseif(k.eq.3.or.k.eq.4)then
912
 
c     Weights
913
 
                     appl_w0(k) = wgtwreal(k) + wgtwdeg(k)
914
 
                     appl_wF(k) = wgtwdegmuf(k)
915
 
c     Flavour map
916
 
                     appl_flavmap(k) = flavour_map(nFKSprocess_used)
917
 
                  endif
918
 
               endif
919
 
            endif
920
 
c     Multiply the weights by the number of equivalent final states that
921
 
c     share both the same matrix element and the same initial luminosity
922
 
            final_state_rescaling = dble(iproc_save(nFKSprocess_used))
923
 
     1           / dble(appl_nproc(flavour_map(nFKSprocess_used)))
924
 
c     Multiply by Vegas weight before filling the APPLgrid histograms,
925
 
c     as done for the NLO plots in the reweighting routine, and by the
926
 
c     final rescaling factor.
927
 
            appl_w0(k) = appl_w0(k) * final_state_rescaling * vegas_weight 
928
 
            appl_wR(k) = appl_wR(k) * final_state_rescaling * vegas_weight
929
 
            appl_wF(k) = appl_wF(k) * final_state_rescaling * vegas_weight
930
 
            appl_wB(k) = appl_wB(k) * final_state_rescaling * vegas_weight
931
 
         enddo
932
 
c     Save the value of the event weight (hadronic differential cross section).
933
 
c     Computed with the reference PDF sets for checks of the amcatnlo interface.
934
 
         appl_event_weight = compute_rwgt_wgt_NLO
935
 
c     Save also the vegas weight
936
 
         appl_vegaswgt = vegas_weight
937
 
      endif
938
 
      return
939
 
      end
940
 
 
941
 
 
942
320
      function compute_rwgt_wgt_Hev(xmuR_over_ref,xmuF1_over_ref,
943
321
     #                              xmuF2_over_ref,xQES_over_ref,
944
322
     #                              kwgtinfo)
1625
1003
      end
1626
1004
 
1627
1005
 
1628
 
      subroutine check_rwgt_wgt(idstring)
1629
 
      implicit none
1630
 
      include "genps.inc"
1631
 
      include "nexternal.inc"
1632
 
      include 'run.inc'
1633
 
      include "reweight.inc"
1634
 
      character*3 idstring
1635
 
      double precision compute_rwgt_wgt_NLO,compute_rwgt_wgt_Hev,
1636
 
     #                 compute_rwgt_wgt_Sev,compute_rwgt_wgt_Sev_nbody
1637
 
      double precision wgtnew,tiny
1638
 
      parameter (tiny=1.d-2)
1639
 
c
1640
 
      INTEGER              IPROC
1641
 
      DOUBLE PRECISION PD(0:MAXPROC)
1642
 
      COMMON /SUBPROC/ PD, IPROC
1643
 
      integer iproc_save
1644
 
      save iproc_save
1645
 
c
1646
 
      iproc_save=iproc
1647
 
      if(idstring.eq."NLO")then
1648
 
        wgtnew=compute_rwgt_wgt_NLO(muR_over_ref,muF1_over_ref,
1649
 
     #                              muF2_over_ref,QES_over_ref,
1650
 
     #                              iwgtinfo)
1651
 
      elseif(idstring.eq."Hev")then
1652
 
        wgtnew=compute_rwgt_wgt_Hev(muR_over_ref,muF1_over_ref,
1653
 
     #                              muF2_over_ref,QES_over_ref,
1654
 
     #                              iwgtinfo)
1655
 
      elseif(idstring.eq."Sev")then
1656
 
        wgtnew=compute_rwgt_wgt_Sev(muR_over_ref,muF1_over_ref,
1657
 
     #                              muF2_over_ref,QES_over_ref,
1658
 
     #                              iwgtinfo)
1659
 
      elseif(idstring.eq."nbd")then
1660
 
        wgtnew=compute_rwgt_wgt_Sev_nbody(muR_over_ref,muF1_over_ref,
1661
 
     #                              muF2_over_ref,QES_over_ref,
1662
 
     #                              iwgtinfo)
1663
 
      else
1664
 
        write(*,*)'Error in check_rwgt_wgt'
1665
 
        write(*,*)' Unknown function: ',idstring
1666
 
        stop
1667
 
      endif
1668
 
c
1669
 
      if (idstring.eq."nbd") then
1670
 
      if( (abs(wgtref_nbody).ge.1.d0 .and.
1671
 
     #     abs(1.d0-wgtnew/wgtref_nbody).gt.tiny) .or.
1672
 
     #    (abs(wgtref_nbody).lt.1.d0 .and.
1673
 
     #     abs(wgtnew-wgtref_nbody).gt.tiny) )then
1674
 
         write(*,*)'Error in check_rwgt_wgt: ',idstring,wgtref_nbody
1675
 
     &        ,wgtnew
1676
 
        stop
1677
 
      endif
1678
 
      else
1679
 
      if( (abs(wgtref).ge.1.d0 .and.
1680
 
     #     abs(1.d0-wgtnew/wgtref).gt.tiny) .or.
1681
 
     #    (abs(wgtref).lt.1.d0 .and.
1682
 
     #     abs(wgtnew-wgtref).gt.tiny) )then
1683
 
        write(*,*)'Error in check_rwgt_wgt: ',idstring,wgtref,wgtnew
1684
 
        stop
1685
 
      endif
1686
 
      endif
1687
 
      iproc=iproc_save
1688
 
c
1689
 
      return
1690
 
      end
1691
 
 
1692
 
 
1693
 
      subroutine fill_rwgt_NLOplot()
1694
 
      implicit none
1695
 
      include "genps.inc"
1696
 
      include "nexternal.inc"
1697
 
      include 'coupl.inc'
1698
 
      include 'run.inc'
1699
 
      include "reweight.inc"
1700
 
      include "reweightNLO.inc"
1701
 
      integer izero
1702
 
      parameter (izero=0)
1703
 
      integer kr,kf,n
1704
 
      double precision pr_muR_over_ref,pr_muF1_over_ref,
1705
 
     # pr_muF2_over_ref,dummy,compute_rwgt_wgt_NLO
1706
 
c
1707
 
      dummy=compute_rwgt_wgt_NLO(ymuR_over_ref,ymuF1_over_ref,
1708
 
     #                           ymuF2_over_ref,yQES_over_ref,
1709
 
     #                           iwgtinfo)
1710
 
      wgtrefNLO11=wgtNLO11
1711
 
      wgtrefNLO12=wgtNLO12
1712
 
      wgtrefNLO20=wgtNLO20
1713
 
c
1714
 
      if(do_rwgt_scale)then
1715
 
        do kr=1,numscales
1716
 
          do kf=1,numscales
1717
 
            pr_muR_over_ref=ymuR_over_ref*yfactR(kr)
1718
 
            pr_muF1_over_ref=ymuF1_over_ref*yfactF(kf)
1719
 
            pr_muF2_over_ref=pr_muF1_over_ref
1720
 
            dummy=compute_rwgt_wgt_NLO(pr_muR_over_ref,pr_muF1_over_ref,
1721
 
     #                                 pr_muF2_over_ref,yQES_over_ref,
1722
 
     #                                 iwgtinfo)
1723
 
            wgtNLOxsecmu(1,kr,kf)=wgtNLO11
1724
 
            wgtNLOxsecmu(2,kr,kf)=wgtNLO12
1725
 
            wgtNLOxsecmu(3,kr,kf)=wgtNLO20
1726
 
          enddo
1727
 
        enddo
1728
 
      endif
1729
 
c
1730
 
      if(do_rwgt_pdf)then
1731
 
        do n=0,numPDFs-1
1732
 
          call InitPDF(n)
1733
 
          dummy=compute_rwgt_wgt_NLO(ymuR_over_ref,ymuF1_over_ref,
1734
 
     #                               ymuF2_over_ref,yQES_over_ref,
1735
 
     #                               iwgtinfo)
1736
 
          wgtNLOxsecPDF(1,n)=wgtNLO11
1737
 
          wgtNLOxsecPDF(2,n)=wgtNLO12
1738
 
          wgtNLOxsecPDF(3,n)=wgtNLO20
1739
 
        enddo
1740
 
c Restore default PDFs
1741
 
        call InitPDF(izero)
1742
 
      endif
1743
 
c
1744
 
      return
1745
 
      end
1746
 
 
1747
1006
 
1748
1007
      subroutine setup_fill_rwgt_NLOplot()
1749
1008
      implicit none
1811
1070
      implicit none
1812
1071
      double precision scale,vev,mbmb,apimuR,apimZ,apimb,mbmuR,alphas,pi
1813
1072
      parameter (pi=3.14159265358979323846d0)
 
1073
      include "nexternal.inc"
1814
1074
      include "genps.inc"
1815
1075
      include "reweight.inc"
1816
1076
      include "coupl.inc"