← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.3

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_v2_hepmc.f

  • Committer: olivier-mattelaer
  • Date: 2017-03-08 12:31:17 UTC
  • Revision ID: olivier-mattelaer-20170308123117-h0zkqjyh9sihsc61
empty version to have an effective freeze of the code

Show diffs side-by-side

added added

removed removed

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