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

      SUBROUTINE PIGEN
C
C*********************************************************************
C  DESIGN  : D. HECK   IK3  FZK KARLSRUHE
C  DATE    : JUL  31, 1989
C*********************************************************************
C  THIS SUBROUTINE STEERS THE PHOTONUCLEAR REACTION:
C    FOR PRODUCTION OF 1 PION, PIGEN1 IS CALLED.
C    FOR PRODUCTION OF 2 PIONS, PIGEN2 IS CALLED.
C    FOR PRODUCTION OF MORE PARTICLES, SDPM IS CALLED.
C*********************************************************************
      DOUBLE PRECISION PEIG,REGPAR,REGGEN,REGLVL
      DOUBLE PRECISION ENERN
      DIMENSION REGPAR(12)
*KEEP,BUFFS.
      COMMON /BUFFS/   RUNH,RUNE,EVTH,EVTE,DATAB,LH
      INTEGER          MAXBUF,MAXLEN
      PARAMETER        (MAXBUF=39*7)
      PARAMETER        (MAXLEN=12)
      REAL             RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
     *                 RUNE(MAXBUF),DATAB(MAXBUF)
      INTEGER          LH
      CHARACTER*4      CRUNH,CRUNE,CEVTH,CEVTE
      EQUIVALENCE      (RUNH(1),CRUNH), (RUNE(1),CRUNE)
      EQUIVALENCE      (EVTH(1),CEVTH), (EVTE(1),CEVTE)
*KEEP,GENER.
      COMMON /GENER/   GEN,ALEVEL
      DOUBLE PRECISION GEN,ALEVEL
*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,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,STACKE.
      COMMON/STACKE/   E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
      DOUBLE PRECISION E(60),TIME(60)
      REAL             X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
      INTEGER          IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
*KEND.
      COMMON/ACLOCK/NCLOCK,JCLOCK
C_____IF (NCLOCK.GT.JCLOCK) THEN
C______WRITE(MDEBUG,* )' PIGEN: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
C______CALL AUSGB2
C_____END IF
      IF(DEBUG)WRITE(MDEBUG,*)'PIGEN : E=',E(NP)
C***  INCREASE AGE, WE HAVE HADRONIC INTERACTION
      IGEN(NP)=IGEN(NP)+1
      SECPAR(9)=IGEN(NP)
      SECPAR(10)=-Z(NP)
      PEIG=E(NP)
      IF (U(NP)**2+V(NP)**2.GT.3.E-38) THEN
       ANGLEX = -ATAN2(V(NP),U(NP))
      ELSE
       ANGLEX = 0.
      END IF
C ***  SUBTRACT EM SUBSHOWER FROM NKG CALCULATION
      IF ( FNKG ) THEN
       SECPAR(3) = W(NP)
       SECPAR(4) = ANGLEX
       SECPAR(5) = -Z(NP)
       ENERN = -PEIG*1.D-3
       CALL NKG(ENERN)
      ENDIF
      CALL RMMAR(RNNO90,1,2)
      IF (RNNO90.GT.(PEIG-400.D0)/1000.D0) THEN
C ***  FOR ENERGIES BETWEEN 400 MEV AND 1400 MEV (=1000+400) DECIDE
C ***  BY CHANCE WHETHER ONE OR TWO PIONS ARE GENERATED
C ***  PIGEN1 TREATES THE PRODUCTION OF 1 PION
       CALL PIGEN1
      ELSE IF(RNNO90.GT.(PEIG-2000.D0)/1000.D0) THEN
C ***  FOR ENERGIES BETWEEN 2000MEV AND 3000MEV (=1000+2000) DECIDE
C ***  BY CHANCE WHETHER 2 (PIGEN2) OR MORE PIONS (DPM) ARE GENERATED
C ***  PIGEN2 TREATES THE PRODUCTION OF 2 PIONS
       CALL PIGEN2
      ELSE
C ***  AT HIGHER ENERGIES MORE THAN 2 PIONS ARE GENERATED BY
C ***  THE DUAL PARTON MODEL, BY VENUS, BY SIBYLL, BY QGS, OR BY DPMJET
        DO 191 K=1,MAXLEN
C  ***  SAVE CURPAR PARTICLE INTO REGISTER REGPAR
        REGPAR(K)=CURPAR(K)
191    CONTINUE
192    CONTINUE
       REGGEN = GEN
       REGLVL = ALEVEL
C ***  FILL CURRENT EGS4-PARTICLE INTO CURPAR
       ITYPE=1
       CURPAR(1)=1.D0
       CURPAR(2)=PEIG*1.D-3
       CURPAR(3)=W(NP)
       CURPAR(4)=ANGLEX
       CURPAR(5)=-Z(NP)
       CURPAR(6)=TIME(NP)
       CURPAR(7)=X(NP)
       CURPAR(8)=-Y(NP)
       CURPAR(9)=0.D0
       CURPAR(10)=1.D0
       CURPAR(12)=SQRT(PAMA(14)*(PAMA(14)+PEIG*2.D-3))
       CURPAR(11)=(PEIG*1.D-3+PAMA(14))/CURPAR(12)
       GEN = IGEN(NP)
       ALEVEL = -Z(NP)
C ***  ELIMINATE GAMMA FROM EGS-STACK
       NP=NP-1
C ***  HDPM, VENUS, SIBYLL, QGS, DPMJET GIVE ALL PARTICLES TO SECPAR
       CALL TSTINI
       CALL SDPM
       CALL TSTEND
        DO 201 K=1,MAXLEN
C  ***  RESTORE CURPAR PARTICLE FROM REGPAR
        CURPAR(K)=REGPAR(K)
201    CONTINUE
202    CONTINUE
       GEN = REGGEN
       ALEVEL = REGLVL
C ***  END OF MANY PION GENERATION
      END IF
      RETURN
      END