SUBROUTINE MUBREM C----------------------------------------------------------------------- C MU(ON) BREM(SSTRAHLUNG) C C TREATES MUON BREMSSTRAHLUNG (WITHOUT POLARISATION EFFECTS) C IN ANALOGY WITH SUBROUTINE GBREMM FROM GEANT WRITTEN BY L. URBAN C EXPLANATIONS SEE CERN PROGRM LIBRARY LONG WRITEUP W5013 C THIS SUBROUTINE IS CALLED FROM MUTRAC 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,ELABCT. COMMON /ELABCT/ ELCUT DOUBLE PRECISION ELCUT(4) *KEEP,GENER. COMMON /GENER/ GEN,ALEVEL DOUBLE PRECISION GEN,ALEVEL *KEEP,MUPART. COMMON /MUPART/ AMUPAR,BCUT,CMUON,FMUBRM,FMUORG DOUBLE PRECISION AMUPAR(14),BCUT,CMUON(11) LOGICAL FMUBRM,FMUORG *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,REST. COMMON /REST/ CONTNE,TAR,LT DOUBLE PRECISION CONTNE(3),TAR INTEGER LT *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 ALFA1,AUXIL,BETA1,COSTH3,COSTH4,CREJ,D,F1, * EKIN,EMUON,PHI3,SCREJ,SINTH3,THETA3,U,UMAX, * V,VC,VM,V1,W1,Z INTEGER I DATA ALFA1/0.625D0/ C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9) 444 FORMAT(' MUBREM: CURPAR=',1P,9E10.3) C COPY COORDINATES TO SECPAR DO 11 I = 5,8 SECPAR(I) = CURPAR(I) 11 CONTINUE SECPAR( 9) = GEN C TOTAL AND KINETIC ENERGY OF MUON EMUON = PAMA(5) * GAMMA EKIN = EMUON - PAMA(5) IF ( EKIN .LE. BCUT ) THEN C MUON ENERGY IS TOO LOW TO PRODUCE BREMSSTRAHLUNG SECPAR(2) = CURPAR(2) GOTO 900 ENDIF VC = BCUT/EMUON VM = 1.D0 - CMUON(6+LT)/EMUON IF ( VM .LE. 0.D0 ) THEN C MAXIMUM OF BREMSSTRAHLUNG SPECTRUM IS NEGATIVE, NO BREMSSTRAHLUNG SECPAR(2) = CURPAR(2) GOTO 900 ENDIF CREJ = CMUON(3+LT)/EMUON 50 CALL RMMAR(RD,2,1) V = VC*(VM/VC)**RD(1) V1 = 1.D0 - V C COMPUTE REJECTION FUNCTION F1 = CMUON(LT) - LOG(1.D0 + CREJ*V/V1) SCREJ = (V1 + 0.75D0*V*V)*F1/CMUON(LT) IF ( RD(2) .GT. SCREJ ) GOTO 50 C PHOTON ENERGY SECPAR(2) = EMUON * V C RADIATED GAMMA BELOW CUT? THEN REDUCE ENERGY OF MUON IF ( SECPAR(2) .LE. ELCUT(4) ) THEN GO TO 800 ENDIF C SET MATERIAL CONSTANTS CMUON(.) ACCORDING TO C TARGET INDEX LT (1=N, 2=O, 3=AR) WHICH HAS BEEN SET IN BOX2 IF ( LT .EQ. 1 ) THEN Z = 7.D0 ELSEIF ( LT .EQ. 2 ) THEN Z = 8.D0 ELSE Z = 18.D0 ENDIF C GENERATE EMITTED PHOTON ANGLES WITH RESPECT TO MUON DIRECTION C PHI IS GENERATED ISOTROPICALLY AND THETA IS ASSIGNED A UNIVERSAL C ANGULAR DISTRIBUTION WITH D=D(Z,E,V) C THIS FUNCTION APPROXIMATES THE REAL DISTRIBUTION FUNCTION WHICH CAN C BE FOUND IN: YUNG-SU TSAI, REV. MOD. PHYS. 46(1974)815 C +ERRATUM: REV. MOD. PHYS. 49(1977)421 D = 0.13D0 *(0.8D0 + 1.3D0/Z) * (100.D0 + 1.D0/EMUON) * (1.D0 + V) W1 = 9.D0 / (9.D0 + D) UMAX = EMUON * PI / PAMA(5) 10 CALL RMMAR(RD,3,1) IF ( RD(1) .LE. W1 ) THEN BETA1 = ALFA1 ELSE BETA1 = 3.D0 * ALFA1 ENDIF U = -( LOG(RD(2) * RD(3)) ) / BETA1 C CUT: THETA SHOULD BE .LE. PI ! C THIS CONDITION DEPENDS ON E IN THE CASE OF D=CONST TOO! IF ( U .GE. UMAX ) GOTO 10 THETA3 = U * PAMA(ITYPE) / EMUON COSTH3 = COS( THETA3 ) SINTH3 = SIN( THETA3 ) CALL RMMAR(RD,1,1) PHI3 = PI2 * RD(1) CALL ADDANG( COSTHE,PHI, COSTH3,PHI3, SECPAR(3),SECPAR(4)) IF ( SECPAR(3) .GT. C(29) ) THEN C WRITE BREMSSTRAHLUNG PHOTON TO STACK SECPAR( 1) = 1.D0 SECPAR(10) = H CALL TSTACK ENDIF C REDUCE ENERGY OF MUON 800 CONTINUE EMUON = EMUON * V1 SECPAR(2) = EMUON/PAMA(5) 900 CONTINUE C WRITE MUON TO STACK SECPAR( 1) = CURPAR(1) SECPAR( 3) = CURPAR(3) SECPAR( 4) = CURPAR(4) SECPAR(10) = ALEVEL CALL TSTACK RETURN END