source: trunk/MagicSoft/Simulation/Corsika/Mmcs/vapor.f

      SUBROUTINE VAPOR(MAPROJ,INEW,JFIN,ITYP,PFRX,PFRY)
 
C-----------------------------------------------------------------------
C  (E)VAPOR(ATION OF NUCLEONS AND ALPHA PARTICLES FROM FRAGMENT)
C
C  TREATES THE REMAINING UNFRAGMENTED NUCLEUS
C  EVAPORATION FOLLOWING CAMPI APPROXIMATION
C  SEE: X. CAMPI AND J. HUEFNER, PHYS.REV. C24 (1981) 2199
C  AND  J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990)
C  THIS SUBROUTINE IS CALLED FROM SDPM AND VSTORE
C
C  ARGUMENTS INPUT:
C   MAPROJ       = NUMBER OF NUCLEONS OF PROJECTILE
C   INEW         = PARTICLE TYPE OF SPECTATOR FRAGMENT
C  ARGUMENTS OUTPUT:
C   JFIN         = NUMBER OF FRAGMENTS
C   ITYP(1:JFIN) = NATURE (PARTICLE CODE) OF FRAGMENTS (GEANT)
C   PFRX(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN X-DIRECTION
C   PFRY(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN Y-DIRECTION
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,DPMFLG.
      COMMON /DPMFLG/  NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
      INTEGER          NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
*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 PFR(60),PFRX(60),PFRY(60)
      DOUBLE PRECISION AFIN,AGLH,APRF,BGLH,EEX,PHIFR,RANNOR,SPFRX,SPFRY
      INTEGER          ITYP(60),IARM,INEW,ITYPRM,INRM,IS,IZRM,JC,JFIN,
     *                 K,L,LS,MAPROJ,MF,NFIN,NINTA,NNUC,NPRF,NSTEP
      EXTERNAL         RANNOR
C-----------------------------------------------------------------------
 
      IF (DEBUG) WRITE(MDEBUG,*) 'VAPOR : MAPROJ,INEW = ',MAPROJ,INEW
 
      ITYPRM = INEW
      NPRF   = INEW/100
      NINTA  = MAPROJ - NPRF
      IF ( NINTA .EQ. 0 ) THEN
C  NO NUCLEON HAS INTERACTED
        JFIN    = 1
        PFR(1)  = 0.D0
        ITYP(1) = INEW
        IF (DEBUG) WRITE(MDEBUG,*) 'VAPOR : JFIN,NINTA= ',JFIN,NINTA
        RETURN
      ENDIF
 
C  EXCITATION ENERGY EEX OF PREFRAGMENT
C  SEE: J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990), CHPT. 4.2
      EEX = 0.D0
      CALL RMMAR(RD,2*NINTA,1)
      DO  22  L = 1,NINTA
        IF ( RD(NINTA+L) .LT. RD(L) ) RD(L) = 1. - RD(L)
        EEX = EEX + RD(L)
 22   CONTINUE
C  DEPTH OF WOODS-SAXON POTENTIAL TO FERMI SURFACE IS 0.040 GEV
      IF (DEBUG) WRITE(MDEBUG,*)'VAPOR : EEX = ',SNGL(EEX*0.04D0),' GEV'
C  EVAPORATION: EACH EVAPORATION STEP NEEDS ABOUT 0.020 GEV, THEREFORE
C  NSTEP IS EEX * 0.04/0.02 = EEX * 2.
      NSTEP = INT(EEX*2.D0)
 
      IF ( NSTEP .LE. 0 ) THEN
C  EXCITATION ENERGY TOO SMALL, NO EVAPORATION
        JFIN = 1
        PFR(1)  = 0.D0
        ITYP(1) = INEW
        IF (DEBUG) WRITE(MDEBUG,*) 'VAPOR : JFIN,EEX = ',JFIN,SNGL(EEX)
        RETURN
      ENDIF
 
C  AFIN IS ATOMIC NUMBER OF FINAL NUCLEUS
      APRF = FLOAT(NPRF)
      AFIN = APRF - 1.6D0 * FLOAT(NSTEP)
      NFIN = MAX( INT(AFIN+0.5D0), 0 )
C  CORRESPONDS TO DEFINITION; FRAGMENTATION-EVAPORATION
C  CONVOLUTION EMU07 /MODEL ABRASION EVAPORATION (JNC FZK APRIL 94)
C  NNUC IS NUMBER OF EVAPORATING NUCLEONS
      NNUC = NPRF - NFIN
      IF (DEBUG) WRITE(MDEBUG,*) 'VAPOR : NFIN,NNUC = ',NFIN,NNUC
      JC = 0
 
      IF     ( NNUC .LE. 0 ) THEN
C  NO EVAPORATION
        JFIN = 1
        PFR(1)  = 0.D0
        ITYP(1) = INEW
        RETURN
 
      ELSEIF ( NNUC .GE. 4 ) THEN
C  EVAPORATION WITH FORMATION OF ALPHA PARTICLES POSSIBLE
C  IARM, IZRM, INRM ARE NUMBER OF NUCLEONS, PROTONS, NEUTRONS OF
C  REMAINDER
        DO  31  LS = 1,NSTEP
          IARM = ITYPRM/100
          IF ( IARM .LE. 0 ) GOTO 100
          IZRM = MOD(ITYPRM,100)
          INRM = IARM - IZRM
          JC = JC + 1
          CALL RMMAR(RD,2,1)
          IF ( RD(1).LT.0.2 .AND. IZRM.GE.2 .AND. INRM.GE.2 ) THEN
            ITYP(JC) = 402
            NNUC   = NNUC - 4
            ITYPRM = ITYPRM - 402
          ELSE
            IF ( RD(2)*(IZRM+INRM) .LT. IZRM ) THEN
              ITYP(JC) = 14
              ITYPRM   = ITYPRM - 101
            ELSE
              ITYP(JC) = 13
              ITYPRM   = ITYPRM - 100
            ENDIF
            NNUC = NNUC - 1
          ENDIF
          IF ( NNUC .LE. 0 ) GOTO 50
 31     CONTINUE
      ENDIF
 
      IF ( NNUC .LT. 4 ) THEN
C  EVAPORATION WITHOUT FORMATION OF ALPHA PARTICLES
        CALL RMMAR(RD,NNUC,1)
        DO  32  IS = 1,NNUC
          IARM = ITYPRM/100
          IF ( IARM .LE. 0 ) GOTO 100
          IZRM = MOD(ITYPRM,100)
          JC = JC + 1
          IF ( RD(IS)*IARM .LT. IZRM ) THEN
            ITYP(JC) = 14
            ITYPRM   = ITYPRM - 101
          ELSE
            ITYP(JC) = 13
            ITYPRM   = ITYPRM - 100
          ENDIF
 32     CONTINUE
      ENDIF
 
 50   CONTINUE
      JC = JC + 1
      IF     ( ITYPRM .GT. 101 ) THEN
        ITYP(JC) = ITYPRM
      ELSEIF ( ITYPRM .EQ. 101 ) THEN
        ITYP(JC) = 14
      ELSEIF ( ITYPRM .EQ. 100 ) THEN
        ITYP(JC) = 13
      ELSE
        JC = JC - 1
        IF ( ITYPRM .NE. 0 ) WRITE(MONIOU,*)
     *                  'VAPOR : ILLEGAL PARTICLE ITYPRM =',ITYPRM
      ENDIF
 
 100  JFIN = JC
      IF (DEBUG) WRITE(MDEBUG,*) 'VAPOR : NO   ITYP           PFR'
      IF     ( NFRAGM .EQ. 2 ) THEN
C  EVAPORATION WITH PT AFTER PARAMETRIZED JACEE DATA
        DO  150  MF = 1,JFIN
          PFR(MF) = RANNOR(0.088D0,0.044D0)
          IF (DEBUG) WRITE(MDEBUG,*) MF,ITYP(MF),SNGL(PFR(MF))
 150    CONTINUE
      ELSEIF ( NFRAGM .EQ. 3 ) THEN
C  EVAPORATION WITH PT AFTER GOLDHABER'S MODEL (PHYS.LETT.53B(1974)306)
        DO  160  MF = 1,JFIN
          K = MAX( 1, ITYP(MF)/100 )
          BGLH  = K * (MAPROJ - K) / FLOAT(MAPROJ-1)
C  THE VALUE 0.103 [GEV] IS SIGMA(0)=P(FERMI)/SQRT(5.)
*         AGLH = 0.103D0 * SQRT( BGLH )
C  THE VALUE 0.090 [GEV] IS EXPERIMENTALLY DETERMINED SIGMA(0)
          AGLH = 0.090D0 * SQRT( BGLH )
          PFR(MF) = RANNOR(0.D0,AGLH)
          IF (DEBUG) WRITE(MDEBUG,*) MF,ITYP(MF),SNGL(PFR(MF))
 160    CONTINUE
      ELSE
C  EVAPORATION WITHOUT TRANSVERSE MOMENTUM
        DO  165  MF = 1,JFIN
          PFR(MF) = 0.D0
          IF (DEBUG) WRITE(MDEBUG,*) MF,ITYP(MF),SNGL(PFR(MF))
 165    CONTINUE
      ENDIF
C  CALCULATE RESIDUAL TRANSVERSE MOMENTUM
      SPFRX = 0.D0
      SPFRY = 0.D0
      CALL RMMAR(RD,JFIN,1)
      DO  170  MF = 1,JFIN
        PHIFR = PI * RD(MF)
        PFRX(MF) = PFR(MF) * COS(PHIFR)
        PFRY(MF) = PFR(MF) * SIN(PHIFR)
        SPFRY = SPFRY + PFRY(MF)
        SPFRX = SPFRX + PFRX(MF)
 170  CONTINUE
C  CORRECT ALL TRANSVERSE MOMENTA FOR MOMENTUM CONSERVATION
      SPFRX = SPFRX / JFIN
      SPFRY = SPFRY / JFIN
      DO  180  MF = 1,JFIN
        PFRX(MF) = PFRX(MF) - SPFRX
        PFRY(MF) = PFRY(MF) - SPFRY
 180  CONTINUE
 
      IF (DEBUG) WRITE(MDEBUG,*) 'VAPOR : NINTA,JFIN= ',NINTA,JFIN
      RETURN
      END