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
7 C
8  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
9  * ,ampiz,ampi,amro,gamro,ama1,gama1
10  * ,amk,amkz,amkst,gamkst
11 C
12  amrop=1.1
13  gamrop=0.36
14  amom=.782
15  gamom=0.0084
16 C 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
126 C
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
143 C
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
185 CAM normalised to e nu nutau channel
186 CAM enu munu pinu rhonu A1nu Knu K*nu pi
187 CAM 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 *
200 C
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,
218 C 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)
403 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
404 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
405 C CONTENTS
406  COMMON / taukle / bra1,brk0,brk0b,brks
407  real*4 bra1,brk0,brk0b,brks
408  COMMON / taubra / gamprt(30),jlist(30),nchan
409 
410 C Subroutine TAURDF is disabled
411  RETURN
412 
413 
414  IF (kto.EQ.1) THEN
415 C ==================
416 C 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
422 C ====
423 C 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
429 C =====
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
448 C ----------------------------------------------------------------------
449 C INITIALISATION OF MASSES
450 C
451 C called by : KORALZ
452 C ----------------------------------------------------------------------
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
460 C
461 C IN-COMING / OUT-GOING FERMION MASSES
462  amtau = 1.7842
463 C --- 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
483 C
484 C
485 C IN-COMING / OUT-GOING FERMION MASSES
486 !! AMNUTA = PKORB(1,2)
487 !! AMNUE = PKORB(1,4)
488 !! AMNUMU = PKORB(1,6)
489 C
490 C 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)
501 C
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)
507 C take effective beam which is less massive, it should be irrelevant
508 C but in case HEPEVT is particulary dirty may help.
509 C this routine calculate reduced system transver and cosine of scattering
510 C 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
525 C 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
542 C 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)
551 C this function calculates probability for the helicity +1 +1 configuration
552 C 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
556 C >>>>> IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
557 C >>>>> 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 
568 C write(*,*) 'svar= ', svar
569 C write(*,*) 'COSTHE=', COSTHE
570 C write(*,*) ide,' ',idf
571 C write(*,*) 'PLZAP0=', PLZAP0
572 C COSTHE=0.999
573 C write(*,*) 'TBORN+=', T_BORN(0,SVAR,COSTHE,1D0,1D0)
574 
575 C COSTHE=-0.999
576 C write(*,*) 'TBORN-=', T_BORN(0,SVAR,COSTHE,-1D0,-1D0)
577 C 100 format (A,E8.4)
578 
579 ! PLZAP0=0.5
580  END
581  FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
582 C ----------------------------------------------------------------------
583 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
584 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
585 C EXAMPLE OF THE METHOD APPLIED THERE
586 C INPUT PARAMETERS ARE: SVAR -- transfer
587 C COSTHE -- cosine of angle between tau+ and 1st beam
588 C TA,TB -- helicity states of tau+ tau-
589 C
590 C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
591 C WARNING: temporarily tau-mass terms inactivated, Dec 10, 2019
592 C ----------------------------------------------------------------------
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
603 C=====================================================================
604  COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
605  real*8 swsq,amw,amz,amh,amtop,gammz
606 C SWSQ = sin2 (theta Weinberg)
607 C AMW,AMZ = W & Z boson masses respectively
608 C AMH = the Higgs mass
609 C AMTOP = the top mass
610 C 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/
624 C
625 C MEMORIZATION =========================================================
626  IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
627  $ .OR.ide0.NE.ide)THEN
628 C
629  keygsw=1
630 C ** 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
638 C 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))
643 C 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) ! poor man's alpha running
650  ! corr.: alpha(m_z)/alpha(0)
651  ! t_born is called
652  ! in tau_reweight_lib.cxx
653  ! not perfect choice
654  ! for sig of small SVAR's,
655  ! but spin correlations OK.
656  propz =1d0/dcmplx(svar-amz**2,amz*gammz)
657  ! zbw Dec 10, 2019 : residuum of Z adjusted now to effective born
658  gfermi = 1.166389e-5
659  alphainv=128.86674175 !137.03604D0
660  zetvpi = gfermi *amz**2 *alphainv /(dsqrt(2.d0)*8.d0*pi)
661  $ *(swsq*(1d0-swsq)) *16d0
662  propz =propz *zetvpi
663  propz =propz *(137.03604d0/128.86674175d0)
664 ! PROPA=PROPA*1.037
665 ! write(*,*) 'GF=',gfermi,' ',ALPHAINV,' ',swsq,' ',amz,' == ',ZetVPi
666  IF (keygsw.EQ.0) propz=0.d0
667  DO 50 i=1,2
668  DO 50 j=1,2
669  regula= (3-2*i)*(3-2*j) + costhe
670  regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
671  aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
672  azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
673  aborn(i,j)=aphot(i,j)+azett(i,j)
674  aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
675  azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
676  abornm(i,j)=aphotm(i,j)+azettm(i,j)
677  50 CONTINUE
678  ENDIF
679 C
680 C******************
681 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
682 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
683 C* HELICITY CONSERVATION EXPLICITLY OBEYED
684  polar1= (seps1)
685  polar2= (-seps2)
686  born=0d0
687  DO 150 i=1,2
688  helic= 3-2*i
689  DO 150 j=1,2
690  helit=3-2*j
691  factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
692  factom=factor*(1+helit*ta)*(1-helit*tb)
693  factor=factor*(1+helit*ta)*(1+helit*tb)
694 
695  born=born+cdabs(aborn(i,j))**2*factor
696 C MASS TERM IN BORN
697  IF (mode.GE.1) THEN
698  born=born+cdabs(abornm(i,j))**2*factom
699  ENDIF
700 
701  150 CONTINUE
702 C************
703  funt=born
704  IF(funt.LT.0.d0) funt=born
705 
706 C
707  IF (svar.GT.4d0*amfin**2) THEN
708 C PHASE SPACE THRESHOLD FACTOR
709  thresh=sqrt(1-4d0*amfin**2/svar)
710 ! desactivated, necessary for not anymore recommended solution !
711 ! THRESH=1.D0 ! zbw 10.12.2019 necessary kill for tauspinner electroweak
712  t_born= funt*svar**2*thresh
713  ELSE
714  thresh=0.d0
715  t_born=0.d0
716  ENDIF
717  END
718 
719  SUBROUTINE initwk(IDEX,IDFX,SVAR)
720 ! initialization routine coupling masses etc.
721  IMPLICIT REAL*8 (a-h,o-z)
722  COMMON / t_beampm / ene ,amin,amfin,ide,idf
723  real*8 ene ,amin,amfin
724  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
725  & ,xupgi ,xupzi ,xupgf ,xupzf
726  & ,ndiag0,ndiaga,keya,keyz
727  & ,itce,jtce,itcf,jtcf,kolor
728  real*8 ss,poln,t3e,qe,t3f,qf
729  & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
730  COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
731  real*8 swsq,amw,amz,amh,amtop,gammz
732 C SWSQ = sin2 (theta Weinberg)
733 C AMW,AMZ = W & Z boson masses respectively
734 C AMH = the Higgs mass
735 C AMTOP = the top mass
736 C GAMMZ = Z0 width
737 C
738  ene=sqrt(svar)/2
739  amin=0.511d-3
740  swsq=0.2315200d0 ! 0.23147
741  amz=91.18870000000d0! 91.1882
742  gammz=2.4952d0
743  IF (idfx.EQ. 15) then
744  idf=2 ! denotes tau +2 tau-
745  amfin=1.77703 !this mass is irrelevant if small, used in ME only
746  ELSEIF (idfx.EQ.-15) then
747  idf=-2 ! denotes tau -2 tau-
748  amfin=1.77703 !this mass is irrelevant if small, used in ME only
749  ELSE
750  WRITE(*,*) 'INITWK: WRONG IDFX'
751  stop
752  ENDIF
753 
754  IF (idex.EQ. 11) then !electron
755  ide= 2
756  amin=0.511d-3
757  ELSEIF (idex.EQ.-11) then !positron
758  ide=-2
759  amin=0.511d-3
760  ELSEIF (idex.EQ. 13) then !mu+
761  ide= 2
762  amin=0.105659
763  ELSEIF (idex.EQ.-13) then !mu-
764  ide=-2
765  amin=0.105659
766  ELSEIF (idex.EQ. 1) then !d
767  ide= 4
768  amin=0.05
769  ELSEIF (idex.EQ.- 1) then !d~
770  ide=-4
771  amin=0.05
772  ELSEIF (idex.EQ. 2) then !u
773  ide= 3
774  amin=0.02
775  ELSEIF (idex.EQ.- 2) then !u~
776  ide=-3
777  amin=0.02
778  ELSEIF (idex.EQ. 3) then !s
779  ide= 4
780  amin=0.3
781  ELSEIF (idex.EQ.- 3) then !s~
782  ide=-4
783  amin=0.3
784  ELSEIF (idex.EQ. 4) then !c
785  ide= 3
786  amin=1.3
787  ELSEIF (idex.EQ.- 4) then !c~
788  ide=-3
789  amin=1.3
790  ELSEIF (idex.EQ. 5) then !b
791  ide= 4
792  amin=4.5
793  ELSEIF (idex.EQ.- 5) then !b~
794  ide=-4
795  amin=4.5
796  ELSEIF (idex.EQ. 12) then !nu_e
797  ide= 1
798  amin=0.1d-3
799  ELSEIF (idex.EQ.- 12) then !nu_e~
800  ide=-1
801  amin=0.1d-3
802  ELSEIF (idex.EQ. 14) then !nu_mu
803  ide= 1
804  amin=0.1d-3
805  ELSEIF (idex.EQ.- 14) then !nu_mu~
806  ide=-1
807  amin=0.1d-3
808  ELSEIF (idex.EQ. 16) then !nu_tau
809  ide= 1
810  amin=0.1d-3
811  ELSEIF (idex.EQ.- 16) then !nu_tau~
812  ide=-1
813  amin=0.1d-3
814 
815  ELSE
816  WRITE(*,*) 'INITWK: WRONG IDEX'
817  stop
818  ENDIF
819 
820 C ----------------------------------------------------------------------
821 C
822 C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
823 C
824 C called by : KORALZ
825 C ----------------------------------------------------------------------
826  itce=ide/iabs(ide)
827  jtce=(1-itce)/2
828  itcf=idf/iabs(idf)
829  jtcf=(1-itcf)/2
830  CALL t_givizo( ide, 1,aizor,qe,kdumm)
831  CALL t_givizo( ide,-1,aizol,qe,kdumm)
832  xupgi(1)=qe
833  xupgi(2)=qe
834  t3e = aizol+aizor
835  xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
836  xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
837  CALL t_givizo( idf, 1,aizor,qf,kolor)
838  CALL t_givizo( idf,-1,aizol,qf,kolor)
839  xupgf(1)=qf
840  xupgf(2)=qf
841  t3f = aizol+aizor
842  xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
843  xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
844 C
845  ndiag0=2
846  ndiaga=11
847  keya = 1
848  keyz = 1
849 C
850 C
851  RETURN
852  END
853 
854  SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
855 C ----------------------------------------------------------------------
856 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
857 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
858 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
859 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
860 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
861 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
862 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
863 C
864 C called by : EVENTE, EVENTM, FUNTIH, .....
865 C ----------------------------------------------------------------------
866  IMPLICIT REAL*8(a-h,o-z)
867 C
868  IF(idferm.EQ.0.OR.iabs(idferm).GT.4) GOTO 901
869  IF(iabs(ihelic).NE.1) GOTO 901
870  ih =ihelic
871  idtype =iabs(idferm)
872  ic =idferm/idtype
873  lepqua=int(idtype*0.4999999d0)
874  iupdow=idtype-2*lepqua-1
875  charge =(-iupdow+2d0/3d0*lepqua)*ic
876  sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
877  kolor=1+2*lepqua
878 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
879 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
880  RETURN
881  901 print *,' STOP IN GIVIZO: WRONG PARAMS.'
882  stop
883  END
884  SUBROUTINE phyfix(NSTOP,NSTART)
885  common/lujets/n,k(4000,5),p(4000,5),v(4000,5)
886  SAVE /lujets/
887 C NSTOP NSTART : when PHYTIA history ends and event starts.
888  nstop=0
889  nstart=1
890  DO i=1, n
891  IF(k(i,1).NE.21) THEN
892  nstop = i-1
893  nstart= i
894  GOTO 500
895  ENDIF
896  ENDDO
897  500 CONTINUE
898  END
899 
900 
901  SUBROUTINE taupi0(PI0,K)
902 C no initialization required. Must be called once after every:
903 C 1) CALL DEKAY(1+10,...)
904 C 2) CALL DEKAY(2+10,...)
905 C 3) CALL DEXAY(1,...)
906 C 4) CALL DEXAY(2,...)
907 C subroutine to decay originating from TAUOLA's taus:
908 C 1) etas (with CALL TAUETA(JAK))
909 C 2) later pi0's from taus.
910 C 3) extensions to other applications possible.
911 C this routine belongs to >tauola universal interface<, but uses
912 C routines from >tauola< utilities as well. 25.08.2005
913 C this is the hepevt class in old style. No d_h_ class pre-name
914 
915 C position of taus, must be defined by host program:
916  COMMON /taupos/ np1,np2
917 c
918  REAL PHOT1(4),PHOT2(4)
919  real*8 r,x(4),y(4),pi0(4)
920  INTEGER JEZELI(3),K
921  DATA jezeli /0,0,0/
922  SAVE jezeli
923 
924 ! random 3 vector on the sphere, masless
925  r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
926  CALL spherd(r,x)
927  x(4)=r
928  y(4)=r
929 
930  y(1)=-x(1)
931  y(2)=-x(2)
932  y(3)=-x(3)
933 ! boost to lab and to real*4
934  CALL bostdq(-1,pi0,x,x)
935  CALL bostdq(-1,pi0,y,y)
936  DO l=1,4
937  phot1(l)=x(l)
938  phot2(l)=y(l)
939  ENDDO
940 C to hepevt
941  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
942  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
943 
944 C
945  END
946  SUBROUTINE taueta(PETA,K)
947 C subroutine to decay etas's from taus.
948 C this routine belongs to tauola universal interface, but uses
949 C routines from tauola utilities. Just flat phase space, but 4 channels.
950 C it is called at the beginning of SUBR. TAUPI0(JAK)
951 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
952 C this is the hepevt class in old style. No d_h_ class pre-name
953 
954  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
955  * ,ampiz,ampi,amro,gamro,ama1,gama1
956  * ,amk,amkz,amkst,gamkst
957 *
958  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
959  * ,ampiz,ampi,amro,gamro,ama1,gama1
960  * ,amk,amkz,amkst,gamkst
961 
962 C position of taus, must be defined by host program:
963 c
964  REAL RRR(1),BRSUM(3), RR(2)
965  REAL PHOT1(4),PHOT2(4),PHOT3(4)
966  real*8 x(4), y(4), z(4)
967  REAL YM1,YM2,YM3
968  real*8 r,ru,peta(4),xm1,xm2,xm3,xlam
969  real*8 a,b,c
970  INTEGER K
971  xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
972 C position of decaying particle:
973 C DO L=1,4
974 C PETA(L)= phep(L,K) ! eta 4 momentum
975 C ENDDO
976 C eta cumulated branching ratios:
977  brsum(1)=0.389 ! gamma gamma
978  brsum(2)=brsum(1)+0.319 ! 3 pi0
979  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
980  CALL ranmar(rrr,1)
981 
982  IF (rrr(1).LT.brsum(1)) THEN ! gamma gamma channel exactly like pi0
983 ! random 3 vector on the sphere, masless
984  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
985  CALL spherd(r,x)
986  x(4)=r
987  y(4)=r
988 
989  y(1)=-x(1)
990  y(2)=-x(2)
991  y(3)=-x(3)
992 ! boost to lab and to real*4
993  CALL bostdq(-1,peta,x,x)
994  CALL bostdq(-1,peta,y,y)
995  DO l=1,4
996  phot1(l)=x(l)
997  phot2(l)=y(l)
998  ENDDO
999 C to hepevt
1000  CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1001  CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1002  ELSE ! 3 body channels
1003  IF(rrr(1).LT.brsum(2)) THEN ! 3 pi0
1004  id1= 111
1005  id2= 111
1006  id3= 111
1007  xm1=ampiz ! masses
1008  xm2=ampiz
1009  xm3=ampiz
1010  ELSEIF(rrr(1).LT.brsum(3)) THEN ! pi+ pi- pi0
1011  id1= 211
1012  id2=-211
1013  id3= 111
1014  xm1=ampi ! masses
1015  xm2=ampi
1016  xm3=ampiz
1017  ELSE ! pi+ pi- gamma
1018  id1= 211
1019  id2=-211
1020  id3= 22
1021  xm1=ampi ! masses
1022  xm2=ampi
1023  xm3=0.0
1024  ENDIF
1025  7 CONTINUE ! we generate mass of the first pair:
1026  CALL ranmar(rr,2)
1027  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1028  amin=xm1+xm2
1029  amax=r-xm3
1030  am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1031 C weight for flat phase space
1032  wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1033  & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1034  & /r**2 /am2**2
1035  IF (rr(2).GT.wt) GOTO 7
1036 
1037  ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2 ! momenta of the
1038  ! first two products
1039  ! in the rest frame of that pair
1040  CALL spherd(ru,x)
1041  x(4)=sqrt(ru**2+xm1**2)
1042  y(4)=sqrt(ru**2+xm2**2)
1043 
1044  y(1)=-x(1)
1045  y(2)=-x(2)
1046  y(3)=-x(3)
1047 C generate momentum of that pair in rest frame of eta:
1048  ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1049  CALL spherd(ru,z)
1050  z(4)=sqrt(ru**2+am2**2)
1051 C and boost first two decay products to rest frame of eta.
1052  CALL bostdq(-1,z,x,x)
1053  CALL bostdq(-1,z,y,y)
1054 C redefine Z(4) to 4-momentum of the last decay product:
1055  z(1)=-z(1)
1056  z(2)=-z(2)
1057  z(3)=-z(3)
1058  z(4)=sqrt(ru**2+xm3**2)
1059 C boost all to lab and move to real*4; also masses
1060  CALL bostdq(-1,peta,x,x)
1061  CALL bostdq(-1,peta,y,y)
1062  CALL bostdq(-1,peta,z,z)
1063  DO l=1,4
1064  phot1(l)=x(l)
1065  phot2(l)=y(l)
1066  phot3(l)=z(l)
1067  ENDDO
1068  ym1=xm1
1069  ym2=xm2
1070  ym3=xm3
1071 C to hepevt
1072  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1073  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1074  CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1075  ENDIF
1076 
1077 
1078 C
1079  END
1080  SUBROUTINE tauk0s(PETA,K)
1081 C subroutine to decay K0S's from taus.
1082 C this routine belongs to tauola universal interface, but uses
1083 C routines from tauola utilities. Just flat phase space, but 4 channels.
1084 C it is called at the beginning of SUBR. TAUPI0(JAK)
1085 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1086 
1087  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1088  * ,ampiz,ampi,amro,gamro,ama1,gama1
1089  * ,amk,amkz,amkst,gamkst
1090 *
1091  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1092  * ,ampiz,ampi,amro,gamro,ama1,gama1
1093  * ,amk,amkz,amkst,gamkst
1094 
1095 C position of taus, must be defined by host program:
1096  COMMON /taupos/ np1,np2
1097 c
1098  REAL RRR(1),BRSUM(3)
1099  REAL PHOT1(4),PHOT2(4)
1100  real*8 x(4), y(4)
1101  REAL YM1,YM2
1102  real*8 r,peta(4),xm1,xm2,xlam
1103  real*8 a,b,c
1104  INTEGER K
1105  xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1106 C position of decaying particle:
1107 
1108 
1109 ! DO L=1,4
1110 ! PETA(L)= phep(L,K) ! K0S 4 momentum (this is cloned from eta decay)
1111 ! ENDDO
1112 C K0S cumulated branching ratios:
1113  brsum(1)=0.313 ! 2 PI0
1114  brsum(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI-
1115  brsum(3)=brsum(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma
1116  CALL ranmar(rrr,1)
1117 
1118  IF(rrr(1).LT.brsum(1)) THEN ! 2 pi0
1119  id1= 111
1120  id2= 111
1121  xm1=ampiz ! masses
1122  xm2=ampiz
1123  ELSEIF(rrr(1).LT.brsum(2)) THEN ! pi+ pi-
1124  id1= 211
1125  id2=-211
1126  xm1=ampi ! masses
1127  xm2=ampi
1128  ELSE ! gamma gamma unused !!!
1129  id1= 22
1130  id2= 22
1131  xm1= 0.0 ! masses
1132  xm2= 0.0
1133  ENDIF
1134 
1135 ! random 3 vector on the sphere, of equal mass !!
1136  r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1137  r4=r
1138  r=sqrt(abs(r**2-xm1**2))
1139  CALL spherd(r,x)
1140  x(4)=r4
1141  y(4)=r4
1142 
1143  y(1)=-x(1)
1144  y(2)=-x(2)
1145  y(3)=-x(3)
1146 ! boost to lab and to real*4
1147  CALL bostdq(-1,peta,x,x)
1148  CALL bostdq(-1,peta,y,y)
1149  DO l=1,4
1150  phot1(l)=x(l)
1151  phot2(l)=y(l)
1152  ENDDO
1153 
1154  ym1=xm1
1155  ym2=xm2
1156 C to hepevt
1157  CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1158  CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1159 
1160 C
1161  END
1162 
1163  subroutine bostdq(idir,vv,pp,q)
1164 * *******************************
1165 c Boost along arbitrary vector v (see eg. J.D. Jacson, Classical
1166 c Electrodynamics).
1167 c Four-vector pp is boosted from an actual frame to the rest frame
1168 c of the four-vector v (for idir=1) or back (for idir=-1).
1169 c q is a resulting four-vector.
1170 c Note: v must be time-like, pp may be arbitrary.
1171 c
1172 c Written by: Wieslaw Placzek date: 22.07.1994
1173 c Last update: 3/29/95 by: M.S.
1174 c
1175  implicit DOUBLE PRECISION (a-h,o-z)
1176  parameter(nout=6)
1177  DOUBLE PRECISION v(4),p(4),q(4),pp(4),vv(4)
1178  save
1179 !
1180  do 1 i=1,4
1181  v(i)=vv(i)
1182  1 p(i)=pp(i)
1183  amv=(v(4)**2-v(1)**2-v(2)**2-v(3)**2)
1184  if (amv.le.0d0) then
1185  write(6,*) 'bosstv: warning amv**2=',amv
1186  endif
1187  amv=sqrt(abs(amv))
1188  if (idir.eq.-1) then
1189  q(4)=( p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4))/amv
1190  wsp =(q(4)+p(4))/(v(4)+amv)
1191  elseif (idir.eq.1) then
1192  q(4)=(-p(1)*v(1)-p(2)*v(2)-p(3)*v(3)+p(4)*v(4))/amv
1193  wsp =-(q(4)+p(4))/(v(4)+amv)
1194  else
1195  write(nout,*)' >>> boostv: wrong value of idir = ',idir
1196  endif
1197  q(1)=p(1)+wsp*v(1)
1198  q(2)=p(2)+wsp*v(2)
1199  q(3)=p(3)+wsp*v(3)
1200  end
1201 
1202 
1203 
1204 
Generated by doxygen 1.9.1