45 COMMON / taubra / gamprt(30),jlist(30),nchan
52 IF(nchan.LE.0.OR.nchan.GT.30)
GOTO 902
59 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
64 9020
FORMAT(
' ----- JAKER: WRONG NCHAN')
67 SUBROUTINE dekay(KTO,HX)
99 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
105 COMMON /taupos/ np1,np2
106 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
107 REAL*4 GAMPMC ,GAMPER
108 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
110 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
112 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
115 CHARACTER NAMES(NMODE)*31
116 COMMON / inout / inut,iout
117 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4)
131 IF (iwarm.EQ.1) x=5/(iwarm-1)
133 WRITE(iout,7001) jak1,jak2
137 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN 138 CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
139 CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
140 CALL dadmpi(-1,idum,pdum,pdum1,pdum2)
141 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
142 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
143 CALL dadmkk(-1,idum,pdum,pdum1,pdum2)
144 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
145 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
151 ELSEIF(kto.EQ.1)
THEN 155 IF(iwarm.EQ.0)
GOTO 902
158 #elif defined (ALEPH) 163 CALL dekay1(0,h,isgn)
164 ELSEIF(kto.EQ.2)
THEN 168 IF(iwarm.EQ.0)
GOTO 902
171 #elif defined (ALEPH) 176 CALL dekay2(0,h,isgn)
177 ELSEIF(kto.EQ.11)
THEN 182 CALL dekay1(1,h,isgn)
183 ELSEIF(kto.EQ.12)
THEN 188 CALL dekay2(1,h,isgn)
189 ELSEIF(kto.EQ.100)
THEN 191 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN 192 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
193 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
194 CALL dadmpi( 1,idum,pdum,pdum1,pdum2)
195 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
196 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
197 CALL dadmkk( 1,idum,pdum,pdum1,pdum2)
198 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
199 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
200 WRITE(iout,7010) nev1,nev2,nevtot
201 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
203 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
214 7001
FORMAT(///1x,15(5h*****)
216 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
217 $ /,
' *', 25x,
'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
218 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
219 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
220 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
221 $ /,
' *', 25x,
'Physics initialization by ALEPH collab ',9x,1h*,
222 $ /,
' *', 25x,
'it is suggested to use this version ',9x,1h*,
223 $ /,
' *', 25x,
' with the help of the collab. advice ',9x,1h*,
224 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
225 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
227 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
228 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
229 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
230 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
231 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
232 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
233 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
234 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
235 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
236 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
237 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
238 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
240 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
241 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
242 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
243 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
244 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
245 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
246 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
247 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
248 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
249 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
251 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
252 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
253 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
255 7010
FORMAT(///1x,15(5h*****)
257 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
258 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
260 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
261 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
263 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
264 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
265 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
266 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
267 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
268 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
269 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
270 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
271 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
272 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
273 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
275 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
277 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
278 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
279 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
280 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
281 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
282 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
283 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
285 $ ,i10,2f12.7,a31 ,8x,1h*)
287 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
288 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
291 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
294 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
297 SUBROUTINE dekay1(IMOD,HH,ISGN)
301 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
302 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
303 REAL*4 GAMPMC ,GAMPER
304 COMMON / decp4 / pp1(4),pp2(4),kff1,kff2
308 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
309 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
310 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
311 REAL*4 GAMPMC ,GAMPER
314 REAL HV(4),PNU(4),PPI(4)
315 REAL PWB(4),PMU(4),PNM(4)
316 REAL PRHO(4),PIC(4),PIZ(4)
317 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
324 IF(jak1.EQ.-1)
RETURN 329 IF(jak1.EQ.0)
CALL jaker(jak)
331 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
332 ELSEIF(jak.EQ.2)
THEN 333 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
334 ELSEIF(jak.EQ.3)
THEN 335 CALL dadmpi(0, isgn,hv,ppi,pnu)
336 ELSEIF(jak.EQ.4)
THEN 337 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
338 ELSEIF(jak.EQ.5)
THEN 339 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
340 ELSEIF(jak.EQ.6)
THEN 341 CALL dadmkk(0, isgn,hv,pkk,pnu)
342 ELSEIF(jak.EQ.7)
THEN 343 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
345 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
351 ELSEIF(imd.EQ.1)
THEN 355 nevdec(jak)=nevdec(jak)+1
360 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
361 CALL dwrph(ktom,phot)
365 ELSEIF(jak.EQ.2)
THEN 366 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
367 CALL dwrph(ktom,phot)
371 ELSEIF(jak.EQ.3)
THEN 372 CALL dwlupi(1,isgn,ppi,pnu)
376 ELSEIF(jak.EQ.4)
THEN 377 CALL dwluro(1,isgn,pnu,prho,pic,piz)
381 ELSEIF(jak.EQ.5)
THEN 382 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
385 ELSEIF(jak.EQ.6)
THEN 386 CALL dwlukk(1,isgn,pkk,pnu)
389 ELSEIF(jak.EQ.7)
THEN 390 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
395 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
403 SUBROUTINE dekay2(IMOD,HH,ISGN)
407 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
408 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
409 REAL*4 GAMPMC ,GAMPER
410 COMMON / decp4 / pp1(4),pp2(4),kff1,kff2
414 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
415 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
416 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
417 REAL*4 GAMPMC ,GAMPER
420 REAL HV(4),PNU(4),PPI(4)
421 REAL PWB(4),PMU(4),PNM(4)
422 REAL PRHO(4),PIC(4),PIZ(4)
423 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
430 IF(jak2.EQ.-1)
RETURN 435 IF(jak2.EQ.0)
CALL jaker(jak)
437 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
438 ELSEIF(jak.EQ.2)
THEN 439 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
440 ELSEIF(jak.EQ.3)
THEN 441 CALL dadmpi(0, isgn,hv,ppi,pnu)
442 ELSEIF(jak.EQ.4)
THEN 443 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
444 ELSEIF(jak.EQ.5)
THEN 445 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
446 ELSEIF(jak.EQ.6)
THEN 447 CALL dadmkk(0, isgn,hv,pkk,pnu)
448 ELSEIF(jak.EQ.7)
THEN 449 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
451 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
456 ELSEIF(imd.EQ.1)
THEN 460 nevdec(jak)=nevdec(jak)+1
465 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
466 CALL dwrph(ktom,phot)
470 ELSEIF(jak.EQ.2)
THEN 471 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
472 CALL dwrph(ktom,phot)
476 ELSEIF(jak.EQ.3)
THEN 477 CALL dwlupi(2,isgn,ppi,pnu)
481 ELSEIF(jak.EQ.4)
THEN 482 CALL dwluro(2,isgn,pnu,prho,pic,piz)
486 ELSEIF(jak.EQ.5)
THEN 487 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
490 ELSEIF(jak.EQ.6)
THEN 491 CALL dwlukk(2,isgn,pkk,pnu)
494 ELSEIF(jak.EQ.7)
THEN 495 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
500 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
508 SUBROUTINE dexay(KTO,POL)
522 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
523 REAL*4 GAMPMC ,GAMPER
524 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
526 COMMON /taupos/ np1,np2
527 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
529 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
531 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
534 CHARACTER NAMES(NMODE)*31
535 COMMON / inout / inut,iout
537 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
551 WRITE(iout, 7001) jak1,jak2
555 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN 556 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
557 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
558 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
559 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
560 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
561 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
562 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
563 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
569 ELSEIF(kto.EQ.1)
THEN 574 IF(iwarm.EQ.0)
GOTO 902
577 CALL dexay1(kto,jak1,jakp,pol,isgn)
578 ELSEIF(kto.EQ.2)
THEN 583 IF(iwarm.EQ.0)
GOTO 902
584 isgn=-idff/iabs(idff)
586 CALL dexay1(kto,jak2,jakm,pol,isgn)
587 ELSEIF(kto.EQ.100)
THEN 589 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN 590 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
591 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
592 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
593 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
594 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
595 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
596 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
597 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
598 WRITE(iout,7010) nev1,nev2,nevtot
599 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
601 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
608 7001
FORMAT(///1x,15(5h*****)
610 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
611 $ /,
' *', 25x,
'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
612 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
613 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
614 $ /,
' *', 25x,
'Physics initialization by ALEPH collab ',9x,1h*,
615 $ /,
' *', 25x,
'it is suggested to use this version ',9x,1h*,
616 $ /,
' *', 25x,
' with the help of the collab. advice ',9x,1h*,
617 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
618 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
620 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
621 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
622 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
623 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
624 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
625 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
626 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
627 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
628 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
629 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
630 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
631 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
633 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
634 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
635 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
636 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
637 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
638 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
639 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
640 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
642 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
643 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
644 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
645 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
646 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
650 7010
FORMAT(///1x,15(5h*****)
652 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
653 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
655 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
656 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
658 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
659 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
660 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
661 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
662 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
663 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
664 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
665 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
666 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
667 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
668 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
670 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
672 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
673 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
674 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
675 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
676 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
677 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
678 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
680 $ ,i10,2f12.7,a31 ,8x,1h*)
682 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
683 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
685 902
WRITE(iout, 9020)
686 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
688 910
WRITE(iout, 9100)
689 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
692 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
698 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
699 REAL*4 GAMPMC ,GAMPER
700 COMMON / inout / inut,iout
703 REAL PRHO(4),PIC(4),PIZ(4)
704 REAL PWB(4),PMU(4),PNM(4)
705 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
711 IF(jakin.EQ.-1)
RETURN 718 IF(jak.EQ.0)
CALL jaker(jak)
721 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
722 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
723 CALL dwrph(kto,phot )
724 ELSEIF(jak.EQ.2)
THEN 725 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
726 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
727 CALL dwrph(kto,phot )
728 ELSEIF(jak.EQ.3)
THEN 729 CALL dexpi(0, isgn,polar,ppi,pnu)
730 CALL dwlupi(kto,isgn,ppi,pnu)
731 ELSEIF(jak.EQ.4)
THEN 732 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
733 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
734 ELSEIF(jak.EQ.5)
THEN 735 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
736 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
737 ELSEIF(jak.EQ.6)
THEN 738 CALL dexkk(0, isgn,polar,pkk,pnu)
739 CALL dwlukk(kto,isgn,pkk,pnu)
740 ELSEIF(jak.EQ.7)
THEN 741 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
742 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
745 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
746 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
748 nevdec(jak)=nevdec(jak)+1
750 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
757 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
763 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
766 ELSEIF(mode.EQ. 0)
THEN 769 IF(iwarm.EQ.0)
GOTO 902
770 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
771 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
774 IF(rn(1).GT.wt)
GOTO 300
776 ELSEIF(mode.EQ. 1)
THEN 778 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
784 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
787 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
796 COMMON / inout / inut,iout
797 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
803 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
806 ELSEIF(mode.EQ. 0)
THEN 809 IF(iwarm.EQ.0)
GOTO 902
810 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
811 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
814 IF(rn(1).GT.wt)
GOTO 300
816 ELSEIF(mode.EQ. 1)
THEN 818 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
823 902
WRITE(iout, 9020)
824 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
827 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
832 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
833 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
834 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
835 * ,ampiz,ampi,amro,gamro,ama1,gama1
836 * ,amk,amkz,amkst,gamkst
838 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
839 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
840 * ,AMK,AMKZ,AMKST,GAMKST
841 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
842 REAL*4 GAMPMC ,GAMPER
847 COMMON / inout / inut,iout
852 REAL HHV(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4)
853 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
856 DATA pi /3.141592653589793238462643/
869 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
870 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
874 ELSEIF(mode.EQ. 0)
THEN 877 IF(iwarm.EQ.0)
GOTO 902
879 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
885 IF(wt.GT.wtmax) nevovr=nevovr+1
886 IF(rn*wtmax.GT.wt)
GOTO 300
893 CALL rotor2(thet,pnu,pnu)
894 CALL rotor3( phi,pnu,pnu)
895 CALL rotor2(thet,pwb,pwb)
896 CALL rotor3( phi,pwb,pwb)
897 CALL rotor2(thet,q1,q1)
898 CALL rotor3( phi,q1,q1)
899 CALL rotor2(thet,q2,q2)
900 CALL rotor3( phi,q2,q2)
901 CALL rotor2(thet,hv,hv)
902 CALL rotor3( phi,hv,hv)
903 CALL rotor2(thet,phx,phx)
904 CALL rotor3( phi,phx,phx)
906 44 hhv(i)=-isgn*hv(i)
909 ELSEIF(mode.EQ. 1)
THEN 911 IF(nevraw.EQ.0)
RETURN 912 pargam=swt/float(nevraw+1)
914 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
916 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
924 7010
FORMAT(///1x,15(5h*****)
925 $ /,
' *', 25x,
'******** DADMEL FINAL REPORT ******** ',9x,1h*
926 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
927 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
928 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
929 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
930 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
931 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
932 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
933 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
935 902
WRITE(iout, 9020)
936 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
939 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
941 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
942 * ,ampiz,ampi,amro,gamro,ama1,gama1
943 * ,amk,amkz,amkst,gamkst
945 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
946 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
947 * ,AMK,AMKZ,AMKST,GAMKST
948 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
949 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
950 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
951 REAL*4 GAMPMC ,GAMPER
952 COMMON / inout / inut,iout
954 REAL HHV(4),HV(4),PNU(4),PWB(4),Q1(4),Q2(4)
955 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
958 DATA pi /3.141592653589793238462643/
971 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
972 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
976 ELSEIF(mode.EQ. 0)
THEN 979 IF(iwarm.EQ.0)
GOTO 902
981 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
987 IF(wt.GT.wtmax) nevovr=nevovr+1
988 IF(rn*wtmax.GT.wt)
GOTO 300
993 CALL rotor2(thet,pnu,pnu)
994 CALL rotor3( phi,pnu,pnu)
995 CALL rotor2(thet,pwb,pwb)
996 CALL rotor3( phi,pwb,pwb)
997 CALL rotor2(thet,q1,q1)
998 CALL rotor3( phi,q1,q1)
999 CALL rotor2(thet,q2,q2)
1000 CALL rotor3( phi,q2,q2)
1001 CALL rotor2(thet,hv,hv)
1002 CALL rotor3( phi,hv,hv)
1003 CALL rotor2(thet,phx,phx)
1004 CALL rotor3( phi,phx,phx)
1006 44 hhv(i)=-isgn*hv(i)
1009 ELSEIF(mode.EQ. 1)
THEN 1011 IF(nevraw.EQ.0)
RETURN 1012 pargam=swt/float(nevraw+1)
1014 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1016 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1024 7010
FORMAT(///1x,15(5h*****)
1025 $ /,
' *', 25x,
'******** DADMMU FINAL REPORT ******** ',9x,1h*
1026 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
1027 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
1028 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1029 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
1030 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1031 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1032 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
1033 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
1034 $ /,1x,15(5h*****)/)
1035 902
WRITE(iout, 9020)
1036 9020
FORMAT(
' ----- DADMMU: LACK OF INITIALISATION')
1039 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1045 REAL*4 HVX(4),PAAX(4),XAX(4),QPX(4),XNX(4)
1046 REAL*8 HV(4),PH(4),PAA(4),XA(4),QP(4),XN(4)
1049 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1060 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1066 REAL*4 HVX(4),PAAX(4),XAX(4),QPX(4),XNX(4)
1067 REAL*8 HV(4),PH(4),PAA(4),XA(4),QP(4),XN(4)
1070 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1081 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
1082 IMPLICIT REAL*8 (a-h,o-z)
1087 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1088 * ,ampiz,ampi,amro,gamro,ama1,gama1
1089 * ,amk,amkz,amkst,gamkst
1091 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1092 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1093 * ,AMK,AMKZ,AMKST,GAMKST
1094 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1095 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1097 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1098 REAL*4 GAMPMC ,GAMPER
1100 COMMON / inout / inut,iout
1101 COMMON / taurad / xk0dec,itdkrc
1103 REAL*8 HV(4),PT(4),PH(4),PAA(4),XA(4),QP(4),XN(4)
1107 DATA pi /3.141592653589793238462643d0/
1111 xlam(x,y,z)=sqrt((x-y-z)**2-4.0*y*z)
1117 phspac=1./2**17/pi**8
1127 IF (ielmu.EQ.1)
THEN 1134 IF ( itdkrc.EQ.0) prhard=0d0
1136 IF(prsoft.LT.0.1)
THEN 1137 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1142 ihard=(rr5.GT.prsoft)
1146 ams1=(amu+amnuta)**2
1149 xl1=log(xk1/2/xk0dec)
1154 phspac=phspac*ams2*xl1*xk
1155 phspac=phspac/prhard
1158 phspac=phspac*2**6*pi**3
1159 phspac=phspac/prsoft
1168 am2sq=ams1+ rr2*(ams2-ams1)
1170 phspac=phspac*(ams2-ams1)
1172 enq1=(am2sq+amnuta**2)/(2*am2)
1173 enq2=(am2sq-amnuta**2)/(2*am2)
1174 ppi= enq1**2-amnuta**2
1175 pppi=sqrt(abs(enq1**2-amnuta**2))
1176 phspac=phspac*(4*pi)*(2*pppi/am2)
1178 CALL spherd(pppi,xn)
1188 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1189 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1190 ppi = pr(4)**2-am2**2
1194 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1196 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1198 exe=(pr(4)+pr(3))/am2
1199 CALL bostd3(exe,xn,xn)
1200 CALL bostd3(exe,xa,xa)
1204 eps=4*(amu/amtax)**2
1205 xl1=log((2+eps)/eps)
1207 eta =exp(xl1*rr3+xl0)
1210 phspac=phspac*xl1/2*eta
1212 CALL rotpox(thet,phi,xn)
1213 CALL rotpox(thet,phi,xa)
1214 CALL rotpox(thet,phi,qp)
1215 CALL rotpox(thet,phi,pr)
1221 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1222 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1223 ppi = paa(4)**2-am3**2
1224 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1232 exe=(paa(4)+paa(3))/am3
1233 CALL bostd3(exe,xn,xn)
1234 CALL bostd3(exe,xa,xa)
1235 CALL bostd3(exe,qp,qp)
1236 CALL bostd3(exe,pr,pr)
1238 thet =acos(-1.+2*rr3)
1240 CALL rotpox(thet,phi,xn)
1241 CALL rotpox(thet,phi,xa)
1242 CALL rotpox(thet,phi,qp)
1243 CALL rotpox(thet,phi,pr)
1258 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1259 dgamt=1/(2.*amtax)*amplit*phspac
1261 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1262 IMPLICIT REAL*8 (a-h,o-z)
1268 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1269 * ,ampiz,ampi,amro,gamro,ama1,gama1
1270 * ,amk,amkz,amkst,gamkst
1272 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1273 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1274 * ,AMK,AMKZ,AMKST,GAMKST
1275 REAL*8 HV(4),QP(4),XN(4),XA(4),XK(4)
1279 IF(xk(4).LT.0.1d0*ak0)
THEN 1280 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1282 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1286 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1299 IMPLICIT REAL*8(a-h,o-z)
1300 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1301 * ,ampiz,ampi,amro,gamro,ama1,gama1
1302 * ,amk,amkz,amkst,gamkst
1304 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1305 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1306 * ,AMK,AMKZ,AMKST,GAMKST
1307 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1308 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1309 COMMON / qedprm /alfinv,alfpi,xk0
1310 REAL*8 ALFINV,ALFPI,XK0
1311 REAL*8 QP(4),XN(4),XA(4),XK(4)
1314 REAL*8 S0(3),RXA(3),RXK(3),RQP(3)
1315 DATA pi /3.141592653589793238462643d0/
1321 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1329 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1331 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1332 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1334 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1335 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1336 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1338 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1339 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1349 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1351 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1352 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1353 const4=256*pi/alphai*gf**2
1354 IF (itdkrc.EQ.0) const4=0d0
1357 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1359 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1360 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1361 5 hv(i)=s0(i)/s1-1.d0
1364 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1378 IMPLICIT REAL*8(a-h,o-z)
1379 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1380 * ,ampiz,ampi,amro,gamro,ama1,gama1
1381 * ,amk,amkz,amkst,gamkst
1383 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1384 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1385 * ,AMK,AMKZ,AMKST,GAMKST
1386 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1387 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1388 COMMON / qedprm /alfinv,alfpi,xk0
1389 REAL*8 ALFINV,ALFPI,XK0
1390 dimension qp(4),xn(4),xa(4)
1393 REAL*8 RXA(3),RXN(3),RQP(3)
1394 REAL*8 BORNPL(3),AM3POL(3),XM3POL(3)
1395 DATA pi /3.141592653589793238462643d0/
1408 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1409 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1411 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1415 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1417 w0=(xn(4)+xa(4))/tmass
1429 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1430 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1431 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1432 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1433 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1434 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1437 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1438 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1439 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1444 const3=1/(2*alphai*pi)*64*gf**2
1445 IF (itdkrc.EQ.0) const3=0d0
1446 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1447 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1450 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1451 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1452 born= 32*(gfermi**2/2.)*brak
1454 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1455 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1456 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1457 am3pol(i)=xm3pol(i)*const3
1460 & (gv+ga)**2*tmass*xnxa*qp(i)
1461 & -(gv-ga)**2*tmass*qpxn*xa(i)
1462 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1463 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1465 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1467 IF (thb/born.LT.0.1d0)
THEN 1468 print *,
'ERROR IN THB, THB/BORN=',thb/born
1477 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1484 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
1488 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1491 ELSEIF(mode.EQ. 0)
THEN 1494 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1495 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1498 IF(rn(1).GT.wt)
GOTO 300
1500 ELSEIF(mode.EQ. 1)
THEN 1502 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1508 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1510 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1511 * ,ampiz,ampi,amro,gamro,ama1,gama1
1512 * ,amk,amkz,amkst,gamkst
1514 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1515 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1516 * ,AMK,AMKZ,AMKST,GAMKST
1517 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1518 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1519 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1520 REAL*4 GAMPMC ,GAMPER
1521 COMMON / inout / inut,iout
1522 REAL PPI(4),PNU(4),HV(4)
1523 DATA pi /3.141592653589793238462643/
1528 ELSEIF(mode.EQ. 0)
THEN 1531 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1532 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1533 xpi= sqrt(epi**2-ampi**2)
1535 CALL sphera(xpi,ppi)
1543 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1544 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1545 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1547 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1550 ELSEIF(mode.EQ. 1)
THEN 1552 IF(nevtot.EQ.0)
RETURN 1558 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1560 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1561 $ -4*ampi**2*amnuta**2 )/amtau**2
1564 WRITE(iout, 7010) nevtot,gamm,rat,error
1571 7010
FORMAT(///1x,15(5h*****)
1572 $ /,
' *', 25x,
'******** DADMPI FINAL REPORT ******** ',9x,1h*
1573 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1574 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1575 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1576 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1577 $ /,1x,15(5h*****)/)
1579 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1588 COMMON / inout / inut,iout
1589 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
1595 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1599 ELSEIF(mode.EQ. 0)
THEN 1602 IF(iwarm.EQ.0)
GOTO 902
1603 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1604 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1609 IF(rn(1).GT.wt)
GOTO 300
1611 ELSEIF(mode.EQ. 1)
THEN 1613 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1619 902
WRITE(iout, 9020)
1620 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1623 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1625 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1626 * ,ampiz,ampi,amro,gamro,ama1,gama1
1627 * ,amk,amkz,amkst,gamkst
1629 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1630 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1631 * ,AMK,AMKZ,AMKST,GAMKST
1632 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1633 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1634 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1635 REAL*4 GAMPMC ,GAMPER
1636 COMMON / inout / inut,iout
1638 REAL HV(4),PRO(4),PNU(4),PIC(4),PIZ(4)
1639 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
1642 DATA pi /3.141592653589793238462643/
1655 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1656 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1661 ELSEIF(mode.EQ. 0)
THEN 1664 IF(iwarm.EQ.0)
GOTO 902
1665 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1672 IF(wt.GT.wtmax) nevovr=nevovr+1
1673 IF(rn*wtmax.GT.wt)
GOTO 300
1675 costhe=-1.+2.*rrr(2)
1678 CALL rotor2(thet,pnu,pnu)
1679 CALL rotor3( phi,pnu,pnu)
1680 CALL rotor2(thet,pro,pro)
1681 CALL rotor3( phi,pro,pro)
1682 CALL rotor2(thet,pic,pic)
1683 CALL rotor3( phi,pic,pic)
1684 CALL rotor2(thet,piz,piz)
1685 CALL rotor3( phi,piz,piz)
1686 CALL rotor2(thet,hv,hv)
1687 CALL rotor3( phi,hv,hv)
1689 44 hhv(i)=-isgn*hv(i)
1692 ELSEIF(mode.EQ. 1)
THEN 1694 IF(nevraw.EQ.0)
RETURN 1695 pargam=swt/float(nevraw+1)
1697 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1699 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1707 7003
FORMAT(///1x,15(5h*****)
1708 $ /,
' *', 25x,
'******** DADMRO INITIALISATION ********',9x,1h*
1709 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1710 $ /,1x,15(5h*****)/)
1711 7010
FORMAT(///1x,15(5h*****)
1712 $ /,
' *', 25x,
'******** DADMRO FINAL REPORT ******** ',9x,1h*
1713 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1714 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1715 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1716 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1717 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1718 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1719 $ /,1x,15(5h*****)/)
1720 902
WRITE(iout, 9020)
1721 9020
FORMAT(
' ----- DADMRO: LACK OF INITIALISATION')
1724 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1729 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1730 * ,ampiz,ampi,amro,gamro,ama1,gama1
1731 * ,amk,amkz,amkst,gamkst
1733 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1734 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1735 * ,AMK,AMKZ,AMKST,GAMKST
1736 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1737 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1738 REAL HV(4),PT(4),PN(4),PR(4),PIC(4),PIZ(4),QQ(4),RR1(1)
1739 DATA pi /3.141592653589793238462643/
1743 phspac=1./2**11/pi**5
1750 ams1=(ampi+ampiz)**2
1751 ams2=(amtau-amnuta)**2
1761 alp1=atan((ams1-amro**2)/amro/gamro)
1762 alp2=atan((ams2-amro**2)/amro/gamro)
1766 alp=alp1+rr1(1)*(alp2-alp1)
1767 amx2=amro**2+amro*gamro*tan(alp)
1769 IF(amx.LT.2.*ampi)
GO TO 100
1771 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1772 phspac=phspac*(alp2-alp1)
1777 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1778 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1782 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1784 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1787 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1788 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1789 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1790 phspac=phspac*(4*pi)*(2*pppi/amx)
1792 CALL sphera(pppi,pic)
1798 exe=(pr(4)+pr(3))/amx
1800 CALL bostr3(exe,pic,pic)
1801 CALL bostr3(exe,piz,piz)
1803 30 qq(i)=pic(i)-piz(i)
1806 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
1808 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
1809 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
1810 & +(gv**2-ga**2)*amtau*amnuta*qq2
1811 amplit=(gfermi*ccabib)**2*brak*2*fpirho(amx)
1812 dgamt=1/(2.*amtau)*amplit*phspac
1814 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1817 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1828 COMMON / inout / inut,iout
1829 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
1835 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1838 ELSEIF(mode.EQ. 0)
THEN 1841 IF(iwarm.EQ.0)
GOTO 902
1842 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1843 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1846 IF(rn(1).GT.wt)
GOTO 300
1848 ELSEIF(mode.EQ. 1)
THEN 1850 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1855 902
WRITE(iout, 9020)
1856 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
1859 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1863 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1864 * ,ampiz,ampi,amro,gamro,ama1,gama1
1865 * ,amk,amkz,amkst,gamkst
1867 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1868 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1869 * ,AMK,AMKZ,AMKST,GAMKST
1870 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1871 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
1872 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1873 REAL*4 GAMPMC ,GAMPER
1874 COMMON / inout / inut,iout
1876 REAL HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4)
1877 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
1880 DATA pi /3.141592653589793238462643/
1893 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1894 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1898 ELSEIF(mode.EQ. 0)
THEN 1901 IF(iwarm.EQ.0)
GOTO 902
1902 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1911 sswt=sswt+dble(wt)**2
1916 IF(wt.GT.wtmax) nevovr=nevovr+1
1917 IF(rn*wtmax.GT.wt)
GOTO 300
1919 costhe=-1.+2.*rrr(2)
1922 CALL rotpol(thet,phi,pnu)
1923 CALL rotpol(thet,phi,paa)
1924 CALL rotpol(thet,phi,pim1)
1925 CALL rotpol(thet,phi,pim2)
1926 CALL rotpol(thet,phi,pipl)
1927 CALL rotpol(thet,phi,hv)
1929 44 hhv(i)=-isgn*hv(i)
1932 ELSEIF(mode.EQ. 1)
THEN 1934 IF(nevraw.EQ.0)
RETURN 1935 pargam=swt/float(nevraw+1)
1937 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1939 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1947 7003
FORMAT(///1x,15(5h*****)
1948 $ /,
' *', 25x,
'******** DADMAA INITIALISATION ********',9x,1h*
1949 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1950 $ /,1x,15(5h*****)/)
1951 7010
FORMAT(///1x,15(5h*****)
1952 $ /,
' *', 25x,
'******** DADMAA FINAL REPORT ******** ',9x,1h*
1953 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
1954 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
1955 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1956 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
1957 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1958 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1959 $ /,1x,15(5h*****)/)
1960 902
WRITE(iout, 9020)
1961 9020
FORMAT(
' ----- DADMAA: LACK OF INITIALISATION')
1964 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1969 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1970 * ,ampiz,ampi,amro,gamro,ama1,gama1
1971 * ,amk,amkz,amkst,gamkst
1973 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1974 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1975 * ,AMK,AMKZ,AMKST,GAMKST
1976 COMMON / taukle / bra1,brk0,brk0b,brks
1977 REAL*4 BRA1,BRK0,BRK0B,BRKS
1978 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
1988 IF (rmod.LT.bra1)
THEN 2000 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2002 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
2009 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
2013 CALL dadmkk(-1,isgn,hv,pkk,pnu)
2016 ELSEIF(mode.EQ. 0)
THEN 2019 CALL dadmkk( 0,isgn,hv,pkk,pnu)
2020 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2023 IF(rn(1).GT.wt)
GOTO 300
2025 ELSEIF(mode.EQ. 1)
THEN 2027 CALL dadmkk( 1,isgn,hv,pkk,pnu)
2033 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
2038 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2039 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
2041 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2042 * ,ampiz,ampi,amro,gamro,ama1,gama1
2043 * ,amk,amkz,amkst,gamkst
2045 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
2046 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
2047 * ,AMK,AMKZ,AMKST,GAMKST
2049 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2050 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
2053 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2054 REAL*4 GAMPMC ,GAMPER
2055 COMMON / inout / inut,iout
2056 REAL PKK(4),PNU(4),HV(4)
2057 DATA pi /3.141592653589793238462643/
2062 ELSEIF(mode.EQ. 0)
THEN 2065 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
2066 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
2067 xkk= sqrt(ekk**2-amk**2)
2069 CALL sphera(xkk,pkk)
2077 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
2078 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
2079 & +(gv**2-ga**2)*amtau*amnuta*amk**2
2081 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
2084 ELSEIF(mode.EQ. 1)
THEN 2086 IF(nevtot.EQ.0)
RETURN 2093 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2095 $ sqrt((amtau**2-amk**2-amnuta**2)**2
2096 $ -4*amk**2*amnuta**2 )/amtau**2
2101 WRITE(iout, 7010) nevtot,gamm,rat,error
2108 7010
FORMAT(///1x,15(5h*****)
2109 $ /,
' *', 25x,
'******** DADMKK FINAL REPORT ********',9x,1h*
2110 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
2111 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
2112 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2113 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2114 $ /,1x,15(5h*****)/)
2116 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2128 COMMON / inout / inut,iout
2129 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
2136 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2140 ELSEIF(mode.EQ. 0)
THEN 2143 IF(iwarm.EQ.0)
GOTO 902
2144 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2145 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2150 IF(rn(1).GT.wt)
GOTO 300
2152 ELSEIF(mode.EQ. 1)
THEN 2154 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2160 902
WRITE(iout, 9020)
2161 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2164 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2166 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2167 * ,ampiz,ampi,amro,gamro,ama1,gama1
2168 * ,amk,amkz,amkst,gamkst
2170 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
2171 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
2172 * ,AMK,AMKZ,AMKST,GAMKST
2173 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2174 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
2175 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2176 REAL*4 GAMPMC ,GAMPER
2177 COMMON / taukle / bra1,brk0,brk0b,brks
2178 REAL*4 BRA1,BRK0,BRK0B,BRKS
2179 COMMON / inout / inut,iout
2181 REAL HV(4),PKS(4),PNU(4),PKK(4),PPI(4)
2182 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
2183 REAL*4 RRR(3),RMOD(1)
2185 DATA pi /3.141592653589793238462643/
2200 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2201 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2206 ELSEIF(mode.EQ. 0)
THEN 2208 IF(iwarm.EQ.0)
GOTO 902
2214 IF(rmod(1).LT.dec1)
THEN 2219 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2222 IF(wt.GT.wtmax) nevovr=nevovr+1
2226 IF(rn*wtmax.GT.wt)
GOTO 400
2228 costhe=-1.+2.*rrr(2)
2231 CALL rotor2(thet,pnu,pnu)
2232 CALL rotor3( phi,pnu,pnu)
2233 CALL rotor2(thet,pks,pks)
2234 CALL rotor3( phi,pks,pks)
2235 CALL rotor2(thet,pkk,pkk)
2236 CALL rotor3(phi,pkk,pkk)
2237 CALL rotor2(thet,ppi,ppi)
2238 CALL rotor3( phi,ppi,ppi)
2239 CALL rotor2(thet,hv,hv)
2240 CALL rotor3( phi,hv,hv)
2242 44 hhv(i)=-isgn*hv(i)
2245 ELSEIF(mode.EQ. 1)
THEN 2247 IF(nevraw.EQ.0)
RETURN 2248 pargam=swt/float(nevraw+1)
2250 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2252 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2260 7003
FORMAT(///1x,15(5h*****)
2261 $ /,
' *', 25x,
'******** DADMKS INITIALISATION ********',9x,1h*
2262 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2263 $ /,1x,15(5h*****)/)
2264 7010
FORMAT(///1x,15(5h*****)
2265 $ /,
' *', 25x,
'******** DADMKS FINAL REPORT ********',9x,1h*
2266 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2267 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2268 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2269 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2270 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2271 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2272 $ /,1x,15(5h*****)/)
2273 902
WRITE(iout, 9020)
2274 9020
FORMAT(
' ----- DADMKS: LACK OF INITIALISATION')
2277 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2286 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2287 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
2289 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2290 * ,ampiz,ampi,amro,gamro,ama1,gama1
2291 * ,amk,amkz,amkst,gamkst
2293 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
2294 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
2295 * ,AMK,AMKZ,AMKST,GAMKST
2297 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2298 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
2301 REAL HV(4),PT(4),PN(4),PKS(4),PKK(4),PPI(4),QQ(4),RR1(1)
2308 DATA pi /3.141592653589793238462643/
2312 phspac=1./2**11/pi**5
2324 ams2=(amtau-amnuta)**2
2330 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2331 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2332 alp=alp1+rr1(1)*(alp2-alp1)
2333 amx2=amkst**2+amkst*gamkst*tan(alp)
2335 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2337 phspac=phspac*(alp2-alp1)
2342 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2343 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2348 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2350 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2353 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2354 pppi=sqrt((enpi-ampi)*(enpi+ampi))
2355 phspac=phspac*(4*pi)*(2*pppi/amx)
2357 CALL sphera(pppi,ppi)
2362 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2363 exe=(pks(4)+pks(3))/amx
2365 CALL bostr3(exe,ppi,ppi)
2366 CALL bostr3(exe,pkk,pkk)
2368 30 qq(i)=ppi(i)-pkk(i)
2370 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2371 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2373 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
2376 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2378 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2379 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2380 & +(gv**2-ga**2)*amtau*amnuta*qq2
2384 fks=cabs(bwigm(amx2,amkst,gamkst,ampi,amkz))**2
2386 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2388 amplit=(gfermi*scabib)**2*brak*2*fks
2389 dgamt=1/(2.*amtau)*amplit*phspac
2391 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2394 ELSEIF(jkst.EQ.20)
THEN 2398 ams2=(amtau-amnuta)**2
2408 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2409 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2410 alp=alp1+rr1(1)*(alp2-alp1)
2411 amx2=amkst**2+amkst*gamkst*tan(alp)
2413 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2415 phspac=phspac*(alp2-alp1)
2420 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2421 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2425 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2427 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2430 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2431 pppi=sqrt((enpi-ampiz)*(enpi+ampiz))
2432 phspac=phspac*(4*pi)*(2*pppi/amx)
2434 CALL sphera(pppi,ppi)
2439 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2440 exe=(pks(4)+pks(3))/amx
2442 CALL bostr3(exe,ppi,ppi)
2443 CALL bostr3(exe,pkk,pkk)
2445 60 qq(i)=pkk(i)-ppi(i)
2447 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2448 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2450 61 qq(i)=qq(i)-pks(i)*qqpks/pksd
2453 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2455 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2456 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2457 & +(gv**2-ga**2)*amtau*amnuta*qq2
2461 fks=cabs(bwigm(amx2,amkst,gamkst,amk,ampiz))**2
2463 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2465 amplit=(gfermi*scabib)**2*brak*2*fks
2466 dgamt=1/(2.*amtau)*amplit*phspac
2468 70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2476 SUBROUTINE dphnpi(DGAMT,HV,PN,PR,PPI,JNPI)
2478 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2484 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2485 * ,ampiz,ampi,amro,gamro,ama1,gama1
2486 * ,amk,amkz,amkst,gamkst
2488 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
2489 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
2490 * ,AMK,AMKZ,AMKST,GAMKST
2491 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2492 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
2493 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2495 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2497 CHARACTER NAMES(NMODE)*31
2500 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2502 CHARACTER NAMES(NMODE)*31
2507 REAL PN(4),PR(4),PPI(4,9),HV(4)
2508 REAL PV(5,9),PT(4),UE(3),BE(3)
2509 REAL*4 RRR(9),RORD(9),RR1(1)
2512 DATA pi /3.141592653589793238462643/
2513 DATA dpar/2.,5.,15.,60.,250.,1500.,1.2e4,1.2e5/
2516 pawt(a,b,c)=sqrt(max(0.,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.*a)
2518 REAL*8 PN(4),PR(4),PPI(4,9),HV(4)
2519 REAL*4 PNX(4),PRX(4),PPIX(4,9),HVX(4)
2520 REAL*8 PV(5,9),PT(4),UE(3),BE(3)
2521 REAL*8 PAWT,AMX,AMS1,AMS2,PA,PHS,PHSMAX,PMIN,PMAX
2525 REAL*8 GAM,BEP,PHI,A,B,C
2527 REAL*4 RRR(9),RRX(2),RN(1),RR2(1)
2529 DATA pi /3.141592653589793238462643/
2530 DATA wetmax /20*1d-15/
2535 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2538 ampik(i,j)=dcdmas(idffin(i,j))
2543 IF ((jnpi.LE.0).OR.jnpi.GT.20)
THEN 2544 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2562 phspac = 1./2.**5 /pi**2
2564 4 ps =ps+ampik(i,jnpi)
2571 ams2=(amtau-amnuta)**2
2575 amx2=ams1+ rr1(1)*(ams2-ams1)
2577 amx2=ams1+ rr2(1)*(ams2-ams1)
2581 phspac=phspac * (ams2-ams1)
2586 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2587 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2591 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2593 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2600 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2606 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2607 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2610 amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2611 dgamt=1./(2.*amtau)*amplit*phspac
2615 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2616 phsmax = 1./dpar(nd-2)
2623 pv(5,nd)=ampik(nd,jnpi)
2625 pmax=amw-ps+ampik(nd,jnpi)
2628 pmax=pmax+ampik(il,jnpi)
2629 pmin=pmin+ampik(il+1,jnpi)
2631 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))
2633 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2635 CALL ranmar(rrr,nd-1)
2639 IF(rsav.LE.rord(jl))
GOTO 260
2640 250 rord(jl+1)=rord(jl)
2645 pv(5,il)=pv(5,il+1)+ampik(il,jnpi)
2646 & +(rord(il)-rord(il+1))*(pv(5,1)-ps)
2647 270 phs=phs*pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2649 IF(phs.LT.rn*phsmax)
GOTO 240
2651 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2658 223 ams1=ams1+ampik(jl,jnpi)
2660 amx =(amx-ampik(il,jnpi))
2662 phsmax=phsmax * (ams2-ams1)
2669 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2671 CALL ranmar(rrr,nd-2)
2675 231 ams1=ams1+ampik(jl,jnpi)
2677 ams2=(amx-ampik(il,jnpi))**2
2679 amx2=ams1+ rr1*(ams2-ams1)
2682 phspac=phspac * (ams2-ams1)
2684 phs=phs* (ams2-ams1)
2685 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2686 phs =phs *pa/pv(5,il)
2688 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2689 phs =phs *pa/pv(5,nd-1)
2691 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2692 IF (ncont.EQ.500 000)
THEN 2695 xnpi=xnpi+ampik(kk,jnpi)
2697 WRITE(6,*)
'ROUNDING INSTABILITY IN DPHNPI ?' 2698 WRITE(6,*)
'AMW=',amw,
'XNPI=',xnpi
2699 WRITE(6,*)
'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT' 2700 WRITE(6,*)
'PHS=',phs,
'PHSMAX=',phsmax
2703 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs)
GO TO 100
2706 280
DO 300 il=1,nd-1
2707 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2712 ue(1)=sqrt(1.-ue(3)**2)*cos(phi)
2713 ue(2)=sqrt(1.-ue(3)**2)*sin(phi)
2718 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2719 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2723 290 pv(j,il+1)=-pa*ue(j)
2724 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2725 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2726 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2730 310 ppi(j,nd)=pv(j,nd)
2733 320 be(j)=pv(j,il)/pv(4,il)
2734 gam=pv(4,il)/pv(5,il)
2736 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2739 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.+gam)+ppi(4,i))*be(j)
2741 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2743 ppi(4,i)=gam*(ppi(4,i)+bep)
2763 FUNCTION sigee(Q2,JNP)
2782 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2783 * ,ampiz,ampi,amro,gamro,ama1,gama1
2784 * ,amk,amkz,amkst,gamkst
2786 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
2787 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
2788 * ,AMK,AMKZ,AMKST,GAMKST
2792 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2793 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2794 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2795 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2798 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2801 DATA pi /3.141592653589793238462643/
2811 IF (ampi.LT.0.139) ampi = 0.1395675
2817 datsig(i,2) = datsig(i,2)/2.
2818 datsig(i,1) = datsig(i,1) + datsig(i,2)
2824 IF(t . gt. s-ampi )
GO TO 201
2826 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2827 fact = fact * (datsig(j,1)+datsig(j+1,1))
2828 200 datsig(i,3) = datsig(i,3) + fact
2829 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2830 datsig(i,4) = datsig(i,3)
2831 datsig(i,6) = datsig(i,5)
2834 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2840 sigee=datsig(1,jnpi)+
2841 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2842 ELSEIF(q.LT.1.8)
THEN 2845 IF(q.LT.qmax)
GO TO 2
2848 2 sigee=datsig(i,jnpi)+
2849 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2850 ELSEIF(q.GT.1.8)
THEN 2851 sigee=datsig(17,jnpi)+
2852 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2854 IF(sigee.LT..0) sigee=0.
2856 sigee = sigee/(6.*pi**2*sig0)
2861 FUNCTION sigold(Q2,JNPI)
2880 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2881 * ,ampiz,ampi,amro,gamro,ama1,gama1
2882 * ,amk,amkz,amkst,gamkst
2884 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
2885 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
2886 * ,AMK,AMKZ,AMKST,GAMKST
2890 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2891 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2892 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2893 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2895 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2897 DATA pi /3.141592653589793238462643/
2905 datsig(i,2) = datsig(i,2)/2.
2906 datsig(i,1) = datsig(i,1) + datsig(i,2)
2912 IF(t . gt. s-ampi )
GO TO 201
2914 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2915 fact = fact * (datsig(j,1)+datsig(j+1,1))
2916 200 datsig(i,3) = datsig(i,3) + fact
2917 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2920 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2926 sigee=datsig(1,jnpi)+
2927 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2928 ELSEIF(q.LT.1.8)
THEN 2931 IF(q.LT.qmax)
GO TO 2
2934 2 sigee=datsig(i,jnpi)+
2935 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2936 ELSEIF(q.GT.1.8)
THEN 2937 sigee=datsig(17,jnpi)+
2938 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2940 IF(sigee.LT..0) sigee=0.
2942 sigee = sigee/(6.*pi**2*sig0)
2947 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2952 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2954 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2956 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2959 CHARACTER NAMES(NMODE)*31
2961 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PNPI(4,9)
2968 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
2969 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
2970 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
2972 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2984 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
3001 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3002 * ,ampiz,ampi,amro,gamro,ama1,gama1
3003 * ,amk,amkz,amkst,gamkst
3005 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3006 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
3007 * ,AMK,AMKZ,AMKST,GAMKST
3008 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3009 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3010 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
3013 DATA pi /3.141592653589793238462643/
3015 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3020 phspac=1./2**17/pi**8
3030 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
3031 $ amrx,gamrx,amra,gamra,amrb,gamrb)
3032 IF (ichan.EQ.1)
THEN 3035 ELSEIF (ichan.EQ.2)
THEN 3044 ams1=(amp1+amp2+amp3)**2
3045 ams2=(amtau-amnuta)**2
3051 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3052 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3053 alp=alp1+rr1*(alp2-alp1)
3054 am3sq =amrx**2+amrx*gamrx*tan(alp)
3056 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3057 phspac=phspac*(alp2-alp1)
3062 IF (ichan.LE.2)
THEN 3068 alp1=atan((ams1-amra**2)/amra/gamra)
3069 alp2=atan((ams2-amra**2)/amra/gamra)
3070 alp=alp1+rr2*(alp2-alp1)
3071 am2sq =amra**2+amra*gamra*tan(alp)
3083 am2sq=ams1+ rr2*(ams2-ams1)
3092 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
3093 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
3094 ppi= enq1**2-amp3**2
3095 pppi=sqrt(abs(enq1**2-amp3**2))
3097 phf1=(4*pi)*(2*pppi/am2)
3103 CALL sphera(pppi,pipl)
3121 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
3122 pr(3)= sqrt(abs(pr(4)**2-am2**2))
3123 ppi = pr(4)**2-am2**2
3127 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
3129 phf2=(4*pi)*(2*pr(3)/am3)
3135 exe=(pr(4)+pr(3))/am2
3136 CALL bostr3(exe,pipl,pipl)
3137 CALL bostr3(exe,pim1,pim1)
3144 thet =acos(-1.+2*rr3)
3146 CALL rotpol(thet,phi,pipl)
3147 CALL rotpol(thet,phi,pim1)
3148 CALL rotpol(thet,phi,pim2)
3149 CALL rotpol(thet,phi,pr)
3155 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
3156 paa(3)= sqrt(abs(paa(4)**2-am3**2))
3157 ppi = paa(4)**2-am3**2
3158 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3162 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
3168 alp1=atan((ams1-amra**2)/amra/gamra)
3169 alp2=atan((ams2-amra**2)/amra/gamra)
3170 xpro = (pim1(3)+pipl(3))**2
3171 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
3172 am2sq=-xpro+(pim1(4)+pipl(4))**2
3174 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3175 ff1 =ff1 *(alp2-alp1)
3177 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
3179 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
3180 xjaje=gg1*(ams2-ams1)
3184 alp1=atan((ams1-amrb**2)/amrb/gamrb)
3185 alp2=atan((ams2-amrb**2)/amrb/gamrb)
3186 xpro = (pim2(3)+pipl(3))**2
3187 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
3188 am2sq=-xpro+(pim2(4)+pipl(4))**2
3189 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
3190 ff2 =ff2 *(alp2-alp1)
3191 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
3192 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
3193 xjadw=gg2*(ams2-ams1)
3200 IF (ichan.EQ.2)
THEN 3205 IF (xjac1.NE.0.0) a1=prob1/xjac1
3206 IF (xjac2.NE.0.0) a2=prob2/xjac2
3207 IF (xjac3.NE.0.0) a3=prob3/xjac3
3209 IF (a1+a2+a3.NE.0.0)
THEN 3210 phspac=phspac/(a1+a2+a3)
3222 exe=(paa(4)+paa(3))/am3
3223 CALL bostr3(exe,pipl,pipl)
3224 CALL bostr3(exe,pim1,pim1)
3225 CALL bostr3(exe,pim2,pim2)
3226 CALL bostr3(exe,pr,pr)
3229 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3233 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3235 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN 3246 dgamt=1/(2.*amtau)*amplit*phspac
3248 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3258 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3259 * ,ampiz,ampi,amro,gamro,ama1,gama1
3260 * ,amk,amkz,amkst,gamkst
3262 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3263 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
3264 * ,AMK,AMKZ,AMKST,GAMKST
3265 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3266 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3267 COMMON /testa1/ keya1
3268 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3269 REAL PAA(4),VEC1(4),VEC2(4)
3270 REAL PIVEC(4),PIAKS(4),HVM(4)
3271 COMPLEX BWIGN,HADCUR(4),FPIK
3278 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3282 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3284 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3285 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3286 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3287 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3288 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3290 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3291 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3292 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3293 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3295 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3296 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3298 IF (keya1.EQ.1)
THEN 3302 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3304 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3305 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3306 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3309 fnorm=2.0*sqrt(2.)/3.0/fpi
3310 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3312 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3313 $ *(cmplx(vec1(i))*fpik(xmro1)
3314 $ +cmplx(vec2(i))*fpik(xmro2))
3319 CALL clvec(hadcur,pn,pivec)
3320 CALL claxi(hadcur,pn,piaks)
3321 CALL clnut(hadcur,brakm,hvm)
3323 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3324 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3325 amplit=(gfermi*ccabib)**2*brak/2.
3330 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3331 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3341 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3342 * ,ampiz,ampi,amro,gamro,ama1,gama1
3343 * ,amk,amkz,amkst,gamkst
3345 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3346 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
3347 * ,AMK,AMKZ,AMKST,GAMKST
3349 IF (qkwa.LT.(amro+ampi)**2)
THEN 3350 gfun=4.1*(qkwa-9*ampiz**2)**3
3351 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3353 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3356 COMPLEX FUNCTION bwigs(S,M,G)
3361 REAL PI,PIM,QS,QM,W,GS,MK
3369 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3370 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3387 IF (s.GT.(pim+mk)**2)
THEN 3389 qs=p(s,pim**2,mk**2)
3390 qm=p(m**2,pim**2,mk**2)
3392 gs=g*(m/w)*(qs/qm)**3
3394 bw=m**2/cmplx(m**2-s,-m*gs)
3395 qpm=p(pm**2,pim**2,mk**2)
3396 g1=pg*(pm/w)*(qs/qpm)**3
3397 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3398 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3399 #elif defined (ALEPH) 3403 bwigs=m**2/cmplx(m**2-s,-m*gs)
3405 bwigs=m**2/cmplx(m**2-s,-m*gs)
3409 COMPLEX FUNCTION bwig(S,M,G)
3414 REAL PI,PIM,QS,QM,W,GS
3423 IF (s.GT.4.*pim**2)
THEN 3424 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3425 qm=sqrt(m**2/4.-pim**2)
3427 gs=g*(m/w)*(qs/qm)**3
3431 bwig=m**2/cmplx(m**2-s,-m*gs)
3434 COMPLEX FUNCTION fpik(W)
3439 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
3444 IF (init.EQ.0 )
THEN 3464 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3473 fpirho=cabs(fpik(w))**2
3475 SUBROUTINE clvec(HJ,PN,PIV)
3485 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3486 hh=
REAL(HJ(4)*CONJG(HJ(4))-HJ(3)*CONJG(HJ(3))
3487 $ -HJ(2)*CONJG(HJ(2))-HJ(1)*CONJG(HJ(1)))
3489 10 piv(i)=4.*
REAL(HN*CONJG(HJ(I)))-2.*HH*PN(I)
3492 SUBROUTINE claxi(HJ,PN,PIA)
3499 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3500 COMMON / idfc / idff
3502 COMPLEX HJ(4),HJC(4)
3505 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3510 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN 3511 sign= idff/abs(idff)
3512 ELSEIF (ktom.EQ.2)
THEN 3513 sign=-idff/abs(idff)
3515 print *,
'STOP IN CLAXI: KTOM=',ktom
3520 10 hjc(i)=conjg(hj(i))
3521 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3522 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3523 pia(3)= 2.*pn(4)*det2(1,2)
3524 pia(4)= 2.*pn(3)*det2(1,2)
3527 20 pia(i)=pia(i)*sign
3529 SUBROUTINE clnut(HJ,B,HV)
3541 b=
REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3))
& - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) 3544 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3558 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3559 * ,ampiz,ampi,amro,gamro,ama1,gama1
3560 * ,amk,amkz,amkst,gamkst
3562 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3563 * ,ampiz,ampi,amro,gamro,ama1,gama1
3564 * ,amk,amkz,amkst,gamkst
3565 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3566 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3567 COMMON /testa1/ keya1
3568 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3569 REAL PAA(4),VEC1(4),VEC2(4)
3570 REAL PIVEC(4),PIAKS(4),HVM(4)
3571 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3577 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3585 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3588 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3589 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3590 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3591 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3593 prod1 =vec1(1)*pipl(1)
3594 prod2 =vec2(2)*pipl(2)
3595 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3596 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3597 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3598 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3599 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3600 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3602 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3604 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3606 vec1(i)= vec1(i)/gnorm
3608 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3609 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3610 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3611 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3612 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3613 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3614 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3615 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3616 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3617 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3618 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3620 fnorm=formom(xmaa,xmom)
3625 hadcur(i) = fnorm *(
3626 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3627 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3628 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3630 hadcur(i) = fnorm *(
3631 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3632 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3633 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3638 CALL clvec(hadcur,pn,pivec)
3639 CALL claxi(hadcur,pn,piaks)
3640 CALL clnut(hadcur,brakm,hvm)
3642 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3643 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3645 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3646 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3650 amplit=(gfermi*ccabib)**2*brak/2.
3659 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3673 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3674 * ,ampiz,ampi,amro,gamro,ama1,gama1
3675 * ,amk,amkz,amkst,gamkst
3677 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3678 * ,ampiz,ampi,amro,gamro,ama1,gama1
3679 * ,amk,amkz,amkst,gamkst
3680 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3681 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3682 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3683 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3684 REAL PIVEC(4),PIAKS(4),HVM(4)
3685 REAL FNORM(0:7),COEF(1:5,0:7)
3686 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3688 COMPLEX F1,F2,F3,F4,F5
3690 EXTERNAL form1,form2,form3,form4,form5
3691 DATA pi /3.141592653589793238462643/
3695 IF (icont.EQ.0)
THEN 3703 fnorm(4)=scabib/fpi/dwapi0
3708 coef(1,0)= 2.0*sqrt(2.)/3.0
3709 coef(2,0)=-2.0*sqrt(2.)/3.0
3712 coef(3,0)= 2.0*sqrt(2.)/3.0
3719 coef(1,1)=-sqrt(2.)/3.0
3720 coef(2,1)= sqrt(2.)/3.0
3725 coef(1,2)=-sqrt(2.)/3.0
3726 coef(2,2)= sqrt(2.)/3.0
3744 coef(1,4)= 1.0/sqrt(2.)/3.0
3745 coef(2,4)=-1.0/sqrt(2.)/3.0
3750 coef(1,5)=-sqrt(2.)/3.0
3751 coef(2,5)= sqrt(2.)/3.0
3773 coef(5,7)=-sqrt(2.0/3.0)
3778 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3779 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3780 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3781 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3782 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3783 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3784 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3785 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3787 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3788 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3789 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3790 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3791 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3792 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3794 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3795 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3796 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3797 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3798 CALL prod5(pim1,pim2,pim3,vec5)
3803 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3804 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3805 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3807 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3808 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3809 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3812 hadcur(i)= cmplx(fnorm(mnum)) * (
3813 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3814 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3818 hadcur(i)= cmplx(fnorm(mnum)) * (
3819 $cmplx(vec1(i)*coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)+
3820 $cmplx(vec2(i)*coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)+
3821 $cmplx(vec3(i)*coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)+
3823 $cmplx(vec4(i)*coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,
3824 $ xmro2**2,xmro3**2) +
3825 $(-1.0)*uroj/4.0/pi**2/fpi**2*
3826 $cmplx(vec5(i)*coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2))
3831 CALL clvec(hadcur,pn,pivec)
3832 CALL claxi(hadcur,pn,piaks)
3833 CALL clnut(hadcur,brakm,hvm)
3835 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3836 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3837 amplit=(gfermi)**2*brak/2.
3839 print *,
'MNUM=',mnum
3845 IF (k.EQ.4) znak=1.0
3846 xm1=znak*pim1(k)**2+xm1
3847 xm2=znak*pim2(k)**2+xm2
3848 xm3=znak*pim3(k)**2+xm3
3849 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3850 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3851 print *,
'************************************************' 3855 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3856 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3861 SUBROUTINE prod5(P1,P2,P3,PIA)
3867 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3868 COMMON / idfc / idff
3869 REAL PIA(4),P1(4),P2(4),P3(4)
3870 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3872 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN 3873 sign= idff/abs(idff)
3874 ELSEIF (ktom.EQ.2)
THEN 3875 sign=-idff/abs(idff)
3877 print *,
'STOP IN PROD5: KTOM=',ktom
3883 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3884 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3885 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3886 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3889 20 pia(i)=pia(i)*sign
3892 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3909 COMMON / inout / inut,iout
3910 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3916 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3923 ELSEIF(mode.EQ. 0)
THEN 3926 IF(iwarm.EQ.0)
GOTO 902
3927 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3928 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3931 IF(rn(1).GT.wt)
GOTO 300
3933 ELSEIF(mode.EQ. 1)
THEN 3935 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3940 902
WRITE(iout, 9020)
3941 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3944 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3946 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3947 * ,ampiz,ampi,amro,gamro,ama1,gama1
3948 * ,amk,amkz,amkst,gamkst
3950 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3951 * ,ampiz,ampi,amro,gamro,ama1,gama1
3952 * ,amk,amkz,amkst,gamkst
3953 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3954 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3955 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3956 REAL*4 GAMPMC ,GAMPER
3959 COMMON / inout / inut,iout
3961 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3963 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3965 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3968 CHARACTER NAMES(NMODE)*31
3970 COMMON / inout / inut,iout
3973 REAL*4 PNU(4),PWB(4),PNPI(4,9),HV(4),HHV(4)
3974 REAL*4 PDUM1(4),PDUM2(4),PDUMI(4,9)
3977 REAL*8 SWT(NMODE),SSWT(NMODE)
3978 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3980 DATA pi /3.141592653589793238462643/
3996 #if defined (CePeCe) 3998 #elif defined (ALEPH) 4004 IF (jnpi.LE.nm4)
THEN 4005 a = pkorb(3,37+jnpi)
4017 ELSEIF(jnpi.LE.nm4)
THEN 4018 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4019 ELSEIF(jnpi.LE.nm4+nm5)
THEN 4020 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4021 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN 4022 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4023 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN 4024 inum=jnpi-nm4-nm5-nm6
4025 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
4026 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN 4027 inum=jnpi-nm4-nm5-nm6-nm3
4028 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
4032 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
4034 #if defined (CePeCe) 4035 #elif defined (ALEPH) 4044 ELSEIF(mode.EQ. 0)
THEN 4046 IF(iwarm.EQ.0)
GOTO 902
4051 ELSEIF(jnpi.LE.nm4)
THEN 4052 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4053 ELSEIF(jnpi.LE.nm4+nm5)
THEN 4054 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4055 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN 4056 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
4057 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN 4058 inum=jnpi-nm4-nm5-nm6
4059 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
4060 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN 4061 inum=jnpi-nm4-nm5-nm6-nm3
4062 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
4070 nevraw(jnpi)=nevraw(jnpi)+1
4071 swt(jnpi)=swt(jnpi)+wt
4073 sswt(jnpi)=sswt(jnpi)+wt**2
4077 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
4082 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
4083 IF(rn*wtmax(jnpi).GT.wt)
GOTO 300
4085 costhe=-1.+2.*rrr(2)
4088 CALL rotor2(thet,pnu,pnu)
4089 CALL rotor3( phi,pnu,pnu)
4090 CALL rotor2(thet,pwb,pwb)
4091 CALL rotor3( phi,pwb,pwb)
4092 CALL rotor2(thet,hv,hv)
4093 CALL rotor3( phi,hv,hv)
4096 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
4097 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
4099 nevacc(jnpi)=nevacc(jnpi)+1
4101 ELSEIF(mode.EQ. 1)
THEN 4104 IF(nevraw(jnpi).EQ.0)
GOTO 500
4105 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
4107 IF(nevraw(jnpi).NE.0)
4108 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
4110 WRITE(iout, 7010) names(jnpi),
4111 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
4113 gampmc(8+jnpi-1)=rat
4114 gamper(8+jnpi-1)=error
4120 7003
FORMAT(///1x,15(5h*****)
4121 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
4123 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
4125 $ /,1x,15(5h*****)/)
4126 7010
FORMAT(///1x,15(5h*****)
4127 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
4128 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
4129 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
4130 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
4131 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
4132 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
4133 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
4134 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
4135 $ /,1x,15(5h*****)/)
4136 902
WRITE(iout, 9020)
4137 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
4139 903
WRITE(iout, 9030) jnpi,mode
4140 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
4145 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4155 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4156 * ,ampiz,ampi,amro,gamro,ama1,gama1
4157 * ,amk,amkz,amkst,gamkst
4159 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4160 * ,ampiz,ampi,amro,gamro,ama1,gama1
4161 * ,amk,amkz,amkst,gamkst
4162 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4163 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4165 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4166 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4168 CHARACTER NAMES(NMODE)*31
4171 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
4174 REAL*8 UU,FF,FF1,FF2,FF3,FF4,GG1,GG2,GG3,GG4,RR
4175 DATA pi /3.141592653589793238462643/
4177 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
4182 phspac=1./2**23/pi**11
4186 amp1=dcdmas(idffin(1,jnpi))
4187 amp2=dcdmas(idffin(2,jnpi))
4188 amp3=dcdmas(idffin(3,jnpi))
4189 amp4=dcdmas(idffin(4,jnpi))
4193 #if defined (CePeCe) 4200 #elif defined (CLEO) 4235 CALL choice(100+jnpi,rrdum,ichan,prob1,prob2,prob3,
4238 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
4240 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4254 ams1=(amp1+amp2+amp3+amp4)**2
4255 ams2=(amtau-amnuta)**2
4256 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4257 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4258 alp=alp1+rr1*(alp2-alp1)
4259 am4sq =amrop**2+amrop*gamrop*tan(alp)
4262 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4263 phspac=phspac*(alp2-alp1)
4267 ams1=(amp2+amp3+amp4)**2
4269 IF (rrr(9).GT.prez)
THEN 4270 am3sq=ams1+ rr1*(ams2-ams1)
4276 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4277 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4278 alp=alp1+rr1*(alp2-alp1)
4279 am3sq =amrx**2+amrx*gamrx*tan(alp)
4282 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4290 am2sq=ams1+ rr2*(ams2-ams1)
4296 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4297 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4298 ppi= enq1**2-amp3**2
4299 pppi=sqrt(abs(enq1**2-amp3**2))
4301 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4302 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4303 ppi= enq1**2-amp4**2
4304 pppi=sqrt(abs(enq1**2-amp4**2))
4306 phspac=phspac*(4*pi)*(2*pppi/am2)
4312 CALL sphera(pppi,piz)
4330 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4331 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4332 ppi = pr(4)**2-am2**2
4340 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4343 ff3=(4*pi)*(2*pr(3)/am3)
4345 exe=(pr(4)+pr(3))/am2
4346 CALL bostr3(exe,piz,piz)
4347 CALL bostr3(exe,pipl,pipl)
4350 thet =acos(-1.+2*rr3)
4352 CALL rotpol(thet,phi,pipl)
4353 CALL rotpol(thet,phi,pim1)
4354 CALL rotpol(thet,phi,piz)
4355 CALL rotpol(thet,phi,pr)
4365 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4366 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4367 ppi = pr(4)**2-am3**2
4375 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4378 ff4=(4*pi)*(2*pr(3)/am4)
4380 exe=(pr(4)+pr(3))/am3
4381 CALL bostr3(exe,piz,piz)
4382 CALL bostr3(exe,pipl,pipl)
4383 CALL bostr3(exe,pim1,pim1)
4386 thet =acos(-1.+2*rr3)
4388 CALL rotpol(thet,phi,pipl)
4389 CALL rotpol(thet,phi,pim1)
4390 CALL rotpol(thet,phi,pim2)
4391 CALL rotpol(thet,phi,piz)
4392 CALL rotpol(thet,phi,pr)
4398 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4399 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4400 ppi = paa(4)**2-am4**2
4401 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4402 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4406 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4409 IF(rrr(9).LE.0.5*prez)
THEN 4418 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4419 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4421 ams1=(amp1+amp3+amp4)**2
4424 ams2=(sqrt(am3sq)-amp1)**2
4426 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4427 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4430 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4431 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4433 ams1=(amp1+amp3+amp4)**2
4434 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4435 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4436 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4439 ams2=(sqrt(am3sq)-amp1)**2
4441 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4442 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4445 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4446 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4448 ams1=(amp2+amp3+amp4)**2
4449 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4450 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4451 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4454 ams2=(sqrt(am3sq)-amp2)**2
4456 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4457 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4460 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN 4461 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4489 exe=(paa(4)+paa(3))/am4
4490 CALL bostr3(exe,piz,piz)
4491 CALL bostr3(exe,pipl,pipl)
4492 CALL bostr3(exe,pim1,pim1)
4493 CALL bostr3(exe,pim2,pim2)
4494 CALL bostr3(exe,pr,pr)
4499 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4506 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4507 ELSEIF (jnpi.EQ.2)
THEN 4508 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4511 dgamt=1/(2.*amtau)*amplit*phspac
4526 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4540 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4541 * ,ampiz,ampi,amro,gamro,ama1,gama1
4542 * ,amk,amkz,amkst,gamkst
4544 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4545 * ,ampiz,ampi,amro,gamro,ama1,gama1
4546 * ,amk,amkz,amkst,gamkst
4547 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4548 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4549 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4550 REAL PIVEC(4),PIAKS(4),HVM(4)
4551 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4552 EXTERNAL form1,form2,form3,form4,form5
4553 DATA pi /3.141592653589793238462643/
4557 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4559 CALL curr(mnum,pim1,pim2,pim3,pim4,hadcur)
4563 CALL clvec(hadcur,pn,pivec)
4564 CALL claxi(hadcur,pn,piaks)
4565 CALL clnut(hadcur,brakm,hvm)
4567 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4568 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4569 amplit=(ccabib*gfermi)**2*brak/2.
4572 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4573 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4579 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4584 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4585 * ,ampiz,ampi,amro,gamro,ama1,gama1
4586 * ,amk,amkz,amkst,gamkst
4588 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4589 * ,ampiz,ampi,amro,gamro,ama1,gama1
4592 * ,amk,amkz,amkst,gamkst
4593 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4594 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4595 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4597 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4599 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4602 CHARACTER NAMES(NMODE)*31
4603 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4604 REAL*4 PR(4),PI1(4),PI2(4),PI3(4),PI4(4),PI5(4)
4605 REAL*8 AMP1,AMP2,AMP3,AMP4,AMP5,ams1,ams2,amom,gamom
4606 REAL*8 AM5SQ,AM4SQ,AM3SQ,AM2SQ,AM5,AM4,AM3
4608 REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4617 DATA pi /3.141592653589793238462643/
4624 bwign(xm,am,gammab)=xm**2/cmplx(xm**2-am**2,gammab*am)
4626 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4635 phspac=1./2**29/pi**14
4638 amp1=dcdmas(idffin(1,jnpi))
4639 amp2=dcdmas(idffin(2,jnpi))
4640 amp3=dcdmas(idffin(3,jnpi))
4641 amp4=dcdmas(idffin(4,jnpi))
4642 amp5=dcdmas(idffin(5,jnpi))
4657 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4658 ams2=(amtau-amnuta)**2
4659 am5sq=ams1+ rr1*(ams2-ams1)
4661 phspac=phspac*(ams2-ams1)
4666 ams1=(amp2+amp3+amp4+amp5)**2
4668 am4sq=ams1+ rr1*(ams2-ams1)
4675 ams1=(amp2+amp3+amp4)**2
4677 alp1=atan((ams1-amom**2)/amom/gamom)
4678 alp2=atan((ams2-amom**2)/amom/gamom)
4679 alp=alp1+rr1*(alp2-alp1)
4680 am3sq =amom**2+amom*gamom*tan(alp)
4683 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4696 am2sq=ams1+ rr2*(ams2-ams1)
4702 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4703 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4704 ppi= enq1**2-amp3**2
4705 pppi=sqrt(abs(enq1**2-amp3**2))
4706 ff1=(4*pi)*(2*pppi/am2)
4708 call sphera(pppi,pi3)
4719 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4720 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4721 ppi = pr(4)**2-am2**2
4725 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4728 ff2=(4*pi)*(2*pr(3)/am3)
4730 exe=(pr(4)+pr(3))/am2
4731 call bostr3(exe,pi3,pi3)
4732 call bostr3(exe,pi4,pi4)
4735 thet =acos(-1.+2*rr3)
4737 call rotpol(thet,phi,pi2)
4738 call rotpol(thet,phi,pi3)
4739 call rotpol(thet,phi,pi4)
4745 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4746 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4747 ppi = pr(4)**2-am3**2
4751 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4754 ff3=(4*pi)*(2*pr(3)/am4)
4756 exe=(pr(4)+pr(3))/am3
4757 call bostr3(exe,pi2,pi2)
4758 call bostr3(exe,pi3,pi3)
4759 call bostr3(exe,pi4,pi4)
4762 thet =acos(-1.+2*rr3)
4764 call rotpol(thet,phi,pi2)
4765 call rotpol(thet,phi,pi3)
4766 call rotpol(thet,phi,pi4)
4767 call rotpol(thet,phi,pi5)
4773 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4774 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4775 ppi = pr(4)**2-am4**2
4779 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4782 ff4=(4*pi)*(2*pr(3)/am5)
4784 exe=(pr(4)+pr(3))/am4
4785 call bostr3(exe,pi2,pi2)
4786 call bostr3(exe,pi3,pi3)
4787 call bostr3(exe,pi4,pi4)
4788 call bostr3(exe,pi5,pi5)
4791 thet =acos(-1.+2*rr3)
4793 call rotpol(thet,phi,pi1)
4794 call rotpol(thet,phi,pi2)
4795 call rotpol(thet,phi,pi3)
4796 call rotpol(thet,phi,pi4)
4797 call rotpol(thet,phi,pi5)
4806 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4807 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4808 ppi = paa(4)**2-am5sq
4809 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4813 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4816 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4820 exe=(paa(4)+paa(3))/am5
4821 call bostr3(exe,pi1,pi1)
4822 call bostr3(exe,pi2,pi2)
4823 call bostr3(exe,pi3,pi3)
4824 call bostr3(exe,pi4,pi4)
4825 call bostr3(exe,pi5,pi5)
4833 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4834 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4835 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4836 fompp = cabs(bwign(am3,amom,gamom))**2
4841 amplit=ccabib**2*gfermi**2/2. * brak
4842 amplit = amplit * fompp * fnorm
4845 dgamt=1/(2.*amtau)*amplit*phspac
4868 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4874 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4875 * ,ampiz,ampi,amro,gamro,ama1,gama1
4876 * ,amk,amkz,amkst,gamkst
4878 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4879 * ,ampiz,ampi,amro,gamro,ama1,gama1
4880 * ,amk,amkz,amkst,gamkst
4881 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4882 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4883 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4889 DATA pi /3.141592653589793238462643/
4893 phspac=1./2**11/pi**5
4901 ams2=(amtau-amnuta)**2
4904 amx2=ams1+ rr1(1)*(ams2-ams1)
4924 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4925 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4929 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4931 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4934 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4935 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4936 pppi=sqrt((enq1-amk)*(enq1+amk))
4937 phspac=phspac*(4*pi)*(2*pppi/amx)
4939 CALL sphera(pppi,pkc)
4945 exe=(pr(4)+pr(3))/amx
4947 CALL bostr3(exe,pkc,pkc)
4948 CALL bostr3(exe,pkz,pkz)
4950 30 qq(i)=pkc(i)-pkz(i)
4952 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4953 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4955 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4958 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4960 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4961 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4962 & +(gv**2-ga**2)*amtau*amnuta*qq2
4963 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4964 dgamt=1/(2.*amtau)*amplit*phspac
4966 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4977 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4978 * ,ampiz,ampi,amro,gamro,ama1,gama1
4979 * ,amk,amkz,amkst,gamkst
4981 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4982 * ,ampiz,ampi,amro,gamro,ama1,gama1
4983 * ,amk,amkz,amkst,gamkst
4986 fpirk=cabs(fpikm(w,amk,amkz))**2
4989 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4994 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4999 IF (init.EQ.0 )
THEN 5012 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
5021 parameter(nmxhep=2000)
5022 common/hepevtx/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
5023 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
5028 SUBROUTINE dwrph(KTO,PHX)
5032 IMPLICIT REAL*8 (a-h,o-z)
5043 IF (qhot(4).GT.1.e-5)
CALL dwluph(kto,qhot)
5046 SUBROUTINE dwluph(KTO,PHOT)
5057 COMMON /taupos/ np1,np2
5063 COMMON /taupos/ np1,np2
5067 IF (phot(4).LE.0.0)
RETURN 5070 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN 5077 IF(ktos.GT.10) ktos=ktos-10
5079 CALL tralo4(ktos,phot,phot,am)
5080 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
5085 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
5096 COMMON /taupos/ np1,np2
5099 REAL PNU(4),PWB(4),PEL(4),PNE(4)
5102 COMMON /taupos/ np1,np2
5113 CALL tralo4(kto,pnu,pnu,am)
5114 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5117 CALL tralo4(kto,pwb,pwb,am)
5121 CALL tralo4(kto,pel,pel,am)
5122 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
5125 CALL tralo4(kto,pne,pne,am)
5126 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
5130 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
5141 COMMON /taupos/ np1,np2
5144 REAL PNU(4),PWB(4),PMU(4),PNM(4)
5147 COMMON /taupos/ np1,np2
5158 CALL tralo4(kto,pnu,pnu,am)
5159 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5162 CALL tralo4(kto,pwb,pwb,am)
5166 CALL tralo4(kto,pmu,pmu,am)
5167 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5170 CALL tralo4(kto,pnm,pnm,am)
5171 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5175 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5186 COMMON /taupos/ np1,np2
5196 CALL tralo4(kto,pnu,pnu,am)
5197 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5200 CALL tralo4(kto,ppi,ppi,am)
5201 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5205 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5216 COMMON /taupos/ np1,np2
5219 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5222 COMMON /taupos/ np1,np2
5233 CALL tralo4(kto,pnu,pnu,am)
5234 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5237 CALL tralo4(kto,prho,prho,am)
5238 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5241 CALL tralo4(kto,pic,pic,am)
5242 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5245 CALL tralo4(kto,piz,piz,am)
5246 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5250 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5262 COMMON /taupos/ np1,np2
5265 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5268 COMMON /taupos/ np1,np2
5279 CALL tralo4(kto,pnu,pnu,am)
5280 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5283 CALL tralo4(kto,paa,paa,am)
5284 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5292 CALL tralo4(kto,pim2,pim2,am)
5293 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5296 CALL tralo4(kto,pim1,pim1,am)
5297 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5300 CALL tralo4(kto,pipl,pipl,am)
5301 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5303 ELSE IF (jaa.EQ.2)
THEN 5308 CALL tralo4(kto,pim2,pim2,am)
5309 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5312 CALL tralo4(kto,pim1,pim1,am)
5313 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5316 CALL tralo4(kto,pipl,pipl,am)
5317 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5323 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5333 COMMON /taupos/ np1,np2
5347 CALL tralo4 (kto,pnu,pnu,am)
5348 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5351 CALL tralo4 (kto,pkk,pkk,am)
5352 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5356 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5357 COMMON / taukle / bra1,brk0,brk0b,brks
5358 REAL*4 BRA1,BRK0,BRK0B,BRKS
5360 COMMON /taupos/ np1,np2
5373 REAL PNU(4),PKS(4),PKK(4),PPI(4)
5375 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5376 COMMON /taupos/ np1,np2
5387 CALL tralo4(kto,pnu,pnu,am)
5388 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5391 CALL tralo4(kto,pks,pks,am)
5392 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5400 CALL tralo4(kto,ppi,ppi,am)
5401 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5404 IF (isgn.EQ.-1) bran=brk0
5407 IF(xio(1).GT.bran)
THEN 5413 CALL tralo4(kto,pkk,pkk,am)
5414 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5416 ELSE IF(jkst.EQ.20)
THEN 5421 CALL tralo4(kto,ppi,ppi,am)
5422 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5425 CALL tralo4(kto,pkk,pkk,am)
5426 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5432 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5442 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5444 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5446 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5449 COMMON /taupos/ np1,np2
5450 CHARACTER NAMES(NMODE)*31
5451 REAL PNU(4),PWB(4),PNPI(4,9)
5463 CALL tralo4(kto,pnu,pnu,am)
5464 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5467 CALL tralo4(kto,pwb,pwb,am)
5468 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5477 kfpi=lunpik(idffin(i,jnpi), isgn)
5479 kfpi=lunpik(idffin(i,jnpi),-isgn)
5486 CALL tralo4(kto,ppi,ppi,am)
5487 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5492 #if defined (CePeCe) 5501 IMPLICIT REAL*8 (a-h,o-z)
5503 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5505 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5517 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5518 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5527 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5528 DATA pi /3.141592653589793238462643d0/
5530 IF(abs(y).LT.abs(x))
THEN 5532 IF(x.LE.0d0) the=pi-the
5534 the=acos(x/sqrt(x**2+y**2))
5545 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5546 DATA pi /3.141592653589793238462643d0/
5548 IF(abs(y).LT.abs(x))
THEN 5550 IF(x.LE.0d0) the=pi-the
5552 the=acos(x/sqrt(x**2+y**2))
5554 IF(y.LT.0d0) the=2d0*pi-the
5557 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5562 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5563 dimension pvec(4),qvec(4),rvec(4)
5571 qvec(2)= cs*rvec(2)-sn*rvec(3)
5572 qvec(3)= sn*rvec(2)+cs*rvec(3)
5576 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5581 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5582 dimension pvec(4),qvec(4),rvec(4)
5589 qvec(1)= cs*rvec(1)+sn*rvec(3)
5591 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5595 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5600 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5602 dimension pvec(4),qvec(4),rvec(4)
5608 qvec(1)= cs*rvec(1)-sn*rvec(2)
5609 qvec(2)= sn*rvec(1)+cs*rvec(2)
5613 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5619 REAL*4 PVEC(4),QVEC(4),RVEC(4)
5632 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5638 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5639 dimension pvec(4),qvec(4),rvec(4)
5653 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5658 REAL*4 PVEC(4),QVEC(4),RVEC(4)
5666 qvec(2)= cs*rvec(2)-sn*rvec(3)
5667 qvec(3)= sn*rvec(2)+cs*rvec(3)
5670 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5675 IMPLICIT REAL*4(a-h,o-z)
5676 REAL*4 PVEC(4),QVEC(4),RVEC(4)
5683 qvec(1)= cs*rvec(1)+sn*rvec(3)
5685 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5688 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5693 REAL*4 PVEC(4),QVEC(4),RVEC(4)
5699 qvec(1)= cs*rvec(1)-sn*rvec(2)
5700 qvec(2)= sn*rvec(1)+cs*rvec(2)
5704 SUBROUTINE spherd(R,X)
5709 REAL*8 R,X(4),PI,COSTH,SINTH
5711 DATA pi /3.141592653589793238462643d0/
5715 sinth=sqrt(1 -costh**2)
5716 x(1)=r*sinth*cos(2*pi*rrr(2))
5717 x(2)=r*sinth*sin(2*pi*rrr(2))
5721 SUBROUTINE rotpox(THET,PHI,PP)
5722 IMPLICIT REAL*8 (a-h,o-z)
5732 CALL rotod2(thet,pp,pp)
5733 CALL rotod3( phi,pp,pp)
5736 SUBROUTINE sphera(R,X)
5744 DATA pi /3.141592653589793238462643/
5748 sinth=sqrt(1.-costh**2)
5749 x(1)=r*sinth*cos(2*pi*rrr(2))
5750 x(2)=r*sinth*sin(2*pi*rrr(2))
5754 SUBROUTINE rotpol(THET,PHI,PP)
5761 CALL rotor2(thet,pp,pp)
5762 CALL rotor3( phi,pp,pp)
5765 #include "../randg/tauola-random.h" 5768 IMPLICIT REAL*8(a-h,o-z)
5771 IF(x .LT.-1.0)
GO TO 1
5772 IF(x .LE. 0.5)
GO TO 2
5773 IF(x .EQ. 1.0)
GO TO 3
5774 IF(x .LE. 2.0)
GO TO 4
5778 z=z-0.5* log(abs(x))**2
5784 3 dilogt=1.64493406684822
5788 z=1.64493406684822 - log(x)* log(abs(t))
5789 5 y=2.66666666666666 *t+0.66666666666666
5790 b= 0.00000 00000 00001
5791 a=y*b +0.00000 00000 00004
5792 b=y*a-b+0.00000 00000 00011
5793 a=y*b-a+0.00000 00000 00037
5794 b=y*a-b+0.00000 00000 00121
5795 a=y*b-a+0.00000 00000 00398
5796 b=y*a-b+0.00000 00000 01312
5797 a=y*b-a+0.00000 00000 04342
5798 b=y*a-b+0.00000 00000 14437
5799 a=y*b-a+0.00000 00000 48274
5800 b=y*a-b+0.00000 00001 62421
5801 a=y*b-a+0.00000 00005 50291
5802 b=y*a-b+0.00000 00018 79117
5803 a=y*b-a+0.00000 00064 74338
5804 b=y*a-b+0.00000 00225 36705
5805 a=y*b-a+0.00000 00793 87055
5806 b=y*a-b+0.00000 02835 75385
5807 a=y*b-a+0.00000 10299 04264
5808 b=y*a-b+0.00000 38163 29463
5809 a=y*b-a+0.00001 44963 00557
5810 b=y*a-b+0.00005 68178 22718
5811 a=y*b-a+0.00023 20021 96094
5812 b=y*a-b+0.00100 16274 96164
5813 a=y*b-a+0.00468 63619 59447
5814 b=y*a-b+0.02487 93229 24228
5815 a=y*b-a+0.16607 30329 27855
5816 a=y*a-b+1.93506 43008 6996