IAP GITLAB

DT_GLBINI.f 8.21 KB
Newer Older
Tanguy Pierog's avatar
Tanguy Pierog committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

      SUBROUTINE DT_GLBINI(What)
 
C***********************************************************************
C Pre-initialization of profile function                               *
C This version dated 28.11.00 is written by S. Roesler.                *
C***********************************************************************
 
      IMPLICIT NONE
      DOUBLE PRECISION debin , e , ecm , ecmini , ehi , elab , elo , 
     &                 ONE , plab , q2i , TINY14 , What , xi , ZERO
      INTEGER i , i0 , iasav , ibsav , idx , ie , ij , ijpini , ioffst , 
     &        iproj , itarg , j , jpeach , jpstep , k , KBACC , kproj , 
     &        MAXMSS , MAXOFF , nasav
      INTEGER nbsav , nebin , nlines , nproj , ntarg
      SAVE 
 
      INCLUDE 'inc/dtflka'
 
      PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY14=1.D-14)
 
      LOGICAL lcms
 
C particle properties (BAMJET index convention)
      INCLUDE 'inc/dtpart'
C properties of interacting particles
      INCLUDE 'inc/dtprta'
C emulsion treatment
      INCLUDE 'inc/dtcomp'
C Glauber formalism: flags and parameters for statistics
      INCLUDE 'inc/dtglgp'
32 33 34 35 36 37 38

#ifdef FOR_CORSIKA
cdh  datadir for path to the data sets to be read in by dpmjet/phojet
      COMMON /DATADIR/ DATADIR
      CHARACTER*132    DATADIR
#endif

Tanguy Pierog's avatar
Tanguy Pierog committed
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
C number of data sets other than protons and nuclei
C at the moment = 2 (pions and kaons)
      PARAMETER (MAXOFF=2)
      DIMENSION ijpini(5) , ioffst(25)
      DATA ijpini/13 , 15 , 0 , 0 , 0/
C Glauber data-set to be used for hadron projectiles
C (0=proton, 1=pion, 2=kaon)
      DATA (ioffst(k),k=1,25)/0 , 0 , -1 , -1 , -1 , -1 , -1 , 0 , 0 , 
     &      -1 , -1 , 2 , 1 , 1 , 2 , 2 , 0 , 0 , 2 , 0 , 0 , 0 , 1 , 
     &      2 , 2/
C Acceptance interval for target nucleus mass
      PARAMETER (KBACC=6)
 
      PARAMETER (MAXMSS=100)
      DIMENSION iasav(MAXMSS) , ibsav(MAXMSS)
      DIMENSION What(6)
 
      DATA jpeach , jpstep/18 , 5/
 
C temporary patch until fix has been implemented in phojet:
C  maximum energy for pion projectile
C       DATA ECMXPI / 100000.0D0 /
C
C--------------------------------------------------------------------------
C general initializations
C
C  steps in projectile mass number for initialization
      IF ( What(4).GT.ZERO ) jpeach = INT(What(4))
      IF ( What(5).GT.ZERO ) jpstep = INT(What(5))
C
C  energy range and binning
      elo = ABS(What(1))
      ehi = ABS(What(2))
      IF ( elo.GT.ehi ) elo = ehi
      nebin = MAX(INT(What(3)),1)
      IF ( elo.EQ.ehi ) nebin = 0
      lcms = (What(1).LT.ZERO) .OR. (What(2).LT.ZERO)
      IF ( lcms ) THEN
         ecmini = ehi
      ELSE
         ecmini = SQRT(AAM(IJProj)**2+AAM(IJTarg)**2+2.0D0*AAM(IJTarg)
     &            *ehi)
      END IF
C
C  default arguments for Glauber-routine
      xi = ZERO
      q2i = ZERO
C
C  initialize nuclear parameters, etc.
 
Cc    CALL BERTTP
Cc    CALL INCINI
 
C
C  open Glauber-data output file
      idx = INDEX(CGLb,' ')
      k = 8
      IF ( idx.GT.1 ) k = idx - 1
97
#ifndef FOR_CORSIKA
Tanguy Pierog's avatar
Tanguy Pierog committed
98
      OPEN (LDAt,FILE=CGLb(1:k)//'.glb',STATUS='UNKNOWN')
99 100 101 102 103
#else
c  modification for use with corsika using path to data file in DATADIR
      OPEN(LDAT,STATUS='UNKNOWN',
     &  FILE=DATADIR(1:INDEX(DATADIR,' ')-1)//CGLB(1:K)//'.glb')
#endif
Tanguy Pierog's avatar
Tanguy Pierog committed
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
C
C--------------------------------------------------------------------------
C Glauber-initialization for proton and nuclei projectiles
C
C  initialize phojet for proton-proton interactions
      elab = ZERO
      plab = ZERO
      CALL DT_LTINI(IJProj,IJTarg,elab,plab,ecmini,1)
      CALL DT_PHOINI
C
C  record projectile masses
      nasav = 0
      nproj = MIN(IP,jpeach)
      DO kproj = 1 , nproj
         nasav = nasav + 1
         IF ( nasav.GT.MAXMSS ) STOP ' GLBINI: NASAV > MAXMSS ! '
         iasav(nasav) = kproj
      END DO
      IF ( IP.GT.jpeach ) THEN
         nproj = DBLE(IP-jpeach)/DBLE(jpstep)
         IF ( nproj.EQ.0 ) THEN
            nasav = nasav + 1
            IF ( nasav.GT.MAXMSS ) STOP ' GLBINI: NASAV > MAXMSS ! '
            iasav(nasav) = IP
         ELSE
            DO iproj = 1 , nproj
               kproj = jpeach + iproj*jpstep
               nasav = nasav + 1
               IF ( nasav.GT.MAXMSS ) STOP ' GLBINI: NASAV > MAXMSS ! '
               iasav(nasav) = kproj
            END DO
            IF ( kproj.LT.IP ) THEN
               nasav = nasav + 1
               IF ( nasav.GT.MAXMSS ) STOP ' GLBINI: NASAV > MAXMSS ! '
               iasav(nasav) = IP
            END IF
         END IF
      END IF
C
C  record target masses
      nbsav = 0
      ntarg = 1
      IF ( NCOmpo.GT.0 ) ntarg = NCOmpo
      DO itarg = 1 , ntarg
         nbsav = nbsav + 1
         IF ( nbsav.GT.MAXMSS ) STOP ' GLBINI: NBSAV > MAXMSS ! '
         IF ( NCOmpo.GT.0 ) THEN
            ibsav(nbsav) = IEMuma(itarg)
         ELSE
            ibsav(nbsav) = IT
         END IF
      END DO
C
C  print masses
      WRITE (LDAt,99010) nebin , ': ' , SIGN(elo,What(1)) , 
     &                   SIGN(ehi,What(2))
99010 FORMAT (I4,A,1P,2E13.5)
      nlines = DBLE(nasav)/18.0D0
      IF ( nlines.GT.0 ) THEN
         DO i = 1 , nlines
            IF ( i.EQ.1 ) THEN
               WRITE (LDAt,'(I4,A,18I4)') nasav , ': ' , 
     &                (iasav(j),j=1,18)
            ELSE
               WRITE (LDAt,'(6X,18I4)') (iasav(j),j=18*i-17,18*i)
            END IF
         END DO
      END IF
      i0 = 18*nlines + 1
      IF ( i0.LE.nasav ) THEN
         IF ( i0.EQ.1 ) THEN
            WRITE (LDAt,'(I4,A,18I4)') nasav , ': ' , 
     &             (iasav(j),j=i0,nasav)
         ELSE
            WRITE (LDAt,'(6X,18I4)') (iasav(j),j=i0,nasav)
         END IF
      END IF
      nlines = DBLE(nbsav)/18.0D0
      IF ( nlines.GT.0 ) THEN
         DO i = 1 , nlines
            IF ( i.EQ.1 ) THEN
               WRITE (LDAt,'(I4,A,18I4)') nbsav , ': ' , 
     &                (ibsav(j),j=1,18)
            ELSE
               WRITE (LDAt,'(6X,18I4)') (ibsav(j),j=18*i-17,18*i)
            END IF
         END DO
      END IF
      i0 = 18*nlines + 1
      IF ( i0.LE.nbsav ) THEN
         IF ( i0.EQ.1 ) THEN
            WRITE (LDAt,'(I4,A,18I4)') nbsav , ': ' , 
     &             (ibsav(j),j=i0,nbsav)
         ELSE
            WRITE (LDAt,'(6X,18I4)') (ibsav(j),j=i0,nbsav)
         END IF
      END IF
C
C  calculate Glauber-data for each energy and mass combination
C
C   loop over energy bins
      elo = LOG10(elo)
      ehi = LOG10(ehi)
      debin = (ehi-elo)/MAX(DBLE(nebin),ONE)
      DO ie = 1 , nebin + 1
         e = elo + DBLE(ie-1)*debin
         e = 10**e
         IF ( lcms ) THEN
            e = MAX(2.0D0*AAM(IJProj)+0.1D0,e)
            ecm = e
         ELSE
            plab = ZERO
            ecm = ZERO
            e = MAX(AAM(IJProj)+0.1D0,e)
            CALL DT_LTINI(IJProj,IJTarg,e,plab,ecm,0)
         END IF
C
C   loop over projectile and target masses
         DO itarg = 1 , nbsav
            DO iproj = 1 , nasav
               CALL DT_XSGLAU(iasav(iproj),ibsav(itarg),IJProj,xi,q2i,
     &                        ecm,1,1,-1)
            END DO
         END DO
C
      END DO
C
C--------------------------------------------------------------------------
C Glauber-initialization for pion, kaon, ... projectiles
C
      DO ij = 1 , MAXOFF
C
C  initialize phojet for this interaction
         elab = ZERO
         plab = ZERO
         IJProj = ijpini(ij)
         IP = 1
         IPZ = 1
 
C* newer PHOJET versions initialize new proj/targ combinations dynamically
C* no need to call the initialization again
         CALL DT_LTINI(IJProj,IJTarg,elab,plab,ecmini,1)
C
C  calculate Glauber-data for each energy and mass combination
C
C   loop over energy bins
         DO ie = 1 , nebin + 1
            e = elo + DBLE(ie-1)*debin
            e = 10**e
            IF ( lcms ) THEN
               e = MAX(2.0D0*AAM(IJProj)+TINY14,e)
               ecm = e
            ELSE
               plab = ZERO
               ecm = ZERO
               e = MAX(AAM(IJProj)+TINY14,e)
               CALL DT_LTINI(IJProj,IJTarg,e,plab,ecm,0)
            END IF
C
C   loop over projectile and target masses
            DO itarg = 1 , nbsav
               CALL DT_XSGLAU(1,ibsav(itarg),IJProj,xi,q2i,ecm,1,1,-1)
            END DO
C
         END DO
C
      END DO
 
C--------------------------------------------------------------------------
C close output unit(s), etc.
C
      CLOSE (LDAt)
 
      END SUBROUTINE