← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_tj_hepmc.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

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