1 /* copyright(c) 1991-2018 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public 6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public 15 license along with the gnu c library;
if not, see
16 <http://www.gnu.org/licenses/>. */
19 /* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27 /* glibc
's intent is to support the IEC 559 math functionality, real 28 and complex. If the GCC (4.9 and later) predefined macros 29 specifying compiler intent are available, use them to determine 30 whether the overall intent is to support these features; otherwise, 31 presume an older compiler has intent to support these features and 32 define these macros by default. */ 36 /* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is 37 synchronized with ISO/IEC 10646:2017, fifth edition, plus 38 the following additions from Amendment 1 to the fifth edition: 41 - 3 additional Zanabazar Square characters */ 43 SUBROUTINE TAUOLA(MODE,KEYPOL) 44 C ************************************* 45 C general tauola interface, should work in every case until 46 C hepevt is OK, does not check if hepevt is 'clean
' 47 C in particular will decay decayed taus... 48 C only longitudinal spin effects are included. 49 C in W decay v-a vertex is assumed 50 C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001 51 C this is the hepevt class in old style. No d_h_ class pre-name 53 PARAMETER (NMXHEP=10000) 54 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 55 INTEGER nevhep,nhep,isthep,idhep,jmohep, 58 $ nevhep, ! serial number 59 $ nhep, ! number of particles 60 $ isthep(nmxhep), ! status code 61 $ idhep(nmxhep), ! particle ident KF 62 $ jmohep(2,nmxhep), ! parent particles 63 $ jdahep(2,nmxhep), ! childreen particles 64 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 65 $ vhep(4,nmxhep) ! vertex [mm] 66 * ---------------------------------------------------------------------- 69 $ qedrad(nmxhep) ! Photos flag 70 * ---------------------------------------------------------------------- 72 COMMON /TAUPOS/ NP1, NP2 73 REAL*4 PHOI(4),PHOF(4) 74 double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4) 75 COMMON / MOMDEC / Q1,Q2,P1,P2,P3,P4 76 * tauola, photos and jetset overall switches 77 COMMON /LIBRA/ JAK1,JAK2,ITDKRC,IFPHOT,IFHADM,IFHADP 80 common /pseudocoup/ csc,ssc 83 COMMON / INOUT / INUT,IOUT 85 C to switch tau polarization OFF in taus 86 DIMENSION POL1(4), POL2(4) 87 double precision POL1x(4), POL2x(4) 89 DATA POL1 /0.0,0.0,0.0,0.0/ 90 DATA POL2 /0.0,0.0,0.0,0.0/ 91 DATA PI /3.141592653589793238462643D0/ 93 C store decay vertexes 94 DIMENSION IMOTHER (20) 97 C store daughter pointers 99 COMMON /ISONS_TAU/ISON(2) 103 C *********************** 105 JAK1 = 0 ! decay mode first tau 106 JAK2 = 0 ! decay mode second tau 107 ITDKRC=1.0 ! switch of radiative corrections in decay 108 IFPHOT=1.0 ! PHOTOS switch 111 POL=1.0 ! tau polarization dipswitch must be 1 or 0 121 C couplings of the 'pseudoscalar higgs
' as in CERN-TH/2003-166 123 xmtau=1.777 ! tau mass 124 xmh=120 ! higgs boson mass 125 betah=sqrt(1d0-4*xmtau**2/xmh**2) 128 C write(*,*) ' scalar component=
',csc,' pseudo-scalar component=
',ssc 129 .EQ.
IF (IFPHOT1) CALL PHOINI ! this if PHOTOS was not initialized earlier 130 CALL INIETC(JAK1,JAK2,ITDKRC,IFPHOT) 134 C activation of pi0 and eta decays: (1) means on, (0) off 138 CALL TAUPI0(-1,1,ION) 140 WRITE(IOUT,7001) pol,psi,ION(1),ION(2),ION(3) 141 .EQ.
ELSEIF(MODE0) THEN 142 C *********************** 144 C..... find tau-s and fill common block /TAUPOS/ 145 C this is to avoid LUND history fillings. This call is optional 146 CALL PHYFIX(NSTOP,NSTART) 147 C clear mothers of the previous event 155 C ... and to find mothers giving taus. 157 C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0) 159 .EQ..AND..EQ..AND.
IF(ABS(IDHEP(I))KFTAUISTHEP(I)1 160 .GE..OR..LT.
$ (ISTHEP(I)125ISTHEP(I)120)) THEN 162 .EQ.
DO WHILE (ABS(IDHEP(IMOTH))KFTAU) ! KEEP WALKING UP 163 IMOTH=JMOHEP(1,IMOTH) 165 .EQ..OR.
IF (ISTHEP(IMOTH)3 166 .GE..AND..LE.
$ (ISTHEP(IMOTH)120ISTHEP(IMOTH)125)) THEN 167 DO J=NSTART,NHEP ! WE HAVE WALKED INTO HARD RECORD 168 .EQ..AND.
IF (IDHEP(J)IDHEP(IMOTH) 169 .EQ..AND.
$ JMOHEP(1,J)IMOTH 170 .EQ.
$ ISTHEP(J)2) THEN 180 .EQ.
IF(JMOTHIMOTHER(II)) GOTO 9999 189 C ... taus of every mother are treated in this main loop 198 C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT.. 200 .EQ.
IF (IDHEP(JMOHEP(1,IM0))IDHEP(IM0)) IM0=JMOHEP(1,IM0) 203 .EQ..OR.
IF (ISTHEP(I)3 204 .GE..AND..LE.
$ (ISTHEP(I)120ISTHEP(I)125)) THEN ! HARD RECORD 208 .EQ..OR..EQ.
DO WHILE (IDHEP(IMOTH)IDHEP(I)ABS(IDHEP(IMOTH))KFTAU) ! KEEP WALKING UP 209 IMOTH=JMOHEP(1,IMOTH) 211 .EQ..OR..EQ..AND..EQ.
IF ((IMOTHIM0IMOTHIM)ISEL-1) THEN 214 .EQ..OR..EQ..AND..EQ.
ELSEIF ((IMOTHIM0IMOTHIM)ISEL0) THEN 216 .NE..AND..NE..AND..EQ.
ELSEIF ((IMOTHIM0IMOTHIM)ISEL0) THEN 226 C ... we correct HEPEVT (fix developped with Catherine BISCARAT) 227 .EQ.
c IF (JDAHEP(2,IM)0) THEN ! ID of second daughter was missing 229 c DO I=JDAHEP(1,IM)+1,NHEP ! OK lets look for it 230 .EQ..AND..EQ.
c IF (JMOHEP(1,I)IMISECU1) THEN ! we have found one 232 .EQ..AND..NE.
c ELSEIF (JMOHEP(1,I)IMISECU1) THEN ! we have found one after there 233 c JDAHEP(2,IM)=0 ! was something else, lets kill game 235 .NE.
c IF (JMOHEP(1,I)IM) ISECU=0 ! other stuff starts 239 C ... we check whether there are just two or more tau-likes 241 .EQ..OR..EQ.
IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NCOUNT=NCOUNT+1 242 .EQ..OR..EQ.
IF(IDHEP(I) KFTAUIDHEP(I) KFNUE) NCOUNT=NCOUNT+1 245 C ... if there will be more we will come here again 249 DO I=MAX(NP1+1,ISON(1)),ISON(2) 251 .EQ..OR..EQ.
IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NP1=I 254 DO I=MAX(NP2+1,ISON(1)),ISON(2) 256 .EQ..OR..EQ.
IF(IDHEP(I) KFTAUIDHEP(I) KFNUE) NP2=I 259 P1(I)= PHEP(I,NP1) !momentum of tau+ 260 P2(I)= PHEP(I,NP2) !momentum of tau- 267 .EQ.
IF(KEYPOL1) THEN 268 c.....include polarisation effect 271 .EQ..OR..EQ..OR.
IF(IDHEP(IM)KFHIGGS(1)IDHEP(IM)KFHIGGS(2) 272 .EQ.
$ IDHEP(IM)KFHIGGS(3)) THEN ! case of Higgs 273 .LT.
IF(RRR(1)0.5) THEN 280 .EQ..OR..EQ.
ELSEIF((IDHEP(IM)KFZ0)(IDHEP(IM)KFGAM)) THEN ! case of gamma/Z 281 C there is no angular dependence in gamma/Z polarization 282 C there is no s-dependence in gamma/Z polarization at all 283 C there is even no Z polarization in any form 284 C main reason is that nobody asked ... 285 C but it is prepared and longitudinal correlations 286 C can be included up to KORALZ standards 288 POLZ0=PLZAPX(.true.,IM,NP1,NP2) 289 .LT.
IF(RRR(1)POLZ0) THEN 296 .EQ.
ELSEIF(IDHEP(NP1)-IDHEP(NP2))THEN ! undef orig. only s-dep poss. 297 POLZ0=PLZAPX(.true.,IM,NP1,NP2) 298 .LT.
IF(RRR(1)POLZ0) THEN 305 .EQ.
ELSEIF(ABS(IDHEP(IM))KFHIGCH) THEN ! case of charged Higgs 308 ELSE ! case of W+ or W- 312 c.....include polarisation effect 315 .EQ..OR..EQ..OR.
IF(IDHEP(IM)KFHIGGS(1)IDHEP(IM)KFHIGGS(2) 316 .EQ.
$ IDHEP(IM)KFHIGGS(3)) THEN 317 .EQ..AND.
IF(IDHEP(NP1)-KFTAU 318 .LE..OR..GT..AND.
$ (JDAHEP(1,NP1)NP1JDAHEP(1,NP1)NHEP) 319 .EQ..AND.
$ IDHEP(NP2) KFTAU 320 .LE..OR..GT.
$ (JDAHEP(1,NP2)NP2JDAHEP(1,NP2)NHEP) 322 .EQ.
IF (IDHEP(IM)KFHIGGS(1)) THEN 324 .EQ.
ELSEIF (IDHEP(IM)KFHIGGS(2)) THEN 326 .EQ.
ELSEIF (IDHEP(IM)KFHIGGS(3)) THEN 329 WRITE(*,*) 'warning from tauola:
' 330 WRITE(*,*) 'i stop this run, wrong idhep(im)=
',IDHEP(IM) 333 CALL SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2) 334 .EQ.
IF (IFPHOT1) CALL PHOTOS(IM) ! Bremsstrahlung in Higgs decay 335 ! AFTER adding taus !! 340 .EQ..AND.
IF(IDHEP(NP1)-KFTAU 341 .LE..OR..GT.
$ (JDAHEP(1,NP1)NP1JDAHEP(1,NP1)NHEP)) THEN 342 C here check on if NP1 was not decayed should be verified 344 .EQ.
IF (IFPHOT1) CALL PHOTOS(NP1) 348 .EQ..AND.
IF(IDHEP(NP2) KFTAU 349 .LE..OR..GT.
$ (JDAHEP(1,NP2)NP2JDAHEP(1,NP2)NHEP)) THEN 350 C here check on if NP2 was not decayed should be added 352 .EQ.
IF (IFPHOT1) CALL PHOTOS(NP2) 357 .GT.
IF (NCOUNT0) GOTO 666 360 .EQ.
ELSEIF(MODE1) THEN 361 C *********************** 364 CALL DEKAY(100,POL1x) 368 7001 FORMAT(///1X,15(5H*****) 369 $ /,' *
', 25X,'*****tauola universal interface: ******
',9X,1H*, 370 $ /,' *
', 25X,'*****version 1.22, april 2009 (gfort)**
',9X,1H*, 371 $ /,' *
', 25X,'**authors: p. golonka, b. kersevan, ***
',9X,1H*, 372 $ /,' *
', 25X,'**t. pierzchala, e. richter-was, ******
',9X,1H*, 373 $ /,' *
', 25X,'****** z. was, m. worek ***************
',9X,1H*, 374 $ /,' *
', 25X,'**useful discussions, in particular ***
',9X,1H*, 375 $ /,' *
', 25X,'*with c. biscarat and s. slabospitzky**
',9X,1H*, 376 $ /,' *
', 25X,'****** are warmly acknowledged ********
',9X,1H*, 377 $ /,' *
', 25X,' ',9X,1H*, 378 $ /,' *
', 25X,'********** initialization ************
',9X,1H*, 379 $ /,' *
',F20.5,5X,'tau polarization switch must be 1 or 0
',9X,1H*, 380 $ /,' *
',F20.5,5X,'higs scalar/pseudo mix cern-th/2003-166
',9X,1H*, 381 $ /,' *
',I10, 15X,'pi0 decay switch must be 1 or 0
',9X,1H*, 382 $ /,' *
',I10, 15X,'eta decay switch must be 1 or 0
',9X,1H*, 383 $ /,' *
',I10, 15X,'k0s decay switch must be 1 or 0
',9X,1H*, 386 7002 FORMAT(///1X,15(5H*****) 387 $ /,' *
', 25X,'*****tauola universal interface: ******
',9X,1H*, 388 $ /,' *
', 25X,'*****version 1.22, april 2009 (gfort)**
',9X,1H*, 389 $ /,' *
', 25X,'**authors: p. golonka, b. kersevan, ***
',9X,1H*, 390 $ /,' *
', 25X,'**t. pierzchala, e. richter-was, ******
',9X,1H*, 391 $ /,' *
', 25X,'****** z. was, m. worek ***************
',9X,1H*, 392 $ /,' *
', 25X,'**useful discussions, in particular ***
',9X,1H*, 393 $ /,' *
', 25X,'*with c. biscarat and s. slabospitzky**
',9X,1H*, 394 $ /,' *
', 25X,'****** are warmly acknowledged ********
',9X,1H*, 395 $ /,' *
', 25X,'****** end of
MODULE operation ********
',9X,1H*, 400 SUBROUTINE SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2) 402 REAL*8 HH1,HH2,wthiggs 403 DIMENSION POL1(4), POL2(4),HH1(4),HH2(4), RRR(1) 404 ! CALL DEXAY(1,POL1) ! Kept for tests 405 ! CALL DEXAY(2,POL2) ! Kept for tests 411 wt=wthiggs(IFPSEUDO,HH1,HH2) 412 .GT.
IF (RRR(1)WT) GOTO 10 418 FUNCTION wthiggs(IFPSEUDO,HH1,HH2) 420 common /pseudocoup/ csc,ssc 423 REAL*8 HH1(4),HH2(4),R(4,4),wthiggs 431 R(4,4)= 1D0 ! unpolarized part 432 R(3,3)=-1D0 ! longitudinal 437 R(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2) 438 R(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2) 439 R(1,2)=2*csc*ssc/(csc**2+ssc**2) 440 R(2,1)=-2*csc*ssc/(csc**2+ssc**2) 450 WTHIGGS=WTHIGGS+R(K,L)*HH1(K)*HH2(L) 456 FUNCTION PLZAPX(HOPEin,IM0,NP1,NP2) 457 C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block. 458 C the purpose of this routine is to calculate polarization of Z along tau direction. 459 C this is highly non-trivial due to necessity of reading infromation from hard process 460 C history in HEPEVT, which is often written not up to the gramatic rules. 461 REAL*8 PLZAP0,SVAR,COSTHE,sini,sfin,ZPROP2, 462 $ P1(4),P2(4),Q1(4),Q2(4),QQ(4),PH(4),PD1(4),PD2(4), 463 $ PQ1(4),PQ2(4),PB(4),PA(4) 467 C this is the hepevt class in old style. No d_h_ class pre-name 469 PARAMETER (NMXHEP=10000) 470 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 471 INTEGER nevhep,nhep,isthep,idhep,jmohep, 474 $ nevhep, ! serial number 475 $ nhep, ! number of particles 476 $ isthep(nmxhep), ! status code 477 $ idhep(nmxhep), ! particle ident KF 478 $ jmohep(2,nmxhep), ! parent particles 479 $ jdahep(2,nmxhep), ! childreen particles 480 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 481 $ vhep(4,nmxhep) ! vertex [mm] 482 * ---------------------------------------------------------------------- 485 $ qedrad(nmxhep) ! Photos flag 486 * ---------------------------------------------------------------------- 489 C(BPK)--> BROTHERS OF TAU ALREADY FOUND 491 COMMON /ISONS_TAU/ISON(2) 494 C >> STEP 1: find where are particles in hepevent and pick them up 497 C sometimes shade Z of Z is its mother ... 500 C to protect against check on mother of beam particles. 502 .EQ.
IF (IDHEP(IM0)IDHEP(IM00)) IM=JMOHEP(1,IM0) 505 C find (host generator-level) incoming beam-bare-particles which form Z and co. 509 C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS 511 .EQ.
IF (ISTHEP(IM00)120) THEN 518 C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first
' in hepevt. 519 .EQ..AND..EQ.
IF (IMO10IMO20) THEN 522 .EQ.
IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES 526 .EQ.
IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES 529 C case when it was like qq --> tau+tau- gammas and qq were NOT 'first
' in hepevt. 531 .NE..AND..NE.
ELSEIF (IDHEP(IM)22IDHEP(IM)23) THEN 534 .EQ.
IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES 538 .EQ.
IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES 544 C and check if it really happened 545 .EQ.
IF (IMO10) HOPE=.FALSE. 546 .EQ.
IF (IMO20) HOPE=.FALSE. 547 .EQ.
IF (IMO1IMO2) HOPE=.FALSE. 551 Q1(I)= PHEP(I,NP1) !momentum of tau+ 552 Q2(I)= PHEP(I,NP2) !momentum of tau- 555 C corrections due to possible differences in 4-momentum of shadow vs true Z. 557 .EQ..AND.
IF (IMJMOHEP(1,IM0) 558 .EQ..OR..EQ.
$ (IDHEP(IM)22IDHEP(IM)23)) THEN 565 CALL BOSTDQ( 1,PA, Q1, Q1) 566 CALL BOSTDQ( 1,PA, Q2, Q2) 567 CALL BOSTDQ(-1,PB, Q1, Q1) 568 CALL BOSTDQ(-1,PB, Q2, Q2) 573 QQ(I)= Q1(I)+Q2(I) !momentum of Z 574 IF (HOPE) P1(I)=PHEP(I,IMO1) !momentum of beam1 575 IF (HOPE) P2(I)=PHEP(I,IMO2) !momentum of beam2 580 ! These momenta correspond to quarks, gluons photons or taus 583 IF (HOPE) IDFP1=IDHEP(IMO1) 584 IF (HOPE) IDFP2=IDHEP(IMO2) 586 SVAR=QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2 588 C options which may be useful in some cases of two heavy boson production 589 C need individual considerations. To be developed. 590 C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam 591 C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z 592 PLZAPX=0.5 ! pure gamma 596 C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles 598 C let us define beginning and end of particles which are produced in parallel to Z 599 C in principle following should work 601 C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS 605 INBR=IM ! OK, HARD RECORD Z/GAMMA 606 .EQ.
IF (IFFULL1) INBR=NP1 ! OK, NO Z/GAMMA 607 .EQ.
IF (IDHEP(JMOHEP(1,INBR))IDHEP(INBR)) INBR=JMOHEP(1,INBR) ! FORCE HARD RECORD 609 .EQ..OR..EQ.
IF(NX10NX20) THEN 613 .EQ.
IF(JMOHEP(1,INBR-K)JMOHEP(1,INBR)) THEN 622 .EQ.
IF(JMOHEP(1,K)JMOHEP(1,INBR)) THEN 631 C case of annihilation of two bosons is hopeles 632 .GE..AND..GE.
IF (ABS(IDFP1)20ABS(IDFP2)20) HOPE=.FALSE. 633 C case of annihilation of non-matching flavors is hopeless 634 .LE..AND..LE..AND..NE.
IF (ABS(IDFP1)20ABS(IDFP2)20IDFP1+IDFP20) 637 C options which may be useful in some cases of two heavy boson production 638 C need individual considerations. To be developed. 639 C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam 640 C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z 641 PLZAPX=0.5 ! pure gamma 644 .LT.
IF (ABS(IDFP1)20) IDE= IDFP1 645 .LT.
IF (ABS(IDFP2)20) IDE=-IDFP2 649 C >> STEP 3 we combine gluons, photons into incoming effective beams 652 C in the following we will ignore the possibility of photon emission from taus 653 C however at certain moment it will be necessary to take care of 665 IFLAV=min(ABS(IDFP1),ABS(IDFP2)) 667 *-------------------------------------------------------------------------- 668 c IFLAV=min(ABS(IDFP1),ABS(IDFP2)) 669 c that means that always quark or lepton i.e. process like 670 c f g(gamma) --> f Z0 --> tau tau 671 c we glue fermions to effective beams that is f f --> Z0 --> tau tau 672 c with gamma/g emission from initial fermion. 673 *--------------------------------------------------------------------------- 675 .GE.
IF (ABS(IDFP1)20) THEN 678 .EQ.
IF (ABS(IDP)IFLAV) THEN 686 .GE.
IF (ABS(IDFP2)20) THEN 689 .EQ.
IF (ABS(IDP)IFLAV) THEN 697 C if first beam was boson: gluing 699 .GE.
IF (ABS(IDFP1)20) THEN 703 xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2 704 $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2) 705 xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2 706 $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2) 707 .LT.
IF (XM1XM2) THEN 718 C if second beam was boson: gluing 721 .GE.
IF (ABS(IDFP2)20) THEN 725 xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2 726 $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2) 727 xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2 728 $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2) 729 .LT.
IF (XM1XM2) THEN 745 .EQ.
IF (IDHEP(JMOHEP(1,NP1))IDHEP(NP1)) NPH1=JMOHEP(1,NP1) ! SHOULD PUT US IN HARD REC. 746 .EQ.
IF (IDHEP(JMOHEP(1,NP2))IDHEP(NP2)) NPH2=JMOHEP(1,NP2) ! SHOULD PUT US IN HARD REC. 750 .NE..AND..NE..AND.
IF (ABS(IDHEP(K))IFLAVKIM 752 .NE..AND..NE.
$ KNPH1KNPH2) THEN 754 .EQ..AND..EQ.
IF(IDHEP(K)22IFFULL1) THEN 758 xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2 759 $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2) 760 xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2 761 $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2) 762 xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2 763 $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2) 764 xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2 765 $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2) 768 sini=abs((PD1(4)+PD2(4)-PH(4))**2-(PD1(3)+PD2(3)-PH(3))**2 769 $ -(PD1(2)+PD2(2)-PH(2))**2-(PD1(1)+PD2(1)-PH(1))**2) 770 sfin=abs((PD1(4)+PD2(4) )**2-(PD1(3)+PD2(3) )**2 771 $ -(PD1(2)+PD2(2) )**2-(PD1(1)+PD2(1) )**2) 781 XM=MIN(XM1,XM2,XM3,XM4) 786 .EQ.
ELSEIF (XM2XM) THEN 790 .EQ.
ELSEIF (XM3XM) THEN 803 xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2 804 $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2) 805 xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2 806 $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2) 807 .LT.
IF (XM1XM2) THEN 822 C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in 823 c >> effective outcoming taus 825 C let us define beginning and end of particles which are produced in 830 C find outcoming particles which come from Z 835 C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER 836 .EQ..OR..EQ.
IF (ABS(IDHEP(IM0))22abs(IDHEP(IM0))23) THEN 838 .EQ.
IF(ABS(IDHEP(K))22) THEN 845 xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2 846 $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2) 847 xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2 848 $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2) 867 *------------------------------------------------------------------------ 870 C out of effective momenta we calculate COSTHE and later polarization 871 CALL ANGULU(PD1,PD2,Q1,Q2,COSTHE) 873 PLZAPX=PLZAP0(IDE,IDFQ1,SVAR,COSTHE) 876 SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE) 877 REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4) 878 C take effective beam which is less massive, it should be irrelevant 879 C but in case HEPEVT is particulary dirty may help. 880 C this routine calculate reduced system transver and cosine of scattering 883 XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2) 884 XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2) 885 .LT.
IF (XM1XM2) THEN 896 C calculate space like part of P (in Z restframe) 902 XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2) 904 QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1) 906 QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2 909 PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1) 911 P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2 914 PXP =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2) 915 QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2) 916 PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4) 917 COSTHE=PXQT/PXP/QTXQT 921 FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0) 922 C this function calculates probability for the helicity +1 +1 configuration 923 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle 924 REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN 927 .LT.
C >>>>> IF (IDE*IDF0) COSTHE=-COSTH0 ! this is probably not needed ID 928 C >>>>> of first beam is used by T_GIVIZ0 including sign 931 CALL INITWK(IDE,IDF,SVAR) 933 CALL INITWK(-IDE,-IDF,SVAR) 935 PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0) 936 $ /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0)) 940 FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB) 941 C ---------------------------------------------------------------------- 942 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME 943 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER 944 C EXAMPLE OF THE METHOD APPLIED THERE 945 C INPUT PARAMETERS ARE: SVAR -- transfer 946 C COSTHE -- cosine of angle between tau+ and 1st beam 947 C TA,TB -- helicity states of tau+ tau- 949 C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT 950 C ---------------------------------------------------------------------- 951 IMPLICIT REAL*8(A-H,O-Z) 952 COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF 953 REAL*8 ENE ,AMIN,AMFIN 954 COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF 955 & ,XUPGI ,XUPZI ,XUPGF ,XUPZF 956 & ,NDIAG0,NDIAGA,KEYA,KEYZ 957 & ,ITCE,JTCE,ITCF,JTCF,KOLOR 958 REAL*8 SS,POLN,T3E,QE,T3F,QF 959 & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2) 961 C===================================================================== 962 COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ 963 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ 964 C SWSQ = sin2 (theta Weinberg) 965 C AMW,AMZ = W & Z boson masses respectively 966 C AMH = the Higgs mass 967 C AMTOP = the top mass 969 COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2) 970 COMPLEX*16 XUPZFP(2),XUPZIP(2) 971 COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2) 972 COMPLEX*16 PROPA,PROPZ 974 COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF 975 COMPLEX*16 XTHING,XVE,XVF,XVEF 976 DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/ 979 DATA SVAR0,COST0 /-5.D0,-6.D0/ 980 DATA PI /3.141592653589793238462643D0/ 981 DATA SEPS1,SEPS2 /0D0,0D0/ 983 C MEMORIZATION ========================================================= 984 .NE..OR..NE..OR..NE.
IF ( MODEMODE0SVARSVAR0COSTHECOST0 985 .OR..NE.
$ IDE0IDE)THEN 993 SINTHE=SQRT(1.D0-COSTHE**2) 994 BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR)) 995 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR. 996 XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2)) 997 XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2)) 998 XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2)) 999 XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2)) 1000 C FINAL STATE VECTOR COUPLING 1001 XUPF =0.5D0*(XUPZF(1)+XUPZF(2)) 1002 XUPI =0.5D0*(XUPZI(1)+XUPZI(2)) 1006 PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ) 1007 .EQ.
IF (KEYGSW0) PROPZ=0.D0 1010 REGULA= (3-2*I)*(3-2*J) + COSTHE 1011 REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR) 1012 APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA) 1013 AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA 1014 ABORN(I,J)=APHOT(I,J)+AZETT(I,J) 1015 APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM 1016 AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM 1017 ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J) 1022 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS 1023 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.) 1024 C* HELICITY CONSERVATION EXPLICITLY OBEYED 1032 FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0 1033 FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB) 1034 FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB) 1036 BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR 1038 .GE.
IF (MODE1) THEN 1039 BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM 1045 .LT.
IF(FUNT0.D0) FUNT=BORN 1048 .GT.
IF (SVAR4D0*AMFIN**2) THEN 1049 C PHASE SPACE THRESHOLD FACTOR 1050 THRESH=SQRT(1-4D0*AMFIN**2/SVAR) 1051 T_BORN= FUNT*SVAR**2*THRESH 1056 C ZW HERE WAS AN ERROR 19. 05. 1989 1057 ! write(*,*) 'kkkk
',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF 1058 ! write(*,*) 'kkkk x
',svar,costhe,TA,TB,T_BORN 1061 SUBROUTINE INITWK(IDEX,IDFX,SVAR) 1062 ! initialization routine coupling masses etc. 1063 IMPLICIT REAL*8 (A-H,O-Z) 1064 COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF 1065 REAL*8 ENE ,AMIN,AMFIN 1066 COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF 1067 & ,XUPGI ,XUPZI ,XUPGF ,XUPZF 1068 & ,NDIAG0,NDIAGA,KEYA,KEYZ 1069 & ,ITCE,JTCE,ITCF,JTCF,KOLOR 1070 REAL*8 SS,POLN,T3E,QE,T3F,QF 1071 & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2) 1072 COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ 1073 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ 1074 C SWSQ = sin2 (theta Weinberg) 1075 C AMW,AMZ = W & Z boson masses respectively 1076 C AMH = the Higgs mass 1077 C AMTOP = the top mass 1085 .EQ.
IF (IDFX 15) then 1086 IDF=2 ! denotes tau +2 tau- 1087 AMFIN=1.77703 !this mass is irrelevant if small, used in ME only 1088 .EQ.
ELSEIF (IDFX-15) then 1089 IDF=-2 ! denotes tau -2 tau- 1090 AMFIN=1.77703 !this mass is irrelevant if small, used in ME only 1092 WRITE(*,*) 'initwk: wrong idfx
' 1096 .EQ.
IF (IDEX 11) then !electron 1099 .EQ.
ELSEIF (IDEX-11) then !positron 1102 .EQ.
ELSEIF (IDEX 13) then !mu+ 1105 .EQ.
ELSEIF (IDEX-13) then !mu- 1108 .EQ.
ELSEIF (IDEX 1) then !d 1111 .EQ.
ELSEIF (IDEX- 1) then !d~ 1114 .EQ.
ELSEIF (IDEX 2) then !u 1117 .EQ.
ELSEIF (IDEX- 2) then !u~ 1120 .EQ.
ELSEIF (IDEX 3) then !s 1123 .EQ.
ELSEIF (IDEX- 3) then !s~ 1126 .EQ.
ELSEIF (IDEX 4) then !c 1129 .EQ.
ELSEIF (IDEX- 4) then !c~ 1132 .EQ.
ELSEIF (IDEX 5) then !b 1135 .EQ.
ELSEIF (IDEX- 5) then !b~ 1138 .EQ.
ELSEIF (IDEX 12) then !nu_e 1141 .EQ.
ELSEIF (IDEX- 12) then !nu_e~ 1144 .EQ.
ELSEIF (IDEX 14) then !nu_mu 1147 .EQ.
ELSEIF (IDEX- 14) then !nu_mu~ 1150 .EQ.
ELSEIF (IDEX 16) then !nu_tau 1153 .EQ.
ELSEIF (IDEX- 16) then !nu_tau~ 1158 WRITE(*,*) 'initwk: wrong idex
' 1162 C ---------------------------------------------------------------------- 1164 C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX 1166 C called by : KORALZ 1167 C ---------------------------------------------------------------------- 1172 CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM) 1173 CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM) 1177 XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) 1178 XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) 1179 CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR) 1180 CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR) 1184 XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ)) 1185 XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ)) 1196 SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR) 1197 C ---------------------------------------------------------------------- 1198 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION 1199 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK 1200 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE 1201 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY) 1202 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF) 1203 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE 1204 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS 1206 C called by : EVENTE, EVENTM, FUNTIH, ..... 1207 C ---------------------------------------------------------------------- 1208 IMPLICIT REAL*8(A-H,O-Z) 1210 .EQ..OR..GT.
IF(IDFERM0IABS(IDFERM)4) GOTO 901 1211 .NE.
IF(IABS(IHELIC)1) GOTO 901 1213 IDTYPE =IABS(IDFERM) 1215 LEPQUA=INT(IDTYPE*0.4999999D0) 1216 IUPDOW=IDTYPE-2*LEPQUA-1 1217 CHARGE =(-IUPDOW+2D0/3D0*LEPQUA)*IC 1218 SIZO3 =0.25D0*(IC-IH)*(1-2*IUPDOW) 1220 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS 1221 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) 1223 901 PRINT *,' stop in givizo: wrong params.
' 1226 SUBROUTINE PHYFIX(NSTOP,NSTART) 1227 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) 1229 C NSTOP NSTART : when PHYTIA history ends and event starts. 1233 .NE.
IF(K(I,1)21) THEN 1241 SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG) 1242 C ---------------------------------------------------------------------- 1243 C this subroutine fills one entry into the HEPEVT common 1244 C and updates the information for affected mother entries 1246 C written by Martin W. Gruenewald (91/01/28) 1248 C called by : ZTOHEP,BTOHEP,DWLUxy 1249 C ---------------------------------------------------------------------- 1251 C this is the hepevt class in old style. No d_h_ class pre-name 1252 C this is the hepevt class in old style. No d_h_ class pre-name 1254 PARAMETER (NMXHEP=10000) 1255 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 1256 INTEGER nevhep,nhep,isthep,idhep,jmohep, 1259 $ nevhep, ! serial number 1260 $ nhep, ! number of particles 1261 $ isthep(nmxhep), ! status code 1262 $ idhep(nmxhep), ! particle ident KF 1263 $ jmohep(2,nmxhep), ! parent particles 1264 $ jdahep(2,nmxhep), ! childreen particles 1265 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 1266 $ vhep(4,nmxhep) ! vertex [mm] 1267 * ---------------------------------------------------------------------- 1270 $ qedrad(nmxhep) ! Photos flag 1271 * ---------------------------------------------------------------------- 1277 C check address mode 1282 .GT.
ELSE IF (N0) THEN 1293 .LE..OR..GT.
IF ((IHEP0)(IHEPNMXHEP)) RETURN 1300 .LT.
IF(JMO10)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP 1302 .LT.
IF(JMO20)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP 1309 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations 1313 C FLAG FOR PHOTOS... 1317 DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP) 1320 C if there is a daughter at IHEP, mother entry at IP has decayed 1321 .EQ.
IF(ISTHEP(IP)1)ISTHEP(IP)=2 1323 C and daughter pointers of mother entry must be updated 1324 .EQ.
IF(JDAHEP(1,IP)0)THEN 1328 JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP)) 1337 FUNCTION IHEPDIM(DUM) 1338 C this is the hepevt class in old style. No d_h_ class pre-name 1339 C this is the hepevt class in old style. No d_h_ class pre-name 1341 PARAMETER (NMXHEP=10000) 1342 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 1343 INTEGER nevhep,nhep,isthep,idhep,jmohep, 1346 $ nevhep, ! serial number 1347 $ nhep, ! number of particles 1348 $ isthep(nmxhep), ! status code 1349 $ idhep(nmxhep), ! particle ident KF 1350 $ jmohep(2,nmxhep), ! parent particles 1351 $ jdahep(2,nmxhep), ! childreen particles 1352 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 1353 $ vhep(4,nmxhep) ! vertex [mm] 1354 * ---------------------------------------------------------------------- 1357 $ qedrad(nmxhep) ! Photos flag 1358 * ---------------------------------------------------------------------- 1363 IMPLICIT REAL*8(A-H,O-Z) 1364 COMPLEX*16 CPRZ0,CPRZ0M 1367 CPRZ0=DCMPLX((S-AMZ**2),S/AMZ*GAMMZ) 1369 ZPROP2=(ABS(CPRZ0M))**2 1372 SUBROUTINE TAUPI0(MODE,JAK,ION) 1373 C no initialization required. Must be called once after every: 1374 C 1) CALL DEKAY(1+10,...) 1375 C 2) CALL DEKAY(2+10,...) 1376 C 3) CALL DEXAY(1,...) 1377 C 4) CALL DEXAY(2,...) 1378 C subroutine to decay originating from TAUOLA's taus:
1379 c 1) etas(with
CALL taueta(jak))
1380 c 2) later pi0
's from taus. 1381 C 3) extensions to other applications possible. 1382 C this routine belongs to >tauola universal interface<, but uses 1383 C routines from >tauola< utilities as well. 25.08.2005 1384 C this is the hepevt class in old style. No d_h_ class pre-name 1386 PARAMETER (NMXHEP=10000) 1387 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 1388 INTEGER nevhep,nhep,isthep,idhep,jmohep, 1391 $ nevhep, ! serial number 1392 $ nhep, ! number of particles 1393 $ isthep(nmxhep), ! status code 1394 $ idhep(nmxhep), ! particle ident KF 1395 $ jmohep(2,nmxhep), ! parent particles 1396 $ jdahep(2,nmxhep), ! childreen particles 1397 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 1398 $ vhep(4,nmxhep) ! vertex [mm] 1399 * ---------------------------------------------------------------------- 1402 $ qedrad(nmxhep) ! Photos flag 1403 * ---------------------------------------------------------------------- 1406 C position of taus, must be defined by host program: 1407 COMMON /TAUPOS/ NP1,NP2 1409 REAL PHOT1(4),PHOT2(4) 1410 REAL*8 R,X(4),Y(4),PI0(4) 1411 INTEGER JEZELI(3),ION(3) 1414 .EQ.
IF (MODE-1) THEN 1420 .EQ.
IF (JEZELI(1)0) RETURN 1421 .EQ.
IF (JEZELI(2)1) CALL TAUETA(JAK) 1422 .EQ.
IF (JEZELI(3)1) CALL TAUK0S(JAK) 1423 C position of decaying particle: 1424 .EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN 1429 nhepM=nhep ! to avoid infinite loop 1430 DO K=JDAHEP(1,NPS),nhepM ! we search for pi0's from tau till eor.
1431 IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k)
THEN 1436 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1445 CALL bostdq(-1,pi0,x,x)
1446 CALL bostdq(-1,pi0,y,y)
1452 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1453 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1458 SUBROUTINE taueta(JAK)
1459 c
subroutine to decay etas
's from taus. 1460 C this routine belongs to tauola universal interface, but uses 1461 C routines from tauola utilities. Just flat phase space, but 4 channels. 1462 C it is called at the beginning of SUBR. TAUPI0(JAK) 1463 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005 1464 C this is the hepevt class in old style. No d_h_ class pre-name 1466 PARAMETER (NMXHEP=10000) 1467 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 1468 INTEGER nevhep,nhep,isthep,idhep,jmohep, 1471 $ nevhep, ! serial number 1472 $ nhep, ! number of particles 1473 $ isthep(nmxhep), ! status code 1474 $ idhep(nmxhep), ! particle ident KF 1475 $ jmohep(2,nmxhep), ! parent particles 1476 $ jdahep(2,nmxhep), ! childreen particles 1477 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 1478 $ vhep(4,nmxhep) ! vertex [mm] 1479 * ---------------------------------------------------------------------- 1482 $ qedrad(nmxhep) ! Photos flag 1483 * ---------------------------------------------------------------------- 1485 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU 1486 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 1487 * ,AMK,AMKZ,AMKST,GAMKST 1489 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU 1490 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 1491 * ,AMK,AMKZ,AMKST,GAMKST 1493 C position of taus, must be defined by host program: 1494 COMMON /TAUPOS/ NP1,NP2 1496 REAL RRR(1),BRSUM(3), RR(2) 1497 REAL PHOT1(4),PHOT2(4),PHOT3(4) 1498 REAL*8 X(4), Y(4), Z(4) 1500 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM 1502 XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c)) 1503 C position of decaying particle: 1504 .EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN 1509 nhepM=nhep ! to avoid infinite loop 1510 DO K=JDAHEP(1,NPS),nhepM ! we search for etas's from tau till eor.
1511 IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k)
THEN 1515 c eta cumulated branching ratios:
1517 brsum(2)=brsum(1)+0.319
1518 brsum(3)=brsum(2)+0.237
1521 IF (rrr(1).LT.brsum(1))
THEN 1523 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1532 CALL bostdq(-1,peta,x,x)
1533 CALL bostdq(-1,peta,y,y)
1539 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1540 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1542 IF(rrr(1).LT.brsum(2))
THEN 1549 ELSEIF(rrr(1).LT.brsum(3))
THEN 1566 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1569 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1570 c weight for flat phase space
1571 wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1572 & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1574 IF (rr(2).GT.wt)
GOTO 7
1576 ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2
1580 x(4)=sqrt(ru**2+xm1**2)
1581 y(4)=sqrt(ru**2+xm2**2)
1586 c generate momentum of that pair in rest frame of eta:
1587 ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1589 z(4)=sqrt(ru**2+am2**2)
1590 c and boost first two decay products to rest frame of eta.
1591 CALL bostdq(-1,z,x,x)
1592 CALL bostdq(-1,z,y,y)
1593 c redefine z(4) to 4-momentum of the last decay product:
1597 z(4)=sqrt(ru**2+xm3**2)
1598 c boost all to lab and move to real*4; also masses
1599 CALL bostdq(-1,peta,x,x)
1600 CALL bostdq(-1,peta,y,y)
1601 CALL bostdq(-1,peta,z,z)
1611 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1612 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1613 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1620 SUBROUTINE tauk0s(JAK)
1621 c
subroutine to decay k0s
's from taus. 1622 C this routine belongs to tauola universal interface, but uses 1623 C routines from tauola utilities. Just flat phase space, but 4 channels. 1624 C it is called at the beginning of SUBR. TAUPI0(JAK) 1625 C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005 1626 C this is the hepevt class in old style. No d_h_ class pre-name 1628 PARAMETER (NMXHEP=10000) 1629 REAL*8 phep, vhep ! to be real*4/ *8 depending on host 1630 INTEGER nevhep,nhep,isthep,idhep,jmohep, 1633 $ nevhep, ! serial number 1634 $ nhep, ! number of particles 1635 $ isthep(nmxhep), ! status code 1636 $ idhep(nmxhep), ! particle ident KF 1637 $ jmohep(2,nmxhep), ! parent particles 1638 $ jdahep(2,nmxhep), ! childreen particles 1639 $ phep(5,nmxhep), ! four-momentum, mass [GeV] 1640 $ vhep(4,nmxhep) ! vertex [mm] 1641 * ---------------------------------------------------------------------- 1644 $ qedrad(nmxhep) ! Photos flag 1645 * ---------------------------------------------------------------------- 1648 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU 1649 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 1650 * ,AMK,AMKZ,AMKST,GAMKST 1652 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU 1653 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 1654 * ,AMK,AMKZ,AMKST,GAMKST 1656 C position of taus, must be defined by host program: 1657 COMMON /TAUPOS/ NP1,NP2 1659 REAL RRR(1),BRSUM(3), RR(2) 1660 REAL PHOT1(4),PHOT2(4),PHOT3(4) 1661 REAL*8 X(4), Y(4), Z(4) 1663 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM 1665 XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c)) 1666 C position of decaying particle: 1667 .EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN 1672 nhepM=nhep ! to avoid infinite loop 1673 DO K=JDAHEP(1,NPS),nhepM ! we search for K0S's from tau till eor.
1674 IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k)
THEN 1680 c k0s cumulated branching ratios:
1683 brsum(3)=brsum(2)+0.237
1686 IF(rrr(1).LT.brsum(1))
THEN 1691 ELSEIF(rrr(1).LT.brsum(2))
THEN 1704 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1706 r=sqrt(abs(r**2-xm1**2))
1715 CALL bostdq(-1,peta,x,x)
1716 CALL bostdq(-1,peta,y,y)
1725 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1726 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)