« back to all changes in this revision
Viewing changes to MG4_DECAY/decay_mom.f
- browse files at revision 36.21.85
- compare with another revision
- download diff
- download tarball
- view history from revision 36.21.85
- Committer: Chia-Hsien Shen
- Date: 2012-06-30 05:53:02 UTC
- mto: (36.21.97 decay_calculator_205) (238.3.11 1.6.0)
- mto: This revision was merged to the branch mainline in revision 249.
- Revision ID: anakin8120@gmail.com-20120630055302-5azmh3tx4vv2hlp1
move ./decay to ./mg5decay; resolve unit tests (n.b. __init__ does not check keys of input dictionaries, followed last revision)
- files added:
- DECAY
- mg5decay
- files removed:
- MG4_DECAY
- files modified:
- tests/input_files/full_sm_UFO/model.pkl
added
removed
MG4_DECAY/decay_mom.f
1
c***********************************************************************
2
c Routines for generating momentum of decay products
3
c***********************************************************************
4
5
SUBROUTINE PHASESPACE(PD,PSWGT)
6
C************************************************************************
7
C WRAPPER FOR CALLING THE PHASE SPACE ROUTINES
8
C************************************************************************
9
IMPLICIT NONE
10
include 'decay.inc'
11
C
12
C ARGUMENTS
13
C
14
REAL*8 PD(0:3,MaxDecPart)
15
REAL*8 PSWGT
16
C
17
C LOCAL
18
C
19
LOGICAL DONE
20
INTEGER IMAX
21
C
22
C----------
23
C START
24
C----------
25
C
26
PSWGT=0d0
27
IMAX=0
28
29
IF(ND.EQ.2) THEN
30
CALL TWOMOM(PD(0,1),PD(0,2),PD(0,3),PSWGT)
31
ELSEIF(ND.EQ.3) THEN
32
CALL THREEMOM(PD(0,1),PD(0,2),PD(0,3),PD(0,4),PSWGT)
33
ELSEIF(ND.EQ.4) THEN
34
CALL FOURMOM(PD(0,1),PD(0,2),PD(0,3),PD(0,4),PD(0,5),PSWGT)
35
ELSE
36
WRITE(*,*) 'NO PS AVAILABLE '
37
STOP
38
ENDIF
39
40
41
RETURN
42
END
43
44
45
SUBROUTINE TWOMOM(P1,P2,P3,PSWGT)
46
C************************************************************************
47
C GENERIC TWO BODY PHACE SPACE FOR THE DECAYS P1->P2+P3
48
C
49
C GIVEN P1 MOMENTUM SETS UP THE MOMENTA P2,P3
50
C ALSO SETS UP THE APPROPRIATE PHASESPACE WEIGTH PSWGT
51
C
52
C X(1) = COSTHETA OF FIRST DECAY
53
C X(2) = PHI OF FIRST DECAY
54
C
55
C************************************************************************
56
IMPLICIT NONE
57
include 'decay.inc'
58
C
59
C ARGUMENTS
60
C
61
REAL*8 P1(0:3),P2(0:3),P3(0:3)
62
REAL*8 PSWGT
63
C
64
C LOCAL
65
C
66
REAL*8 COSTH,PHI,JACPOLE
67
REAL*8 RS
68
REAL*8 XA2,XB2,S,JAC
69
REAL*8 a
70
integer i
71
C
72
C EXTERNAL
73
C
74
REAL*8 LAMBDA,DOT
75
C
76
C----------
77
C BEGIN
78
C----------
79
C
80
C CMS ENERGY = (P1)^2
81
C
82
S=DOT(P1,P1)
83
RS=SQRT(S)
84
C
85
C PICK VALUE OF THETA AND PHI
86
C
87
COSTH=-1.D0+2.D0*X(1)
88
PHI = 2D0*PI*X(2)
89
C
90
C DETERMINE JACOBIAN FOR THETA AND PHI
91
C
92
JAC = 4D0*PI
93
C
94
C CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
95
C OF THE FINAL STATE MASSLESS PARTONS IN THEIR CM FRAME
96
C
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
100
C
101
JAC = JAC /(2D0*PI)**2 !FOR THE (2 PI)^3N-4 OUT FRONT
102
C
103
C CALCULATE PHASE SPACE FACTOR DSIG/DCOS(THETA)
104
C
105
PSWGT = 1.D0/(2.D0*RS) !FLUX FACTOR
106
XA2 = M2*M2/S
107
XB2 = M3*M3/S
108
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI)
109
110
PSWGT=PSWGT*JAC
111
112
c write(*,*) 'from decay_mom',x(1),x(2),pswgt
113
RETURN
114
END
115
116
117
118
SUBROUTINE THREEMOM(P1,P2,P3,P4,PSWGT)
119
C************************************************************************
120
C This routine generates the momenta for the decay products of
121
C
122
C 1 -> 2 W -> 2 3 4
123
C
124
C when the W can be on-shell or not.
125
C It ALSO SETS UP THE APPROPRIATE PHASESPACE WEIGTH PSWGT
126
C
127
C X(1) = COSTHETA OF FIRST DECAY TOP -> W*+B
128
C X(2) = PHI OF FIRST DECAY TOP -> W*+B
129
C X(3) = M_W*
130
C X(4) = COSTHETA OF DECAY W* -> E+VE
131
C X(5) = PHI OF DECAY W* -> E+VE
132
C
133
C************************************************************************
134
IMPLICIT NONE
135
include 'decay.inc'
136
C
137
C ARGUMENTS
138
C
139
REAL*8 P1(0:3),P2(0:3),P3(0:3),P4(0:3),PW(0:3)
140
REAL*8 PSWGT
141
C
142
C LOCAL
143
C
144
REAL*8 COSTH,PHI,JACPOLE
145
REAL*8 RS
146
REAL*8 COSTH1,PHI1,S1,RS1
147
REAL*8 COSTH2,PHI2,S2,RS2
148
REAL*8 XA2,XB2,S,JAC
149
REAL*8 a
150
INTEGER I
151
C
152
C EXTERNAL
153
C
154
REAL*8 LAMBDA,DOT
155
C
156
C LOGICAL
157
C
158
LOGICAL FIRSTTIME
159
DATA FIRSTTIME /.TRUE./
160
C
161
C----------
162
C BEGIN
163
C----------
164
C
165
C CMS ENERGY
166
C
167
168
S=DOT(P1,P1)
169
RS=SQRT(S)
170
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)
173
S2 = S2*S
174
RS2 = SQRT(S2)
175
JAC=JACPOLE*S
176
ELSE ! W is off-shell
177
S2=(S-M2*M2)*X(3)+M2*M2 ! generated uniformily
178
RS2 = SQRT(S2)
179
JAC=(S-M2*M2)
180
ENDIF
181
C
182
C check total energy conservation
183
C
184
IF(RS2.GE.(RS-M2).or.RS2.LT.(M3+M4)) THEN
185
PSWGT=0D0
186
RETURN
187
ENDIF
188
C
189
C PICK VALUE OF THETA AND PHI FOR FIRST DECAY TOP -> W*+B
190
C
191
COSTH=-1.D0+2.D0*X(1)
192
PHI = 2D0*PI*X(2)
193
C
194
C DETERMINE JACOBIAN FOR THETA AND PHI
195
C
196
JAC = JAC * 4D0*PI
197
C
198
C CHOOSE COS(THETA) AND PHI FOR DECAY W*-> E+VE
199
C
200
COSTH1 = -1 + 2.*X(4)
201
PHI1 = 2D0*PI*X(5)
202
C
203
C JACOBIAN FACTOR FOR DECAY OF W*-> E+VE
204
C
205
JAC=JAC*4D0*PI
206
C
207
C CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
208
C OF THE FINAL STATE MASSLESS PARTONS IN THEIR CM FRAME
209
C
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
216
217
JAC = JAC /(2D0*PI)**5 !FOR THE (2 PI)^3N-4 OUT FRONT
218
C
219
C CALCULATE PHASE SPACE FACTOR DSIG/DCOS(THETA)
220
C
221
222
PSWGT = 1.D0/(2.D0*RS) !FLUX FACTOR
223
XA2 = S2/S
224
XB2 = M2*M2/S
225
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !T->W*+B
226
XA2 = M3*M3/S
227
XB2 = M4*M4/S
228
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !W*->E+VE
229
230
231
PSWGT=PSWGT*JAC
232
233
RETURN
234
END
235
236
237
SUBROUTINE FOURMOM(P1,P2,P3,P4,P5,PSWGT)
238
C************************************************************************
239
C This routine generates the momenta for the decay products of
240
C
241
C 1 -> v v'* -> 2 3 4 5
242
C
243
C where the v'* can be on-shell or not.
244
c
245
C GIVEN X(8) AND THE P1 MOMENTUM SETS UP THE MOMENTA
246
C ALSO SETS UP THE APPROPRIATE PHASESPACE WEIGTH PSWGT
247
C
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)
250
C X(3) = M_(2,3)
251
C X(4) = 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
256
C
257
C************************************************************************
258
IMPLICIT NONE
259
include "decay.inc"
260
C
261
C ARGUMENTS
262
C
263
REAL*8 P1(0:3),P2(0:3),P3(0:3),P4(0:3),P5(0:3)
264
REAL*8 PSWGT
265
C
266
C LOCAL
267
C
268
REAL*8 JAC,JACPOLE
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
272
REAL*8 XA2,XB2,SMIN
273
REAL*8 TEMP
274
INTEGER I
275
C
276
C EXTERNAL
277
C
278
REAL*8 LAMBDA,DOT
279
real xran1
280
integer iseed
281
data iseed/1/ !no effect if already called
282
C
283
C LOGICAL
284
C
285
LOGICAL FIRSTTIME
286
DATA FIRSTTIME /.TRUE./
287
C
288
C----------
289
C BEGIN
290
C----------
291
c
292
JAC=1D0
293
C
294
C CMS ENERGY = (HIGGS MASS)^2
295
C
296
SH=DOT(P1,P1)
297
RSH=SQRT(SH)
298
C
299
C PICK VALUE OF THETA AND PHI FOR FIRST DECAY SH->M3+S1
300
C
301
COSTH=-1.D0+2.D0*X(1)
302
PHI = 2D0*PI*X(2)
303
C
304
C DETERMINE JACOBIAN FOR THETA AND PHI
305
C
306
JAC = JAC * 4D0*PI
307
C
308
C MASSES OF THE DECAY 1->V V
309
C
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)
312
S1 = S1*SH
313
RS1 = SQRT(S1)
314
JAC=JACPOLE*SH*JAC
315
ELSE ! NO W'S CAN BE ON SHELL
316
S1=SH*X(3) ! generated uniformily
317
RS1 = SQRT(S1)
318
JAC=SH*JAC
319
ENDIF
320
321
IF(RS1.GE.RSH) THEN ! ENERGY CONSERVATION
322
PSWGT=0D0
323
RETURN
324
ENDIF
325
C
326
C DETERMINE INVARIENT MASS OF SYSTEM M4+M5
327
C
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)
330
S2 = S2*SH
331
RS2 = SQRT(S2)
332
JAC=JACPOLE*SH*JAC
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
338
RS2 = SQRT(S2)
339
JAC = JAC*SH*(X4MAX-X4MIN)
340
ENDIF
341
342
IF(RS2+RS1.GT.RSH) THEN ! ENERGY CONSERVATION
343
PSWGT=0D0
344
RETURN
345
ENDIF
346
c
347
c switch order to ensure no bias
348
c
349
IF(xran1(iseed).lt.0.5e0)then
350
temp=rs1
351
rs1=rs2
352
rs2=temp
353
temp=s1
354
s1=s2
355
s2=temp
356
ENDIF
357
358
IF(RS1 .LT. (M2+M3)) THEN ! MINIMUM inv mass
359
PSWGT=0D0
360
RETURN
361
ENDIF
362
363
IF(RS2 .LT. (M4+M5)) THEN ! MINIMUM inv mass
364
PSWGT=0D0
365
RETURN
366
ENDIF
367
368
369
C
370
C CHOOSE COS(THETA) AND PHI FOR DECAY S1->M3+M4
371
C
372
COSTH1 = -1 + 2.*X(5)
373
PHI1 = 2D0*PI*X(6)
374
C
375
C CHOOSE COS(THETA) AND PHI FOR DECAY S2->M4+M5
376
C
377
COSTH2 = -1 + 2.*X(7)
378
PHI2 = 2D0*PI*X(8)
379
C
380
C JACOBIAN FACTOR FOR DECAY OF S1->M3+M4
381
C
382
JAC=JAC*4D0*PI
383
C
384
C JACOBIAN FACTOR FOR DECAY OF S2->M4+M5
385
C
386
JAC=JAC*4D0*PI
387
C
388
C CALCULATE COMPONENTS OF MOMENTUM FOUR-VECTORS
389
C OF THE FINAL STATE
390
C
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
400
C
401
C
402
JAC = JAC /(2D0*PI)**8 !FOR THE (2 PI)^3N-4 OUT FRONT
403
C
404
C
405
C CALCULATE PHASE SPACE FACTOR DSIG/DCOS(THETA)
406
C
407
408
PSWGT = 1.D0/(2.D0*RSH) !FLUX FACTOR
409
XA2 = S1/SH
410
XB2 = S2/SH
411
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !S->S1+S2
412
XA2 = M2*M2/S1
413
XB2 = M3*M3/S1
414
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !S1->M2+M3
415
XA2 = M4*M4/S2
416
XB2 = M5*M5/S2
417
PSWGT = PSWGT*.5D0*PI*SQRT(LAMBDA(ONE,XA2,XB2))/(4.D0*PI) !S2->M4+M5
418
419
PSWGT=PSWGT*JAC
420
421
422
423
RETURN
424
END
425
426
427
REAL*8 FUNCTION LAMBDA(S,MA2,MB2)
428
IMPLICIT NONE
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****************************************************************************
434
REAL*8 MA2,MB2,S
435
LAMBDA=S**2+MA2**2+MB2**2-2.*S*MA2-2.*MA2*MB2-2.*S*MB2
436
RETURN
437
END
438
439
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**********************************************************************
446
IMPLICIT NONE
447
C
448
C ARGUMENTS
449
C
450
REAL*8 POLE,WIDTH,Y,JAC
451
REAL*4 X
452
C
453
C LOCAL
454
C
455
REAL*8 Z,ZMIN,ZMAX
456
c-----
457
c Begin Code
458
c-----
459
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)
465
RETURN
466
END
467
468
Loggerhead is a web-based interface for Breezy
Version: 2.0.1