← 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_ttx_v2.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_ttx_v2.f
 
1
c
 
2
c Example analysis for "p p > t t~ [QCD]" process.
 
3
c
 
4
C----------------------------------------------------------------------
 
5
      SUBROUTINE RCLOS()
 
6
C     DUMMY IF HBOOK IS USED
 
7
C----------------------------------------------------------------------
 
8
      END
 
9
 
 
10
 
 
11
C----------------------------------------------------------------------
 
12
      SUBROUTINE HWABEG
 
13
C     USER'S ROUTINE FOR INITIALIZATION
 
14
C----------------------------------------------------------------------
 
15
      INCLUDE 'HERWIG65.INC'
 
16
      include 'reweight0.inc'
 
17
      REAL*8 pi
 
18
      integer j,kk,l,i
 
19
      PARAMETER (PI=3.14159265358979312D0)
 
20
      character*5 cc(2)
 
21
      data cc/'     ','     '/
 
22
      integer nwgt,max_weight,nwgt_analysis
 
23
      common/cnwgt/nwgt
 
24
      common/c_analysis/nwgt_analysis
 
25
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
26
      character*15 weights_info(max_weight)
 
27
      common/cwgtsinfo/weights_info
 
28
c
 
29
      call inihist
 
30
      nwgt_analysis=nwgt
 
31
      do i=1,1
 
32
      do kk=1,nwgt_analysis
 
33
        l=(kk-1)*16+(i-1)*8
 
34
        call mbook(l+ 1,'total rate    '//weights_info(kk)//cc(i),
 
35
     &       1.0d0,0.5d0,5.5d0)
 
36
        call mbook(l+ 2,'t rap         '//weights_info(kk)//cc(i),
 
37
     &       0.2d0,-5d0,5d0)
 
38
        call mbook(l+ 3,'tx rap        '//weights_info(kk)//cc(i),
 
39
     &       0.2d0,-5d0,5d0)
 
40
        call mbook(l+ 4,'t-tx pair rap '//weights_info(kk)//cc(i),
 
41
     &       0.1d0,-3d0,3d0)
 
42
        call mbook(l+ 5,'m t-tx        '//weights_info(kk)//cc(i),
 
43
     &       10d0,0d0,1000d0)
 
44
        call mbook(l+ 6,'pt t          '//weights_info(kk)//cc(i),
 
45
     &       4d0,0d0,400d0)
 
46
        call mbook(l+ 7,'pt tx         '//weights_info(kk)//cc(i),
 
47
     &       4d0,0d0,400d0)
 
48
        call mbook(l+ 8,'pt t-tx       '//weights_info(kk)//cc(i),
 
49
     &       2d0,0d0,200d0)
 
50
      enddo
 
51
      enddo
 
52
      END
 
53
 
 
54
 
 
55
C----------------------------------------------------------------------
 
56
      SUBROUTINE HWAEND
 
57
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
58
C----------------------------------------------------------------------
 
59
      INCLUDE 'HERWIG65.INC'
 
60
      REAL*8 XNORM
 
61
      INTEGER I,J,KK,l,nwgt_analysis
 
62
      integer NPL
 
63
      parameter(NPL=15000)
 
64
      common/c_analysis/nwgt_analysis
 
65
      OPEN(UNIT=99,FILE='HERQQ.TOP',STATUS='UNKNOWN')
 
66
C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG 
 
67
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
 
68
      XNORM=1.D3/DFLOAT(NEVHEP)
 
69
      DO I=1,NPL
 
70
        CALL MFINAL3(I)             
 
71
        CALL MCOPY(I,I+NPL)
 
72
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
73
        CALL MFINAL3(I+NPL)             
 
74
      ENDDO                          
 
75
C
 
76
      do i=1,1
 
77
      do kk=1,nwgt_analysis
 
78
         l=(kk-1)*16+(i-1)*8
 
79
        call multitop(NPL+l+ 1,NPL-1,3,2,'total rate   ',' ','LIN')
 
80
        call multitop(NPL+l+ 2,NPL-1,3,2,'t rap        ',' ','LOG')
 
81
        call multitop(NPL+l+ 3,NPL-1,3,2,'tx rap       ',' ','LOG')
 
82
        call multitop(NPL+l+ 4,NPL-1,3,2,'t-tx pair rap',' ','LOG')
 
83
        call multitop(NPL+l+ 5,NPL-1,3,2,'m t-tx       ',' ','LOG')
 
84
        call multitop(NPL+l+ 6,NPL-1,3,2,'pt t         ',' ','LOG')
 
85
        call multitop(NPL+l+ 7,NPL-1,3,2,'pt tx        ',' ','LOG')
 
86
        call multitop(NPL+l+ 8,NPL-1,3,2,'pt t-tx      ',' ','LOG')
 
87
      enddo
 
88
      enddo
 
89
      CLOSE(99)
 
90
      END
 
91
 
 
92
 
 
93
C----------------------------------------------------------------------
 
94
      SUBROUTINE HWANAL
 
95
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
 
96
C----------------------------------------------------------------------
 
97
      INCLUDE 'HERWIG65.INC'
 
98
      include 'reweight0.inc'
 
99
      DOUBLE PRECISION HWVDOT,PSUM(4)
 
100
      INTEGER ICHSUM,ICHINI,IHEP
 
101
      LOGICAL DIDSOF,flcuts,siq1flag,siq2flag,ddflag
 
102
      INTEGER ID,ID1,IST,IQ1,IQ2,IT1,IT2,ILP,INU,IBQ,ILM,INB,IBB,IJ
 
103
      DOUBLE PRECISION YCUT,PTCUT,ptlp,ylp,getrapidity,ptnu,ynu,
 
104
     # ptbq,ybq,ptlm,ylm,ptnb,ynb,ptbb,ybb,ptbqbb,dphibqbb,
 
105
     # getdelphi,xmbqbb,getinvm,ptlplm,dphilplm,xmlplm,ptbqlm,
 
106
     # dphibqlm,xmbqlm,ptbblp,dphibblp,xmbblp,ptbqnb,dphibqnb,
 
107
     # xmbqnb,ptbbnu,dphibbnu,xmbbnu,ptq1,ptq2,ptg,yq1,yq2,
 
108
     # etaq1,getpseudorap,etaq2,azi,azinorm,qqm,dr,yqq
 
109
      DOUBLE PRECISION XPTQ(5),XPTB(5),XPLP(5),XPNU(5),XPBQ(5),XPLM(5),
 
110
     # XPNB(5),XPBB(5),p_t(0:3),p_tx(0:3),pttx(0:3),
 
111
     # mtt,pt_t,pt_tx,pt_ttx,yt,ytx,yttx,var 
 
112
      DOUBLE PRECISION YPBQBB(4),YPLPLM(4),YPBQLM(4),YPBBLP(4),
 
113
     # YPBQNB(4),YPBBNU(4),YPTQTB(4)
 
114
      REAL*8 PI
 
115
      PARAMETER (PI=3.14159265358979312D0)
 
116
      REAL*8 WWW0
 
117
      INTEGER KK,IVLEP1,IVLEP2,i,l
 
118
      COMMON/VVLIN/IVLEP1,IVLEP2
 
119
      integer nwgt_analysis,max_weight
 
120
      common/c_analysis/nwgt_analysis
 
121
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
122
      double precision ww(max_weight),www(max_weight)
 
123
      common/cww/ww
 
124
c
 
125
      IF (IERROR.NE.0) RETURN
 
126
      IF (WW(1).EQ.0D0) THEN
 
127
         WRITE(*,*)'WW(1) = 0. Stopping'
 
128
         STOP
 
129
      ENDIF
 
130
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
 
131
C EFFECT, SO THROW THE EVENT AWAY
 
132
 
 
133
      IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN
 
134
        CALL HWWARN('HWANAL',111)
 
135
        GOTO 999
 
136
      ENDIF
 
137
      DO I=1,nwgt_analysis
 
138
         WWW(I)=EVWGT*ww(i)/ww(1)
 
139
      ENDDO
 
140
      CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PSUM)
 
141
      CALL HWVSCA(4,-1D0,PSUM,PSUM)
 
142
      ICHSUM=0
 
143
      ICHINI=ICHRG(IDHW(1))+ICHRG(IDHW(2))
 
144
      DIDSOF=.FALSE.
 
145
      IQ1=0
 
146
      IQ2=0
 
147
      DO 100 IHEP=1,NHEP
 
148
C UNCOMMENT THE FOLLOWING WHEN REMOVING THE CHECK ON MOMENTUM 
 
149
C        IF(IQ1*IQ2.EQ.1) GOTO 11
 
150
        IF (IDHW(IHEP).EQ.16) DIDSOF=.TRUE.
 
151
        IF (ISTHEP(IHEP).EQ.1) THEN
 
152
          CALL HWVSUM(4,PHEP(1,IHEP),PSUM,PSUM)
 
153
          ICHSUM=ICHSUM+ICHRG(IDHW(IHEP))
 
154
        ENDIF
 
155
        IST=ISTHEP(IHEP)      
 
156
        ID=IDHW(IHEP)
 
157
        ID1=IDHEP(IHEP)
 
158
        IF(IST.EQ.155.AND.ID1.EQ.6)THEN
 
159
C FOUND A TOP; KEEP ONLY THE FIRST ON RECORD
 
160
          IQ1=IQ1+1
 
161
          IF(IQ1.EQ.1)IT1=IHEP
 
162
        ELSEIF(IST.EQ.155.AND.ID1.EQ.-6)THEN
 
163
C FOUND AN ANTITOP; KEEP ONLY THE FIRST ON RECORD
 
164
          IQ2=IQ2+1
 
165
          IF(IQ2.EQ.1)IT2=IHEP
 
166
        ENDIF
 
167
  100 CONTINUE
 
168
      IF(IQ1*IQ2.EQ.0.AND.IERROR.EQ.0)CALL HWWARN('HWANAL',501)
 
169
C CHECK MOMENTUM AND CHARGE CONSERVATION
 
170
      IF (HWVDOT(3,PSUM,PSUM).GT.1.E-4*PHEP(4,1)**2) THEN
 
171
         CALL HWUEPR
 
172
         CALL HWWARN('HWANAL',112)
 
173
         GOTO 999
 
174
      ENDIF
 
175
      IF (ICHSUM.NE.ICHINI) THEN
 
176
         CALL HWUEPR
 
177
         CALL HWWARN('HWANAL',113)
 
178
         GOTO 999
 
179
      ENDIF
 
180
C FILL THE FOUR-MOMENTA
 
181
      DO IJ=1,5
 
182
         XPTQ(IJ)=PHEP(IJ,IT1)
 
183
         XPTB(IJ)=PHEP(IJ,IT2)
 
184
      ENDDO
 
185
      DO IJ=1,4
 
186
         p_t(MOD(IJ,4))=XPTQ(IJ)
 
187
         p_tx(MOD(IJ,4))=XPTB(IJ)
 
188
         pttx(MOD(IJ,4))=XPTQ(IJ)+XPTB(IJ)
 
189
      ENDDO
 
190
      mtt    = getinvm(pttx(0),pttx(1),pttx(2),pttx(3))
 
191
      pt_t   = dsqrt(p_t(1)**2 + p_t(2)**2)
 
192
      pt_tx  = dsqrt(p_tx(1)**2 + p_tx(2)**2)
 
193
      pt_ttx = dsqrt(pttx(1)**2 + pttx(2)**2)
 
194
      yt  = getrapidity(p_t(0), p_t(3))
 
195
      ytx = getrapidity(p_tx(0), p_tx(3))
 
196
      yttx= getrapidity(pttx(0), pttx(3))
 
197
      var=1.d0
 
198
      do i=1,1
 
199
         do kk=1,nwgt_analysis
 
200
            l=(kk-1)*16+(i-1)*8
 
201
            call mfill(l+1,var,www(kk))
 
202
            call mfill(l+2,yt,www(kk))
 
203
            call mfill(l+3,ytx,www(kk))
 
204
            call mfill(l+4,yttx,www(kk))
 
205
            call mfill(l+5,mtt,www(kk))
 
206
            call mfill(l+6,pt_t,www(kk))
 
207
            call mfill(l+7,pt_tx,www(kk))
 
208
            call mfill(l+8,pt_ttx,www(kk))
 
209
         enddo
 
210
      enddo
 
211
c
 
212
 999  return
 
213
      end
 
214
 
 
215
 
 
216
      function getrapidity(en,pl)
 
217
      implicit none
 
218
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
219
      parameter (tiny=1.d-8)
 
220
c
 
221
      xplus=en+pl
 
222
      xminus=en-pl
 
223
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
224
        if( (xplus/xminus).gt.tiny )then
 
225
          y=0.5d0*log( xplus/xminus )
 
226
        else
 
227
          y=sign(1.d0,pl)*1.d8
 
228
        endif
 
229
      else
 
230
        y=sign(1.d0,pl)*1.d8
 
231
      endif
 
232
      getrapidity=y
 
233
      return
 
234
      end
 
235
 
 
236
 
 
237
      function getpseudorap(en,ptx,pty,pl)
 
238
      implicit none
 
239
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
240
      parameter (tiny=1.d-5)
 
241
c
 
242
      pt=sqrt(ptx**2+pty**2)
 
243
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
244
        eta=sign(1.d0,pl)*1.d8
 
245
      else
 
246
        th=atan2(pt,pl)
 
247
        eta=-log(tan(th/2.d0))
 
248
      endif
 
249
      getpseudorap=eta
 
250
      return
 
251
      end
 
252
 
 
253
 
 
254
      function getinvm(en,ptx,pty,pl)
 
255
      implicit none
 
256
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
257
      parameter (tiny=1.d-5)
 
258
c
 
259
      tmp=en**2-ptx**2-pty**2-pl**2
 
260
      if(tmp.gt.0.d0)then
 
261
        tmp=sqrt(tmp)
 
262
      elseif(tmp.gt.-tiny)then
 
263
        tmp=0.d0
 
264
      else
 
265
        write(*,*)'Attempt to compute a negative mass'
 
266
        stop
 
267
      endif
 
268
      getinvm=tmp
 
269
      return
 
270
      end
 
271
 
 
272
 
 
273
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
274
      implicit none
 
275
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
276
      parameter (tiny=1.d-5)
 
277
c
 
278
      pt1=sqrt(ptx1**2+pty1**2)
 
279
      pt2=sqrt(ptx2**2+pty2**2)
 
280
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
281
        tmp=ptx1*ptx2+pty1*pty2
 
282
        tmp=tmp/(pt1*pt2)
 
283
        if(abs(tmp).gt.1.d0+tiny)then
 
284
          write(*,*)'Cosine larger than 1'
 
285
          stop
 
286
        elseif(abs(tmp).ge.1.d0)then
 
287
          tmp=sign(1.d0,tmp)
 
288
        endif
 
289
        tmp=acos(tmp)
 
290
      else
 
291
        tmp=1.d8
 
292
      endif
 
293
      getdelphi=tmp
 
294
      return
 
295
      end