← Back to branch summary

~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3

« back to all changes in this revision

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

  • Committer: Marco Zaro
  • Date: 2016-05-24 12:48:05 UTC
  • mfrom: (78.337.150 2.4.1)
  • Revision ID: marco.zaro@gmail.com-20160524124805-d0y3ij808wdroloh
merged with 2.4.1 rev 410

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/reweight_xsec_events.f
6
6
c and PDF variations
7
7
c Compile with makefile_rwgt
8
8
      implicit none
9
 
      include "nexternal.inc"
10
 
      include "genps.inc"
11
 
      include "nFKSconfigs.inc"
 
9
      include "run.inc"
12
10
      include "reweight_all.inc"
13
 
      include "run.inc"
14
 
      character*7 pdlabel,epa_label
15
 
      integer lhaid
 
11
      integer i,ii,jj,kk,isave,idpdf(0:maxPDFs),itmp,lef,ifile,maxevt
 
12
     $     ,iSorH_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
 
13
     $     ,kwgtinfo,kexternal,jwgtnumpartn,ofile,kf,kr,n,nn
 
14
      double precision yfactR(maxscales),yfactF(maxscales),value(20)
 
15
     $     ,scale1_lhe,scale2_lhe,wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
 
16
     $     ,wgtpdfmax,saved_weight,xsecPDFr_acc(0:maxPDFs,maxPDFsets)
 
17
     $     ,xsecScale_acc(maxscales,maxscales,maxdynscales)
 
18
      logical AddInfoLHE,unweighted
 
19
      character*9 ch1
 
20
      character*10 MonteCarlo
 
21
      character*20 parm(20)
 
22
      character*80 event_file,fname1
 
23
      character*140 buff
 
24
c Parameters
 
25
      integer    izero
 
26
      parameter (izero=0)
 
27
c Common blocks
 
28
      character*7         pdlabel,epa_label
 
29
      integer       lhaid
16
30
      common/to_pdf/lhaid,pdlabel,epa_label
17
 
      integer maxevt,ifile,ofile,i,jj,isave,ii
18
 
      double precision saved_weight
19
 
      logical unweighted
 
31
c Les Houches Event File info:
20
32
      integer IDBMUP(2),PDFGUP(2),PDFSUP(2),IDWTUP,NPRUP,LPRUP
21
33
      double precision EBMUP(2),XSECUP,XERRUP,XMAXUP
22
34
      INTEGER MAXNUP
23
35
      PARAMETER (MAXNUP=500)
24
 
      INTEGER NUP,IDPRUP,IDUP(MAXNUP),ISTUP(MAXNUP),
25
 
     # MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP)
26
 
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
27
 
     # PUP(5,MAXNUP),VTIMUP(MAXNUP),SPINUP(MAXNUP)
28
 
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
29
 
      double precision scale1_lhe,scale2_lhe,percentage
30
 
      integer kwgtinfo,kexternal,jwgtnumpartn
31
 
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
32
 
      double precision xmuR_over_ref,xmuF1_over_ref,xmuF2_over_ref,
33
 
     # xQES_over_ref,pr_muR_over_ref,pr_muF1_over_ref,pr_muF2_over_ref,
34
 
     # tmp,yfactR(maxscales),yfactF(maxscales),xsecPDFr(0:maxPDFs)
35
 
      double precision xsecPDFr_acc(0:maxPDFs),xsecScale_acc(maxscales
36
 
     $     ,maxscales)
37
 
      double precision compute_rwgt_wgt_Sev,compute_rwgt_wgt_Sev_nbody
38
 
     &     ,compute_rwgt_wgt_Hev
39
 
      integer kr,kf,n,nng,nps,npairs,nsets,izero,itmp,idpdf(0:maxPDFs)
40
 
      parameter (izero=0)
41
 
      integer lef
42
 
      character*80 event_file,fname1
43
 
      character*140 buff
44
 
      character*10 MonteCarlo
45
 
      character*9 ch1
46
 
      character*20 parm(20)
47
 
      double precision value(20)
48
 
      logical AddInfoLHE
49
 
      external compute_rwgt_wgt_Sev,compute_rwgt_wgt_Sev_nbody
50
 
     &     ,compute_rwgt_wgt_Hev
51
 
      integer i_process
52
 
      common/c_i_process/i_process
 
36
      INTEGER NUP,IDPRUP,IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP)
 
37
     $     ,ICOLUP(2,MAXNUP)
 
38
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP(5,MAXNUP)
 
39
     $     ,VTIMUP(MAXNUP),SPINUP(MAXNUP)
53
40
c
54
41
      call setrun                !Sets up run parameters
55
42
 
56
 
 
57
43
      write(*,*) 'Enter event file name'
58
44
      read(*,*) event_file
59
45
 
61
47
      write(*,*)'      0 otherwise'
62
48
      read(*,*)isave
63
49
      if(isave.eq.1)then
64
 
        isave=9
 
50
        if (store_rwgt_info)then
 
51
           isave = -9
 
52
        else
 
53
           isave=9
 
54
        endif
65
55
      else
66
56
        isave=0
67
57
      endif
68
58
 
69
 
 
70
 
      xQES_over_ref=QES_over_ref
71
 
      xmuR_over_ref=muR_over_ref
72
 
      xmuF1_over_ref=muF1_over_ref
73
 
      xmuF2_over_ref=muF2_over_ref
74
 
      write(*,*) 'Using:'  
75
 
      write(*,*) 'QES_over_ref: ', xQES_over_ref
76
 
      write(*,*) 'muR_over_ref: ', xmuR_over_ref
77
 
      write(*,*) 'muF1_over_ref: ', xmuF1_over_ref
78
 
      write(*,*) 'muF2_over_ref: ', xmuF2_over_ref
79
 
      if (xmuF1_over_ref .ne. xmuF2_over_ref) then
80
 
          write(*,*) "The variables muF1_over_ref and muF2_over_ref" //
81
 
     1     " have to be set equal in the run_card.dat." //
82
 
     1     " Run cannot continue, quitting..."
83
 
          stop
84
 
      endif
85
 
 
86
59
      if(do_rwgt_scale)then
87
 
        yfactR(1)=1.d0
88
 
        yfactR(2)=rw_Rscale_up
89
 
        yfactR(3)=rw_Rscale_down
90
 
        yfactF(1)=1.d0
91
 
        yfactF(2)=rw_Fscale_up
92
 
        yfactF(3)=rw_Fscale_down
93
 
        write(*,*) 'Doing scale reweight:'
94
 
        write(*,*) rw_Fscale_down, ' < mu_F < ', rw_Fscale_up
95
 
        write(*,*) rw_Rscale_down, ' < mu_R < ', rw_Rscale_up
96
 
        numscales=3
97
 
      else
98
 
        numscales=0
 
60
         do nn=1,dyn_scale(0)
 
61
            if (lscalevar(nn)) then
 
62
               write (*,*) "Including scale uncertainties for"/
 
63
     $              /" dynamical_scale_choice",dyn_scale(nn)
 
64
               write (*,*) "  - renormalisation scale variation:"
 
65
     $              ,(scalevarR(i),i=1,nint(scalevarR(0)))
 
66
               write (*,*) "  - factorisation scale variation:"
 
67
     $              ,(scalevarF(i),i=1,nint(scalevarF(0)))
 
68
            else
 
69
               write (*,*) "Including central value for"/
 
70
     $              /" dynamical_scale_choice",dyn_scale(nn)
 
71
            endif
 
72
         enddo
99
73
      endif
100
74
 
101
 
c Note: when ipdf#0, the central PDF set will be used also as a reference
102
 
c for the scale uncertainty
103
75
      if(do_rwgt_pdf)then
104
 
 
105
 
        idpdf(0)=lhaid
106
 
        idpdf(1)=pdf_set_min
107
 
        itmp=pdf_set_max
108
 
        nsets=itmp-idpdf(1)+1
109
 
        write(*,*) 'Doing PDF reweight:'
110
 
        write(*,*) 'Central set id: ', idpdf(0)
111
 
        write(*,*) 'Min error set id: ', idpdf(1)
112
 
        write(*,*) 'Max error set id: ', itmp
113
 
        if(mod(nsets,2).ne.0)then
114
 
          write(*,*)'The number of error sets must be even',nsets
115
 
          stop
116
 
        else
117
 
          npairs=nsets/2
118
 
        endif
119
 
        do i=2,nsets
120
 
          idpdf(i)=idpdf(1)+i-1
121
 
        enddo
122
 
        if(nsets.gt.maxPDFs)then
123
 
          write(*,*)'Too many PDFs: increase maxPDFs in reweight0.inc'
124
 
          stop
125
 
        endif
126
 
c
127
 
        value(1)=idpdf(0)
128
 
        parm(1)='DEFAULT'
129
 
        call pdfset(parm,value)
130
 
c
131
 
        numPDFpairs=npairs
132
 
      else
133
 
        numPDFpairs=0
 
76
         do nn=1,lhaPDFid(0)
 
77
            if (lpdfvar(nn)) then
 
78
               write (*,*) "Including central PDF with "/
 
79
     $              /"uncertainties for "//lhaPDFsetname(nn)
 
80
            else
 
81
               write (*,*) "Including central PDF for "
 
82
     $              //lhaPDFsetname(nn)
 
83
            endif
 
84
c Load all the PDF sets (the 1st one has already by loaded by the call
 
85
c to "setrun")
 
86
            if (nn.gt.1) then
 
87
               call initpdfsetbynamem(nn,lhaPDFsetname(nn))
 
88
               if (lpdfvar(nn)) then
 
89
                  call numberPDFm(nn,nmemPDF(nn))
 
90
               else
 
91
                  nmemPDF(nn)=0
 
92
               endif
 
93
            endif
 
94
            if(nmemPDF(nn)+1.gt.maxPDFs)then
 
95
               write(*,*)'Too many PDFs: increase maxPDFs in '/
 
96
     $              /'reweight0.inc to ',numPDFs+1
 
97
               stop
 
98
            endif
 
99
         enddo
 
100
c start with central member of the first set
 
101
         call InitPDFm(1,0)
134
102
      endif
135
103
 
136
 
c$$$      call fk88strcat(event_file,'.rwgt',fname1)
137
104
      lef=index(event_file,' ')-1
138
105
      fname1=event_file(1:lef)//'.rwgt'
139
106
 
147
114
     &     XSECUP,XERRUP,XMAXUP,LPRUP)
148
115
 
149
116
      do i=1,min(10,maxevt)
150
 
        call read_lhef_event(ifile,
151
 
     &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
152
 
     &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
153
 
 
154
 
        if(buff(1:1).ne.'#')then
155
 
          write(*,*)'This event file cannot be reweighted [1]',i
156
 
          stop
157
 
        endif
158
 
        read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
159
 
     #                    fksfather_lhe,ipartner_lhe,
160
 
     #                    scale1_lhe,scale2_lhe,
161
 
     #                    kwgtinfo,kexternal,jwgtnumpartn,
162
 
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
163
 
        if((kwgtinfo.lt.1.or.kwgtinfo.gt.5) .and. kwgtinfo.ne.-5)then
164
 
          write(*,*)'This event file cannot be reweighted [2]',i
165
 
          write(*,*)kwgtinfo
166
 
          stop 1
 
117
         call read_lhef_event(ifile,
 
118
     &        NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
119
     &        IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
120
         if(buff(1:1).ne.'#')then
 
121
            write (*,*) 'This event file cannot be reweighted [1]',i
 
122
            stop
 
123
         endif
 
124
         read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,fksfather_lhe
 
125
     $        ,ipartner_lhe,scale1_lhe,scale2_lhe,kwgtinfo,kexternal
 
126
     $        ,jwgtnumpartn,wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
 
127
     $        ,wgtpdfmax
 
128
        if(kwgtinfo.ne.-5)then
 
129
           write (*,*) 'This event file cannot be reweighted [2]',i
 
130
           write (*,*) kwgtinfo
 
131
           stop 1
167
132
        endif
168
133
        if(i.eq.1)then
169
134
          saved_weight=abs(XWGTUP)
198
163
     &     XSECUP,XERRUP,XMAXUP,LPRUP)
199
164
 
200
165
c To keep track of the accumulated results:
201
 
      do ii=1,numscales
202
 
         do jj=1,numscales
203
 
            xsecScale_acc(jj,ii)=0d0
204
 
         enddo
205
 
      enddo
206
 
      do n=0,nsets
207
 
         xsecPDFr_acc(n)=0d0
208
 
      enddo
209
 
 
210
 
       nScontributions=1
 
166
      do kk=1,dyn_scale(0)
 
167
         if (lscalevar(kk)) then
 
168
            do ii=1,nint(scalevarF(0))
 
169
               do jj=1,nint(scalevarR(0))
 
170
                  xsecScale_acc(jj,ii,kk)=0d0
 
171
               enddo
 
172
            enddo
 
173
         else
 
174
            xsecScale_acc(1,1,kk)=0d0
 
175
         endif
 
176
      enddo
 
177
      do nn=1,lhaPDFid(0)
 
178
         if (lpdfvar(nn)) then
 
179
            do n=0,nmemPDF(nn)
 
180
               xsecPDFr_acc(n,nn)=0d0
 
181
            enddo
 
182
         else
 
183
            xsecPDFr_acc(0,nn)=0d0
 
184
         endif
 
185
      enddo
 
186
      nScontributions=1
211
187
 
212
188
c Determine the flavor map between the NLO and Born
213
189
      call find_iproc_map()
214
 
 
215
 
      
216
 
 
217
190
      do i=1,maxevt
218
 
 
219
 
         
220
191
         call read_lhef_event(ifile,
221
192
     &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
222
193
     &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
223
 
 
224
194
         if(buff(1:1).ne.'#')then
225
195
            write(*,*)'This event file cannot be reweighted [3]',i
226
196
            stop
230
200
     $        ,jwgtnumpartn,wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
231
201
     $        ,wgtpdfmax
232
202
 
233
 
         if (kwgtinfo.ne.-5) then
234
 
         call reweight_fill_extra_inverse()
235
 
 
236
 
        if(kwgtinfo.lt.1.or.kwgtinfo.gt.5)then
237
 
          write(*,*)'This event file cannot be reweighted [4]',i
238
 
          write(*,*)kwgtinfo
239
 
          stop
240
 
        endif
241
 
        if(wgtcentral.ne.0.d0.or.wgtmumin.ne.0.d0.or.
242
 
     #     wgtmumax.ne.0.d0.or.wgtpdfmin.ne.0.d0.or.
243
 
     #     wgtpdfmax.ne.0.d0)then
244
 
          write(*,*)'This event file was already reweighted',i
245
 
          write(*,*)wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
246
 
          stop
247
 
        endif
248
 
 
249
 
        if (kwgtinfo.eq.5) call reweight_settozero()
250
 
 
251
 
        if(do_rwgt_scale)then
252
 
 
253
 
          wgtmumin=1.d40
254
 
          wgtmumax=-1.d40
255
 
 
256
 
          do kr=1,3
257
 
            do kf=1,3
258
 
              wgtref=0d0
259
 
              pr_muR_over_ref=xmuR_over_ref*yfactR(kr)
260
 
              pr_muF1_over_ref=xmuF1_over_ref*yfactF(kf)
261
 
              pr_muF2_over_ref=pr_muF1_over_ref
262
 
              wgtxsecmu(kr,kf)=0d0
263
 
              if(iSorH_lhe.eq.1)then
264
 
c The nbody contributions
265
 
                 if (kwgtinfo.eq.5) then
266
 
                    call fill_reweight0inc_nbody(i_process)
267
 
                    wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)
268
 
     &                   +compute_rwgt_wgt_Sev_nbody(pr_muR_over_ref
269
 
     &                   ,pr_muF1_over_ref, pr_muF2_over_ref
270
 
     &                   ,xQES_over_ref, kwgtinfo)
271
 
                    call reweight_settozero()
272
 
                 endif
273
 
                 do ii=1,nScontributions
274
 
                    nFKSprocess_used=nFKSprocess_reweight(ii)
275
 
                    if (kwgtinfo.eq.5)
276
 
     &                   call fill_reweight0inc(nFKSprocess_used*2-1
277
 
     $                   ,i_process)
278
 
                    wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)+
279
 
     &                   compute_rwgt_wgt_Sev(pr_muR_over_ref
280
 
     &                   ,pr_muF1_over_ref, pr_muF2_over_ref
281
 
     &                   ,xQES_over_ref, kwgtinfo)
282
 
                    if (kwgtinfo.eq.5) call reweight_settozero()
283
 
                 enddo
284
 
              elseif(iSorH_lhe.eq.2)then
285
 
                 if (kwgtinfo.eq.5)
286
 
     &                call fill_reweight0inc(nFKSprocess_used*2,
287
 
     &                i_process)
288
 
                 wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)+
289
 
     &                compute_rwgt_wgt_Hev(pr_muR_over_ref
290
 
     &                ,pr_muF1_over_ref, pr_muF2_over_ref
291
 
     &                ,xQES_over_ref, kwgtinfo)
292
 
                 if (kwgtinfo.eq.5) call reweight_settozero()
293
 
              else
294
 
                 write(*,*)'Invalid value of iSorH_lhe',iSorH_lhe
295
 
                 stop
296
 
              endif
297
 
c
298
 
              tmp=wgtxsecmu(kr,kf)
299
 
              if(tmp.lt.wgtmumin)wgtmumin=tmp
300
 
              if(tmp.gt.wgtmumax)wgtmumax=tmp
301
 
            enddo
302
 
          enddo
303
 
 
304
 
          if (kwgtinfo.eq.5) then
305
 
             if (iSorH_lhe.eq.1) then
306
 
                wgtref=wgtref_nbody_all(i_process)
307
 
                do ii=1,nScontributions
308
 
                   wgtref=wgtref+ wgtref_all(nFKSprocess_reweight(ii)*2
309
 
     $                  -1,i_process)
310
 
                enddo
311
 
             else
312
 
                wgtref=wgtref_all(nFKSprocess_used*2,i_process)
313
 
             endif
314
 
          endif
315
 
 
316
 
          if(unweighted)then
317
 
            wgtcentral=wgtxsecmu(1,1)/wgtref
318
 
            wgtmumin=wgtmumin/wgtref
319
 
            wgtmumax=wgtmumax/wgtref
320
 
          else
321
 
            wgtcentral=wgtxsecmu(1,1)
322
 
          endif
323
 
 
324
 
        endif
325
 
 
326
 
        if(do_rwgt_pdf)then
327
 
 
328
 
          do n=0,nsets
329
 
             wgtref=0d0
330
 
             call InitPDF(n)
331
 
             wgtxsecPDF(n)=0d0
332
 
 
333
 
             if(iSorH_lhe.eq.1)then
334
 
c The nbody contributions
335
 
                if (kwgtinfo.eq.5) then
336
 
                   call fill_reweight0inc_nbody(i_process)
337
 
                   wgtxsecPDF(n)=wgtxsecPDF(n)
338
 
     $                  +compute_rwgt_wgt_Sev_nbody(xmuR_over_ref
339
 
     $                  ,xmuF1_over_ref, xmuF2_over_ref ,xQES_over_ref,
340
 
     $                  kwgtinfo)
341
 
                   call reweight_settozero()
342
 
                endif
343
 
                do ii=1,nScontributions
344
 
                   nFKSprocess_used=nFKSprocess_reweight(ii)
345
 
                   if (kwgtinfo.eq.5)
346
 
     &                  call fill_reweight0inc(nFKSprocess_used*2-1
347
 
     $                   ,i_process)
348
 
                   wgtxsecPDF(n)=wgtxsecPDF(n)+
349
 
     &                  compute_rwgt_wgt_Sev(xmuR_over_ref
350
 
     &                  ,xmuF1_over_ref, xmuF2_over_ref ,xQES_over_ref,
351
 
     &                  kwgtinfo)
352
 
                   if (kwgtinfo.eq.5) call reweight_settozero()
353
 
                enddo
354
 
             elseif(iSorH_lhe.eq.2)then
355
 
                if (kwgtinfo.eq.5)
356
 
     &               call fill_reweight0inc(nFKSprocess_used*2,
357
 
     &               i_process)
358
 
                wgtxsecPDF(n)=wgtxsecPDF(n)+
359
 
     &               compute_rwgt_wgt_Hev(xmuR_over_ref ,xmuF1_over_ref,
360
 
     &               xmuF2_over_ref ,xQES_over_ref, kwgtinfo)
361
 
                if (kwgtinfo.eq.5) call reweight_settozero()
362
 
             else
363
 
                write(*,*)'Invalid value of iSorH_lhe',iSorH_lhe
364
 
                stop
365
 
             endif
366
 
c
367
 
             if (kwgtinfo.eq.5) then
368
 
                if (iSorH_lhe.eq.1) then
369
 
                   wgtref=wgtref_nbody_all(i_process)
370
 
                   do ii=1,nScontributions
371
 
                      wgtref=wgtref+wgtref_all(nFKSprocess_reweight(ii)
372
 
     $                     *2-1,i_process)
373
 
                   enddo
374
 
                else
375
 
                   wgtref=wgtref_all(nFKSprocess_used*2,i_process)
376
 
                endif
377
 
             endif
378
 
             
379
 
             if(unweighted)then
380
 
                xsecPDFr(n)=wgtxsecPDF(n)/wgtref
381
 
             else
382
 
                xsecPDFr(n)=wgtxsecPDF(n)
383
 
             endif
384
 
          enddo
385
 
 
386
 
          if(do_rwgt_scale)then
387
 
            if(abs(xsecPDFr(0)/wgtcentral-1.d0).gt.1.d-6)then
388
 
              write(*,*)'Central valued computed with mu and PDF differ'
389
 
              write(*,*)xsecPDFr(0),wgtcentral
390
 
              stop
391
 
            endif
392
 
          else
393
 
            wgtcentral=xsecPDFr(0)
394
 
c The following serves to write on tape the reference cross section
395
 
c computed with the new parameters
396
 
            wgtxsecmu(1,1)=wgtxsecPDF(0)
397
 
          endif
398
 
 
399
 
          wgtpdfmin=0.d0
400
 
          wgtpdfmax=0.d0
401
 
 
402
 
          do n=1,npairs
403
 
            nps=2*n-1
404
 
            nng=2*n
405
 
 
406
 
            wgtpdfmin=wgtpdfmin+
407
 
     #                ( max(0.d0,
408
 
     #                      xsecPDFr(0)-xsecPDFr(nps),
409
 
     #                      xsecPDFr(0)-xsecPDFr(nng)) )**2
410
 
            wgtpdfmax=wgtpdfmax+
411
 
     #                ( max(0.d0,
412
 
     #                      xsecPDFr(nps)-xsecPDFr(0),
413
 
     #                      xsecPDFr(nng)-xsecPDFr(0)) )**2
414
 
          enddo
415
 
          wgtpdfmin=wgtcentral-sqrt(wgtpdfmin)
416
 
          wgtpdfmax=wgtcentral+sqrt(wgtpdfmax)
417
 
 
418
 
c Restore default PDFs
419
 
          call InitPDF(izero)
420
 
 
421
 
        endif
422
 
 
423
 
        else                      ! kwgtinfo.eq.-5
424
 
         
425
 
        call fill_wgt_info_from_rwgt_lines
426
 
        if (do_rwgt_scale)call reweight_scale_ext(yfactR,yfactF)
427
 
        if (do_rwgt_pdf)  call reweight_pdf_ext
428
 
        call fill_rwgt_arrays
429
 
 
430
 
        endif
431
 
 
432
 
        write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe,jfks_lhe,
433
 
     $       fksfather_lhe,ipartner_lhe, scale1_lhe,scale2_lhe, isave
434
 
     $       ,izero,izero, wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
435
 
     $       ,wgtpdfmax
436
 
 
 
203
c Do the actual reweighting.
 
204
         call fill_wgt_info_from_rwgt_lines
 
205
         if (do_rwgt_scale)call reweight_scale_ext
 
206
         if (do_rwgt_pdf)  call reweight_pdf_ext
 
207
         call fill_rwgt_arrays
 
208
 
 
209
         write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe,jfks_lhe,
 
210
     $        fksfather_lhe,ipartner_lhe, scale1_lhe,scale2_lhe, isave
 
211
     $        ,mexternal,izero, wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
 
212
     $        ,wgtpdfmax
437
213
 
438
214
c renormalize all the scale & PDF weights to have the same normalization
439
215
c as XWGTUP
440
 
        if(do_rwgt_scale)then
441
 
           do kr=1,3
442
 
              do kf=1,3
443
 
                 wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)/wgtref*XWGTUP
444
 
              enddo
445
 
           enddo
446
 
        endif
447
 
        if(do_rwgt_pdf)then
448
 
           do n=0,nsets
449
 
              wgtxsecPDF(n)=wgtxsecPDF(n)/wgtref*XWGTUP
450
 
           enddo
451
 
        endif
 
216
         if(do_rwgt_scale)then
 
217
            do kk=1,dyn_scale(0)
 
218
               if (lscalevar(kk)) then
 
219
                  do ii=1,nint(scalevarF(0))
 
220
                     do jj=1,nint(scalevarR(0))
 
221
                        wgtxsecmu(jj,ii,kk)=
 
222
     &                       wgtxsecmu(jj,ii,kk)/wgtref*XWGTUP
 
223
                     enddo
 
224
                  enddo
 
225
               else
 
226
                  wgtxsecmu(1,1,kk)=wgtxsecmu(1,1,kk)/wgtref*XWGTUP
 
227
               endif
 
228
            enddo
 
229
         endif
 
230
         if (do_rwgt_pdf) then
 
231
            do nn=1,lhaPDFid(0)
 
232
               if (lpdfvar(nn)) then
 
233
                  do n=0,nmemPDF(nn)
 
234
                     wgtxsecPDF(n,nn)=wgtxsecPDF(n,nn)/wgtref*XWGTUP
 
235
                  enddo
 
236
               else
 
237
                  wgtxsecPDF(0,nn)=wgtxsecPDF(0,nn)/wgtref*XWGTUP
 
238
               endif
 
239
            enddo
 
240
         endif
452
241
 
453
242
c Keep track of the accumulated results:
454
 
        if (numscales.gt.0) then
455
 
           do ii=1,numscales
456
 
              do jj=1,numscales
457
 
                 xsecScale_acc(ii,jj)=xsecScale_acc(ii,jj)+wgtxsecmu(ii
458
 
     $                ,jj)
459
 
              enddo
460
 
           enddo
461
 
        endif
462
 
        if (nsets.gt.0) then
463
 
           do n=0,nsets
464
 
              xsecPDFr_acc(n)=xsecPDFr_acc(n)+wgtxsecPDF(n)
465
 
           enddo
466
 
        endif
467
 
 
 
243
         if (do_rwgt_scale) then
 
244
            do kk=1,dyn_scale(0)
 
245
               if (lscalevar(kk)) then
 
246
                  do ii=1,nint(scalevarF(0))
 
247
                     do jj=1,nint(scalevarR(0))
 
248
                        xsecScale_acc(jj,ii,kk)=xsecScale_acc(jj,ii,kk)
 
249
     $                       +wgtxsecmu(jj,ii,kk)
 
250
                        
 
251
                     enddo
 
252
                  enddo
 
253
               else
 
254
                  xsecScale_acc(1,1,kk)=xsecScale_acc(1,1,kk)
 
255
     $                 +wgtxsecmu(1,1,kk)
 
256
               endif
 
257
            enddo
 
258
         endif
 
259
         if (do_rwgt_pdf) then
 
260
            do nn=1,lhaPDFid(0)
 
261
               if (lpdfvar(nn)) then
 
262
                  do n=0,nmemPDF(nn)
 
263
                     xsecPDFr_acc(n,nn)=xsecPDFr_acc(n,nn)+wgtxsecPDF(n
 
264
     $                    ,nn)
 
265
                  enddo
 
266
               else
 
267
                  xsecPDFr_acc(0,nn)=xsecPDFr_acc(0,nn)+wgtxsecPDF(0,nn)
 
268
               endif
 
269
            enddo
 
270
         endif
468
271
c Write event to disk:
469
 
        call write_lhef_event(ofile,
470
 
     &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
471
 
     &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
472
 
 
 
272
         call write_lhef_event(ofile,
 
273
     &        NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
274
     &        IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
275
         
473
276
      enddo
474
277
 
475
278
      write(ofile,'(a)')'</LesHouchesEvents>'
478
281
 
479
282
c Write the accumulated results to a file
480
283
      open (unit=34,file='scale_pdf_dependence.dat',status='unknown')
481
 
      write (34,*) numscales**2
482
 
      if (numscales.gt.0) then
483
 
         write (34,*) ((xsecScale_acc(ii,jj),ii=1,numscales),jj=1
484
 
     $        ,numscales)
485
 
      else
486
 
         write (34,*) ''
 
284
      if (do_rwgt_scale) then
 
285
         write (34,*) "scale variations:"
 
286
         do kk=1,dyn_scale(0)
 
287
            if (lscalevar(kk)) then
 
288
               write (34,*) dyn_scale(kk),nint(scalevarR(0))
 
289
     $              ,nint(scalevarF(0))
 
290
               write (34,*) ((xsecScale_acc(jj,ii,kk),jj=1
 
291
     $              ,nint(scalevarR(0))),ii=1,nint(scalevarF(0)))
 
292
            else
 
293
               write (34,*) dyn_scale(kk),1,1
 
294
               write (34,*) xsecScale_acc(1,1,kk)
 
295
            endif
 
296
         enddo
487
297
      endif
488
 
      if (nsets.gt.0) then
489
 
         write (34,*) nsets + 1
490
 
         write (34,*) (xsecPDFr_acc(n),n=0,nsets)
491
 
      else
492
 
         write(34,*) nsets
493
 
         write (34,*) ''
 
298
      if (do_rwgt_pdf) then
 
299
         write (34,*) "pdf variations:"
 
300
         do nn=1,lhaPDFid(0)
 
301
            if (lpdfvar(nn)) then
 
302
               write (34,*) trim(adjustl(lhaPDFsetname(nn))),
 
303
     $              nmemPDF(nn)+1
 
304
               write (34,*) (xsecPDFr_acc(n,nn),n=0,nmemPDF(nn))
 
305
            else
 
306
               write(34,*) lhaPDFsetname(nn),nmemPDF(nn) + 1
 
307
               write (34,*) xsecPDFr_acc(0,nn)
 
308
            endif
 
309
         enddo
494
310
      endif
495
311
      close(34)
496
312
 
497
 
 201  format(a9,1x,i1,4(1x,i2),2(1x,e14.8),1x,i1,2(1x,i2),5(1x,e14.8))
 
313
 201  format(a9,1x,i1,4(1x,i2),2(1x,e14.8),1x,i2,2(1x,i2),5(1x,e14.8))
498
314
 
499
315
      end
500
316
 
501
 
 
502
317
c Dummy subroutine (normally used with vegas/mint when resuming plots)
503
318
      subroutine resume()
504
319
      end
561
376
      end
562
377
 
563
378
 
564
 
 
565
 
 
566
 
 
567
 
 
568
 
 
569
 
 
570
 
 
571
 
      
572
379
      subroutine fill_wgt_info_from_rwgt_lines
573
380
      implicit none
574
381
      include 'nexternal.inc'
575
382
      include 'c_weight.inc'
576
383
      include 'reweight0.inc'
577
 
      integer i,idum,j,k,momenta_conf
 
384
      integer i,idum,j,k,momenta_conf(2),ii
578
385
      icontr=n_ctr_found
579
386
      iwgt=1
580
387
      do i=1,icontr
581
 
         read(n_ctr_str(i),*)(wgt(j,i),j=1,3),idum,(pdg(j,i),j=1
582
 
     &        ,nexternal),QCDpower(i),(bjx(j,i),j=1,2),(scales2(j,i),j=1
583
 
     &        ,3),momenta_conf,itype(i),nFKS(i),wgts(1,i)
584
 
         do j=1,nexternal
585
 
            do k=0,3
586
 
               momenta(k,j,i)=momenta_str(k,j,momenta_conf)
 
388
         read(n_ctr_str(i),*)(wgt(j,i),j=1,3),(wgt_ME_tree(j,i),j=1,2)
 
389
     &        ,idum,(pdg(j,i),j=1,nexternal),QCDpower(i),(bjx(j,i),j=1
 
390
     &        ,2),(scales2(j,i),j=1,3),g_strong(i),(momenta_conf(j),j=1
 
391
     &        ,2),itype(i),nFKS(i),idum,idum,idum,wgts(1,i)
 
392
         do ii=1,2
 
393
            do j=1,nexternal
 
394
               do k=0,3
 
395
                  if (momenta_conf(ii).gt.0) then
 
396
                     momenta_m(k,j,ii,i)=momenta_str(k,j
 
397
     $                                               ,momenta_conf(ii))
 
398
                  else
 
399
                     momenta_m(k,j,ii,i)=-99d0
 
400
                     exit
 
401
                  endif
 
402
               enddo
587
403
            enddo
588
404
         enddo
589
405
      enddo
590
406
      end
591
407
      
592
 
      subroutine reweight_scale_ext(yfactR,yfactF)
 
408
      subroutine reweight_scale_ext
593
409
      implicit none
594
410
      include 'nexternal.inc'
595
411
      include 'c_weight.inc'
596
412
      include 'run.inc'
597
 
      integer i,pd,lp,iwgt_save,kr,kf
598
 
      double precision yfactR(3),yfactF(3),mu2_f(3),mu2_r(3),xlum(3)
599
 
     &     ,pdg2pdf,mu2_q,rwgt_muR_dep_fac,g(3),alphas,pi
 
413
      include 'reweight0.inc'
 
414
      integer i,pd,lp,iwgt_save,kr,kf,dd
 
415
      double precision mu2_f(maxscales),mu2_r(maxscales),xlum(maxscales)
 
416
     $     ,pdg2pdf,mu2_q,rwgt_muR_dep_fac,g(maxscales),alphas,pi
 
417
     $     ,c_mu2_r,c_mu2_f
600
418
      parameter (pi=3.14159265358979323846d0)
601
419
      external pdg2pdf,rwgt_muR_dep_fac,alphas
602
420
      iwgt_save=iwgt
603
421
      do i=1,icontr
604
422
         iwgt=iwgt_save
605
423
         mu2_q=scales2(1,i)
606
 
         do kr=1,3
607
 
            mu2_r(kr)=scales2(2,i)*yfactR(kr)**2
 
424
         do dd=1,dyn_scale(0)
 
425
            call set_mu_central(i,dd,c_mu2_r,c_mu2_f)
 
426
            do kr=1,nint(scalevarR(0))
 
427
               if ((.not. lscalevar(dd)) .and. kr.ne.1) exit
 
428
               mu2_r(kr)=c_mu2_r*scalevarR(kr)**2
608
429
c Update the strong coupling
609
 
            g(kr)=sqrt(4d0*pi*alphas(sqrt(mu2_r(kr))))
610
 
         enddo
611
 
         do kf=1,3
612
 
            mu2_f(kf)=scales2(3,i)*yfactF(kf)**2
 
430
               g(kr)=sqrt(4d0*pi*alphas(sqrt(mu2_r(kr))))
 
431
            enddo
 
432
            do kf=1,nint(scalevarF(0))
 
433
               if ((.not. lscalevar(dd)) .and. kf.ne.1) exit
 
434
               mu2_f(kf)=c_mu2_f*scalevarF(kf)**2
613
435
c call the PDFs
614
 
            xlum(kf)=1d0
615
 
            LP=SIGN(1,LPP(1))
616
 
            pd=pdg(1,i)
617
 
            if (pd.eq.21) pd=0
618
 
            xlum(kf)=xlum(kf)*PDG2PDF(ABS(LPP(1)),pd*LP,bjx(1,i)
619
 
     &           ,DSQRT(mu2_f(kf)))
620
 
            LP=SIGN(1,LPP(2))
621
 
            pd=pdg(2,i)
622
 
            if (pd.eq.21) pd=0
623
 
            xlum(kf)=xlum(kf)*PDG2PDF(ABS(LPP(2)),pd*LP,bjx(2,i)
624
 
     &           ,DSQRT(mu2_f(kf)))
625
 
         enddo
626
 
         do kr=1,3
627
 
            do kf=1,3
628
 
               iwgt=iwgt+1 ! increment the iwgt for the wgts() array
629
 
               if (iwgt.gt.max_wgt) then
630
 
                  write (*,*) 'ERROR too many weights in reweight_scale'
631
 
     &                 ,iwgt,max_wgt
632
 
                  stop 1
633
 
               endif
 
436
               xlum(kf)=1d0
 
437
               LP=SIGN(1,LPP(1))
 
438
               pd=pdg(1,i)
 
439
               if (pd.eq.21) pd=0
 
440
               xlum(kf)=xlum(kf)*PDG2PDF(ABS(LPP(1)),pd*LP,bjx(1,i)
 
441
     &              ,DSQRT(mu2_f(kf)))
 
442
               LP=SIGN(1,LPP(2))
 
443
               pd=pdg(2,i)
 
444
               if (pd.eq.21) pd=0
 
445
               xlum(kf)=xlum(kf)*PDG2PDF(ABS(LPP(2)),pd*LP,bjx(2,i)
 
446
     &              ,DSQRT(mu2_f(kf)))
 
447
            enddo
 
448
            do kf=1,nint(scalevarF(0))
 
449
               if ((.not. lscalevar(dd)) .and. kf.ne.1) exit
 
450
               do kr=1,nint(scalevarR(0))
 
451
                  if ((.not. lscalevar(dd)) .and. kr.ne.1) exit
 
452
                  iwgt=iwgt+1   ! increment the iwgt for the wgts() array
 
453
                  if (iwgt.gt.max_wgt) then
 
454
                     write (*,*) 'ERROR too many weights in '/
 
455
     $                    /'reweight_scale',iwgt,max_wgt
 
456
                     stop 1
 
457
                  endif
634
458
c add the weights to the array
635
 
               wgts(iwgt,i)=xlum(kf) * (wgt(1,i)+wgt(2,i)*log(mu2_r(kr)
636
 
     &              /mu2_q)+wgt(3,i)*log(mu2_f(kf)/mu2_q))*g(kr)
637
 
     &              **QCDpower(i)
638
 
               wgts(iwgt,i)=wgts(iwgt,i)
639
 
     &              *rwgt_muR_dep_fac(sqrt(mu2_r(kr)),sqrt(mu2_r(1)))
 
459
                  wgts(iwgt,i)=xlum(kf) * (wgt(1,i)+wgt(2,i)
 
460
     $                 *log(mu2_r(kr)/mu2_q)+wgt(3,i)*log(mu2_f(kf)
 
461
     $                 /mu2_q))*g(kr)**QCDpower(i)
 
462
                  wgts(iwgt,i)=wgts(iwgt,i)*rwgt_muR_dep_fac(
 
463
     &                          sqrt(mu2_r(kr)),sqrt(scales2(2,i)))
 
464
               enddo
640
465
            enddo
641
466
         enddo
642
467
      enddo
650
475
      include 'c_weight.inc'
651
476
      include 'run.inc'
652
477
      include 'reweight0.inc'
653
 
      integer i,pd,lp,iwgt_save,izero,n
 
478
      integer i,pd,lp,iwgt_save,izero,n,nn,iset,imem
654
479
      parameter (izero=0)
655
480
      double precision mu2_f,mu2_r,pdg2pdf,mu2_q,rwgt_muR_dep_fac
656
481
     &     ,xlum,alphas,g,pi
657
482
      parameter (pi=3.14159265358979323846d0)
658
483
      external pdg2pdf,rwgt_muR_dep_fac,alphas
659
 
      do n=0,numPDFpairs*2
660
 
         iwgt=iwgt+1
661
 
         if (iwgt.gt.max_wgt) then
662
 
            write (*,*) 'ERROR too many weights in reweight_pdf',iwgt
663
 
     &           ,max_wgt
664
 
            stop 1
665
 
         endif
666
 
         call InitPDF(n)
667
 
         do i=1,icontr
668
 
            mu2_q=scales2(1,i)
669
 
            mu2_r=scales2(2,i)
670
 
            mu2_f=scales2(3,i)
 
484
      do nn=1,lhaPDFid(0)
 
485
         do n=0,nmemPDF(nn)
 
486
            if ((.not. lpdfvar(nn)) .and. n.ne.0) exit
 
487
            iwgt=iwgt+1
 
488
            if (iwgt.gt.max_wgt) then
 
489
               write (*,*) 'ERROR too many weights in reweight_pdf',iwgt
 
490
     &              ,max_wgt
 
491
               stop 1
 
492
            endif
 
493
            call InitPDFm(nn,n)
 
494
            do i=1,icontr
 
495
               mu2_q=scales2(1,i)
 
496
               mu2_r=scales2(2,i)
 
497
               mu2_f=scales2(3,i)
671
498
c alpha_s
672
 
            g=sqrt(4d0*pi*alphas(sqrt(mu2_r)))
 
499
               g=sqrt(4d0*pi*alphas(sqrt(mu2_r)))
673
500
c call the PDFs
674
 
            xlum=1d0
675
 
            LP=SIGN(1,LPP(1))
676
 
            pd=pdg(1,i)
677
 
            if (pd.eq.21) pd=0
678
 
            xlum=xlum*PDG2PDF(ABS(LPP(1)),pd*LP,bjx(1,i),DSQRT(mu2_f))
679
 
            LP=SIGN(1,LPP(2))
680
 
            pd=pdg(2,i)
681
 
            if (pd.eq.21) pd=0
682
 
            xlum=xlum*PDG2PDF(ABS(LPP(2)),pd*LP,bjx(2,i),DSQRT(mu2_f))
 
501
               xlum=1d0
 
502
               LP=SIGN(1,LPP(1))
 
503
               pd=pdg(1,i)
 
504
               if (pd.eq.21) pd=0
 
505
               xlum=xlum*
 
506
     &            PDG2PDF(ABS(LPP(1)),pd*LP,bjx(1,i),DSQRT(mu2_f))
 
507
               LP=SIGN(1,LPP(2))
 
508
               pd=pdg(2,i)
 
509
               if (pd.eq.21) pd=0
 
510
               xlum=xlum*
 
511
     &             PDG2PDF(ABS(LPP(2)),pd*LP,bjx(2,i),DSQRT(mu2_f))
683
512
c add the weights to the array
684
 
            wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q) +
685
 
     &           wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
686
 
            wgts(iwgt,i)=wgts(iwgt,i)
687
 
     &           *rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r))
 
513
               wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q)
 
514
     &              +wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
 
515
               wgts(iwgt,i)=wgts(iwgt,i)
 
516
     &              *rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r))
 
517
            enddo
688
518
         enddo
689
519
      enddo
690
 
      call InitPDF(izero)
 
520
c reset to the 0th member of the 1st set
 
521
      call InitPDFm(1,0)
691
522
      return
692
523
      end
693
524
      
697
528
      include 'nexternal.inc'
698
529
      include 'c_weight.inc'
699
530
      include 'reweight0.inc'
700
 
      integer kr,kf,n,iw,i
701
 
      do kr=1,numscales
702
 
         do kf=1,numscales
703
 
            wgtxsecmu(kr,kf)=0d0
704
 
         enddo
 
531
      include 'run.inc'
 
532
      integer ii,jj,kk,nn,n,iw,i
 
533
      do kk=1,dyn_scale(0)
 
534
         if (lscalevar(kk)) then
 
535
            do ii=1,nint(scalevarF(0))
 
536
               do jj=1,nint(scalevarR(0))
 
537
                  wgtxsecmu(jj,ii,kk)=0d0
 
538
               enddo
 
539
            enddo
 
540
         else
 
541
            wgtxsecmu(1,1,kk)=0d0
 
542
         endif
705
543
      enddo
706
 
      do n=0,numPDFpairs*2
707
 
         wgtxsecPDF(n)=0d0
 
544
      do nn=1,lhaPDFid(0)
 
545
         if (lpdfvar(nn)) then
 
546
            do n=0,nmemPDF(nn)
 
547
               wgtxsecPDF(n,nn)=0d0
 
548
            enddo
 
549
         else
 
550
            wgtxsecPDF(0,nn)=0d0
 
551
         endif
708
552
      enddo
709
553
      do i=1,icontr
710
554
         iw=2
711
 
         do kr=1,numscales
712
 
            do kf=1,numscales
713
 
               wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)+wgts(iw,i)
714
 
               iw=iw+1
715
 
            enddo
716
 
         enddo
717
 
         do n=0,numPDFpairs*2
718
 
            wgtxsecPDF(n)=wgtxsecPDF(n)+wgts(iw,i)
719
 
            iw=iw+1
720
 
         enddo
 
555
         if (do_rwgt_scale) then
 
556
            do kk=1,dyn_scale(0)
 
557
               if (lscalevar(kk)) then
 
558
                  do ii=1,nint(scalevarF(0))
 
559
                     do jj=1,nint(scalevarR(0))
 
560
                        wgtxsecmu(jj,ii,kk)=wgtxsecmu(jj,ii,kk)+wgts(iw,i)
 
561
                        iw=iw+1
 
562
                     enddo
 
563
                  enddo
 
564
               else
 
565
                  wgtxsecmu(1,1,kk)=wgtxsecmu(1,1,kk)+wgts(iw,i)
 
566
                  iw=iw+1
 
567
               endif
 
568
            enddo
 
569
         endif
 
570
         if (do_rwgt_pdf) then
 
571
            do nn=1,lhaPDFid(0)
 
572
               if (lpdfvar(nn)) then
 
573
                  do n=0,nmemPDF(nn)
 
574
                     wgtxsecPDF(n,nn)=wgtxsecPDF(n,nn)+wgts(iw,i)
 
575
                     iw=iw+1
 
576
                  enddo
 
577
               else
 
578
                  wgtxsecPDF(0,nn)=wgtxsecPDF(0,nn)+wgts(iw,i)
 
579
                  iw=iw+1
 
580
               endif
 
581
            enddo
 
582
         endif
721
583
      enddo
722
 
      if (numscales.eq.0) then
723
 
         wgtxsecmu(1,1)=wgtxsecPDF(0)
724
 
      endif
725
584
      return
726
585
      end
727
586
 
728
587
      
 
588
      subroutine set_mu_central(ic,dd,c_mu2_r,c_mu2_f)
 
589
      implicit none
 
590
      include 'nexternal.inc'
 
591
      include 'c_weight.inc'
 
592
      include 'reweight0.inc'
 
593
      include 'run.inc'
 
594
      integer ic,dd,i,j
 
595
      double precision c_mu2_r,c_mu2_f,muR,muF,pp(0:3,nexternal)
 
596
      if (dd.eq.1) then
 
597
         c_mu2_r=scales2(2,ic)
 
598
         c_mu2_f=scales2(3,ic)
 
599
      else
 
600
c need to recompute the scales using the momenta
 
601
         dynamical_scale_choice=dyn_scale(dd)
 
602
         do i=1,nexternal
 
603
            do j=0,3
 
604
               pp(j,i)=momenta(j,i,ic)
 
605
            enddo
 
606
         enddo
 
607
         call set_ren_scale(pp,muR)
 
608
         c_mu2_r=muR**2
 
609
         call set_fac_scale(pp,muF)
 
610
         c_mu2_f=muF**2
 
611
c     reset the default dynamical_scale_choice
 
612
         dynamical_scale_choice=dyn_scale(1)
 
613
      endif
 
614
      return
 
615
      end