Higgs: rewrite.f

File rewrite.f, 4.1 KB (added by trac, 7 years ago)
Line 
1      program rw
2      implicit none
3      integer i,j
4      double precision p(0:4,10),wgt,aqcd,aqed,scale
5      integer lin,lun,ievent,ic(7,4),nexternal
6      character*140 buff2,buff,infile,outfile
7      logical done
8
9
10
11      lun=13
12      lin=14
13
14      write (*,*) 'enter input file name:'
15      read(*,*) infile
16
17      write (*,*) 'output file name is out.lhe.'
18
19
20      open(unit=lun,file=infile)
21      open(unit=lin,file='out.lhe',status='unknown')
22
23      done=.false.
24
25
26      do while (.not.done)
27         call read_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
28         
29         if (done)return
30
31         nexternal=3
32         ic(6,3)=1
33         do j=0,3
34            p(j,3)=p(j,1)+p(j,2)
35         enddo
36         p(4,3)=dsqrt(p(0,3)**2-p(1,3)**2-p(2,3)**2-p(3,3)**2)
37
38         call write_event(lin,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff)
39      enddo
40
41      write (lin,*)"</LesHouchesEvents>"
42      write (*,*) ""
43
44      close(lun)
45      close(lin)
46
47
48      return
49      end
50
51      subroutine read_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
52c********************************************************************
53c     Reads one event from data file #lun
54c     ic(*,1) = Particle ID
55c     ic(*,2) = Mothup(1)
56c     ic(*,3) = Mothup(2)
57c     ic(*,4) = ICOLUP(1)
58c     ic(*,5) = ICOLUP(2)
59c     ic(*,6) = ISTUP   -1=initial state +1=final  +2=decayed
60c     ic(*,7) = Helicity
61c********************************************************************
62      implicit none
63
64      double precision pi
65      parameter (pi = 3.1415926d0)
66c
67c     Arguments
68c     
69      integer lun
70      integer nexternal, ic(7,10)
71      logical done
72      double precision P(0:4,10),wgt,aqcd,aqed,scale
73      integer ievent
74      character*140 buff2
75c
76c     Local
77c
78      integer i,j,k
79      character*(256) buff
80      double precision xdum1,xdum2
81c
82c     Global
83c
84      logical banner_open
85      integer lun_ban
86      common/to_banner/banner_open, lun_ban
87
88      data lun_ban/37/
89      data banner_open/.false./
90c-----
91c  Begin Code
92c-----     
93      done=.false.
94      if (.not. banner_open) then
95         open (unit=lun_ban, status='scratch')
96         banner_open=.true.
97      endif
98 11   read(lun,'(a256)',end=99,err=99) buff
99      do while(index(buff,"<event") .eq. 0)
100         write(14,'(a)') buff(1:80)
101         read(lun,'(a256)',end=99,err=99) buff
102      enddo
103      read(lun,*,err=11, end=11) nexternal,k,wgt,scale,aqed,aqcd
104      do i=1,nexternal-2  !****Put here -2 because we don't care about the photons.***
105         read(lun,*,err=99,end=99) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
106     $     (p(j,i),j=1,3),p(0,i),p(4,i),xdum1,xdum2
107         ic(7,i)=int(xdum2)
108      enddo
109      do while(index(buff,"</event") .eq. 0)
110         read(lun,'(a256)',end=99,err=99) buff
111         if(buff(1:1).eq.'#') buff2=buff(1:140)
112      enddo
113      return
114 99   done=.true.
115      return
116 55   format(i3,5e19.11)         
117      end
118
119      subroutine write_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff)
120c********************************************************************
121c     Writes one event from data file #lun according to LesHouches
122c     ic(1,*) = Particle ID
123c     ic(2.*) = Mothup(1)
124c     ic(3,*) = Mothup(2)
125c     ic(4,*) = ICOLUP(1)
126c     ic(5,*) = ICOLUP(2)
127c     ic(6,*) = ISTUP   -1=initial state +1=final  +2=decayed
128c     ic(7,*) = Helicity
129c********************************************************************
130      implicit none
131c
132c     parameters
133c
134      double precision pi
135      parameter (pi = 3.1415926d0)
136c
137c     Arguments
138c     
139      integer lun, ievent
140      integer nexternal, ic(7,*)
141      double precision P(0:4,*),wgt
142      double precision aqcd, aqed, scale
143      character*140 buff
144c
145c     Local
146c
147      integer i,j,k
148c
149c     Global
150c
151
152c-----
153c  Begin Code
154c-----     
155      write(lun,'(a)') '<event>'
156      write(lun,'(i2,i4,4e15.7)') nexternal,100,wgt,scale,aqed,aqcd
157      do i=1,nexternal
158         write(lun,51) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
159     $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
160      enddo
161      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))
162      write(lun,'(a)') '</event>'
163      return
164 51   format(i9,5i5,5e19.11,f3.0,f4.0)
165      end
166