1
c***********************************************************************
2
c Routines for generating momentum of decay products
3
c***********************************************************************
5
SUBROUTINE PHASESPACE(PD,PSWGT)
6
C************************************************************************
7
C WRAPPER FOR CALLING THE PHASE SPACE ROUTINES
8
C************************************************************************
14
REAL*8 PD(0:3,MaxDecPart)
30
CALL TWOMOM(PD(0,1),PD(0,2),PD(0,3),PSWGT)
32
CALL THREEMOM(PD(0,1),PD(0,2),PD(0,3),PD(0,4),PSWGT)
34
CALL FOURMOM(PD(0,1),PD(0,2),PD(0,3),PD(0,4),PD(0,5),PSWGT)
36
WRITE(*,*) 'NO PS AVAILABLE '
45
SUBROUTINE TWOMOM(P1,P2,P3,PSWGT)
46
C************************************************************************
47
C GENERIC TWO BODY PHACE SPACE FOR THE DECAYS P1->P2+P3
49
C GIVEN P1 MOMENTUM SETS UP THE MOMENTA P2,P3
50
C ALSO SETS UP THE APPROPRIATE PHASESPACE WEIGTH PSWGT
52
C X(1) = COSTHETA OF FIRST DECAY
53
C X(2) = PHI OF FIRST DECAY
55
C************************************************************************
61
REAL*8 P1(0:3),P2(0:3),P3(0:3)
66
REAL*8 COSTH,PHI,JACPOLE
85
C PICK VALUE OF THETA AND PHI
90
C DETERMINE JACOBIAN FOR THETA AND PHI
94
C CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
95
C OF THE FINAL STATE MASSLESS PARTONS IN THEIR CM FRAME
97
CALL MOM2CX(RS,M2,M3,COSTH,PHI,P2,P3) !DECAY P1
98
CALL BOOSTX(P2,P1,P2) !BOOST DECAY TO LAB
99
CALL BOOSTX(P3,P1,P3) !BOOST DECAY TO LAB
101
JAC = JAC /(2D0*PI)**2 !FOR THE (2 PI)^3N-4 OUT FRONT
103
C CALCULATE PHASE SPACE FACTOR DSIG/DCOS(THETA)
105
PSWGT = 1.D0/(2.D0*RS) !FLUX FACTOR
108
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI)
112
c write(*,*) 'from decay_mom',x(1),x(2),pswgt
118
SUBROUTINE THREEMOM(P1,P2,P3,P4,PSWGT)
119
C************************************************************************
120
C This routine generates the momenta for the decay products of
124
C when the W can be on-shell or not.
125
C It ALSO SETS UP THE APPROPRIATE PHASESPACE WEIGTH PSWGT
127
C X(1) = COSTHETA OF FIRST DECAY TOP -> W*+B
128
C X(2) = PHI OF FIRST DECAY TOP -> W*+B
130
C X(4) = COSTHETA OF DECAY W* -> E+VE
131
C X(5) = PHI OF DECAY W* -> E+VE
133
C************************************************************************
139
REAL*8 P1(0:3),P2(0:3),P3(0:3),P4(0:3),PW(0:3)
144
REAL*8 COSTH,PHI,JACPOLE
146
REAL*8 COSTH1,PHI1,S1,RS1
147
REAL*8 COSTH2,PHI2,S2,RS2
159
DATA FIRSTTIME /.TRUE./
171
IF(RS.GT.(M2+MV)) THEN ! W CAN BE ON SHELL
172
CALL TRANSPOLE(MV**2/S,MV*GV/S,X(3),S2,JACPOLE)
176
ELSE ! W is off-shell
177
S2=(S-M2*M2)*X(3)+M2*M2 ! generated uniformily
182
C check total energy conservation
184
IF(RS2.GE.(RS-M2).or.RS2.LT.(M3+M4)) THEN
189
C PICK VALUE OF THETA AND PHI FOR FIRST DECAY TOP -> W*+B
191
COSTH=-1.D0+2.D0*X(1)
194
C DETERMINE JACOBIAN FOR THETA AND PHI
198
C CHOOSE COS(THETA) AND PHI FOR DECAY W*-> E+VE
200
COSTH1 = -1 + 2.*X(4)
203
C JACOBIAN FACTOR FOR DECAY OF W*-> E+VE
207
C CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
208
C OF THE FINAL STATE MASSLESS PARTONS IN THEIR CM FRAME
210
CALL MOM2CX(RS,RS2,M2,COSTH ,PHI ,PW,P2) !DECAY 1
211
CALL MOM2CX(RS2,M3,M4,COSTH1,PHI1,P3,P4) !DECAY W*
212
CALL BOOSTX(P2,P1,P2) !BOOST DECAY TO LAB
213
CALL BOOSTX(PW,P1,PW) !BOOST DECAY TO LAB
214
CALL BOOSTX(P3,PW,P3) !BOOST DECAY TO LAB
215
CALL BOOSTX(P4,PW,P4) !BOOST DECAY TO LAB
217
JAC = JAC /(2D0*PI)**5 !FOR THE (2 PI)^3N-4 OUT FRONT
219
C CALCULATE PHASE SPACE FACTOR DSIG/DCOS(THETA)
222
PSWGT = 1.D0/(2.D0*RS) !FLUX FACTOR
225
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !T->W*+B
228
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !W*->E+VE
237
SUBROUTINE FOURMOM(P1,P2,P3,P4,P5,PSWGT)
238
C************************************************************************
239
C This routine generates the momenta for the decay products of
241
C 1 -> v v'* -> 2 3 4 5
243
C where the v'* can be on-shell or not.
245
C GIVEN X(8) AND THE P1 MOMENTUM SETS UP THE MOMENTA
246
C ALSO SETS UP THE APPROPRIATE PHASESPACE WEIGTH PSWGT
248
C X(1) = COSTHETA OF FIRST DECAY SH -> M_(2+3)+M_(4,5)
249
C X(2) = PHI OF FIRST DECAY SH -> M_(2+3)+M_(4,5)
252
C X(5) = COSTHETA OF DECAY M_(2,3) -> M2+M3
253
C X(6) = PHI OF FIRST DECAY M_(2,3) -> M2+M3
254
C X(7) = COSTHETA OF DECAY M_(4,5) -> M4+M5
255
C X(8) = PHI OF FIRST DECAY M_(4,5) -> M4+M5
257
C************************************************************************
263
REAL*8 P1(0:3),P2(0:3),P3(0:3),P4(0:3),P5(0:3)
269
REAL*8 RSH,SH,COSINIT,PHIINIT,COSTH,PHI
270
REAL*8 COSTH1,PHI1,S1,PTEMP1(0:3),X3MIN,X3MAX,RS1
271
REAL*8 COSTH2,PHI2,S2,PTEMP2(0:3),X4MIN,X4MAX,RS2
281
data iseed/1/ !no effect if already called
286
DATA FIRSTTIME /.TRUE./
294
C CMS ENERGY = (HIGGS MASS)^2
299
C PICK VALUE OF THETA AND PHI FOR FIRST DECAY SH->M3+S1
301
COSTH=-1.D0+2.D0*X(1)
304
C DETERMINE JACOBIAN FOR THETA AND PHI
308
C MASSES OF THE DECAY 1->V V
310
IF(RSH.GT.MV) THEN ! AT LEAST A W CAN BE ON SHELL
311
CALL TRANSPOLE(MV**2/SH,MV*GV/SH,X(3),S1,JACPOLE)
315
ELSE ! NO W'S CAN BE ON SHELL
316
S1=SH*X(3) ! generated uniformily
321
IF(RS1.GE.RSH) THEN ! ENERGY CONSERVATION
326
C DETERMINE INVARIENT MASS OF SYSTEM M4+M5
328
IF(RSH-RS1.GT.MV) THEN ! THE SECOND W CAN BE ON-SHELL
329
CALL TRANSPOLE(MV**2/SH,MV*GV/SH,X(4),S2,JACPOLE)
333
else ! the second w cannot be on-shell
334
c X4MIN = (M4+M5)**2/SH
335
X4MIN = 0d0 ! less efficient but symmetric in m2,m3 <-> m4,m5
336
X4MAX = (1D0-RS1/RSH)**2
337
S2 = ((X4MAX-X4MIN)*X(4)+X4MIN)*SH
339
JAC = JAC*SH*(X4MAX-X4MIN)
342
IF(RS2+RS1.GT.RSH) THEN ! ENERGY CONSERVATION
347
c switch order to ensure no bias
349
IF(xran1(iseed).lt.0.5e0)then
358
IF(RS1 .LT. (M2+M3)) THEN ! MINIMUM inv mass
363
IF(RS2 .LT. (M4+M5)) THEN ! MINIMUM inv mass
370
C CHOOSE COS(THETA) AND PHI FOR DECAY S1->M3+M4
372
COSTH1 = -1 + 2.*X(5)
375
C CHOOSE COS(THETA) AND PHI FOR DECAY S2->M4+M5
377
COSTH2 = -1 + 2.*X(7)
380
C JACOBIAN FACTOR FOR DECAY OF S1->M3+M4
384
C JACOBIAN FACTOR FOR DECAY OF S2->M4+M5
388
C CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
391
CALL MOM2CX(RSH,RS1,RS2,COSTH,PHI,PTEMP1,PTEMP2)
392
CALL MOM2CX(RS1,M2,M3,COSTH1,PHI1,P2,P3) !DECAY PTEMP1
393
CALL MOM2CX(RS2,M4,M5,COSTH2,PHI2,P4,P5) !DECAY PTEMP2
394
CALL BOOSTX(PTEMP1,P1,PTEMP1) !BOOST DECAY TO CM
395
CALL BOOSTX(PTEMP2,P1,PTEMP2) !BOOST DECAY TO CM
396
CALL BOOSTX(P2,PTEMP1,P2) !BOOST DECAY TO CM
397
CALL BOOSTX(P3,PTEMP1,P3) !BOOST DECAY TO CM
398
CALL BOOSTX(P4,PTEMP2,P4) !BOOST DECAY TO CM
399
CALL BOOSTX(P5,PTEMP2,P5) !BOOST DECAY TO CM
402
JAC = JAC /(2D0*PI)**8 !FOR THE (2 PI)^3N-4 OUT FRONT
405
C CALCULATE PHASE SPACE FACTOR DSIG/DCOS(THETA)
408
PSWGT = 1.D0/(2.D0*RSH) !FLUX FACTOR
411
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !S->S1+S2
414
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !S1->M2+M3
417
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !S2->M4+M5
427
REAL*8 FUNCTION LAMBDA(S,MA2,MB2)
429
C****************************************************************************
430
C THIS IS THE LAMBDA FUNCTION FROM VERNONS BOOK COLLIDER PHYSICS P 662
431
C MA2 AND MB2 ARE THE MASS SQUARED OF THE FINAL STATE PARTICLES
432
C 2-D PHASE SPACE = .5*PI*SQRT(1.,MA2/S^2,MB2/S^2)*(D(OMEGA)/4PI)
433
C****************************************************************************
435
LAMBDA=S**2+MA2**2+MB2**2-2.*S*MA2-2.*MA2*MB2-2.*S*MB2
440
SUBROUTINE TRANSPOLE(POLE,WIDTH,X,Y,JAC)
441
C**********************************************************************
442
C THIS ROUTINE TRANSFERS EVENLY SPACED X VALUES BETWEEN 0 AND 1
443
C TO Y VALUES WITH A POLE AT Y=POLE WITH WIDTH WIDTH AND RETURNS
444
C THE APPROPRIATE JACOBIAN FOR THIS
445
C**********************************************************************
450
REAL*8 POLE,WIDTH,Y,JAC
460
ZMIN = ATAN((-POLE)/WIDTH)/WIDTH
461
ZMAX = ATAN((1.-POLE)/WIDTH)/WIDTH
462
Z = ZMIN+(ZMAX-ZMIN)*X
463
Y = POLE+WIDTH*TAN(WIDTH*Z)
464
JAC= (WIDTH/COS(WIDTH*Z))**2*(ZMAX-ZMIN)