← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/madfks_plot.uncer.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/madfks_plot.uncer.f
 
1
c
 
2
c
 
3
c Plotting routines
 
4
c
 
5
c
 
6
      subroutine initplot
 
7
      implicit none
 
8
c Book histograms in this routine. Use mbook or bookup. The entries
 
9
c of these routines are real*8
 
10
c
 
11
c begin reweight
 
12
      include 'reweight0.inc'
 
13
      include 'reweightNLO.inc'
 
14
      character*4 eesc(maxscales,maxscales)
 
15
      character*4 eepdf(0:maxPDFs)
 
16
      common/ceestrings/eesc,eepdf
 
17
c The following are useful, but not crucial. One has maxset replicas
 
18
c (corresponding eg to no-cuts, cuts1, cuts2,...) of up to maxplot plots
 
19
c (corresponding eg to pt1, eta1, pt2, ...), for each PDF and scale choices.
 
20
c The integers istartscales and istartPDFs are defined below
 
21
      integer maxset,maxplot,istartscales,istartPDFs
 
22
      common/cnumplots/maxset,maxplot,istartscales,istartPDFs
 
23
      integer isc,jsc,ipdf
 
24
c end reweight
 
25
      include 'run.inc'
 
26
      integer i,kk,kk1
 
27
      character*4 ddstr
 
28
      character*6 cc(2)
 
29
      data cc/' NLO  ',' Born '/
 
30
c
 
31
c resets histograms
 
32
      call inihist
 
33
c begin reweight
 
34
      call plotunc_setup()
 
35
      maxset=2
 
36
      maxplot=50
 
37
      istartscales=maxset*maxplot
 
38
      istartPDFs=(1+numscales*numscales)*maxset*maxplot
 
39
c end reweight
 
40
c
 
41
      do i=1,maxset
 
42
c default
 
43
        kk=(i-1)*maxplot
 
44
        ddstr='def'
 
45
        call bookup(kk+1,'total rate'//cc(i)//ddstr,
 
46
     #              1.0d0,0.5d0,5.5d0)
 
47
c begin reweight
 
48
        if(.not.doNLOscaleunc)goto 111
 
49
        do isc=1,numscales
 
50
          do jsc=1,numscales
 
51
            ddstr=eesc(isc,jsc)
 
52
            kk1=kk+istartscales+
 
53
     #          maxset*maxplot*((isc-1)*numscales+(jsc-1))
 
54
            call bookup(kk1+1,'total rate'//cc(i)//ddstr,
 
55
     #                  1.0d0,0.5d0,5.5d0)
 
56
          enddo
 
57
        enddo
 
58
 111    continue
 
59
        if(.not.doNLOPDFunc)goto 222
 
60
        do ipdf=0,numPDFs-1
 
61
          ddstr=eepdf(ipdf)
 
62
          kk1=kk+istartPDFs+maxset*maxplot*ipdf
 
63
          call bookup(kk1+1,'total rate'//cc(i)//ddstr,
 
64
     #                1.0d0,0.5d0,5.5d0)
 
65
        enddo
 
66
 222    continue
 
67
c end reweight
 
68
      enddo
 
69
      return
 
70
      end
 
71
 
 
72
 
 
73
      subroutine topout
 
74
      implicit none
 
75
      character*14 ytit
 
76
      logical usexinteg,mint
 
77
      common/cusexinteg/usexinteg,mint
 
78
      integer itmax,ncall
 
79
      common/citmax/itmax,ncall
 
80
      real*8 xnorm1,xnorm2
 
81
      logical unwgt
 
82
      double precision evtsgn
 
83
      common /c_unwgt/evtsgn,unwgt
 
84
      integer i,kk,kk1
 
85
      include 'dbook.inc'
 
86
c begin reweight
 
87
      include 'reweight0.inc'
 
88
      include 'reweightNLO.inc'
 
89
      character*4 eesc(maxscales,maxscales)
 
90
      character*4 eepdf(0:maxPDFs)
 
91
      common/ceestrings/eesc,eepdf
 
92
      integer isc,jsc,ipdf
 
93
c The following are defined in initplot
 
94
      integer maxset,maxplot,istartscales,istartPDFs
 
95
      common/cnumplots/maxset,maxplot,istartscales,istartPDFs
 
96
c end reweight
 
97
c
 
98
      if (unwgt) then
 
99
         ytit='events per bin'
 
100
      else
 
101
         ytit='sigma per bin '
 
102
      endif
 
103
      xnorm1=1.d0/float(itmax)
 
104
      xnorm2=1.d0/float(ncall*itmax)
 
105
      do i=1,NPLOTS
 
106
        if(usexinteg.and..not.mint) then
 
107
           call mopera(i,'+',i,i,xnorm1,0.d0)
 
108
        elseif(mint) then
 
109
           call mopera(i,'+',i,i,xnorm2,0.d0)
 
110
        endif
 
111
        call mfinal(i)
 
112
      enddo
 
113
c default
 
114
      do i=1,maxset
 
115
        kk=(i-1)*maxplot
 
116
        call multitop(kk+1,3,2,'total rate',ytit,'LIN')
 
117
      enddo
 
118
c begin reweight
 
119
      if(.not.doNLOscaleunc)goto 111
 
120
      do isc=1,numscales
 
121
        do jsc=1,numscales
 
122
c within these do loops, repeat what done for default
 
123
          do i=1,maxset
 
124
            kk=(i-1)*maxplot
 
125
            kk1=kk+istartscales+
 
126
     #          maxset*maxplot*((isc-1)*numscales+(jsc-1))
 
127
            call multitop(kk1+1,3,2,'total rate',ytit,'LIN')
 
128
          enddo
 
129
c end repeat
 
130
        enddo
 
131
      enddo
 
132
 111  continue
 
133
      if(.not.doNLOPDFunc)goto 222
 
134
      do ipdf=0,numPDFs-1
 
135
c within this do loop, repeat what done for default
 
136
        do i=1,maxset
 
137
          kk=(i-1)*maxplot
 
138
          kk1=kk+istartPDFs+maxset*maxplot*ipdf
 
139
          call multitop(kk1+1,3,2,'total rate',ytit,'LIN')
 
140
        enddo
 
141
c end repeat
 
142
      enddo
 
143
 222  continue
 
144
c end reweight
 
145
      return                
 
146
      end
 
147
 
 
148
 
 
149
      subroutine outfun(pp,ybst_til_tolab,www,itype)
 
150
C
 
151
C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
 
152
C
 
153
C In MadFKS, the momenta PP given in input to this function are in the
 
154
C reduced parton c.m. frame. If need be, boost them to the lab frame.
 
155
C The rapidity of this boost is
 
156
C
 
157
C       YBST_TIL_TOLAB
 
158
C
 
159
C also given in input
 
160
C
 
161
C This is the rapidity that enters in the arguments of the sinh() and
 
162
C cosh() of the boost, in such a way that
 
163
C       ylab = ycm - ybst_til_tolab
 
164
C where ylab is the rapidity in the lab frame and ycm the rapidity
 
165
C in the center-of-momentum frame.
 
166
C
 
167
C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
 
168
      implicit none
 
169
      include 'nexternal.inc'
 
170
      real*8 pp(0:3,nexternal),ybst_til_tolab,www
 
171
      integer itype
 
172
      real*8 var
 
173
      integer kk,kk1
 
174
c begin reweight
 
175
      include 'reweight.inc'
 
176
      include 'reweightNLO.inc'
 
177
      integer irwgt,isc,jsc,ipdf
 
178
      double precision wgtden,wwwrsc
 
179
c end reweight
 
180
c The following are defined in initplot
 
181
      integer maxset,maxplot,istartscales,istartPDFs
 
182
      common/cnumplots/maxset,maxplot,istartscales,istartPDFs
 
183
c masses
 
184
      double precision pmass(nexternal)
 
185
      common/to_mass/pmass
 
186
c
 
187
      kk=0
 
188
      if(itype.eq.20)kk=50
 
189
c begin reweight
 
190
      if(itype.eq.11)then
 
191
        irwgt=1
 
192
        wgtden=wgtrefNLO11
 
193
      elseif(itype.eq.12)then
 
194
        irwgt=2
 
195
        wgtden=wgtrefNLO12
 
196
      elseif(itype.eq.20)then
 
197
        irwgt=3
 
198
        wgtden=wgtrefNLO20
 
199
      else
 
200
        write(*,*)'Error in outfun: unknown itype',itype
 
201
        stop
 
202
      endif
 
203
c end reweight
 
204
c
 
205
      var=1.d0
 
206
      call mfill(kk+1,var,www)
 
207
c begin reweight
 
208
      if(.not.doNLOscaleunc.or.wgtden.eq.0.d0)goto 111
 
209
      do isc=1,numscales
 
210
        do jsc=1,numscales
 
211
          kk1=kk+istartscales+
 
212
     #        maxset*maxplot*((isc-1)*numscales+(jsc-1))
 
213
          wwwrsc=www * wgtNLOxsecmu(irwgt,isc,jsc)/wgtden
 
214
          call mfill(kk1+1,var,wwwrsc)
 
215
        enddo
 
216
      enddo
 
217
 111  continue
 
218
      if(.not.doNLOPDFunc.or.wgtden.eq.0.d0)goto 222
 
219
      do ipdf=0,numPDFs-1
 
220
        kk1=kk+istartPDFs+maxset*maxplot*ipdf
 
221
        wwwrsc=www * wgtNLOxsecPDF(irwgt,ipdf)/wgtden
 
222
        call mfill(kk1+1,var,wwwrsc)
 
223
      enddo
 
224
 222  continue
 
225
c end reweight
 
226
 999  return      
 
227
      end
 
228
 
 
229
 
 
230
 
 
231
c
 
232
c
 
233
c Utilities
 
234
c
 
235
c
 
236
      subroutine plotunc_setup()
 
237
      implicit none
 
238
      include 'reweight0.inc'
 
239
      include 'reweightNLO.inc'
 
240
      integer i,j
 
241
      character*4 eesc(maxscales,maxscales)
 
242
      character*4 eepdf(0:maxPDFs)
 
243
      common/ceestrings/eesc,eepdf
 
244
c
 
245
      if(.not.doNLOscaleunc)goto 111
 
246
c numscales*numscales tags sij
 
247
      if(numscales.gt.9)then
 
248
        write(*,*)'Number of scales too large: wrong sij tags',
 
249
     #            numscales
 
250
        stop
 
251
      endif
 
252
      do i=1,numscales
 
253
        do j=1,numscales
 
254
          eesc(i,j)='s'
 
255
          call fk88strnum(eesc(i,j),i)
 
256
          call fk88strnum(eesc(i,j),j)
 
257
        enddo
 
258
      enddo
 
259
 111  continue
 
260
      if(.not.doNLOPDFunc)goto 222
 
261
c numPDFs tags pxyz
 
262
      if(numPDFs.gt.999)then
 
263
        write(*,*)'Number of PDFs too large: wrong pxyz tags',
 
264
     #            numPDFs
 
265
        stop
 
266
      endif
 
267
      do i=0,numPDFs-1
 
268
        if(i.lt.10)then
 
269
          eepdf(i)='p00'
 
270
        elseif(i.lt.100)then
 
271
          eepdf(i)='p0'
 
272
        else
 
273
          eepdf(i)='p'
 
274
        endif
 
275
        call fk88strnum(eepdf(i),i)
 
276
      enddo
 
277
 222  continue
 
278
      return
 
279
      end
 
280
 
 
281
 
 
282
      subroutine fk88strnum(string,num)
 
283
c- writes the number num on the string string starting at the blank
 
284
c- following the last non-blank character
 
285
      character * (*) string
 
286
      character * 20 tmp
 
287
      l = len(string)
 
288
      write(tmp,'(i15)')num
 
289
      j=1
 
290
      dowhile(tmp(j:j).eq.' ')
 
291
        j=j+1
 
292
      enddo
 
293
      ipos = ifk88istrl(string)
 
294
      ito = ipos+1+(15-j)
 
295
      if(ito.gt.l) then
 
296
         write(*,*)'error, string too short'
 
297
         write(*,*) string
 
298
         stop
 
299
      endif
 
300
 
 
301
 
 
302
      string(ipos+1:ito)=tmp(j:)
 
303
      end
 
304
      function ifk88istrl(string)
 
305
c returns the position of the last non-blank character in string
 
306
      character * (*) string
 
307
      i = len(string)
 
308
      dowhile(i.gt.0.and.string(i:i).eq.' ')
 
309
         i=i-1
 
310
      enddo
 
311
      ifk88istrl = i
 
312
      end