SUBROUTINE DECAY6(AM0,AM3,AM4,AM5,PARAMA,PARAMB,PARAMC,AMPMX,MODE) C----------------------------------------------------------------------- C DECAY (INTO 3 PARTICLES) C C TREATES DECAY INTO 3 PARTICLES; FULLY CONSERVING ENERGY AND MOMENTA C KINEMATIC RANGE PARAMETRISATION SEE PHYS. LETT. 204B (1988) 90-91 C FOR LEPTONIC KAON DACAY: THE POLARIZATION OF THE MUON AND C THE NEUTRINO PRODUCTION IS INCLUDED. C THIS SUBROUTINE IS CALLED FROM ETADEC, KDECAY, AND PI0DEC C ARGUMENTS: C AM0 = MASS OF DECAYING PARTICLE C AM3, AM4, AM5 = MASSES OF RESULTING PARTICLES C PARAMA = DALITZ AMPLITUDE PARAMETER (SEE BELOW) C PARAMB = DALITZ AMPLITUDE PARAMETER (SEE BELOW) C PARAMC = DALITZ AMPLITUDE PARAMETER (SEE BELOW) C AMPMX = MAXIMUM AMPLITUDE OF DALITZ PLOT C MODE = 1 FOR DECAY KAON ----> 3 PIONS C = 2 FOR DECAY ETA ----> 3 PIONS OR 2 PIONS + GAMMA C FOR DECAY PI(0) ----> ELECTRON + POSITRON + GAMMA C = 3 FOR DECAY KAON ----> PION + MUON + NEUTRINO C = 4 FOR DECAY KAON ----> PION + ELECTRON + NEUTRINO C C AMPLITUDE PARAMETERS PARAMA, PARAMB, PARAMC ARE DEPENDENT ON MODE: C FOR MODE=1: PARAMA = G DALITZ AMPLITUDE PARAMETRISATION SEE C PARAMB = H PHYS. LETT. 204B (1988) 181 - 193 C PARAMC = K C C FOR MODE=2: PARAMA = A DALITZ AMPLITUDE PARAMETRISATION SEE C PARAMB = DUMMY PHYS. LETT. 204B (1988) 173 - 175; C PARAMC = DUMMY J.G. LAYTER ET.AL. PHYS.REV.D7(1973)2565 C C FOR MODE>2: PARAMA = LAMBDA-PLUS DALITZ AMPLITUDE PARAMETRISATION C PARAMB = LAMBDA-ZERO SEE PHYS. LETT. 204B (1988) C PARAMC = DUMMY 182 - 194 C C DESIGN : D. HECK IK3 FZK KARLSRUHE C----------------------------------------------------------------------- IMPLICIT NONE *KEEP,CONST. COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER *KEEP,DECAY. COMMON /DECAY/ GAM345,COS345,PHI345 DOUBLE PRECISION GAM345(3),COS345(3),PHI345(3) *KEEP,PAM. COMMON /PAM/ PAMA,SIGNUM DOUBLE PRECISION PAMA(6000),SIGNUM(6000) *KEEP,PARPAR. COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C, * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14), * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH INTEGER ITYPE,LEVL *KEEP,PARPAE. DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE), * (CURPAR(4), PHI ), (CURPAR(5), H ), * (CURPAR(6), T ), (CURPAR(7), X ), * (CURPAR(8), Y ), (CURPAR(9), CHI ), * (CURPAR(10),BETA), (CURPAR(11),GCM ), * (CURPAR(12),ECM ) *KEEP,POLAR. COMMON /POLAR/ POLART,POLARF DOUBLE PRECISION POLART,POLARF *KEEP,RANDPA. COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR DOUBLE PRECISION FAC,U1,U2 REAL RD(3000) INTEGER ISEED(103,10),NSEQ LOGICAL KNOR *KEEP,RUNPAR. COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB, * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN, * MONIOU,MDEBUG,NUCNUC, * CETAPE, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE, * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG COMMON /RUNPAC/ DSN,HOST,USER DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB REAL STEPFC INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE INTEGER CETAPE CHARACTER*79 DSN CHARACTER*20 HOST,USER LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG *KEND. DOUBLE PRECISION ABYM,AMPLI,AMPMX,AM0,AM3,AM34I,AM34SQ,AM35SQ, * AM4,AM5,APARAL,APERPN,AUXA,AUXB,AUX1,AUX2,AUX2A, * AUX3,AUX4,AUX4A,AUX5,AUX6, * AUX7,AUX8,AUX10,AUX12,AUX14,BBYM,BOFQ, * CM0SQ,CM3SQ,CM3SQI,CM4SQ,CM5SQ,COSALF,COSBET, * COSFI4,COSOME,COSPHI,COSPSI,COS3CM,COS4CM,COS5CM, * DISCR,EPIPRM,E3CM,E3STAR,E4CM,E5CM,E5STAR,FACT, * GRLAMD,OMEGA,PA,PARAMA,PARAMB,PARAMC,PB,PC,PSI, * P3CM,P3SQ,P4CM,P4SQ,P5CM,P5SQ,ROOT1,ROOT2, * SINALF,SINBET,SINFI4,SINFI5,SINOMG,SINPHI,SINPSI, * SINT4,SINT4I,SINT5I,SIN3CM,S0,TBYMSS,XIT,XI0 INTEGER MODE C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,444) AM0,AM3,AM4,AM5 444 FORMAT(' DECAY6: AM0',1P,E10.3,' AM3',E10.3,' AM4',E10.3, * ' AM5',E10.3) C CALCULATE AUXILIARY QUANTITIES CM0SQ = AM0**2 CM3SQ = AM3**2 CM4SQ = AM4**2 CM5SQ = AM5**2 AUX1 = (AM3 + AM4)**2 AUX2A = (AM0 - AM5)**2 AUX2 = AUX2A - AUX1 AUX3 = (AM3 + AM5)**2 AUX4A = (AM0 - AM4)**2 AUX4 = AUX4A - AUX3 AUX5 = CM3SQ - CM4SQ AUX6 = CM0SQ - CM5SQ AUX7 = 0.5D0 / AM0 IF ( MODE .EQ. 1 ) THEN AUX8 = (AM0 - AM3)**2 S0 = OB3 * ( CM0SQ + CM3SQ + CM4SQ + CM5SQ ) AUX10 = 1.D0 / PAMA(8)**2 ELSEIF ( MODE .EQ. 2 ) THEN AUX14 = 1.D0 / (AM0 - AM3 - AM4 - AM5) ELSEIF ( MODE .EQ. 3 .OR. MODE .EQ. 4 ) THEN CM3SQI = 1.D0 / CM3SQ AUX12 = (CM0SQ + CM3SQ - CM4SQ) * AUX7 C XI0 IS XI(0); GRLAMD IS GREAT LAMBDA XI0 = ( CM0SQ - CM3SQ) * CM3SQI * (PARAMB - PARAMA) GRLAMD = -XI0 * PARAMA ELSE WRITE(MONIOU,*) 'DECAY6: UNEXPECTED MODE =',MODE RETURN ENDIF 100 CALL RMMAR( RD,3,1 ) C ARE INVARIANT MASS SQUARES INSIDE BOUNDARY OF DALITZ PLOT? AM34SQ = AUX2 * RD(1) + AUX1 AM35SQ = AUX4 * RD(2) + AUX3 AM34I = 0.5D0 / SQRT(AM34SQ) E3STAR = (AUX5 + AM34SQ) * AM34I E5STAR = (AUX6 - AM34SQ) * AM34I ROOT1 = SQRT(E3STAR**2 - CM3SQ ) ROOT2 = SQRT(E5STAR**2 - CM5SQ ) DISCR = AM35SQ - (E3STAR + E5STAR)**2 C REJECT RANDOM NUMBERS, IF OUTSIDE KINEMATIC BOUNDARY OF DALITZ PLOT IF ( DISCR .GT. -(ROOT1 - ROOT2)**2 ) GOTO 100 IF ( DISCR .LT. -(ROOT1 + ROOT2)**2 ) GOTO 100 C E3CM, E4CM, E5CM ARE ENERGIES IN THE C. M. SYSTEM E4CM = (CM0SQ + CM4SQ - AM35SQ) * AUX7 E5CM = (CM0SQ + CM5SQ - AM34SQ) * AUX7 E3CM = AM0 - E4CM - E5CM IF ( MODE .EQ. 1 ) THEN FACT = AUX10 * (AUX2A - 2.D0*AM0*(E5CM-AM5) - S0) C AMPLITUDE OF SQUARED MATRIX ELEMENT (SEE PHYS. LETT. B204 (1988) 181) AMPLI = 1.D0 + PARAMA*FACT + PARAMB*FACT**2 + PARAMC*( AUX10 * * ( AUX4A -AUX8 -2.D0*(E4CM-AM4-E3CM+AM3)*AM0 ) )**2 ELSEIF ( MODE .EQ. 2 ) THEN C AMPLITUDE OF SQUARED MATRIX ELEMENT (SEE PHYS. LETT. B204 (1988) 173) C REF: J. G. LAYTER ET AL., PHYS. REV. D7 (1973) 2565 AMPLI = 1.D0 + PARAMA * ( 3.D0 * (E5CM - AM5) * AUX14 - 1.D0 ) ELSE C EPIPRM IS (ENERGY OF PION)PRIMED EPIPRM = AUX12 - E3CM C PA, PB, AND PC ARE THE A, B, AND C PARAMETERS PA = AM0 * ( 2.D0 * E4CM * E5CM - AM0 * EPIPRM ) * + CM4SQ * ( 0.25D0 * EPIPRM - E5CM ) PB = CM4SQ * ( E5CM - 0.5D0 * EPIPRM ) PC = CM4SQ * EPIPRM * 0.25D0 C TBYMSS IS T DIVIDED BY MASS SQUARE OF PION TBYMSS = (CM0SQ + CM3SQ - 2.D0 * AM0 * E3CM) * CM3SQI C XIT IS XI(T) XIT = XI0 + GRLAMD*TBYMSS C AMPLITUDE OF SQUARED MATRIX ELEMENT (PHYS. LETT. B204 (1988) 183) AMPLI = (1.D0 + PARAMA*TBYMSS)**2 * ( PA + XIT*PB + XIT**2 *PC ) ENDIF C REJECT RANDOM NUMBERS, IF RD(3) IS LARGER THAN DALITZ PLOT AMPLITUDE IF ( RD(3)*AMPMX .GT. AMPLI ) GOTO 100 IF (DEBUG) WRITE(MDEBUG,*)'DECAY6: E3CM,E4CM,E5CM=', * SNGL(E3CM),SNGL(E4CM),SNGL(E5CM) C P3CM, P4CM, P5CM ARE MOMENTA IN THE C.M. SYSTEM C P3SQ, P4SQ, P5SQ ARE SQUARED MOMENTA IN C.M. SYSTEM P5SQ = E5CM**2 - CM5SQ P5CM = SQRT(P5SQ) P4SQ = E4CM**2 - CM4SQ P4CM = SQRT(P4SQ) P3SQ = E3CM**2 - CM3SQ P3CM = SQRT(P3SQ) C ANGLE ALFA AND BETA ARE BETWEEN PARTICLE 3 AND 4 RSP. 3 AND 5 COSALF = (P5SQ - P3SQ - P4SQ) / (2.D0 * P3CM * P4CM) SINALF = -SQRT(1.D0 - COSALF**2) COSBET = (P4SQ - P3SQ - P5SQ) / (2.D0 * P3CM * P5CM) SINBET = SQRT(1.D0 - COSBET**2) C NOW SELECT RANDOM NUMBERS FOR THREE INDEPENDENT ANGLES IN CM-SYSTEM C COS3CM AND PHI ARE ANGLES OF PARTICLE 3 RELATIVE TO DECAYING PARTICLE CALL RMMAR( RD,3,1 ) COS3CM = 2.D0*RD(1) - 1.D0 SIN3CM = SQRT(1.D0 - COS3CM**2) PHI345(1) = PI2 * RD(2) COSPHI = COS( PHI345(1) ) SINPHI = SIN( PHI345(1) ) C ANGLE PSI GIVES ROTATION OF PLANE (3,4,5) RELATIVE TO PLANE (1,3) PSI = PI2 * RD(3) COSPSI = COS(PSI) SINPSI = SIN(PSI) C CALCULATE ALL NEEDED POLAR AND AZIMUTHAL ANGLES IN THE CM-SYSTEM COS4CM = COS3CM * COSALF - SIN3CM * COSPSI * SINALF IF ( ABS(COS4CM) .LT. 1.D0 ) THEN SINT4 = SQRT(1.D0 - COS4CM**2) SINT4I = 1.D0 / SINT4 AUXA = COS3CM * COSPSI * SINALF + SIN3CM * COSALF COSFI4 = (COSPHI*AUXA-SINPHI*SINPSI*SINALF) * SINT4I PHI345(2) = ACOS( MAX( -1.D0, MIN( 1.D0, COSFI4 ) ) ) SINFI4 = (SINPHI*AUXA+COSPHI*SINPSI*SINALF) * SINT4I IF ( SINFI4 .LE. 0.D0 ) PHI345(2) = PI2 - PHI345(2) ELSE PHI345(2) = 0.D0 ENDIF C CALCULATE GAMMA FACTORS AND POLAR ANGLES IN LABORATORY SYSTEM GAM345(1) = GAMMA * (E3CM + BETA * P3CM * COS3CM) / AM3 COS345(1) = MIN( 1.D0, (BETA * E3CM + P3CM * COS3CM) * GAMMA * / ( AM3 * SQRT(GAM345(1)**2 - 1.D0) ) ) GAM345(2) = GAMMA * (E4CM + BETA * P4CM * COS4CM) / AM4 COS345(2) = MIN( 1.D0, (BETA * E4CM + P4CM * COS4CM) * GAMMA * / ( AM4 * SQRT(GAM345(2)**2 - 1.D0) ) ) C CALCULATE PARAMETERS OF PARTICLE 5, IF NEEDED IF ( MODE .LE. 2 ) THEN COS5CM = COS3CM * COSBET - SIN3CM * COSPSI * SINBET IF ( ABS(COS5CM) .LT. 1.D0 ) THEN SINT5I = 1.D0 / SQRT(1.D0 - COS5CM**2) AUXB = COS3CM * COSPSI * SINBET + SIN3CM * COSBET PHI345(3) = ACOS((COSPHI*AUXB-SINPHI*SINPSI*SINBET) * SINT5I) SINFI5 = (SINPHI*AUXB+COSPHI*SINPSI*SINBET) * SINT5I IF ( SINFI5 .LE. 0.D0 ) PHI345(3) = PI2 - PHI345(3) ELSE PHI345(3) = 0.D0 ENDIF IF ( AM5 .NE. 0.D0 ) THEN GAM345(3) = GAMMA * (E5CM + BETA * P5CM * COS5CM) / AM5 COS345(3) = MIN( 1.D0, (BETA * E5CM + P5CM * COS5CM) * GAMMA * / ( AM5 * SQRT(GAM345(3)**2 - 1.D0) ) ) ELSE C IF PARTICLE 5 IS GAMMA RAY OR NEUTRINO, THEN GAM345(3) IS THE ENERGY GAM345(3) = GAMMA * (E5CM + BETA * P5CM * COS5CM) COS345(3) = MIN( 1.D0, (BETA * E5CM + P5CM * COS5CM) * GAMMA * / GAM345(3) ) ENDIF ENDIF IF ( MODE .EQ. 3 ) THEN C CALCULATION OF MUON POLARIZATION. WE FOLLOW THE DESCRIPTION OF C L. JAUNEAU, IN: METHODS IN SUBNUCLEAR PHYSICS, VOL. 3, M. NIKOLIC ED. C (GORDON + BREACH, NEW YORK, 1969), P. 123 C SEE ALSO: L.M. CHOUNET ET AL., PHYS. REP. 4 (1972) 199, APPENDIX 1. C SEE ALSO: N. CABBIBO, A. MAKSYMOWICZ, PHYS. LETT. 9 (1964) 352 C (CORRECTIONS IN: PHYS. LETT. 11 (1964) 360; 14 (1965) 72) C WE DEFINE BOFQ (READ: B OF Q), WHICH IS -B(Q**2)*4 BOFQ = 1.D0 - XIT C ABYM AND BBYM (READ A BY M; B BY M) ARE THE QUANTITIES A/M AND B/M ABYM = AM0 * ( BOFQ * EPIPRM - 2.D0 * E5CM ) BBYM = CM0SQ + 0.25D0 * CM4SQ * BOFQ**2 - BOFQ * AM0 * E4CM C NOW CALCULATE THE COMPONENTS APARAL (PARALLEL TO MU DIRECTION) AND C APERPN (PERPENDICULAR TO MU DIRECTION) USING QUANTITIES DEFINED IN C KAON REST SYSTEM. NOTE OUR DEFINITION OF SINALF (ALWAYS WITH NEGATIVE C SIGN) OPPOSITE TO CABBIBO'S SIN(PSI) AND JAUNEAU'S SIN(THETA) APARAL = -P3CM*AM4*BBYM*COSALF - P4CM * ( AM0*ABYM - BBYM * * ( P3CM*SINALF*(E4CM-AM4)/P4CM + AM0 - E3CM ) ) APERPN = P3CM*AM4*BBYM*SINALF C NOW NORMALIZE THE PARALLEL COMPONENT OF POLARIZATION; POLART IS C COSINE OF THE ANGLE BETWEEN MUON MOMENTUM AND POLARISATION POLART = APARAL / SQRT(APARAL**2 + APERPN**2) C THE POLARIZATION VECTOR LIES IN THE PLANE OF MOMENTA (PION,MUON). C OMEGA IS THE ANGLE BY WHICH THE DECAY PLANE (PION,MUON) IS ROTATET C AROUND THE DIRECTION OF MUON RELATIVE TO THE PLANE (KAON,MUON) IF ( ABS(COS4CM) .LT. 1.D0 .AND. SINALF .NE. 0.D0 ) THEN COSOME = (COS4CM*COSALF - COS3CM)*SINT4I/SINALF OMEGA = ACOS( MAX( -1.D0, MIN( 1.D0, COSOME ) ) ) IF ( SINFI4 .NE. 0.D0 ) THEN SINOMG = ( COSFI4 * ( COSALF - COS3CM*COS4CM ) * SINT4I * - SIN3CM * COSPHI ) / (SINALF*SINFI4) IF ( SINOMG .LT. 0.D0 ) OMEGA = PI2 - OMEGA ENDIF ELSE OMEGA = 0.D0 ENDIF POLARF = OMEGA ENDIF RETURN END