C++ Interface to Tauola
tauola_extras.f
1
2 SUBROUTINE choice(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
3 $ AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
4 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
5 * ,ampiz,ampi,amro,gamro,ama1,gama1
6 * ,amk,amkz,amkst,gamkst
7C
8 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
9 * ,ampiz,ampi,amro,gamro,ama1,gama1
10 * ,amk,amkz,amkst,gamkst
11C
12 amrop=1.1
13 gamrop=0.36
14 amom=.782
15 gamom=0.0084
16C XXXXA CORRESPOND TO S2 CHANNEL !
17 IF(mnum.EQ.0) THEN
18 prob1=0.5
19 prob2=0.5
20 amrx =ama1
21 gamrx=gama1
22 amra =amro
23 gamra=gamro
24 amrb =amro
25 gamrb=gamro
26 ELSEIF(mnum.EQ.1) THEN
27 prob1=0.5
28 prob2=0.5
29 amrx =1.57
30 gamrx=0.9
31 amrb =amkst
32 gamrb=gamkst
33 amra =amro
34 gamra=gamro
35 ELSEIF(mnum.EQ.2) THEN
36 prob1=0.5
37 prob2=0.5
38 amrx =1.57
39 gamrx=0.9
40 amrb =amkst
41 gamrb=gamkst
42 amra =amro
43 gamra=gamro
44 ELSEIF(mnum.EQ.3) THEN
45 prob1=0.5
46 prob2=0.5
47 amrx =1.27
48 gamrx=0.3
49 amra =amkst
50 gamra=gamkst
51 amrb =amkst
52 gamrb=gamkst
53 ELSEIF(mnum.EQ.4) THEN
54 prob1=0.5
55 prob2=0.5
56 amrx =1.27
57 gamrx=0.3
58 amra =amkst
59 gamra=gamkst
60 amrb =amkst
61 gamrb=gamkst
62 ELSEIF(mnum.EQ.5) THEN
63 prob1=0.5
64 prob2=0.5
65 amrx =1.27
66 gamrx=0.3
67 amra =amkst
68 gamra=gamkst
69 amrb =amro
70 gamrb=gamro
71 ELSEIF(mnum.EQ.6) THEN
72 prob1=0.4
73 prob2=0.4
74 amrx =1.27
75 gamrx=0.3
76 amra =amro
77 gamra=gamro
78 amrb =amkst
79 gamrb=gamkst
80 ELSEIF(mnum.EQ.7) THEN
81 prob1=0.0
82 prob2=1.0
83 amrx =1.27
84 gamrx=0.9
85 amra =amro
86 gamra=gamro
87 amrb =amro
88 gamrb=gamro
89 ELSEIF(mnum.EQ.8) THEN
90 prob1=0.0
91 prob2=1.0
92 amrx =amrop
93 gamrx=gamrop
94 amrb =amom
95 gamrb=gamom
96 amra =amro
97 gamra=gamro
98 ELSEIF(mnum.EQ.101) THEN
99 prob1=.35
100 prob2=.35
101 amrx =1.2
102 gamrx=.46
103 amrb =amom
104 gamrb=gamom
105 amra =amom
106 gamra=gamom
107 ELSEIF(mnum.EQ.102) THEN
108 prob1=0.0
109 prob2=0.0
110 amrx =1.4
111 gamrx=.6
112 amrb =amom
113 gamrb=gamom
114 amra =amom
115 gamra=gamom
116 ELSE
117 prob1=0.0
118 prob2=0.0
119 amrx =ama1
120 gamrx=gama1
121 amra =amro
122 gamra=gamro
123 amrb =amro
124 gamrb=gamro
125 ENDIF
126C
127 IF (rr.LE.prob1) THEN
128 ichan=1
129 ELSEIF(rr.LE.(prob1+prob2)) THEN
130 ichan=2
131 ax =amra
132 gx =gamra
133 amra =amrb
134 gamra=gamrb
135 amrb =ax
136 gamrb=gx
137 px =prob1
138 prob1=prob2
139 prob2=px
140 ELSE
141 ichan=3
142 ENDIF
143C
144 prob3=1.0-prob1-prob2
145 END
146 SUBROUTINE initdk
147* ----------------------------------------------------------------------
148* INITIALISATION OF TAU DECAY PARAMETERS and routines
149*
150* called by : KORALZ
151* ----------------------------------------------------------------------
152
153 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
154 real*4 gfermi,gv,ga,ccabib,scabib,gamel
155 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
156 * ,ampiz,ampi,amro,gamro,ama1,gama1
157 * ,amk,amkz,amkst,gamkst
158*
159 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
160 * ,ampiz,ampi,amro,gamro,ama1,gama1
161 * ,amk,amkz,amkst,gamkst
162 COMMON / taubra / gamprt(30),jlist(30),nchan
163 COMMON / taukle / bra1,brk0,brk0b,brks
164 real*4 bra1,brk0,brk0b,brks
165
166
167
168
169
170
171 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
172 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
173 & ,names
174 CHARACTER NAMES(NMODE)*31
175
176 CHARACTER OLDNAMES(7)*31
177 CHARACTER*80 bxINIT
178 parameter(
179 $ bxinit ='(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
180 $ )
181 real*4 pi
182*
183*
184* LIST OF BRANCHING RATIOS
185CAM normalised to e nu nutau channel
186CAM enu munu pinu rhonu A1nu Knu K*nu pi
187CAM DATA JLIST / 1, 2, 3, 4, 5, 6, 7,
188*AM DATA GAMPRT /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,O.O811,0.616
189*AM
190*AM multipion decays
191*
192* conventions of particles names
193* K-,P-,K+, K0,P-,KB, K-,P0,K0
194* 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
195* P0,P0,K-, K-,P-,P+, P-,KB,P0
196* 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
197* ET,P-,P0 P-,P0,GM
198* 9, 1, 2 , 1, 2, 8
199*
200C
201 dimension nopik(6,nmode),npik(nmode)
202*AM outgoing multiplicity and flavors of multi-pion /multi-K modes
203 DATA npik / 4, 4,
204 1 5, 5,
205 2 6, 6,
206 3 3, 3,
207 4 3, 3,
208 5 3, 3,
209 6 3, 3,
210 7 2 /
211 DATA nopik / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
212 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
213 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
214 3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
215 4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
216 5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
217 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
218C AJWMOD fix sign bug, 2/22/99
219 7 -3,-4, 0, 0, 0, 0 /
220* LIST OF BRANCHING RATIOS
221 nchan = nmode + 7
222 DO 1 i = 1,30
223 IF (i.LE.nchan) THEN
224 jlist(i) = i
225 IF(i.EQ. 1) gamprt(i) =0.1800
226 IF(i.EQ. 2) gamprt(i) =0.1751
227 IF(i.EQ. 3) gamprt(i) =0.1110
228 IF(i.EQ. 4) gamprt(i) =0.2515
229 IF(i.EQ. 5) gamprt(i) =0.1790
230 IF(i.EQ. 6) gamprt(i) =0.0071
231 IF(i.EQ. 7) gamprt(i) =0.0134
232 IF(i.EQ. 8) gamprt(i) =0.0450
233 IF(i.EQ. 9) gamprt(i) =0.0100
234 IF(i.EQ.10) gamprt(i) =0.0009
235 IF(i.EQ.11) gamprt(i) =0.0004
236 IF(i.EQ.12) gamprt(i) =0.0003
237 IF(i.EQ.13) gamprt(i) =0.0005
238 IF(i.EQ.14) gamprt(i) =0.0015
239 IF(i.EQ.15) gamprt(i) =0.0015
240 IF(i.EQ.16) gamprt(i) =0.0015
241 IF(i.EQ.17) gamprt(i) =0.0005
242 IF(i.EQ.18) gamprt(i) =0.0050
243 IF(i.EQ.19) gamprt(i) =0.0055
244 IF(i.EQ.20) gamprt(i) =0.0017
245 IF(i.EQ.21) gamprt(i) =0.0013
246 IF(i.EQ.22) gamprt(i) =0.0010
247 IF(i.EQ. 1) oldnames(i)=' TAU- --> E- '
248 IF(i.EQ. 2) oldnames(i)=' TAU- --> MU- '
249 IF(i.EQ. 3) oldnames(i)=' TAU- --> PI- '
250 IF(i.EQ. 4) oldnames(i)=' TAU- --> PI-, PI0 '
251 IF(i.EQ. 5) oldnames(i)=' TAU- --> A1- (two subch) '
252 IF(i.EQ. 6) oldnames(i)=' TAU- --> K- '
253 IF(i.EQ. 7) oldnames(i)=' TAU- --> K*- (two subch) '
254 IF(i.EQ. 8) names(i-7)=' TAU- --> 2PI-, PI0, PI+ '
255 IF(i.EQ. 9) names(i-7)=' TAU- --> 3PI0, PI- '
256 IF(i.EQ.10) names(i-7)=' TAU- --> 2PI-, PI+, 2PI0 '
257 IF(i.EQ.11) names(i-7)=' TAU- --> 3PI-, 2PI+, '
258 IF(i.EQ.12) names(i-7)=' TAU- --> 3PI-, 2PI+, PI0 '
259 IF(i.EQ.13) names(i-7)=' TAU- --> 2PI-, PI+, 3PI0 '
260 IF(i.EQ.14) names(i-7)=' TAU- --> K-, PI-, K+ '
261 IF(i.EQ.15) names(i-7)=' TAU- --> K0, PI-, K0B '
262 IF(i.EQ.16) names(i-7)=' TAU- --> K-, K0, PI0 '
263 IF(i.EQ.17) names(i-7)=' TAU- --> PI0 PI0 K- '
264 IF(i.EQ.18) names(i-7)=' TAU- --> K- PI- PI+ '
265 IF(i.EQ.19) names(i-7)=' TAU- --> PI- K0B PI0 '
266 IF(i.EQ.20) names(i-7)=' TAU- --> ETA PI- PI0 '
267 IF(i.EQ.21) names(i-7)=' TAU- --> PI- PI0 GAM '
268 IF(i.EQ.22) names(i-7)=' TAU- --> K- K0 '
269 ELSE
270 jlist(i) = 0
271 gamprt(i) = 0.
272 ENDIF
273 1 CONTINUE
274 DO i=1,nmode
275 mulpik(i)=npik(i)
276 DO j=1,mulpik(i)
277 idffin(j,i)=nopik(j,i)
278 ENDDO
279 ENDDO
280*
281*
282* --- COEFFICIENTS TO FIX RATIO OF:
283* --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
284* --- PROBABILITY OF K0 TO BE KS
285* --- PROBABILITY OF K0B TO BE KS
286* --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
287* --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
288* --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
289* --- NEGLECTS MASS-PHASE SPACE EFFECTS
290 bra1=0.5
291 brk0=0.5
292 brk0b=0.5
293 brks=0.6667
294*
295
296 gfermi = 1.16637e-5
297 ccabib = 0.975
298 gv = 1.0
299 ga =-1.0
300
301
302
303* ZW 13.04.89 HERE WAS AN ERROR
304 scabib = sqrt(1.-ccabib**2)
305 pi =4.*atan(1.)
306 gamel = gfermi**2*amtau**5/(192*pi**3)
307*
308* CALL DEXAY(-1,pol1)
309*
310 RETURN
311 END
312 FUNCTION dcdmas(IDENT)
313 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
314 * ,ampiz,ampi,amro,gamro,ama1,gama1
315 * ,amk,amkz,amkst,gamkst
316*
317 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
318 * ,ampiz,ampi,amro,gamro,ama1,gama1
319 * ,amk,amkz,amkst,gamkst
320 IF (ident.EQ. 1) THEN
321 apkmas=ampi
322 ELSEIF (ident.EQ.-1) THEN
323 apkmas=ampi
324 ELSEIF (ident.EQ. 2) THEN
325 apkmas=ampiz
326 ELSEIF (ident.EQ.-2) THEN
327 apkmas=ampiz
328 ELSEIF (ident.EQ. 3) THEN
329 apkmas=amk
330 ELSEIF (ident.EQ.-3) THEN
331 apkmas=amk
332 ELSEIF (ident.EQ. 4) THEN
333 apkmas=amkz
334 ELSEIF (ident.EQ.-4) THEN
335 apkmas=amkz
336 ELSEIF (ident.EQ. 8) THEN
337 apkmas=0.0001
338 ELSEIF (ident.EQ.-8) THEN
339 apkmas=0.0001
340 ELSEIF (ident.EQ. 9) THEN
341 apkmas=0.5488
342 ELSEIF (ident.EQ.-9) THEN
343 apkmas=0.5488
344 ELSE
345 print *, 'STOP IN APKMAS, WRONG IDENT=',ident
346 stop
347 ENDIF
348 dcdmas=apkmas
349 END
350 FUNCTION lunpik(ID,ISGN)
351 COMMON / taukle / bra1,brk0,brk0b,brks
352 real*4 bra1,brk0,brk0b,brks
353 real*4 xio(1)
354 ident=id*isgn
355 IF (ident.EQ. 1) THEN
356 ipkdef=-211
357 ELSEIF (ident.EQ.-1) THEN
358 ipkdef= 211
359 ELSEIF (ident.EQ. 2) THEN
360 ipkdef=111
361 ELSEIF (ident.EQ.-2) THEN
362 ipkdef=111
363 ELSEIF (ident.EQ. 3) THEN
364 ipkdef=-321
365 ELSEIF (ident.EQ.-3) THEN
366 ipkdef= 321
367 ELSEIF (ident.EQ. 4) THEN
368*
369* K0 --> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
370 CALL ranmar(xio,1)
371 IF (xio(1).GT.brk0) THEN
372 ipkdef= 130
373 ELSE
374 ipkdef= 310
375 ENDIF
376 ELSEIF (ident.EQ.-4) THEN
377*
378* K0B--> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
379 CALL ranmar(xio,1)
380 IF (xio(1).GT.brk0b) THEN
381 ipkdef= 130
382 ELSE
383 ipkdef= 310
384 ENDIF
385 ELSEIF (ident.EQ. 8) THEN
386 ipkdef= 22
387 ELSEIF (ident.EQ.-8) THEN
388 ipkdef= 22
389 ELSEIF (ident.EQ. 9) THEN
390 ipkdef= 221
391 ELSEIF (ident.EQ.-9) THEN
392 ipkdef= 221
393 ELSE
394 print *, 'STOP IN IPKDEF, WRONG IDENT=',ident
395 stop
396 ENDIF
397 lunpik=ipkdef
398 END
399
400
401
402 SUBROUTINE taurdf(KTO)
403C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
404C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
405C CONTENTS
406 COMMON / taukle / bra1,brk0,brk0b,brks
407 real*4 bra1,brk0,brk0b,brks
408 COMMON / taubra / gamprt(30),jlist(30),nchan
409
410C Subroutine TAURDF is disabled
411 RETURN
412
413
414 IF (kto.EQ.1) THEN
415C ==================
416C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
417 bra1 = pkorb(4,1)
418 brks = pkorb(4,3)
419 brk0 = pkorb(4,5)
420 brk0b = pkorb(4,6)
421 ELSE
422C ====
423C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
424 bra1 = pkorb(4,2)
425 brks = pkorb(4,4)
426 brk0 = pkorb(4,5)
427 brk0b = pkorb(4,6)
428 ENDIF
429C =====
430 END
431
432 SUBROUTINE iniphy(XK00)
433* ----------------------------------------------------------------------
434* INITIALISATION OF PARAMETERS
435* USED IN QED and/or GSW ROUTINES
436* ----------------------------------------------------------------------
437 COMMON / qedprm /alfinv,alfpi,xk0
438 real*8 alfinv,alfpi,xk0
439 real*8 pi8,xk00
440*
441 pi8 = 4.d0*datan(1.d0)
442 alfinv = 137.03604d0
443 alfpi = 1d0/(alfinv*pi8)
444 xk0=xk00
445 END
446
447 SUBROUTINE inimas
448C ----------------------------------------------------------------------
449C INITIALISATION OF MASSES
450C
451C called by : KORALZ
452C ----------------------------------------------------------------------
453 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
454 * ,ampiz,ampi,amro,gamro,ama1,gama1
455 * ,amk,amkz,amkst,gamkst
456*
457 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
458 * ,ampiz,ampi,amro,gamro,ama1,gama1
459 * ,amk,amkz,amkst,gamkst
460C
461C IN-COMING / OUT-GOING FERMION MASSES
462 amtau = 1.7842
463C --- let us update tau mass ...
464 amtau = 1.777
465 amnuta = 0.010
466 amel = 0.0005111
467 amnue = 0.0
468 ammu = 0.105659
469 amnumu = 0.0
470*
471* MASSES USED IN TAU DECAYS
472 ampiz = 0.134964
473 ampi = 0.139568
474 amro = 0.773
475 gamro = 0.145
476*C GAMRO = 0.666
477 ama1 = 1.251
478 gama1 = 0.599
479 amk = 0.493667
480 amkz = 0.49772
481 amkst = 0.8921
482 gamkst = 0.0513
483C
484C
485C IN-COMING / OUT-GOING FERMION MASSES
486!! AMNUTA = PKORB(1,2)
487!! AMNUE = PKORB(1,4)
488!! AMNUMU = PKORB(1,6)
489C
490C MASSES USED IN TAU DECAYS Cleo settings
491!! AMPIZ = PKORB(1,7)
492!! AMPI = PKORB(1,8)
493!! AMRO = PKORB(1,9)
494!! GAMRO = PKORB(2,9)
495 ama1 = 1.275 !! PKORB(1,10)
496 gama1 = 0.615 !! PKORB(2,10)
497!! AMK = PKORB(1,11)
498!! AMKZ = PKORB(1,12)
499!! AMKST = PKORB(1,13)
500!! GAMKST = PKORB(2,13)
501C
502
503 RETURN
504 END
505 SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
506 real*8 pd1(4),pd2(4),q1(4),q2(4),costhe,p(4),qq(4),qt(4)
507C take effective beam which is less massive, it should be irrelevant
508C but in case HEPEVT is particulary dirty may help.
509C this routine calculate reduced system transver and cosine of scattering
510C angle.
511
512 xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
513 xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
514 IF (xm1.LT.xm2) THEN
515 sign=1d0
516 DO k=1,4
517 p(k)=pd1(k)
518 ENDDO
519 ELSE
520 sign=-1d0
521 DO k=1,4
522 p(k)=pd2(k)
523 ENDDO
524 ENDIF
525C calculate space like part of P (in Z restframe)
526 DO k=1,4
527 qq(k)=q1(k)+q2(k)
528 qt(k)=q1(k)-q2(k)
529 ENDDO
530
531 xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
532
533 qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
534 DO k=1,4
535 qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
536 ENDDO
537
538 pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
539 DO k=1,4
540 p(k)=p(k)-qq(k)*pxqq/xmqq**2
541 ENDDO
542C calculate costhe
543 pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
544 qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
545 pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
546 costhe=pxqt/pxp/qtxqt
547 costhe=costhe*sign
548 END
549
550 FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
551C this function calculates probability for the helicity +1 +1 configuration
552C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
553 real*8 plzap0,svar,costhe,costh0,t_born
554
555 costhe=costh0
556C >>>>> IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
557C >>>>> of first beam is used by T_GIVIZ0 including sign
558
559 IF (idf.GT.0) THEN
560 CALL initwk(ide,idf,svar)
561 ELSE
562 CALL initwk(-ide,-idf,svar)
563 ENDIF
564 plzap0=t_born(0,svar,costhe,1d0,1d0)
565 $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
566
567
568C write(*,*) 'svar= ', svar
569C write(*,*) 'COSTHE=', COSTHE
570C write(*,*) ide,' ',idf
571C write(*,*) 'PLZAP0=', PLZAP0
572C COSTHE=0.999
573C write(*,*) 'TBORN+=', T_BORN(0,SVAR,COSTHE,1D0,1D0)
574
575C COSTHE=-0.999
576C write(*,*) 'TBORN-=', T_BORN(0,SVAR,COSTHE,-1D0,-1D0)
577C 100 format (A,E8.4)
578
579! PLZAP0=0.5
580 END
581 FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
582C ----------------------------------------------------------------------
583C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
584C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
585C EXAMPLE OF THE METHOD APPLIED THERE
586C INPUT PARAMETERS ARE: SVAR -- transfer
587C COSTHE -- cosine of angle between tau+ and 1st beam
588C TA,TB -- helicity states of tau+ tau-
589C
590C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
591C WARNING: temporarily tau-mass terms inactivated, Dec 10, 2019
592C ----------------------------------------------------------------------
593 IMPLICIT REAL*8(a-h,o-z)
594 COMMON / t_beampm / ene ,amin,amfin,ide,idf
595 real*8 ene ,amin,amfin
596 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
597 & ,xupgi ,xupzi ,xupgf ,xupzf
598 & ,ndiag0,ndiaga,keya,keyz
599 & ,itce,jtce,itcf,jtcf,kolor
600 real*8 ss,poln,t3e,qe,t3f,qf
601 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
602 real*8 seps1,seps2
603C=====================================================================
604 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
605 real*8 swsq,amw,amz,amh,amtop,gammz
606C SWSQ = sin2 (theta Weinberg)
607C AMW,AMZ = W & Z boson masses respectively
608C AMH = the Higgs mass
609C AMTOP = the top mass
610C GAMMZ = Z0 width
611 COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
612 COMPLEX*16 XUPZFP(2),XUPZIP(2)
613 COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
614 COMPLEX*16 PROPA,PROPZ
615 COMPLEX*16 XR,XI
616 COMPLEX*16 XUPF,XUPI
617 COMPLEX*16 XTHING
618 DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
619 DATA mode0 /-5/
620 DATA ide0 /-55/
621 DATA svar0,cost0 /-5.d0,-6.d0/
622 DATA pi /3.141592653589793238462643d0/
623 DATA seps1,seps2 /0d0,0d0/
624C
625C MEMORIZATION =========================================================
626 IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
627 $ .OR.ide0.NE.ide)THEN
628C
629 keygsw=1
630C ** PROPAGATORS
631 ide0=ide
632 mode0=mode
633 svar0=svar
634 cost0=costhe
635 sinthe=sqrt(1.d0-costhe**2)
636 beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
637 beta=1.0 ! zbw 10.12.2019 necessary kill for tauspinner electroweak
638C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
639 xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5d0*beta*(xupzf(1)-xupzf(2))
640 xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5d0*beta*(xupzf(1)-xupzf(2))
641 xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5d0*(xupzi(1)-xupzi(2))
642 xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5d0*(xupzi(1)-xupzi(2))
643C FINAL STATE VECTOR COUPLING
644 xupf =0.5d0*(xupzf(1)+xupzf(2))
645 xupi =0.5d0*(xupzi(1)+xupzi(2))
646 xthing =0d0
647
648 propa =1d0/svar
649 propa =propa *(137.03604d0/128.86674175d0) ! alpha(m_z)/alpha(0)
650 ! note t_born call
651 ! in tau_reweight_lib.cxx
652 propz =1d0/dcmplx(svar-amz**2,amz*gammz)
653 ! zbw Dec 10, 2019 : residuum of Z adjusted now to effective born
654 gfermi = 1.166389e-5
655 alphainv=137.03604d0
656 zetvpi = gfermi *amz**2 *alphainv /(dsqrt(2.d0)*8.d0*pi)
657 $ *(swsq*(1d0-swsq)) *16d0
658 propz =propz *zetvpi
659! PROPA=PROPA*1.037
660! write(*,*) 'GF=',gfermi,' ',ALPHAINV,' ',swsq,' ',amz,' == ',ZetVPi
661 IF (keygsw.EQ.0) propz=0.d0
662 DO 50 i=1,2
663 DO 50 j=1,2
664 regula= (3-2*i)*(3-2*j) + costhe
665 regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
666 aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
667 azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
668 aborn(i,j)=aphot(i,j)+azett(i,j)
669 aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
670 azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
671 abornm(i,j)=aphotm(i,j)+azettm(i,j)
672 50 CONTINUE
673 ENDIF
674C
675C******************
676C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
677C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
678C* HELICITY CONSERVATION EXPLICITLY OBEYED
679 polar1= (seps1)
680 polar2= (-seps2)
681 born=0d0
682 DO 150 i=1,2
683 helic= 3-2*i
684 DO 150 j=1,2
685 helit=3-2*j
686 factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
687 factom=factor*(1+helit*ta)*(1-helit*tb)
688 factor=factor*(1+helit*ta)*(1+helit*tb)
689
690 born=born+cdabs(aborn(i,j))**2*factor
691C MASS TERM IN BORN
692 IF (mode.GE.1) THEN
693 born=born+cdabs(abornm(i,j))**2*factom
694 ENDIF
695
696 150 CONTINUE
697C************
698 funt=born
699 IF(funt.LT.0.d0) funt=born
700
701C
702 IF (svar.GT.4d0*amfin**2) THEN
703C PHASE SPACE THRESHOLD FACTOR
704 thresh=sqrt(1-4d0*amfin**2/svar)
705 thresh=1.d0 ! zbw 10.12.2019 necessary kill for tauspinner electroweak
706 t_born= funt*svar**2*thresh
707 ELSE
708 thresh=0.d0
709 t_born=0.d0
710 ENDIF
711 END
712
713 SUBROUTINE initwk(IDEX,IDFX,SVAR)
714! initialization routine coupling masses etc.
715 IMPLICIT REAL*8 (a-h,o-z)
716 COMMON / t_beampm / ene ,amin,amfin,ide,idf
717 real*8 ene ,amin,amfin
718 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
719 & ,xupgi ,xupzi ,xupgf ,xupzf
720 & ,ndiag0,ndiaga,keya,keyz
721 & ,itce,jtce,itcf,jtcf,kolor
722 real*8 ss,poln,t3e,qe,t3f,qf
723 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
724 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
725 real*8 swsq,amw,amz,amh,amtop,gammz
726C SWSQ = sin2 (theta Weinberg)
727C AMW,AMZ = W & Z boson masses respectively
728C AMH = the Higgs mass
729C AMTOP = the top mass
730C GAMMZ = Z0 width
731C
732 ene=sqrt(svar)/2
733 amin=0.511d-3
734 swsq=0.2315200d0 ! 0.23147
735 amz=91.18870000000d0! 91.1882
736 gammz=2.4952d0
737 IF (idfx.EQ. 15) then
738 idf=2 ! denotes tau +2 tau-
739 amfin=1.77703 !this mass is irrelevant if small, used in ME only
740 ELSEIF (idfx.EQ.-15) then
741 idf=-2 ! denotes tau -2 tau-
742 amfin=1.77703 !this mass is irrelevant if small, used in ME only
743 ELSE
744 WRITE(*,*) 'INITWK: WRONG IDFX'
745 stop
746 ENDIF
747
748 IF (idex.EQ. 11) then !electron
749 ide= 2
750 amin=0.511d-3
751 ELSEIF (idex.EQ.-11) then !positron
752 ide=-2
753 amin=0.511d-3
754 ELSEIF (idex.EQ. 13) then !mu+
755 ide= 2
756 amin=0.105659
757 ELSEIF (idex.EQ.-13) then !mu-
758 ide=-2
759 amin=0.105659
760 ELSEIF (idex.EQ. 1) then !d
761 ide= 4
762 amin=0.05
763 ELSEIF (idex.EQ.- 1) then !d~
764 ide=-4
765 amin=0.05
766 ELSEIF (idex.EQ. 2) then !u
767 ide= 3
768 amin=0.02
769 ELSEIF (idex.EQ.- 2) then !u~
770 ide=-3
771 amin=0.02
772 ELSEIF (idex.EQ. 3) then !s
773 ide= 4
774 amin=0.3
775 ELSEIF (idex.EQ.- 3) then !s~
776 ide=-4
777 amin=0.3
778 ELSEIF (idex.EQ. 4) then !c
779 ide= 3
780 amin=1.3
781 ELSEIF (idex.EQ.- 4) then !c~
782 ide=-3
783 amin=1.3
784 ELSEIF (idex.EQ. 5) then !b
785 ide= 4
786 amin=4.5
787 ELSEIF (idex.EQ.- 5) then !b~
788 ide=-4
789 amin=4.5
790 ELSEIF (idex.EQ. 12) then !nu_e
791 ide= 1
792 amin=0.1d-3
793 ELSEIF (idex.EQ.- 12) then !nu_e~
794 ide=-1
795 amin=0.1d-3
796 ELSEIF (idex.EQ. 14) then !nu_mu
797 ide= 1
798 amin=0.1d-3
799 ELSEIF (idex.EQ.- 14) then !nu_mu~
800 ide=-1
801 amin=0.1d-3
802 ELSEIF (idex.EQ. 16) then !nu_tau
803 ide= 1
804 amin=0.1d-3
805 ELSEIF (idex.EQ.- 16) then !nu_tau~
806 ide=-1
807 amin=0.1d-3
808
809 ELSE
810 WRITE(*,*) 'INITWK: WRONG IDEX'
811 stop
812 ENDIF
813
814C ----------------------------------------------------------------------
815C
816C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
817C
818C called by : KORALZ
819C ----------------------------------------------------------------------
820 itce=ide/iabs(ide)
821 jtce=(1-itce)/2
822 itcf=idf/iabs(idf)
823 jtcf=(1-itcf)/2
824 CALL t_givizo( ide, 1,aizor,qe,kdumm)
825 CALL t_givizo( ide,-1,aizol,qe,kdumm)
826 xupgi(1)=qe
827 xupgi(2)=qe
828 t3e = aizol+aizor
829 xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
830 xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
831 CALL t_givizo( idf, 1,aizor,qf,kolor)
832 CALL t_givizo( idf,-1,aizol,qf,kolor)
833 xupgf(1)=qf
834 xupgf(2)=qf
835 t3f = aizol+aizor
836 xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
837 xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
838C
839 ndiag0=2
840 ndiaga=11
841 keya = 1
842 keyz = 1
843C
844C
845 RETURN
846 END
847
848 SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
849C ----------------------------------------------------------------------
850C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
851C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
852C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
853C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
854C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
855C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
856C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
857C
858C called by : EVENTE, EVENTM, FUNTIH, .....
859C ----------------------------------------------------------------------
860 IMPLICIT REAL*8(a-h,o-z)
861C
862 IF(idferm.EQ.0.OR.iabs(idferm).GT.4) GOTO 901
863 IF(iabs(ihelic).NE.1) GOTO 901
864 ih =ihelic
865 idtype =iabs(idferm)
866 ic =idferm/idtype
867 lepqua=int(idtype*0.4999999d0)
868 iupdow=idtype-2*lepqua-1
869 charge =(-iupdow+2d0/3d0*lepqua)*ic
870 sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
871 kolor=1+2*lepqua
872C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
873C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
874 RETURN
875 901 print *,' STOP IN GIVIZO: WRONG PARAMS.'
876 stop
877 END
878 SUBROUTINE phyfix(NSTOP,NSTART)
879 common/lujets/n,k(4000,5),p(4000,5),v(4000,5)
880 SAVE /lujets/
881C NSTOP NSTART : when PHYTIA history ends and event starts.
882 nstop=0
883 nstart=1
884 DO i=1, n
885 IF(k(i,1).NE.21) THEN
886 nstop = i-1
887 nstart= i
888 GOTO 500
889 ENDIF
890 ENDDO
891 500 CONTINUE
892 END
893
894
895 SUBROUTINE taupi0(PI0,K)
896C no initialization required. Must be called once after every:
897C 1) CALL DEKAY(1+10,...)
898C 2) CALL DEKAY(2+10,...)
899C 3) CALL DEXAY(1,...)
900C 4) CALL DEXAY(2,...)
901C subroutine to decay originating from TAUOLA's taus:
902C 1) etas (with CALL TAUETA(JAK))
903C 2) later pi0's from taus.
904C 3) extensions to other applications possible.
905C this routine belongs to >tauola universal interface<, but uses
906C routines from >tauola< utilities as well. 25.08.2005
907C this is the hepevt class in old style. No d_h_ class pre-name
908
909C position of taus, must be defined by host program:
910 COMMON /taupos/ np1,np2
911c
912 REAL PHOT1(4),PHOT2(4)
913 real*8 r,x(4),y(4),pi0(4)
914 INTEGER JEZELI(3),K
915 DATA jezeli /0,0,0/
916 SAVE jezeli
917
918! random 3 vector on the sphere, masless
919 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
920 CALL spherd(r,x)
921 x(4)=r
922 y(4)=r
923
924 y(1)=-x(1)
925 y(2)=-x(2)
926 y(3)=-x(3)
927! boost to lab and to real*4
928 CALL bostdq(-1,pi0,x,x)
929 CALL bostdq(-1,pi0,y,y)
930 DO l=1,4
931 phot1(l)=x(l)
932 phot2(l)=y(l)
933 ENDDO
934C to hepevt
935 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
936 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
937
938C
939 END
940 SUBROUTINE taueta(PETA,K)
941C subroutine to decay etas's from taus.
942C this routine belongs to tauola universal interface, but uses
943C routines from tauola utilities. Just flat phase space, but 4 channels.
944C it is called at the beginning of SUBR. TAUPI0(JAK)
945C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
946C this is the hepevt class in old style. No d_h_ class pre-name
947
948 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
949 * ,ampiz,ampi,amro,gamro,ama1,gama1
950 * ,amk,amkz,amkst,gamkst
951*
952 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
953 * ,ampiz,ampi,amro,gamro,ama1,gama1
954 * ,amk,amkz,amkst,gamkst
955
956C position of taus, must be defined by host program:
957c
958 REAL RRR(1),BRSUM(3), RR(2)
959 REAL PHOT1(4),PHOT2(4),PHOT3(4)
960 real*8 x(4), y(4), z(4)
961 REAL YM1,YM2,YM3
962 real*8 r,ru,peta(4),xm1,xm2,xm3,xlam
963 real*8 a,b,c
964 INTEGER K
965 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
966C position of decaying particle:
967C DO L=1,4
968C PETA(L)= phep(L,K) ! eta 4 momentum
969C ENDDO
970C eta cumulated branching ratios:
971 brsum(1)=0.389 ! gamma gamma
972 brsum(2)=brsum(1)+0.319 ! 3 pi0
973 brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
974 CALL ranmar(rrr,1)
975
976 IF (rrr(1).LT.brsum(1)) THEN ! gamma gamma channel exactly like pi0
977! random 3 vector on the sphere, masless
978 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
979 CALL spherd(r,x)
980 x(4)=r
981 y(4)=r
982
983 y(1)=-x(1)
984 y(2)=-x(2)
985 y(3)=-x(3)
986! boost to lab and to real*4
987 CALL bostdq(-1,peta,x,x)
988 CALL bostdq(-1,peta,y,y)
989 DO l=1,4
990 phot1(l)=x(l)
991 phot2(l)=y(l)
992 ENDDO
993C to hepevt
994 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
995 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
996 ELSE ! 3 body channels
997 IF(rrr(1).LT.brsum(2)) THEN ! 3 pi0
998 id1= 111
999 id2= 111
1000 id3= 111
1001 xm1=ampiz ! masses
1002 xm2=ampiz
1003 xm3=ampiz
1004 ELSEIF(rrr(1).LT.brsum(3)) THEN ! pi+ pi- pi0
1005 id1= 211
1006 id2=-211
1007 id3= 111
1008 xm1=ampi ! masses
1009 xm2=ampi
1010 xm3=ampiz
1011 ELSE ! pi+ pi- gamma
1012 id1= 211
1013 id2=-211
1014 id3= 22
1015 xm1=ampi ! masses
1016 xm2=ampi
1017 xm3=0.0
1018 ENDIF
1019 7 CONTINUE ! we generate mass of the first pair:
1020 CALL ranmar(rr,2)
1021 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1022 amin=xm1+xm2
1023 amax=r-xm3
1024 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1025C weight for flat phase space
1026 wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1027 & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1028 & /r**2 /am2**2
1029 IF (rr(2).GT.wt) GOTO 7
1030
1031 ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2 ! momenta of the
1032 ! first two products
1033 ! in the rest frame of that pair
1034 CALL spherd(ru,x)
1035 x(4)=sqrt(ru**2+xm1**2)
1036 y(4)=sqrt(ru**2+xm2**2)
1037
1038 y(1)=-x(1)
1039 y(2)=-x(2)
1040 y(3)=-x(3)
1041C generate momentum of that pair in rest frame of eta:
1042 ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1043 CALL spherd(ru,z)
1044 z(4)=sqrt(ru**2+am2**2)
1045C and boost first two decay products to rest frame of eta.
1046 CALL bostdq(-1,z,x,x)
1047 CALL bostdq(-1,z,y,y)
1048C redefine Z(4) to 4-momentum of the last decay product:
1049 z(1)=-z(1)
1050 z(2)=-z(2)
1051 z(3)=-z(3)
1052 z(4)=sqrt(ru**2+xm3**2)
1053C boost all to lab and move to real*4; also masses
1054 CALL bostdq(-1,peta,x,x)
1055 CALL bostdq(-1,peta,y,y)
1056 CALL bostdq(-1,peta,z,z)
1057 DO l=1,4
1058 phot1(l)=x(l)
1059 phot2(l)=y(l)
1060 phot3(l)=z(l)
1061 ENDDO
1062 ym1=xm1
1063 ym2=xm2
1064 ym3=xm3
1065C to hepevt
1066 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1067 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1068 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1069 ENDIF
1070
1071
1072C
1073 END
1074 SUBROUTINE tauk0s(PETA,K)
1075C subroutine to decay K0S's from taus.
1076C this routine belongs to tauola universal interface, but uses
1077C routines from tauola utilities. Just flat phase space, but 4 channels.
1078C it is called at the beginning of SUBR. TAUPI0(JAK)
1079C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1080
1081 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1082 * ,ampiz,ampi,amro,gamro,ama1,gama1
1083 * ,amk,amkz,amkst,gamkst
1084*
1085 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1086 * ,ampiz,ampi,amro,gamro,ama1,gama1
1087 * ,amk,amkz,amkst,gamkst
1088
1089C position of taus, must be defined by host program:
1090 COMMON /taupos/ np1,np2
1091c
1092 REAL RRR(1),BRSUM(3)
1093 REAL PHOT1(4),PHOT2(4)
1094 real*8 x(4), y(4)
1095 REAL YM1,YM2
1096 real*8 r,peta(4),xm1,xm2,xlam
1097 real*8 a,b,c
1098 INTEGER K
1099 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1100C position of decaying particle:
1101
1102
1103! DO L=1,4
1104! PETA(L)= phep(L,K) ! K0S 4 momentum (this is cloned from eta decay)
1105! ENDDO
1106C K0S cumulated branching ratios:
1107 brsum(1)=0.313 ! 2 PI0
1108 brsum(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI-
1109 brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1110 CALL ranmar(rrr,1)
1111
1112 IF(rrr(1).LT.brsum(1)) THEN ! 2 pi0
1113 id1= 111
1114 id2= 111
1115 xm1=ampiz ! masses
1116 xm2=ampiz
1117 ELSEIF(rrr(1).LT.brsum(2)) THEN ! pi+ pi-
1118 id1= 211
1119 id2=-211
1120 xm1=ampi ! masses
1121 xm2=ampi
1122 ELSE ! gamma gamma unused !!!
1123 id1= 22
1124 id2= 22
1125 xm1= 0.0 ! masses
1126 xm2= 0.0
1127 ENDIF
1128
1129! random 3 vector on the sphere, of equal mass !!
1130 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1131 r4=r
1132 r=sqrt(abs(r**2-xm1**2))
1133 CALL spherd(r,x)
1134 x(4)=r4
1135 y(4)=r4
1136
1137 y(1)=-x(1)
1138 y(2)=-x(2)
1139 y(3)=-x(3)
1140! boost to lab and to real*4
1141 CALL bostdq(-1,peta,x,x)
1142 CALL bostdq(-1,peta,y,y)
1143 DO l=1,4
1144 phot1(l)=x(l)
1145 phot2(l)=y(l)
1146 ENDDO
1147
1148 ym1=xm1
1149 ym2=xm2
1150C to hepevt
1151 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1152 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1153
1154C
1155 END
1156
1157 subroutine bostdq(idir,vv,pp,q)
1158* *******************************
1159c Boost along arbitrary vector v (see eg. J.D. Jacson, Classical
1160c Electrodynamics).
1161c Four-vector pp is boosted from an actual frame to the rest frame
1162c of the four-vector v (for idir=1) or back (for idir=-1).
1163c q is a resulting four-vector.
1164c Note: v must be time-like, pp may be arbitrary.
1165c
1166c Written by: Wieslaw Placzek date: 22.07.1994
1167c Last update: 3/29/95 by: M.S.
1168c
1169 implicit DOUBLE PRECISION (a-h,o-z)
1170 parameter(nout=6)
1171 DOUBLE PRECISION v(4),p(4),q(4),pp(4),vv(4)
1172 save
1173!
1174 do 1 i=1,4
1175 v(i)=vv(i)
1176 1 p(i)=pp(i)
1177 amv=(v(4)**2-v(1)**2-v(2)**2-v(3)**2)
1178 if (amv.le.0d0) then
1179 write(6,*) 'bosstv: warning amv**2=',amv
1180 endif
1181 amv=sqrt(abs(amv))
1182 if (idir.eq.-1) then
1183 q(4)=( p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4))/amv
1184 wsp =(q(4)+p(4))/(v(4)+amv)
1185 elseif (idir.eq.1) then
1186 q(4)=(-p(1)*v(1)-p(2)*v(2)-p(3)*v(3)+p(4)*v(4))/amv
1187 wsp =-(q(4)+p(4))/(v(4)+amv)
1188 else
1189 write(nout,*)' >>> boostv: wrong value of idir = ',idir
1190 endif
1191 q(1)=p(1)+wsp*v(1)
1192 q(2)=p(2)+wsp*v(2)
1193 q(3)=p(3)+wsp*v(3)
1194 end
1195
1196
1197
1198