← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

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

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_tj.f
 
1
c
 
2
c Example analysis for "p p > t j [QCD]" process.
 
3
c Example analysis for "p p > t j $$ w+ w- [QCD]" process.
 
4
c Example analysis for "p p > w+ > t j  [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 HWABEG
 
15
C     USER'S ROUTINE FOR INITIALIZATION
 
16
C----------------------------------------------------------------------
 
17
      INCLUDE 'HERWIG65.INC'
 
18
      include 'reweight0.inc'
 
19
      integer j,kk,l
 
20
      character*5 cc(2)
 
21
      data cc/'     ','     '/
 
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
      call inihist
 
29
      nwgt_analysis=nwgt
 
30
      do j=1,2
 
31
      do kk=1,nwgt_analysis
 
32
      l=(kk-1)*48+(j-1)*24
 
33
      call mbook(l+ 1,'t pt        '//weights_info(kk)//cc(j)
 
34
     &     ,5d0,0d0,200d0)
 
35
      call mbook(l+ 2,'t log pt    '//weights_info(kk)//cc(j)
 
36
     &     ,0.05d0,0d0,5d0)
 
37
      call mbook(l+ 3,'t y         '//weights_info(kk)//cc(j)
 
38
     &     ,0.25d0,-6d0,6d0)
 
39
      call mbook(l+ 4,'t eta       '//weights_info(kk)//cc(j)
 
40
     &     ,0.25d0,-6d0,6d0)
 
41
c
 
42
      call mbook(l+ 5,'j1 pt       '//weights_info(kk)//cc(j)
 
43
     &     ,5d0,0d0,200d0)
 
44
      call mbook(l+ 6,'j1 log pt   '//weights_info(kk)//cc(j)
 
45
     &     ,0.05d0,0d0,5d0)
 
46
      call mbook(l+ 7,'j1 y        '//weights_info(kk)//cc(j)
 
47
     &     ,0.25d0,-6d0,6d0)
 
48
      call mbook(l+ 8,'j1 eta      '//weights_info(kk)//cc(j)
 
49
     &     ,0.25d0,-6d0,6d0)
 
50
c
 
51
      call mbook(l+ 9,'j2 pt       '//weights_info(kk)//cc(j)
 
52
     &     ,5d0,0d0,200d0)
 
53
      call mbook(l+10,'j2 log pt   '//weights_info(kk)//cc(j)
 
54
     &     ,0.05d0,0d0,5d0)
 
55
      call mbook(l+11,'j2 y        '//weights_info(kk)//cc(j)
 
56
     &     ,0.25d0,-6d0,6d0)
 
57
      call mbook(l+12,'j2 eta      '//weights_info(kk)//cc(j)
 
58
     &     ,0.25d0,-6d0,6d0)
 
59
c
 
60
      call mbook(l+13,'bj1 pt      '//weights_info(kk)//cc(j)
 
61
     &     ,5d0,0d0,200d0)
 
62
      call mbook(l+14,'bj1 log pt  '//weights_info(kk)//cc(j)
 
63
     &     ,0.05d0,0d0,5d0)
 
64
      call mbook(l+15,'bj1 y       '//weights_info(kk)//cc(j)
 
65
     &     ,0.25d0,-6d0,6d0)
 
66
      call mbook(l+16,'bj1 eta     '//weights_info(kk)//cc(j)
 
67
     &     ,0.25d0,-6d0,6d0)
 
68
c
 
69
      call mbook(l+17,'bj2 pt      '//weights_info(kk)//cc(j)
 
70
     &     ,5d0,0d0,200d0)
 
71
      call mbook(l+18,'bj2 log pt  '//weights_info(kk)//cc(j)
 
72
     &     ,0.05d0,0d0,5d0)
 
73
      call mbook(l+19,'bj2 y       '//weights_info(kk)//cc(j)
 
74
     &     ,0.25d0,-6d0,6d0)
 
75
      call mbook(l+20,'bj2 eta     '//weights_info(kk)//cc(j)
 
76
     &     ,0.25d0,-6d0,6d0)
 
77
c
 
78
      call mbook(l+21,'syst pt     '//weights_info(kk)//cc(j)
 
79
     &     ,5d0,0d0,200d0)
 
80
      call mbook(l+22,'syst log pt '//weights_info(kk)//cc(j)
 
81
     &     ,0.05d0,0d0,5d0)
 
82
      call mbook(l+23,'syst y      '//weights_info(kk)//cc(j)
 
83
     &     ,0.25d0,-6d0,6d0)
 
84
      call mbook(l+24,'syst eta    '//weights_info(kk)//cc(j)
 
85
     &     ,0.25d0,-6d0,6d0)
 
86
c
 
87
      enddo
 
88
      enddo
 
89
      END
 
90
 
 
91
 
 
92
C----------------------------------------------------------------------
 
93
      SUBROUTINE HWAEND
 
94
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
95
C----------------------------------------------------------------------
 
96
      INCLUDE 'HERWIG65.INC'
 
97
      REAL*8 XNORM
 
98
      INTEGER I,J,KK,l,nwgt_analysis
 
99
      integer NPL
 
100
      parameter(NPL=15000)
 
101
      common/c_analysis/nwgt_analysis
 
102
      OPEN(UNIT=99,FILE='HERST.top',STATUS='UNKNOWN')
 
103
C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG 
 
104
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
 
105
      XNORM=1.D3/DFLOAT(NEVHEP)
 
106
      DO I=1,NPL
 
107
        CALL MFINAL3(I)             
 
108
        CALL MCOPY(I,I+NPL)
 
109
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
110
        CALL MFINAL3(I+NPL)
 
111
      ENDDO
 
112
      do i=1,1
 
113
      do kk=1,nwgt_analysis
 
114
      l=(kk-1)*48+(i-1)*24
 
115
      call multitop(NPL+l+ 1,NPL-1,2,3,'t pt',' ','LOG')
 
116
      call multitop(NPL+l+ 2,NPL-1,2,3,'t log pt',' ','LOG')
 
117
      call multitop(NPL+l+ 3,NPL-1,2,3,'t y',' ','LOG')
 
118
      call multitop(NPL+l+ 4,NPL-1,2,3,'t eta',' ','LOG')
 
119
c
 
120
      call multitop(NPL+l+ 5,NPL-1,2,3,'j1 pt',' ','LOG')
 
121
      call multitop(NPL+l+ 6,NPL-1,2,3,'j1 log pt',' ','LOG')
 
122
      call multitop(NPL+l+ 7,NPL-1,2,3,'j1 y',' ','LOG')
 
123
      call multitop(NPL+l+ 8,NPL-1,2,3,'j1 eta',' ','LOG')
 
124
c
 
125
      call multitop(NPL+l+ 9,NPL-1,2,3,'j2 pt',' ','LOG')
 
126
      call multitop(NPL+l+10,NPL-1,2,3,'j2 log pt',' ','LOG')
 
127
      call multitop(NPL+l+11,NPL-1,2,3,'j2 y',' ','LOG')
 
128
      call multitop(NPL+l+12,NPL-1,2,3,'j2 eta',' ','LOG')
 
129
c
 
130
      call multitop(NPL+l+13,NPL-1,2,3,'bj1 pt',' ','LOG')
 
131
      call multitop(NPL+l+14,NPL-1,2,3,'bj1 log pt',' ','LOG')
 
132
      call multitop(NPL+l+15,NPL-1,2,3,'bj1 y',' ','LOG')
 
133
      call multitop(NPL+l+16,NPL-1,2,3,'bj1 eta',' ','LOG')
 
134
c
 
135
      call multitop(NPL+l+17,NPL-1,2,3,'bj2 pt',' ','LOG')
 
136
      call multitop(NPL+l+18,NPL-1,2,3,'bj2 log pt',' ','LOG')
 
137
      call multitop(NPL+l+19,NPL-1,2,3,'bj2 y',' ','LOG')
 
138
      call multitop(NPL+l+20,NPL-1,2,3,'bj2 eta',' ','LOG')
 
139
c
 
140
      call multitop(NPL+l+21,NPL-1,2,3,'syst pt',' ','LOG')
 
141
      call multitop(NPL+l+22,NPL-1,2,3,'syst log pt',' ','LOG')
 
142
      call multitop(NPL+l+23,NPL-1,2,3,'syst y',' ','LOG')
 
143
      call multitop(NPL+l+24,NPL-1,2,3,'syst eta',' ','LOG')
 
144
c
 
145
      enddo
 
146
      enddo
 
147
      CLOSE(99)
 
148
      END
 
149
 
 
150
 
 
151
C----------------------------------------------------------------------
 
152
      SUBROUTINE HWANAL
 
153
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
 
154
C     BASED ON AN ANALYSIS FILE WRITTEN BY E.RE
 
155
C----------------------------------------------------------------------
 
156
      INCLUDE 'HERWIG65.INC'
 
157
      include 'reweight0.inc'
 
158
      DOUBLE PRECISION HWVDOT,PSUM(4)
 
159
      REAL*8 PI
 
160
      PARAMETER (PI=3.14159265358979312D0)
 
161
      REAL*8 WWW0
 
162
      INTEGER kk,mu,jpart,i,j,ihep,ichsum,nt,nb,nbjet,ist,id,
 
163
     &njet,id1,i1,i2,ibmatch,ichini,k,ihadr,count_j,count_bj
 
164
      integer maxtrack,maxjet,maxnum,l
 
165
      parameter (maxtrack=2048,maxjet=2048,maxnum=30)
 
166
      integer ntracks,jetvec(maxtrack),ib1
 
167
      double precision pttop,etatop,ytop,ptj1,etaj1,yj1,ptbj1,etabj1,
 
168
     &ybj1,ptbj2,etabj2,ybj2,jet_ktradius,jet_ktptmin,palg,pttemp_spec,
 
169
     &pttemp_bjet,pttemp,tmp,getrapidity,getpseudorap,getpt,
 
170
     &pjet(4,maxtrack),ptrack(4,maxtrack),p_top(4,maxnum),p_b(4,maxnum),
 
171
     &p_bjet(4,maxnum),psyst(4),ptsyst,ysyst,etasyst,ptj2,yj2,etaj2
 
172
      logical is_b_jet(maxnum)
 
173
      integer btrack(maxnum),ib(maxnum)
 
174
      integer nwgt_analysis,max_weight
 
175
      common/c_analysis/nwgt_analysis
 
176
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
177
      double precision ww(max_weight),www(max_weight)
 
178
      common/cww/ww
 
179
c
 
180
      IF (IERROR.NE.0) RETURN
 
181
      IF (WW(1).EQ.0D0) THEN
 
182
         WRITE(*,*)'WW(1) = 0. Stopping'
 
183
         STOP
 
184
      ENDIF
 
185
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
 
186
C EFFECT, SO THROW THE EVENT AWAY
 
187
      IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN
 
188
        CALL HWWARN('HWANAL',111)
 
189
        GOTO 999
 
190
      ENDIF
 
191
      DO I=1,nwgt_analysis
 
192
         WWW(I)=EVWGT*ww(i)/ww(1)
 
193
      ENDDO
 
194
      CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PSUM)
 
195
      CALL HWVSCA(4,-1D0,PSUM,PSUM)
 
196
      ICHSUM=0
 
197
      ICHINI=ICHRG(IDHW(1))+ICHRG(IDHW(2))
 
198
 
 
199
C INITIALIZE
 
200
      NT=0
 
201
      NB=0
 
202
      NBJET=0
 
203
      NTRACKS=0
 
204
      NJET=0
 
205
 
 
206
      DO IHEP=1,NHEP
 
207
         IF (ISTHEP(IHEP).EQ.1) THEN
 
208
            CALL HWVSUM(4,PHEP(1,IHEP),PSUM,PSUM)
 
209
            ICHSUM=ICHSUM+ICHRG(IDHW(IHEP))
 
210
         ENDIF
 
211
         IST=ISTHEP(IHEP)      
 
212
         ID=IDHEP(IHEP)
 
213
         ID1=IHADR(ID) ! equal to the PDG of the massive quark in hadron
 
214
C TOP
 
215
        IF(IST.EQ.155.AND.ABS(ID).EQ.6) THEN
 
216
           DO MU=1,4
 
217
              P_TOP(MU,1)=PHEP(MU,IHEP)
 
218
           ENDDO
 
219
        ENDIF
 
220
c Define particles that go into jet. 
 
221
        IF (IST.EQ.1.AND.ABS(ID).GE.100)THEN
 
222
           NTRACKS=NTRACKS+1
 
223
           if (abs(id1).eq.5) THEN
 
224
c FOUND A stable B-FLAVOURED HADRON.
 
225
              NB=NB+1
 
226
              IB(NB)=IHEP
 
227
              DO MU=1,4
 
228
                 P_B(MU,NB)=PHEP(MU,IHEP)
 
229
              ENDDO
 
230
              BTRACK(NB)=NTRACKS
 
231
           endif
 
232
           DO MU=1,4
 
233
              PTRACK(MU,NTRACKS)=PHEP(MU,IHEP)
 
234
           ENDDO
 
235
           IF(NTRACKS.EQ.MAXTRACK) THEN
 
236
              WRITE(*,*)'HWANAL: TOO MANY PARTICLES, INCREASE MAXTRACK'
 
237
              STOP
 
238
           ENDIF
 
239
        ENDIF
 
240
      ENDDO
 
241
C END OF LOOP OVER IHEP
 
242
 
 
243
C CHECK MOMENTUM AND CHARGE CONSERVATION
 
244
      IF (HWVDOT(3,PSUM,PSUM).GT.1.E-4*PHEP(4,1)**2) THEN
 
245
         CALL HWUEPR
 
246
         CALL HWWARN('HWANAL',112)
 
247
         GOTO 999
 
248
      ENDIF
 
249
      IF (ICHSUM.NE.ICHINI) THEN
 
250
         CALL HWUEPR
 
251
         CALL HWWARN('HWANAL',113)
 
252
         GOTO 999
 
253
      ENDIF
 
254
 
 
255
      IF (NTRACKS.EQ.0) THEN
 
256
         WRITE(*,*) 'NO TRACKS FOUND, DROP ANALYSIS OF THIS EVENT'
 
257
         GOTO 999
 
258
      ENDIF
 
259
         
 
260
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
261
C KT ALGORITHM, FASTJET IMPLEMENTATION
 
262
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
263
      NJET=0
 
264
      JET_KTRADIUS = 0.7D0          
 
265
      JET_KTPTMIN  = 5D0
 
266
      PALG=1D0
 
267
      CALL fastjetppgenkt(PTRACK,NTRACKS,JET_KTRADIUS,JET_KTPTMIN,PALG,
 
268
     $     PJET,NJET,JETVEC)
 
269
 
 
270
c Check that jets are ordered in pt
 
271
      do i=1,njet-1
 
272
         if (getpt(pjet(1,i)).lt.getpt(pjet(1,i+1)) ) then
 
273
            write (*,*) 'ERROR jets not ordered'
 
274
            stop
 
275
         endif
 
276
      enddo
 
277
         
 
278
C b-jet 
 
279
      do i=1,njet
 
280
         is_b_jet(i)=.false.
 
281
         do j=1,NB
 
282
            if (JETVEC(BTRACK(j)).eq.i) then
 
283
               is_b_jet(i)=.true.
 
284
               exit
 
285
            endif
 
286
         enddo
 
287
      enddo
 
288
 
 
289
      pttop = getpt(p_top(1,1))
 
290
      ytop  = getrapidity(p_top(1,1))
 
291
      etatop= getpseudorap(p_top(1,1))
 
292
 
 
293
      count_j=0
 
294
      count_bj=0
 
295
      do i=1,njet
 
296
         if(.not.is_b_jet(i))then
 
297
            count_j=count_j+1
 
298
            if(count_j.eq.1)then
 
299
               ptj1 = getpt(pjet(1,i))
 
300
               yj1  = getrapidity(pjet(1,i))
 
301
               etaj1= getpseudorap(pjet(1,i))
 
302
               do mu=1,4
 
303
                  psyst(mu)=p_top(mu,1)+pjet(mu,i)
 
304
               enddo
 
305
               ptsyst = getpt(psyst)
 
306
               ysyst  = getrapidity(psyst)
 
307
               etasyst= getpseudorap(psyst)
 
308
            elseif(count_j.eq.2)then
 
309
               ptj2 = getpt(pjet(1,i))
 
310
               yj2  = getrapidity(pjet(1,i))
 
311
               etaj2= getpseudorap(pjet(1,i))
 
312
            endif
 
313
         elseif(is_b_jet(i))then
 
314
            count_bj=count_bj+1
 
315
            if (count_bj.eq.1) then
 
316
               ptbj1 = getpt(pjet(1,i))
 
317
               ybj1  = getrapidity(pjet(1,i))
 
318
               etabj1= getpseudorap(pjet(1,i))
 
319
            elseif (count_bj.eq.2) then
 
320
               ptbj2 = getpt(pjet(1,i))
 
321
               ybj2  = getrapidity(pjet(1,i))
 
322
               etabj2= getpseudorap(pjet(1,i))
 
323
            endif
 
324
         endif
 
325
      enddo
 
326
      nbjet=count_bj
 
327
c fill the histograms
 
328
      do i=1,2
 
329
         do kk=1,nwgt_analysis
 
330
            l=(kk-1)*48+(i-1)*24
 
331
            call mfill(l+1,pttop,www(kk))
 
332
            if(pttop.gt.0d0) call mfill(l+2,log10(pttop),www(kk))
 
333
            call mfill(l+3,ytop,www(kk))
 
334
            call mfill(l+4,etatop,www(kk))
 
335
            if(njet.ge.1)then
 
336
               call mfill(l+5,ptj1,www(kk))
 
337
               if (ptj1.gt.0d0) call mfill(l+6,log10(ptj1),www(kk))
 
338
               call mfill(l+7,yj1,www(kk))
 
339
               call mfill(l+8,etaj1,www(kk))
 
340
               call mfill(l+21,ptsyst,www(kk))
 
341
               if(ptsyst.gt.0d0) call mfill(l+22,log10(ptsyst),www(kk))
 
342
               call mfill(l+23,ysyst,www(kk))
 
343
               call mfill(l+24,etasyst,www(kk))
 
344
            endif
 
345
            if(njet.ge.2)then
 
346
               call mfill(l+9,ptj2,www(kk))
 
347
               if(ptj2.gt.0d0) call mfill(l+10,log10(ptj2),www(kk))
 
348
               call mfill(l+11,yj2,www(kk))
 
349
               call mfill(l+12,etaj2,www(kk))
 
350
            endif
 
351
            if(nbjet.ge.1)then
 
352
               call mfill(l+13,ptbj1,www(kk))
 
353
               if (ptbj1.gt.0d0) call mfill(l+14,log10(ptbj1),www(kk))
 
354
               call mfill(l+15,ybj1,www(kk))
 
355
               call mfill(l+16,etabj1,www(kk))
 
356
            endif
 
357
            if(nbjet.ge.2)then
 
358
               call mfill(l+17,ptbj2,www(kk))
 
359
               if (ptbj2.gt.0d0) call mfill(l+18,log10(ptbj2),www(kk))
 
360
               call mfill(l+19,ybj2,www(kk))
 
361
               call mfill(l+20,etabj2,www(kk))
 
362
            endif
 
363
         enddo
 
364
      enddo
 
365
      
 
366
 999  RETURN
 
367
      END
 
368
 
 
369
 
 
370
      function getrapidity(p)
 
371
      implicit none
 
372
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y,p(4)
 
373
      parameter (tiny=1.d-5)
 
374
c
 
375
      en=p(4)
 
376
      pl=p(3)
 
377
      xplus=en+pl
 
378
      xminus=en-pl
 
379
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
380
        if( (xplus/xminus).gt.tiny )then
 
381
          y=0.5d0*log( xplus/xminus )
 
382
        else
 
383
          y=sign(1.d0,pl)*1.d8
 
384
        endif
 
385
      else
 
386
        y=sign(1.d0,pl)*1.d8
 
387
      endif
 
388
      getrapidity=y
 
389
      return
 
390
      end
 
391
 
 
392
 
 
393
      function getpseudorap(p)
 
394
      implicit none
 
395
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th,p(4)
 
396
      parameter (tiny=1.d-5)
 
397
c
 
398
      en=p(4)
 
399
      ptx=p(1)
 
400
      pty=p(2)
 
401
      pl=p(3)
 
402
      pt=sqrt(ptx**2+pty**2)
 
403
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
404
        eta=sign(1.d0,pl)*1.d8
 
405
      else
 
406
        th=atan2(pt,pl)
 
407
        eta=-log(tan(th/2.d0))
 
408
      endif
 
409
      getpseudorap=eta
 
410
      return
 
411
      end
 
412
 
 
413
      function getpt(p)
 
414
      implicit none
 
415
      real*8 getpt,p(4)
 
416
      getpt=dsqrt(p(1)**2+p(2)**2)
 
417
      return
 
418
      end
 
419
 
 
420
      function getinvm(en,ptx,pty,pl)
 
421
      implicit none
 
422
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
423
      parameter (tiny=1.d-5)
 
424
c
 
425
      tmp=en**2-ptx**2-pty**2-pl**2
 
426
      if(tmp.gt.0.d0)then
 
427
        tmp=sqrt(tmp)
 
428
      elseif(tmp.gt.-tiny)then
 
429
        tmp=0.d0
 
430
      else
 
431
        write(*,*)'Attempt to compute a negative mass'
 
432
        stop
 
433
      endif
 
434
      getinvm=tmp
 
435
      return
 
436
      end
 
437
 
 
438
 
 
439
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
440
      implicit none
 
441
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
442
      parameter (tiny=1.d-5)
 
443
c
 
444
      pt1=sqrt(ptx1**2+pty1**2)
 
445
      pt2=sqrt(ptx2**2+pty2**2)
 
446
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
447
        tmp=ptx1*ptx2+pty1*pty2
 
448
        tmp=tmp/(pt1*pt2)
 
449
        if(abs(tmp).gt.1.d0+tiny)then
 
450
          write(*,*)'Cosine larger than 1'
 
451
          stop
 
452
        elseif(abs(tmp).ge.1.d0)then
 
453
          tmp=sign(1.d0,tmp)
 
454
        endif
 
455
        tmp=acos(tmp)
 
456
      else
 
457
        tmp=1.d8
 
458
      endif
 
459
      getdelphi=tmp
 
460
      return
 
461
      end
 
462
 
 
463
 
 
464
      FUNCTION IHADR(ID)
 
465
c Returns the PDG code of the heavier quark in the hadron of PDG code ID
 
466
      IMPLICIT NONE
 
467
      INTEGER IHADR,ID,ID1
 
468
C
 
469
      IF(ID.NE.0)THEN
 
470
        ID1=ABS(ID)
 
471
        IF(ID1.GT.10000)ID1=ID1-1000*INT(ID1/1000)
 
472
        IHADR=ID1/(10**INT(LOG10(DFLOAT(ID1))))
 
473
      ELSE
 
474
        IHADR=0
 
475
      ENDIF
 
476
      RETURN
 
477
      END
 
478
 
 
479