1 /* copyright(c) 1991-2012 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public
6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public
15 license along with the gnu c library;
if not, see
16 <http://www.gnu.org/licenses/>. */
19 /* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27 /* we
do support the iec 559 math functionality,
real and complex. */
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
32 /* we
do not support c11 <threads.h>. */
35 c *********************
37 c **********************************************************************
39 c *********tauola library: version 2.7 ******** *
40 c **************august 1995****************** *
41 c ** authors: s.jadach, z.was ***** *
42 c ** r. decker, m. jezabek, j.h.kuehn, ***** *
43 c ********available from: wasm at cernvm ****** *
44 c *******published in comp. phys. comm.******** *
45 c *** preprint cern-th-5856 september 1990 **** *
46 c *** preprint cern-th-6195 october 1991 **** *
47 c *** preprint cern-th-6793 november 1992 **** *
48 c **********************************************************************
50 c ----------------------------------------------------------------------
52 c chooses decay mode according to list of branching ratios
63 c ----------------------------------------------------------------------
64 COMMON / taubra / gamprt(30),jlist(30),nchan
68 IF(nchan.LE.0.OR.nchan.GT.30) goto 902
75 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
80 9020
FORMAT(
' ----- JAKER: WRONG NCHAN')
83 SUBROUTINE dekay(KTO,HX)
84 c ***********************
85 c this dekay is in spirit of the
'DECAY' which
86 c was included in koral-b program, comp. phys. commun.
87 c vol. 36 (1985) 191, see comments on general philosophy there.
88 c kto=0 initialisation(obligatory)
89 c kto=1,11 denotes tau+ and kto=2,12 tau-
90 c dekay(1,h) and dekay(2,h) is called internally by mc generator.
91 c h denotes the polarimetric vector, used by the host
PROGRAM for
92 c calculation of the spin weight.
93 c user may optionally CALL dekay(11,h) dekay(12,h) in order
94 c to transform decay products to cms and
WRITE lund record in /lujets/.
95 c kto=100, print final report(
OPTIONAL).
97 c jak=1 electron decay
105 c jak=0 inclusive: jak=1,2,3,4,5,6,7,8
108 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
110 COMMON /taupos/ np1,np2
111 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
112 REAL*4 gampmc ,gamper
113 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
114 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
116 CHARACTER names(nmode)*31
117 COMMON / inout / inut,iout
118 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4),hdum(4)
124 c initialisation or reinitialisation
125 c first or second tau positions in hepevt as in koralb/z
129 IF (iwarm.EQ.1) x=5/(iwarm-1)
131 WRITE(iout,7001) jak1,jak2
135 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
136 CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
137 CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
138 CALL dadmpi(-1,idum,hdum,pdum1,pdum2)
139 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
140 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
141 CALL dadmkk(-1,idum,hdum,pdum1,pdum2)
142 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
143 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
149 ELSEIF(kto.EQ.1)
THEN
150 c =====================
151 c decay of tau+ in the tau rest frame
153 IF(iwarm.EQ.0) goto 902
155 c ajwmod to change brs depending on sign:
157 CALL dekay1(0,h,isgn)
158 ELSEIF(kto.EQ.2)
THEN
159 c =================================
160 c decay of tau- in the tau rest frame
162 IF(iwarm.EQ.0) goto 902
164 c ajwmod to change brs depending on sign:
166 CALL dekay2(0,h,isgn)
167 ELSEIF(kto.EQ.11)
THEN
168 c ======================
169 c rest of decay procedure for accepted tau+ decay
172 CALL dekay1(1,h,isgn)
173 ELSEIF(kto.EQ.12)
THEN
174 c ======================
175 c rest of decay procedure for accepted tau- decay
178 CALL dekay2(1,h,isgn)
179 ELSEIF(kto.EQ.100)
THEN
180 c =======================
181 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
182 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
183 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
184 CALL dadmpi( 1,idum,hdum,pdum1,pdum2)
185 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
186 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
187 CALL dadmkk( 1,idum,hdum,pdum1,pdum2)
188 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
189 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
190 WRITE(iout,7010) nev1,nev2,nevtot
191 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
193 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
204 7001
FORMAT(///1x,15(5h*****)
205 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
206 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
207 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
208 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
209 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
210 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
211 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
212 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
213 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
214 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
215 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
216 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
217 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
218 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
219 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
221 7010
FORMAT(///1x,15(5h*****)
222 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
223 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
224 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
225 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
226 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
227 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
228 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
229 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
230 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
231 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
232 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
233 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
234 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
236 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
238 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
239 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
240 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
241 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
242 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
243 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
244 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
246 $ ,i10,2f12.7,a31 ,8x,1h*)
248 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
249 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
252 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
255 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
258 SUBROUTINE dekay1(IMOD,HH,ISGN)
259 c *******************************
260 c this routine simulates tau+ decay
261 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
262 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
263 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
264 REAL*4 gampmc ,gamper
266 REAL hv(4),pnu(4),ppi(4)
267 REAL pwb(4),pmu(4),pnm(4)
268 REAL prho(4),pic(4),piz(4)
269 REAL paa(4),pim1(4),pim2(4),pipl(4)
276 IF(jak1.EQ.-1)
RETURN
281 IF(jak1.EQ.0) CALL jaker(jak)
283 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
284 ELSEIF(jak.EQ.2)
THEN
285 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
286 ELSEIF(jak.EQ.3)
THEN
287 CALL dadmpi(0, isgn,hv,ppi,pnu)
288 ELSEIF(jak.EQ.4)
THEN
289 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
290 ELSEIF(jak.EQ.5)
THEN
291 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
292 ELSEIF(jak.EQ.6)
THEN
293 CALL dadmkk(0, isgn,hv,pkk,pnu)
294 ELSEIF(jak.EQ.7)
THEN
295 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
297 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
303 ELSEIF(imd.EQ.1)
THEN
304 c =====================
307 nevdec(jak)=nevdec(jak)+1
312 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
313 CALL dwrph(ktom,phot)
317 ELSEIF(jak.EQ.2)
THEN
318 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
319 CALL dwrph(ktom,phot)
323 ELSEIF(jak.EQ.3)
THEN
324 CALL dwlupi(1,isgn,ppi,pnu)
328 ELSEIF(jak.EQ.4)
THEN
329 CALL dwluro(1,isgn,pnu,prho,pic,piz)
333 ELSEIF(jak.EQ.5)
THEN
334 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
337 ELSEIF(jak.EQ.6)
THEN
338 CALL dwlukk(1,isgn,pkk,pnu)
341 ELSEIF(jak.EQ.7)
THEN
342 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
347 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
355 SUBROUTINE dekay2(IMOD,HH,ISGN)
356 c *******************************
357 c this routine simulates tau- decay
358 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
359 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
360 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
361 REAL*4 gampmc ,gamper
363 REAL hv(4),pnu(4),ppi(4)
364 REAL pwb(4),pmu(4),pnm(4)
365 REAL prho(4),pic(4),piz(4)
366 REAL paa(4),pim1(4),pim2(4),pipl(4)
373 IF(jak2.EQ.-1)
RETURN
378 IF(jak2.EQ.0) CALL jaker(jak)
380 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
381 ELSEIF(jak.EQ.2)
THEN
382 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
383 ELSEIF(jak.EQ.3)
THEN
384 CALL dadmpi(0, isgn,hv,ppi,pnu)
385 ELSEIF(jak.EQ.4)
THEN
386 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
387 ELSEIF(jak.EQ.5)
THEN
388 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
389 ELSEIF(jak.EQ.6)
THEN
390 CALL dadmkk(0, isgn,hv,pkk,pnu)
391 ELSEIF(jak.EQ.7)
THEN
392 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
394 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
399 ELSEIF(imd.EQ.1)
THEN
400 c =====================
403 nevdec(jak)=nevdec(jak)+1
408 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
409 CALL dwrph(ktom,phot)
413 ELSEIF(jak.EQ.2)
THEN
414 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
415 CALL dwrph(ktom,phot)
419 ELSEIF(jak.EQ.3)
THEN
420 CALL dwlupi(2,isgn,ppi,pnu)
424 ELSEIF(jak.EQ.4)
THEN
425 CALL dwluro(2,isgn,pnu,prho,pic,piz)
429 ELSEIF(jak.EQ.5)
THEN
430 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
433 ELSEIF(jak.EQ.6)
THEN
434 CALL dwlukk(2,isgn,pkk,pnu)
437 ELSEIF(jak.EQ.7)
THEN
438 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
443 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
451 SUBROUTINE dexay(KTO,POL)
452 c ----------------------------------------------------------------------
453 c this
'DEXAY' is a routine which generates decay of the single
454 c polarized tau, pol is a polarization vector(not a polarimeter
455 c vector as in dekay) of the tau and it is an input
PARAMETER.
456 c kto=0 initialisation(obligatory)
457 c kto=1 denotes tau+ and kto=2 tau-
458 c dexay(1,pol) and dexay(2,pol) are called internally by mc generator.
459 c decay products are transformed readily
460 c to cms and writen in the lund record in /lujets/
461 c kto=100, print final report(
OPTIONAL).
464 c ----------------------------------------------------------------------
465 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
466 REAL*4 gampmc ,gamper
467 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
469 COMMON /taupos/ np1,np2
470 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
471 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
473 CHARACTER names(nmode)*31
474 COMMON / inout / inut,iout
476 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
485 c initialisation or reinitialisation
486 c first or second tau positions in hepevt as in koralb/z
490 WRITE(iout, 7001) jak1,jak2
494 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
495 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
496 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
497 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
498 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
499 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
500 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
501 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
502 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
508 ELSEIF(kto.EQ.1)
THEN
509 c =====================
510 c decay of tau+ in the tau rest frame
513 IF(iwarm.EQ.0) goto 902
515 cam CALL dexay1(pol,isgn)
516 CALL dexay1(kto,jak1,jakp,pol,isgn)
517 ELSEIF(kto.EQ.2)
THEN
518 c =================================
519 c decay of tau- in the tau rest frame
522 IF(iwarm.EQ.0) goto 902
523 isgn=-idff/iabs(idff)
524 cam CALL dexay2(pol,isgn)
525 CALL dexay1(kto,jak2,jakm,pol,isgn)
526 ELSEIF(kto.EQ.100)
THEN
527 c =======================
528 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
529 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
530 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
531 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
532 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
533 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
534 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
535 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
536 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
537 WRITE(iout,7010) nev1,nev2,nevtot
538 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
540 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
547 7001
FORMAT(///1x,15(5h*****)
548 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
549 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
550 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
551 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
552 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
553 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
554 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
555 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
556 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
557 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
558 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
559 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
560 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
561 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
562 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
563 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
564 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
566 chbu
format 7010 had more than 19 continuation lines
568 7010
FORMAT(///1x,15(5h*****)
569 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
570 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
571 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
572 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
573 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
574 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
575 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
576 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
577 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
578 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
579 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
580 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
581 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
583 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
585 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
586 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
587 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
588 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
589 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
590 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
591 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
593 $ ,i10,2f12.7,a31 ,8x,1h*)
595 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
596 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
598 902
WRITE(iout, 9020)
599 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
601 910
WRITE(iout, 9100)
602 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
605 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
606 c ---------------------------------------------------------------------
607 c this routine simulates tau+- decay
610 c ---------------------------------------------------------------------
611 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
612 REAL*4 gampmc ,gamper
613 COMMON / inout / inut,iout
616 REAL prho(4),pic(4),piz(4)
617 REAL pwb(4),pmu(4),pnm(4)
618 REAL paa(4),pim1(4),pim2(4),pipl(4)
624 IF(jakin.EQ.-1)
RETURN
631 IF(jak.EQ.0) CALL jaker(jak)
634 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
635 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
636 CALL dwrph(kto,phot )
637 ELSEIF(jak.EQ.2)
THEN
638 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
639 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
640 CALL dwrph(kto,phot )
641 ELSEIF(jak.EQ.3)
THEN
642 CALL dexpi(0, isgn,polar,ppi,pnu)
643 CALL dwlupi(kto,isgn,ppi,pnu)
644 ELSEIF(jak.EQ.4)
THEN
645 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
646 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
647 ELSEIF(jak.EQ.5)
THEN
648 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
649 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
650 ELSEIF(jak.EQ.6)
THEN
651 CALL dexkk(0, isgn,polar,pkk,pnu)
652 CALL dwlukk(kto,isgn,pkk,pnu)
653 ELSEIF(jak.EQ.7)
THEN
654 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
655 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
658 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
659 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
661 nevdec(jak)=nevdec(jak)+1
663 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
664 c ----------------------------------------------------------------------
665 c this simulates tau decay in tau rest frame
666 c into electron and two neutrinos
668 c called by : dexay,dexay1
669 c ----------------------------------------------------------------------
670 REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
674 c ===================
676 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
677 cc CALL hbook1(813,
'WEIGHT DISTRIBUTION DEXEL $',100,0,2)
679 ELSEIF(mode.EQ. 0)
THEN
680 c =======================
682 IF(iwarm.EQ.0) goto 902
683 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
684 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
685 cc CALL hfill(813,wt)
687 IF(rn(1).GT.wt) goto 300
689 ELSEIF(mode.EQ. 1)
THEN
690 c =======================
691 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
697 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
700 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
701 c ----------------------------------------------------------------------
702 c this simulates tau decay in its rest frame
703 c into muon and two neutrinos
704 c output four momenta: pnu tauneutrino,
708 c ----------------------------------------------------------------------
709 COMMON / inout / inut,iout
710 REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
714 c ===================
716 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
717 cc CALL hbook1(814,
'WEIGHT DISTRIBUTION DEXMU $',100,0,2)
719 ELSEIF(mode.EQ. 0)
THEN
720 c =======================
722 IF(iwarm.EQ.0) goto 902
723 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
724 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
725 cc CALL hfill(814,wt)
727 IF(rn(1).GT.wt) goto 300
729 ELSEIF(mode.EQ. 1)
THEN
730 c =======================
731 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
736 902
WRITE(iout, 9020)
737 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
740 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
741 c ----------------------------------------------------------------------
743 c called by : dexel,(dekay,dekay1)
744 c ----------------------------------------------------------------------
745 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
746 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
747 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
748 * ,ampiz,ampi,amro,gamro,ama1,gama1
749 * ,amk,amkz,amkst,gamkst
751 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
752 * ,ampiz,ampi,amro,gamro,ama1,gama1
753 * ,amk,amkz,amkst,gamkst
754 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
755 REAL*4 gampmc ,gamper
757 COMMON / inout / inut,iout
758 REAL hhv(4),hv(4),pwb(4),pnu(4),q1(4),q2(4)
759 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
762 DATA pi /3.141592653589793238462643/
766 c ===================
775 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
776 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
778 cc CALL hbook1(803,
'WEIGHT DISTRIBUTION DADMEL $',100,0,2)
780 ELSEIF(mode.EQ. 0)
THEN
781 c =======================
783 IF(iwarm.EQ.0) goto 902
785 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
786 cc CALL hfill(803,wt/wtmax)
791 IF(wt.GT.wtmax) nevovr=nevovr+1
792 IF(rn*wtmax.GT.wt) goto 300
793 c rotations to basic tau rest frame
799 CALL rotor2(thet,pnu,pnu)
800 CALL rotor3( phi,pnu,pnu)
801 CALL rotor2(thet,pwb,pwb)
802 CALL rotor3( phi,pwb,pwb)
803 CALL rotor2(thet,q1,q1)
804 CALL rotor3( phi,q1,q1)
805 CALL rotor2(thet,q2,q2)
806 CALL rotor3( phi,q2,q2)
807 CALL rotor2(thet,hv,hv)
808 CALL rotor3( phi,hv,hv)
809 CALL rotor2(thet,phx,phx)
810 CALL rotor3( phi,phx,phx)
812 44 hhv(i)=-isgn*hv(i)
815 ELSEIF(mode.EQ. 1)
THEN
816 c =======================
817 IF(nevraw.EQ.0)
RETURN
818 pargam=swt/float(nevraw+1)
820 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
822 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
830 7010
FORMAT(///1x,15(5h*****)
831 $ /,
' *', 25x,
'******** DADMEL FINAL REPORT ******** ',9x,1h*
832 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
833 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
834 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
835 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
836 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
837 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
838 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
839 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
841 902
WRITE(iout, 9020)
842 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
845 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
846 c ----------------------------------------------------------------------
847 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
848 * ,ampiz,ampi,amro,gamro,ama1,gama1
849 * ,amk,amkz,amkst,gamkst
851 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
852 * ,ampiz,ampi,amro,gamro,ama1,gama1
853 * ,amk,amkz,amkst,gamkst
854 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
855 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
856 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
857 REAL*4 gampmc ,gamper
858 COMMON / inout / inut,iout
860 REAL hhv(4),hv(4),pnu(4),pwb(4),q1(4),q2(4)
861 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
864 DATA pi /3.141592653589793238462643/
868 c ===================
877 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
878 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
880 cc CALL hbook1(802,
'WEIGHT DISTRIBUTION DADMMU $',100,0,2)
882 ELSEIF(mode.EQ. 0)
THEN
883 c =======================
885 IF(iwarm.EQ.0) goto 902
887 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
888 cc CALL hfill(802,wt/wtmax)
893 IF(wt.GT.wtmax) nevovr=nevovr+1
894 IF(rn*wtmax.GT.wt) goto 300
895 c rotations to basic tau rest frame
899 CALL rotor2(thet,pnu,pnu)
900 CALL rotor3( phi,pnu,pnu)
901 CALL rotor2(thet,pwb,pwb)
902 CALL rotor3( phi,pwb,pwb)
903 CALL rotor2(thet,q1,q1)
904 CALL rotor3( phi,q1,q1)
905 CALL rotor2(thet,q2,q2)
906 CALL rotor3( phi,q2,q2)
907 CALL rotor2(thet,hv,hv)
908 CALL rotor3( phi,hv,hv)
909 CALL rotor2(thet,phx,phx)
910 CALL rotor3( phi,phx,phx)
912 44 hhv(i)=-isgn*hv(i)
915 ELSEIF(mode.EQ. 1)
THEN
916 c =======================
917 IF(nevraw.EQ.0)
RETURN
918 pargam=swt/float(nevraw+1)
920 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
922 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
930 7010
FORMAT(///1x,15(5h*****)
931 $ /,
' *', 25x,
'******** DADMMU FINAL REPORT ******** ',9x,1h*
932 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
933 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
934 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
935 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
936 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
937 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
938 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
939 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
941 902
WRITE(iout, 9020)
942 9020
FORMAT(
' ----- DADMMU: LACK OF INITIALISATION')
945 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
946 c xnx,xna was flipped in parameters of dphsel and dphsmu
947 c *********************************************************************
948 c * electron decay mode *
949 c *********************************************************************
951 REAL*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
952 REAL*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
955 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
966 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
967 c xnx,xna was flipped in parameters of dphsel and dphsmu
968 c *********************************************************************
969 c * muon decay mode *
970 c *********************************************************************
972 REAL*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
973 REAL*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
976 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
987 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
988 IMPLICIT REAL*8 (a-h,o-z)
989 c ----------------------------------------------------------------------
990 * it simulates e,mu channels of tau decay in its rest frame with
991 * qed order alpha corrections
992 c ----------------------------------------------------------------------
993 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
994 * ,ampiz,ampi,amro,gamro,ama1,gama1
995 * ,amk,amkz,amkst,gamkst
997 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
998 * ,ampiz,ampi,amro,gamro,ama1,gama1
999 * ,amk,amkz,amkst,gamkst
1000 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1001 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1002 COMMON / inout / inut,iout
1003 COMMON / taurad / xk0dec,itdkrc
1005 REAL*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1009 DATA pi /3.141592653589793238462643d0/
1010 c ajwmod to satisfy compiler, comment out this unused function.
1011 c amro, gamro is only a
PARAMETER for geting hight efficiency
1013 c three body phase space normalised as in bjorken-drell
1014 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
1015 phspac=1./2**17/pi**8
1025 IF (ielmu.EQ.1)
THEN
1032 IF ( itdkrc.EQ.0) prhard=0d0
1034 IF(prsoft.LT.0.1)
THEN
1035 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1040 ihard=(rr5.GT.prsoft)
1042 c tau decay to
'TAU+photon'
1044 ams1=(amu+amnuta)**2
1047 xl1=log(xk1/2/xk0dec)
1052 phspac=phspac*ams2*xl1*xk
1053 phspac=phspac/prhard
1056 phspac=phspac*2**6*pi**3
1057 phspac=phspac/prsoft
1059 c mass of neutrina system
1066 am2sq=ams1+ rr2*(ams2-ams1)
1068 phspac=phspac*(ams2-ams1)
1069 * neutrina rest frame, define xn and xa
1070 enq1=(am2sq+amnuta**2)/(2*am2)
1071 enq2=(am2sq-amnuta**2)/(2*am2)
1072 ppi= enq1**2-amnuta**2
1073 pppi=sqrt(abs(enq1**2-amnuta**2))
1074 phspac=phspac*(4*pi)*(2*pppi/am2)
1075 * nu tau in nunu rest frame
1076 CALL spherd(pppi,xn)
1078 * nu light in nunu rest frame
1082 * tau-prim rest frame, define qp(muon
1086 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1087 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1088 ppi = pr(4)**2-am2**2
1092 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1094 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1095 * neutrina boosted from their frame to tau-prim rest frame
1096 exe=(pr(4)+pr(3))/am2
1097 CALL bostd3(exe,xn,xn)
1098 CALL bostd3(exe,xa,xa)
1102 eps=4*(amu/amtax)**2
1103 xl1=log((2+eps)/eps)
1105 eta =exp(xl1*rr3+xl0)
1108 phspac=phspac*xl1/2*eta
1110 CALL rotpox(thet,phi,xn)
1111 CALL rotpox(thet,phi,xa)
1112 CALL rotpox(thet,phi,qp)
1113 CALL rotpox(thet,phi,pr)
1115 * now to the tau rest frame, define tau-prim and gamma momenta
1119 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1120 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1121 ppi = paa(4)**2-am3**2
1122 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1128 * all momenta boosted from tau-prim rest frame to tau rest frame
1129 * z-axis antiparallel to photon momentum
1130 exe=(paa(4)+paa(3))/am3
1131 CALL bostd3(exe,xn,xn)
1132 CALL bostd3(exe,xa,xa)
1133 CALL bostd3(exe,qp,qp)
1134 CALL bostd3(exe,pr,pr)
1136 thet =acos(-1.+2*rr3)
1138 CALL rotpox(thet,phi,xn)
1139 CALL rotpox(thet,phi,xa)
1140 CALL rotpox(thet,phi,qp)
1141 CALL rotpox(thet,phi,pr)
1143 * now to the tau rest frame, define tau-prim and gamma momenta
1155 c partial width consists of phase space and amplitude
1156 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1157 dgamt=1/(2.*amtax)*amplit*phspac
1159 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1160 IMPLICIT REAL*8 (a-h,o-z)
1161 c ----------------------------------------------------------------------
1162 c it calculates matrix element for the
1163 c tau --> mu(e) nu nubar decay mode
1164 c including complete order alpha qed corrections.
1165 c ----------------------------------------------------------------------
1166 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1167 * ,ampiz,ampi,amro,gamro,ama1,gama1
1168 * ,amk,amkz,amkst,gamkst
1170 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1171 * ,ampiz,ampi,amro,gamro,ama1,gama1
1172 * ,amk,amkz,amkst,gamkst
1173 REAL*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1177 IF(xk(4).LT.0.1d0*ak0)
THEN
1178 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1180 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1184 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1186 c **********************************************************************
1187 c
REAL photon matrix element squared *
1189 c hv- polarimetric four-vector of tau *
1190 c qp,xn,xa,xk - 4-momenta of electron(muon), nu, nubar and photon *
1191 c all four-vectors in tau rest frame(in gev) *
1192 c ak0 - infrared cutoff, minimal energy of hard photons(gev) *
1193 c sqm2 - value for s=0 *
1194 c see eqs. (2.9)-(2.10) from cjk( nucl.phys.b(1991) ) *
1195 c **********************************************************************
1197 IMPLICIT REAL*8(a-h,o-z)
1198 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1199 * ,ampiz,ampi,amro,gamro,ama1,gama1
1200 * ,amk,amkz,amkst,gamkst
1202 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1203 * ,ampiz,ampi,amro,gamro,ama1,gama1
1204 * ,amk,amkz,amkst,gamkst
1205 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1206 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1207 COMMON / qedprm /alfinv,alfpi,xk0
1208 REAL*8 alfinv,alfpi,xk0
1209 REAL*8 qp(4),xn(4),xa(4),xk(4)
1212 REAL*8 s0(3),rxa(3),rxk(3),rqp(3)
1213 DATA pi /3.141592653589793238462643d0/
1219 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1221 c scalar products of four-momenta
1227 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1228 c rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1229 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1230 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1232 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1233 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1234 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1235 c xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1236 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1237 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1247 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1249 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1250 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1251 const4=256*pi/alphai*gf**2
1252 IF (itdkrc.EQ.0) const4=0d0
1255 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1257 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1258 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1259 5 hv(i)=s0(i)/s1-1.d0
1262 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1264 c **********************************************************************
1265 c born +virtual+soft photon matrix element**2 o(alpha) *
1267 c hv- polarimetric four-vector of tau *
1268 c qp,xn,xa - four-momenta of electron(muon), nu and nubar in gev *
1269 c all four-vectors in tau rest frame *
1270 c ak0 - infrared cutoff, minimal energy of hard photons *
1271 c thb - value for s=0 *
1272 c see eqs. (2.2),(2.4)-(2.5) from cjk(nucl.phys.b351(1991)70 *
1273 c and(c.2) from jk(nucl.phys.b320(1991)20 ) *
1274 c **********************************************************************
1276 IMPLICIT REAL*8(a-h,o-z)
1277 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1278 * ,ampiz,ampi,amro,gamro,ama1,gama1
1279 * ,amk,amkz,amkst,gamkst
1281 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1282 * ,ampiz,ampi,amro,gamro,ama1,gama1
1283 * ,amk,amkz,amkst,gamkst
1284 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1285 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1286 COMMON / qedprm /alfinv,alfpi,xk0
1287 REAL*8 alfinv,alfpi,xk0
1288 dimension qp(4),xn(4),xa(4)
1291 REAL*8 rxa(3),rxn(3),rqp(3)
1292 REAL*8 bornpl(3),am3pol(3),xm3pol(3)
1293 DATA pi /3.141592653589793238462643d0/
1306 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1307 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1308 c rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1309 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1311 c quasi two-body variables
1313 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1315 w0=(xn(4)+xa(4))/tmass
1327 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1328 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1329 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1330 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1331 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1332 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1334 c scalar products of four-momenta
1335 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1336 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1337 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1341 c decay differential width without and with polarization
1342 const3=1/(2*alphai*pi)*64*gf**2
1343 IF (itdkrc.EQ.0) const3=0d0
1344 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1345 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1347 c v-a and v+a couplings, but in the born part only
1348 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1349 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1350 born= 32*(gfermi**2/2.)*brak
1352 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1353 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1354 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1355 am3pol(i)=xm3pol(i)*const3
1356 c v-a and v+a couplings, but in the born part only
1358 & (gv+ga)**2*tmass*xnxa*qp(i)
1359 & -(gv-ga)**2*tmass*qpxn*xa(i)
1360 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1361 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1363 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1365 IF (thb/born.LT.0.1d0)
THEN
1366 print *,
'ERROR IN THB, THB/BORN=',thb/born
1371 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1372 c ----------------------------------------------------------------------
1373 c tau decay into pion and tau-neutrino
1375 c output four momenta: pnu tauneutrino,
1377 c ----------------------------------------------------------------------
1378 REAL pol(4),hv(4),pnu(4),ppi(4),rn(1)
1381 c ===================
1382 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1383 cc CALL hbook1(815,
'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1385 ELSEIF(mode.EQ. 0)
THEN
1386 c =======================
1388 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1389 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1390 cc CALL hfill(815,wt)
1392 IF(rn(1).GT.wt) goto 300
1394 ELSEIF(mode.EQ. 1)
THEN
1395 c =======================
1396 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1402 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1403 c ----------------------------------------------------------------------
1404 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1405 * ,ampiz,ampi,amro,gamro,ama1,gama1
1406 * ,amk,amkz,amkst,gamkst
1408 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1409 * ,ampiz,ampi,amro,gamro,ama1,gama1
1410 * ,amk,amkz,amkst,gamkst
1411 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1412 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1413 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1414 REAL*4 gampmc ,gamper
1415 COMMON / inout / inut,iout
1416 REAL ppi(4),pnu(4),hv(4)
1417 DATA pi /3.141592653589793238462643/
1420 c ===================
1422 ELSEIF(mode.EQ. 0)
THEN
1423 c =======================
1425 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1426 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1427 xpi= sqrt(epi**2-ampi**2)
1429 CALL sphera(xpi,ppi)
1431 c tau-neutrino momentum
1437 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1438 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1439 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1441 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1444 ELSEIF(mode.EQ. 1)
THEN
1445 c =======================
1446 IF(nevtot.EQ.0)
RETURN
1448 c gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1449 c * (brak/amtau**4)**2
1450 czw 7.02.93 here was an error affecting non standard model
1451 c configurations only
1452 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1454 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1455 $ -4*ampi**2*amnuta**2 )/amtau**2
1458 WRITE(iout, 7010) nevtot,gamm,rat,error
1461 cam nevdec(3)=nevtot
1465 7010
FORMAT(///1x,15(5h*****)
1466 $ /,
' *', 25x,
'******** DADMPI FINAL REPORT ******** ',9x,1h*
1467 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1468 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1469 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1470 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1471 $ /,1x,15(5h*****)/)
1473 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1474 c ----------------------------------------------------------------------
1475 c this simulates tau decay in tau rest frame
1476 c into nu rho, next rho decays into pion pair.
1477 c output four momenta: pnu tauneutrino,
1481 c ----------------------------------------------------------------------
1482 COMMON / inout / inut,iout
1483 REAL pol(4),hv(4),pro(4),pnu(4),pic(4),piz(4),rn(1)
1487 c ===================
1489 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1490 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXRO $',100,0,2)
1491 cc CALL hbook1(916,
'ABS2 OF HV IN ROUTINE DEXRO $',100,0,2)
1493 ELSEIF(mode.EQ. 0)
THEN
1494 c =======================
1496 IF(iwarm.EQ.0) goto 902
1497 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1498 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1499 cc CALL hfill(816,wt)
1500 cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
1501 cc CALL hfill(916,xhelp)
1503 IF(rn(1).GT.wt) goto 300
1505 ELSEIF(mode.EQ. 1)
THEN
1506 c =======================
1507 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1513 902
WRITE(iout, 9020)
1514 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1517 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1518 c ----------------------------------------------------------------------
1519 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1520 * ,ampiz,ampi,amro,gamro,ama1,gama1
1521 * ,amk,amkz,amkst,gamkst
1523 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1524 * ,ampiz,ampi,amro,gamro,ama1,gama1
1525 * ,amk,amkz,amkst,gamkst
1526 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1527 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1528 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1529 REAL*4 gampmc ,gamper
1530 COMMON / inout / inut,iout
1532 REAL hv(4),pro(4),pnu(4),pic(4),piz(4)
1533 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
1536 DATA pi /3.141592653589793238462643/
1540 c ===================
1549 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1550 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1552 cc CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMRO $',100,0,2)
1555 ELSEIF(mode.EQ. 0)
THEN
1556 c =======================
1558 IF(iwarm.EQ.0) goto 902
1559 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1560 cc CALL hfill(801,wt/wtmax)
1566 IF(wt.GT.wtmax) nevovr=nevovr+1
1567 IF(rn*wtmax.GT.wt) goto 300
1568 c rotations to basic tau rest frame
1569 costhe=-1.+2.*rrr(2)
1572 CALL rotor2(thet,pnu,pnu)
1573 CALL rotor3( phi,pnu,pnu)
1574 CALL rotor2(thet,pro,pro)
1575 CALL rotor3( phi,pro,pro)
1576 CALL rotor2(thet,pic,pic)
1577 CALL rotor3( phi,pic,pic)
1578 CALL rotor2(thet,piz,piz)
1579 CALL rotor3( phi,piz,piz)
1580 CALL rotor2(thet,hv,hv)
1581 CALL rotor3( phi,hv,hv)
1583 44 hhv(i)=-isgn*hv(i)
1586 ELSEIF(mode.EQ. 1)
THEN
1587 c =======================
1588 IF(nevraw.EQ.0)
RETURN
1589 pargam=swt/float(nevraw+1)
1591 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1593 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1597 cam nevdec(4)=nevacc
1601 7003
FORMAT(///1x,15(5h*****)
1602 $ /,
' *', 25x,
'******** DADMRO INITIALISATION ********',9x,1h*
1603 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1604 $ /,1x,15(5h*****)/)
1605 7010
FORMAT(///1x,15(5h*****)
1606 $ /,
' *', 25x,
'******** DADMRO FINAL REPORT ******** ',9x,1h*
1607 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1608 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1609 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1610 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1611 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1612 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1613 $ /,1x,15(5h*****)/)
1614 902
WRITE(iout, 9020)
1615 9020
FORMAT(
' ----- DADMRO: LACK OF INITIALISATION')
1618 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1619 c ----------------------------------------------------------------------
1620 c it simulates rho decay in tau rest frame with
1621 c z-axis along rho momentum
1622 c ----------------------------------------------------------------------
1623 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1624 * ,ampiz,ampi,amro,gamro,ama1,gama1
1625 * ,amk,amkz,amkst,gamkst
1627 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1628 * ,ampiz,ampi,amro,gamro,ama1,gama1
1629 * ,amk,amkz,amkst,gamkst
1630 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1631 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1632 REAL hv(4),pt(4),pn(4),pr(4),pic(4),piz(4),qq(4),rr1(1)
1633 DATA pi /3.141592653589793238462643/
1636 c three body phase space normalised as in bjorken-drell
1637 phspac=1./2**11/pi**5
1643 c mass of(real/virtual) rho
1644 ams1=(ampi+ampiz)**2
1645 ams2=(amtau-amnuta)**2
1647 c amx2=ams1+ rr1*(ams2-ams1)
1649 c phspac=phspac*(ams2-ams1)
1650 c phase space with sampling for rho resonance
1651 alp1=atan((ams1-amro**2)/amro/gamro)
1652 alp2=atan((ams2-amro**2)/amro/gamro)
1656 alp=alp1+rr1(1)*(alp2-alp1)
1657 amx2=amro**2+amro*gamro*tan(alp)
1659 IF(amx.LT.2.*ampi) go to 100
1661 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1662 phspac=phspac*(alp2-alp1)
1664 c tau-neutrino momentum
1667 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1668 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1672 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1674 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1677 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1678 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1679 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1680 phspac=phspac*(4*pi)*(2*pppi/amx)
1681 c charged pi momentum in rho rest frame
1682 CALL sphera(pppi,pic)
1684 c neutral pi momentum in rho rest frame
1688 exe=(pr(4)+pr(3))/amx
1689 c pions boosted from rho rest frame to tau rest frame
1690 CALL bostr3(exe,pic,pic)
1691 CALL bostr3(exe,piz,piz)
1693 30 qq(i)=pic(i)-piz(i)
1696 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
1698 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
1699 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
1700 & +(gv**2-ga**2)*amtau*amnuta*qq2
1701 amplit=(gfermi*ccabib)**2*brak*2*fpirho(amx)
1702 dgamt=1/(2.*amtau)*amplit*phspac
1704 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1707 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1708 c ----------------------------------------------------------------------
1709 * this simulates tau decay in tau rest frame
1710 * into nu a1, next a1 decays into rho pi and finally rho into pi pi.
1711 * output four momenta: pnu tauneutrino,
1713 * pim1 pion minus(or pi0) 1 (for tau minus)
1714 * pim2 pion minus(or pi0) 2
1715 * pipl pion plus(or pi-)
1716 * (pipl,pim1) form a rho
1717 c ----------------------------------------------------------------------
1718 COMMON / inout / inut,iout
1719 REAL pol(4),hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4),rn(1)
1723 c ===================
1725 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1726 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
1728 ELSEIF(mode.EQ. 0)
THEN
1729 * =======================
1731 IF(iwarm.EQ.0) goto 902
1732 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1733 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1734 cc CALL hfill(816,wt)
1736 IF(rn(1).GT.wt) goto 300
1738 ELSEIF(mode.EQ. 1)
THEN
1739 * =======================
1740 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1745 902
WRITE(iout, 9020)
1746 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
1749 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1750 c ----------------------------------------------------------------------
1751 * a1 decay unweighted events
1752 c ----------------------------------------------------------------------
1753 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1754 * ,ampiz,ampi,amro,gamro,ama1,gama1
1755 * ,amk,amkz,amkst,gamkst
1757 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1758 * ,ampiz,ampi,amro,gamro,ama1,gama1
1759 * ,amk,amkz,amkst,gamkst
1760 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1761 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1762 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1763 REAL*4 gampmc ,gamper
1764 COMMON / inout / inut,iout
1766 REAL hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4)
1767 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
1770 DATA pi /3.141592653589793238462643/
1774 c ===================
1783 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1784 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1786 cc CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMAA $',100,0,2)
1788 ELSEIF(mode.EQ. 0)
THEN
1789 c =======================
1791 IF(iwarm.EQ.0) goto 902
1792 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1793 cc CALL hfill(801,wt/wtmax)
1798 sswt=sswt+dble(wt)**2
1802 IF(wt.GT.wtmax) nevovr=nevovr+1
1803 IF(rn*wtmax.GT.wt) goto 300
1804 c rotations to basic tau rest frame
1805 costhe=-1.+2.*rrr(2)
1808 CALL rotpol(thet,phi,pnu)
1809 CALL rotpol(thet,phi,paa)
1810 CALL rotpol(thet,phi,pim1)
1811 CALL rotpol(thet,phi,pim2)
1812 CALL rotpol(thet,phi,pipl)
1813 CALL rotpol(thet,phi,hv)
1815 44 hhv(i)=-isgn*hv(i)
1818 ELSEIF(mode.EQ. 1)
THEN
1819 c =======================
1820 IF(nevraw.EQ.0)
RETURN
1821 pargam=swt/float(nevraw+1)
1823 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1825 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1829 cam nevdec(5)=nevacc
1833 7003
FORMAT(///1x,15(5h*****)
1834 $ /,
' *', 25x,
'******** DADMAA INITIALISATION ********',9x,1h*
1835 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1836 $ /,1x,15(5h*****)/)
1837 7010
FORMAT(///1x,15(5h*****)
1838 $ /,
' *', 25x,
'******** DADMAA FINAL REPORT ******** ',9x,1h*
1839 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
1840 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
1841 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1842 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
1843 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1844 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1845 $ /,1x,15(5h*****)/)
1846 902
WRITE(iout, 9020)
1847 9020
FORMAT(
' ----- DADMAA: LACK OF INITIALISATION')
1850 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1851 c ----------------------------------------------------------------------
1852 * it simulates a1 decay in tau rest frame with
1853 * z-axis along a1 momentum
1854 c ----------------------------------------------------------------------
1855 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1856 * ,ampiz,ampi,amro,gamro,ama1,gama1
1857 * ,amk,amkz,amkst,gamkst
1859 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1860 * ,ampiz,ampi,amro,gamro,ama1,gama1
1861 * ,amk,amkz,amkst,gamkst
1862 COMMON / taukle / bra1,brk0,brk0b,brks
1863 REAL*4 bra1,brk0,brk0b,brks
1864 REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4)
1868 c matrix element number:
1870 c
TYPE of the generation:
1874 IF (rmod.LT.bra1)
THEN
1886 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
1888 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
1889 c ----------------------------------------------------------------------
1890 c tau decay into kaon and tau-neutrino
1892 c output four momenta: pnu tauneutrino,
1894 c ----------------------------------------------------------------------
1895 REAL pol(4),hv(4),pnu(4),pkk(4),rn(1)
1898 c ===================
1899 CALL dadmkk(-1,isgn,hv,pkk,pnu)
1900 cc CALL hbook1(815,
'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1902 ELSEIF(mode.EQ. 0)
THEN
1903 c =======================
1905 CALL dadmkk( 0,isgn,hv,pkk,pnu)
1906 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1907 cc CALL hfill(815,wt)
1909 IF(rn(1).GT.wt) goto 300
1911 ELSEIF(mode.EQ. 1)
THEN
1912 c =======================
1913 CALL dadmkk( 1,isgn,hv,pkk,pnu)
1919 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
1920 c ----------------------------------------------------------------------
1922 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1923 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1924 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1925 * ,ampiz,ampi,amro,gamro,ama1,gama1
1926 * ,amk,amkz,amkst,gamkst
1928 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1929 * ,ampiz,ampi,amro,gamro,ama1,gama1
1930 * ,amk,amkz,amkst,gamkst
1931 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1932 REAL*4 gampmc ,gamper
1933 COMMON / inout / inut,iout
1934 REAL pkk(4),pnu(4),hv(4)
1935 DATA pi /3.141592653589793238462643/
1938 c ===================
1940 ELSEIF(mode.EQ. 0)
THEN
1941 c =======================
1943 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
1944 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
1945 xkk= sqrt(ekk**2-amk**2)
1947 CALL sphera(xkk,pkk)
1949 c tau-neutrino momentum
1955 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
1956 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
1957 & +(gv**2-ga**2)*amtau*amnuta*amk**2
1959 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
1962 ELSEIF(mode.EQ. 1)
THEN
1963 c =======================
1964 IF(nevtot.EQ.0)
RETURN
1966 cfz there was brak/amtau**4 before
1967 c gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
1968 c * (brak/amtau**4)**2
1969 czw 7.02.93 here was an error affecting non standard model
1970 c configurations only
1971 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
1973 $ sqrt((amtau**2-amk**2-amnuta**2)**2
1974 $ -4*amk**2*amnuta**2 )/amtau**2
1979 WRITE(iout, 7010) nevtot,gamm,rat,error
1982 cam nevdec(6)=nevtot
1986 7010
FORMAT(///1x,15(5h*****)
1987 $ /,
' *', 25x,
'******** DADMKK FINAL REPORT ********',9x,1h*
1988 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
1989 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
1990 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1991 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1992 $ /,1x,15(5h*****)/)
1994 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
1995 c ----------------------------------------------------------------------
1996 c this simulates tau decay in tau rest frame
1997 c into nu k*,
THEN k* decays into pi0,k+-(jkst=20)
1998 c or pi+-,k0(jkst=10).
1999 c output four momenta: pnu tauneutrino,
2005 c ----------------------------------------------------------------------
2006 COMMON / inout / inut,iout
2007 REAL pol(4),hv(4),pks(4),pnu(4),pkk(4),ppi(4),rn(1)
2011 c ===================
2013 cfz initialisation done with the gharged pion neutral kaon mode(jkst=10
2014 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2015 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXKS $',100,0,2)
2016 cc CALL hbook1(916,
'ABS2 OF HV IN ROUTINE DEXKS $',100,0,2)
2018 ELSEIF(mode.EQ. 0)
THEN
2019 c =======================
2021 IF(iwarm.EQ.0) goto 902
2022 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2023 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2024 cc CALL hfill(816,wt)
2025 cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
2026 cc CALL hfill(916,xhelp)
2028 IF(rn(1).GT.wt) goto 300
2030 ELSEIF(mode.EQ. 1)
THEN
2031 c ======================================
2032 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2038 902
WRITE(iout, 9020)
2039 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2042 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2043 c ----------------------------------------------------------------------
2044 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2045 * ,ampiz,ampi,amro,gamro,ama1,gama1
2046 * ,amk,amkz,amkst,gamkst
2048 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2049 * ,ampiz,ampi,amro,gamro,ama1,gama1
2050 * ,amk,amkz,amkst,gamkst
2051 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2052 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2053 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2054 REAL*4 gampmc ,gamper
2055 COMMON / taukle / bra1,brk0,brk0b,brks
2056 REAL*4 bra1,brk0,brk0b,brks
2057 COMMON / inout / inut,iout
2059 REAL hv(4),pks(4),pnu(4),pkk(4),ppi(4)
2060 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
2061 REAL*4 rrr(3),rmod(1)
2063 DATA pi /3.141592653589793238462643/
2067 c ===================
2076 c the initialisation is done with the 66.7% MODE
2078 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2079 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2081 cc CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMKS $',100,0,2)
2083 cc CALL hbook1(112,
'-------- K* MASS -------- $',100,0.,2.)
2084 ELSEIF(mode.EQ. 0)
THEN
2085 c =====================================
2086 IF(iwarm.EQ.0) goto 902
2087 c here we choose randomly between k0 pi+_(66.7%)
2088 c and k+_ pi0(33.3%)
2092 IF(rmod(1).LT.dec1)
THEN
2097 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2100 IF(wt.GT.wtmax) nevovr=nevovr+1
2104 IF(rn*wtmax.GT.wt) goto 400
2105 c rotations to basic tau rest frame
2106 costhe=-1.+2.*rrr(2)
2109 CALL rotor2(thet,pnu,pnu)
2110 CALL rotor3( phi,pnu,pnu)
2111 CALL rotor2(thet,pks,pks)
2112 CALL rotor3( phi,pks,pks)
2113 CALL rotor2(thet,pkk,pkk)
2114 CALL rotor3(phi,pkk,pkk)
2115 CALL rotor2(thet,ppi,ppi)
2116 CALL rotor3( phi,ppi,ppi)
2117 CALL rotor2(thet,hv,hv)
2118 CALL rotor3( phi,hv,hv)
2120 44 hhv(i)=-isgn*hv(i)
2123 ELSEIF(mode.EQ. 1)
THEN
2124 c =======================
2125 IF(nevraw.EQ.0)
RETURN
2126 pargam=swt/float(nevraw+1)
2128 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2130 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2134 cam nevdec(7)=nevacc
2138 7003
FORMAT(///1x,15(5h*****)
2139 $ /,
' *', 25x,
'******** DADMKS INITIALISATION ********',9x,1h*
2140 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2141 $ /,1x,15(5h*****)/)
2142 7010
FORMAT(///1x,15(5h*****)
2143 $ /,
' *', 25x,
'******** DADMKS FINAL REPORT ********',9x,1h*
2144 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2145 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2146 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2147 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2148 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2149 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2150 $ /,1x,15(5h*****)/)
2151 902
WRITE(iout, 9020)
2152 9020
FORMAT(
' ----- DADMKS: LACK OF INITIALISATION')
2155 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2156 c ----------------------------------------------------------------------
2157 c it simulates kaon* decay in tau rest frame with
2158 c z-axis along kaon* momentum
2159 c jkst=10 for k* --->k0 + pi+-
2160 c jkst=20 for k* --->k+- + pi0
2161 c ----------------------------------------------------------------------
2162 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2163 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2164 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2165 * ,ampiz,ampi,amro,gamro,ama1,gama1
2166 * ,amk,amkz,amkst,gamkst
2168 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2169 * ,ampiz,ampi,amro,gamro,ama1,gama1
2170 * ,amk,amkz,amkst,gamkst
2171 REAL hv(4),pt(4),pn(4),pks(4),pkk(4),ppi(4),qq(4),rr1(1)
2173 DATA pi /3.141592653589793238462643/
2176 c three body phase space normalised as in bjorken-drell
2177 phspac=1./2**11/pi**5
2184 c here begin the k0,pi+_ decay
2186 c ==================
2187 c mass of(real/virtual) k*
2189 ams2=(amtau-amnuta)**2
2191 c amx2=ams1+ rr1(1)*(ams2-ams1)
2193 c phspac=phspac*(ams2-ams1)
2194 c phase space with sampling for k* resonance
2195 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2196 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2197 alp=alp1+rr1(1)*(alp2-alp1)
2198 amx2=amkst**2+amkst*gamkst*tan(alp)
2200 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2202 phspac=phspac*(alp2-alp1)
2204 c tau-neutrino momentum
2207 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2208 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2213 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2215 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2218 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2219 pppi=sqrt((enpi-ampi)*(enpi+ampi))
2220 phspac=phspac*(4*pi)*(2*pppi/amx)
2221 c charged pi momentum in kaon* rest frame
2222 CALL sphera(pppi,ppi)
2224 c neutral kaon momentum in k* rest frame
2227 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2228 exe=(pks(4)+pks(3))/amx
2229 c pion and k boosted from k* rest frame to tau rest frame
2230 CALL bostr3(exe,ppi,ppi)
2231 CALL bostr3(exe,pkk,pkk)
2233 30 qq(i)=ppi(i)-pkk(i)
2234 c qq transverse to pks
2235 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2236 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2238 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
2241 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2243 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2244 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2245 & +(gv**2-ga**2)*amtau*amnuta*qq2
2246 c a simple breit-wigner is chosen for k* resonance
2247 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2248 amplit=(gfermi*scabib)**2*brak*2*fks
2249 dgamt=1/(2.*amtau)*amplit*phspac
2251 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2253 c here begin the k+-,pi0 decay
2254 ELSEIF(jkst.EQ.20)
THEN
2255 c ======================
2256 c mass of(real/virtual) k*
2258 ams2=(amtau-amnuta)**2
2260 c amx2=ams1+ rr1*(ams2-ams1)
2262 c phspac=phspac*(ams2-ams1)
2263 c phase space with sampling for k* resonance
2264 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2265 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2266 alp=alp1+rr1(1)*(alp2-alp1)
2267 amx2=amkst**2+amkst*gamkst*tan(alp)
2269 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2271 phspac=phspac*(alp2-alp1)
2273 c tau-neutrino momentum
2276 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2277 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2281 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2283 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2286 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2287 pppi=sqrt((enpi-ampiz)*(enpi+ampiz))
2288 phspac=phspac*(4*pi)*(2*pppi/amx)
2289 c neutral pi momentum in k* rest frame
2290 CALL sphera(pppi,ppi)
2292 c charged kaon momentum in k* rest frame
2295 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2296 exe=(pks(4)+pks(3))/amx
2297 c pion and k boosted from k* rest frame to tau rest frame
2298 CALL bostr3(exe,ppi,ppi)
2299 CALL bostr3(exe,pkk,pkk)
2301 60 qq(i)=pkk(i)-ppi(i)
2302 c qq transverse to pks
2303 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2304 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2306 61 qq(i)=qq(i)-pks(i)*qqpks/pksd
2309 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2311 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2312 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2313 & +(gv**2-ga**2)*amtau*amnuta*qq2
2314 c a simple breit-wigner is chosen for the k* resonance
2315 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2316 amplit=(gfermi*scabib)**2*brak*2*fks
2317 dgamt=1/(2.*amtau)*amplit*phspac
2319 70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2326 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2327 c ----------------------------------------------------------------------
2328 c it simulates multipi decay in tau rest frame with
2329 c z-axis opposite to neutrino momentum
2330 c ----------------------------------------------------------------------
2331 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2332 * ,ampiz,ampi,amro,gamro,ama1,gama1
2333 * ,amk,amkz,amkst,gamkst
2335 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2336 * ,ampiz,ampi,amro,gamro,ama1,gama1
2337 * ,amk,amkz,amkst,gamkst
2338 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2339 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2340 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2341 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2343 CHARACTER names(nmode)*31
2346 REAL*8 pn(4),pr(4),ppi(4,9),hv(4)
2347 REAL*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2348 REAL*8 pv(5,9),pt(4),ue(3),be(3)
2349 REAL*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2353 REAL*8 gam,bep,phi,a,b,c
2355 REAL*4 rrr(9),rrx(2),rn(1),rr2(1)
2357 DATA pi /3.141592653589793238462643/
2358 DATA wetmax /20*1d-15/
2360 cc-- pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2.*a)
2363 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2365 ampik(i,j)=dcdmas(idffin(i,j))
2368 IF ((jnpi.LE.0).OR.jnpi.GT.20)
THEN
2369 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2383 phspac = 1./2.**5 /pi**2
2385 4 ps =ps+ampik(i,jnpi)
2388 ams2=(amtau-amnuta)**2
2391 amx2=ams1+ rr2(1)*(ams2-ams1)
2394 phspac=phspac * (ams2-ams1)
2396 c tau-neutrino momentum
2399 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2400 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2404 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2406 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2408 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
2409 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
2413 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2414 c here was an error. 20.10.91 (zw)
2415 c brak=2*(gv**2+ga**2)*(2*pxq*pxn+amx2*qxn)
2416 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2417 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2418 cam assume neutrino mass=0. and sum over final polarisation
2419 c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
2420 amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2421 dgamt=1./(2.*amtau)*amplit*phspac
2423 c isotropic w decay in w rest frame
2428 pv(5,nd)=ampik(nd,jnpi)
2429 c compute max. phase space factor
2430 pmax=amw-ps+ampik(nd,jnpi)
2433 pmax=pmax+ampik(il,jnpi)
2434 pmin=pmin+ampik(il+1,jnpi)
2435 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2437 c --- 2.02.94 zw 9 lines
2442 223 ams1=ams1+ampik(jl,jnpi)
2444 amx =(amx-ampik(il,jnpi))
2446 phsmax=phsmax * (ams2-ams1)
2451 cam generate nd-2 effective masses
2453 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2455 CALL ranmar(rrr,nd-2)
2459 231 ams1=ams1+ampik(jl,jnpi)
2461 ams2=(amx-ampik(il,jnpi))**2
2463 amx2=ams1+ rr1*(ams2-ams1)
2466 phspac=phspac * (ams2-ams1)
2467 c --- 2.02.94 zw 1 line
2468 phs=phs* (ams2-ams1)
2469 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2470 phs =phs *pa/pv(5,il)
2472 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2473 phs =phs *pa/pv(5,nd-1)
2475 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2476 IF (ncont.EQ.500 000)
THEN
2479 xnpi=xnpi+ampik(kk,jnpi)
2481 WRITE(6,*)
'ROUNDING INSTABILITY IN DPHNPI ?'
2482 WRITE(6,*)
'AMW=',amw,
'XNPI=',xnpi
2483 WRITE(6,*)
'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2484 WRITE(6,*)
'PHS=',phs,
'PHSMAX=',phsmax
2487 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs) go to 100
2488 c...perform successive two-particle decays in respective cm frame
2489 280
DO 300 il=1,nd-1
2490 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2494 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2495 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2498 290 pv(j,il+1)=-pa*ue(j)
2499 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2500 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2501 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2503 c...lorentz transform decay products to tau frame
2505 310 ppi(j,nd)=pv(j,nd)
2508 320 be(j)=pv(j,il)/pv(4,il)
2509 gam=pv(4,il)/pv(5,il)
2511 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2513 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2514 ppi(4,i)=gam*(ppi(4,i)+bep)
2531 FUNCTION sigee(Q2,JNP)
2532 c ----------------------------------------------------------------------
2533 c e+e- cross section in the(1.gev2,amtau**2) region
2534 c normalised to sig0 = 4/3 pi alfa2
2535 c used in matrix element for multipion tau decays
2536 c cf ys.tsai phys.rev d4 ,2821(1971)
2537 c f.gilman et al phys.rev d17,1846(1978)
2538 c c.kiesling, to be pub. in high energy e+e- physics(1988)
2539 c datsig(*,1) = e+e- -> pi+pi-2pi0
2540 c datsig(*,2) = e+e- -> 2pi+2pi-
2541 c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2542 c(phys lett 78b,623(1978)
2543 c datsig(*,5) = e+e- -> 6pi
2545 c 4- and 6-pion cross sections from
data
2546 c 5-pion contribution related to 4-pion cross section
2549 c ----------------------------------------------------------------------
2550 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2551 * ,ampiz,ampi,amro,gamro,ama1,gama1
2552 * ,amk,amkz,amkst,gamkst
2554 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2555 * ,ampiz,ampi,amro,gamro,ama1,gama1
2556 * ,amk,amkz,amkst,gamkst
2560 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2561 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2562 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2563 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2566 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2569 DATA pi /3.141592653589793238462643/
2577 c ajwmod: initialize
if called from outside qq:
2578 IF (ampi.LT.0.139) ampi = 0.1395675
2582 datsig(i,2) = datsig(i,2)/2.
2583 datsig(i,1) = datsig(i,1) + datsig(i,2)
2589 IF(t . gt. s-ampi ) go to 201
2591 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2592 fact = fact * (datsig(j,1)+datsig(j+1,1))
2593 200 datsig(i,3) = datsig(i,3) + fact
2594 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2595 datsig(i,4) = datsig(i,3)
2596 datsig(i,6) = datsig(i,5)
2598 c
WRITE(6,1000) datsig
2599 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2605 sigee=datsig(1,jnpi)+
2606 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2607 ELSEIF(q.LT.1.8)
THEN
2610 IF(q.LT.qmax) go to 2
2613 2 sigee=datsig(i,jnpi)+
2614 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2615 ELSEIF(q.GT.1.8)
THEN
2616 sigee=datsig(17,jnpi)+
2617 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2619 IF(sigee.LT..0) sigee=0.
2621 sigee = sigee/(6.*pi**2*sig0)
2626 FUNCTION sigold(Q2,JNPI)
2627 c ----------------------------------------------------------------------
2628 c e+e- cross section in the(1.gev2,amtau**2) region
2629 c normalised to sig0 = 4/3 pi alfa2
2630 c used in matrix element for multipion tau decays
2631 c cf ys.tsai phys.rev d4 ,2821(1971)
2632 c f.gilman et al phys.rev d17,1846(1978)
2633 c c.kiesling, to be pub. in high energy e+e- physics(1988)
2634 c datsig(*,1) = e+e- -> pi+pi-2pi0
2635 c datsig(*,2) = e+e- -> 2pi+2pi-
2636 c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2637 c(phys lett 78b,623(1978)
2638 c datsig(*,4) = e+e- -> 6pi
2640 c 4- and 6-pion cross sections from
data
2641 c 5-pion contribution related to 4-pion cross section
2644 c ----------------------------------------------------------------------
2645 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2646 * ,ampiz,ampi,amro,gamro,ama1,gama1
2647 * ,amk,amkz,amkst,gamkst
2649 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2650 * ,ampiz,ampi,amro,gamro,ama1,gama1
2651 * ,amk,amkz,amkst,gamkst
2655 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2656 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2657 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2658 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2660 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2662 DATA pi /3.141592653589793238462643/
2670 datsig(i,2) = datsig(i,2)/2.
2671 datsig(i,1) = datsig(i,1) + datsig(i,2)
2677 IF(t . gt. s-ampi ) go to 201
2679 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2680 fact = fact * (datsig(j,1)+datsig(j+1,1))
2681 200 datsig(i,3) = datsig(i,3) + fact
2682 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2684 c
WRITE(6,1000) datsig
2685 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2691 sigee=datsig(1,jnpi)+
2692 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2693 ELSEIF(q.LT.1.8)
THEN
2696 IF(q.LT.qmax) go to 2
2699 2 sigee=datsig(i,jnpi)+
2700 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2701 ELSEIF(q.GT.1.8)
THEN
2702 sigee=datsig(17,jnpi)+
2703 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2705 IF(sigee.LT..0) sigee=0.
2707 sigee = sigee/(6.*pi**2*sig0)
2712 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2713 c ----------------------------------------------------------------------
2714 * it simulates three pi(k) decay in the tau rest frame
2715 * z-axis along hadronic system
2716 c ----------------------------------------------------------------------
2717 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2718 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2720 CHARACTER names(nmode)*31
2722 REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pnpi(4,9)
2723 c matrix element number:
2725 c
TYPE of the generation:
2728 c --- masses of the decay products
2729 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
2730 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
2731 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
2733 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2745 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
2746 c ----------------------------------------------------------------------
2747 * it simulates a1 decay in tau rest frame with
2748 * z-axis along a1 momentum
2749 * it can be also used to generate k k pi and k pi pi tau decays.
2751 * keyt - algorithm controlling switch
2752 * 2 - flat phase space pim1 pim2 symmetrized statistical factor 1/2
2753 * 1 - like 1 but peaked around a1 and rho(two channels) masses.
2754 * 3 - peaked around omega, all particles different
2755 * other- flat phase space, all particles different
2756 * amp1 - mass of first pi, etc. (1-3)
2757 * mnum - matrix element
type
2758 * 0 - a1 matrix element
2759 * 1-6 - matrix element for k pi pi, k k pi decay modes
2760 * 7 - pi- pi0 gamma matrix element
2761 c ----------------------------------------------------------------------
2762 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2763 * ,ampiz,ampi,amro,gamro,ama1,gama1
2764 * ,amk,amkz,amkst,gamkst
2766 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2767 * ,ampiz,ampi,amro,gamro,ama1,gama1
2768 * ,amk,amkz,amkst,gamkst
2769 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2770 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2771 REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4)
2774 DATA pi /3.141592653589793238462643/
2776 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
2777 c amro, gamro is only a
PARAMETER for geting hight efficiency
2779 c three body phase space normalised as in bjorken-drell
2780 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
2781 phspac=1./2**17/pi**8
2791 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
2792 $ amrx,gamrx,amra,gamra,amrb,gamrb)
2793 IF (ichan.EQ.1)
THEN
2796 ELSEIF (ichan.EQ.2)
THEN
2805 ams1=(amp1+amp2+amp3)**2
2806 ams2=(amtau-amnuta)**2
2807 * phase space with sampling for a1 resonance
2808 alp1=atan((ams1-amrx**2)/amrx/gamrx)
2809 alp2=atan((ams2-amrx**2)/amrx/gamrx)
2810 alp=alp1+rr1*(alp2-alp1)
2811 am3sq =amrx**2+amrx*gamrx*tan(alp)
2813 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
2814 phspac=phspac*(alp2-alp1)
2815 c mass of(real/virtual) rho -
2819 IF (ichan.LE.2)
THEN
2820 * phase space with sampling for rho resonance,
2821 alp1=atan((ams1-amra**2)/amra/gamra)
2822 alp2=atan((ams2-amra**2)/amra/gamra)
2823 alp=alp1+rr2*(alp2-alp1)
2824 am2sq =amra**2+amra*gamra*tan(alp)
2826 c --- this part of the jacobian will be recovered later ---------------
2827 c phspac=phspac*(alp2-alp1)
2828 c phspac=phspac*((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
2829 c----------------------------------------------------------------------
2832 am2sq=ams1+ rr2*(ams2-ams1)
2836 * rho restframe, define pipl and pim1
2837 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
2838 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
2839 ppi= enq1**2-amp3**2
2840 pppi=sqrt(abs(enq1**2-amp3**2))
2841 c --- this part of jacobian will be recovered later
2842 phf1=(4*pi)*(2*pppi/am2)
2843 * pi minus momentum in rho rest frame
2844 CALL sphera(pppi,pipl)
2846 * pi0 1 momentum in rho rest frame
2850 * a1 rest frame, define pim2
2854 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
2855 pr(3)= sqrt(abs(pr(4)**2-am2**2))
2856 ppi = pr(4)**2-am2**2
2860 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
2862 phf2=(4*pi)*(2*pr(3)/am3)
2863 * old pions boosted from rho rest frame to a1 rest frame
2864 exe=(pr(4)+pr(3))/am2
2865 CALL bostr3(exe,pipl,pipl)
2866 CALL bostr3(exe,pim1,pim1)
2870 thet =acos(-1.+2*rr3)
2872 CALL rotpol(thet,phi,pipl)
2873 CALL rotpol(thet,phi,pim1)
2874 CALL rotpol(thet,phi,pim2)
2875 CALL rotpol(thet,phi,pr)
2877 * now to the tau rest frame, define a1 and neutrino momenta
2881 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
2882 paa(3)= sqrt(abs(paa(4)**2-am3**2))
2883 ppi = paa(4)**2-am3**2
2884 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
2885 * tau-neutrino momentum
2888 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
2890 c here we correct for the jacobians of the two chains
2891 c ---first channel ------- pim1+pipl
2894 alp1=atan((ams1-amra**2)/amra/gamra)
2895 alp2=atan((ams2-amra**2)/amra/gamra)
2896 xpro = (pim1(3)+pipl(3))**2
2897 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
2898 am2sq=-xpro+(pim1(4)+pipl(4))**2
2899 c jacobian of speeding
2900 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
2901 ff1 =ff1 *(alp2-alp1)
2902 c lambda of rho decay
2903 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
2904 c lambda of a1 decay
2905 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
2906 xjaje=gg1*(ams2-ams1)
2907 c ---second channel ------ pim2+pipl
2910 alp1=atan((ams1-amrb**2)/amrb/gamrb)
2911 alp2=atan((ams2-amrb**2)/amrb/gamrb)
2912 xpro = (pim2(3)+pipl(3))**2
2913 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
2914 am2sq=-xpro+(pim2(4)+pipl(4))**2
2915 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
2916 ff2 =ff2 *(alp2-alp1)
2917 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
2918 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
2919 xjadw=gg2*(ams2-ams1)
2926 IF (ichan.EQ.2)
THEN
2931 IF (xjac1.NE.0.0) a1=prob1/xjac1
2932 IF (xjac2.NE.0.0) a2=prob2/xjac2
2933 IF (xjac3.NE.0.0) a3=prob3/xjac3
2935 IF (a1+a2+a3.NE.0.0)
THEN
2936 phspac=phspac/(a1+a2+a3)
2946 * all pions boosted from a1 rest frame to tau rest frame
2947 * z-axis antiparallel to neutrino momentum
2948 exe=(paa(4)+paa(3))/am3
2949 CALL bostr3(exe,pipl,pipl)
2950 CALL bostr3(exe,pim1,pim1)
2951 CALL bostr3(exe,pim2,pim2)
2952 CALL bostr3(exe,pr,pr)
2953 c partial width consists of phase space and amplitude
2955 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
2956 c
ELSEIF (mnum.EQ.0)
THEN
2957 c CALL dampaa(pt,pn,pim1,pim2,pipl,amplit,hv)
2959 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
2961 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN
2962 c the statistical factor for identical pi-s is cancelled with
2963 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
2967 dgamt=1/(2.*amtau)*amplit*phspac
2969 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
2970 c ----------------------------------------------------------------------
2971 * calculates differential cross section and polarimeter vector
2972 * for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
2973 * all spin effects in the full decay chain are taken into account.
2974 * calculations done in tau rest frame with z-axis along neutrino moment
2975 * the routine is writen for zero neutrino mass.
2977 c called by : dphsaa
2978 c ----------------------------------------------------------------------
2979 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2980 * ,ampiz,ampi,amro,gamro,ama1,gama1
2981 * ,amk,amkz,amkst,gamkst
2983 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2984 * ,ampiz,ampi,amro,gamro,ama1,gama1
2985 * ,amk,amkz,amkst,gamkst
2986 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2987 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2988 COMMON /testa1/ keya1
2989 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
2990 REAL paa(4),vec1(4),vec2(4)
2991 REAL pivec(4),piaks(4),hvm(4)
2992 COMPLEX bwign,hadcur(4),fpik
2995 * f constants for a1, a1-rho-pi, and rho-pi-pi
2998 * this inline funct. calculates the scalar part of the propagator
2999 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3001 * four momentum of a1
3003 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3004 * masses of a1, and of two pi-pairs which may form rho
3005 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3006 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3007 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3008 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3009 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3010 * elements of hadron current
3011 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3012 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3013 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3014 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3016 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3017 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3018 * hadron current saturated with a1 and rho resonances
3019 IF (keya1.EQ.1)
THEN
3023 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3025 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3026 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3027 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3030 fnorm=2.0*sqrt(2.)/3.0/fpi
3031 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3033 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3034 $ *(cmplx(vec1(i))*fpik(xmro1)
3035 $ +cmplx(vec2(i))*fpik(xmro2))
3039 * calculate pi-vectors: vector and axial
3040 CALL clvec(hadcur,pn,pivec)
3041 CALL claxi(hadcur,pn,piaks)
3042 CALL clnut(hadcur,brakm,hvm)
3043 * spin independent part of decay diff-cross-sect. in tau rest frame
3044 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3045 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3046 amplit=(gfermi*ccabib)**2*brak/2.
3047 c the statistical factor for identical pi-s was cancelled with
3048 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3049 c polarimeter vector in tau rest frame
3051 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3052 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3053 c hv is defined for tau- with gamma=b+hv*pol
3059 c ****************************************************************
3060 c g-
FUNCTION used to inroduce energy dependence in a1 width
3061 c ****************************************************************
3062 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3063 * ,ampiz,ampi,amro,gamro,ama1,gama1
3064 * ,amk,amkz,amkst,gamkst
3066 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3067 * ,ampiz,ampi,amro,gamro,ama1,gama1
3068 * ,amk,amkz,amkst,gamkst
3070 IF (qkwa.LT.(amro+ampi)**2)
THEN
3071 gfun=4.1*(qkwa-9*ampiz**2)**3
3072 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3074 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3077 COMPLEX FUNCTION bwigs(S,M,G)
3078 c **********************************************************
3079 c p-wave breit-wigner for k*
3080 c **********************************************************
3082 REAL pi,pim,qs,qm,w,gs,mk
3083 c ajw: add k*-prim possibility:
3087 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3088 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3089 c ------------ parameters --------------------
3095 c ajw: add k*-prim possibility:
3099 c ------- breit-wigner -----------------------
3101 qs=p(s,pim**2,mk**2)
3102 qm=p(m**2,pim**2,mk**2)
3104 gs=g*(m/w)*(qs/qm)**3
3105 bw=m**2/cmplx(m**2-s,-m*gs)
3106 qpm=p(pm**2,pim**2,mk**2)
3107 g1=pg*(pm/w)*(qs/qpm)**3
3108 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3109 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3112 COMPLEX FUNCTION bwig(S,M,G)
3113 c **********************************************************
3114 c p-wave breit-wigner for rho
3115 c **********************************************************
3117 REAL pi,pim,qs,qm,w,gs
3119 c ------------ parameters --------------------
3124 c ------- breit-wigner -----------------------
3126 IF (s.GT.4.*pim**2)
THEN
3127 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3128 qm=sqrt(m**2/4.-pim**2)
3130 gs=g*(m/w)*(qs/qm)**3
3134 bwig=m**2/cmplx(m**2-s,-m*gs)
3137 COMPLEX FUNCTION fpik(W)
3138 c **********************************************************
3140 c **********************************************************
3142 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
3146 c ------------ parameters --------------------
3147 IF (init.EQ.0 )
THEN
3157 c -----------------------------------------------
3159 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3164 c **********************************************************
3165 c square of pion form factor
3166 c **********************************************************
3168 fpirho=cabs(fpik(w))**2
3170 SUBROUTINE clvec(HJ,PN,PIV)
3171 c ----------------------------------------------------------------------
3172 * calculates the
"VECTOR TYPE" pi-vector piv
3173 * note that the neutrino mom. pn is assumed to be along z-axis
3175 c called by : dampaa
3176 c ----------------------------------------------------------------------
3180 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3181 hh=
REAL(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3182 $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3184 10 piv(i)=4.*
REAL(hn*conjg(hj(i)))-2.*hh*pn(i)
3187 SUBROUTINE claxi(HJ,PN,PIA)
3188 c ----------------------------------------------------------------------
3189 * calculates the
"AXIAL TYPE" pi-vector pia
3190 * note that the neutrino mom. pn is assumed to be along z-axis
3191 c sign is chosen +/- for decay of tau +/- respectively
3192 c called by : dampaa, clnut
3193 c ----------------------------------------------------------------------
3194 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3195 COMMON / idfc / idff
3197 COMPLEX hj(4),hjc(4)
3198 c det2(i,j)=aimag(hj(i)*hjc(j)-hj(j)*hjc(i))
3199 c -- here was an error(zw, 21.11.1991)
3200 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3201 c -- it was affecting sign of a_lr asymmetry in a1 decay.
3202 c -- note also collision of notation of gamma_va as defined in
3203 c -- tauola paper and j.h. kuhn and santamaria z. phys c 48 (1990) 445
3204 * -----------------------------------
3205 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3206 sign= idff/abs(idff)
3207 ELSEIF (ktom.EQ.2)
THEN
3208 sign=-idff/abs(idff)
3210 print *,
'STOP IN CLAXI: KTOM=',ktom
3215 10 hjc(i)=conjg(hj(i))
3216 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3217 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3218 pia(3)= 2.*pn(4)*det2(1,2)
3219 pia(4)= 2.*pn(3)*det2(1,2)
3220 c all four indices are up so pia(3) and pia(4) have same sign
3222 20 pia(i)=pia(i)*sign
3224 SUBROUTINE clnut(HJ,B,HV)
3225 c ----------------------------------------------------------------------
3226 * calculates the contribution by neutrino mass
3227 * note the tau is assumed to be at rest
3229 c called by : dampaa
3230 c ----------------------------------------------------------------------
3236 b=
REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3))
& - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) )
3239 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3240 c ----------------------------------------------------------------------
3241 * calculates differential cross section and polarimeter vector
3242 * for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
3243 * all spin effects in the full decay chain are taken into account.
3244 * calculations done in tau rest frame with z-axis along neutrino moment
3245 * the routine is writen for zero neutrino mass.
3247 c called by : dphsaa
3248 c ----------------------------------------------------------------------
3249 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3250 * ,ampiz,ampi,amro,gamro,ama1,gama1
3251 * ,amk,amkz,amkst,gamkst
3253 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3254 * ,ampiz,ampi,amro,gamro,ama1,gama1
3255 * ,amk,amkz,amkst,gamkst
3256 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3257 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3258 COMMON /testa1/ keya1
3259 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
3260 REAL paa(4),vec1(4),vec2(4)
3261 REAL pivec(4),piaks(4),hvm(4)
3262 COMPLEX bwign,hadcur(4),fnorm,formom
3264 * this inline funct. calculates the scalar part of the propagator
3265 c ajwmod to satisfy compiler, comment out this unused function.
3267 * four momentum of a1
3272 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3274 * masses of a1, and of two pi-pairs which may form rho
3275 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3276 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3277 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3278 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3279 * elements of hadron current
3280 prod1 =vec1(1)*pipl(1)
3281 prod2 =vec2(2)*pipl(2)
3282 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3283 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3284 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3285 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3286 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3287 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3289 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3291 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3293 vec1(i)= vec1(i)/gnorm
3295 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3296 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3297 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3298 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3299 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3300 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3301 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3302 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3303 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3304 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3305 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3307 fnorm=formom(xmaa,xmom)
3312 hadcur(i) = fnorm *(
3313 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3314 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3315 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3317 hadcur(i) = fnorm *(
3318 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3319 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3320 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3324 * calculate pi-vectors: vector and axial
3325 CALL clvec(hadcur,pn,pivec)
3326 CALL claxi(hadcur,pn,piaks)
3327 CALL clnut(hadcur,brakm,hvm)
3328 * spin independent part of decay diff-cross-sect. in tau rest frame
3329 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3330 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3332 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3333 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3335 c hv is defined for tau- with gamma=b+hv*pol
3337 amplit=(gfermi*ccabib)**2*brak/2.
3338 c the statistical factor for identical pi-s was cancelled with
3339 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3340 c polarimeter vector in tau rest frame
3346 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3347 c ----------------------------------------------------------------------
3348 * calculates differential cross section and polarimeter vector
3349 * for tau decay into k k pi, k pi pi.
3350 * all spin effects in the full decay chain are taken into account.
3351 * calculations done in tau rest frame with z-axis along neutrino moment
3352 c mnum decay mode identifier.
3354 c called by : dphsaa
3355 c ----------------------------------------------------------------------
3356 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3357 * ,ampiz,ampi,amro,gamro,ama1,gama1
3358 * ,amk,amkz,amkst,gamkst
3360 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3361 * ,ampiz,ampi,amro,gamro,ama1,gama1
3362 * ,amk,amkz,amkst,gamkst
3363 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3364 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3365 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4)
3366 REAL paa(4),vec1(4),vec2(4),vec3(4),vec4(4),vec5(4)
3367 REAL pivec(4),piaks(4),hvm(4)
3368 REAL fnorm(0:7),coef(1:5,0:7)
3369 COMPLEX hadcur(4),form1,form2,form3,form4,form5,uroj
3370 COMPLEX f1,f2,f3,f4,f5
3371 EXTERNAL form1,form2,form3,form4,form5
3372 DATA pi /3.141592653589793238462643/
3376 IF (icont.EQ.0)
THEN
3384 fnorm(4)=scabib/fpi/dwapi0
3389 coef(1,0)= 2.0*sqrt(2.)/3.0
3390 coef(2,0)=-2.0*sqrt(2.)/3.0
3391 c ajw 2/98: add in the d-wave and i=0 3pi substructure:
3392 coef(3,0)= 2.0*sqrt(2.)/3.0
3396 coef(1,1)=-sqrt(2.)/3.0
3397 coef(2,1)= sqrt(2.)/3.0
3402 coef(1,2)=-sqrt(2.)/3.0
3403 coef(2,2)= sqrt(2.)/3.0
3408 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3415 coef(1,4)= 1.0/sqrt(2.)/3.0
3416 coef(2,4)=-1.0/sqrt(2.)/3.0
3421 coef(1,5)=-sqrt(2.)/3.0
3422 coef(2,5)= sqrt(2.)/3.0
3427 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3438 coef(5,7)=-sqrt(2.0/3.0)
3443 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3444 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3445 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3446 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3447 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3448 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3449 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3450 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3451 * elements of hadron current
3452 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3453 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3454 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3455 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3456 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3457 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3459 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3460 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3461 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3462 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3463 CALL prod5(pim1,pim2,pim3,vec5)
3465 c be aware that sign of vec2 is opposite to sign of vec1 in a1
case
3466 c rationalize this code:
3467 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3468 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3469 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3471 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3472 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3473 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3476 hadcur(i)= cmplx(fnorm(mnum)) * (
3477 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3478 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3481 * calculate pi-vectors: vector and axial
3482 CALL clvec(hadcur,pn,pivec)
3483 CALL claxi(hadcur,pn,piaks)
3484 CALL clnut(hadcur,brakm,hvm)
3485 * spin independent part of decay diff-cross-sect. in tau rest frame
3486 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3487 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3488 amplit=(gfermi)**2*brak/2.
3490 print *,
'MNUM=',mnum
3496 IF (k.EQ.4) znak=1.0
3497 xm1=znak*pim1(k)**2+xm1
3498 xm2=znak*pim2(k)**2+xm2
3499 xm3=znak*pim3(k)**2+xm3
3500 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3501 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3502 print *,
'************************************************'
3504 c polarimeter vector in tau rest frame
3506 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3507 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3508 c hv is defined for tau- with gamma=b+hv*pol
3512 SUBROUTINE prod5(P1,P2,P3,PIA)
3513 c ----------------------------------------------------------------------
3514 c
external product of p1, p2, p3 4-momenta.
3515 c sign is chosen +/- for decay of tau +/- respectively
3516 c called by : dampaa, clnut
3517 c ----------------------------------------------------------------------
3518 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3519 COMMON / idfc / idff
3520 REAL pia(4),p1(4),p2(4),p3(4)
3521 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3522 * -----------------------------------
3523 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3524 sign= idff/abs(idff)
3525 ELSEIF (ktom.EQ.2)
THEN
3526 sign=-idff/abs(idff)
3528 print *,
'STOP IN PROD5: KTOM=',ktom
3532 c epsilon( p1(1), p2(2), p3(3), (4) ) = 1
3534 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3535 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3536 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3537 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3538 c all four indices are up so pia(3) and pia(4) have same sign
3540 20 pia(i)=pia(i)*sign
3543 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3544 c ----------------------------------------------------------------------
3545 * this simulates tau decay in tau rest frame
3546 * into nu a1, next a1 decays into rho pi and finally rho into pi pi.
3547 * output four momenta: pnu tauneutrino,
3549 * pim1 pion minus(or pi0) 1 (for tau minus)
3550 * pim2 pion minus(or pi0) 2
3551 * pipl pion plus(or pi-)
3552 * (pipl,pim1) form a rho
3553 c ----------------------------------------------------------------------
3554 COMMON / inout / inut,iout
3555 REAL pol(4),hv(4),paa(4),pnu(4),pnpi(4,9),rn(1)
3559 c ===================
3561 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3562 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3564 ELSEIF(mode.EQ. 0)
THEN
3565 * =======================
3567 IF(iwarm.EQ.0) goto 902
3568 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3569 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3570 cc CALL hfill(816,wt)
3572 IF(rn(1).GT.wt) goto 300
3574 ELSEIF(mode.EQ. 1)
THEN
3575 * =======================
3576 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3581 902
WRITE(iout, 9020)
3582 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3585 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3586 c ----------------------------------------------------------------------
3587 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3588 * ,ampiz,ampi,amro,gamro,ama1,gama1
3589 * ,amk,amkz,amkst,gamkst
3591 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3592 * ,ampiz,ampi,amro,gamro,ama1,gama1
3593 * ,amk,amkz,amkst,gamkst
3594 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3595 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3596 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3597 REAL*4 gampmc ,gamper
3598 COMMON / inout / inut,iout
3599 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3600 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3602 CHARACTER names(nmode)*31
3604 REAL*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3605 REAL*4 pdum1(4),pdum2(4),pdumi(4,9)
3608 REAL*8 swt(nmode),sswt(nmode)
3609 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3611 DATA pi /3.141592653589793238462643/
3615 c ===================
3616 c -- at the moment only two decay modes of multipions have m. elem
3627 c for 4pi phase space, need lots more trials at initialization,
3628 c or
use the wtmax determined with many trials for default model:
3630 IF (jnpi.LE.nm4)
THEN
3631 a = pkorb(3,37+jnpi)
3642 ELSEIF(jnpi.LE.nm4)
THEN
3643 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3644 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3645 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3646 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3647 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3648 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3649 inum=jnpi-nm4-nm5-nm6
3650 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
3651 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3652 inum=jnpi-nm4-nm5-nm6-nm3
3653 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
3657 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
3659 c print *,
' DADNEW JNPI,NTRIALS,WTMAX =',jnpi,ntrials,wtmax(jnpi)
3660 c CALL hbook1(801,
'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
3661 c print 7004,wtmax(jnpi)
3665 ELSEIF(mode.EQ. 0)
THEN
3666 c =======================
3667 IF(iwarm.EQ.0) goto 902
3672 ELSEIF(jnpi.LE.nm4)
THEN
3673 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3674 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3675 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3676 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3677 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
3678 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3679 inum=jnpi-nm4-nm5-nm6
3680 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
3681 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3682 inum=jnpi-nm4-nm5-nm6-nm3
3683 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
3690 c CALL hfill(801,wt/wtmax(jnpi))
3691 nevraw(jnpi)=nevraw(jnpi)+1
3692 swt(jnpi)=swt(jnpi)+wt
3694 cc sswt(jnpi)=sswt(jnpi)+wt**2
3695 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
3699 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
3700 IF(rn*wtmax(jnpi).GT.wt) goto 300
3701 c rotations to basic tau rest frame
3702 costhe=-1.+2.*rrr(2)
3705 CALL rotor2(thet,pnu,pnu)
3706 CALL rotor3( phi,pnu,pnu)
3707 CALL rotor2(thet,pwb,pwb)
3708 CALL rotor3( phi,pwb,pwb)
3709 CALL rotor2(thet,hv,hv)
3710 CALL rotor3( phi,hv,hv)
3713 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
3714 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
3716 nevacc(jnpi)=nevacc(jnpi)+1
3718 ELSEIF(mode.EQ. 1)
THEN
3719 c =======================
3721 IF(nevraw(jnpi).EQ.0) goto 500
3722 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
3724 IF(nevraw(jnpi).NE.0)
3725 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
3727 WRITE(iout, 7010) names(jnpi),
3728 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
3730 gampmc(8+jnpi-1)=rat
3731 gamper(8+jnpi-1)=error
3732 cam nevdec(8+jnpi-1)=nevacc(jnpi)
3737 7003
FORMAT(///1x,15(5h*****)
3738 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
3740 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
3742 $ /,1x,15(5h*****)/)
3743 7010
FORMAT(///1x,15(5h*****)
3744 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
3745 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
3746 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
3747 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
3748 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
3749 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
3750 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
3751 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
3752 $ /,1x,15(5h*****)/)
3753 902
WRITE(iout, 9020)
3754 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
3756 903
WRITE(iout, 9030) jnpi,mode
3757 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
3762 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
3763 c ----------------------------------------------------------------------
3764 * it simulates a1 decay in tau rest frame with
3765 * z-axis along a1 momentum
3766 c ----------------------------------------------------------------------
3767 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3768 * ,ampiz,ampi,amro,gamro,ama1,gama1
3769 * ,amk,amkz,amkst,gamkst
3771 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3772 * ,ampiz,ampi,amro,gamro,ama1,gama1
3773 * ,amk,amkz,amkst,gamkst
3774 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3775 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3776 REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pmult(4,9)
3779 REAL*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
3780 DATA pi /3.141592653589793238462643/
3782 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3783 c amro, gamro is only a
PARAMETER for geting hight efficiency
3785 c three body phase space normalised as in bjorken-drell
3786 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
3787 phspac=1./2**23/pi**11
3797 c ajw: cant simply change amrop, etc, here
3798 c choice is a by-hand tuning/optimization, no simple relationship
3799 c to actual resonance masses(accd to z.was).
3800 c what matters in the
end is what you put in formf/curr .
3816 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
3817 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
3827 * masses of 4, 3 and 2 pi systems
3828 c 3 pi with sampling for resonance
3831 ams1=(amp1+amp2+amp3+amp4)**2
3832 ams2=(amtau-amnuta)**2
3833 alp1=atan((ams1-amrop**2)/amrop/gamrop)
3834 alp2=atan((ams2-amrop**2)/amrop/gamrop)
3835 alp=alp1+rr1*(alp2-alp1)
3836 am4sq =amrop**2+amrop*gamrop*tan(alp)
3839 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
3840 phspac=phspac*(alp2-alp1)
3844 ams1=(amp2+amp3+amp4)**2
3846 IF (rrr(9).GT.prez)
THEN
3847 am3sq=ams1+ rr1*(ams2-ams1)
3849 c --- this part of jacobian will be recovered later
3852 * phase space with sampling for omega resonance,
3853 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3854 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3855 alp=alp1+rr1*(alp2-alp1)
3856 am3sq =amrx**2+amrx*gamrx*tan(alp)
3858 c --- this part of the jacobian will be recovered later ---------------
3859 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3867 am2sq=ams1+ rr2*(ams2-ams1)
3869 c --- this part of jacobian will be recovered later
3871 * 2 restframe, define piz and pipl
3872 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
3873 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
3874 ppi= enq1**2-amp4**2
3875 pppi=sqrt(abs(enq1**2-amp4**2))
3876 phspac=phspac*(4*pi)*(2*pppi/am2)
3877 * piz momentum in 2 rest frame
3878 CALL sphera(pppi,piz)
3880 * pipl momentum in 2 rest frame
3884 * 3 rest frame, define pim1
3888 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
3889 pr(3)= sqrt(abs(pr(4)**2-am2**2))
3890 ppi = pr(4)**2-am2**2
3894 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
3896 c --- this part of jacobian will be recovered later
3897 ff3=(4*pi)*(2*pr(3)/am3)
3898 * old pions boosted from 2 rest frame to 3 rest frame
3899 exe=(pr(4)+pr(3))/am2
3900 CALL bostr3(exe,piz,piz)
3901 CALL bostr3(exe,pipl,pipl)
3904 thet =acos(-1.+2*rr3)
3906 CALL rotpol(thet,phi,pipl)
3907 CALL rotpol(thet,phi,pim1)
3908 CALL rotpol(thet,phi,piz)
3909 CALL rotpol(thet,phi,pr)
3910 * 4 rest frame, define pim2
3914 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
3915 pr(3)= sqrt(abs(pr(4)**2-am3**2))
3916 ppi = pr(4)**2-am3**2
3920 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
3922 c --- this part of jacobian will be recovered later
3923 ff4=(4*pi)*(2*pr(3)/am4)
3924 * old pions boosted from 3 rest frame to 4 rest frame
3925 exe=(pr(4)+pr(3))/am3
3926 CALL bostr3(exe,piz,piz)
3927 CALL bostr3(exe,pipl,pipl)
3928 CALL bostr3(exe,pim1,pim1)
3931 thet =acos(-1.+2*rr3)
3933 CALL rotpol(thet,phi,pipl)
3934 CALL rotpol(thet,phi,pim1)
3935 CALL rotpol(thet,phi,pim2)
3936 CALL rotpol(thet,phi,piz)
3937 CALL rotpol(thet,phi,pr)
3939 * now to the tau rest frame, define paa and neutrino momenta
3943 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
3944 paa(3)= sqrt(abs(paa(4)**2-am4**2))
3945 ppi = paa(4)**2-am4**2
3946 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3947 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
3948 * tau-neutrino momentum
3951 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
3953 c zbw 20.12.2002 bug fix
3954 IF(rrr(9).LE.0.5*prez)
THEN
3961 c we include remaining part of the jacobian
3963 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
3964 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
3966 ams1=(amp1+amp3+amp4)**2
3969 ams2=(sqrt(am3sq)-amp1)**2
3971 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
3972 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
3975 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
3976 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
3978 ams1=(amp1+amp3+amp4)**2
3979 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3980 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3981 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3984 ams2=(sqrt(am3sq)-amp1)**2
3986 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
3987 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
3989 c --- second channel
3990 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
3991 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
3993 ams1=(amp2+amp3+amp4)**2
3994 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3995 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3996 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3999 ams2=(sqrt(am3sq)-amp2)**2
4001 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4002 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4004 c --- jacobian averaged over the two
4005 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4006 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4011 * momenta of the two pi-minus are randomly symmetrised
4022 c momenta of pi0-s are generated uniformly only
IF prez=0.0
4032 * all pions boosted from 4 rest frame to tau rest frame
4033 * z-axis antiparallel to neutrino momentum
4034 exe=(paa(4)+paa(3))/am4
4035 CALL bostr3(exe,piz,piz)
4036 CALL bostr3(exe,pipl,pipl)
4037 CALL bostr3(exe,pim1,pim1)
4038 CALL bostr3(exe,pim2,pim2)
4039 CALL bostr3(exe,pr,pr)
4040 c partial width consists of phase space and amplitude
4041 c check on consistency with dadnpi,
THEN, code breakes uniform pion
4042 c distribution in hadronic system
4043 cam assume neutrino mass=0. and sum over final polarisation
4045 c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
4046 c amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,1)
4048 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4049 ELSEIF (jnpi.EQ.2)
THEN
4050 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4052 dgamt=1/(2.*amtau)*amplit*phspac
4062 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4063 c ----------------------------------------------------------------------
4064 * calculates differential cross section and polarimeter vector
4065 * for tau decay into 4 pi modes
4066 * all spin effects in the full decay chain are taken into account.
4067 * calculations done in tau rest frame with z-axis along neutrino moment
4068 c mnum decay mode identifier.
4070 c called by : dphsaa
4071 c ----------------------------------------------------------------------
4072 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4073 * ,ampiz,ampi,amro,gamro,ama1,gama1
4074 * ,amk,amkz,amkst,gamkst
4076 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4077 * ,ampiz,ampi,amro,gamro,ama1,gama1
4078 * ,amk,amkz,amkst,gamkst
4079 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4080 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4081 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4),pim4(4)
4082 REAL pivec(4),piaks(4),hvm(4)
4083 COMPLEX hadcur(4),form1,form2,form3,form4,form5
4084 EXTERNAL form1,form2,form3,form4,form5
4085 DATA pi /3.141592653589793238462643/
4088 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4090 * calculate pi-vectors: vector and axial
4091 CALL clvec(hadcur,pn,pivec)
4092 CALL claxi(hadcur,pn,piaks)
4093 CALL clnut(hadcur,brakm,hvm)
4094 * spin independent part of decay diff-cross-sect. in tau rest frame
4095 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4096 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4097 amplit=(ccabib*gfermi)**2*brak/2.
4098 c polarimeter vector in tau rest frame
4100 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4101 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4102 c hv is defined for tau- with gamma=b+hv*pol
4107 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4108 c ----------------------------------------------------------------------
4109 * it simulates 5pi decay in tau rest frame with
4110 * z-axis along 5pi momentum
4111 c ----------------------------------------------------------------------
4112 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4113 * ,ampiz,ampi,amro,gamro,ama1,gama1
4114 * ,amk,amkz,amkst,gamkst
4116 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4117 * ,ampiz,ampi,amro,gamro,ama1,gama1
4120 * ,amk,amkz,amkst,gamkst
4121 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4122 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4123 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4124 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4126 CHARACTER names(nmode)*31
4127 REAL hv(4),pt(4),pn(4),paa(4),pmult(4,9)
4128 REAL*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4129 REAL*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4130 REAL*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4132 REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4137 DATA pi /3.141592653589793238462643/
4143 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4149 c 6 body phase space normalised as in bjorken-drell
4150 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4151 phspac=1./2**29/pi**14
4152 c phspac=1./2**5/pi**2
4153 c init 5pi decay mode(jnpi)
4154 amp1=dcdmas(idffin(1,jnpi))
4155 amp2=dcdmas(idffin(2,jnpi))
4156 amp3=dcdmas(idffin(3,jnpi))
4157 amp4=dcdmas(idffin(4,jnpi))
4158 amp5=dcdmas(idffin(5,jnpi))
4168 c masses of 5, 4, 3 and 2 pi systems
4169 c 3 pi with sampling for omega resonance
4173 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4174 ams2=(amtau-amnuta)**2
4175 am5sq=ams1+ rr1*(ams2-ams1)
4177 phspac=phspac*(ams2-ams1)
4182 ams1=(amp2+amp3+amp4+amp5)**2
4184 am4sq=ams1+ rr1*(ams2-ams1)
4189 c phase space with sampling for omega resonance
4191 ams1=(amp2+amp3+amp4)**2
4193 alp1=atan((ams1-amom**2)/amom/gamom)
4194 alp2=atan((ams2-amom**2)/amom/gamom)
4195 alp=alp1+rr1*(alp2-alp1)
4196 am3sq =amom**2+amom*gamom*tan(alp)
4198 c --- this part of the jacobian will be recovered later ---------------
4199 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4202 c am3sq=ams1+ rr1*(ams2-ams1)
4204 c --- this part of jacobian will be recovered later
4212 am2sq=ams1+ rr2*(ams2-ams1)
4214 c --- this part of jacobian will be recovered later
4217 c(34) restframe, define pi3 and pi4
4218 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4219 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4220 ppi= enq1**2-amp3**2
4221 pppi=sqrt(abs(enq1**2-amp3**2))
4222 ff1=(4*pi)*(2*pppi/am2)
4223 c pi3 momentum in(34) rest frame
4224 call sphera(pppi,pi3)
4226 c pi4 momentum in(34) rest frame
4231 c(234) rest frame, define pi2
4235 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4236 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4237 ppi = pr(4)**2-am2**2
4241 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4243 c --- this part of jacobian will be recovered later
4244 ff2=(4*pi)*(2*pr(3)/am3)
4245 c old pions boosted from 2 rest frame to 3 rest frame
4246 exe=(pr(4)+pr(3))/am2
4247 call bostr3(exe,pi3,pi3)
4248 call bostr3(exe,pi4,pi4)
4251 thet =acos(-1.+2*rr3)
4253 call rotpol(thet,phi,pi2)
4254 call rotpol(thet,phi,pi3)
4255 call rotpol(thet,phi,pi4)
4257 c(2345) rest frame, define pi5
4261 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4262 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4263 ppi = pr(4)**2-am3**2
4267 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4269 c --- this part of jacobian will be recovered later
4270 ff3=(4*pi)*(2*pr(3)/am4)
4271 c old pions boosted from 3 rest frame to 4 rest frame
4272 exe=(pr(4)+pr(3))/am3
4273 call bostr3(exe,pi2,pi2)
4274 call bostr3(exe,pi3,pi3)
4275 call bostr3(exe,pi4,pi4)
4278 thet =acos(-1.+2*rr3)
4280 call rotpol(thet,phi,pi2)
4281 call rotpol(thet,phi,pi3)
4282 call rotpol(thet,phi,pi4)
4283 call rotpol(thet,phi,pi5)
4285 c(12345) rest frame, define pi1
4289 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4290 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4291 ppi = pr(4)**2-am4**2
4295 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4297 c --- this part of jacobian will be recovered later
4298 ff4=(4*pi)*(2*pr(3)/am5)
4299 c old pions boosted from 4 rest frame to 5 rest frame
4300 exe=(pr(4)+pr(3))/am4
4301 call bostr3(exe,pi2,pi2)
4302 call bostr3(exe,pi3,pi3)
4303 call bostr3(exe,pi4,pi4)
4304 call bostr3(exe,pi5,pi5)
4307 thet =acos(-1.+2*rr3)
4309 call rotpol(thet,phi,pi1)
4310 call rotpol(thet,phi,pi2)
4311 call rotpol(thet,phi,pi3)
4312 call rotpol(thet,phi,pi4)
4313 call rotpol(thet,phi,pi5)
4315 * now to the tau rest frame, define paa and neutrino momenta
4319 c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4320 c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4321 c ppi = paa(4)**2-am5**2
4322 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4323 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4324 ppi = paa(4)**2-am5sq
4325 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4326 * tau-neutrino momentum
4329 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4332 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4334 c all pions boosted from 5 rest frame to tau rest frame
4335 c z-axis antiparallel to neutrino momentum
4336 exe=(paa(4)+paa(3))/am5
4337 call bostr3(exe,pi1,pi1)
4338 call bostr3(exe,pi2,pi2)
4339 call bostr3(exe,pi3,pi3)
4340 call bostr3(exe,pi4,pi4)
4341 call bostr3(exe,pi5,pi5)
4343 c partial width consists of phase space and amplitude
4344 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
4345 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
4349 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4350 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4351 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4352 fompp = cabs(bwign(am3,amom,gamom))**2
4353 c normalisation factor(to some numerical undimensioned factor;
4354 c cf r.fischer et al zphys c3, 313 (1980))
4356 c amplit=ccabib**2*gfermi**2/2. * brak * am5sq*sigee(am5sq,jnpi)
4357 amplit=ccabib**2*gfermi**2/2. * brak
4358 amplit = amplit * fompp * fnorm
4360 c amplit = amplit * fnorm
4361 dgamt=1/(2.*amtau)*amplit*phspac
4374 c missing: transposition of identical particles, startistical factors
4375 c for identical matrices, polarimetric vector. matrix element rather naive.
4376 c flat phase space in pion system + with breit wigner for omega
4377 c anyway it is better than nothing, and code is improvable.
4379 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4380 c ----------------------------------------------------------------------
4381 c it simulates rho decay in tau rest frame with
4382 c z-axis along rho momentum
4383 c rho decays to k kbar
4384 c ----------------------------------------------------------------------
4385 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4386 * ,ampiz,ampi,amro,gamro,ama1,gama1
4387 * ,amk,amkz,amkst,gamkst
4389 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4390 * ,ampiz,ampi,amro,gamro,ama1,gama1
4391 * ,amk,amkz,amkst,gamkst
4392 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4393 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4394 REAL hv(4),pt(4),pn(4),pr(4),pkc(4),pkz(4),qq(4),pmult(4,9)
4396 DATA pi /3.141592653589793238462643/
4399 c three body phase space normalised as in bjorken-drell
4400 phspac=1./2**11/pi**5
4406 c mass of(real/virtual) rho
4408 ams2=(amtau-amnuta)**2
4411 amx2=ams1+ rr1(1)*(ams2-ams1)
4413 phspac=phspac*(ams2-ams1)
4414 c phase space with sampling for rho resonance
4415 c alp1=atan((ams1-amro**2)/amro/gamro)
4416 c alp2=atan((ams2-amro**2)/amro/gamro)
4419 c CALL ranmar(rr1,1)
4420 c alp=alp1+rr1(1)*(alp2-alp1)
4421 c amx2=amro**2+amro*gamro*tan(alp)
4423 c
IF(amx.LT.(amk+amkz)) go to 100
4425 c phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
4426 c phspac=phspac*(alp2-alp1)
4428 c tau-neutrino momentum
4431 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4432 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4436 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4438 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4441 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4442 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4443 pppi=sqrt((enq1-amk)*(enq1+amk))
4444 phspac=phspac*(4*pi)*(2*pppi/amx)
4445 c charged pi momentum in rho rest frame
4446 CALL sphera(pppi,pkc)
4448 c neutral pi momentum in rho rest frame
4452 exe=(pr(4)+pr(3))/amx
4453 c pions boosted from rho rest frame to tau rest frame
4454 CALL bostr3(exe,pkc,pkc)
4455 CALL bostr3(exe,pkz,pkz)
4457 30 qq(i)=pkc(i)-pkz(i)
4458 c qq transverse to pr
4459 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4460 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4462 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4465 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4467 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4468 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4469 & +(gv**2-ga**2)*amtau*amnuta*qq2
4470 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4471 dgamt=1/(2.*amtau)*amplit*phspac
4473 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4481 c ----------------------------------------------------------
4482 c square of pion form factor
4483 c ----------------------------------------------------------
4484 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4485 * ,ampiz,ampi,amro,gamro,ama1,gama1
4486 * ,amk,amkz,amkst,gamkst
4488 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4489 * ,ampiz,ampi,amro,gamro,ama1,gama1
4490 * ,amk,amkz,amkst,gamkst
4493 fpirk=cabs(fpikm(w,amk,amkz))**2
4494 c fpirk=cabs(fpikmk(w,amk,amkz))**2
4496 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4497 c **********************************************************
4499 c **********************************************************
4501 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
4505 c ------------ parameters --------------------
4506 IF (init.EQ.0 )
THEN
4517 c -----------------------------------------------
4519 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4525 c initialize lund
COMMON
4528 SUBROUTINE dwrph(KTO,PHX)
4530 c -------------------------
4532 IMPLICIT REAL*8 (a-h,o-z)
4539 c
CASE of tau radiative decays.
4540 c filling of the lund
COMMON block.
4543 IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
4546 SUBROUTINE dwluph(KTO,PHOT)
4547 c---------------------------------------------------------------------
4548 c lorentz transformation to cmsystem and
4549 c updating of hepevt record
4551 c called by : dexay1,(dekay1,dekay2)
4553 c used when radiative corrections in decays are generated
4554 c---------------------------------------------------------------------
4557 COMMON /taupos/ np1,np2
4560 IF (phot(4).LE.0.0)
RETURN
4562 c position of decaying particle:
4563 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
4570 IF(ktos.GT.10) ktos=ktos-10
4571 c boost and append photon(gamma is 22)
4572 CALL tralo4(ktos,phot,phot,am)
4573 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4578 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4579 c ----------------------------------------------------------------------
4580 c lorentz transformation to cmsystem and
4581 c updating of hepevt record
4583 c isgn = 1/-1 for tau-/tau+
4585 c called by : dexay,(dekay1,dekay2)
4586 c ----------------------------------------------------------------------
4588 REAL pnu(4),pwb(4),pel(4),pne(4)
4589 COMMON /taupos/ np1,np2
4591 c position of decaying particle:
4598 c tau neutrino(nu_tau is 16)
4599 CALL tralo4(kto,pnu,pnu,am)
4600 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4603 CALL tralo4(kto,pwb,pwb,am)
4604 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4606 c electron(e- is 11)
4607 CALL tralo4(kto,pel,pel,am)
4608 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4610 c anti electron neutrino(nu_e is 12)
4611 CALL tralo4(kto,pne,pne,am)
4612 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4616 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4617 c ----------------------------------------------------------------------
4618 c lorentz transformation to cmsystem and
4619 c updating of hepevt record
4621 c isgn = 1/-1 for tau-/tau+
4623 c called by : dexay,(dekay1,dekay2)
4624 c ----------------------------------------------------------------------
4626 REAL pnu(4),pwb(4),pmu(4),pnm(4)
4627 COMMON /taupos/ np1,np2
4629 c position of decaying particle:
4636 c tau neutrino(nu_tau is 16)
4637 CALL tralo4(kto,pnu,pnu,am)
4638 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4641 CALL tralo4(kto,pwb,pwb,am)
4642 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4645 CALL tralo4(kto,pmu,pmu,am)
4646 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
4648 c anti muon neutrino(nu_mu is 14)
4649 CALL tralo4(kto,pnm,pnm,am)
4650 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
4654 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
4655 c ----------------------------------------------------------------------
4656 c lorentz transformation to cmsystem and
4657 c updating of hepevt record
4659 c isgn = 1/-1 for tau-/tau+
4661 c called by : dexay,(dekay1,dekay2)
4662 c ----------------------------------------------------------------------
4665 COMMON /taupos/ np1,np2
4667 c position of decaying particle:
4674 c tau neutrino(nu_tau is 16)
4675 CALL tralo4(kto,pnu,pnu,am)
4676 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4678 c charged pi meson(pi+ is 211)
4679 CALL tralo4(kto,ppi,ppi,am)
4680 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
4684 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
4685 c ----------------------------------------------------------------------
4686 c lorentz transformation to cmsystem and
4687 c updating of hepevt record
4689 c isgn = 1/-1 for tau-/tau+
4691 c called by : dexay,(dekay1,dekay2)
4692 c ----------------------------------------------------------------------
4694 REAL pnu(4),prho(4),pic(4),piz(4)
4695 COMMON /taupos/ np1,np2
4697 c position of decaying particle:
4704 c tau neutrino(nu_tau is 16)
4705 CALL tralo4(kto,pnu,pnu,am)
4706 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4708 c charged rho meson(rho+ is 213)
4709 CALL tralo4(kto,prho,prho,am)
4710 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
4712 c charged pi meson(pi+ is 211)
4713 CALL tralo4(kto,pic,pic,am)
4714 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
4716 c pi0 meson(pi0 is 111)
4717 CALL tralo4(kto,piz,piz,am)
4718 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
4722 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
4723 c ----------------------------------------------------------------------
4724 c lorentz transformation to cmsystem and
4725 c updating of hepevt record
4727 c isgn = 1/-1 for tau-/tau+
4728 c jaa = 1 (2) for a_1- decay to pi+ 2pi- (pi- 2pi0)
4730 c called by : dexay,(dekay1,dekay2)
4731 c ----------------------------------------------------------------------
4733 REAL pnu(4),paa(4),pim1(4),pim2(4),pipl(4)
4734 COMMON /taupos/ np1,np2
4736 c position of decaying particle:
4743 c tau neutrino(nu_tau is 16)
4744 CALL tralo4(kto,pnu,pnu,am)
4745 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4747 c charged a_1 meson(a_1+ is 20213)
4748 CALL tralo4(kto,paa,paa,am)
4749 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
4751 c two possible decays of the charged a1 meson
4754 c a1 --> pi+ pi- pi- (or charged conjugate)
4756 c pi minus(or c.c.) (pi+ is 211)
4757 CALL tralo4(kto,pim2,pim2,am)
4758 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
4760 c pi minus(or c.c.) (pi+ is 211)
4761 CALL tralo4(kto,pim1,pim1,am)
4762 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
4764 c pi plus(or c.c.) (pi+ is 211)
4765 CALL tralo4(kto,pipl,pipl,am)
4766 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
4768 ELSE IF (jaa.EQ.2)
THEN
4770 c a1 --> pi- pi0 pi0(or charged conjugate)
4772 c pi zero(pi0 is 111)
4773 CALL tralo4(kto,pim2,pim2,am)
4774 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
4776 c pi zero(pi0 is 111)
4777 CALL tralo4(kto,pim1,pim1,am)
4778 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
4780 c pi minus(or c.c.) (pi+ is 211)
4781 CALL tralo4(kto,pipl,pipl,am)
4782 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
4788 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
4789 c ----------------------------------------------------------------------
4790 c lorentz transformation to cmsystem and
4791 c updating of hepevt record
4793 c isgn = 1/-1 for tau-/tau+
4795 c ----------------------------------------------------------------------
4798 COMMON /taupos/ np1,np2
4800 c position of decaying particle
4807 c tau neutrino(nu_tau is 16)
4808 CALL tralo4(kto,pnu,pnu,am)
4809 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4811 c k meson(k+ is 321)
4812 CALL tralo4(kto,pkk,pkk,am)
4813 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
4817 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
4818 COMMON / taukle / bra1,brk0,brk0b,brks
4819 REAL*4 bra1,brk0,brk0b,brks
4820 c ----------------------------------------------------------------------
4821 c lorentz transformation to cmsystem and
4822 c updating of hepevt record
4824 c isgn = 1/-1 for tau-/tau+
4825 c jkst=10 (20) corresponds to k0b pi- (k- pi0) decay
4827 c ----------------------------------------------------------------------
4829 REAL pnu(4),pks(4),pkk(4),ppi(4),xio(1)
4830 COMMON /taupos/ np1,np2
4832 c position of decaying particle
4839 c tau neutrino(nu_tau is 16)
4840 CALL tralo4(kto,pnu,pnu,am)
4841 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4843 c charged k* meson(k*+ is 323)
4844 CALL tralo4(kto,pks,pks,am)
4845 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
4847 c two possible decay modes of charged k*
4850 c k*- --> pi- k0b(or charged conjugate)
4852 c charged pi meson(pi+ is 211)
4853 CALL tralo4(kto,ppi,ppi,am)
4854 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
4857 IF (isgn.EQ.-1) bran=brk0
4858 c k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
4860 IF(xio(1).GT.bran)
THEN
4866 CALL tralo4(kto,pkk,pkk,am)
4867 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
4869 ELSE IF(jkst.EQ.20)
THEN
4873 c pi zero(pi0 is 111)
4874 CALL tralo4(kto,ppi,ppi,am)
4875 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
4877 c charged k meson(k+ is 321)
4878 CALL tralo4(kto,pkk,pkk,am)
4879 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
4885 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
4886 c ----------------------------------------------------------------------
4887 c lorentz transformation to cmsystem and
4888 c updating of hepevt record
4890 c isgn = 1/-1 for tau-/tau+
4892 c called by : dexay,(dekay1,dekay2)
4893 c ----------------------------------------------------------------------
4895 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4896 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4898 COMMON /taupos/ np1,np2
4899 CHARACTER names(nmode)*31
4900 REAL pnu(4),pwb(4),pnpi(4,9)
4904 c position of decaying particle
4911 c tau neutrino(nu_tau is 16)
4912 CALL tralo4(kto,pnu,pnu,am)
4913 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4916 CALL tralo4(kto,pwb,pwb,am)
4917 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4919 c multi pi mode jnpi
4921 c get multiplicity of mode jnpi
4924 kfpi=lunpik(idffin(i,jnpi),-isgn)
4925 c for charged conjugate
case, change charged pions only
4926 c
IF(kfpi.NE.111)kfpi=kfpi*isgn
4930 CALL tralo4(kto,ppi,ppi,am)
4931 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
4937 c ----------------------------------------------------------------------
4938 c calculates mass of pp(double precision)
4941 c ----------------------------------------------------------------------
4942 IMPLICIT REAL*8 (a-h,o-z)
4944 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
4946 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
4951 c ******************
4952 c ----------------------------------------------------------------------
4953 c calculates mass of pp
4956 c ----------------------------------------------------------------------
4958 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
4959 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
4964 c ----------------------------------------------------------------------
4966 c used by : koralz radkor
4967 c ----------------------------------------------------------------------
4968 IMPLICIT DOUBLE PRECISION (a-h,o-z)
4969 DATA pi /3.141592653589793238462643d0/
4971 IF(abs(y).LT.abs(x))
THEN
4973 IF(x.LE.0d0) the=pi-the
4975 the=acos(x/sqrt(x**2+y**2))
4981 c ----------------------------------------------------------------------
4982 * calculates angle in(0,2*pi) range out of x-y
4984 c used by : koralz radkor
4985 c ----------------------------------------------------------------------
4986 IMPLICIT DOUBLE PRECISION (a-h,o-z)
4987 DATA pi /3.141592653589793238462643d0/
4989 IF(abs(y).LT.abs(x))
THEN
4991 IF(x.LE.0d0) the=pi-the
4993 the=acos(x/sqrt(x**2+y**2))
4995 IF(y.LT.0d0) the=2d0*pi-the
4998 SUBROUTINE rotod1(PH1,PVEC,QVEC)
4999 c ----------------------------------------------------------------------
5002 c ----------------------------------------------------------------------
5003 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5004 dimension pvec(4),qvec(4),rvec(4)
5012 qvec(2)= cs*rvec(2)-sn*rvec(3)
5013 qvec(3)= sn*rvec(2)+cs*rvec(3)
5017 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5018 c ----------------------------------------------------------------------
5020 c used by : koralz radkor
5021 c ----------------------------------------------------------------------
5022 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5023 dimension pvec(4),qvec(4),rvec(4)
5030 qvec(1)= cs*rvec(1)+sn*rvec(3)
5032 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5036 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5037 c ----------------------------------------------------------------------
5039 c used by : koralz radkor
5040 c ----------------------------------------------------------------------
5041 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5043 dimension pvec(4),qvec(4),rvec(4)
5049 qvec(1)= cs*rvec(1)-sn*rvec(2)
5050 qvec(2)= sn*rvec(1)+cs*rvec(2)
5054 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5055 c ----------------------------------------------------------------------
5056 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5058 c used by : tauola koralz(?)
5059 c ----------------------------------------------------------------------
5060 REAL*4 pvec(4),qvec(4),rvec(4)
5073 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5074 c ----------------------------------------------------------------------
5075 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5077 c used by : koralz radkor
5078 c ----------------------------------------------------------------------
5079 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5080 dimension pvec(4),qvec(4),rvec(4)
5094 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5095 c ----------------------------------------------------------------------
5098 c ----------------------------------------------------------------------
5099 REAL*4 pvec(4),qvec(4),rvec(4)
5107 qvec(2)= cs*rvec(2)-sn*rvec(3)
5108 qvec(3)= sn*rvec(2)+cs*rvec(3)
5111 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5112 c ----------------------------------------------------------------------
5115 c ----------------------------------------------------------------------
5116 IMPLICIT REAL*4(a-h,o-z)
5117 REAL*4 pvec(4),qvec(4),rvec(4)
5124 qvec(1)= cs*rvec(1)+sn*rvec(3)
5126 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5129 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5130 c ----------------------------------------------------------------------
5133 c ----------------------------------------------------------------------
5134 REAL*4 pvec(4),qvec(4),rvec(4)
5140 qvec(1)= cs*rvec(1)-sn*rvec(2)
5141 qvec(2)= sn*rvec(1)+cs*rvec(2)
5145 SUBROUTINE spherd(R,X)
5146 c ----------------------------------------------------------------------
5147 c generates uniformly three-vector x on sphere of radius r
5148 c double precison version of sphera
5149 c ----------------------------------------------------------------------
5150 REAL*8 r,x(4),pi,costh,sinth
5152 DATA pi /3.141592653589793238462643d0/
5156 sinth=sqrt(1 -costh**2)
5157 x(1)=r*sinth*cos(2*pi*rrr(2))
5158 x(2)=r*sinth*sin(2*pi*rrr(2))
5162 SUBROUTINE rotpox(THET,PHI,PP)
5163 IMPLICIT REAL*8 (a-h,o-z)
5164 c ----------------------------------------------------------------------
5166 c ----------------------------------------------------------------------
5169 CALL rotod2(thet,pp,pp)
5170 CALL rotod3( phi,pp,pp)
5173 SUBROUTINE sphera(R,X)
5174 c ----------------------------------------------------------------------
5175 c generates uniformly three-vector x on sphere of radius r
5177 c called by : dphsxx,dadmpi,dadmkk
5178 c ----------------------------------------------------------------------
5181 DATA pi /3.141592653589793238462643/
5185 sinth=sqrt(1.-costh**2)
5186 x(1)=r*sinth*cos(2*pi*rrr(2))
5187 x(2)=r*sinth*sin(2*pi*rrr(2))
5191 SUBROUTINE rotpol(THET,PHI,PP)
5192 c ----------------------------------------------------------------------
5194 c called by : dadmaa,dphsaa
5195 c ----------------------------------------------------------------------
5198 CALL rotor2(thet,pp,pp)
5199 CALL rotor3( phi,pp,pp)
5202 SUBROUTINE ranmar(RVEC,LENV)
5203 c ----------------------------------------------------------------------
5204 c<<<<<
FUNCTION ranmar(IDUMM)
5205 c cernlib v113, version with automatic default initialization
5206 c transformed to
SUBROUTINE to be as in cernlib
5207 c am.lutz november 1988, feb. 1989
5210 c in report fsu-scri-87-50
5211 c modified by f. james, 1988 and 1989, to generate a vector
5212 c of pseudorandom numbers rvec of length lenv, and to put in
5213 c the
COMMON block everything needed to specify currrent state,
5214 c and to add input and output entry points rmarin, rmarut.
5216 c unique random number used in the
program
5217 c ----------------------------------------------------------------------
5218 COMMON / inout / inut,iout
5220 common/raset1/u(97),c,i97,j97
5221 parameter(modcns=1000000000)
5222 DATA ntot,ntot2,ijkl/-1,0,0/
5224 IF (ntot .GE. 0) go to 50
5226 c default initialization. user has called ranmar without rmarin.
5233 entry rmarin(ijklin, ntotin,ntot2n)
5234 c initializing routine for ranmar, may be called before
5235 c generating pseudorandom numbers with ranmar. the input
5236 c values should be in the ranges: 0<=ijklin<=900 ooo ooo
5237 c 0<=ntotin<=999 999 999
5238 c 0<=ntot2n<<999 999 999
5239 c to get the standard values in marsaglia-s paper, ijklin=54217137
5242 ntot = max(ntotin,0)
5243 ntot2= max(ntot2n,0)
5245 c always come here to initialize
5248 kl = ijkl - 30082*ij
5249 i = mod(ij/177, 177) + 2
5250 j = mod(ij, 177) + 2
5251 k = mod(kl/169, 178) + 1
5253 WRITE(iout,201) ijkl,ntot,ntot2
5254 201
FORMAT(1x,
' RANMAR INITIALIZED: ',i10,2x,2i10)
5259 m = mod(mod(i*j,179)*k, 179)
5263 l = mod(53*l+1, 169)
5264 IF (mod(l*m,64) .GE. 32) s = s+t
5269 4 twom24 = 0.5*twom24
5271 cd = 7654321.*twom24
5272 cm = 16777213.*twom24
5275 c complete initialization by skipping
5276 c(ntot2*modcns + ntot) random numbers
5277 DO 45 loop2= 1, ntot2+1
5279 IF (loop2 .EQ. ntot2+1) now=ntot
5280 IF (now .GT. 0)
THEN
5281 WRITE (iout,
'(A,I15)')
' RMARIN SKIPPING OVER ',now
5282 DO 40 idum = 1, ntot
5284 IF (uni .LT. 0.) uni=uni+1.
5287 IF (i97 .EQ. 0) i97=97
5289 IF (j97 .EQ. 0) j97=97
5291 IF (c .LT. 0.) c=c+cm
5295 IF (kalled .EQ. 1)
RETURN
5297 c normal entry to generate lenv random numbers
5299 DO 100 ivec= 1, lenv
5301 IF (uni .LT. 0.) uni=uni+1.
5304 IF (i97 .EQ. 0) i97=97
5306 IF (j97 .EQ. 0) j97=97
5308 IF (c .LT. 0.) c=c+cm
5310 IF (uni .LT. 0.) uni=uni+1.
5311 c replace exact zeroes by uniform distr. *2**-24
5312 IF (uni .EQ. 0.)
THEN
5314 c an exact zero here is very unlikely, but lets be safe.
5315 IF (uni .EQ. 0.) uni= twom24*twom24
5320 IF (ntot .GE. modcns)
THEN
5322 ntot = ntot - modcns
5325 c entry to output current status
5326 entry rmarut(ijklut,ntotut,ntot2t)
5334 IMPLICIT REAL*8(a-h,o-z)
5335 cern c304 version 29/07/71 dilog 59 c
5337 IF(x .LT.-1.0) go to 1
5338 IF(x .LE. 0.5) go to 2
5339 IF(x .EQ. 1.0) go to 3
5340 IF(x .LE. 2.0) go to 4
5344 z=z-0.5* log(abs(x))**2
5350 3 dilogt=1.64493406684822
5354 z=1.64493406684822 - log(x)* log(abs(t))
5355 5 y=2.66666666666666 *t+0.66666666666666
5356 b= 0.00000 00000 00001
5357 a=y*b +0.00000 00000 00004
5358 b=y*a-b+0.00000 00000 00011
5359 a=y*b-a+0.00000 00000 00037
5360 b=y*a-b+0.00000 00000 00121
5361 a=y*b-a+0.00000 00000 00398
5362 b=y*a-b+0.00000 00000 01312
5363 a=y*b-a+0.00000 00000 04342
5364 b=y*a-b+0.00000 00000 14437
5365 a=y*b-a+0.00000 00000 48274
5366 b=y*a-b+0.00000 00001 62421
5367 a=y*b-a+0.00000 00005 50291
5368 b=y*a-b+0.00000 00018 79117
5369 a=y*b-a+0.00000 00064 74338
5370 b=y*a-b+0.00000 00225 36705
5371 a=y*b-a+0.00000 00793 87055
5372 b=y*a-b+0.00000 02835 75385
5373 a=y*b-a+0.00000 10299 04264
5374 b=y*a-b+0.00000 38163 29463
5375 a=y*b-a+0.00001 44963 00557
5376 b=y*a-b+0.00005 68178 22718
5377 a=y*b-a+0.00023 20021 96094
5378 b=y*a-b+0.00100 16274 96164
5379 a=y*b-a+0.00468 63619 59447
5380 b=y*a-b+0.02487 93229 24228
5381 a=y*b-a+0.16607 30329 27855
5382 a=y*a-b+1.93506 43008 6996
5385 c=======================================================================
5386 c===================
END of cpc part ====================================
5387 c=======================================================================
5390 c-----------------------------------------------------------------------
5391 c initialize rchl currents
5392 c dummy routine for compatibility with new updates of tauola
5395 c-----------------------------------------------------------------------
5396 SUBROUTINE inirchl(FLAG)
5401 25
FORMAT(1x,
"TAUOLA IniRChL: Fatal error, FLAG=",i2,
" but RChL currents missing")
5402 WRITE(*,*)
" in loaded version of TAUOLA-FORTRAN library."