DPMJET/user3.0-4.f

   1 *
   2 *===program dpmjet=====================================================*
   3 *
   4       PROGRAM DPMJET
   5
   6       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   7       SAVE
   8
   9       EXTERNAL BDANDI, BDEMFL, BDEVAP, BDMULS, BDNOPT, BDPHFL, BDPWXS,
  10      &         BDRNDM, BDSQZI, BLKDT1, BLKDT2, BLKDT3, BLKDT4, BLKDT5,
  11      &         BLKDT6, BLKDT7, COUNTR, ELPROI, HADINI, BDPREE
  12
  13 C     EXTERNAL PYDATA
  14
  15 * event flag
  16       COMMON /DTEVNO/ NEVENT,ICASCA
  17
  18 *-----------------------------------------------------------------------
  19 * initialization
  20
  21 *   the following statement provides a call to DT_USRHIS(MODE=1) for
  22 *   histogram initialization etc.
  23       CALL DT_DTUINI(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IEMU)
  24
  25 *-----------------------------------------------------------------------
  26 * generation of events
  27
  28       DO 1 IEVT=1,NEVTS
  29
  30 *   some defaults, do not change!
  31          NEVENT = IEVT
  32          KKMAT  = -1
  33          ELAB   = EPN
  34 *   uncomment if dpmjet3 is linked to particle transport code
  35 C        ICASCA = 1
  36
  37 ************************************************************************
  38 * The following lines show how to select the target nucleus for runs
  39 * with composite targets (and fixed projectile and energy!).
  40 *
  41 *   Sampling of the target nucleus (mass number NTMASS, charge NTCHAR)
  42 *   according to the fractions defined with EMULSION input-cards.
  43 *   The different nuclei are numbered as KKMAT = 1,2,3,...  according to
  44 *   their appearance in the input-file.
  45          IF (IEMU.GT.0) THEN
  46 *   Replace this selection by your own one if needed.
  47             CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
  48 *   Kkmat has to be negative for composite targets!
  49             KKMAT = -KKMAT
  50          ENDIF
  51 ************************************************************************
  52
  53 ************************************************************************
  54 * The following lines show how to define projectile, target and energy
  55 * for this event in runs with Glauber-data file pre-initialized for a
  56 * certain range of projectiles, targets and energies. The definitions
  57 * have to be within the pre-initialized parameter range.
  58 *
  59 *   projectile-id (for hadron projectiles)
  60 C        IDP    = 1
  61 *   projectile mass and charge numbers
  62 C        NPMASS = 12
  63 C        NPCHAR = 6
  64 *   target mass and charge numbers
  65 C        NTMASS = 16
  66 C        NTCHAR = 8
  67 *   lab energy
  68 C        ELAB = 200.0D0
  69 ************************************************************************
  70
  71 ************************************************************************
  72 * If an energy-range has been defined with the ENERGY input-card the
  73 * laboratory energy ELAB can be set to any value within that range. For
  74 * example:
  75 C        ELO  = 10.0D0
  76 C        EHI  = 1000.0D0
  77 C        ELAB = DT_RNDM(ELAB)*(EHI-ELO)+ELO
  78 ************************************************************************
  79
  80 *   sampling of one event
  81          CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,ELAB,KKMAT,IREJ)
  82          IF (IREJ.NE.0) GOTO 1
  83
  84 *   the following statement provides a call to DT_USRHIS(MODE=2) from
  85 *   where the final state particles can be obtained
  86
  87          CALL PHO_PHIST(2000,DUM)
  88
  89     1 CONTINUE
  90
  91 *-----------------------------------------------------------------------
  92 * output, statistics etc.
  93
  94 *   the following statement provides a call to DT_USRHIS(MODE=3) in
  95 *   order to calculate histograms etc.
  96       CALL DT_DTUOUT
  97
  98       END
  99 *
 100 *===usrhis=============================================================*
 101 *
 102 CDECK  ID>, DT_USRHIS
 103       SUBROUTINE DT_USRHIS(MODE)
 104
 105       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 106       SAVE
 107 *
 108 * COMMON /DTEVT1/ :
 109 *                   NHKK         number of entries in common block
 110 *                   NEVHKK       number of the event
 111 *                   ISTHKK(i)    status code for entry i
 112 *                   IDHKK(i)     identifier for the entry
 113 *                                (for particles: identifier according
 114 *                                 to the PDG numbering scheme)
 115 *                   JMOHKK(1,i)  pointer to the entry of the first mother
 116 *                                of entry i
 117 *                   JMOHKK(2,i)  pointer to the entry of the second mother
 118 *                                of entry i
 119 *                   JDAHKK(1,i)  pointer to the entry of the first daughter
 120 *                                of entry i
 121 *                   JDAHKK(2,i)  pointer to the entry of the second daughter
 122 *                                of entry i
 123 *                   PHKK(1..3,i) 3-momentum
 124 *                   PHKK(4,i)    energy
 125 *                   PHKK(5,i)    mass
 126 *
 127 * event history
 128
 129       PARAMETER (NMXHKK=200000)
 130
 131       COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
 132      &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
 133      &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
 134 * extended event history
 135       COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
 136      &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
 137      &                IHIST(2,NMXHKK)
 138
 139       GOTO (1,2,3) MODE
 140
 141 *------------------------------------------------------------------
 142 *
 143     1 CONTINUE
 144 *
 145 * initializations
 146 *
 147 *  Called with MODE=1 once at the beginning of the run.
 148 *
 149       RETURN
 150 *
 151 *------------------------------------------------------------------
 152 *
 153     2 CONTINUE
 154 *
 155 * scoring of the present event
 156 *
 157 *  Called with MODE=2 every time one event has been finished.
 158 *
 159 *  The final state particles from the actual event (number NEVHKK)
 160 *  can be found in DTEVT1 and identified by their status:
 161 *
 162 *     ISTHKK(i) = 1    final state particle produced in
 163 *                      photon-/hadron-/nucleon-nucleon collisions or
 164 *                      in intranuclear cascade processes
 165 *                -1    nucleons, deuterons, H-3, He-3, He-4 evaporated
 166 *                      from excited nucleus and
 167 *                      photons produced in nuclear deexcitation processes
 168 *                1001  residual nucleus (ground state)
 169 *
 170 *  The types of these particles/nuclei are given in IDHKK as follows
 171 *
 172 *     all final state part. except nuclei :
 173 *       IDHKK(i)=particle identifier according to PDG numbering scheme
 174 *     nuclei (evaporation products, and residual nucleus) :
 175 *       IDHKK(i)=80000, IDRES(i)=mass number, IDXRES(i)=charge number
 176 *
 177 *  The 4-momenta and masses can be found in PHKK (target nucleus rest frame):
 178 *                   PHKK(1..3,i) 3-momentum (p_x,p_y,p_z)
 179 *                   PHKK(4,i)    energy
 180 *                   PHKK(5,i)    mass
 181 *
 182 *
 183 *
 184 *  Pick out the final state particles from DTEVT1 in each event for
 185 *  instance by the following loop (NHKK=number of entries in the present
 186 *  event) and fill your histograms
 187       NPART = 0
 188       DO 19 I=1,NHKK
 189 C         IF (ABS(ISTHKK(I)).NE. -1 .AND. ABS(ISTHKK(I)) .NE. 1001 ) THEN
 190          IF (ABS(ISTHKK(I)).EQ. 1) THEN
 191             NPART=NPART+1
 192          ENDIF
 193  19   CONTINUE
 194       write(16,*) npart
 195       DO 20 I=1,NHKK
 196 C         IF (ABS(ISTHKK(I)).NE. -1 .AND. ABS(ISTHKK(I)) .NE. 1001 ) THEN
 197          IF (ABS(ISTHKK(I)).EQ. 1) THEN
 198             WRITE(16,100)
 199      &           i, isthkk(i), idhkk(i),
 200      &           phkk(1,i), phkk(2,i), phkk(3,i), phkk(4,i), phkk(5,i),
 201      &           jmohkk(1,i), jmohkk(2,i), jdahkk(1,i), jdahkk(2,i)
 202          ELSEIF (ABS(ISTHKK(I)).EQ.1001) THEN
 203          ENDIF
 204  20   CONTINUE
 205  100  FORMAT(3i6, 5e13.5, 4i7)
 206 *  At any time during the run a list of the actual entries in DTEVT1 and
 207 *  DTEVT2 can be obtained (output unit 6) by the following statement:
 208 C     CALL DT_EVTOUT(4)
 209
 210       RETURN
 211 *
 212 *------------------------------------------------------------------
 213 *
 214     3 CONTINUE
 215 *
 216 * output/statistics/histograms etc.
 217 *
 218 *  Called with MODE=3 once after all events have been sampled.
 219 *
 220       RETURN
 221
 222       END