1/* copyright(c) 1991-2023 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 <https://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/* glibc
's intent is to support the IEC 559 math functionality, real
28 and complex. If the GCC (4.9 and later) predefined macros
29 specifying compiler intent are available, use them to determine
30 whether the overall intent is to support these features; otherwise,
31 presume an older compiler has intent to support these features and
32 define these macros by default. */
36/* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is
37 synchronized with ISO/IEC 10646:2017, fifth edition, plus
38 the following additions from Amendment 1 to the fifth edition:
41 - 3 additional Zanabazar Square characters */
44C *********************
46C **********************************************************************
48C *********TAUOLA LIBRARY: VERSION 2.7 ******** *
49C **************August 1995****************** *
50C ** AUTHORS: S.JADACH, Z.WAS ***** *
51C ** R. DECKER, M. JEZABEK, J.H.KUEHN, ***** *
52C ********AVAILABLE FROM: WASM AT CERNVM ****** *
53C *******PUBLISHED IN COMP. PHYS. COMM.******** *
54C *** PREPRINT CERN-TH-5856 SEPTEMBER 1990 **** *
55C *** PREPRINT CERN-TH-6195 OCTOBER 1991 **** *
56C *** PREPRINT CERN-TH-6793 NOVEMBER 1992 **** *
57C **********************************************************************
59C ----------------------------------------------------------------------
61C CHOOSES DECAY MODE ACCORDING TO LIST OF BRANCHING RATIOS
72C ----------------------------------------------------------------------
73 COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
77.LE..OR..GT.
IF(NCHAN0NCHAN30) GOTO 902
84.LT.
IF(RRR(1)CUMUL(I)/CUMUL(NCHAN)) JI=I
89 9020 FORMAT(' ----- jaker: wrong nchan
')
92 SUBROUTINE DEKAY(KTO,HX)
93C ***********************
94C THIS DEKAY IS IN SPIRIT OF THE 'decay
' WHICH
95C WAS INCLUDED IN KORAL-B PROGRAM, COMP. PHYS. COMMUN.
96C VOL. 36 (1985) 191, SEE COMMENTS ON GENERAL PHILOSOPHY THERE.
97C KTO=0 INITIALISATION (OBLIGATORY)
98C KTO=1,11 DENOTES TAU+ AND KTO=2,12 TAU-
99C DEKAY(1,H) AND DEKAY(2,H) IS CALLED INTERNALLY BY MC GENERATOR.
100C H DENOTES THE POLARIMETRIC VECTOR, USED BY THE HOST PROGRAM FOR
101C CALCULATION OF THE SPIN WEIGHT.
102C USER MAY OPTIONALLY CALL DEKAY(11,H) DEKAY(12,H) IN ORDER
103C TO TRANSFORM DECAY PRODUCTS TO CMS AND WRITE LUND RECORD IN /LUJETS/.
104C KTO=100, PRINT FINAL REPORT (OPTIONAL).
106C JAK=1 ELECTRON DECAY
114C JAK=0 INCLUSIVE: JAK=1,2,3,4,5,6,7,8
117 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
119 COMMON /TAUPOS/ NP1,NP2
120 COMMON / TAUBMC / GAMPMC(30),GAMPER(30),NEVDEC(30)
121 REAL*4 GAMPMC ,GAMPER
122 PARAMETER (NMODE=15,NM1=0,NM2=1,NM3=8,NM4=2,NM5=1,NM6=3)
123 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
125 CHARACTER NAMES(NMODE)*31
126 COMMON / INOUT / INUT,IOUT
127 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4)
133C INITIALISATION OR REINITIALISATION
134C first or second tau positions in HEPEVT as in KORALB/Z
138.EQ.
IF (IWARM1) X=5/(IWARM-1)
140 WRITE(IOUT,7001) JAK1,JAK2
144.NE..OR..NE.
IF(JAK1-1JAK2-1) THEN
145 CALL DADMEL(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
146 CALL DADMMU(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
147 CALL DADMPI(-1,IDUM,HDUM,PDUM1,PDUM2)
148 CALL DADMRO(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4)
149 CALL DADMAA(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JDUM)
150 CALL DADMKK(-1,IDUM,HDUM,PDUM1,PDUM2)
151 CALL DADMKS(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,JDUM)
152 CALL DADNEW(-1,IDUM,HDUM,PDUM1,PDUM2,PDUMX,JDUM)
158.EQ.
ELSEIF(KTO1) THEN
159C =====================
160C DECAY OF TAU+ IN THE TAU REST FRAME
162.EQ.
IF(IWARM0) GOTO 902
164C AJWMOD to change BRs depending on sign:
166 CALL DEKAY1(0,H,ISGN)
167.EQ.
ELSEIF(KTO2) THEN
168C =================================
169C DECAY OF TAU- IN THE TAU REST FRAME
171.EQ.
IF(IWARM0) GOTO 902
173C AJWMOD to change BRs depending on sign:
175 CALL DEKAY2(0,H,ISGN)
176.EQ.
ELSEIF(KTO11) THEN
177C ======================
178C REST OF DECAY PROCEDURE FOR ACCEPTED TAU+ DECAY
181 CALL DEKAY1(1,H,ISGN)
182.EQ.
ELSEIF(KTO12) THEN
183C ======================
184C REST OF DECAY PROCEDURE FOR ACCEPTED TAU- DECAY
187 CALL DEKAY2(1,H,ISGN)
188.EQ.
ELSEIF(KTO100) THEN
189C =======================
190.NE..OR..NE.
IF(JAK1-1JAK2-1) THEN
191 CALL DADMEL( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
192 CALL DADMMU( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
193 CALL DADMPI( 1,IDUM,HDUM,PDUM1,PDUM2)
194 CALL DADMRO( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4)
195 CALL DADMAA( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JDUM)
196 CALL DADMKK( 1,IDUM,HDUM,PDUM1,PDUM2)
197 CALL DADMKS( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,JDUM)
198 CALL DADNEW( 1,IDUM,HDUM,PDUM1,PDUM2,PDUMX,JDUM)
199 WRITE(IOUT,7010) NEV1,NEV2,NEVTOT
200 WRITE(IOUT,7011) (NEVDEC(I),GAMPMC(I),GAMPER(I),I= 1,7)
202 $ (NEVDEC(I),GAMPMC(I),GAMPER(I),NAMES(I-7),I=8,7+NMODE)
213 7001 FORMAT(///1X,15(5H*****)
214 $ /,' *
', 25X,'*****tauola library: version 2.7 ******
',9X,1H*,
215 $ /,' *
', 25X,'***********august 1995***************
',9X,1H*,
216 $ /,' *
', 25X,'**authors: s.jadach, z.was*************
',9X,1H*,
217 $ /,' *
', 25X,'**r. decker, m. jezabek, j.h.kuehn*****
',9X,1H*,
218 $ /,' *
', 25X,'**available from: wasm at cernvm ******
',9X,1H*,
219 $ /,' *
', 25X,'***** published in comp. phys. comm.***
',9X,1H*,
220 $ /,' *
', 25X,' physics initialization by cleo collab
',9X,1H*,
221 $ /,' *
', 25X,' see alain weinstein www home page:
',9X,1H*,
222 $ /,' *
', 25X,'http://www.cithep.caltech.edu/~ajw/
',9X,1H*,
223 $ /,' *
', 25X,'/korb_doc.html
#files ',9X,1H*,
224 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
225 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
226 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
227 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
228 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
230 7010
FORMAT(///1x,15(5h*****)
231 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
232 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
233 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
234 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
235 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
236 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
237 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
238 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
239 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
240 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
241 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
242 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
243 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
245 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
247 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
248 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
249 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
250 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
251 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
252 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
253 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
255 $ ,i10,2f12.7,a31 ,8x,1h*)
257 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
258 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
261 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
264 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
267 SUBROUTINE dekay1(IMOD,HH,ISGN)
268c *******************************
269c this routine simulates tau+ decay
270 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
271 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
272 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
273 real*4 gampmc ,gamper
275 REAL HV(4),PNU(4),PPI(4)
276 REAL PWB(4),PMU(4),PNM(4)
277 REAL PRHO(4),PIC(4),PIZ(4)
278 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
285 IF(jak1.EQ.-1)
RETURN
290 IF(jak1.EQ.0)
CALL jaker(jak)
292 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
293 ELSEIF(jak.EQ.2)
THEN
294 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
295 ELSEIF(jak.EQ.3)
THEN
296 CALL dadmpi(0, isgn,hv,ppi,pnu)
297 ELSEIF(jak.EQ.4)
THEN
298 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
299 ELSEIF(jak.EQ.5)
THEN
300 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
301 ELSEIF(jak.EQ.6)
THEN
302 CALL dadmkk(0, isgn,hv,pkk,pnu)
303 ELSEIF(jak.EQ.7)
THEN
304 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
306 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
312 ELSEIF(imd.EQ.1)
THEN
313c =====================
316 nevdec(jak)=nevdec(jak)+1
321 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
322 CALL dwrph(ktom,phot)
326 ELSEIF(jak.EQ.2)
THEN
327 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
328 CALL dwrph(ktom,phot)
332 ELSEIF(jak.EQ.3)
THEN
333 CALL dwlupi(1,isgn,ppi,pnu)
337 ELSEIF(jak.EQ.4)
THEN
338 CALL dwluro(1,isgn,pnu,prho,pic,piz)
342 ELSEIF(jak.EQ.5)
THEN
343 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
346 ELSEIF(jak.EQ.6)
THEN
347 CALL dwlukk(1,isgn,pkk,pnu)
350 ELSEIF(jak.EQ.7)
THEN
351 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
356 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
364 SUBROUTINE dekay2(IMOD,HH,ISGN)
365c *******************************
366c this routine simulates tau- decay
367 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
368 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
369 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
370 real*4 gampmc ,gamper
372 REAL HV(4),PNU(4),PPI(4)
373 REAL PWB(4),PMU(4),PNM(4)
374 REAL PRHO(4),PIC(4),PIZ(4)
375 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
382 IF(jak2.EQ.-1)
RETURN
387 IF(jak2.EQ.0)
CALL jaker(jak)
389 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
390 ELSEIF(jak.EQ.2)
THEN
391 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
392 ELSEIF(jak.EQ.3)
THEN
393 CALL dadmpi(0, isgn,hv,ppi,pnu)
394 ELSEIF(jak.EQ.4)
THEN
395 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
396 ELSEIF(jak.EQ.5)
THEN
397 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
398 ELSEIF(jak.EQ.6)
THEN
399 CALL dadmkk(0, isgn,hv,pkk,pnu)
400 ELSEIF(jak.EQ.7)
THEN
401 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
403 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
408 ELSEIF(imd.EQ.1)
THEN
409c =====================
412 nevdec(jak)=nevdec(jak)+1
417 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
418 CALL dwrph(ktom,phot)
422 ELSEIF(jak.EQ.2)
THEN
423 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
424 CALL dwrph(ktom,phot)
428 ELSEIF(jak.EQ.3)
THEN
429 CALL dwlupi(2,isgn,ppi,pnu)
433 ELSEIF(jak.EQ.4)
THEN
434 CALL dwluro(2,isgn,pnu,prho,pic,piz)
438 ELSEIF(jak.EQ.5)
THEN
439 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
442 ELSEIF(jak.EQ.6)
THEN
443 CALL dwlukk(2,isgn,pkk,pnu)
446 ELSEIF(jak.EQ.7)
THEN
447 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
452 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
460 SUBROUTINE dexay(KTO,POL)
461c ----------------------------------------------------------------------
462c this
'DEXAY' is a routine which generates decay of the single
463c polarized tau, pol is a polarization vector(not a polarimeter
464c vector as in dekay) of the tau and it is an input
PARAMETER.
465c kto=0 initialisation(obligatory)
466c kto=1 denotes tau+ and kto=2 tau-
467c dexay(1,pol) and dexay(2,pol) are called internally by mc generator.
468c decay products are transformed readily
469c to cms and writen in the lund record in /lujets/
470c kto=100, print final report(
OPTIONAL).
473c ----------------------------------------------------------------------
474 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
475 real*4 gampmc ,gamper
476 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
478 COMMON /taupos/ np1,np2
479 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
480 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
482 CHARACTER NAMES(NMODE)*31
483 COMMON / inout / inut,iout
485 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
494c initialisation or reinitialisation
495c first or second tau positions in hepevt as in koralb/z
499 WRITE(iout, 7001) jak1,jak2
503 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
504 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
505 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
506 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
507 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
508 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
509 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
510 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
511 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
517 ELSEIF(kto.EQ.1)
THEN
518c =====================
519c decay of tau+ in the tau rest frame
522 IF(iwarm.EQ.0)
GOTO 902
524cam
CALL dexay1(pol,isgn)
525 CALL dexay1(kto,jak1,jakp,pol,isgn)
526 ELSEIF(kto.EQ.2)
THEN
527c =================================
528c decay of tau- in the tau rest frame
531 IF(iwarm.EQ.0)
GOTO 902
532 isgn=-idff/iabs(idff)
533cam
CALL dexay2(pol,isgn)
534 CALL dexay1(kto,jak2,jakm,pol,isgn)
535 ELSEIF(kto.EQ.100)
THEN
536c =======================
537 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
538 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
539 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
540 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
541 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
542 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
543 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
544 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
545 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
546 WRITE(iout,7010) nev1,nev2,nevtot
547 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
549 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
556 7001
FORMAT(///1x,15(5h*****)
557 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
558 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
559 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
560 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
561 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
562 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
563 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
564 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
565 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
566 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
567 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
568 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
569 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
570 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
571 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
572 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
573 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
575chbu
format 7010 had more than 19 continuation lines
577 7010
FORMAT(///1x,15(5h*****)
578 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
579 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
580 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
581 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
582 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
583 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
584 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
585 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
586 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
587 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
588 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
589 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
590 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
592 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
594 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
595 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
596 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
597 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
598 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
599 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
600 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
602 $ ,i10,2f12.7,a31 ,8x,1h*)
604 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
605 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
607 902
WRITE(iout, 9020)
608 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
610 910
WRITE(iout, 9100)
611 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
614 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
615c ---------------------------------------------------------------------
616c this routine simulates tau+- decay
619c ---------------------------------------------------------------------
620 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
621 real*4 gampmc ,gamper
622 COMMON / inout / inut,iout
625 REAL PRHO(4),PIC(4),PIZ(4)
626 REAL PWB(4),PMU(4),PNM(4)
627 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
633 IF(jakin.EQ.-1)
RETURN
640 IF(jak.EQ.0)
CALL jaker(jak)
643 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
644 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
645 CALL dwrph(kto,phot )
646 ELSEIF(jak.EQ.2)
THEN
647 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
648 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
649 CALL dwrph(kto,phot )
650 ELSEIF(jak.EQ.3)
THEN
651 CALL dexpi(0, isgn,polar,ppi,pnu)
652 CALL dwlupi(kto,isgn,ppi,pnu)
653 ELSEIF(jak.EQ.4)
THEN
654 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
655 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
656 ELSEIF(jak.EQ.5)
THEN
657 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
658 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
659 ELSEIF(jak.EQ.6)
THEN
660 CALL dexkk(0, isgn,polar,pkk,pnu)
661 CALL dwlukk(kto,isgn,pkk,pnu)
662 ELSEIF(jak.EQ.7)
THEN
663 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
664 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
667 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
668 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
670 nevdec(jak)=nevdec(jak)+1
672 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
673c ----------------------------------------------------------------------
674c this simulates tau decay in tau rest frame
675c into electron and two neutrinos
677c called by : dexay,dexay1
678c ----------------------------------------------------------------------
679 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
685 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
686cc
CALL hbook1(813,
'WEIGHT DISTRIBUTION DEXEL $',100,0,2)
688 ELSEIF(mode.EQ. 0)
THEN
689c =======================
691 IF(iwarm.EQ.0)
GOTO 902
692 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
693 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
696 IF(rn(1).GT.wt)
GOTO 300
698 ELSEIF(mode.EQ. 1)
THEN
699c =======================
700 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
706 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
709 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
710c ----------------------------------------------------------------------
711c this simulates tau decay in its rest frame
712c into muon and two neutrinos
713c output four momenta: pnu tauneutrino,
717c ----------------------------------------------------------------------
718 COMMON / inout / inut,iout
719 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
725 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
726cc
CALL hbook1(814,
'WEIGHT DISTRIBUTION DEXMU $',100,0,2)
728 ELSEIF(mode.EQ. 0)
THEN
729c =======================
731 IF(iwarm.EQ.0)
GOTO 902
732 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
733 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
736 IF(rn(1).GT.wt)
GOTO 300
738 ELSEIF(mode.EQ. 1)
THEN
739c =======================
740 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
745 902
WRITE(iout, 9020)
746 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
749 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
750c ----------------------------------------------------------------------
752c called by : dexel,(dekay,dekay1)
753c ----------------------------------------------------------------------
754 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
755 real*4 gfermi,gv,ga,ccabib,scabib,gamel
756 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
757 * ,ampiz,ampi,amro,gamro,ama1,gama1
758 * ,amk,amkz,amkst,gamkst
760 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
761 * ,ampiz,ampi,amro,gamro,ama1,gama1
762 * ,amk,amkz,amkst,gamkst
763 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
764 real*4 gampmc ,gamper
766 COMMON / inout / inut,iout
767 REAL HHV(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4)
768 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
771 DATA pi /3.141592653589793238462643/
784 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
785 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
787cc
CALL hbook1(803,
'WEIGHT DISTRIBUTION DADMEL $',100,0,2)
789 ELSEIF(mode.EQ. 0)
THEN
790c =======================
792 IF(iwarm.EQ.0)
GOTO 902
794 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
795cc
CALL hfill(803,wt/wtmax)
800 IF(wt.GT.wtmax) nevovr=nevovr+1
801 IF(rn*wtmax.GT.wt)
GOTO 300
802c rotations to basic tau rest frame
808 CALL rotor2(thet,pnu,pnu)
809 CALL rotor3( phi,pnu,pnu)
810 CALL rotor2(thet,pwb,pwb)
811 CALL rotor3( phi,pwb,pwb)
812 CALL rotor2(thet,q1,q1)
813 CALL rotor3( phi,q1,q1)
814 CALL rotor2(thet,q2,q2)
815 CALL rotor3( phi,q2,q2)
816 CALL rotor2(thet,hv,hv)
817 CALL rotor3( phi,hv,hv)
818 CALL rotor2(thet,phx,phx)
819 CALL rotor3( phi,phx,phx)
821 44 hhv(i)=-isgn*hv(i)
824 ELSEIF(mode.EQ. 1)
THEN
825c =======================
826 IF(nevraw.EQ.0)
RETURN
827 pargam=swt/float(nevraw+1)
829 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
831 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
839 7010
FORMAT(///1x,15(5h*****)
840 $ /,
' *', 25x,
'******** DADMEL FINAL REPORT ******** ',9x,1h*
841 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
842 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
843 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
844 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
845 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
846 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
847 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
848 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
850 902
WRITE(iout, 9020)
851 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
854 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
855c ----------------------------------------------------------------------
856 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
857 * ,ampiz,ampi,amro,gamro,ama1,gama1
858 * ,amk,amkz,amkst,gamkst
860 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
861 * ,ampiz,ampi,amro,gamro,ama1,gama1
862 * ,amk,amkz,amkst,gamkst
863 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
864 real*4 gfermi,gv,ga,ccabib,scabib,gamel
865 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
866 real*4 gampmc ,gamper
867 COMMON / inout / inut,iout
869 REAL HHV(4),HV(4),PNU(4),PWB(4),Q1(4),Q2(4)
870 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
873 DATA pi /3.141592653589793238462643/
886 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
887 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
889cc
CALL hbook1(802,
'WEIGHT DISTRIBUTION DADMMU $',100,0,2)
891 ELSEIF(mode.EQ. 0)
THEN
892c =======================
894 IF(iwarm.EQ.0)
GOTO 902
896 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
897cc
CALL hfill(802,wt/wtmax)
902 IF(wt.GT.wtmax) nevovr=nevovr+1
903 IF(rn*wtmax.GT.wt)
GOTO 300
904c rotations to basic tau rest frame
908 CALL rotor2(thet,pnu,pnu)
909 CALL rotor3( phi,pnu,pnu)
910 CALL rotor2(thet,pwb,pwb)
911 CALL rotor3( phi,pwb,pwb)
912 CALL rotor2(thet,q1,q1)
913 CALL rotor3( phi,q1,q1)
914 CALL rotor2(thet,q2,q2)
915 CALL rotor3( phi,q2,q2)
916 CALL rotor2(thet,hv,hv)
917 CALL rotor3( phi,hv,hv)
918 CALL rotor2(thet,phx,phx)
919 CALL rotor3( phi,phx,phx)
921 44 hhv(i)=-isgn*hv(i)
924 ELSEIF(mode.EQ. 1)
THEN
925c =======================
926 IF(nevraw.EQ.0)
RETURN
927 pargam=swt/float(nevraw+1)
929 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
931 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
939 7010
FORMAT(///1x,15(5h*****)
940 $ /,
' *', 25x,
'******** DADMMU FINAL REPORT ******** ',9x,1h*
941 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
942 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
943 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
944 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
945 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
946 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
947 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
948 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
950 902
WRITE(iout, 9020)
951 9020
FORMAT(
' ----- DADMMU: LACK OF INITIALISATION')
954 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
955c xnx,xna was flipped in parameters of dphsel and dphsmu
956c *********************************************************************
957c * electron decay mode *
958c *********************************************************************
960 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
961 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
964 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
975 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
976c xnx,xna was flipped in parameters of dphsel and dphsmu
977c *********************************************************************
979c *********************************************************************
981 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
982 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
985 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
996 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
997 IMPLICIT REAL*8 (a-h,o-z)
998c ----------------------------------------------------------------------
999* it simulates e,mu channels of tau decay in its rest frame with
1000* qed order alpha corrections
1001c ----------------------------------------------------------------------
1002 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1003 * ,ampiz,ampi,amro,gamro,ama1,gama1
1004 * ,amk,amkz,amkst,gamkst
1006 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1007 * ,ampiz,ampi,amro,gamro,ama1,gama1
1008 * ,amk,amkz,amkst,gamkst
1009 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1010 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1011 COMMON / inout / inut,iout
1012 COMMON / taurad / xk0dec,itdkrc
1014 real*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1018 DATA pi /3.141592653589793238462643d0/
1019c ajwmod to satisfy compiler, comment out this unused function.
1020c amro, gamro is only a
PARAMETER for geting hight efficiency
1022c three body phase space normalised as in bjorken-drell
1023c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
1024 phspac=1./2**17/pi**8
1034 IF (ielmu.EQ.1)
THEN
1041 IF ( itdkrc.EQ.0) prhard=0d0
1043 IF(prsoft.LT.0.1)
THEN
1044 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1049 ihard=(rr5.GT.prsoft)
1051c tau decay to
'TAU+photon'
1053 ams1=(amu+amnuta)**2
1056 xl1=log(xk1/2/xk0dec)
1061 phspac=phspac*ams2*xl1*xk
1062 phspac=phspac/prhard
1065 phspac=phspac*2**6*pi**3
1066 phspac=phspac/prsoft
1068c mass of neutrina system
1075 am2sq=ams1+ rr2*(ams2-ams1)
1077 phspac=phspac*(ams2-ams1)
1078* neutrina rest frame, define xn and xa
1079 enq1=(am2sq+amnuta**2)/(2*am2)
1080 enq2=(am2sq-amnuta**2)/(2*am2)
1081 ppi= enq1**2-amnuta**2
1082 pppi=sqrt(abs(enq1**2-amnuta**2))
1083 phspac=phspac*(4*pi)*(2*pppi/am2)
1084* nu tau in nunu rest frame
1085 CALL spherd(pppi,xn)
1087* nu light in nunu rest frame
1091* tau-prim rest frame, define qp(muon
1095 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1096 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1097 ppi = pr(4)**2-am2**2
1101 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1103 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1104* neutrina boosted from their frame to tau-prim rest frame
1105 exe=(pr(4)+pr(3))/am2
1106 CALL bostd3(exe,xn,xn)
1107 CALL bostd3(exe,xa,xa)
1111 eps=4*(amu/amtax)**2
1112 xl1=log((2+eps)/eps)
1114 eta =exp(xl1*rr3+xl0)
1117 phspac=phspac*xl1/2*eta
1119 CALL rotpox(thet,phi,xn)
1120 CALL rotpox(thet,phi,xa)
1121 CALL rotpox(thet,phi,qp)
1122 CALL rotpox(thet,phi,pr)
1124* now to the tau rest frame, define tau-prim and gamma momenta
1128 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1129 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1130 ppi = paa(4)**2-am3**2
1131 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1137* all momenta boosted from tau-prim rest frame to tau rest frame
1138* z-axis antiparallel to photon momentum
1139 exe=(paa(4)+paa(3))/am3
1140 CALL bostd3(exe,xn,xn)
1141 CALL bostd3(exe,xa,xa)
1142 CALL bostd3(exe,qp,qp)
1143 CALL bostd3(exe,pr,pr)
1145 thet =acos(-1.+2*rr3)
1147 CALL rotpox(thet,phi,xn)
1148 CALL rotpox(thet,phi,xa)
1149 CALL rotpox(thet,phi,qp)
1150 CALL rotpox(thet,phi,pr)
1152* now to the tau rest frame, define tau-prim and gamma momenta
1164c partial width consists of phase space and amplitude
1165 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1166 dgamt=1/(2.*amtax)*amplit*phspac
1168 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1169 IMPLICIT REAL*8 (a-h,o-z)
1170c ----------------------------------------------------------------------
1171c it calculates matrix element for the
1172c tau --> mu(e) nu nubar decay mode
1173c including complete order alpha qed corrections.
1174c ----------------------------------------------------------------------
1175 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1176 * ,ampiz,ampi,amro,gamro,ama1,gama1
1177 * ,amk,amkz,amkst,gamkst
1179 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1180 * ,ampiz,ampi,amro,gamro,ama1,gama1
1181 * ,amk,amkz,amkst,gamkst
1182 real*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1186 IF(xk(4).LT.0.1d0*ak0)
THEN
1187 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1189 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1193 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1195c **********************************************************************
1196c real photon matrix element squared *
1198c hv- polarimetric four-vector of tau *
1199c qp,xn,xa,xk - 4-momenta of electron(muon), nu, nubar and photon *
1200c all four-vectors in tau rest frame(in gev) *
1201c ak0 - infrared cutoff, minimal energy of hard photons(gev) *
1202c sqm2 -
value for s=0 *
1203c see eqs. (2.9)-(2.10) from cjk( nucl.phys.b(1991) ) *
1204c **********************************************************************
1206 IMPLICIT REAL*8(a-h,o-z)
1207 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1208 * ,ampiz,ampi,amro,gamro,ama1,gama1
1209 * ,amk,amkz,amkst,gamkst
1211 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1212 * ,ampiz,ampi,amro,gamro,ama1,gama1
1213 * ,amk,amkz,amkst,gamkst
1214 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1215 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1216 COMMON / qedprm /alfinv,alfpi,xk0
1217 real*8 alfinv,alfpi,xk0
1218 real*8 qp(4),xn(4),xa(4),xk(4)
1221 real*8 s0(3),rxa(3),rxk(3),rqp(3)
1222 DATA pi /3.141592653589793238462643d0/
1228 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1230c scalar products of four-momenta
1236 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1237c rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1238 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1239 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1241 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1242 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1243 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1244c xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1245 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1246 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1256 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1258 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1259 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1260 const4=256*pi/alphai*gf**2
1261 IF (itdkrc.EQ.0) const4=0d0
1264 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1266 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1267 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1268 5 hv(i)=s0(i)/s1-1.d0
1271 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1273c **********************************************************************
1274c born +virtual+soft photon matrix element**2 o(alpha) *
1276c hv- polarimetric four-vector of tau *
1277c qp,xn,xa - four-momenta of electron(muon), nu and nubar in gev *
1278c all four-vectors in tau rest frame *
1279c ak0 - infrared cutoff, minimal energy of hard photons *
1280c thb -
VALUE for s=0 *
1281c see eqs. (2.2),(2.4)-(2.5) from cjk(nucl.phys.b351(1991)70 *
1282c and(c.2) from jk(nucl.phys.b320(1991)20 ) *
1283c **********************************************************************
1285 IMPLICIT REAL*8(a-h,o-z)
1286 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1287 * ,ampiz,ampi,amro,gamro,ama1,gama1
1288 * ,amk,amkz,amkst,gamkst
1290 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1291 * ,ampiz,ampi,amro,gamro,ama1,gama1
1292 * ,amk,amkz,amkst,gamkst
1293 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1294 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1295 COMMON / qedprm /alfinv,alfpi,xk0
1296 real*8 alfinv,alfpi,xk0
1297 dimension qp(4),xn(4),xa(4)
1300 real*8 rxa(3),rxn(3),rqp(3)
1301 real*8 bornpl(3),am3pol(3),xm3pol(3)
1302 DATA pi /3.141592653589793238462643d0/
1315 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1316 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1317c rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1318 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1320c quasi two-body variables
1322 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1324 w0=(xn(4)+xa(4))/tmass
1336 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1337 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1338 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1339 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1340 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1341 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1343c scalar products of four-momenta
1344 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1345 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1346 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1350c decay differential width without and with polarization
1351 const3=1/(2*alphai*pi)*64*gf**2
1352 IF (itdkrc.EQ.0) const3=0d0
1353 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1354 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1356c v-a and v+a couplings, but in the born part only
1357 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1358 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1359 born= 32*(gfermi**2/2.)*brak
1361 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1362 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1363 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1364 am3pol(i)=xm3pol(i)*const3
1365c v-a and v+a couplings, but in the born part only
1367 & (gv+ga)**2*tmass*xnxa*qp(i)
1368 & -(gv-ga)**2*tmass*qpxn*xa(i)
1369 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1370 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1372 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1374 IF (thb/born.LT.0.1d0)
THEN
1375 print *,
'ERROR IN THB, THB/BORN=',thb/born
1380 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1381c ----------------------------------------------------------------------
1382c tau decay into pion and tau-neutrino
1384c output four momenta: pnu tauneutrino,
1386c ----------------------------------------------------------------------
1387 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
1390c ===================
1391 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1392cc
CALL hbook1(815,
'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1394 ELSEIF(mode.EQ. 0)
THEN
1395c =======================
1397 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1398 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1399cc
CALL hfill(815,wt)
1401 IF(rn(1).GT.wt)
GOTO 300
1403 ELSEIF(mode.EQ. 1)
THEN
1404c =======================
1405 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1411 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1412c ----------------------------------------------------------------------
1413 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1414 * ,ampiz,ampi,amro,gamro,ama1,gama1
1415 * ,amk,amkz,amkst,gamkst
1417 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1418 * ,ampiz,ampi,amro,gamro,ama1,gama1
1419 * ,amk,amkz,amkst,gamkst
1420 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1421 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1422 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1423 real*4 gampmc ,gamper
1424 COMMON / inout / inut,iout
1425 REAL PPI(4),PNU(4),HV(4)
1426 DATA pi /3.141592653589793238462643/
1429c ===================
1431 ELSEIF(mode.EQ. 0)
THEN
1432c =======================
1434 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1435 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1436 xpi= sqrt(epi**2-ampi**2)
1438 CALL sphera(xpi,ppi)
1440c tau-neutrino momentum
1446 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1447 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1448 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
145040 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1453 ELSEIF(mode.EQ. 1)
THEN
1454c =======================
1455 IF(nevtot.EQ.0)
RETURN
1457c gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1458c * (brak/amtau**4)**2
1459czw 7.02.93 here was an error affecting non standard model
1460c configurations only
1461 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1463 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1464 $ -4*ampi**2*amnuta**2 )/amtau**2
1467 WRITE(iout, 7010) nevtot,gamm,rat,error
1474 7010
FORMAT(///1x,15(5h*****)
1475 $ /,
' *', 25x,
'******** DADMPI FINAL REPORT ******** ',9x,1h*
1476 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1477 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1478 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1479 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1480 $ /,1x,15(5h*****)/)
1482 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1483c ----------------------------------------------------------------------
1484c this simulates tau decay in tau rest frame
1485c into nu rho, next rho decays into pion pair.
1486c output four momenta: pnu tauneutrino,
1490c ----------------------------------------------------------------------
1491 COMMON / inout / inut,iout
1492 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
1496c ===================
1498 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1499cc
CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXRO $',100,0,2)
1500cc
CALL hbook1(916,
'ABS2 OF HV IN ROUTINE DEXRO $',100,0,2)
1502 ELSEIF(mode.EQ. 0)
THEN
1503c =======================
1505 IF(iwarm.EQ.0)
GOTO 902
1506 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1507 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1508cc
CALL hfill(816,wt)
1509cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
1510cc
CALL hfill(916,xhelp)
1512 IF(rn(1).GT.wt)
GOTO 300
1514 ELSEIF(mode.EQ. 1)
THEN
1515c =======================
1516 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1522 902
WRITE(iout, 9020)
1523 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1526 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1527c ----------------------------------------------------------------------
1528 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1529 * ,ampiz,ampi,amro,gamro,ama1,gama1
1530 * ,amk,amkz,amkst,gamkst
1532 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1533 * ,ampiz,ampi,amro,gamro,ama1,gama1
1534 * ,amk,amkz,amkst,gamkst
1535 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1536 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1537 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1538 real*4 gampmc ,gamper
1539 COMMON / inout / inut,iout
1541 REAL HV(4),PRO(4),PNU(4),PIC(4),PIZ(4)
1542 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
1545 DATA pi /3.141592653589793238462643/
1549c ===================
1558 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1559 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1561cc
CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMRO $',100,0,2)
1564 ELSEIF(mode.EQ. 0)
THEN
1565c =======================
1567 IF(iwarm.EQ.0)
GOTO 902
1568 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1569cc
CALL hfill(801,wt/wtmax)
1575 IF(wt.GT.wtmax) nevovr=nevovr+1
1576 IF(rn*wtmax.GT.wt)
GOTO 300
1577c rotations to basic tau rest frame
1578 costhe=-1.+2.*rrr(2)
1581 CALL rotor2(thet,pnu,pnu)
1582 CALL rotor3( phi,pnu,pnu)
1583 CALL rotor2(thet,pro,pro)
1584 CALL rotor3( phi,pro,pro)
1585 CALL rotor2(thet,pic,pic)
1586 CALL rotor3( phi,pic,pic)
1587 CALL rotor2(thet,piz,piz)
1588 CALL rotor3( phi,piz,piz)
1589 CALL rotor2(thet,hv,hv)
1590 CALL rotor3( phi,hv,hv)
1592 44 hhv(i)=-isgn*hv(i)
1595 ELSEIF(mode.EQ. 1)
THEN
1596c =======================
1597 IF(nevraw.EQ.0)
RETURN
1598 pargam=swt/float(nevraw+1)
1600 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1602 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1610 7003
FORMAT(///1x,15(5h*****)
1611 $ /,
' *', 25x,
'******** DADMRO INITIALISATION ********',9x,1h*
1612 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1613 $ /,1x,15(5h*****)/)
1614 7010
FORMAT(///1x,15(5h*****)
1615 $ /,
' *', 25x,
'******** DADMRO FINAL REPORT ******** ',9x,1h*
1616 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1617 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1618 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1619 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1620 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1621 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1622 $ /,1x,15(5h*****)/)
1623 902
WRITE(iout, 9020)
1624 9020
FORMAT(
' ----- DADMRO: LACK OF INITIALISATION')
1627 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1628c ----------------------------------------------------------------------
1629c it simulates rho decay in tau rest frame with
1630c z-axis along rho momentum
1631c ----------------------------------------------------------------------
1632 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1633 * ,ampiz,ampi,amro,gamro,ama1,gama1
1634 * ,amk,amkz,amkst,gamkst
1636 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1637 * ,ampiz,ampi,amro,gamro,ama1,gama1
1638 * ,amk,amkz,amkst,gamkst
1639 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1640 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1641 REAL HV(4),PT(4),PN(4),PR(4),PIC(4),PIZ(4),QQ(4),RR1(1)
1642 DATA pi /3.141592653589793238462643/
1645c three body phase space normalised as in bjorken-drell
1646 phspac=1./2**11/pi**5
1652c mass of(real/virtual) rho
1653 ams1=(ampi+ampiz)**2
1654 ams2=(amtau-amnuta)**2
1656c amx2=ams1+ rr1*(ams2-ams1)
1658c phspac=phspac*(ams2-ams1)
1659c phase space with sampling for rho resonance
1660 alp1=atan((ams1-amro**2)/amro/gamro)
1661 alp2=atan((ams2-amro**2)/amro/gamro)
1665 alp=alp1+rr1(1)*(alp2-alp1)
1666 amx2=amro**2+amro*gamro*tan(alp)
1668 IF(amx.LT.2.*ampi)
GO TO 100
1670 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1671 phspac=phspac*(alp2-alp1)
1673c tau-neutrino momentum
1676 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1677 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1681 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1683 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1686 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1687 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1688 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1689 phspac=phspac*(4*pi)*(2*pppi/amx)
1690c charged pi momentum in rho rest frame
1691 CALL sphera(pppi,pic)
1693c neutral pi momentum in rho rest frame
1697 exe=(pr(4)+pr(3))/amx
1698c pions boosted from rho rest frame to tau rest frame
1699 CALL bostr3(exe,pic,pic)
1700 CALL bostr3(exe,piz,piz)
170230 qq(i)=pic(i)-piz(i)
1705 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
1707 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
1708 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
1709 & +(gv**2-ga**2)*amtau*amnuta*qq2
1710 amplit=(gfermi*ccabib)**2*brak*2*fpirho(amx)
1711 dgamt=1/(2.*amtau)*amplit*phspac
1713 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1716 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1717c ----------------------------------------------------------------------
1718* this simulates tau decay in tau rest frame
1719* into nu a1, next a1 decays into rho pi and finally rho into pi pi.
1720* output four momenta: pnu tauneutrino,
1722* pim1 pion minus(or pi0) 1 (for tau minus)
1723* pim2 pion minus(or pi0) 2
1724* pipl pion plus(or pi-)
1725* (pipl,pim1) form a rho
1726c ----------------------------------------------------------------------
1727 COMMON / inout / inut,iout
1728 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
1732c ===================
1734 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1735cc
CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
1737 ELSEIF(mode.EQ. 0)
THEN
1738* =======================
1740 IF(iwarm.EQ.0)
GOTO 902
1741 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1742 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1743cc
CALL hfill(816,wt)
1745 IF(rn(1).GT.wt)
GOTO 300
1747 ELSEIF(mode.EQ. 1)
THEN
1748* =======================
1749 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1754 902
WRITE(iout, 9020)
1755 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
1758 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1759c ----------------------------------------------------------------------
1760* a1 decay unweighted events
1761c ----------------------------------------------------------------------
1762 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1763 * ,ampiz,ampi,amro,gamro,ama1,gama1
1764 * ,amk,amkz,amkst,gamkst
1766 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1767 * ,ampiz,ampi,amro,gamro,ama1,gama1
1768 * ,amk,amkz,amkst,gamkst
1769 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1770 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1771 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1772 real*4 gampmc ,gamper
1773 COMMON / inout / inut,iout
1775 REAL HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4)
1776 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
1779 DATA pi /3.141592653589793238462643/
1783c ===================
1792 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1793 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1795cc
CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMAA $',100,0,2)
1797 ELSEIF(mode.EQ. 0)
THEN
1798c =======================
1800 IF(iwarm.EQ.0)
GOTO 902
1801 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1802cc
CALL hfill(801,wt/wtmax)
1807 sswt=sswt+dble(wt)**2
1811 IF(wt.GT.wtmax) nevovr=nevovr+1
1812 IF(rn*wtmax.GT.wt)
GOTO 300
1813c rotations to basic tau rest frame
1814 costhe=-1.+2.*rrr(2)
1817 CALL rotpol(thet,phi,pnu)
1818 CALL rotpol(thet,phi,paa)
1819 CALL rotpol(thet,phi,pim1)
1820 CALL rotpol(thet,phi,pim2)
1821 CALL rotpol(thet,phi,pipl)
1822 CALL rotpol(thet,phi,hv)
1824 44 hhv(i)=-isgn*hv(i)
1827 ELSEIF(mode.EQ. 1)
THEN
1828c =======================
1829 IF(nevraw.EQ.0)
RETURN
1830 pargam=swt/float(nevraw+1)
1832 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1834 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1842 7003
FORMAT(///1x,15(5h*****)
1843 $ /,
' *', 25x,
'******** DADMAA INITIALISATION ********',9x,1h*
1844 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1845 $ /,1x,15(5h*****)/)
1846 7010
FORMAT(///1x,15(5h*****)
1847 $ /,
' *', 25x,
'******** DADMAA FINAL REPORT ******** ',9x,1h*
1848 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
1849 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
1850 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1851 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
1852 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1853 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1854 $ /,1x,15(5h*****)/)
1855 902
WRITE(iout, 9020)
1856 9020
FORMAT(
' ----- DADMAA: LACK OF INITIALISATION')
1859 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1860c ----------------------------------------------------------------------
1861* it simulates a1 decay in tau rest frame with
1862* z-axis along a1 momentum
1863c ----------------------------------------------------------------------
1864 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1865 * ,ampiz,ampi,amro,gamro,ama1,gama1
1866 * ,amk,amkz,amkst,gamkst
1868 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1869 * ,ampiz,ampi,amro,gamro,ama1,gama1
1870 * ,amk,amkz,amkst,gamkst
1871 COMMON / taukle / bra1,brk0,brk0b,brks
1872 real*4 bra1,brk0,brk0b,brks
1873 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
1877c matrix element number:
1879c
TYPE of the generation:
1883 IF (rmod.LT.bra1)
THEN
1895 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
1897 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
1898c ----------------------------------------------------------------------
1899c tau decay into kaon and tau-neutrino
1901c output four momenta: pnu tauneutrino,
1903c ----------------------------------------------------------------------
1904 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
1907c ===================
1908 CALL dadmkk(-1,isgn,hv,pkk,pnu)
1909cc
CALL hbook1(815,
'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1911 ELSEIF(mode.EQ. 0)
THEN
1912c =======================
1914 CALL dadmkk( 0,isgn,hv,pkk,pnu)
1915 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1916cc
CALL hfill(815,wt)
1918 IF(rn(1).GT.wt)
GOTO 300
1920 ELSEIF(mode.EQ. 1)
THEN
1921c =======================
1922 CALL dadmkk( 1,isgn,hv,pkk,pnu)
1928 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
1929c ----------------------------------------------------------------------
1931 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1932 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1933 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1934 * ,ampiz,ampi,amro,gamro,ama1,gama1
1935 * ,amk,amkz,amkst,gamkst
1937 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1938 * ,ampiz,ampi,amro,gamro,ama1,gama1
1939 * ,amk,amkz,amkst,gamkst
1940 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1941 real*4 gampmc ,gamper
1942 COMMON / inout / inut,iout
1943 REAL PKK(4),PNU(4),HV(4)
1944 DATA pi /3.141592653589793238462643/
1947c ===================
1949 ELSEIF(mode.EQ. 0)
THEN
1950c =======================
1952 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
1953 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
1954 xkk= sqrt(ekk**2-amk**2)
1956 CALL sphera(xkk,pkk)
1958c tau-neutrino momentum
1964 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
1965 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
1966 & +(gv**2-ga**2)*amtau*amnuta*amk**2
196840 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
1971 ELSEIF(mode.EQ. 1)
THEN
1972c =======================
1973 IF(nevtot.EQ.0)
RETURN
1975cfz there was brak/amtau**4 before
1976c gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
1977c * (brak/amtau**4)**2
1978czw 7.02.93 here was an error affecting non standard model
1979c configurations only
1980 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
1982 $ sqrt((amtau**2-amk**2-amnuta**2)**2
1983 $ -4*amk**2*amnuta**2 )/amtau**2
1988 WRITE(iout, 7010) nevtot,gamm,rat,error
1995 7010
FORMAT(///1x,15(5h*****)
1996 $ /,
' *', 25x,
'******** DADMKK FINAL REPORT ********',9x,1h*
1997 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
1998 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
1999 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2000 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2001 $ /,1x,15(5h*****)/)
2003 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2004c ----------------------------------------------------------------------
2005c this simulates tau decay in tau rest frame
2006c into nu k*,
THEN k* decays into pi0,k+-(jkst=20)
2007c or pi+-,k0(jkst=10).
2008c output four momenta: pnu tauneutrino,
2014c ----------------------------------------------------------------------
2015 COMMON / inout / inut,iout
2016 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
2020c ===================
2022cfz initialisation done with the gharged pion neutral kaon mode(jkst=10
2023 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2024cc
CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXKS $',100,0,2)
2025cc
CALL hbook1(916,
'ABS2 OF HV IN ROUTINE DEXKS $',100,0,2)
2027 ELSEIF(mode.EQ. 0)
THEN
2028c =======================
2030 IF(iwarm.EQ.0)
GOTO 902
2031 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2032 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2033cc
CALL hfill(816,wt)
2034cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
2035cc
CALL hfill(916,xhelp)
2037 IF(rn(1).GT.wt)
GOTO 300
2039 ELSEIF(mode.EQ. 1)
THEN
2040c ======================================
2041 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2047 902
WRITE(iout, 9020)
2048 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2051 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2052c ----------------------------------------------------------------------
2053 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2054 * ,ampiz,ampi,amro,gamro,ama1,gama1
2055 * ,amk,amkz,amkst,gamkst
2057 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2058 * ,ampiz,ampi,amro,gamro,ama1,gama1
2059 * ,amk,amkz,amkst,gamkst
2060 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2061 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2062 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2063 real*4 gampmc ,gamper
2064 COMMON / taukle / bra1,brk0,brk0b,brks
2065 real*4 bra1,brk0,brk0b,brks
2066 COMMON / inout / inut,iout
2068 REAL HV(4),PKS(4),PNU(4),PKK(4),PPI(4)
2069 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
2070 real*4 rrr(3),rmod(1)
2072 DATA pi /3.141592653589793238462643/
2076c ===================
2085c the initialisation is done with the 66.7% MODE
2087 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2088 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2090cc
CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMKS $',100,0,2)
2092cc
CALL hbook1(112,
'-------- K* MASS -------- $',100,0.,2.)
2093 ELSEIF(mode.EQ. 0)
THEN
2094c =====================================
2095 IF(iwarm.EQ.0)
GOTO 902
2096c here we choose randomly between k0 pi+_(66.7%)
2101 IF(rmod(1).LT.dec1)
THEN
2106 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2109 IF(wt.GT.wtmax) nevovr=nevovr+1
2113 IF(rn*wtmax.GT.wt)
GOTO 400
2114c rotations to basic tau rest frame
2115 costhe=-1.+2.*rrr(2)
2118 CALL rotor2(thet,pnu,pnu)
2119 CALL rotor3( phi,pnu,pnu)
2120 CALL rotor2(thet,pks,pks)
2121 CALL rotor3( phi,pks,pks)
2122 CALL rotor2(thet,pkk,pkk)
2123 CALL rotor3(phi,pkk,pkk)
2124 CALL rotor2(thet,ppi,ppi)
2125 CALL rotor3( phi,ppi,ppi)
2126 CALL rotor2(thet,hv,hv)
2127 CALL rotor3( phi,hv,hv)
2129 44 hhv(i)=-isgn*hv(i)
2132 ELSEIF(mode.EQ. 1)
THEN
2133c =======================
2134 IF(nevraw.EQ.0)
RETURN
2135 pargam=swt/float(nevraw+1)
2137 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2139 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2147 7003
FORMAT(///1x,15(5h*****)
2148 $ /,
' *', 25x,
'******** DADMKS INITIALISATION ********',9x,1h*
2149 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2150 $ /,1x,15(5h*****)/)
2151 7010
FORMAT(///1x,15(5h*****)
2152 $ /,
' *', 25x,
'******** DADMKS FINAL REPORT ********',9x,1h*
2153 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2154 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2155 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2156 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2157 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2158 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2159 $ /,1x,15(5h*****)/)
2160 902
WRITE(iout, 9020)
2161 9020
FORMAT(
' ----- DADMKS: LACK OF INITIALISATION')
2164 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2165c ----------------------------------------------------------------------
2166c it simulates kaon* decay in tau rest frame with
2167c z-axis along kaon* momentum
2168c jkst=10 for k* --->k0 + pi+-
2169c jkst=20 for k* --->k+- + pi0
2170c ----------------------------------------------------------------------
2171 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2172 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2173 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2174 * ,ampiz,ampi,amro,gamro,ama1,gama1
2175 * ,amk,amkz,amkst,gamkst
2177 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2178 * ,ampiz,ampi,amro,gamro,ama1,gama1
2179 * ,amk,amkz,amkst,gamkst
2180 REAL HV(4),PT(4),PN(4),PKS(4),PKK(4),PPI(4),QQ(4),RR1(1)
2182 DATA pi /3.141592653589793238462643/
2185c three body phase space normalised as in bjorken-drell
2186 phspac=1./2**11/pi**5
2193c here begin the k0,pi+_ decay
2196c mass of(real/virtual) k*
2198 ams2=(amtau-amnuta)**2
2200c amx2=ams1+ rr1(1)*(ams2-ams1)
2202c phspac=phspac*(ams2-ams1)
2203c phase space with sampling for k* resonance
2204 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2205 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2206 alp=alp1+rr1(1)*(alp2-alp1)
2207 amx2=amkst**2+amkst*gamkst*tan(alp)
2209 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2211 phspac=phspac*(alp2-alp1)
2213c tau-neutrino momentum
2216 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2217 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2222 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2224 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2227 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2228 pppi=sqrt((enpi-ampi)*(enpi+ampi))
2229 phspac=phspac*(4*pi)*(2*pppi/amx)
2230c charged pi momentum in kaon* rest frame
2231 CALL sphera(pppi,ppi)
2233c neutral kaon momentum in k* rest frame
2236 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2237 exe=(pks(4)+pks(3))/amx
2238c pion and k boosted from k* rest frame to tau rest frame
2239 CALL bostr3(exe,ppi,ppi)
2240 CALL bostr3(exe,pkk,pkk)
224230 qq(i)=ppi(i)-pkk(i)
2243c qq transverse to pks
2244 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2245 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
224731 qq(i)=qq(i)-pks(i)*qqpks/pksd
2250 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2252 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2253 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2254 & +(gv**2-ga**2)*amtau*amnuta*qq2
2255c a simple breit-wigner is chosen for k* resonance
2256 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2257 amplit=(gfermi*scabib)**2*brak*2*fks
2258 dgamt=1/(2.*amtau)*amplit*phspac
2260 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2262c here begin the k+-,pi0 decay
2263 ELSEIF(jkst.EQ.20)
THEN
2264c ======================
2265c mass of(real/virtual) k*
2267 ams2=(amtau-amnuta)**2
2269c amx2=ams1+ rr1*(ams2-ams1)
2271c phspac=phspac*(ams2-ams1)
2272c phase space with sampling for k* resonance
2273 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2274 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2275 alp=alp1+rr1(1)*(alp2-alp1)
2276 amx2=amkst**2+amkst*gamkst*tan(alp)
2278 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2280 phspac=phspac*(alp2-alp1)
2282c tau-neutrino momentum
2285 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2286 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2290 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2292 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2295 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2296 pppi=sqrt((enpi-ampiz)*(enpi+ampiz))
2297 phspac=phspac*(4*pi)*(2*pppi/amx)
2298c neutral pi momentum in k* rest frame
2299 CALL sphera(pppi,ppi)
2301c charged kaon momentum in k* rest frame
2304 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2305 exe=(pks(4)+pks(3))/amx
2306c pion and k boosted from k* rest frame to tau rest frame
2307 CALL bostr3(exe,ppi,ppi)
2308 CALL bostr3(exe,pkk,pkk)
231060 qq(i)=pkk(i)-ppi(i)
2311c qq transverse to pks
2312 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2313 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
231561 qq(i)=qq(i)-pks(i)*qqpks/pksd
2318 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2320 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2321 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2322 & +(gv**2-ga**2)*amtau*amnuta*qq2
2323c a simple breit-wigner is chosen for the k* resonance
2324 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2325 amplit=(gfermi*scabib)**2*brak*2*fks
2326 dgamt=1/(2.*amtau)*amplit*phspac
2328 70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2335 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2336c ----------------------------------------------------------------------
2337c it simulates multipi decay in tau rest frame with
2338c z-axis opposite to neutrino momentum
2339c ----------------------------------------------------------------------
2340 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2341 * ,ampiz,ampi,amro,gamro,ama1,gama1
2342 * ,amk,amkz,amkst,gamkst
2344 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2345 * ,ampiz,ampi,amro,gamro,ama1,gama1
2346 * ,amk,amkz,amkst,gamkst
2347 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2348 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2349 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2350 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2352 CHARACTER NAMES(NMODE)*31
2355 real*8 pn(4),pr(4),ppi(4,9),hv(4)
2356 real*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2357 real*8 pv(5,9),pt(4),ue(3),be(3)
2358 real*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2362 real*8 gam,bep,phi,a,b,c
2364 real*4 rrr(9),rrx(2),rn(1),rr2(1)
2366 DATA pi /3.141592653589793238462643/
2367 DATA wetmax /20*1d-15/
2369cc-- pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2.*a)
2372 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2374 ampik(i,j)=dcdmas(idffin(i,j))
2377 IF ((jnpi.LE.0).OR.jnpi.GT.20)
THEN
2378 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2392 phspac = 1./2.**5 /pi**2
23944 ps =ps+ampik(i,jnpi)
2397 ams2=(amtau-amnuta)**2
2400 amx2=ams1+ rr2(1)*(ams2-ams1)
2403 phspac=phspac * (ams2-ams1)
2405c tau-neutrino momentum
2408 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2409 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2413 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2415 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2417c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
2418c or f.gilman sh.rhie phys.rev.d31,1066(1985)
2422 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2423c here was an error. 20.10.91 (zw)
2424c brak=2*(gv**2+ga**2)*(2*pxq*pxn+amx2*qxn)
2425 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2426 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2427cam assume neutrino mass=0. and sum over final polarisation
2428c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
2429 amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2430 dgamt=1./(2.*amtau)*amplit*phspac
2432c isotropic w decay in w rest frame
2437 pv(5,nd)=ampik(nd,jnpi)
2438c compute max. phase space factor
2439 pmax=amw-ps+ampik(nd,jnpi)
2442 pmax=pmax+ampik(il,jnpi)
2443 pmin=pmin+ampik(il+1,jnpi)
2444 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2446c --- 2.02.94 zw 9 lines
2451 223 ams1=ams1+ampik(jl,jnpi)
2453 amx =(amx-ampik(il,jnpi))
2455 phsmax=phsmax * (ams2-ams1)
2460cam generate nd-2 effective masses
2462 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2464 CALL ranmar(rrr,nd-2)
2468 231 ams1=ams1+ampik(jl,jnpi)
2470 ams2=(amx-ampik(il,jnpi))**2
2472 amx2=ams1+ rr1*(ams2-ams1)
2475 phspac=phspac * (ams2-ams1)
2476c --- 2.02.94 zw 1 line
2477 phs=phs* (ams2-ams1)
2478 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2479 phs =phs *pa/pv(5,il)
2481 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2482 phs =phs *pa/pv(5,nd-1)
2484 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2485 IF (ncont.EQ.500 000)
THEN
2488 xnpi=xnpi+ampik(kk,jnpi)
2490 WRITE(6,*)
'ROUNDING INSTABILITY IN DPHNPI ?'
2491 WRITE(6,*)
'AMW=',amw,
'XNPI=',xnpi
2492 WRITE(6,*)
'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2493 WRITE(6,*)
'PHS=',phs,
'PHSMAX=',phsmax
2496 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs)
GO TO 100
2497c...perform successive two-particle decays in respective cm frame
2498 280
DO 300 il=1,nd-1
2499 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2503 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2504 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2507 290 pv(j,il+1)=-pa*ue(j)
2508 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2509 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2510 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2512c...lorentz transform decay products to tau frame
2514 310 ppi(j,nd)=pv(j,nd)
2517 320 be(j)=pv(j,il)/pv(4,il)
2518 gam=pv(4,il)/pv(5,il)
2520 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2522 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2523 ppi(4,i)=gam*(ppi(4,i)+bep)
2540 FUNCTION sigee(Q2,JNP)
2541c ----------------------------------------------------------------------
2542c e+e- cross section in the(1.gev2,amtau**2) region
2543c normalised to sig0 = 4/3 pi alfa2
2544c used in matrix element for multipion tau decays
2545c cf ys.tsai phys.rev d4 ,2821(1971)
2546c f.gilman et al phys.rev d17,1846(1978)
2547c c.kiesling, to be pub. in high energy e+e- physics(1988)
2548c datsig(*,1) = e+e- -> pi+pi-2pi0
2549c datsig(*,2) = e+e- -> 2pi+2pi-
2550c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2551c(phys lett 78b,623(1978)
2552c datsig(*,5) = e+e- -> 6pi
2554c 4- and 6-pion cross sections from
data
2555c 5-pion contribution related to 4-pion cross section
2558c ----------------------------------------------------------------------
2559 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2560 * ,ampiz,ampi,amro,gamro,ama1,gama1
2561 * ,amk,amkz,amkst,gamkst
2563 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2564 * ,ampiz,ampi,amro,gamro,ama1,gama1
2565 * ,amk,amkz,amkst,gamkst
2569 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2570 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2571 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2572 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2575 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2578 DATA pi /3.141592653589793238462643/
2586c ajwmod: initialize
if called from outside qq:
2587 IF (ampi.LT.0.139) ampi = 0.1395675
2591 datsig(i,2) = datsig(i,2)/2.
2592 datsig(i,1) = datsig(i,1) + datsig(i,2)
2598 IF(t . gt. s-ampi )
GO TO 201
2600 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2601 fact = fact * (datsig(j,1)+datsig(j+1,1))
2602 200 datsig(i,3) = datsig(i,3) + fact
2603 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2604 datsig(i,4) = datsig(i,3)
2605 datsig(i,6) = datsig(i,5)
2607c
WRITE(6,1000) datsig
2608 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2614 sigee=datsig(1,jnpi)+
2615 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2616 ELSEIF(q.LT.1.8)
THEN
2619 IF(q.LT.qmax)
GO TO 2
2622 2 sigee=datsig(i,jnpi)+
2623 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2624 ELSEIF(q.GT.1.8)
THEN
2625 sigee=datsig(17,jnpi)+
2626 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2628 IF(sigee.LT..0) sigee=0.
2630 sigee = sigee/(6.*pi**2*sig0)
2635 FUNCTION sigold(Q2,JNPI)
2636c ----------------------------------------------------------------------
2637c e+e- cross section in the(1.gev2,amtau**2) region
2638c normalised to sig0 = 4/3 pi alfa2
2639c used in matrix element for multipion tau decays
2640c cf ys.tsai phys.rev d4 ,2821(1971)
2641c f.gilman et al phys.rev d17,1846(1978)
2642c c.kiesling, to be pub. in high energy e+e- physics(1988)
2643c datsig(*,1) = e+e- -> pi+pi-2pi0
2644c datsig(*,2) = e+e- -> 2pi+2pi-
2645c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2646c(phys lett 78b,623(1978)
2647c datsig(*,4) = e+e- -> 6pi
2649c 4- and 6-pion cross sections from
data
2650c 5-pion contribution related to 4-pion cross section
2653c ----------------------------------------------------------------------
2654 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2655 * ,ampiz,ampi,amro,gamro,ama1,gama1
2656 * ,amk,amkz,amkst,gamkst
2658 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2659 * ,ampiz,ampi,amro,gamro,ama1,gama1
2660 * ,amk,amkz,amkst,gamkst
2664 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2665 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2666 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2667 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2669 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2671 DATA pi /3.141592653589793238462643/
2679 datsig(i,2) = datsig(i,2)/2.
2680 datsig(i,1) = datsig(i,1) + datsig(i,2)
2686 IF(t . gt. s-ampi )
GO TO 201
2688 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2689 fact = fact * (datsig(j,1)+datsig(j+1,1))
2690 200 datsig(i,3) = datsig(i,3) + fact
2691 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2693c
WRITE(6,1000) datsig
2694 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2700 sigee=datsig(1,jnpi)+
2701 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2702 ELSEIF(q.LT.1.8)
THEN
2705 IF(q.LT.qmax)
GO TO 2
2708 2 sigee=datsig(i,jnpi)+
2709 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2710 ELSEIF(q.GT.1.8)
THEN
2711 sigee=datsig(17,jnpi)+
2712 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2714 IF(sigee.LT..0) sigee=0.
2716 sigee = sigee/(6.*pi**2*sig0)
2721 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2722c ----------------------------------------------------------------------
2723* it simulates three pi(k) decay in the tau rest frame
2724* z-axis along hadronic system
2725c ----------------------------------------------------------------------
2726 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2727 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2729 CHARACTER NAMES(NMODE)*31
2731 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PNPI(4,9)
2732c matrix element number:
2734c
TYPE of the generation:
2737c --- masses of the decay products
2738 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
2739 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
2740 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
2742 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2754 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
2755c ----------------------------------------------------------------------
2756* it simulates a1 decay in tau rest frame with
2757* z-axis along a1 momentum
2758* it can be also used to generate k k pi and k pi pi tau decays.
2760* keyt - algorithm controlling switch
2761* 2 - flat phase space pim1 pim2 symmetrized statistical factor 1/2
2762* 1 - like 1 but peaked around a1 and rho(two channels) masses.
2763* 3 - peaked around omega, all particles different
2764* other- flat phase space, all particles different
2765* amp1 - mass of first pi, etc. (1-3)
2766* mnum - matrix element
type
2767* 0 - a1 matrix element
2768* 1-6 - matrix element for k pi pi, k k pi decay modes
2769* 7 - pi- pi0 gamma matrix element
2770c ----------------------------------------------------------------------
2771 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2772 * ,ampiz,ampi,amro,gamro,ama1,gama1
2773 * ,amk,amkz,amkst,gamkst
2775 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2776 * ,ampiz,ampi,amro,gamro,ama1,gama1
2777 * ,amk,amkz,amkst,gamkst
2778 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2779 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2780 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
2783 DATA pi /3.141592653589793238462643/
2785 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
2786c amro, gamro is only a
PARAMETER for geting hight efficiency
2788c three body phase space normalised as in bjorken-drell
2789c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
2790 phspac=1./2**17/pi**8
2800 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
2801 $ amrx,gamrx,amra,gamra,amrb,gamrb)
2802 IF (ichan.EQ.1)
THEN
2805 ELSEIF (ichan.EQ.2)
THEN
2814 ams1=(amp1+amp2+amp3)**2
2815 ams2=(amtau-amnuta)**2
2816* phase space with sampling for a1 resonance
2817 alp1=atan((ams1-amrx**2)/amrx/gamrx)
2818 alp2=atan((ams2-amrx**2)/amrx/gamrx)
2819 alp=alp1+rr1*(alp2-alp1)
2820 am3sq =amrx**2+amrx*gamrx*tan(alp)
2822 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
2823 phspac=phspac*(alp2-alp1)
2824c mass of(real/virtual) rho -
2828 IF (ichan.LE.2)
THEN
2829* phase space with sampling for rho resonance,
2830 alp1=atan((ams1-amra**2)/amra/gamra)
2831 alp2=atan((ams2-amra**2)/amra/gamra)
2832 alp=alp1+rr2*(alp2-alp1)
2833 am2sq =amra**2+amra*gamra*tan(alp)
2835c --- this part of the jacobian will be recovered later ---------------
2836c phspac=phspac*(alp2-alp1)
2837c phspac=phspac*((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
2838c----------------------------------------------------------------------
2841 am2sq=ams1+ rr2*(ams2-ams1)
2845* rho restframe, define pipl and pim1
2846 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
2847 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
2848 ppi= enq1**2-amp3**2
2849 pppi=sqrt(abs(enq1**2-amp3**2))
2850c --- this part of jacobian will be recovered later
2851 phf1=(4*pi)*(2*pppi/am2)
2852* pi minus momentum in rho rest frame
2853 CALL sphera(pppi,pipl)
2855* pi0 1 momentum in rho rest frame
2859* a1 rest frame, define pim2
2863 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
2864 pr(3)= sqrt(abs(pr(4)**2-am2**2))
2865 ppi = pr(4)**2-am2**2
2869 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
2871 phf2=(4*pi)*(2*pr(3)/am3)
2872* old pions boosted from rho rest frame to a1 rest frame
2873 exe=(pr(4)+pr(3))/am2
2874 CALL bostr3(exe,pipl,pipl)
2875 CALL bostr3(exe,pim1,pim1)
2879 thet =acos(-1.+2*rr3)
2881 CALL rotpol(thet,phi,pipl)
2882 CALL rotpol(thet,phi,pim1)
2883 CALL rotpol(thet,phi,pim2)
2884 CALL rotpol(thet,phi,pr)
2886* now to the tau rest frame, define a1 and neutrino momenta
2890 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
2891 paa(3)= sqrt(abs(paa(4)**2-am3**2))
2892 ppi = paa(4)**2-am3**2
2893 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
2894* tau-neutrino momentum
2897 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
2899c here we correct for the jacobians of the two chains
2900c ---first channel ------- pim1+pipl
2903 alp1=atan((ams1-amra**2)/amra/gamra)
2904 alp2=atan((ams2-amra**2)/amra/gamra)
2905 xpro = (pim1(3)+pipl(3))**2
2906 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
2907 am2sq=-xpro+(pim1(4)+pipl(4))**2
2908c jacobian of speeding
2909 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
2910 ff1 =ff1 *(alp2-alp1)
2911c lambda of rho decay
2912 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
2914 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
2915 xjaje=gg1*(ams2-ams1)
2916c ---second channel ------ pim2+pipl
2919 alp1=atan((ams1-amrb**2)/amrb/gamrb)
2920 alp2=atan((ams2-amrb**2)/amrb/gamrb)
2921 xpro = (pim2(3)+pipl(3))**2
2922 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
2923 am2sq=-xpro+(pim2(4)+pipl(4))**2
2924 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
2925 ff2 =ff2 *(alp2-alp1)
2926 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
2927 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
2928 xjadw=gg2*(ams2-ams1)
2935 IF (ichan.EQ.2)
THEN
2940 IF (xjac1.NE.0.0) a1=prob1/xjac1
2941 IF (xjac2.NE.0.0) a2=prob2/xjac2
2942 IF (xjac3.NE.0.0) a3=prob3/xjac3
2944 IF (a1+a2+a3.NE.0.0)
THEN
2945 phspac=phspac/(a1+a2+a3)
2955* all pions boosted from a1 rest frame to tau rest frame
2956* z-axis antiparallel to neutrino momentum
2957 exe=(paa(4)+paa(3))/am3
2958 CALL bostr3(exe,pipl,pipl)
2959 CALL bostr3(exe,pim1,pim1)
2960 CALL bostr3(exe,pim2,pim2)
2961 CALL bostr3(exe,pr,pr)
2962c partial width consists of phase space and amplitude
2964 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
2965c
ELSEIF (mnum.EQ.0)
THEN
2966c
CALL dampaa(pt,pn,pim1,pim2,pipl,amplit,hv)
2968 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
2970 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN
2971c the statistical factor for identical pi-s is cancelled with
2972c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
2976 dgamt=1/(2.*amtau)*amplit*phspac
2978 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
2979c ----------------------------------------------------------------------
2980* calculates differential cross section and polarimeter vector
2981* for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
2982* all spin effects in the full decay chain are taken into account.
2983* calculations done in tau rest frame with z-axis along neutrino moment
2984* the routine is writen for zero neutrino mass.
2987c ----------------------------------------------------------------------
2988 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2989 * ,ampiz,ampi,amro,gamro,ama1,gama1
2990 * ,amk,amkz,amkst,gamkst
2992 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2993 * ,ampiz,ampi,amro,gamro,ama1,gama1
2994 * ,amk,amkz,amkst,gamkst
2995 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2996 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2997 COMMON /testa1/ keya1
2998 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
2999 REAL PAA(4),VEC1(4),VEC2(4)
3000 REAL PIVEC(4),PIAKS(4),HVM(4)
3001 COMPLEX BWIGN,HADCUR(4),FPIK
3004* f constants for a1, a1-rho-pi, and rho-pi-pi
3007* this inline funct. calculates the scalar part of the propagator
3008 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3010* four momentum of a1
3012 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3013* masses of a1, and of two pi-pairs which may form rho
3014 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3015 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3016 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3017 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3018 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3019* elements of hadron current
3020 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3021 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3022 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3023 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3025 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3026 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3027* hadron current saturated with a1 and rho resonances
3028 IF (keya1.EQ.1)
THEN
3032 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3034 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3035 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3036 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3039 fnorm=2.0*sqrt(2.)/3.0/fpi
3040 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3042 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3043 $ *(cmplx(vec1(i))*fpik(xmro1)
3044 $ +cmplx(vec2(i))*fpik(xmro2))
3048* calculate pi-vectors: vector and axial
3049 CALL clvec(hadcur,pn,pivec)
3050 CALL claxi(hadcur,pn,piaks)
3051 CALL clnut(hadcur,brakm,hvm)
3052* spin independent part of decay diff-cross-sect. in tau rest frame
3053 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3054 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3055 amplit=(gfermi*ccabib)**2*brak/2.
3056c the statistical factor for identical pi-s was cancelled with
3057c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3058c polarimeter vector in tau rest frame
3060 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3061 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3062c hv is defined for tau- with gamma=b+hv*pol
3068c ****************************************************************
3069c g-
FUNCTION used to inroduce energy dependence in a1 width
3070c ****************************************************************
3071 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3072 * ,ampiz,ampi,amro,gamro,ama1,gama1
3073 * ,amk,amkz,amkst,gamkst
3075 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3076 * ,ampiz,ampi,amro,gamro,ama1,gama1
3077 * ,amk,amkz,amkst,gamkst
3079 IF (qkwa.LT.(amro+ampi)**2)
THEN
3080 gfun=4.1*(qkwa-9*ampiz**2)**3
3081 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3083 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3086 COMPLEX FUNCTION bwigs(S,M,G)
3087c **********************************************************
3088c p-wave breit-wigner for k*
3089c **********************************************************
3091 REAL PI,PIM,QS,QM,W,GS,MK
3092c ajw: add k*-prim possibility:
3096 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3097 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3098c ------------ parameters --------------------
3104c ajw: add k*-prim possibility:
3108c ------- breit-wigner -----------------------
3110 qs=p(s,pim**2,mk**2)
3111 qm=p(m**2,pim**2,mk**2)
3113 gs=g*(m/w)*(qs/qm)**3
3114 bw=m**2/cmplx(m**2-s,-m*gs)
3115 qpm=p(pm**2,pim**2,mk**2)
3116 g1=pg*(pm/w)*(qs/qpm)**3
3117 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3118 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3121 COMPLEX FUNCTION bwig(S,M,G)
3122c **********************************************************
3123c p-wave breit-wigner for rho
3124c **********************************************************
3126 REAL PI,PIM,QS,QM,W,GS
3128c ------------ parameters --------------------
3133c ------- breit-wigner -----------------------
3135 IF (s.GT.4.*pim**2)
THEN
3136 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3137 qm=sqrt(m**2/4.-pim**2)
3139 gs=g*(m/w)*(qs/qm)**3
3143 bwig=m**2/cmplx(m**2-s,-m*gs)
3146 COMPLEX FUNCTION fpik(W)
3147c **********************************************************
3149c **********************************************************
3151 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
3155c ------------ parameters --------------------
3156 IF (init.EQ.0 )
THEN
3166c -----------------------------------------------
3168 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3173c **********************************************************
3174c square of pion form factor
3175c **********************************************************
3177 fpirho=cabs(fpik(w))**2
3179 SUBROUTINE clvec(HJ,PN,PIV)
3180c ----------------------------------------------------------------------
3181* calculates the
"VECTOR TYPE" pi-vector piv
3182* note that the neutrino mom. pn is assumed to be along z-axis
3185c ----------------------------------------------------------------------
3189 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3190 hh= real(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3191 $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3193 10 piv(i)=4.*real(hn*conjg(hj(i)))-2.*hh*pn(i)
3196 SUBROUTINE claxi(HJ,PN,PIA)
3197c ----------------------------------------------------------------------
3198* calculates the
"AXIAL TYPE" pi-vector pia
3199* note that the neutrino mom. pn is assumed to be along z-axis
3200c sign is chosen +/- for decay of tau +/- respectively
3201c called by : dampaa, clnut
3202c ----------------------------------------------------------------------
3203 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3204 COMMON / idfc / idff
3206 COMPLEX HJ(4),HJC(4)
3207c det2(i,j)=aimag(hj(i)*hjc(j)-hj(j)*hjc(i))
3208c -- here was an error(zw, 21.11.1991)
3209 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3210c -- it was affecting sign of a_lr asymmetry in a1 decay.
3211c -- note also collision of notation of gamma_va as defined in
3212c -- tauola paper and j.h. kuhn and santamaria z. phys c 48 (1990) 445
3213* -----------------------------------
3214 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3215 sign= idff/abs(idff)
3216 ELSEIF (ktom.EQ.2)
THEN
3217 sign=-idff/abs(idff)
3219 print *,
'STOP IN CLAXI: KTOM=',ktom
3224 10 hjc(i)=conjg(hj(i))
3225 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3226 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3227 pia(3)= 2.*pn(4)*det2(1,2)
3228 pia(4)= 2.*pn(3)*det2(1,2)
3229c all four indices are up so pia(3) and pia(4) have same sign
3231 20 pia(i)=pia(i)*sign
3233 SUBROUTINE clnut(HJ,B,HV)
3234c ----------------------------------------------------------------------
3235* calculates the contribution by neutrino mass
3236* note the tau is assumed to be at rest
3239c ----------------------------------------------------------------------
3245 b=real( hj(4)*aimag(hj(4)) - hj(3)*aimag(hj(3))
3246 & - hj(2)*aimag(hj(2)) - hj(1)*aimag(hj(1)) )
3249 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3250c ----------------------------------------------------------------------
3251* calculates differential cross section and polarimeter vector
3252* for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
3253* all spin effects in the full decay chain are taken into account.
3254* calculations done in tau rest frame with z-axis along neutrino moment
3255* the routine is writen for zero neutrino mass.
3258c ----------------------------------------------------------------------
3259 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3260 * ,ampiz,ampi,amro,gamro,ama1,gama1
3261 * ,amk,amkz,amkst,gamkst
3263 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3264 * ,ampiz,ampi,amro,gamro,ama1,gama1
3265 * ,amk,amkz,amkst,gamkst
3266 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3267 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3268 COMMON /testa1/ keya1
3269 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3270 REAL PAA(4),VEC1(4),VEC2(4)
3271 REAL PIVEC(4),PIAKS(4),HVM(4)
3272 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3274* this inline funct. calculates the scalar part of the propagator
3275c ajwmod to satisfy compiler, comment out this unused function.
3277* four momentum of a1
3282 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3284* masses of a1, and of two pi-pairs which may form rho
3285 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3286 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3287 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3288 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3289* elements of hadron current
3290 prod1 =vec1(1)*pipl(1)
3291 prod2 =vec2(2)*pipl(2)
3292 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3293 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3294 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3295 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3296 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3297 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3299 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3301 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3303 vec1(i)= vec1(i)/gnorm
3305 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3306 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3307 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3308 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3309 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3310 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3311 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3312 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3313 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3314 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3315 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3317 fnorm=formom(xmaa,xmom)
3322 hadcur(i) = fnorm *(
3323 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3324 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3325 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3327 hadcur(i) = fnorm *(
3328 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3329 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3330 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3334* calculate pi-vectors: vector and axial
3335 CALL clvec(hadcur,pn,pivec)
3336 CALL claxi(hadcur,pn,piaks)
3337 CALL clnut(hadcur,brakm,hvm)
3338* spin independent part of decay diff-cross-sect. in tau rest frame
3339 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3340 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3342 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3343 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3345c hv is defined for tau- with gamma=b+hv*pol
3347 amplit=(gfermi*ccabib)**2*brak/2.
3348c the statistical factor for identical pi-s was cancelled with
3349c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3350c polarimeter vector in tau rest frame
3356 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3357c ----------------------------------------------------------------------
3358* calculates differential cross section and polarimeter vector
3359* for tau decay into k k pi, k pi pi.
3360* all spin effects in the full decay chain are taken into account.
3361* calculations done in tau rest frame with z-axis along neutrino moment
3362c mnum decay mode identifier.
3365c ----------------------------------------------------------------------
3366 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3367 * ,ampiz,ampi,amro,gamro,ama1,gama1
3368 * ,amk,amkz,amkst,gamkst
3370 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3371 * ,ampiz,ampi,amro,gamro,ama1,gama1
3372 * ,amk,amkz,amkst,gamkst
3373 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3374 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3375 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3376 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3377 REAL PIVEC(4),PIAKS(4),HVM(4)
3378 REAL FNORM(0:7),COEF(1:5,0:7)
3379 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3380 COMPLEX F1,F2,F3,F4,F5
3381 EXTERNAL form1,form2,form3,form4,form5
3382 DATA pi /3.141592653589793238462643/
3386 IF (icont.EQ.0)
THEN
3394 fnorm(4)=scabib/fpi/dwapi0
3399 coef(1,0)= 2.0*sqrt(2.)/3.0
3400 coef(2,0)=-2.0*sqrt(2.)/3.0
3401c ajw 2/98: add in the d-wave and i=0 3pi substructure:
3402 coef(3,0)= 2.0*sqrt(2.)/3.0
3406 coef(1,1)=-sqrt(2.)/3.0
3407 coef(2,1)= sqrt(2.)/3.0
3412 coef(1,2)=-sqrt(2.)/3.0
3413 coef(2,2)= sqrt(2.)/3.0
3418c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3425 coef(1,4)= 1.0/sqrt(2.)/3.0
3426 coef(2,4)=-1.0/sqrt(2.)/3.0
3431 coef(1,5)=-sqrt(2.)/3.0
3432 coef(2,5)= sqrt(2.)/3.0
3437c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3448 coef(5,7)=-sqrt(2.0/3.0)
3453 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3454 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3455 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3456 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3457 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3458 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3459 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3460 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3461* elements of hadron current
3462 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3463 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3464 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3465 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3466 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3467 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3469 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3470 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3471 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3472 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3473 CALL prod5(pim1,pim2,pim3,vec5)
3475c be aware that sign of vec2 is opposite to sign of vec1 in a1
case
3476c rationalize this code:
3477 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3478 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3479 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3481 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3482 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3483 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3486 hadcur(i)= cmplx(fnorm(mnum)) * (
3487 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3488 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3491* calculate pi-vectors: vector and axial
3492 CALL clvec(hadcur,pn,pivec)
3493 CALL claxi(hadcur,pn,piaks)
3494 CALL clnut(hadcur,brakm,hvm)
3495* spin independent part of decay diff-cross-sect. in tau rest frame
3496 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3497 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3498 amplit=(gfermi)**2*brak/2.
3500 print *,
'MNUM=',mnum
3506 IF (k.EQ.4) znak=1.0
3507 xm1=znak*pim1(k)**2+xm1
3508 xm2=znak*pim2(k)**2+xm2
3509 xm3=znak*pim3(k)**2+xm3
3510 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3511 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3512 print *,
'************************************************'
3514c polarimeter vector in tau rest frame
3516 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3517 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3518c hv is defined for tau- with gamma=b+hv*pol
3522 SUBROUTINE prod5(P1,P2,P3,PIA)
3523c ----------------------------------------------------------------------
3524c
external product of P1, P2, P3 4-momenta.
3525c sign is chosen +/- for decay of tau +/- respectively
3526c called by : dampaa, clnut
3527c ----------------------------------------------------------------------
3528 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3529 COMMON / idfc / idff
3530 REAL PIA(4),P1(4),P2(4),P3(4)
3531 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3532* -----------------------------------
3533 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3534 sign= idff/abs(idff)
3535 ELSEIF (ktom.EQ.2)
THEN
3536 sign=-idff/abs(idff)
3538 print *,
'STOP IN PROD5: KTOM=',ktom
3542c epsilon( p1(1), p2(2), p3(3), (4) ) = 1
3544 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3545 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3546 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3547 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3548c all four indices are up so pia(3) and pia(4) have same sign
3550 20 pia(i)=pia(i)*sign
3553 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3554c ----------------------------------------------------------------------
3555* this simulates tau decay in tau rest frame
3556* into nu a1, next a1 decays into rho pi and finally rho into pi pi.
3557* output four momenta: pnu tauneutrino,
3559* pim1 pion minus(or pi0) 1 (for tau minus)
3560* pim2 pion minus(or pi0) 2
3561* pipl pion plus(or pi-)
3562* (pipl,pim1) form a rho
3563c ----------------------------------------------------------------------
3564 COMMON / inout / inut,iout
3565 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3569c ===================
3571 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3572cc
CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3574 ELSEIF(mode.EQ. 0)
THEN
3575* =======================
3577 IF(iwarm.EQ.0)
GOTO 902
3578 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3579 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3580cc
CALL hfill(816,wt)
3582 IF(rn(1).GT.wt)
GOTO 300
3584 ELSEIF(mode.EQ. 1)
THEN
3585* =======================
3586 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3591 902
WRITE(iout, 9020)
3592 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3595 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3596c ----------------------------------------------------------------------
3597 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3598 * ,ampiz,ampi,amro,gamro,ama1,gama1
3599 * ,amk,amkz,amkst,gamkst
3601 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3602 * ,ampiz,ampi,amro,gamro,ama1,gama1
3603 * ,amk,amkz,amkst,gamkst
3604 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3605 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3606 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3607 real*4 gampmc ,gamper
3608 COMMON / inout / inut,iout
3609 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3610 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3612 CHARACTER NAMES(NMODE)*31
3614 real*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3615 real*4 pdum1(4),pdum2(4),pdumi(4,9)
3618 real*8 swt(nmode),sswt(nmode)
3619 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3621 DATA pi /3.141592653589793238462643/
3625c ===================
3626c -- at the moment only two decay modes of multipions have m. elem
3637c for 4pi phase space, need lots more trials at initialization,
3638c or
use the wtmax determined with many trials for default model:
3640 IF (jnpi.LE.nm4)
THEN
3641 a = pkorb(3,37+jnpi)
3652 ELSEIF(jnpi.LE.nm4)
THEN
3653 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3654 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3655 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3656 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3657 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3658 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3659 inum=jnpi-nm4-nm5-nm6
3660 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
3661 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3662 inum=jnpi-nm4-nm5-nm6-nm3
3663 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
3667 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
3669c print *,
' DADNEW JNPI,NTRIALS,WTMAX =',jnpi,ntrials,wtmax(jnpi)
3670c
CALL hbook1(801,
'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
3671c print 7004,wtmax(jnpi)
3675 ELSEIF(mode.EQ. 0)
THEN
3676c =======================
3677 IF(iwarm.EQ.0)
GOTO 902
3682 ELSEIF(jnpi.LE.nm4)
THEN
3683 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3684 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3685 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3686 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3687 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
3688 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3689 inum=jnpi-nm4-nm5-nm6
3690 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
3691 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3692 inum=jnpi-nm4-nm5-nm6-nm3
3693 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
3700c
CALL hfill(801,wt/wtmax(jnpi))
3701 nevraw(jnpi)=nevraw(jnpi)+1
3702 swt(jnpi)=swt(jnpi)+wt
3704cc sswt(jnpi)=sswt(jnpi)+wt**2
3705 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
3709 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
3710 IF(rn*wtmax(jnpi).GT.wt)
GOTO 300
3711c rotations to basic tau rest frame
3712 costhe=-1.+2.*rrr(2)
3715 CALL rotor2(thet,pnu,pnu)
3716 CALL rotor3( phi,pnu,pnu)
3717 CALL rotor2(thet,pwb,pwb)
3718 CALL rotor3( phi,pwb,pwb)
3719 CALL rotor2(thet,hv,hv)
3720 CALL rotor3( phi,hv,hv)
3723 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
3724 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
3726 nevacc(jnpi)=nevacc(jnpi)+1
3728 ELSEIF(mode.EQ. 1)
THEN
3729c =======================
3731 IF(nevraw(jnpi).EQ.0)
GOTO 500
3732 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
3734 IF(nevraw(jnpi).NE.0)
3735 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
3737 WRITE(iout, 7010) names(jnpi),
3738 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
3740 gampmc(8+jnpi-1)=rat
3741 gamper(8+jnpi-1)=error
3742cam nevdec(8+jnpi-1)=nevacc(jnpi)
3747 7003
FORMAT(///1x,15(5h*****)
3748 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
3750 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
3752 $ /,1x,15(5h*****)/)
3753 7010
FORMAT(///1x,15(5h*****)
3754 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
3755 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
3756 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
3757 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
3758 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
3759 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
3760 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
3761 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
3762 $ /,1x,15(5h*****)/)
3763 902
WRITE(iout, 9020)
3764 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
3766 903
WRITE(iout, 9030) jnpi,mode
3767 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
3772 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
3773c ----------------------------------------------------------------------
3774* it simulates a1 decay in tau rest frame with
3775* z-axis along a1 momentum
3776c ----------------------------------------------------------------------
3777 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3778 * ,ampiz,ampi,amro,gamro,ama1,gama1
3779 * ,amk,amkz,amkst,gamkst
3781 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3782 * ,ampiz,ampi,amro,gamro,ama1,gama1
3783 * ,amk,amkz,amkst,gamkst
3784 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3785 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3786 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
3789 real*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
3790 DATA pi /3.141592653589793238462643/
3792 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3793c amro, gamro is only a
PARAMETER for geting hight efficiency
3795c three body phase space normalised as in bjorken-drell
3796c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
3797 phspac=1./2**23/pi**11
3807c ajw: cant simply change amrop, etc, here
3808c choice is a by-hand tuning/optimization, no simple relationship
3809c to actual resonance masses(accd to z.was).
3810c what matters in the
end is what you put in formf/curr .
3826 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
3827 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
3837* masses of 4, 3 and 2 pi systems
3838c 3 pi with sampling for resonance
3841 ams1=(amp1+amp2+amp3+amp4)**2
3842 ams2=(amtau-amnuta)**2
3843 alp1=atan((ams1-amrop**2)/amrop/gamrop)
3844 alp2=atan((ams2-amrop**2)/amrop/gamrop)
3845 alp=alp1+rr1*(alp2-alp1)
3846 am4sq =amrop**2+amrop*gamrop*tan(alp)
3849 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
3850 phspac=phspac*(alp2-alp1)
3854 ams1=(amp2+amp3+amp4)**2
3856 IF (rrr(9).GT.prez)
THEN
3857 am3sq=ams1+ rr1*(ams2-ams1)
3859c --- this part of jacobian will be recovered later
3862* phase space with sampling for omega resonance,
3863 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3864 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3865 alp=alp1+rr1*(alp2-alp1)
3866 am3sq =amrx**2+amrx*gamrx*tan(alp)
3868c --- this part of the jacobian will be recovered later ---------------
3869 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3877 am2sq=ams1+ rr2*(ams2-ams1)
3879c --- this part of jacobian will be recovered later
3881* 2 restframe, define piz and pipl
3882 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
3883 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
3884 ppi= enq1**2-amp4**2
3885 pppi=sqrt(abs(enq1**2-amp4**2))
3886 phspac=phspac*(4*pi)*(2*pppi/am2)
3887* piz momentum in 2 rest frame
3888 CALL sphera(pppi,piz)
3890* pipl momentum in 2 rest frame
3894* 3 rest frame, define pim1
3898 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
3899 pr(3)= sqrt(abs(pr(4)**2-am2**2))
3900 ppi = pr(4)**2-am2**2
3904 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
3906c --- this part of jacobian will be recovered later
3907 ff3=(4*pi)*(2*pr(3)/am3)
3908* old pions boosted from 2 rest frame to 3 rest frame
3909 exe=(pr(4)+pr(3))/am2
3910 CALL bostr3(exe,piz,piz)
3911 CALL bostr3(exe,pipl,pipl)
3914 thet =acos(-1.+2*rr3)
3916 CALL rotpol(thet,phi,pipl)
3917 CALL rotpol(thet,phi,pim1)
3918 CALL rotpol(thet,phi,piz)
3919 CALL rotpol(thet,phi,pr)
3920* 4 rest frame, define pim2
3924 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
3925 pr(3)= sqrt(abs(pr(4)**2-am3**2))
3926 ppi = pr(4)**2-am3**2
3930 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
3932c --- this part of jacobian will be recovered later
3933 ff4=(4*pi)*(2*pr(3)/am4)
3934* old pions boosted from 3 rest frame to 4 rest frame
3935 exe=(pr(4)+pr(3))/am3
3936 CALL bostr3(exe,piz,piz)
3937 CALL bostr3(exe,pipl,pipl)
3938 CALL bostr3(exe,pim1,pim1)
3941 thet =acos(-1.+2*rr3)
3943 CALL rotpol(thet,phi,pipl)
3944 CALL rotpol(thet,phi,pim1)
3945 CALL rotpol(thet,phi,pim2)
3946 CALL rotpol(thet,phi,piz)
3947 CALL rotpol(thet,phi,pr)
3949* now to the tau rest frame, define paa and neutrino momenta
3953 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
3954 paa(3)= sqrt(abs(paa(4)**2-am4**2))
3955 ppi = paa(4)**2-am4**2
3956 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3957 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
3958* tau-neutrino momentum
3961 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
3963c zbw 20.12.2002 bug fix
3964 IF(rrr(9).LE.0.5*prez)
THEN
3971c we include remaining part of the jacobian
3973 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
3974 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
3976 ams1=(amp1+amp3+amp4)**2
3979 ams2=(sqrt(am3sq)-amp1)**2
3981 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
3982 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
3985 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
3986 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
3988 ams1=(amp1+amp3+amp4)**2
3989 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3990 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3991 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3994 ams2=(sqrt(am3sq)-amp1)**2
3996 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
3997 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4000 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4001 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4003 ams1=(amp2+amp3+amp4)**2
4004 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4005 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4006 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4009 ams2=(sqrt(am3sq)-amp2)**2
4011 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4012 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4014c --- jacobian averaged over the two
4015 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4016 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4021* momenta of the two pi-minus are randomly symmetrised
4032c momenta of pi0-s are generated uniformly only
IF prez=0.0
4042* all pions boosted from 4 rest frame to tau rest frame
4043* z-axis antiparallel to neutrino momentum
4044 exe=(paa(4)+paa(3))/am4
4045 CALL bostr3(exe,piz,piz)
4046 CALL bostr3(exe,pipl,pipl)
4047 CALL bostr3(exe,pim1,pim1)
4048 CALL bostr3(exe,pim2,pim2)
4049 CALL bostr3(exe,pr,pr)
4050c partial width consists of phase space and amplitude
4051c check on consistency with dadnpi,
THEN, code breakes uniform pion
4052c distribution in hadronic system
4053cam assume neutrino mass=0. and sum over final polarisation
4055c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
4056c amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,1)
4058 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4059 ELSEIF (jnpi.EQ.2)
THEN
4060 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4062 dgamt=1/(2.*amtau)*amplit*phspac
4072 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4073c ----------------------------------------------------------------------
4074* calculates differential cross section and polarimeter vector
4075* for tau decay into 4 pi modes
4076* all spin effects in the full decay chain are taken into account.
4077* calculations done in tau rest frame with z-axis along neutrino moment
4078c mnum decay mode identifier.
4081c ----------------------------------------------------------------------
4082 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4083 * ,ampiz,ampi,amro,gamro,ama1,gama1
4084 * ,amk,amkz,amkst,gamkst
4086 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4087 * ,ampiz,ampi,amro,gamro,ama1,gama1
4088 * ,amk,amkz,amkst,gamkst
4089 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4090 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4091 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4092 REAL PIVEC(4),PIAKS(4),HVM(4)
4093 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4094 EXTERNAL form1,form2,form3,form4,form5
4095 DATA pi /3.141592653589793238462643/
4098 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4100* calculate pi-vectors: vector and axial
4101 CALL clvec(hadcur,pn,pivec)
4102 CALL claxi(hadcur,pn,piaks)
4103 CALL clnut(hadcur,brakm,hvm)
4104* spin independent part of decay diff-cross-sect. in tau rest frame
4105 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4106 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4107 amplit=(ccabib*gfermi)**2*brak/2.
4108c polarimeter vector in tau rest frame
4110 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4111 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4112c hv is defined for tau- with gamma=b+hv*pol
4117 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4118c ----------------------------------------------------------------------
4119* it simulates 5pi decay in tau rest frame with
4120* z-axis along 5pi momentum
4121c ----------------------------------------------------------------------
4122 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4123 * ,ampiz,ampi,amro,gamro,ama1,gama1
4124 * ,amk,amkz,amkst,gamkst
4126 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4127 * ,ampiz,ampi,amro,gamro,ama1,gama1
4130 * ,amk,amkz,amkst,gamkst
4131 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4132 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4133 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4134 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4136 CHARACTER NAMES(NMODE)*31
4137 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4138 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4139 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4140 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4142 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4147 DATA pi /3.141592653589793238462643/
4153 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4159c 6 body phase space normalised as in bjorken-drell
4160c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4161 phspac=1./2**29/pi**14
4162c phspac=1./2**5/pi**2
4163c init 5pi decay mode(jnpi)
4164 amp1=dcdmas(idffin(1,jnpi))
4165 amp2=dcdmas(idffin(2,jnpi))
4166 amp3=dcdmas(idffin(3,jnpi))
4167 amp4=dcdmas(idffin(4,jnpi))
4168 amp5=dcdmas(idffin(5,jnpi))
4178c masses of 5, 4, 3 and 2 pi systems
4179c 3 pi with sampling for omega resonance
4183 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4184 ams2=(amtau-amnuta)**2
4185 am5sq=ams1+ rr1*(ams2-ams1)
4187 phspac=phspac*(ams2-ams1)
4192 ams1=(amp2+amp3+amp4+amp5)**2
4194 am4sq=ams1+ rr1*(ams2-ams1)
4199c phase space with sampling for omega resonance
4201 ams1=(amp2+amp3+amp4)**2
4203 alp1=atan((ams1-amom**2)/amom/gamom)
4204 alp2=atan((ams2-amom**2)/amom/gamom)
4205 alp=alp1+rr1*(alp2-alp1)
4206 am3sq =amom**2+amom*gamom*tan(alp)
4208c --- this part of the jacobian will be recovered later ---------------
4209 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4212c am3sq=ams1+ rr1*(ams2-ams1)
4214c --- this part of jacobian will be recovered later
4222 am2sq=ams1+ rr2*(ams2-ams1)
4224c --- this part of jacobian will be recovered later
4227c(34) restframe, define pi3 and pi4
4228 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4229 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4230 ppi= enq1**2-amp3**2
4231 pppi=sqrt(abs(enq1**2-amp3**2))
4232 ff1=(4*pi)*(2*pppi/am2)
4233c pi3 momentum in(34) rest frame
4234 call sphera(pppi,pi3)
4236c pi4 momentum in(34) rest frame
4241c(234) rest frame, define pi2
4245 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4246 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4247 ppi = pr(4)**2-am2**2
4251 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4253c --- this part of jacobian will be recovered later
4254 ff2=(4*pi)*(2*pr(3)/am3)
4255c old pions boosted from 2 rest frame to 3 rest frame
4256 exe=(pr(4)+pr(3))/am2
4257 call bostr3(exe,pi3,pi3)
4258 call bostr3(exe,pi4,pi4)
4261 thet =acos(-1.+2*rr3)
4263 call rotpol(thet,phi,pi2)
4264 call rotpol(thet,phi,pi3)
4265 call rotpol(thet,phi,pi4)
4267c(2345) rest frame, define pi5
4271 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4272 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4273 ppi = pr(4)**2-am3**2
4277 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4279c --- this part of jacobian will be recovered later
4280 ff3=(4*pi)*(2*pr(3)/am4)
4281c old pions boosted from 3 rest frame to 4 rest frame
4282 exe=(pr(4)+pr(3))/am3
4283 call bostr3(exe,pi2,pi2)
4284 call bostr3(exe,pi3,pi3)
4285 call bostr3(exe,pi4,pi4)
4288 thet =acos(-1.+2*rr3)
4290 call rotpol(thet,phi,pi2)
4291 call rotpol(thet,phi,pi3)
4292 call rotpol(thet,phi,pi4)
4293 call rotpol(thet,phi,pi5)
4295c(12345) rest frame, define pi1
4299 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4300 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4301 ppi = pr(4)**2-am4**2
4305 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4307c --- this part of jacobian will be recovered later
4308 ff4=(4*pi)*(2*pr(3)/am5)
4309c old pions boosted from 4 rest frame to 5 rest frame
4310 exe=(pr(4)+pr(3))/am4
4311 call bostr3(exe,pi2,pi2)
4312 call bostr3(exe,pi3,pi3)
4313 call bostr3(exe,pi4,pi4)
4314 call bostr3(exe,pi5,pi5)
4317 thet =acos(-1.+2*rr3)
4319 call rotpol(thet,phi,pi1)
4320 call rotpol(thet,phi,pi2)
4321 call rotpol(thet,phi,pi3)
4322 call rotpol(thet,phi,pi4)
4323 call rotpol(thet,phi,pi5)
4325* now to the tau rest frame, define paa and neutrino momenta
4329c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4330c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4331c ppi = paa(4)**2-am5**2
4332 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4333 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4334 ppi = paa(4)**2-am5sq
4335 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4336* tau-neutrino momentum
4339 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4342 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4344c all pions boosted from 5 rest frame to tau rest frame
4345c z-axis antiparallel to neutrino momentum
4346 exe=(paa(4)+paa(3))/am5
4347 call bostr3(exe,pi1,pi1)
4348 call bostr3(exe,pi2,pi2)
4349 call bostr3(exe,pi3,pi3)
4350 call bostr3(exe,pi4,pi4)
4351 call bostr3(exe,pi5,pi5)
4353c partial width consists of phase space and amplitude
4354c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
4355c or f.gilman sh.rhie phys.rev.d31,1066(1985)
4359 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4360 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4361 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4362 fompp = cabs(bwign(am3,amom,gamom))**2
4363c normalisation factor(to some numerical undimensioned factor;
4364c cf r.fischer et al zphys c3, 313 (1980))
4366c amplit=ccabib**2*gfermi**2/2. * brak * am5sq*sigee(am5sq,jnpi)
4367 amplit=ccabib**2*gfermi**2/2. * brak
4368 amplit = amplit * fompp * fnorm
4370c amplit = amplit * fnorm
4371 dgamt=1/(2.*amtau)*amplit*phspac
4384c missing: transposition of identical particles, startistical factors
4385c for identical matrices, polarimetric vector. matrix element rather naive.
4386c flat phase space in pion system + with breit wigner for omega
4387c anyway it is better than nothing, and code is improvable.
4389 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4390c ----------------------------------------------------------------------
4391c it simulates rho decay in tau rest frame with
4392c z-axis along rho momentum
4393c rho decays to k kbar
4394c ----------------------------------------------------------------------
4395 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4396 * ,ampiz,ampi,amro,gamro,ama1,gama1
4397 * ,amk,amkz,amkst,gamkst
4399 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4400 * ,ampiz,ampi,amro,gamro,ama1,gama1
4401 * ,amk,amkz,amkst,gamkst
4402 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4403 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4404 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4406 DATA pi /3.141592653589793238462643/
4409c three body phase space normalised as in bjorken-drell
4410 phspac=1./2**11/pi**5
4416c mass of(real/virtual) rho
4418 ams2=(amtau-amnuta)**2
4421 amx2=ams1+ rr1(1)*(ams2-ams1)
4423 phspac=phspac*(ams2-ams1)
4424c phase space with sampling for rho resonance
4425c alp1=atan((ams1-amro**2)/amro/gamro)
4426c alp2=atan((ams2-amro**2)/amro/gamro)
4430c alp=alp1+rr1(1)*(alp2-alp1)
4431c amx2=amro**2+amro*gamro*tan(alp)
4433c
IF(amx.LT.(amk+amkz))
GO TO 100
4435c phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
4436c phspac=phspac*(alp2-alp1)
4438c tau-neutrino momentum
4441 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4442 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4446 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4448 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4451 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4452 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4453 pppi=sqrt((enq1-amk)*(enq1+amk))
4454 phspac=phspac*(4*pi)*(2*pppi/amx)
4455c charged pi momentum in rho rest frame
4456 CALL sphera(pppi,pkc)
4458c neutral pi momentum in rho rest frame
4462 exe=(pr(4)+pr(3))/amx
4463c pions boosted from rho rest frame to tau rest frame
4464 CALL bostr3(exe,pkc,pkc)
4465 CALL bostr3(exe,pkz,pkz)
4467 30 qq(i)=pkc(i)-pkz(i)
4468c qq transverse to pr
4469 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4470 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
447231 qq(i)=qq(i)-pr(i)*qqpks/pksd
4475 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4477 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4478 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4479 & +(gv**2-ga**2)*amtau*amnuta*qq2
4480 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4481 dgamt=1/(2.*amtau)*amplit*phspac
4483 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4491c ----------------------------------------------------------
4492c square of pion form factor
4493c ----------------------------------------------------------
4494 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4495 * ,ampiz,ampi,amro,gamro,ama1,gama1
4496 * ,amk,amkz,amkst,gamkst
4498 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4499 * ,ampiz,ampi,amro,gamro,ama1,gama1
4500 * ,amk,amkz,amkst,gamkst
4503 fpirk=cabs(fpikm(w,amk,amkz))**2
4504c fpirk=cabs(fpikmk(w,amk,amkz))**2
4506 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4507c **********************************************************
4509c **********************************************************
4511 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4515c ------------ parameters --------------------
4516 IF (init.EQ.0 )
THEN
4527c -----------------------------------------------
4529 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4535c initialize lund
COMMON
4538 SUBROUTINE dwrph(KTO,PHX)
4540c -------------------------
4542 IMPLICIT REAL*8 (a-h,o-z)
4549c
CASE of tau radiative decays.
4550c filling of the lund
COMMON block.
4553 IF (qhot(4).GT.1.e-5)
CALL dwluph(kto,qhot)
4556 SUBROUTINE dwluph(KTO,PHOT)
4557c---------------------------------------------------------------------
4558c lorentz transformation to cmsystem and
4559c updating of hepevt record
4561c called by : dexay1,(dekay1,dekay2)
4563c used when radiative corrections in decays are generated
4564c---------------------------------------------------------------------
4567 COMMON /taupos/ np1,np2
4570 IF (phot(4).LE.0.0)
RETURN
4572c position of decaying particle:
4573 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
4580 IF(ktos.GT.10) ktos=ktos-10
4581c boost and append photon(gamma is 22)
4582 CALL tralo4(ktos,phot,phot,am)
4583 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4588 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4589c ----------------------------------------------------------------------
4590c lorentz transformation to cmsystem and
4591c updating of hepevt record
4593c isgn = 1/-1 for tau-/tau+
4595c called by : dexay,(dekay1,dekay2)
4596c ----------------------------------------------------------------------
4598 REAL PNU(4),PWB(4),PEL(4),PNE(4)
4599 COMMON /taupos/ np1,np2
4601c position of decaying particle:
4608c tau neutrino(nu_tau is 16)
4609 CALL tralo4(kto,pnu,pnu,am)
4610 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4613 CALL tralo4(kto,pwb,pwb,am)
4614c
CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4617 CALL tralo4(kto,pel,pel,am)
4618 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4620c anti electron neutrino(nu_e is 12)
4621 CALL tralo4(kto,pne,pne,am)
4622 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4626 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4627c ----------------------------------------------------------------------
4628c lorentz transformation to cmsystem and
4629c updating of hepevt record
4631c isgn = 1/-1 for tau-/tau+
4633c called by : dexay,(dekay1,dekay2)
4634c ----------------------------------------------------------------------
4636 REAL PNU(4),PWB(4),PMU(4),PNM(4)
4637 COMMON /taupos/ np1,np2
4639c position of decaying particle:
4646c tau neutrino(nu_tau is 16)
4647 CALL tralo4(kto,pnu,pnu,am)
4648 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4651 CALL tralo4(kto,pwb,pwb,am)
4652c
CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4655 CALL tralo4(kto,pmu,pmu,am)
4656 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
4658c anti muon neutrino(nu_mu is 14)
4659 CALL tralo4(kto,pnm,pnm,am)
4660 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
4664 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
4665c ----------------------------------------------------------------------
4666c lorentz transformation to cmsystem and
4667c updating of hepevt record
4669c isgn = 1/-1 for tau-/tau+
4671c called by : dexay,(dekay1,dekay2)
4672c ----------------------------------------------------------------------
4675 COMMON /taupos/ np1,np2
4677c position of decaying particle:
4684c tau neutrino(nu_tau is 16)
4685 CALL tralo4(kto,pnu,pnu,am)
4686 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4688c charged pi meson(pi+ is 211)
4689 CALL tralo4(kto,ppi,ppi,am)
4690 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
4694 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
4695c ----------------------------------------------------------------------
4696c lorentz transformation to cmsystem and
4697c updating of hepevt record
4699c isgn = 1/-1 for tau-/tau+
4701c called by : dexay,(dekay1,dekay2)
4702c ----------------------------------------------------------------------
4704 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
4705 COMMON /taupos/ np1,np2
4707c position of decaying particle:
4714c tau neutrino(nu_tau is 16)
4715 CALL tralo4(kto,pnu,pnu,am)
4716 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4718c charged rho meson(rho+ is 213)
4719 CALL tralo4(kto,prho,prho,am)
4720 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
4722c charged pi meson(pi+ is 211)
4723 CALL tralo4(kto,pic,pic,am)
4724 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
4726c pi0 meson(pi0 is 111)
4727 CALL tralo4(kto,piz,piz,am)
4728 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
4732 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
4733c ----------------------------------------------------------------------
4734c lorentz transformation to cmsystem and
4735c updating of hepevt record
4737c isgn = 1/-1 for tau-/tau+
4738c jaa = 1 (2) for a_1- decay to pi+ 2pi- (pi- 2pi0)
4740c called by : dexay,(dekay1,dekay2)
4741c ----------------------------------------------------------------------
4743 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
4744 COMMON /taupos/ np1,np2
4746c position of decaying particle:
4753c tau neutrino(nu_tau is 16)
4754 CALL tralo4(kto,pnu,pnu,am)
4755 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4757c charged a_1 meson(a_1+ is 20213)
4758 CALL tralo4(kto,paa,paa,am)
4759 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
4761c two possible decays of the charged a1 meson
4764c a1 --> pi+ pi- pi- (or charged conjugate)
4766c pi minus(or c.c.) (pi+ is 211)
4767 CALL tralo4(kto,pim2,pim2,am)
4768 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
4770c pi minus(or c.c.) (pi+ is 211)
4771 CALL tralo4(kto,pim1,pim1,am)
4772 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
4774c pi plus(or c.c.) (pi+ is 211)
4775 CALL tralo4(kto,pipl,pipl,am)
4776 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
4778 ELSE IF (jaa.EQ.2)
THEN
4780c a1 --> pi- pi0 pi0(or charged conjugate)
4782c pi zero(pi0 is 111)
4783 CALL tralo4(kto,pim2,pim2,am)
4784 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
4786c pi zero(pi0 is 111)
4787 CALL tralo4(kto,pim1,pim1,am)
4788 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
4790c pi minus(or c.c.) (pi+ is 211)
4791 CALL tralo4(kto,pipl,pipl,am)
4792 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
4798 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
4799c ----------------------------------------------------------------------
4800c lorentz transformation to cmsystem and
4801c updating of hepevt record
4803c isgn = 1/-1 for tau-/tau+
4805c ----------------------------------------------------------------------
4808 COMMON /taupos/ np1,np2
4810c position of decaying particle
4817c tau neutrino(nu_tau is 16)
4818 CALL tralo4 (kto,pnu,pnu,am)
4819 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4822 CALL tralo4 (kto,pkk,pkk,am)
4823 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
4827 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
4828 COMMON / taukle / bra1,brk0,brk0b,brks
4829 real*4 bra1,brk0,brk0b,brks
4830c ----------------------------------------------------------------------
4831c lorentz transformation to cmsystem and
4832c updating of hepevt record
4834c isgn = 1/-1 for tau-/tau+
4835c jkst=10 (20) corresponds to k0b pi- (k- pi0) decay
4837c ----------------------------------------------------------------------
4839 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
4840 COMMON /taupos/ np1,np2
4842c position of decaying particle
4849c tau neutrino(nu_tau is 16)
4850 CALL tralo4(kto,pnu,pnu,am)
4851 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4853c charged k* meson(k*+ is 323)
4854 CALL tralo4(kto,pks,pks,am)
4855 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
4857c two possible decay modes of charged k*
4860c k*- --> pi- k0b(or charged conjugate)
4862c charged pi meson(pi+ is 211)
4863 CALL tralo4(kto,ppi,ppi,am)
4864 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
4867 IF (isgn.EQ.-1) bran=brk0
4868c k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
4870 IF(xio(1).GT.bran)
THEN
4876 CALL tralo4(kto,pkk,pkk,am)
4877 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
4879 ELSE IF(jkst.EQ.20)
THEN
4883c pi zero(pi0 is 111)
4884 CALL tralo4(kto,ppi,ppi,am)
4885 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
4887c charged k meson(k+ is 321)
4888 CALL tralo4(kto,pkk,pkk,am)
4889 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
4895 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
4896c ----------------------------------------------------------------------
4897c lorentz transformation to cmsystem and
4898c updating of hepevt record
4900c isgn = 1/-1 for tau-/tau+
4902c called by : dexay,(dekay1,dekay2)
4903c ----------------------------------------------------------------------
4905 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4906 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4908 COMMON /taupos/ np1,np2
4909 CHARACTER NAMES(NMODE)*31
4910 REAL PNU(4),PWB(4),PNPI(4,9)
4914c position of decaying particle
4921c tau neutrino(nu_tau is 16)
4922 CALL tralo4(kto,pnu,pnu,am)
4923 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4926 CALL tralo4(kto,pwb,pwb,am)
4927 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4931c get multiplicity of mode jnpi
4934 kfpi=lunpik(idffin(i,jnpi),-isgn)
4935c for charged conjugate
case, change charged pions only
4936c
IF(kfpi.NE.111)kfpi=kfpi*isgn
4940 CALL tralo4(kto,ppi,ppi,am)
4941 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
4947c ----------------------------------------------------------------------
4948c calculates mass of pp(double precision)
4951c ----------------------------------------------------------------------
4952 IMPLICIT REAL*8 (a-h,o-z)
4954 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
4956 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
4962c ----------------------------------------------------------------------
4963c calculates mass of pp
4966c ----------------------------------------------------------------------
4968 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
4969 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
4974c ----------------------------------------------------------------------
4976c used by : koralz radkor
4977c ----------------------------------------------------------------------
4978 IMPLICIT DOUBLE PRECISION (a-h,o-z)
4979 DATA pi /3.141592653589793238462643d0/
4981 IF(abs(y).LT.abs(x))
THEN
4983 IF(x.LE.0d0) the=pi-the
4985 the=acos(x/sqrt(x**2+y**2))
4991c ----------------------------------------------------------------------
4992* calculates angle in(0,2*pi) range out of x-y
4994c used by : koralz radkor
4995c ----------------------------------------------------------------------
4996 IMPLICIT DOUBLE PRECISION (a-h,o-z)
4997 DATA pi /3.141592653589793238462643d0/
4999 IF(abs(y).LT.abs(x))
THEN
5001 IF(x.LE.0d0) the=pi-the
5003 the=acos(x/sqrt(x**2+y**2))
5005 IF(y.LT.0d0) the=2d0*pi-the
5008 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5009c ----------------------------------------------------------------------
5012c ----------------------------------------------------------------------
5013 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5014 dimension pvec(4),qvec(4),rvec(4)
5022 qvec(2)= cs*rvec(2)-sn*rvec(3)
5023 qvec(3)= sn*rvec(2)+cs*rvec(3)
5027 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5028c ----------------------------------------------------------------------
5030c used by : koralz radkor
5031c ----------------------------------------------------------------------
5032 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5033 dimension pvec(4),qvec(4),rvec(4)
5040 qvec(1)= cs*rvec(1)+sn*rvec(3)
5042 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5046 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5047c ----------------------------------------------------------------------
5049c used by : koralz radkor
5050c ----------------------------------------------------------------------
5051 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5053 dimension pvec(4),qvec(4),rvec(4)
5059 qvec(1)= cs*rvec(1)-sn*rvec(2)
5060 qvec(2)= sn*rvec(1)+cs*rvec(2)
5064 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5065c ----------------------------------------------------------------------
5066c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5068c used by : tauola koralz(?)
5069c ----------------------------------------------------------------------
5070 real*4 pvec(4),qvec(4),rvec(4)
5083 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5084c ----------------------------------------------------------------------
5085c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5087c used by : koralz radkor
5088c ----------------------------------------------------------------------
5089 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5090 dimension pvec(4),qvec(4),rvec(4)
5104 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5105c ----------------------------------------------------------------------
5108c ----------------------------------------------------------------------
5109 real*4 pvec(4),qvec(4),rvec(4)
5117 qvec(2)= cs*rvec(2)-sn*rvec(3)
5118 qvec(3)= sn*rvec(2)+cs*rvec(3)
5121 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5122c ----------------------------------------------------------------------
5125c ----------------------------------------------------------------------
5126 IMPLICIT REAL*4(a-h,o-z)
5127 real*4 pvec(4),qvec(4),rvec(4)
5134 qvec(1)= cs*rvec(1)+sn*rvec(3)
5136 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5139 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5140c ----------------------------------------------------------------------
5143c ----------------------------------------------------------------------
5144 real*4 pvec(4),qvec(4),rvec(4)
5150 qvec(1)= cs*rvec(1)-sn*rvec(2)
5151 qvec(2)= sn*rvec(1)+cs*rvec(2)
5155 SUBROUTINE spherd(R,X)
5156c ----------------------------------------------------------------------
5157c generates uniformly three-vector x on sphere of radius r
5158c double precison version of sphera
5159c ----------------------------------------------------------------------
5160 real*8 r,x(4),pi,costh,sinth
5162 DATA pi /3.141592653589793238462643d0/
5166 sinth=sqrt(1 -costh**2)
5167 x(1)=r*sinth*cos(2*pi*rrr(2))
5168 x(2)=r*sinth*sin(2*pi*rrr(2))
5172 SUBROUTINE rotpox(THET,PHI,PP)
5173 IMPLICIT REAL*8 (a-h,o-z)
5174c ----------------------------------------------------------------------
5176c ----------------------------------------------------------------------
5179 CALL rotod2(thet,pp,pp)
5180 CALL rotod3( phi,pp,pp)
5183 SUBROUTINE sphera(R,X)
5184c ----------------------------------------------------------------------
5185c generates uniformly three-vector x on sphere of radius r
5187c called by : dphsxx,dadmpi,dadmkk
5188c ----------------------------------------------------------------------
5191 DATA pi /3.141592653589793238462643/
5195 sinth=sqrt(1.-costh**2)
5196 x(1)=r*sinth*cos(2*pi*rrr(2))
5197 x(2)=r*sinth*sin(2*pi*rrr(2))
5201 SUBROUTINE rotpol(THET,PHI,PP)
5202c ----------------------------------------------------------------------
5204c called by : dadmaa,dphsaa
5205c ----------------------------------------------------------------------
5208 CALL rotor2(thet,pp,pp)
5209 CALL rotor3( phi,pp,pp)
5212 SUBROUTINE ranmar(RVEC,LENV)
5213c ----------------------------------------------------------------------
5214c<<<<<
FUNCTION ranmar(IDUMM)
5215c cernlib v113, version with automatic default initialization
5216c transformed to
SUBROUTINE to be as in cernlib
5217c am.lutz november 1988, feb. 1989
5220c in report fsu-scri-87-50
5221c modified by f. james, 1988 and 1989, to generate a vector
5222c of pseudorandom numbers rvec of length lenv, and to put in
5223c the
COMMON block everything needed to specify currrent state,
5224c and to add input and output entry points rmarin, rmarut.
5226c unique random number used in the
program
5227c ----------------------------------------------------------------------
5228 COMMON / inout / inut,iout
5230 common/raset1/u(97),c,i97,j97
5231 parameter(modcns=1000000000)
5232 DATA ntot,ntot2,ijkl/-1,0,0/
5234 IF (ntot .GE. 0)
GO TO 50
5236c default initialization. user has called ranmar without rmarin.
5243 entry rmarin(ijklin, ntotin,ntot2n)
5244c initializing routine for ranmar, may be called before
5245c generating pseudorandom numbers with ranmar. the input
5246c values should be in the ranges: 0<=ijklin<=900 ooo ooo
5247c 0<=ntotin<=999 999 999
5248c 0<=ntot2n<<999 999 999
5249c to get the standard values in marsaglia-s paper, ijklin=54217137
5252 ntot = max(ntotin,0)
5253 ntot2= max(ntot2n,0)
5255c always come here to initialize
5258 kl = ijkl - 30082*ij
5259 i = mod(ij/177, 177) + 2
5260 j = mod(ij, 177) + 2
5261 k = mod(kl/169, 178) + 1
5263 WRITE(iout,201) ijkl,ntot,ntot2
5264 201
FORMAT(1x,
' RANMAR INITIALIZED: ',i10,2x,2i10)
5269 m = mod(mod(i*j,179)*k, 179)
5273 l = mod(53*l+1, 169)
5274 IF (mod(l*m,64) .GE. 32) s = s+t
5279 4 twom24 = 0.5*twom24
5281 cd = 7654321.*twom24
5282 cm = 16777213.*twom24
5285c complete initialization by skipping
5286c(ntot2*modcns + ntot) random numbers
5287 DO 45 loop2= 1, ntot2+1
5289 IF (loop2 .EQ. ntot2+1) now=ntot
5290 IF (now .GT. 0)
THEN
5291 WRITE (iout,
'(A,I15)')
' RMARIN SKIPPING OVER ',now
5292 DO 40 idum = 1, ntot
5294 IF (uni .LT. 0.) uni=uni+1.
5297 IF (i97 .EQ. 0) i97=97
5299 IF (j97 .EQ. 0) j97=97
5301 IF (c .LT. 0.) c=c+cm
5305 IF (kalled .EQ. 1)
RETURN
5307c normal entry to generate lenv random numbers
5309 DO 100 ivec= 1, lenv
5311 IF (uni .LT. 0.) uni=uni+1.
5314 IF (i97 .EQ. 0) i97=97
5316 IF (j97 .EQ. 0) j97=97
5318 IF (c .LT. 0.) c=c+cm
5320 IF (uni .LT. 0.) uni=uni+1.
5321c replace exact zeroes by uniform distr. *2**-24
5322 IF (uni .EQ. 0.)
THEN
5324c an exact zero here is very unlikely, but lets be safe.
5325 IF (uni .EQ. 0.) uni= twom24*twom24
5330 IF (ntot .GE. modcns)
THEN
5332 ntot = ntot - modcns
5335c entry to output current status
5336 entry rmarut(ijklut,ntotut,ntot2t)
5344 IMPLICIT REAL*8(a-h,o-z)
5345cern c304 version 29/07/71 dilog 59 c
5347 IF(x .LT.-1.0)
GO TO 1
5348 IF(x .LE. 0.5)
GO TO 2
5349 IF(x .EQ. 1.0)
GO TO 3
5350 IF(x .LE. 2.0)
GO TO 4
5354 z=z-0.5* log(abs(x))**2
5360 3 dilogt=1.64493406684822
5364 z=1.64493406684822 - log(x)* log(abs(t))
5365 5 y=2.66666666666666 *t+0.66666666666666
5366 b= 0.00000 00000 00001
5367 a=y*b +0.00000 00000 00004
5368 b=y*a-b+0.00000 00000 00011
5369 a=y*b-a+0.00000 00000 00037
5370 b=y*a-b+0.00000 00000 00121
5371 a=y*b-a+0.00000 00000 00398
5372 b=y*a-b+0.00000 00000 01312
5373 a=y*b-a+0.00000 00000 04342
5374 b=y*a-b+0.00000 00000 14437
5375 a=y*b-a+0.00000 00000 48274
5376 b=y*a-b+0.00000 00001 62421
5377 a=y*b-a+0.00000 00005 50291
5378 b=y*a-b+0.00000 00018 79117
5379 a=y*b-a+0.00000 00064 74338
5380 b=y*a-b+0.00000 00225 36705
5381 a=y*b-a+0.00000 00793 87055
5382 b=y*a-b+0.00000 02835 75385
5383 a=y*b-a+0.00000 10299 04264
5384 b=y*a-b+0.00000 38163 29463
5385 a=y*b-a+0.00001 44963 00557
5386 b=y*a-b+0.00005 68178 22718
5387 a=y*b-a+0.00023 20021 96094
5388 b=y*a-b+0.00100 16274 96164
5389 a=y*b-a+0.00468 63619 59447
5390 b=y*a-b+0.02487 93229 24228
5391 a=y*b-a+0.16607 30329 27855
5392 a=y*a-b+1.93506 43008 6996
5395c=======================================================================
5396c===================
END OF CPC PART ====================================
5397c=======================================================================
5400c-----------------------------------------------------------------------
5401c initialize rchl currents
5402c dummy routine for compatibility with new updates of tauola
5405c-----------------------------------------------------------------------
5406 SUBROUTINE inirchl(FLAG)
5411 25
FORMAT(1x,
"TAUOLA IniRChL: Fatal error, FLAG=",i2,
" but RChL currents missing")
5412 WRITE(*,*)
" in loaded version of TAUOLA-FORTRAN library."