← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_V.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_V.f
 
1
c
 
2
c Example analysis for "p p > w+ [QCD]" process.
 
3
c Example analysis for "p p > w- [QCD]" process.
 
4
c Example analysis for "p p > z [QCD]" process.
 
5
c
 
6
C----------------------------------------------------------------------
 
7
      SUBROUTINE RCLOS()
 
8
C     DUMMY IF HBOOK IS USED
 
9
C----------------------------------------------------------------------
 
10
      END
 
11
 
 
12
 
 
13
C----------------------------------------------------------------------
 
14
      SUBROUTINE PYABEG
 
15
C     USER'S ROUTINE FOR INITIALIZATION
 
16
C----------------------------------------------------------------------
 
17
      implicit double precision(a-h, o-z)
 
18
      implicit integer(i-n)
 
19
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
 
20
      include 'reweight0.inc'
 
21
      real * 8 xm0,gam,xmlow,xmupp,bin
 
22
      real * 8 xmi,xms,pi
 
23
      PARAMETER (PI=3.14159265358979312D0)
 
24
      integer j,kk,l,jpr
 
25
      character*5 cc(2)
 
26
      data cc/'     ','     '/
 
27
      integer nwgt,max_weight,nwgt_analysis
 
28
      common/cnwgt/nwgt
 
29
      common/c_analysis/nwgt_analysis
 
30
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
31
      character*15 weights_info(max_weight)
 
32
      common/cwgtsinfo/weights_info
 
33
c
 
34
      call inihist
 
35
      nwgt_analysis=nwgt
 
36
      xmi=40.d0
 
37
      xms=120.d0
 
38
      bin=1.0d0
 
39
      do j=1,1
 
40
      do kk=1,nwgt_analysis
 
41
      l=(kk-1)*10+(j-1)*5
 
42
      call mbook(l+ 1,'V pt     '//weights_info(kk)//cc(j)
 
43
     &     ,2.d0,0.d0,200.d0)
 
44
      call mbook(l+ 2,'V log pt '//weights_info(kk)//cc(j)
 
45
     &     ,0.05d0,0.d0,5.d0)
 
46
      call mbook(l+ 3,'V y      '//weights_info(kk)//cc(j)
 
47
     &     ,0.25d0,-9.d0,9.d0)
 
48
      call mbook(l+ 4,'V eta    '//weights_info(kk)//cc(j)
 
49
     &     ,0.25d0,-9.d0,9.d0)
 
50
      call mbook(l+ 5,'mV       '//weights_info(kk)//cc(j)
 
51
     &     ,bin,xmi,xms)
 
52
      enddo
 
53
      enddo
 
54
 999  END
 
55
C----------------------------------------------------------------------
 
56
      SUBROUTINE PYAEND(IEVT)
 
57
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
58
C----------------------------------------------------------------------
 
59
      REAL*8 XNORM
 
60
      INTEGER I,J,KK,l,nwgt_analysis
 
61
      integer NPL
 
62
      parameter(NPL=15000)
 
63
      common/c_analysis/nwgt_analysis
 
64
      OPEN(UNIT=99,FILE='PYTSVB.TOP',STATUS='UNKNOWN')
 
65
      XNORM=1.D0/IEVT
 
66
      DO I=1,NPL              
 
67
        CALL MFINAL3(I)             
 
68
        CALL MCOPY(I,I+NPL)
 
69
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
70
        CALL MFINAL3(I+NPL)             
 
71
      ENDDO                          
 
72
C
 
73
      do i=1,1
 
74
      do kk=1,nwgt_analysis
 
75
      l=(kk-1)*10+(i-1)*5
 
76
      call multitop(NPL+l+ 1,NPL-1,3,2,'V pt     ',' ','LOG')
 
77
      call multitop(NPL+l+ 2,NPL-1,3,2,'V log pt ',' ','LOG')
 
78
      call multitop(NPL+l+ 3,NPL-1,3,2,'V y      ',' ','LOG')
 
79
      call multitop(NPL+l+ 4,NPL-1,3,2,'V eta    ',' ','LOG')
 
80
      call multitop(NPL+l+ 5,NPL-1,3,2,'mV       ',' ','LOG')
 
81
      enddo
 
82
      enddo
 
83
      CLOSE(99)
 
84
      END
 
85
 
 
86
C----------------------------------------------------------------------
 
87
      SUBROUTINE PYANAL
 
88
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
 
89
C----------------------------------------------------------------------
 
90
      implicit double precision(a-h, o-z)
 
91
      implicit integer(i-n)
 
92
      include 'reweight0.inc'
 
93
      DOUBLE PRECISION PSUM(4),XME,PPV(5),PPE(5),PPNU(5),
 
94
     # PPDCE(5),PPDCNU(5),WT,ETAEMIN(2),ETAEMAX(2),PTEMIN(2),
 
95
     # XMV,PTV,YV,GETRAPIDITY,PTE,THE,ETAE,PTNU,THNU,ETANU,
 
96
     # PTDCE,THDCE,ETADCE,PTDCNU,THDCNU,ETADCNU,ETAV,GETPSEUDORAP
 
97
      INTEGER ICHSUM,ICHINI,IHEP,JPR,IDENT,IFV,IST,ID,ID1,IHRD,IV,
 
98
     # IJ,IE,INU,J
 
99
      integer pychge
 
100
      external pydata
 
101
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
 
102
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
 
103
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
 
104
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
 
105
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
 
106
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
 
107
      LOGICAL DIDSOF,TEST1,TEST2,FLAG
 
108
      REAL*8 PI
 
109
      PARAMETER (PI=3.14159265358979312D0)
 
110
      REAL*8 TINY,p1(4),p2(4),pihep(4)
 
111
      INTEGER KK,i,l
 
112
      DATA TINY/.1D-5/
 
113
      DATA XME/5.11D-4/
 
114
      integer nwgt_analysis,max_weight
 
115
      common/c_analysis/nwgt_analysis
 
116
      DOUBLE PRECISION EVWEIGHT
 
117
      COMMON/CEVWEIGHT/EVWEIGHT
 
118
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
119
      double precision ww(max_weight),www(max_weight)
 
120
      common/cww/ww
 
121
      INTEGER IFAIL
 
122
      COMMON/CIFAIL/IFAIL
 
123
      SAVE INOBOSON
 
124
c
 
125
      IF(IFAIL.EQ.1)RETURN
 
126
      IF (WW(1).EQ.0D0) THEN
 
127
         WRITE(*,*)'WW(1) = 0. Stopping'
 
128
         STOP
 
129
      ENDIF
 
130
c
 
131
c CHOOSE IDENT=24 FOR W+, IDENT=-24 FOR W-, IDENT=23 FOR Z0
 
132
      IDENT=24
 
133
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
 
134
C EFFECT, SO THROW THE EVENT AWAY
 
135
      IF(SIGN(1.D0,P(3,3)).EQ.SIGN(1.D0,P(4,3)))THEN
 
136
         WRITE(*,*)'WARNING 502 IN PYANAL'
 
137
         GOTO 999
 
138
      ENDIF
 
139
      DO I=1,nwgt_analysis
 
140
         WWW(I)=EVWEIGHT*ww(i)/ww(1)
 
141
      ENDDO
 
142
      do i=1,4
 
143
         p1(i)=p(1,i)
 
144
         p2(i)=p(2,i)
 
145
      enddo
 
146
      CALL VVSUM(4,P1,P2,PSUM)
 
147
      CALL VSCA(4,-1D0,PSUM,PSUM)
 
148
      ICHSUM=0
 
149
      kf1=k(1,2)
 
150
      kf2=k(2,2)
 
151
      ICHINI=pychge(kf1)+pychge(kf2)
 
152
      IFV=0
 
153
      DO 100 IHEP=1,N
 
154
        do j=1,4
 
155
          pihep(j)=p(ihep,j)
 
156
        enddo
 
157
        IST=K(IHEP,1)      
 
158
        ID1=K(IHEP,2)
 
159
        IORI=K(IHEP,3)
 
160
        IF (IST.LE.10) THEN
 
161
          CALL VVSUM(4,PIHEP,PSUM,PSUM)
 
162
          ICHSUM=ICHSUM+pychge(ID1)
 
163
        ENDIF
 
164
        TEST1=IORI.EQ.0
 
165
        TEST2=ID1.EQ.IDENT
 
166
        IF(TEST1.AND.TEST2)IV0=IHEP
 
167
        TEST1=IORI.EQ.IV0
 
168
        IF(TEST1.AND.TEST2)THEN
 
169
          IV=IHEP
 
170
          IFV=IFV+1
 
171
          DO IJ=1,5
 
172
             PPV(IJ)=P(IHEP,ij)
 
173
          ENDDO
 
174
        ENDIF
 
175
  100 CONTINUE
 
176
      IF(IFV.EQ.0) THEN
 
177
         INOBOSON=INOBOSON+1
 
178
         WRITE(*,*)'WARNING 503 IN PYANAL: NO WEAK BOSON ',INOBOSON
 
179
         GOTO 999
 
180
      ENDIF
 
181
C CHECK MOMENTUM AND CHARGE CONSERVATION
 
182
      IF (VDOT(3,PSUM,PSUM).GT.1.E-4*P(1,4)**2) THEN
 
183
         WRITE(*,*)'WARNING 112 IN PYANAL'
 
184
         GOTO 999
 
185
      ENDIF
 
186
      IF (ICHSUM.NE.ICHINI) THEN
 
187
         WRITE(*,*)'WARNING 113 IN PYANAL'
 
188
         GOTO 999
 
189
      ENDIF
 
190
      IF(IFV.GT.1) THEN
 
191
         WRITE(*,*)'WARNING 55 IN PYANAL'
 
192
         GOTO 999
 
193
      ENDIF
 
194
C FILL THE HISTOS
 
195
C Variables of the vector boson
 
196
      xmv=ppv(5)
 
197
      ptv=sqrt(ppv(1)**2+ppv(2)**2)
 
198
      yv=getrapidity(ppv(4),ppv(3))
 
199
      etav=getpseudorap(ppv(4),ppv(1),ppv(2),ppv(3))
 
200
C
 
201
      do i=1,1
 
202
         do kk=1,nwgt_analysis
 
203
            l=(kk-1)*10+(i-1)*5
 
204
            call mfill(l+1,ptv,WWW(kk))
 
205
            if(ptv.gt.0) call mfill(l+2,log10(ptv),WWW(kk))
 
206
            call mfill(l+3,yv,WWW(kk))
 
207
            call mfill(l+4,etav,WWW(kk))
 
208
            call mfill(l+5,xmv,WWW(kk))
 
209
         enddo
 
210
      enddo
 
211
C
 
212
 999  END
 
213
 
 
214
 
 
215
      function getrapidity(en,pl)
 
216
      implicit none
 
217
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
218
      parameter (tiny=1.d-5)
 
219
c
 
220
      xplus=en+pl
 
221
      xminus=en-pl
 
222
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
223
        if( (xplus/xminus).gt.tiny )then
 
224
          y=0.5d0*log( xplus/xminus )
 
225
        else
 
226
          y=sign(1.d0,pl)*1.d8
 
227
        endif
 
228
      else
 
229
        y=sign(1.d0,pl)*1.d8
 
230
      endif
 
231
      getrapidity=y
 
232
      return
 
233
      end
 
234
 
 
235
      function getpseudorap(en,ptx,pty,pl)
 
236
      implicit none
 
237
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
238
      parameter (tiny=1.d-5)
 
239
c
 
240
      pt=sqrt(ptx**2+pty**2)
 
241
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
242
        eta=sign(1.d0,pl)*1.d8
 
243
      else
 
244
        th=atan2(pt,pl)
 
245
        eta=-log(tan(th/2.d0))
 
246
      endif
 
247
      getpseudorap=eta
 
248
      return
 
249
      end
 
250
 
 
251
 
 
252
C-----------------------------------------------------------------------
 
253
      SUBROUTINE VVSUM(N,P,Q,R)
 
254
C-----------------------------------------------------------------------
 
255
C    VECTOR SUM
 
256
C-----------------------------------------------------------------------
 
257
      IMPLICIT NONE
 
258
      INTEGER N,I
 
259
      DOUBLE PRECISION P(N),Q(N),R(N)
 
260
      DO 10 I=1,N
 
261
   10 R(I)=P(I)+Q(I)
 
262
      END
 
263
 
 
264
 
 
265
 
 
266
C-----------------------------------------------------------------------
 
267
      SUBROUTINE VSCA(N,C,P,Q)
 
268
C-----------------------------------------------------------------------
 
269
C     VECTOR TIMES SCALAR
 
270
C-----------------------------------------------------------------------
 
271
      IMPLICIT NONE
 
272
      INTEGER N,I
 
273
      DOUBLE PRECISION C,P(N),Q(N)
 
274
      DO 10 I=1,N
 
275
   10 Q(I)=C*P(I)
 
276
      END
 
277
 
 
278
 
 
279
 
 
280
C-----------------------------------------------------------------------
 
281
      FUNCTION VDOT(N,P,Q)
 
282
C-----------------------------------------------------------------------
 
283
C     VECTOR DOT PRODUCT
 
284
C-----------------------------------------------------------------------
 
285
      IMPLICIT NONE
 
286
      INTEGER N,I
 
287
      DOUBLE PRECISION VDOT,PQ,P(N),Q(N)
 
288
      PQ=0.
 
289
      DO 10 I=1,N
 
290
   10 PQ=PQ+P(I)*Q(I)
 
291
      VDOT=PQ
 
292
      END
 
293