← 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_leptons_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_leptons_hepmc.f
 
1
c
 
2
c Sample analysis for leptonic final states
 
3
c
 
4
C----------------------------------------------------------------------
 
5
      SUBROUTINE RCLOS()
 
6
C     DUMMY IF HBOOK IS USED
 
7
C----------------------------------------------------------------------
 
8
      END
 
9
 
 
10
 
 
11
C----------------------------------------------------------------------
 
12
      SUBROUTINE HWABEG
 
13
C     USER''S ROUTINE FOR INITIALIZATION
 
14
C----------------------------------------------------------------------
 
15
      INCLUDE 'HEPMC.INC'
 
16
      include 'reweight0.inc'
 
17
      integer nwgt,max_weight,nwgt_analysis
 
18
      common/cnwgt/nwgt
 
19
      common/c_analysis/nwgt_analysis
 
20
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
21
      character*15 weights_info(max_weight)
 
22
      common/cwgtsinfo/weights_info
 
23
      integer nsingle,ncorr,nlepton,nplots,ncuts
 
24
      common/cplots/nsingle,ncorr,nlepton,nplots,ncuts
 
25
      integer MAXELM,MAXELP,MAXMUM,MAXMUP
 
26
      common/cmaxlep/MAXELM,MAXELP,MAXMUM,MAXMUP
 
27
      character*60 tmpstr1,tmpstr2,tmpstr3
 
28
      integer maxcuts
 
29
      parameter (maxcuts=2)
 
30
      character*5 cc(maxcuts)
 
31
      data cc/'     ',' cuts'/
 
32
      real*8 pi
 
33
      PARAMETER (PI=3.14159265358979312D0)
 
34
      real*8 sbin(100),smin(100),smax(100)
 
35
      real*8 cbin(100),cmin(100),cmax(100)
 
36
      integer i,kk,l,icuts
 
37
      integer l0,ilep1,ilep2,io,ipair
 
38
c      
 
39
      call inihist
 
40
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
41
c To be changed !!
 
42
      nwgt=1
 
43
      weights_info(nwgt)="central value  "
 
44
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
45
      nwgt_analysis=nwgt
 
46
c The analysis will consider:
 
47
c  MAXELM e-
 
48
c  MAXELP e+
 
49
c  MAXMUM mu-
 
50
c  MAXMUP mu+
 
51
c Set these variables here and only here
 
52
      MAXELM=0
 
53
      MAXELP=0
 
54
      MAXMUM=1
 
55
      MAXMUP=1
 
56
      nlepton=MAXELM+MAXELP+MAXMUM+MAXMUP
 
57
c For each weight and cut configuration, there will be:
 
58
c   nsingle single-inclusive plots (e.g., pt and y)
 
59
c   ncorr correlation plots (e.g., invM, ptpair, dphi, Deltay)
 
60
c to be repeated for each of the nlepton leptons and 
 
61
c nlepton*(nlepton-1)/2 lepton pairs respectively
 
62
      ncuts=1
 
63
      if(ncuts.gt.maxcuts)then
 
64
        WRITE(*,*)ncuts,maxcuts
 
65
        CALL HWWARN('HWABEG',500)
 
66
      endif
 
67
      nsingle=2
 
68
      ncorr=4
 
69
      nplots=nlepton * nsingle +
 
70
     #       nlepton*(nlepton-1)/2 * ncorr
 
71
      do i=1,100
 
72
        sbin(i)=0.d0
 
73
        smin(i)=0.d0
 
74
        smax(i)=0.d0
 
75
        cbin(i)=0.d0
 
76
        cmin(i)=0.d0
 
77
        cmax(i)=0.d0
 
78
      enddo
 
79
c pt
 
80
      sbin(1)=2.d0
 
81
      smin(1)=0.d0
 
82
      smax(1)=200.d0
 
83
c y
 
84
      sbin(2)=0.1d0
 
85
      smin(2)=-5.d0
 
86
      smax(2)=5.d0
 
87
c inv M
 
88
      cbin(1)=5.d0
 
89
      cmin(1)=0.d0
 
90
      cmax(1)=500.d0
 
91
c ptpair
 
92
      cbin(2)=2.d0
 
93
      cmin(2)=0.d0
 
94
      cmax(2)=200.d0
 
95
c dphi
 
96
      cbin(3)=pi/40.d0
 
97
      cmin(3)=0.d0
 
98
      cmax(3)=pi
 
99
c Delta y
 
100
      cbin(4)=0.1d0
 
101
      cmin(4)=-5.d0
 
102
      cmax(4)=5.d0
 
103
c
 
104
      do kk=1,nwgt_analysis
 
105
      do icuts=1,ncuts
 
106
        l0=(kk-1)*ncuts*nplots+(icuts-1)*nplots
 
107
        do ilep1=1,nlepton
 
108
        do io=1,nsingle
 
109
          l=l0+nsingle*(ilep1-1)+io
 
110
          write(tmpstr1,'(i3)')ilep1
 
111
          write(tmpstr3,'(i3)')io
 
112
          call mbook(l,'sing '//tmpstr1(1:3)//
 
113
     &    ' '//tmpstr3(1:3)//' '
 
114
     &       //weights_info(kk)//cc(icuts),sbin(io),smin(io),smax(io))
 
115
        enddo
 
116
        enddo
 
117
        ipair=0
 
118
        do ilep1=1,nlepton-1
 
119
        do ilep2=ilep1+1,nlepton
 
120
        ipair=ipair+1
 
121
        do io=1,ncorr
 
122
          l=l0+nlepton*nsingle+ncorr*(ipair-1)+io
 
123
          write(tmpstr1,'(i3)')ilep1
 
124
          write(tmpstr2,'(i3)')ilep2
 
125
          write(tmpstr3,'(i3)')io
 
126
          call mbook(l,'corr '//tmpstr1(1:3)//' '//
 
127
     &       tmpstr2(1:3)//' '//tmpstr3(1:3)//
 
128
     &     ' '//weights_info(kk)//cc(icuts),cbin(io),cmin(io),cmax(io))
 
129
        enddo
 
130
        enddo
 
131
        enddo
 
132
      enddo
 
133
      enddo
 
134
 999  END
 
135
 
 
136
C----------------------------------------------------------------------
 
137
      SUBROUTINE HWAEND
 
138
C     USER''S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
 
139
C----------------------------------------------------------------------
 
140
      INCLUDE 'HEPMC.INC'
 
141
      REAL*8 XNORM
 
142
      integer nsingle,ncorr,nlepton,nplots,ncuts
 
143
      common/cplots/nsingle,ncorr,nlepton,nplots,ncuts
 
144
      INTEGER I,J,kk,l,nwgt_analysis
 
145
      integer l0,ilep1,ilep2,io,ipair,icuts
 
146
      integer NPL
 
147
      parameter(NPL=15000)
 
148
      common/c_analysis/nwgt_analysis
 
149
      OPEN(UNIT=99,FILE='HERWIG.top',STATUS='UNKNOWN')
 
150
C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG 
 
151
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
 
152
      XNORM=1.D3/DFLOAT(NEVHEP)
 
153
      DO I=1,NPL              
 
154
        CALL MFINAL3(I)             
 
155
        CALL MCOPY(I,I+NPL)
 
156
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
157
        CALL MFINAL3(I+NPL)             
 
158
      ENDDO                          
 
159
C
 
160
      do kk=1,nwgt_analysis
 
161
      do icuts=1,ncuts
 
162
        l0=(kk-1)*ncuts*nplots+(icuts-1)*nplots
 
163
        do ilep1=1,nlepton
 
164
        do io=1,nsingle
 
165
          l=l0+nsingle*(ilep1-1)+io
 
166
          call multitop(NPL+l,NPL-1,3,2,'sing ',' ','LOG')
 
167
        enddo
 
168
        enddo
 
169
        ipair=0
 
170
        do ilep1=1,nlepton-1
 
171
        do ilep2=ilep1+1,nlepton
 
172
        ipair=ipair+1
 
173
        do io=1,ncorr
 
174
          l=l0+nlepton*nsingle+ncorr*(ipair-1)+io
 
175
          call multitop(NPL+l,NPL-1,3,2,'corr ',' ','LOG')
 
176
        enddo
 
177
        enddo
 
178
        enddo
 
179
      enddo
 
180
      enddo
 
181
      CLOSE(99)
 
182
      END
 
183
 
 
184
C----------------------------------------------------------------------
 
185
      SUBROUTINE HWANAL
 
186
C     USER''S ROUTINE TO ANALYSE DATA FROM EVENT
 
187
C----------------------------------------------------------------------
 
188
      INCLUDE 'HEPMC.INC'
 
189
      include 'reweight0.inc'
 
190
      DOUBLE PRECISION HWVDOT,PSUM(4)
 
191
      INTEGER ICHSUM,ICHINI,IHEP,IST,ID
 
192
      LOGICAL DIDSOF
 
193
      integer nwgt_analysis,max_weight
 
194
      common/c_analysis/nwgt_analysis
 
195
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
196
      double precision ww(max_weight),www(max_weight)
 
197
      common/cww/ww
 
198
c
 
199
      integer nsingle,ncorr,nlepton,nplots,ncuts
 
200
      common/cplots/nsingle,ncorr,nlepton,nplots,ncuts
 
201
      integer MAXELM,MAXELP,MAXMUM,MAXMUP
 
202
      common/cmaxlep/MAXELM,MAXELP,MAXMUM,MAXMUP
 
203
c
 
204
      integer i,j,ilep,ilep1,ilep2,io,NELP,NELM,NMUP,NMUM,
 
205
     # IPAIR,kk,L,L0,ICUTS
 
206
      real*8 GETPTV4,GETRAPIDITYV4,GETINVMV4,GETDELPHIV4,PPAIR(4),
 
207
     # PELM(4,25),PELP(4,25),PMUM(4,25),PMUP(4,25),PLEP(4,25),
 
208
     # obs(100)
 
209
c
 
210
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
211
c To be changed !!
 
212
      ww(1)=1d0
 
213
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
214
      IF(MOD(NEVHEP,10000).EQ.0)RETURN
 
215
      IF (WW(1).EQ.0D0) THEN
 
216
         WRITE(*,*)'WW(1) = 0. Stopping'
 
217
         STOP
 
218
      ENDIF
 
219
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT''S A POWER-SUPPRESSED
 
220
C EFFECT, SO THROW THE EVENT AWAY
 
221
      IF(SIGN(1.D0,PHEP(3,1)).EQ.SIGN(1.D0,PHEP(3,2)))THEN
 
222
        CALL HWWARN('HWANAL',111)
 
223
        GOTO 999
 
224
      ENDIF
 
225
      DO I=1,nwgt_analysis
 
226
         WWW(I)=EVWGT*ww(i)/ww(1)
 
227
      ENDDO
 
228
      ICHSUM=0
 
229
      NELP=0
 
230
      NELM=0
 
231
      NMUP=0
 
232
      NMUM=0
 
233
      DO 100 IHEP=1,NHEP
 
234
        IST=ISTHEP(IHEP)
 
235
        ID=IDHEP(IHEP)
 
236
c Find electrons
 
237
        IF(ID.EQ.11.AND.IST.EQ.1) THEN
 
238
           NELM=NELM+1
 
239
           DO I=1,4
 
240
              PELM(I,NELM)=PHEP(I,IHEP)
 
241
           ENDDO
 
242
c Find positrons
 
243
        ELSEIF(ID.EQ.-11.AND.IST.EQ.1) THEN
 
244
           NELP=NELP+1
 
245
           DO I=1,4
 
246
              PELP(I,NELP)=PHEP(I,IHEP)
 
247
           ENDDO
 
248
c Find mu-
 
249
        ELSEIF(ID.EQ.13.AND.IST.EQ.1) THEN
 
250
           NMUM=NMUM+1
 
251
           DO I=1,4
 
252
              PMUM(I,NMUM)=PHEP(I,IHEP)
 
253
           ENDDO
 
254
c Find mu+
 
255
        ELSEIF(ID.EQ.-13.AND.IST.EQ.1) THEN
 
256
           NMUP=NMUP+1
 
257
           DO I=1,4
 
258
              PMUP(I,NMUP)=PHEP(I,IHEP)
 
259
           ENDDO
 
260
        ENDIF
 
261
  100 CONTINUE
 
262
      IF( NELM.LT.MAXELM .OR. NELP.LT.MAXELP .OR.
 
263
     #    NMUM.LT.MAXMUM .OR. NMUP.LT.MAXMUP )THEN
 
264
        CALL HWUEPR
 
265
        WRITE(*,*)NELM,NELP,NMUM,NMUP
 
266
        WRITE(*,*)MAXELM,MAXELP,MAXMUM,MAXMUP
 
267
        CALL HWWARN('HWANAL',500)
 
268
      ENDIF
 
269
c Keep the first MAX** leptons of the species **
 
270
      ILEP=0
 
271
      DO J=1,MAXELM
 
272
        ILEP=ILEP+1
 
273
        DO I=1,4
 
274
           PLEP(I,ILEP)=PELM(I,J)
 
275
        ENDDO
 
276
      ENDDO
 
277
      DO J=1,MAXELP
 
278
        ILEP=ILEP+1
 
279
        DO I=1,4
 
280
           PLEP(I,ILEP)=PELP(I,J)
 
281
        ENDDO
 
282
      ENDDO
 
283
      DO J=1,MAXMUM
 
284
        ILEP=ILEP+1
 
285
        DO I=1,4
 
286
           PLEP(I,ILEP)=PMUM(I,J)
 
287
        ENDDO
 
288
      ENDDO
 
289
      DO J=1,MAXMUP
 
290
        ILEP=ILEP+1
 
291
        DO I=1,4
 
292
           PLEP(I,ILEP)=PMUP(I,J)
 
293
        ENDDO
 
294
      ENDDO
 
295
      IF( ILEP.NE.nlepton )THEN
 
296
        CALL HWUEPR
 
297
        WRITE(*,*)ILEP,nlepton
 
298
        CALL HWWARN('HWANAL',501)
 
299
      ENDIF
 
300
c
 
301
      DO ILEP1=1,NLEPTON
 
302
        OBS(1)=GETPTV4(PLEP(1,ILEP1))
 
303
        OBS(2)=GETRAPIDITYV4(PLEP(1,ILEP1))
 
304
        DO KK=1,NWGT_ANALYSIS
 
305
        DO ICUTS=1,NCUTS
 
306
          L0=(KK-1)*NCUTS*NPLOTS+(ICUTS-1)*NPLOTS
 
307
          DO IO=1,NSINGLE
 
308
            L=L0+NSINGLE*(ILEP1-1)+IO
 
309
            CALL MFILL(L,OBS(IO),WWW(KK))
 
310
          ENDDO
 
311
        ENDDO
 
312
        ENDDO
 
313
      ENDDO
 
314
c
 
315
      IPAIR=0
 
316
      DO ILEP1=1,NLEPTON-1
 
317
        DO ILEP2=ILEP1+1,NLEPTON
 
318
          IPAIR=IPAIR+1
 
319
          DO I=1,4
 
320
             PPAIR(I)=PLEP(I,ILEP1)+PLEP(I,ILEP2)
 
321
          ENDDO
 
322
          OBS(1)=GETINVMV4(PPAIR)
 
323
          OBS(2)=GETPTV4(PPAIR)
 
324
          OBS(3)=GETDELPHIV4(PLEP(1,ILEP1),PLEP(1,ILEP2))
 
325
          OBS(4)=GETRAPIDITYV4(PLEP(1,ILEP1))-
 
326
     #           GETRAPIDITYV4(PLEP(1,ILEP2))
 
327
          DO KK=1,NWGT_ANALYSIS
 
328
          DO ICUTS=1,NCUTS
 
329
            L0=(KK-1)*NCUTS*NPLOTS+(ICUTS-1)*NPLOTS
 
330
            DO IO=1,NCORR
 
331
              L=L0+NLEPTON*NSINGLE+NCORR*(IPAIR-1)+IO
 
332
              CALL MFILL(L,OBS(IO),WWW(KK))
 
333
            ENDDO
 
334
          ENDDO
 
335
          ENDDO
 
336
        ENDDO
 
337
      ENDDO
 
338
 999  END
 
339
 
 
340
 
 
341
      FUNCTION RANDA(SEED)
 
342
*     -----------------
 
343
* Ref.: K. Park and K.W. Miller, Comm. of the ACM 31 (1988) p.1192
 
344
* Use seed = 1 as first value.
 
345
*
 
346
      IMPLICIT INTEGER(A-Z)
 
347
      DOUBLE PRECISION MINV,RANDA
 
348
      SAVE
 
349
      PARAMETER(M=2147483647,A=16807,Q=127773,R=2836)
 
350
      PARAMETER(MINV=0.46566128752458d-09)
 
351
      HI = SEED/Q
 
352
      LO = MOD(SEED,Q)
 
353
      SEED = A*LO - R*HI
 
354
      IF(SEED.LE.0) SEED = SEED + M
 
355
      RANDA = SEED*MINV
 
356
      END
 
357
 
 
358
 
 
359
 
 
360
 
 
361
      function getrapidity(en,pl)
 
362
      implicit none
 
363
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
364
      parameter (tiny=1.d-8)
 
365
c
 
366
      xplus=en+pl
 
367
      xminus=en-pl
 
368
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
369
        if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
 
370
          y=0.5d0*log( xplus/xminus )
 
371
        else
 
372
          y=sign(1.d0,pl)*1.d8
 
373
        endif
 
374
      else
 
375
        y=sign(1.d0,pl)*1.d8
 
376
      endif
 
377
      getrapidity=y
 
378
      return
 
379
      end
 
380
 
 
381
 
 
382
      function getpseudorap(en,ptx,pty,pl)
 
383
      implicit none
 
384
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
385
      parameter (tiny=1.d-5)
 
386
c
 
387
      pt=sqrt(ptx**2+pty**2)
 
388
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
389
        eta=sign(1.d0,pl)*1.d8
 
390
      else
 
391
        th=atan2(pt,pl)
 
392
        eta=-log(tan(th/2.d0))
 
393
      endif
 
394
      getpseudorap=eta
 
395
      return
 
396
      end
 
397
 
 
398
 
 
399
      function getinvm(en,ptx,pty,pl)
 
400
      implicit none
 
401
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
402
      parameter (tiny=1.d-5)
 
403
c
 
404
      tmp=en**2-ptx**2-pty**2-pl**2
 
405
      if(tmp.gt.0.d0)then
 
406
        tmp=sqrt(tmp)
 
407
      elseif(tmp.gt.-tiny)then
 
408
        tmp=0.d0
 
409
      else
 
410
        write(*,*)'Attempt to compute a negative mass'
 
411
        stop
 
412
      endif
 
413
      getinvm=tmp
 
414
      return
 
415
      end
 
416
 
 
417
 
 
418
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
419
      implicit none
 
420
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
421
      parameter (tiny=1.d-5)
 
422
c
 
423
      pt1=sqrt(ptx1**2+pty1**2)
 
424
      pt2=sqrt(ptx2**2+pty2**2)
 
425
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
426
        tmp=ptx1*ptx2+pty1*pty2
 
427
        tmp=tmp/(pt1*pt2)
 
428
        if(abs(tmp).gt.1.d0+tiny)then
 
429
          write(*,*)'Cosine larger than 1'
 
430
          stop
 
431
        elseif(abs(tmp).ge.1.d0)then
 
432
          tmp=sign(1.d0,tmp)
 
433
        endif
 
434
        tmp=acos(tmp)
 
435
      else
 
436
        tmp=1.d8
 
437
      endif
 
438
      getdelphi=tmp
 
439
      return
 
440
      end
 
441
 
 
442
 
 
443
      function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
 
444
      implicit none
 
445
      real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
 
446
     # getpseudorap,getdelphi
 
447
c
 
448
      deta=getpseudorap(en1,ptx1,pty1,pl1)-
 
449
     #     getpseudorap(en2,ptx2,pty2,pl2)
 
450
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
 
451
      getdr=sqrt(dphi**2+deta**2)
 
452
      return
 
453
      end
 
454
 
 
455
 
 
456
      function getdry(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
 
457
      implicit none
 
458
      real*8 getdry,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
 
459
     # getrapidity,getdelphi
 
460
c
 
461
      deta=getrapidity(en1,pl1)-
 
462
     #     getrapidity(en2,pl2)
 
463
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
 
464
      getdry=sqrt(dphi**2+deta**2)
 
465
      return
 
466
      end
 
467
 
 
468
 
 
469
      function getptv(p)
 
470
      implicit none
 
471
      real*8 getptv,p(5)
 
472
c
 
473
      getptv=sqrt(p(1)**2+p(2)**2)
 
474
      return
 
475
      end
 
476
 
 
477
 
 
478
      function getpseudorapv(p)
 
479
      implicit none
 
480
      real*8 getpseudorapv,p(5)
 
481
      real*8 getpseudorap
 
482
c
 
483
      getpseudorapv=getpseudorap(p(4),p(1),p(2),p(3))
 
484
      return
 
485
      end
 
486
 
 
487
 
 
488
      function getrapidityv(p)
 
489
      implicit none
 
490
      real*8 getrapidityv,p(5)
 
491
      real*8 getrapidity
 
492
c
 
493
      getrapidityv=getrapidity(p(4),p(3))
 
494
      return
 
495
      end
 
496
 
 
497
 
 
498
      function getdrv(p1,p2)
 
499
      implicit none
 
500
      real*8 getdrv,p1(5),p2(5)
 
501
      real*8 getdr
 
502
c
 
503
      getdrv=getdr(p1(4),p1(1),p1(2),p1(3),
 
504
     #             p2(4),p2(1),p2(2),p2(3))
 
505
      return
 
506
      end
 
507
 
 
508
 
 
509
      function getinvmv(p)
 
510
      implicit none
 
511
      real*8 getinvmv,p(5)
 
512
      real*8 getinvm
 
513
c
 
514
      getinvmv=getinvm(p(4),p(1),p(2),p(3))
 
515
      return
 
516
      end
 
517
 
 
518
 
 
519
      function getdelphiv(p1,p2)
 
520
      implicit none
 
521
      real*8 getdelphiv,p1(5),p2(5)
 
522
      real*8 getdelphi
 
523
c
 
524
      getdelphiv=getdelphi(p1(1),p1(2),
 
525
     #                     p2(1),p2(2))
 
526
      return
 
527
      end
 
528
 
 
529
 
 
530
      function getptv4(p)
 
531
      implicit none
 
532
      real*8 getptv4,p(4)
 
533
c
 
534
      getptv4=sqrt(p(1)**2+p(2)**2)
 
535
      return
 
536
      end
 
537
 
 
538
 
 
539
      function getpseudorapv4(p)
 
540
      implicit none
 
541
      real*8 getpseudorapv4,p(4)
 
542
      real*8 getpseudorap
 
543
c
 
544
      getpseudorapv4=getpseudorap(p(4),p(1),p(2),p(3))
 
545
      return
 
546
      end
 
547
 
 
548
 
 
549
      function getrapidityv4(p)
 
550
      implicit none
 
551
      real*8 getrapidityv4,p(4)
 
552
      real*8 getrapidity
 
553
c
 
554
      getrapidityv4=getrapidity(p(4),p(3))
 
555
      return
 
556
      end
 
557
 
 
558
 
 
559
      function getdrv4(p1,p2)
 
560
      implicit none
 
561
      real*8 getdrv4,p1(4),p2(4)
 
562
      real*8 getdr
 
563
c
 
564
      getdrv4=getdr(p1(4),p1(1),p1(2),p1(3),
 
565
     #              p2(4),p2(1),p2(2),p2(3))
 
566
      return
 
567
      end
 
568
 
 
569
 
 
570
      function getinvmv4(p)
 
571
      implicit none
 
572
      real*8 getinvmv4,p(4)
 
573
      real*8 getinvm
 
574
c
 
575
      getinvmv4=getinvm(p(4),p(1),p(2),p(3))
 
576
      return
 
577
      end
 
578
 
 
579
 
 
580
      function getdelphiv4(p1,p2)
 
581
      implicit none
 
582
      real*8 getdelphiv4,p1(4),p2(4)
 
583
      real*8 getdelphi
 
584
c
 
585
      getdelphiv4=getdelphi(p1(1),p1(2),
 
586
     #                      p2(1),p2(2))
 
587
      return
 
588
      end
 
589
 
 
590
 
 
591
      function getcosv4(q1,q2)
 
592
      implicit none
 
593
      real*8 getcosv4,q1(4),q2(4)
 
594
      real*8 xnorm1,xnorm2,tmp
 
595
c
 
596
      if(q1(4).lt.0.d0.or.q2(4).lt.0.d0)then
 
597
        getcosv4=-1.d10
 
598
        return
 
599
      endif
 
600
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
 
601
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
 
602
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
 
603
        tmp=-1.d10
 
604
      else
 
605
        tmp=q1(1)*q2(1)+q1(2)*q2(2)+q1(3)*q2(3)
 
606
        tmp=tmp/(xnorm1*xnorm2)
 
607
        if(abs(tmp).gt.1.d0.and.abs(tmp).le.1.001d0)then
 
608
          tmp=sign(1.d0,tmp)
 
609
        elseif(abs(tmp).gt.1.001d0)then
 
610
          write(*,*)'Error in getcosv4',tmp
 
611
          stop
 
612
        endif
 
613
      endif
 
614
      getcosv4=tmp
 
615
      return
 
616
      end
 
617
 
 
618
 
 
619
 
 
620
      function getmod(p)
 
621
      implicit none
 
622
      double precision p(4),getmod
 
623
 
 
624
      getmod=sqrt(p(1)**2+p(2)**2+p(3)**2)
 
625
 
 
626
      return
 
627
      end
 
628
 
 
629
 
 
630
 
 
631
      subroutine getperpenv4(q1,q2,qperp)
 
632
c Normal to the plane defined by \vec{q1},\vec{q2}
 
633
      implicit none
 
634
      real*8 q1(4),q2(4),qperp(4)
 
635
      real*8 xnorm1,xnorm2
 
636
      integer i
 
637
c
 
638
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
 
639
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
 
640
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
 
641
        do i=1,4
 
642
          qperp(i)=-1.d10
 
643
        enddo
 
644
      else
 
645
        qperp(1)=q1(2)*q2(3)-q1(3)*q2(2)
 
646
        qperp(2)=q1(3)*q2(1)-q1(1)*q2(3)
 
647
        qperp(3)=q1(1)*q2(2)-q1(2)*q2(1)
 
648
        do i=1,3
 
649
          qperp(i)=qperp(i)/(xnorm1*xnorm2)
 
650
        enddo
 
651
        qperp(4)=1.d0
 
652
      endif
 
653
      return
 
654
      end
 
655
 
 
656
 
 
657
 
 
658
 
 
659
 
 
660
      subroutine boostwdir2(chybst,shybst,chybstmo,xd,xin,xout)
 
661
c chybstmo = chybst-1; if it can be computed analytically it improves
 
662
c the numerical accuracy
 
663
      implicit none
 
664
      real*8 chybst,shybst,chybstmo,xd(1:3),xin(0:3),xout(0:3)
 
665
      real*8 tmp,en,pz
 
666
      integer i
 
667
c
 
668
      if(abs(xd(1)**2+xd(2)**2+xd(3)**2-1).gt.1.d-6)then
 
669
        write(*,*)'Error #1 in boostwdir2',xd
 
670
        stop
 
671
      endif
 
672
c
 
673
      en=xin(0)
 
674
      pz=xin(1)*xd(1)+xin(2)*xd(2)+xin(3)*xd(3)
 
675
      xout(0)=en*chybst-pz*shybst
 
676
      do i=1,3
 
677
        xout(i)=xin(i)+xd(i)*(pz*chybstmo-en*shybst)
 
678
      enddo
 
679
c
 
680
      return
 
681
      end
 
682
 
 
683
 
 
684
 
 
685
 
 
686
      subroutine boostwdir3(chybst,shybst,chybstmo,xd,xxin,xxout)
 
687
      implicit none
 
688
      real*8 chybst,shybst,chybstmo,xd(1:3),xxin(4),xxout(4)
 
689
      real*8 xin(0:3),xout(0:3)
 
690
      integer i
 
691
c
 
692
      do i=1,4
 
693
         xin(mod(i,4))=xxin(i)
 
694
      enddo
 
695
      call boostwdir2(chybst,shybst,chybstmo,xd,xin,xout)
 
696
      do i=1,4
 
697
         xxout(i)=xout(mod(i,4))
 
698
      enddo
 
699
c
 
700
      return
 
701
      end
 
702
 
 
703
 
 
704
 
 
705
 
 
706
      subroutine getwedge(p1,p2,pout)
 
707
      implicit none
 
708
      real*8 p1(4),p2(4),pout(4)
 
709
 
 
710
      pout(1)=p1(2)*p2(3)-p1(3)*p2(2)
 
711
      pout(2)=p1(3)*p2(1)-p1(1)*p2(3)
 
712
      pout(3)=p1(1)*p2(2)-p1(2)*p2(1)
 
713
      pout(4)=0d0
 
714
 
 
715
      return
 
716
      end
 
717
 
 
718
C-----------------------------------------------------------------------
 
719
      SUBROUTINE HWWARN(SUBRTN,ICODE)
 
720
C-----------------------------------------------------------------------
 
721
C     DEALS WITH ERRORS DURING EXECUTION
 
722
C     SUBRTN = NAME OF CALLING SUBROUTINE
 
723
C     ICODE  = ERROR CODE:    - -1 NONFATAL, KILL EVENT & PRINT NOTHING
 
724
C                            0- 49 NONFATAL, PRINT WARNING & CONTINUE
 
725
C                           50- 99 NONFATAL, PRINT WARNING & JUMP
 
726
C                          100-199 NONFATAL, DUMP & KILL EVENT
 
727
C                          200-299    FATAL, TERMINATE RUN
 
728
C                          300-399    FATAL, DUMP EVENT & TERMINATE RUN
 
729
C                          400-499    FATAL, DUMP EVENT & STOP DEAD
 
730
C                          500-       FATAL, STOP DEAD WITH NO DUMP
 
731
C-----------------------------------------------------------------------
 
732
      INCLUDE 'HEPMC.INC'
 
733
      INTEGER ICODE,NRN,IERROR
 
734
      CHARACTER*6 SUBRTN
 
735
      IF (ICODE.GE.0) WRITE (6,10) SUBRTN,ICODE
 
736
   10 FORMAT(/' HWWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4)
 
737
      IF (ICODE.LT.0) THEN
 
738
         IERROR=ICODE
 
739
         RETURN
 
740
      ELSEIF (ICODE.LT.100) THEN
 
741
         WRITE (6,20) NEVHEP,NRN,EVWGT
 
742
   20    FORMAT(' EVENT',I8,':   SEEDS =',I11,' &',I11,
 
743
     &'  WEIGHT =',E11.4/' EVENT SURVIVES. EXECUTION CONTINUES')
 
744
         IF (ICODE.GT.49) RETURN
 
745
      ELSEIF (ICODE.LT.200) THEN
 
746
         WRITE (6,30) NEVHEP,NRN,EVWGT
 
747
   30    FORMAT(' EVENT',I8,':   SEEDS =',I11,' &',I11,
 
748
     &'  WEIGHT =',E11.4/' EVENT KILLED.   EXECUTION CONTINUES')
 
749
         IERROR=ICODE
 
750
         RETURN
 
751
      ELSEIF (ICODE.LT.300) THEN
 
752
         WRITE (6,40)
 
753
   40    FORMAT(' EVENT SURVIVES.  RUN ENDS GRACEFULLY')
 
754
c$$$         CALL HWEFIN
 
755
c$$$         CALL HWAEND
 
756
         STOP
 
757
      ELSEIF (ICODE.LT.400) THEN
 
758
         WRITE (6,50)
 
759
   50    FORMAT(' EVENT KILLED: DUMP FOLLOWS.  RUN ENDS GRACEFULLY')
 
760
         IERROR=ICODE
 
761
c$$$         CALL HWUEPR
 
762
c$$$         CALL HWUBPR
 
763
c$$$         CALL HWEFIN
 
764
c$$$         CALL HWAEND
 
765
         STOP
 
766
      ELSEIF (ICODE.LT.500) THEN
 
767
         WRITE (6,60)
 
768
   60    FORMAT(' EVENT KILLED: DUMP FOLLOWS.  RUN STOPS DEAD')
 
769
         IERROR=ICODE
 
770
c$$$         CALL HWUEPR
 
771
c$$$         CALL HWUBPR
 
772
         STOP
 
773
      ELSE
 
774
         WRITE (6,70)
 
775
   70    FORMAT(' RUN CANNOT CONTINUE')
 
776
         STOP
 
777
      ENDIF
 
778
      END
 
779
 
 
780
      subroutine HWUEPR
 
781
      INCLUDE 'HEPMC.INC'
 
782
      integer ip,i
 
783
      PRINT *,' EVENT ',NEVHEP
 
784
      DO IP=1,NHEP
 
785
         PRINT '(I4,I8,I4,4I4,1P,5D11.3)',IP,IDHEP(IP),ISTHEP(IP),
 
786
     &        JMOHEP(1,IP),JMOHEP(2,IP),JDAHEP(1,IP),JDAHEP(2,IP),
 
787
     &        (PHEP(I,IP),I=1,5)
 
788
      ENDDO
 
789
      return
 
790
      end