← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

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

  • Committer: andrew.lifson at lu
  • Date: 2021-09-01 15:34:39 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210901153439-7fasjhav4cp4m88r
testing a new repository of a madgraph folder

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_taptam.f
 
1
c
 
2
c Example analysis for "p p > ta+ ta- [QCD]" process.
 
3
c
 
4
c It features the HwU format for histogram booking and output.
 
5
c The details of how to process/manipulate the resulting .HwU file,
 
6
c in particular how to plot it using gnuplot, I refer the reader to this
 
7
c FAQ:
 
8
c
 
9
c      https://answers.launchpad.net/mg5amcnlo/+faq/2671
 
10
c
 
11
c It mostly relies on using the following madgraph5 module in standalone
 
12
c
 
13
c  <MG5_aMC_install_dir>/madgraph/various/histograms.py
 
14
c
 
15
c You can learn about how to run it and what options are available with
 
16
c
 
17
c  python <MG5_aMC_install_dir>/madgraph/various/histograms.py --help
 
18
c
 
19
C----------------------------------------------------------------------
 
20
      SUBROUTINE RCLOS()
 
21
C     DUMMY IF HBOOK IS USED
 
22
C----------------------------------------------------------------------
 
23
      END
 
24
 
 
25
 
 
26
C----------------------------------------------------------------------
 
27
      SUBROUTINE HWABEG
 
28
C     USER''S ROUTINE FOR INITIALIZATION
 
29
C----------------------------------------------------------------------
 
30
      use HwU_wgts_info_len
 
31
      INCLUDE 'HERWIG65.INC'
 
32
      include 'reweight0.inc'
 
33
c
 
34
c     The type suffix of the histogram title, with syntax 
 
35
c     |T@<type_name> is semantic in the HwU format. It allows for
 
36
c     various filtering when using the histogram.py module
 
37
c     (see comment at the beginning of this file).
 
38
c     It is in general a good idea to keep the same title for the
 
39
c     same observable (if they use the same range) and differentiate
 
40
c     them only using the type suffix.
 
41
c
 
42
      character*8 HwUtype(2)
 
43
      data HwUtype/'|T@NOCUT','|T@CUT  '/
 
44
      integer j,jpr,i
 
45
      character*5 cc(2)
 
46
      data cc/'     ','Born '/
 
47
      integer nwgt,max_weight,nwgt_analysis,kk,l
 
48
      common/cnwgt/nwgt
 
49
      common/c_analysis/nwgt_analysis
 
50
      character*(wgts_info_len) weights_info(max_weight_shower)
 
51
      common/cwgtsinfo/weights_info
 
52
c Initialize histograms
 
53
      call HwU_inithist(nwgt,weights_info)
 
54
c Set method for error estimation to '0', i.e., use Poisson statistics
 
55
c for the uncertainty estimate
 
56
      call set_error_estimation(0)
 
57
      nwgt_analysis=nwgt
 
58
      do i=1,1
 
59
        l=(i-1)*4
 
60
        call HwU_book(l+1,'total rate  '//HwUtype(i),5,0.5d0,5.5d0)
 
61
        call HwU_book(l+2,'ta+ta- mass '//HwUtype(i),40,0d0,200d0)
 
62
        call HwU_book(l+3,'ta+ta- rap  '//HwUtype(i),40,-5d0,5d0)
 
63
        call HwU_book(l+4,'ta+ta- pt   '//HwUtype(i),20,0d0,400d0)
 
64
      enddo
 
65
 
 
66
      END
 
67
 
 
68
C----------------------------------------------------------------------
 
69
      SUBROUTINE HWAEND
 
70
C     USER''S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
71
C----------------------------------------------------------------------
 
72
      INCLUDE 'HERWIG65.INC'
 
73
      REAL*8 XNORM
 
74
      INTEGER I,J,KK,l,nwgt_analysis
 
75
      integer NPL
 
76
      parameter(NPL=15000)
 
77
      common/c_analysis/nwgt_analysis
 
78
c Convert from nb to pb using xnorm. This assumes that event weights
 
79
c *average* to the total cross section, so no extra weight needed
 
80
      xnorm=1d3
 
81
c Collect accumulated results
 
82
      call finalize_histograms(nevhep)
 
83
c Write the histograms to disk. 
 
84
      open (unit=99,file='MADatNLO.HwU',status='unknown')
 
85
      call HwU_output(99,xnorm)
 
86
      close (99)
 
87
      return
 
88
      END
 
89
 
 
90
C----------------------------------------------------------------------
 
91
      SUBROUTINE HWANAL
 
92
C     USER''S ROUTINE TO ANALYSE DATA FROM EVENT
 
93
C----------------------------------------------------------------------
 
94
      INCLUDE 'HERWIG65.INC'
 
95
      include 'reweight0.inc'
 
96
      DOUBLE PRECISION HWVDOT,PSUM(4),PH(5),YCUT,XMH,PTH,YH,THV,ETAV,
 
97
     #  PPL(5),PPLB(5),PTL,YL,THL,ETAL,PLL,ENL,PTLB,YLB,THLB,ETALB,
 
98
     #  PLLB,ENLB,PTPAIR,DLL,CLL,AZI,AZINORM,XMLL,DETALLB,VAR
 
99
      INTEGER ICHSUM,ICHINI,IHEP,IV,IFV,IST,ID,IJ,ID1,JPR,IDENT,
 
100
     #  ILL,ILLB,IHRD,ILL0,ILLB0
 
101
      LOGICAL DIDSOF,TEST1,TEST2,flag,ISLP,ISLM,FOUNDP,FOUNDM
 
102
      REAL*8 PI,wmass,wgamma,bwcutoff,getinvm,getdelphi,getrapidity,
 
103
     &getpseudorap
 
104
      PARAMETER (PI=3.14159265358979312D0)
 
105
      REAL*8 WWW0,TINY
 
106
      INTEGER KK,i,l
 
107
      DATA TINY/.1D-5/
 
108
      integer nwgt_analysis,max_weight
 
109
      common/c_analysis/nwgt_analysis
 
110
      integer maxRWGT
 
111
      parameter (maxRWGT=100)
 
112
      double precision wgtxsecRWGT(maxRWGT)
 
113
      parameter (max_weight=maxscales*maxscales+maxpdfs+maxRWGT+1)
 
114
      double precision ww(max_weight),www(max_weight)
 
115
      common/cww/ww
 
116
      IF (IERROR.NE.0) RETURN
 
117
      IF (WW(1).EQ.0D0) THEN
 
118
         WRITE(*,*)'WW(1) = 0. Stopping'
 
119
         STOP
 
120
      ENDIF
 
121
      IDENT=15
 
122
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT''S A POWER-SUPPRESSED
 
123
C EFFECT, SO THROW THE EVENT AWAY
 
124
      IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN
 
125
        CALL HWWARN('HWANAL',111)
 
126
        GOTO 999
 
127
      ENDIF
 
128
      DO I=1,nwgt_analysis
 
129
         WWW(I)=EVWGT*ww(i)/ww(1)
 
130
      ENDDO
 
131
      CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PSUM)
 
132
      CALL HWVSCA(4,-1D0,PSUM,PSUM)
 
133
      ICHSUM=0
 
134
      ICHINI=ICHRG(IDHW(1))+ICHRG(IDHW(2))
 
135
      DIDSOF=.FALSE.
 
136
      FOUNDP=.FALSE.
 
137
      FOUNDM=.FALSE.
 
138
      DO 100 IHEP=1,NHEP
 
139
         IF (IDHW(IHEP).EQ.16) DIDSOF=.TRUE.
 
140
         IF (ISTHEP(IHEP).EQ.1) THEN
 
141
            CALL HWVSUM(4,PHEP(1,IHEP),PSUM,PSUM)
 
142
            ICHSUM=ICHSUM+ICHRG(IDHW(IHEP))
 
143
         ENDIF
 
144
         IST=ISTHEP(IHEP)      
 
145
         ID=IDHW(IHEP)
 
146
         ID1=IDHEP(IHEP)
 
147
         ISLP=ID1.EQ.IDENT
 
148
         ISLM=ID1.EQ.-IDENT
 
149
         IF(((IST.GE.120.AND.IST.LE.125).OR.IST.EQ.1.OR.IST.EQ.198)
 
150
     &       .AND.ISLM.AND..NOT.FOUNDM)THEN
 
151
            ILL=IHEP
 
152
            FOUNDM=.TRUE.
 
153
         ENDIF
 
154
         IF(((IST.GE.120.AND.IST.LE.125).OR.IST.EQ.1.OR.IST.EQ.198)
 
155
     &       .AND.ISLP.AND..NOT.FOUNDP)THEN
 
156
            ILLB=IHEP
 
157
            FOUNDP=.TRUE.
 
158
         ENDIF
 
159
 100  CONTINUE
 
160
      IF(.NOT.FOUNDP.OR..NOT.FOUNDM)THEN
 
161
         WRITE(*,*)'NO TAUS FOUND.'
 
162
         STOP
 
163
      ENDIF
 
164
C CHECK MOMENTUM AND CHARGE CONSERVATION
 
165
      IF (HWVDOT(3,PSUM,PSUM).GT.1.E-4*PHEP(4,1)**2) THEN
 
166
         CALL HWUEPR
 
167
         CALL HWWARN('HWANAL',112)
 
168
         GOTO 999
 
169
      ENDIF
 
170
      IF (ICHSUM.NE.ICHINI) THEN
 
171
         CALL HWUEPR
 
172
         CALL HWWARN('HWANAL',113)
 
173
         GOTO 999
 
174
      ENDIF
 
175
C FIND THE LEPTONS
 
176
      DO IJ=1,5
 
177
        PPL(IJ)=PHEP(IJ,ILL)
 
178
        PPLB(IJ)=PHEP(IJ,ILLB)
 
179
        PH(IJ)=PPL(IJ)+PPLB(IJ)
 
180
      ENDDO
 
181
      xmh = getinvm(ph(4),ph(1),ph(2),ph(3))
 
182
      yh  = getrapidity(ph(4),ph(3))
 
183
      pth = sqrt(max(ph(1)**2+ph(2)**2,0d0))
 
184
      var = 1.d0
 
185
      do i=1,1
 
186
         l=(i-1)*4
 
187
         call HwU_fill(l+1,var,WWW)
 
188
         call HwU_fill(l+2,xmh,WWW)
 
189
         call HwU_fill(l+3,yh,WWW)
 
190
         call HwU_fill(l+4,pth,WWW)
 
191
      enddo
 
192
      call HwU_add_points
 
193
C
 
194
 999  END
 
195
 
 
196
 
 
197
      function getrapidity(en,pl)
 
198
      implicit none
 
199
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
200
      parameter (tiny=1.d-8)
 
201
      xplus=en+pl
 
202
      xminus=en-pl
 
203
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
204
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
205
            y=0.5d0*log( xplus/xminus  )
 
206
         else
 
207
            y=sign(1.d0,pl)*1.d8
 
208
         endif
 
209
      else 
 
210
         y=sign(1.d0,pl)*1.d8
 
211
      endif
 
212
      getrapidity=y
 
213
      return
 
214
      end
 
215
 
 
216
 
 
217
      function getpseudorap(en,ptx,pty,pl)
 
218
      implicit none
 
219
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
220
      parameter (tiny=1.d-5)
 
221
c
 
222
      pt=sqrt(ptx**2+pty**2)
 
223
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
224
        eta=sign(1.d0,pl)*1.d8
 
225
      else
 
226
        th=atan2(pt,pl)
 
227
        eta=-log(tan(th/2.d0))
 
228
      endif
 
229
      getpseudorap=eta
 
230
      return
 
231
      end
 
232
 
 
233
 
 
234
      function getinvm(en,ptx,pty,pl)
 
235
      implicit none
 
236
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
237
      parameter (tiny=1.d-5)
 
238
c
 
239
      tmp=en**2-ptx**2-pty**2-pl**2
 
240
      if(tmp.gt.0.d0)then
 
241
        tmp=sqrt(tmp)
 
242
      elseif(tmp.gt.-tiny)then
 
243
        tmp=0.d0
 
244
      else
 
245
        write(*,*)'Attempt to compute a negative mass'
 
246
        stop
 
247
      endif
 
248
      getinvm=tmp
 
249
      return
 
250
      end
 
251
 
 
252
 
 
253
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
254
      implicit none
 
255
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
256
      parameter (tiny=1.d-5)
 
257
c
 
258
      pt1=sqrt(ptx1**2+pty1**2)
 
259
      pt2=sqrt(ptx2**2+pty2**2)
 
260
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
261
        tmp=ptx1*ptx2+pty1*pty2
 
262
        tmp=tmp/(pt1*pt2)
 
263
        if(abs(tmp).gt.1.d0+tiny)then
 
264
          write(*,*)'Cosine larger than 1'
 
265
          stop
 
266
        elseif(abs(tmp).ge.1.d0)then
 
267
          tmp=sign(1.d0,tmp)
 
268
        endif
 
269
        tmp=acos(tmp)
 
270
      else
 
271
        tmp=1.d8
 
272
      endif
 
273
      getdelphi=tmp
 
274
      return
 
275
      end