← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/Utilities/plot_events_lhef.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/Utilities/plot_events_lhef.f
 
1
c Same as plot_events.f, but relevant to LH showered files (ntuples).
 
2
c Compile with
 
3
c   f77 -g -ffixed-line-length-132 -fno-automatic 
 
4
c      -o plot_events_lhef plot_events_lhef.f madfks_plot*.f madfks_dbook.f 
 
5
c         handling_lhe_events.f any-dependencies-in-madfksplot
 
6
c Copy in this directory
 
7
c  reweight0.inc dbook.inc
 
8
      program plot_events
 
9
      implicit none
 
10
      integer maxevt,ifile,i,npart,itype
 
11
      integer IDBMUP(2),PDFGUP(2),PDFSUP(2),IDWTUP,NPRUP,LPRUP
 
12
      double precision EBMUP(2),XSECUP,XERRUP,XMAXUP
 
13
      INTEGER MAXNUP
 
14
      PARAMETER (MAXNUP=500)
 
15
      INTEGER NUP,IDPRUP,IDUP(MAXNUP),ISTUP(MAXNUP),
 
16
     # MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP)
 
17
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
18
     # PUP(5,MAXNUP),VTIMUP(MAXNUP),SPINUP(MAXNUP)
 
19
      double precision sum_wgt,sum_wgt_resc,rat,www
 
20
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
 
21
      double precision scale1_lhe,scale2_lhe
 
22
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
23
      double precision zero
 
24
      parameter (zero=0.d0)
 
25
      integer ifk88seed
 
26
      real*8 rndec(10),fk88random
 
27
      common/crndec/rndec
 
28
      character*80 event_file
 
29
      character*140 buff
 
30
      character*9 ch1
 
31
      logical AddInfoLHE
 
32
      logical keepevent,rescale,doscale
 
33
      integer kr,kf,npdfset
 
34
      logical usexinteg,mint
 
35
      common/cusexinteg/usexinteg,mint
 
36
      integer itmax,ncall
 
37
      common/citmax/itmax,ncall
 
38
      logical unwgt
 
39
      double precision evtsgn
 
40
      common /c_unwgt/evtsgn,unwgt
 
41
      include 'reweight0.inc'
 
42
 
 
43
      integer nexternal
 
44
      parameter (nexternal=20)
 
45
      integer j,k
 
46
      real*8 ecm,xmass(3*nexternal),xmom(0:3,3*nexternal)
 
47
      character*10 MonteCarlo
 
48
      integer isc,ipdf
 
49
 
 
50
      usexinteg=.false.
 
51
      mint=.false.
 
52
      keepevent=.true.
 
53
      rescale=.true.
 
54
c evtsgn and unwgt used by plotting routines. Sort of fake values here
 
55
      evtsgn=0.d0
 
56
      unwgt=.false.
 
57
      itmax=1
 
58
      ncall=1
 
59
 
 
60
      ifk88seed=1
 
61
      do i=1,10
 
62
        rndec(i)=fk88random(ifk88seed)
 
63
      enddo
 
64
 
 
65
      write(*,*)'Enter event file name'
 
66
      read(*,*)event_file
 
67
 
 
68
      write(*,*)'Enter (kr,kf) -- central result is (1,1)'
 
69
      write(*,*)' Enter (0,0) to skip'
 
70
      read(*,*)kr,kf
 
71
      if(kr.eq.0.or.kf.eq.0)then
 
72
        write(*,*)'Enter number of PDF error set (e.g. 1,..,40)'
 
73
        write(*,*)' Enter 0 to skip'
 
74
        read(*,*)npdfset
 
75
      else
 
76
        npdfset=0
 
77
      endif
 
78
      if(kr.eq.0.and.kf.eq.0.and.npdfset.eq.0)then
 
79
        write(*,*)'Will not rescale weights'
 
80
        rescale=.false.
 
81
      else
 
82
        write(*,*)'Will rescale weights'
 
83
        if(npdfset.eq.0)then
 
84
          doscale=.true.
 
85
        else
 
86
          doscale=.false.
 
87
        endif
 
88
      endif
 
89
 
 
90
 
 
91
      ifile=34
 
92
      open (unit=ifile,file=event_file,status='old')
 
93
      AddInfoLHE=.false.
 
94
 
 
95
      call read_lhef_header_full(ifile,maxevt,isc,ipdf,MonteCarlo)
 
96
      numscales=int(sqrt(dble(isc)))
 
97
      numPDFpairs=ipdf/2
 
98
c Showered LH files have maxevt<0; in that case, it is not the number of
 
99
c events, but its upper bound
 
100
      if(maxevt.gt.0)then
 
101
        write(*,*)'This is unlikely to be a showered LH file'
 
102
        stop
 
103
      endif
 
104
      maxevt=abs(maxevt)
 
105
      call read_lhef_init(ifile,
 
106
     &     IDBMUP,EBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,
 
107
     &     XSECUP,XERRUP,XMAXUP,LPRUP)
 
108
 
 
109
      itype=12
 
110
      sum_wgt=0d0
 
111
      sum_wgt_resc=0d0
 
112
      call initplot
 
113
      i=0
 
114
      dowhile(i.lt.maxevt.and.keepevent)
 
115
         call read_lhef_event_catch(ifile,
 
116
     &        NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
117
     &        IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
118
 
 
119
         if(index(buff,'endoffile').ne.0)then
 
120
           keepevent=.false.
 
121
           goto 111
 
122
         endif
 
123
 
 
124
         i=i+1
 
125
         sum_wgt=sum_wgt+XWGTUP
 
126
         if(rescale)then
 
127
           if(doscale)then
 
128
             rat=wgtxsecmu(kr,kf)/wgtref
 
129
           else
 
130
             rat=wgtxsecPDF(npdfset)/wgtref
 
131
           endif
 
132
           www=XWGTUP*rat
 
133
         else
 
134
           www=XWGTUP
 
135
         endif
 
136
 
 
137
         sum_wgt_resc=sum_wgt_resc+www
 
138
 
 
139
         if(i.eq.1.and.buff(1:1).eq.'#')AddInfoLHE=.true.
 
140
         if(AddInfoLHE)then
 
141
           if(buff(1:1).ne.'#')then
 
142
             write(*,*)'Inconsistency in event file',i,' ',buff
 
143
             stop
 
144
           endif
 
145
           read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
 
146
     #                       fksfather_lhe,ipartner_lhe,
 
147
     #                       scale1_lhe,scale2_lhe,
 
148
     #                       jwgtinfo,mexternal,iwgtnumpartn,
 
149
     #            wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
150
         endif
 
151
 
 
152
         npart=0
 
153
         do k=1,nup
 
154
           if(abs(ISTUP(k)).eq.1)then
 
155
             npart=npart+1
 
156
             xmass(npart)=pup(5,k)
 
157
             do j=1,4
 
158
               xmom(mod(j,4),npart)=pup(j,k)
 
159
             enddo
 
160
           endif
 
161
         enddo
 
162
 
 
163
c Showered LH files only contain final-state particles.
 
164
c Don't check momentum conservation
 
165
 
 
166
         call outfun(xmom,zero,www,itype)
 
167
 
 
168
 111     continue
 
169
 
 
170
      enddo
 
171
      close(34)
 
172
 
 
173
      call mclear
 
174
      open(unit=99,file='events.top',status='unknown')
 
175
      call topout
 
176
      close(99)
 
177
 
 
178
      write (*,*) 'The total number of events is:',i
 
179
      write (*,*) 'The sum of weights is:',sum_wgt
 
180
      if(rescale)
 
181
     #  write (*,*) 'The sum of rescaled weights is:',sum_wgt_resc
 
182
 
 
183
      end
 
184
 
 
185
 
 
186
c Dummy subroutine (normally used with vegas when resuming plots)
 
187
      subroutine resume()
 
188
      end
 
189
 
 
190
 
 
191
      FUNCTION FK88RANDOM(SEED)
 
192
*     -----------------
 
193
* Ref.: K. Park and K.W. Miller, Comm. of the ACM 31 (1988) p.1192
 
194
* Use seed = 1 as first value.
 
195
*
 
196
      IMPLICIT INTEGER(A-Z)
 
197
      REAL*8 MINV,FK88RANDOM
 
198
      SAVE
 
199
      PARAMETER(M=2147483647,A=16807,Q=127773,R=2836)
 
200
      PARAMETER(MINV=0.46566128752458d-09)
 
201
      HI = SEED/Q
 
202
      LO = MOD(SEED,Q)
 
203
      SEED = A*LO - R*HI
 
204
      IF(SEED.LE.0) SEED = SEED + M
 
205
      FK88RANDOM = SEED*MINV
 
206
      END
 
207
 
 
208
 
 
209
      subroutine boostwdir2(chybst,shybst,chybstmo,xd,xin,xout)
 
210
c chybstmo = chybst-1; if it can be computed analytically it improves
 
211
c the numerical accuracy
 
212
      implicit none
 
213
      real*8 chybst,shybst,chybstmo,xd(1:3),xin(0:3),xout(0:3)
 
214
      real*8 tmp,en,pz
 
215
      integer i
 
216
c
 
217
      if(abs(xd(1)**2+xd(2)**2+xd(3)**2-1).gt.1.d-6)then
 
218
        write(*,*)'Error #1 in boostwdir2',xd
 
219
        stop
 
220
      endif
 
221
c
 
222
      en=xin(0)
 
223
      pz=xin(1)*xd(1)+xin(2)*xd(2)+xin(3)*xd(3)
 
224
      xout(0)=en*chybst-pz*shybst
 
225
      do i=1,3
 
226
        xout(i)=xin(i)+xd(i)*(pz*chybstmo-en*shybst)
 
227
      enddo
 
228
c
 
229
      return
 
230
      end