← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin_standalone.f

  • Committer: Marco Zaro
  • Date: 2014-01-27 16:54:10 UTC
  • mfrom: (78.124.55 MG5_aMC_2.1)
  • Revision ID: marco.zaro@gmail.com-20140127165410-5lma8c2hzbzm426j
merged with lp:~maddevelopers/madgraph5/MG5_aMC_2.1 r 267

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin_standalone.f
1
 
C----------------------------------------------------------------------
2
 
      SUBROUTINE UPEVNT
3
 
C----------------------------------------------------------------------
4
 
C  Reads MC@NLO input files and fills Les Houches event common HEPEUP
5
 
C  Event file is written by MadFKS
6
 
C----------------------------------------------------------------------
7
 
      INCLUDE 'HERWIG65.INC'
8
 
C---Les Houches Event Common Block
9
 
      INTEGER MAXNUP
10
 
      PARAMETER (MAXNUP=500)
11
 
      INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
12
 
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
13
 
      COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
14
 
     &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
15
 
     &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
16
 
     &              SPINUP(MAXNUP)
17
 
      INTEGER ISORH_LHE,IFKS_LHE,JFKS_LHE,FKSFATHER_LHE,IPARTNER_LHE
18
 
      DOUBLE PRECISION SCALE1_LHE,SCALE2_LHE
19
 
      INTEGER MQQ
20
 
      COMMON/cMQQ/MQQ
21
 
      INTEGER IUNIT
22
 
      PARAMETER (IUNIT=61)
23
 
      CHARACTER*80 STRING
24
 
      CHARACTER*140 BUFF
25
 
      CHARACTER*9 CH1
26
 
      INTEGER I
27
 
c
28
 
      DOUBLE PRECISION PCM(5),PTR,XSCALE
29
 
      INTEGER ISCALE
30
 
      LOGICAL REMIT
31
 
C
32
 
      IF (IERROR.NE.0) RETURN
33
 
c
34
 
      ISORH_LHE=0
35
 
      read(iunit,'(a)')string
36
 
      if(INDEX(STRING,'<event>').eq.0)then
37
 
        CALL HWWARN('UPEVNT',500)
38
 
      endif
39
 
      read(iunit,503)NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
40
 
C---Les Houches expects mean weight to be the cross section in pb
41
 
      XWGTUP=XWGTUP*MQQ
42
 
      do i=1,nup
43
 
        read(iunit,504)IDUP(I),ISTUP(I),MOTHUP(1,I),MOTHUP(2,I),
44
 
     #                 ICOLUP(1,I),ICOLUP(2,I),
45
 
     #                 PUP(1,I),PUP(2,I),PUP(3,I),PUP(4,I),PUP(5,I),
46
 
     #                 VTIMUP(I),SPINUP(I)
47
 
      enddo
48
 
      read(iunit,'(a)')buff
49
 
      if(buff(1:1).eq.'#')then
50
 
        read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
51
 
     #                    fksfather_lhe,ipartner_lhe,
52
 
     #                    scale1_lhe,scale2_lhe
53
 
        read(iunit,'(a)')string
54
 
      else
55
 
        string=buff(1:len_trim(buff))
56
 
      endif
57
 
      if(INDEX(STRING,'</event>').eq.0)then
58
 
        CALL HWWARN('UPEVNT',501)
59
 
      endif
60
 
c Do as done for ttbar in standard MC@NLO
61
 
      REMIT=ISORH_LHE.EQ.2
62
 
      CALL HWVSUM(4,PUP(1,1),PUP(1,2),PCM)
63
 
      CALL HWUMAS(PCM)
64
 
      XSCALE=SCALUP
65
 
      SCALUP=PCM(5)
66
 
c if shape is used, uncomment the following
67
 
c$$$      IF (XSCALE.GT.0D0.AND.XSCALE.LT.PCM(5)) SCALUP=XSCALE
68
 
      ISCALE=0
69
 
      IF (REMIT) THEN
70
 
c Extra parton is #5 here, not #3 as in standalone MC@NLO
71
 
         IF (ISCALE.EQ.0) THEN
72
 
            PTR=SQRT(PUP(1,5)**2+PUP(2,5)**2)
73
 
            SCALUP=PCM(5)-2.*PTR
74
 
         ELSE
75
 
            CALL HWWARN('UPEVNT',501)
76
 
         ENDIF
77
 
      ENDIF
78
 
c Modify what follows to set scale of H or S events in a different way
79
 
c$$$      IF(ISORH_LHE.EQ.2)THEN
80
 
c$$$c H events
81
 
c$$$        IF(SCALE2_LHE.GT.0.D0)SCALUP=SCALE2_LHE
82
 
c$$$      ENDIF
83
 
 503  format(1x,i2,1x,i6,4(1x,d14.8))
84
 
 504  format(1x,i8,1x,i2,4(1x,i4),5(1x,d14.8),2(1x,d10.4))
85
 
      END
86
 
 
87
 
 
88
 
C----------------------------------------------------------------------
89
 
      SUBROUTINE UPINIT
90
 
C----------------------------------------------------------------------
91
 
C  Reads MC@NLO input headers and fills Les Houches run common HEPRUP
92
 
C  Event file is written by MadFKS
93
 
C----------------------------------------------------------------------
94
 
      INCLUDE 'HERWIG65.INC'
95
 
C--Les Houches Common Blocks
96
 
      INTEGER MAXPUP
97
 
      PARAMETER(MAXPUP=100)
98
 
      INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
99
 
      DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
100
 
      COMMON /HEPRUP/ IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
101
 
     &                IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),
102
 
     &                XMAXUP(MAXPUP),LPRUP(MAXPUP)
103
 
      INTEGER MAXNUP
104
 
      PARAMETER (MAXNUP=500)
105
 
      INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
106
 
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
107
 
      COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
108
 
     &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
109
 
     &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
110
 
     &              SPINUP(MAXNUP)
111
 
c Hard event file (to be entered in Herwig driver)
112
 
      CHARACTER*50 QQIN
113
 
      COMMON/VVJIN/QQIN
114
 
      LOGICAL HEADERS
115
 
      CHARACTER*80 STRING
116
 
      INTEGER MQQ
117
 
      COMMON/cMQQ/MQQ
118
 
C
119
 
      IF (IERROR.NE.0) RETURN
120
 
C--SET UP INPUT FILES
121
 
      OPEN(UNIT=61,FILE=QQIN,STATUS='UNKNOWN')
122
 
C--Read (non compulsory) headers here if need be
123
 
      HEADERS=.TRUE.
124
 
      MQQ=-1
125
 
      DO WHILE(HEADERS)
126
 
         READ(61,'(a)')STRING
127
 
         if(index(string,'<header>').ne.0) then
128
 
            READ(61,*) MQQ
129
 
         endif
130
 
        HEADERS=INDEX(STRING,'<init>').eq.0
131
 
      ENDDO
132
 
      if(MQQ.eq.-1) call HWWARN('UPINIT',501)
133
 
C--Read up to </init> in the event file
134
 
      read(61,501)IDBMUP(1),IDBMUP(2),EBMUP(1),EBMUP(2),
135
 
     #            PDFGUP(1),PDFGUP(2),PDFSUP(1),PDFSUP(2),
136
 
     #            IDWTUP,NPRUP
137
 
      read(61,502)XSECUP(1),XERRUP(1),XMAXUP(1),LPRUP(1)
138
 
      read(61,'(a)')string
139
 
      if(INDEX(STRING,'</init>').eq.0)then
140
 
        CALL HWWARN('UPINIT',500)
141
 
      endif
142
 
c
143
 
 501  format(2(1x,i6),2(1x,d14.8),2(1x,i2),2(1x,i6),1x,i2,1x,i3)
144
 
 502  format(3(1x,d14.8),1x,i6)
145
 
      return
146
 
 999  END
147
 
 
148
 
 
149
 
C----------------------------------------------------------------------
150
 
      SUBROUTINE HWURSC(NP,PP)
151
 
C  RESCALES A SET OF NP (<21) 3-MOMENTA PP(1-3,*) IN
152
 
C  THEIR CMF TO PUT PP ON MASS-SHELL AT MASSES PP(5,*) 
153
 
C----------------------------------------------------------------------
154
 
      IMPLICIT NONE
155
 
      INTEGER NP,IP,IT,NT
156
 
      DOUBLE PRECISION PP(5,*),P(5,20),P2(20),M2(20),SP(5),
157
 
     & TINY,FAC,ECM,DCM,EP,STEP,FRT,HWUSQR
158
 
      DATA TINY,NT/1D-9,20/
159
 
      IF (NP.GT.20) CALL HWWARN('HWURSC',300+NP)
160
 
C--COMPUTE CM MOMENTUM
161
 
      CALL HWVZRO(4,SP)
162
 
      DO IP=1,NP
163
 
         CALL HWVSUM(4,PP(1,IP),SP,SP)
164
 
      ENDDO
165
 
      CALL HWUMAS(SP)
166
 
C--BOOST TO CMF
167
 
      DO IP=1,NP
168
 
         CALL HWULOF(SP,PP(1,IP),P(1,IP))
169
 
         P2(IP)=P(1,IP)**2+P(2,IP)**2+P(3,IP)**2
170
 
         M2(IP)=P(5,IP)**2
171
 
      ENDDO
172
 
C--ITERATE RESCALING OF 3-MOMENTA
173
 
      FAC=1D0
174
 
      DO IT=1,NT
175
 
         ECM=0D0
176
 
         DCM=0D0
177
 
         DO IP=1,NP
178
 
            EP=HWUSQR(M2(IP)+FAC*P2(IP))
179
 
            IF (EP.GT.0D0) THEN
180
 
               ECM=ECM+EP
181
 
               DCM=DCM+P2(IP)/EP
182
 
            ENDIF
183
 
         ENDDO
184
 
         IF (DCM.EQ.0D0) CALL HWWARN('HWURSC',390)
185
 
         STEP=2D0*(ECM-SP(5))/DCM
186
 
         FAC=FAC-STEP
187
 
         IF (ABS(STEP).LT.TINY) GOTO 100
188
 
      ENDDO
189
 
C--FAILED TO CONVERGE
190
 
      CALL HWWARN('HWURSC',1)
191
 
C--CONVERGED: RESCALE 3-MOMENTA AND BOOST BACK 
192
 
 100  IF (FAC.LT.0D0) CALL HWWARN('HWURSC',391)
193
 
      FRT=SQRT(FAC)
194
 
      DO IP=1,NP
195
 
         CALL HWVSCA(3,FRT,P(1,IP),P(1,IP))
196
 
         P(4,IP)=SQRT(M2(IP)+FAC*P2(IP))
197
 
         CALL HWULOB(SP,P(1,IP),PP(1,IP))
198
 
      ENDDO
199
 
      END