← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lvl.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lvl.f
 
1
c
 
2
c Example analysis for "p p > e+ ve [QCD]" process.
 
3
c Example analysis for "p p > e- ve~ [QCD]" process.
 
4
c Example analysis for "p p > mu+ vm [QCD]" process.
 
5
c Example analysis for "p p > mu- vm~ [QCD]" process.
 
6
c Example analysis for "p p > ta+ vt [QCD]" process.
 
7
c Example analysis for "p p > ta- vt~ [QCD]" process.
 
8
c
 
9
C----------------------------------------------------------------------
 
10
      SUBROUTINE RCLOS()
 
11
C     DUMMY IF HBOOK IS USED
 
12
C----------------------------------------------------------------------
 
13
      END
 
14
 
 
15
 
 
16
C----------------------------------------------------------------------
 
17
      SUBROUTINE HWABEG
 
18
C     USER'S ROUTINE FOR INITIALIZATION
 
19
C----------------------------------------------------------------------
 
20
      INCLUDE 'HERWIG65.INC'
 
21
      include 'reweight0.inc'
 
22
      integer j,kk,l,i
 
23
      character*5 cc(2)
 
24
      data cc/'     ','     '/
 
25
      integer nwgt,max_weight,nwgt_analysis
 
26
      common/cnwgt/nwgt
 
27
      common/c_analysis/nwgt_analysis
 
28
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
29
      character*15 weights_info(max_weight)
 
30
      common/cwgtsinfo/weights_info
 
31
c
 
32
      call inihist
 
33
      nwgt_analysis=nwgt
 
34
      do i=1,1
 
35
      do kk=1,nwgt_analysis
 
36
        l=(kk-1)*16+(i-1)*8
 
37
        call mbook(l+1,'total rate '//weights_info(kk)//cc(i),
 
38
     &       1.0d0,0.5d0,5.5d0)
 
39
        call mbook(l+2,'lep rapidity '//weights_info(kk)//cc(i),
 
40
     &       0.5d0,-5d0,5d0)
 
41
        call mbook(l+3,'lep pt '//weights_info(kk)//cc(i),
 
42
     &       10d0,0d0,200d0)
 
43
        call mbook(l+4,'et miss '//weights_info(kk)//cc(i),
 
44
     &       10d0,0d0,200d0)
 
45
        call mbook(l+5,'trans. mass '//weights_info(kk)//cc(i),
 
46
     &       5d0,0d0,200d0)
 
47
        call mbook(l+6,'w rapidity '//weights_info(kk)//cc(i),
 
48
     &       0.5d0,-5d0,5d0)
 
49
        call mbook(l+7,'w pt '//weights_info(kk)//cc(i),
 
50
     &       10d00,0d0,200d0)
 
51
        call mbook(l+8,'cphi[l,vl] '//weights_info(kk)//cc(i),
 
52
     &       0.05d0,-1d0,1d0)
 
53
      enddo
 
54
      enddo
 
55
 999  END
 
56
 
 
57
 
 
58
C----------------------------------------------------------------------
 
59
      SUBROUTINE HWAEND
 
60
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
61
C----------------------------------------------------------------------
 
62
      INCLUDE 'HERWIG65.INC'
 
63
      REAL*8 XNORM
 
64
      INTEGER I,J,KK,l,nwgt_analysis
 
65
      integer NPL
 
66
      parameter(NPL=15000)
 
67
      common/c_analysis/nwgt_analysis
 
68
      OPEN(UNIT=99,FILE='HERLL.TOP',STATUS='UNKNOWN')
 
69
C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG 
 
70
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
 
71
      XNORM=1.D3/DFLOAT(NEVHEP)
 
72
      DO I=1,NPL              
 
73
        CALL MFINAL3(I)             
 
74
        CALL MCOPY(I,I+NPL)
 
75
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
76
        CALL MFINAL3(I+NPL)             
 
77
      ENDDO                          
 
78
C
 
79
      do i=1,1
 
80
      do kk=1,nwgt_analysis
 
81
         l=(kk-1)*16+(i-1)*8
 
82
         call multitop(NPL+l+1,NPL-1,3,2,'total rate   ',' ','LIN')
 
83
         call multitop(NPL+l+2,NPL-1,3,2,'lep rapidity ',' ','LIN')
 
84
         call multitop(NPL+l+3,NPL-1,3,2,'lep pt       ',' ','LOG')
 
85
         call multitop(NPL+l+4,NPL-1,3,2,'et miss      ',' ','LOG')
 
86
         call multitop(NPL+l+5,NPL-1,3,2,'trans. mass  ',' ','LOG')
 
87
         call multitop(NPL+l+6,NPL-1,3,2,'w rapidity   ',' ','LIN')
 
88
         call multitop(NPL+l+7,NPL-1,3,2,'w pt         ',' ','LOG')
 
89
         call multitop(NPL+l+8,NPL-1,3,2,'cphi[l,vl]   ',' ','LOG')
 
90
      enddo
 
91
      enddo
 
92
      CLOSE(99)
 
93
      END
 
94
 
 
95
C----------------------------------------------------------------------
 
96
      SUBROUTINE HWANAL
 
97
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
 
98
C----------------------------------------------------------------------
 
99
      INCLUDE 'HERWIG65.INC'
 
100
      include 'reweight0.inc'
 
101
      DOUBLE PRECISION HWVDOT,PSUM(4),PPV(5),PTW,YW,YE,PPL(5),PPLB(5),
 
102
     & PTE,PLL,PTLB,PLLB,var,mtr,etmiss,cphi
 
103
      INTEGER ICHSUM,ICHINI,IHEP,IV,IFV,IST,ID,IJ,ID1,JPR,IDENT,
 
104
     #  ILL,ILLB,IHRD
 
105
      LOGICAL DIDSOF,FOUNDL,FOUNDN,ISL,ISN
 
106
      REAL*8 PI,getrapidity
 
107
      PARAMETER (PI=3.14159265358979312D0)
 
108
      REAL*8 WWW0,TINY,SIGNL,SIGNN
 
109
      INTEGER KK,I,L,IL,IN
 
110
      DATA TINY/.1D-5/
 
111
      integer nwgt_analysis,max_weight
 
112
      common/c_analysis/nwgt_analysis
 
113
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
114
      double precision ww(max_weight),www(max_weight)
 
115
      common/cww/ww
 
116
c
 
117
      IF (IERROR.NE.0) RETURN
 
118
      IF (WW(1).EQ.0D0) THEN
 
119
         WRITE(*,*)'WW(1) = 0. Stopping'
 
120
         STOP
 
121
      ENDIF
 
122
C CHOOSE IDENT = 11 FOR E - NU_E
 
123
C        IDENT = 13 FOR MU - NU_MU
 
124
C        IDENT = 15 FOR TAU - NU_TAU
 
125
      IDENT=11
 
126
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
 
127
C EFFECT, SO THROW THE EVENT AWAY
 
128
      IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN
 
129
        CALL HWWARN('HWANAL',111)
 
130
        GOTO 999
 
131
      ENDIF
 
132
      DO I=1,nwgt_analysis
 
133
         WWW(I)=EVWGT*ww(i)/ww(1)
 
134
      ENDDO
 
135
      CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PSUM)
 
136
      CALL HWVSCA(4,-1D0,PSUM,PSUM)
 
137
      ICHSUM=0
 
138
      ICHINI=ICHRG(IDHW(1))+ICHRG(IDHW(2))
 
139
      DIDSOF=.FALSE.
 
140
      FOUNDL=.FALSE.
 
141
      FOUNDN=.FALSE.
 
142
      DO 100 IHEP=1,NHEP
 
143
        IF (IDHW(IHEP).EQ.16) DIDSOF=.TRUE.
 
144
        IF (ISTHEP(IHEP).EQ.1) THEN
 
145
          CALL HWVSUM(4,PHEP(1,IHEP),PSUM,PSUM)
 
146
          ICHSUM=ICHSUM+ICHRG(IDHW(IHEP))
 
147
        ENDIF
 
148
        IST=ISTHEP(IHEP)      
 
149
        ID=IDHW(IHEP)
 
150
        ID1=IDHEP(IHEP)
 
151
        ISL=ABS(ID1).EQ.IDENT
 
152
        ISN=ABS(ID1).EQ.IDENT+1
 
153
        IF(((IST.GE.120.AND.IST.LE.125).OR.IST.EQ.1.OR.IST.EQ.198)
 
154
     &      .AND.ISL.AND..NOT.FOUNDL)THEN
 
155
            IL=IHEP
 
156
            FOUNDL=.TRUE.
 
157
            SIGNL=SIGN(1D0,DBLE(ID1))
 
158
         ENDIF
 
159
         IF(((IST.GE.120.AND.IST.LE.125).OR.IST.EQ.1)
 
160
     &       .AND.ISN.AND..NOT.FOUNDN)THEN
 
161
            IN=IHEP
 
162
            FOUNDN=.TRUE.
 
163
            SIGNN=SIGN(1D0,DBLE(ID1))
 
164
         ENDIF
 
165
  100 CONTINUE
 
166
      IF(.NOT.FOUNDL.OR..NOT.FOUNDN)THEN
 
167
         WRITE(*,*)'NO LEPTONS FOUND.'
 
168
         WRITE(*,*)'CURRENTLY THIS ANALYSIS LOOKS FOR'
 
169
         IF(IDENT.EQ.11)WRITE(*,*)'E - NU_E'
 
170
         IF(IDENT.EQ.13)WRITE(*,*)'MU - NU_MU'
 
171
         IF(IDENT.EQ.15)WRITE(*,*)'TAU - NU_TAU'
 
172
         WRITE(*,*)'IF THIS IS NOT MEANT,'
 
173
         WRITE(*,*)'PLEASE CHANGE THE VALUE OF IDENT IN THIS FILE.'
 
174
         STOP
 
175
      ENDIF
 
176
      IF(SIGNN.EQ.SIGNL)THEN
 
177
         WRITE(*,*)'TWO SAME SIGN LEPTONS!'
 
178
         WRITE(*,*)IL,IN,SIGNL,SIGNN
 
179
         STOP
 
180
      ENDIF
 
181
C CHECK MOMENTUM AND CHARGE CONSERVATION
 
182
      IF (HWVDOT(3,PSUM,PSUM).GT.1.E-4*PHEP(4,1)**2) THEN
 
183
         CALL HWUEPR
 
184
         CALL HWWARN('HWANAL',112)
 
185
         GOTO 999
 
186
      ENDIF
 
187
      IF (ICHSUM.NE.ICHINI) THEN
 
188
         CALL HWUEPR
 
189
         CALL HWWARN('HWANAL',113)
 
190
         GOTO 999
 
191
      ENDIF
 
192
      DO IJ=1,5
 
193
        PPL(IJ)=PHEP(IJ,IN)
 
194
        PPLB(IJ)=PHEP(IJ,IL)
 
195
        PPV(IJ)=PPL(IJ)+PPLB(IJ)
 
196
      ENDDO
 
197
      ye     = getrapidity(pplb(4), pplb(3))
 
198
      yw     = getrapidity(ppv(4), ppv(3))
 
199
      pte    = dsqrt(pplb(1)**2 + pplb(2)**2)
 
200
      ptw    = dsqrt(ppv(1)**2+ppv(2)**2)
 
201
      etmiss = dsqrt(ppl(1)**2 + ppl(2)**2)
 
202
      mtr    = dsqrt(2d0*pte*etmiss-2d0*ppl(1)*pplb(1)-2d0*ppl(2)*pplb(2))
 
203
      cphi   = (ppl(1)*pplb(1)+ppl(2)*pplb(2))/pte/etmiss
 
204
      var    = 1.d0
 
205
      do i=1,1
 
206
         do kk=1,nwgt_analysis
 
207
            l=(kk-1)*16+(i-1)*8
 
208
            call mfill(l+1,var,www(kk))
 
209
            call mfill(l+2,ye,www(kk))
 
210
            call mfill(l+3,pte,www(kk))
 
211
            call mfill(l+4,etmiss,www(kk))
 
212
            call mfill(l+5,mtr,www(kk))
 
213
            call mfill(l+6,yw,www(kk))
 
214
            call mfill(l+7,ptw,www(kk))
 
215
            call mfill(l+8,cphi,www(kk))
 
216
         enddo
 
217
      enddo
 
218
 999  END
 
219
      
 
220
      function getrapidity(en,pl)
 
221
      implicit none
 
222
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
223
      parameter (tiny=1.d-8)
 
224
      xplus=en+pl
 
225
      xminus=en-pl
 
226
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
227
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
228
            y=0.5d0*log( xplus/xminus  )
 
229
         else
 
230
            y=sign(1.d0,pl)*1.d8
 
231
         endif
 
232
      else 
 
233
         y=sign(1.d0,pl)*1.d8
 
234
      endif
 
235
      getrapidity=y
 
236
      return
 
237
      end