← 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

merge with 2.0.0beta4 revision 216

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/handling_lhe_events.f
10
10
c Hence, one may need an explicit copy of one onto the other
11
11
c
12
12
 
 
13
      block data
 
14
      integer event_id
 
15
      common /c_event_id/ event_id
 
16
      data event_id /-99/
 
17
      logical rwgt_skip
 
18
      common /crwgt_skip/ rwgt_skip
 
19
      data rwgt_skip /.false./
 
20
      end
 
21
 
13
22
      subroutine write_lhef_header(ifile,nevents,MonteCarlo)
14
23
      implicit none 
15
24
      integer ifile,nevents
35
44
 
36
45
      subroutine write_lhef_header_banner(ifile,nevents,MonteCarlo,path)
37
46
      implicit none 
38
 
      integer ifile,nevents,iseed,i
39
 
      double precision mcmass(-5:21)
 
47
      integer ifile,nevents,iseed,i,pdf_set_min,pdf_set_max,idwgt
 
48
      double precision mcmass(-5:21),rw_Rscale_down,rw_Rscale_up
 
49
     $     ,rw_Fscale_down,rw_Fscale_up
40
50
      character*10 MonteCarlo
41
51
      character*100 path
42
 
      character*72 buffer,buffer2
 
52
      character*72 buffer,buffer_lc,buffer2
 
53
      logical rwgt_skip
 
54
      common /crwgt_skip/ rwgt_skip
 
55
      logical rwgt_skip_pdf,rwgt_skip_scales
 
56
      integer event_id
 
57
      common /c_event_id/ event_id
 
58
      include 'reweight_all.inc'
 
59
c     Set the event_id to 0. If 0 or positive, this value will be update
 
60
c     in write_lhe_event. It is set to -99 through a block data
 
61
c     statement.
 
62
      event_id=0
43
63
c
44
64
      write(ifile,'(a)') '<LesHouchesEvents version="1.0">'
45
65
      write(ifile,'(a)') '  <header>'
64
84
      write(ifile,'(a)') '  <MGRunCard>'
65
85
      open (unit=92,file=path(1:index(path," ")-1)//'run_card.dat'
66
86
     &     ,err=97)
 
87
      rwgt_skip_pdf=.true.
 
88
      rwgt_skip_scales=.true.
 
89
      rwgt_skip=.true.
 
90
      pdf_set_min=-1
 
91
      pdf_set_max=-1
 
92
      numscales=0
67
93
      do
68
94
         read(92,'(a)',err=87,end=87) buffer
 
95
         buffer_lc=buffer
 
96
         call case_trap3(72,buffer_lc)
69
97
c Replace the random number seed with the one used...
70
 
         if (index(buffer,'iseed').ne.0 .and. buffer(1:1).ne.'#') then
 
98
         if (index(buffer_lc,'iseed').ne.0 .and. buffer(1:1).ne.'#') then
71
99
            open (unit=93,file="randinit",status="old",err=96)
72
100
            read(93,'(a)') buffer2
73
101
            if (index(buffer2,'=').eq.0) goto 96
76
104
            close(93)
77
105
            write(buffer,'(i11,a)')iseed,' =  iseed'
78
106
c Update the number of events
79
 
         elseif (index(buffer,'nevents').ne.0 .and.
 
107
         elseif (index(buffer_lc,'nevents').ne.0 .and.
80
108
     &           buffer(1:1).ne.'#') then
81
109
            write(buffer,'(i11,a)')nevents,' = nevents'
 
110
         elseif (index(buffer_lc,'reweight_pdf').ne.0 .and.
 
111
     $           index(buffer_lc,'.true.').ne.0 .and.
 
112
     $           buffer(1:1).ne.'#') then
 
113
            rwgt_skip=.false.
 
114
            rwgt_skip_pdf=.false.
 
115
         elseif (index(buffer_lc,'pdf_set_min').ne.0 .and.
 
116
     &           buffer(1:1).ne.'#') then
 
117
            read(buffer(1:index(buffer,'=')-1),*) pdf_set_min
 
118
         elseif (index(buffer_lc,'pdf_set_max').ne.0 .and.
 
119
     &           buffer(1:1).ne.'#') then
 
120
            read(buffer(1:index(buffer,'=')-1),*) pdf_set_max
 
121
         elseif (index(buffer_lc,'reweight_scale').ne.0 .and.
 
122
     $           index(buffer_lc,'.true.').ne.0 .and.
 
123
     $           buffer(1:1).ne.'#') then
 
124
            rwgt_skip=.false.
 
125
            rwgt_skip_scales=.false.
 
126
            numscales=3
 
127
         elseif (index(buffer_lc,'rw_rscale_down').ne.0 .and.
 
128
     &           buffer(1:1).ne.'#') then
 
129
            read(buffer(1:index(buffer,'=')-1),*) rw_Rscale_down
 
130
         elseif (index(buffer_lc,'rw_rscale_up').ne.0 .and.
 
131
     &           buffer(1:1).ne.'#') then
 
132
            read(buffer(1:index(buffer,'=')-1),*) rw_Rscale_up
 
133
         elseif (index(buffer_lc,'rw_fscale_down').ne.0 .and.
 
134
     &           buffer(1:1).ne.'#') then
 
135
            read(buffer(1:index(buffer,'=')-1),*) rw_Fscale_down
 
136
         elseif (index(buffer_lc,'rw_fscale_up').ne.0 .and.
 
137
     &           buffer(1:1).ne.'#') then
 
138
            read(buffer(1:index(buffer,'=')-1),*) rw_Fscale_up
82
139
         endif
83
140
         goto 95
84
141
 96      write (*,*) '"randinit" file not found in write_lhef_header_'/
86
143
 95      write(ifile,'(a)') buffer
87
144
      enddo
88
145
 87   close(92)
 
146
      if ( (pdf_set_min.ne.-1 .and. pdf_set_max.eq.-1) .or.
 
147
     &     (pdf_set_min.eq.-1 .and. pdf_set_max.ne.-1)) then
 
148
         write (*,*) 'Not consistent PDF reweigthing parameters'/
 
149
     $        /' found in the run_card.',pdf_set_min,pdf_set_max
 
150
         stop
 
151
      endif
 
152
      if (.not.rwgt_skip_pdf) numPDFpairs=(pdf_set_max-pdf_set_min+1)/2
89
153
      write(ifile,'(a)') '  </MGRunCard>'
90
154
      write(ifile,'(a)') '  <MonteCarloMasses>'
91
155
      call fill_MC_mshell_wrap(MonteCarlo,mcmass)
94
158
      enddo
95
159
      write (ifile,'(2x,i6,3x,e12.6)')21,mcmass(21)
96
160
      write(ifile,'(a)') '  </MonteCarloMasses>'
 
161
c Write here the reweight information if need be
 
162
      if (.not.rwgt_skip) then
 
163
         write(ifile,'(a)') '  <initrwgt>'
 
164
         if (numscales.ne.0) then
 
165
            if (numscales.ne.3) then
 
166
               write (*,*) 'Error in handling_lhe_events.f:'
 
167
               write (*,*) 'number of scale not equal to three'
 
168
     $              ,numscales
 
169
               stop
 
170
            endif
 
171
            write(ifile,'(a)') "    <weightgroup "/
 
172
     &           /"type='scale_variation' combine='envelope'>"
 
173
            write(ifile,602) "      <weight id='1001'>"/
 
174
     $           /" muR=",1d0," muF=",1d0," </weight>"
 
175
            write(ifile,602) "      <weight id='1002'>"/
 
176
     $           /" muR=",1d0," muF=",rw_Fscale_up," </weight>"
 
177
            write(ifile,602) "      <weight id='1003'>"/
 
178
     $           /" muR=",1d0," muF=",rw_Fscale_down," </weight>"
 
179
            write(ifile,602) "      <weight id='1004'>"/
 
180
     $           /" muR=",rw_Rscale_up," muF=",1d0," </weight>"
 
181
            write(ifile,602) "      <weight id='1005'> muR="
 
182
     $           ,rw_Rscale_up," muF=",rw_Fscale_up," </weight>"
 
183
            write(ifile,602) "      <weight id='1006'> muR="
 
184
     $           ,rw_Rscale_up," muF=",rw_Fscale_down," </weight>"
 
185
            write(ifile,602) "      <weight id='1007'> muR="
 
186
     $           ,rw_Rscale_down," muF=",1d0," </weight>"
 
187
            write(ifile,602) "      <weight id='1008'> muR="
 
188
     $           ,rw_Rscale_down," muF=",rw_Fscale_up," </weight>"
 
189
            write(ifile,602) "      <weight id='1009'> muR="
 
190
     $           ,rw_Rscale_down," muF=",rw_Fscale_down," </weight>"
 
191
            write(ifile,'(a)') "    </weightgroup>"
 
192
         endif
 
193
         if (numPDFpairs.ne.0) then
 
194
            if (pdf_set_min.lt.90000) then    ! MSTW & CTEQ
 
195
               write(ifile,'(a)') "    <weightgroup "/
 
196
     &              /"type='PDF_variation' combine='hessian'>"
 
197
            else                              ! NNPDF
 
198
               write(ifile,'(a)') "    <weightgroup "/
 
199
     &              /"type='PDF_variation' combine='gaussian'>"
 
200
            endif
 
201
            do i=1,numPDFpairs*2
 
202
               idwgt=2000+i
 
203
               write(ifile,'(a,i4,a,i6,a)') "      <weight id='",idwgt,
 
204
     $              "'> pdfset=",pdf_set_min+(i-1)," </weight>"
 
205
            enddo
 
206
            write(ifile,'(a)') "    </weightgroup>"
 
207
         endif
 
208
         write(ifile,'(a)') '  </initrwgt>'
 
209
      endif
97
210
      write(ifile,'(a)') '  </header>'
98
211
 250  format(1x,i8)
99
212
      return
109
222
     &     /' run_card.dat not found   :',path(1:index(path," ")-1)
110
223
     &     //'run_card.dat'
111
224
      stop
 
225
 602  format(a,e11.5,a,e11.5,a)
112
226
      end
113
227
 
114
228
 
256
370
      return
257
371
      end
258
372
 
259
 
 
260
373
      subroutine write_lhef_event(ifile,
261
374
     # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
262
375
     # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
269
382
      character*1 ch1
270
383
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
271
384
      double precision scale1_lhe,scale2_lhe
272
 
      integer ii,j,nps,nng,iFKS
 
385
      integer ii,j,nps,nng,iFKS,idwgt
273
386
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
387
      logical rwgt_skip
 
388
      common /crwgt_skip/ rwgt_skip
 
389
      integer event_id
 
390
      common /c_event_id/ event_id
274
391
      include 'reweight_all.inc'
 
392
c     if event_id is zero or positive (that means that there was a call
 
393
c     to write_lhef_header_banner) update it and write it
 
394
c RF: don't use the event_id:
 
395
      event_id = -99
275
396
c
276
 
      write(ifile,'(a)')
277
 
     # '  <event>'
 
397
      if (event_id.ge.0) then
 
398
         event_id=event_id+1
 
399
         if (event_id.le.9) then
 
400
            write(ifile,'(a,i1,a)') "  <event id='",event_id,"'>"
 
401
         elseif(event_id.le.99) then
 
402
            write(ifile,'(a,i2,a)') "  <event id='",event_id,"'>"
 
403
         elseif(event_id.le.999) then
 
404
            write(ifile,'(a,i3,a)') "  <event id='",event_id,"'>"
 
405
         elseif(event_id.le.9999) then
 
406
            write(ifile,'(a,i4,a)') "  <event id='",event_id,"'>"
 
407
         elseif(event_id.le.99999) then
 
408
            write(ifile,'(a,i5,a)') "  <event id='",event_id,"'>"
 
409
         elseif(event_id.le.999999) then
 
410
            write(ifile,'(a,i6,a)') "  <event id='",event_id,"'>"
 
411
         elseif(event_id.le.9999999) then
 
412
            write(ifile,'(a,i7,a)') "  <event id='",event_id,"'>"
 
413
         elseif(event_id.le.99999999) then
 
414
            write(ifile,'(a,i8,a)') "  <event id='",event_id,"'>"
 
415
         elseif(event_id.le.999999999) then
 
416
            write(ifile,'(a,i9,a)') "  <event id='",event_id,"'>"
 
417
         else
 
418
            write (ifile,*) "ERROR: EVENT ID TOO LARGE",event_id
 
419
            write (*,*) "ERROR: EVENT ID TOO LARGE",event_id
 
420
            stop
 
421
         endif
 
422
      else
 
423
         write(ifile,'(a)') '  <event>'
 
424
      endif
278
425
      write(ifile,503)NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
279
426
      do i=1,nup
280
427
        write(ifile,504)IDUP(I),ISTUP(I),MOTHUP(1,I),MOTHUP(2,I),
282
429
     #                  PUP(1,I),PUP(2,I),PUP(3,I),PUP(4,I),PUP(5,I),
283
430
     #                  VTIMUP(I),SPINUP(I)
284
431
      enddo
285
 
      if(buff(1:1).eq.'#') then
 
432
      if(buff(1:1).eq.'#' .and. .not.rwgt_skip) then
286
433
        write(ifile,'(a)') buff(1:len_trim(buff))
287
434
        read(buff,200)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
288
435
     #                    fksfather_lhe,ipartner_lhe,
290
437
     #                    jwgtinfo,mexternal,iwgtnumpartn,
291
438
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
292
439
        if(jwgtinfo.ge.1.and.jwgtinfo.le.4)then
293
 
          write(ifile,'(a)')
294
 
     # '  <rwgt>'
 
440
           write(ifile,'(a)') '  <rwgt>'
295
441
          write(ifile,401)wgtref,wgtqes2(2)
296
442
          write(ifile,402)wgtxbj(1,1),wgtxbj(2,1),
297
443
     #                    wgtxbj(1,2),wgtxbj(2,2),
327
473
          enddo
328
474
          if(jwgtinfo.eq.4) write(ifile,'(1x,e14.8,1x,i4,1x,i4)')
329
475
     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
330
 
          write(ifile,'(a)')
331
 
     # '  </rwgt>'
 
476
          write(ifile,'(a)') '  </rwgt>'
332
477
         elseif(jwgtinfo.eq.5) then
333
478
           write(ifile,'(a)')'  <rwgt>'
334
479
           if (iSorH_lhe.eq.1) then ! S-event
400
545
           write(ifile,'(a)')'  </rwgt>'
401
546
 
402
547
        elseif(jwgtinfo.eq.8)then
403
 
          write(ifile,'(a)')
404
 
     # '  <rwgt>'
 
548
           write(ifile,'(a)') '  <rwgt>'
405
549
          write(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
406
550
          do i=1,numscales
407
551
            write(ifile,404)(wgtxsecmu(i,j),j=1,numscales)
411
555
            nng=2*i
412
556
            write(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
413
557
          enddo
414
 
          write(ifile,'(a)')
415
 
     # '  </rwgt>'
 
558
          write(ifile,'(a)') '  </rwgt>'
 
559
 
 
560
        elseif(jwgtinfo.eq.9)then
 
561
           write(ifile,'(a)') '  <rwgt>'
 
562
           do i=1,numscales
 
563
              do j=1,numscales
 
564
                 idwgt=1000+(i-1)*numscales+j
 
565
                 write(ifile,601) "   <wgt id='",idwgt,"'>",wgtxsecmu(i
 
566
     $                ,j)," </wgt>"
 
567
              enddo
 
568
           enddo
 
569
           do i=1,2*numPDFpairs
 
570
              idwgt=2000+i
 
571
              write(ifile,601) "   <wgt id='",idwgt,"'>",wgtxsecPDF(i)
 
572
     $             ," </wgt>"
 
573
           enddo
 
574
           write(ifile,'(a)') '  </rwgt>'
416
575
        endif
417
576
      endif
418
 
      write(ifile,'(a)')
419
 
     # '  </event>'
 
577
      write(ifile,'(a)') '  </event>'
420
578
 200  format(1a,1x,i1,4(1x,i2),2(1x,e14.8),1x,i1,2(1x,i2),5(1x,e14.8))
421
579
 401  format(2(1x,e14.8))
422
580
 402  format(8(1x,e14.8))
428
586
 442  format(1x,e16.10,2(1x,e14.8))
429
587
 503  format(1x,i2,1x,i6,4(1x,e14.8))
430
588
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,e14.8),2(1x,e10.4))
 
589
 601  format(a12,i4,a2,1x,e11.5,a7)
431
590
c
432
591
      return
433
592
      end
443
602
      integer ifile,i
444
603
      character*140 buff
445
604
      character*80 string
 
605
      character*12 dummy12
 
606
      character*2 dummy2
446
607
      character*1 ch1
447
608
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
448
609
      double precision scale1_lhe,scale2_lhe
449
 
      integer ii,j,nps,nng,iFKS
 
610
      integer ii,j,nps,nng,iFKS,idwgt
450
611
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
451
612
      include 'reweight_all.inc'
452
613
c
587
748
            read(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
588
749
          enddo
589
750
          read(ifile,'(a)')string
 
751
        elseif(jwgtinfo.eq.9)then
 
752
           read(ifile,'(a)')string
 
753
           wgtref=XWGTUP
 
754
           do i=1,numscales
 
755
              do j=1,numscales
 
756
                 read(ifile,601) dummy12,idwgt,dummy2,wgtxsecmu(i,j)
 
757
              enddo
 
758
           enddo
 
759
           do i=1,2*numPDFpairs
 
760
              read(ifile,601) dummy12,idwgt,dummy2,wgtxsecPDF(i)
 
761
           enddo
 
762
           if (numscales.eq.0 .and. numPDFpairs.ne.0) then
 
763
              wgtxsecmu(1,1)=XWGTUP
 
764
           endif
 
765
           read(ifile,'(a)')string
590
766
        endif
591
767
        read(ifile,'(a)')string
592
768
      else
604
780
 442  format(1x,e16.10,2(1x,e14.8))
605
781
 503  format(1x,i2,1x,i6,4(1x,e14.8))
606
782
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,e14.8),2(1x,e10.4))
 
783
 601  format(a12,i4,a2,1x,e11.5,a7)
607
784
c
608
785
      return
609
786
      end
620
797
      integer ifile,i
621
798
      character*140 buff
622
799
      character*80 string
 
800
      character*12 dummy12
 
801
      character*2 dummy2
623
802
      character*1 ch1
624
803
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
625
804
      double precision scale1_lhe,scale2_lhe
626
 
      integer ii,j,nps,nng,iFKS
 
805
      integer ii,j,nps,nng,iFKS,idwgt
627
806
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
628
807
      include 'reweight_all.inc'
629
808
c
630
809
      read(ifile,'(a)')string
631
 
      if(index(string,'<event>').eq.0)then
 
810
      if(index(string,'<event').eq.0)then
632
811
        if(index(string,'</LesHouchesEvents>').ne.0)then
633
812
          buff='endoffile'
634
813
          return
774
953
            read(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
775
954
          enddo
776
955
          read(ifile,'(a)')string
 
956
        elseif(jwgtinfo.eq.9)then
 
957
           read(ifile,'(a)')string
 
958
           wgtref=XWGTUP
 
959
           do i=1,numscales
 
960
              do j=1,numscales
 
961
                 read(ifile,601) dummy12,idwgt,dummy2,wgtxsecmu(i,j)
 
962
              enddo
 
963
           enddo
 
964
           do i=1,2*numPDFpairs
 
965
              read(ifile,601) dummy12,idwgt,dummy2,wgtxsecPDF(i)
 
966
           enddo
 
967
           if (numscales.eq.0 .and. numPDFpairs.ne.0) then
 
968
              wgtxsecmu(1,1)=XWGTUP
 
969
           endif
 
970
           read(ifile,'(a)')string
777
971
        endif
778
972
        read(ifile,'(a)')string
779
973
      else
791
985
 442  format(1x,e16.10,2(1x,e14.8))
792
986
 503  format(1x,i2,1x,i6,4(1x,e14.8))
793
987
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,e14.8),2(1x,e10.4))
 
988
 601  format(a12,i4,a2,1x,e11.5,a7)
794
989
c
795
990
      return
796
991
      end
852
1047
 
853
1048
      return
854
1049
      end
 
1050
 
 
1051
 
 
1052
      subroutine case_trap3(ilength,name)
 
1053
c**********************************************************    
 
1054
c change the string to lowercase if the input is not
 
1055
c**********************************************************
 
1056
      implicit none
 
1057
c
 
1058
c     ARGUMENT
 
1059
c      
 
1060
      character*(*) name
 
1061
c
 
1062
c     LOCAL
 
1063
c
 
1064
      integer i,k,ilength
 
1065
 
 
1066
      do i=1,ilength
 
1067
         k=ichar(name(i:i))
 
1068
         if(k.ge.65.and.k.le.90) then  !upper case A-Z
 
1069
            k=ichar(name(i:i))+32   
 
1070
            name(i:i)=char(k)        
 
1071
         endif
 
1072
      enddo
 
1073
 
 
1074
      return
 
1075
      end
 
1076
 
 
1077