source: trunk/MagicSoft/Simulation/Corsika/Mmcs/hmeson.f@ 1118

      SUBROUTINE HMESON( E,AMASS,ASMASS )
 
C-----------------------------------------------------------------------
C  H(EAVY) MESON
C
C  HANDLES PION INITIATED HEAVY MESON AND ITS DECAY IN UP TO 3 PIONS
C  HEAVY MESON EMITTED FORWARD
C  THIS SUBROUTINE IS CALLED FROM BOX67 AND BOX69
C  ARGUMENTS:
C   E      = AVAILABLE ENERGY IN CM
C   AMASS  = MASS OF HEAVY MESON
C   ASMASS = MASS TO BE LEFT OVER FOR OTHER PARTICLES
C-----------------------------------------------------------------------
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*KEEP,BAL.
      COMMON /BAL/     EBAL
      DOUBLE PRECISION EBAL(10)
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*KEEP,ELASTY.
      COMMON /ELASTY/  ELAST,IELIS,IELHM,IELNU,IELPI
      DOUBLE PRECISION ELAST
      INTEGER          IELIS(20),IELHM(20),IELNU(20),IELPI(20)
*KEEP,MULT.
      COMMON /MULT/    EKINL,MSMM,MULTMA,MULTOT
      DOUBLE PRECISION EKINL
      INTEGER          MSMM,MULTMA(37,13),MULTOT(37,13)
*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,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
*KEEP,VKIN.
      COMMON /VKIN/    BETACM
      DOUBLE PRECISION BETACM
*KEND.
 
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,*) 'HMESON: E,AMASS,ASMASS=',
     *         SNGL(E),SNGL(AMASS),SNGL(ASMASS)
 
      IPI   = 0
      EDHM  = 0.D0
      PACC  = 0.D0
      W     = 0.6D0
 
C  GAMMA AND BETA OF HEAVY MESON IN CM AND LAB
C  E > AMASS + ASMASS  TO KEEP GHMCM > 1.
      GHMCM  = ( E**2+AMASS**2-ASMASS**2 ) / ( 2.D0*E*AMASS )
      BHMCM  = SQRT(GHMCM**2 - 1.D0) / GHMCM
      GHMLAB = GCM * GHMCM * (1.D0 + BETACM * BHMCM)
      BHMLAB = SQRT(GHMLAB**2 - 1.D0) / GHMLAB
 
C  DECAY OF HEAVY MESON
    7 CONTINUE
      IPI  = IPI + 1
C  CHOSE TRANSVERSE MOMENTUM RANDOMLY
      PTPI = PTRANS(DUMMY)
C  CHOSE LONGITUDINAL MOMENTUM RANDOMLY
      IF ( IPI .LT. 3 ) THEN
        P  = PCL(C(40),W)
      ELSE
        P2 = RESTE**2 - PAMA(8)**2 - PTPI**2
        P  = SQRT(MAX( P2, 0.D0 ))
      ENDIF
 
      PTPI  = PTRANS(DUMMY)
      GPIHM = SQRT( P**2 / PAMA(8)**2 + 1.D0 )
      BPIHM = SQRT( GPIHM**2-1.D0 ) / GPIHM
      EDHM  = EDHM + SQRT( PAMA(8)**2 + P**2 + PTPI**2 )
      RESTE = AMASS - EDHM
 
C  FOR FIRST 2 PARTICLES CHOSE RANDOMLY WHETHER FORWARD OR BACKWARD
C  FOR 3. PARTICLE DECIDE ACCORDING TO ACCULMULATED P
      CALL RMMAR( RD,3,1 )
      IF ( IPI .EQ. 3 ) THEN
        IF ( PACC .LE. 0.D0 ) THEN
          RD(1) = 0.
        ELSE
          RD(1) = 1.
        ENDIF
      ENDIF
 
      IF ( RD(1) .GT. 0.5 ) THEN
C  BACKWARD PION
        GPILAB = GHMLAB*GPIHM*(1.D0-BHMLAB*BPIHM)
        PACC   = PACC - P
      ELSE
C  FORWARD PION
        GPILAB = GHMLAB*GPIHM*(1.D0+BHMLAB*BPIHM)
        PACC   = PACC + P
      ENDIF
C  CORRECTIVE ACTION IF GPILAB LESS OR EQUAL TO 1.0
      GPILAB   = MAX( GPILAB, 1.D0 )
 
C  GET NEW DIRECTION
      PLLAB2 = PAMA(8)**2 *(GPILAB**2 - 1.D0)
      CTHETA = SQRT( PLLAB2 / (PTPI**2+PLLAB2) )
      IF ( CTHETA .GE. C(27) ) THEN
        CALL ADDANG( COSTHE,PHI, CTHETA,RD(2)*PI2, SECPAR(3),SECPAR(4) )
        IF ( SECPAR(3) .GE. C(29) ) THEN
          SECPAR(2) = GPILAB
 
C  RANDOM CHARGE ASSIGNMENT FOR PIONS
          IF     ( RD(3) .LE. OB3 ) THEN
            SECPAR(1) = 7.D0
          ELSEIF ( RD(3) .LE. TB3 ) THEN
            SECPAR(1) = 8.D0
          ELSE
            SECPAR(1) = 9.D0
          ENDIF
 
          DO  4  J = 5,8
            SECPAR(J) = CURPAR(J)
    4     CONTINUE
 
          CALL TSTACK
        ENDIF
      ENDIF
      IF ( IPI .LT. 3  .AND.  RESTE .GT. PAMA(8) ) GOTO 7
 
C  STATISTICS ON ENERGY BALANCE, MULTIPLICITY AND ELASTICITY
      EBAL(4) = EBAL(4) + RESTE
      MSMM    = MSMM + IPI
 
C  INELASTICITY STATISTICS
      IN = 1.D0 + SECPAR(2) / GAMMA * 20.D0
      IN = MIN( IN, 20 )
      IELHM(IN) = IELHM(IN) + 1
 
      RETURN
      END