← 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_lvl.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_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 PYABEG
 
18
C     USER'S ROUTINE FOR INITIALIZATION
 
19
C----------------------------------------------------------------------
 
20
      implicit none
 
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 PYAEND(IEVT)
 
60
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
61
C----------------------------------------------------------------------
 
62
      REAL*8 XNORM
 
63
      INTEGER I,J,KK,l,nwgt_analysis
 
64
      integer NPL
 
65
      parameter(NPL=15000)
 
66
      common/c_analysis/nwgt_analysis
 
67
      OPEN(UNIT=99,FILE='PYTLL.TOP',STATUS='UNKNOWN')
 
68
      XNORM=1.D0/IEVT
 
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,'lep rapidity ',' ','LIN')
 
81
         call multitop(NPL+l+3,NPL-1,3,2,'lep pt       ',' ','LOG')
 
82
         call multitop(NPL+l+4,NPL-1,3,2,'et miss      ',' ','LOG')
 
83
         call multitop(NPL+l+5,NPL-1,3,2,'trans. mass  ',' ','LOG')
 
84
         call multitop(NPL+l+6,NPL-1,3,2,'w rapidity   ',' ','LIN')
 
85
         call multitop(NPL+l+7,NPL-1,3,2,'w pt         ',' ','LOG')
 
86
         call multitop(NPL+l+8,NPL-1,3,2,'cphi[l,vl]   ',' ','LOG')
 
87
      enddo
 
88
      enddo
 
89
      CLOSE(99)
 
90
      END
 
91
 
 
92
C----------------------------------------------------------------------
 
93
      SUBROUTINE PYANAL
 
94
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
 
95
C----------------------------------------------------------------------
 
96
      implicit double precision(a-h, o-z)
 
97
      implicit integer(i-n)
 
98
      include 'reweight0.inc'
 
99
      DOUBLE PRECISION HWVDOT,PSUM(4),PPV(5),PTW,YW,YE,PPL(5),PPLB(5),
 
100
     & PTE,PLL,PTLB,PLLB,var,mtr,etmiss,cphi
 
101
      INTEGER ICHSUM,ICHINI,IHEP,IV,IFV,IST,ID,IJ,ID1,JPR,IDENT,
 
102
     #  ILL,ILLB,IHRD
 
103
      integer pychge
 
104
      double precision p1(4),p2(4),pihep(4)
 
105
      external pydata
 
106
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
 
107
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
 
108
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
 
109
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
 
110
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
 
111
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
 
112
      LOGICAL DIDSOF,FOUNDL,FOUNDN,ISL,ISN
 
113
      REAL*8 PI,getrapidity
 
114
      PARAMETER (PI=3.14159265358979312D0)
 
115
      REAL*8 WWW0,TINY,SIGNL,SIGNN
 
116
      INTEGER KK,I,L,IL,IN
 
117
      DATA TINY/.1D-5/
 
118
      integer nwgt_analysis,max_weight
 
119
      common/c_analysis/nwgt_analysis
 
120
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
121
      double precision ww(max_weight),www(max_weight)
 
122
      common/cww/ww
 
123
      DOUBLE PRECISION EVWEIGHT
 
124
      COMMON/CEVWEIGHT/EVWEIGHT
 
125
      INTEGER IFAIL
 
126
      COMMON/CIFAIL/IFAIL
 
127
c
 
128
      IF(IFAIL.EQ.1)RETURN
 
129
      IF (WW(1).EQ.0D0) THEN
 
130
         WRITE(*,*)'WW(1) = 0. Stopping'
 
131
         STOP
 
132
      ENDIF
 
133
C CHOOSE IDENT = 11 FOR E - NU_E
 
134
C        IDENT = 13 FOR MU - NU_MU
 
135
C        IDENT = 15 FOR TAU - NU_TAU
 
136
      IDENT=11
 
137
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
 
138
C EFFECT, SO THROW THE EVENT AWAY
 
139
      IF(SIGN(1.D0,P(3,3)).EQ.SIGN(1.D0,P(4,3)))THEN
 
140
         WRITE(*,*)'WARNING 111 IN PYANAL'
 
141
         GOTO 999
 
142
      ENDIF
 
143
      DO I=1,nwgt_analysis
 
144
         WWW(I)=EVWEIGHT*ww(i)/ww(1)
 
145
      ENDDO
 
146
      DO I=1,4
 
147
         P1(I)=P(1,I)
 
148
         P2(I)=P(2,I)
 
149
      ENDDO
 
150
      CALL VVSUM(4,P1,P2,PSUM)
 
151
      CALL VSCA(4,-1D0,PSUM,PSUM)
 
152
      ICHSUM=0
 
153
      KF1=K(1,2)
 
154
      KF2=K(2,2)
 
155
      ICHINI=PYCHGE(KF1)+PYCHGE(KF2)
 
156
C
 
157
      DO 100 IHEP=1,N
 
158
        DO J=1,4
 
159
          PIHEP(J)=P(IHEP,J)
 
160
        ENDDO
 
161
        IST =K(IHEP,1)
 
162
        ID1 =K(IHEP,2)
 
163
        IORI=K(IHEP,3)
 
164
        IF (IST.LE.10) THEN
 
165
          CALL VVSUM(4,PIHEP,PSUM,PSUM)
 
166
          ICHSUM=ICHSUM+PYCHGE(ID1)
 
167
        ENDIF
 
168
        ISL=ABS(ID1).EQ.IDENT
 
169
        ISN=ABS(ID1).EQ.IDENT+1
 
170
        IF(IORI.LE.10..AND.ISL)THEN
 
171
            IL=IHEP
 
172
            FOUNDL=.TRUE.
 
173
            SIGNL=SIGN(1D0,DBLE(ID1))
 
174
         ENDIF
 
175
         IF(IORI.LE.10.AND.ISN)THEN
 
176
            IN=IHEP
 
177
            FOUNDN=.TRUE.
 
178
            SIGNN=SIGN(1D0,DBLE(ID1))
 
179
         ENDIF
 
180
  100 CONTINUE
 
181
      IF(.NOT.FOUNDL.OR..NOT.FOUNDN)THEN
 
182
         WRITE(*,*)'NO LEPTONS FOUND.'
 
183
         WRITE(*,*)'CURRENTLY THIS ANALYSIS LOOKS FOR'
 
184
         IF(IDENT.EQ.11)WRITE(*,*)'E - NU_E'
 
185
         IF(IDENT.EQ.13)WRITE(*,*)'MU - NU_MU'
 
186
         IF(IDENT.EQ.15)WRITE(*,*)'TAU - NU_TAU'
 
187
         WRITE(*,*)'IF THIS IS NOT MEANT,'
 
188
         WRITE(*,*)'PLEASE CHANGE THE VALUE OF IDENT IN THIS FILE.'
 
189
         STOP
 
190
      ENDIF
 
191
      IF(SIGNN.EQ.SIGNL)THEN
 
192
         WRITE(*,*)'TWO SAME SIGN LEPTONS!'
 
193
         WRITE(*,*)IL,IN,SIGNL,SIGNN
 
194
         STOP
 
195
      ENDIF
 
196
C CHECK MOMENTUM AND CHARGE CONSERVATION
 
197
      IF (VDOT(3,PSUM,PSUM).GT.1.E-4*P(1,4)**2) THEN
 
198
         WRITE(*,*)'WARNING 112 IN PYANAL'
 
199
         GOTO 999
 
200
      ENDIF
 
201
      IF (ICHSUM.NE.ICHINI) THEN
 
202
         WRITE(*,*)'WARNING 113 IN PYANAL'
 
203
         GOTO 999
 
204
      ENDIF
 
205
      DO IJ=1,5
 
206
        PPL(IJ) =P(IN, IJ)
 
207
        PPLB(IJ)=P(IL,IJ)
 
208
        PPV(IJ)=PPL(IJ)+PPLB(IJ)
 
209
      ENDDO
 
210
      ye     = getrapidity(pplb(4), pplb(3))
 
211
      yw     = getrapidity(ppv(4), ppv(3))
 
212
      pte    = dsqrt(pplb(1)**2 + pplb(2)**2)
 
213
      ptw    = dsqrt(ppv(1)**2+ppv(2)**2)
 
214
      etmiss = dsqrt(ppl(1)**2 + ppl(2)**2)
 
215
      mtr    = dsqrt(2d0*pte*etmiss-2d0*ppl(1)*pplb(1)-2d0*ppl(2)*pplb(2))
 
216
      cphi   = (ppl(1)*pplb(1)+ppl(2)*pplb(2))/pte/etmiss
 
217
      var    = 1.d0
 
218
      do i=1,1
 
219
         do kk=1,nwgt_analysis
 
220
            l=(kk-1)*16+(i-1)*8
 
221
            call mfill(l+1,var,www(kk))
 
222
            call mfill(l+2,ye,www(kk))
 
223
            call mfill(l+3,pte,www(kk))
 
224
            call mfill(l+4,etmiss,www(kk))
 
225
            call mfill(l+5,mtr,www(kk))
 
226
            call mfill(l+6,yw,www(kk))
 
227
            call mfill(l+7,ptw,www(kk))
 
228
            call mfill(l+8,cphi,www(kk))
 
229
         enddo
 
230
      enddo
 
231
 999  END
 
232
      
 
233
 
 
234
C-----------------------------------------------------------------------
 
235
      SUBROUTINE VVSUM(N,P,Q,R)
 
236
C-----------------------------------------------------------------------
 
237
C    VECTOR SUM
 
238
C-----------------------------------------------------------------------
 
239
      IMPLICIT NONE
 
240
      INTEGER N,I
 
241
      DOUBLE PRECISION P(N),Q(N),R(N)
 
242
      DO 10 I=1,N
 
243
   10 R(I)=P(I)+Q(I)
 
244
      END
 
245
 
 
246
 
 
247
 
 
248
C-----------------------------------------------------------------------
 
249
      SUBROUTINE VSCA(N,C,P,Q)
 
250
C-----------------------------------------------------------------------
 
251
C     VECTOR TIMES SCALAR
 
252
C-----------------------------------------------------------------------
 
253
      IMPLICIT NONE
 
254
      INTEGER N,I
 
255
      DOUBLE PRECISION C,P(N),Q(N)
 
256
      DO 10 I=1,N
 
257
   10 Q(I)=C*P(I)
 
258
      END
 
259
 
 
260
 
 
261
 
 
262
C-----------------------------------------------------------------------
 
263
      FUNCTION VDOT(N,P,Q)
 
264
C-----------------------------------------------------------------------
 
265
C     VECTOR DOT PRODUCT
 
266
C-----------------------------------------------------------------------
 
267
      IMPLICIT NONE
 
268
      INTEGER N,I
 
269
      DOUBLE PRECISION VDOT,PQ,P(N),Q(N)
 
270
      PQ=0.
 
271
      DO 10 I=1,N
 
272
   10 PQ=PQ+P(I)*Q(I)
 
273
      VDOT=PQ
 
274
      END
 
275
 
 
276
      function getrapidity(en,pl)
 
277
      implicit none
 
278
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
279
      parameter (tiny=1.d-8)
 
280
      xplus=en+pl
 
281
      xminus=en-pl
 
282
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
283
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
284
            y=0.5d0*log( xplus/xminus  )
 
285
         else
 
286
            y=sign(1.d0,pl)*1.d8
 
287
         endif
 
288
      else 
 
289
         y=sign(1.d0,pl)*1.d8
 
290
      endif
 
291
      getrapidity=y
 
292
      return
 
293
      end