← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

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

  • Committer: Marco Zaro
  • Date: 2012-08-28 21:06:34 UTC
  • mto: (78.35.14 AutoMint)
  • mto: This revision was merged to the branch mainline in revision 249.
  • Revision ID: marco.zaro@gmail.com-20120828210634-5a06yvda3hplw8ur
doing some renaming:
 Template/FKS-born -> Template/NLO
 fks_born.py -> fks_base.py
 fks_born_helas_objects.py -> fks_helas_objects.py
 export_fks_born.py -> export_fks.py

also functions/classes and tests renamed
all unit tests ok, exporting ok

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/handling_lhe_events.f
 
1
c Utility routines for LHEF. Originally taken from collect_events.f
 
2
c
 
3
c Note: the routines read_lhef_event and write_lhef_event use the common
 
4
c blocks in reweight0.inc, relevant to reweight information. This is
 
5
c independent of the process, and in particular of process-related
 
6
c parameters such as nexternal, which is replaced here by (its supposed)
 
7
c upper bound maxparticles. The arrays which have one dimension defined
 
8
c by maxparticles may have a correspondence with process-specific ones,
 
9
c and the dimensions of the latter are typically defined by nexternal.
 
10
c Hence, one may need an explicit copy of one onto the other
 
11
c
 
12
 
 
13
      subroutine write_lhef_header(ifile,nevents,MonteCarlo)
 
14
      implicit none 
 
15
      integer ifile,nevents
 
16
      character*10 MonteCarlo
 
17
c
 
18
      write(ifile,'(a)')
 
19
     #     '<LesHouchesEvents version="1.0">'
 
20
      write(ifile,'(a)')
 
21
     #     '  <!--'
 
22
      write(ifile,'(a)')
 
23
     #     MonteCarlo
 
24
      write(ifile,'(a)')
 
25
     #     '  -->'
 
26
      write(ifile,'(a)')
 
27
     #     '  <header>'
 
28
      write(ifile,250)nevents
 
29
      write(ifile,'(a)')
 
30
     #     '  </header>'
 
31
 250  format(1x,i8)
 
32
      return
 
33
      end
 
34
 
 
35
 
 
36
 
 
37
      subroutine write_lhef_header_string(ifile,string,MonteCarlo)
 
38
      implicit none 
 
39
      integer ifile
 
40
      character*10 MonteCarlo,string
 
41
c
 
42
      write(ifile,'(a)')
 
43
     #     '<LesHouchesEvents version="1.0">'
 
44
      write(ifile,'(a)')
 
45
     #     '  <!--'
 
46
      write(ifile,'(a)')
 
47
     #     MonteCarlo
 
48
      write(ifile,'(a)')
 
49
     #     '  -->'
 
50
      write(ifile,'(a)')
 
51
     #     '  <header>'
 
52
      write(ifile,'(a)')string
 
53
      write(ifile,'(a)')
 
54
     #     '  </header>'
 
55
      return
 
56
      end
 
57
 
 
58
 
 
59
 
 
60
      subroutine write_lhef_header_banner(ifile,nevents,MonteCarlo,path)
 
61
      implicit none 
 
62
      integer ifile,nevents
 
63
      character*10 MonteCarlo
 
64
      character*100 path
 
65
      character*72 buffer
 
66
c
 
67
      write(ifile,'(a)') '<LesHouchesEvents version="1.0">'
 
68
      write(ifile,'(a)') '  <!--'
 
69
      write(ifile,'(a)') MonteCarlo
 
70
      write(ifile,'(a)') '  -->'
 
71
      write(ifile,'(a)') '  <header>'
 
72
      write(ifile,250) nevents
 
73
      write(ifile,'(a)') '  <MG5ProcCard>'
 
74
      open (unit=71,file=path(1:index(path," ")-1)//'proc_card_mg5.dat'
 
75
     &     ,err=99)
 
76
      do
 
77
         read(71,'(a)',err=89,end=89) buffer
 
78
         write(ifile,'(a)') buffer
 
79
      enddo
 
80
 89   close(71)
 
81
      write(ifile,'(a)') '  </MG5ProcCard>'
 
82
      write(ifile,'(a)') '  <slha>'
 
83
      open (unit=71,file=path(1:index(path," ")-1)//'param_card.dat'
 
84
     &     ,err=98)
 
85
      do
 
86
         read(71,'(a)',err=88,end=88) buffer
 
87
         write(ifile,'(a)') buffer
 
88
      enddo
 
89
 88   close(71)
 
90
      write(ifile,'(a)') '  </slha>'
 
91
      write(ifile,'(a)') '  <MG5RunCard>'
 
92
      open (unit=71,file=path(1:index(path," ")-1)//'run_card.dat'
 
93
     &     ,err=97)
 
94
      do
 
95
         read(71,'(a)',err=87,end=87) buffer
 
96
         write(ifile,'(a)') buffer
 
97
      enddo
 
98
 87   close(71)
 
99
      write(ifile,'(a)') '  </MG5RunCard>'
 
100
      write(ifile,250) nevents
 
101
      write(ifile,'(a)') '  </header>'
 
102
 250  format(1x,i8)
 
103
      return
 
104
 99   write (*,*) 'ERROR in write_lhef_header_banner: '/
 
105
     &     /' proc_card_mg5.dat not found   :',path(1:index(path," ")-1)
 
106
     &     //'proc_card_mg5.dat'
 
107
      stop
 
108
 98   write (*,*) 'ERROR in write_lhef_header_banner: '/
 
109
     &     /' param_card.dat not found   :',path(1:index(path," ")-1)
 
110
     &     //'param_card.dat'
 
111
      stop
 
112
 97   write (*,*) 'ERROR in write_lhef_header_banner: '/
 
113
     &     /' run_card.dat not found   :',path(1:index(path," ")-1)
 
114
     &     //'run_card.dat'
 
115
      stop
 
116
      end
 
117
 
 
118
 
 
119
 
 
120
      subroutine read_lhef_header(ifile,nevents,MonteCarlo)
 
121
      implicit none 
 
122
      integer ifile,nevents
 
123
      character*10 MonteCarlo
 
124
      character*80 string,string0
 
125
c
 
126
      string='  '
 
127
      dowhile(string.ne.'  -->')
 
128
        string0=string
 
129
        read(ifile,'(a)')string
 
130
      enddo
 
131
c Works only if the name of the MC is the last line of the comments
 
132
      MonteCarlo=string0(1:10)
 
133
c Here we are at the end of (user-defined) comments. Now go to end
 
134
c of headers
 
135
      dowhile(string.ne.'  </header>')
 
136
        string0=string
 
137
        read(ifile,'(a)')string
 
138
      enddo
 
139
      read(string0,250)nevents
 
140
 250  format(1x,i8)
 
141
      return
 
142
      end
 
143
 
 
144
 
 
145
      subroutine write_lhef_init(ifile,
 
146
     #  IDBMUP,EBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,
 
147
     #  XSECUP,XERRUP,XMAXUP,LPRUP)
 
148
      implicit none
 
149
      integer ifile,IDBMUP(2),PDFGUP(2),PDFSUP(2),IDWTUP,NPRUP,LPRUP
 
150
      double precision EBMUP(2),XSECUP,XERRUP,XMAXUP
 
151
c
 
152
      write(ifile,'(a)')
 
153
     # '  <init>'
 
154
      write(ifile,501)IDBMUP(1),IDBMUP(2),EBMUP(1),EBMUP(2),
 
155
     #                PDFGUP(1),PDFGUP(2),PDFSUP(1),PDFSUP(2),
 
156
     #                IDWTUP,NPRUP
 
157
      write(ifile,502)XSECUP,XERRUP,XMAXUP,LPRUP
 
158
      write(ifile,'(a)')
 
159
     # '  </init>'
 
160
 501  format(2(1x,i6),2(1x,d14.8),2(1x,i2),2(1x,i6),1x,i2,1x,i3)
 
161
 502  format(3(1x,d14.8),1x,i6)
 
162
c
 
163
      return
 
164
      end
 
165
 
 
166
 
 
167
      subroutine read_lhef_init(ifile,
 
168
     #  IDBMUP,EBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,
 
169
     #  XSECUP,XERRUP,XMAXUP,LPRUP)
 
170
      implicit none
 
171
      integer ifile,IDBMUP(2),PDFGUP(2),PDFSUP(2),IDWTUP,NPRUP,LPRUP
 
172
      double precision EBMUP(2),XSECUP,XERRUP,XMAXUP
 
173
      character*80 string
 
174
c
 
175
      read(ifile,'(a)')string
 
176
      read(ifile,501)IDBMUP(1),IDBMUP(2),EBMUP(1),EBMUP(2),
 
177
     #                PDFGUP(1),PDFGUP(2),PDFSUP(1),PDFSUP(2),
 
178
     #                IDWTUP,NPRUP
 
179
      read(ifile,502)XSECUP,XERRUP,XMAXUP,LPRUP
 
180
      read(ifile,'(a)')string
 
181
 501  format(2(1x,i6),2(1x,d14.8),2(1x,i2),2(1x,i6),1x,i2,1x,i3)
 
182
 502  format(3(1x,d14.8),1x,i6)
 
183
c
 
184
      return
 
185
      end
 
186
 
 
187
 
 
188
      subroutine write_lhef_event(ifile,
 
189
     # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
190
     # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
191
      implicit none
 
192
      INTEGER NUP,IDPRUP,IDUP(*),ISTUP(*),MOTHUP(2,*),ICOLUP(2,*)
 
193
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
194
     # PUP(5,*),VTIMUP(*),SPINUP(*)
 
195
      character*140 buff
 
196
      integer ifile,i
 
197
      character*1 ch1
 
198
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
 
199
      double precision scale1_lhe,scale2_lhe
 
200
      integer ii,j,nps,nng,iFKS
 
201
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
202
      include 'reweight_all.inc'
 
203
c
 
204
      write(ifile,'(a)')
 
205
     # '  <event>'
 
206
      write(ifile,503)NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
 
207
      do i=1,nup
 
208
        write(ifile,504)IDUP(I),ISTUP(I),MOTHUP(1,I),MOTHUP(2,I),
 
209
     #                  ICOLUP(1,I),ICOLUP(2,I),
 
210
     #                  PUP(1,I),PUP(2,I),PUP(3,I),PUP(4,I),PUP(5,I),
 
211
     #                  VTIMUP(I),SPINUP(I)
 
212
      enddo
 
213
      if(buff(1:1).eq.'#') then
 
214
        write(ifile,'(a)') buff(1:len_trim(buff))
 
215
        read(buff,200)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
 
216
     #                    fksfather_lhe,ipartner_lhe,
 
217
     #                    scale1_lhe,scale2_lhe,
 
218
     #                    jwgtinfo,mexternal,iwgtnumpartn,
 
219
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
220
        if(jwgtinfo.ge.1.and.jwgtinfo.le.4)then
 
221
          write(ifile,'(a)')
 
222
     # '  <rwgt>'
 
223
          write(ifile,401)wgtref,wgtqes2(2)
 
224
          write(ifile,402)wgtxbj(1,1),wgtxbj(2,1),
 
225
     #                    wgtxbj(1,2),wgtxbj(2,2),
 
226
     #                    wgtxbj(1,3),wgtxbj(2,3),
 
227
     #                    wgtxbj(1,4),wgtxbj(2,4)
 
228
          if(jwgtinfo.eq.1)then
 
229
            write(ifile,403)wgtmuR2(1),wgtmuF12(1),wgtmuF22(1),
 
230
     #                      wgtmuR2(2),wgtmuF12(2),wgtmuF22(2)
 
231
          elseif(jwgtinfo.eq.2)then
 
232
            ii=iSorH_lhe+1
 
233
            if(ii.eq.3)ii=1
 
234
            write(ifile,404)wgtmuR2(ii),wgtmuF12(ii),wgtmuF22(ii)
 
235
            do i=1,mexternal
 
236
              write(ifile,405)(wgtkinE(j,i,iSorH_lhe),j=0,3)
 
237
            enddo
 
238
          elseif(jwgtinfo.eq.3 .or. jwgtinfo.eq.4)then
 
239
            do i=1,mexternal
 
240
              write(ifile,405)(wgtkinE(j,i,1),j=0,3)
 
241
            enddo
 
242
            do i=1,mexternal
 
243
              write(ifile,405)(wgtkinE(j,i,2),j=0,3)
 
244
            enddo
 
245
          endif
 
246
          write(ifile,441)wgtwreal(1),wgtwreal(2),
 
247
     #                    wgtwreal(3),wgtwreal(4)
 
248
          write(ifile,441)wgtwdeg(3),wgtwdeg(4),
 
249
     #                    wgtwdegmuf(3),wgtwdegmuf(4)
 
250
          write(ifile,405)wgtwborn(2),wgtwns(2),
 
251
     #                    wgtwnsmuf(2),wgtwnsmur(2)
 
252
          do i=1,iwgtnumpartn
 
253
            write(ifile,442)wgtwmcxsecE(i),
 
254
     #                      wgtmcxbjE(1,i),wgtmcxbjE(2,i)
 
255
          enddo
 
256
          if(jwgtinfo.eq.4) write(ifile,'(1x,d14.8,1x,i4,1x,i4)')
 
257
     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
 
258
          write(ifile,'(a)')
 
259
     # '  </rwgt>'
 
260
         elseif(jwgtinfo.eq.5) then
 
261
           write(ifile,'(a)')'  <rwgt>'
 
262
           if (iSorH_lhe.eq.1) then ! S-event
 
263
              write(ifile,'(1x,d14.8,i4)') wgtbpower,nScontributions
 
264
              write(ifile,'(1x,i4,1x,d14.8)') nFKSprocess_used_born
 
265
     &             ,wgtref_nbody
 
266
              do i=1,mexternal
 
267
                 write(ifile,405)(wgtkin_all(j,i,2,0),j=0,3)
 
268
              enddo
 
269
              write(ifile,402) wgtxbj_all(1,2,0),wgtxbj_all(2,2,0)
 
270
              write(ifile,'(1x,d14.8)') wgtqes2_all(2,0)
 
271
              write(ifile,405)wgtwborn_all,wgtwns_all,
 
272
     &             wgtwnsmuf_all,wgtwnsmur_all
 
273
              
 
274
              do ii=1,nScontributions
 
275
                 write(ifile,'(1x,i4)') nFKSprocess_reweight(ii)
 
276
                 iFKS=nFKSprocess_reweight(ii)*2-1
 
277
                 write(ifile,'(1x,d14.8,1x,i4)')
 
278
     &                wgtref_all(iFKS),iwgtnumpartn_all(iFKS)
 
279
                 do i=1,mexternal
 
280
                    write(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
 
281
                 enddo
 
282
c$$$                 do i=1,mexternal
 
283
c$$$                    write(ifile,405)(wgtkin_all(j,i,2,iFKS),j=0,3)
 
284
c$$$                 enddo
 
285
                 write(ifile,402)
 
286
     &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
 
287
     &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
 
288
     &                wgtxbj_all(1,3,iFKS),wgtxbj_all(2,3,iFKS),
 
289
     &                wgtxbj_all(1,4,iFKS),wgtxbj_all(2,4,iFKS)
 
290
                 write(ifile,'(1x,d14.8)') wgtqes2_all(2,iFKS)
 
291
                 write(ifile,441)wgtwreal_all(1,iFKS),wgtwreal_all(2
 
292
     &                ,iFKS),wgtwreal_all(3,iFKS),wgtwreal_all(4,iFKS)
 
293
                 write(ifile,441)wgtwdeg_all(3,iFKS),wgtwdeg_all(4,iFKS)
 
294
     &                ,wgtwdegmuf_all(3,iFKS),wgtwdegmuf_all(4,iFKS)
 
295
                 do i=1,iwgtnumpartn_all(iFKS)
 
296
                    write(ifile,442)wgtwmcxsec_all(i,iFKS),
 
297
     &                   wgtmcxbj_all(1,i,iFKS),wgtmcxbj_all(2,i,iFKS)
 
298
                 enddo
 
299
                 
 
300
              enddo
 
301
           elseif (iSorH_lhe.eq.2) then ! H-event
 
302
              write(ifile,'(1x,d14.8)') wgtbpower
 
303
              iFKS=nFKSprocess_used*2
 
304
              write(ifile,'(1x,i4)') nFKSprocess_used
 
305
              write(ifile,'(1x,d14.8,1x,i4)') wgtref_all(iFKS)
 
306
     &             ,iwgtnumpartn_all(iFKS)
 
307
              do i=1,mexternal
 
308
                 write(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
 
309
              enddo
 
310
              do i=1,mexternal
 
311
                 write(ifile,405)(wgtkin_all(j,i,2,iFKS),j=0,3)
 
312
              enddo
 
313
              write(ifile,402)
 
314
     &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
 
315
     &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
 
316
     &                wgtxbj_all(1,3,iFKS),wgtxbj_all(2,3,iFKS),
 
317
     &                wgtxbj_all(1,4,iFKS),wgtxbj_all(2,4,iFKS)
 
318
                 write(ifile,441)wgtwreal_all(1,iFKS),wgtwreal_all(2
 
319
     &                ,iFKS),wgtwreal_all(3,iFKS),wgtwreal_all(4,iFKS)
 
320
                 do i=1,iwgtnumpartn_all(iFKS)
 
321
                    write(ifile,442)wgtwmcxsec_all(i,iFKS),
 
322
     &                   wgtmcxbj_all(1,i,iFKS),wgtmcxbj_all(2,i,iFKS)
 
323
                 enddo
 
324
           else
 
325
              write (*,*) 'Not an S- or H-event in write_lhef_event'
 
326
              stop
 
327
           endif
 
328
           write(ifile,'(a)')'  </rwgt>'
 
329
 
 
330
        elseif(jwgtinfo.eq.8)then
 
331
          write(ifile,'(a)')
 
332
     # '  <rwgt>'
 
333
          write(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
 
334
          do i=1,numscales
 
335
            write(ifile,404)(wgtxsecmu(i,j),j=1,numscales)
 
336
          enddo
 
337
          do i=1,numPDFpairs
 
338
            nps=2*i-1
 
339
            nng=2*i
 
340
            write(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
 
341
          enddo
 
342
          write(ifile,'(a)')
 
343
     # '  </rwgt>'
 
344
        endif
 
345
      endif
 
346
      write(ifile,'(a)')
 
347
     # '  </event>'
 
348
 200  format(1a,1x,i1,4(1x,i2),2(1x,d14.8),1x,i1,2(1x,i2),5(1x,d14.8))
 
349
 401  format(2(1x,d14.8))
 
350
 402  format(8(1x,d14.8))
 
351
 403  format(6(1x,d14.8))
 
352
 404  format(3(1x,d14.8))
 
353
 405  format(4(1x,d14.8))
 
354
 406  format(2(1x,d14.8),2(1x,i3))
 
355
 441  format(4(1x,d16.10))
 
356
 442  format(1x,d16.10,2(1x,d14.8))
 
357
 503  format(1x,i2,1x,i6,4(1x,d14.8))
 
358
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,d14.8),2(1x,d10.4))
 
359
c
 
360
      return
 
361
      end
 
362
 
 
363
 
 
364
      subroutine read_lhef_event(ifile,
 
365
     # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
366
     # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
367
      implicit none
 
368
      INTEGER NUP,IDPRUP,IDUP(*),ISTUP(*),MOTHUP(2,*),ICOLUP(2,*)
 
369
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
370
     # PUP(5,*),VTIMUP(*),SPINUP(*)
 
371
      integer ifile,i
 
372
      character*140 buff
 
373
      character*80 string
 
374
      character*1 ch1
 
375
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
 
376
      double precision scale1_lhe,scale2_lhe
 
377
      integer ii,j,nps,nng,iFKS
 
378
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
379
      include 'reweight_all.inc'
 
380
c
 
381
      read(ifile,'(a)')string
 
382
      read(ifile,*)NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
 
383
      do i=1,nup
 
384
        read(ifile,*)IDUP(I),ISTUP(I),MOTHUP(1,I),MOTHUP(2,I),
 
385
     #                  ICOLUP(1,I),ICOLUP(2,I),
 
386
     #                  PUP(1,I),PUP(2,I),PUP(3,I),PUP(4,I),PUP(5,I),
 
387
     #                  VTIMUP(I),SPINUP(I)
 
388
      enddo
 
389
      read(ifile,'(a)')buff
 
390
      if(buff(1:1).eq.'#')then
 
391
        read(buff,200)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
 
392
     #                    fksfather_lhe,ipartner_lhe,
 
393
     #                    scale1_lhe,scale2_lhe,
 
394
     #                    jwgtinfo,mexternal,iwgtnumpartn,
 
395
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
396
        if(jwgtinfo.ge.1.and.jwgtinfo.le.4)then
 
397
          read(ifile,'(a)')string
 
398
          read(ifile,401)wgtref,wgtqes2(2)
 
399
          read(ifile,402)wgtxbj(1,1),wgtxbj(2,1),
 
400
     #                   wgtxbj(1,2),wgtxbj(2,2),
 
401
     #                   wgtxbj(1,3),wgtxbj(2,3),
 
402
     #                   wgtxbj(1,4),wgtxbj(2,4)
 
403
          if(jwgtinfo.eq.1)then
 
404
            read(ifile,403)wgtmuR2(1),wgtmuF12(1),wgtmuF22(1),
 
405
     #                     wgtmuR2(2),wgtmuF12(2),wgtmuF22(2)
 
406
          elseif(jwgtinfo.eq.2)then
 
407
            ii=iSorH_lhe+1
 
408
            if(ii.eq.3)ii=1
 
409
            read(ifile,404)wgtmuR2(ii),wgtmuF12(ii),wgtmuF22(ii)
 
410
            do i=1,mexternal
 
411
              read(ifile,405)(wgtkinE(j,i,iSorH_lhe),j=0,3)
 
412
            enddo
 
413
          elseif(jwgtinfo.eq.3 .or. jwgtinfo.eq.4)then
 
414
            do i=1,mexternal
 
415
              read(ifile,405)(wgtkinE(j,i,1),j=0,3)
 
416
            enddo
 
417
            do i=1,mexternal
 
418
              read(ifile,405)(wgtkinE(j,i,2),j=0,3)
 
419
            enddo
 
420
          endif
 
421
          read(ifile,441)wgtwreal(1),wgtwreal(2),
 
422
     #                   wgtwreal(3),wgtwreal(4)
 
423
          read(ifile,441)wgtwdeg(3),wgtwdeg(4),
 
424
     #                   wgtwdegmuf(3),wgtwdegmuf(4)
 
425
          read(ifile,405)wgtwborn(2),wgtwns(2),
 
426
     #                    wgtwnsmuf(2),wgtwnsmur(2)
 
427
          do i=1,iwgtnumpartn
 
428
            read(ifile,442)wgtwmcxsecE(i),
 
429
     #                     wgtmcxbjE(1,i),wgtmcxbjE(2,i)
 
430
          enddo
 
431
          if(jwgtinfo.eq.4) read(ifile,'(1x,d14.8,1x,i4,1x,i4)')
 
432
     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
 
433
          read(ifile,'(a)')string
 
434
        elseif(jwgtinfo.eq.5) then
 
435
           read(ifile,'(a)')string
 
436
           if (iSorH_lhe.eq.1) then ! S-event
 
437
              read(ifile,'(1x,d14.8,i4)') wgtbpower,nScontributions
 
438
              read(ifile,'(1x,i4,1x,d14.8)') nFKSprocess_used_born
 
439
     &             ,wgtref_nbody
 
440
              do i=1,mexternal
 
441
                 read(ifile,405)(wgtkin_all(j,i,2,0),j=0,3)
 
442
              enddo
 
443
              read(ifile,402) wgtxbj_all(1,2,0),wgtxbj_all(2,2,0)
 
444
              read(ifile,'(1x,d14.8)') wgtqes2_all(2,0)
 
445
              read(ifile,405)wgtwborn_all,wgtwns_all,
 
446
     &             wgtwnsmuf_all,wgtwnsmur_all
 
447
              
 
448
              do ii=1,nScontributions
 
449
                 read(ifile,'(1x,i4)') nFKSprocess_reweight(ii)
 
450
                 iFKS=nFKSprocess_reweight(ii)*2-1
 
451
                 read(ifile,'(1x,d14.8,1x,i4)') wgtref_all(iFKS)
 
452
     &                ,iwgtnumpartn_all(iFKS)
 
453
                 do i=1,mexternal
 
454
                    read(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
 
455
                 enddo
 
456
                 do i=1,mexternal
 
457
                    do j=0,3
 
458
                       wgtkin_all(j,i,2,iFKS)=wgtkin_all(j,i,2,0)
 
459
                    enddo
 
460
                 enddo
 
461
                 read(ifile,402)
 
462
     &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
 
463
     &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
 
464
     &                wgtxbj_all(1,3,iFKS),wgtxbj_all(2,3,iFKS),
 
465
     &                wgtxbj_all(1,4,iFKS),wgtxbj_all(2,4,iFKS)
 
466
                 read(ifile,'(1x,d14.8)') wgtqes2_all(2,iFKS)
 
467
                 read(ifile,441)wgtwreal_all(1,iFKS),wgtwreal_all(2
 
468
     &                ,iFKS),wgtwreal_all(3,iFKS),wgtwreal_all(4,iFKS)
 
469
                 read(ifile,441)wgtwdeg_all(3,iFKS),wgtwdeg_all(4,iFKS)
 
470
     &                ,wgtwdegmuf_all(3,iFKS),wgtwdegmuf_all(4,iFKS)
 
471
                 do i=1,iwgtnumpartn_all(iFKS)
 
472
                    read(ifile,442)wgtwmcxsec_all(i,iFKS),
 
473
     &                   wgtmcxbj_all(1,i,iFKS),wgtmcxbj_all(2,i,iFKS)
 
474
                 enddo
 
475
              
 
476
              enddo
 
477
           elseif (iSorH_lhe.eq.2) then ! H-event
 
478
              read(ifile,'(1x,d14.8)') wgtbpower
 
479
              read(ifile,'(1x,i4)') nFKSprocess_used
 
480
              iFKS=nFKSprocess_used*2
 
481
              read(ifile,'(1x,d14.8,1x,i4)') wgtref_all(iFKS)
 
482
     &             ,iwgtnumpartn_all(iFKS)
 
483
              do i=1,mexternal
 
484
                 read(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
 
485
              enddo
 
486
              do i=1,mexternal
 
487
                 read(ifile,405)(wgtkin_all(j,i,2,iFKS),j=0,3)
 
488
              enddo
 
489
              read(ifile,402)
 
490
     &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
 
491
     &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
 
492
     &                wgtxbj_all(1,3,iFKS),wgtxbj_all(2,3,iFKS),
 
493
     &                wgtxbj_all(1,4,iFKS),wgtxbj_all(2,4,iFKS)
 
494
                 read(ifile,441)wgtwreal_all(1,iFKS),wgtwreal_all(2
 
495
     &                ,iFKS),wgtwreal_all(3,iFKS),wgtwreal_all(4,iFKS)
 
496
                 do i=1,iwgtnumpartn_all(iFKS)
 
497
                    read(ifile,442)wgtwmcxsec_all(i,iFKS),
 
498
     &                   wgtmcxbj_all(1,i,iFKS),wgtmcxbj_all(2,i,iFKS)
 
499
                 enddo
 
500
           else
 
501
              write (*,*) 'Not an S- or H-event in write_lhef_event'
 
502
              stop
 
503
           endif
 
504
           read(ifile,'(a)')string
 
505
 
 
506
        elseif(jwgtinfo.eq.8)then
 
507
          read(ifile,'(a)')string
 
508
          read(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
 
509
          do i=1,numscales
 
510
            read(ifile,404)(wgtxsecmu(i,j),j=1,numscales)
 
511
          enddo
 
512
          do i=1,numPDFpairs
 
513
            nps=2*i-1
 
514
            nng=2*i
 
515
            read(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
 
516
          enddo
 
517
          read(ifile,'(a)')string
 
518
        endif
 
519
        read(ifile,'(a)')string
 
520
      else
 
521
        string=buff(1:len_trim(buff))
 
522
        buff=' '
 
523
      endif
 
524
 200  format(1a,1x,i1,4(1x,i2),2(1x,d14.8),1x,i1,2(1x,i2),5(1x,d14.8))
 
525
 401  format(2(1x,d14.8))
 
526
 402  format(8(1x,d14.8))
 
527
 403  format(6(1x,d14.8))
 
528
 404  format(3(1x,d14.8))
 
529
 405  format(4(1x,d14.8))
 
530
 406  format(2(1x,d14.8),2(1x,i3))
 
531
 441  format(4(1x,d16.10))
 
532
 442  format(1x,d16.10,2(1x,d14.8))
 
533
 503  format(1x,i2,1x,i6,4(1x,d14.8))
 
534
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,d14.8),2(1x,d10.4))
 
535
c
 
536
      return
 
537
      end
 
538
 
 
539
 
 
540
c Same as read_lhef_event, except for the end-of-file catch
 
541
      subroutine read_lhef_event_catch(ifile,
 
542
     # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
543
     # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
544
      implicit none
 
545
      INTEGER NUP,IDPRUP,IDUP(*),ISTUP(*),MOTHUP(2,*),ICOLUP(2,*)
 
546
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
547
     # PUP(5,*),VTIMUP(*),SPINUP(*)
 
548
      integer ifile,i
 
549
      character*140 buff
 
550
      character*80 string
 
551
      character*1 ch1
 
552
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
 
553
      double precision scale1_lhe,scale2_lhe
 
554
      integer ii,j,nps,nng,iFKS
 
555
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
556
      include 'reweight_all.inc'
 
557
c
 
558
      read(ifile,'(a)')string
 
559
      if(index(string,'<event>').eq.0)then
 
560
        if(index(string,'</LesHouchesEvents>').ne.0)then
 
561
          buff='endoffile'
 
562
          return
 
563
        else
 
564
          write(*,*)'Unknown structure in read_lhef_event_catch:'
 
565
          write(*,*)string(1:len_trim(string))
 
566
          stop
 
567
        endif
 
568
      endif
 
569
      read(ifile,*)NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
 
570
      do i=1,nup
 
571
        read(ifile,*)IDUP(I),ISTUP(I),MOTHUP(1,I),MOTHUP(2,I),
 
572
     #                  ICOLUP(1,I),ICOLUP(2,I),
 
573
     #                  PUP(1,I),PUP(2,I),PUP(3,I),PUP(4,I),PUP(5,I),
 
574
     #                  VTIMUP(I),SPINUP(I)
 
575
      enddo
 
576
      read(ifile,'(a)')buff
 
577
      if(buff(1:1).eq.'#')then
 
578
        read(buff,200)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
 
579
     #                    fksfather_lhe,ipartner_lhe,
 
580
     #                    scale1_lhe,scale2_lhe,
 
581
     #                    jwgtinfo,mexternal,iwgtnumpartn,
 
582
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
583
        if(jwgtinfo.ge.1.and.jwgtinfo.le.4)then
 
584
          read(ifile,'(a)')string
 
585
          read(ifile,401)wgtref,wgtqes2(2)
 
586
          read(ifile,402)wgtxbj(1,1),wgtxbj(2,1),
 
587
     #                   wgtxbj(1,2),wgtxbj(2,2),
 
588
     #                   wgtxbj(1,3),wgtxbj(2,3),
 
589
     #                   wgtxbj(1,4),wgtxbj(2,4)
 
590
          if(jwgtinfo.eq.1)then
 
591
            read(ifile,403)wgtmuR2(1),wgtmuF12(1),wgtmuF22(1),
 
592
     #                     wgtmuR2(2),wgtmuF12(2),wgtmuF22(2)
 
593
          elseif(jwgtinfo.eq.2)then
 
594
            ii=iSorH_lhe+1
 
595
            if(ii.eq.3)ii=1
 
596
            read(ifile,404)wgtmuR2(ii),wgtmuF12(ii),wgtmuF22(ii)
 
597
            do i=1,mexternal
 
598
              read(ifile,405)(wgtkinE(j,i,iSorH_lhe),j=0,3)
 
599
            enddo
 
600
          elseif(jwgtinfo.eq.3 .or. jwgtinfo.eq.4)then
 
601
            do i=1,mexternal
 
602
              read(ifile,405)(wgtkinE(j,i,1),j=0,3)
 
603
            enddo
 
604
            do i=1,mexternal
 
605
              read(ifile,405)(wgtkinE(j,i,2),j=0,3)
 
606
            enddo
 
607
          endif
 
608
          read(ifile,441)wgtwreal(1),wgtwreal(2),
 
609
     #                   wgtwreal(3),wgtwreal(4)
 
610
          read(ifile,441)wgtwdeg(3),wgtwdeg(4),
 
611
     #                   wgtwdegmuf(3),wgtwdegmuf(4)
 
612
          read(ifile,405)wgtwborn(2),wgtwns(2),
 
613
     #                    wgtwnsmuf(2),wgtwnsmur(2)
 
614
          do i=1,iwgtnumpartn
 
615
            read(ifile,442)wgtwmcxsecE(i),
 
616
     #                     wgtmcxbjE(1,i),wgtmcxbjE(2,i)
 
617
          enddo
 
618
          if(jwgtinfo.eq.4) read(ifile,'(1x,d14.8,1x,i4,1x,i4)')
 
619
     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
 
620
          read(ifile,'(a)')string
 
621
        elseif(jwgtinfo.eq.5) then
 
622
           read(ifile,'(a)')string
 
623
           if (iSorH_lhe.eq.1) then ! S-event
 
624
              read(ifile,'(1x,d14.8,i4)') wgtbpower,nScontributions
 
625
              read(ifile,'(1x,i4,1x,d14.8)') nFKSprocess_used_born
 
626
     &             ,wgtref_nbody
 
627
              do i=1,mexternal
 
628
                 read(ifile,405)(wgtkin_all(j,i,2,0),j=0,3)
 
629
              enddo
 
630
              read(ifile,402) wgtxbj_all(1,2,0),wgtxbj_all(2,2,0)
 
631
              read(ifile,'(1x,d14.8)') wgtqes2_all(2,0)
 
632
              read(ifile,405)wgtwborn_all,wgtwns_all,
 
633
     &             wgtwnsmuf_all,wgtwnsmur_all
 
634
              
 
635
              do ii=1,nScontributions
 
636
                 read(ifile,'(1x,i4)') nFKSprocess_reweight(ii)
 
637
                 iFKS=nFKSprocess_reweight(ii)*2-1
 
638
                 read(ifile,'(1x,d14.8,1x,i4)') wgtref_all(iFKS)
 
639
     &                ,iwgtnumpartn_all(iFKS)
 
640
                 do i=1,mexternal
 
641
                    read(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
 
642
                 enddo
 
643
                 do i=1,mexternal
 
644
                    do j=0,3
 
645
                       wgtkin_all(j,i,2,iFKS)=wgtkin_all(j,i,2,0)
 
646
                    enddo
 
647
                 enddo
 
648
                 read(ifile,402)
 
649
     &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
 
650
     &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
 
651
     &                wgtxbj_all(1,3,iFKS),wgtxbj_all(2,3,iFKS),
 
652
     &                wgtxbj_all(1,4,iFKS),wgtxbj_all(2,4,iFKS)
 
653
                 read(ifile,'(1x,d14.8)') wgtqes2_all(2,iFKS)
 
654
                 read(ifile,441)wgtwreal_all(1,iFKS),wgtwreal_all(2
 
655
     &                ,iFKS),wgtwreal_all(3,iFKS),wgtwreal_all(4,iFKS)
 
656
                 read(ifile,441)wgtwdeg_all(3,iFKS),wgtwdeg_all(4,iFKS)
 
657
     &                ,wgtwdegmuf_all(3,iFKS),wgtwdegmuf_all(4,iFKS)
 
658
                 do i=1,iwgtnumpartn_all(iFKS)
 
659
                    read(ifile,442)wgtwmcxsec_all(i,iFKS),
 
660
     &                   wgtmcxbj_all(1,i,iFKS),wgtmcxbj_all(2,i,iFKS)
 
661
                 enddo
 
662
                 
 
663
              enddo
 
664
           elseif (iSorH_lhe.eq.2) then ! H-event
 
665
              read(ifile,'(1x,d14.8)') wgtbpower
 
666
              read(ifile,'(1x,i4)') nFKSprocess_used
 
667
              iFKS=nFKSprocess_used*2
 
668
              read(ifile,'(1x,d14.8,1x,i4)') wgtref_all(iFKS)
 
669
     &             ,iwgtnumpartn_all(iFKS)
 
670
              do i=1,mexternal
 
671
                 read(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
 
672
              enddo
 
673
              do i=1,mexternal
 
674
                 read(ifile,405)(wgtkin_all(j,i,2,iFKS),j=0,3)
 
675
              enddo
 
676
              read(ifile,402)
 
677
     &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
 
678
     &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
 
679
     &                wgtxbj_all(1,3,iFKS),wgtxbj_all(2,3,iFKS),
 
680
     &                wgtxbj_all(1,4,iFKS),wgtxbj_all(2,4,iFKS)
 
681
                 read(ifile,441)wgtwreal_all(1,iFKS),wgtwreal_all(2
 
682
     &                ,iFKS),wgtwreal_all(3,iFKS),wgtwreal_all(4,iFKS)
 
683
                 do i=1,iwgtnumpartn_all(iFKS)
 
684
                    read(ifile,442)wgtwmcxsec_all(i,iFKS),
 
685
     &                   wgtmcxbj_all(1,i,iFKS),wgtmcxbj_all(2,i,iFKS)
 
686
                 enddo
 
687
           else
 
688
              write (*,*) 'Not an S- or H-event in write_lhef_event'
 
689
              stop
 
690
           endif
 
691
           read(ifile,'(a)')string
 
692
 
 
693
        elseif(jwgtinfo.eq.8)then
 
694
          read(ifile,'(a)')string
 
695
          read(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
 
696
          do i=1,numscales
 
697
            read(ifile,404)(wgtxsecmu(i,j),j=1,numscales)
 
698
          enddo
 
699
          do i=1,numPDFpairs
 
700
            nps=2*i-1
 
701
            nng=2*i
 
702
            read(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
 
703
          enddo
 
704
          read(ifile,'(a)')string
 
705
        endif
 
706
        read(ifile,'(a)')string
 
707
      else
 
708
        string=buff(1:len_trim(buff))
 
709
        buff=' '
 
710
      endif
 
711
 200  format(1a,1x,i1,4(1x,i2),2(1x,d14.8),1x,i1,2(1x,i2),5(1x,d14.8))
 
712
 401  format(2(1x,d14.8))
 
713
 402  format(8(1x,d14.8))
 
714
 403  format(6(1x,d14.8))
 
715
 404  format(3(1x,d14.8))
 
716
 405  format(4(1x,d14.8))
 
717
 406  format(2(1x,d14.8),2(1x,i3))
 
718
 441  format(4(1x,d16.10))
 
719
 442  format(1x,d16.10,2(1x,d14.8))
 
720
 503  format(1x,i2,1x,i6,4(1x,d14.8))
 
721
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,d14.8),2(1x,d10.4))
 
722
c
 
723
      return
 
724
      end