source: trunk/MagicSoft/Simulation/Corsika/Mmcs/difrac.f@ 10099

      SUBROUTINE DIFRAC( NRETFL )
 
C-----------------------------------------------------------------------
C  (SINGLE) DIF(F)RAC(TION)
C
C  SETS PARAMETERS FOR HDPM IN CASE OF SINGLE DIFFRACTION
C  THIS SUBROUTINE IS CALLED FROM HDPM
C  ARGUMENT:
C   NRETFL = 0  CORRECT ENDING OF SUBROUTINE
C          = 1  INCORRECT ENDING OF SUBROUTINE
C-----------------------------------------------------------------------
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*KEEP,DPMFLG.
      COMMON /DPMFLG/  NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
      INTEGER          NFLAIN,NFLDIF,NFLPI0,NFLCHE,NFLPIF,NFRAGM
*KEEP,INTER.
      COMMON /INTER/   AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
     *                 FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
     *                 RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
     *                 WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN,
     *                 IDIF,ITAR
      DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
     *                 FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
     *                 RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
     *                 WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN
      INTEGER          IDIF,ITAR
*KEEP,LEPAR.
      COMMON /LEPAR/   LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
      INTEGER          LEPAR1,LEPAR2,LASTPI,NRESPC,NRESPN,NCPLUS
*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,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.
 
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,*) 'DIFRAC:'
 
C  DECIDE FIRST, WHETHER PROJECTILE OR TARGET DIFFRACTION
      CALL RMMAR( RD,1,1 )
      IF ( RD(1) .LE. 0.5 ) THEN
C  PROJECTILE DIFFRACTON, TARGET DIFFRACTION FLAG IS NOT SET
        NFTARD = 0
C  MASS OF INCOMING PARTICLE AND PI(0) MASS
C  PI(0) IS MINIMAL OUTCOME OF SECONDARIES IN DIFRAC
        XM0    = ( PAMA(LEPAR1) + PAMA(7) )**2
      ELSE
C  TARGET DIFFRACTON, SET TARGET DIFFRACTION FLAG
        NFTARD = 1
C  MASS OF NUCLEON AND PI(0) MASS
C  PI(0) IS MINIMAL OUTCOME OF SECONDARIES IN DIFRAC
        XM0    = ( PAMA(LEPAR2) + PAMA(7) )**2
      ENDIF
C  MAXIMAL DIFFRACTIVE MASS, FACTOR 0.15 GIVEN BY COHERENCE CONDITION
      XMX    = 0.15D0 * S
 
C  THROW MAXIMAL 200 TIMES TO GET A GOOD DIFFRACTIVE MASS
      NCDIFL = 0
  7   CONTINUE
C  GET DIFFRACTIVE MASS
      CALL RMMAR( RD,2,1 )
C  GET S (=ECM**2) (WHY THIS WAY OF THROWING ???)
      SDIF = (XMX/XM0)**RD(1) * XM0
 
      IF ( SDIF .LE. XM0 ) THEN
        IF ( NCDIFL .LE. 200 ) THEN
          NCDIFL = NCDIFL + 1
          GOTO 7
        ELSE
C  SET RETURN FLAG TO ERROR
          NRETFL = 1
          RETURN
        ENDIF
      ENDIF
 
C  DISTRIBUTION OF DIFFRACTIVE MASS FLATTENS OFF FOR DIFFRACTIVE
C  MASS SQUARED .LE. 2 GEV
      IF ( SDIF .LE. 2.D0 ) THEN
C----- SO GEHT DAS NICHT!!   16.12.91 D.H.
        SDIF = RD(2) * (2.D0 - XM0) + XM0
      ENDIF
C  SQRT(S) IS ECM
      ECMDIF = SQRT(SDIF)
C  LOG(S), LOG(S)**2
      DLOG   = LOG(SDIF)
      DLOGSQ = DLOG**2
      IF ( DEBUG ) WRITE(MDEBUG,*) 'DIFRAC: SDIF,ECMDIF,NFTARD=',
     *               SNGL(SDIF),SNGL(ECMDIF),NFTARD
 
C  RAPIDITY IN CMS OF DIFFRACTIVE SYSTEM
C  TO CALCULATE DMLOG, SUBTRACT SUM OF MASS SQUARES FROM SDIF
C  PI(0) MASS SQUARED IS 0.0182.
      IF ( NFTARD .EQ. 0 ) THEN
        YY0   =  LOG(ECMDPM/ECMDIF)
        DMLOG =  LOG(SDIF - 0.0182D0 - PAMA(LEPAR1)**2)
      ELSE
        YY0   = -LOG(ECMDPM/ECMDIF)
        DMLOG =  LOG(SDIF - 0.0182D0 - PAMA(LEPAR2)**2)
      ENDIF
      IF ( DEBUG ) WRITE(MDEBUG,*) 'DIFRAC: YY0,DMLOG=',
     *                           SNGL(YY0),SNGL(DMLOG)
C  CENTRAL RAPIDITY DENSITY IN CMS OF DIFFRACTIVE SYSTEM
C  PARAMETRISATION SEE CAPDEVIELLE,J.PHYS.G:NUCL.PHYS.16(1990)1539 EQ.7
C  WE USE ONLY THE LOW-ENERGY PART OF THE PARAMETRISATION, AS SDIF DOES
C  NOT REACH THE HIGHER VALUES
      DC0 = 0.82D0 * (SDIF**0.107D0)
 
C  THERE ARE 3 ENERGY DEPENDENT FORMULAS FOR AVERAGE CHARGED
C  MULTIPLICITY ( AVCH1 );
C  PARAMETRISATIONS SEE CAPDEVIELLE,J.PHYS.G:NUCL.PHYS.16(1990)1539 EQ.8
      IF     ( ECMDIF .LE. 187.5D0 ) THEN
C  CHARGED MULTIPLICITY (M**2 IN PLACE OF S)
        AVCH1 = 0.57D0 + 0.584D0*DLOG + 0.127D0*DLOGSQ
      ELSEIF ( ECMDIF .LE. 945.5D0 ) THEN
        AVCH1 = -6.55D0 + 6.89D0 * SDIF**0.131D0
      ELSE
        AVCH1 = 3.4D0 * SDIF**0.17D0
      ENDIF
C  PARAMETRISATION IS BASED ON COLLIDER DATA WHERE PROTON AND ANTIPROTON
C  ARE INCLUDED. LOWER LIMIT FOR AVERAGE CHARGED MULTIPLICITY IS 1.
      AVCH1 = MAX( 1.D0, AVCH1 )
 
C  CENTER OF GAUSSIAN 1ST+2ND STRING OF FRAGMENTATION SYSTEM
      POSC2 = 0.146D0 * DMLOG + 0.072D0
C  WIDTH  OF GAUSSIAN 1ST+2ND STRING OF FRAGMENTATION SYSTEM
      WIDC2 = 0.120D0 * DMLOG + 0.180D0
C  INTERACTION FACTOR GNU FOR INTERACTION WITH NUCLEUS;
      IF ( NFLAIN .EQ. 0 ) THEN
        GNU   = 1.D0
        AVCH3 = 0.D0
        POSC3 = 0.D0
        WIDC3 = 1.D0
      ELSE
C  NEW PARAMETRIZATION OF J.N.CAPDEVIELLE (MARCH 93)
        GNU   = (0.4826D0 + 3.522D-2 * DLOG) * TAR**0.31D0
C  CENTER OF GAUSSIAN FOR 3RD STRING (FROM TARGET)
        POSC3 = +3.D0 - 2.575D0 * EXP( -0.081756452D0 * GNU )
C  WIDTH  OF GAUSSIAN FOR 3RD STRING (FROM TARGET)
        WIDC3 = 1.2338466D0 + 0.078969916D0 * LOG(GNU)
        IF ( ECMDIF .LE. 137.D0 ) THEN
          AVCH3 = 0.57D0 * AVCH1 * (GNU-1.D0)
        ELSE
          AVCH3 = 0.5D0  * AVCH1 * (GNU-1.D0)
        ENDIF
      ENDIF
      IF ( DEBUG ) WRITE(MDEBUG,100)
     *          SNGL(POSC2),SNGL(WIDC2),SNGL(POSC3),SNGL(WIDC3)
 100  FORMAT(' DIFRAC: POSC2,WIDC2,POSC3,WIDC3=',4F12.7)
C  AVERAGE CHARGED, INCLUDING THOSE FROM TARGET
      AVCH = AVCH1 + AVCH3
C  THE FOLOWING PROCEDURE IS TO PRODUCE PHOTONS FROM UNKNOWN NEUTRAL
C  DECAYS FOLLOWING CORRELATION WITH CHARGED PARTICLES BASED ON
C  PHOTON EXCESS AT COLLIDER EXPERIMENTS. SEUGP IS <N_PHOTON>
C  PROBLEM OF THE RISE OF THE UNKNOWN ETA PRODUCTION CROSS SECTION
C  IS SOLVED WITH PARAMETRISATION OF UA5 (Z. PHYS. C43 (1989) 75)
      IF ( ECMDIF .LE. 103.D0 ) THEN
        SEUGP = -1.27D0 + 0.52D0 * DLOG + 0.148D0 * DLOGSQ
      ELSE
C  AT HIGH DIFFRACTIVE MASS USE PARAMETRISATION OF THOUW ????
        SEUGP = -18.7D0 + 11.55D0 * SDIF**0.1195D0
      ENDIF
      SEUGP = MAX( 0.5D0, SEUGP )
      IF ( DEBUG ) WRITE(MDEBUG,110)
     *   SNGL(DC0),SNGL(AVCH1),SNGL(AVCH3),SNGL(AVCH),SNGL(SEUGP)
 110  FORMAT(' DIFRAC: DC0,AVCH1,AVCH3,AVCH,SEUGP=',5F12.6)
 
C  SET RETURN FLAG TO OK
      NRETFL = 0
      RETURN
      END