← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2-alpha0

« back to all changes in this revision

Viewing changes to Template/Source/kin_functions.f

Added Template and HELAS into bzr

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/Source/kin_functions.f
 
1
c************************************************************************
 
2
c  THIS FILE CONTAINS THE DEFINITIONS OF USEFUL FUNCTIONS OF MOMENTA:
 
3
c  
 
4
c  DOT(p1,p2)         : 4-Vector Dot product
 
5
c  R2(p1,p2)          : distance in eta,phi between two particles
 
6
c  SumDot(P1,P2,dsign): invariant mass of 2 particles
 
7
c  rap(p)             : rapidity of particle in the lab frame (p in CM frame)
 
8
C  RAP2(P)            : rapidity of particle in the lab frame (p in lab frame)
 
9
c  DELTA_PHI(P1, P2)  : separation in phi of two particles 
 
10
c  ET(p)              : transverse energy of particle
 
11
c  PT(p)              : transverse momentum of particle
 
12
c  DJ(p1,p2)          : y*S (Durham) value for two partons
 
13
c  DJB(p1,p2)         : mT^2=m^2+pT^2 for one particle
 
14
c  PYJB(p1,p2)        : The Pythia ISR pT^2=(1-x)*Q^2
 
15
c  DJ2(p1,p2)         : scalar product squared
 
16
c  threedot(p1,p2)    : 3-vector Dot product (accept 4 vector in entry)
 
17
c  rho                : |p| in lab frame
 
18
c  eta                : pseudo-rapidity
 
19
c  phi                : phi
 
20
c  four_momentum      : (theta,phi,rho,mass)-> 4 vector
 
21
c  four_momentum_set2 : (pt,eta,phi,mass--> 4 vector
 
22
c
 
23
c************************************************************************
 
24
 
 
25
      DOUBLE PRECISION FUNCTION R2(P1,P2)
 
26
c************************************************************************
 
27
c     Distance in eta,phi between two particles.
 
28
c************************************************************************
 
29
      IMPLICIT NONE
 
30
c
 
31
c     Arguments
 
32
c
 
33
      double precision p1(0:3),p2(0:3)
 
34
c
 
35
c     External
 
36
c
 
37
      double precision rap,DELTA_PHI
 
38
      external rap,delta_phi
 
39
c-----
 
40
c  Begin Code
 
41
c-----
 
42
      R2 = (DELTA_PHI(P1,P2))**2+(rap(p1)-rap(p2))**2
 
43
      RETURN
 
44
      END
 
45
 
 
46
      DOUBLE PRECISION FUNCTION SumDot(P1,P2,dsign)
 
47
c************************************************************************
 
48
c     Invarient mass of 2 particles
 
49
c************************************************************************
 
50
      IMPLICIT NONE
 
51
c
 
52
c     Arguments
 
53
c
 
54
      double precision p1(0:3),p2(0:3),dsign
 
55
c
 
56
c     Local
 
57
c      
 
58
      integer i
 
59
      double precision ptot(0:3)
 
60
c
 
61
c     External
 
62
c
 
63
      double precision dot
 
64
      external dot
 
65
c-----
 
66
c  Begin Code
 
67
c-----
 
68
 
 
69
      do i=0,3
 
70
         ptot(i)=p1(i)+dsign*p2(i)
 
71
      enddo
 
72
      SumDot = dot(ptot,ptot)
 
73
      RETURN
 
74
      END
 
75
 
 
76
      DOUBLE PRECISION  FUNCTION rap(p)
 
77
c************************************************************************
 
78
c     Returns rapidity of particle in the lab frame
 
79
c************************************************************************
 
80
      IMPLICIT NONE
 
81
c
 
82
c     Arguments
 
83
c
 
84
      double precision  p(0:3)
 
85
c
 
86
c     Local
 
87
c
 
88
      double precision pm
 
89
c
 
90
c     Global
 
91
c
 
92
      include 'run.inc'
 
93
c-----
 
94
c  Begin Code
 
95
c-----
 
96
c      pm=dsqrt(p(1)**2+p(2)**2+p(3)**2)
 
97
      pm = p(0)
 
98
      rap = .5d0*dlog((pm+p(3))/(pm-p(3)))+
 
99
     $     .5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
 
100
      end
 
101
      DOUBLE PRECISION  FUNCTION rap2(p)
 
102
c************************************************************************
 
103
c     Returns rapidity of particle in the lab frame
 
104
c************************************************************************
 
105
      IMPLICIT NONE
 
106
c
 
107
c     Arguments
 
108
c
 
109
      double precision  p(0:3)
 
110
c
 
111
c     Local
 
112
c
 
113
      double precision pm
 
114
c
 
115
c     Global
 
116
c
 
117
      include 'run.inc'
 
118
c-----
 
119
c  Begin Code
 
120
c-----
 
121
c      pm=dsqrt(p(1)**2+p(2)**2+p(3)**2)
 
122
      pm = p(0)
 
123
      rap2 = .5d0*dlog((pm+p(3))/(pm-p(3)))
 
124
      end
 
125
 
 
126
      DOUBLE PRECISION FUNCTION DELTA_PHI(P1, P2)
 
127
c************************************************************************
 
128
c     Returns separation in phi of two particles p1,p2
 
129
c************************************************************************
 
130
      IMPLICIT NONE
 
131
c
 
132
c     Arguments
 
133
c
 
134
      double precision p1(0:3),p2(0:3)
 
135
c
 
136
c     Local
 
137
c
 
138
      REAL*8 DENOM, TEMP
 
139
c-----
 
140
c  Begin Code
 
141
c-----
 
142
      DENOM = SQRT(P1(1)**2 + P1(2)**2) * SQRT(P2(1)**2 + P2(2)**2)
 
143
      TEMP = MAX(-0.99999999D0, (P1(1)*P2(1) + P1(2)*P2(2)) / DENOM)
 
144
      TEMP = MIN( 0.99999999D0, TEMP)
 
145
      DELTA_PHI = ACOS(TEMP)
 
146
      END
 
147
 
 
148
 
 
149
 
 
150
      double precision function et(p)
 
151
c************************************************************************
 
152
c     Returns transverse energy of particle
 
153
c************************************************************************
 
154
      IMPLICIT NONE
 
155
c
 
156
c     Arguments
 
157
c
 
158
      double precision p(0:3)
 
159
c
 
160
c     Local
 
161
c
 
162
      double precision pt
 
163
c-----
 
164
c  Begin Code
 
165
c-----
 
166
      pt = dsqrt(p(1)**2+p(2)**2)
 
167
      if (pt .gt. 0d0) then
 
168
         et = p(0)*pt/dsqrt(pt**2+p(3)**2)
 
169
      else
 
170
         et = 0d0
 
171
      endif
 
172
      end
 
173
 
 
174
      double precision function pt(p)
 
175
c************************************************************************
 
176
c     Returns transverse momentum of particle
 
177
c************************************************************************
 
178
      IMPLICIT NONE
 
179
c
 
180
c     Arguments
 
181
c
 
182
      double precision p(0:3)
 
183
c-----
 
184
c  Begin Code
 
185
c-----
 
186
 
 
187
      pt = dsqrt(p(1)**2+p(2)**2)
 
188
 
 
189
      return
 
190
      end
 
191
 
 
192
      double precision function DJ(p1,p2)
 
193
c***************************************************************************
 
194
c     Uses Durham algorythm to calculate the y value for two partons
 
195
c     If collision type is hh, hadronic jet measure is used
 
196
c       y_{ij} = 2min[p_{i,\perp}^2,p_{j,\perp}^2]/S
 
197
c                  (cosh(\eta_i-\eta_j)-cos(\phi_1-\phi_2))
 
198
c***************************************************************************
 
199
      implicit none
 
200
c
 
201
c     Arguments
 
202
c
 
203
      double precision p1(0:4),p2(0:4) ! 4 is mass**2
 
204
c
 
205
c     Global
 
206
c
 
207
      double precision D
 
208
      common/to_dj/D
 
209
c
 
210
c     Local
 
211
c
 
212
 
 
213
      include 'run.inc'
 
214
 
 
215
      double precision pt1,pt2,ptm1,ptm2,eta1,eta2,phi1,phi2,p1a,p2a,costh,sumdot
 
216
      integer j
 
217
c-----
 
218
c  Begin Code
 
219
c-----
 
220
      if ((lpp(1).eq.0).and.(lpp(2).eq.0)) then
 
221
         p1a = dsqrt(p1(1)**2+p1(2)**2+p1(3)**2)
 
222
         p2a = dsqrt(p2(1)**2+p2(2)**2+p2(3)**2)
 
223
         if (p1a*p2a .ne. 0d0) then
 
224
            costh = (p1(1)*p2(1)+p1(2)*p2(2)+p1(3)*p2(3))/(p1a*p2a)
 
225
            dj = 2d0*min(p1(0)**2,p2(0)**2)*(1d0-costh) !Durham
 
226
c            dj = 2d0*p1(0)*p2(0)*(1d0-costh)    !JADE
 
227
         else
 
228
            print*,'Warning 0 momentum in Durham algorythm'
 
229
            write(*,'(4e15.5)') (p1(j),j=0,3)
 
230
            write(*,'(4e15.5)') (p2(j),j=0,3)
 
231
            dj = 0d0
 
232
         endif
 
233
      else
 
234
        pt1 = p1(1)**2+p1(2)**2
 
235
        pt2 = p2(1)**2+p2(2)**2
 
236
        p1a = dsqrt(pt1+p1(3)**2)
 
237
        p2a = dsqrt(pt2+p2(3)**2)
 
238
        eta1 = 0.5d0*log((p1a+p1(3))/(p1a-p1(3)))
 
239
        eta2 = 0.5d0*log((p2a+p2(3))/(p2a-p2(3)))
 
240
c        ptm1 = max((p1(0)-p1(3))*(p1(0)+p1(3)),0d0)
 
241
c        ptm2 = max((p2(0)-p2(3))*(p2(0)+p2(3)),0d0)
 
242
        dj = max(p1(4),p2(4))+min(pt1,pt2)*2d0*(cosh(eta1-eta2)-
 
243
     &     (p1(1)*p2(1)+p1(2)*p2(2))/dsqrt(pt1*pt2))/D**2
 
244
 
 
245
c        write(*,'(a,5e16.4)')'Mom(1): ',(p1(j),j=1,3),p1(0),p1(4)
 
246
c        write(*,'(a,5e16.4)')'Mom(2): ',(p2(j),j=1,3),p2(0),p2(4)
 
247
c       print *,'pT1: ',sqrt(pt1),' pT2: ',sqrt(pt2)
 
248
c       print *,'deltaR: ',sqrt(2d0*(cosh(eta1-eta2)-
 
249
c     &     (p1(1)*p2(1)+p1(2)*p2(2))/dsqrt(pt1*pt2))/D**2),
 
250
c     $      ' m: ',sqrt(SumDot(p1,p2,1d0))
 
251
c     write(*,*) 'p1  = ',p1(0),',',p1(1),',',p1(2),',',p1(3)
 
252
c     write(*,*) 'pm1 = ',pm1,', p1a = ',p1a,'eta1 = ',eta1
 
253
c     write(*,*) 'p2  = ',p2(0),',',p2(1),',',p2(2),',',p2(3)
 
254
c     write(*,*) 'pm2 = ',pm2,', p2a = ',p2a,'eta2 = ',eta2
 
255
c     write(*,*) 'dj = ',dj
 
256
      endif
 
257
      end
 
258
      
 
259
      double precision function PYDJ(p1,p2)
 
260
c***************************************************************************
 
261
c     Uses Durham algorythm to calculate the y value for two partons
 
262
c     If collision type is hh, hadronic jet measure is used
 
263
c       y_{ij} = 2min[p_{i,\perp}^2,p_{j,\perp}^2]/S
 
264
c                  (cosh(\eta_i-\eta_j)-cos(\phi_1-\phi_2))
 
265
c***************************************************************************
 
266
      implicit none
 
267
c
 
268
c     Arguments
 
269
c
 
270
      double precision p1(0:4),p2(0:4) ! 4 is mass**2
 
271
c
 
272
c     Global
 
273
c
 
274
      double precision D
 
275
      common/to_dj/D
 
276
c
 
277
c     Local
 
278
c
 
279
 
 
280
      double precision SumDot
 
281
      external SumDot
 
282
c-----
 
283
c  Begin Code
 
284
c-----
 
285
 
 
286
      pydj = p1(0)*p2(0)/(p1(0)+p2(0))**2*SumDot(p1,p2,1d0)
 
287
 
 
288
      end
 
289
      
 
290
      double precision function DJ1(p1,p2)
 
291
c***************************************************************************
 
292
c     Uses single-sided Durham algorythm to calculate the y value for 
 
293
c     parton radiated off non-parton
 
294
c     If collision type is hh, hadronic jet measure is used
 
295
c       y_{ij} = 2min[p_{i,\perp}^2,p_{j,\perp}^2]/S
 
296
c                  (cosh(\eta_i-\eta_j)-cos(\phi_1-\phi_2))
 
297
c***************************************************************************
 
298
      implicit none
 
299
c
 
300
c     Arguments
 
301
c
 
302
      double precision p1(0:3),p2(0:3)
 
303
c
 
304
c     Local
 
305
c
 
306
 
 
307
      include 'run.inc'
 
308
 
 
309
      double precision pt1,pt2,ptm1,eta1,eta2,phi1,phi2,p1a,p2a,costh
 
310
      integer j
 
311
c-----
 
312
c  Begin Code
 
313
c-----
 
314
      if ((lpp(1).eq.0).and.(lpp(2).eq.0)) then
 
315
      p1a = dsqrt(p1(1)**2+p1(2)**2+p1(3)**2)
 
316
      p2a = dsqrt(p2(1)**2+p2(2)**2+p2(3)**2)
 
317
      if (p1a*p2a .ne. 0d0) then
 
318
         costh = (p1(1)*p2(1)+p1(2)*p2(2)+p1(3)*p2(3))/(p1a*p2a)
 
319
         dj1 = 2d0*p1(0)**2*(1d0-costh)   !Durham
 
320
c         dj = 2d0*p1(0)*p2(0)*(1d0-costh)    !JADE
 
321
      else
 
322
         print*,'Warning 0 momentum in Durham algorythm'
 
323
         write(*,'(4e15.5)') (p1(j),j=0,3)
 
324
         write(*,'(4e15.5)') (p2(j),j=0,3)
 
325
         dj1 = 0d0
 
326
      endif
 
327
      else
 
328
        pt1 = p1(1)**2+p1(2)**2
 
329
        pt2 = p2(1)**2+p2(2)**2
 
330
        p1a = dsqrt(pt1+p1(3)**2)
 
331
        p2a = dsqrt(pt2+p2(3)**2)
 
332
        eta1 = 0.5d0*log((p1a+p1(3))/(p1a-p1(3)))
 
333
        eta2 = 0.5d0*log((p2a+p2(3))/(p2a-p2(3)))
 
334
        ptm1 = max((p1(0)-p1(3))*(p1(0)+p1(3)),0d0)
 
335
        dj1 = 2d0*ptm1*(cosh(eta1-eta2)-
 
336
     &     (p1(1)*p2(1)+p1(2)*p2(2))/dsqrt(pt1*pt2))
 
337
c     write(*,*) 'p1  = ',p1(0),',',p1(1),',',p1(2),',',p1(3)
 
338
c     write(*,*) 'pm1 = ',pm1,', p1a = ',p1a,'eta1 = ',eta1
 
339
c     write(*,*) 'p2  = ',p2(0),',',p2(1),',',p2(2),',',p2(3)
 
340
c     write(*,*) 'pm2 = ',pm2,', p2a = ',p2a,'eta2 = ',eta2
 
341
c     write(*,*) 'dj = ',dj
 
342
      endif
 
343
      end
 
344
      
 
345
      double precision function DJB(p1)
 
346
c***************************************************************************
 
347
c     Uses kt algorythm to calculate the y value for one parton
 
348
c       y_i    = p_{i,\perp}^2/S
 
349
c***************************************************************************
 
350
      implicit none
 
351
c
 
352
c     Arguments
 
353
c
 
354
      double precision p1(0:4)  ! 4 is mass**2
 
355
c
 
356
c     Local
 
357
c
 
358
      double precision pm1
 
359
      include 'run.inc'
 
360
 
 
361
c-----
 
362
c  Begin Code
 
363
c-----
 
364
c      pm1=max(p1(0)**2-p1(1)**2-p1(2)**2-p1(3)**2,0d0)
 
365
      if ((lpp(1).eq.0).and.(lpp(2).eq.0)) then
 
366
c        write(*,*) 'kin_functions.f: Error. No jet measure w.r.t. beam.'
 
367
c        djb = 0d0
 
368
         djb=max(p1(0),0d0)
 
369
      else
 
370
        djb = (p1(0)-p1(3))*(p1(0)+p1(3)) ! p1(1)**2+p1(2)**2+pm1
 
371
c        djb = p1(1)**2+p1(2)**2+p1(4)
 
372
      endif
 
373
      end
 
374
 
 
375
      double precision function PYJB(p2,p1,ppart,z)
 
376
c***************************************************************************
 
377
c     Calculate the Pythia ISR evolution pT2
 
378
c       pTE2 = (1-z)(Q2+m2), Q2=-(p1-p2)**2, z=sred/sprev
 
379
c     Note! p1 and p2 must have mass**2 component!
 
380
c***************************************************************************
 
381
      implicit none
 
382
c
 
383
c     Arguments
 
384
c
 
385
      double precision p1(0:4),p2(0:4),ppart(0:3),z
 
386
c
 
387
c     Local
 
388
c
 
389
      double precision sred,sprev,Q2,pstar(0:3),pm2
 
390
      integer i
 
391
      double precision dot,SumDot,PT
 
392
 
 
393
c-----
 
394
c  Begin Code
 
395
c-----
 
396
      pm2=0
 
397
 
 
398
      if(p1(4).gt.0.or.p2(4).gt.0.and..not.
 
399
     $   (p1(4).gt.0.and.p2(4).gt.0)) pm2=max(p1(4),p2(4))
 
400
      do i=0,3
 
401
        pstar(i)=p1(i)-p2(i)
 
402
      enddo
 
403
      Q2=-dot(pstar,pstar)+pm2
 
404
      if(Q2.lt.0)then
 
405
c        print *,'Error in PYJB: Q2 = ',Q2
 
406
        PYJB=1d30
 
407
        return
 
408
      endif
 
409
      sprev=SumDot(p1,ppart,1d0)
 
410
      sred=SumDot(pstar,ppart,1d0)
 
411
 
 
412
      if(sred.lt.1d0)then
 
413
        PYJB=1d20
 
414
        z=0d0
 
415
        return
 
416
      endif
 
417
 
 
418
      z=sred/sprev
 
419
      if(z.gt.1.or.z.lt.0)then
 
420
        print *,'Error in PYJB: z = ',z,', sprev = ',sprev,
 
421
     $     ', sred = ',sred,', Q2 = ',Q2
 
422
        stop
 
423
      endif
 
424
      PYJB=(1d0-z)*Q2
 
425
      end
 
426
 
 
427
      double precision function zclus(p2,p1,ppart)
 
428
c***************************************************************************
 
429
c     Calculate the Pythia ISR evolution pT2
 
430
c     z=sred/sprev
 
431
c***************************************************************************
 
432
      implicit none
 
433
c
 
434
c     Arguments
 
435
c
 
436
      double precision p1(0:3),p2(0:3),ppart(0:3)
 
437
c
 
438
c     Local
 
439
c
 
440
      double precision sred,sprev,pstar(0:3)
 
441
      integer i, nerr
 
442
      data nerr/0/
 
443
      double precision dot,SumDot
 
444
 
 
445
c-----
 
446
c  Begin Code
 
447
c-----
 
448
      do i=0,3
 
449
        pstar(i)=p1(i)-p2(i)
 
450
      enddo
 
451
      sprev=SumDot(p1,ppart,1d0)
 
452
      sred=SumDot(pstar,ppart,1d0)
 
453
 
 
454
      if(sred.lt.1d0)then
 
455
        zclus=0d0
 
456
        return
 
457
      endif
 
458
 
 
459
      zclus=sred/sprev
 
460
      if((zclus.gt.1.or.zclus.lt.0).and.nerr.le.10)then
 
461
        print *,'Error in zclus: zclus = ',zclus,', sprev = ',sprev,
 
462
     $     ', sred = ',sred
 
463
        nerr=nerr+1
 
464
        if(nerr.eq.10)
 
465
     $       print *,'No more zclus errors will be printed'
 
466
      endif
 
467
 
 
468
      return
 
469
      end
 
470
 
 
471
      double precision function DJ2(p1,p2)
 
472
c***************************************************************************
 
473
c     Uses Lorentz
 
474
c***************************************************************************
 
475
      implicit none
 
476
c
 
477
c     Arguments
 
478
c
 
479
      double precision p1(0:3),p2(0:3)
 
480
c
 
481
c     Local
 
482
c
 
483
      integer j
 
484
c
 
485
c     External
 
486
c
 
487
      double precision dot
 
488
c-----
 
489
c  Begin Code
 
490
c-----
 
491
      dj2 = dot(p1,p1)+2d0*dot(p1,p2)+dot(p2,p2)
 
492
      return
 
493
      end
 
494
 
 
495
      subroutine switchmom(p1,p,ic,jc,nexternal)
 
496
c**************************************************************************
 
497
c     Changes stuff for crossings
 
498
c**************************************************************************
 
499
      implicit none
 
500
      integer nexternal
 
501
      integer jc(nexternal),ic(nexternal)
 
502
      real*8 p1(0:3,nexternal),p(0:3,nexternal)
 
503
      integer i,j
 
504
c-----
 
505
c Begin Code
 
506
c-----
 
507
      do i=1,nexternal
 
508
         do j=0,3
 
509
            p(j,ic(i))=p1(j,i)
 
510
         enddo
 
511
      enddo
 
512
      do i=1,nexternal
 
513
         jc(i)=1
 
514
      enddo
 
515
      jc(ic(1))=-1
 
516
      jc(ic(2))=-1
 
517
      end
 
518
 
 
519
      double precision function dot(p1,p2)
 
520
C****************************************************************************
 
521
C     4-Vector Dot product
 
522
C****************************************************************************
 
523
      implicit none
 
524
      double precision p1(0:3),p2(0:3)
 
525
      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
 
526
 
 
527
      if(dabs(dot).lt.1d-6)then ! solve numerical problem 
 
528
         dot=0d0
 
529
      endif
 
530
 
 
531
      end
 
532
C*****************************************************************************
 
533
C*****************************************************************************
 
534
C                      MadWeight function
 
535
C*****************************************************************************
 
536
C*****************************************************************************
 
537
 
 
538
      double precision function threedot(p1,p2)
 
539
C****************************************************************************
 
540
C     3-Vector  product
 
541
C****************************************************************************
 
542
      implicit none
 
543
      double precision p1(0:3),p2(0:3)
 
544
      threedot=p1(1)*p2(1)+p1(2)*p2(2)+p1(3)*p2(3)
 
545
 
 
546
      end
 
547
 
 
548
 
 
549
      double precision function rho(p1)
 
550
C****************************************************************************
 
551
C     computes rho(p)=dsqrt (p(1)**2+p(2)**2+p(3)**2)
 
552
C****************************************************************************
 
553
      implicit none
 
554
      double precision p1(0:3)
 
555
      double precision  threedot
 
556
      external  threedot
 
557
c
 
558
      rho=dsqrt(threedot(p1,p1))
 
559
 
 
560
      end
 
561
 
 
562
      double precision function theta(p)
 
563
c************************************************************************
 
564
c     Returns polar angle of particle
 
565
c************************************************************************
 
566
      IMPLICIT NONE
 
567
c
 
568
c     Arguments
 
569
c
 
570
      double precision p(0:3)
 
571
c-----
 
572
c  Begin Code
 
573
c-----
 
574
 
 
575
      theta=dacos(p(3)/dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
576
 
 
577
      return
 
578
      end
 
579
 
 
580
      double precision function eta(p)
 
581
c************************************************************************
 
582
c     Returns pseudo rapidity of particle
 
583
c************************************************************************
 
584
      IMPLICIT NONE
 
585
c
 
586
c     Arguments
 
587
c
 
588
      double precision p(0:3)
 
589
c
 
590
c     external
 
591
c
 
592
      double precision theta
 
593
      external theta
 
594
c-----
 
595
c  Begin Code
 
596
c-----
 
597
 
 
598
      eta=-dlog(dtan(theta(p)/2))
 
599
 
 
600
      return
 
601
      end
 
602
 
 
603
      subroutine four_momentum(theta,phi,rho,m,p)
 
604
c****************************************************************************
 
605
c     modif 3/07/07 : this subroutine defines 4-momentum from theta,phi,rho,m
 
606
c     with rho=px**2+py**2+pz**2
 
607
c****************************************************************************
 
608
c
 
609
c     argument
 
610
c
 
611
      double precision theta,phi,rho,m,p(0:3)
 
612
c
 
613
      P(1)=rho*dsin(theta)*dcos(phi)
 
614
      P(2)=rho*dsin(theta)*dsin(phi)
 
615
      P(3)=rho*dcos(theta)
 
616
      P(0)=dsqrt(rho**2+m**2)
 
617
 
 
618
      return
 
619
      end
 
620
      subroutine four_momentum_set2(eta,phi,PT,m,p)
 
621
c****************************************************************************
 
622
c     modif 16/11/06 : this subroutine defines 4-momentum from PT,eta,phi,m
 
623
c****************************************************************************
 
624
c
 
625
c     argument
 
626
c
 
627
      double precision PT,eta,phi,m,p(0:3)
 
628
c
 
629
c
 
630
c
 
631
      P(1)=PT*dcos(phi)
 
632
      P(2)=PT*dsin(phi)
 
633
      P(3)=PT*dsinh(eta)
 
634
      P(0)=dsqrt(p(1)**2+p(2)**2+p(3)**2+m**2)  
 
635
      return
 
636
      end
 
637
 
 
638
 
 
639
 
 
640
      DOUBLE PRECISION  FUNCTION phi(p)
 
641
c************************************************************************
 
642
c     MODIF 16/11/06 : this subroutine defines phi angle
 
643
c                      phi is defined from 0 to 2 pi
 
644
c************************************************************************
 
645
      IMPLICIT NONE
 
646
c
 
647
c     Arguments
 
648
c
 
649
      double precision  p(0:3)
 
650
c
 
651
c     Parameter
 
652
c
 
653
 
 
654
      double precision pi,zero
 
655
      parameter (pi=3.141592654d0,zero=0d0)
 
656
c-----
 
657
c  Begin Code
 
658
c-----
 
659
c 
 
660
      if(p(1).gt.zero) then
 
661
      phi=datan(p(2)/p(1))
 
662
      else if(p(1).lt.zero) then
 
663
      phi=datan(p(2)/p(1))+pi
 
664
      else if(p(2).GE.zero) then !remind that p(1)=0
 
665
      phi=pi/2d0
 
666
      else if(p(2).lt.zero) then !remind that p(1)=0
 
667
      phi=-pi/2d0
 
668
      endif
 
669
      if(phi.lt.zero) phi=phi+2*pi
 
670
      return
 
671
      end
 
672