← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

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

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lvl_hepmc.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 HWABEG
 
18
C     USER'S ROUTINE FOR INITIALIZATION
 
19
C----------------------------------------------------------------------
 
20
      INCLUDE 'HEPMC.INC'
 
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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
34
c To be changed !!
 
35
      nwgt=1
 
36
      weights_info(nwgt)="central value  "
 
37
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
38
      nwgt_analysis=nwgt
 
39
      do i=1,1
 
40
      do kk=1,nwgt_analysis
 
41
        l=(kk-1)*16+(i-1)*8
 
42
        call mbook(l+1,'total rate '//weights_info(kk)//cc(i),
 
43
     &       1.0d0,0.5d0,5.5d0)
 
44
        call mbook(l+2,'lep rapidity '//weights_info(kk)//cc(i),
 
45
     &       0.5d0,-5d0,5d0)
 
46
        call mbook(l+3,'lep pt '//weights_info(kk)//cc(i),
 
47
     &       10d0,0d0,200d0)
 
48
        call mbook(l+4,'et miss '//weights_info(kk)//cc(i),
 
49
     &       10d0,0d0,200d0)
 
50
        call mbook(l+5,'trans. mass '//weights_info(kk)//cc(i),
 
51
     &       5d0,0d0,200d0)
 
52
        call mbook(l+6,'w rapidity '//weights_info(kk)//cc(i),
 
53
     &       0.5d0,-5d0,5d0)
 
54
        call mbook(l+7,'w pt '//weights_info(kk)//cc(i),
 
55
     &       10d00,0d0,200d0)
 
56
        call mbook(l+8,'cphi[l,vl] '//weights_info(kk)//cc(i),
 
57
     &       0.05d0,-1d0,1d0)
 
58
      enddo
 
59
      enddo
 
60
 999  END
 
61
 
 
62
 
 
63
C----------------------------------------------------------------------
 
64
      SUBROUTINE HWAEND
 
65
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
66
C----------------------------------------------------------------------
 
67
      INCLUDE 'HEPMC.INC'
 
68
      REAL*8 XNORM
 
69
      INTEGER I,J,KK,l,nwgt_analysis
 
70
      integer NPL
 
71
      parameter(NPL=15000)
 
72
      common/c_analysis/nwgt_analysis
 
73
      OPEN(UNIT=99,FILE='HERLL.TOP',STATUS='UNKNOWN')
 
74
C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG 
 
75
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
 
76
      XNORM=1.D3/DFLOAT(NEVHEP)
 
77
      DO I=1,NPL
 
78
        CALL MFINAL3(I)             
 
79
        CALL MCOPY(I,I+NPL)
 
80
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
81
        CALL MFINAL3(I+NPL)             
 
82
      ENDDO                          
 
83
C
 
84
      do i=1,1
 
85
      do kk=1,nwgt_analysis
 
86
         l=(kk-1)*16+(i-1)*8
 
87
         call multitop(NPL+l+1,NPL-1,3,2,'total rate   ',' ','LIN')
 
88
         call multitop(NPL+l+2,NPL-1,3,2,'lep rapidity ',' ','LIN')
 
89
         call multitop(NPL+l+3,NPL-1,3,2,'lep pt       ',' ','LOG')
 
90
         call multitop(NPL+l+4,NPL-1,3,2,'et miss      ',' ','LOG')
 
91
         call multitop(NPL+l+5,NPL-1,3,2,'trans. mass  ',' ','LOG')
 
92
         call multitop(NPL+l+6,NPL-1,3,2,'w rapidity   ',' ','LIN')
 
93
         call multitop(NPL+l+7,NPL-1,3,2,'w pt         ',' ','LOG')
 
94
         call multitop(NPL+l+8,NPL-1,3,2,'cphi[l,vl]   ',' ','LOG')
 
95
      enddo
 
96
      enddo
 
97
      CLOSE(99)
 
98
      END
 
99
 
 
100
C----------------------------------------------------------------------
 
101
      SUBROUTINE HWANAL
 
102
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
 
103
C----------------------------------------------------------------------
 
104
      INCLUDE 'HEPMC.INC'
 
105
      include 'reweight0.inc'
 
106
      DOUBLE PRECISION HWVDOT,PSUM(4),PPV(5),PTW,YW,YE,PPL(5),PPLB(5),
 
107
     & PTE,PLL,PTLB,PLLB,var,mtr,etmiss,cphi
 
108
      INTEGER ICHSUM,ICHINI,IHEP,IV,IFV,IST,ID,IJ,ID1,JPR,IDENT,
 
109
     #  ILL,ILLB,IHRD,NL,NN
 
110
      LOGICAL DIDSOF,FOUNDL,FOUNDN,ISL,ISN
 
111
      REAL*8 PI,getrapidity
 
112
      PARAMETER (PI=3.14159265358979312D0)
 
113
      REAL*8 WWW0,TINY,SIGNL,SIGNN
 
114
      INTEGER KK,I,L,IL,IN
 
115
      DATA TINY/.1D-5/
 
116
      integer nwgt_analysis,max_weight
 
117
      common/c_analysis/nwgt_analysis
 
118
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
119
      double precision ww(max_weight),www(max_weight)
 
120
      common/cww/ww
 
121
c
 
122
      IF(MOD(NEVHEP,10000).EQ.0)RETURN
 
123
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
124
c To be changed !!
 
125
      ww(1)=1d0
 
126
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
127
      IF (WW(1).EQ.0D0) THEN
 
128
         WRITE(*,*)'WW(1) = 0. Stopping'
 
129
         STOP
 
130
      ENDIF
 
131
C CHOOSE IDENT = 11 FOR E - NU_E
 
132
C        IDENT = 13 FOR MU - NU_MU
 
133
C        IDENT = 15 FOR TAU - NU_TAU
 
134
      IDENT=11
 
135
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
 
136
C EFFECT, SO THROW THE EVENT AWAY
 
137
      IF(SIGN(1.D0,PHEP(3,1)).EQ.SIGN(1.D0,PHEP(3,2)))THEN
 
138
        CALL HWWARN('HWANAL',111)
 
139
        GOTO 999
 
140
      ENDIF
 
141
      DO I=1,nwgt_analysis
 
142
         WWW(I)=EVWGT*ww(i)/ww(1)
 
143
      ENDDO
 
144
      ICHSUM=0
 
145
      DIDSOF=.FALSE.
 
146
      FOUNDL=.FALSE.
 
147
      FOUNDN=.FALSE.
 
148
      NL=0
 
149
      NN=0
 
150
      DO 100 IHEP=1,NHEP
 
151
        IST=ISTHEP(IHEP)      
 
152
        ID1=IDHEP(IHEP)
 
153
        ISL=ABS(ID1).EQ.IDENT
 
154
        ISN=ABS(ID1).EQ.IDENT+1
 
155
        IF(IST.EQ.1.AND.NL.EQ.0.AND.ISL)THEN
 
156
            NL=NL+1
 
157
            IL=IHEP
 
158
            FOUNDL=.TRUE.
 
159
            SIGNL=SIGN(1D0,DBLE(ID1))
 
160
         ENDIF
 
161
         IF(IST.EQ.1.AND.NN.EQ.0.AND.ISN)THEN
 
162
            NN=NN+1
 
163
            IN=IHEP
 
164
            FOUNDN=.TRUE.
 
165
            SIGNN=SIGN(1D0,DBLE(ID1))
 
166
         ENDIF
 
167
  100 CONTINUE
 
168
      IF(.NOT.FOUNDL.OR..NOT.FOUNDN)THEN
 
169
         WRITE(*,*)'NO LEPTONS FOUND.'
 
170
         WRITE(*,*)'CURRENTLY THIS ANALYSIS LOOKS FOR'
 
171
         IF(IDENT.EQ.11)WRITE(*,*)'E - NU_E'
 
172
         IF(IDENT.EQ.13)WRITE(*,*)'MU - NU_MU'
 
173
         IF(IDENT.EQ.15)WRITE(*,*)'TAU - NU_TAU'
 
174
         WRITE(*,*)'IF THIS IS NOT MEANT,'
 
175
         WRITE(*,*)'PLEASE CHANGE THE VALUE OF IDENT IN THIS FILE.'
 
176
         STOP
 
177
      ENDIF
 
178
      IF(SIGNN.EQ.SIGNL)THEN
 
179
         WRITE(*,*)'TWO SAME SIGN LEPTONS!'
 
180
         WRITE(*,*)IL,IN,SIGNL,SIGNN
 
181
         STOP
 
182
      ENDIF
 
183
      DO IJ=1,5
 
184
        PPL(IJ)=PHEP(IJ,IN)
 
185
        PPLB(IJ)=PHEP(IJ,IL)
 
186
        PPV(IJ)=PPL(IJ)+PPLB(IJ)
 
187
      ENDDO
 
188
      ye     = getrapidity(pplb(4), pplb(3))
 
189
      yw     = getrapidity(ppv(4), ppv(3))
 
190
      pte    = dsqrt(pplb(1)**2 + pplb(2)**2)
 
191
      ptw    = dsqrt(ppv(1)**2+ppv(2)**2)
 
192
      etmiss = dsqrt(ppl(1)**2 + ppl(2)**2)
 
193
      mtr    = dsqrt(2d0*pte*etmiss-2d0*ppl(1)*pplb(1)-2d0*ppl(2)*pplb(2))
 
194
      cphi   = (ppl(1)*pplb(1)+ppl(2)*pplb(2))/pte/etmiss
 
195
      var    = 1.d0
 
196
      do i=1,1
 
197
         do kk=1,nwgt_analysis
 
198
            l=(kk-1)*16+(i-1)*8
 
199
            call mfill(l+1,var,www(kk))
 
200
            call mfill(l+2,ye,www(kk))
 
201
            call mfill(l+3,pte,www(kk))
 
202
            call mfill(l+4,etmiss,www(kk))
 
203
            call mfill(l+5,mtr,www(kk))
 
204
            call mfill(l+6,yw,www(kk))
 
205
            call mfill(l+7,ptw,www(kk))
 
206
            call mfill(l+8,cphi,www(kk))
 
207
         enddo
 
208
      enddo
 
209
 999  END
 
210
 
 
211
 
 
212
C-----------------------------------------------------------------------
 
213
      SUBROUTINE HWWARN(SUBRTN,ICODE)
 
214
C-----------------------------------------------------------------------
 
215
C     DEALS WITH ERRORS DURING EXECUTION
 
216
C     SUBRTN = NAME OF CALLING SUBROUTINE
 
217
C     ICODE  = ERROR CODE:    - -1 NONFATAL, KILL EVENT & PRINT NOTHING
 
218
C                            0- 49 NONFATAL, PRINT WARNING & CONTINUE
 
219
C                           50- 99 NONFATAL, PRINT WARNING & JUMP
 
220
C                          100-199 NONFATAL, DUMP & KILL EVENT
 
221
C                          200-299    FATAL, TERMINATE RUN
 
222
C                          300-399    FATAL, DUMP EVENT & TERMINATE RUN
 
223
C                          400-499    FATAL, DUMP EVENT & STOP DEAD
 
224
C                          500-       FATAL, STOP DEAD WITH NO DUMP
 
225
C-----------------------------------------------------------------------
 
226
      INCLUDE 'HEPMC.INC'
 
227
      INTEGER ICODE,NRN,IERROR
 
228
      CHARACTER*6 SUBRTN
 
229
      IF (ICODE.GE.0) WRITE (6,10) SUBRTN,ICODE
 
230
   10 FORMAT(/' HWWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4)
 
231
      IF (ICODE.LT.0) THEN
 
232
         IERROR=ICODE
 
233
         RETURN
 
234
      ELSEIF (ICODE.LT.100) THEN
 
235
         WRITE (6,20) NEVHEP,NRN,EVWGT
 
236
   20    FORMAT(' EVENT',I8,':   SEEDS =',I11,' &',I11,
 
237
     &'  WEIGHT =',E11.4/' EVENT SURVIVES. EXECUTION CONTINUES')
 
238
         IF (ICODE.GT.49) RETURN
 
239
      ELSEIF (ICODE.LT.200) THEN
 
240
         WRITE (6,30) NEVHEP,NRN,EVWGT
 
241
   30    FORMAT(' EVENT',I8,':   SEEDS =',I11,' &',I11,
 
242
     &'  WEIGHT =',E11.4/' EVENT KILLED.   EXECUTION CONTINUES')
 
243
         IERROR=ICODE
 
244
         RETURN
 
245
      ELSEIF (ICODE.LT.300) THEN
 
246
         WRITE (6,40)
 
247
   40    FORMAT(' EVENT SURVIVES.  RUN ENDS GRACEFULLY')
 
248
c$$$         CALL HWEFIN
 
249
c$$$         CALL HWAEND
 
250
         STOP
 
251
      ELSEIF (ICODE.LT.400) THEN
 
252
         WRITE (6,50)
 
253
   50    FORMAT(' EVENT KILLED: DUMP FOLLOWS.  RUN ENDS GRACEFULLY')
 
254
         IERROR=ICODE
 
255
c$$$         CALL HWUEPR
 
256
c$$$         CALL HWUBPR
 
257
c$$$         CALL HWEFIN
 
258
c$$$         CALL HWAEND
 
259
         STOP
 
260
      ELSEIF (ICODE.LT.500) THEN
 
261
         WRITE (6,60)
 
262
   60    FORMAT(' EVENT KILLED: DUMP FOLLOWS.  RUN STOPS DEAD')
 
263
         IERROR=ICODE
 
264
c$$$         CALL HWUEPR
 
265
c$$$         CALL HWUBPR
 
266
         STOP
 
267
      ELSE
 
268
         WRITE (6,70)
 
269
   70    FORMAT(' RUN CANNOT CONTINUE')
 
270
         STOP
 
271
      ENDIF
 
272
      END
 
273
 
 
274
 
 
275
      subroutine HWUEPR
 
276
      INCLUDE 'HEPMC.INC'
 
277
      integer ip,i
 
278
      PRINT *,' EVENT ',NEVHEP
 
279
      DO IP=1,NHEP
 
280
         PRINT '(I4,I8,I4,4I4,1P,5D11.3)',IP,IDHEP(IP),ISTHEP(IP),
 
281
     &        JMOHEP(1,IP),JMOHEP(2,IP),JDAHEP(1,IP),JDAHEP(2,IP),
 
282
     &        (PHEP(I,IP),I=1,5)
 
283
      ENDDO
 
284
      return
 
285
      end
 
286
 
 
287
      function getrapidity(en,pl)
 
288
      implicit none
 
289
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
290
      parameter (tiny=1.d-8)
 
291
      xplus=en+pl
 
292
      xminus=en-pl
 
293
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
294
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
295
            y=0.5d0*log( xplus/xminus  )
 
296
         else
 
297
            y=sign(1.d0,pl)*1.d8
 
298
         endif
 
299
      else 
 
300
         y=sign(1.d0,pl)*1.d8
 
301
      endif
 
302
      getrapidity=y
 
303
      return
 
304
      end
 
305
 
 
306
 
 
307
      function getpseudorap(en,ptx,pty,pl)
 
308
      implicit none
 
309
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
310
      parameter (tiny=1.d-5)
 
311
c
 
312
      pt=sqrt(ptx**2+pty**2)
 
313
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
314
        eta=sign(1.d0,pl)*1.d8
 
315
      else
 
316
        th=atan2(pt,pl)
 
317
        eta=-log(tan(th/2.d0))
 
318
      endif
 
319
      getpseudorap=eta
 
320
      return
 
321
      end
 
322
 
 
323
 
 
324
      function getinvm(en,ptx,pty,pl)
 
325
      implicit none
 
326
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
327
      parameter (tiny=1.d-5)
 
328
c
 
329
      tmp=en**2-ptx**2-pty**2-pl**2
 
330
      if(tmp.gt.0.d0)then
 
331
        tmp=sqrt(tmp)
 
332
      elseif(tmp.gt.-tiny)then
 
333
        tmp=0.d0
 
334
      else
 
335
        write(*,*)'Attempt to compute a negative mass'
 
336
        stop
 
337
      endif
 
338
      getinvm=tmp
 
339
      return
 
340
      end
 
341
 
 
342
 
 
343
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
344
      implicit none
 
345
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
346
      parameter (tiny=1.d-5)
 
347
c
 
348
      pt1=sqrt(ptx1**2+pty1**2)
 
349
      pt2=sqrt(ptx2**2+pty2**2)
 
350
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
351
        tmp=ptx1*ptx2+pty1*pty2
 
352
        tmp=tmp/(pt1*pt2)
 
353
        if(abs(tmp).gt.1.d0+tiny)then
 
354
          write(*,*)'Cosine larger than 1'
 
355
          stop
 
356
        elseif(abs(tmp).ge.1.d0)then
 
357
          tmp=sign(1.d0,tmp)
 
358
        endif
 
359
        tmp=acos(tmp)
 
360
      else
 
361
        tmp=1.d8
 
362
      endif
 
363
      getdelphi=tmp
 
364
      return
 
365
      end