00001
00002 SUBROUTINE JAKER(JAK)
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 COMMON / TAUBRA / GAMPRT(500),JLIST(500),NCHAN
00037
00038
00039
00040 REAL CUMUL(500),RRR(1)
00041
00042 IF(NCHAN.LE.0.OR.NCHAN.GT.500) GOTO 902
00043 CALL RANMAR(RRR,1)
00044 SUM=0
00045 DO 20 I=1,NCHAN
00046 SUM=SUM+GAMPRT(I)
00047 20 CUMUL(I)=SUM
00048 DO 25 I=NCHAN,1,-1
00049 IF(RRR(1).LT.CUMUL(I)/CUMUL(NCHAN)) JI=I
00050 25 CONTINUE
00051 JAK=JLIST(JI)
00052 RETURN
00053 902 PRINT 9020
00054 9020 FORMAT(' ----- JAKER: WRONG NCHAN')
00055 STOP
00056 END
00057 SUBROUTINE DEKAY(KTO,HX)
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 REAL H(4)
00083 REAL*8 HX(4)
00084 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
00085
00086 COMMON / IDFC / IDF
00087
00088 COMMON /TAUPOS/ NP1,NP2
00089 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00090 REAL*4 GAMPMC ,GAMPER
00091 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
00092
00093 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
00094
00095 & ,NAMES
00096 CHARACTER NAMES(NMODE)*31
00097 COMMON / INOUT / INUT,IOUT
00098 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4)
00099 REAL PDUMX(4,9)
00100 DATA IWARM/0/
00101 KTOM=KTO
00102
00103 IF(KTO.EQ.-1) THEN
00104
00105
00106
00107 NP1=3
00108 NP2=4
00109 KTOM=1
00110 IF (IWARM.EQ.1) X=5/(IWARM-1)
00111 IWARM=1
00112 WRITE(IOUT,7001) JAK1,JAK2
00113 NEVTOT=0
00114 NEV1=0
00115 NEV2=0
00116 IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN
00117 CALL DADMEL(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00118 CALL DADMMU(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00119 CALL DADMPI(-1,IDUM,PDUM,PDUM1,PDUM2)
00120 CALL DADMRO(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4)
00121 CALL DADMAA(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JDUM)
00122 CALL DADMKK(-1,IDUM,PDUM,PDUM1,PDUM2)
00123 CALL DADMKS(-1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,JDUM)
00124 CALL DADNEW(-1,IDUM,HDUM,PDUM1,PDUM2,PDUMX,JDUM)
00125 ENDIF
00126 DO 21 I=1,500
00127 NEVDEC(I)=0
00128 GAMPMC(I)=0
00129 21 GAMPER(I)=0
00130 ELSEIF(KTO.EQ.1) THEN
00131
00132
00133 NEVTOT=NEVTOT+1
00134 IF(IWARM.EQ.0) GOTO 902
00135 ISGN= IDF/IABS(IDF)
00136
00137
00138 CALL TAURDF(KTO)
00139
00140 CALL DEKAY1(0,H,ISGN)
00141 ELSEIF(KTO.EQ.2) THEN
00142
00143
00144 NEVTOT=NEVTOT+1
00145 IF(IWARM.EQ.0) GOTO 902
00146 ISGN=-IDF/IABS(IDF)
00147
00148
00149 CALL TAURDF(KTO)
00150
00151 CALL DEKAY2(0,H,ISGN)
00152 ELSEIF(KTO.EQ.11) THEN
00153
00154
00155 NEV1=NEV1+1
00156 ISGN= IDF/IABS(IDF)
00157 CALL DEKAY1(1,H,ISGN)
00158 ELSEIF(KTO.EQ.12) THEN
00159
00160
00161 NEV2=NEV2+1
00162 ISGN=-IDF/IABS(IDF)
00163 CALL DEKAY2(1,H,ISGN)
00164 ELSEIF(KTO.EQ.100) THEN
00165
00166 IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN
00167 CALL DADMEL( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00168 CALL DADMMU( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00169 CALL DADMPI( 1,IDUM,PDUM,PDUM1,PDUM2)
00170 CALL DADMRO( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4)
00171 CALL DADMAA( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JDUM)
00172 CALL DADMKK( 1,IDUM,PDUM,PDUM1,PDUM2)
00173 CALL DADMKS( 1,IDUM,HDUM,PDUM1,PDUM2,PDUM3,PDUM4,JDUM)
00174 CALL DADNEW( 1,IDUM,HDUM,PDUM1,PDUM2,PDUMX,JDUM)
00175 WRITE(IOUT,7010) NEV1,NEV2,NEVTOT
00176 WRITE(IOUT,7011) (NEVDEC(I),GAMPMC(I),GAMPER(I),I= 1,7)
00177 WRITE(IOUT,7012)
00178 $ (NEVDEC(I),GAMPMC(I),GAMPER(I),NAMES(I-7),I=8,7+NMODE)
00179 WRITE(IOUT,7013)
00180 ENDIF
00181 ELSE
00182
00183 GOTO 910
00184 ENDIF
00185
00186 DO 78 K=1,4
00187 78 HX(K)=H(K)
00188 RETURN
00189 7001 FORMAT(///1X,15(5H*****)
00190
00191 $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9X,1H*,
00192 $ /,' *', 25X,'***********Sep 2005***************',9X,1H*,
00193 $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*,
00194 $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*,
00195 $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*,
00196 $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*,
00197 $ /,' *', 25X,' Physics initialization by CLEO collab ',9X,1H*,
00198 $ /,' *', 25X,' see Alain Weinstein www home page: ',9X,1H*,
00199 $ /,' *', 25X,'http://www.cithep.caltech.edu/~ajw/ ',9X,1H*,
00200 $ /,' *', 25X,'/korb_doc.html#files ',9X,1H*,
00201 $ /,' *', 25X,'*******CERN TH-6793 NOVEMBER 1992*****',9X,1H*,
00202 $ /,' *', 25X,'**5 or more pi dec.: precision limited ',9X,1H*,
00203
00204 $ /,' *', 25X,'****DEKAY ROUTINE: INITIALIZATION******',9X,1H*,
00205 $ /,' *',I20 ,5X,'JAK1 = DECAY MODE TAU+ ',9X,1H*,
00206 $ /,' *',I20 ,5X,'JAK2 = DECAY MODE TAU- ',9X,1H*,
00207 $ /,1X,15(5H*****)/)
00208 7010 FORMAT(///1X,15(5H*****)
00209
00210 $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9X,1H*,
00211 $ /,' *', 25X,'***********Sep 2005***************',9X,1H*,
00212
00213 $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*,
00214 $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*,
00215 $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*,
00216 $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*,
00217 $ /,' *', 25X,'*******CERN-TH-5856 SEPTEMBER 1990*****',9X,1H*,
00218 $ /,' *', 25X,'*******CERN-TH-6195 SEPTEMBER 1991*****',9X,1H*,
00219 $ /,' *', 25X,'*******CERN TH-6793 NOVEMBER 1992*****',9X,1H*,
00220 $ /,' *', 25X,'*****DEKAY ROUTINE: FINAL REPORT*******',9X,1H*,
00221 $ /,' *',I20 ,5X,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9X,1H*,
00222 $ /,' *',I20 ,5X,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9X,1H*,
00223 $ /,' *',I20 ,5X,'NEVTOT = SUM ',9X,1H*,
00224 $ /,' *',' NOEVTS ',
00225 $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9X,1H*)
00226 7011 FORMAT(1X,'*'
00227 $ ,I10,2F12.7 ,' DADMEL ELECTRON ',9X,1H*
00228 $ /,' *',I10,2F12.7 ,' DADMMU MUON ',9X,1H*
00229 $ /,' *',I10,2F12.7 ,' DADMPI PION ',9X,1H*
00230 $ /,' *',I10,2F12.7, ' DADMRO RHO (->2PI) ',9X,1H*
00231 $ /,' *',I10,2F12.7, ' DADMAA A1 (->3PI) ',9X,1H*
00232 $ /,' *',I10,2F12.7, ' DADMKK KAON ',9X,1H*
00233 $ /,' *',I10,2F12.7, ' DADMKS K* ',9X,1H*)
00234 7012 FORMAT(1X,'*'
00235 $ ,I10,2F12.7,A31 ,8X,1H*)
00236 7013 FORMAT(1X,'*'
00237 $ ,20X,'THE ERROR IS RELATIVE AND PART.WIDTH ',10X,1H*
00238 $ /,' *',20X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10X,1H*
00239 $ /,1X,15(5H*****)/)
00240 902 PRINT 9020
00241 9020 FORMAT(' ----- DEKAY: LACK OF INITIALISATION')
00242 STOP
00243 910 PRINT 9100
00244 9100 FORMAT(' ----- DEKAY: WRONG VALUE OF KTO ')
00245 STOP
00246 END
00247 SUBROUTINE DEKAY1(IMOD,HH,ISGN)
00248
00249
00250
00251 COMMON / DECP4 / PP1(4),PP2(4),KF1,KF2
00252 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
00253 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00254 REAL*4 GAMPMC ,GAMPER
00255
00256 REAL HH(4)
00257 REAL HV(4),PNU(4),PPI(4)
00258 REAL PWB(4),PMU(4),PNM(4)
00259 REAL PRHO(4),PIC(4),PIZ(4)
00260 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
00261 REAL PKK(4),PKS(4)
00262 REAL PNPI(4,9)
00263 REAL PHOT(4)
00264 REAL PDUM(4)
00265 DATA NEV,NPRIN/0,10/
00266 KTO=1
00267 IF(JAK1.EQ.-1) RETURN
00268 IMD=IMOD
00269 IF(IMD.EQ.0) THEN
00270
00271 JAK=JAK1
00272 IF(JAK1.EQ.0) CALL JAKER(JAK)
00273 IF(JAK.EQ.1) THEN
00274 CALL DADMEL(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT)
00275 ELSEIF(JAK.EQ.2) THEN
00276 CALL DADMMU(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT)
00277 ELSEIF(JAK.EQ.3) THEN
00278 CALL DADMPI(0, ISGN,HV,PPI,PNU)
00279 ELSEIF(JAK.EQ.4) THEN
00280 CALL DADMRO(0, ISGN,HV,PNU,PRHO,PIC,PIZ)
00281 ELSEIF(JAK.EQ.5) THEN
00282 CALL DADMAA(0, ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
00283 ELSEIF(JAK.EQ.6) THEN
00284 CALL DADMKK(0, ISGN,HV,PKK,PNU)
00285 ELSEIF(JAK.EQ.7) THEN
00286 CALL DADMKS(0, ISGN,HV,PNU,PKS ,PKK,PPI,JKST)
00287 ELSE
00288 CALL DADNEW(0, ISGN,HV,PNU,PWB,PNPI,JAK-7)
00289 ENDIF
00290 DO 33 I=1,3
00291 33 HH(I)=HV(I)
00292 HH(4)=1.0
00293
00294 ELSEIF(IMD.EQ.1) THEN
00295
00296 NEV=NEV+1
00297 IF (JAK.LT.501) THEN
00298 NEVDEC(JAK)=NEVDEC(JAK)+1
00299 ENDIF
00300 DO 34 I=1,4
00301 34 PDUM(I)=.0
00302 IF(JAK.EQ.1) THEN
00303 CALL DWLUEL(1,ISGN,PNU,PWB,PMU,PNM)
00304 CALL DWRPH(KTOM,PHOT)
00305 DO 10 I=1,4
00306 10 PP1(I)=PMU(I)
00307
00308 ELSEIF(JAK.EQ.2) THEN
00309 CALL DWLUMU(1,ISGN,PNU,PWB,PMU,PNM)
00310 CALL DWRPH(KTOM,PHOT)
00311 DO 20 I=1,4
00312 20 PP1(I)=PMU(I)
00313
00314 ELSEIF(JAK.EQ.3) THEN
00315 CALL DWLUPI(1,ISGN,PPI,PNU)
00316 DO 30 I=1,4
00317 30 PP1(I)=PPI(I)
00318
00319 ELSEIF(JAK.EQ.4) THEN
00320 CALL DWLURO(1,ISGN,PNU,PRHO,PIC,PIZ)
00321 DO 40 I=1,4
00322 40 PP1(I)=PRHO(I)
00323
00324 ELSEIF(JAK.EQ.5) THEN
00325 CALL DWLUAA(1,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
00326 DO 50 I=1,4
00327 50 PP1(I)=PAA(I)
00328 ELSEIF(JAK.EQ.6) THEN
00329 CALL DWLUKK(1,ISGN,PKK,PNU)
00330 DO 60 I=1,4
00331 60 PP1(I)=PKK(I)
00332 ELSEIF(JAK.EQ.7) THEN
00333 CALL DWLUKS(1,ISGN,PNU,PKS,PKK,PPI,JKST)
00334 DO 70 I=1,4
00335 70 PP1(I)=PKS(I)
00336 ELSE
00337
00338 CALL DWLNEW(1,ISGN,PNU,PWB,PNPI,JAK)
00339 DO 80 I=1,4
00340 80 PP1(I)=PWB(I)
00341 ENDIF
00342
00343 ENDIF
00344
00345 END
00346 SUBROUTINE DEKAY2(IMOD,HH,ISGN)
00347
00348
00349
00350 COMMON / DECP4 / PP1(4),PP2(4),KF1,KF2
00351 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
00352 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00353 REAL*4 GAMPMC ,GAMPER
00354
00355 REAL HH(4)
00356 REAL HV(4),PNU(4),PPI(4)
00357 REAL PWB(4),PMU(4),PNM(4)
00358 REAL PRHO(4),PIC(4),PIZ(4)
00359 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
00360 REAL PKK(4),PKS(4)
00361 REAL PNPI(4,9)
00362 REAL PHOT(4)
00363 REAL PDUM(4)
00364 DATA NEV,NPRIN/0,10/
00365 KTO=2
00366 IF(JAK2.EQ.-1) RETURN
00367 IMD=IMOD
00368 IF(IMD.EQ.0) THEN
00369
00370 JAK=JAK2
00371 IF(JAK2.EQ.0) CALL JAKER(JAK)
00372 IF(JAK.EQ.1) THEN
00373 CALL DADMEL(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT)
00374 ELSEIF(JAK.EQ.2) THEN
00375 CALL DADMMU(0, ISGN,HV,PNU,PWB,PMU,PNM,PHOT)
00376 ELSEIF(JAK.EQ.3) THEN
00377 CALL DADMPI(0, ISGN,HV,PPI,PNU)
00378 ELSEIF(JAK.EQ.4) THEN
00379 CALL DADMRO(0, ISGN,HV,PNU,PRHO,PIC,PIZ)
00380 ELSEIF(JAK.EQ.5) THEN
00381 CALL DADMAA(0, ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
00382 ELSEIF(JAK.EQ.6) THEN
00383 CALL DADMKK(0, ISGN,HV,PKK,PNU)
00384 ELSEIF(JAK.EQ.7) THEN
00385 CALL DADMKS(0, ISGN,HV,PNU,PKS ,PKK,PPI,JKST)
00386 ELSE
00387 CALL DADNEW(0, ISGN,HV,PNU,PWB,PNPI,JAK-7)
00388 ENDIF
00389 DO 33 I=1,3
00390 33 HH(I)=HV(I)
00391 HH(4)=1.0
00392 ELSEIF(IMD.EQ.1) THEN
00393
00394 NEV=NEV+1
00395 IF (JAK.LT.501) THEN
00396 NEVDEC(JAK)=NEVDEC(JAK)+1
00397 ENDIF
00398 DO 34 I=1,4
00399 34 PDUM(I)=.0
00400 IF(JAK.EQ.1) THEN
00401 CALL DWLUEL(2,ISGN,PNU,PWB,PMU,PNM)
00402 CALL DWRPH(KTOM,PHOT)
00403 DO 10 I=1,4
00404 10 PP2(I)=PMU(I)
00405
00406 ELSEIF(JAK.EQ.2) THEN
00407 CALL DWLUMU(2,ISGN,PNU,PWB,PMU,PNM)
00408 CALL DWRPH(KTOM,PHOT)
00409 DO 20 I=1,4
00410 20 PP2(I)=PMU(I)
00411
00412 ELSEIF(JAK.EQ.3) THEN
00413 CALL DWLUPI(2,ISGN,PPI,PNU)
00414 DO 30 I=1,4
00415 30 PP2(I)=PPI(I)
00416
00417 ELSEIF(JAK.EQ.4) THEN
00418 CALL DWLURO(2,ISGN,PNU,PRHO,PIC,PIZ)
00419 DO 40 I=1,4
00420 40 PP2(I)=PRHO(I)
00421
00422 ELSEIF(JAK.EQ.5) THEN
00423 CALL DWLUAA(2,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
00424 DO 50 I=1,4
00425 50 PP2(I)=PAA(I)
00426 ELSEIF(JAK.EQ.6) THEN
00427 CALL DWLUKK(2,ISGN,PKK,PNU)
00428 DO 60 I=1,4
00429 60 PP1(I)=PKK(I)
00430 ELSEIF(JAK.EQ.7) THEN
00431 CALL DWLUKS(2,ISGN,PNU,PKS,PKK,PPI,JKST)
00432 DO 70 I=1,4
00433 70 PP1(I)=PKS(I)
00434 ELSE
00435
00436 CALL DWLNEW(2,ISGN,PNU,PWB,PNPI,JAK)
00437 DO 80 I=1,4
00438 80 PP1(I)=PWB(I)
00439 ENDIF
00440
00441 ENDIF
00442
00443 END
00444 SUBROUTINE DEXAY(KTO,POL)
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00459 REAL*4 GAMPMC ,GAMPER
00460 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
00461 COMMON / IDFC / IDFF
00462 COMMON /TAUPOS/ NP1,NP2
00463 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
00464
00465 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
00466
00467 & ,NAMES
00468 CHARACTER NAMES(NMODE)*31
00469 COMMON / INOUT / INUT,IOUT
00470 REAL POL(4)
00471 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
00472 REAL PDUM(4)
00473 REAL PDUMI(4,9)
00474 DATA IWARM/0/
00475 KTOM=KTO
00476
00477 IF(KTO.EQ.-1) THEN
00478
00479
00480
00481
00482 NP1=3
00483 NP2=4
00484 IWARM=1
00485 WRITE(IOUT, 7001) JAK1,JAK2
00486 NEVTOT=0
00487 NEV1=0
00488 NEV2=0
00489 IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN
00490 CALL DEXEL(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00491 CALL DEXMU(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00492 CALL DEXPI(-1,IDUM,PDUM,PDUM1,PDUM2)
00493 CALL DEXRO(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4)
00494 CALL DEXAA(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,IDUM)
00495 CALL DEXKK(-1,IDUM,PDUM,PDUM1,PDUM2)
00496 CALL DEXKS(-1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,IDUM)
00497 CALL DEXNEW(-1,IDUM,PDUM,PDUM1,PDUM2,PDUMI,IDUM)
00498 ENDIF
00499 DO 21 I=1,500
00500 NEVDEC(I)=0
00501 GAMPMC(I)=0
00502 21 GAMPER(I)=0
00503 ELSEIF(KTO.EQ.1) THEN
00504
00505
00506 NEVTOT=NEVTOT+1
00507 NEV1=NEV1+1
00508 IF(IWARM.EQ.0) GOTO 902
00509 ISGN=IDFF/IABS(IDFF)
00510
00511 CALL DEXAY1(KTO,JAK1,JAKP,POL,ISGN)
00512 ELSEIF(KTO.EQ.2) THEN
00513
00514
00515 NEVTOT=NEVTOT+1
00516 NEV2=NEV2+1
00517 IF(IWARM.EQ.0) GOTO 902
00518 ISGN=-IDFF/IABS(IDFF)
00519
00520 CALL DEXAY1(KTO,JAK2,JAKM,POL,ISGN)
00521 ELSEIF(KTO.EQ.100) THEN
00522
00523 IF(JAK1.NE.-1.OR.JAK2.NE.-1) THEN
00524 CALL DEXEL( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00525 CALL DEXMU( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00526 CALL DEXPI( 1,IDUM,PDUM,PDUM1,PDUM2)
00527 CALL DEXRO( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4)
00528 CALL DEXAA( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,IDUM)
00529 CALL DEXKK( 1,IDUM,PDUM,PDUM1,PDUM2)
00530 CALL DEXKS( 1,IDUM,PDUM,PDUM1,PDUM2,PDUM3,PDUM4,IDUM)
00531 CALL DEXNEW( 1,IDUM,PDUM,PDUM1,PDUM2,PDUMI,IDUM)
00532 WRITE(IOUT,7010) NEV1,NEV2,NEVTOT
00533 WRITE(IOUT,7011) (NEVDEC(I),GAMPMC(I),GAMPER(I),I= 1,7)
00534 WRITE(IOUT,7012)
00535 $ (NEVDEC(I),GAMPMC(I),GAMPER(I),NAMES(I-7),I=8,7+NMODE)
00536 WRITE(IOUT,7013)
00537 ENDIF
00538 ELSE
00539 GOTO 910
00540 ENDIF
00541 RETURN
00542 7001 FORMAT(///1X,15(5H*****)
00543
00544 $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9X,1H*,
00545 $ /,' *', 25X,'*********** Sep 2005***************',9X,1H*,
00546 $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*,
00547 $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*,
00548 $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*,
00549 $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*,
00550 $ /,' *', 25X,' Physics initialization by CLEO collab ',9X,1H*,
00551 $ /,' *', 25X,' see Alain Weinstein www home page: ',9X,1H*,
00552 $ /,' *', 25X,'http://www.cithep.caltech.edu/~ajw/ ',9X,1H*,
00553 $ /,' *', 25X,'/korb_doc.html#files ',9X,1H*,
00554 $ /,' *', 25X,'*******CERN TH-6793 NOVEMBER 1992*****',9X,1H*,
00555 $ /,' *', 25X,'**5 or more pi dec.: precision limited ',9X,1H*,
00556
00557 $ /,' *', 25X,'*******CERN-TH-6793 NOVEMBER 1992*****',9X,1H*,
00558 $ /,' *', 25X,'**5 or more pi dec.: precision limited ',9X,1H*,
00559 $ /,' *', 25X,'******DEXAY ROUTINE: INITIALIZATION****',9X,1H*
00560 $ /,' *',I20 ,5X,'JAK1 = DECAY MODE FERMION1 (TAU+) ',9X,1H*
00561 $ /,' *',I20 ,5X,'JAK2 = DECAY MODE FERMION2 (TAU-) ',9X,1H*
00562 $ /,1X,15(5H*****)/)
00563
00564
00565 7010 FORMAT(///1X,15(5H*****)
00566
00567 $ /,' *', 25X,'*****TAUOLA LIBRARY: VERSION 2.9 ******',9X,1H*,
00568 $ /,' *', 25X,'***********Sep 2005***************',9X,1H*,
00569 $ /,' *', 25X,'**AUTHORS: S.JADACH, Z.WAS*************',9X,1H*,
00570 $ /,' *', 25X,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9X,1H*,
00571 $ /,' *', 25X,'**AVAILABLE FROM: WASM AT CERNVM ******',9X,1H*,
00572 $ /,' *', 25X,'***** PUBLISHED IN COMP. PHYS. COMM.***',9X,1H*,
00573 $ /,' *', 25X,'*******CERN-TH-5856 SEPTEMBER 1990*****',9X,1H*,
00574 $ /,' *', 25X,'*******CERN-TH-6195 SEPTEMBER 1991*****',9X,1H*,
00575 $ /,' *', 25X,'*******CERN-TH-6793 NOVEMBER 1992*****',9X,1H*,
00576 $ /,' *', 25X,'******DEXAY ROUTINE: FINAL REPORT******',9X,1H*
00577 $ /,' *',I20 ,5X,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9X,1H*
00578 $ /,' *',I20 ,5X,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9X,1H*
00579 $ /,' *',I20 ,5X,'NEVTOT = SUM ',9X,1H*
00580 $ /,' *',' NOEVTS ',
00581 $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9X,1H*)
00582 7011 FORMAT(1X,'*'
00583 $ ,I10,2F12.7 ,' DADMEL ELECTRON ',9X,1H*
00584 $ /,' *',I10,2F12.7 ,' DADMMU MUON ',9X,1H*
00585 $ /,' *',I10,2F12.7 ,' DADMPI PION ',9X,1H*
00586 $ /,' *',I10,2F12.7, ' DADMRO RHO (->2PI) ',9X,1H*
00587 $ /,' *',I10,2F12.7, ' DADMAA A1 (->3PI) ',9X,1H*
00588 $ /,' *',I10,2F12.7, ' DADMKK KAON ',9X,1H*
00589 $ /,' *',I10,2F12.7, ' DADMKS K* ',9X,1H*)
00590 7012 FORMAT(1X,'*'
00591 $ ,I10,2F12.7,A31 ,8X,1H*)
00592 7013 FORMAT(1X,'*'
00593 $ ,20X,'THE ERROR IS RELATIVE AND PART.WIDTH ',10X,1H*
00594 $ /,' *',20X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10X,1H*
00595 $ /,1X,15(5H*****)/)
00596 902 WRITE(IOUT, 9020)
00597 9020 FORMAT(' ----- DEXAY: LACK OF INITIALISATION')
00598 STOP
00599 910 WRITE(IOUT, 9100)
00600 9100 FORMAT(' ----- DEXAY: WRONG VALUE OF KTO ')
00601 STOP
00602 END
00603 SUBROUTINE DEXAY1(KTO,JAKIN,JAK,POL,ISGN)
00604
00605
00606
00607
00608
00609 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00610 REAL*4 GAMPMC ,GAMPER
00611 COMMON / INOUT / INUT,IOUT
00612 REAL POL(4),POLAR(4)
00613 REAL PNU(4),PPI(4)
00614 REAL PRHO(4),PIC(4),PIZ(4)
00615 REAL PWB(4),PMU(4),PNM(4)
00616 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
00617 REAL PKK(4),PKS(4)
00618 REAL PNPI(4,9)
00619 REAL PHOT(4)
00620 REAL PDUM(4)
00621
00622 IF(JAKIN.EQ.-1) RETURN
00623 DO 33 I=1,3
00624 33 POLAR(I)=POL(I)
00625 POLAR(4)=0.
00626 DO 34 I=1,4
00627 34 PDUM(I)=.0
00628 JAK=JAKIN
00629 IF(JAK.EQ.0) CALL JAKER(JAK)
00630
00631 IF(JAK.EQ.1) THEN
00632 CALL DEXEL(0, ISGN,POLAR,PNU,PWB,PMU,PNM,PHOT)
00633 CALL DWLUEL(KTO,ISGN,PNU,PWB,PMU,PNM)
00634 CALL DWRPH(KTO,PHOT )
00635 ELSEIF(JAK.EQ.2) THEN
00636 CALL DEXMU(0, ISGN,POLAR,PNU,PWB,PMU,PNM,PHOT)
00637 CALL DWLUMU(KTO,ISGN,PNU,PWB,PMU,PNM)
00638 CALL DWRPH(KTO,PHOT )
00639 ELSEIF(JAK.EQ.3) THEN
00640 CALL DEXPI(0, ISGN,POLAR,PPI,PNU)
00641 CALL DWLUPI(KTO,ISGN,PPI,PNU)
00642 ELSEIF(JAK.EQ.4) THEN
00643 CALL DEXRO(0, ISGN,POLAR,PNU,PRHO,PIC,PIZ)
00644 CALL DWLURO(KTO,ISGN,PNU,PRHO,PIC,PIZ)
00645 ELSEIF(JAK.EQ.5) THEN
00646 CALL DEXAA(0, ISGN,POLAR,PNU,PAA,PIM1,PIM2,PIPL,JAA)
00647 CALL DWLUAA(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
00648 ELSEIF(JAK.EQ.6) THEN
00649 CALL DEXKK(0, ISGN,POLAR,PKK,PNU)
00650 CALL DWLUKK(KTO,ISGN,PKK,PNU)
00651 ELSEIF(JAK.EQ.7) THEN
00652 CALL DEXKS(0, ISGN,POLAR,PNU,PKS,PKK,PPI,JKST)
00653 CALL DWLUKS(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
00654 ELSE
00655 JNPI=JAK-7
00656 CALL DEXNEW(0, ISGN,POLAR,PNU,PWB,PNPI,JNPI)
00657 CALL DWLNEW(KTO,ISGN,PNU,PWB,PNPI,JAK)
00658 ENDIF
00659 NEVDEC(JAK)=NEVDEC(JAK)+1
00660 END
00661 SUBROUTINE DEXEL(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
00662
00663
00664
00665
00666
00667
00668 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
00669 DATA IWARM/0/
00670
00671 IF(MODE.EQ.-1) THEN
00672
00673 IWARM=1
00674 CALL DADMEL( -1,ISGN,HV,PNU,PWB,Q1,Q2,PH)
00675
00676
00677 ELSEIF(MODE.EQ. 0) THEN
00678
00679 300 CONTINUE
00680 IF(IWARM.EQ.0) GOTO 902
00681 CALL DADMEL( 0,ISGN,HV,PNU,PWB,Q1,Q2,PH)
00682 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
00683
00684 CALL RANMAR(RN,1)
00685 IF(RN(1).GT.WT) GOTO 300
00686
00687 ELSEIF(MODE.EQ. 1) THEN
00688
00689 CALL DADMEL( 1,ISGN,HV,PNU,PWB,Q1,Q2,PH)
00690
00691 ENDIF
00692
00693 RETURN
00694 902 PRINT 9020
00695 9020 FORMAT(' ----- DEXEL: LACK OF INITIALISATION')
00696 STOP
00697 END
00698 SUBROUTINE DEXMU(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
00699
00700
00701
00702
00703
00704
00705
00706
00707 COMMON / INOUT / INUT,IOUT
00708 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
00709 DATA IWARM/0/
00710
00711 IF(MODE.EQ.-1) THEN
00712
00713 IWARM=1
00714 CALL DADMMU( -1,ISGN,HV,PNU,PWB,Q1,Q2,PH)
00715
00716
00717 ELSEIF(MODE.EQ. 0) THEN
00718
00719 300 CONTINUE
00720 IF(IWARM.EQ.0) GOTO 902
00721 CALL DADMMU( 0,ISGN,HV,PNU,PWB,Q1,Q2,PH)
00722 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
00723
00724 CALL RANMAR(RN,1)
00725 IF(RN(1).GT.WT) GOTO 300
00726
00727 ELSEIF(MODE.EQ. 1) THEN
00728
00729 CALL DADMMU( 1,ISGN,HV,PNU,PWB,Q1,Q2,PH)
00730
00731 ENDIF
00732
00733 RETURN
00734 902 WRITE(IOUT, 9020)
00735 9020 FORMAT(' ----- DEXMU: LACK OF INITIALISATION')
00736 STOP
00737 END
00738 SUBROUTINE DADMEL(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
00739
00740
00741
00742
00743 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
00744 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
00745 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00746 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00747 * ,AMK,AMKZ,AMKST,GAMKST
00748
00749 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00750 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00751 * ,AMK,AMKZ,AMKST,GAMKST
00752 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00753 REAL*4 GAMPMC ,GAMPER
00754
00755 REAL*4 PHX(4)
00756
00757 COMMON / INOUT / INUT,IOUT
00758
00759
00760 REAL HHV(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4)
00761 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
00762 REAL*4 RRR(3)
00763 REAL*8 SWT, SSWT
00764 DATA PI /3.141592653589793238462643/
00765 DATA IWARM/0/
00766
00767 IF(MODE.EQ.-1) THEN
00768
00769 IWARM=1
00770 NEVRAW=0
00771 NEVACC=0
00772 NEVOVR=0
00773 SWT=0
00774 SSWT=0
00775 WTMAX=1E-20
00776 DO 15 I=1,500
00777 CALL DPHSEL(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00778 IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2
00779 15 CONTINUE
00780
00781
00782 ELSEIF(MODE.EQ. 0) THEN
00783
00784 300 CONTINUE
00785 IF(IWARM.EQ.0) GOTO 902
00786 NEVRAW=NEVRAW+1
00787 CALL DPHSEL(WT,HV,PNU,PWB,Q1,Q2,PHX)
00788
00789 SWT=SWT+WT
00790 SSWT=SSWT+WT**2
00791 CALL RANMAR(RRR,3)
00792 RN=RRR(1)
00793 IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1
00794 IF(RN*WTMAX.GT.WT) GOTO 300
00795
00796 RR2=RRR(2)
00797 COSTHE=-1.+2.*RR2
00798 THET=ACOS(COSTHE)
00799 RR3=RRR(3)
00800 PHI =2*PI*RR3
00801 CALL ROTOR2(THET,PNU,PNU)
00802 CALL ROTOR3( PHI,PNU,PNU)
00803 CALL ROTOR2(THET,PWB,PWB)
00804 CALL ROTOR3( PHI,PWB,PWB)
00805 CALL ROTOR2(THET,Q1,Q1)
00806 CALL ROTOR3( PHI,Q1,Q1)
00807 CALL ROTOR2(THET,Q2,Q2)
00808 CALL ROTOR3( PHI,Q2,Q2)
00809 CALL ROTOR2(THET,HV,HV)
00810 CALL ROTOR3( PHI,HV,HV)
00811 CALL ROTOR2(THET,PHX,PHX)
00812 CALL ROTOR3( PHI,PHX,PHX)
00813 DO 44,I=1,3
00814 44 HHV(I)=-ISGN*HV(I)
00815 NEVACC=NEVACC+1
00816
00817 ELSEIF(MODE.EQ. 1) THEN
00818
00819 IF(NEVRAW.EQ.0) RETURN
00820 PARGAM=SWT/FLOAT(NEVRAW+1)
00821 ERROR=0
00822 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW))
00823 RAT=PARGAM/GAMEL
00824 WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR
00825
00826 GAMPMC(1)=RAT
00827 GAMPER(1)=ERROR
00828
00829 ENDIF
00830
00831 RETURN
00832 7010 FORMAT(///1X,15(5H*****)
00833 $ /,' *', 25X,'******** DADMEL FINAL REPORT ******** ',9X,1H*
00834 $ /,' *',I20 ,5X,'NEVRAW = NO. OF EL DECAYS TOTAL ',9X,1H*
00835 $ /,' *',I20 ,5X,'NEVACC = NO. OF EL DECS. ACCEPTED ',9X,1H*
00836 $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H*
00837 $ /,' *',E20.5,5X,'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9X,1H*
00838 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
00839 $ /,' *',F20.9,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H*
00840 $ /,' *',25X, 'COMPLETE QED CORRECTIONS INCLUDED ',9X,1H*
00841 $ /,' *',25X, 'BUT ONLY V-A CUPLINGS ',9X,1H*
00842 $ /,1X,15(5H*****)/)
00843 902 WRITE(IOUT, 9020)
00844 9020 FORMAT(' ----- DADMEL: LACK OF INITIALISATION')
00845 STOP
00846 END
00847 SUBROUTINE DADMMU(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
00848
00849 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00850 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00851 * ,AMK,AMKZ,AMKST,GAMKST
00852
00853 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00854 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00855 * ,AMK,AMKZ,AMKST,GAMKST
00856 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
00857 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
00858 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
00859 REAL*4 GAMPMC ,GAMPER
00860 COMMON / INOUT / INUT,IOUT
00861 REAL*4 PHX(4)
00862 REAL HHV(4),HV(4),PNU(4),PWB(4),Q1(4),Q2(4)
00863 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
00864 REAL*4 RRR(3)
00865 REAL*8 SWT, SSWT
00866 DATA PI /3.141592653589793238462643/
00867 DATA IWARM /0/
00868
00869 IF(MODE.EQ.-1) THEN
00870
00871 IWARM=1
00872 NEVRAW=0
00873 NEVACC=0
00874 NEVOVR=0
00875 SWT=0
00876 SSWT=0
00877 WTMAX=1E-20
00878 DO 15 I=1,500
00879 CALL DPHSMU(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5)
00880 IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2
00881 15 CONTINUE
00882
00883
00884 ELSEIF(MODE.EQ. 0) THEN
00885
00886 300 CONTINUE
00887 IF(IWARM.EQ.0) GOTO 902
00888 NEVRAW=NEVRAW+1
00889 CALL DPHSMU(WT,HV,PNU,PWB,Q1,Q2,PHX)
00890
00891 SWT=SWT+WT
00892 SSWT=SSWT+WT**2
00893 CALL RANMAR(RRR,3)
00894 RN=RRR(1)
00895 IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1
00896 IF(RN*WTMAX.GT.WT) GOTO 300
00897
00898 COSTHE=-1.+2.*RRR(2)
00899 THET=ACOS(COSTHE)
00900 PHI =2*PI*RRR(3)
00901 CALL ROTOR2(THET,PNU,PNU)
00902 CALL ROTOR3( PHI,PNU,PNU)
00903 CALL ROTOR2(THET,PWB,PWB)
00904 CALL ROTOR3( PHI,PWB,PWB)
00905 CALL ROTOR2(THET,Q1,Q1)
00906 CALL ROTOR3( PHI,Q1,Q1)
00907 CALL ROTOR2(THET,Q2,Q2)
00908 CALL ROTOR3( PHI,Q2,Q2)
00909 CALL ROTOR2(THET,HV,HV)
00910 CALL ROTOR3( PHI,HV,HV)
00911 CALL ROTOR2(THET,PHX,PHX)
00912 CALL ROTOR3( PHI,PHX,PHX)
00913 DO 44,I=1,3
00914 44 HHV(I)=-ISGN*HV(I)
00915 NEVACC=NEVACC+1
00916
00917 ELSEIF(MODE.EQ. 1) THEN
00918
00919 IF(NEVRAW.EQ.0) RETURN
00920 PARGAM=SWT/FLOAT(NEVRAW+1)
00921 ERROR=0
00922 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW))
00923 RAT=PARGAM/GAMEL
00924 WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR
00925
00926 GAMPMC(2)=RAT
00927 GAMPER(2)=ERROR
00928
00929 ENDIF
00930
00931 RETURN
00932 7010 FORMAT(///1X,15(5H*****)
00933 $ /,' *', 25X,'******** DADMMU FINAL REPORT ******** ',9X,1H*
00934 $ /,' *',I20 ,5X,'NEVRAW = NO. OF MU DECAYS TOTAL ',9X,1H*
00935 $ /,' *',I20 ,5X,'NEVACC = NO. OF MU DECS. ACCEPTED ',9X,1H*
00936 $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H*
00937 $ /,' *',E20.5,5X,'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9X,1H*
00938 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
00939 $ /,' *',F20.9,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H*
00940 $ /,' *',25X, 'COMPLETE QED CORRECTIONS INCLUDED ',9X,1H*
00941 $ /,' *',25X, 'BUT ONLY V-A CUPLINGS ',9X,1H*
00942 $ /,1X,15(5H*****)/)
00943 902 WRITE(IOUT, 9020)
00944 9020 FORMAT(' ----- DADMMU: LACK OF INITIALISATION')
00945 STOP
00946 END
00947 SUBROUTINE DPHSEL(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
00948
00949
00950
00951
00952 REAL*4 PHX(4)
00953 REAL*4 HVX(4),PAAX(4),XAX(4),QPX(4),XNX(4)
00954 REAL*8 HV(4),PH(4),PAA(4),XA(4),QP(4),XN(4)
00955 REAL*8 DGAMT
00956 IELMU=1
00957 CALL DRCMU(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
00958 DO 7 K=1,4
00959 HVX(K)=HV(K)
00960 PHX(K)=PH(K)
00961 PAAX(K)=PAA(K)
00962 XAX(K)=XA(K)
00963 QPX(K)=QP(K)
00964 XNX(K)=XN(K)
00965 7 CONTINUE
00966 DGAMX=DGAMT
00967 END
00968 SUBROUTINE DPHSMU(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
00969
00970
00971
00972
00973 REAL*4 PHX(4)
00974 REAL*4 HVX(4),PAAX(4),XAX(4),QPX(4),XNX(4)
00975 REAL*8 HV(4),PH(4),PAA(4),XA(4),QP(4),XN(4)
00976 REAL*8 DGAMT
00977 IELMU=2
00978 CALL DRCMU(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
00979 DO 7 K=1,4
00980 HVX(K)=HV(K)
00981 PHX(K)=PH(K)
00982 PAAX(K)=PAA(K)
00983 XAX(K)=XA(K)
00984 QPX(K)=QP(K)
00985 XNX(K)=XN(K)
00986 7 CONTINUE
00987 DGAMX=DGAMT
00988 END
00989 SUBROUTINE DRCMU(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
00990 IMPLICIT REAL*8 (A-H,O-Z)
00991
00992
00993
00994
00995 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
00996 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
00997 * ,AMK,AMKZ,AMKST,GAMKST
00998
00999 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01000 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01001 * ,AMK,AMKZ,AMKST,GAMKST
01002 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01003 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01004
01005 COMMON / INOUT / INUT,IOUT
01006 COMMON / TAURAD / XK0DEC,ITDKRC
01007 REAL*8 XK0DEC
01008 REAL*8 HV(4),PT(4),PH(4),PAA(4),XA(4),QP(4),XN(4)
01009 REAL*8 PR(4)
01010 REAL*4 RRR(6)
01011 LOGICAL IHARD
01012 DATA PI /3.141592653589793238462643D0/
01013
01014
01015
01016
01017
01018
01019
01020 PHSPAC=1./2**17/PI**8
01021 AMTAX=AMTAU
01022
01023 PT(1)=0.D0
01024 PT(2)=0.D0
01025 PT(3)=0.D0
01026 PT(4)=AMTAX
01027
01028 CALL RANMAR(RRR,6)
01029
01030 IF (IELMU.EQ.1) THEN
01031 AMU=AMEL
01032 ELSE
01033 AMU=AMMU
01034 ENDIF
01035
01036 PRHARD=0.30D0
01037 IF ( ITDKRC.EQ.0) PRHARD=0D0
01038 PRSOFT=1.-PRHARD
01039 IF(PRSOFT.LT.0.1) THEN
01040 PRINT *, 'ERROR IN DRCMU; PRSOFT=',PRSOFT
01041 STOP
01042 ENDIF
01043
01044 RR5=RRR(5)
01045 IHARD=(RR5.GT.PRSOFT)
01046 IF (IHARD) THEN
01047
01048 RR1=RRR(1)
01049 AMS1=(AMU+AMNUTA)**2
01050 AMS2=(AMTAX)**2
01051 XK1=1-AMS1/AMS2
01052 XL1=LOG(XK1/2/XK0DEC)
01053 XL0=LOG(2*XK0DEC)
01054 XK=EXP(XL1*RR1+XL0)
01055 AM3SQ=(1-XK)*AMS2
01056 AM3 =SQRT(AM3SQ)
01057 PHSPAC=PHSPAC*AMS2*XL1*XK
01058 PHSPAC=PHSPAC/PRHARD
01059 ELSE
01060 AM3=AMTAX
01061 PHSPAC=PHSPAC*2**6*PI**3
01062 PHSPAC=PHSPAC/PRSOFT
01063 ENDIF
01064
01065 RR2=RRR(2)
01066 AMS1=(AMNUTA)**2
01067 AMS2=(AM3-AMU)**2
01068
01069
01070
01071 AM2SQ=AMS1+ RR2*(AMS2-AMS1)
01072 AM2 =SQRT(AM2SQ)
01073 PHSPAC=PHSPAC*(AMS2-AMS1)
01074
01075 ENQ1=(AM2SQ+AMNUTA**2)/(2*AM2)
01076 ENQ2=(AM2SQ-AMNUTA**2)/(2*AM2)
01077 PPI= ENQ1**2-AMNUTA**2
01078 PPPI=SQRT(ABS(ENQ1**2-AMNUTA**2))
01079 PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AM2)
01080
01081 CALL SPHERD(PPPI,XN)
01082 XN(4)=ENQ1
01083
01084 DO 30 I=1,3
01085 30 XA(I)=-XN(I)
01086 XA(4)=ENQ2
01087
01088
01089 PR(1)=0
01090 PR(2)=0
01091 PR(4)=1.D0/(2*AM3)*(AM3**2+AM2**2-AMU**2)
01092 PR(3)= SQRT(ABS(PR(4)**2-AM2**2))
01093 PPI = PR(4)**2-AM2**2
01094
01095 QP(1)=0
01096 QP(2)=0
01097 QP(4)=1.D0/(2*AM3)*(AM3**2-AM2**2+AMU**2)
01098 QP(3)=-PR(3)
01099 PHSPAC=PHSPAC*(4*PI)*(2*PR(3)/AM3)
01100
01101 EXE=(PR(4)+PR(3))/AM2
01102 CALL BOSTD3(EXE,XN,XN)
01103 CALL BOSTD3(EXE,XA,XA)
01104 RR3=RRR(3)
01105 RR4=RRR(4)
01106 IF (IHARD) THEN
01107 EPS=4*(AMU/AMTAX)**2
01108 XL1=LOG((2+EPS)/EPS)
01109 XL0=LOG(EPS)
01110 ETA =EXP(XL1*RR3+XL0)
01111 CTHET=1+EPS-ETA
01112 THET =ACOS(CTHET)
01113 PHSPAC=PHSPAC*XL1/2*ETA
01114 PHI = 2*PI*RR4
01115 CALL ROTPOX(THET,PHI,XN)
01116 CALL ROTPOX(THET,PHI,XA)
01117 CALL ROTPOX(THET,PHI,QP)
01118 CALL ROTPOX(THET,PHI,PR)
01119
01120
01121
01122 PAA(1)=0
01123 PAA(2)=0
01124 PAA(4)=1/(2*AMTAX)*(AMTAX**2+AM3**2)
01125 PAA(3)= SQRT(ABS(PAA(4)**2-AM3**2))
01126 PPI = PAA(4)**2-AM3**2
01127 PHSPAC=PHSPAC*(4*PI)*(2*PAA(3)/AMTAX)
01128
01129 PH(1)=0
01130 PH(2)=0
01131 PH(4)=PAA(3)
01132 PH(3)=-PAA(3)
01133
01134
01135 EXE=(PAA(4)+PAA(3))/AM3
01136 CALL BOSTD3(EXE,XN,XN)
01137 CALL BOSTD3(EXE,XA,XA)
01138 CALL BOSTD3(EXE,QP,QP)
01139 CALL BOSTD3(EXE,PR,PR)
01140 ELSE
01141 THET =ACOS(-1.+2*RR3)
01142 PHI = 2*PI*RR4
01143 CALL ROTPOX(THET,PHI,XN)
01144 CALL ROTPOX(THET,PHI,XA)
01145 CALL ROTPOX(THET,PHI,QP)
01146 CALL ROTPOX(THET,PHI,PR)
01147
01148
01149
01150 PAA(1)=0
01151 PAA(2)=0
01152 PAA(4)=AMTAX
01153 PAA(3)=0
01154
01155 PH(1)=0
01156 PH(2)=0
01157 PH(4)=0
01158 PH(3)=0
01159 ENDIF
01160
01161 CALL DAMPRY(ITDKRC,XK0DEC,PH,XA,QP,XN,AMPLIT,HV)
01162 DGAMT=1/(2.*AMTAX)*AMPLIT*PHSPAC
01163 END
01164 SUBROUTINE DAMPRY(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
01165 IMPLICIT REAL*8 (A-H,O-Z)
01166
01167
01168
01169
01170
01171 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01172 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01173 * ,AMK,AMKZ,AMKST,GAMKST
01174
01175 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01176 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01177 * ,AMK,AMKZ,AMKST,GAMKST
01178 REAL*8 HV(4),QP(4),XN(4),XA(4),XK(4)
01179
01180 HV(4)=1.D0
01181 AK0=XK0DEC*AMTAU
01182 IF(XK(4).LT.0.1D0*AK0) THEN
01183 AMPLIT=THB(ITDKRC,QP,XN,XA,AK0,HV)
01184 ELSE
01185 AMPLIT=SQM2(ITDKRC,QP,XN,XA,XK,AK0,HV)
01186 ENDIF
01187 RETURN
01188 END
01189 FUNCTION SQM2(ITDKRC,QP,XN,XA,XK,AK0,HV)
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202 IMPLICIT REAL*8(A-H,O-Z)
01203 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01204 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01205 * ,AMK,AMKZ,AMKST,GAMKST
01206
01207 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01208 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01209 * ,AMK,AMKZ,AMKST,GAMKST
01210 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01211 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01212 COMMON / QEDPRM /ALFINV,ALFPI,XK0
01213 REAL*8 ALFINV,ALFPI,XK0
01214 REAL*8 QP(4),XN(4),XA(4),XK(4)
01215 REAL*8 R(4)
01216 REAL*8 HV(4)
01217 REAL*8 S0(3),RXA(3),RXK(3),RQP(3)
01218 DATA PI /3.141592653589793238462643D0/
01219
01220 TMASS=AMTAU
01221 GF=GFERMI
01222 ALPHAI=ALFINV
01223 TMASS2=TMASS**2
01224 EMASS2=QP(4)**2-QP(1)**2-QP(2)**2-QP(3)**2
01225 R(4)=TMASS
01226
01227 DO 7 I=1,3
01228 R(1)=0.D0
01229 R(2)=0.D0
01230 R(3)=0.D0
01231 R(I)=TMASS
01232 RXA(I)=R(4)*XA(4)-R(1)*XA(1)-R(2)*XA(2)-R(3)*XA(3)
01233
01234 RXK(I)=R(4)*XK(4)-R(1)*XK(1)-R(2)*XK(2)-R(3)*XK(3)
01235 RQP(I)=R(4)*QP(4)-R(1)*QP(1)-R(2)*QP(2)-R(3)*QP(3)
01236 7 CONTINUE
01237 QPXN=QP(4)*XN(4)-QP(1)*XN(1)-QP(2)*XN(2)-QP(3)*XN(3)
01238 QPXA=QP(4)*XA(4)-QP(1)*XA(1)-QP(2)*XA(2)-QP(3)*XA(3)
01239 QPXK=QP(4)*XK(4)-QP(1)*XK(1)-QP(2)*XK(2)-QP(3)*XK(3)
01240
01241 XNXK=XN(4)*XK(4)-XN(1)*XK(1)-XN(2)*XK(2)-XN(3)*XK(3)
01242 XAXK=XA(4)*XK(4)-XA(1)*XK(1)-XA(2)*XK(2)-XA(3)*XK(3)
01243 TXN=TMASS*XN(4)
01244 TXA=TMASS*XA(4)
01245 TQP=TMASS*QP(4)
01246 TXK=TMASS*XK(4)
01247
01248 X= XNXK/QPXN
01249 Z= TXK/TQP
01250 A= 1+X
01251 B= 1+ X*(1+Z)/2+Z/2
01252 S1= QPXN*TXA*( -EMASS2/QPXK**2*A + 2*TQP/(QPXK*TXK)*B-
01253 $TMASS2/TXK**2) +
01254 $QPXN/TXK**2* ( TMASS2*XAXK - TXA*TXK+ XAXK*TXK) -
01255 $TXA*TXN/TXK - QPXN/(QPXK*TXK)* (TQP*XAXK-TXK*QPXA)
01256 CONST4=256*PI/ALPHAI*GF**2
01257 IF (ITDKRC.EQ.0) CONST4=0D0
01258 SQM2=S1*CONST4
01259 DO 5 I=1,3
01260 S0(I) = QPXN*RXA(I)*(-EMASS2/QPXK**2*A + 2*TQP/(QPXK*TXK)*B-
01261 $ TMASS2/TXK**2) +
01262 $ QPXN/TXK**2* (TMASS2*XAXK - TXA*RXK(I)+ XAXK*RXK(I))-
01263 $ RXA(I)*TXN/TXK - QPXN/(QPXK*TXK)*(RQP(I)*XAXK- RXK(I)*QPXA)
01264 5 HV(I)=S0(I)/S1-1.D0
01265 RETURN
01266 END
01267 FUNCTION THB(ITDKRC,QP,XN,XA,AK0,HV)
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281 IMPLICIT REAL*8(A-H,O-Z)
01282 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01283 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01284 * ,AMK,AMKZ,AMKST,GAMKST
01285
01286 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01287 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01288 * ,AMK,AMKZ,AMKST,GAMKST
01289 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01290 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01291 COMMON / QEDPRM /ALFINV,ALFPI,XK0
01292 REAL*8 ALFINV,ALFPI,XK0
01293 DIMENSION QP(4),XN(4),XA(4)
01294 REAL*8 HV(4)
01295 DIMENSION R(4)
01296 REAL*8 RXA(3),RXN(3),RQP(3)
01297 REAL*8 BORNPL(3),AM3POL(3),XM3POL(3)
01298 DATA PI /3.141592653589793238462643D0/
01299
01300 TMASS=AMTAU
01301 GF=GFERMI
01302 ALPHAI=ALFINV
01303
01304 TMASS2=TMASS**2
01305 R(4)=TMASS
01306 DO 7 I=1,3
01307 R(1)=0.D0
01308 R(2)=0.D0
01309 R(3)=0.D0
01310 R(I)=TMASS
01311 RXA(I)=R(4)*XA(4)-R(1)*XA(1)-R(2)*XA(2)-R(3)*XA(3)
01312 RXN(I)=R(4)*XN(4)-R(1)*XN(1)-R(2)*XN(2)-R(3)*XN(3)
01313
01314 RQP(I)=R(4)*QP(4)-R(1)*QP(1)-R(2)*QP(2)-R(3)*QP(3)
01315 7 CONTINUE
01316
01317 U0=QP(4)/TMASS
01318 U3=SQRT(QP(1)**2+QP(2)**2+QP(3)**2)/TMASS
01319 W3=U3
01320 W0=(XN(4)+XA(4))/TMASS
01321 UP=U0+U3
01322 UM=U0-U3
01323 WP=W0+W3
01324 WM=W0-W3
01325 YU=LOG(UP/UM)/2
01326 YW=LOG(WP/WM)/2
01327 EPS2=U0**2-U3**2
01328 EPS=SQRT(EPS2)
01329 Y=W0**2-W3**2
01330 AL=AK0/TMASS
01331
01332 F0=2*U0/U3*( DILOGT(1-(UM*WM/(UP*WP)))- DILOGT(1-WM/WP) +
01333 $DILOGT(1-UM/UP) -2*YU+ 2*LOG(UP)*(YW+YU) ) +
01334 $1/Y* ( 2*U3*YU + (1-EPS2- 2*Y)*LOG(EPS) ) +
01335 $ 2 - 4*(U0/U3*YU -1)* LOG(2*AL)
01336 FP= YU/(2*U3)*(1 + (1-EPS2)/Y ) + LOG(EPS)/Y
01337 FM= YU/(2*U3)*(1 - (1-EPS2)/Y ) - LOG(EPS)/Y
01338 F3= EPS2*(FP+FM)/2
01339
01340 QPXN=QP(4)*XN(4)-QP(1)*XN(1)-QP(2)*XN(2)-QP(3)*XN(3)
01341 QPXA=QP(4)*XA(4)-QP(1)*XA(1)-QP(2)*XA(2)-QP(3)*XA(3)
01342 XNXA=XN(4)*XA(4)-XN(1)*XA(1)-XN(2)*XA(2)-XN(3)*XA(3)
01343 TXN=TMASS*XN(4)
01344 TXA=TMASS*XA(4)
01345 TQP=TMASS*QP(4)
01346
01347 CONST3=1/(2*ALPHAI*PI)*64*GF**2
01348 IF (ITDKRC.EQ.0) CONST3=0D0
01349 XM3= -( F0* QPXN*TXA + FP*EPS2* TXN*TXA +
01350 $FM* QPXN*QPXA + F3* TMASS2*XNXA )
01351 AM3=XM3*CONST3
01352
01353 BRAK= (GV+GA)**2*TQP*XNXA+(GV-GA)**2*TXA*QPXN
01354 & -(GV**2-GA**2)*TMASS*AMNUTA*QPXA
01355 BORN= 32*(GFERMI**2/2.)*BRAK
01356 DO 5 I=1,3
01357 XM3POL(I)= -( F0* QPXN*RXA(I) + FP*EPS2* TXN*RXA(I) +
01358 $ FM* QPXN* (QPXA + (RXA(I)*TQP-TXA*RQP(I))/TMASS2 ) +
01359 $ F3* (TMASS2*XNXA +TXN*RXA(I) -RXN(I)*TXA) )
01360 AM3POL(I)=XM3POL(I)*CONST3
01361
01362 BORNPL(I)=BORN+(
01363 & (GV+GA)**2*TMASS*XNXA*QP(I)
01364 & -(GV-GA)**2*TMASS*QPXN*XA(I)
01365 & +(GV**2-GA**2)*AMNUTA*TXA*QP(I)
01366 & -(GV**2-GA**2)*AMNUTA*TQP*XA(I) )*
01367 & 32*(GFERMI**2/2.)
01368 5 HV(I)=(BORNPL(I)+AM3POL(I))/(BORN+AM3)-1.D0
01369 THB=BORN+AM3
01370 IF (THB/BORN.LT.0.1D0) THEN
01371 PRINT *, 'ERROR IN THB, THB/BORN=',THB/BORN
01372
01373 THB=0.D0
01374
01375 ENDIF
01376 RETURN
01377 END
01378 SUBROUTINE DEXPI(MODE,ISGN,POL,PPI,PNU)
01379
01380
01381
01382
01383
01384
01385 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
01386
01387 IF(MODE.EQ.-1) THEN
01388
01389 CALL DADMPI(-1,ISGN,HV,PPI,PNU)
01390
01391
01392 ELSEIF(MODE.EQ. 0) THEN
01393
01394 300 CONTINUE
01395 CALL DADMPI( 0,ISGN,HV,PPI,PNU)
01396 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
01397
01398 CALL RANMAR(RN,1)
01399 IF(RN(1).GT.WT) GOTO 300
01400
01401 ELSEIF(MODE.EQ. 1) THEN
01402
01403 CALL DADMPI( 1,ISGN,HV,PPI,PNU)
01404
01405 ENDIF
01406
01407 RETURN
01408 END
01409 SUBROUTINE DADMPI(MODE,ISGN,HV,PPI,PNU)
01410
01411 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01412 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01413 * ,AMK,AMKZ,AMKST,GAMKST
01414
01415 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01416 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01417 * ,AMK,AMKZ,AMKST,GAMKST
01418 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01419 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01420 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
01421 REAL*4 GAMPMC ,GAMPER
01422 COMMON / INOUT / INUT,IOUT
01423 REAL PPI(4),PNU(4),HV(4)
01424 DATA PI /3.141592653589793238462643/
01425
01426 IF(MODE.EQ.-1) THEN
01427
01428 NEVTOT=0
01429 ELSEIF(MODE.EQ. 0) THEN
01430
01431 NEVTOT=NEVTOT+1
01432 EPI= (AMTAU**2+AMPI**2-AMNUTA**2)/(2*AMTAU)
01433 ENU= (AMTAU**2-AMPI**2+AMNUTA**2)/(2*AMTAU)
01434 XPI= SQRT(EPI**2-AMPI**2)
01435
01436 CALL SPHERA(XPI,PPI)
01437 PPI(4)=EPI
01438
01439 DO 30 I=1,3
01440 30 PNU(I)=-PPI(I)
01441 PNU(4)=ENU
01442 PXQ=AMTAU*EPI
01443 PXN=AMTAU*ENU
01444 QXN=PPI(4)*PNU(4)-PPI(1)*PNU(1)-PPI(2)*PNU(2)-PPI(3)*PNU(3)
01445 BRAK=(GV**2+GA**2)*(2*PXQ*QXN-AMPI**2*PXN)
01446 & +(GV**2-GA**2)*AMTAU*AMNUTA*AMPI**2
01447 DO 40 I=1,3
01448 40 HV(I)=-ISGN*2*GA*GV*AMTAU*(2*PPI(I)*QXN-PNU(I)*AMPI**2)/BRAK
01449 HV(4)=1
01450
01451 ELSEIF(MODE.EQ. 1) THEN
01452
01453 IF(NEVTOT.EQ.0) RETURN
01454 FPI=0.1284
01455
01456
01457
01458
01459 GAMM=(GFERMI*FPI)**2/(16.*PI)*AMTAU**3*
01460 $ (BRAK/AMTAU**4)*
01461 $ SQRT((AMTAU**2-AMPI**2-AMNUTA**2)**2
01462 $ -4*AMPI**2*AMNUTA**2 )/AMTAU**2
01463 ERROR=0
01464 RAT=GAMM/GAMEL
01465 WRITE(IOUT, 7010) NEVTOT,GAMM,RAT,ERROR
01466 GAMPMC(3)=RAT
01467 GAMPER(3)=ERROR
01468
01469 ENDIF
01470
01471 RETURN
01472 7010 FORMAT(///1X,15(5H*****)
01473 $ /,' *', 25X,'******** DADMPI FINAL REPORT ******** ',9X,1H*
01474 $ /,' *',I20 ,5X,'NEVTOT = NO. OF PI DECAYS TOTAL ',9X,1H*
01475 $ /,' *',E20.5,5X,'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9X,1H*
01476 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
01477 $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9X,1H*
01478 $ /,1X,15(5H*****)/)
01479 END
01480 SUBROUTINE DEXRO(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
01481
01482
01483
01484
01485
01486
01487
01488
01489 COMMON / INOUT / INUT,IOUT
01490 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
01491 DATA IWARM/0/
01492
01493 IF(MODE.EQ.-1) THEN
01494
01495 IWARM=1
01496 CALL DADMRO( -1,ISGN,HV,PNU,PRO,PIC,PIZ)
01497
01498
01499
01500 ELSEIF(MODE.EQ. 0) THEN
01501
01502 300 CONTINUE
01503 IF(IWARM.EQ.0) GOTO 902
01504 CALL DADMRO( 0,ISGN,HV,PNU,PRO,PIC,PIZ)
01505 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
01506
01507
01508
01509 CALL RANMAR(RN,1)
01510 IF(RN(1).GT.WT) GOTO 300
01511
01512 ELSEIF(MODE.EQ. 1) THEN
01513
01514 CALL DADMRO( 1,ISGN,HV,PNU,PRO,PIC,PIZ)
01515
01516
01517 ENDIF
01518
01519 RETURN
01520 902 WRITE(IOUT, 9020)
01521 9020 FORMAT(' ----- DEXRO: LACK OF INITIALISATION')
01522 STOP
01523 END
01524 SUBROUTINE DADMRO(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
01525
01526 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01527 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01528 * ,AMK,AMKZ,AMKST,GAMKST
01529
01530 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01531 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01532 * ,AMK,AMKZ,AMKST,GAMKST
01533 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01534 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01535 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
01536 REAL*4 GAMPMC ,GAMPER
01537 COMMON / INOUT / INUT,IOUT
01538 REAL HHV(4)
01539 REAL HV(4),PRO(4),PNU(4),PIC(4),PIZ(4)
01540 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
01541 REAL*4 RRR(3)
01542 REAL*8 SWT, SSWT
01543 DATA PI /3.141592653589793238462643/
01544 DATA IWARM/0/
01545
01546 IF(MODE.EQ.-1) THEN
01547
01548 IWARM=1
01549 NEVRAW=0
01550 NEVACC=0
01551 NEVOVR=0
01552 SWT=0
01553 SSWT=0
01554 WTMAX=1E-20
01555 DO 15 I=1,500
01556 CALL DPHSRO(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4)
01557 IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2
01558 15 CONTINUE
01559
01560
01561
01562 ELSEIF(MODE.EQ. 0) THEN
01563
01564 300 CONTINUE
01565 IF(IWARM.EQ.0) GOTO 902
01566 CALL DPHSRO(WT,HV,PNU,PRO,PIC,PIZ)
01567
01568 NEVRAW=NEVRAW+1
01569 SWT=SWT+WT
01570 SSWT=SSWT+WT**2
01571 CALL RANMAR(RRR,3)
01572 RN=RRR(1)
01573 IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1
01574 IF(RN*WTMAX.GT.WT) GOTO 300
01575
01576 COSTHE=-1.+2.*RRR(2)
01577 THET=ACOS(COSTHE)
01578 PHI =2*PI*RRR(3)
01579 CALL ROTOR2(THET,PNU,PNU)
01580 CALL ROTOR3( PHI,PNU,PNU)
01581 CALL ROTOR2(THET,PRO,PRO)
01582 CALL ROTOR3( PHI,PRO,PRO)
01583 CALL ROTOR2(THET,PIC,PIC)
01584 CALL ROTOR3( PHI,PIC,PIC)
01585 CALL ROTOR2(THET,PIZ,PIZ)
01586 CALL ROTOR3( PHI,PIZ,PIZ)
01587 CALL ROTOR2(THET,HV,HV)
01588 CALL ROTOR3( PHI,HV,HV)
01589 DO 44 I=1,3
01590 44 HHV(I)=-ISGN*HV(I)
01591 NEVACC=NEVACC+1
01592
01593 ELSEIF(MODE.EQ. 1) THEN
01594
01595 IF(NEVRAW.EQ.0) RETURN
01596 PARGAM=SWT/FLOAT(NEVRAW+1)
01597 ERROR=0
01598 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW))
01599 RAT=PARGAM/GAMEL
01600 WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR
01601
01602 GAMPMC(4)=RAT
01603 GAMPER(4)=ERROR
01604
01605 ENDIF
01606
01607 RETURN
01608 7003 FORMAT(///1X,15(5H*****)
01609 $ /,' *', 25X,'******** DADMRO INITIALISATION ********',9X,1H*
01610 $ /,' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H*
01611 $ /,1X,15(5H*****)/)
01612 7010 FORMAT(///1X,15(5H*****)
01613 $ /,' *', 25X,'******** DADMRO FINAL REPORT ******** ',9X,1H*
01614 $ /,' *',I20 ,5X,'NEVRAW = NO. OF RHO DECAYS TOTAL ',9X,1H*
01615 $ /,' *',I20 ,5X,'NEVACC = NO. OF RHO DECS. ACCEPTED ',9X,1H*
01616 $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H*
01617 $ /,' *',E20.5,5X,'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9X,1H*
01618 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
01619 $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H*
01620 $ /,1X,15(5H*****)/)
01621 902 WRITE(IOUT, 9020)
01622 9020 FORMAT(' ----- DADMRO: LACK OF INITIALISATION')
01623 STOP
01624 END
01625 SUBROUTINE DPHSRO(DGAMT,HV,PN,PR,PIC,PIZ)
01626
01627
01628
01629
01630 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01631 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01632 * ,AMK,AMKZ,AMKST,GAMKST
01633
01634 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01635 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01636 * ,AMK,AMKZ,AMKST,GAMKST
01637 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01638 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01639 REAL HV(4),PT(4),PN(4),PR(4),PIC(4),PIZ(4),QQ(4),RR1(1)
01640 DATA PI /3.141592653589793238462643/
01641 DATA ICONT /0/
01642
01643
01644 PHSPAC=1./2**11/PI**5
01645
01646 PT(1)=0.
01647 PT(2)=0.
01648 PT(3)=0.
01649 PT(4)=AMTAU
01650
01651 AMS1=(AMPI+AMPIZ)**2
01652 AMS2=(AMTAU-AMNUTA)**2
01653
01654
01655
01656
01657
01658
01659
01660 ALP1=ATAN((AMS1-AMRO**2)/AMRO/GAMRO)
01661 ALP2=ATAN((AMS2-AMRO**2)/AMRO/GAMRO)
01662
01663 100 CONTINUE
01664 CALL RANMAR(RR1,1)
01665 ALP=ALP1+RR1(1)*(ALP2-ALP1)
01666 AMX2=AMRO**2+AMRO*GAMRO*TAN(ALP)
01667 AMX=SQRT(AMX2)
01668 IF(AMX.LT.2.*AMPI) GO TO 100
01669
01670 PHSPAC=PHSPAC*((AMX2-AMRO**2)**2+(AMRO*GAMRO)**2)/(AMRO*GAMRO)
01671 PHSPAC=PHSPAC*(ALP2-ALP1)
01672
01673
01674 PN(1)=0
01675 PN(2)=0
01676 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AMX**2)
01677 PN(3)=-SQRT((PN(4)-AMNUTA)*(PN(4)+AMNUTA))
01678
01679 PR(1)=0
01680 PR(2)=0
01681 PR(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AMX**2)
01682 PR(3)=-PN(3)
01683 PHSPAC=PHSPAC*(4*PI)*(2*PR(3)/AMTAU)
01684
01685
01686 ENQ1=(AMX2+AMPI**2-AMPIZ**2)/(2.*AMX)
01687 ENQ2=(AMX2-AMPI**2+AMPIZ**2)/(2.*AMX)
01688 PPPI=SQRT((ENQ1-AMPI)*(ENQ1+AMPI))
01689 PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AMX)
01690
01691 CALL SPHERA(PPPI,PIC)
01692 PIC(4)=ENQ1
01693
01694 DO 20 I=1,3
01695 20 PIZ(I)=-PIC(I)
01696 PIZ(4)=ENQ2
01697 EXE=(PR(4)+PR(3))/AMX
01698
01699 CALL BOSTR3(EXE,PIC,PIC)
01700 CALL BOSTR3(EXE,PIZ,PIZ)
01701 DO 30 I=1,4
01702 30 QQ(I)=PIC(I)-PIZ(I)
01703
01704 PRODPQ=PT(4)*QQ(4)
01705 PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
01706 PRODPN=PT(4)*PN(4)
01707 QQ2= QQ(4)**2-QQ(1)**2-QQ(2)**2-QQ(3)**2
01708 BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
01709 & +(GV**2-GA**2)*AMTAU*AMNUTA*QQ2
01710 AMPLIT=(GFERMI*CCABIB)**2*BRAK*2*FPIRHO(AMX)
01711 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
01712 DO 40 I=1,3
01713 40 HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
01714 RETURN
01715 END
01716 SUBROUTINE DEXAA(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727 COMMON / INOUT / INUT,IOUT
01728 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
01729 DATA IWARM/0/
01730
01731 IF(MODE.EQ.-1) THEN
01732
01733 IWARM=1
01734 CALL DADMAA( -1,ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
01735
01736
01737 ELSEIF(MODE.EQ. 0) THEN
01738
01739 300 CONTINUE
01740 IF(IWARM.EQ.0) GOTO 902
01741 CALL DADMAA( 0,ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
01742 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
01743
01744 CALL RANMAR(RN,1)
01745 IF(RN(1).GT.WT) GOTO 300
01746
01747 ELSEIF(MODE.EQ. 1) THEN
01748
01749 CALL DADMAA( 1,ISGN,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
01750
01751 ENDIF
01752
01753 RETURN
01754 902 WRITE(IOUT, 9020)
01755 9020 FORMAT(' ----- DEXAA: LACK OF INITIALISATION')
01756 STOP
01757 END
01758 SUBROUTINE DADMAA(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
01759
01760
01761
01762 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01763 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01764 * ,AMK,AMKZ,AMKST,GAMKST
01765
01766 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01767 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01768 * ,AMK,AMKZ,AMKST,GAMKST
01769 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01770 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01771 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
01772 REAL*4 GAMPMC ,GAMPER
01773 COMMON / INOUT / INUT,IOUT
01774 REAL HHV(4)
01775 REAL HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4)
01776 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
01777 REAL*4 RRR(3)
01778 REAL*8 SWT, SSWT
01779 DATA PI /3.141592653589793238462643/
01780 DATA IWARM/0/
01781
01782 IF(MODE.EQ.-1) THEN
01783
01784 IWARM=1
01785 NEVRAW=0
01786 NEVACC=0
01787 NEVOVR=0
01788 SWT=0
01789 SSWT=0
01790 WTMAX=1E-20
01791 DO 15 I=1,500
01792 CALL DPHSAA(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,PDUM5,JAA)
01793 IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2
01794 15 CONTINUE
01795
01796
01797 ELSEIF(MODE.EQ. 0) THEN
01798
01799 300 CONTINUE
01800 IF(IWARM.EQ.0) GOTO 902
01801 CALL DPHSAA(WT,HV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
01802
01803 NEVRAW=NEVRAW+1
01804 SWT=SWT+WT
01805
01806
01807
01808 SSWT=SSWT+dble(WT)**2
01809
01810
01811 CALL RANMAR(RRR,3)
01812 RN=RRR(1)
01813 IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1
01814 IF(RN*WTMAX.GT.WT) GOTO 300
01815
01816 COSTHE=-1.+2.*RRR(2)
01817 THET=ACOS(COSTHE)
01818 PHI =2*PI*RRR(3)
01819 CALL ROTPOL(THET,PHI,PNU)
01820 CALL ROTPOL(THET,PHI,PAA)
01821 CALL ROTPOL(THET,PHI,PIM1)
01822 CALL ROTPOL(THET,PHI,PIM2)
01823 CALL ROTPOL(THET,PHI,PIPL)
01824 CALL ROTPOL(THET,PHI,HV)
01825 DO 44 I=1,3
01826 44 HHV(I)=-ISGN*HV(I)
01827 NEVACC=NEVACC+1
01828
01829 ELSEIF(MODE.EQ. 1) THEN
01830
01831 IF(NEVRAW.EQ.0) RETURN
01832 PARGAM=SWT/FLOAT(NEVRAW+1)
01833 ERROR=0
01834 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW))
01835 RAT=PARGAM/GAMEL
01836 WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR
01837
01838 GAMPMC(5)=RAT
01839 GAMPER(5)=ERROR
01840
01841 ENDIF
01842
01843 RETURN
01844 7003 FORMAT(///1X,15(5H*****)
01845 $ /,' *', 25X,'******** DADMAA INITIALISATION ********',9X,1H*
01846 $ /,' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H*
01847 $ /,1X,15(5H*****)/)
01848 7010 FORMAT(///1X,15(5H*****)
01849 $ /,' *', 25X,'******** DADMAA FINAL REPORT ******** ',9X,1H*
01850 $ /,' *',I20 ,5X,'NEVRAW = NO. OF A1 DECAYS TOTAL ',9X,1H*
01851 $ /,' *',I20 ,5X,'NEVACC = NO. OF A1 DECS. ACCEPTED ',9X,1H*
01852 $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H*
01853 $ /,' *',E20.5,5X,'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9X,1H*
01854 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
01855 $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H*
01856 $ /,1X,15(5H*****)/)
01857 902 WRITE(IOUT, 9020)
01858 9020 FORMAT(' ----- DADMAA: LACK OF INITIALISATION')
01859 STOP
01860 END
01861 SUBROUTINE DPHSAA(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
01862
01863
01864
01865
01866 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01867 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01868 * ,AMK,AMKZ,AMKST,GAMKST
01869
01870 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01871 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01872 * ,AMK,AMKZ,AMKST,GAMKST
01873 COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
01874 REAL*4 BRA1,BRK0,BRK0B,BRKS
01875 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
01876
01877
01878 REAL*4 RRR(1)
01879
01880 MNUM=0
01881
01882 KEYT=1
01883 CALL RANMAR(RRR,1)
01884 RMOD=RRR(1)
01885 IF (RMOD.LT.BRA1) THEN
01886 JAA=1
01887 AMP1=AMPI
01888 AMP2=AMPI
01889 AMP3=AMPI
01890 ELSE
01891 JAA=2
01892 AMP1=AMPIZ
01893 AMP2=AMPIZ
01894 AMP3=AMPI
01895 ENDIF
01896 CALL
01897 $ DPHTRE(DGAMT,HV,PN,PAA,PIM1,AMP1,PIM2,AMP2,PIPL,AMP3,KEYT,MNUM)
01898 END
01899 SUBROUTINE DEXKK(MODE,ISGN,POL,PKK,PNU)
01900
01901
01902
01903
01904
01905
01906 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
01907
01908 IF(MODE.EQ.-1) THEN
01909
01910 CALL DADMKK(-1,ISGN,HV,PKK,PNU)
01911
01912
01913 ELSEIF(MODE.EQ. 0) THEN
01914
01915 300 CONTINUE
01916 CALL DADMKK( 0,ISGN,HV,PKK,PNU)
01917 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
01918
01919 CALL RANMAR(RN,1)
01920 IF(RN(1).GT.WT) GOTO 300
01921
01922 ELSEIF(MODE.EQ. 1) THEN
01923
01924 CALL DADMKK( 1,ISGN,HV,PKK,PNU)
01925
01926 ENDIF
01927
01928 RETURN
01929 END
01930 SUBROUTINE DADMKK(MODE,ISGN,HV,PKK,PNU)
01931
01932
01933
01934 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01935 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
01936
01937 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01938 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01939 * ,AMK,AMKZ,AMKST,GAMKST
01940
01941 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01942 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01943 * ,AMK,AMKZ,AMKST,GAMKST
01944
01945
01946 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
01947 REAL*4 GAMPMC ,GAMPER
01948 COMMON / INOUT / INUT,IOUT
01949 REAL PKK(4),PNU(4),HV(4)
01950 DATA PI /3.141592653589793238462643/
01951
01952 IF(MODE.EQ.-1) THEN
01953
01954 NEVTOT=0
01955 ELSEIF(MODE.EQ. 0) THEN
01956
01957 NEVTOT=NEVTOT+1
01958 EKK= (AMTAU**2+AMK**2-AMNUTA**2)/(2*AMTAU)
01959 ENU= (AMTAU**2-AMK**2+AMNUTA**2)/(2*AMTAU)
01960 XKK= SQRT(EKK**2-AMK**2)
01961
01962 CALL SPHERA(XKK,PKK)
01963 PKK(4)=EKK
01964
01965 DO 30 I=1,3
01966 30 PNU(I)=-PKK(I)
01967 PNU(4)=ENU
01968 PXQ=AMTAU*EKK
01969 PXN=AMTAU*ENU
01970 QXN=PKK(4)*PNU(4)-PKK(1)*PNU(1)-PKK(2)*PNU(2)-PKK(3)*PNU(3)
01971 BRAK=(GV**2+GA**2)*(2*PXQ*QXN-AMK**2*PXN)
01972 & +(GV**2-GA**2)*AMTAU*AMNUTA*AMK**2
01973 DO 40 I=1,3
01974 40 HV(I)=-ISGN*2*GA*GV*AMTAU*(2*PKK(I)*QXN-PNU(I)*AMK**2)/BRAK
01975 HV(4)=1
01976
01977 ELSEIF(MODE.EQ. 1) THEN
01978
01979 IF(NEVTOT.EQ.0) RETURN
01980 FKK=0.0354
01981
01982
01983
01984
01985
01986 GAMM=(GFERMI*FKK)**2/(16.*PI)*AMTAU**3*
01987 $ (BRAK/AMTAU**4)*
01988 $ SQRT((AMTAU**2-AMK**2-AMNUTA**2)**2
01989 $ -4*AMK**2*AMNUTA**2 )/AMTAU**2
01990 ERROR=0
01991
01992 ERROR=0
01993 RAT=GAMM/GAMEL
01994 WRITE(IOUT, 7010) NEVTOT,GAMM,RAT,ERROR
01995 GAMPMC(6)=RAT
01996 GAMPER(6)=ERROR
01997
01998 ENDIF
01999
02000 RETURN
02001 7010 FORMAT(///1X,15(5H*****)
02002 $ /,' *', 25X,'******** DADMKK FINAL REPORT ********',9X,1H*
02003 $ /,' *',I20 ,5X,'NEVTOT = NO. OF K DECAYS TOTAL ',9X,1H*,
02004 $ /,' *',E20.5,5X,'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9X,1H*,
02005 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
02006 $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9X,1H*
02007 $ /,1X,15(5H*****)/)
02008 END
02009 SUBROUTINE DEXKS(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021 COMMON / INOUT / INUT,IOUT
02022 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
02023 DATA IWARM/0/
02024
02025 IF(MODE.EQ.-1) THEN
02026
02027 IWARM=1
02028
02029 CALL DADMKS( -1,ISGN,HV,PNU,PKS,PKK,PPI,JKST)
02030
02031
02032
02033 ELSEIF(MODE.EQ. 0) THEN
02034
02035 300 CONTINUE
02036 IF(IWARM.EQ.0) GOTO 902
02037 CALL DADMKS( 0,ISGN,HV,PNU,PKS,PKK,PPI,JKST)
02038 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
02039
02040
02041
02042 CALL RANMAR(RN,1)
02043 IF(RN(1).GT.WT) GOTO 300
02044
02045 ELSEIF(MODE.EQ. 1) THEN
02046
02047 CALL DADMKS( 1,ISGN,HV,PNU,PKS,PKK,PPI,JKST)
02048
02049
02050 ENDIF
02051
02052 RETURN
02053 902 WRITE(IOUT, 9020)
02054 9020 FORMAT(' ----- DEXKS: LACK OF INITIALISATION')
02055 STOP
02056 END
02057 SUBROUTINE DADMKS(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
02058
02059 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02060 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02061 * ,AMK,AMKZ,AMKST,GAMKST
02062
02063 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02064 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02065 * ,AMK,AMKZ,AMKST,GAMKST
02066 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02067 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02068 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
02069 REAL*4 GAMPMC ,GAMPER
02070 COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
02071 REAL*4 BRA1,BRK0,BRK0B,BRKS
02072 COMMON / INOUT / INUT,IOUT
02073 REAL HHV(4)
02074 REAL HV(4),PKS(4),PNU(4),PKK(4),PPI(4)
02075 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
02076 REAL*4 RRR(3),RMOD(1)
02077 REAL*8 SWT, SSWT
02078 DATA PI /3.141592653589793238462643/
02079 DATA IWARM/0/
02080
02081 IF(MODE.EQ.-1) THEN
02082
02083 IWARM=1
02084 NEVRAW=0
02085 NEVACC=0
02086 NEVOVR=0
02087 SWT=0
02088 SSWT=0
02089 WTMAX=1E-20
02090 DO 15 I=1,500
02091
02092 JKST=10
02093 CALL DPHSKS(WT,HV,PDUM1,PDUM2,PDUM3,PDUM4,JKST)
02094 IF(WT.GT.WTMAX/1.2) WTMAX=WT*1.2
02095 15 CONTINUE
02096
02097
02098
02099 ELSEIF(MODE.EQ. 0) THEN
02100
02101 IF(IWARM.EQ.0) GOTO 902
02102
02103
02104 DEC1=BRKS
02105 400 CONTINUE
02106 CALL RANMAR(RMOD,1)
02107 IF(RMOD(1).LT.DEC1) THEN
02108 JKST=10
02109 ELSE
02110 JKST=20
02111 ENDIF
02112 CALL DPHSKS(WT,HV,PNU,PKS,PKK,PPI,JKST)
02113 CALL RANMAR(RRR,3)
02114 RN=RRR(1)
02115 IF(WT.GT.WTMAX) NEVOVR=NEVOVR+1
02116 NEVRAW=NEVRAW+1
02117 SWT=SWT+WT
02118 SSWT=SSWT+WT**2
02119 IF(RN*WTMAX.GT.WT) GOTO 400
02120
02121 COSTHE=-1.+2.*RRR(2)
02122 THET=ACOS(COSTHE)
02123 PHI =2*PI*RRR(3)
02124 CALL ROTOR2(THET,PNU,PNU)
02125 CALL ROTOR3( PHI,PNU,PNU)
02126 CALL ROTOR2(THET,PKS,PKS)
02127 CALL ROTOR3( PHI,PKS,PKS)
02128 CALL ROTOR2(THET,PKK,PKK)
02129 CALL ROTOR3(PHI,PKK,PKK)
02130 CALL ROTOR2(THET,PPI,PPI)
02131 CALL ROTOR3( PHI,PPI,PPI)
02132 CALL ROTOR2(THET,HV,HV)
02133 CALL ROTOR3( PHI,HV,HV)
02134 DO 44 I=1,3
02135 44 HHV(I)=-ISGN*HV(I)
02136 NEVACC=NEVACC+1
02137
02138 ELSEIF(MODE.EQ. 1) THEN
02139
02140 IF(NEVRAW.EQ.0) RETURN
02141 PARGAM=SWT/FLOAT(NEVRAW+1)
02142 ERROR=0
02143 IF(NEVRAW.NE.0) ERROR=SQRT(SSWT/SWT**2-1./FLOAT(NEVRAW))
02144 RAT=PARGAM/GAMEL
02145 WRITE(IOUT, 7010) NEVRAW,NEVACC,NEVOVR,PARGAM,RAT,ERROR
02146
02147 GAMPMC(7)=RAT
02148 GAMPER(7)=ERROR
02149
02150 ENDIF
02151
02152 RETURN
02153 7003 FORMAT(///1X,15(5H*****)
02154 $ /,' *', 25X,'******** DADMKS INITIALISATION ********',9X,1H*
02155 $ /,' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H*
02156 $ /,1X,15(5H*****)/)
02157 7010 FORMAT(///1X,15(5H*****)
02158 $ /,' *', 25X,'******** DADMKS FINAL REPORT ********',9X,1H*
02159 $ /,' *',I20 ,5X,'NEVRAW = NO. OF K* DECAYS TOTAL ',9X,1H*,
02160 $ /,' *',I20 ,5X,'NEVACC = NO. OF K* DECS. ACCEPTED ',9X,1H*,
02161 $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H*
02162 $ /,' *',E20.5,5X,'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9X,1H*,
02163 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
02164 $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H*
02165 $ /,1X,15(5H*****)/)
02166 902 WRITE(IOUT, 9020)
02167 9020 FORMAT(' ----- DADMKS: LACK OF INITIALISATION')
02168 STOP
02169 END
02170 SUBROUTINE DPHSKS(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
02171
02172
02173
02174
02175
02176
02177
02178 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02179 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02180
02181 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02182 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02183 * ,AMK,AMKZ,AMKST,GAMKST
02184
02185 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02186 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02187 * ,AMK,AMKZ,AMKST,GAMKST
02188
02189
02190 REAL HV(4),PT(4),PN(4),PKS(4),PKK(4),PPI(4),QQ(4),RR1(1)
02191
02192 COMPLEX BWIGS
02193
02194 DATA PI /3.141592653589793238462643/
02195
02196 DATA ICONT /0/
02197
02198 PHSPAC=1./2**11/PI**5
02199
02200 PT(1)=0.
02201 PT(2)=0.
02202 PT(3)=0.
02203 PT(4)=AMTAU
02204 CALL RANMAR(RR1,1)
02205
02206 IF(JKST.EQ.10)THEN
02207
02208
02209 AMS1=(AMPI+AMKZ)**2
02210 AMS2=(AMTAU-AMNUTA)**2
02211
02212
02213
02214
02215
02216 ALP1=ATAN((AMS1-AMKST**2)/AMKST/GAMKST)
02217 ALP2=ATAN((AMS2-AMKST**2)/AMKST/GAMKST)
02218 ALP=ALP1+RR1(1)*(ALP2-ALP1)
02219 AMX2=AMKST**2+AMKST*GAMKST*TAN(ALP)
02220 AMX=SQRT(AMX2)
02221 PHSPAC=PHSPAC*((AMX2-AMKST**2)**2+(AMKST*GAMKST)**2)
02222 & /(AMKST*GAMKST)
02223 PHSPAC=PHSPAC*(ALP2-ALP1)
02224
02225
02226 PN(1)=0
02227 PN(2)=0
02228 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AMX**2)
02229 PN(3)=-SQRT((PN(4)-AMNUTA)*(PN(4)+AMNUTA))
02230
02231
02232 PKS(1)=0
02233 PKS(2)=0
02234 PKS(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AMX**2)
02235 PKS(3)=-PN(3)
02236 PHSPAC=PHSPAC*(4*PI)*(2*PKS(3)/AMTAU)
02237
02238
02239 ENPI=( AMX**2+AMPI**2-AMKZ**2 ) / ( 2*AMX )
02240 PPPI=SQRT((ENPI-AMPI)*(ENPI+AMPI))
02241 PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AMX)
02242
02243 CALL SPHERA(PPPI,PPI)
02244 PPI(4)=ENPI
02245
02246 DO 20 I=1,3
02247 20 PKK(I)=-PPI(I)
02248 PKK(4)=( AMX**2+AMKZ**2-AMPI**2 ) / ( 2*AMX )
02249 EXE=(PKS(4)+PKS(3))/AMX
02250
02251 CALL BOSTR3(EXE,PPI,PPI)
02252 CALL BOSTR3(EXE,PKK,PKK)
02253 DO 30 I=1,4
02254 30 QQ(I)=PPI(I)-PKK(I)
02255
02256 PKSD =PKS(4)*PKS(4)-PKS(3)*PKS(3)-PKS(2)*PKS(2)-PKS(1)*PKS(1)
02257 QQPKS=PKS(4)* QQ(4)-PKS(3)* QQ(3)-PKS(2)* QQ(2)-PKS(1)* QQ(1)
02258 DO 31 I=1,4
02259 31 QQ(I)=QQ(I)-PKS(I)*QQPKS/PKSD
02260
02261 PRODPQ=PT(4)*QQ(4)
02262 PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
02263 PRODPN=PT(4)*PN(4)
02264 QQ2= QQ(4)**2-QQ(1)**2-QQ(2)**2-QQ(3)**2
02265 BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
02266 & +(GV**2-GA**2)*AMTAU*AMNUTA*QQ2
02267
02268
02269 FKS=CABS(BWIGS(AMX2,AMKST,GAMKST))**2
02270
02271 AMPLIT=(GFERMI*SCABIB)**2*BRAK*2*FKS
02272 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
02273 DO 40 I=1,3
02274 40 HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
02275
02276
02277 ELSEIF(JKST.EQ.20)THEN
02278
02279
02280 AMS1=(AMPIZ+AMK)**2
02281 AMS2=(AMTAU-AMNUTA)**2
02282
02283
02284
02285
02286
02287
02288
02289 ALP1=ATAN((AMS1-AMKST**2)/AMKST/GAMKST)
02290 ALP2=ATAN((AMS2-AMKST**2)/AMKST/GAMKST)
02291 ALP=ALP1+RR1(1)*(ALP2-ALP1)
02292 AMX2=AMKST**2+AMKST*GAMKST*TAN(ALP)
02293 AMX=SQRT(AMX2)
02294 PHSPAC=PHSPAC*((AMX2-AMKST**2)**2+(AMKST*GAMKST)**2)
02295 & /(AMKST*GAMKST)
02296 PHSPAC=PHSPAC*(ALP2-ALP1)
02297
02298
02299 PN(1)=0
02300 PN(2)=0
02301 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AMX**2)
02302 PN(3)=-SQRT((PN(4)-AMNUTA)*(PN(4)+AMNUTA))
02303
02304 PKS(1)=0
02305 PKS(2)=0
02306 PKS(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AMX**2)
02307 PKS(3)=-PN(3)
02308 PHSPAC=PHSPAC*(4*PI)*(2*PKS(3)/AMTAU)
02309
02310
02311 ENPI=( AMX**2+AMPIZ**2-AMK**2 ) / ( 2*AMX )
02312 PPPI=SQRT((ENPI-AMPIZ)*(ENPI+AMPIZ))
02313 PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AMX)
02314
02315 CALL SPHERA(PPPI,PPI)
02316 PPI(4)=ENPI
02317
02318 DO 50 I=1,3
02319 50 PKK(I)=-PPI(I)
02320 PKK(4)=( AMX**2+AMK**2-AMPIZ**2 ) / ( 2*AMX )
02321 EXE=(PKS(4)+PKS(3))/AMX
02322
02323 CALL BOSTR3(EXE,PPI,PPI)
02324 CALL BOSTR3(EXE,PKK,PKK)
02325 DO 60 I=1,4
02326 60 QQ(I)=PKK(I)-PPI(I)
02327
02328 PKSD =PKS(4)*PKS(4)-PKS(3)*PKS(3)-PKS(2)*PKS(2)-PKS(1)*PKS(1)
02329 QQPKS=PKS(4)* QQ(4)-PKS(3)* QQ(3)-PKS(2)* QQ(2)-PKS(1)* QQ(1)
02330 DO 61 I=1,4
02331 61 QQ(I)=QQ(I)-PKS(I)*QQPKS/PKSD
02332
02333 PRODPQ=PT(4)*QQ(4)
02334 PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
02335 PRODPN=PT(4)*PN(4)
02336 QQ2= QQ(4)**2-QQ(1)**2-QQ(2)**2-QQ(3)**2
02337 BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
02338 & +(GV**2-GA**2)*AMTAU*AMNUTA*QQ2
02339
02340
02341 FKS=CABS(BWIGS(AMX2,AMKST,GAMKST))**2
02342
02343 AMPLIT=(GFERMI*SCABIB)**2*BRAK*2*FKS
02344 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
02345 DO 70 I=1,3
02346 70 HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
02347 ENDIF
02348 RETURN
02349 END
02350
02351
02352
02353
02354 SUBROUTINE DPHNPI(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
02355
02356
02357
02358
02359
02360 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02361 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02362 * ,AMK,AMKZ,AMKST,GAMKST
02363
02364 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02365 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02366 * ,AMK,AMKZ,AMKST,GAMKST
02367 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02368 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02369 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
02370
02371 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
02372 & ,NAMES
02373 CHARACTER NAMES(NMODE)*31
02374 REAL*8 WETMAX(500)
02375
02376
02377
02378 REAL*8 PN(4),PR(4),PPI(4,9),HV(4)
02379 REAL*4 PNX(4),PRX(4),PPIX(4,9),HVX(4)
02380 REAL*8 PV(5,9),PT(4),UE(3),BE(3)
02381 REAL*8 PAWT,AMX,AMS1,AMS2,PA,PHS,PHSMAX,PMIN,PMAX
02382
02383 REAL*8 PHSPAC
02384
02385 REAL*8 GAM,BEP,PHI,A,B,C
02386 REAL*8 AMPIK
02387 REAL*4 RRR(9),RRX(2),RN(1),RR2(1)
02388
02389 DATA PI /3.141592653589793238462643/
02390 DATA WETMAX /500*1D-15/
02391
02392
02393
02394 PAWT(A,B,C)=
02395 $ SQRT(MAX(0.D0,(A**2-(B+C)**2)*(A**2-(B-C)**2)))/(2.D0*A)
02396
02397
02398 AMPIK(I,J)=DCDMAS(IDFFIN(I,J))
02399
02400
02401
02402 IF ((JNPI.LE.0).OR.JNPI.GT.100) THEN
02403 WRITE(6,*) 'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',JNPI
02404 STOP
02405 ENDIF
02406
02407
02408
02409 PT(1)=0.
02410 PT(2)=0.
02411 PT(3)=0.
02412 PT(4)=AMTAU
02413
02414
02415 500 CONTINUE
02416
02417
02418 ND=MULPIK(JNPI)
02419 PS=0.
02420 PHSPAC = 1./2.**5 /PI**2
02421 DO 4 I=1,ND
02422 4 PS =PS+AMPIK(I,JNPI)
02423
02424 CALL RANMAR(RR2,1)
02425
02426 AMS1=PS**2
02427 AMS2=(AMTAU-AMNUTA)**2
02428
02429
02430
02431 AMX2=AMS1+ RR2(1)*(AMS2-AMS1)
02432
02433 AMX =SQRT(AMX2)
02434 AMW =AMX
02435 PHSPAC=PHSPAC * (AMS2-AMS1)
02436
02437
02438 PN(1)=0
02439 PN(2)=0
02440 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AMX2)
02441 PN(3)=-SQRT((PN(4)-AMNUTA)*(PN(4)+AMNUTA))
02442
02443 PR(1)=0
02444 PR(2)=0
02445 PR(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AMX2)
02446 PR(3)=-PN(3)
02447 PHSPAC=PHSPAC * (4.*PI) * (2.*PR(3)/AMTAU)
02448
02449
02450
02451
02452 PXQ=AMTAU*PR(4)
02453 PXN=AMTAU*PN(4)
02454 QXN=PR(4)*PN(4)-PR(1)*PN(1)-PR(2)*PN(2)-PR(3)*PN(3)
02455
02456
02457
02458
02459 BRAK=2*(GV**2+GA**2)*(2*PXQ*QXN+AMX2*PXN)
02460 & -6*(GV**2-GA**2)*AMTAU*AMNUTA*AMX2
02461
02462
02463 jn=JNPI-nm4-nm5+3
02464
02465
02466 AMPLIT=CCABIB**2*GFERMI**2/2.* BRAK*AMX2*SIGEE(AMX2,JN)
02467 DGAMT=1./(2.*AMTAU)*AMPLIT*PHSPAC
02468
02469
02470
02471 PHSMAX = 1.
02472
02473 DO 200 I=1,4
02474 200 PV(I,1)=PR(I)
02475 PV(5,1)=AMW
02476 PV(5,ND)=AMPIK(ND,JNPI)
02477
02478 PMAX=AMW-PS+AMPIK(ND,JNPI)
02479 PMIN=.0
02480 DO 220 IL=ND-1,1,-1
02481 PMAX=PMAX+AMPIK(IL,JNPI)
02482 PMIN=PMIN+AMPIK(IL+1,JNPI)
02483
02484 220 PHSMAX=PHSMAX*PAWT(PMAX,PMIN,AMPIK(IL,JNPI))/PMAX
02485
02486
02487 AMX=AMW
02488 DO 222 IL=1,ND-2
02489 AMS1=.0
02490 DO 223 JL=IL+1,ND
02491 223 AMS1=AMS1+AMPIK(JL,JNPI)
02492 AMS1=AMS1**2
02493 AMX =(AMX-AMPIK(IL,JNPI))
02494 AMS2=(AMX)**2
02495 PHSMAX=PHSMAX * (AMS2-AMS1)
02496 222 CONTINUE
02497 NCONT=0
02498 100 CONTINUE
02499 NCONT=NCONT+1
02500
02501 PHS=1.D0
02502 PHSPAC = 1./2.**(6*ND-7) /PI**(3*ND-4)
02503 AMX=AMW
02504 CALL RANMAR(RRR,ND-2)
02505 DO 230 IL=1,ND-2
02506 AMS1=.0D0
02507 DO 231 JL=IL+1,ND
02508 231 AMS1=AMS1+AMPIK(JL,JNPI)
02509 AMS1=AMS1**2
02510 AMS2=(AMX-AMPIK(IL,JNPI))**2
02511 RR1=RRR(IL)
02512 AMX2=AMS1+ RR1*(AMS2-AMS1)
02513 AMX=SQRT(AMX2)
02514 PV(5,IL+1)=AMX
02515 PHSPAC=PHSPAC * (AMS2-AMS1)
02516
02517 PHS=PHS* (AMS2-AMS1)
02518 PA=PAWT(PV(5,IL),PV(5,IL+1),AMPIK(IL,JNPI))
02519 PHS =PHS *PA/PV(5,IL)
02520 230 CONTINUE
02521 PA=PAWT(PV(5,ND-1),AMPIK(ND-1,JNPI),AMPIK(ND,JNPI))
02522 PHS =PHS *PA/PV(5,ND-1)
02523 CALL RANMAR(RN,1)
02524 WETMAX(JNPI)=1.2D0*MAX(WETMAX(JNPI)/1.2D0,PHS/PHSMAX)
02525 IF (NCONT.EQ.500 000) THEN
02526 XNPI=0.0
02527 DO KK=1,ND
02528 XNPI=XNPI+AMPIK(KK,JNPI)
02529 ENDDO
02530 WRITE(6,*) 'ROUNDING INSTABILITY IN DPHNPI ?'
02531 WRITE(6,*) 'AMW=',AMW,'XNPI=',XNPI
02532 WRITE(6,*) 'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
02533 WRITE(6,*) 'PHS=',PHS,'PHSMAX=',PHSMAX
02534 GOTO 500
02535 ENDIF
02536 IF(RN(1)*PHSMAX*WETMAX(JNPI).GT.PHS) GO TO 100
02537
02538
02539 280 DO 300 IL=1,ND-1
02540 PA=PAWT(PV(5,IL),PV(5,IL+1),AMPIK(IL,JNPI))
02541
02542 CALL RANMAR(RRX,2)
02543 UE(3)=2.*RRX(1)-1.
02544 PHI=2.*PI*RRX(2)
02545 UE(1)=SQRT(1.D0-UE(3)**2)*COS(PHI)
02546 UE(2)=SQRT(1.D0-UE(3)**2)*SIN(PHI)
02547
02548 DO 290 J=1,3
02549 PPI(J,IL)=PA*UE(J)
02550 290 PV(J,IL+1)=-PA*UE(J)
02551 PPI(4,IL)=SQRT(PA**2+AMPIK(IL,JNPI)**2)
02552 PV(4,IL+1)=SQRT(PA**2+PV(5,IL+1)**2)
02553 PHSPAC=PHSPAC *(4.*PI)*(2.*PA/PV(5,IL))
02554 300 CONTINUE
02555
02556 DO 310 J=1,4
02557 310 PPI(J,ND)=PV(J,ND)
02558 DO 340 IL=ND-1,1,-1
02559 DO 320 J=1,3
02560 320 BE(J)=PV(J,IL)/PV(4,IL)
02561 GAM=PV(4,IL)/PV(5,IL)
02562 DO 340 I=IL,ND
02563 BEP=BE(1)*PPI(1,I)+BE(2)*PPI(2,I)+BE(3)*PPI(3,I)
02564 DO 330 J=1,3
02565
02566 330 PPI(J,I)=PPI(J,I)+GAM*(GAM*BEP/(1.D0+GAM)+PPI(4,I))*BE(J)
02567
02568 PPI(4,I)=GAM*(PPI(4,I)+BEP)
02569 340 CONTINUE
02570
02571 HV(4)=1.
02572 HV(3)=0.
02573 HV(2)=0.
02574 HV(1)=0.
02575
02576 DO K=1,4
02577 PNX(K)=PN(K)
02578 PRX(K)=PR(K)
02579 HVX(K)=HV(K)
02580 DO L=1,ND
02581 PPIX(K,L)=PPI(K,L)
02582 ENDDO
02583 ENDDO
02584
02585 RETURN
02586 END
02587 FUNCTION SIGEE(Q2,JNP)
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597
02598
02599
02600
02601
02602
02603
02604
02605
02606 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02607 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02608 * ,AMK,AMKZ,AMKST,GAMKST
02609
02610 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02611 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02612 * ,AMK,AMKZ,AMKST,GAMKST
02613 REAL*4 DATSIG(17,6)
02614
02615 DATA DATSIG/
02616 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
02617 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
02618 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
02619 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
02620 5 17*.0,
02621 6 17*.0,
02622 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
02623 8 17*.0/
02624 DATA SIG0 / 86.8 /
02625 DATA PI /3.141592653589793238462643/
02626 DATA INIT / 0 /
02627
02628
02629 JNPI=JNP
02630 IF(JNPI.GT.6) JNPI=6
02631
02632 IF(JNP.EQ.4) JNPI=3
02633 IF(JNP.EQ.3) JNPI=4
02634 IF(INIT.EQ.0) THEN
02635 INIT=1
02636
02637
02638
02639
02640 AMPI2=AMPI**2
02641 FPI = .943*AMPI
02642 DO 100 I=1,17
02643 DATSIG(I,2) = DATSIG(I,2)/2.
02644 DATSIG(I,1) = DATSIG(I,1) + DATSIG(I,2)
02645 S = 1.025+(I-1)*.05
02646 FACT=0.
02647 S2=S**2
02648 DO 200 J=1,17
02649 T= 1.025+(J-1)*.05
02650 IF(T . GT. S-AMPI ) GO TO 201
02651 T2=T**2
02652 FACT=(T2/S2)**2*SQRT((S2-T2-AMPI2)**2-4.*T2*AMPI2)/S2 *2.*T*.05
02653 FACT = FACT * (DATSIG(J,1)+DATSIG(J+1,1))
02654 200 DATSIG(I,3) = DATSIG(I,3) + FACT
02655 201 DATSIG(I,3) = DATSIG(I,3) /(2*PI*FPI)**2
02656 DATSIG(I,4) = DATSIG(I,3)
02657 DATSIG(I,6) = DATSIG(I,5)
02658 100 CONTINUE
02659
02660 1000 FORMAT(///1X,' EE SIGMA USED IN MULTIPI DECAYS'/
02661 % (17F7.2/))
02662 ENDIF
02663 Q=SQRT(Q2)
02664 QMIN=1.
02665 IF(Q.LT.QMIN) THEN
02666 SIGEE=DATSIG(1,JNPI)+
02667 & (DATSIG(2,JNPI)-DATSIG(1,JNPI))*(Q-1.)/.05
02668 ELSEIF(Q.LT.1.8) THEN
02669 DO 1 I=1,16
02670 QMAX = QMIN + .05
02671 IF(Q.LT.QMAX) GO TO 2
02672 QMIN = QMIN + .05
02673 1 CONTINUE
02674 2 SIGEE=DATSIG(I,JNPI)+
02675 & (DATSIG(I+1,JNPI)-DATSIG(I,JNPI)) * (Q-QMIN)/.05
02676 ELSEIF(Q.GT.1.8) THEN
02677 SIGEE=DATSIG(17,JNPI)+
02678 & (DATSIG(17,JNPI)-DATSIG(16,JNPI)) * (Q-1.8)/.05
02679 ENDIF
02680 IF(SIGEE.LT..0) SIGEE=0.
02681
02682 SIGEE = SIGEE/(6.*PI**2*SIG0)
02683
02684 RETURN
02685 END
02686
02687 FUNCTION SIGOLD(Q2,JNPI)
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02707 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02708 * ,AMK,AMKZ,AMKST,GAMKST
02709
02710 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02711 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02712 * ,AMK,AMKZ,AMKST,GAMKST
02713 REAL*4 DATSIG(17,4)
02714
02715 DATA DATSIG/
02716 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
02717 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
02718 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
02719 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
02720 5 17*.0,
02721 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
02722 DATA SIG0 / 86.8 /
02723 DATA PI /3.141592653589793238462643/
02724 DATA INIT / 0 /
02725
02726 IF(INIT.EQ.0) THEN
02727 INIT=1
02728 AMPI2=AMPI**2
02729 FPI = .943*AMPI
02730 DO 100 I=1,17
02731 DATSIG(I,2) = DATSIG(I,2)/2.
02732 DATSIG(I,1) = DATSIG(I,1) + DATSIG(I,2)
02733 S = 1.025+(I-1)*.05
02734 FACT=0.
02735 S2=S**2
02736 DO 200 J=1,17
02737 T= 1.025+(J-1)*.05
02738 IF(T . GT. S-AMPI ) GO TO 201
02739 T2=T**2
02740 FACT=(T2/S2)**2*SQRT((S2-T2-AMPI2)**2-4.*T2*AMPI2)/S2 *2.*T*.05
02741 FACT = FACT * (DATSIG(J,1)+DATSIG(J+1,1))
02742 200 DATSIG(I,3) = DATSIG(I,3) + FACT
02743 201 DATSIG(I,3) = DATSIG(I,3) /(2*PI*FPI)**2
02744 100 CONTINUE
02745
02746 1000 FORMAT(///1X,' EE SIGMA USED IN MULTIPI DECAYS'/
02747 % (17F7.2/))
02748 ENDIF
02749 Q=SQRT(Q2)
02750 QMIN=1.
02751 IF(Q.LT.QMIN) THEN
02752 SIGEE=DATSIG(1,JNPI)+
02753 & (DATSIG(2,JNPI)-DATSIG(1,JNPI))*(Q-1.)/.05
02754 ELSEIF(Q.LT.1.8) THEN
02755 DO 1 I=1,16
02756 QMAX = QMIN + .05
02757 IF(Q.LT.QMAX) GO TO 2
02758 QMIN = QMIN + .05
02759 1 CONTINUE
02760 2 SIGEE=DATSIG(I,JNPI)+
02761 & (DATSIG(I+1,JNPI)-DATSIG(I,JNPI)) * (Q-QMIN)/.05
02762 ELSEIF(Q.GT.1.8) THEN
02763 SIGEE=DATSIG(17,JNPI)+
02764 & (DATSIG(17,JNPI)-DATSIG(16,JNPI)) * (Q-1.8)/.05
02765 ENDIF
02766 IF(SIGEE.LT..0) SIGEE=0.
02767
02768 SIGEE = SIGEE/(6.*PI**2*SIG0)
02769 SIGOLD=SIGEE
02770
02771 RETURN
02772 END
02773 SUBROUTINE DPHSPK(DGAMT,HV,PN,PAA,PNPI,JAA)
02774
02775
02776
02777
02778 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
02779
02780 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
02781
02782 & ,NAMES
02783 CHARACTER NAMES(NMODE)*31
02784
02785 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PNPI(4,9)
02786
02787 MNUM=JAA
02788
02789 KEYT=4
02790 IF(JAA.EQ.7) KEYT=3
02791 IF(JAA.EQ.9) KEYT=1
02792
02793 AMP1=DCDMAS(IDFFIN(1,JAA+NM4+NM5+NM6))
02794 AMP2=DCDMAS(IDFFIN(2,JAA+NM4+NM5+NM6))
02795 AMP3=DCDMAS(IDFFIN(3,JAA+NM4+NM5+NM6))
02796 CALL
02797 $ DPHTRE(DGAMT,HV,PN,PAA,PIM1,AMP1,PIM2,AMP2,PIPL,AMP3,KEYT,MNUM)
02798 DO I=1,4
02799 PNPI(I,1)=PIM1(I)
02800 PNPI(I,2)=PIM2(I)
02801 PNPI(I,3)=PIPL(I)
02802 ENDDO
02803 END
02804
02805
02806
02807
02808 SUBROUTINE
02809 $ DPHTRE(DGAMT,HV,PN,PAA,PIM1,AMPA,PIM2,AMPB,PIPL,AMP3,KEYT,MNUM)
02810
02811
02812
02813
02814
02815
02816
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02827 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02828 * ,AMK,AMKZ,AMKST,GAMKST
02829
02830 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
02831 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
02832 * ,AMK,AMKZ,AMKST,GAMKST
02833 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02834 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
02835 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
02836 REAL PR(4)
02837 REAL*4 RRR(5)
02838 DATA PI /3.141592653589793238462643/
02839 DATA ICONT /0/
02840 XLAM(X,Y,Z)=SQRT(ABS((X-Y-Z)**2-4.0*Y*Z))
02841
02842
02843
02844
02845 PHSPAC=1./2**17/PI**8
02846
02847 PT(1)=0.
02848 PT(2)=0.
02849 PT(3)=0.
02850 PT(4)=AMTAU
02851
02852 CALL RANMAR(RRR,5)
02853 RR=RRR(5)
02854
02855 CALL CHOICE(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
02856 $ AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
02857 IF (ICHAN.EQ.1) THEN
02858 AMP1=AMPB
02859 AMP2=AMPA
02860 ELSEIF (ICHAN.EQ.2) THEN
02861 AMP1=AMPA
02862 AMP2=AMPB
02863 ELSE
02864 AMP1=AMPB
02865 AMP2=AMPA
02866 ENDIF
02867
02868 RR1=RRR(1)
02869 AMS1=(AMP1+AMP2+AMP3)**2
02870 AMS2=(AMTAU-AMNUTA)**2
02871
02872
02873
02874 ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX)
02875 ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX)
02876 ALP=ALP1+RR1*(ALP2-ALP1)
02877 AM3SQ =AMRX**2+AMRX*GAMRX*TAN(ALP)
02878 AM3 =SQRT(AM3SQ)
02879 PHSPAC=PHSPAC*((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX)
02880 PHSPAC=PHSPAC*(ALP2-ALP1)
02881
02882 RR2=RRR(2)
02883 AMS1=(AMP2+AMP3)**2
02884 AMS2=(AM3-AMP1)**2
02885 IF (ICHAN.LE.2) THEN
02886
02887
02888
02889 ALP1=ATAN((AMS1-AMRA**2)/AMRA/GAMRA)
02890 ALP2=ATAN((AMS2-AMRA**2)/AMRA/GAMRA)
02891 ALP=ALP1+RR2*(ALP2-ALP1)
02892 AM2SQ =AMRA**2+AMRA*GAMRA*TAN(ALP)
02893 AM2 =SQRT(AM2SQ)
02894
02895
02896
02897
02898 ELSE
02899
02900
02901
02902 AM2SQ=AMS1+ RR2*(AMS2-AMS1)
02903 AM2 =SQRT(AM2SQ)
02904 PHF0=(AMS2-AMS1)
02905 ENDIF
02906
02907
02908
02909 ENQ1=(AM2SQ-AMP2**2+AMP3**2)/(2*AM2)
02910 ENQ2=(AM2SQ+AMP2**2-AMP3**2)/(2*AM2)
02911 PPI= ENQ1**2-AMP3**2
02912 PPPI=SQRT(ABS(ENQ1**2-AMP3**2))
02913
02914 PHF1=(4*PI)*(2*PPPI/AM2)
02915
02916
02917
02918 CALL SPHERA(PPPI,PIPL)
02919 PIPL(4)=ENQ1
02920
02921
02922
02923 DO 30 I=1,3
02924 30 PIM1(I)=-PIPL(I)
02925 PIM1(4)=ENQ2
02926
02927
02928
02929
02930 PR(1)=0
02931 PR(2)=0
02932 PR(4)=1./(2*AM3)*(AM3**2+AM2**2-AMP1**2)
02933 PR(3)= SQRT(ABS(PR(4)**2-AM2**2))
02934 PPI = PR(4)**2-AM2**2
02935
02936 PIM2(1)=0
02937 PIM2(2)=0
02938 PIM2(4)=1./(2*AM3)*(AM3**2-AM2**2+AMP1**2)
02939 PIM2(3)=-PR(3)
02940 PHF2=(4*PI)*(2*PR(3)/AM3)
02941
02942
02943
02944 EXE=(PR(4)+PR(3))/AM2
02945 CALL BOSTR3(EXE,PIPL,PIPL)
02946 CALL BOSTR3(EXE,PIM1,PIM1)
02947 RR3=RRR(3)
02948 RR4=RRR(4)
02949
02950
02951
02952 THET =ACOS(-1.+2*RR3)
02953 PHI = 2*PI*RR4
02954 CALL ROTPOL(THET,PHI,PIPL)
02955 CALL ROTPOL(THET,PHI,PIM1)
02956 CALL ROTPOL(THET,PHI,PIM2)
02957 CALL ROTPOL(THET,PHI,PR)
02958
02959
02960
02961 PAA(1)=0
02962 PAA(2)=0
02963 PAA(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AM3**2)
02964 PAA(3)= SQRT(ABS(PAA(4)**2-AM3**2))
02965 PPI = PAA(4)**2-AM3**2
02966 PHSPAC=PHSPAC*(4*PI)*(2*PAA(3)/AMTAU)
02967
02968 PN(1)=0
02969 PN(2)=0
02970 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AM3**2)
02971 PN(3)=-PAA(3)
02972
02973
02974 AMS1=(AMP2+AMP3)**2
02975 AMS2=(AM3-AMP1)**2
02976 ALP1=ATAN((AMS1-AMRA**2)/AMRA/GAMRA)
02977 ALP2=ATAN((AMS2-AMRA**2)/AMRA/GAMRA)
02978 XPRO = (PIM1(3)+PIPL(3))**2
02979 $ +(PIM1(2)+PIPL(2))**2+(PIM1(1)+PIPL(1))**2
02980 AM2SQ=-XPRO+(PIM1(4)+PIPL(4))**2
02981
02982 FF1 = ((AM2SQ-AMRA**2)**2+(AMRA*GAMRA)**2)/(AMRA*GAMRA)
02983 FF1 =FF1 *(ALP2-ALP1)
02984
02985 GG1 = (4*PI)*(XLAM(AM2SQ,AMP2**2,AMP3**2)/AM2SQ)
02986
02987 GG1 =GG1 *(4*PI)*SQRT(4*XPRO/AM3SQ)
02988 XJAJE=GG1*(AMS2-AMS1)
02989
02990 AMS1=(AMP1+AMP3)**2
02991 AMS2=(AM3-AMP2)**2
02992 ALP1=ATAN((AMS1-AMRB**2)/AMRB/GAMRB)
02993 ALP2=ATAN((AMS2-AMRB**2)/AMRB/GAMRB)
02994 XPRO = (PIM2(3)+PIPL(3))**2
02995 $ +(PIM2(2)+PIPL(2))**2+(PIM2(1)+PIPL(1))**2
02996 AM2SQ=-XPRO+(PIM2(4)+PIPL(4))**2
02997 FF2 = ((AM2SQ-AMRB**2)**2+(AMRB*GAMRB)**2)/(AMRB*GAMRB)
02998 FF2 =FF2 *(ALP2-ALP1)
02999 GG2 = (4*PI)*(XLAM(AM2SQ,AMP1**2,AMP3**2)/AM2SQ)
03000 GG2 =GG2 *(4*PI)*SQRT(4*XPRO/AM3SQ)
03001 XJADW=GG2*(AMS2-AMS1)
03002
03003 A1=0.0
03004 A2=0.0
03005 A3=0.0
03006 XJAC1=FF1*GG1
03007 XJAC2=FF2*GG2
03008 IF (ICHAN.EQ.2) THEN
03009 XJAC3=XJADW
03010 ELSE
03011 XJAC3=XJAJE
03012 ENDIF
03013 IF (XJAC1.NE.0.0) A1=PROB1/XJAC1
03014 IF (XJAC2.NE.0.0) A2=PROB2/XJAC2
03015 IF (XJAC3.NE.0.0) A3=PROB3/XJAC3
03016
03017 IF (A1+A2+A3.NE.0.0) THEN
03018 PHSPAC=PHSPAC/(A1+A2+A3)
03019 ELSE
03020 PHSPAC=0.0
03021 ENDIF
03022 IF(ICHAN.EQ.2) THEN
03023 DO 70 I=1,4
03024 X=PIM1(I)
03025 PIM1(I)=PIM2(I)
03026 70 PIM2(I)=X
03027 ENDIF
03028
03029
03030 EXE=(PAA(4)+PAA(3))/AM3
03031 CALL BOSTR3(EXE,PIPL,PIPL)
03032 CALL BOSTR3(EXE,PIM1,PIM1)
03033 CALL BOSTR3(EXE,PIM2,PIM2)
03034 CALL BOSTR3(EXE,PR,PR)
03035
03036 IF (MNUM.EQ.8) THEN
03037 CALL DAMPOG(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
03038
03039
03040 ELSE
03041
03042 CALL DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
03043
03044
03045 ENDIF
03046
03047 IF (KEYT.EQ.1.OR.KEYT.EQ.2) THEN
03048
03049
03050
03051 PHSPAC=PHSPAC/2.
03052
03053 ENDIF
03054 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
03055 END
03056 SUBROUTINE DAMPAA(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
03057
03058
03059
03060
03061
03062
03063
03064
03065
03066 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03067 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03068 * ,AMK,AMKZ,AMKST,GAMKST
03069
03070 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03071 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03072 * ,AMK,AMKZ,AMKST,GAMKST
03073 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03074 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03075 COMMON /TESTA1/ KEYA1
03076 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
03077 REAL PAA(4),VEC1(4),VEC2(4)
03078 REAL PIVEC(4),PIAKS(4),HVM(4)
03079 COMPLEX BWIGN,HADCUR(4),FPIK
03080 DATA ICONT /1/
03081
03082
03083
03084 DATA FPI /93.3E-3/
03085
03086 BWIGN(XM,AM,GAMMA)=1./CMPLX(XM**2-AM**2,GAMMA*AM)
03087
03088
03089 DO 10 I=1,4
03090 10 PAA(I)=PIM1(I)+PIM2(I)+PIPL(I)
03091
03092 XMAA =SQRT(ABS(PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2))
03093 XMRO1 =SQRT(ABS((PIPL(4)+PIM1(4))**2-(PIPL(1)+PIM1(1))**2
03094 $ -(PIPL(2)+PIM1(2))**2-(PIPL(3)+PIM1(3))**2))
03095 XMRO2 =SQRT(ABS((PIPL(4)+PIM2(4))**2-(PIPL(1)+PIM2(1))**2
03096 $ -(PIPL(2)+PIM2(2))**2-(PIPL(3)+PIM2(3))**2))
03097
03098 PROD1 =PAA(4)*(PIM1(4)-PIPL(4))-PAA(1)*(PIM1(1)-PIPL(1))
03099 $ -PAA(2)*(PIM1(2)-PIPL(2))-PAA(3)*(PIM1(3)-PIPL(3))
03100 PROD2 =PAA(4)*(PIM2(4)-PIPL(4))-PAA(1)*(PIM2(1)-PIPL(1))
03101 $ -PAA(2)*(PIM2(2)-PIPL(2))-PAA(3)*(PIM2(3)-PIPL(3))
03102 DO 40 I=1,4
03103 VEC1(I)= PIM1(I)-PIPL(I) -PAA(I)*PROD1/XMAA**2
03104 40 VEC2(I)= PIM2(I)-PIPL(I) -PAA(I)*PROD2/XMAA**2
03105
03106 IF (KEYA1.EQ.1) THEN
03107 FA1=9.87
03108 FAROPI=1.0
03109 FRO2PI=1.0
03110 FNORM=FA1/SQRT(2.)*FAROPI*FRO2PI
03111 DO 45 I=1,4
03112 HADCUR(I)= CMPLX(FNORM) *AMA1**2*BWIGN(XMAA,AMA1,GAMA1)
03113 $ *(CMPLX(VEC1(I))*AMRO**2*BWIGN(XMRO1,AMRO,GAMRO)
03114 $ +CMPLX(VEC2(I))*AMRO**2*BWIGN(XMRO2,AMRO,GAMRO))
03115 45 CONTINUE
03116 ELSE
03117 FNORM=2.0*SQRT(2.)/3.0/FPI
03118 GAMAX=GAMA1*GFUN(XMAA**2)/GFUN(AMA1**2)
03119 DO 46 I=1,4
03120 HADCUR(I)= CMPLX(FNORM) *AMA1**2*BWIGN(XMAA,AMA1,GAMAX)
03121 $ *(CMPLX(VEC1(I))*FPIK(XMRO1)
03122 $ +CMPLX(VEC2(I))*FPIK(XMRO2))
03123 46 CONTINUE
03124 ENDIF
03125
03126
03127 CALL CLVEC(HADCUR,PN,PIVEC)
03128 CALL CLAXI(HADCUR,PN,PIAKS)
03129 CALL CLNUT(HADCUR,BRAKM,HVM)
03130
03131 BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4)
03132 & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM
03133 AMPLIT=(GFERMI*CCABIB)**2*BRAK/2.
03134
03135
03136
03137 DO 90 I=1,3
03138 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I)))
03139 & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I)
03140
03141 HV(I)=-HV(I)/BRAK
03142 90 CONTINUE
03143 END
03144
03145 FUNCTION GFUN(QKWA)
03146
03147
03148
03149 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03150 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03151 * ,AMK,AMKZ,AMKST,GAMKST
03152
03153 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03154 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03155 * ,AMK,AMKZ,AMKST,GAMKST
03156
03157 IF (QKWA.LT.(AMRO+AMPI)**2) THEN
03158 GFUN=4.1*(QKWA-9*AMPIZ**2)**3
03159 $ *(1.-3.3*(QKWA-9*AMPIZ**2)+5.8*(QKWA-9*AMPIZ**2)**2)
03160 ELSE
03161 GFUN=QKWA*(1.623+10.38/QKWA-9.32/QKWA**2+0.65/QKWA**3)
03162 ENDIF
03163 END
03164 COMPLEX FUNCTION BWIGS(S,M,G)
03165
03166
03167
03168 REAL S,M,G
03169 REAL PI,PIM,QS,QM,W,GS,MK
03170
03171
03172 REAL PM, PG, PBETA
03173 COMPLEX BW,BWP
03174
03175 DATA INIT /0/
03176 P(A,B,C)=SQRT(ABS(ABS(((A+B-C)**2-4.*A*B)/4./A)
03177 $ +(((A+B-C)**2-4.*A*B)/4./A))/2.0)
03178
03179 IF (INIT.EQ.0) THEN
03180 INIT=1
03181 PI=3.141592654
03182 PIM=.139
03183 MK=.493667
03184
03185
03186 PM = PKORB(1,16)
03187 PG = PKORB(2,16)
03188 PBETA = PKORB(3,16)
03189
03190
03191 ENDIF
03192
03193 QS=P(S,PIM**2,MK**2)
03194 QM=P(M**2,PIM**2,MK**2)
03195 W=SQRT(S)
03196 GS=G*(M/W)*(QS/QM)**3
03197
03198 BW=M**2/CMPLX(M**2-S,-M*GS)
03199 QPM=P(PM**2,PIM**2,MK**2)
03200 G1=PG*(PM/W)*(QS/QPM)**3
03201 BWP=PM**2/CMPLX(PM**2-S,-PM*G1)
03202 BWIGS= (BW+PBETA*BWP)/(1+PBETA)
03203
03204 RETURN
03205 END
03206 COMPLEX FUNCTION BWIG(S,M,G)
03207
03208
03209
03210 REAL S,M,G
03211 REAL PI,PIM,QS,QM,W,GS
03212 DATA INIT /0/
03213
03214 IF (INIT.EQ.0) THEN
03215 INIT=1
03216 PI=3.141592654
03217 PIM=.139
03218
03219 ENDIF
03220 IF (S.GT.4.*PIM**2) THEN
03221 QS=SQRT(ABS(ABS(S/4.-PIM**2)+(S/4.-PIM**2))/2.0)
03222 QM=SQRT(M**2/4.-PIM**2)
03223 W=SQRT(S)
03224 GS=G*(M/W)*(QS/QM)**3
03225 ELSE
03226 GS=0.0
03227 ENDIF
03228 BWIG=M**2/CMPLX(M**2-S,-M*GS)
03229 RETURN
03230 END
03231 COMPLEX FUNCTION FPIK(W)
03232
03233
03234
03235 COMPLEX BWIG
03236 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
03237 EXTERNAL BWIG
03238 DATA INIT /0/
03239
03240
03241 IF (INIT.EQ.0 ) THEN
03242 INIT=1
03243 PI=3.141592654
03244 PIM=.140
03245
03246 ROM=PKORB(1,9)
03247 ROG=PKORB(2,9)
03248 ROM1=PKORB(1,15)
03249 ROG1=PKORB(2,15)
03250 BETA1=PKORB(3,15)
03251
03252 ENDIF
03253
03254 S=W**2
03255 FPIK= (BWIG(S,ROM,ROG)+BETA1*BWIG(S,ROM1,ROG1))
03256 & /(1+BETA1)
03257 RETURN
03258 END
03259 FUNCTION FPIRHO(W)
03260
03261
03262
03263 COMPLEX FPIK
03264 FPIRHO=CABS(FPIK(W))**2
03265 END
03266 SUBROUTINE CLVEC(HJ,PN,PIV)
03267
03268
03269
03270
03271
03272
03273 REAL PIV(4),PN(4)
03274 COMPLEX HJ(4),HN
03275
03276 HN= HJ(4)*CMPLX(PN(4))-HJ(3)*CMPLX(PN(3))
03277 HH= REAL(HJ(4)*CONJG(HJ(4))-HJ(3)*CONJG(HJ(3))
03278 $ -HJ(2)*CONJG(HJ(2))-HJ(1)*CONJG(HJ(1)))
03279 DO 10 I=1,4
03280 10 PIV(I)=4.*REAL(HN*CONJG(HJ(I)))-2.*HH*PN(I)
03281 RETURN
03282 END
03283 SUBROUTINE CLAXI(HJ,PN,PIA)
03284
03285
03286
03287
03288
03289
03290 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
03291 COMMON / IDFC / IDFF
03292 REAL PIA(4),PN(4)
03293 COMPLEX HJ(4),HJC(4)
03294
03295
03296 DET2(I,J)=AIMAG(HJC(I)*HJ(J)-HJC(J)*HJ(I))
03297
03298
03299
03300
03301 IF (KTOM.EQ.1.OR.KTOM.EQ.-1) THEN
03302 SIGN= IDFF/ABS(IDFF)
03303 ELSEIF (KTOM.EQ.2) THEN
03304 SIGN=-IDFF/ABS(IDFF)
03305 ELSE
03306 PRINT *, 'STOP IN CLAXI: KTOM=',KTOM
03307 STOP
03308 ENDIF
03309
03310 DO 10 I=1,4
03311 10 HJC(I)=CONJG(HJ(I))
03312 PIA(1)= -2.*PN(3)*DET2(2,4)+2.*PN(4)*DET2(2,3)
03313 PIA(2)= -2.*PN(4)*DET2(1,3)+2.*PN(3)*DET2(1,4)
03314 PIA(3)= 2.*PN(4)*DET2(1,2)
03315 PIA(4)= 2.*PN(3)*DET2(1,2)
03316
03317 DO 20 I=1,4
03318 20 PIA(I)=PIA(I)*SIGN
03319 END
03320 SUBROUTINE CLNUT(HJ,B,HV)
03321
03322
03323
03324
03325
03326
03327 COMPLEX HJ(4)
03328 REAL HV(4),P(4)
03329 DATA P /3*0.,1.0/
03330
03331 CALL CLAXI(HJ,P,HV)
03332 B=REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3))
03333 & - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) )
03334 RETURN
03335 END
03336 SUBROUTINE DAMPOG(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
03337
03338
03339
03340
03341
03342
03343
03344
03345
03346
03347
03348 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03349 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03350 * ,AMK,AMKZ,AMKST,GAMKST
03351
03352 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03353 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03354 * ,AMK,AMKZ,AMKST,GAMKST
03355 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03356 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03357 COMMON /TESTA1/ KEYA1
03358 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
03359 REAL PAA(4),VEC1(4),VEC2(4)
03360 REAL PIVEC(4),PIAKS(4),HVM(4)
03361 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
03362 DATA ICONT /1/
03363
03364
03365
03366
03367
03368
03369 DO 10 I=1,4
03370 VEC1(I)=0.0
03371 VEC2(I)=0.0
03372 HV(I) =0.0
03373 10 PAA(I)=PIM1(I)+PIM2(I)+PIPL(I)
03374 VEC1(1)=1.0
03375
03376 XMAA =SQRT(ABS(PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2))
03377 XMOM =SQRT(ABS( (PIM2(4)+PIPL(4))**2-(PIM2(3)+PIPL(3))**2
03378 $ -(PIM2(2)+PIPL(2))**2-(PIM2(1)+PIPL(1))**2 ))
03379 XMRO2 =(PIPL(1))**2 +(PIPL(2))**2 +(PIPL(3))**2
03380
03381 PROD1 =VEC1(1)*PIPL(1)
03382 PROD2 =VEC2(2)*PIPL(2)
03383 P12 =PIM1(4)*PIM2(4)-PIM1(1)*PIM2(1)
03384 $ -PIM1(2)*PIM2(2)-PIM1(3)*PIM2(3)
03385 P1PL =PIM1(4)*PIPL(4)-PIM1(1)*PIPL(1)
03386 $ -PIM1(2)*PIPL(2)-PIM1(3)*PIPL(3)
03387 P2PL =PIPL(4)*PIM2(4)-PIPL(1)*PIM2(1)
03388 $ -PIPL(2)*PIM2(2)-PIPL(3)*PIM2(3)
03389 DO 40 I=1,3
03390 VEC1(I)= (VEC1(I)-PROD1/XMRO2*PIPL(I))
03391 40 CONTINUE
03392 GNORM=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
03393 DO 41 I=1,3
03394 VEC1(I)= VEC1(I)/GNORM
03395 41 CONTINUE
03396 VEC2(1)=(VEC1(2)*PIPL(3)-VEC1(3)*PIPL(2))/SQRT(XMRO2)
03397 VEC2(2)=(VEC1(3)*PIPL(1)-VEC1(1)*PIPL(3))/SQRT(XMRO2)
03398 VEC2(3)=(VEC1(1)*PIPL(2)-VEC1(2)*PIPL(1))/SQRT(XMRO2)
03399 P1VEC1 =PIM1(4)*VEC1(4)-PIM1(1)*VEC1(1)
03400 $ -PIM1(2)*VEC1(2)-PIM1(3)*VEC1(3)
03401 P2VEC1 =VEC1(4)*PIM2(4)-VEC1(1)*PIM2(1)
03402 $ -VEC1(2)*PIM2(2)-VEC1(3)*PIM2(3)
03403 P1VEC2 =PIM1(4)*VEC2(4)-PIM1(1)*VEC2(1)
03404 $ -PIM1(2)*VEC2(2)-PIM1(3)*VEC2(3)
03405 P2VEC2 =VEC2(4)*PIM2(4)-VEC2(1)*PIM2(1)
03406 $ -VEC2(2)*PIM2(2)-VEC2(3)*PIM2(3)
03407
03408 FNORM=FORMOM(XMAA,XMOM)
03409 BRAK=0.0
03410 DO 120 JJ=1,2
03411 DO 45 I=1,4
03412 IF (JJ.EQ.1) THEN
03413 HADCUR(I) = FNORM *(
03414 $ VEC1(I)*(AMPI**2*P1PL-P2PL*(P12-P1PL))
03415 $ -PIM2(I)*(P2VEC1*P1PL-P1VEC1*P2PL)
03416 $ +PIPL(I)*(P2VEC1*P12 -P1VEC1*(AMPI**2+P2PL)) )
03417 ELSE
03418 HADCUR(I) = FNORM *(
03419 $ VEC2(I)*(AMPI**2*P1PL-P2PL*(P12-P1PL))
03420 $ -PIM2(I)*(P2VEC2*P1PL-P1VEC2*P2PL)
03421 $ +PIPL(I)*(P2VEC2*P12 -P1VEC2*(AMPI**2+P2PL)) )
03422 ENDIF
03423 45 CONTINUE
03424
03425
03426 CALL CLVEC(HADCUR,PN,PIVEC)
03427 CALL CLAXI(HADCUR,PN,PIAKS)
03428 CALL CLNUT(HADCUR,BRAKM,HVM)
03429
03430 BRAK=BRAK+(GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4)
03431 & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM
03432 DO 90 I=1,3
03433 HV(I)=HV(I)-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I)))
03434 & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I)
03435 90 CONTINUE
03436
03437 120 CONTINUE
03438 AMPLIT=(GFERMI*CCABIB)**2*BRAK/2.
03439
03440
03441
03442 DO 91 I=1,3
03443 HV(I)=-HV(I)/BRAK
03444 91 CONTINUE
03445
03446 END
03447 SUBROUTINE DAMPPK(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
03448
03449
03450
03451
03452
03453
03454
03455
03456
03457
03458
03459 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03460 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03461 * ,AMK,AMKZ,AMKST,GAMKST
03462
03463 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03464 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03465 * ,AMK,AMKZ,AMKST,GAMKST
03466 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03467 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03468 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
03469 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
03470 REAL PIVEC(4),PIAKS(4),HVM(4)
03471 REAL FNORM(0:19),COEF(1:5,0:19)
03472 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
03473
03474 COMPLEX F1,F2,F3,F4,F5
03475
03476 EXTERNAL FORM1,FORM2,FORM3,FORM4,FORM5
03477 DATA PI /3.141592653589793238462643/
03478 DATA ICONT /0/
03479
03480 DATA FPI /93.3E-3/
03481 IF (ICONT.EQ.0) THEN
03482 ICONT=1
03483 UROJ=CMPLX(0.0,1.0)
03484 DWAPI0=SQRT(2.0)
03485 FNORM(0)=CCABIB/FPI
03486 FNORM(1)=CCABIB/FPI
03487 FNORM(2)=CCABIB/FPI
03488 FNORM(3)=CCABIB/FPI
03489 FNORM(4)=SCABIB/FPI/DWAPI0
03490 FNORM(5)=SCABIB/FPI
03491 FNORM(6)=SCABIB/FPI
03492 FNORM(7)=CCABIB/FPI
03493 FNORM(8)=0.0
03494 FNORM(9)=CCABIB/FPI
03495 DO K=10,19
03496 FNORM(K)=FNORM(9)
03497 ENDDO
03498
03499 COEF(1,0)= 2.0*SQRT(2.)/3.0
03500 COEF(2,0)=-2.0*SQRT(2.)/3.0
03501
03502
03503 COEF(3,0)= 2.0*SQRT(2.)/3.0
03504
03505 COEF(4,0)= FPI
03506 COEF(5,0)= 0.0
03507
03508 COEF(1,1)=-SQRT(2.)/3.0
03509 COEF(2,1)= SQRT(2.)/3.0
03510 COEF(3,1)= 0.0
03511 COEF(4,1)= FPI
03512 COEF(5,1)= SQRT(2.)
03513
03514 COEF(1,2)=-SQRT(2.)/3.0
03515 COEF(2,2)= SQRT(2.)/3.0
03516 COEF(3,2)= 0.0
03517 COEF(4,2)= 0.0
03518 COEF(5,2)=-SQRT(2.)
03519
03520
03521
03522 COEF(1,3)= 1./3.
03523 COEF(2,3)=-2./3.
03524 COEF(3,3)= 2./3.
03525
03526 COEF(4,3)= 0.0
03527 COEF(5,3)= 0.0
03528
03529 COEF(1,4)= 1.0/SQRT(2.)/3.0
03530 COEF(2,4)=-1.0/SQRT(2.)/3.0
03531 COEF(3,4)= 0.0
03532 COEF(4,4)= 0.0
03533 COEF(5,4)= 0.0
03534
03535 COEF(1,5)=-SQRT(2.)/3.0
03536 COEF(2,5)= SQRT(2.)/3.0
03537 COEF(3,5)= 0.0
03538 COEF(4,5)= 0.0
03539 COEF(5,5)=-SQRT(2.)
03540
03541
03542
03543 COEF(1,6)= 1./3.
03544 COEF(2,6)=-2./3.
03545 COEF(3,6)= 2./3.
03546
03547 COEF(4,6)= 0.0
03548 COEF(5,6)=-2.0
03549
03550 COEF(1,7)= 0.0
03551 COEF(2,7)= 0.0
03552 COEF(3,7)= 0.0
03553 COEF(4,7)= 0.0
03554 COEF(5,7)=-SQRT(2.0/3.0)
03555
03556 COEF(1,8)= 0.0
03557 COEF(2,8)= 0.0
03558 COEF(3,8)= 0.0
03559 COEF(4,8)= 0.0
03560 COEF(5,8)= 0.0
03561
03562
03563 COEF(1,9)= 2.0*SQRT(2.)/3.0
03564 COEF(2,9)=-2.0*SQRT(2.)/3.0
03565 COEF(3,9)= 2.0*SQRT(2.)/3.0
03566
03567 COEF(4,9)= FPI
03568 COEF(5,9)= 0.0
03569
03570 DO K=10,19
03571 COEF(1,K)= 2.0*SQRT(2.)/3.0
03572 COEF(2,K)=-2.0*SQRT(2.)/3.0
03573 COEF(3,K)= 2.0*SQRT(2.)/3.0
03574
03575 COEF(4,K)= FPI
03576 COEF(5,K)= 0.0
03577
03578 ENDDO
03579
03580
03581 ENDIF
03582
03583 DO 10 I=1,4
03584 10 PAA(I)=PIM1(I)+PIM2(I)+PIM3(I)
03585 XMAA =SQRT(ABS(PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2))
03586 XMRO1 =SQRT(ABS((PIM3(4)+PIM2(4))**2-(PIM3(1)+PIM2(1))**2
03587 $ -(PIM3(2)+PIM2(2))**2-(PIM3(3)+PIM2(3))**2))
03588 XMRO2 =SQRT(ABS((PIM3(4)+PIM1(4))**2-(PIM3(1)+PIM1(1))**2
03589 $ -(PIM3(2)+PIM1(2))**2-(PIM3(3)+PIM1(3))**2))
03590 XMRO3 =SQRT(ABS((PIM1(4)+PIM2(4))**2-(PIM1(1)+PIM2(1))**2
03591 $ -(PIM1(2)+PIM2(2))**2-(PIM1(3)+PIM2(3))**2))
03592
03593 PROD1 =PAA(4)*(PIM2(4)-PIM3(4))-PAA(1)*(PIM2(1)-PIM3(1))
03594 $ -PAA(2)*(PIM2(2)-PIM3(2))-PAA(3)*(PIM2(3)-PIM3(3))
03595 PROD2 =PAA(4)*(PIM3(4)-PIM1(4))-PAA(1)*(PIM3(1)-PIM1(1))
03596 $ -PAA(2)*(PIM3(2)-PIM1(2))-PAA(3)*(PIM3(3)-PIM1(3))
03597 PROD3 =PAA(4)*(PIM1(4)-PIM2(4))-PAA(1)*(PIM1(1)-PIM2(1))
03598 $ -PAA(2)*(PIM1(2)-PIM2(2))-PAA(3)*(PIM1(3)-PIM2(3))
03599 DO 40 I=1,4
03600 VEC1(I)= PIM2(I)-PIM3(I) -PAA(I)*PROD1/XMAA**2
03601 VEC2(I)= PIM3(I)-PIM1(I) -PAA(I)*PROD2/XMAA**2
03602 VEC3(I)= PIM1(I)-PIM2(I) -PAA(I)*PROD3/XMAA**2
03603 40 VEC4(I)= PIM1(I)+PIM2(I)+PIM3(I)
03604 CALL PROD5(PIM1,PIM2,PIM3,VEC5)
03605
03606
03607
03608
03609 F1 = CMPLX(COEF(1,MNUM))*FORM1(MNUM,XMAA**2,XMRO1**2,XMRO2**2)
03610 F2 = CMPLX(COEF(2,MNUM))*FORM2(MNUM,XMAA**2,XMRO2**2,XMRO1**2)
03611 F3 = CMPLX(COEF(3,MNUM))*FORM3(MNUM,XMAA**2,XMRO3**2,XMRO1**2)
03612 F4 = (-1.0*UROJ)*
03613 $CMPLX(COEF(4,MNUM))*FORM4(MNUM,XMAA**2,XMRO1**2,XMRO2**2,XMRO3**2)
03614 F5 = (-1.0)*UROJ/4.0/PI**2/FPI**2*
03615 $ CMPLX(COEF(5,MNUM))*FORM5(MNUM,XMAA**2,XMRO1**2,XMRO2**2)
03616
03617 DO 45 I=1,4
03618 HADCUR(I)= CMPLX(FNORM(MNUM)) * (
03619 $ CMPLX(VEC1(I))*F1+CMPLX(VEC2(I))*F2+CMPLX(VEC3(I))*F3+
03620 $ CMPLX(VEC4(I))*F4+CMPLX(VEC5(I))*F5)
03621 45 CONTINUE
03622
03623
03624
03625 CALL CLVEC(HADCUR,PN,PIVEC)
03626 CALL CLAXI(HADCUR,PN,PIAKS)
03627 CALL CLNUT(HADCUR,BRAKM,HVM)
03628
03629 BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4)
03630 & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM
03631 AMPLIT=(GFERMI)**2*BRAK/2.
03632 IF (MNUM.GE.20) THEN
03633 PRINT *, 'MNUM=',MNUM
03634 ZNAK=-1.0
03635 XM1=0.0
03636 XM2=0.0
03637 XM3=0.0
03638 DO 77 K=1,4
03639 IF (K.EQ.4) ZNAK=1.0
03640 XM1=ZNAK*PIM1(K)**2+XM1
03641 XM2=ZNAK*PIM2(K)**2+XM2
03642 XM3=ZNAK*PIM3(K)**2+XM3
03643 77 PRINT *, 'PIM1=',PIM1(K),'PIM2=',PIM2(K),'PIM3=',PIM3(K)
03644 PRINT *, 'XM1=',SQRT(XM1),'XM2=',SQRT(XM2),'XM3=',SQRT(XM3)
03645 PRINT *, '************************************************'
03646 ENDIF
03647
03648 DO 90 I=1,3
03649 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I)))
03650 & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I)
03651
03652 HV(I)=-HV(I)/BRAK
03653 90 CONTINUE
03654 END
03655 SUBROUTINE PROD5(P1,P2,P3,PIA)
03656
03657
03658
03659
03660
03661 COMMON / JAKI / JAK1,JAK2,JAKP,JAKM,KTOM
03662 COMMON / IDFC / IDFF
03663 REAL PIA(4),P1(4),P2(4),P3(4)
03664 DET2(I,J)=P1(I)*P2(J)-P2(I)*P1(J)
03665
03666 IF (KTOM.EQ.1.OR.KTOM.EQ.-1) THEN
03667 SIGN= IDFF/ABS(IDFF)
03668 ELSEIF (KTOM.EQ.2) THEN
03669 SIGN=-IDFF/ABS(IDFF)
03670 ELSE
03671 PRINT *, 'STOP IN PROD5: KTOM=',KTOM
03672 STOP
03673 ENDIF
03674
03675
03676
03677 PIA(1)= -P3(3)*DET2(2,4)+P3(4)*DET2(2,3)+P3(2)*DET2(3,4)
03678 PIA(2)= -P3(4)*DET2(1,3)+P3(3)*DET2(1,4)-P3(1)*DET2(3,4)
03679 PIA(3)= P3(4)*DET2(1,2)-P3(2)*DET2(1,4)+P3(1)*DET2(2,4)
03680 PIA(4)= P3(3)*DET2(1,2)-P3(2)*DET2(1,3)+P3(1)*DET2(2,3)
03681
03682 DO 20 I=1,4
03683 20 PIA(I)=PIA(I)*SIGN
03684 END
03685
03686 SUBROUTINE DEXNEW(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
03687
03688
03689
03690
03691
03692
03693
03694
03695
03696
03697
03698
03699 COMMON / INOUT / INUT,IOUT
03700 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
03701 DATA IWARM/0/
03702
03703 IF(MODE.EQ.-1) THEN
03704
03705 IWARM=1
03706 CALL DADNEW( -1,ISGN,HV,PNU,PAA,PNPI,JDUMM)
03707
03708
03709
03710
03711 ELSEIF(MODE.EQ. 0) THEN
03712
03713 300 CONTINUE
03714 IF(IWARM.EQ.0) GOTO 902
03715 CALL DADNEW( 0,ISGN,HV,PNU,PAA,PNPI,JNPI)
03716 WT=(1+POL(1)*HV(1)+POL(2)*HV(2)+POL(3)*HV(3))/2.
03717
03718 CALL RANMAR(RN,1)
03719 IF(RN(1).GT.WT) GOTO 300
03720
03721 ELSEIF(MODE.EQ. 1) THEN
03722
03723 CALL DADNEW( 1,ISGN,HV,PNU,PAA,PNPI,JDUMM)
03724
03725 ENDIF
03726
03727 RETURN
03728 902 WRITE(IOUT, 9020)
03729 9020 FORMAT(' ----- DEXNEW: LACK OF INITIALISATION')
03730 STOP
03731 END
03732 SUBROUTINE DADNEW(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
03733
03734 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03735 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03736 * ,AMK,AMKZ,AMKST,GAMKST
03737
03738 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03739 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03740 * ,AMK,AMKZ,AMKST,GAMKST
03741 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03742 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03743 COMMON / TAUBMC / GAMPMC(500),GAMPER(500),NEVDEC(500)
03744 REAL*4 GAMPMC ,GAMPER
03745
03746 COMMON / INOUT / INUT,IOUT
03747
03748 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
03749
03750 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
03751
03752 & ,NAMES
03753 CHARACTER NAMES(NMODE)*31
03754
03755
03756 REAL*4 PNU(4),PWB(4),PNPI(4,9),HV(4),HHV(4)
03757 REAL*4 PDUM1(4),PDUM2(4),PDUMI(4,9)
03758 REAL*4 RRR(3)
03759 REAL*4 WTMAX(NMODE)
03760 REAL*8 SWT(NMODE),SSWT(NMODE)
03761 DIMENSION NEVRAW(NMODE),NEVOVR(NMODE),NEVACC(NMODE)
03762
03763 DATA PI /3.141592653589793238462643/
03764 DATA IWARM/0/
03765
03766 IF(MODE.EQ.-1) THEN
03767
03768
03769 NMOD=NMODE
03770 IWARM=1
03771
03772 DO 1 JNPI=1,NMOD
03773 NEVRAW(JNPI)=0
03774 NEVACC(JNPI)=0
03775 NEVOVR(JNPI)=0
03776 SWT(JNPI)=0
03777 SSWT(JNPI)=0
03778 WTMAX(JNPI)=-1.
03779
03780
03781
03782 NTRIALS = 45000
03783 IF (JNPI.LE.2) THEN
03784 A = PKORB(3,37+JNPI)
03785 IF (A.LT.0.) THEN
03786 NTRIALS = 30000
03787 ELSE
03788 NTRIALS = 0
03789 WTMAX(JNPI)=A
03790 END IF
03791 END IF
03792
03793 DO I=1,NTRIALS
03794
03795 IF (JNPI.LE.0) THEN
03796 GOTO 903
03797 ELSEIF(JNPI.LE.NM4) THEN
03798 CALL DPH4PI(WT,HV,PDUM1,PDUM2,PDUMI,JNPI)
03799
03800 ELSEIF(JNPI.LE.NM4+NM5) THEN
03801 CALL DPH5PI(WT,HV,PDUM1,PDUM2,PDUMI,JNPI)
03802 ELSEIF(JNPI.LE.NM4+NM5+NM6) THEN
03803 CALL DPHNPI(WT,HV,PDUM1,PDUM2,PDUMI,JNPI)
03804 ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3) THEN
03805 INUM=JNPI-NM4-NM5-NM6
03806 CALL DPHSPK(WT,HV,PDUM1,PDUM2,PDUMI,INUM)
03807 ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3+NM2) THEN
03808 INUM=JNPI-NM4-NM5-NM6-NM3
03809 CALL DPHSRK(WT,HV,PDUM1,PDUM2,PDUMI,INUM)
03810 ELSE
03811 GOTO 903
03812 ENDIF
03813 IF(WT.GT.WTMAX(JNPI)/1.2) WTMAX(JNPI)=WT*1.2
03814 ENDDO
03815
03816
03817
03818
03819
03820 1 CONTINUE
03821 WRITE(IOUT,7005)
03822
03823 ELSEIF(MODE.EQ. 0) THEN
03824
03825 IF(IWARM.EQ.0) GOTO 902
03826
03827 300 CONTINUE
03828 IF (JNPI.LE.0) THEN
03829 GOTO 903
03830 ELSEIF(JNPI.LE.NM4) THEN
03831 CALL DPH4PI(WT,HHV,PNU,PWB,PNPI,JNPI)
03832 ELSEIF(JNPI.LE.NM4+NM5) THEN
03833 CALL DPH5PI(WT,HHV,PNU,PWB,PNPI,JNPI)
03834 ELSEIF(JNPI.LE.NM4+NM5+NM6) THEN
03835 CALL DPHNPI(WT,HHV,PNU,PWB,PNPI,JNPI)
03836 ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3) THEN
03837 INUM=JNPI-NM4-NM5-NM6
03838 CALL DPHSPK(WT,HHV,PNU,PWB,PNPI,INUM)
03839 ELSEIF(JNPI.LE.NM4+NM5+NM6+NM3+NM2) THEN
03840 INUM=JNPI-NM4-NM5-NM6-NM3
03841 CALL DPHSRK(WT,HHV,PNU,PWB,PNPI,INUM)
03842 ELSE
03843 GOTO 903
03844 ENDIF
03845 DO I=1,4
03846 HV(I)=-ISGN*HHV(I)
03847 ENDDO
03848
03849 NEVRAW(JNPI)=NEVRAW(JNPI)+1
03850 SWT(JNPI)=SWT(JNPI)+WT
03851
03852
03853
03854 SSWT(JNPI)=SSWT(JNPI)+dble(WT)**2
03855
03856
03857 CALL RANMAR(RRR,3)
03858 RN=RRR(1)
03859 IF(WT.GT.WTMAX(JNPI)) NEVOVR(JNPI)=NEVOVR(JNPI)+1
03860 IF(RN*WTMAX(JNPI).GT.WT) GOTO 300
03861
03862 COSTHE=-1.+2.*RRR(2)
03863 THET=ACOS(COSTHE)
03864 PHI =2*PI*RRR(3)
03865 CALL ROTOR2(THET,PNU,PNU)
03866 CALL ROTOR3( PHI,PNU,PNU)
03867 CALL ROTOR2(THET,PWB,PWB)
03868 CALL ROTOR3( PHI,PWB,PWB)
03869 CALL ROTOR2(THET,HV,HV)
03870 CALL ROTOR3( PHI,HV,HV)
03871 ND=MULPIK(JNPI)
03872 DO 301 I=1,ND
03873 CALL ROTOR2(THET,PNPI(1,I),PNPI(1,I))
03874 CALL ROTOR3( PHI,PNPI(1,I),PNPI(1,I))
03875 301 CONTINUE
03876 NEVACC(JNPI)=NEVACC(JNPI)+1
03877
03878 ELSEIF(MODE.EQ. 1) THEN
03879
03880 DO 500 JNPI=1,NMOD
03881 IF(NEVRAW(JNPI).EQ.0) GOTO 500
03882 PARGAM=SWT(JNPI)/FLOAT(NEVRAW(JNPI)+1)
03883 ERROR=0
03884 IF(NEVRAW(JNPI).NE.0)
03885 & ERROR=SQRT(SSWT(JNPI)/SWT(JNPI)**2-1./FLOAT(NEVRAW(JNPI)))
03886 RAT=PARGAM/GAMEL
03887 WRITE(IOUT, 7010) NAMES(JNPI),
03888 & NEVRAW(JNPI),NEVACC(JNPI),NEVOVR(JNPI),PARGAM,RAT,ERROR
03889
03890 GAMPMC(8+JNPI-1)=RAT
03891 GAMPER(8+JNPI-1)=ERROR
03892
03893 500 CONTINUE
03894 ENDIF
03895
03896 RETURN
03897 7003 FORMAT(///1X,15(5H*****)
03898 $ /,' *', 25X,'******** DADNEW INITIALISATION ********',9X,1H*
03899 $ )
03900 7004 FORMAT(' *',E20.5,5X,'WTMAX = MAXIMUM WEIGHT ',9X,1H*/)
03901 7005 FORMAT(
03902 $ /,1X,15(5H*****)/)
03903 7010 FORMAT(///1X,15(5H*****)
03904 $ /,' *', 25X,'******** DADNEW FINAL REPORT ******** ',9X,1H*
03905 $ /,' *', 25X,'CHANNEL:',A31 ,9X,1H*
03906 $ /,' *',I20 ,5X,'NEVRAW = NO. OF DECAYS TOTAL ',9X,1H*
03907 $ /,' *',I20 ,5X,'NEVACC = NO. OF DECAYS ACCEPTED ',9X,1H*
03908 $ /,' *',I20 ,5X,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9X,1H*
03909 $ /,' *',E20.5,5X,'PARTIAL WTDTH IN GEV UNITS ',9X,1H*
03910 $ /,' *',F20.9,5X,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9X,1H*
03911 $ /,' *',F20.8,5X,'RELATIVE ERROR OF PARTIAL WIDTH ',9X,1H*
03912 $ /,1X,15(5H*****)/)
03913 902 WRITE(IOUT, 9020)
03914 9020 FORMAT(' ----- DADNEW: LACK OF INITIALISATION')
03915 STOP
03916 903 WRITE(IOUT, 9030) JNPI,MODE
03917 9030 FORMAT(' ----- DADNEW: WRONG JNPI',2I5)
03918 STOP
03919 END
03920
03921
03922 SUBROUTINE DPH4PI(DGAMT,HV,PN,PAA,PMULT,JNPI)
03923
03924
03925
03926
03927
03928
03929 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03930 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03931 * ,AMK,AMKZ,AMKST,GAMKST
03932
03933 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
03934 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
03935 * ,AMK,AMKZ,AMKST,GAMKST
03936 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03937 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
03938
03939
03940 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
03941 REAL PR(4),PIZ(4)
03942 REAL*4 RRR(9)
03943 REAL*8 UU,FF,FF1,FF2,FF3,FF4,GG1,GG2,GG3,GG4,RR
03944 DATA PI /3.141592653589793238462643/
03945 DATA ICONT /0/
03946 XLAM(X,Y,Z)=SQRT(ABS((X-Y-Z)**2-4.0*Y*Z))
03947
03948
03949
03950
03951 PHSPAC=1./2**23/PI**11
03952 PHSP=1./2**5/PI**2
03953
03954 IF (JNPI.EQ.1) THEN
03955 PREZ=0.7
03956
03957 AMP1=AMPI
03958 AMP2=AMPI
03959 AMP3=AMPI
03960 AMP4=AMPIZ
03961 AMRX=PKORB(1,14)
03962 GAMRX=PKORB(2,14)
03963
03964
03965
03966
03967
03968 AMROP =1.2
03969 GAMROP=.46
03970 ELSE
03971 PREZ=0.0
03972
03973 AMP1=AMPIZ
03974 AMP2=AMPIZ
03975 AMP3=AMPIZ
03976 AMP4=AMPI
03977
03978 AMRX=1.4
03979 GAMRX=.6
03980 AMROP =AMRX
03981 GAMROP=GAMRX
03982
03983 ENDIF
03984
03985 RRB=0.3
03986 CALL CHOICE(100+JNPI,RRB,ICHAN,PROB1,PROB2,PROB3,
03987
03988 $ AMROP,GAMROP,AMRX,GAMRX,AMRB,GAMRB)
03989 PREZ=PROB1+PROB2
03990
03991 PT(1)=0.
03992 PT(2)=0.
03993 PT(3)=0.
03994 PT(4)=AMTAU
03995
03996 CALL RANMAR(RRR,9)
03997
03998
03999
04000
04001 RR1=RRR(6)
04002 AMS1=(AMP1+AMP2+AMP3+AMP4)**2
04003 AMS2=(AMTAU-AMNUTA)**2
04004 ALP1=ATAN((AMS1-AMROP**2)/AMROP/GAMROP)
04005 ALP2=ATAN((AMS2-AMROP**2)/AMROP/GAMROP)
04006 ALP=ALP1+RR1*(ALP2-ALP1)
04007 AM4SQ =AMROP**2+AMROP*GAMROP*TAN(ALP)
04008 AM4 =SQRT(AM4SQ)
04009 PHSPAC=PHSPAC*
04010 $ ((AM4SQ-AMROP**2)**2+(AMROP*GAMROP)**2)/(AMROP*GAMROP)
04011 PHSPAC=PHSPAC*(ALP2-ALP1)
04012
04013
04014 RR1=RRR(1)
04015 AMS1=(AMP2+AMP3+AMP4)**2
04016 AMS2=(AM4-AMP1)**2
04017 IF (RRR(9).GT.PREZ) THEN
04018 AM3SQ=AMS1+ RR1*(AMS2-AMS1)
04019 AM3 =SQRT(AM3SQ)
04020
04021 FF1=AMS2-AMS1
04022 ELSE
04023
04024 ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX)
04025 ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX)
04026 ALP=ALP1+RR1*(ALP2-ALP1)
04027 AM3SQ =AMRX**2+AMRX*GAMRX*TAN(ALP)
04028 AM3 =SQRT(AM3SQ)
04029
04030 FF1=((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX)
04031 FF1=FF1*(ALP2-ALP1)
04032 ENDIF
04033
04034 RR2=RRR(2)
04035 AMS1=(AMP3+AMP4)**2
04036 AMS2=(AM3-AMP2)**2
04037
04038 AM2SQ=AMS1+ RR2*(AMS2-AMS1)
04039 AM2 =SQRT(AM2SQ)
04040
04041 FF2=(AMS2-AMS1)
04042
04043
04044 ENQ1=(AM2SQ-AMP3**2+AMP4**2)/(2*AM2)
04045 ENQ2=(AM2SQ+AMP3**2-AMP4**2)/(2*AM2)
04046 PPI= ENQ1**2-AMP4**2
04047 PPPI=SQRT(ABS(ENQ1**2-AMP4**2))
04048
04049 PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AM2)
04050
04051
04052
04053 CALL SPHERA(PPPI,PIZ)
04054 PIZ(4)=ENQ1
04055
04056
04057
04058 DO 30 I=1,3
04059 30 PIPL(I)=-PIZ(I)
04060 PIPL(4)=ENQ2
04061
04062
04063
04064
04065 PR(1)=0
04066 PR(2)=0
04067 PR(4)=1./(2*AM3)*(AM3**2+AM2**2-AMP2**2)
04068 PR(3)= SQRT(ABS(PR(4)**2-AM2**2))
04069 PPI = PR(4)**2-AM2**2
04070
04071
04072
04073 PIM1(1)=0
04074 PIM1(2)=0
04075 PIM1(4)=1./(2*AM3)*(AM3**2-AM2**2+AMP2**2)
04076 PIM1(3)=-PR(3)
04077
04078 FF3=(4*PI)*(2*PR(3)/AM3)
04079
04080 EXE=(PR(4)+PR(3))/AM2
04081 CALL BOSTR3(EXE,PIZ,PIZ)
04082 CALL BOSTR3(EXE,PIPL,PIPL)
04083 RR3=RRR(3)
04084 RR4=RRR(4)
04085 THET =ACOS(-1.+2*RR3)
04086 PHI = 2*PI*RR4
04087 CALL ROTPOL(THET,PHI,PIPL)
04088 CALL ROTPOL(THET,PHI,PIM1)
04089 CALL ROTPOL(THET,PHI,PIZ)
04090 CALL ROTPOL(THET,PHI,PR)
04091
04092
04093
04094
04095 PR(1)=0
04096 PR(2)=0
04097 PR(4)=1./(2*AM4)*(AM4**2+AM3**2-AMP1**2)
04098 PR(3)= SQRT(ABS(PR(4)**2-AM3**2))
04099 PPI = PR(4)**2-AM3**2
04100
04101
04102
04103 PIM2(1)=0
04104 PIM2(2)=0
04105 PIM2(4)=1./(2*AM4)*(AM4**2-AM3**2+AMP1**2)
04106 PIM2(3)=-PR(3)
04107
04108 FF4=(4*PI)*(2*PR(3)/AM4)
04109
04110 EXE=(PR(4)+PR(3))/AM3
04111 CALL BOSTR3(EXE,PIZ,PIZ)
04112 CALL BOSTR3(EXE,PIPL,PIPL)
04113 CALL BOSTR3(EXE,PIM1,PIM1)
04114 RR3=RRR(7)
04115 RR4=RRR(8)
04116 THET =ACOS(-1.+2*RR3)
04117 PHI = 2*PI*RR4
04118 CALL ROTPOL(THET,PHI,PIPL)
04119 CALL ROTPOL(THET,PHI,PIM1)
04120 CALL ROTPOL(THET,PHI,PIM2)
04121 CALL ROTPOL(THET,PHI,PIZ)
04122 CALL ROTPOL(THET,PHI,PR)
04123
04124
04125
04126 PAA(1)=0
04127 PAA(2)=0
04128 PAA(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AM4**2)
04129 PAA(3)= SQRT(ABS(PAA(4)**2-AM4**2))
04130 PPI = PAA(4)**2-AM4**2
04131 PHSPAC=PHSPAC*(4*PI)*(2*PAA(3)/AMTAU)
04132 PHSP=PHSP*(4*PI)*(2*PAA(3)/AMTAU)
04133
04134 PN(1)=0
04135 PN(2)=0
04136 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AM4**2)
04137 PN(3)=-PAA(3)
04138
04139 IF(RRR(9).LE.0.5*PREZ) THEN
04140 DO 72 I=1,4
04141 X=PIM1(I)
04142 PIM1(I)=PIM2(I)
04143 72 PIM2(I)=X
04144 ENDIF
04145
04146
04147
04148 AM3SQ=(PIM1(4)+PIZ(4)+PIPL(4))**2-(PIM1(3)+PIZ(3)+PIPL(3))**2
04149 $ -(PIM1(2)+PIZ(2)+PIPL(2))**2-(PIM1(1)+PIZ(1)+PIPL(1))**2
04150 AMS2=(AM4-AMP2)**2
04151 AMS1=(AMP1+AMP3+AMP4)**2
04152 FF1=(AMS2-AMS1)
04153 AMS1=(AMP3+AMP4)**2
04154 AMS2=(SQRT(AM3SQ)-AMP1)**2
04155 FF2=AMS2-AMS1
04156 FF3=(4*PI)*(XLAM(AM2**2,AMP1**2,AM3SQ)/AM3SQ)
04157 FF4=(4*PI)*(XLAM(AM3SQ,AMP2**2,AM4**2)/AM4**2)
04158 UU=FF1*FF2*FF3*FF4
04159
04160 AM3SQ=(PIM1(4)+PIZ(4)+PIPL(4))**2-(PIM1(3)+PIZ(3)+PIPL(3))**2
04161 $ -(PIM1(2)+PIZ(2)+PIPL(2))**2-(PIM1(1)+PIZ(1)+PIPL(1))**2
04162 AMS2=(AM4-AMP2)**2
04163 AMS1=(AMP1+AMP3+AMP4)**2
04164 ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX)
04165 ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX)
04166 FF1=((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX)
04167 FF1=FF1*(ALP2-ALP1)
04168 AMS1=(AMP3+AMP4)**2
04169 AMS2=(SQRT(AM3SQ)-AMP1)**2
04170 FF2=AMS2-AMS1
04171 FF3=(4*PI)*(XLAM(AM2**2,AMP1**2,AM3SQ)/AM3SQ)
04172 FF4=(4*PI)*(XLAM(AM3SQ,AMP2**2,AM4**2)/AM4**2)
04173 FF=FF1*FF2*FF3*FF4
04174
04175 AM3SQ=(PIM2(4)+PIZ(4)+PIPL(4))**2-(PIM2(3)+PIZ(3)+PIPL(3))**2
04176 $ -(PIM2(2)+PIZ(2)+PIPL(2))**2-(PIM2(1)+PIZ(1)+PIPL(1))**2
04177 AMS2=(AM4-AMP1)**2
04178 AMS1=(AMP2+AMP3+AMP4)**2
04179 ALP1=ATAN((AMS1-AMRX**2)/AMRX/GAMRX)
04180 ALP2=ATAN((AMS2-AMRX**2)/AMRX/GAMRX)
04181 GG1=((AM3SQ-AMRX**2)**2+(AMRX*GAMRX)**2)/(AMRX*GAMRX)
04182 GG1=GG1*(ALP2-ALP1)
04183 AMS1=(AMP3+AMP4)**2
04184 AMS2=(SQRT(AM3SQ)-AMP2)**2
04185 GG2=AMS2-AMS1
04186 GG3=(4*PI)*(XLAM(AM2**2,AMP2**2,AM3SQ)/AM3SQ)
04187 GG4=(4*PI)*(XLAM(AM3SQ,AMP1**2,AM4**2)/AM4**2)
04188 GG=GG1*GG2*GG3*GG4
04189
04190 IF ( ( (FF+GG)*UU+FF*GG ).GT.0.0D0) THEN
04191 RR=FF*GG*UU/(0.5*PREZ*(FF+GG)*UU+(1.0-PREZ)*FF*GG)
04192 PHSPAC=PHSPAC*RR
04193 ELSE
04194 PHSPAC=0.0
04195 ENDIF
04196
04197 IF (JNPI.EQ.1) THEN
04198 RR5= RRR(5)
04199 IF(RR5.LE.0.5) THEN
04200 DO 70 I=1,4
04201 X=PIM1(I)
04202 PIM1(I)=PIM2(I)
04203 70 PIM2(I)=X
04204 ENDIF
04205 PHSPAC=PHSPAC/2.
04206 ELSE
04207
04208 RR5= RRR(5)
04209 IF(RR5.LE.0.5) THEN
04210 DO 71 I=1,4
04211 X=PIM1(I)
04212 PIM1(I)=PIM2(I)
04213 71 PIM2(I)=X
04214 ENDIF
04215 PHSPAC=PHSPAC/6.
04216 ENDIF
04217
04218
04219 EXE=(PAA(4)+PAA(3))/AM4
04220 CALL BOSTR3(EXE,PIZ,PIZ)
04221 CALL BOSTR3(EXE,PIPL,PIPL)
04222 CALL BOSTR3(EXE,PIM1,PIM1)
04223 CALL BOSTR3(EXE,PIM2,PIM2)
04224 CALL BOSTR3(EXE,PR,PR)
04225
04226
04227
04228
04229
04230
04231
04232
04233 IF (JNPI.EQ.1) THEN
04234 CALL DAM4PI(JNPI,PT,PN,PIM1,PIM2,PIZ,PIPL,AMPLIT,HV)
04235 ELSEIF (JNPI.EQ.2) THEN
04236 CALL DAM4PI(JNPI,PT,PN,PIM1,PIM2,PIPL,PIZ,AMPLIT,HV)
04237 ELSE
04238 CALL DAM4PI(JNPI,PT,PN,PIM1,PIM2,PIPL,PIZ,AMPLIT,HV)
04239 ENDIF
04240
04241 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
04242
04243
04244 DO 77 K=1,4
04245 PMULT(K,1)=PIM1(K)
04246 PMULT(K,2)=PIM2(K)
04247
04248 PMULT(K,3)=PIPL(K)
04249 PMULT(K,4)=PIZ (K)
04250
04251 77 CONTINUE
04252 END
04253 SUBROUTINE DAM4PI(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
04254
04255
04256
04257
04258
04259
04260
04261
04262
04263
04264
04265 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04266 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04267 * ,AMK,AMKZ,AMKST,GAMKST
04268
04269 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04270 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04271 * ,AMK,AMKZ,AMKST,GAMKST
04272 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04273 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04274 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
04275 REAL PIVEC(4),PIAKS(4),HVM(4)
04276 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
04277 EXTERNAL FORM1,FORM2,FORM3,FORM4,FORM5
04278 DATA PI /3.141592653589793238462643/
04279 DATA ICONT /0/
04280
04281
04282 CALL CURR_CLEO(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
04283
04284
04285
04286 CALL CLVEC(HADCUR,PN,PIVEC)
04287 CALL CLAXI(HADCUR,PN,PIAKS)
04288 CALL CLNUT(HADCUR,BRAKM,HVM)
04289
04290 BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4)
04291 & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM
04292 AMPLIT=(CCABIB*GFERMI)**2*BRAK/2.
04293
04294 DO 90 I=1,3
04295 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I)))
04296 & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I)
04297
04298 IF (BRAK.NE.0.0)
04299 &HV(I)=-HV(I)/BRAK
04300 90 CONTINUE
04301 END
04302 SUBROUTINE DAM5PI(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,PIM5,AMPLIT,HV)
04303
04304
04305
04306
04307
04308
04309
04310
04311
04312
04313
04314 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04315 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04316 * ,AMK,AMKZ,AMKST,GAMKST
04317
04318 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04319 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04320 * ,AMK,AMKZ,AMKST,GAMKST
04321 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04322 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04323 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4),PIM5(4)
04324 REAL PIVEC(4),PIAKS(4),HVM(4)
04325 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
04326 EXTERNAL FORM1,FORM2,FORM3,FORM4,FORM5
04327 DATA PI /3.141592653589793238462643/
04328 DATA ICONT /0/
04329
04330
04331 CALL CURR5(MNUM,PIM1,PIM2,PIM3,PIM4,PIM5,HADCUR)
04332
04333
04334
04335 CALL CLVEC(HADCUR,PN,PIVEC)
04336 CALL CLAXI(HADCUR,PN,PIAKS)
04337 CALL CLNUT(HADCUR,BRAKM,HVM)
04338
04339 BRAK= (GV**2+GA**2)*PT(4)*PIVEC(4) +2.*GV*GA*PT(4)*PIAKS(4)
04340 & +2.*(GV**2-GA**2)*AMNUTA*AMTAU*BRAKM
04341 AMPLIT=(CCABIB*GFERMI)**2*BRAK/2.
04342
04343 DO 90 I=1,3
04344 HV(I)=-(AMTAU*((GV**2+GA**2)*PIAKS(I)+2.*GV*GA*PIVEC(I)))
04345 & +(GV**2-GA**2)*AMNUTA*AMTAU*HVM(I)
04346
04347 IF (BRAK.NE.0.0)
04348 &HV(I)=-HV(I)/BRAK
04349 90 CONTINUE
04350 END
04351 SUBROUTINE DPH5PI(DGAMT,HV,PN,PAA,PMULT,JNPI)
04352
04353
04354
04355
04356 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04357 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04358 * ,AMK,AMKZ,AMKST,GAMKST
04359
04360 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04361 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04362
04363
04364 * ,AMK,AMKZ,AMKST,GAMKST
04365 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04366 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04367 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
04368
04369 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
04370
04371 & ,NAMES
04372 CHARACTER NAMES(NMODE)*31
04373 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
04374 REAL*4 PR(4),PI1(4),PI2(4),PI3(4),PI4(4),PI5(4)
04375 REAL*8 AMP1,AMP2,AMP3,AMP4,AMP5,ams1,ams2,amom,gamom
04376 REAL*8 AM5SQ,AM4SQ,AM3SQ,AM2SQ,AM5,AM4,AM3
04377 REAL*4 RRR(12)
04378 REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
04379
04380 REAL*8 XM,AM,GAMMA
04381
04382 real*8 phspac
04383
04384
04385 DATA PI /3.141592653589793238462643/
04386 DATA ICONT /0/
04387 data fpi /93.3e-3/
04388
04389 COMPLEX BWIGN
04390
04391
04392 BWIGN(XM,AM,GAMMA)=XM**2/CMPLX(XM**2-AM**2,GAMMA*AM)
04393
04394
04395 PROBa2=0.7
04396 PROBOM=0.7
04397 ama2=1.260
04398 gama2=0.400
04399
04400 AMOM=.782
04401 GAMOM=.7
04402
04403 IF (JNPI.EQ.23.OR.JNPI.EQ.24) GAMOM= 0.0085
04404
04405
04406
04407
04408 PHSPAC=1./2**29/PI**14
04409
04410
04411 AMP1=DCDMAS(IDFFIN(1,JNPI))
04412 AMP2=DCDMAS(IDFFIN(2,JNPI))
04413 AMP3=DCDMAS(IDFFIN(3,JNPI))
04414 AMP4=DCDMAS(IDFFIN(4,JNPI))
04415 AMP5=DCDMAS(IDFFIN(5,JNPI))
04416
04417
04418 PT(1)=0.
04419 PT(2)=0.
04420 PT(3)=0.
04421 PT(4)=AMTAU
04422
04423 CALL RANMAR(RRR,12)
04424
04425
04426
04427
04428
04429 IF (RRR(11).GT.PROBa2) THEN
04430
04431 rr1=rrr(10)
04432 ams1=(amp1+amp2+amp3+amp4+amp5)**2
04433 ams2=(amtau-amnuta)**2
04434 alp1=atan((ams1-ama2**2)/ama2/gama2)
04435 alp2=atan((ams2-ama2**2)/ama2/gama2)
04436 am5sq=ams1+ rr1*(ams2-ams1)
04437 am5 =sqrt(am5sq)
04438
04439
04440 ELSE
04441 rr1=rrr(10)
04442 ams1=(amp1+amp2+amp3+amp4+amp5)**2
04443 ams2=(amtau-amnuta)**2
04444 alp1=atan((ams1-ama2**2)/ama2/gama2)
04445 alp2=atan((ams2-ama2**2)/ama2/gama2)
04446 alp=alp1+rr1*(alp2-alp1)
04447 am5sq =ama2**2+ama2*gama2*tan(alp)
04448 am5 =sqrt(am5sq)
04449 ENDIF
04450
04451 gg5=((am5sq-ama2**2)**2+(ama2*gama2)**2)/(ama2*gama2)
04452 gg5=gg5*(alp2-alp1)
04453 phspac=phspac/(PROBa2/gg5+(1D0-PROBa2)/(ams2-ams1) )
04454
04455
04456
04457 rr1=rrr(9)
04458 ams1=(amp2+amp3+amp4+amp5)**2
04459 ams2=(am5-amp1)**2
04460 am4sq=ams1+ rr1*(ams2-ams1)
04461 am4 =sqrt(am4sq)
04462 gg1=ams2-ams1
04463
04464
04465
04466 IF (RRR(12).LT.PROBom) THEN
04467
04468 rr1=rrr(1)
04469 ams1=(amp2+amp3+amp4)**2
04470 ams2=(am4-amp5)**2
04471 alp1=atan((ams1-amom**2)/amom/gamom)
04472 alp2=atan((ams2-amom**2)/amom/gamom)
04473 alp=alp1+rr1*(alp2-alp1)
04474 am3sq =amom**2+amom*gamom*tan(alp)
04475 am3 =sqrt(am3sq)
04476 ELSE
04477
04478 rr1=rrr(1)
04479 ams1=(amp2+amp3+amp4)**2
04480 ams2=(am4-amp5)**2
04481 alp1=atan((ams1-amom**2)/amom/gamom)
04482 alp2=atan((ams2-amom**2)/amom/gamom)
04483
04484 am3sq=ams1+ rr1*(ams2-ams1)
04485 am3 =sqrt(am3sq)
04486
04487 ENDIF
04488
04489 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
04490 gg2=gg2*(alp2-alp1)
04491 gg2=1D0/(PROBOM/gg2+(1D0-PROBOM)/(ams2-ams1))
04492
04493
04494 rr2=rrr(2)
04495 ams1=(amp3+amp4)**2
04496 ams2=(am3-amp2)**2
04497
04498 am2sq=ams1+ rr2*(ams2-ams1)
04499 am2 =sqrt(am2sq)
04500
04501 gg3=ams2-ams1
04502
04503
04504 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
04505 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
04506 ppi= enq1**2-amp3**2
04507 pppi=sqrt(abs(enq1**2-amp3**2))
04508 ff1=(4*pi)*(2*pppi/am2)
04509
04510 call sphera(pppi,pi3)
04511 pi3(4)=enq1
04512
04513 do 30 i=1,3
04514 30 pi4(i)=-pi3(i)
04515 pi4(4)=enq2
04516
04517
04518
04519 pr(1)=0
04520 pr(2)=0
04521 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
04522 pr(3)= sqrt(abs(pr(4)**2-am2**2))
04523 ppi = pr(4)**2-am2**2
04524
04525 pi2(1)=0
04526 pi2(2)=0
04527 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
04528 pi2(3)=-pr(3)
04529
04530 ff2=(4*pi)*(2*pr(3)/am3)
04531
04532 exe=(pr(4)+pr(3))/am2
04533 call bostr3(exe,pi3,pi3)
04534 call bostr3(exe,pi4,pi4)
04535 rr3=rrr(3)
04536 rr4=rrr(4)
04537 thet =acos(-1.+2*rr3)
04538 phi = 2*pi*rr4
04539 call rotpol(thet,phi,pi2)
04540 call rotpol(thet,phi,pi3)
04541 call rotpol(thet,phi,pi4)
04542
04543
04544
04545 pr(1)=0
04546 pr(2)=0
04547 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
04548 pr(3)= sqrt(abs(pr(4)**2-am3**2))
04549 ppi = pr(4)**2-am3**2
04550
04551 pi5(1)=0
04552 pi5(2)=0
04553 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
04554 pi5(3)=-pr(3)
04555
04556 ff3=(4*pi)*(2*pr(3)/am4)
04557
04558 exe=(pr(4)+pr(3))/am3
04559 call bostr3(exe,pi2,pi2)
04560 call bostr3(exe,pi3,pi3)
04561 call bostr3(exe,pi4,pi4)
04562 rr3=rrr(5)
04563 rr4=rrr(6)
04564 thet =acos(-1.+2*rr3)
04565 phi = 2*pi*rr4
04566 call rotpol(thet,phi,pi2)
04567 call rotpol(thet,phi,pi3)
04568 call rotpol(thet,phi,pi4)
04569 call rotpol(thet,phi,pi5)
04570
04571
04572
04573 pr(1)=0
04574 pr(2)=0
04575 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
04576 pr(3)= sqrt(abs(pr(4)**2-am4**2))
04577 ppi = pr(4)**2-am4**2
04578
04579 pi1(1)=0
04580 pi1(2)=0
04581 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
04582 pi1(3)=-pr(3)
04583
04584 ff4=(4*pi)*(2*pr(3)/am5)
04585
04586 exe=(pr(4)+pr(3))/am4
04587 call bostr3(exe,pi2,pi2)
04588 call bostr3(exe,pi3,pi3)
04589 call bostr3(exe,pi4,pi4)
04590 call bostr3(exe,pi5,pi5)
04591 rr3=rrr(7)
04592 rr4=rrr(8)
04593 thet =acos(-1.+2*rr3)
04594 phi = 2*pi*rr4
04595 call rotpol(thet,phi,pi1)
04596 call rotpol(thet,phi,pi2)
04597 call rotpol(thet,phi,pi3)
04598 call rotpol(thet,phi,pi4)
04599 call rotpol(thet,phi,pi5)
04600
04601
04602
04603 paa(1)=0
04604 paa(2)=0
04605
04606
04607
04608 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
04609 paa(3)= sqrt(abs(paa(4)**2-am5sq))
04610 ppi = paa(4)**2-am5sq
04611 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
04612
04613 pn(1)=0
04614 pn(2)=0
04615 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
04616 pn(3)=-paa(3)
04617
04618 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
04619
04620
04621
04622 exe=(paa(4)+paa(3))/am5
04623 call bostr3(exe,pi1,pi1)
04624 call bostr3(exe,pi2,pi2)
04625 call bostr3(exe,pi3,pi3)
04626 call bostr3(exe,pi4,pi4)
04627 call bostr3(exe,pi5,pi5)
04628
04629
04630
04631
04632
04633 PXQ=AMTAU*PAA(4)
04634 PXN=AMTAU*PN(4)
04635 QXN=PAA(4)*PN(4)-PAA(1)*PN(1)-PAA(2)*PN(2)-PAA(3)*PN(3)
04636 BRAK=2*(GV**2+GA**2)*(2*PXQ*QXN+AM5SQ*PXN)
04637 & -6*(GV**2-GA**2)*AMTAU*AMNUTA*AM5SQ
04638 fompp = cabs(bwign(am3,amom,gamom))**2
04639
04640
04641 fnorm = 1/fpi**6
04642
04643 AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK
04644 amplit = amplit * fompp * fnorm
04645
04646
04647
04648
04649
04650 DO 40 I=1,3
04651 40 HV(I)=0.
04652
04653
04654
04655
04656 if (jnpi.gt.23)
04657 $ CALL DAM5PI(jnpi,PT,PN,PI1,PI2,PI3,PI4,PI5,AMPLIT,HV)
04658 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
04659
04660 do 77 k=1,4
04661 pmult(k,1)=pi1(k)
04662 pmult(k,2)=pi2(k)
04663 pmult(k,3)=pi3(k)
04664 pmult(k,4)=pi4(k)
04665 pmult(k,5)=pi5(k)
04666 77 continue
04667 return
04668
04669
04670
04671
04672
04673
04674 end
04675 SUBROUTINE DPHSRK(DGAMT,HV,PN,PR,PMULT,INUM)
04676
04677
04678
04679
04680
04681 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04682 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04683 * ,AMK,AMKZ,AMKST,GAMKST
04684
04685 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04686 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04687 * ,AMK,AMKZ,AMKST,GAMKST
04688 COMMON / DECPAR / GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04689 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
04690 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
04691
04692 REAL RR1(1)
04693
04694 DATA PI /3.141592653589793238462643/
04695 DATA ICONT /0/
04696
04697
04698 PHSPAC=1./2**11/PI**5
04699
04700 PT(1)=0.
04701 PT(2)=0.
04702 PT(3)=0.
04703 PT(4)=AMTAU
04704
04705 AMS1=(AMK+AMKZ)**2
04706 AMS2=(AMTAU-AMNUTA)**2
04707
04708 CALL RANMAR(RR1,1)
04709 AMX2=AMS1+ RR1(1)*(AMS2-AMS1)
04710 AMX=SQRT(AMX2)
04711 PHSPAC=PHSPAC*(AMS2-AMS1)
04712
04713
04714
04715
04716 100 CONTINUE
04717
04718
04719
04720
04721
04722
04723
04724
04725
04726
04727 PN(1)=0
04728 PN(2)=0
04729 PN(4)=1./(2*AMTAU)*(AMTAU**2+AMNUTA**2-AMX**2)
04730 PN(3)=-SQRT((PN(4)-AMNUTA)*(PN(4)+AMNUTA))
04731
04732 PR(1)=0
04733 PR(2)=0
04734 PR(4)=1./(2*AMTAU)*(AMTAU**2-AMNUTA**2+AMX**2)
04735 PR(3)=-PN(3)
04736 PHSPAC=PHSPAC*(4*PI)*(2*PR(3)/AMTAU)
04737
04738
04739 ENQ1=(AMX2+AMK**2-AMKZ**2)/(2.*AMX)
04740 ENQ2=(AMX2-AMK**2+AMKZ**2)/(2.*AMX)
04741 PPPI=SQRT((ENQ1-AMK)*(ENQ1+AMK))
04742 PHSPAC=PHSPAC*(4*PI)*(2*PPPI/AMX)
04743
04744 CALL SPHERA(PPPI,PKC)
04745 PKC(4)=ENQ1
04746
04747 DO 20 I=1,3
04748 20 PKZ(I)=-PKC(I)
04749 PKZ(4)=ENQ2
04750 EXE=(PR(4)+PR(3))/AMX
04751
04752 CALL BOSTR3(EXE,PKC,PKC)
04753 CALL BOSTR3(EXE,PKZ,PKZ)
04754 DO 30 I=1,4
04755 30 QQ(I)=PKC(I)-PKZ(I)
04756
04757 PKSD =PR(4)*PR(4)-PR(3)*PR(3)-PR(2)*PR(2)-PR(1)*PR(1)
04758 QQPKS=PR(4)* QQ(4)-PR(3)* QQ(3)-PR(2)* QQ(2)-PR(1)* QQ(1)
04759 DO 31 I=1,4
04760 31 QQ(I)=QQ(I)-PR(I)*QQPKS/PKSD
04761
04762 PRODPQ=PT(4)*QQ(4)
04763 PRODNQ=PN(4)*QQ(4)-PN(1)*QQ(1)-PN(2)*QQ(2)-PN(3)*QQ(3)
04764 PRODPN=PT(4)*PN(4)
04765 QQ2= QQ(4)**2-QQ(1)**2-QQ(2)**2-QQ(3)**2
04766 BRAK=(GV**2+GA**2)*(2*PRODPQ*PRODNQ-PRODPN*QQ2)
04767 & +(GV**2-GA**2)*AMTAU*AMNUTA*QQ2
04768 AMPLIT=(GFERMI*CCABIB)**2*BRAK*2*FPIRK(AMX)
04769 DGAMT=1/(2.*AMTAU)*AMPLIT*PHSPAC
04770 DO 40 I=1,3
04771 40 HV(I)=2*GV*GA*AMTAU*(2*PRODNQ*QQ(I)-QQ2*PN(I))/BRAK
04772 do 77 k=1,4
04773 pmult(k,1)=pkc(k)
04774 pmult(k,2)=pkz(k)
04775 77 continue
04776 RETURN
04777 END
04778 FUNCTION FPIRK(W)
04779
04780
04781
04782 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04783 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04784 * ,AMK,AMKZ,AMKST,GAMKST
04785
04786 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
04787 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
04788 * ,AMK,AMKZ,AMKST,GAMKST
04789
04790 COMPLEX FPIKM
04791 FPIRK=CABS(FPIKM(W,AMK,AMKZ))**2
04792
04793 END
04794 COMPLEX FUNCTION FPIKMK(W,XM1,XM2)
04795
04796
04797
04798 COMPLEX BWIGM
04799 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
04800 EXTERNAL BWIG
04801 DATA INIT /0/
04802
04803
04804 IF (INIT.EQ.0 ) THEN
04805 INIT=1
04806 PI=3.141592654
04807 PIM=.140
04808 ROM=0.773
04809 ROG=0.145
04810 ROM1=1.570
04811 ROG1=0.510
04812
04813 BETA1=-0.221
04814 ENDIF
04815
04816 S=W**2
04817 FPIKMK=(BWIGM(S,ROM,ROG,XM1,XM2)+BETA1*BWIGM(S,ROM1,ROG1,XM1,XM2))
04818 & /(1+BETA1)
04819 RETURN
04820 END
04821 SUBROUTINE RESLUX
04822
04823
04824
04825
04826 NHEP=0
04827 END
04828 SUBROUTINE DWRPH(KTO,PHX)
04829
04830
04831
04832 IMPLICIT REAL*8 (A-H,O-Z)
04833 REAL*4 PHX(4)
04834 REAL*4 QHOT(4)
04835
04836 DO 9 K=1,4
04837 QHOT(K) =0.0
04838 9 CONTINUE
04839
04840
04841 DO 1002 I=1,4
04842 1002 QHOT(I)=PHX(I)
04843 IF (QHOT(4).GT.1.E-5) CALL DWLUPH(KTO,QHOT)
04844 RETURN
04845 END
04846 SUBROUTINE DWLUPH(KTO,PHOT)
04847
04848
04849
04850
04851
04852
04853
04854
04855
04856
04857
04858 REAL PHOT(4)
04859
04860 COMMON /TAUPOS/ NP1,NP2
04861
04862
04863
04864 IF (PHOT(4).LE.0.0) RETURN
04865
04866
04867 IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
04868 NPS=NP1
04869 ELSE
04870 NPS=NP2
04871 ENDIF
04872
04873 KTOS=KTO
04874 IF(KTOS.GT.10) KTOS=KTOS-10
04875
04876 CALL TRALO4(KTOS,PHOT,PHOT,AM)
04877 CALL FILHEP(0,1,22,NPS,NPS,0,0,PHOT,0.0,.TRUE.)
04878
04879 RETURN
04880 END
04881
04882 SUBROUTINE DWLUEL(KTO,ISGN,PNU,PWB,PEL,PNE)
04883
04884
04885
04886
04887
04888
04889
04890
04891
04892
04893
04894 REAL PNU(4),PWB(4),PEL(4),PNE(4)
04895
04896 COMMON /TAUPOS/ NP1,NP2
04897
04898
04899
04900 IF(KTO.EQ. 1) THEN
04901 NPS=NP1
04902 ELSE
04903 NPS=NP2
04904 ENDIF
04905
04906
04907 CALL TRALO4(KTO,PNU,PNU,AM)
04908 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
04909
04910
04911 CALL TRALO4(KTO,PWB,PWB,AM)
04912
04913
04914
04915 CALL TRALO4(KTO,PEL,PEL,AM)
04916 CALL FILHEP(0,1,11*ISGN,NPS,NPS,0,0,PEL,AM,.FALSE.)
04917
04918
04919 CALL TRALO4(KTO,PNE,PNE,AM)
04920 CALL FILHEP(0,1,-12*ISGN,NPS,NPS,0,0,PNE,AM,.TRUE.)
04921
04922 RETURN
04923 END
04924 SUBROUTINE DWLUMU(KTO,ISGN,PNU,PWB,PMU,PNM)
04925
04926
04927
04928
04929
04930
04931
04932
04933
04934
04935
04936 REAL PNU(4),PWB(4),PMU(4),PNM(4)
04937
04938 COMMON /TAUPOS/ NP1,NP2
04939
04940
04941
04942 IF(KTO.EQ. 1) THEN
04943 NPS=NP1
04944 ELSE
04945 NPS=NP2
04946 ENDIF
04947
04948
04949 CALL TRALO4(KTO,PNU,PNU,AM)
04950 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
04951
04952
04953 CALL TRALO4(KTO,PWB,PWB,AM)
04954
04955
04956
04957 CALL TRALO4(KTO,PMU,PMU,AM)
04958 CALL FILHEP(0,1,13*ISGN,NPS,NPS,0,0,PMU,AM,.FALSE.)
04959
04960
04961 CALL TRALO4(KTO,PNM,PNM,AM)
04962 CALL FILHEP(0,1,-14*ISGN,NPS,NPS,0,0,PNM,AM,.TRUE.)
04963
04964 RETURN
04965 END
04966 SUBROUTINE DWLUPI(KTO,ISGN,PPI,PNU)
04967
04968
04969
04970
04971
04972
04973
04974
04975
04976 REAL PNU(4),PPI(4)
04977 COMMON /TAUPOS/ NP1,NP2
04978
04979
04980 IF(KTO.EQ. 1) THEN
04981 NPS=NP1
04982 ELSE
04983 NPS=NP2
04984 ENDIF
04985
04986
04987 CALL TRALO4(KTO,PNU,PNU,AM)
04988 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
04989
04990
04991 CALL TRALO4(KTO,PPI,PPI,AM)
04992 CALL FILHEP(0,1,-211*ISGN,NPS,NPS,0,0,PPI,AM,.TRUE.)
04993
04994 RETURN
04995 END
04996 SUBROUTINE DWLURO(KTO,ISGN,PNU,PRHO,PIC,PIZ)
04997
04998
04999
05000
05001
05002
05003
05004
05005
05006
05007
05008 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
05009
05010 COMMON /TAUPOS/ NP1,NP2
05011
05012
05013
05014 IF(KTO.EQ. 1) THEN
05015 NPS=NP1
05016 ELSE
05017 NPS=NP2
05018 ENDIF
05019
05020
05021 CALL TRALO4(KTO,PNU,PNU,AM)
05022 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
05023
05024
05025 CALL TRALO4(KTO,PRHO,PRHO,AM)
05026 CALL FILHEP(0,2,-213*ISGN,NPS,NPS,0,0,PRHO,AM,.TRUE.)
05027
05028
05029 CALL TRALO4(KTO,PIC,PIC,AM)
05030 CALL FILHEP(0,1,-211*ISGN,-1,-1,0,0,PIC,AM,.TRUE.)
05031
05032
05033 CALL TRALO4(KTO,PIZ,PIZ,AM)
05034 CALL FILHEP(0,1,111,-2,-2,0,0,PIZ,AM,.TRUE.)
05035
05036 RETURN
05037 END
05038 SUBROUTINE DWLUAA(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
05039
05040
05041
05042
05043
05044
05045
05046
05047
05048
05049
05050
05051 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
05052
05053 COMMON /TAUPOS/ NP1,NP2
05054
05055
05056
05057 IF(KTO.EQ. 1) THEN
05058 NPS=NP1
05059 ELSE
05060 NPS=NP2
05061 ENDIF
05062
05063
05064 CALL TRALO4(KTO,PNU,PNU,AM)
05065 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
05066
05067
05068 CALL TRALO4(KTO,PAA,PAA,AM)
05069 CALL FILHEP(0,1,-20213*ISGN,NPS,NPS,0,0,PAA,AM,.TRUE.)
05070
05071
05072 IF(JAA.EQ.1) THEN
05073
05074
05075
05076
05077 CALL TRALO4(KTO,PIM2,PIM2,AM)
05078 CALL FILHEP(0,1,-211*ISGN,-1,-1,0,0,PIM2,AM,.TRUE.)
05079
05080
05081 CALL TRALO4(KTO,PIM1,PIM1,AM)
05082 CALL FILHEP(0,1,-211*ISGN,-2,-2,0,0,PIM1,AM,.TRUE.)
05083
05084
05085 CALL TRALO4(KTO,PIPL,PIPL,AM)
05086 CALL FILHEP(0,1, 211*ISGN,-3,-3,0,0,PIPL,AM,.TRUE.)
05087
05088 ELSE IF (JAA.EQ.2) THEN
05089
05090
05091
05092
05093 CALL TRALO4(KTO,PIM2,PIM2,AM)
05094 CALL FILHEP(0,1,111,-1,-1,0,0,PIM2,AM,.TRUE.)
05095
05096
05097 CALL TRALO4(KTO,PIM1,PIM1,AM)
05098 CALL FILHEP(0,1,111,-2,-2,0,0,PIM1,AM,.TRUE.)
05099
05100
05101 CALL TRALO4(KTO,PIPL,PIPL,AM)
05102 CALL FILHEP(0,1,-211*ISGN,-3,-3,0,0,PIPL,AM,.TRUE.)
05103
05104 ENDIF
05105
05106 RETURN
05107 END
05108 SUBROUTINE DWLUKK (KTO,ISGN,PKK,PNU)
05109
05110
05111
05112
05113
05114
05115
05116
05117 REAL PKK(4),PNU(4)
05118 COMMON /TAUPOS/ NP1,NP2
05119
05120
05121
05122 IF (KTO.EQ.1) THEN
05123
05124 NPS=NP1
05125 ELSE
05126 NPS=NP2
05127 ENDIF
05128
05129
05130 CALL TRALO4 (KTO,PNU,PNU,AM)
05131 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
05132
05133
05134 CALL TRALO4 (KTO,PKK,PKK,AM)
05135 CALL FILHEP(0,1,-321*ISGN,NPS,NPS,0,0,PKK,AM,.TRUE.)
05136
05137 RETURN
05138 END
05139 SUBROUTINE DWLUKS(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
05140 COMMON / TAUKLE / BRA1,BRK0,BRK0B,BRKS
05141 REAL*4 BRA1,BRK0,BRK0B,BRKS
05142
05143
05144
05145
05146
05147
05148
05149
05150
05151
05152
05153 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
05154 COMMON /TAUPOS/ NP1,NP2
05155
05156
05157
05158 IF(KTO.EQ. 1) THEN
05159 NPS=NP1
05160 ELSE
05161 NPS=NP2
05162 ENDIF
05163
05164
05165 CALL TRALO4(KTO,PNU,PNU,AM)
05166 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
05167
05168
05169 CALL TRALO4(KTO,PKS,PKS,AM)
05170 CALL FILHEP(0,1,-323*ISGN,NPS,NPS,0,0,PKS,AM,.TRUE.)
05171
05172
05173 IF(JKST.EQ.10) THEN
05174
05175
05176
05177
05178 CALL TRALO4(KTO,PPI,PPI,AM)
05179 CALL FILHEP(0,1,-211*ISGN,-1,-1,0,0,PPI,AM,.TRUE.)
05180
05181 BRAN=BRK0B
05182 IF (ISGN.EQ.-1) BRAN=BRK0
05183
05184 CALL RANMAR(XIO,1)
05185 IF(XIO(1).GT.BRAN) THEN
05186 K0TYPE = 130
05187 ELSE
05188 K0TYPE = 310
05189 ENDIF
05190
05191 CALL TRALO4(KTO,PKK,PKK,AM)
05192 CALL FILHEP(0,1,K0TYPE,-2,-2,0,0,PKK,AM,.TRUE.)
05193
05194 ELSE IF(JKST.EQ.20) THEN
05195
05196
05197
05198
05199 CALL TRALO4(KTO,PPI,PPI,AM)
05200 CALL FILHEP(0,1,111,-1,-1,0,0,PPI,AM,.TRUE.)
05201
05202
05203 CALL TRALO4(KTO,PKK,PKK,AM)
05204 CALL FILHEP(0,1,-321*ISGN,-2,-2,0,0,PKK,AM,.TRUE.)
05205
05206 ENDIF
05207
05208 RETURN
05209 END
05210 SUBROUTINE DWLNEW(KTO,ISGN,PNU,PWB,PNPI,MODE)
05211
05212
05213
05214
05215
05216
05217
05218
05219
05220 PARAMETER (NMODE=86,NM1=0,NM2=11,NM3=19,NM4=22,NM5=21,NM6=13)
05221
05222 COMMON / TAUDCD /IDFFIN(9,NMODE),MULPIK(NMODE)
05223
05224 & ,NAMES
05225 COMMON /TAUPOS/ NP1,NP2
05226 CHARACTER NAMES(NMODE)*31
05227 REAL PNU(4),PWB(4),PNPI(4,9)
05228 REAL PPI(4)
05229
05230 JNPI=MODE-7
05231
05232 IF(KTO.EQ. 1) THEN
05233 NPS=NP1
05234 ELSE
05235 NPS=NP2
05236 ENDIF
05237
05238
05239 CALL TRALO4(KTO,PNU,PNU,AM)
05240 CALL FILHEP(0,1,16*ISGN,NPS,NPS,0,0,PNU,AM,.TRUE.)
05241
05242
05243 CALL TRALO4(KTO,PWB,PWB,AM)
05244 CALL FILHEP(0,1,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
05245
05246
05247
05248
05249 ND=MULPIK(JNPI)
05250 DO I=1,ND
05251
05252 KFPI=LUNPIK(IDFFIN(I,JNPI),-ISGN)
05253
05254
05255
05256 DO J=1,4
05257 PPI(J)=PNPI(J,I)
05258 END DO
05259 CALL TRALO4(KTO,PPI,PPI,AM)
05260 CALL FILHEP(0,1,KFPI,-I,-I,0,0,PPI,AM,.TRUE.)
05261 END DO
05262
05263 RETURN
05264 END
05265
05266
05267 FUNCTION AMAST(PP)
05268
05269
05270
05271
05272
05273 IMPLICIT REAL*8 (A-H,O-Z)
05274 REAL*8 PP(4)
05275 AAA=PP(4)**2-PP(3)**2-PP(2)**2-PP(1)**2
05276
05277 IF(AAA.NE.0.0) AAA=AAA/SQRT(ABS(AAA))
05278 AMAST=AAA
05279 RETURN
05280 END
05281 FUNCTION AMAS4(PP)
05282
05283
05284
05285
05286
05287
05288 REAL PP(4)
05289 AAA=PP(4)**2-PP(3)**2-PP(2)**2-PP(1)**2
05290 IF(AAA.NE.0.0) AAA=AAA/SQRT(ABS(AAA))
05291 AMAS4=AAA
05292 RETURN
05293 END
05294 FUNCTION ANGXY(X,Y)
05295
05296
05297
05298
05299 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
05300 DATA PI /3.141592653589793238462643D0/
05301
05302 IF(ABS(Y).LT.ABS(X)) THEN
05303 THE=ATAN(ABS(Y/X))
05304 IF(X.LE.0D0) THE=PI-THE
05305 ELSE
05306 THE=ACOS(X/SQRT(X**2+Y**2))
05307 ENDIF
05308 ANGXY=THE
05309 RETURN
05310 END
05311 FUNCTION ANGFI(X,Y)
05312
05313
05314
05315
05316
05317 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
05318 DATA PI /3.141592653589793238462643D0/
05319
05320 IF(ABS(Y).LT.ABS(X)) THEN
05321 THE=ATAN(ABS(Y/X))
05322 IF(X.LE.0D0) THE=PI-THE
05323 ELSE
05324 THE=ACOS(X/SQRT(X**2+Y**2))
05325 ENDIF
05326 IF(Y.LT.0D0) THE=2D0*PI-THE
05327 ANGFI=THE
05328 END
05329 SUBROUTINE ROTOD1(PH1,PVEC,QVEC)
05330
05331
05332
05333
05334 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
05335 DIMENSION PVEC(4),QVEC(4),RVEC(4)
05336
05337 PHI=PH1
05338 CS=COS(PHI)
05339 SN=SIN(PHI)
05340 DO 10 I=1,4
05341 10 RVEC(I)=PVEC(I)
05342 QVEC(1)=RVEC(1)
05343 QVEC(2)= CS*RVEC(2)-SN*RVEC(3)
05344 QVEC(3)= SN*RVEC(2)+CS*RVEC(3)
05345 QVEC(4)=RVEC(4)
05346 RETURN
05347 END
05348 SUBROUTINE ROTOD2(PH1,PVEC,QVEC)
05349
05350
05351
05352
05353 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
05354 DIMENSION PVEC(4),QVEC(4),RVEC(4)
05355
05356 PHI=PH1
05357 CS=COS(PHI)
05358 SN=SIN(PHI)
05359 DO 10 I=1,4
05360 10 RVEC(I)=PVEC(I)
05361 QVEC(1)= CS*RVEC(1)+SN*RVEC(3)
05362 QVEC(2)=RVEC(2)
05363 QVEC(3)=-SN*RVEC(1)+CS*RVEC(3)
05364 QVEC(4)=RVEC(4)
05365 RETURN
05366 END
05367 SUBROUTINE ROTOD3(PH1,PVEC,QVEC)
05368
05369
05370
05371
05372 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
05373
05374 DIMENSION PVEC(4),QVEC(4),RVEC(4)
05375 PHI=PH1
05376 CS=COS(PHI)
05377 SN=SIN(PHI)
05378 DO 10 I=1,4
05379 10 RVEC(I)=PVEC(I)
05380 QVEC(1)= CS*RVEC(1)-SN*RVEC(2)
05381 QVEC(2)= SN*RVEC(1)+CS*RVEC(2)
05382 QVEC(3)=RVEC(3)
05383 QVEC(4)=RVEC(4)
05384 END
05385 SUBROUTINE BOSTR3(EXE,PVEC,QVEC)
05386
05387
05388
05389
05390
05391 REAL*4 PVEC(4),QVEC(4),RVEC(4)
05392
05393 DO 10 I=1,4
05394 10 RVEC(I)=PVEC(I)
05395 RPL=RVEC(4)+RVEC(3)
05396 RMI=RVEC(4)-RVEC(3)
05397 QPL=RPL*EXE
05398 QMI=RMI/EXE
05399 QVEC(1)=RVEC(1)
05400 QVEC(2)=RVEC(2)
05401 QVEC(3)=(QPL-QMI)/2
05402 QVEC(4)=(QPL+QMI)/2
05403 END
05404 SUBROUTINE BOSTD3(EXE,PVEC,QVEC)
05405
05406
05407
05408
05409
05410 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
05411 DIMENSION PVEC(4),QVEC(4),RVEC(4)
05412
05413 DO 10 I=1,4
05414 10 RVEC(I)=PVEC(I)
05415 RPL=RVEC(4)+RVEC(3)
05416 RMI=RVEC(4)-RVEC(3)
05417 QPL=RPL*EXE
05418 QMI=RMI/EXE
05419 QVEC(1)=RVEC(1)
05420 QVEC(2)=RVEC(2)
05421 QVEC(3)=(QPL-QMI)/2
05422 QVEC(4)=(QPL+QMI)/2
05423 RETURN
05424 END
05425 SUBROUTINE ROTOR1(PH1,PVEC,QVEC)
05426
05427
05428
05429
05430 REAL*4 PVEC(4),QVEC(4),RVEC(4)
05431
05432 PHI=PH1
05433 CS=COS(PHI)
05434 SN=SIN(PHI)
05435 DO 10 I=1,4
05436 10 RVEC(I)=PVEC(I)
05437 QVEC(1)=RVEC(1)
05438 QVEC(2)= CS*RVEC(2)-SN*RVEC(3)
05439 QVEC(3)= SN*RVEC(2)+CS*RVEC(3)
05440 QVEC(4)=RVEC(4)
05441 END
05442 SUBROUTINE ROTOR2(PH1,PVEC,QVEC)
05443
05444
05445
05446
05447 IMPLICIT REAL*4(A-H,O-Z)
05448 REAL*4 PVEC(4),QVEC(4),RVEC(4)
05449
05450 PHI=PH1
05451 CS=COS(PHI)
05452 SN=SIN(PHI)
05453 DO 10 I=1,4
05454 10 RVEC(I)=PVEC(I)
05455 QVEC(1)= CS*RVEC(1)+SN*RVEC(3)
05456 QVEC(2)=RVEC(2)
05457 QVEC(3)=-SN*RVEC(1)+CS*RVEC(3)
05458 QVEC(4)=RVEC(4)
05459 END
05460 SUBROUTINE ROTOR3(PHI,PVEC,QVEC)
05461
05462
05463
05464
05465 REAL*4 PVEC(4),QVEC(4),RVEC(4)
05466
05467 CS=COS(PHI)
05468 SN=SIN(PHI)
05469 DO 10 I=1,4
05470 10 RVEC(I)=PVEC(I)
05471 QVEC(1)= CS*RVEC(1)-SN*RVEC(2)
05472 QVEC(2)= SN*RVEC(1)+CS*RVEC(2)
05473 QVEC(3)=RVEC(3)
05474 QVEC(4)=RVEC(4)
05475 END
05476 SUBROUTINE SPHERD(R,X)
05477
05478
05479
05480
05481 REAL*8 R,X(4),PI,COSTH,SINTH
05482 REAL*4 RRR(2)
05483 DATA PI /3.141592653589793238462643D0/
05484
05485 CALL RANMAR(RRR,2)
05486 COSTH=-1+2*RRR(1)
05487 SINTH=SQRT(1 -COSTH**2)
05488 X(1)=R*SINTH*COS(2*PI*RRR(2))
05489 X(2)=R*SINTH*SIN(2*PI*RRR(2))
05490 X(3)=R*COSTH
05491 RETURN
05492 END
05493 SUBROUTINE ROTPOX(THET,PHI,PP)
05494 IMPLICIT REAL*8 (A-H,O-Z)
05495
05496
05497
05498
05499
05500 DIMENSION PP(4)
05501
05502 CALL ROTOD2(THET,PP,PP)
05503 CALL ROTOD3( PHI,PP,PP)
05504 RETURN
05505 END
05506 SUBROUTINE SPHERA(R,X)
05507
05508
05509
05510
05511
05512 REAL X(4)
05513 REAL*4 RRR(2)
05514 DATA PI /3.141592653589793238462643/
05515
05516 CALL RANMAR(RRR,2)
05517 COSTH=-1.+2.*RRR(1)
05518 SINTH=SQRT(1.-COSTH**2)
05519 X(1)=R*SINTH*COS(2*PI*RRR(2))
05520 X(2)=R*SINTH*SIN(2*PI*RRR(2))
05521 X(3)=R*COSTH
05522 RETURN
05523 END
05524 SUBROUTINE ROTPOL(THET,PHI,PP)
05525
05526
05527
05528
05529 REAL PP(4)
05530
05531 CALL ROTOR2(THET,PP,PP)
05532 CALL ROTOR3( PHI,PP,PP)
05533 RETURN
05534 END
05535 SUBROUTINE RANMAR(RVEC,LENV)
05536
05537
05538
05539
05540
05541
05542
05543
05544
05545
05546
05547
05548
05549
05550
05551 COMMON / INOUT / INUT,IOUT
05552 DIMENSION RVEC(*)
05553 COMMON/RASET1/U(97),C,I97,J97
05554 PARAMETER (MODCNS=1000000000)
05555 DATA NTOT,NTOT2,IJKL/-1,0,0/
05556
05557 IF (NTOT .GE. 0) GO TO 50
05558
05559
05560 IJKL = 54217137
05561 NTOT = 0
05562 NTOT2 = 0
05563 KALLED = 0
05564 GO TO 1
05565
05566 ENTRY RMARIN(IJKLIN, NTOTIN,NTOT2N)
05567
05568
05569
05570
05571
05572
05573
05574 IJKL = IJKLIN
05575 NTOT = MAX(NTOTIN,0)
05576 NTOT2= MAX(NTOT2N,0)
05577 KALLED = 1
05578
05579 1 CONTINUE
05580 IJ = IJKL/30082
05581 KL = IJKL - 30082*IJ
05582 I = MOD(IJ/177, 177) + 2
05583 J = MOD(IJ, 177) + 2
05584 K = MOD(KL/169, 178) + 1
05585 L = MOD(KL, 169)
05586 WRITE(IOUT,201) IJKL,NTOT,NTOT2
05587 201 FORMAT(1X,' RANMAR INITIALIZED: ',I10,2X,2I10)
05588 DO 2 II= 1, 97
05589 S = 0.
05590 T = .5
05591 DO 3 JJ= 1, 24
05592 M = MOD(MOD(I*J,179)*K, 179)
05593 I = J
05594 J = K
05595 K = M
05596 L = MOD(53*L+1, 169)
05597 IF (MOD(L*M,64) .GE. 32) S = S+T
05598 3 T = 0.5*T
05599 2 U(II) = S
05600 TWOM24 = 1.0
05601 DO 4 I24= 1, 24
05602 4 TWOM24 = 0.5*TWOM24
05603 C = 362436.*TWOM24
05604 CD = 7654321.*TWOM24
05605 CM = 16777213.*TWOM24
05606 I97 = 97
05607 J97 = 33
05608
05609
05610 DO 45 LOOP2= 1, NTOT2+1
05611 NOW = MODCNS
05612 IF (LOOP2 .EQ. NTOT2+1) NOW=NTOT
05613 IF (NOW .GT. 0) THEN
05614 WRITE (IOUT,'(A,I15)') ' RMARIN SKIPPING OVER ',NOW
05615 DO 40 IDUM = 1, NTOT
05616 UNI = U(I97)-U(J97)
05617 IF (UNI .LT. 0.) UNI=UNI+1.
05618 U(I97) = UNI
05619 I97 = I97-1
05620 IF (I97 .EQ. 0) I97=97
05621 J97 = J97-1
05622 IF (J97 .EQ. 0) J97=97
05623 C = C - CD
05624 IF (C .LT. 0.) C=C+CM
05625 40 CONTINUE
05626 ENDIF
05627 45 CONTINUE
05628 IF (KALLED .EQ. 1) RETURN
05629
05630
05631 50 CONTINUE
05632 DO 100 IVEC= 1, LENV
05633 UNI = U(I97)-U(J97)
05634 IF (UNI .LT. 0.) UNI=UNI+1.
05635 U(I97) = UNI
05636 I97 = I97-1
05637 IF (I97 .EQ. 0) I97=97
05638 J97 = J97-1
05639 IF (J97 .EQ. 0) J97=97
05640 C = C - CD
05641 IF (C .LT. 0.) C=C+CM
05642 UNI = UNI-C
05643 IF (UNI .LT. 0.) UNI=UNI+1.
05644
05645 IF (UNI .EQ. 0.) THEN
05646 UNI = TWOM24*U(2)
05647
05648 IF (UNI .EQ. 0.) UNI= TWOM24*TWOM24
05649 ENDIF
05650 RVEC(IVEC) = UNI
05651 100 CONTINUE
05652 NTOT = NTOT + LENV
05653 IF (NTOT .GE. MODCNS) THEN
05654 NTOT2 = NTOT2 + 1
05655 NTOT = NTOT - MODCNS
05656 ENDIF
05657 RETURN
05658
05659 ENTRY RMARUT(IJKLUT,NTOTUT,NTOT2T)
05660 IJKLUT = IJKL
05661 NTOTUT = NTOT
05662 NTOT2T = NTOT2
05663 RETURN
05664 END
05665 FUNCTION DILOGT(X)
05666
05667 IMPLICIT REAL*8(A-H,O-Z)
05668
05669 Z=-1.64493406684822
05670 IF(X .LT.-1.0) GO TO 1
05671 IF(X .LE. 0.5) GO TO 2
05672 IF(X .EQ. 1.0) GO TO 3
05673 IF(X .LE. 2.0) GO TO 4
05674 Z=3.2898681336964
05675 1 T=1.0/X
05676 S=-0.5
05677 Z=Z-0.5* LOG(ABS(X))**2
05678 GO TO 5
05679 2 T=X
05680 S=0.5
05681 Z=0.
05682 GO TO 5
05683 3 DILOGT=1.64493406684822
05684 RETURN
05685 4 T=1.0-X
05686 S=-0.5
05687 Z=1.64493406684822 - LOG(X)* LOG(ABS(T))
05688 5 Y=2.66666666666666 *T+0.66666666666666
05689 B= 0.00000 00000 00001
05690 A=Y*B +0.00000 00000 00004
05691 B=Y*A-B+0.00000 00000 00011
05692 A=Y*B-A+0.00000 00000 00037
05693 B=Y*A-B+0.00000 00000 00121
05694 A=Y*B-A+0.00000 00000 00398
05695 B=Y*A-B+0.00000 00000 01312
05696 A=Y*B-A+0.00000 00000 04342
05697 B=Y*A-B+0.00000 00000 14437
05698 A=Y*B-A+0.00000 00000 48274
05699 B=Y*A-B+0.00000 00001 62421
05700 A=Y*B-A+0.00000 00005 50291
05701 B=Y*A-B+0.00000 00018 79117
05702 A=Y*B-A+0.00000 00064 74338
05703 B=Y*A-B+0.00000 00225 36705
05704 A=Y*B-A+0.00000 00793 87055
05705 B=Y*A-B+0.00000 02835 75385
05706 A=Y*B-A+0.00000 10299 04264
05707 B=Y*A-B+0.00000 38163 29463
05708 A=Y*B-A+0.00001 44963 00557
05709 B=Y*A-B+0.00005 68178 22718
05710 A=Y*B-A+0.00023 20021 96094
05711 B=Y*A-B+0.00100 16274 96164
05712 A=Y*B-A+0.00468 63619 59447
05713 B=Y*A-B+0.02487 93229 24228
05714 A=Y*B-A+0.16607 30329 27855
05715 A=Y*A-B+1.93506 43008 6996
05716 DILOGT=S*T*(A-B)+Z
05717 RETURN
05718
05719
05720
05721 END
